Hi everybody,
I’m trying to solve a linear elasticity problem (see MWE below) with Crouzeix-Raviart non-conforming elements (note the stabilization term penalizing the jump of the displacement across faces of two adjacent elements) on a periodic domain. Everything works well in serial but I obtain the error below when I run it in parallel. Everything works well also (in serial and parallel) when I use conforming Lagrange elements, or when I use the Crouzeix-Raviart non-conforming elements for a non-periodic problem. Would you have any idea of what is going on and how to fix that? I am using FEniCS 2018 on Docker Toolbox.
Thank you very much.
Victor
Error message:
Traceback (most recent call last):
File "elasticity_periodic.py", line 71, in <module>
solve(lhs(F)==rhs(F), u, bc1)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 247, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** [fenics-support@googlegroups.com](mailto:fenics-support@googlegroups.com)
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** Process: 1
***
*** DOLFIN version: 2018.1.0
*** Git changeset: 948dc42cc4e06ed9227d0201ad50f94ac94cbf9f
*** -------------------------------------------------------------------------
MWE:
#!/usr/bin/env python
# coding: utf-
from __future__ import print_function
from dolfin import *
import numpy as np
comm = MPI.comm_world
comm_rank = MPI.rank(comm)
set_log_level(30)
parameters["ghost_mode"] = "shared_facet"
mesh = UnitSquareMesh(20,20)
def corner(x):
return ((near(x[0],0.)) and (near(x[1],0.)))
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
left=bool(near(x[0],0.))
bottom=bool(near(x[1],0.))
right=bool(near(x[0],1.))
top=bool(near(x[1],1.))
return bool((left or bottom) and (not ((left and top) or (bottom and right))) and on_boundary)
def map(self, x, y):
if near(x[0], 1.) and near(x[1], 1.):
y[0] = x[0] - 1.
y[1] = x[1] - 1.
elif near(x[0], 1.):
y[0] = x[0] - 1.
y[1] = x[1]
elif near(x[1], 1.):
y[0] = x[0]
y[1] = x[1] - 1.
else:
y[0] = -1.
y[1] = -1.
V = VectorFunctionSpace(mesh, "CR", 1, constrained_domain=PeriodicBoundary())
u = Function(V)
v = TrialFunction(V)
dv = TestFunction(V)
mu=1
lmbda=5
f = Expression(("4*(l + 2*m)*pi*pi*cos(2*pi*x[0])","4*(l + 2*m)*pi*pi*sin(2*pi*x[1])"), l=lmbda, m=mu, degree=5)
bc1 = DirichletBC(V, Expression(("1","0"),degree=1), corner, method="pointwise")
def eps(v):
return sym(grad(v))
F = (lmbda*tr(eps(v))*tr(eps(dv))+2*mu*inner(eps(v),eps(dv)))*dx - inner(f, dv)*dx
h=FacetArea(mesh)
h_avg = (h('+') + h('-'))/2
F += (1.0/h_avg)*dot(jump(v),jump(dv))*dS
solve(lhs(F)==rhs(F), u, bc1)