gradient of finite element function inaccurately integrated

Main Forum Old-Forum - do not post here
1

Hello,

I am using NGsolve to solve Poisson problem over a square domain using pure Neumann boundary conditions with mixed H1 and number space. On the right hand side I have a term (grad(p_exact)nq) on the boundary of the domain (p is my unknown). I define p_exact to be a finite element function to be able to take the gradient. When I solve my problem using 1st order elements I get convergence of order 1 (while I need to get order 2). Then I have realized when I change the grad(p_exact) with its real value I actually get convergence of order 2. Below I attached a MWE, I am not sure if there is a bug in NGsolve.

 #############################
import time
from ngsolve.internal import *
from ngsolve import *
from xfem import *
from xfem.lsetcurv import *
import numpy as np
from netgen.geom2d import unit_square, SplineGeometry
import math

# Geometry approximation order
geo_order = 1
order=1
nu = 1.0
alpha = 10*order**2
geo = SplineGeometry()
# Add channel
pi = math.pi
geo.AddRectangle(p1=(-1,-1), p2=(1,1), bcs=("wall", "outlet", "wall", "inlet"))
Lx=0
maxhh=(1/4)*(2**(-Lx))
mesh = Mesh(geo.GenerateMesh(maxh=maxhh))
mesh.Curve(geo_order)
Lt=0
deltatt=(1/16)*(2**(-Lt))



def solve(mesh):

  p_ex = CoefficientFunction(x*x)
  f = CoefficientFunction((-4))
  Q = H1(mesh, order=order, dgjumps = True)

  
  R = FESpace("number", mesh) 
  QR = FESpace([Q, R], dgjumps = True)
  p, r =  QR.TrialFunction()
  q, s =  QR.TestFunction()
  p1 = GridFunction(Q)
  p1_ = GridFunction(QR)
  
  Q_ex = H1(mesh, order=4)
  p_ex_h = GridFunction(Q_ex) 
  p_ex_h.Set(p_ex)
  
  
  p_exact = GridFunction(Q)
  p_exact.Set(p_ex)
   
  # Mesh size + normal
  h = specialcf.mesh_size
  n = specialcf.normal(2)
    
  ah_p = BilinearForm(QR)
  ah_p += SymbolicBFI(nu*InnerProduct(grad(p),grad(q)))
  ah_p += SymbolicBFI(p*s+q*r)  
                                 
  # RHS for step1 
  lh_p = LinearForm(QR)
  lh_p += SymbolicLFI((f*q))
  lh_p += SymbolicLFI(InnerProduct(p_ex, s))
  #using lines below I get 1st order 
  lh_p += SymbolicLFI(InnerProduct(grad(p_exact)*n,q), VOL_or_BND=BND, skeleton=True, definedon=mesh.Boundaries("inlet|outlet|wall"))
  #switching to line below I get 2nd order (I discard the wall from the boundary integral since it needs to give me a 0 p_analutical is x^2.)
  #lh_p += SymbolicLFI((2*q), VOL_or_BND=BND, skeleton=True, definedon=mesh.Boundaries("inlet|outlet")) 
                                                                                                                                
  ah_p.Assemble()
  p_mat_inv  = ah_p.mat.Inverse(QR.FreeDofs())

  lh_p.Assemble() 
  p1_.vec.data = p_mat_inv*lh_p.vec
  p1.Set(p1_.components[0])

  l2_error_p = sqrt(Integrate(InnerProduct(p1-p_ex_h , p1-p_ex_h), mesh))
  print ("refinement level = %d   L2-error (p) L_2norm = %e" % (ref_level, l2_error_p))

  Redraw(blocking=True)
  l2_error_p = 0
# Solve problem
ref_level_max =4
ref_level = 0

for ref_level in range(ref_level_max+1):    
    if ref_level > 0:
        mesh.Curve(geo_order)
        mesh.Refine()
    solve(mesh)
2

Hello,

you define “p_exact” to be a first order GridFunction, thus the gradient “grad(p_exact)” is a piece-wise constant discontinuous function. Your right hand side then does not have enough regularity to obtain the highest possible order.

  1. Use CoefficientFunctions as long as possible since you loose accuracy when you interpolate your function with a GridFunction (using the Set() command).

You could define the gradient ad CoefficientFunction:

gradp = CoefficientFunction((2*x,0)) lh_p += SymbolicLFI(InnerProduct(gradp*n,q), VOL_or_BND=BND, skeleton=True, definedon=mesh.Boundaries("inlet|outlet|wall"))
This gives me second order convergence.

  1. Use another space with higher order (I tried an p=2 H1 space) and calculate the gradient as you did.

  2. You could also use sympy to calculate the gradient and then convert those expressions to a NGSolve CoefficientFunction. I just used it for scalar function, but you can manage the dimension yourself (see the attached file).

Best regards,
Christoph

Attachment: gradient.py

3

Dear Christoph,

I would like to make sure I am computing my errors in the correct way. You wrote in the previous message to use coefficient functions as long as possible because interpolation causes a loss in accuracy. The attached file solves the Poisson equation with mixed boundary conditions. The finite element space has linear basis functions:

fes = H1(mesh, order = 1, dirichlet = 'left|bottom|right')

I have defined the coefficient function uex = (1 - y) * sin(pi * x)
and interpolated it:

gfuex = GridFunction(fes) 
gfuex.Set(uex)

I get quite different values when I compute

sqrt(Integrate((gfu - uex) ** 2, mesh)
sqrt(Integrate((gfu - gfuex) ** 2, mesh)

I know that the default order for the Integrate method is 5. Should I set this value to align with the order of the finite element space?

Thank you,

Barry

Attachment: primal.py