Tried to fix #54.

Problems related to T-intersection, with interseciton through fracture
vertexes.

examples/example3/test_simple_networks.py

@@ -245,3 +245,21 @@ def test_T_intersection_one_boundary_one_inside(**kwargs):
     domain = {'xmin':-1, 'xmax': 2, 'ymin': -2, 'ymax': 2, 'zmin': -1,
               'zmax':2}
     grids = meshing.simplex_grid([p1, p2], domain)
+
+def test_issue_54():
+    mesh_kwargs = {}
+    mesh_size = 5e-1
+    mesh_kwargs['mesh_size'] = {'mode': 'constant',
+    'value': mesh_size, 'bound_value': 1.2*mesh_size}
+    mesh_kwargs['file_name'] = 'bounding_box_test'
+    domain = {'xmin': -1, 'xmax': 1,
+              'ymin': 0, 'ymax': 1,
+              'zmin': 0, 'zmax': 1}
+    f_1 = np.array([[.5,.5,.5,.5],
+                    [.4,.5,.5,.4],
+                    [0.2,0.2,.8,.8]])
+    f_2 = np.array([[0,.8,.8,0],
+                    [.5,.5,.5,.5],
+                    [0.2,0.2,.8,.8]])
+    grids = meshing.simplex_grid([f_1, f_2], domain,
+                                 ensure_matching_face_cell=False)

src/porepy/fracs/fractures.py

@@ -457,11 +457,11 @@ def point_on_segment(ip, poly):
                 return int_points, on_boundary_self, on_boundary_other
 
         # Else, we have boundary intersection, and need to process them
-        bound_pt_self, self_segment, _, self_cuts_through = \
+        bound_pt_self, self_segment, is_vertex_self, self_cuts_through = \
                 self._process_segment_isect(bound_sect_self_other, self.p, tol,
                                             other.p)
 
-        bound_pt_other, other_segment, _, other_cuts_through = \
+        bound_pt_other, other_segment, is_vertex_other, other_cuts_through = \
                 self._process_segment_isect(bound_sect_other_self, other.p,
                                             tol, self.p)
 
@@ -503,28 +503,45 @@ def point_on_segment(ip, poly):
         if self_segment:
             assert other_segment  # This should be reflexive
             assert bound_pt_self.shape[1] == 2
-            assert np.allclose(bound_pt_self, bound_pt_other)
+
+            # Depending on whether the intersection points are also vertexes
+            # of the polygons, bound_pt_self and other need not be identical.
+            # Return the uniquified combination of the two
+            bound_pt = np.hstack((bound_pt_other, bound_pt_self))
+            bound_pt, _, _ = setmembership.unique_columns_tol(bound_pt,
+                                                               tol=tol)
+
             on_boundary_self = [True, True]
             on_boundary_other = [True, True]
-            return bound_pt_self, on_boundary_self, on_boundary_other
+            return bound_pt, on_boundary_self, on_boundary_other
 
         # Case where one boundary segment of one fracture cuts through two
         # boundary segment (thus the surface) of the other fracture. Gives a
         # T-type intersection
         if self_cuts_through:
-            assert np.allclose(bound_pt_self, bound_pt_other)
             on_boundary_self = True
             on_boundary_other = False
-            return bound_pt_self, on_boundary_self, on_boundary_other
+            # Depending on whether the intersection points are also vertexes
+            # of the polygons, bound_pt_self and other need not be identical.
+            # Return the uniquified combination of the two
+            bound_pt = np.hstack((bound_pt_other, bound_pt_self))
+            bound_pt, _, _ = setmembership.unique_columns_tol(bound_pt,
+                                                               tol=tol)
+
+            return bound_pt, on_boundary_self, on_boundary_other
         elif other_cuts_through:
-            assert np.allclose(bound_pt_self, bound_pt_other)
             on_boundary_self = False
             on_boundary_other = True
-            return bound_pt_self, on_boundary_self, on_boundary_other
 
-        # By now, there should be a single member of bound_pt
-        assert bound_pt_self.shape[1] == 1 or bound_pt_self.shape[1] == 2
-        assert bound_pt_other.shape[1] == 1 or bound_pt_other.shape[1] == 2
+            # Depending on whether the intersection points are also vertexes
+            # of the polygons, bound_pt_self and other need not be identical.
+            # Return the uniquified combination of the two
+            bound_pt = np.hstack((bound_pt_other, bound_pt_self))
+            bound_pt, _, _ = setmembership.unique_columns_tol(bound_pt,
+                                                               tol=tol)
+
+            return bound_pt, on_boundary_self, on_boundary_other
+
         bound_pt = np.hstack((int_points, bound_pt_self))
 
         # Check if there are two boundary points on the same segment. If so,
@@ -617,8 +634,8 @@ def _process_segment_isect(self, isect_bound, poly, tol, other_poly):
             elif ip_unique.shape[1] == 1:
                 # Either a vertex or single intersection point.
                 poly_ext, _, _ = setmembership.unique_columns_tol(
-                    np.hstack((self.p, ip_unique)), tol=tol)
-                if poly_ext.shape[1] == self.p.shape[1]:
+                    np.hstack((poly, ip_unique)), tol=tol)
+                if poly_ext.shape[1] == poly.shape[1]:
                     # This is a vertex, we skip it
                     pass
                 else: