# Evaluation of cell sizes

**URL:** <https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381>\
**Category:** dolfinx\
**Created:** [August 11, 2021, 7:21am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381 "2021-08-11T07:21:19Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Mario-B](https://avatars.discourse-cdn.com/v4/letter/m/e47774/32.png) [@Mario-B](https://fenicsproject.discourse.group/u/Mario-B)\
**Post date:** [August 11, 2021, 7:21am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/1 "2021-08-11T07:21:19Z")

</div>

Hi everybody,

Currently, I translate my ‘classical’ dolfincode to dolfinx. Among others, I face the problem of evaluating all cell sizes by doing:

```auto
elemsize = [c.h() for c in cells(mesh)]

```

But that doesn’t work under dolfinx. As far as I noticed, the ‘hmin’ attribute (see [dolfinx.cpp.mesh — DOLFINX documentation](https://fenicsproject.org/olddocs/dolfinx/dev/python/generated/dolfinx.cpp.mesh.html#dolfinx.cpp.mesh.Mesh)) should give me the option the get the smallest cell size (and via ‘hmax’ the maximal cell size). But unfortunately, it doesn’t work.

I should mention that I read xdmf-meshes into dolfinx that are created by converting msh files into xdmf. My python function doing this reads

```auto
import sys

# converts msh to xdmf
# essential for dolfinx
# see
# https://fenicsproject.discourse.group/t/transitioning-from-mesh-xml-to-mesh-xdmf-from-dolfin-convert-to-meshio/412/147
model_name = sys.argv[1] 
dim = sys.argv[2] # 2d or 3d

import meshio # python3 -m pip install meshio
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

# Read in mesh
msh = meshio.read(f"{model_name}.msh")

# Create and save one file for the mesh, and one file for the facets 
if dim == '2d':
    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write(f"{model_name}_mesh.xdmf", triangle_mesh)
    meshio.write(f"{model_name}_mt.xdmf", line_mesh)    
elif dim == '3d':
    tetra_mesh = create_mesh(msh, "tetra")
    triangle_mesh = create_mesh(msh, "triangle")
    meshio.write(f"{model_name}_mesh.xdmf", tetra_mesh)
    meshio.write(f"{model_name}_mt.xdmf", triangle_mesh)       

```

Has anybody an idea how I can get the whole array of cell sizes for doing some ‘statistics’?

Thanks for the help!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 11, 2021, 8:09am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/2 "2021-08-11T08:09:48Z")

</div>

Consider:

```python
import dolfinx
from mpi4py import MPI

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 5)
tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
h = dolfinx.cpp.mesh.h(mesh, tdim, range(num_cells))
print(h)

```

This code is a simplified version of [dolfinx/test\_mesh.py at 1087856ac3594a1c3c7d99c27a26e65f98da1fa2 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/1087856ac3594a1c3c7d99c27a26e65f98da1fa2/python/test/unit/mesh/test_mesh.py#L278-L284)

---

<div class="post-metadata">

**Author:** ![Mario-B](https://avatars.discourse-cdn.com/v4/letter/m/e47774/32.png) [@Mario-B](https://fenicsproject.discourse.group/u/Mario-B)\
**Post date:** [August 11, 2021, 8:34am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/3 "2021-08-11T08:34:13Z")

</div>

@dokken, thank you very much!

---

<div class="post-metadata">

**Author:** ![hawkspar](https://avatars.discourse-cdn.com/v4/letter/h/b487fb/32.png) [@hawkspar](https://fenicsproject.discourse.group/u/hawkspar)\
**Post date:** [August 29, 2022, 2:47pm UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/4 "2022-08-29T14:47:07Z")

</div>

Thank you for that helpful snippet !

I would like to extend it into a cell-valued `dolfinx.Function` so I could use it as an `Argument` to a `ufl.Form`. My current toy problem goes :

```python
import ufl
import dolfinx as dfx
from mpi4py.MPI import COMM_WORLD as comm

mesh = dfx.UnitSquareMesh(comm, 10, 5)
tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
FE_constant=ufl.FiniteElement("DG",mesh.ufl_cell(),0)
W = dfx.FunctionSpace(mesh,FE_constant)
h = dfx.cpp.mesh.h(mesh, tdim, range(num_cells))
hf = dfx.Function(W)
hf.x.array[:]=h

```

This works in serial but fails in parallel. It’s troubling to me because you seem to have done some effort in writing a parallel snippet and using `Function` then `x.array` directly in parallel usually works well for me. Traceback goes :

```bash
Traceback (most recent call last):
  File "/home/shared/src/sanity_check.py", line 13, in <module>
    hf.x.array[:]=h
ValueError: could not broadcast input array from shape (10,) into shape (17,)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 29, 2022, 3:16pm UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/5 "2022-08-29T15:16:08Z")

</div>

> [@hawkspar](#):
>
> ```auto
> num_cells = mesh.topology.index_map(tdim).size_local
> 
> ```

Change this to

```python
num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts

```

---

<div class="post-metadata">

**Author:** ![hawkspar](https://avatars.discourse-cdn.com/v4/letter/h/b487fb/32.png) [@hawkspar](https://fenicsproject.discourse.group/u/hawkspar)\
**Post date:** [August 30, 2022, 9:25am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/6 "2022-08-30T09:25:18Z")

</div>

Thank you ! It works beautifully.

---

<div class="post-metadata">

**Author:** ![mcjbo](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mcjbo/32/7655_2.png) [@mcjbo](https://fenicsproject.discourse.group/u/mcjbo)\
**Post date:** [June 12, 2025, 8:37am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/7 "2025-06-12T08:37:21Z")

</div>

Hi dokken,  
apparently, this does not work in dolfinx 0.9.0 causing

```auto
TypeError: h(): incompatible function arguments. The following argument types are supported:
    1. h(mesh: dolfinx.cpp.mesh.Mesh_float32, dim: int, entities: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> numpy.ndarray[dtype=float32,]
    2. h(mesh: dolfinx.cpp.mesh.Mesh_float64, dim: int, entities: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> numpy.ndarray[dtype=float64,]

Invoked with types: dolfinx.mesh.Mesh, int, range

```

In general I was wondering that `range(num_cells)` will always start from 0 on each process. However, index 0 should be on one process only if executed in parallel.

I have worked on a different solution, which in turn results in PETSC ERROR.

```auto
import dolfinx
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 5)

tdim = mesh.topology.dim
local_range = mesh.topology.index_map(tdim).local_range
h = mesh.h(tdim, np.arange(local_range[0],local_range[1]))
print(f"rank {MPI.COMM_WORLD.Get_rank()} h = \n{h}")

```

It works in serial, but fails in parallel on all ranks \>0:

```auto
[1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[1]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[1]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[1]PETSC ERROR: to get more information on the crash.
[1]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 1 (rank 1 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 1

```

Another option I tried is using

```auto
h = mesh.h(tdim, mesh.topology.original_cell_index)

```

causing the same error on rank 0.

I am using dolfinx 0.9.0 with petsc 3.22.0 and petsc4py 3.22.0 from conda-forge.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 12, 2025, 9:08am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/8 "2025-06-12T09:08:34Z")

</div>

> [@mcjbo](#):
>
> In general I was wondering that `range(num_cells)` will always start from 0 on each process. However, index 0 should be on one process only if executed in parallel.

`mesh.h` takes in the local cell numbering (within a process), not the global cell numbering.  
See for instance Figure 8 in [https://zenodo.org/records/10447666:](https://zenodo.org/records/10447666:)

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/99cebea86cf2a72130a7347e640307791fda3f12.png)  
To fix the first error you had, change

> [@hawkspar](#):
>
> `range(num_cells)`

`np.arange(num_cells, dtype=np.int32)`

---

<div class="post-metadata">

**Author:** ![mcjbo](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mcjbo/32/7655_2.png) [@mcjbo](https://fenicsproject.discourse.group/u/mcjbo)\
**Post date:** [June 12, 2025, 9:39am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/9 "2025-06-12T09:39:37Z")

</div>

> [@dokken](#):
>
> `mesh.h` takes in the local cell numbering (within a process), not the global cell numbering.  
> See for instance Figure 8 in [https://zenodo.org/records/10447666:](https://zenodo.org/records/10447666:)

Perfect, that’s what I needed. However,

```auto
h = dolfinx.cpp.mesh.h(mesh, tdim, np.arange(num_cells, dtype=np.int32))

```

still produces `TypeError`. Only calling `h` from the mesh directly works for me:

```auto
h = mesh.h(tdim, np.arange(num_cells,dtype=np.int32))

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 12, 2025, 9:50am UTC](https://fenicsproject.discourse.group/t/evaluation-of-cell-sizes/6381/10 "2025-06-12T09:50:30Z")

</div>

> [@mcjbo](#):
>
> `= dolfinx.cpp.mesh.h(mesh, tdim, np.arange(num_cells, dtype=np.int32))`

That is because this is a C++ binding, and expecting `mesh._cpp_object` rather than `mesh`.

> [@mcjbo](#):
>
> `h = mesh.h(tdim, np.arange(num_cells,dtype=np.int32))`

This is the preferred wrapper that takes care of the handling of the C++ objects