Boundary conditions on edges or nodes

1

Hi All,

I am a novice Fenics user and am playing with some demos to get going.

I would like to run an analysis on a hyperelastic material I formulated, using a cube under hydrostatic pressure. To avoid rigid body motions, I need to apply some Dirichlet boundary conditions, say constrain a node in 3 DOFs, its neighbour in 2 and a suitably chosen third one in 1 (6DOF in total on 3 nodes to remove rigid body motions). Maybe there is a more clever way to do it in Fenics, would also be nice to know.

I am though struggling to apply BCs on nodes.
Here below I have a minimal example, essentially one hyperelasticity demo with modified BCs. It covers a cube, on whose right face a Dirichlet boundary conditions is applied, while on the left face only the edges are constrained.

from dolfin import *
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
tol = 1E-14
def right(x, on_boundary):
    return on_boundary and (near(x[0], 1, tol))
def leftplus1(x, on_boundary):
    return on_boundary and (near(x[0], 0, tol)) and (near(x[1], 1, tol)) 
def leftplus2(x, on_boundary):
    return on_boundary and (near(x[0], 0, tol)) and (near(x[1], 0, tol)) 
def leftplus3(x, on_boundary):
    return on_boundary and (near(x[0], 0, tol)) and (near(x[2], 0, tol)) 
def leftplus4(x, on_boundary):
    return on_boundary and (near(x[0], 0, tol)) and (near(x[2], 1, tol)) 
c = Constant((0.0, 0.0, 0.0))
r = Constant((1.0, 0.00, 0.0))
bcl = DirichletBC(V, c, leftplus1)
bc2 = DirichletBC(V, c, leftplus2)
bc3 = DirichletBC(V, c, leftplus3)
bc4 = DirichletBC(V, c, leftplus4)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bc2,bc3,bc4, bcr]

du = TrialFunction(V)           
v  = TestFunction(V)           
u  = Function(V)                 
B  = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T  = Constant((0.0,  0.0, 0.0))  # Traction force on the boundary
d = len(u)
I = Identity(d)            
F = I + grad(u)             
C = F.T*F                 
Ic = tr (C)
J  = det(F)
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
F = derivative(Pi, u, v)

J = derivative(F, u, du)
solve(F == 0, u, bcs, J=J)

The edge constraint is given, for example

 def leftplus1(x, on_boundary):
     return on_boundary and (near(x[0], 0, tol)) and (near(x[1], 1, tol))

So I am getting points on the boundary AND points with x[0] = 0 AND points with x[1] = 0. The intersection of two planes and the boundary should give the edge of interest, or have I totally misunderstood how to define domains?
Yet there is no displacement in the model. I checked for silly mistakes and anything trivial but could not find any, but I could not find any. ,The BC on the right side is still applied, so how come no displacements are evident? I started wondering if anything is wrong with the function definition.

Thanks a lot

2 Likes
2

Please format your code appropriately using ``` to encapsulate it.
Also, remove all unused imports and all lines that are commented out.
Make sure that the code is runnable if you do a copy-paste into an editor (that indentation is correct), and that it reproduces your error.
Follow all the instructions in: Read before posting: How do I get my question answered? - #3

3

My apologies, I certainly did not intend to leave the post in the abysmal original state, could see only the typing window and had to post it to see what state it was appearing in, sorry about that.

4

Instead of constraining the dofs, you can set the nullspace operator, see for instance: Singular Poisson — DOLFIN documentation
and to define the nullspace operator: Bitbucket

You could also have a look at: Apply a DirichletBC to a single node - #2 by klunkean

5

Thank you very much. The null space approach is elegant and I am looking into it, and will use it. I would though also be interested to understand why the simplest approach fails in my case. I checked the answer by klunkean and now implemented my boundary conditions as

def right(x, on_boundary):
    return  (near(x[0], 1, tol))
class Gammathree(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, DOLFIN_EPS) and near(x[1], 0, DOLFIN_EPS) and near(x[2], 0, DOLFIN_EPS)
class Gammatwo(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, DOLFIN_EPS) and near(x[1], 0, DOLFIN_EPS) and near(x[2], 1, DOLFIN_EPS)
class Gammaone(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, DOLFIN_EPS) and near(x[1], 1, DOLFIN_EPS) and near(x[2], 0, DOLFIN_EPS)
point3 = Gammathree()
point2 = Gammatwo()
point1 = Gammaone()

# Define Dirichlet boundary (x = 0 or x = 1)
c = Constant((0.0, 0.0, 0.0))
r = Constant((1.0, 0.00, 0.0))
bc1= DirichletBC(V, c, point3)
bc21 = DirichletBC(V.sub(0), Constant(0.0), point2)
bc22 = DirichletBC(V.sub(1), Constant(0.0), point2)
bc3 = DirichletBC(V.sub(0), Constant(0.0), point1)
bcr = DirichletBC(V, r, right)
bcs = [bc1,bc21,bc22,bc3,bcr]

Point 1 is the point at (0,0,0), constrained in all 3 DOFs. Point 2 is the point at (0,0,1), constrained in the x- and y- directions (hence I use, if I understood, V.sub(0) and V.sub(1)).
Point 3 is the point at (0,1,0), constrained in the x-direction only.
Then I apply a translation in the x-direction to the right face.
I do not get any errors but there is no displacement in the results.

Thanks a lot

6

First of all, you are missing one thing, which is that in DirichletBC in the other post, the argument pointwise is used.

Secondly,
You can inspect your Dirichlet conditions with the following:

bc1= DirichletBC(V, c, point3, method="pointwise")
bc21 = DirichletBC(V.sub(0), Constant(0.0), point2, method="pointwise")
bc22 = DirichletBC(V.sub(1), Constant(0.0), point2, method="pointwise")
bc3 = DirichletBC(V.sub(0), Constant(0.0), point1, method="pointwise")
bcr = DirichletBC(V, r, right)
u_ = Function(V)
u_.vector()[:] = 0.5

bcs = [bc1,bc21,bc22,bc3,bcr]
[bc.apply(u_.vector()) for bc in bcs]
File("u.pvd") <<u_

and inspect u.pvd in for instance paraview (I suggest using the point visualization there).
If you use the modified BCs, I get a result by using for instance a 4x4x4 cube.
Displaced mesh is showed below:

elasticity_fixed
elasticity_fixed1280Ă—720 114 KB

To me, it seems you have not gotten your elasticity parameters right.

7

Hi there, thank you very much for your reply, indeed I was careless and forgot the “pointwise” method. It now does what is was supposed to do, although the test is admittedly rather artificial.
I tried to understand what is that your snippet does, I searched and experimented but could not make sense of it. Running my script, with the correct boundary conditions, I get the following result

bc_pointwise
bc_pointwise2250Ă—874 70.6 KB
, which is what I expect.
Running your snippet you added yo your last post, I get this (I visualised it using points, but it is not that intellegible on the post, so I use “surfaces” here)

u_def
u_def2250Ă—874 16.8 KB

So, what is does your snippet does?


bcs = [bc1,bc21,bc22,bc3,bcr]
[bc.apply(u_.vector()) for bc in bcs]
File("u.pvd") <<u_```

what is the 0.5 standing for? and what is the purpose of the iteration over the bc ?
I understand it must be something obvious, sorry I am totally missing it and the fact I get that funny sheared shape puzzles me. 

Thanks a lot really much appreciated
8

Setting it to 0.5 is just so that it is possible to identify zero boundary conditions. Simply set all dof values to some value that is not in any of the boundary conditions. They apply function applies the bc to each dof. Thus looking at it in paraview can identify which dofs that there has been applied to.

9

Sorry I am being really thick here.
You saying that you are displacing all the DOFs, except the constrained one, by 0.5?
Then I cannot make sense of your plot. apart from the 3 constrained nodes, the deformation should be affine, which does not seem to be. Moreover the max magnitude in the plot equals 1.0, not 0.5. Also I do not understand why my plot is different.
Thanks so much for your patience…

10

So the figured I showed in my post is not u_, but the solution to your artificial problem.

u_ in my case would be perturbed 0.5 in each direction for nodes except those that have a Dirichlet BC applied to them.

11

I see, then it all makes perfect sense, thanks a lot