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=0.5)).Curve(3)

\[V^{\text{Lame}} \, j \;=\; (\tfrac{1}{2}I + K^{\text{Lame}}) u_0\]

with $$K^{\text{Lame}}u = Ku - V (M(u)) + 2\mu V^{\text{Lame}} (M(u))$$

fesL2 = VectorValued(SurfaceL2(mesh, order=3, dual_mapping=False))
fesH1 = VectorH1(mesh, order=4)

print (fesL2.ndof)
print (fesH1.ndof)

u,v = fesL2.TnT()
uH1, vH1 = fesH1.TnT()

3360
2694

n = specialcf.normal(3)
def MM(u) : 
    g = Grad(u).Trace()
    return g.trans * n - Trace(g) * n

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

E, nu = 210, 0.2
mu = E/(2*(1+nu))
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()
u_0 = GridFunction(fesH1)
u_0.Set(u0, definedon=mesh.Boundaries(".*"))
Draw (u_0, 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
    M = BilinearForm( uH1 * v * ds(bonus_intorder=3)).Assemble()
    K = (LaplaceDL(uH1 * ds)*v*ds).mat - (LaplaceSL(MM(uH1)*ds) * v *ds).mat + 2*mu*(LameSL(MM(uH1)*ds, E, nu) * v *ds).mat    
    rhs = ( (0.5 * M.mat + K)*u_0.vec).Evaluate()
    CG(mat = V.mat, pre=pre, rhs = rhs, 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 SingularMLExpansionwe know what we do - evaluateDeriv not implemented for dipoles in SingularMLExpansion

CG iteration 1, residual = 0.05306471866672728

CG iteration 2, residual = 0.008197460079730573

CG iteration 3, residual = 0.0010597906553412024

CG iteration 4, residual = 0.00031908317255985134

CG iteration 5, residual = 0.00033143392843269976

CG iteration 6, residual = 0.0001264907107817599

CG iteration 7, residual = 7.432163205289407e-05

CG iteration 8, residual = 7.10101249538099e-05

CG iteration 9, residual = 2.9950398912906156e-05

CG iteration 10, residual = 1.0244901208426782e-05

CG iteration 11, residual = 9.073440784449057e-06

CG iteration 12, residual = 8.785461169428362e-06

CG iteration 13, residual = 4.8245544453714716e-06

CG iteration 14, residual = 2.521016765925809e-06

CG iteration 15, residual = 2.0338754710516108e-06

CG iteration 16, residual = 2.0708766058917018e-06

CG iteration 17, residual = 1.59989736297026e-06

CG iteration 18, residual = 1.1814602709432336e-06

CG iteration 19, residual = 1.687005843547034e-06

CG iteration 20, residual = 8.505522028934996e-07

CG iteration 21, residual = 5.050487782905122e-07

CG iteration 22, residual = 3.197185594622128e-07

CG iteration 23, residual = 2.2532733695494e-07

CG iteration 24, residual = 2.440796447436819e-07

CG iteration 25, residual = 1.986463312238126e-07

CG iteration 26, residual = 1.0523869769803827e-07

CG iteration 27, residual = 8.661581936683917e-08

CG iteration 28, residual = 7.831678383893163e-08

CG iteration 29, residual = 4.6610722091371704e-08

CG iteration 30, residual = 3.854842398782723e-08

CG iteration 31, residual = 2.8912765166362366e-08

CG iteration 32, residual = 3.847227901706019e-08

CG iteration 33, residual = 2.360850844329942e-08

CG iteration 34, residual = 1.3409715607542192e-08

CG iteration 35, residual = 7.754756851070537e-09

CG iteration 36, residual = 1.1396264654174743e-08

CG iteration 37, residual = 5.17109731977978e-09

CG iteration 38, residual = 5.609865885165493e-09

CG iteration 39, residual = 7.429290065281057e-09

CG iteration 40, residual = 3.2870000517268482e-09

CG iteration 41, residual = 2.0991591437649257e-09

CG iteration 42, residual = 1.6467201531983622e-09

CG iteration 43, residual = 1.3262785703190753e-09

CG iteration 44, residual = 9.539383855912127e-10

CG iteration 45, residual = 1.0254522377990573e-09

CG iteration 46, residual = 8.69488178719896e-10

CG iteration 47, residual = 7.339292827494729e-10

CG iteration 48, residual = 4.2051407792875237e-10

Draw (j, order=3);

lam = E * nu / ((1+nu) * (1-2*nu))
grad_u = CF( (u0.Diff(x), u0.Diff(y), u0.Diff(z)) )
eps_u = 1/2 * (grad_u.Reshape((3,3)) + grad_u.Reshape((3,3)).trans) 
sigma_1 = lam * (u0.Diff(x)[0] + u0.Diff(y)[1] + u0.Diff(z)[2])
sigma_1 = CF((sigma_1, 0, 0, 0, sigma_1, 0, 0, 0, sigma_1)).Reshape((3,3))
sigma = sigma_1 + 2 * mu * eps_u
uexa = sigma * n
Draw (uexa, mesh, draw_vol=False, order=3);
Draw( j - uexa, mesh, draw_vol=False, order=3)

BaseWebGuiScene

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)*v*ds).GetPotential(j)
# DL = (LaplaceDL(uH1 * ds)*v*ds).GetPotential(u_0) - (LaplaceSL(MM(uH1)*ds)*v*ds).GetPotential(u_0) + 2 * mu * (LameSL(MM(uH1)*ds, E, nu)*v*ds).GetPotential(u_0)    
DL = LaplaceDL(uH1 * ds)(u_0) - LaplaceSL(MM(uH1)*ds)(u_0) + 2 * mu * LameSL(MM(uH1)*ds, E, nu)(u_0)    
repformula = SL-DL 
sol.Set (repformula, definedon=vismesh.Boundaries(".*"))
Draw (sol, vismesh, order=3);
Draw (u0, vismesh, order=3);
Draw (sol-u0, vismesh, order=3);