Volume of an element

mesh
1

Is there a way to find the volume of each element / of a given element in the domain, in DOLFINx? I would like to do this for a tetrahedral mesh of ~100k elements.

I’ve seen the topic Mesh connectivity in Fenicsx - #2 by dokken , but the suggested approach (see MWE) doesn’t seem to run:

from mpi4py import MPI
from dolfinx import mesh, fem
import ufl

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

V = fem.FunctionSpace(domain,('DG',0))
fun = fem.Function(V)
vol = ufl.CellVolume(domain)
fun.interpolate(vol)
2

If you read a bit further in the post, there is an approach that works around this issue: Mesh connectivity in Fenicsx - #5 by dokken

3

I’ve tried the code snippet in the post you reference, modified as my actual domain is 3D:

V = fem.FunctionSpace(domain,('V',0))
elmVol_ufl = ufl.CellVolume(domain)
elmVol = fem.Function(V,name='Element Volume')
num_nodes = len(V.dofmap.list)
elmVol_expr = fem.Expression(elmVol_ufl, np.array([num_nodes,3]))

q_degree = 3
vol2 = fem.Function(V)
fem.petsc.assemble_vector(vol2.vector, fem.form(
    1*ufl.TestFunction(V)*ufl.dx(metadata={"quadrature_degree": q_degree})))
vol2.x.scatter_forward()
print(vol2.x.array - elmVol.x.array)
assert np.allclose(vol2.x.array, elmVol.x.array)
vol2.name = "ElementVolume (assembly)"

but get an assertion error at the fem.Expression line. It’s hard to debug this as I don’t understand the code snippet.

I’m not sure why I can’t define this fem.Expression, and I don’t understand what vol2 represents.

4

That is because it is not supoprted at the moment, and the syntax:

is not correct.
Why don’t you use:

Ve = fem.FunctionSpace(domain, ("DG", 0))
q_degree = 3
vol = fem.Function(Ve)
fem.petsc.assemble_vector(vol.vector, fem.form(
    1*ufl.TestFunction(Ve)*ufl.dx(metadata={"quadrature_degree": q_degree})))
vol.x.scatter_forward()

which amounts to assembling 1 over each element, i.e. getting the volume of it with linear complexity.

2 Likes