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 polygons, 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()

Generate the mesh, the faces of the geometry give the materials of the mesh:

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

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

Piece-wise constant coefficients in sub-domains:

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("chip"))

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

gfu = GridFunction(fes)

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.006210901701303183

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.005434323839293221
ndof = 583  toterr = 0.0034474161165093395
ndof = 997  toterr = 0.002477627244137173
ndof = 1492  toterr = 0.001764188010920619
ndof = 2158  toterr = 0.0011758942561252656

ndof = 2950  toterr = 0.0007720684985710706
ndof = 3826  toterr = 0.0005050543890735538
ndof = 4651  toterr = 0.00033457395108488805
ndof = 5422  toterr = 0.00022627521645150598
ndof = 6121  toterr = 0.00015855835364824745

ndof = 6586  toterr = 0.00012787180844763262
ndof = 7321  toterr = 9.64400518726508e-05
ndof = 8200  toterr = 7.79102661322021e-05
ndof = 9184  toterr = 6.082706314040289e-05
ndof = 9895  toterr = 5.020406565163159e-05

ndof = 11209  toterr = 3.913525053079808e-05
ndof = 12994  toterr = 2.805370426289392e-05
ndof = 14365  toterr = 2.327350691219258e-05
ndof = 16252  toterr = 1.8486313859438926e-05
ndof = 18613  toterr = 1.4268462247050254e-05

ndof = 21763  toterr = 1.0776534037956513e-05
ndof = 25099  toterr = 8.414666975753977e-06
ndof = 28927  toterr = 6.561001565597816e-06
ndof = 33445  toterr = 5.037684819305905e-06
ndof = 38899  toterr = 3.8443944840868665e-06

ndof = 45985  toterr = 2.9011683599063498e-06
ndof = 53179  toterr = 2.277809617359276e-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();