Relation of midpoints to node ID with 2nd order meshes

mesh
1

I simulate FSI problems on a single 2nd order mesh in three-dimension. I need to retrieve the node coordinates and connectivity of vertices and midpoints on each fluid-structure interface facet, respectively. Say if there are n triangular facets on the surface of interest, I want to retrieve two nx3 arrays vertex_indices and midpoint_indices (a facet with 3 vertices + 3 midpoints). So that on ith facet, I have connection of vertices and midpoints by fsi_mesh.geometry.x[vertex_indices[i]] and fsi_mesh.geometry.x[midpoint_indices[i]].

The attached MWE and figure show my process to retrieve connection of vertices and midpoints on the surface of a spherical mesh. For vertices, the connection can always be found by Facet (topology) -> Vertices (topology) -> Node (geometry) as

# Facet (topology) -> Vertices (topology) -> Node (geometry) 
vertices = mesh.compute_incident_entities(topology, np.array([facet], dtype=np.int32), fdim, 0)
nodes = mesh.entities_to_geometry(_mesh, 0, np.array(vertices, dtype=np.int32)).reshape(-1)
vertex_coords = _mesh.geometry.x[nodes]

However, the relation of midpoints to node ID with 2nd order meshes is not explicit. I first create a mesh-matched 2nd order Lagrange functionspace and use the found 6 DOF on a facet to match which are vertices and which are midpoints (this is must not the optimal one) as

# For each facet, find the 3 vertices by matching the coordinates of vertices, 
# then the remaining 3 entries are the midpoints.
vertex_index = []
for i in range(3):       # 3 vertices
    for j in range(6):   # total 6 DOFs: 3 vertices + 3 midpoints
        if np.linalg.norm(vertex_coords[i] - V_coord[dof_at_facet[j]]) < 1e-6:
            vertex_index.append(dof_at_facet[j])
            break

assert len(vertex_index) == 3, f"len(vertex_ind  ex) = {len(vertex_index)}"    
vertex_indices.append(vertex_index)
midpoint_indices.append(np.setdiff1d(dof_at_facet, vertex_index))

Questions:

  1. Any recommended ways to retrieve the connectivity of midpoints on 2nd order meshes?
  2. Any potential issues related to my process as the found midpoint_indices is in terms of topology, as there is a missing conversion from topology to geometry. (e.g. the order of a Lagrange functionspace DOF to node order of an independent mesh is matched, however the submesh does not follow the same order.)
# mwe.py
# Run $ python mwe.py with DOLFINx v0.9.0

import gmsh
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, io 

comm = MPI.COMM_WORLD
mesh_order = 2
name = "sphere"

# Gmsh to build the mesh
gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 2)
model = gmsh.model()
model.add(name)
model.setCurrent(name)
sphere = model.occ.addSphere(0, 0, 0, 1, tag=1)
model.occ.synchronize()
model.add_physical_group(dim=3, tags=[sphere])
model.mesh.generate(dim=3)
model.mesh.setOrder(mesh_order)

# Load to DOLFINx
_mesh, _, _ = io.gmshio.model_to_mesh(model, comm, rank=0)

# Find facets on boundary of interest
topology = _mesh.topology
tdim = topology.dim
fdim = tdim - 1
topology.create_connectivity(fdim, tdim)
topology.create_connectivity(0, tdim)
boundary_facets = mesh.exterior_facet_indices(topology)

# Build function space and use DOF to find the coordinates
V = fem.functionspace(_mesh, ("Lagrange", mesh_order))
vertex_indices = []
midpoint_indices = []

V_coord = V.tabulate_dof_coordinates()

for facet in boundary_facets:
    dof_at_facet = fem.locate_dofs_topological(V, fdim, [facet])

    assert len(dof_at_facet) == 3 * mesh_order, f"len(dof_at_facet) = {len(dof_at_facet)}"
    
    # Facet (topology) -> Vertices (topology) -> Node (geometry) 
    vertices = mesh.compute_incident_entities(topology, np.array([facet], dtype=np.int32), fdim, 0)
    nodes = mesh.entities_to_geometry(_mesh, 0, np.array(vertices, dtype=np.int32)).reshape(-1)
    vertex_coords = _mesh.geometry.x[nodes]

    # For each facet, find the 3 vertices by matching the coordinates of vertices, 
    # then the remaining 3 entries are the midpoints.
    vertex_index = []
    for i in range(3):
        for j in range(6):
            if np.linalg.norm(vertex_coords[i] - V_coord[dof_at_facet[j]]) < 1e-6:
                vertex_index.append(dof_at_facet[j])
                break
            
    assert len(vertex_index) == 3, f"len(vertex_index) = {len(vertex_index)}"
    
    vertex_indices.append(vertex_index)
    midpoint_indices.append(np.setdiff1d(dof_at_facet, vertex_index))

# Load to numpy array with shape n x 3 
vertex_indices = np.array(vertex_indices)
midpoint_indices = np.array(midpoint_indices)

# Validation
assert np.all(np.isclose(_mesh.geometry.x[vertex_indices], V_coord[vertex_indices]))
assert np.all(np.isclose(_mesh.geometry.x[midpoint_indices], V_coord[midpoint_indices]))

# Visualization
import matplotlib.pyplot as plt
def plot_surface(points, triangles):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # Extract x, y, z coordinates from points
    x = points[:, 0]
    y = points[:, 1]
    z = points[:, 2]

    # Plot the surface using triangles
    ax.plot_trisurf(x, y, z, triangles=triangles, edgecolor='black', alpha=0.5)

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    plt.show()
    
plot_surface(V_coord, vertex_indices)  # plot triangles by vertices
# plot_surface(V_coord, midpoint_indices)  # plot triangles by midpoints

image
image1409×639 81.3 KB