normal vector computation after refine

Hi, I hope this is the right place to report a bug. I'm new in using dolfinx and tried to transfer some of my FEM codes to dolfinx. One feature I do need is an adaptive refinement. Dolfinx offers a refinement based on markers. After refining a 3D surface triangle mesh I found, that the vertex ordering of the refined triangles is not preserved. As consequence the sign of normal vectors depends on the ordering of the vertices. Since I'm new in using dolfinx I don't know if I made a mistake. Below you can find a small piece of code demonstrating my observation.

Best, Markus

import numpy as np
import numba

import ufl
import dolfinx as dfx
from dolfinx import cpp as _cpp

from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner, FiniteElement, Cell, Mesh, VectorElement

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import math
import meshio

def extract_geometricial_data(msh, dim, entities):
    """For a set of entities in a mesh, return the coordinates of the
    vertices"""
    mesh_nodes = []
    geom = msh.geometry
    g_indices = _cpp.mesh.entities_to_geometry(msh, dim, np.array(entities, dtype=np.int32), False)
    for cell in g_indices:
        nodes = np.zeros((len(cell), 3), dtype=np.float64)
        for j, entity in enumerate(cell):
            nodes[j] = geom.x[entity]
        mesh_nodes.append(nodes)
    return mesh_nodes

@numba.njit(fastmath=True)
def mark_all(x):
    return np.ones(x.shape[1])

@numba.njit(fastmath=True)
def get_marked_edges(dim, edges, h, hmax):
    marked_edges = []
    dist_min = 1.0e12
    dist_max = -1.0e12
    for e in edges:
        dist = h[e]
        if dist > dist_max:
            dist_max = dist
        if dist < dist_min:
            dist_min = dist    
        if  dist > hmax:
            marked_edges.append(e)  
    return dist_min, dist_max, np.array(marked_edges)  

def mark_large_edges(msh, dim, hmax):
    """For a set of edge entities in a mesh, return the list of the
    edges largen than hmax"""

    edges = mesh.locate_entities(msh, dim-1, mark_all)
    h = _cpp.mesh.h(msh, dim-1, edges)
    return get_marked_edges(dim, edges, h, hmax)

def get_node_normals(cell_indices, cell_normals, num_nodes):
    node_normals = np.zeros((num_nodes,3))
    node_counts  = np.zeros(num_nodes)

    for cell_id, cell in enumerate(cell_indices):
        for node_id in cell:
            node_counts[node_id]  += 1
            node_normals[node_id] += cell_normals[cell_id]
    for node in range(num_nodes):
        node_normals[node] /= node_counts[node]

    return node_normals

# set up mesh

l0 = 1 
points = [
    [-l0, -l0, 0.0],
    [l0, -l0, 0.0],
    [l0, l0, 0.0],
    [-l0, l0, 0.0]
]
triangles = [("triangle", [[0, 1, 2], [0, 2, 3]])]
mshio = meshio.Mesh(points, triangles)

shape = "triangle"
degree = 1
cell = Cell(shape)
domain = Mesh(VectorElement("Lagrange", cell, degree))
msh = mesh.create_mesh(MPI.COMM_WORLD, mshio.cells_dict['triangle'], mshio.points, domain)
cell_dim = msh.topology.dim
msh.topology.create_entities(cell_dim)

# get some properties and calculate cell normals for first mesh

cell_dim = msh.topology.dim
cells = dfx.mesh.locate_entities(msh, cell_dim, mark_all)
edges = mesh.locate_entities(msh, cell_dim-1, mark_all)
nodes = mesh.locate_entities(msh, cell_dim-2, mark_all)
print("number of cells ", cells.size, ", number of edges ", edges.size, ", number of nodes ", nodes.size )
cell_normals = _cpp.mesh.cell_normals(msh, cell_dim, cells)
print(cell_normals)

cell_nodes = extract_geometricial_data(msh,cell_dim,cells)
facet_nodes = extract_geometricial_data(msh,cell_dim-1,edges)

# refine some times and calculate cell normals for refined meshes

hmin = 0.8

dist_min, dist_max, marked = mark_large_edges(msh,cell_dim,hmin) 

while marked.size > 0:
    dist_min, dist_max, marked = mark_large_edges(msh,cell_dim,hmin) 
    print("dist_min ", dist_min, ", dist_max", dist_max, ", number of edges to be refined: ", marked.size )
    if marked.size > 0:
        msh.topology.create_entities(cell_dim)
        msh = mesh.refine(msh, marked, redistribute=False)
        cells = dfx.mesh.locate_entities(msh, cell_dim, mark_all)
        edges = mesh.locate_entities(msh, cell_dim-1, mark_all)
        nodes = mesh.locate_entities(msh, cell_dim-2, mark_all)
        num_cells = cells.size
        num_edges = edges.size
        num_nodes = nodes.size
        msh.topology.create_entities(cell_dim)
        cell_indices = _cpp.mesh.entities_to_geometry(msh, cell_dim, np.array(cells, dtype=np.int32), False)
        print("number of cells ", cells.size, ", number of edges ", edges.size, ", number of nodes ", nodes.size )
        cell_normals = _cpp.mesh.cell_normals(msh, cell_dim, cells)
        node_normals = get_node_normals(cell_indices, cell_normals, num_nodes)
        print(cell_normals)

FEniCS / dolfinx

normal vector computation after refine #2612