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

[mnotazio]

Hello everyone! I just started using dolfinx recently, so I apologize in advance if the question turns out to be silly!

I am using dolfinx version 0.5.2 and am creating a mesh from a list of points and a cell list specifying the connection between the points.
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.

Taking the simple test as an example (from Create mesh in fenicsx from point and element array):

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

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]

for ii in range(pp.shape[0]):
print(f “pp[{ii}, :] = {pp[ii, :]}\t\tdofs[{ii}, :] = {dofs[ii, :]}”
f “mesh.geometry.x[{ii}, :] = {mesh.geometry.x[ii, :]}”)

The sorting of dofs results differently than the sorting of pp points.

Could someone help me with this? Thanks in advance!

[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).

[mnotazio]

Thank you so much for your answer! It works fine!

[ivanmaione]

I wonder if there is a similar approach to find the order of the mesh cells with respect to the initial order given in ee. Indeed, the final assertion fails when running the following code (I am aware that the code is valid only for CG1).
I am asking because I have a similar problem in which, in addition to the connectivity information (i.e. pp and ee), I have also a marker array that I would like to use as a MeshTag for further operations. Unfortunately, due to the change in the cells’ order, I am not able to directly pass the marker array in the meshtag_from_entities function.

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([[1, 4, 5], [1, 2, 4], [0, 1, 5], [2, 3, 4]])
markers = np.array([1,1,2,3])

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, :]=}")

cells_conn = mesh.topology.connectivity(mesh.topology.dim,0)

# The following line will not produce the expected cell meshtag since markers are created using ee order
ct = dolfinx.mesh.meshtags_from_entities(mesh, mesh.topology.dim, cells_conn, markers)

conn1 = []
for ii in range(len(cells_conn.array)):
    conn1.append(geom_loc_to_org[top_to_geom[cells_conn.array[ii]]])

conn1 = np.array(conn1).reshape(-1,3)
for ii in range(len(conn1)):
    assert np.allclose(conn1[ii], ee[ii])

I am also using dolfinx version 0.5.2.

[dokken]

@ivanmaione

I wonder if there is a similar approach to find the order of the mesh cells with respect to the initial order given in ee.

You can use:

o_cell_idx =  mesh.topology.original_cell_index
for ii in range(len(conn1)):
    assert np.allclose(conn1[ii], ee[o_cell_idx[ii]])

[ivanmaione]

Works like a charm!

Thanks a lot, @dokken.

[mnotazio]

Taking the example from @ivanmaione:

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([[1, 4, 5], [1, 2, 4], [0, 1, 5], [2, 3, 4]])
markers = np.array([1,1,2,3])

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, :]=}")

cells_conn = mesh.topology.connectivity(mesh.topology.dim,0)
conn1 = []
for ii in range(len(cells_conn.array)):
    conn1.append(geom_loc_to_org[top_to_geom[cells_conn.array[ii]]])

o_cell_idx =  mesh.topology.original_cell_index
for ii in range(len(conn1)):
    assert np.allclose(conn1[ii], ee[o_cell_idx[ii]])

ct = dolfinx.mesh.meshtags_from_entities(
        mesh,
        mesh.topology.dim,
        cells_conn,
        markers[o_cell_idx].astype(np.int32),
    )

In order to have a complete import from any mesh, given arrays of nodes, elements, and cell markers, a dolfinx mesh would also need facet tags (at the moment ignored). Would it be possible to define facet tags based on cell tags to distinguish the various domains within the mesh?

[dokken]

To me it is unclear what you mean by

@mnotazio

Would it be possible to define facet tags based on cell tags to distinguish the various domains within the mesh?

Do you mean that you want to make facet tags based on the cell domains (if so, what should be the criterion for each marker), or do you mean make facet markers based on the vertex numbers, say a cell with vertices [0,1,2] has the facets [0,1], [1,2],[2,0].

[mnotazio]

Yes, I would like to make facet tags based on the cell domains in order to define the facets of the cells which define the edge of each subdomain. The facet marker would have the same value of cell marker for a certain cell.

Something similar has been discussed here:

Get boundary points from custom cell tags

Sorry, I’m not expert about it so I’ll try to explain my self better. I need the two vertexes for each cell marked with 1 (in yellow) which are part of the boundary for subdomain 1 (In order to define the boundary, for example, of a more complex domain) → (the red points in the figure) [a]

(Referring to the picture of the linked question, I would like to identify the facets “corresponding” to the points with the red cross. If the yellow subdomain is defined by cell tags with value 1, the drawn facets would have value 1).

The question is: is there a sort of built-in method in dolfinx or do I have to apply the same logic for each subdomain? In this case I need facets in order to use locate_dofs_topological method to identify dofs on the edge of a certain subdomain.

I hope I have clarified the concept. Thank you for your help!

[dokken]

@mnotazio

The facet marker would have the same value of cell marker for a certain cell.

We again then come to the issue: What if a facet is connected to two cells, and each of the cells have different markers? You have not explained what should happen if the domain marked with one has an exterior facet.

Here is a simple script to illustrate how to uniquely identify all facets of a mesh.

import dolfinx
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

def some_cells(x):
    return (x[0]<0.2) | ((x[0]>0.5) & (x[0]<0.8) & (x[1]>0.5) & (x[1]<0.8))
all_cells = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)
values = np.zeros_like(all_cells, dtype=np.int32)
values[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, some_cells)] = 1

ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, all_cells, values)
with dolfinx.io.XDMFFile(mesh.comm, "ct.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ct, mesh.geometry)


mesh.topology.create_entities(mesh.topology.dim-1)
all_facets = np.arange(mesh.topology.index_map(mesh.topology.dim-1).size_local, dtype=np.int32)
facet_values = np.zeros_like(all_facets)


zero_facets = dolfinx.mesh.compute_incident_entities(mesh.topology, ct.find(0), mesh.topology.dim, mesh.topology.dim-1)

ones_facets = dolfinx.mesh.compute_incident_entities(mesh.topology, ct.find(1), mesh.topology.dim, mesh.topology.dim-1)

interface_facets = np.intersect1d(zero_facets, ones_facets)

facet_values[zero_facets] = 5
facet_values[ones_facets] = 10
facet_values[interface_facets] = 15
ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, all_facets, facet_values)

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with dolfinx.io.XDMFFile(mesh.comm, "ft.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ft, mesh.geometry)

[mnotazio]

Thank you @dokken, it works perfectly for my purpose!