Getting coordinates of facets in dolfinx

Hi, i am trying to get the coordinates of facets, cells or surfaces in dolfinx. Some of the functions I use have a dependency on the coordinates of my system. As an example: I am trying to find the lowest point of this mesh (Taken from the dolfinx tutorial).

gmsh.initialize()
proc = MPI.COMM_WORLD.rank
top_marker = 2
bottom_marker = 1
left_marker = 1
if proc == 0:
    # We create one rectangle for each subdomain
    gmsh.model.occ.addRectangle(0, 0, 0, 1, 0.5, tag=1)
    gmsh.model.occ.addRectangle(0, 0.5, 0, 1, 0.5, tag=2)
    # We fuse the two rectangles and keep the interface between them
    gmsh.model.occ.fragment([(2, 1)], [(2, 2)])
    gmsh.model.occ.synchronize()

    # Mark the top (2) and bottom (1) rectangle
    top, bottom = None, None
    for surface in gmsh.model.getEntities(dim=2):
        com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
        if np.allclose(com, [0.5, 0.25, 0]):
            bottom = surface[1]
        else:
            top = surface[1]
    gmsh.model.addPhysicalGroup(2, [bottom], bottom_marker)
    gmsh.model.addPhysicalGroup(2, [top], top_marker)
    # Tag the left boundary
    left = []
    for line in gmsh.model.getEntities(dim=1):
        com = gmsh.model.occ.getCenterOfMass(line[0], line[1])
        if np.isclose(com[0], 0):
            left.append(line[1])
    gmsh.model.addPhysicalGroup(1, left, left_marker)
    gmsh.model.mesh.generate(2)
    gmsh.write("mesh.msh")
gmsh.finalize()

I can’t find an implementation similar to this post (Finding cell facet coordinates) in dolfinx, unfortunately. The method cells(), I can’t find in the python API. Did I miss something or is there a better implementation for this?
I import the data, using the approach of Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken
Many thanks in advance

Consider the following to get the coordinates for all entities of your mesh (can be restricted by changing np.arange(num_facets_owned_by_proc, dtype=np.int32) to the facets of interest)

import dolfinx

from mpi4py import MPI
import numpy as np
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 2, 2)
fdim = mesh.topology.dim - 1
mesh.topology.create_connectivity(fdim, 0)

num_facets_owned_by_proc = mesh.topology.index_map(fdim).size_local
geometry_entitites = dolfinx.cpp.mesh.entities_to_geometry(mesh, fdim, np.arange(num_facets_owned_by_proc, dtype=np.int32), False)
points = mesh.geometry.x
for e, entity in enumerate(geometry_entitites):
    print(e, points[entity])
3 Likes

Thank you very much :slight_smile:

1 Like