# Normal vector computation after refine

Hello,
I want to use 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 and changes it sign with refinements. Below you can find a small piece of code demonstrating my observation. From my opinion the cell normals should always point into the same direction since the surface is just a part of the xy-plane.

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)

@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)
``````