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.
Recently we looked closer into a similar use case for topologically one dimensional meshes (networks). To obtain ordering information on the final mesh one needs to track possible reordering (or reorientation) throughout mesh construction by hand and account for it when one depends on it later.
The bigger question to me is, what do you actually want to do with this “consistently oriented” mesh? Do you want to use it in coupling, or are you just looking to post-process something?
I’m just focusing on the orientation of (internal) boundaries of a mesh, meaning the normal vector should always point “outwards” (or at least in the same direction for each part of the boundary). Also for quad and hex grids the normal vector at their boundaries should have a unique orientation or not? It should just depend in which order the node indices for each facet is stored.
For interior facets there’s of course no general orientation, but for boundary facets you have an outward pointing normal and the orientation of the boundary edge/triangle/quad is just a question of node ordering. I’m not sure, if this is taken into account, when the topology of the mesh is generated, because the facets includes both boundary and internal facets.
If I understand the implementation of
right, when you create a mesh with gmsh, only the topology of the highest topological dimension (codim=0) is used for creating the mesh in dolfinx.
This means you loose topological information about entities with a lower topological dimension, e.g. the orientation of boundary facets (codim=1), which is provided by gmsh, get lost.
For that reason the normal vectors at the boundaries can point in any direction as you can see running the code below with reorient_facets=False.
One has to do some reorientation after the mesh construction like the check_orientation function in my code or include topological information of the boundaries (codim=1) from gmsh in the mesh construction.
Regarding the permute-flag in dolfinx.mesh.entities_to_geometry: I’ve noted it, but it doesn’t do what I want.
If I understand it right, it change local node ordering between dolfinx and gmsh for some element types like quadrilateral, but for edges and triangles it should not make any difference. It does not change the node order of boundary edges and triangles. You can test it with my code below, where you can set permute=True and still have some wrong oriented normal vectors on the outer surface as well as on the interior surface.
@dokken: You asked about the intention of this question. To give you a broader overview, I’m actually considering a fluid-structure interaction (FSI) problem. Therefore, I use a Stokes solver together with an ALE method (additional Poisson equation for the mesh-velocity) to move my mesh. Here you can see 2D and 3D examples for a droplet in a channel.
Since the mesh gets strongly distorted after a few time steps, remeshing (in the fluid domain) is necessary. Therefore, I’ve extracted the mesh of the inner and outer boundaries and remeshed them with gmsh. In the code below, you can see how I do that in 2D and 3D.
First I thought, this would require to have a consistently oriented connectivity array for the edges/triangles. As I’ve found out that actually not necessary, because in 2D you can
set reorient=True in gmsh.model.geo.add_curve_loop, so the orientation of boundary edges doesn’t matter, and in 3D the orientation of the triangles doesn’t matter at all in the way I’ve implemented the remeshing.
So for the remeshing the orientation does not matter. The bigger question here is: How to generalize the remeshing to work with several MPI-processes? Probably one needs to move the whole geometry and topology information of the boundary to a single process and do the actual remeshing only on that process.
There’s another reason why I wanted a consistent boundary orientation in FSI example above. In future I would like to include interfacial mechanics on the interior surface, e.g. surface tension, bending forces, elastic shearing forces, etc. For that it’s also necessary that the normals on the internal surface point outwards everywhere.
Here’s the example code, that does the remeshing and checks the orientation of the normal vectors on boundaries (just works for a single MPI-process).
It’s a bit lengthy, but I wanted to keep the remeshing as general as possible and it should work in 2D and 3D.
import gmsh
import dolfinx
from mpi4py import MPI
import numpy as np
import pyvista as pv
from typing import Sequence
import numpy.typing as npt
def compute_normals(pts: npt.NDArray[np.float64],
f2n: npt.NDArray[np.int32]) -> npt.NDArray[npt.float64]:
"""Compute the normals along each facet.
Args:
pts: Points of the mesh of as array of shape (num_nodes, 3).
f2n: Facet to node mapping as array of shape (num_facets, dim) with dim = 2 or 3. Each row defines the node indices
of an edge (dim = 2) or a triangle (dim = 3).
Returns:
The normals for each facet as array of shape (num_facets, 3).
"""
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: dolfinx.mehs.Mesh,
facets: npt.NDArray[np.int32],
cell_markers: npt.NDArray[np.int32] | None = None,
inner_domain_id: int | None = None,
reorient_faces: bool = False,
permute: bool = True) -> npt.NDArray[np.int32]:
"""Checks the orientation for some facets in a 2D or 3D mesh.
Args:
mesh: Mesh to check.
facets: IDs of the facets to check. Array of shape (num_facets, ).
cell_markers: The tags for each cell inside the mesh.
If not provided all facets must lie on the outer boundary of the domain.
inner_domain_id: The tag corresponding to cell_markers, that is used to identify the inner side of each facet.
If not provided all facets must lie on the outer boundary of the domain.
reorient_faces: Whether to reorient the facets.
permute: Set permute, when calling dolfinx.mesh.entities_to_geometry.
Returns:
An array of integers of shape (facets.shape[0], mesh.topology.dim) containing the connectivity for all edges
if mesh.topology.dim==2 or all triangles if mesh.topology.dim==3. The columns for each edges/triangles will be reoriented, if
reorient_facets is True.
"""
tdim = mesh.topology.dim
fdim = tdim - 1
assert tdim == 2 or tdim == 3
if permute:
mesh.topology.create_entity_permutations()
f2n = dolfinx.mesh.entities_to_geometry(mesh, fdim, facets, permute=permute)
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: npt.NDArray[np.float64],
f2n: npt.NDArray[np.int32],
normals: npt.NDArray[np.float64]) -> None:
"""Plot facets with normals using pyvista
Args:
pts: An array of shape (num_nodes, 3) containing all points of the mesh
f2n: Integer array of shape (num_facets, dim) containing the node indices for all edges (dim = 2) or triangles (dim=3)
normals: An array of shape (num_facets, 3) containing the normals for each facet.
"""
# 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 plot_mesh(mesh: dolfinx.mesh.Mesh,
tags: npt.NDArray[np.int32] | None = None) -> None:
"""Plot the mesh (dim=2) or a slice of the mesh (dim=3).
Args:
mesh: The mesh to plot
tags: An array with tags for each cell, that will be used for coloring, if provided.
"""
gdim = mesh.geometry.dim
tdim = mesh.topology.dim
pv_mesh = pv.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh, gdim))
if tags is not None:
pv_mesh.cell_data["tags"] = tags
pl = pv.Plotter()
scalars = None if tags is None else "tags"
if tdim == 2:
pl.add_mesh(pv_mesh, show_edges=True, scalars=scalars, color="blue", edge_color="yellow", cmap="coolwarm")
pl.view_xy()
else:
slice = pv_mesh.slice(origin=[0.5, 0.5, 0.5], normal=[0.0, 0.0, 1.0])
pl.add_mesh(slice, show_edges=True, color="blue", edge_color="yellow", scalars=scalars, cmap="coolwarm")
pl.view_xy()
pl.show()
def create_mesh(dim: int,
radius: flaot = 1.0,
length: float = 4.0,
mesh_width: float | None = None) -> tuple[dolfinx.io.gmsh.MeshData, dict[str, int]]:
""" Create a mesh for a circle in a rectangle for dim = 2 and a sphere in a box for dim = 3 using gmsh.
Args:
dim: The dimension of the mesh
radius: The radius of the circle/sphere.
length: The length of the rectangle/box.
mesh_width: The mesh width.
Returns:
The mesh data (mesh and mesh tags) and a dictionary mapping each facet and subdomain to its tag.
"""
# 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:
tag_info = {"inner_boundary": 1,
"left": 2, "right": 3,
"bottom": 4, "top": 5,
"inner_domain": 10, "outer_domain": 11}
circle_id = 1
outer_id = 2
circle = model.occ.add_disk(0, 0, 0, radius, radius, tag=circle_id)
rectangle = model.occ.add_rectangle(-0.5 * length, -0.5 * length, 0, length, length)
outer_domain_dim_tags, _ = model.occ.cut([(dim, rectangle)], [(dim, circle)], removeTool=False)
model.occ.synchronize()
# Group tags of boundaries
outer_domain_boundaries = model.get_boundary(outer_domain_dim_tags)
circle_boundaries = model.get_boundary([(dim, circle)])
outer_boundaries = [tag for dim2, tag in outer_domain_boundaries if
dim2 == 1 and (dim2, tag) not in circle_boundaries]
com_outer_boundaries = np.array([model.occ.get_center_of_mass(1, tag) for tag in outer_boundaries])
left, right, top, bottom = -1, -1, -1, -1
for com, tag in zip(com_outer_boundaries, outer_boundaries):
if np.allclose(com, [-0.5 * length, 0.0, 0.0]):
left = tag
elif np.allclose(com, [0.5 * length, 0.0, 0.0]):
right = tag
elif np.allclose(com, [0.0, -0.5 * length, 0.0]):
bottom = tag
elif np.allclose(com, [0.0, 0.5 * length, 0.0]):
top = tag
else:
raise Exception(f"Boundary with com = {com}, tag = {tag} cannot be identified.")
assert left > 0 and right > 0 and top > 0 and bottom > 0
# Physical groups for edges
model.add_physical_group(1, [left], tag=tag_info["left"])
model.add_physical_group(1, [right], tag=tag_info["right"])
model.add_physical_group(1, [bottom], tag=tag_info["bottom"])
model.add_physical_group(1, [top], tag=tag_info["top"])
model.add_physical_group(1, [tag for dim2, tag in circle_boundaries if dim2 == 1], tag=tag_info["inner_boundary"])
# Physical groups for surfaces
model.add_physical_group(2, [circle], tag=tag_info["inner_domain"])
model.add_physical_group(2, [tag for dim2, tag in outer_domain_dim_tags if dim2 == 2], tag=tag_info["outer_domain"])
else:
tag_info = {"inner_boundary": 1,
"left": 2, "right": 3,
"bottom": 4, "top": 5,
"back": 6, "front": 7,
"inner_domain": 10, "outer_domain": 11}
sphere = model.occ.add_sphere(0, 0, 0, 1.0, tag=tag_info["inner_domain"])
box = model.occ.add_box(-0.5 * length, -0.5 * length, -0.5 * length, length, length, length)
outer_domain_dim_tags, _ = model.occ.cut([(dim, box)], [(dim, sphere)], removeTool=False)
model.occ.synchronize()
# Group tags of boundaries
outer_domain_boundaries = model.get_boundary(outer_domain_dim_tags)
sphere_boundaries = model.get_boundary([(dim, sphere)])
outer_boundaries = [tag for dim2, tag in outer_domain_boundaries if
dim2 == 2 and (dim2, tag) not in sphere_boundaries]
com_outer_boundaries = np.array([model.occ.get_center_of_mass(2, tag) for tag in outer_boundaries])
left, right, top, bottom, back, front = -1, -1, -1, -1, -1, -1
for com, tag in zip(com_outer_boundaries, outer_boundaries):
if np.allclose(com, [-0.5 * length, 0.0, 0.0]):
left = tag
elif np.allclose(com, [0.5 * length, 0.0, 0.0]):
right = tag
elif np.allclose(com, [0.0, -0.5 * length, 0.0]):
bottom = tag
elif np.allclose(com, [0.0, 0.5 * length, 0.0]):
top = tag
elif np.allclose(com, [0.0, 0.0, -0.5 * length]):
back = tag
elif np.allclose(com, [0.0, 0.0, 0.5 * length]):
front = tag
else:
raise Exception(f"Boundary with com = {com}, tag = {tag} cannot be identified.")
assert left > 0 and right > 0 and top > 0 and bottom > 0 and back > 0 and front > 0
# Physical groups for surfaces
model.add_physical_group(2, [left], tag=tag_info["left"])
model.add_physical_group(2, [right], tag=tag_info["right"])
model.add_physical_group(2, [bottom], tag=tag_info["bottom"])
model.add_physical_group(2, [top], tag=tag_info["top"])
model.add_physical_group(2, [back], tag=tag_info["back"])
model.add_physical_group(2, [front], tag=tag_info["front"])
model.add_physical_group(2, [tag for dim2, tag in sphere_boundaries if dim2 == 2], tag=tag_info["inner_boundary"])
# Physical groups for volumes
model.add_physical_group(3, [sphere], tag=tag_info["inner_domain"])
model.add_physical_group(3, [tag for dim2, tag in outer_domain_dim_tags if dim2 == 3], tag=tag_info["outer_domain"])
if mesh_width is not None:
gmsh.option.set_number("Mesh.MeshSizeMin", mesh_width)
gmsh.option.set_number("Mesh.MeshSizeMax", mesh_width)
model.mesh.generate(dim=dim)
mesh_data = dolfinx.io.gmsh.model_to_mesh(model, MPI.COMM_WORLD, rank=0, gdim=dim)
gmsh.finalize()
return mesh_data, tag_info
def remesh_2d_mesh(mesh: dolfinx.mesh.Mesh,
facets_list: Sequence[npt.NDArray[np.int32]],
facet_ids: Sequence[int],
curve_loop_facet_indices: Sequence[Sequence[int]],
surfaces: Sequence[Sequence[int]],
domain_ids: Sequence[int],
mesh_width: float = 0.0,
permute: bool = False) -> dolfinx.io.gmsh.MeshData:
""" Remesh the surface of two-dimensional mesh providing all facets and curves, that are external and internal boundaries,
which should be kept fixed.
Args:
mesh: The 2D-mesh to be remeshed.
facets_list: A list containing different parts of the boundary. Each entry contains an integer array marking the facet IDs
of this part of the boundary.
facet_ids: Tags for each facet in facets_list.
curve_loop_facet_indices: A list of curve loops. Each curve loop is defined by a subset of facets, e.g. curve_loop_facets_indices[0] == [0, 2],
would mean the first and third facet in facets_list form the first curve loop. The orientation of the facets doesn't matter
surfaces: List of surfaces defined by the curve loops, they are bounded by, e.g. surfaces[0] == [0, 2] means, the first surface is bounded
by the first and the third curve loop. The first curve loop should mark the outer boundary of the surface.
domain_ids: Tags for each surface in surfaces.
mesh_width: Mesh width, that is used for mesh generation.
permute: Set permute, when calling dolfinx.mesh.entities_to_geometry.
Returns:
The mesh data (mesh and mesh tags)
"""
tdim = mesh.topology.dim
assert tdim == 2
assert len(facet_ids) == len(facets_list)
assert len(surfaces) == len(domain_ids)
fdim = tdim - 1
gmsh.initialize()
gmsh.option.set_number("General.Terminal", 0)
model = gmsh.model
name = "remesh_mesh"
model.add(name)
model.set_current(name)
all_facets = np.hstack(facets_list)
boundary_node_ids = dolfinx.mesh.compute_incident_entities(mesh.topology, all_facets, fdim, 0)
boundary_points = mesh.geometry.x[boundary_node_ids].ravel()
# gmsh tags must be positive
boundary_node_ids += 1
# Add nodes
for i, node_id in enumerate(boundary_node_ids):
x, y, z = boundary_points[3 * i:3 * i + 3]
model.geo.add_point(x, y, z, tag=node_id, meshSize=mesh_width)
all_line_tags = []
# Add lines
for facets in facets_list:
if permute:
mesh.topology.create_entity_permutations()
f2n = dolfinx.mesh.entities_to_geometry(mesh, fdim, facets, permute=permute)
line_tags = facets + 1
node_tags_lines = f2n + 1
for line_tag, node_tags_line in zip(line_tags, node_tags_lines, strict=True):
model.geo.add_line(*node_tags_line, tag=line_tag)
all_line_tags.append(line_tags)
# Add curve loops
curve_loops = []
for relative_facet_indices in curve_loop_facet_indices:
line_tags_curve_loop = []
for facet_id in relative_facet_indices:
line_tags_curve_loop.extend(all_line_tags[facet_id])
curve_loop = model.geo.add_curve_loop(line_tags_curve_loop, reorient=True)
curve_loops.append(curve_loop)
# Add surfaces
for curve_loop_indices, domain_id in zip(surfaces, domain_ids, strict=True):
curve_loops_surface = [curve_loops[i] for i in curve_loop_indices]
surface = model.geo.add_plane_surface(curve_loops_surface, tag=domain_id)
model.geo.synchronize()
# Add physical groups for facets (lines)
for i, facet_id in enumerate(facet_ids):
name = f"Facet_{facet_id}"
model.add_physical_group(fdim, all_line_tags[i], name=name, tag=facet_id)
# Add physical groups for subdomains
for domain_id in domain_ids:
model.add_physical_group(tdim, [domain_id], name=f"Domain_{domain_id}", tag=domain_id)
# Generate the mesh and convert it to dolfinx
model.mesh.generate(tdim)
mesh_data = dolfinx.io.gmsh.model_to_mesh(model, mesh.comm, rank=0, gdim=mesh.geometry.dim)
gmsh.finalize()
return mesh_data
def remesh_3d_mesh(
mesh: dolfinx.mesh.Mesh,
facets_list: Sequence[npt.NDArray[np.int32]],
facet_ids: Sequence[int],
surface_loop_facet_indices: Sequence[Sequence[int]],
volumes: Sequence[Sequence[int]],
domain_ids: Sequence[int],
mesh_width: float | None = None,
permute: bool = False
) -> dolfinx.io.gmsh.MeshData:
""" Remesh the volume of three-dimensional mesh providing all facets and surface loops, that are external and internal boundaries,
which should be kept fixed.
Args:
mesh: The 3D-mesh to be remeshed.
facets_list: A list containing different parts of the boundary. Each entry contains an integer array marking the facet IDs
of this part of the boundary.
facet_ids: Tags for each facet in facets_list.
surface_loop_facet_indices: A list of surface loops. Each surface loop is defined by a subset of facets, e.g.
surface_loop_facets_indices[0] == [0, 2], would mean the first and third facet in facets_list form the
first surface loop. The orientation of the facets doesn't matter.
volumes: List of volumes defined by the surface loops, they are bounded by, e.g. volumes[0] == [0, 2] means, the first volumes is bounded
by the first and the third surface loop. The first surface loop should mark the outer boundary of the volume.
domain_ids: Tags for each volume in volumes.
mesh_width: Mesh width, that is used for mesh generation.
permute: Set permute, when calling dolfinx.mesh.entities_to_geometry.
Returns:
The mesh data (mesh and mesh tags)
"""
tdim = mesh.topology.dim
assert tdim == 3
assert len(facet_ids) == len(facets_list)
assert len(volumes) == len(domain_ids)
fdim = tdim - 1
gmsh.initialize()
gmsh.option.set_number("General.Terminal", 0)
model = gmsh.model
name = "remesh_mesh"
model.add(name)
model.set_current(name)
all_facets = np.hstack(facets_list)
boundary_node_ids = dolfinx.mesh.compute_incident_entities(mesh.topology, all_facets, fdim, 0)
boundary_points = mesh.geometry.x[boundary_node_ids].ravel()
# gmsh tags must be positive
boundary_node_ids += 1
# Add nodes
nodes_tag = 1
model.add_discrete_entity(dim=0, tag=nodes_tag)
model.mesh.add_nodes(dim=0, tag=nodes_tag, nodeTags=boundary_node_ids, coord=boundary_points)
# Add facets
for facets, facet_id in zip(facets_list, facet_ids, strict=True):
if permute:
mesh.topology.create_entity_permutations()
f2n = dolfinx.mesh.entities_to_geometry(mesh, fdim, facets, permute=permute)
model.add_discrete_entity(dim=fdim, tag=facet_id)
element_tags = facets + 1
node_tags_element = f2n.ravel() + 1
triangle_id = 2 # ID for a 3-node-triangle in gmsh
model.mesh.add_elements_by_type(tag=facet_id, elementType=triangle_id, elementTags=element_tags,
nodeTags=node_tags_element)
# Add surface loops
surface_loops = []
for relative_facet_indices in surface_loop_facet_indices:
facet_tags_surface_loop = [facet_ids[i] for i in relative_facet_indices]
surface_loop = model.geo.add_surface_loop(facet_tags_surface_loop)
surface_loops.append(surface_loop)
# Add volumes
for surface_loop_indices, domain_id in zip(volumes, domain_ids, strict=True):
surface_loops_volume = [surface_loops[i] for i in surface_loop_indices]
volume = model.geo.add_volume(surface_loops_volume, tag=domain_id)
model.geo.synchronize()
# Add physical groups for facets (surfaces)
for facet_id in facet_ids:
name = f"Facet_{facet_id}"
model.add_physical_group(fdim, [facet_id], name=name, tag=facet_id)
# Add physical groups for domains
for domain_id in domain_ids:
model.add_physical_group(tdim, [domain_id], name=f"Domain_{domain_id}", tag=domain_id)
if mesh_width is not None:
gmsh.option.set_number("Mesh.MeshSizeMin", mesh_width)
gmsh.option.set_number("Mesh.MeshSizeMax", mesh_width)
# Generate the mesh and convert it to dolfinx
model.mesh.generate(tdim)
mesh_data = dolfinx.io.gmsh.model_to_mesh(model, mesh.comm, rank=0, gdim=mesh.geometry.dim)
gmsh.finalize()
return mesh_data
def main():
"""Main function"""
dim = 3 # Dimension of the mesh
permute = False # Set permute, when calling dolfinx.mesh.entities_to_geometry?
reorient_facets = False # Reorient the facets/normals?
# Properties of the mesh
radius = 1.0
length = 3.0
mesh_width = 0.5
# Generate the mesh
mesh_data, tag_info = create_mesh(dim, radius, length, mesh_width)
mesh = mesh_data.mesh
cell_tags = mesh_data.cell_tags
facet_tags = mesh_data.facet_tags
tdim = mesh.topology.dim
fdim = tdim - 1
# Plot the initial mesh
plot_mesh(mesh, tags=cell_tags.values)
# Identify different boundaries of the mesh
left_facets = facet_tags.find(tag_info["left"])
right_facets = facet_tags.find(tag_info["right"])
bottom_facets = facet_tags.find(tag_info["bottom"])
top_facets = facet_tags.find(tag_info["top"])
if tdim == 2:
boundary_facets = np.hstack([left_facets, right_facets, bottom_facets, top_facets])
else:
back_facets = facet_tags.find(tag_info["back"])
front_facets = facet_tags.find(tag_info["front"])
boundary_facets = np.hstack([left_facets, right_facets, bottom_facets, top_facets, back_facets, front_facets])
# Compute and plot normals for the outer and the internal surface
f2n_outer = check_orientation(mesh, facets=boundary_facets, permute=permute, reorient_faces=reorient_facets)
normals_outer = compute_normals(mesh.geometry.x, f2n_outer)
plot_normals(mesh.geometry.x, f2n_outer, normals_outer)
assert (cell_tags.indices == np.arange(cell_tags.indices.shape[0])).all()
inner_facets = facet_tags.find(tag_info["inner_boundary"])
f2n_inner = check_orientation(mesh, facets=inner_facets, cell_markers=cell_tags.values,
inner_domain_id=tag_info["inner_domain"], permute=permute,
reorient_faces=reorient_facets)
normals_inner = compute_normals(mesh.geometry.x, f2n_inner)
plot_normals(mesh.geometry.x, f2n_inner, normals_inner)
# Remeshing of the domain
if dim == 2:
facets_list = [left_facets, right_facets, bottom_facets, top_facets, inner_facets]
facet_ids = [tag_info["left"], tag_info["right"], tag_info["bottom"], tag_info["top"], tag_info["inner_boundary"]]
# Each entry of the list corresponds to a single curve loop. Each entry define the indices of the facets from
# facet_list above a curve_loop consists of.
curve_loops = [[0, 1, 2, 3], [4]] # Outer rectangle and circle boundary
# Each entry defines a subdomain. Each entry define the indices of the curve loops, which are the boundaries of
# the subdomain
domains = [[0, 1], [1]] # Rectangle minus circle and circle
domain_ids = [tag_info["outer_domain"], tag_info["inner_domain"]]
mesh_data_remeshed = remesh_2d_mesh(mesh, facets_list, facet_ids, curve_loops, domains, domain_ids,
permute=permute, mesh_width=mesh_width)
elif dim == 3:
facets_list = [left_facets, right_facets, bottom_facets, top_facets, back_facets, front_facets, inner_facets]
facet_ids = [tag_info["left"], tag_info["right"], tag_info["bottom"], tag_info["top"],
tag_info["back"], tag_info["front"], tag_info["inner_boundary"]]
# Each entry of the list corresponds to a single surface loop. Each entry define the indices of the facets from
# facet_list above a surface_loop consists of.
surface_loops = [[0, 1, 2, 3, 4, 5], [6]] # Outer box and sphere surface
# Each entry defines a subdomain. Each entry define the indices of the curve loops, which are the boundaries of
# the subdomain
domains = [[0, 1], [1]] # Box minus sphere and sphere
domain_ids = [tag_info["outer_domain"], tag_info["inner_domain"]]
mesh_data_remeshed = remesh_3d_mesh(mesh, facets_list, facet_ids, surface_loops, domains, domain_ids,
permute=permute, mesh_width=mesh_width)
else:
assert False
remeshed_mesh = mesh_data_remeshed.mesh
cell_tags = mesh_data_remeshed.cell_tags
assert (cell_tags.indices == np.arange(cell_tags.indices.shape[0])).all()
cell_tags = cell_tags.values
plot_mesh(remeshed_mesh, cell_tags)
if __name__ == '__main__':
main()
ufl.FacetNormal will always point outwards of a cell, so therefore I am curious as why you need this extra extraction of entities. You can use ufl.Expression with facet entities to tabulate these, as shown in: FacetNormal vector components and director vector - #9 by dokken
There are several examples of this on the forum and in scifem.
The reorient function in legacy FEniCS relates to the CellNormal, as well as the orientation of facets shared between two cells. The issue I am mentioning above is illustrated in figure 4. Of v3 (arxiv) of https://arxiv.org/pdf/2102.11901v3
For the remeshing, you would have to gather the boundary on a single process. This is due to the fact that very few mesh generation tools run with MPI.
