Hi everyone,
I’m trying to extract the (oriented) connectivity array for a surface of my mesh, which can lie on the boundary or be an internal surface as well. I know you can simply get the connectivity array calling f2n = dolfinx.mesh.entities_to_geometry(mesh, fdim, facet_ids) with fdim=1 (2D) or fdim=2 (3D) and facet_ids marking all facets (edges/triangles) of the surface in consideration.
The problem is, that this connectivity array does not yield any information about the orientation of the facets.
In 2D I would expect, if \mathbf{x}_0 and \mathbf{x}_1 are the nodes of an edge (ordered according to the columns of f2n above), and \mathbf{e} = \mathbf{x}_1 - \mathbf{x}_0 be the edge vector, then the normal to this edge is given by \mathbf{n} = (e_2, -e_1)^T / ||\mathbf{e}||_2.
Similar in 3D, the normal vector of a triangle with nodes \mathbf{x}_0, \mathbf{x}_1 and \mathbf{x}_2 (again ordered according to the columns of f2n above) the normal vector should be \mathbf{n} = \frac{(\mathbf{x}_1 - \mathbf{x}_0) \times (\mathbf{x}_2 - \mathbf{x}_0)}{||(\mathbf{x}_1 - \mathbf{x}_0) \times (\mathbf{x}_2 - \mathbf{x}_0)||_2}.
Is there a function in dolfinx to get the f2n array with sorted columns according to the orientation of the edges/triangles? I know in general it might be not possible to get a unique orientation for interior surfaces, but I would at least expect that you can reorient the outer boundary consistently. What does determine the order of the node indices, if I call f2n = dolfinx.mesh.entities_to_geometry(mesh, fdim, facet_ids) as above? Especially, if the mesh was created with gmsh, curves and surfaces have a consistent orientation, which can be destroyed, when calling dolfinx.io.gmsh.model_to_mesh (probably due to the mesh-partitioner?).
I’ve also seen the issues Computing cell orientation in FEniCSx and Wrong FacetNormal vector on internal boundaries regarding orientation problems, but I found no intrinsic dolfinx function to reorient the surface. In the post Computing cell orientation in FEniCSx it was mentioned in legacy fenics there was a function mesh.init_cell_orientation() and I’m searching for something similar in FEniCSx.
Here’s an example code, that demonstrate that in general the direction of the normal vectors are wrong/inconsistent (setting reorient_faces = False). It should work in 2D and 3D for the outer boundary of the domain (interior_surface = False) as well as for an interior surface (interior_surface = True) on a single process.
import pyvista as pv
import dolfinx
from mpi4py import MPI
import ufl
import numpy as np
import gmsh
def compute_normals(pts, f2n):
num_facets, dim = f2n.shape
# First and second node of the edges/triangles
x0 = pts[f2n[:, 0]]
x1 = pts[f2n[:, 1]]
# Compute the normals
if dim == 2:
e = x1 - x0
normals = np.zeros((num_facets, 3), dtype=x0.dtype)
normals[:, 0] = e[:, 1]
normals[:, 1] = -e[:, 0]
else:
# Third node of the triangles
x2 = pts[f2n[:, 2]]
normals = np.cross(x1 - x0, x2 - x0)
normals /= np.linalg.norm(normals, axis=1).reshape(-1, 1)
return normals
def check_orientation(mesh, facets, cell_markers=None, inner_domain_id=None, reorient_faces=False):
tdim = mesh.topology.dim
fdim = tdim - 1
assert tdim == 2 or tdim == 3
f2n = dolfinx.mesh.entities_to_geometry(mesh, fdim, facets)
f2c = mesh.topology.connectivity(fdim, tdim)
c2n = mesh.geometry.dofmaps[0]
num_facets = facets.shape[0]
# Nodes of the edges or triangles with shape (num_facets, fdim+1, 3)
# The second dimension corresponds to nodes of the edges/triangles
nodes_facet = mesh.geometry.x[f2n[:]]
# The third/fourth node of the triangle/tetrahedron neigbored to the boundary edge/triangle
nodes_cell = np.empty((num_facets, 3), dtype=nodes_facet.dtype)
is_interior_facet = cell_markers is not None and inner_domain_id is not None
# Find the fourth node in the tetrahedron neigbored to this triangle
for i, facet in enumerate(facets):
cells = f2c.links(facet)
if is_interior_facet:
# We search for the cell with cell_marker equal to inner_domain_id
ids_inner_cells = np.flatnonzero(cell_markers[cells] == inner_domain_id)
# There can only be one interior cell neighbored to the facet
assert(len(ids_inner_cells) == 1)
cell_id = cells[ids_inner_cells[0]]
else:
# A boundary facet can only be neighbored to one cell
assert(len(cells) == 1)
cell_id = cells[0]
cell_node_ids = c2n[cell_id]
# Find the index in cell_node_ids, that is not contained in f2n[i, :]
neq_facet_id = cell_node_ids != f2n[i, 0]
for j in range(1, f2n.shape[1]):
neq_facet_id = np.logical_and(neq_facet_id, cell_node_ids != f2n[i, j])
local_id_remaining_node = np.argwhere(neq_facet_id)
assert local_id_remaining_node.size == 1
nodes_cell[i] = mesh.geometry.x[cell_node_ids[local_id_remaining_node[0, 0]]]
normals = compute_normals(mesh.geometry.x, f2n)
# Check the orientation of the normal vector computing the scalar product
# between the normal vector and the vector conecting the facet center with the remaining node
# of the cell owning this boundary facet
facet_centers = np.mean(nodes_facet, axis=1)
sign_facets = np.sign(np.sum(normals * (facet_centers - nodes_cell), axis=1))
orientation_wrong = sign_facets < 0
ids_wrong_facets = np.argwhere(orientation_wrong)
if ids_wrong_facets.shape[0] > 0:
id_wrong_0 = ids_wrong_facets[0, 0]
print("There are wrong oriented facets. The first one has")
print(f" Index = {id_wrong_0}")
print(f" Normal = {normals[id_wrong_0]}")
print(f" Center = {facet_centers[id_wrong_0]}")
else:
print("All facets have the right orientation.")
if reorient_faces:
# Reverse the order of the node indices for wrong oriented facets
f2n_oriented = f2n.copy()
f2n_oriented[orientation_wrong] = np.fliplr(f2n[orientation_wrong])
return f2n_oriented
else:
return f2n
def plot_normals(pts, f2n, normals):
# Plot surface with normals in pyvista
num_facets, dim = f2n.shape
padding = np.full((num_facets, 1), dim, dtype=f2n.dtype)
cells = np.hstack((padding, f2n)).ravel()
cell_type = pv.CellType.LINE if dim == 2 else pv.CellType.TRIANGLE
cell_types = np.full(num_facets, cell_type, dtype=np.uint8)
grid = pv.UnstructuredGrid(cells, cell_types, pts)
grid.cell_data["normals"] = normals
nomal_glyphs = grid.glyph(orient="normals", scale=True, factor=0.3)
plotter = pv.Plotter()
color = "black" if dim == 2 else "lightblue"
plotter.add_mesh(grid, color=color, line_width=5, show_edges=True, opacity=0.5)
plotter.add_mesh(nomal_glyphs, color="blue")
if dim == 2:
plotter.view_xy()
plotter.add_axes()
plotter.show()
def create_mesh(dim):
# Create a mesh for a circle inside a rectangle in 2D
# and a mesh for a sphere in a box in 3D
gmsh.initialize()
gmsh.option.set_number("General.Terminal", 0)
model = gmsh.model()
name = "circle_in_rect" if dim == 2 else "sphere_in_box"
model.add(name)
model.setCurrent(name)
if dim == 2:
circle_id = 1
outer_id = 2
circle = model.occ.add_disk(0, 0, 0, 1.0, 1.0, tag=circle_id)
rectangle = model.occ.add_rectangle(-2.0, -2.0, 0, 4.0, 4.0)
outer_domain_dim_tags, _ = model.occ.cut([(dim, rectangle)], [(dim, circle)], removeTool=False)
model.occ.synchronize()
# Physical groups for edges
edges = gmsh.model.get_entities(1)
for dim2, tag in edges:
gmsh.model.add_physical_group(dim2, [tag], tag)
# Physical groups for surfaces
gmsh.model.add_physical_group(2, [circle], tag=circle_id)
gmsh.model.add_physical_group(2, [tag for dim2, tag in outer_domain_dim_tags if dim2 == 2], tag=outer_id)
model.mesh.generate(dim=2)
else:
sphere_id = 1
outer_id = 2
sphere = model.occ.add_sphere(0, 0, 0, 1.0, tag=sphere_id)
box = model.occ.add_box(-2.0, -2.0, -2.0, 4.0, 4.0, 4.0)
outer_domain_dim_tags, _ = model.occ.cut([(dim, box)], [(dim, sphere)], removeTool=False)
model.occ.synchronize()
# Physical groups for surfaces
surfaces = gmsh.model.get_entities(2)
volumes = gmsh.model.get_entities(3)
for dim2, tag in surfaces:
gmsh.model.add_physical_group(dim2, [tag], tag)
# Physical groups for volumes
gmsh.model.add_physical_group(3, [sphere], tag=sphere_id)
gmsh.model.add_physical_group(3, [tag for dim2, tag in outer_domain_dim_tags if dim2 == 3], tag=outer_id)
model.mesh.generate(dim=3)
mesh_data = dolfinx.io.gmsh.model_to_mesh(model, MPI.COMM_WORLD, rank=0, gdim=dim)
gmsh.finalize()
return mesh_data
def main():
dim = 2 # Dimension of the mesh
interior_surface = True # If true an interior surface (here x=0.5), if false the boundary of the mesh
reorient_faces = False # Reoriente the faces if set to true. Otherwise just plot the (partially wrong oriented) normals
assert dim == 2 or dim == 3
mesh_data = create_mesh(dim)
mesh = mesh_data.mesh
facet_tags = mesh_data.facet_tags
cell_tags = mesh_data.cell_tags
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
if interior_surface:
# 1 is the tag for the surface of the circle/sphere
facets = facet_tags.find(1)
inner_domain_id = 1
num_cells = cell_tags.indices.shape[0]
assert (cell_tags.indices == np.arange(num_cells, dtype=cell_tags.indices.dtype)).all()
cell_markers = cell_tags.values
else:
facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
cell_markers = None
inner_domain_id = None
f2n = check_orientation(mesh, facets, cell_markers=cell_markers, inner_domain_id=inner_domain_id, reorient_faces=reorient_faces)
normals = compute_normals(mesh.geometry.x, f2n)
plot_normals(mesh.geometry.x, f2n, normals)
if __name__ == "__main__":
main()
The check_orientation function actually contains a workaround for my problem (using a similar approach as in Wrong FacetNormal vector on internal boundaries - #2 by dokken, set reorient_faces=True to illustrate it), but it might become inefficient for high mesh resolutions due to the for loop in python and it does not work in parallel yet. So I’m searching for a better solution.
General background of the question: I want to do remeshing with gmsh (reading the contours/surfaces and create a new meshes for the areas/volumes). Therefore, I need to read the boundary points and boundary facets. The latter requires an oriented f2n array, because otherwise gmsh throw errors due to wrong oriented edges/triangles.
Thanks for your help.