Exact DOFs for RT2 spaces

[cuiyu-he1]

Hello! I have a question for the exact definition of DOFs in RT2 space (with 8 DOFs in each element in 2D). My understanding is that the basis functions on the edges (2D) should have zero moments in the element, i.e., if multiplying one edge basis by any test function in VectorFunctionSpace(mesh,‘DG’,0), its integration should be 0. But in fact, that is not true based on my computation. I put a short code as follows. the Any help is highly appreciated! Thanks.
P.S. I am using the older version 1.5.0 of Fenics.

from fenics import *
mesh = UnitSquareMesh(2,2,‘crossed’)
mesh.init()
RT2 = FunctionSpace(mesh,‘RT’,2)
dofmap = RT2.dofmap()

tau = Function(RT2)
tau.vector()[:] = 0.0

%Set all boundary DOF = 1.0 and all interior DOF = 0.0.
for cell in cells(mesh):
cell_dof = dofmap.cell_dofs(cell.index())
for l1, l2 in enumerate(cell.entities(1)):
for index in dofmap.tabulate_facet_dofs(l1):
tau.vector()[index] = 1.0

DG0 = VectorFunctionSpace(mesh,‘DG’,0)
w = TestFunction(DG0)
L = inner(tau,w)*dx
b = assemble(L).array().reshape(-1,2)
print(b)

---

%I should get b[:]==0. But I get the following:
array([[ 0.0625 , -0.02083333],
[ 0.01041667, -0.03125 ],
[ 0.01041667, 0.03125 ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ]])

[mscroggs]

Hello.

In FIAT, the RT2 basis functions are defined using point evaluations rather than integral moments, so these integrals will not be 0. You can however create the integral moment version by passing the variant "integral" fo FIAT when creating the element (although I can remember how to do this from dolfin).

In Basix (the replacement of FIAT in FEniCSx), the integral moment version of RT is used by default.

[dokken]

To add to @mscroggs answer, here is a minimal example using dolfinx:

import dolfinx
from mpi4py import MPI
import ufl
import numpy as np
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD,2,2)
RT2 = dolfinx.FunctionSpace(mesh,("RT",2))
def boundary(x):
    return np.full(x.shape[1], True)
        
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, boundary)
boundary_dofs = dolfinx.fem.locate_dofs_topological(RT2, mesh.topology.dim-1, boundary_facets)

tau = dolfinx.Function(RT2)
with tau.vector.localForm() as loc:
    loc[boundary_dofs] = np.ones(len(boundary_dofs))

DG0 = dolfinx.VectorFunctionSpace(mesh,("DG",0))
w = ufl.TestFunction(DG0)
L = ufl.inner(tau,w)*ufl.dx
b = dolfinx.fem.assemble_vector(L).array.reshape(-1,2)
print(b)

which yields:

[[-1.16226473e-16 -7.63278329e-17]
 [ 7.63278329e-17  1.16226473e-16]
 [ 3.46944695e-17  6.93889390e-17]
 [ 0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00]
 [-6.93889390e-17 -3.46944695e-17]
 [ 1.45716772e-16  1.45716772e-16]
 [-1.45716772e-16 -1.45716772e-16]]

[cuiyu-he1]

Dear Jorgen,

Thanks for the clarification. The code looks nice in FenicsX! It seems FenicsX indeed has changed a lot and in the better directions from my point of view. When Fenics will produce a stable version for users?

Unfortunately, I have to use Fenics 1.5.0 for my CutFEM method at this moment. For one my my project, I need to locally construct a flux in RT2 space. It is fine that boundary DOFs do not zero moments in the interior element. Could you then kindly tell me what exactly is the interior DOFs for RT_k (k >1) elements?

[mscroggs]

The DOFs of an RT2 space in FIAT are:

On edges

0. Point evaluation at (2/3, 1/3) in outward normal direction
1. Point evaluation at (1/3, 2/3) in outward normal direction
2. Point evaluation at (0, 1/3) in inward normal direction
3. Point evaluation at (0, 2/3) in inward normal direction
4. Point evaluation at (1/3, 0) in outward normal direction
5. Point evaluation at (1/3, 0) in outward normal direction

On interior

6. Point evaluation at (1/3, 1/3) in direction (1, 0)
7. Point evaluation at (1/3, 1/3) in direction (0, 1)

In general, the DOFs of in order k space will have k normal point evaluations on each edge and a point evaluation in the two directions at a triangular array of k(k-1)/2 points in the interior.

[dokken]

@cuiyu-he1

Thanks for the clarification. The code looks nice in FenicsX! It seems FenicsX indeed has changed a lot and in the better directions from my point of view. When Fenics will produce a stable version for users?

We are working towards that and I hope in a month or so there could be a stable release.

[cuiyu-he1]

Thanks, Matthew.

However, the DOFs in the interior is probably not correct.

from fenics import *
import numpy as np
mesh = UnitSquareMesh(2,2,'crossed')
mesh.init()
RT2 = FunctionSpace(mesh,'RT',2)
dofmap = RT2.dofmap()
v = dofmap.tabulate_all_coordinates(mesh).reshape(-1,2)
tau = Function(RT2)
tau.vector()[:] = 0.0
for cell in cells(mesh):
    cell_dof = dofmap.cell_dofs(cell.index())
    tau.vector()[cell_dof[6]]= 1.0
    tau.vector()[cell_dof[7]]= 0.0
    value = tau(v[cell_dof[6]])
    #I should expect the first value =1 and second value = 0, right?
    print('on cell', cell.index(),'1st part=',np.inner(value,[1,0]),'2nd-part=',np.inner(value,[0,1]))

Below is what I got:

# ('on cell', 0, '1st part=', 4.0, '2nd-part=', 2.2204460492503091e-15)

# ('on cell', 1, '1st part=', -2.2204460492503091e-15, '2nd-part=', -4.0)

# ('on cell', 2, '1st part=', -2.2204460492503225e-15, '2nd-part=', 3.9999999999999876)

# ('on cell', 3, '1st part=', -3.9999999999999876, '2nd-part=', 2.2204460492503225e-15)

# ('on cell', 4, '1st part=', 4.0, '2nd-part=', 7.1054273576009892e-15)

# ('on cell', 5, '1st part=', -2.2204460492502945e-15, '2nd-part=', -4.0000000000000124)

# ('on cell', 6, '1st part=', -2.2204460492503387e-15, '2nd-part=', 3.9999999999999734)

# ('on cell', 7, '1st part=', -3.9999999999999876, '2nd-part=', 7.1054273576010318e-15)

# ('on cell', 8, '1st part=', 4.0000000000000124, '2nd-part=', 2.2204460492502945e-15)

# ('on cell', 9, '1st part=', -7.1054273576009892e-15, '2nd-part=', -4.0)

# ('on cell', 10, '1st part=', -7.1054273576010318e-15, '2nd-part=', 3.9999999999999876)

# ('on cell', 11, '1st part=', -3.9999999999999734, '2nd-part=', 2.2204460492503387e-15)

# ('on cell', 12, '1st part=', 4.0000000000000124, '2nd-part=', 7.1054273576009419e-15)

# ('on cell', 13, '1st part=', -7.1054273576009419e-15, '2nd-part=', -4.0000000000000124)

# ('on cell', 14, '1st part=', -7.1054273576010839e-15, '2nd-part=', 3.9999999999999734)

# ('on cell', 15, '1st part=', -3.9999999999999734, '2nd-part=', 7.1054273576010839e-15)

I don’t understand what I have missed. Thanks again for your help!

[mscroggs]

Looks like I missed a detail when parsing the FIAT source code: the DOFs are scaled by the area of the triangle, which explains why you’re getting fours, not ones.

The DOFs are in the directions (1, 0) and (0, 1) on the reference element. The RT elements use a contravariant Piola mapping to map from the reference to each cell. This changes the directions each DOF is associated with. (The definitions of these mapping are written down here.)

[cuiyu-he1]

Thanks a lot for the clarification and timely reply, Matthew. The DOFs for RT spaces are indeed more complex than I thought in Fenics. My original difficulty was to manually set up DOFs for RT2 spaces. Unfortunately, I am still not quite sure how to set up them. I will post a separate question with more precise information. Many thanks again!