Floating point exception with cpp.mesh.entities_to_geometry

dolfinx
,
1

Hi,
I am experiencing some issues when using cpp.mesh.entities_to_geometry in dolfinx 0.9.0 (installed via conda install -c conda-forge adios4dolfinx fenics-dolfinx). It was working without problems in dolfinx 0.7.3.
In the following a MWE:

from mpi4py import MPI
from dolfinx import fem, mesh, cpp
import numpy as np

comm = MPI.COMM_WORLD

L = 7.0
mesh_size = 0.4
N = int(round(L/mesh_size))
if N % 2 == 0:
    N += 1
volume_tag = 1
surface_tag = 10


domain = mesh.create_box(comm=comm, points=[[0.0, 0.0, 0.0], [L, L, L]],n=[N, N, N], cell_type=mesh.CellType.tetrahedron)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
domain.topology.create_connectivity(tdim, 0)
domain.topology.create_connectivity(fdim, 0)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
facet_tags = mesh.meshtags(domain, fdim, mesh.exterior_facet_indices(domain.topology), surface_tag)
cell_tags = mesh.meshtags(domain, tdim, np.array(range(domain.topology.index_map(tdim).size_local)), volume_tag)

V_test = fem.functionspace(domain, ("CG",1))
print(f"Rank {comm.rank}: Function space created successfully.", flush=True)

n_nodes = domain.geometry.index_map().size_local
bnd_facets = facet_tags.find(surface_tag)
local_bnd_facets = [f for f in bnd_facets if f<domain.topology.index_map(fdim).size_local]
local_bnd_facets_array = np.array(local_bnd_facets, dtype=np.int32)
local_bnd_facets_array.flags.writeable = False       
domain.topology.create_entity_permutations()


bnd_f_to_n = cpp.mesh.entities_to_geometry(domain._cpp_object, fdim, local_bnd_facets_array, True)

print(f"Rank {comm.rank}: here.", flush=True)

bnd_nodes = list(set(bnd_f_to_n.reshape(-1)))
        

if comm.rank ==0:
    print("END")

I create a box and perform some operations on the mesh. Observed behavior:

  • The code runs fast and without failing when requesting up to 22 cores
  • If I increase the number of cores (30, for example), the code runs very slowly and eventually fails when using cpp.mesh.entities_to_geometry with a floating point exception

[13]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

I added the following checks on local_bnd_facets_array and it looks like:

Rank 13: local_bnd_facets_array =

I assume this happens because in parallel, not all ranks own boundary facets.
However, shouldn’t cpp.mesh.entities_to_geometry be able to handle this gracefully as I assume it was the case in previous versions of dolfinx? Or is there a reason why this is not the case?
Also, the fact that the code runs very slowly before failing is an indication that something is already corrupted ahead of this function?

Please note that I need to use cpp.mesh.entities_to_geometry (not the corresponding function from topology) to map the entities to their geometric representation. In the sequel of my code, facets must be oriented.

Many thanks for your help!

# 1. Check type
if not isinstance(local_bnd_facets_array, np.ndarray):
    raise TypeError("Input must be a NumPy array.")
if local_bnd_facets_array.dtype != np.int32:
    raise TypeError(f"Array must be of type np.int32, got {local_bnd_facets_array.dtype}.")
if local_bnd_facets_array.ndim != 1:
    raise ValueError("Input array must be 1-dimensional.")

# 2. Check index range
size_local = domain.topology.index_map(fdim).size_local
if not np.all((local_bnd_facets_array >= 0) & (local_bnd_facets_array < size_local)):
    bad = local_bnd_facets_array[
        (local_bnd_facets_array < 0) | (local_bnd_facets_array >= size_local)
    ]
    raise ValueError(
        f"Rank {comm.rank}: Found invalid facet indices {bad}. "
        f"Valid range is 0 to {size_local-1}."
    )

# 3. Check for duplicates (optional, warning only)
if len(local_bnd_facets_array) != len(np.unique(local_bnd_facets_array)):
    print(f"Rank {comm.rank}: Warning: Duplicate facet indices detected.", flush=True)

# 4. Print for debugging
print(f"Rank {comm.rank}: local_bnd_facets_array = {local_bnd_facets_array}", flush=True)
print(f"Rank {comm.rank}: Valid index range is 0 to {size_local-1}", flush=True)

print(f"Rank {comm.rank}: Facet array validation passed.", flush=True)
2

This was fixed in:

and will be part of the 0.10 release.