Extra dofs selected when adjacent faces are combined in boundary conditions with 3D BDM (N2F) elements

dolfinx
1

Extra facets are inadvertently selected when combining faces in locate_entities_boundary (indepedent of element type), and subsequently extra dofs in locate_dofs_topological (on e.g. 3D BDM (N2F) elements, but not e.g. 3D P elements). Separating them out into seperate boundary conditions results in consistency with dolfin.

from ufl import *
from dolfinx import *
import numpy as np
from mpi4py import MPI
from slepc4py import SLEPc
mesh=UnitCubeMesh(MPI.COMM_WORLD,1,1,1)
BDM=FiniteElement("BDM",mesh.ufl_cell(),1)
W=FunctionSpace(mesh,BDM)
def noslip_boundary(x): return np.isclose(x[0],0.0)
facdim=2
facets=dolfinx.mesh.locate_entities_boundary(mesh,facdim,noslip_boundary)
print(facets)
dofs=dolfinx.fem.locate_dofs_topological(W,facdim,facets)
print(dofs)
bc0=DirichletBC(Function(W),dofs)
def noslip_boundary(x): return np.isclose(x[1],1.0)
facets=dolfinx.mesh.locate_entities_boundary(mesh,facdim,noslip_boundary)
print(facets)
dofs=dolfinx.fem.locate_dofs_topological(W,facdim,facets)
print(dofs)
bc1=DirichletBC(Function(W),dofs)
def noslip_boundary(x): return np.logical_or(np.isclose(x[0],0.0),np.isclose(x[1],1.0))
facets=dolfinx.mesh.locate_entities_boundary(mesh,facdim,noslip_boundary)
print(facets)
dofs=dolfinx.fem.locate_dofs_topological(W,facdim,facets)
print(dofs)
bc=DirichletBC(Function(W),dofs)
us=TrialFunction(W)
vs=TestFunction(W)
a=inner(us,vs)*dx
A = dolfinx.fem.assemble_matrix(a,[bc])
A.assemble()
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)
eigensolver.solve()
vr, vi = A.createVecs()
l = eigensolver.getEigenpair(0 ,vr, vi)
print(l.real)
A = dolfinx.fem.assemble_matrix(a,[bc0,bc1])
A.assemble()
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)
eigensolver.solve()
vr, vi = A.createVecs()
l = eigensolver.getEigenpair(0 ,vr, vi)
print(l.real)

Results (just printing facets and dofs)

[8 4]
[[ 0]
[ 1]
[ 2]
[21]
[22]
[23]]
[15 14]
[[27]
[28]
[29]
[42]
[43]
[44]]
[15 8 17 4 3 14] <- facets 17 and 3 are also inadvertently selected
[[ 0]
[ 1]
[ 2]
[ 9] <- extra dof selected
[10] <- extra dof selected
[11] <- extra dof selected
[21]
[22]
[23]
[27] <- extra dof selected
[28] <- extra dof selected
[29] <- extra dof selected
[37]
[38]
[39]
[42]
[43]
[44]]
1.7746788776946878 (eigenvalue trivial system doing everything in one go w.r.t. boundary condtion bc)
1.8034741465722868 (eigenvalue trivial system separating out boundary conditions bc0 and bc1)

(This does not happen on all adjacent faces)

Just following the Stokes equations with Taylor-Hood elements demo here with respect to the boundary conditions.
Versioning: latest kali linux with apt get fenics fenicx + latest git ffcx and dolfinx (3 days ago)

2

How about using locate_dofs_geometrical?

Note that encapsulating code with ``` helps on formatting and readability

3

I have not tried alternative techniques. But I think this should be fixed and looked into, as it is the same methodology as in a reference demo. It is also not obvious to the user. The workaround is in the post: separate out the faces into separate boundary conditions.

4

This is the exact reason for having both locate_dofs_geometrical and locate_dofs_topological.
The function you use when you use locate_entities_boundary will check if all vertices of a given facet satisfies the function that you are sending in. Which is does for more facets when you take the union of these two sets. Locate dofs geometrical checks the dof coordinates of each dof, and will return the correct number of dofs for your approach:

from ufl import *
from dolfinx import *
from mpi4py import MPI
import numpy as np

mesh=UnitCubeMesh(MPI.COMM_WORLD,1,1,1)
BDM=FiniteElement("BDM",mesh.ufl_cell(),1)
W=FunctionSpace(mesh,BDM)


def noslip_boundary(x):
    return np.isclose(x[0],0.0)

dofs0 = dolfinx.fem.locate_dofs_geometrical(W, noslip_boundary)
print(len(dofs0), dofs0.T)

def noslip_boundary(x):
    return np.isclose(x[1],1.0)
dofs1 = dolfinx.fem.locate_dofs_geometrical(W, noslip_boundary)
print(len(dofs1), dofs1.T)


def noslip_boundary(x):
    return np.logical_or(np.isclose(x[0],0.0),np.isclose(x[1],1.0))
dofs = dolfinx.fem.locate_dofs_geometrical(W, noslip_boundary)
print(len(dofs), dofs.T)
print(np.isin(dofs, dofs0).T, np.isin(dofs, dofs1).T)

which returns:

6 [[ 0  1  2 21 22 23]]
6 [[27 28 29 42 43 44]]
12 [[ 0  1  2 21 22 23 27 28 29 42 43 44]]
[[ True  True  True  True  True  True False False False False False False]] [[False False False False False False  True  True  True  True  True  True]]
5

Thank you for your alternative, and more direct method. But to come back on what goes wrong originally: Preliminary answer: Apparently, only the vertices are checked. So then also internal facets are chosen. Why doesn’t there appear to be a problem in 2D but there is in 3D? Due to the name of the function, locate_entities_boundary should check that facets are only connected to a single element. (moreover, it is not immediately clear how locate_dof_geometric works properly with mixed dofs, as it appears to return pairs of dofs)