Hello everyone,

I’m trying to solve a problem where I need to numerically verify **Green’s formula** in 2D for piecewise linear functions. The formula I’m trying to implement is as follows:

Where:

- ( \Gamma = \partial \Omega )
- ( \nu ) is the outer unit normal vector
- ( a ) is a constant 2D vector

For ease, I define the following integrals:

The exact formula should satisfy:

My task is to implement this in **Python** using **NGSolve** and evaluate the error:

And verify that the error is on the order of machine precision (( \approx 10^{-15} )).

**Code Attempted:**

I’m new to both **NGSolve** and **Python**, and although I’ve tried several approaches, I haven’t been able to solve the problem. Below is the code I’m working with:

```
# Importing necessary libraries from NGSolve
import ngsolve as ng
from ngsolve import Mesh, H1, GridFunction, CoefficientFunction, grad, InnerProduct
from netgen.geom2d import unit_square
import numpy as np
# Create mesh for the unit square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
# Define the finite element space (H1 space with piecewise linear functions)
fes = H1(mesh, order=1)
# Trial and test functions
u = fes.TrialFunction()
v = fes.TestFunction()
# Define the constant vector 'a'
a = ng.CoefficientFunction((1, 1)) # Example: a = (1, 1)
# Define some functions for u and v in H1 (choose different examples to test)
u_exact = ng.sin(ng.pi * ng.x) * ng.sin(ng.pi * ng.y)
v_exact = ng.cos(ng.pi * ng.x) * ng.cos(ng.pi * ng.y)
# Create GridFunction for u and v in the finite element space
u_gf = GridFunction(fes)
v_gf = GridFunction(fes)
u_gf.Set(u_exact)
v_gf.Set(v_exact)
# Calculate the integrals I1, I2, I3
# I1 = Integral of (a · ∇v) u dx over Omega
I1 = ng.Integrate(InnerProduct(a, grad(v)) * u_gf, mesh)
# I2 = Integral of (ν · a) u v dS over the boundary Gamma
n = ng.specialcf.normal(2) # normal vector
I2 = ng.Integrate(InnerProduct(n, a) * u_gf * v_gf, mesh.BoundaryCF())
# I3 = Integral of (a · ∇u) v dx over Omega
I3 = ng.Integrate(InnerProduct(a, grad(u_gf)) * v_gf, mesh)
# Error defined as: ERR = I1 - I2 + I3
ERR = I1 - I2 + I3
# Output the results
print(f"I1: {I1}")
print(f"I2: {I2}")
print(f"I3: {I3}")
print(f"Error (ERR): {ERR}")
# Check if error is close to machine precision (~10^-15)
assert np.abs(ERR) < 1e-14, "Error is not within machine precision!"
```

**Problem:**

Despite following this approach, I can’t verify that the error is on the order of ( 10^{-15} ), and the code doesn’t seem to generate the expected results. Could anyone help me figure out what I might be doing wrong or suggest improvements to the code?

Thank you in advance for your help!