OrientedBoundary are not correct with from_meshio() when lines are given as shapely.Polygon.exterior

duarte-jfs commented 3 months ago

I have the code below to illustrate the issue. One box I create a OrientedBoundary with mesh.facets_around_elements() and the other is generated by giving the exterior of a shapely.Polygon.

from skfem.io.meshio import from_meshio
from femwell.mesh import mesh_from_OrderedDict
import shapely

from skfem import ElementTriP1, Basis

import matplotlib.pyplot as plt
import numpy as np

%matplotlib widget

polygons = {'poly_box': shapely.box(0.2,0.2,0.4,0.4),
            'poly_exterior': shapely.box(0.6, 0.6, 0.8, 0.8).exterior,
            'background': shapely.box(0, 0, 1, 1)}

mesh = from_meshio(
    mesh_from_OrderedDict(polygons, {}, default_resolution_max=0.1, filename="mesh.msh")
)
basis = Basis(mesh, ElementTriP1())

boundary1 = mesh.facets_around(mesh.subdomains['poly_box'])
boundary2 = mesh.boundaries['poly_exterior']

facets1 = mesh.facets[:, boundary1]
facets2 = mesh.facets[:, boundary2]

facet_basis1 = basis.boundary(boundary1)
facet_basis2 = basis.boundary(boundary2)

plt.figure()
plt.scatter(*mesh.p, marker = 'o', facecolor = 'None', edgecolor = 'black', s = 15)

for facets, facet_basis in zip([facets1, facets2],[facet_basis1, facet_basis2]):
    colors = plt.cm.jet(np.linspace(0,1,facets.shape[1]))
    for i in range(facets.shape[1]):
        x = [mesh.p[0,facets[0,i]], mesh.p[0,facets[1,i]]]
        y = [mesh.p[1,facets[0,i]], mesh.p[1,facets[1,i]]]

        norm = np.asarray([facet_basis.normals[0,i,0], 
                           facet_basis.normals[1,i,0]])

        norm = norm/np.sqrt(np.sum(np.abs(norm)**2)) * 0.1

        plt.arrow(np.mean(x), np.mean(y), dx = norm[0], dy = norm[1], color = colors[i], 
                  head_width = 0.03, 
                  head_length = 0.01,
                  length_includes_head = True)

        plt.plot(x, y, color = colors[i])
        plt.text(np.mean(x), np.mean(y), f"{i}")

I'm unsure how to solve this, but my guess is that it originates on the from_meshio(). This is a very important feature when considering line integrals

@HelgeGehring

kinnala commented 3 months ago

What is mesh_from_OrderedDict and how is the orientation information represented in its output?

duarte-jfs commented 3 months ago

mesh_from_OrderedDict comes from femwell (https://github.com/HelgeGehring/femwell/blob/main/femwell/mesh/mesh.py).

how is the orientation information represented in its output

That I don't quite follow. What should I be looking for? I believe the orientation is still correctly passed from mesh_from_OrderedDict to from_meshio() because when following the tangents of the boundary from:

tangents = (mtmp.p[:, facets[1]] - mtmp.p[:, facets[0]])

They follow the rectangle in an anticlockwise manner. The problem was on the t1, _ = mtmp.f2t[:, boundaries[k]]. The t1 elements returned for the poly_exterior were outside the rectangle boundary whereas the elements represented by t1 for the boundary surrounding the poly_box (which the mesh returns as a boundary too) returned the elements inside the polygon boundary. But I must admit, I debugging beyond that proved quite difficult to me

The following code allows us to see that:

def from_meshio(m,
                out=None,
                int_data_to_sets=False,
                force_meshio_type=None,
                ignore_orientation=False,
                ignore_interior_facets=False):

    if ignore_interior_facets:
        logger.warning("kwarg ignore_interior_facets is unused.")

    cells = m.cells_dict
    meshio_type = None

    if force_meshio_type is None:
        # detect 3D
        for k in cells:
            if k in {'tetra',
                     'hexahedron',
                     'tetra10',
                     'hexahedron27',
                     'wedge'}:
                meshio_type = k
                break

        if meshio_type is None:
            # detect 2D
            for k in cells:
                if k in {'triangle',
                         'quad',
                         'triangle6',
                         'quad9'}:
                    meshio_type = k
                    break

        if meshio_type is None:
            # detect 1D
            for k in cells:
                if k == 'line':
                    meshio_type = k
                    break
    else:
        meshio_type = force_meshio_type

    if meshio_type is None:
        raise NotImplementedError("Mesh type(s) not supported "
                                  "in import: {}.".format(cells.keys()))

    mesh_type = MESH_TYPE_MAPPING[meshio_type]

    # create p and t
    p = np.ascontiguousarray(mesh_type.strip_extra_coordinates(m.points).T)
    t = np.ascontiguousarray(cells[meshio_type].T)

    # reorder t if needed
    if meshio_type == 'hexahedron':
        t = t[INV_HEX_MAPPING[:8]]
    elif meshio_type == 'hexahedron27':
        t = t[INV_HEX_MAPPING]

    if int_data_to_sets:
        m.int_data_to_sets()

    subdomains = {}
    boundaries = {}

    # parse any subdomains from cell_sets
    if m.cell_sets:
        subdomains = {k: v[meshio_type].astype(np.int64)
                      for k, v in m.cell_sets_dict.items()
                      if meshio_type in v}

    # create temporary mesh for matching boundary elements
    mtmp = mesh_type(p, t, validate=False)
    bnd_type = BOUNDARY_TYPE_MAPPING[meshio_type]
    #############################################
    mesh = mtmp

    plt.figure()
    plt.scatter(*mesh.p, marker = 'o', facecolor = 'None', edgecolor = 'black', s = 15)

     #################################################   
    # parse boundaries from cell_sets
    # print('inside meshio:', bnd_type)
    if m.cell_sets and bnd_type in m.cells_dict:
        p2f = mtmp.p2f
        for k, v in m.cell_sets_dict.items():
            if bnd_type in v and k.split(":")[0] != "gmsh":
                facets = m.cells_dict[bnd_type][v[bnd_type]].T
                sorted_facets = np.sort(facets, axis=0)
                ind = p2f[:, sorted_facets[0]]
                for itr in range(sorted_facets.shape[0] - 1):
                    ind = ind.multiply(p2f[:, sorted_facets[itr + 1]])
                boundaries[k] = np.nonzero(ind)[0]

                if not ignore_orientation:
                    try:
                        print(boundaries[k])
                        ori = np.zeros_like(boundaries[k], dtype=np.float64)
                        t1, _ = mtmp.f2t[:, boundaries[k]]
                        if facets.shape[0] == 2:
                            tangents = (mtmp.p[:, facets[1]]
                                        - mtmp.p[:, facets[0]])
                            normals = np.array([-tangents[1], tangents[0]])

                            tangents
                        elif facets.shape[0] == 3:
                            tangents1 = (mtmp.p[:, facets[1]]
                                         - mtmp.p[:, facets[0]])
                            tangents2 = (mtmp.p[:, facets[2]]
                                         - mtmp.p[:, facets[0]])
                            normals = -np.cross(tangents1.T, tangents2.T).T
                        else:
                            raise NotImplementedError

                        for r in range(len(t1)):
                            plt.plot(*mtmp.p[:, mtmp.t[(0,1), t1[r]]], color = 'blue')
                            plt.plot(*mtmp.p[:, mtmp.t[(1,2), t1[r]]], color = 'blue')
                            plt.plot(*mtmp.p[:, mtmp.t[(0,2), t1[r]]], color = 'blue')

                        for itr in range(mtmp.t.shape[0]): #iterate over the points of the triangles that contain the boundary

                            plt.scatter(*mtmp.p[:, facets[0]], edgecolor = 'red', facecolor = 'none', s = 20)
                            plt.scatter(*mtmp.p[:, mtmp.t[itr, t1]], edgecolor = 'blue', facecolor = 'none', s = 25)
                            for g in range(len(facets[0])):
                                mesh = mtmp
                                x = [mesh.p[0,facets[0,g]], mesh.p[0,facets[1,g]]]
                                y = [mesh.p[1,facets[0,g]], mesh.p[1,facets[1,g]]]
                                plt.text(np.mean(x), np.mean(y), f"{g}")

                            # for l in range(len(facets[0])):
                            #     plt.arrow(*mtmp.p[:, facets[0][l]],*(mtmp.p[:, facets[0][l]]- 
                            #                                          mtmp.p[:, mtmp.t[itr, t1][l]]))

                            # print('##########################################')
                            ori += np.sum(normals
                                          * (mtmp.p[:, facets[0]] #This is every point of the line
                                             - mtmp.p[:, mtmp.t[itr, t1]]),
                                          axis=0)
                        ori = 1 * (ori > 0)
                        # print(ori)
                        boundaries[k] = OrientedBoundary(boundaries[k],
                                                         ori)
                    except Exception as e:
                        print(e)
                        logger.warning("Failure to orient a boundary.")
    # print(boundaries)
    # MSH 2.2 tag parsing
    if len(boundaries) == 0 and m.cell_data and m.field_data:
        try:
            elements_tag = m.cell_data_dict['gmsh:physical'][meshio_type]
            subdomains = {}
            tags = np.unique(elements_tag)

            def find_tagname(tag):
                for key in m.field_data:
                    if m.field_data[key][0] == tag:
                        return key
                return None

            for tag in tags:
                t_set = np.nonzero(tag == elements_tag)[0]
                subdomains[find_tagname(tag)] = t_set

            # find tagged boundaries
            if bnd_type in m.cell_data_dict['gmsh:physical']:
                facets = m.cells_dict[bnd_type]
                facets_tag = m.cell_data_dict['gmsh:physical'][bnd_type]

            # put meshio facets to dict
            dic = {tuple(np.sort(facets[i])): facets_tag[i]
                   for i in range(facets.shape[0])}

            # get index of corresponding Mesh.facets for each meshio
            # facet found in the dict
            index = np.array([[dic[tuple(np.sort(mtmp.facets[:, i]))], i]
                              for i in mtmp.boundary_facets()
                              if tuple(np.sort(mtmp.facets[:, i])) in dic])

            # read meshio tag numbers and names
            tags = index[:, 0]
            boundaries = {}
            for tag in np.unique(tags):
                tagindex = np.nonzero(tags == tag)[0]
                boundaries[find_tagname(tag)] = index[tagindex, 1]

        except Exception:
            logger.warning("Failure to parse tags from meshio.")

    # attempt parsing skfem tags
    if m.cell_data:
        _boundaries, _subdomains = mtmp._decode_cell_data(m.cell_data)
        boundaries.update(_boundaries)
        subdomains.update(_subdomains)

    # export mesh data
    if out is not None and isinstance(out, list):
        for i, field in enumerate(out):
            out[i] = getattr(m, field)

    return mesh_type(
        p,
        t,
        None if len(boundaries) == 0 else boundaries,
        None if len(subdomains) == 0 else subdomains,
    )

duarte-jfs commented 3 months ago

One possible fix that can work generally for any 2D mesh is by making use of the correct calculation of the tangents and normals, and extend it to the third dimension, then performing a cross product:

def from_meshio(m,
                out=None,
                int_data_to_sets=False,
                force_meshio_type=None,
                ignore_orientation=False,
                ignore_interior_facets=False):

    if ignore_interior_facets:
        logger.warning("kwarg ignore_interior_facets is unused.")

    cells = m.cells_dict
    meshio_type = None

    if force_meshio_type is None:
        # detect 3D
        for k in cells:
            if k in {'tetra',
                     'hexahedron',
                     'tetra10',
                     'hexahedron27',
                     'wedge'}:
                meshio_type = k
                break

        if meshio_type is None:
            # detect 2D
            for k in cells:
                if k in {'triangle',
                         'quad',
                         'triangle6',
                         'quad9'}:
                    meshio_type = k
                    break

        if meshio_type is None:
            # detect 1D
            for k in cells:
                if k == 'line':
                    meshio_type = k
                    break
    else:
        meshio_type = force_meshio_type

    if meshio_type is None:
        raise NotImplementedError("Mesh type(s) not supported "
                                  "in import: {}.".format(cells.keys()))

    mesh_type = MESH_TYPE_MAPPING[meshio_type]

    # create p and t
    p = np.ascontiguousarray(mesh_type.strip_extra_coordinates(m.points).T)
    t = np.ascontiguousarray(cells[meshio_type].T)

    # reorder t if needed
    if meshio_type == 'hexahedron':
        t = t[INV_HEX_MAPPING[:8]]
    elif meshio_type == 'hexahedron27':
        t = t[INV_HEX_MAPPING]

    if int_data_to_sets:
        m.int_data_to_sets()

    subdomains = {}
    boundaries = {}

    # parse any subdomains from cell_sets
    if m.cell_sets:
        subdomains = {k: v[meshio_type].astype(np.int64)
                      for k, v in m.cell_sets_dict.items()
                      if meshio_type in v}

    # create temporary mesh for matching boundary elements
    mtmp = mesh_type(p, t, validate=False)
    bnd_type = BOUNDARY_TYPE_MAPPING[meshio_type]

    # parse boundaries from cell_sets

    if m.cell_sets and bnd_type in m.cells_dict:
        p2f = mtmp.p2f
        for k, v in m.cell_sets_dict.items():
            if bnd_type in v and k.split(":")[0] != "gmsh":
                facets = m.cells_dict[bnd_type][v[bnd_type]].T
                sorted_facets = np.sort(facets, axis=0)
                ind = p2f[:, sorted_facets[0]]
                for itr in range(sorted_facets.shape[0] - 1):
                    ind = ind.multiply(p2f[:, sorted_facets[itr + 1]])
                boundaries[k] = np.nonzero(ind)[0]

                if not ignore_orientation:
                    try:
                        ori = np.zeros_like(boundaries[k], dtype=np.float64)
                        t1, _ = mtmp.f2t[:, boundaries[k]]
                        if facets.shape[0] == 2:
                            tangents = (mtmp.p[:, facets[1]]
                                        - mtmp.p[:, facets[0]])
                            normals = np.array([-tangents[1], tangents[0]])
                            print(tangents)
                            tangents = np.vstack((tangents, np.zeros(facets.shape[1])))
                            normals = np.vstack((normals, np.zeros(facets.shape[1])))

                            ori = (np.sign(-np.cross(tangents.T, normals.T).T[-1])-1)/2
                        elif facets.shape[0] == 3:
                            tangents1 = (mtmp.p[:, facets[1]]
                                         - mtmp.p[:, facets[0]])
                            tangents2 = (mtmp.p[:, facets[2]]
                                         - mtmp.p[:, facets[0]])
                            normals = -np.cross(tangents1.T, tangents2.T).T

                            for itr in range(mtmp.t.shape[0]): 
                                ori += np.sum(normals
                                              * (mtmp.p[:, facets[0]] #This is every point of the line
                                                 - mtmp.p[:, mtmp.t[itr, t1]]),
                                              axis=0)
                            ori = 1 * (ori > 0)
                        else:
                            raise NotImplementedError

                        boundaries[k] = OrientedBoundary(boundaries[k],
                                                         ori)
                    except Exception as e:
                        logger.warning("Failure to orient a boundary.")
    # print(boundaries)
    # MSH 2.2 tag parsing
    if len(boundaries) == 0 and m.cell_data and m.field_data:
        try:
            elements_tag = m.cell_data_dict['gmsh:physical'][meshio_type]
            subdomains = {}
            tags = np.unique(elements_tag)

            def find_tagname(tag):
                for key in m.field_data:
                    if m.field_data[key][0] == tag:
                        return key
                return None

            for tag in tags:
                t_set = np.nonzero(tag == elements_tag)[0]
                subdomains[find_tagname(tag)] = t_set

            # find tagged boundaries
            if bnd_type in m.cell_data_dict['gmsh:physical']:
                facets = m.cells_dict[bnd_type]
                facets_tag = m.cell_data_dict['gmsh:physical'][bnd_type]

            # put meshio facets to dict
            dic = {tuple(np.sort(facets[i])): facets_tag[i]
                   for i in range(facets.shape[0])}

            # get index of corresponding Mesh.facets for each meshio
            # facet found in the dict
            index = np.array([[dic[tuple(np.sort(mtmp.facets[:, i]))], i]
                              for i in mtmp.boundary_facets()
                              if tuple(np.sort(mtmp.facets[:, i])) in dic])

            # read meshio tag numbers and names
            tags = index[:, 0]
            boundaries = {}
            for tag in np.unique(tags):
                tagindex = np.nonzero(tags == tag)[0]
                boundaries[find_tagname(tag)] = index[tagindex, 1]

        except Exception:
            logger.warning("Failure to parse tags from meshio.")

    # attempt parsing skfem tags
    if m.cell_data:
        _boundaries, _subdomains = mtmp._decode_cell_data(m.cell_data)
        boundaries.update(_boundaries)
        subdomains.update(_subdomains)

    # export mesh data
    if out is not None and isinstance(out, list):
        for i, field in enumerate(out):
            out[i] = getattr(m, field)

    return mesh_type(
        p,
        t,
        None if len(boundaries) == 0 else boundaries,
        None if len(subdomains) == 0 else subdomains,
    )

I don't know if the issue will remain for 3D, or even if this is a good permanent fix, but this is as far as I got. Hope it helps

kinnala commented 3 months ago

I believe the orientation part of from_meshio was created while testing against Gmsh generated meshes, loaded from file via meshio, and if I recall correctly, in Gmsh meshes the boundary/interface facets were listed in a specific (probably anticlockwise) order. In scikit-fem this information was made more explicit by introducing the concept of OrientedBoundary, so that a consistent normal vector could be used.

However, I’ll need to study where these meshes are originating from and how do they differ with the Gmsh ones, or if there is another bug I am not aware of.

Does the problem persist if you do a save-load cycle with the meshio.Mesh object? Can you save one of those meshio objects to file and post it here so I can test it more easily?

duarte-jfs commented 3 months ago

Does the problem persist if you do a save-load cycle with the meshio.Mesh object?

I have done the following:


mesh = from_meshio(
    mesh_from_OrderedDict(polygons, {}, default_resolution_max=0.1, filename="mesh.msh")
)

mesh = from_file('mesh.msh', out = None)

and the problem persists. Is this what you meant?

Can you save one of those meshio objects to file and post it here so I can test it more easily?

Do you mean the .msh file? mesh.zip

kinnala commented 3 months ago

Thank you. I believe this issue originates from the optimizations done in #1117.

In fact, it seems that

                        for itr in range(mtmp.t.shape[0]):
                            ori += np.sum(normals
                                          * (mtmp.p[:, facets[0]]
                                             - mtmp.p[:, mtmp.t[itr, t1]]),
                                          axis=0)

which is supposed to check the dot product between the proposed normal and the edges of the elements is completely incorrect.

More precisely, it seems that the order of facets in facets and the order of the corresponding triangles in t1 is not equal. I will continue looking into it and see if there is a simple fix for all dimensions.

kinnala commented 3 months ago

I believe this is fixed by the PR #1152

kinnala commented 3 months ago

I have released this in 10.0.1.

kinnala / scikit-fem

OrientedBoundary are not correct with from_meshio() when lines are given as shapely.Polygon.exterior #1150