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.053064718562771604
CG iteration 2, residual = 0.008197459938925584
CG iteration 3, residual = 0.001059790755147785
CG iteration 4, residual = 0.00031908945107483503
CG iteration 5, residual = 0.0003314502538995826
CG iteration 6, residual = 0.00012648958237663536
CG iteration 7, residual = 7.432180667730537e-05
CG iteration 8, residual = 7.101062535294532e-05
CG iteration 9, residual = 2.995040025633676e-05
CG iteration 10, residual = 1.0244926294702609e-05
CG iteration 11, residual = 9.073944831220187e-06
CG iteration 12, residual = 8.786094194486355e-06
CG iteration 13, residual = 4.824430160896618e-06
CG iteration 14, residual = 2.521056797991032e-06
CG iteration 15, residual = 2.0339274689865908e-06
CG iteration 16, residual = 2.0709553713592717e-06
CG iteration 17, residual = 1.599929653405409e-06
CG iteration 18, residual = 1.1814795338186915e-06
CG iteration 19, residual = 1.6870073190569184e-06
CG iteration 20, residual = 8.505760157445821e-07
CG iteration 21, residual = 5.050607861306973e-07
CG iteration 22, residual = 3.1972756542118314e-07
CG iteration 23, residual = 2.2533378701115862e-07
CG iteration 24, residual = 2.4409446136331244e-07
CG iteration 25, residual = 1.9883895346575447e-07
CG iteration 26, residual = 1.0521840971210262e-07
CG iteration 27, residual = 8.664970218396085e-08
CG iteration 28, residual = 7.832077789844356e-08
CG iteration 29, residual = 4.6612060647979265e-08
CG iteration 30, residual = 3.8553465487636226e-08
CG iteration 31, residual = 2.891534459519822e-08
CG iteration 32, residual = 3.842567012584932e-08
CG iteration 33, residual = 2.3626551691786042e-08
CG iteration 34, residual = 1.3413057095215874e-08
CG iteration 35, residual = 7.757017837721011e-09
CG iteration 36, residual = 1.1400883883291996e-08
CG iteration 37, residual = 5.1732430514244025e-09
CG iteration 38, residual = 5.610948350576802e-09
CG iteration 39, residual = 7.432106898759783e-09
CG iteration 40, residual = 3.288223606603221e-09
CG iteration 41, residual = 2.100047887912177e-09
CG iteration 42, residual = 1.647558243968189e-09
CG iteration 43, residual = 1.326994674264469e-09
CG iteration 44, residual = 9.545597600628558e-10
CG iteration 45, residual = 1.0299636401805815e-09
CG iteration 46, residual = 8.698290982577107e-10
CG iteration 47, residual = 7.330897248422344e-10
CG iteration 48, residual = 4.20744243787917e-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);
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[9], line 4
1 vismesh = (WorkPlane().RectangleC(4,4).Face()*Sphere((0,0,0),1)).GenerateMesh(maxh=0.1).Curve(4)
2 sol = GridFunction(VectorH1(vismesh,order=3))
3
----> 4 SL = (LameSL(u*ds(bonus_intorder=4), E, nu)*v*ds).GetPotential(j)
5 # 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)
6 DL = LaplaceDL(uH1 * ds)(u_0) - LaplaceSL(MM(uH1)*ds)(u_0) + 2 * mu * LameSL(MM(uH1)*ds, E, nu)(u_0)
7 repformula = SL-DL
AttributeError: 'ngsolve.bem.IntegralOperator' object has no attribute 'GetPotential'