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