Dirichlet Lame Indirect Method

Dirichlet Lame Indirect Method#

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *
from ngsolve import Projector, Preconditioner
from ngsolve.krylovspace import CG

Define the geometry \(\Omega \subset \mathbb R^3\) and create a mesh:

sp = Glue ( Sphere( (0,0,0), 1).faces)
mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=1)).Curve(3)
Draw( mesh )

WebGLScene

Create test and trial function finite element spaces for \(H^{-\frac12}(\Gamma)\) according to the given mesh:

fesL2 = VectorValued(SurfaceL2(mesh, order=3, dual_mapping=False))
# fesL2 = VectorH1(mesh, order=3)
print (fesL2.ndof)
u,v = fesL2.TnT()

Define Dirichlet data \(u_0\) and compute the right hand side vector:

p0 = CF( (2,2,2) )
X = CF( (x,y,z) )

E, nu = 210, 0.2
alpha = (1+nu)/((1-nu)*2*E)
norm = Norm(X-p0)
lapkernel = alpha/(4*pi) * 1/norm

u0 = (3-4*nu) * lapkernel * CF( (1,0,0) ) \
   + lapkernel/norm**2 * (X-p0) * (X-p0)[0] 
u0 *= 1000
rhs = LinearForm (u0*v.Trace()*ds(bonus_intorder=3)).Assemble()

Draw (u0, mesh)

WebGLScene

j = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()

with TaskManager():
    V = LameSL(u*ds,E,nu) *v*ds
    CG(mat = V.mat, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, maxsteps=100, initialize=False, printrates=True)

we know what we do - evaluateDeriv not implemented for dipoles in SingularMLExpansion
we know what we do - evaluateDeriv not implemented for dipoles in SingularMLExpansion
CG iteration 1, residual = 0.7606077059115475     

CG iteration 2, residual = 0.09170406274108531

CG iteration 3, residual = 0.019994607736604556

CG iteration 4, residual = 0.0073918965713169315

CG iteration 5, residual = 0.0013106864458220478

CG iteration 6, residual = 0.0006070847932383691

CG iteration 7, residual = 0.0004835676718434552

CG iteration 8, residual = 0.0006162054439700274

CG iteration 9, residual = 0.00019498405945363067

CG iteration 10, residual = 6.868587919008121e-05

CG iteration 11, residual = 5.228338761664081e-05

CG iteration 12, residual = 3.8101242313303015e-05

CG iteration 13, residual = 1.5701076228976925e-05

CG iteration 14, residual = 1.7901632405230254e-05

CG iteration 15, residual = 2.079466070865471e-05

CG iteration 16, residual = 1.3638046312363413e-05

CG iteration 17, residual = 9.09326470618203e-06

CG iteration 18, residual = 5.9109820869175995e-06

CG iteration 19, residual = 4.942680979090482e-06

CG iteration 20, residual = 4.311618581513869e-06

CG iteration 21, residual = 3.7651084141081805e-06

CG iteration 22, residual = 3.6058763228439366e-06

CG iteration 23, residual = 1.8771355068176316e-06

CG iteration 24, residual = 1.1860479641320578e-06

CG iteration 25, residual = 6.975842905598553e-07

CG iteration 26, residual = 5.957668501857826e-07

CG iteration 27, residual = 6.731271220900453e-07

CG iteration 28, residual = 4.80139355560749e-07

CG iteration 29, residual = 2.2937246516420627e-07

CG iteration 30, residual = 2.057295835168276e-07

CG iteration 31, residual = 2.1556064840081673e-07

CG iteration 32, residual = 1.3240380460051351e-07

CG iteration 33, residual = 1.6328667169896184e-07

CG iteration 34, residual = 1.106972806898702e-07

CG iteration 35, residual = 6.474813616862117e-08

CG iteration 36, residual = 5.2100626399660794e-08

CG iteration 37, residual = 2.634679032035377e-08

CG iteration 38, residual = 4.467105273836708e-08

CG iteration 39, residual = 2.5093178417819498e-08

CG iteration 40, residual = 1.8486864748934583e-08

CG iteration 41, residual = 1.7078196379884866e-08

CG iteration 42, residual = 1.040384180946171e-08

CG iteration 43, residual = 7.1036707146355344e-09

Draw (j, order=3);

vismesh = (WorkPlane().RectangleC(4,4).Face()*Sphere((0,0,0),1)).GenerateMesh(maxh=0.1).Curve(4)
sol = GridFunction(VectorH1(vismesh,order=3))

SL = LameSL(u*ds(bonus_intorder=4), E, nu)(j)
# SL = (LameSL(u*ds(bonus_intorder=4), E, nu) * v*ds) . GetPotential(j)

sol.Set (SL, definedon=vismesh.Boundaries(".*"))
Draw (sol, vismesh, order=3);
Draw (u0, vismesh, order=3);
Draw (sol-u0, vismesh, order=3);