FiPy: How to find node (vertex) at interface between two mesh

I’ve defined two meshes in FiPy via Gmsh and would like to find the nodes at the interface between the two mesh. Is there a way to do this in FiPy?

siliconGeometry = '''
SetFactory("OpenCASCADE");
//set node spacing
ns = 1e-1;
ns2 = 1e-2;
x1 = 0;
y1 = 0;
x2 = 1;
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("Silicon") = {1};
'''

oxideGeometry = '''
SetFactory("OpenCASCADE");
//set node spacing
ns = 1e-1;
ns2 = 1e-2;
x1 = 0;
y1 = 0.5;
x2 = 1;
y2 = 1;
Point(5) = {x1, y1, 0, ns2};
Point(6) = {x2, y1, 0, ns2};
Point(7) = {x2, y2, 0, ns};
Point(8) = {x1, y2, 0, ns};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Physical Surface("Oxide") = {2};

m0 = Gmsh2D(siliconGeometry)
m1 = Gmsh2D(oxideGeometry)

I’d like to get all the nodes (or lines) at the interface between mesh m0 and m1.

Answer

FiPy has a function called nearest which gets the nearest values for two vectors. So, to get the overlapping IDs from both m0 and m1 then use

from fipy import Gmsh2D
from fipy.tools.numerix import nearest
import numpy as np

siliconGeometry = '''
...
'''

oxideGeometry = '''
...
'''
 
m0 = Gmsh2D(siliconGeometry)
m1 = Gmsh2D(oxideGeometry)

near_ids = nearest(m0.vertexCoords, m1.vertexCoords)
mask = np.all(np.isclose(m0.vertexCoords[:, near_ids], m1.vertexCoords), axis=0)
m1_close_ids = np.arange(len(m1.vertexCoords[0]))[mask]
m0_close_ids = near_ids[mask]

for i in range(len(m0_close_ids)):
    m0_id = m0_close_ids[i]
    m1_id = m1_close_ids[i]
    print()
    print(f'm0: {m0_id}, {m0.vertexCoords[:, m0_id]}')
    print(f'm1: {m1_id}, {m1.vertexCoords[:, m1_id]}')

We’ve used nearest to get the near_ids for m0 and then checked if the nodes are close using np.isclose. Using mask we can then construct the corresponding close IDs between the two meshes.

Note that it is probably much faster to use Scipy’s KDTree to find the nearest values. However, for this small problem it doesn’t matter, but for a very large mesh I imagine it would be important. In fact FiPy should really be using that internally (it isn’t currently).

Also, another tweak here is that you know that the upper and lower faces overlap so you can simply ask for the m0.facesTop and m1.facesBottom masks to reduce the problem to just the overlapping faces, but this relies on prior knowledge of the meshes.