mesh_from_Dict does not handle MultiLineStrings()

[duarte-jfs]

I have noticed an issue when using mesh_from_Dict(), stemming from the following lines:

       # Add overall bounding box
        total = unary_union(list(shapes_dict.values()))
        all_polygons = [total, *list(shapes_dict.values())]
        # Break up all shapes so that plane is tiled with non-overlapping layers, get the maximal number of fragments
        # Equivalent to BooleanFragments
        np.seterr(invalid="ignore")
        listpoly = [a.intersection(b) for a, b in combinations(all_polygons, 2)]
        np.seterr(**initial_settings)
        rings = [
            LineString(list(object.exterior.coords))
            for object in listpoly
            if not (object.geom_type in ["Point", "LineString"] or object.is_empty)
        ]

This is due to the fact that listpoly will have a MultiLineString when intersecting a LinearString with total. The change that could (probably) solve this is:

       # Add overall bounding box
        total = unary_union(list(shapes_dict.values()))
        all_polygons = [total, *list(shapes_dict.values())]
        # Break up all shapes so that plane is tiled with non-overlapping layers, get the maximal number of fragments
        # Equivalent to BooleanFragments
        np.seterr(invalid="ignore")
        listpoly = [a.intersection(b) for a, b in combinations(all_polygons, 2)]
        np.seterr(**initial_settings)
        rings = [
            LineString(list(object.exterior.coords))
            for object in listpoly
            if not (object.geom_type in ["Point", "LineString", "MultiLineString"] or object.is_empty)
        ]

Minimal working example

from collections import OrderedDict

import numpy as np

from shapely.geometry import box, LineString, MultiLineString
from skfem.io.meshio import from_meshio

from femwell.mesh import mesh_from_OrderedDict,mesh_from_Dict

box1 = box(-1,-1,1,1)
box2 = box(-1,5,1,6)

surface = box1.buffer(0.5).exterior
polygons = dict(box1 = box1,
               box2 = box2,
               surface = surface)

resolutions = dict(box1 = {'resolution': 0.2, 'distance':1},
                   box2 =  {'resolution': 0.2, 'distance':1},
                   surface =  {'resolution': 0.2, 'distance':1})

mesh = from_meshio(mesh_from_Dict(polygons, resolutions))

mesh.draw()

[HelgeGehring]

Great catch! Could you make a pull request?

We mostly use mesh_from_OrderedDict, as that way you can define which objects should be in the mesh in case things are overlapping. (I don´t remember what was the advantage of mesh_from_Dict, do you @simbilod ?) Could you make sure that it also works with that one?

[duarte-jfs]

Could you make a pull request?

Once I figure out how, I will XD

We mostly use mesh_from_OrderedDict, as that way you can define which objects should be in the mesh in case things are overlapping. (I don´t remember what was the advantage of mesh_from_Dict, do you @simbilod ?) Could you make sure that it also works with that one?

The same test with OrderedDict worked fine, and I assume is because it filters if it is a Polygon or not, rather than if it is a point or line.

[HelgeGehring]

Yay, sounds great! 🥳

[simbilod]

I think mesh_from_Dict was meant to add the ability to have different polygons be tagged the same way

I should really finish the conversion to meshwell at some point so we have one place to discuss the meshing issues :)