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.0062027218384066805
Adaptive loop:
[7]:
level = 0
while fes.ndof < 50000:
mesh.Refine()
SolveBVP()
CalcError()
level = level+1
if level%5 == 0:
Draw (gfu)
ndof = 328 toterr = 0.005424419894719641
ndof = 583 toterr = 0.0034370983010913294
ndof = 997 toterr = 0.0024617055250623224
ndof = 1492 toterr = 0.0017417208970402825
ndof = 2143 toterr = 0.001159036769265472
ndof = 2917 toterr = 0.0007578273431000435
ndof = 3775 toterr = 0.000493096916783832
ndof = 4651 toterr = 0.0003247245805364616
ndof = 5440 toterr = 0.00021839607818093112
ndof = 6211 toterr = 0.00015253103584489278
ndof = 6910 toterr = 0.00011329794969631048
ndof = 7456 toterr = 9.041817442322814e-05
ndof = 8263 toterr = 7.05016031884768e-05
ndof = 9433 toterr = 5.622462330206224e-05
ndof = 10495 toterr = 4.587268141504037e-05
ndof = 12130 toterr = 3.351750230659081e-05
ndof = 13675 toterr = 2.564899368696983e-05
ndof = 15259 toterr = 2.077431963572974e-05
ndof = 18055 toterr = 1.5345303474890484e-05
ndof = 20347 toterr = 1.2390591822142242e-05
ndof = 23761 toterr = 9.324033982234093e-06
ndof = 27262 toterr = 7.144658685394859e-06
ndof = 31426 toterr = 5.654098735383225e-06
ndof = 37588 toterr = 4.121876481044875e-06
ndof = 43645 toterr = 3.159814850590886e-06
ndof = 50794 toterr = 2.4145846823594603e-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();
