Wrong order of dofs after creating mesh from lists of points and cells

[dokken]

You are making some major assumptions here that simply cannot be held by DOLFINx

@mnotazio

I need the order of dofs (necessarily for a Function Space of type CG1) of the created mesh to match with the order of the specified points.

This is because:

1. Only simple function spaces has the same number of dofs as the mesh geometry
2. The input ordering does not make sense once one wants to run in parallel.

More information about this has been discussed in:

I/O from XDMF/HDF5 files in dolfin-x

I want to share what is now working much faster for me than my old approach via the coordinates. In preprocessing (serial), read in your mesh with the standard XDMF routine, set values to a function, and write a text file using igi = msh.geometry.input_global_indices with the following structure igi1 fx1 fy1 fz1 igi2 fx2 fy2 fz2 ... ign fxn fyn fzn so the input node index followed by the x,y,z values of your vector field (or the one value of a scalar field). Then, in parallel, use def rea…

Consider the following MWE to map a vertex->mesh node->original_input_position

import dolfinx
import matplotlib.pyplot as plt
import numpy as np
import ufl
from mpi4py import MPI

pp = np.array([[0.5, 1], [2, 1], [3, 1.5], [3.5, 2.5], [2.2, 2], [1, 2.2]])
ee = np.array([[0, 1, 5], [1, 4, 5], [1, 2, 4], [2, 3, 4]])
plt.triplot(pp[:, 0], pp[:, 1], ee)
plt.show()

gdim = 2
shape = "triangle"
degree = 1

cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))

cells = np.array(ee, dtype=np.int32)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, pp, domain)

V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
dofs = V.tabulate_dof_coordinates()[:, 0:2]

geom_loc_to_org = mesh.geometry.input_global_indices
vertices = np.arange(mesh.topology.index_map(0).size_local, dtype=np.int32)
top_to_geom = dolfinx.cpp.mesh.entities_to_geometry(mesh, 0, vertices, False).reshape(-1)


for ii in range(len(vertices)):
    geom_node= top_to_geom[vertices[ii]]
    original_pos = geom_loc_to_org[top_to_geom[ii]]
    assert np.allclose(pp[original_pos,:], mesh.geometry.x[geom_node, :mesh.geometry.dim])
    print(f"{pp[original_pos,:]=}, {mesh.geometry.x[geom_node, :]=}")

The reason for
vertex->mesh node is that DOLFINx supports isoparametric mesh elements (i.e. you can have higher order cells that have more nodes than just vertices).