Hello,
I want to apply the following boundary condition on a 2D face : x[0]=0 and no constrain on x[1] and x[2].
I wanna do it on this hyperelasticity example, on the right face.
If i do :
r1=Expression(("0","x[1]","x[2]"),degree=2)
it doesn’t work.
Here is my code :
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
E=1e6
nu=0.3
#optimisation de la compilation
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
# import mesh
mesh = BoxMesh(Point(0,0,0),Point(1,1,0.1),20,20,2)
# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
class PeriodicBoundary(SubDomain):
def __init__(self, tolerance=DOLFIN_EPS):
""" vertices stores the coordinates of the 4 unit cell corners"""
SubDomain.__init__(self, tolerance)
self.tol = tolerance
def inside(self, x, on_boundary):
# return True if on bottom boundary
return bool((near(x[1], 0, self.tol)) and on_boundary)
def map(self, x, y):
y[0] = x[0]
y[1] = x[1] - 1
y[2] = x[2]
PBC = PeriodicBoundary(tolerance=1e-8)
V = VectorFunctionSpace(mesh, "Lagrange", 1,constrained_domain=PBC)
v = TestFunction(V)
du = TrialFunction(V)
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, 0.0, -100000.0)) # Body force per unit volume
T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
r0= Constant((0.,0.,0.))
r1= Constant((0.,0.,0.))
bc_L = DirichletBC(V,r0,left)
bc_R = DirichletBC(V,r1,right)
bcs =[bc_L,bc_R]
# Kinematics
d = len(u)
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Lame parameters
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(F, u, du)
# Solve variational problem
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['maximum_iterations'] = 20
prm['newton_solver']['relaxation_parameter'] = 1.0
solver.solve()
# Save solution in VTK format
file = File("displacement.pvd");
file << u;
Thanks