This page was generated from wta/adaptivity.ipynb.

Adaptivity

[1]:
from ngsolve import *
from ngsolve.webgui import Draw

Define the geometry by 2D Netgen-OpenCascade modeling (new in Netgen/NGSolve 2105):

[2]:
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo

def MakeGeometryOCC():
    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"
    DrawGeo(geo)
    return OCCGeometry(geo, dim=2)

geo = MakeGeometryOCC()

Define the geometry by curves (old-style segments geometry):

[3]:
#   point numbers 0, 1, ... 11
#   sub-domain numbers (1), (2), (3)
#
#
#             7-------------6
#             |             |
#             |     (2)     |
#             |             |
#      3------4-------------5------2
#      |                           |
#      |             11            |
#      |           /   \           |
#      |         10 (3) 9          |
#      |           \   /     (1)   |
#      |             8             |
#      |                           |
#      0---------------------------1
#

def MakeGeometry():
    from netgen.geom2d import SplineGeometry
    geometry = SplineGeometry()

    # point coordinates ...
    pnts = [ (0,0), (1,0), (1,0.6), (0,0.6), \
             (0.2,0.6), (0.8,0.6), (0.8,0.8), (0.2,0.8), \
             (0.5,0.15), (0.65,0.3), (0.5,0.45), (0.35,0.3) ]
    pnums = [geometry.AppendPoint(*p) for p in pnts]

    # start-point, end-point, boundary-condition, left-domain, right-domain:
    lines = [ (0,1,"bot",1,0), (1,2,"outer",1,0), (2,5,"outer",1,0), (5,4,"inner",1,2), (4,3,"outer",1,0), (3,0,"outer",1,0), \
              (5,6,"outer",2,0), (6,7,"outer",2,0), (7,4,"outer",2,0), \
              (8,9,"inner",3,1), (9,10,"inner",3,1), (10,11,"inner",3,1), (11,8,"inner",3,1) ]

    for p1,p2,bc,left,right in lines:
        geometry.Append(["line", pnums[p1], pnums[p2]], bc=bc, leftdomain=left, rightdomain=right)

    geometry.SetMaterial(1,"base")
    geometry.SetMaterial(2,"top")
    geometry.SetMaterial(3,"chip")

    return geometry

# geo = MakeGeometry()
# Draw(geo)

Piece-wise constant coefficients in sub-domains:

[4]:
mesh = Mesh(geo.GenerateMesh(maxh=0.2))

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

lam = CoefficientFunction([1, 1000, 10])
a = BilinearForm(fes)
a += lam*grad(u)*grad(v)*dx

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

c = Preconditioner(a, type="multigrid", inverse="sparsecholesky")

gfu = GridFunction(fes, autoupdate=True)

Assemble and solve problem:

[5]:
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\]
[6]:
l = []    # l = list of estimated total error
space_flux = HDiv(mesh, order=2, autoupdate=True)
gf_flux = GridFunction(space_flux, "flux", autoupdate=True)

def CalcError():

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

    # 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 for refinement:
    maxerr = max(eta2)
    # marking with Python loop:
    # for el in mesh.Elements():
    #    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)

    # marking using numpy:
    mesh.ngmesh.Elements2D().NumPy()["refine"] = \
        eta2.NumPy() > 0.25*maxerr

CalcError()
ndof = 208  toterr = 0.006146664512048994

Adaptive loop:

[7]:
level = 0
while fes.ndof < 50000:
    mesh.Refine()
    SolveBVP()
    CalcError()
    level = level+1
    if level%5 == 0:
        Draw (gfu)
ndof = 355  toterr = 0.005412444050627485
ndof = 610  toterr = 0.0034764703269970164
ndof = 1057  toterr = 0.0024407198787842655
ndof = 1498  toterr = 0.0017725398350735106
ndof = 2176  toterr = 0.0011753184815645348
ndof = 2977  toterr = 0.0007678238856745071
ndof = 3895  toterr = 0.0004992073011966705
ndof = 4711  toterr = 0.0003283496614114705
ndof = 5509  toterr = 0.00022040792539933297
ndof = 6271  toterr = 0.00015347867373531965
ndof = 6934  toterr = 0.00011528702327131788
ndof = 7678  toterr = 9.010822245637461e-05
ndof = 8611  toterr = 6.96416842564254e-05
ndof = 9745  toterr = 5.301219839103373e-05
ndof = 10642  toterr = 4.57838567246872e-05
ndof = 12292  toterr = 3.365421383995064e-05
ndof = 13735  toterr = 2.6279981724029656e-05
ndof = 15721  toterr = 2.0122030884399945e-05
ndof = 17944  toterr = 1.5652788759240043e-05
ndof = 20470  toterr = 1.2331216391261612e-05
ndof = 23848  toterr = 9.222303534562167e-06
ndof = 27451  toterr = 7.2979670533946725e-06
ndof = 32353  toterr = 5.5243333430289305e-06
ndof = 38077  toterr = 4.0931323780630575e-06
ndof = 43858  toterr = 3.172523131672133e-06
ndof = 51115  toterr = 2.433173856683243e-06
[8]:
Draw (gfu);
[9]:
%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();
../../_images/i-tutorials_wta_adaptivity_15_0.png