Mesh connectivity in Fenicsx

dolfinx
1

Hello everyone!
I used to use FreeFem++. But considering the advantages of coding with Python, I decided to migrate to Fenicsx. However, it seems that there is no comprehensive tutorial for Fenicx and it should be learned mostly from examples and demos.
I have questions about working with meshes, some of which I got the answers from the examples, but I don’t know the best way. Thank you for your advice.

Having domain mesh, I want to get:
1- The total number of elements in each subdomain
2- The total number of nodes in each subdomain
3- The relationship between the indices of nodes and the indices of dofs
4- The element including the point x0 and y0
5- The area (volume) of the i-th element
6- the indices of the boundary nodes

My solutions are:

1- domain.geometry.x.shape[0]
2- domain.topology.connectivity(2,0).array.size/3
3- ?
4- ?
5- getting the coordinates of vertices and then using the corresponding formula for calculating area (volume).
6- ?

3 Likes
2

There are alot of more detailed examples at: The FEniCSx tutorial — FEniCSx tutorial

domain.topology.index_map(dim) gives you a map of indices for entities of dimension dim. So say you have a 2D mesh, mesh.topology.index_map(2) will give you a map for all triangles (or quads) in the domain. You can access size_local (number of triangles owned by the process), size_global (total number of triangles) across all mpi ranks, num_ghosts (number of triangles on the process that is owned by another rank, typical at partition boundaries).

This is highly dependent on your function space. For first order Lagrange spaces, you can get a relationship between the dofs and the mesh vertices. However, for spaces such as higher order Lagrange, Nedelec
or RT elements, dofs are associated with edges/facets and not vertices. You can get more information about the layout from the FunctionSpace dofmap: dolfinx.fem — DOLFINx 0.8.0.0 documentation
which is documented in the C++ API DOLFINx: ElementDofLayout Class Reference
and relates to the indices in
dolfinx.fem — DOLFINx 0.8.0.0 documentation

See: Implementation — FEniCSx tutorial

This is highly dependent on what cell type (and order you are using).
You could use ufl.CellVolume(domain) and interpolate this expression into a DG-0 function similar to: Implementation — FEniCSx tutorial

As I’m not sure if you mean vertices or degrees of freedom here;

  1. Vertices: dolfinx.mesh — DOLFINx 0.8.0.0 documentation
  2. Could use 1. And then: dolfinx.fem — DOLFINx 0.8.0.0 documentation
5 Likes
3

Thank you so much @dokken! this really helped me.

4

Hi dokken,

I am trying to calculate cell volume using ufl.Cell_Volume

However, I got a RuntimeError: Not handled: <class ‘ufl.geometry.CellVolume’>

It also shows that “Only know how to compute the cell volume of an affine cell.”

The original codes are here:

domain = mesh.create_rectangle( MPI.COMM_WORLD,
np.array([[0,0],[80,80]]),
[80, 50],
cell_type=mesh.CellType.quadrilateral)
Ve = fem.FunctionSpace(domain, (“DG”, 0))
elmVol_ufl = ufl.CellVolume(domain)
elmVol_expr = fem.Expression(elmVol_ufl, np.array([4131,2]))
elmVol = fem.Function(Ve,name=‘Element Volume’)
elmVol.interpolate(elmVol_expr)

Is there something wrong with my code?

5

I’ve made an issue regarding this at:

A workaround would be to use:

Ve = fem.FunctionSpace(domain, ("DG", 0))
q_degree = 3
vol2 = fem.Function(Ve)
fem.petsc.assemble_vector(vol2.vector, fem.form(
    1*ufl.TestFunction(Ve)*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)"
1 Like
6

Hi dokken, I have some questions regarding this.

What is the meaning of the following command,i.e., what is printed out by this?

print(vol2.x.array - elmVol.x.array)

Where is the calculated volume stored? Is that in vol2? It seems to me that vol2 is defined as a Function. Is it possible to store the volume of all the cells in a np.array? If so, how can I achieve that?

Thank you!

7

The volume is stored in vol2, and you can Get the array by accessing vol2.x.array

8

I noticed this requires dolfinx. Now I am working with dolfin. Does dolfin have the similar functions?

9

Yes, you can do the same in DOLFIN.

from dolfin import *

mesh = UnitSquareMesh(10, 10)

for cell in cells(mesh):
    print(cell.volume())

as shown in: