Find coordinates of cells in a subdomain

General
1

I need to rotate the mesh in a set of subdomains by transforming relevant points in mesh.geometry.x. As a requirement, I need to do it through domain tags defined in gmsh and not geometrical markers. How can we know which points in mesh.geometry.x belong to cells in a specific subdomain?

Below is a simple example of the geometry. I will appreciate if someone could complete it by adding how to determine the indices of mesh.geometry.x given the domain tag.

import gmsh, mpi4py
import dolfinx
from dolfinx.io import gmshio

# disp = True
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)  # disable output messages
gmsh.clear()
gdim = 2
model_rank = 0

occ = gmsh.model.occ

# first we create the outer domain
outer_box = occ.add_rectangle(-1, -1, 0, 2, 2)
occ.synchronize()
air_shell = occ.addDisk(0, 0, 0, 0.5, 0.5)
dom1 = occ.cut([(gdim, outer_box)], [(gdim, air_shell)])  # one domain

# now the inner domain which is the wheel
dom2 = occ.addDisk(0, 0, 0, 0.5, 0.5)
occ.synchronize() # now we have two non-conformal domains

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.refine()
gmsh.model.mesh.refine()

model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(
    gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)

gmsh.finalize()

# TO ADD: how to retrieve indices of points in mesh.geometry.x belonging to subdomain 2?
2

Something like

points = mesh.geometry.x
cells_map = mesh.topology.index_map(mesh.topology.dim)
num_cells = cells_map.size_local + cells_map.num_ghosts
mesh.topology.create_connectivity(mesh.topology.dim, 0)
connectivity_cells_to_vertices = mesh.topology.connectivity(mesh.topology.dim, 0)
vertex_map = {
    topology_index: geometry_index for c in range(num_cells) for (topology_index, geometry_index) in zip(
        connectivity_cells_to_vertices.links(c), mesh.geometry.dofmap[c])
}

subdomain_marker = 1
subdomain_indices = ct.indices[ct.values == subdomain_marker]
for c in subdomain_indices:
    vertices = [vertex_map[v] for v in connectivity_cells_to_vertices.links(c)]
    print(c, points[vertices])

should work.

I’ve run it successfuly, but given who complex your domain is I can’t really verify if the code is correct, because I have no idea what the subdomain with ID 1 is. It would have been much simpler if you had provided a simple mock domain (something like a square split into two parts), so I’ll leave it to you to check my code for correctness.

1 Like
3

entities_to_geometry can give you the geometry i ndices of a cell.
Get the cells with
cells = ct.find(marker) and pass it through
to the aforementioned function. See for instance Dolfinx.cpp.mesh.entities_to_geometry problem - #2 by dokken

1 Like
4

Thanks a lot @francesco-ballarin and @dokken for your quick feedback!