How to obtain surfaces tag in Gmsh Python api?

i am trying to generate geometries and meshes with the Python api of Gmsh, planning to use it in FEniCS.

I started creating my geometry following the steps reported here: https://jsdokken.com/src/tutorial_gmsh.html

The author first create the volume and then retrieve the surfaces with the command:

surfaces = gmsh.model.occ.getEntities(dim=2)

Finally, he is able to relate the surface to the tag simply by finding the center of mass (com). He uses the command gmsh.model.occ.getCenterOfMass(dim,tag) and compares it with the know com position of his surfaces, like this:

inlet_marker, outlet_marker, wall_marker, obstacle_marker = 1, 3, 5, 7
walls = []
obstacles = []
for surface in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
    if np.allclose(com, [0, B/2, H/2]):
        gmsh.model.addPhysicalGroup(surface[0], [surface[1]], inlet_marker)
        inlet = surface[1]
        gmsh.model.setPhysicalName(surface[0], inlet_marker, "Fluid inlet")
    elif np.allclose(com, [L, B/2, H/2]):
        gmsh.model.addPhysicalGroup(surface[0], [surface[1]], outlet_marker)
        gmsh.model.setPhysicalName(surface[0], outlet_marker, "Fluid outlet")
    elif np.isclose(com[2], 0) or np.isclose(com[1], B) or np.isclose(com[2], H) or np.isclose(com[1],0):
        walls.append(surface[1])
    else:
        obstacles.append(surface[1])

Now, my problem is that this cannot work if two or more surfaces share the same com, such as two concentric cylinders.

How can i discriminate between them in such situation? For example in case of an hollow cylinder, i would like to have a tag for each surface in order to apply different boundary conditions in FEniCS.

Thanks in advance!

Solution

You can make use of gmsh.model.getAdjacencies(dim,tag) where dim and tag are the dimension and tag of your entity of interest. This functions returns two lists up, down. The first one gives you the tags of all entities adjacent (=neighbouring) to the entity of interest with dimension dim + 1. The second list gives you the tags of all entities adjacent to the entity of interest with dimension dim - 1.

In 3D (i.e. dim = 3) the up list will be empty because there are no 4D structures in gmsh. The down list will contain all surface tags the boundary of the volume is made of.

Below is an example code. Part 1 is straight forward and in Part 2 I added a functions that sorts the surface tags by their x-coordinate.

import gmsh

gmsh.initialize()

## PART 1:
tag_cylinder_1 = gmsh.model.occ.addCylinder(0, 0, 0, 1, 0, 0, 0.1)
tag_cylinder_2 = gmsh.model.occ.addCylinder(0, 0, 0, 1, 0, 0, 0.2)
gmsh.model.occ.synchronize()

up_cyl_1, down_cyl_1 = gmsh.model.getAdjacencies(3,tag_cylinder_1)
up_cyl_2, down_cyl_2 = gmsh.model.getAdjacencies(3,tag_cylinder_2)

com_1 = gmsh.model.occ.getCenterOfMass(2, down_cyl_1[0])
com_2 = gmsh.model.occ.getCenterOfMass(2, down_cyl_1[1])
com_3 = gmsh.model.occ.getCenterOfMass(2, down_cyl_1[2])


## PART 2:

def calcBoxVolume(box):
    dx = box[3] - box[0]    
    dy = box[4] - box[1] 
    dz = box[5] - box[2] 
    
    return dx*dy*dz

def getOrderedTags(cylTag):
    up, down = gmsh.model.getAdjacencies(3,cylTag)
    surf_COM = []
    for surface in down:
        com = [surface] + list(gmsh.model.occ.getCenterOfMass(2, surface))
        surf_COM.append(com)
        
    orderedSurfaces = sorted(surf_COM,key = lambda x: x[1])
    orderedSurfaceTags = [item[0] for item in orderedSurfaces]
    return orderedSurfaceTags

def setPhysicalTags(name,cylTag):
    orderedSurfaces = getOrderedTags(cylTag)
    
    gmsh.model.addPhysicalGroup(2, [orderedSurfaces[0]],name="inlet_"+name)
    gmsh.model.addPhysicalGroup(2, [orderedSurfaces[1]],name="tube_"+name)
    gmsh.model.addPhysicalGroup(2, [orderedSurfaces[3]],name="outlet_"+name)
    
def setPhysicalTagsCylDiff(name,cylTag):
    orderedSurfaces = getOrderedTags(cylTag)
    
    tag_A = orderedSurfaces[1]
    tag_B = orderedSurfaces[2]
    
    box_tube_A = gmsh.model.getBoundingBox(2,tag_A)
    box_tube_B = gmsh.model.getBoundingBox(2,tag_B)
    
    volBoxA = calcBoxVolume(box_tube_A)
    volBoxB = calcBoxVolume(box_tube_B)
    
    if volBoxA > volBoxB:
        innerTag = tag_B
        outerTag = tag_A
    else:
        innerTag = tag_A
        outerTag = tag_B
    
    gmsh.model.addPhysicalGroup(2, [orderedSurfaces[0]],name="inlet_"+name)
    gmsh.model.addPhysicalGroup(2, [innerTag],name="tube_inner_"+name)
    gmsh.model.addPhysicalGroup(2, [outerTag],name="tube_outer_"+name)
    gmsh.model.addPhysicalGroup(2, [orderedSurfaces[3]],name="outlet_"+name)
    
# setPhysicalTags("Cylinder_1",tag_cylinder_1)
# setPhysicalTags("Cylinder_2",tag_cylinder_2)

outDimTags, outDimTagsMap = gmsh.model.occ.cut([(3,tag_cylinder_2)],[(3,tag_cylinder_1)])
cylDiffTag = outDimTags[0][1]
gmsh.model.occ.synchronize()
setPhysicalTagsCylDiff("CylDiff",cylDiffTag)

gmsh.model.mesh.generate(2)
gmsh.fltk.run()
gmsh.finalize()