How to get a numpy array with tetrahedrons from a mesh

Question

When I build a domain $\Omega$ in dolfin:

from dolfin import *
import numpy as np
omega = UnitCubeMesh(10,10,10)

after this I need a numpy array with the elements (tetrahedra)

elements = # ... I have not idea what to put here

and I want to do this

>>>elements[6]
[1, 24, 49, 56]

therefore, the tetrahedron with the index 6 have the points with the index 1, 24, 49 and 56 this index is over the points array.

Answer 1

Hi, it sounds like you want a dictionary with cell.index() as a key and vertices of the cell
as values. You can build one by

elements = dict((cell.index(), cell.entities(0)) for cell in cells(mesh))

Just for illustration, here is how you'd get physical coordinates of vertices that form cells
in your mesh using the elements dictionary

from dolfin import *

mesh = UnitCubeMesh(2, 2, 2)
coordinates = mesh.coordinates()
elements = dict((cell.index(), cell.entities(0)) for cell in cells(mesh))

for element in range(mesh.num_cells()):
    print "Vertices of element", element
    for vertex in elements[element]:
        print "\t", vertex, coordinates[vertex]


    print

Comments

Thanks

¿What is the difference between cell.entities(0) and cell.entities(1)? ¿in general cell.entities(n)?

---

From the docstring:

"Return number of incident mesh entities of given topological dimension"

cell.entities(0) return incident vertices and cell.entities(1) return incident edges.

The above code could be written with entity iterators instead:

from dolfin import *
mesh = UnitCubeMesh(2, 2, 2)
coordinates = mesh.coordinates()

for cell in cells(mesh):
    print "Vertices of element", cell.index()
    for vertex in vertices(cell):
        print "\t", vertex.index(), coordinates[vertex.index()]
print

---

then cell.entities(2) return indicent faces? thanks

---

Yes, that is correct.