Mark the boundary faces with no-matching boundary mesh

1

Hi

I am trying to mark boundary faces. The boundary mesh is given in a file (say boundary.xml) other than the mesh file (say mesh.xml). The fact is the boundary mesh in boundary.xml will not exactly match what we have in the mesh.xml. Now I want to identify boundary faces of interest in mesh.xml and mark them. By “of interest”, I mean there exists faces in boundary.xml s.t. the face in mesh.xml is aligned with them (close enough and outward normals point to about the same direction).

How can I do this (if feasible)?

Thanks,
Victor

1 Like
2

Hi, if you define “closeness” of facets by vertices, isn’t it then redundant to include normals into the considerations? For the former, the following is a reasonable starting point

from dolfin import *
import numpy as np


def bdry_facets(mesh):
    '''Bdry facets given in therm of vertex indices of the mesh'''
    fdim = mesh.topology().dim() - 1 
    f = MeshFunction('size_t', mesh, fdim, 0)
    DomainBoundary().mark(f, 1)

    mesh.init(fdim, 0)
    f2v = mesh.topology()(fdim, 0)

    idx = [facet.index() for facet in SubsetIterator(f, 1)]

    return zip(map(f2v, idx), idx)


def first(iterable):
    '''First element'''
    return next(iter(iterable))


# --------------------------------------------------------------------


if __name__ == '__main__':
    from scipy.spatial.distance import cdist
    import itertools
    
    mesh1, mesh2 = UnitSquareMesh(5, 5), RectangleMesh(Point(1, 0), Point(2, 1), 2, 5)

    # Get the distance matrix between boundary points
    f1, f2 = bdry_facets(mesh1), bdry_facets(mesh2)
    v1, v2 = np.unique(np.hstack(map(first, f1))), np.unique(map(first, f2))

    distance_mat = cdist(mesh1.coordinates()[v1], mesh2.coordinates()[v2])

    close_enough = 1E-13
    # Points that are close enough
    close_pairs = [p for p, v in 
                   zip(itertools.product(v1, v2), distance_mat.flat)
                   if v < close_enough]
    close1, close2 = zip(*close_pairs)

    # Boundary facets of close points
    f1 = filter(lambda ((v0, v1), idx): v0 in close1 and v1 in close1, f1)
    f2 = filter(lambda ((v0, v1), idx): v0 in close2 and v1 in close2, f2)

    # Change lookup structure for facets of 2nd mesh
    keys, values = zip(*f2)
    f2 = dict(zip([tuple(sorted(k)) for k in keys], values))
    
    vertex_map = dict(zip(close1, close2))
    # Link the facets
    facet_map = {}
    while f1:
        (v0, v1), id1 = f1.pop()
        # Mapped facet
        mv0, mv1 = vertex_map[v0], vertex_map[v1]
        mf = (mv0, mv1) if mv0 < mv1 else (mv1, mv0)

        if mf in f2: facet_map[id1] = f2[mf]

    assert max(Facet(mesh1, f1).midpoint().distance(Facet(mesh2, f2).midpoint())
               for f1, f2 in facet_map.items()) < close_enough
3

Thanks MiroK for the detailed code. Back to the discussion, yes, for simple geometries, normal will be redundant. There are cases tho making this not hold. In 3D geometries, an edge is associated with two boundaries. Let’s say only one of them is your target boundary (and of course you will have the boundary mesh given). Now if you measure distance from triangles in this boundary mesh to the surface triangles from the tet mesh, it is possible that triangles on the other boundary associated to the edge being “close enough” to some triangles from the target boundary if you don’t have a smart way of defining the close_enough, then here is where the normal coming into play that if you find two triangles from different mesh are close, you check if their normal in similar direction to make sure you get the right triangles. This is not perfect, but safer.

Thanks anyways, will get back to you once I get a chance to read your code and see how I can incorporate the idea to my project.

Victor