Create mesh from arrays: How to preserve order of vertices?

dolfinx
1

I want to load an external tetrahedral mesh into dolfinx by loading an array containing the vertex coordinates as well as an array specifying the connectivity. For my application it is important, that the order of the vertices is preserved, while loading them into dolfinx. I.e., that the vertices are stored in mesh.geometry.x in the same order as in the input array.

But unfortunately, dolfinx seems to change the order. Cf. the following minimal working example:

import numpy as np
from mpi4py import MPI
import dolfinx as dfx # type: ignore
import basix.ufl

#Function for generating a dolfinx mesh from connectivity
def mesh_from_connectivity(cells, vertices):
    cell_type = "tetrahedron"  
    degree = 1
    domain = dfx.mesh.ufl.Mesh(basix.ufl.element("Lagrange", cell_type, degree, basix.LagrangeVariant.equispaced, shape=(3, )))
    return dfx.mesh.create_mesh(MPI.COMM_WORLD, cells, vertices, domain)

#Connectivity array of 5-tetrahedra unitcube, see e.g. https://cecas.clemson.edu/cvel/modeling/EMAG/EMAP/emap4/Meshing.html
UC_vertices=np.array([[1,0,0],[1,1,0],[0,0,0],[0,1,0],[1,0,1],[1,1,1],[0,0,1],[0,1,1]], dtype = float)
UC_cells = np.array([[0,1,3,5],[0,5,6,4],[0,3,2,6],[6,7,5,3],[6,5,0,3]])

#Generate dolfinx mesh for 5-tetrahedra unitcube
mesh = mesh_from_connectivity(UC_cells, UC_vertices)

print("Vertex coordinates of dolfinx mesh:")
print(mesh.geometry.x)
print("Original vertex coordinate array:")
print(UC_vertices)

Which produces the output

Vertex coordinates of dolfinx mesh:
[[1. 0. 0.]
 [1. 1. 0.]
 [0. 1. 0.]
 [1. 1. 1.]
 [0. 0. 1.]
 [0. 1. 1.]
 [0. 0. 0.]
 [1. 0. 1.]]
Original vertex coordinate array:
[[1. 0. 0.]
 [1. 1. 0.]
 [0. 0. 0.]
 [0. 1. 0.]
 [1. 0. 1.]
 [1. 1. 1.]
 [0. 0. 1.]
 [0. 1. 1.]]

Is there a way to preserve the order of the vertices, when one loads a dolfinx-mesh from arrays?

2

No, you cannot enforce that, as

  1. The order of vertices only makes sense in serial.
  2. Data locality is important for performance, and therefore it is reordered to ensure this.

A question that comes to mind:
Why do you need to preserve the vertex/node order?

You can get the mapping by looking at mesh.geometry.input_global_indices,
or: How to apply a boundary condition to a range of DOFs? - #7 by dokken

3

Thanks for the fast answer! mesh.geometry.input_global_indices should be enough to solve the issue for me.

A question that comes to mind:
Why do you need to preserve the vertex/node order?

I am using dolfinx within a larger simulation environment, where global vertex indices are used.