Question about numbering of local DoFs

dolfinx
1

When I was trying out some examples, I found that there seems to be some difference in the numbering of DoFs between basix and dolfinx.

If I tabulate the DoFs of a quadrilateral cell using dolfinx routine like this:

domain = mesh.create_rectangle(comm, [[0.0, 0.0], [1.0, 1.0]], [1, 1], mesh.CellType.quadrilateral, dtype=default_scalar_type,  ghost_mode=mesh.GhostMode.shared_facet)
V = fem.functionspace(domain, ("Lagrange", 1))
print(V.dofmap.list)
print(V.tabulate_dof_coordinates())

The results will be:

[[0 1 2 3]]
[[ 2.77555756e-17  5.55111512e-17  0.00000000e+00]
 [-2.77555756e-17  1.00000000e+00  0.00000000e+00]
 [ 1.00000000e+00 -2.77555756e-17  0.00000000e+00]
 [ 1.00000000e+00  1.00000000e+00  0.00000000e+00]]

As shown above, the second DoF is assigned to the vertex on the y-axis.

However, if I tabulate the DoFs (associated basis functions) of a reference cell in basix like this:

lagrange = basix.create_element(
    basix.ElementFamily.P, CellType.quadrilateral, 1, basix.LagrangeVariant.equispaced
)
points = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
tab = lagrange.tabulate(0, points)[0, :, :, 0]
print(tab)

The results will be:

[[ 1.00000000e+00  0.00000000e+00 -2.77555756e-17  5.55111512e-17]
 [ 1.38777878e-16  1.00000000e+00  2.77555756e-17 -5.55111512e-17]
 [ 1.38777878e-16  2.77555756e-17  1.00000000e+00 -5.55111512e-17]
 [-8.32667268e-17  0.00000000e+00  0.00000000e+00  1.00000000e+00]]

As shown above, the second basis funtion (or the second DoF) is the one located at the vertex on the x-axis. And this is exactly the numbering rule defined in DefElement.

I’m wondering how does this difference happen? Is it related to the “orientation” of the cell entities?

I care about this detail because I’m working with a custom assembler. I need to replace the element stiffness matrix with a pre-calculated numpy array, so I have to figure out what is the local numbering rule of dofs in a cell.

2

You can see that the “dofmap” coordinates is logical from the point of view of the mesh geometry`:

from mpi4py import MPI
from dolfinx import mesh, default_scalar_type, fem

comm = MPI.COMM_WORLD

domain = mesh.create_rectangle(comm, [[0.0, 0.0], [1.0, 1.0]], [1, 1], mesh.CellType.quadrilateral, dtype=default_scalar_type,  ghost_mode=mesh.GhostMode.shared_facet)
print(domain.geometry.x)
V = fem.functionspace(domain, ("Lagrange", 1))
print(V.dofmap.list)
print(V.tabulate_dof_coordinates())

returning the same as the dofmap:

[[0. 0. 0.]
 [0. 1. 0.]
 [1. 0. 0.]
 [1. 1. 0.]]
[[0 1 2 3]]
[[ 5.55111512e-17  5.55111512e-17  0.00000000e+00]
 [-5.55111512e-17  1.00000000e+00  0.00000000e+00]
 [ 1.00000000e+00 -5.55111512e-17  0.00000000e+00]
 [ 1.00000000e+00  1.00000000e+00  0.00000000e+00]]

Could I ask why you need to know the orientation of the cell to do a custom assembly?
The local dof ordering in the assembly kernel follows the basix ordering, you just cannot assume that an element is traversed as (x_0, y_0), (x_1, y_0), (x_0, y_1), (x_1, y_1), as a general mesh has no concept of what the logical first direction should be.

Some other notes that might be interesting is

and that in general the mesh sent in is sent through several re-ordering algorithms to ensure data-locality. Of course, in a one cell mesh, there isn’t much data locality going on.

Another note is that you could create the 1-cell mesh you want with:

from mpi4py import MPI
from dolfinx import mesh, fem
import numpy as np
import basix.ufl
comm = MPI.COMM_WORLD

points = np.array([[0,0],[1,0],[0,1],[1,1]], dtype=np.float64)
cells = np.array([[0,1,2,3]], dtype=np.int64)
cmap = basix.ufl.element("Lagrange", "quadrilateral", 1, shape=(points.shape[1],))
domain = mesh.create_mesh(MPI.COMM_WORLD, cells,points,  cmap)
print(domain.geometry.x)
V = fem.functionspace(domain, ("Lagrange", 1))


print(V.dofmap.list)
print(V.tabulate_dof_coordinates())
from dolfinx import io
with io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)

but as said above there is probably a better way of assembling this local matrix of your than assuming something out the coordinate directions.

3

Thank you very much for your detailed reply! Things I want to do are described below.

I can somehow calculate the local stiffness matrix A_local on a quadrilateral reference cell. And this special cell has some ‘‘orientation’’, for example, the DoFs on the left edge is different from the DoFs on the right edge. Now I have to study a mesh consisting of this special quadilateral cells oriented in the same direction (In fact, this mesh is “virtual” and I just want to use the parallel data structures generated by dolfinx, like dofmap, index map, etc.).

What I want to figure out is, if I build this mesh with mesh.create_rectangle(), will the cells be oriented in the same direction as well? For example, for all of the cells, the first reference DoF will be on the bottom left corner, and the second reference DoF will be on the top left corner, and all that.

Meanwhile, your last code snippet inspires me. Does it means that, if I want to build a mesh with all the cells oriented the same way, I can manually set the point sets?

I hope I have made myself clear. Thank you again!

4

You should not do that (I.e. building a hand oriented mesh).

My question is:

  1. As long as you know that there is a local dof layout, why does it matter if the actual cell is rotated in physical space.

For custom assembly there are some illustrative examples in section 8.3 of: DOLFINx: The next generation FEniCS problem solving environment

The important thing for a quadrilateral cell is that the first and last dof of the local ordering is in opposite corners (i.e. that the dofs are neither clockwise or counterclockwise)

5

Because the ‘‘cell’’ actually originates from a reduced order modeling process, so it has some physical meanings and ‘‘directions’’. The ‘‘mesh’’ consists of repeated cells and they are all oriented in the same direction in my case, or, their directions are arbitrary but pre-assigned. A related study is: EIFEM.

So for a regular quadrilateral mesh, I hope to make sure things go like (an example): for dofmap[0] the first entry is the global numbering for the left-bottom corner, the second entry is the global numbering for the left-top corner, and so on. And the same is for dofmap[1], dofmap[2], .... There should not be a dofmap[i] whose first entry is for the left-bottom corner but the second entry is for the right-bottom corner.

This of course destroys the data-locality to some extent, so do I have to realize the underlying data structrue by myself, or is there any way I can realize it in dolfinx? Maybe with some efficiency trade-off.

6

You can figure out the rotation of any regularly ordered input cell, by looking at the cell coordinates that are input to your custom kernel, and rotate your custom kernel likewise.

See fig 12 of the dolfinx paper.

7

Oops, I have ignored such a simple solution. Thanks a lot!

Besides, I have made a script to test the orientation of the regular meshes generated by routines like dolfinx.mesh.create_box(). The results show that the cells in such regular meshes are indeed oriented to the same direction, no matter parallel or not.

from dolfinx import fem, mesh, default_scalar_type
from mpi4py import MPI
comm = MPI.COMM_WORLD

def isOriented(coord):
    if (np.allclose(coord[1]-coord[0], [1., 0., 0.]) and 
        np.allclose(coord[2]-coord[0], [0., 1., 0.]) and 
        np.allclose(coord[3]-coord[0], [1., 1., 0.]) and
        np.allclose(coord[4]-coord[0], [0., 0., 1.]) and
        np.allclose(coord[5]-coord[0], [1., 0., 1.]) and
        np.allclose(coord[6]-coord[0], [0., 1., 1.]) and
        np.allclose(coord[7]-coord[0], [1., 1., 1.])):
        return True
    else:
        return False

nx, ny, nz = [10, 10, 10]
domain = mesh.create_box(comm, [[0.0, 0.0, 0.0], [nx, ny, nz]], [nx, ny, nz], mesh.CellType.hexahedron, dtype=default_scalar_type,  ghost_mode=mesh.GhostMode.shared_facet)
V = fem.functionspace(domain, ("Lagrange", 1))
coord = V.tabulate_dof_coordinates()
num_owned_cells = domain.topology.index_map(domain.topology.dim).size_local
num_ghost_cells = domain.topology.index_map(domain.topology.dim).num_ghosts
local_dof_map = V.dofmap.list

print('On rank', comm.rank)
print('number of owned cells', num_owned_cells)
print('number of ghost cells', num_ghost_cells)
print('are owned cells oriented', np.array([isOriented(coord[local_dof_map[cell]]) for cell in range(num_owned_cells)]).all())
print('are ghost cells oriented', np.array([isOriented(coord[local_dof_map[cell]]) for cell in range(num_owned_cells, num_owned_cells+num_ghost_cells)]).all())