Dirichlet boundary conditions on CR elements

1

Hello,

I am trying to use non-conforming elements for my simulation on FEniCS (still the 2019.1.0 version), the Crouzeix-Raviart element in particular. Given a simple 2D square geometry and triangular mesh, after applying a Dirichlet boundary condition with constant values on all boundaries, the result does not look reasonable. I don’t understand if it is just because of the non-conforming elements, or due to something else that I didn’t catch. CG and DG elements have also been tested, but they all work well as expected.
The first snapshot is the u_x field with CR elements, and the second one is that with CG elements. The u_y field faces the same issue as well.

Any advice or insight would be much appreciated!

from dolfin import *
# Create mesh
L = 0.5e3
H = 0.5e3
Nx = 20
Ny = 20
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed")
# Define functions
V = VectorFunctionSpace(mesh, "CR", 1)  # Function space for u
u = Function(V)                 
# Boundary conditions
boundary = CompiledSubDomain( "on_boundary")
bcb = DirichletBC(V, Constant((2, -2)), boundary, "geometric")
bcs = [bcb]
# Test BC
for bc in bcs:
    bc.apply(u.vector())
file_results = XDMFFile("test_bc.xdmf")
u.rename("u", "displacement field")
file_results.write(u)

cr-x
cr-x1541×796 36.5 KB

CG-x
CG-x1541×796 37.5 KB

2

I believe using XDMFFile.write(u) causes the function to be interpolated into a CG-1 space before it is written to file. If you project your function to a DG-1 space and use write_checkpoint, the BC appears to be applied correctly:

from dolfin import *
# Create mesh
L = 0.5e3
H = 0.5e3
Nx = 20
Ny = 20
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed")
# Define functions
V = VectorFunctionSpace(mesh, "CR", 1)  # Function space for u
u = Function(V)                 
# Boundary conditions
boundary = CompiledSubDomain( "on_boundary")
bcb = DirichletBC(V, Constant((2, -2)), boundary, "geometric")
bcs = [bcb]
# Test BC
for bc in bcs:
    bc.apply(u.vector())

file_results = XDMFFile("test_bc_me.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.parameters["rewrite_function_mesh"] = True
V_output = VectorFunctionSpace(mesh, "DG", 1)
u_output = project(u, V_output)
u_output.rename("u", "displacement field")
file_results.write(mesh)
file_results.write_checkpoint(u_output, "displacement field")
file_results.close()

screenshot
screenshot732×563 16.6 KB

2 Likes
3

Thank you! That makes perfect sense.

1 Like