Dirichlet Lame Direct Method

Dirichlet Lame Direct 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)

WebGLScene

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 SingularMLExpansion

CG iteration 1, residual = 0.05306471871114103

CG iteration 2, residual = 0.008197460078677587

CG iteration 3, residual = 0.0010597906742765516

CG iteration 4, residual = 0.00031908470188886265

CG iteration 5, residual = 0.0003314379110644896

CG iteration 6, residual = 0.00012649044358746224

CG iteration 7, residual = 7.432166148329243e-05

CG iteration 8, residual = 7.101012461880533e-05

CG iteration 9, residual = 2.9950407362877855e-05

CG iteration 10, residual = 1.0244906942084607e-05

CG iteration 11, residual = 9.073574509932192e-06

CG iteration 12, residual = 8.785632013794162e-06

CG iteration 13, residual = 4.824527275208457e-06

CG iteration 14, residual = 2.521035600289753e-06

CG iteration 15, residual = 2.0339005839093778e-06

CG iteration 16, residual = 2.070902132041345e-06

CG iteration 17, residual = 1.599917588196883e-06

CG iteration 18, residual = 1.1814604651451862e-06

CG iteration 19, residual = 1.68687143732224e-06

CG iteration 20, residual = 8.505874310177408e-07

CG iteration 21, residual = 5.050553309259229e-07

CG iteration 22, residual = 3.197233900687517e-07

CG iteration 23, residual = 2.2533146540075478e-07

CG iteration 24, residual = 2.440886874766224e-07

CG iteration 25, residual = 1.9882195254340797e-07

CG iteration 26, residual = 1.0521713216017165e-07

CG iteration 27, residual = 8.661971447450749e-08

CG iteration 28, residual = 7.832094504023174e-08

CG iteration 29, residual = 4.6614024040877463e-08

CG iteration 30, residual = 3.8551873306784166e-08

CG iteration 31, residual = 2.8913493015738247e-08

CG iteration 32, residual = 3.840772323720404e-08

CG iteration 33, residual = 2.3629210302854758e-08

CG iteration 34, residual = 1.3412158739454453e-08

CG iteration 35, residual = 7.75637290363118e-09

CG iteration 36, residual = 1.1398701773654828e-08

CG iteration 37, residual = 5.172314583272166e-09

CG iteration 38, residual = 5.611277531290539e-09

CG iteration 39, residual = 7.430531151615057e-09

CG iteration 40, residual = 3.2880595669911775e-09

CG iteration 41, residual = 2.099849931988543e-09

CG iteration 42, residual = 1.6473646434533858e-09

CG iteration 43, residual = 1.3268244526784796e-09

CG iteration 44, residual = 9.544017850296361e-10

CG iteration 45, residual = 1.0301835094517962e-09

CG iteration 46, residual = 8.694423520075224e-10

CG iteration 47, residual = 7.329716599174737e-10

CG iteration 48, residual = 4.206884257356553e-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)

WebGLScene

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);

called base class apply, type = N6ngsbem30DifferentialOperatorWithFactorE
called base class apply, type = N6ngsbem30DifferentialOperatorWithFactorE
called base class apply, type = N6ngsbem30DifferentialOperatorWithFactorE