Hi everyone,
I was wondering if in dolfinx exists the possibility of computing the size of the elements of a mesh.
I have found that a characteristic size can be found as in Evaluation of cell sizes - #2 by dokken, however I was looking for something a bit different. Assuming that the mesh is composed by triangles, is it possible to find the area of each triangle?
Thank you in advance
Assuming that the mesh is consting of linear triangles, you can get this with the following code:
import dolfinx
from mpi4py import MPI
import ufl
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 1)
DG0 = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
v = ufl.TestFunction(DG0)
cell_area_form = dolfinx.fem.form(v*ufl.dx)
cell_area = dolfinx.fem.assemble_vector(cell_area_form)
print(cell_area.array)
Thank you for the quick answer. A similar procedure should holds also for 3D meshes made of tetrahedron for instance?
Yes, you’re just computing \int_{\kappa} 1 \; \mathrm{d}\vec{x}, for each element \kappa in the mesh. This is particularly useful for high order elements, for example.
Thank you very much! 
Hi there,
I’m trying to use this method to find the area of a small domain defined in GMSH (with its own defined physical group)
I tried a simpler example for a simple 2D mesh. And
len(cell_area.array)
is different from the number of triangular elements on the mesh. So I cannot index to the cell tags.
Any ideas on how to solve this?
Many thanks,
Katie
Please provide a minimal reproducible example.
Hi Dokken, here is a minimal reproducible example: the msh file (copied at the bottom) creates this mesh, which has 26 triangular elements
The code is
import gmsh
from mpi4py import MPI
import meshio
from dolfinx import mesh, io, fem
from ufl import TrialFunction, TestFunction, dx
import meshio
import numpy as np
gmsh.initialize()
current_mesh = "mesh_test_1"
gmsh.open(current_mesh + ".msh")
model = gmsh.model
model.setCurrent(current_mesh)
print('Model ' + model.getCurrent() + ' (' +str(model.getDimension()) + 'D)')
gmshiomsh, cell_markers, facet_markers = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2) # Extract the data
def create_mesh(mesh, cell_type, prune_z=False): # These are all functions for a meshio mesh
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
print(f'There are {len(cell_data)} elements of type {cell_type}')
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
# Read in mesh
msh = meshio.read("mesh_test_1.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mt.xdmf", line_mesh)
with io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
mymesh = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mymesh, name="Grid")
mymesh.topology.create_connectivity(mymesh.topology.dim, mymesh.topology.dim-1)
with io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(mymesh, name="Grid")
V = fem.FunctionSpace(mymesh, ("CG", 1))
v = TestFunction(V)
cell_area_form = fem.form(v*dx)
cell_area = fem.assemble_vector(cell_area_form)
print(f'Cell area array has length {len(cell_area.array)}')
and the msh file is
$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
4 4 1 0
1 0 0 0 0
2 1 0 0 0
3 1 1 0 0
4 0 1 0 0
1 0 0 0 1 0 0 1 1 2 1 -2
2 1 0 0 1 1 0 0 2 2 -3
3 0 1 0 1 1 0 0 2 3 -4
4 0 0 0 0 1 0 0 2 4 -1
1 0 0 0 1 1 0 1 1 4 1 2 3 4
$EndEntities
$Nodes
9 20 1 20
0 1 0 1
1
0 0 0
0 2 0 1
2
1 0 0
0 3 0 1
3
1 1 0
0 4 0 1
4
0 1 0
1 1 0 2
5
6
0.3333333333325025 0 0
0.6666666666657889 0 0
1 2 0 2
7
8
1 0.3333333333325025 0
1 0.6666666666657889 0
1 3 0 2
9
10
0.6666666666675911 1 0
0.3333333333347203 1 0
1 4 0 2
11
12
0 0.6666666666675911 0
0 0.3333333333347203 0
2 1 0 8
13
14
15
16
17
18
19
20
0.7113248654055673 0.4999999999991457 0
0.2886751345942123 0.5000000000011557 0
0.5000000000006117 0.7525600817161773 0
0.4999999999993867 0.2474399182839603 0
0.2423197548524782 0.7576802451481532 0
0.757680245147464 0.2423197548520695 0
0.2423197548507857 0.2423197548513912 0
0.7576802451491019 0.7576802451486099 0
$EndNodes
$Elements
2 29 1 29
1 1 1 3
1 1 5
2 5 6
3 6 2
2 1 2 26
4 13 14 16
5 14 13 15
6 7 8 13
7 11 12 14
8 5 6 16
9 9 10 15
10 10 4 17
11 6 2 18
12 2 7 18
13 4 11 17
14 12 1 19
15 8 3 20
16 1 5 19
17 3 9 20
18 14 12 19
19 13 8 20
20 7 13 18
21 11 14 17
22 16 14 19
23 15 13 20
24 14 15 17
25 13 16 18
26 16 6 18
27 15 10 17
28 5 16 19
29 9 15 20
$EndElements
Many thanks for your time,
Katie
You are using a CG 1 (or better known as a first order lagrange space). This corresponds to each vertex of your mesh, not each cell. For having a value per cell, you need to change the function space to ("DG", 0).
Thank you for clarifying, that makes much more sense!
Many thanks,
Katie