This page was generated from wta/adaptivity.ipynb.

Adaptivity

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

Define the geometry by 2D Netgen-OpenCascade modeling: Define rectangles and polygones, and glue them together to one shape object. OCC maintains the full geometry topology of vertices, edges, faces and solids.

def MakeGeometry():
    base = Rectangle(1, 0.6).Face()
    chip = MoveTo(0.5,0.15).Line(0.15,0.15).Line(-0.15,0.15).Line(-0.15,-0.15).Close().Face()
    top = MoveTo(0.2,0.6).Rectangle(0.6,0.2).Face()
    base -= chip

    base.faces.name="base"
    chip.faces.name="chip"
    chip.faces.col=(1,0,0)
    top.faces.name="top"
    geo = Glue([base,chip,top])
    geo.edges.name="default"
    geo.edges.Min(Y).name="bot"
    Draw(geo)
    return geo

geo = MakeGeometry()

Piece-wise constant coefficients in sub-domains:

mesh = geo.GenerateMesh(maxh=0.2, dim=2)
print (mesh.GetMaterials())

('base', 'chip', 'top')

fes = H1(mesh, order=3, dirichlet="bot")
u, v = fes.TnT()

lam = mesh.MaterialCF( { "base" : 1, "chip" : 1000, "top" : 20 } )
a = BilinearForm(lam*grad(u)*grad(v)*dx)

# heat-source in inner subdomain
f = LinearForm(1*v*dx(definedon="chip"))

c = preconditioners.MultiGrid(a, inverse="sparsecholesky")

gfu = GridFunction(fes)

WARNING: kwarg 'inverse' is an undocumented flags option for class <class 'ngsolve.comp.MultiGridPreconditioner'>, maybe there is a typo?

Assemble and solve problem:

def SolveBVP():
    a.Assemble()
    f.Assemble()
    inv = CGSolver(a.mat, c.mat)
    gfu.vec.data = inv * f.vec

SolveBVP()
Draw (gfu, mesh);

Gradient recovery error estimator: Interpolate finite element flux

\[q_h := I_h (\lambda \nabla u_h)\]

and take difference as element error indicator:

\[\eta_T := \tfrac{1}{\lambda} \| q_h - \lambda \nabla u_h \|_{L_2(T)}^2\]

l = []    # list of dofs and estimated total error
gf_flux = GridFunction(HDiv(mesh, order=2))

def CalcError():
    flux = lam * grad(gfu)   # the FEM-flux
    gf_flux.Set(flux)        # interpolate into H(div)

    # compute estimator:
    err = 1/lam*(flux-gf_flux)*(flux-gf_flux)
    eta2 = Integrate(err, mesh, VOL, element_wise=True)
    l.append ((fes.ndof, sqrt(sum(eta2))))
    print("ndof =", fes.ndof, " toterr =", sqrt(sum(eta2)))

    # mark elements for refinement using numpy vectorization
    mesh.ngmesh.Elements2D().NumPy()["refine"] = eta2.NumPy() > 0.25*max(eta2)

CalcError()

ndof = 208  toterr = 0.006210901699713995

Adaptive loop:

level = 0
while fes.ndof < 50000:
    mesh.Refine()
    SolveBVP()
    CalcError()
    level = level+1
    if level%5 == 0:
        Draw (gfu)

ndof = 328  toterr = 0.005434323836822083
ndof = 583  toterr = 0.003447416114778153
ndof = 997  toterr = 0.002477627243723903
ndof = 1492  toterr = 0.0017641880110845322
ndof = 2158  toterr = 0.0011758942559793512

ndof = 2950  toterr = 0.0007720684984829111
ndof = 3826  toterr = 0.0005050543890373647
ndof = 4651  toterr = 0.00033457395107638915
ndof = 5422  toterr = 0.00022627521645824218
ndof = 6121  toterr = 0.00015855835366346445

ndof = 6586  toterr = 0.00012787180844394836
ndof = 7321  toterr = 9.644005188341731e-05
ndof = 8200  toterr = 7.791026613586427e-05
ndof = 9184  toterr = 6.082706313476574e-05
ndof = 9895  toterr = 5.0204065645794995e-05

ndof = 11209  toterr = 3.913525052794873e-05
ndof = 12994  toterr = 2.8053704250526714e-05
ndof = 14365  toterr = 2.3273506899091517e-05
ndof = 16252  toterr = 1.8486313851777513e-05
ndof = 18613  toterr = 1.4268462235108226e-05

ndof = 21763  toterr = 1.077653403096017e-05
ndof = 25099  toterr = 8.414666972635488e-06
ndof = 28927  toterr = 6.5610015619228555e-06
ndof = 33445  toterr = 5.037684816781611e-06
ndof = 38899  toterr = 3.8443944828715225e-06

ndof = 45985  toterr = 2.9011683578169164e-06
ndof = 53179  toterr = 2.2778096154866847e-06

Draw (gfu);

Plotting convergence:

%matplotlib inline
import matplotlib.pyplot as plt
plt.yscale('log')
plt.xscale('log')
plt.xlabel("ndof")
plt.ylabel("H1 error-estimate")
ndof,err = zip(*l)
plt.plot(ndof,err, "-*")

plt.ion()
plt.show();