Can you create separate meshes in FiPy by reading "physical names" from Gmsh?

Here is a refinement on the brute-force solution I initially offered:

import fipy as fp

mesh = fp.Gmsh2D('''
SetFactory("OpenCASCADE");
//
ns = 1e-1;
ns2 = 1e-2;
//
x1 = 0;
x2 = 1;
y1 = 0;
y2 = 0.5;
Point(1) = {x1, y1, 0, ns};
Point(2) = {x2, y1, 0, ns};
Point(3) = {x2, y2, 0, ns2};
Point(4) = {x1, y2, 0, ns2};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Physical Surface("Oxide") = {1};
//
y1 = 0.5;
y2 = 1.0;
shiftX = 0.4;
Point(51) = {x1, y1, 0, ns2};
Point(61) = {x1+shiftX, y1, 0, ns2};
Point(71) = {x1+shiftX, y2, 0, ns2};
Point(81) = {x1, y2, 0, ns2};
Line(51) = {51, 61};
Line(61) = {61, 71};
Line(71) = {71, 81};
Line(81) = {81, 51};
Curve Loop(21) = {51, 61, 71, 81};
Plane Surface(21) = {21};
Physical Surface("Silicon") = {21};
//
Point(52) = {x2-shiftX, y1, 0, ns2};
Point(62) = {x2, y1, 0, ns2};
Point(72) = {x2, y2, 0, ns2};
Point(82) = {x2-shiftX, y2, 0, ns2};
Line(52) = {52, 62};
Line(62) = {62, 72};
Line(72) = {72, 82};
Line(82) = {82, 52};
Curve Loop(22) = {52, 62, 72, 82};
Plane Surface(22) = {22};
Physical Surface("Silicon") += {22};
//
y1 = 1.0;
y2 = 1.5;
shiftX = 0.4;
shiftY = 0.5;
Point(9) = {x1, y1, 0, ns2};
Point(91) = {x1+shiftX, y1, 0, ns2};
Point(92) = {x1+shiftX, y1-shiftY, 0, ns2};
Point(101) = {x2-shiftX, y1-shiftY, 0, ns2};
Point(102) = {x2-shiftX, y1, 0, ns2};
Point(10) = {x2, y1, 0, ns2};
Point(11) = {x2, y2, 0, ns};
Point(12) = {x1, y2, 0, ns};
Line(9) = {9, 91};
Line(91) = {91, 92};
Line(92) = {92, 101};
Line(101) = {101, 102};
Line(102) = {102, 10};
Line(10) = {10, 11};
Line(11) = {11, 12};
Line(12) = {12, 9};
Curve Loop(3) = {9, 91, 92, 101, 102, 10, 11, 12};
Plane Surface(3) = {3};
Physical Surface("Gas") = {3};
//''')

from fipy.meshes.mesh2D import Mesh2D

def extract_mesh(mesh, mask):
    cellFaceIDs = mesh.cellFaceIDs[..., mask]
    faceIDs = numerix.unique(cellFaceIDs.flatten())
    facemap = numerix.zeros(mesh.faceVertexIDs.shape[1], dtype=int)
    facemap[faceIDs] = faceIDs.argsort()
    
    faceVertexIDs = mesh.faceVertexIDs[..., faceIDs]
    vertIDs = numerix.unique(faceVertexIDs.flatten())
    vertmap = numerix.zeros(mesh.vertexCoords.shape[1], dtype=int)
    vertmap[vertIDs] = vertIDs.argsort()

    return Mesh2D(mesh.vertexCoords[..., vertIDs],
                  vertmap[faceVertexIDs],
                  facemap[cellFaceIDs])
    
mesh_gas = extract_mesh(mesh=mesh, mask=mesh.physicalCells["Gas"])
mesh_silicon = extract_mesh(mesh=mesh, mask=mesh.physicalCells["Silicon"])
mesh_oxide = extract_mesh(mesh=mesh, mask=mesh.physicalCells["Oxide"])

The derived meshes have compact vertex and face lists and reasonable numbers of exterior faces.

nearest still won't do what you want, because the vast majority of the faces (exterior or otherwise) of the gas mesh are not anywhere near the oxide mesh. You only want the faces that are coincident. To see how complicated that is, you need about 98% of this code.

I think a much better answer is to define the boundaries you are interested in using the Gmsh GEO script, with Physical Line definitions. These can then be extracted as mesh.physicalFaces (some further work would need to be done to extract_mesh(), but I think it's tractable).

Thanks jeguyer. I was able to create separate mesh with the physicalCell identifier. In regards to finding the vertices at gas/oxide interface, I was able to trace the interface contour using surfCoords from the following code:

gasFaces = gasMesh.exteriorFaces
gasVertices = numerix.unique(gasMesh.faceVertexIDs[..., gasFaces].flatten())
gasVertexCoords = gasMesh.vertexCoords[..., gasVertices]

oxideFaces = oxideMesh.exteriorFaces
oxideVertices = numerix.unique(oxideMesh.faceVertexIDs[..., oxideFaces].flatten())
oxideVertexCoords = oxideMesh.vertexCoords[..., oxideVertices]

surfVertexID = numerix.nearest(gasVertexCoords, oxideVertexCoords)
surfCoords = gasVertexCoords[..., numerix.unique(surfVertexID)]