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

CG iteration 3, residual = 0.019994607755071177

CG iteration 4, residual = 0.007391896578715234

CG iteration 5, residual = 0.0013106864475630883

CG iteration 6, residual = 0.0006070848742120193

CG iteration 7, residual = 0.0004835732062998249

CG iteration 8, residual = 0.0006161957442112141

CG iteration 9, residual = 0.0001949839984774121

CG iteration 10, residual = 6.868587106813723e-05

CG iteration 11, residual = 5.228322496859438e-05

CG iteration 12, residual = 3.8101308998404785e-05

CG iteration 13, residual = 1.570107233335883e-05

CG iteration 14, residual = 1.7900685381357247e-05

CG iteration 15, residual = 2.0795608113198933e-05

CG iteration 16, residual = 1.3638231740302822e-05

CG iteration 17, residual = 9.093272504424768e-06

CG iteration 18, residual = 5.9109835825745736e-06

CG iteration 19, residual = 4.9426762555611975e-06

CG iteration 20, residual = 4.311585688771188e-06

CG iteration 21, residual = 3.7649233574636032e-06

CG iteration 22, residual = 3.6060191363000934e-06

CG iteration 23, residual = 1.8771388914443682e-06

CG iteration 24, residual = 1.1860463206868464e-06

CG iteration 25, residual = 6.975827388468029e-07

CG iteration 26, residual = 5.957529754290756e-07

CG iteration 27, residual = 6.73123536963093e-07

CG iteration 28, residual = 4.80143903851059e-07

CG iteration 29, residual = 2.2937129261250096e-07

CG iteration 30, residual = 2.0572819853703785e-07

CG iteration 31, residual = 2.1555916292513847e-07

CG iteration 32, residual = 1.324024475673343e-07

CG iteration 33, residual = 1.632800489286422e-07

CG iteration 34, residual = 1.1069816409915897e-07

CG iteration 35, residual = 6.474743374612029e-08

CG iteration 36, residual = 5.209998825254321e-08

CG iteration 37, residual = 2.6346343961126585e-08

CG iteration 38, residual = 4.466882377267805e-08

CG iteration 39, residual = 2.5092824633130478e-08

CG iteration 40, residual = 1.8486471476542753e-08

CG iteration 41, residual = 1.7077647238878858e-08

CG iteration 42, residual = 1.0403591921981068e-08

CG iteration 43, residual = 7.103449785905316e-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);