Imposed Displacement

NGSolve
1

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

i think

gfu.Set(cf, definedon=mesh.Boundaries("bndname"))

is what you want

3

Thanks,
But it still don’t work, It don’t make any displacement :confused:
I don’t understand why it doesn’t work.

4 
gfu.Set(displacement, definedon=mesh.Boundaries("Face4"))

None of your faces is called “FaceDisplacement”, but “Face4”

5

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 :sweat_smile::

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

you need to homogenize your equation. have a look at the Dirichlet itutorial

7

I’ve read it, but I don’t get it :confused:

8

Hi Thomas,

I am facing a similar problem I extended your code showing combinations of homogenous (fix or zero BCs) and non-homogenous Dirichlet BCs (moving BC).

BR Dieter

Dirichlet_BCs.ipynb (8.1 KB)