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)

BaseWebGuiScene

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

CG iteration 2, residual = 0.09170406274182927

CG iteration 3, residual = 0.019994607741467343

CG iteration 4, residual = 0.0073918965728583304

CG iteration 5, residual = 0.0013106864505248172

CG iteration 6, residual = 0.0006070848446565741

CG iteration 7, residual = 0.0004835705170618677

CG iteration 8, residual = 0.0006162005006621329

CG iteration 9, residual = 0.00019498403216602585

CG iteration 10, residual = 6.868588026067958e-05

CG iteration 11, residual = 5.228330941636588e-05

CG iteration 12, residual = 3.810127662244888e-05

CG iteration 13, residual = 1.5701068722123848e-05

CG iteration 14, residual = 1.7900792143821442e-05

CG iteration 15, residual = 2.0795480541679236e-05

CG iteration 16, residual = 1.3638195284053943e-05

CG iteration 17, residual = 9.093264948745973e-06

CG iteration 18, residual = 5.9109813696019045e-06

CG iteration 19, residual = 4.942678129960301e-06

CG iteration 20, residual = 4.311594700980456e-06

CG iteration 21, residual = 3.7649389100084064e-06

CG iteration 22, residual = 3.606011023214073e-06

CG iteration 23, residual = 1.8771402509141382e-06

CG iteration 24, residual = 1.1860484252292992e-06

CG iteration 25, residual = 6.975844671117768e-07

CG iteration 26, residual = 5.957540675122877e-07

CG iteration 27, residual = 6.73125965820208e-07

CG iteration 28, residual = 4.801465775714914e-07

CG iteration 29, residual = 2.2937274358873808e-07

CG iteration 30, residual = 2.0572710920199737e-07

CG iteration 31, residual = 2.1556065646630723e-07

CG iteration 32, residual = 1.32404116502714e-07

CG iteration 33, residual = 1.6327995563678024e-07

CG iteration 34, residual = 1.1069993109915118e-07

CG iteration 35, residual = 6.474870374843146e-08

CG iteration 36, residual = 5.210061973081227e-08

CG iteration 37, residual = 2.6346812909115976e-08

CG iteration 38, residual = 4.4668879836913424e-08

CG iteration 39, residual = 2.509239627055493e-08

CG iteration 40, residual = 1.8487400114362422e-08

CG iteration 41, residual = 1.70780372085372e-08

CG iteration 42, residual = 1.040395132777244e-08

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