Why does tabulate_dof_coordinates give a different result than mesh.coordinates()?

General
1

I use pygalmesh to create a mesh from a binary segmentation mask. I am trying to run an FEM simulation on the mesh but I was getting nans in some of the output. After digging deeper, I noticed some of the dof coordinates are zero, which they are not in the mesh coordinates.

I am using a linear function space (“P”, 1), so I expect the set of dof coordinates to be the same as the set of mesh vertices. However, this is not the case in the code below. Some of the mesh coordinates turn into zero coordinates in the output of tabulate_dof_coordinates.

Minimum working example:

import nibabel as nib
import pygalmesh
import meshio
from mpi4py import MPI
import fenics as fe

def convert_to_fe_mesh(points, cells):
    '''
    Convert from meshio to fenics mesh.
    '''
    mesh_file = 'temp.xdmf'
    meshio.write_points_cells(mesh_file, points, [(cells.type, cells.data)])
    fe_mesh = fe.Mesh()
    with fe.XDMFFile(MPI.COMM_WORLD, mesh_file) as f:
        f.read(fe_mesh)
    return fe_mesh

nifti = nib.load('lung_combined_mask.nii.gz')
max_radius = 5.0 # try 10.0, 20.0

mesh = pygalmesh.generate_from_array(
    nifti.get_fdata().astype(np.uint16),
    voxel_size=nifti.header.get_zooms(),
    max_cell_circumradius=max_radius,
    odt=True
)
fe_mesh = convert_to_fe_mesh(mesh.points, mesh.cells[1])
fe_space = fe.FunctionSpace(fe_mesh, 'P', 1)

# compare mesh coordinates to dof coordinates for linear FEM basis
v1 = [tuple(v) for v in mesh.points]
v2 = [tuple(v) for v in fe_mesh.coordinates()]
v3 = [tuple(v) for v in fe_space.tabulate_dof_coordinates()]

assert set(v1) == set(v2), (set(v1) - set(v2), set(v2) - set(v1)) # these are the same
assert set(v1) == set(v3), (set(v1) - set(v3), set(v3) - set(v1)) # these are not

This is the output:

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[268], line 25
     22 v3 = [tuple(v) for v in fe_space.tabulate_dof_coordinates()]
     24 assert set(v1) == set(v2), (set(v1) - set(v2), set(v2) - set(v1))
---> 25 assert set(v1) == set(v3), (set(v1) - set(v3), set(v3) - set(v1))

AssertionError: ({(169.13596, 130.50793, 130.2164), (97.680405, 161.96782, 28.77537)}, {(0.0, 0.0, 0.0)})

I noticed that the number of coordinates that turn to zero depends on the mesh size (controlled by max_cell_circumradius). If I use a lower resolution mesh, I can get the coordinate sets to be the same. If I increase the resolution, I get one or two coordinates at the end of the array returned by tabulate_dof_coordinates() set to all zeros, which causes nan errors downstream.

How should I understand this behavior?

2

Here is a google drive link to the NIFTI file containing the binary mark:

I am using fenics version 2019.1.0.

3

If you look at the number of points in your mesh:

mesh.points.shape

it yields:

(14010, 3)

while if you look at the number of unique points used in the cells are:

len(np.unique(mesh.cells[1].data.reshape(-1)))

yields

14009

We can then find which node in the mesh that is not part of a cell:

unique_in_cells = np.unique(mesh.cells[1].data.reshape(-1))
num_points = np.arange(mesh.points.shape[0])
print(np.flatnonzero(np.invert(np.isin(num_points, unique_in_cells))))

yielding

array([7853])

For further questions, please note that it would be beneficial to add links to install instructions for cgal to your post, as it is required to reproduce the results.

4

Thank you for the rapid response, this solved my issue.