Help with Numerical Verification of Green’s Formula (NGSolve + Python)

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:

\int_{\Omega} (a \cdot \nabla v) u \, dx = \int_{\Gamma} (\nu \cdot a) u v \, dS(x) - \int_{\Omega} (a \cdot \nabla u) v \, dx, \quad \forall u, v \in H^1(\Omega)

Where:

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

For ease, I define the following integrals:

I_1 := \int_{\Omega} (a \cdot \nabla v) u \, dx, \quad I_2 := \int_{\Gamma} (\nu \cdot a) u v \, dS(x), \quad I_3 := \int_{\Omega} (a \cdot \nabla u) v \, dx

The exact formula should satisfy:

I_1 = I_2 - I_3

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

\text{ERR} := I_1 - I_2 + I_3

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!

Integrating over boundaries:

I2 = ng.Integrate(InnerProduct(n, a) * u_gf * v_gf, mesh.Boundaries(".*"))

in one of the volume terms you are using the v - Proxyfunction instead of the v_gf

best

Dear Christopher,

Thank you for the tip. I’ve corrected what you pointed out, but now I’m encountering the following error:

"I2 = ng.Integrate(InnerProduct(n, a) * u_gf * v_gf, mesh.BoundaryCF())

TypeError: BoundaryCF(): incompatible function arguments. The following argument types are supported:
1. (self: ngsolve.comp.Mesh, values: dict, default: ngsolve.fem.CoefficientFunction = None) → ngsolve.fem.CoefficientFunction"

I’m really sorry for asking these questions, but how can I solve this? Could you possibly show me an example?

Best regards.

This:

:slight_smile:

Thanks christopher, solved my problem!

I’m glad to see someone from my FEM class posting questions here. This was their homework problem! :slight_smile:

Note: I posted an example template code to help with this.

-Shawn