Jb merge endorse

setup.py

@@ -36,7 +36,7 @@ def read(*names, **kwargs):
     license='GPL 3.0',
     description='Multilevel Monte Carlo method.',
     long_description=long_description,
-    author='Jan Brezina, Martin Spetlik, Klara Steklova',
+    author='Jan Brezina, Martin Spetlik',
     author_email='jan.brezina@tul.cz',
     url='https://github.com/GeoMop/MLMC',
     classifiers=[
@@ -65,9 +65,9 @@ def read(*names, **kwargs):
     py_modules=[splitext(basename(path))[0] for path in glob.glob('src/*.py')],
     package_data={
         # If any package contains *.txt or *.rst files, include them:
-        '': ['*.txt', '*.rst'],
+        #'': ['*.txt', '*.rst'],
         # And include any *.msg files found in the 'hello' package, too:
-        'hello': ['*.msg'],
+        #'hello': ['*.msg'],
     },
 
     # include automatically all files in the template MANIFEST.in
@@ -77,9 +77,8 @@ def read(*names, **kwargs):
         # eg: 'aspectlib==1.1.1', 'six>=1.7',
     ],
     extras_require={
-        # eg:
-        #   'rst': ['docutils>=0.11'],
-        #   ':python_version=="2.6"': ['argparse'],
+        # requirements for optional features
+        'gmsh': ['gmsh-tools'],
     }
     # entry_points={
     #     'console_scripts': [

src/frac_geom.py

@@ -5,7 +5,7 @@
 import geomop.format_last as lg
 import geomop.layers_io
 import geomop.geometry
-#from geomop.plot_polygons import plot_polygon_decomposition
+from geomop.plot_polygons import plot_polygon_decomposition
 
 
 
@@ -15,85 +15,86 @@
 
 
 
-def make_frac_mesh(box, mesh_step, fractures, frac_step):
+def make_frac_mesh(root_polygon, mesh_step:float, fractures, frac_step:float, mesh_base="fractured_2d"):
     """
+    :param root_polygon: List[Point2d]
+    :param fractures: List[(Point2d, Point2d)]
+
     Make geometry and mesh for given 2d box and set of fractures.
     :param box: [min_point, max_point]; points are np.arrays
     :param fractures: Array Nx2x2, one row for every fracture given by endpoints: [p0, p1]
     :return: GmshIO object with physical groups:
-        box: 1,
-        fractures: 1000 + i, i = 0, ... , N-1
+        "bulk": 1
+        "side_<n>", n + 1
+        "frac_<n>", 1000 + n
     """
-    regions = make_regions(mesh_step, fractures, frac_step)
-    decomp, reg_map = make_decomposition(box, fractures, regions)
-    geom = fill_lg(decomp, reg_map, regions)
+    regions = []
+    add_reg(regions, "NONE", -1, not_used=True)
+    i_r_bulk = add_reg(regions, "bulk", 2, mesh_step)
+    i_r_side = [
+        add_reg(regions, "side_{}".format(s_id), 1, step=frac_step, bc=True)
+        for s_id in range(len(root_polygon))
+    ]
+    i_r_frac = [
+        add_reg(regions, "frac_{}".format(f_id), 1, step=frac_step)
+        for f_id in range(len(fractures))
+    ]
+    decomp = make_decomposition(root_polygon, fractures, regions, i_r_bulk, i_r_side, i_r_frac, frac_step)
+
+    geom = fill_lg(decomp, regions, mesh_base=mesh_base)
     return make_mesh(geom)
 
 
 def add_reg(regions, name, dim, step=0.0, bc=False, not_used =False):
     reg = lg.Region(dict(name=name, dim=dim, mesh_step=step, boundary=bc, not_used=not_used))
     reg._id = len(regions)
     regions.append(reg)
+    return reg._id
 
-def make_regions(mesh_step, fractures, frac_step):
-    regions = []
-    add_reg(regions, "NONE", -1, not_used=True)
-    add_reg(regions, "bulk_0", 2, mesh_step)
-    add_reg(regions, ".bc_inflow", 1, bc=True)
-    add_reg(regions, ".bc_outflow", 1, bc=True)
-    for f_id in range(len(fractures)):
-        add_reg(regions, "frac_{}".format(f_id), 1, frac_step)
-    return regions
-
-
-def make_decomposition(box, fractures, regions):
-    box_pd = poly.PolygonDecomposition()
-    p00, p11 = box
-    p01 = np.array([p00[0], p11[1]])
-    p10 = np.array([p11[0], p00[1]])
-    box_pd.add_line(p00, p01)
-    seg_outflow, = box_pd.add_line(p01, p11)
-    box_pd.add_line(p11, p10)
-    seg_inflow, = box_pd.add_line(p10, p00)
-
-    decompositions = [box_pd]
-    for p0, p1 in fractures:
-        pd = poly.PolygonDecomposition()
-        pd.add_line(p0, p1)
-        decompositions.append(pd)
-
-    common_decomp, maps = merge.intersect_decompositions(decompositions)
-    #plot_polygon_decomposition(common_decomp)
-    #print(maps)
-
-    # Map common_decomp objects to regions.
-    none_region_id = 0
-    box_reg_id = 1
-    bc_inflow_id = 2
-    bc_outflow_id = 3
-    frac_id_shift = 4
-    decomp_shapes = [common_decomp.points, common_decomp.segments, common_decomp.polygons]
-    reg_map = [{key: regions[none_region_id] for key in decomp_shapes[d].keys()} for d in range(3)]
-    for i_frac, f_map in enumerate(maps[1:]):
-        for id, orig_seg_id in f_map[1].items():
-            reg_map[1][id] = regions[frac_id_shift + i_frac]
-
-    for id, orig_poly_id in maps[0][2].items():
-        if orig_poly_id == 0:
-            continue
-        reg_map[2][id] = regions[box_reg_id]
-
-    for id, orig_seg_id in maps[0][1].items():
-        if orig_seg_id == seg_inflow.id:
-            reg_map[1][id] = regions[bc_inflow_id]
-        if orig_seg_id == seg_outflow.id:
-            reg_map[1][id] = regions[bc_outflow_id]
-
-
-    return common_decomp, reg_map
-
-
-def fill_lg(decomp, reg_map, regions):
+
+
+def make_decomposition(root_polygon_points, fractures, regions, i_r_bulk, i_r_side, i_r_frac, tol):
+    # Create boundary polygon
+    box_pd = poly.PolygonDecomposition([regions[0], regions[0], regions[0]], tol)
+    last_pt = root_polygon_points[-1]
+    side_segments = {}
+    for i_side, pt in enumerate(root_polygon_points):
+        sub_segments = box_pd.add_line(last_pt, pt, attr=regions[i_r_side[i_side]])
+        last_pt = pt
+        assert type(sub_segments) == list and len(sub_segments) == 1
+        seg = sub_segments[0]
+        side_segments[seg.id] = i_side
+    assert len(box_pd.polygons) == 2
+    box_pd.polygons[1].attr = regions[i_r_bulk]
+
+    # Add fractures
+    outer_wire = box_pd.outer_polygon.outer_wire.childs
+    assert len(outer_wire) == 1
+    outer_wire = next(iter(outer_wire))
+    for i_fr, (p0, p1) in enumerate(fractures):
+        box_pd.decomp.check_consistency()
+        print(i_fr, "fr size:", np.linalg.norm(p1 - p0))
+
+        segments = box_pd.add_line(p0, p1, attr=regions[i_r_frac[i_fr]])
+
+        if type(segments) == list:
+            for seg in segments:
+                if seg.wire[0] == seg.wire[1] and seg.wire[0] == outer_wire:
+                    points = seg.vtxs
+                    box_pd.delete_segment(seg)
+                    for pt in points:
+                        if pt.is_free():
+                            box_pd.remove_free_point(pt.id)
+
+        # TODO: remove outer segments
+
+
+    #common_decomp, maps = merge.intersect_decompositions(decompositions)
+    plot_polygon_decomposition(box_pd)
+    return box_pd
+
+
+def fill_lg(decomp, regions, mesh_base="fractured_2d"):
     """
     Create LayerGeometry object.
     """
@@ -112,9 +113,9 @@ def fill_lg(decomp, reg_map, regions):
     layer = lg.FractureLayer(dict(
         name = "layer",
         top = iface_ns,
-        polygon_region_ids = [ reg_map[2][poly.id]._id for poly in decomp.polygons.values() ],
-        segment_region_ids = [ reg_map[1][seg.id]._id for seg in decomp.segments.values() ],
-        node_region_ids = [ reg_map[0][node.id]._id for node in decomp.points.values() ]
+        polygon_region_ids = [ poly.attr._id for poly in decomp.polygons.values() ],
+        segment_region_ids = [ seg.attr._id for seg in decomp.segments.values() ],
+        node_region_ids = [ node.attr._id for node in decomp.points.values() ]
     ))
     geom.layers = [ layer ]
     #geom.surfaces = [ClassFactory(Surface)]
@@ -132,9 +133,9 @@ def fill_lg(decomp, reg_map, regions):
         nodes = nodes
     ))
     geom.node_sets = [ nodeset ]
-    geomop.layers_io.write_geometry("fractured_2d.json", geom)
+    #geomop.layers_io.write_geometry(mesh_base + ".json", geom)
     return geom
 
 
-def make_mesh(geometry):
-    return geomop.geometry.make_geometry(geometry=geometry, layers_file="fractured_2d.json", mesh_step=1.0)
+def make_mesh(geometry, mesh_base="fractured_2d"):
+    return geomop.geometry.make_geometry(geometry=geometry, layers_file=mesh_base + ".json", mesh_step=1.0)

src/geomop/aabb_lookup.py

@@ -61,9 +61,11 @@ def closest_candidates(self, point):
             return []
         inf_dists = np.max(np.maximum(self.boxes[:, 0:2] - point, point - self.boxes[:, 2:4]), axis=1)
         if np.amin(inf_dists) > 0.0:
+            # closest box not containing the point
             i_closest = np.argmin(inf_dists)
             c_boxes = boxes[i_closest:i_closest+1, :]
         else:
+            # point is inside all boxes
             c_boxes = boxes[np.where(inf_dists<=0.0)]
         assert c_boxes.shape[0] != 0
         # Max distance of closest boxes