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)*n*q) 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)
```