Non-conforming Crouzeix-Raviart elements on a periodic domain in parallel

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)