I am trying to bend a rectangular plate by moving the opposite edges towards each other. I am following the example here, but with different boundary conditions. When I apply boundary conditions to the angle (to clamp the edges). I get the following warning “*** FFC warning: evaluate_dof(s) for enriched element not implemented.”
If I apply only displacement boundary conditions, I can see that the boundary conditions are applied correctly. However, I am now trying to apply force boundary conditions and am running into convergence issues – and seeing this warning thousands of times makes it difficult to see what is going on or how to to fix it. My questions are as follows:
- Is this warning a symptom of issues in my code?
- If this warning is not causing any actual issues, is it possible to suppress it?
My code is below:
from dolfin import *
from fenics_shells import *
import mshr
import numpy as np
import matplotlib.pyplot as plt
parameters["form_compiler"]["quadrature_degree"] = 4
length, H = 1, 0.4
t = 2E-3
ts = Constant((t))
meshlength = int(100*length)
meshheight = int(100*H)
mesh = RectangleMesh(Point(0, 0), Point(length, H), meshlength, meshheight)
E = Constant(1.0)
nu = Constant(0.3)
kappa = Constant(5.0/6.0)
A = (E*ts/t**3/(1. - nu**2))*as_tensor([[1., nu, 0.],[nu, 1., 0.],[0., 0., (1. - nu)/2]])
D = (E*ts**3/t**3/(12.*(1. - nu**2)))*as_tensor([[1., nu, 0.],[nu, 1., 0.],[0., 0., (1. - nu)/2]])
S = E*kappa*ts/t**3/(2*(1. + nu))
P1 = FiniteElement("Lagrange", triangle, degree = 1)
P2 = FiniteElement("Lagrange", triangle, degree = 2)
bubble = FiniteElement("B", triangle, degree = 3)
enriched = P1 + bubble
element = MixedElement([VectorElement(P2, dim=2), P2, VectorElement(enriched, dim=2)])
Q = FunctionSpace(mesh, element)
q_, q, q_t = Function(Q), TrialFunction(Q), TestFunction(Q)
v_, w_, theta_ = split(q_)
e = sym(grad(v_)) + 0.5*outer(grad(w_), grad(w_))
class LeftDirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0]) < DOLFIN_EPS and on_boundary
class RightDirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - length) < DOLFIN_EPS and on_boundary
def applyforce(x, on_boundary):
return x[1] >= 0.0 and x[1] <= 0.4 and on_boundary and near(x[0],0)
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(applyforce).mark(boundary_subdomains,2)
dss = ds(subdomain_data=boundary_subdomains)
Traction = Expression(("a","0","0","0","0"),a = 0,degree = 0)
def eps(v,w):
return sym(grad(v)) + 0.5*outer(grad(w), grad(w))
ev = strain_to_voigt(e)
psi_m = 0.5*dot(A*ev, ev)
gamma = grad(w_) - theta_
psi_s = 0.5*dot(S*gamma, gamma)
k_T = as_tensor(Expression((("c/imp","0.0"),("0.0","c*imp")), c=1., imp=.999, degree=0))
k = sym(grad(theta_)) - k_T
kv = strain_to_voigt(k)
psi_b = 0.5*dot(D*kv, kv)
dx_h = dx(metadata={'quadrature_degree': 2})
h = CellDiameter(mesh)
alpha = project(t**2/h**2, FunctionSpace(mesh,'DG',0))
Pi_PSRI = psi_b*dx + alpha*psi_m*dx + alpha*psi_s*dx + (1.0 - alpha)*psi_s*dx_h + (1.0 - alpha)*psi_m*dx_h - inner(Traction, q_) * dss(2)
Pi = Pi_PSRI
dPi = derivative(Pi, q_, q_t)
J = derivative(dPi, q_, q)
#Left and right boundary conditions
u_right = Expression(("a*x[0]","a*x[1]"), a = 0 ,degree=1) #INITIAL STRAIN
u_left = Expression(("a*x[0]+i*0.0005","a*x[1]"), a = 0 ,i = 0,degree=1) #INITIAL STRAIN
bcl = DirichletBC(Q.sub(0), u_left, LeftDirichletBoundary())
bcr = DirichletBC(Q.sub(0), u_right, RightDirichletBoundary())
rightvert = DirichletBC(Q.sub(1),Constant(0),RightDirichletBoundary())
leftvert = DirichletBC(Q.sub(1),Constant(0),LeftDirichletBoundary())
bcltheta = DirichletBC(Q.sub(2), Constant((0,0)), LeftDirichletBoundary())
bcrtheta = DirichletBC(Q.sub(2), Constant((0,0)), RightDirichletBoundary())
bcs = [bcr,rightvert,leftvert,bcltheta,bcrtheta]
init = Function(Q)
q_.assign(init)
problem = NonlinearVariationalProblem(dPi, q_, bcs, J = J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-8
# c_cr = 0.02
# cs = np.linspace(0.0, c_cr, 21)
c_cr = 0.0516
cs = np.linspace(0.0, 1.5*c_cr, 30)
# and we solve as usual::
xdmf_file = XDMFFile("bending-demo.xdmf")
xdmf_file.parameters["flush_output"] = True
xdmf_file.parameters["functions_share_mesh"] = True
for count, i in enumerate(cs):
if count < 3:
k_T.c = i
else:
k_T.c = 0
# u_left.i = count
Traction.a = count*50
solver.solve()
v_h, w_h, theta_h = q_.split(deepcopy=True)
Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
# Calculate strain for plotting
strain = Function(Vsig, name="Strain")
strain.assign(project(eps(v_h,w_h),Vsig))
a = VectorFunctionSpace(mesh,'Lagrange',degree = 2)
b = FunctionSpace(mesh,'CG', 1)
disp = Function(a, name = "displacement")
disp.assign(project(v_h,a))
disp_z = Function(b, name = "displacement_z")
disp_z.assign(project(w_h,b))
xdmf_file.write(strain, i)
xdmf_file.write(disp,i)
xdmf_file.write(disp_z, i)
Here is the final shape when I apply displacement boundary conditions:
