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)

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

CG iteration 3, residual = 0.019994607755119093

CG iteration 4, residual = 0.007391896578447634

CG iteration 5, residual = 0.0013106864471977356

CG iteration 6, residual = 0.000607084872809374

CG iteration 7, residual = 0.0004835731512155814

CG iteration 8, residual = 0.0006161958381162534

CG iteration 9, residual = 0.0001949839988106078

CG iteration 10, residual = 6.868587106207078e-05

CG iteration 11, residual = 5.228322539075288e-05

CG iteration 12, residual = 3.810130883406934e-05

CG iteration 13, residual = 1.5701072152926592e-05

CG iteration 14, residual = 1.7900665184733963e-05

CG iteration 15, residual = 2.0795628210325495e-05

CG iteration 16, residual = 1.3638235015270745e-05

CG iteration 17, residual = 9.093272436834613e-06

CG iteration 18, residual = 5.910983503932787e-06

CG iteration 19, residual = 4.9426761185771955e-06

CG iteration 20, residual = 4.311585053903089e-06

CG iteration 21, residual = 3.7649189193884797e-06

CG iteration 22, residual = 3.6060223914659646e-06

CG iteration 23, residual = 1.8771388866073025e-06

CG iteration 24, residual = 1.1860462287064805e-06

CG iteration 25, residual = 6.975826650395629e-07

CG iteration 26, residual = 5.957526099424508e-07

CG iteration 27, residual = 6.731234161539171e-07

CG iteration 28, residual = 4.801439918264029e-07

CG iteration 29, residual = 2.2937125063302526e-07

CG iteration 30, residual = 2.0572814104877823e-07

CG iteration 31, residual = 2.155591018192311e-07

CG iteration 32, residual = 1.3240240105476483e-07

CG iteration 33, residual = 1.63279880324102e-07

CG iteration 34, residual = 1.1069815574337509e-07

CG iteration 35, residual = 6.47474020696459e-08

CG iteration 36, residual = 5.209995818246054e-08

CG iteration 37, residual = 2.6346325012779743e-08

CG iteration 38, residual = 4.466877185314126e-08

CG iteration 39, residual = 2.5092796222595818e-08

CG iteration 40, residual = 1.848645866116395e-08

CG iteration 41, residual = 1.707762908718787e-08

CG iteration 42, residual = 1.0403581899779125e-08

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