Triangle orientation preserved on load?

Question

Hello,

I have an oriented triangulated surface mesh in $\mathbb{R}^3$ given as a Fenics XML file.
The orientation is given in terms of the ordering of the vertices for each <triangle> tag. In particular, if a triangle is given by
$t=[v_0, v_1, v_2]$ then there is an induced normal direction by this ordering, say, we always take
$n = \mathrm{normalize}(e_1 \times e_2)$ where $e_1 := p_1 - p_0$ and $e_2 := p_2 - p_0$, and $p_i$ are the geometric points corresponding to vertex indices $v_i$.
My orientation of the triangles is consistent in the sense that the normals computed this way all point either "upwards" or "downwards", but don't flip.

I now try to mimic this in Fenics by deriving from Expression:

class MySurfaceNormal(Expression):
    def __init__
        ...
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self._mesh, ufc_cell.index)
        n = cell.cell_normal()
        values[0] = n[0]
        values[1] = n[1]
        values[2] = n[2]

However the resulting normals are not consistently oriented - some point upwards, some downwards, though the normal space they span looks pretty correct. My guess is that this has to do with Dolfin's internal (local?) indexing scheme for vertex coordinates. A quick look at dolfin/io/XMLMesh::read_mesh suggests that at least the import preserves the order of the indices in each cell, so I probably haven't understood completely the local indexing or connectivity vector in MeshTopology, or mixed up something else (maybe with DoF indices?).

Q: How can I get the normals for the surface mesh that come from the initial vertex ordering in each triangle, following the choice of cross product mentioned above?

I am aware of Marie Rognes' nice article on solving PDEs on the Moebius strip and sphere, but there she's using a global normal field to define the orientation of each cell, using init_cell_orientations. Still, I should be able to deduce an orientation by the ordering of the triangles.

Thanks!

Answer 1

Meanwhile I've found out what's giving me a headache: Fenics/Dolfin uses a renumbering scheme of mesh entities. In particular after the mesh is loaded (and constructed via MeshEditor), the MeshEditor::close() triggers this reordering of vertex indices that appear in in each edge and triangle, currently using a lexicographic ordering as prescribed by UFC (see chapter 16.5 in the Fenics book). This is also touched in the answer to
this question posted here.

As a workaround I've created a table which tells me if my original orientation is preserved by lexicographic ordering or not, and multiply the normal with -1 in the latter case, which I have to pass to the MySurfaceNormal expression.
However this sounds a little clumsy and in my opinion it would be easy to keep track of this orientation within the Fenics Mesh class. In particular for surfaces in 3d this would come in handy if the raw mesh data is already oriented.
Is this a reasonable issue for a feature request? Should be fairly easy to implement.

Comments

Hi Konstantin,

Currently I have a similar problem: I need to keep the initial direction of the normals for a mesh read from an XML file (obtained through Gmsh -> dolfin-convert). I am therefore very interested by the method you suggest to use for normal/face orientation housekeeping.
As I am just starting to learn Python as well as Fenics, I was wondering if you could provide some piece of code to do that.
Thanks a lot for any help!

---

For any triangle $t = [i_0, i_1, i_2]$ you basically just have to check if the index ordering is equivalent to lexicographic ordering, which means that it is an even permutation of the lexicographic ordering. For triangles these are just cyclic shifts of the indices. A simple python function could look like this:

def index_parity(t):
    if t[0] < t[1] < t[2] or t[1] < t[2] < t[0] or t[2] < t[0] < t[1]:
        return 1
    else:
        return -1

Then just multiply the normal with either $1$ or $-1$, depending on the ordering.

You still have to read the original XML though, since the Dolfin import destroys the original ordering. To do so just use any XML parser lib you like, for instance etree.ElementTree from Python's xml module. You can get the triangles from the XML with triangles = tree.findall('mesh/cells/triangle') and each triangle's vertex indices with t.get('v0') and so on (here t is an element in 'triangles').