Evaluation of cell sizes

dolfinx
1

Hi everybody,

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

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

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!

2

Consider:

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

3 Likes
3

@dokken, thank you very much!

4

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 :

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 :

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,)
5

Change this to

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

Thank you ! It works beautifully.