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

duarte-jfs · 2024-07-30T15:37:12Z

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

The text was updated successfully, but these errors were encountered:

kinnala · 2024-07-30T16:58:27Z

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

duarte-jfs · 2024-07-30T17:05:51Z

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 · 2024-07-30T17:11:15Z

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 · 2024-07-30T17:25:23Z

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 · 2024-07-30T20:16:59Z

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 · 2024-07-31T07:30:49Z

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 · 2024-07-31T09:23:32Z

I believe this is fixed by the PR #1152

kinnala · 2024-08-06T09:17:26Z

I have released this in 10.0.1.

kinnala mentioned this issue Jul 31, 2024

preserve order of Mesh.boundaries during loading #1152

Merged

duarte-jfs mentioned this issue Aug 5, 2024

Improvements on RF waveguide design tutorial HelgeGehring/femwell#164

Open

kinnala closed this as completed in #1152 Aug 6, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

duarte-jfs commented Jul 30, 2024 •

edited

Loading

kinnala commented Jul 30, 2024

duarte-jfs commented Jul 30, 2024 •

edited

Loading

duarte-jfs commented Jul 30, 2024

kinnala commented Jul 30, 2024

duarte-jfs commented Jul 30, 2024

kinnala commented Jul 31, 2024

kinnala commented Jul 31, 2024

kinnala commented Aug 6, 2024

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

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

Comments

duarte-jfs commented Jul 30, 2024 • edited Loading

kinnala commented Jul 30, 2024

duarte-jfs commented Jul 30, 2024 • edited Loading

duarte-jfs commented Jul 30, 2024

kinnala commented Jul 30, 2024

duarte-jfs commented Jul 30, 2024

kinnala commented Jul 31, 2024

kinnala commented Jul 31, 2024

kinnala commented Aug 6, 2024

duarte-jfs commented Jul 30, 2024 •

edited

Loading

duarte-jfs commented Jul 30, 2024 •

edited

Loading