Identifying cells and their vertices

mesh
1

Hi all,
I have a code in dolfinx that solves a system of PDEs in a rather complex geometry that I mesh using gmsh. The result is a vector field \mathbf m. I would now need to calculate some quantities for each cell of the mesh, for which I need to identify all the vertices that belong to them. Graphically, it would be something like this:

0 -------- 1 -------- 4 
|          |          |
|  cell 0  |  cell 1  |
|          |          |
2 -------- 3 -------- 5

I would need to know which vertices (in this case, 0, 1, 2, 3) belong to cell 0, calculate some stuff with the value of \mathbf m at them, then do the same for cell 1 (with vertices 1, 4, 5, 3), etc.
Can anyone give me a hint on how to achieve this?
Thanks in advance!

Here’s a MWE on how I create the mesh:

import gmsh
model = gmsh.model; geom = gmsh.model.geo

from mpi4py import MPI
mesh_comm = MPI.COMM_WORLD

ms = 0.1
Rbase = 1
gdim = 2

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
model.add("model_acorn")

p = [
    geom.addPoint(Rbase, np.pi, 0, ms),
    geom.addPoint(Rbase, np.pi/2, 0, ms),
    geom.addPoint(0, np.pi/2, 0, ms),
    geom.addPoint(0, np.pi, 0, ms)
]
l = [
    geom.addLine(p[0], p[1]),
    geom.addLine(p[1], p[2]),
    geom.addLine(p[2], p[3]),
    geom.addLine(p[3], p[0])
]

cl = geom.addCurveLoop([l[j] for j in range(len(l))])
s = geom.addPlaneSurface([cl]) 

acorn = model.addPhysicalGroup(gdim, [s]) 
model.setPhysicalName(gdim, acorn, "Flat acorn")

base = model.addPhysicalGroup(gdim-1, [l[0]])
model.setPhysicalName(gdim-1, base, "Base")

cap = model.addPhysicalGroup(gdim-1, [l[1]])
model.setPhysicalName(gdim-1, cap, "Cap")

geom.synchronize()
model.mesh.generate(2)

domain, cell_markers, facet_markers = gmshio.model_to_mesh(model, mesh_comm, 0, gdim=gdim)
2 
import dolfinx
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2, cell_type=dolfinx.mesh.CellType.quadrilateral)
c2v = mesh.topology.connectivity(mesh.topology.dim, 0)
print(c2v)

  0: [0 1 2 3 ]
  1: [1 4 3 5 ]
  2: [2 3 6 7 ]
  3: [3 5 7 8 ]

image

1 Like
3

Thanks a lot, @nate ! Could you please give me some more details on how to access these values and how to evaluate a fem.Function at them?
Let’s say we have the mesh domain that I wrote in my previos post, and there I define

element = ufl.VectorElement("CG", domain.ufl_cell(), 1, 3)
V = fem.FunctionSpace(domain, element)
m = fem.Function(V)

I then solve a non-linear problem (I don’t think this is relevant, bust just in case) to find the value for m at domain. So far so good. But now I want to sum the values of m at all the vertices of a given cell i. How can I evaluate something like m[cell_i[vertex_j]]?
Thanks!

4

See: Application of point forces, mapping vertex indices to corresponding DOFs - #2 by dokken

2 Likes
5

Hi @dokken, thanks for your reply. I have looked into the post you linked and I think I know now how to access the vertices of a cell and operate with them. However, I don’t really understand the difference between vertices and dofs. I find that vertices = c_to_v.links(cell) and dofs = V0.dofmap.cell_dofs(cell) give me identical arrays. Is this always the case?

6

No, that is not always the case.
A vertex is part of the mesh. Every cell has N number of vertices.

A degree of freedom does not necessarily belong to a vertex (this is only the case for first order Lagrange elements). See for instance: DefElement

1 Like
7

Hi again, @nate and @dokken. I realise now that I forgot to say that I actually need to access the nodes of a given cell in order, forming a closed loop (clockwise or anticlockwise is the same as long as I know which one it is). So in my first example,

0 -------- 1 -------- 4 
|          |          |
|  cell 0  |  cell 1  |
|          |          |
2 -------- 3 -------- 5

I would need to operate with the value of the fem.Function m at nodes 0 2 3 1 0 (in that order) of cell 0, then at nodes 1 3 5 4 1 of cell 1, etc. Is there any way to achieve this? As I see from nate’s response, mesh.topology.connectivity(mesh.topology.dim, 0) doesn’t necessarily return the nodes in order.
Thank you for your help!

8

I don’t understand why the way you are traversing the cell matters. You could always compute the orientation by computing the normal vector based on the geometry. See:

You should do something similar to this depending on your application.