Hi,

as the title suggests, I would like to make sure, that the inhomogeneous boundary conditions are preserved during the adaptive mesh refinement scheme. In my current program, the first solution looks correct. However, after the error estimation and the refinement, the solution is 0 (or near 0) on the boundary. I think that the error estimation function is not correct, but I am struggling to find my blunder.

The PDE to be solved is:

[tex]-\Delta u(x) = -(18x^2-6)e^{-1.5(x^2 + y^2)}-(18y^2-6)e^{-1.5(x^2 + y^2)}[/tex]

on the domain

[tex]\Omega: x,y \in [-1,1][/tex]

subjected to the boundary condition

[tex]u(x,y)|_{\partial \Omega} = 2e^{-1.5(x^2 + y^2)}[/tex]

```
import netgen.gui
from ngsolve import *
import netgen.geom2d as geom2d
geo = geom2d.SplineGeometry()
geo.AddRectangle((-1,-1), (1,1), bcs=["bottom", "right", "top", "left"])
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
mesh.GetBoundaries()
fes = H1(mesh, order=2, dirichlet="left|right|top|bottom")
g = 2*exp(-1.5*(x**2 + y**2))
gfu = GridFunction(fes)
gfu.Set(g, definedon=mesh.Boundaries("left|right|top|bottom"))
Draw(gfu)
u, v = fes.TnT()
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
f = LinearForm(fes)
f += (-(18*x*x-6)*exp(-1.5*(x*x + y*y))-(18*y*y-6)*exp(-1.5*(x*x + y*y)))*v*dx
c = Preconditioner(a,"multigrid")
space_flux = HDiv(mesh, order=2)
gf_flux = GridFunction(space_flux, "flux")
def solveStep():
fes.Update()
gfu.Update()
Draw(gfu)
a.Assemble()
f.Assemble()
solvers.BVP(bf=a, lf=f, gf=gfu, pre=c)
gfu.Update()
Redraw(blocking=True)
def estimateError():
space_flux.Update()
gf_flux.Update()
flux = grad(gfu)
gf_flux.Set(flux)
err = (flux-gf_flux)*(flux-gf_flux)
Draw(err, mesh, 'error_representation')
eta2 = Integrate(err, mesh, VOL, element_wise=True)
maxerr = max(eta2)
for el in mesh.Elements():
mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)
solveStep()
estimateError()
mesh.Refine()
solveStep()
Draw(gfu)
```

Attachment: adaptive_and_bc.ipynb