Hi everyone,
I’m trying to impose a displacement on a face of a 3D mesh.
I would like to impose (0,0,100)
In order to do this, I made this:
displacement = CoefficientFunction((0, 0, 100))
gfu.Set(displacement, definedon="FaceDisplacement")
I don’t understand why it doesn’t work.
I’ve tried to use:
- a.AssembleLinearization(gfu.vec)
- inv = CGSolver(a.mat, pre, printrates=‘\r’, tol=1e-6, maxiter=500)
gfu.vec.data += inv * (f.vec - a.Apply(gfu.vec))
But nothing happen.
Do you have any Idea ?
Thanks
Thomas
My code:
from netgen.occ import *
from ngsolve import *
# Variable
E = 210000
nu = 0.3
# Dirichlet Boundaries
Face_Dirichletx = 'Face4'
Face_Dirichlety = 'Face4'
Face_Dirichletz = 'Face4'
# Displacement boundary
FaceDisplacement = 'Face2'
# Degree of Element
EleDeg = 2
# Mesh the geometry and get the boundaries --------------------------------------------------
box = Box((0,-10,-10), (200,10,10))
box.faces.Max(X).name="Face4"
box.faces.Min(X).name="Face2"
geo = OCCGeometry(box)
mesh = Mesh(geo.GenerateMesh(maxh=6)).Curve(3)
boundaries = list(mesh.GetBoundaries())
# Stress Function ----------------------------------------------------------------------------
mu = E / 2 / (1 + nu)
lam = E * nu / ((1 + nu) * (1 - 2 * nu))
def Stress(strain):
return 2 * mu * strain + lam * Trace(strain) * Id(3)
# FE Space ---------------------------------------------------------------------------------
fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichletx, dirichlety=Face_Dirichlety, dirichletz=Face_Dirichletz)
u, v = fes.TnT()
gfu = GridFunction(fes)
# Apply displacement boundary condition -----------------------------------------------------
displacement = CoefficientFunction((0, 0, 100))
gfu.Set(displacement, definedon="FaceDisplacement")
# Assemble a ----------------------------------------------------------------------------
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile() * dx)
pre = Preconditioner(a, "bddc")
a.Assemble()
# Create zero force as there is no force applied -------------------------------------------
f = LinearForm(fes)
f.Assemble()
# Run the Solver --------------------------------------------------------------------------
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-6, maxiter=500)
gfu.vec.data = inv * f.vec
# Print Result --------------------------------------------------------------------------
with TaskManager():
EFstress = MatrixValued(H1(mesh, order=EleDeg), symmetric=True)
VMstress = GridFunction(EFstress)
VMstress.Interpolate(Stress(Sym(Grad(gfu))))
Draw (gfu)
i think
gfu.Set(cf, definedon=mesh.Boundaries("bndname"))
is what you want
Thanks,
But it still don’t work, It don’t make any displacement ![:confused: :confused:](https://forum.ngsolve.org/images/emoji/twitter/confused.png?v=12)
I don’t understand why it doesn’t work.
gfu.Set(displacement, definedon=mesh.Boundaries("Face4"))
None of your faces is called “FaceDisplacement”, but “Face4”
Yes effectively, I made a mistake in the code that I sent you.
However, This still don’t work.
Have I to modifie “a” in the preconditioner ?
It’s very weird.
Thanks for your help.
My code without errors
:
from netgen.occ import *
from ngsolve import *
# Variable
E = 210000
nu = 0.3
# Dirichlet Boundaries
Face_Dirichletx = 'Face4'
Face_Dirichlety = 'Face4'
Face_Dirichletz = 'Face4'
# Displacement boundary
FaceDisplacement = 'Face2'
# Degree of Element
EleDeg = 2
# Mesh the geometry and get the boundaries --------------------------------------------------
box = Box((0,-10,-10), (200,10,10))
box.faces.Max(X).name="Face4"
box.faces.Min(X).name="Face2"
geo = OCCGeometry(box)
mesh = Mesh(geo.GenerateMesh(maxh=6)).Curve(3)
boundaries = list(mesh.GetBoundaries())
# Stress Function ----------------------------------------------------------------------------
mu = E / 2 / (1 + nu)
lam = E * nu / ((1 + nu) * (1 - 2 * nu))
def Stress(strain):
return 2 * mu * strain + lam * Trace(strain) * Id(3)
# FE Space ---------------------------------------------------------------------------------
fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichletx, dirichlety=Face_Dirichlety, dirichletz=Face_Dirichletz)
u, v = fes.TnT()
gfu = GridFunction(fes)
# Apply displacement boundary condition -----------------------------------------------------
displacement = CF((0, 0, 100))
gfu.Set(displacement, definedon=mesh.Boundaries(FaceDisplacement))
# Assemble a ----------------------------------------------------------------------------
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile() * dx)
pre = Preconditioner(a, "bddc")
a.Assemble()
# Create zero force as there is no force applied -------------------------------------------
force = CF((0, 0, 1))
f = LinearForm(fes)
# f += force * v * ds(FaceDisplacement)
f.Assemble()
# Run the Solver --------------------------------------------------------------------------
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-6, maxiter=500)
gfu.vec.data = inv * f.vec
# Print Result --------------------------------------------------------------------------
with TaskManager():
EFstress = MatrixValued(H1(mesh, order=EleDeg), symmetric=True)
VMstress = GridFunction(EFstress)
VMstress.Interpolate(Stress(Sym(Grad(gfu))))
you need to homogenize your equation. have a look at the Dirichlet itutorial
I’ve read it, but I don’t get it ![:confused: :confused:](https://forum.ngsolve.org/images/emoji/twitter/confused.png?v=12)