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.05306471859540084
CG iteration 2, residual = 0.008197459952126156
CG iteration 3, residual = 0.0010597907586978824
CG iteration 4, residual = 0.0003190901724297698
CG iteration 5, residual = 0.00033145213921634536
CG iteration 6, residual = 0.0001264894302717878
CG iteration 7, residual = 7.432182198766852e-05
CG iteration 8, residual = 7.101076582659451e-05
CG iteration 9, residual = 2.9950381238816466e-05
CG iteration 10, residual = 1.0244924415823622e-05
CG iteration 11, residual = 9.073975131814644e-06
CG iteration 12, residual = 8.786133456908317e-06
CG iteration 13, residual = 4.824409555162587e-06
CG iteration 14, residual = 2.521046201657611e-06
CG iteration 15, residual = 2.0339107810054843e-06
CG iteration 16, residual = 2.070926896750771e-06
CG iteration 17, residual = 1.5999120848799151e-06
CG iteration 18, residual = 1.181494205614956e-06
CG iteration 19, residual = 1.6872936901142937e-06
CG iteration 20, residual = 8.505106857553488e-07
CG iteration 21, residual = 5.05049902434909e-07
CG iteration 22, residual = 3.197195810608308e-07
CG iteration 23, residual = 2.2532892973853176e-07
CG iteration 24, residual = 2.440824931909175e-07
CG iteration 25, residual = 1.987356900707483e-07
CG iteration 26, residual = 1.0522741539855747e-07
CG iteration 27, residual = 8.662174740800775e-08
CG iteration 28, residual = 7.831809325892171e-08
CG iteration 29, residual = 4.6611471526802415e-08
CG iteration 30, residual = 3.85497056248148e-08
CG iteration 31, residual = 2.8915498390028458e-08
CG iteration 32, residual = 3.851636137375147e-08
CG iteration 33, residual = 2.359930508211299e-08
CG iteration 34, residual = 1.3410529515049874e-08
CG iteration 35, residual = 7.755382870617966e-09
CG iteration 36, residual = 1.1397563896497992e-08
CG iteration 37, residual = 5.171637928264035e-09
CG iteration 38, residual = 5.6105726870542115e-09
CG iteration 39, residual = 7.4337991667889386e-09
CG iteration 40, residual = 3.286816373540478e-09
CG iteration 41, residual = 2.099357279336055e-09
CG iteration 42, residual = 1.6469140115857947e-09
CG iteration 43, residual = 1.3264481750921197e-09
CG iteration 44, residual = 9.540813832220956e-10
CG iteration 45, residual = 1.0246511606340486e-09
CG iteration 46, residual = 8.697299484048321e-10
CG iteration 47, residual = 7.343262224078423e-10
CG iteration 48, residual = 4.2057498318108603e-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'