### How to find the coordinates of the dofs of a cell given only an element type and the cell coordinates (i.e. without using a function space or a mesh)?

121
views
1
3 months ago by
Given the coordinates of the vertices of a cell and a finite element type (e.g. continuous Lagrange, degree 3) I would need to find the cell dofs coordinates without defining a mesh or a function space (I only access the vertices of only one cell at a time and I cannot/do not want to store the whole mesh).

Basically, I want to do something like:

from dolfin import *
import numpy as np

editor = MeshEditor()
mesh = Mesh()
editor.open(mesh, "triangle", 2, 2)
editor.init_vertices(3)
editor.init_cells(1)
editor.close()

V = FunctionSpace(mesh, "CG", 3)

dof_coor = V.tabulate_dof_coordinates().reshape((-1,2))

cell_dofs = dof_coor[V.dofmap().cell_dofs(0),:]​

without creating a mesh or a function space, given only the vertex coordinates (0,0), (0,1), (1,0) and the fact that I want the dofs of a P3 element. Is there an easy, efficient way to do this (apart from computing the dof coordinates myself), at least for continuous and discontinuous scalar Lagrange elements (of any degree)? Thanks a lot!

Community: FEniCS Project

6
3 months ago by
You can create the FIAT element (e.g., using ffc.fiatinterface.create_element(ufl_element)) from which you can query coordinates on reference cell using FIAT's interface. For affine transformed cells (simplices) you get physical coordinates by affine transformation (matrix multiplication). For quads/hexes the transformation is tensor product.

Alternatively one can use code generated by ffc.jit. For example

ffc.jit(ufl.FiniteElement("P", ufl.triangle, 2))​

returns two integers which are pointers to instances of ufc::finite_element and ufc::dofmap. By some ctypes trickery one should be able to get to generated tabulate_dof_coordinates which does exactly what you want and its calls should within Numba loop, for example, would be fast. But that's quite tricky to make it working. (Moreover this functionality might be completely different in future FEniCS which will probably abandon UFC. That's why first FIAT solution should be preferable in most circumstances.)
Thanks a lot Jan! Thanks to your help and by giving a look at FIAT/lagrange.py I found that the following works:

from dolfin import *
from FIAT import reference_element as re
import numpy as np

dim = 3
degree = 4
sdim = dim + 1

simplex_coor = np.random.randn(sdim, dim)

# NOTE: currently dolfin uses ufc_simplex rather than FIAT default simplex. The dofs order is different if default_simplex is used.
#ref_el = re.default_simplex(dim)
ref_el = re.ufc_simplex(dim)
top = ref_el.get_topology()

nodes = np.vstack([item for dimension in sorted(top) for entity in sorted(top[dimension]) for item in ref_el.make_points(dimension,entity,degree)])
A,b = re.make_affine_mapping(ref_el.vertices, simplex_coor)

mapped_nodes = A.dot(nodes.T).T + b
​

Here is a simple testing routine using mshr:

from dolfin import *
from mshr import generate_mesh, Rectangle, Box
from FIAT import reference_element as re
import numpy as np

for dim in range(1,4):

if dim == 1:
mesh = IntervalMesh(25, -1, 1)
elif dim == 2:
mesh = generate_mesh(Rectangle(Point(-1,-1), Point(1.,1.)), 10)
else:
mesh = generate_mesh(Box(Point(-1,-1, -1), Point(1.,1.,1.)), 10)

for degree in range(1,5):

V = FunctionSpace(mesh, "CG", degree)
dof_coor  = V.tabulate_dof_coordinates().reshape((-1,dim))

# NOTE: currently dolfin uses ufc_simplex rather than FIAT default simplex. The dofs order is different if default_simplex is used.
#ref_el = re.default_simplex(dim)
ref_el = re.ufc_simplex(dim)
top = ref_el.get_topology()

nodes = np.vstack([item for dimension in sorted(top) for entity in sorted(top[dimension]) for item in ref_el.make_points(dimension,entity,degree)])

for c in cells(mesh):
A,b = re.make_affine_mapping(ref_el.vertices, np.array(c.get_vertex_coordinates()).reshape((-1,dim)))
mapped_nodes = A.dot(nodes.T).T + b
assert np.allclose(mapped_nodes, dof_coor[V.dofmap().cell_dofs(c.index()),:])
​

This works for Lagrange elements of any degree apart from P0 (but the P0 case is straightforward).
written 3 months ago by Matteo Croci
1
I thought that you could get the nodes from instance of FIAT.FiniteElement. That works for any element.
import ufl, ffc
import itertools

e = ufl.FiniteElement("RT", ufl.triangle, 2)
fe = ffc.fiatinterface.create_element(e)
nodes = list(itertools.chain(*(iter(node.get_point_dict().keys()) for node in fe.dual_basis())))
written 3 months ago by Jan Blechta