From beaef8b157e59f7ffee8df75d2afc402e34f00c5 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 25 Aug 2017 14:48:18 +0200 Subject: [PATCH 001/118] Bugfix in comp geom method to compute mutual distances in set of segments. --- src/porepy/utils/comp_geom.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index a30c511f16..e9910e7817 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1643,8 +1643,8 @@ def dist_segment_set(start, end): for j in range(i+1, ns): dl, cpi, cpj = dist_two_segments(start[:, i], end[:, i], start[:, j], end[:, j]) - d[i, j] = d - d[j, i] = d + d[i, j] = dl + d[j, i] = dl cp[i, j, :] = cpi cp[j, i, :] = cpj From ab78a6b7d3be9c46b776aefeb7f0c91cfdf8f0a1 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 25 Aug 2017 14:48:58 +0200 Subject: [PATCH 002/118] Method to compute distance between segment and set of segments. In-between version of two segments and segment sets. One of these must be superfluous. --- src/porepy/utils/comp_geom.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index e9910e7817..290e8b4b09 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1652,6 +1652,36 @@ def dist_segment_set(start, end): #------------------------------------------------------------------------------# +def dist_segment_segment_set(start, end, start_set, end_set): + """ Compute distance and closest points between a segment and a set of + segments. + + Parameters: + + """ + if start.size < 4: + start = start.reshape((-1, 1)) + if end.size < 4: + end = end.reshape((-1, 1)) + + nd = start.shape[0] + ns = start_set.shape[1] + + d = np.zeros( ns) + cp_set = np.zeros(( nd, ns)) + cp = np.zeros((nd, ns)) + + for i in range(ns): + dl, cpi, cpj = dist_two_segments(start, end, start_set[:, j], + end_set[:, j]) + d[i] = dl + cp[:, i] = cpi + cp_set[:, i] = cpj + + return d, cp, cp_set + +#------------------------------------------------------------------------------# + def dist_two_segments(s1_start, s1_end, s2_start, s2_end): """ Compute the distance between two line segments. From 330749bfdbf90dbdfe8d12ddaf35bc7aed28e8b2 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 28 Aug 2017 10:26:10 +0200 Subject: [PATCH 003/118] Grid has method to identify cells closest to point of clouds. --- src/porepy/grids/grid.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/porepy/grids/grid.py b/src/porepy/grids/grid.py index c9d0d4a5ef..293cf3446d 100644 --- a/src/porepy/grids/grid.py +++ b/src/porepy/grids/grid.py @@ -693,6 +693,32 @@ def cell_connection_map(self): return c2c + def closest_cell(self, p): + """ For a set of points, find closest cell by cell center. + + Parameters: + p (np.ndarray, 3xn): Point coordinates. If p.shape[0] < 3, + additional points will be treated as zeros. + + Returns: + np.ndarray of ints: For each point, index of the cell with center + closest to the point. + """ + dim_p = p.shape[0] + if p.shape[0] < 3: + z = np.zeros((3 - p.shape[0], p.shape[1])) + p = np.vstack((p, z)) + + def min_dist(pts): + c = self.cell_centers + d = np.sum(np.power(c - pts, 2), axis=0) + return np.argmin(d) + + ci = np.empty(p.shape[1], dtype=np.int) + for i in range(p.shape[1]): + ci[i] = min_dist(p[:, i].reshape((3, -1))) + return ci + def add_face_tag(self, f, tag): self.face_tags[f] = np.bitwise_or(self.face_tags[f], tag) From fa2e4e7077d4ee5171061ca13fc89e7d84a83cd8 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 29 Aug 2017 12:13:05 +0200 Subject: [PATCH 004/118] Minor bugfix in biot discretization --- src/porepy/numerics/fv/biot.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/porepy/numerics/fv/biot.py b/src/porepy/numerics/fv/biot.py index 31a0237b70..edaa3d83d5 100644 --- a/src/porepy/numerics/fv/biot.py +++ b/src/porepy/numerics/fv/biot.py @@ -496,7 +496,7 @@ def slv(b): x = la.spsolve(A, b) return x elif solver == 'factorized': - slv = la.factorized(A.to_csc()) + slv = la.factorized(A.tocsc()) else: raise ValueError('Unknown solver ' + solver) @@ -587,9 +587,10 @@ def compute_stress(self, g, u, data): all stress values on the first face, then the second etc. """ + param = data['param'] stress_discr = data['stress'] bound_stress = data['bound_stress'] - bound_val = data['bc_val_mech'] + bound_val = param.get_bc_val_mechanics() d = self.extractD(g, u, as_vector=True) stress = np.squeeze(stress_discr * d) + (bound_stress * bound_val) return stress From b82bd7b8bafa7b90b89e1e7b8d24ddea55a26da4 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 29 Aug 2017 12:13:35 +0200 Subject: [PATCH 005/118] Minor changes to Parameter class --- src/porepy/params/data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/params/data.py b/src/porepy/params/data.py index 2f073cf1b2..2a4c79e9a7 100644 --- a/src/porepy/params/data.py +++ b/src/porepy/params/data.py @@ -478,7 +478,7 @@ def get_bc_transport(self): else: return BoundaryCondition(self.g) - conductivity = property(get_conductivity) + bc_transport = property(get_bc_transport) def get_bc_mechanics(self): """ Stiffness matrix, defined as fourth order tensor @@ -488,7 +488,7 @@ def get_bc_mechanics(self): else: return BoundaryCondition(self.g) - stiffness = property(get_stiffness) + bc_mechanics = property(get_bc_mechanics) # Boundary value From ad5632dae10d267ce132072bbb07e8b43f169860 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 29 Aug 2017 15:18:53 +0200 Subject: [PATCH 006/118] Distance point-polygon also indicate whether closest point is on polygon boundary. --- src/porepy/utils/comp_geom.py | 10 ++++++---- test/unit/test_distance_computations.py | 6 ++++-- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 290e8b4b09..59d3a8107e 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1925,6 +1925,8 @@ def dist_points_polygon(p, poly, tol=1e-5): np.array (n_pts): Distance from points to polygon np.array (nd x n_pts): For each point, the closest point on the polygon. + np.array (n_pts, bool): True if the point is found in the interior, + false if on a bounding segment. """ @@ -1981,7 +1983,7 @@ def dist_points_polygon(p, poly, tol=1e-5): cp[:, in_poly] = center + irot.dot(cp_inpoly) if np.all(in_poly): - return d, cp + return d, cp, in_poly # Next, points that are outside the extruded polygons. These will have # their closest point among one of the edges @@ -1998,7 +2000,7 @@ def dist_points_polygon(p, poly, tol=1e-5): d[pi] = d_outside[i, mi] cp[:, pi] = p_outside[i, mi, :] - return d, cp + return d, cp, in_poly #----------------------------------------------------------------------------# @@ -2109,8 +2111,8 @@ def dist_segments_polygon(start, end, poly, tol=1e-5): end = orig_end # Distance from endpoints to - d_start_poly, cp_s_p = dist_points_polygon(start, poly) - d_end_poly, cp_e_p = dist_points_polygon(end, poly) + d_start_poly, cp_s_p, s_in_poly = dist_points_polygon(start, poly) + d_end_poly, cp_e_p, e_in_poly = dist_points_polygon(end, poly) # Loop over all segments that did not cross the polygon. The minimum is # found either by the endpoints, or as between two segments. diff --git a/test/unit/test_distance_computations.py b/test/unit/test_distance_computations.py index ea67d760f8..5efb4c0917 100644 --- a/test/unit/test_distance_computations.py +++ b/test/unit/test_distance_computations.py @@ -226,11 +226,12 @@ def test_norot_poly(self): p = np.array([[0.5, 0.5, 0], [0.5, 0, 0], [0.5, 0.5, 1], [0, 0, 1], [0.5, 0, 1], [2, 0.5, 0], [2, 0, 1]]).T - d, cp = cg.dist_points_polygon(p, poly) + d, cp, in_poly = cg.dist_points_polygon(p, poly) known_d = np.array([0, 0, 1, 1, 1, 1, np.sqrt(2)]) known_cp = np.array([[0.5, 0.5, 0], [0.5, 0, 0], [0.5, 0.5, 0], [0, 0, 0], [0.5, 0, 0], [1, 0.5, 0], [1, 0, 0]]).T + known_inp = np.array([1, 0, 0, 0, 0, 0], dtype=np.bool) assert np.allclose(d, known_d) assert np.allclose(cp, known_cp) @@ -241,12 +242,13 @@ def test_rot_poly(self): p = np.array([[0, 0, 0], [0, 0.5, 0.5], [2, 0.5, 0.5], [0, 0, 0.5], [0, -1, 0.5], [1, 0, 0], [1, 0.5, 0.5]]).T - d, cp = cg.dist_points_polygon(p, poly) + d, cp, in_poly = cg.dist_points_polygon(p, poly) known_d = np.array([1, 1, 1, 1, np.sqrt(2), 0, 0]) known_cp = np.array([[1, 0, 0], [1, 0.5, 0.5], [1, 0.5, 0.5], [1, 0, 0.5], [1, 0, 0.5], [1, 0, 0], [1, 0.5, 0.5]]).T + known_inp = np.array([0, 1, 0, 0, 0, 0, 1], dtype=np.bool) assert np.allclose(d, known_d) assert np.allclose(cp, known_cp) From 20b027e8c2ee1666b7b89db3350aea674918b54b Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 4 Sep 2017 17:25:35 +0200 Subject: [PATCH 007/118] Minor update to extrusion module --- src/porepy/fracs/extrusion.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index d63dc7b5fe..fd3e8126a4 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -333,6 +333,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point): # Eliminate points that are on the other side. eliminate = np.where(sgn * right_sign < 0) + main_frac = np.delete(main_frac, eliminate, axis=1) # Add intersection points on the main fracture. One of these may already be # present, as the point of extrusion, but add_point will uniquify the point From 9a6d49684f843525be9ce652dbf952ef12e513f4 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 4 Sep 2017 17:25:54 +0200 Subject: [PATCH 008/118] Minor update to simplex fracture gridding --- src/porepy/fracs/simplex.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/simplex.py b/src/porepy/fracs/simplex.py index 4ffe2c2cbf..0d8aa33e5a 100644 --- a/src/porepy/fracs/simplex.py +++ b/src/porepy/fracs/simplex.py @@ -6,7 +6,6 @@ import numpy as np from meshio import gmsh_io -import porepy from porepy.grids import constants from porepy.grids.gmsh import gmsh_interface, mesh_2_grid from porepy.fracs import fractures, utils @@ -235,8 +234,8 @@ def triangle_grid(fracs, domain, tol=1e-4, **kwargs): # Gridding size if 'mesh_size' in kwargs.keys(): mesh_size, mesh_size_bound, pts_split, lines_split = \ - utils.determine_mesh_size(pts_split, lines_split, - **kwargs['mesh_size']) + utils.determine_mesh_size_2d(pts_split, lines_split, + **kwargs['mesh_size']) else: mesh_size = None mesh_size_bound = None From 93abd3a2137e1ae9ec7eeec26936f1cb4547d604 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 4 Sep 2017 17:26:59 +0200 Subject: [PATCH 009/118] Continued working on fracture network point tolerance in 3d --- src/porepy/fracs/fractures.py | 168 +++++++++++++++++++++++++++++----- 1 file changed, 147 insertions(+), 21 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 0aba8d430e..2bd102e952 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -98,7 +98,7 @@ def points_2_ccw(self): """ Ensure that the points are sorted in a counter-clockwise order. - Implementation note: + mplementation note: For now, the ordering of nodes in based on a simple angle argument. This will not be robust for general point clouds, but we expect the fractures to be regularly shaped in this sense. In particular, we @@ -847,25 +847,25 @@ def get_normal(self, frac): def get_center(self, frac): return self._fractures[frac].center - def get_intersections(self, frac): - """ - Get all intersections in the plane of a given fracture. + def intersections_of_fracture(self, frac): + """ Get all known intersections for a fracture. + + If called before find_intersections(), the returned list will be empty. + + Paremeters: + frac: A fracture in the network + + Returns: + np.array (Intersection): Array of intersections + """ - if frac is None: - frac = np.arange(len(self._fractures)) - if isinstance(frac, np.ndarray): - frac = list(frac) - elif isinstance(frac, int): - frac = [frac] - isect = [] - for fi in frac: - f = self[fi] - isect_loc = [] - for i in self.intersections: - if i.first.index == f.index or i.second.index == f.index: - isect_loc.append(i) - isect.append(isect_loc) - return isect + fi = frac.index + frac_arr = [] + for i in self.intersections: + if i.first == frac or i.second == frac: + ac_arr.append(i) + return np.asarray(arr).astype(np.int) + def find_intersections(self, use_orig_points=False): """ @@ -1019,6 +1019,22 @@ def split_intersections(self): logger.info('Finished fracture splitting after %.5f seconds', time.time() - start_time) + def _fracs_2_edges(self, edges_2_frac): + """ Invert the mapping between edges and fractures. + + Returns: + List of list: For each fracture, index of all edges that points to + it. + """ + f2e = [] + for fi in len(self._fractures): + f_l = [] + for ei, e in enumerate(edges_2_frac): + if fi in e: + f_l.append(ei) + f2e.append(f_l) + return f2e + def _point_and_edge_lists(self): """ Obtain lists of all points and connections necessary to describe @@ -2012,10 +2028,120 @@ def _determine_mesh_size(self, **kwargs): else: raise ValueError('Unknown mesh size mode ' + mode) - def distance_point_segment(): + def compute_distances(self, h_ideal, hmin): + + isect_pt = np.zeros((3, 0), dtype=np.double) + + def dist_p(a, b): + a = a.reshape((-1, 1)) + b = a.reshape((-1, 1)) + d = b - a + return np.sqrt(np.sum(d * d)) + + + edges = self.decomposition['edges'] + + frac_2_edge = self._fracs_2_edges(self.decomposition['edges_2_frac']) + + intersecting_frac = [] + # Loop over all fractures + for f in self._fractures: + + # First compare segments with intersections to this fracture + nfp = f.p.shape[1] + + # Keep track of which other fractures are intersecting - will be + # needed later on + isect_f = [] + for i in self.intersections_of_fracture(f): + if f is i.first: + isect_f.append(i.second.index) + else: + isect_f.append(i.first.index) + + # Assuming the fracture is convex, the closest point for + dist, cp = cg.dist_points_segments(i.coord, f.p, + np.roll(f.p, 1, axis=1)) + # Insert a (candidate) point only at the segment closest to the + # intersection point. If the intersection line runs parallel + # with a segment, this may be insufficient, but we will deal + # with this if necessary. + closest_segment = np.argmin(dist, axis=1) + + for pi, si, di in enumerate(zip(closest_segment, + dist[closest_segment])): + if di < h_ideal: + d_1 = dist_p(cp[:, pi], f.p[:, si]) + d_2 = dist_p(cp[:, pi], f.p[:, (si+1)%nfp]) + if d_1 > hmin and d_2 > hmin: + np.insert(f.p, (si+1)%nfp, cp[:, pi], axis=1) + intersecting_frac.append(isect_f) + + for fi, f in enumerate(self._fractures): + nfp = fi.p.shape[1] + for of in self._fractures: + # Can do some box arguments here to avoid computations + + # First, check if we are intersecting, this is covered already + if of in intersecting_fracs[fi]: + continue + + # Then, compare distance between segments + for si, fs in enumerate(f.segments()): + of_start = of.p + of_end = np.roll(of_start, 1, axis=1) + d, cp_f, _ = cg.dist_segment_segment_set(f_start, f_end, + of_start, of_end) + mi = np.argmin(d) + if d[mi] < h_ideal: + d_1 = dist_p(cp_f[:, mi], f.p[:, si]) + d_2 = dist_p(cp_f[:, mi], f.p[:, (si+1)%nfp]) + if d_1 > hmin and d_2 > hmin: + np.insert(f.p, (si+1)%nfp, cp[:, pi], axis=1) + + # Finally, cover the case where the smallest distance is given + # by a point. Points with false in_point should already be + # covered by the above iteration over segments. + d, cp, in_poly = cg.dist_points_polygon(of.p, f.p) + for di, cpi, ip in zip(d, cp, in_poly): + # Closest point on segment is covered above + if not ip: + continue + if di < h_ideal: + pass + + + + intersecting_frac = np.unique(intersecting_frac) + + f_start = f.p + f_end = np.roll(f_start, 1, axis=1) + + # Loop over all fractures once more. + for of in self._fractures: + # If the fractures are intersecting, no need to do anything + # (covered above). If not, find the closest points by comparing + # segments. + if of.index in intersecting_frac: + continue + + # again information to be picked up - identify the distance + # with the closest points, and the two respective fractures + # / segments + + + + + + + + + + + def distance_point_segment(self): pass - def distance_segment_segment(): + def distance_segment_segment(self): # What happens if two segments are close? pass From db540049e9847ee3440ad864ec76b84ed2efb05b Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 08:59:49 +0200 Subject: [PATCH 010/118] Minor bugfix in fracture extrusion --- src/porepy/fracs/extrusion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 1dfc50dbc3..d0d7e025b8 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -22,8 +22,8 @@ def _intersection_by_num_node(edges, num): num_crosses = crosses.size - edges_of_crosses = np.zeros((num_abut, num), dtype=np.int) - for i, pi in enumerate(abutments): + edges_of_crosses = np.zeros((num_crosses, num), dtype=np.int) + for i, pi in enumerate(crosses): edges_of_crosses[i] = np.where(np.any(edges[:2] == pi, axis=0))[0] return crosses, edges_of_crosses From 999a8b3ce4b912fb4e8a08a8dd76d7211869ca8f Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 09:22:22 +0200 Subject: [PATCH 011/118] Improved sympy representation of fracture polygons. --- src/porepy/fracs/fractures.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index a22a0b9061..771d60fd3d 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -222,9 +222,22 @@ def compute_centroid(self): # Project back again. self.center = rot.transpose().dot(np.append(center, z)) - def as_sp_polygon(self): - sp = [sympy.geometry.Point(self.p[:, i]) - for i in range(self.p.shape[1])] + def as_sp_polygon(self, p=None): + """ Represent polygon as a sympy object. + + Parameters: + p (np.array, nd x npt, optional): Points for the polygon. Defaults + to None, in which case self.p is used. + + Returns: + sympy.geometry.Polygon: Representation of the polygon formed by p. + + """ + if p is None: + p = self.p + + sp = [sympy.geometry.Point(p[:, i]) + for i in range(p.shape[1])] return sympy.geometry.Polygon(*sp) def intersects(self, other, tol, check_point_contact=True): From 74bae6731da3f23f3943922d9fa9b5f964282436 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 09:23:37 +0200 Subject: [PATCH 012/118] Fracture convexity check on 2d representation of polygon. Previous version used 3d, and ran into trouble with sympy --- src/porepy/fracs/fractures.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 771d60fd3d..9cf36f2455 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -191,7 +191,8 @@ def check_convexity(self): boolean, true if the polygon is convex. """ - return self.as_sp_polygon().is_convex() + p_2d = self.plane_coordinates() + return self.as_sp_polygon(p_2d).is_convex() def compute_centroid(self): """ From a51015eb243288dd572a6bcd19733d5031f9cbd3 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 09:24:13 +0200 Subject: [PATCH 013/118] Bugfix in add_point() of Fracture. Sorting of vertices to ccw did not update vertex order. --- src/porepy/fracs/fractures.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 9cf36f2455..76e3c97f5f 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -142,7 +142,9 @@ def points_2_ccw(self): def add_points(self, p, check_convexity=True, tol=1e-4): """ - Add a point to the polygon, and enforce ccw sorting. + Add a point to the polygon with ccw sorting enforced. + + Test for convexity after points are added. Parameters: p (np.ndarray, 3xn): Points to add @@ -157,7 +159,8 @@ def add_points(self, p, check_convexity=True, tol=1e-4): self.p = np.hstack((self.p, p)) self.p, _, _ = setmembership.unique_columns_tol(self.p, tol=tol) - self.points_2_ccw() + # Sort points to ccw + self.p = self.p[:, self.points_2_ccw()] return self.check_convexity() From fff71ab4b15717df6a53f187e5d14b2539f25968 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 11:38:32 +0200 Subject: [PATCH 014/118] Bugfix in fracture extrusion, t-intersection treatment. --- src/porepy/fracs/extrusion.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index d0d7e025b8..e9ab7a1567 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -59,22 +59,28 @@ def t_intersections(edges): num_abut = abutments.size primal_frac = np.zeros(num_abut, dtype=np.int) sec_frac = np.zeros(num_abut, dtype=np.int) - for i, ei in enumerate(frac_num[edges_of_abutments]): - fi, count = np.unique(ei, return_counts=True) + other_point = np.zeros(num_abut, dtype=np.int) + for i, (pi, ei) in enumerate(zip(abutments, edges_of_abutments)): + # Count number of occurences for each fracture associated with this + # intersection. + fi_all = frac_num[edges_of_abutments] + fi, count = np.unique(fi_all, return_counts=True) assert fi.size == 2 + # Find the fracture number associated with main and abutting edge. if count[0] == 1: primal_frac[i] = fi[1] sec_frac[i] = fi[0] else: primal_frac[i] = fi[0] sec_frac[i] = fi[1] - - other_point = np.zeros(num_abut, dtype=np.int) - for i, pi in enumerate(abutments): - if edges[0, sec_frac[i]] == pi: - other_point[i] = edges[1, sec_frac[i]] + # Also find the other point of the abutting edge + ind = np.where(fi_all == sec_frac[i])[1] + ei_abut = ei[ind] + assert ei_abut.size == 1 + if edges[0, ei_abut] == pi: + other_point[i] = edges[1, ei_abut] else: - other_point[i] = edges[0, sec_frac[i]] + other_point[i] = edges[0, ei_abut] return abutments, primal_frac, sec_frac, other_point From 3cf4754016095d18752dc29bd04722e27b7838a2 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 11:38:48 +0200 Subject: [PATCH 015/118] Apparently working method to cut fracture discs by abutting relations. Used in fracture extrusion. --- src/porepy/fracs/extrusion.py | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index e9ab7a1567..a34db156a1 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -284,7 +284,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): main_min = main_frac.p.min(axis=1) main_max = main_frac.p.max(axis=1) - # Equation for the plane through the other fracture, on the form + # Equation for the plane through the other fracture, on the form # n_x(x-c_x) + n_y(y-c_y) + n_z(z-c_z) = 0 n = cg.compute_normal(other_frac.p).reshape((-1, 1)) c = other_frac.center @@ -296,7 +296,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): # Define points in the plane of the second fracture with min and max # coordinates picked from the main fracture. # The below tricks with indices are needed to find a dimension with a - # non-zero gradient of the plane, so that we can divide safely. + # non-zero gradient of the plane, so that we can divide safely. # Implementation note: It might have been possible to do this with a # rotation to the natural plane of the other fracture, but it is not clear # this will really be simpler. @@ -305,7 +305,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): # We should perhaps do a scaling of coordinates. non_zero = np.where(np.abs(n) > 1e-8)[0] if non_zero.size == 0: - raise ValueError('Points of second fracture too close') + raise ValueError('Could not compute normal vector of other fracture') ind = np.setdiff1d(np.arange(3), non_zero[0]) i0 = ind[0] i1 = ind[1] @@ -315,23 +315,23 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): # A for-loop might have been possible here. p[i0, 0] = main_min[i0] p[i1, 0] = main_min[i1] - p[i2, 0] = c[i0] - (n[i0] * (main_min[i0] - c[i0]) - + n[i1] * (main_min[i1] - c[i1])) / n[i0] + p[i2, 0] = c[i2] - (n[i0] * (main_min[i0] - c[i0]) + + n[i1] * (main_min[i1] - c[i1])) / n[i2] p[i0, 1] = main_max[i0] p[i1, 1] = main_min[i1] - p[i2, 1] = c[i0] - (n[i0] * (main_max[i0] - c[i0]) - + n[i1] * (main_min[i1] - c[i1])) / n[i0] + p[i2, 1] = c[i2] - (n[i0] * (main_max[i0] - c[i0]) + + n[i1] * (main_min[i1] - c[i1])) / n[i2] p[i0, 2] = main_max[i0] p[i1, 2] = main_max[i1] - p[i2, 2] = c[i0] - (n[i0] * (main_max[i0] - c[i0]) - + n[i1] * (main_max[i1] - c[i1])) / n[i0] + p[i2, 2] = c[i2] - (n[i0] * (main_max[i0] - c[i0]) + + n[i1] * (main_max[i1] - c[i1])) / n[i2] p[i0, 3] = main_min[i0] p[i1, 3] = main_max[i1] - p[i2, 3] = c[i0] - (n[i0] * (main_min[i0] - c[i0]) - + n[i1] * (main_max[i1] - c[i1])) / n[i0] + p[i2, 3] = c[i2] - (n[i0] * (main_min[i0] - c[i0]) + + n[i1] * (main_max[i1] - c[i1])) / n[i2] # Create an auxiliary fracture that spans the same plane as the other # fracture, and with a larger extension than the main fracture. @@ -342,11 +342,13 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): # Next step is to eliminate points in the main fracture that are on the # wrong side of the other fracture. v = main_frac.p - other_frac.center.reshape((-1, 1)) - sgn = np.sign(v * n) - right_sign = np.sign(reference_point * v) + sgn = np.sign(np.sum(v * n, axis=0)) + ref_v = reference_point - other_frac.center.reshape((-1, 1)) + right_sign = np.sign(np.sum(ref_v * n, axis=0)) - # Eliminate points that are on the other side. - eliminate = np.where(sgn * right_sign < 0) + # Eliminate points that are on the other side. + eliminate = np.where(sgn * right_sign < 0)[0] + main_frac.remove_points(eliminate) # Add intersection points on the main fracture. One of these may already be # present, as the point of extrusion, but add_point will uniquify the point From eeb44f25d17654b8a33e3f944f4967cbbc42fe4e Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 11:53:39 +0200 Subject: [PATCH 016/118] Fracture class has method to check if its polygon is planar. --- src/porepy/fracs/fractures.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 76e3c97f5f..91789ee010 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -58,6 +58,9 @@ def __init__(self, points, index=None): self.index = index + assert self.is_planar(), 'Points define non-planar fracture' + assert self.check_convexity(), 'Points form non-convex polygon' + def set_index(self, i): self.index = i @@ -197,6 +200,21 @@ def check_convexity(self): p_2d = self.plane_coordinates() return self.as_sp_polygon(p_2d).is_convex() + def is_planar(self, tol=1e-4): + """ Check if the points forming this fracture lies in a plane. + + Parameters: + tol (double): Tolerance for non-planarity. Treated as an absolute + quantity (no scaling with fracture extent) + + Returns: + boolean, True if the polygon is planar. False if not. + """ + p = self.p - np.mean(self.p, axis=1).reshape((-1, 1)) + rot = cg.project_plane_matrix(p) + p_2d = rot.dot(p) + return np.max(np.abs(p_2d[2])) < tol + def compute_centroid(self): """ Compute, and redefine, center of the fracture in the form of the From ee0a8cd5ff939ff007e09bc5b84511fb7683f9e1 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 11:54:11 +0200 Subject: [PATCH 017/118] Fracture class has method to remove points. --- src/porepy/fracs/fractures.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 91789ee010..7179902022 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -165,7 +165,20 @@ def add_points(self, p, check_convexity=True, tol=1e-4): # Sort points to ccw self.p = self.p[:, self.points_2_ccw()] - return self.check_convexity() + return self.check_convexity() and self.is_planar(tol) + + def remove_points(self, ind, keep_orig=False): + """ Remove points from the fracture definition + + Parameters: + ind (np array-like): Indices of points to remove. + keep_orig (boolean, optional): Whether to keep the original points + in the attribute orig_p. Defaults to False. + + """ + self.p = np.delete(self.p, ind, axis=1) + if not keep_orig: + self.orig_p = self.p def plane_coordinates(self): """ From 8af5f93ef4d19079e9ab9af84ba50d24de6a0a0f Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 12:47:35 +0200 Subject: [PATCH 018/118] Fracture test of convexity is optional --- src/porepy/fracs/extrusion.py | 6 +++--- src/porepy/fracs/fractures.py | 17 ++++++++++++----- 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index a34db156a1..ea7e86b5ad 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -204,7 +204,7 @@ def discs_from_exposure(pt, edges): for i in range(num_fracs): # The simplest way of distributing points along the disc seems to be to - # create an elliptic fracture, and pick out the points. + # create an elliptic fracture, and pick out the points. f = EllipticFracture(center=center[:, i], major_axis=radius[i], minor_axis=radius[i], dip_angle=np.pi/2, strike_angle=strike_angle[i], major_axis_angle=0) @@ -214,7 +214,7 @@ def discs_from_exposure(pt, edges): # have on the fracture f.add_points(np.vstack((np.hstack((p0[:, i], 0)), np.hstack((p1[:, i], 0)))).T) - fracs.append(Fracture(f.p)) + fracs.append(Fracture(f.p, check_convexity=False)) return fracs @@ -335,7 +335,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): # Create an auxiliary fracture that spans the same plane as the other # fracture, and with a larger extension than the main fracture. - aux_frac = Fracture(p) + aux_frac = Fracture(p, check_convexity=False) isect_pt, _, _ = main_frac.intersects(aux_frac, tol) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 0cedcc9eb0..dde0712837 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -47,7 +47,7 @@ class Fracture(object): - def __init__(self, points, index=None): + def __init__(self, points, index=None, check_convexity=True): self.p = points # Ensure the points are ccw self.points_2_ccw() @@ -59,7 +59,8 @@ def __init__(self, points, index=None): self.index = index assert self.is_planar(), 'Points define non-planar fracture' - assert self.check_convexity(), 'Points form non-convex polygon' + if check_convexity: + assert self.check_convexity(), 'Points form non-convex polygon' def set_index(self, i): self.index = i @@ -147,12 +148,15 @@ def add_points(self, p, check_convexity=True, tol=1e-4): """ Add a point to the polygon with ccw sorting enforced. - Test for convexity after points are added. + Always run a test to check that the points are still planar. By + default, a check of convexity is also performed, however, this can be + turned off to speed up simulations (the test uses sympy, which turns + out to be slow in many cases). Parameters: p (np.ndarray, 3xn): Points to add check_convexity (boolean, optional): Verify that the polygon is - convex. Defaults to true + convex. Defaults to true. tol (double): Tolerance used to check if the point already exists. Return: @@ -165,7 +169,10 @@ def add_points(self, p, check_convexity=True, tol=1e-4): # Sort points to ccw self.p = self.p[:, self.points_2_ccw()] - return self.check_convexity() and self.is_planar(tol) + if check_convexity: + return self.check_convexity() and self.is_planar(tol) + else: + return self.check_convexity() def remove_points(self, ind, keep_orig=False): """ Remove points from the fracture definition From f007f5c8f243b9e9a62532eda2db91ee3d778e5e Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 15:05:11 +0200 Subject: [PATCH 019/118] Minor fixes of segment distance computations in comp_geom --- src/porepy/utils/comp_geom.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index de2e5b878a..63940af430 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1678,10 +1678,8 @@ def dist_segment_segment_set(start, end, start_set, end_set): Parameters: """ - if start.size < 4: - start = start.reshape((-1, 1)) - if end.size < 4: - end = end.reshape((-1, 1)) + start = np.squeeze(start) + end = np.squeeze(end) nd = start.shape[0] ns = start_set.shape[1] @@ -1691,8 +1689,8 @@ def dist_segment_segment_set(start, end, start_set, end_set): cp = np.zeros((nd, ns)) for i in range(ns): - dl, cpi, cpj = dist_two_segments(start, end, start_set[:, j], - end_set[:, j]) + dl, cpi, cpj = dist_two_segments(start, end, start_set[:, i], + end_set[:, i]) d[i] = dl cp[:, i] = cpi cp_set[:, i] = cpj @@ -1746,7 +1744,7 @@ def dist_two_segments(s1_start, s1_end, s2_start, s2_end): dot_2_starts = d2.dot(d_starts) discr = dot_1_1 * dot_2_2 - dot_1_2 ** 2 # Sanity check - assert discr >= 0 + assert discr >= -SMALL_TOLERANCE sc = sN = sD = discr tc = tN = tD = discr From 5c572adc98f214532c0f082b7c3e1dd9e1a09122 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 15:10:04 +0200 Subject: [PATCH 020/118] Made test of fracture planarity relative to fracture size. Still not sure if this is a good option, but it should reduce the number of error messages here. --- src/porepy/fracs/fractures.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index dde0712837..162be6f6ec 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1910,7 +1910,9 @@ def _points_2_plane(self, p_loc, edges_loc, p_ind_loc): rot = cg.project_plane_matrix(p_loc) p_2d = rot.dot(p_loc) - assert np.max(np.abs(p_2d[2])) < 2*self.tol * np.sqrt(3) + extent = p_2d.max(axis=1) - p_2d.min(axis=1) + + assert extent[2] < np.max(extent[:2]) * self.tol # Dump third coordinate p_2d = p_2d[:2] From aa724d5a43c4bd3f3ad616ce788d67e6029d5033 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 15:11:09 +0200 Subject: [PATCH 021/118] Calculate bounding box for fracture network if None is provided. --- src/porepy/fracs/fractures.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 162be6f6ec..12e1bdd91f 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1975,12 +1975,13 @@ def compute_distances(self): return p_2_p, s_2_s - def impose_external_boundary(self, box, truncate_fractures=True): + def impose_external_boundary(self, box=None, truncate_fractures=True): """ Set an external boundary for the fracture set. The boundary takes the form of a 3D box, described by its minimum and - maximum coordinates. + maximum coordinates. If no bounding box is provided, a box will be + fited outside the fracture network. If desired, the fratures will be truncated to lay within the bounding box; that is, Fracture.p will be modified. The orginal coordinates of @@ -1998,6 +1999,20 @@ def impose_external_boundary(self, box, truncate_fractures=True): boundary will be truncated. """ + if box is None: + OVERLAP = 0.05 + cmin = np.ones((3, 1)) * float('inf') + cmax = -np.ones((3, 1)) * float('inf') + for f in self._fractures: + cmin = np.min(np.hstack((cmin, f.p)), axis=1).reshape((-1, 1)) + cmax = np.max(np.hstack((cmax, f.p)), axis=1).reshape((-1, 1)) + cmin = cmin[:, 0] + cmax = cmax[:, 0] + dx = OVERLAP * (cmax - cmin) + box = {'xmin': cmin[0] - dx[0], 'xmax': cmax[0] + dx[0], + 'ymin': cmin[1] - dx[1], 'ymax': cmax[1] + dx[1], + 'zmin': cmin[2] - dx[2], 'zmax': cmax[2] + dx[2]} + # Insert boundary in the form of a box, and kick out (parts of) # fractures outside the box self.domain = box From 2180f2c349c18cab601384bf2a81a59a67958b2b Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 15:12:48 +0200 Subject: [PATCH 022/118] Bugfix in method to find intersecting fractures in a FractureNetwork. --- src/porepy/fracs/fractures.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 12e1bdd91f..5c63f7967f 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -971,9 +971,11 @@ def intersections_of_fracture(self, frac): fi = frac.index frac_arr = [] for i in self.intersections: + if i.coord.size == 0: + continue if i.first == frac or i.second == frac: - ac_arr.append(i) - return np.asarray(arr).astype(np.int) + frac_arr.append(i) + return frac_arr def find_intersections(self, use_orig_points=False): From d96f7f11e5a5b7220e07164d145b375cfa2beb61 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 15:13:37 +0200 Subject: [PATCH 023/118] Several changes to distance computations in fracture network. --- src/porepy/fracs/fractures.py | 62 +++++++++++++---------------------- 1 file changed, 22 insertions(+), 40 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 5c63f7967f..60346dc893 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -2182,12 +2182,7 @@ def dist_p(a, b): d = b - a return np.sqrt(np.sum(d * d)) - - edges = self.decomposition['edges'] - - frac_2_edge = self._fracs_2_edges(self.decomposition['edges_2_frac']) - - intersecting_frac = [] + intersecting_fracs = [] # Loop over all fractures for f in self._fractures: @@ -2211,32 +2206,46 @@ def dist_p(a, b): # with a segment, this may be insufficient, but we will deal # with this if necessary. closest_segment = np.argmin(dist, axis=1) + min_dist = dist[np.arange(i.coord.shape[1]), closest_segment] - for pi, si, di in enumerate(zip(closest_segment, - dist[closest_segment])): + for pi, (si, di) in enumerate(zip(closest_segment, min_dist)): if di < h_ideal: d_1 = dist_p(cp[:, pi], f.p[:, si]) d_2 = dist_p(cp[:, pi], f.p[:, (si+1)%nfp]) + # If the intersection point is not very close to any of + # the points on the segment, we split the segment. if d_1 > hmin and d_2 > hmin: np.insert(f.p, (si+1)%nfp, cp[:, pi], axis=1) - intersecting_frac.append(isect_f) + + # Take note of the intersecting fractures + intersecting_fracs.append(isect_f) for fi, f in enumerate(self._fractures): - nfp = fi.p.shape[1] + nfp = f.p.shape[1] for of in self._fractures: # Can do some box arguments here to avoid computations # First, check if we are intersecting, this is covered already - if of in intersecting_fracs[fi]: + if of.index in intersecting_fracs[fi]: continue + f_start = f.p + f_end = np.roll(f_start, 1, axis=1) # Then, compare distance between segments + # Loop over fracture segments of this fracture for si, fs in enumerate(f.segments()): of_start = of.p of_end = np.roll(of_start, 1, axis=1) - d, cp_f, _ = cg.dist_segment_segment_set(f_start, f_end, - of_start, of_end) + # Compute distance to all fractures of the other fracture, + # and find the closest segment on the other fracture + fs = f_start[:, si].squeeze()#.reshape((-1, 1)) + fe = f_end[:, si].squeeze()#.reshape((-1, 1)) + d, cp_f, _ = cg.dist_segment_segment_set(fs, fe, of_start, + of_end) mi = np.argmin(d) + # If the distance is smaller than ideal length, but the + # closets point is not too close to the segment endpoints, + # we add a new point if d[mi] < h_ideal: d_1 = dist_p(cp_f[:, mi], f.p[:, si]) d_2 = dist_p(cp_f[:, mi], f.p[:, (si+1)%nfp]) @@ -2254,33 +2263,6 @@ def dist_p(a, b): if di < h_ideal: pass - - - intersecting_frac = np.unique(intersecting_frac) - - f_start = f.p - f_end = np.roll(f_start, 1, axis=1) - - # Loop over all fractures once more. - for of in self._fractures: - # If the fractures are intersecting, no need to do anything - # (covered above). If not, find the closest points by comparing - # segments. - if of.index in intersecting_frac: - continue - - # again information to be picked up - identify the distance - # with the closest points, and the two respective fractures - # / segments - - - - - - - - - def distance_point_segment(self): pass From 9a93e2d9570a9442cf7f18e6acb9898d3d871986 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 16:03:07 +0200 Subject: [PATCH 024/118] Minor enhancement of distance computation in point cloud --- src/porepy/utils/comp_geom.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 63940af430..8c574b3bbc 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1903,11 +1903,13 @@ def dist_point_pointset(p, pset, exponent=2): #----------------------------------------------------------------------------# -def dist_pointset(p): +def dist_pointset(p, max_diag=False): """ Compute mutual distance between all points in a point set. Parameters: p (np.ndarray, 3xn): Points + max_diag (boolean, defaults to True): If True, the diagonal values will + are set to a large value, rather than 0. Returns: np.array (nxn): Distance between points. @@ -1922,6 +1924,10 @@ def dist_pointset(p): for i in range(n): d[i] = dist_point_pointset(p[:, i], p) + if max_diag: + for i in range(n): + d[i, i] = 2 * np.max(d[i]) + return d #----------------------------------------------------------------------------# From 515a9b1aa523eeb288841de87134ad8c3b4bce27 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 16:03:45 +0200 Subject: [PATCH 025/118] Gmsh interface handles varying point size better. --- src/porepy/grids/gmsh/gmsh_interface.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/porepy/grids/gmsh/gmsh_interface.py b/src/porepy/grids/gmsh/gmsh_interface.py index ae1cd8e835..ea2131204c 100644 --- a/src/porepy/grids/gmsh/gmsh_interface.py +++ b/src/porepy/grids/gmsh/gmsh_interface.py @@ -18,7 +18,7 @@ class GmshWriter(object): def __init__(self, pts, lines, polygons=None, domain=None, nd=None, mesh_size=None, mesh_size_bound=None, line_type=None, intersection_points=None, tolerance=None, edges_2_frac=None, - meshing_algorithm=None): + meshing_algorithm=None, point_size=None): """ :param pts: np.ndarary, Points @@ -42,6 +42,8 @@ def __init__(self, pts, lines, polygons=None, domain=None, nd=None, if domain is not None: self.domain = domain + self.point_size = point_size + # Points that should be decleared physical (intersections between 3 # fractures) self.intersection_points = intersection_points @@ -207,8 +209,8 @@ def __write_points(self): s += 'p' + str(i) + ' = newp; Point(p' + str(i) + ') = ' s += '{' + str(p[0, i]) + ', ' + str(p[1, i]) + ', '\ + str(p[2, i]) - if self.lchar is not None: - s += ', ' + str(self.lchar[i]) + ' };\n' + if self.point_size is not None: + s += ', ' + str(self.point_size[i]) + ' };\n' else: s += '};\n' From e0e5ff6ee1d1f4b0fe9958fe28dac93de9600dfd Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 16:04:35 +0200 Subject: [PATCH 026/118] Parameter tuning in FractureNetwork class. --- src/porepy/fracs/fractures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 60346dc893..e691a6a47d 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1914,7 +1914,7 @@ def _points_2_plane(self, p_loc, edges_loc, p_ind_loc): extent = p_2d.max(axis=1) - p_2d.min(axis=1) - assert extent[2] < np.max(extent[:2]) * self.tol + assert extent[2] < np.max(extent[:2]) * self.tol * 5 # Dump third coordinate p_2d = p_2d[:2] From f4f90e31c5e764807d3a8102ed3845cebb8f8062 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 16:04:57 +0200 Subject: [PATCH 027/118] Deleted unused function in FractureNetwork. --- src/porepy/fracs/fractures.py | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index e691a6a47d..8f226f8662 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1948,34 +1948,6 @@ def plot_fractures(self, ind=None): f.plot_frame(ax) return fig - def compute_distances(self): - nf = len(self._fractures) - - # Distance between vertexes - p_2_p = np.zeros((nf, nf)) - # Distance between segments - s_2_s = np.zeros((nf, nf)) - - for i, f_1 in enumerate(self._fractures): - for j, f_2 in enumerate(self._fractures[i + 1:]): - d = np.Inf - for p_1 in f_1.points(): - for p_2 in f_2.points(): - d = np.minimum(d, cg.dist_point_pointset(p_1, p_2)[0]) - p_2_p[i, i + j + 1] = d - - d = np.Inf - for s_1 in f_1.segments(): - for s_2 in f_2.segments(): - d = np.minimum(d, - cg.distance_segment_segment(s_1[:, 0], - s_1[:, 1], - s_2[:, 0], - s_2[:, 1])) - s_2_s[i, i + j + 1] = d - - return p_2_p, s_2_s - def impose_external_boundary(self, box=None, truncate_fractures=True): """ From 263ca75c6e9339a844fb0be56f683032efedac86 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Oct 2017 16:05:15 +0200 Subject: [PATCH 028/118] Prototype method for adaptive mesh size in 3d fracture networks. Adaptation relative to fracture geometry. --- src/porepy/fracs/fractures.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 8f226f8662..544ff01a9f 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -2125,7 +2125,7 @@ def _determine_mesh_size(self, **kwargs): See the gmsh manual for further details. """ - mode = kwargs.get('mode', 'constant') + mode = kwargs.get('mode', 'distance') num_pts = self.decomposition['points'].shape[1] @@ -2141,11 +2141,22 @@ def _determine_mesh_size(self, **kwargs): else: mesh_size_bound = None return mesh_size, mesh_size_bound + elif mode == 'distance': + + p = self.decomposition['points'] + dist = cg.dist_pointset(p, max_diag=True) + mesh_size = np.min(dist, axis=1) + mesh_size = np.max(mesh_size, self.hmin * np.ones_like(p)) + mesh_size = np.min(mesh_size, self.h_ideal * np.ones_like(p)) + + mesh_size_bound = self.h_ideal else: raise ValueError('Unknown mesh size mode ' + mode) def compute_distances(self, h_ideal, hmin): + self.h_ideal = h_ideal + self.hmin = hmin isect_pt = np.zeros((3, 0), dtype=np.double) def dist_p(a, b): From 6e6a0bfbc7bda7cd65abe9b674163fe1c652d771 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 9 Oct 2017 09:51:18 +0200 Subject: [PATCH 029/118] Minor fix of FractureNetwork --- src/porepy/fracs/fractures.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 544ff01a9f..d8bd018e68 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -2150,6 +2150,7 @@ def _determine_mesh_size(self, **kwargs): mesh_size = np.min(mesh_size, self.h_ideal * np.ones_like(p)) mesh_size_bound = self.h_ideal + return mesh_size, mehs_size_bound else: raise ValueError('Unknown mesh size mode ' + mode) @@ -2343,8 +2344,7 @@ def to_gmsh(self, file_name, **kwargs): determine_mesh_size(self.decomposition['points'], **kwargs['mesh_size']) else: - mesh_size = None - mesh_size_bound = None + mesh_size, mesh_size_bound = self._determine_mesh_size() # The tolerance applied in gmsh should be consistent with the tolerance # used in the splitting of the fracture network. The documentation of From d3c75b0ca963c097739abfb7a259d19e328fafe4 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 09:44:44 +0200 Subject: [PATCH 030/118] Minor bugfixes in FractureNetwork adaptive mesh size. --- src/porepy/fracs/fractures.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index d8bd018e68..cb235eb5ba 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -2146,11 +2146,11 @@ def _determine_mesh_size(self, **kwargs): p = self.decomposition['points'] dist = cg.dist_pointset(p, max_diag=True) mesh_size = np.min(dist, axis=1) - mesh_size = np.max(mesh_size, self.hmin * np.ones_like(p)) - mesh_size = np.min(mesh_size, self.h_ideal * np.ones_like(p)) + mesh_size = np.maximum(mesh_size, self.hmin * np.ones_like(p)) + mesh_size = np.minimum(mesh_size, self.h_ideal * np.ones_like(p)) mesh_size_bound = self.h_ideal - return mesh_size, mehs_size_bound + return mesh_size, mesh_size_bound else: raise ValueError('Unknown mesh size mode ' + mode) @@ -2162,7 +2162,7 @@ def compute_distances(self, h_ideal, hmin): def dist_p(a, b): a = a.reshape((-1, 1)) - b = a.reshape((-1, 1)) + b = b.reshape((-1, 1)) d = b - a return np.sqrt(np.sum(d * d)) From 89a322205c1745cea553a301a1ce656b18ad8dc7 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 10:23:58 +0200 Subject: [PATCH 031/118] Cleanup of mesh size parameters in gmsh export --- src/porepy/grids/gmsh/gmsh_interface.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/porepy/grids/gmsh/gmsh_interface.py b/src/porepy/grids/gmsh/gmsh_interface.py index ea2131204c..d19ee9fd41 100644 --- a/src/porepy/grids/gmsh/gmsh_interface.py +++ b/src/porepy/grids/gmsh/gmsh_interface.py @@ -18,7 +18,7 @@ class GmshWriter(object): def __init__(self, pts, lines, polygons=None, domain=None, nd=None, mesh_size=None, mesh_size_bound=None, line_type=None, intersection_points=None, tolerance=None, edges_2_frac=None, - meshing_algorithm=None, point_size=None): + meshing_algorithm=None): """ :param pts: np.ndarary, Points @@ -36,13 +36,11 @@ def __init__(self, pts, lines, polygons=None, domain=None, nd=None, else: self.nd = nd - self.lchar = mesh_size - self.lchar_bound = mesh_size_bound - if domain is not None: self.domain = domain - self.point_size = point_size + self.mesh_size = mesh_size + self.mesh_size_bound = mesh_size_bound # Points that should be decleared physical (intersections between 3 # fractures) @@ -154,10 +152,10 @@ def __write_boundary_3d(self): zmax = str(self.domain['zmax']) # + ', ' # Add mesh size on boundary points if these are provided - if self.lchar_bound is not None: + if self.mesh_size_bound is not None: zmin += ', ' zmax += ', ' - h = str(self.lchar_bound) + '};' + h = str(self.mesh_size_bound) + '};' else: h = '};' ls = '\n' @@ -209,8 +207,8 @@ def __write_points(self): s += 'p' + str(i) + ' = newp; Point(p' + str(i) + ') = ' s += '{' + str(p[0, i]) + ', ' + str(p[1, i]) + ', '\ + str(p[2, i]) - if self.point_size is not None: - s += ', ' + str(self.point_size[i]) + ' };\n' + if self.mesh_size is not None: + s += ', ' + str(self.mesh_size[i]) + ' };\n' else: s += '};\n' From 63b44be6aeebc45637306bdbfba04be26535798e Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 10:27:18 +0200 Subject: [PATCH 032/118] Various bugfixes in FractureNetwork adaptive mesh size. --- src/porepy/fracs/fractures.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index cb235eb5ba..fdd83fa361 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -2144,10 +2144,11 @@ def _determine_mesh_size(self, **kwargs): elif mode == 'distance': p = self.decomposition['points'] + num_pts = p.shape[1] dist = cg.dist_pointset(p, max_diag=True) mesh_size = np.min(dist, axis=1) - mesh_size = np.maximum(mesh_size, self.hmin * np.ones_like(p)) - mesh_size = np.minimum(mesh_size, self.h_ideal * np.ones_like(p)) + mesh_size = np.maximum(mesh_size, self.hmin * np.ones(num_pts)) + mesh_size = np.minimum(mesh_size, self.h_ideal * np.ones(num_pts)) mesh_size_bound = self.h_ideal return mesh_size, mesh_size_bound @@ -2164,7 +2165,7 @@ def dist_p(a, b): a = a.reshape((-1, 1)) b = b.reshape((-1, 1)) d = b - a - return np.sqrt(np.sum(d * d)) + return np.sqrt(np.sum(d**2)) intersecting_fracs = [] # Loop over all fractures @@ -2194,12 +2195,12 @@ def dist_p(a, b): for pi, (si, di) in enumerate(zip(closest_segment, min_dist)): if di < h_ideal: - d_1 = dist_p(cp[:, pi], f.p[:, si]) - d_2 = dist_p(cp[:, pi], f.p[:, (si+1)%nfp]) + d_1 = dist_p(cp[pi, si], f.p[:, si]) + d_2 = dist_p(cp[pi, si], f.p[:, (si+1)%nfp]) # If the intersection point is not very close to any of # the points on the segment, we split the segment. if d_1 > hmin and d_2 > hmin: - np.insert(f.p, (si+1)%nfp, cp[:, pi], axis=1) + np.insert(f.p, (si+1)%nfp, cp[pi, si], axis=1) # Take note of the intersecting fractures intersecting_fracs.append(isect_f) From 0eb57f7c6ba76ce1d6a1780e74671e480883e1de Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 10:28:39 +0200 Subject: [PATCH 033/118] Deleted old method for splitting fractures into polygons before gmsh output. Method was no longer needed, fracture intersections should rather be embedded into the fracture planes. --- src/porepy/fracs/fractures.py | 330 ---------------------------------- 1 file changed, 330 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index fdd83fa361..4a17b21fe9 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1089,14 +1089,6 @@ def split_intersections(self): if self.verbose > 1: self._verify_fractures_in_plane(all_p, edges, edges_2_frac) - - # With all edges being non-intersecting, the next step is to split the - # fractures into polygons formed of boundary edges, intersection edges - # and auxiliary edges that connect the boundary and intersections. -# all_p, edges, is_boundary_edge, poly_segments, poly_2_frac\ -# = self._split_into_polygons(all_p, edges, edges_2_frac, -# is_boundary_edge) - # Store the full decomposition. self.decomposition = {'points': all_p, 'edges': edges.astype('int'), @@ -1453,328 +1445,6 @@ def _remove_edge_intersections(self, all_p, edges, edges_2_frac, # way that ensures that no polygon lies on both sides of an # intersection edge. - def _split_into_polygons(self, all_p, edges, edges_2_frac, is_boundary_edge): - """ - Split the fracture surfaces into non-intersecting polygons. - - Starting from a description of the fracture polygons and their - intersection lines (which should be non-intersecting), the fractures - are split into sub-polygons which they may share a boundary, but no - sub-polygon intersect the surface of another. This step is necessary - for the fracture network to be ammenable to gmsh. - - The sub-polygons are defined by drawing auxiliary lines between - fracture points. This expands the global set of edges, but not the - points. The auxiliary lines will be sent to gmsh as constraints on the - gridding algorithm. - - TODO: The restriction that the auxiliary lines can only be drawn - between existing lines can create unfortunate constraints in the - gridding, in terms of sharp angles etc. This can be alleviated by - adding additional points to the global description. However, if a point - is added to an intersection edge, this will require an update of all - fractures sharing that intersection. Not sure what the right thing to - do here is. - - Parameters: - all_p (np.ndarray, 3xn): Coordinates of all points used to describe - the fracture polygons, and their intersections. Should be - unique. - edges (np.ndarray, 2xn): Connections between points, formed either - by a fracture boundary, or a fracture intersection. - edges_2_frac (list): For each edge, index of all fractures that - point to the edge. - is_boundary_edge (np.ndarray of bool, size=num_edges): A flag - telling whether the edge is on the boundary of a fracture. - - Returns: - np.ndarray (3xn_pt): Coordinates of all points used to describe the - fractures. - np.ndarray (3xn_edges): Connections between points, formed by - fracture boundaries, fracture intersections, or auxiliary lines - needed to form sub-polygons. The auxiliary lines will be added - towards the end of the array. - np.ndarray (boolean, 3xn_edges_orig): Flag telling if an edge is on - the boundary of a fracture. NOTE: The list is *not* expanded as - auxiliary lines are added; these are known to not be on the - boundary. The legth of this output, relative to the second, can - thus be used to flag auxiliary lines. - list of np.ndarray: Each item contains a sub-polygon, specified by - the global indices of its vertices. - list of int: For each sub-polygon, index of the fracture it forms a - part of. - - Note that for the moment, the first and third outputs are not - modified compared to respective input. This will change if the - splitting is allowed to introduce new additional points. - - """ - logger.info('Split fractures into non-intersecting polygons') - start_time = time.time() - - # For each polygon, list of segments that make up the polygon - poly_segments = [] - # Which fracture is the polygon part of - poly_2_frac = [] - - for fi, _ in enumerate(self._fractures): - - if self.verbose > 1: - print(' Split fracture no ' + str(fi)) - - # Identify the edges associated with this fracture - # It would have been more convenient to use an inverse - # relationship frac_2_edge here, but that would have made the - # update for new edges (towards the end of this loop) more - # cumbersome. - edges_loc_ind = [] - for ei, e in enumerate(edges_2_frac): - if np.any(e == fi): - edges_loc_ind.append(ei) - - edges_loc = np.vstack((edges[:, edges_loc_ind], - np.array(edges_loc_ind))) - p_ind_loc = np.unique(edges_loc[:2]) - p_loc = all_p[:, p_ind_loc] - - p_2d, edges_2d, _, _ = self._points_2_plane(p_loc, edges_loc, - p_ind_loc) - - - # Add a tag to trace the edges during splitting - edges_2d[2] = edges_loc[2] - - # Flag of whether the constraint is an internal boundary (as - # opposed to boundary of the polygon) - internal_boundary = np.where([not is_boundary_edge[i] - for i in edges_2d[2]])[0] - - # Create a dictionary of input information to triangle - triangle_opts = {'vertices': p_2d.transpose(), - 'segments': edges_2d[:2].transpose(), - 'segment_markers': - edges_2d[2].transpose().astype('int32') - } - # Run a constrained Delaunay triangulation. - t = triangle.triangulate(triangle_opts, 'p') - - # Pull out information on the triangulation - segment_markers = \ - np.squeeze(np.asarray(t['segment_markers']).transpose()) - tri = np.sort(np.asarray(t['triangles']).transpose(), axis=0) - p = np.asarray(t['vertices']).transpose() - segments = np.sort(np.asarray(t['segments']).transpose(), axis=0) - - p = cg.snap_to_grid(p, tol=self.tol) - - # It turns out that Triangle sometime creates (almost) duplicate - # points. Our algorithm are based on the points staying, thus these - # need to go. - # The problem is characterized by output from Triangle having more - # points than the input. - if p.shape[1] > p_2d.shape[1]: - - # Identify duplicate points - # Not sure about tolerance used here - p, _, old_2_new = \ - setmembership.unique_columns_tol(p, tol=np.sqrt(self.tol)) - # Update segment and triangle coordinates to refer to unique - # points - segments = old_2_new[segments] - tri = old_2_new[tri] - - # Remove segments with the same start and endpoint - segments_remove = np.argwhere(np.diff(segments, axis=0)[0] == 0) - segments = np.delete(segments, segments_remove, axis=1) - segment_markers = np.delete(segment_markers, segments_remove) - - # Remove segments that exist in duplicates (a, b) and (b, a) - segments.sort(axis=0) - segments, new_2_old_seg, _ = \ - setmembership.unique_columns_tol(segments) - segment_markers = segment_markers[new_2_old_seg] - - # Look for degenerate triangles, where the same point occurs - # twice - tri.sort(axis=0) - degenerate_tri = np.any(np.diff(tri, axis=0) == 0, axis=0) - tri = tri[:, np.logical_not(degenerate_tri)] - - # Remove duplicate triangles, if any - tri, _, _ = setmembership.unique_columns_tol(tri) - - - # We need a connection map between cells. The simplest option was - # to construct a simplex grid, and use the associated function. - # While we were at it, we could have used this to find cell - # centers, cell-face (segment) mappings etc, but for reason, that - # did not happen. - g = simplex.TriangleGrid(p, tri) - g.compute_geometry() - c2c = g.cell_connection_map() - - def cells_of_segments(cells, s): - # Find the cells neighboring - hit_0 = np.logical_or(cells[0] == s[0], cells[0] == - s[1]).astype('int') - hit_1 = np.logical_or( - cells[1] == s[0], cells[1] == s[1]).astype('int') - hit_2 = np.logical_or( - cells[2] == s[0], cells[2] == s[1]).astype('int') - hit = np.squeeze(np.where(hit_0 + hit_1 + hit_2 == 2)) - return hit - - px = p[0] - py = p[1] - # Cell centers (in 2d coordinates) - cell_c = np.vstack((np.mean(px[tri], axis=0), - np.mean(py[tri], axis=0))) - - # The target sub-polygons should have intersections as an external - # boundary. The sub-polygons are grown from triangles on each side of - # the intersection, denoted cw and ccw. - cw_cells = [] - ccw_cells = [] - - # For each internal boundary, find the cells on each side of the - # boundary; these will be assigned to different polygons. For the - # moment, there should be a single cell along each side of the - # boundary, but this may change in the future. - - if internal_boundary.size == 0: - # For non-intersecting fractures, we just need a single seed - cw_cells = [0] - else: - # Full treatment - for ib in internal_boundary: - segment_match = np.squeeze(np.where(segment_markers == - edges_2d[2, ib])) - loc_segments = segments[:, segment_match] - loc_cells = cells_of_segments(tri, loc_segments) - - p_0 = p_2d[:, edges_2d[0, ib]] - p_1 = p_2d[:, edges_2d[1, ib]] - - cw_loc = [] - ccw_loc = [] - - for ci in np.atleast_1d(loc_cells): - if cg.is_ccw_polyline(p_0, p_1, cell_c[:, ci]): - ccw_loc.append(ci) - else: - cw_loc.append(ci) - cw_cells.append(cw_loc) - ccw_cells.append(ccw_loc) - - polygon = np.zeros(tri.shape[1], dtype='int') - counter = 1 - for c in cw_cells + ccw_cells: - polygon[c] = counter - counter += 1 - - # A cell in the partitioning may border more than one internal - # constraint, and thus be assigned two values (one overwritten) in - # the above loop. This should be fine; the corresponding index - # will have no cells assigned to it, but that is taken care of - # below. - - num_iter = 0 - - # From the seeds, assign a polygon index to all triangles - while np.any(polygon == 0): - - num_occurences = np.bincount(polygon) - # Disregard the 0 (which we know is there, since the while loop - # took another round). - order = np.argsort(num_occurences[1:]) - for pi in order: - vec = 0 * polygon - # We know that the values in polygon are mapped to - vec[np.where(polygon == pi + 1)] = 1 - new_val = c2c * vec - not_found = np.where(np.logical_and(polygon == 0, - new_val > 0)) - polygon[not_found] = pi + 1 - - num_iter += 1 - if num_iter >= tri.shape[1]: - # Sanity check, since we're in a while loop. The code - # should terminate long before this, but I have not - # tried to find an exact maximum number of iterations. - raise ValueError('This must be too many iterations') - - logger.debug('Split fracture %i into %i polygons', fi, - np.unique(polygon).size) - - # For each cell, find its boundary - for poly in np.unique(polygon): - tri_loc = tri[:, np.squeeze(np.where(polygon == poly))] - # Special case of a single triangle - if tri_loc.size == 3: - tri_loc = tri_loc.reshape((-1, 1)) - - seg_loc = np.sort(np.hstack((tri_loc[:2], - tri_loc[1:], - tri_loc[::2])), axis=0) - # Really need something like accumarray here.. - unique_s, _, all_2_unique = setmembership.unique_columns_tol( - seg_loc) - bound_segment_ind = np.squeeze(np.where(np.bincount(all_2_unique) - == 1)) - bound_segment = unique_s[:, bound_segment_ind] - boundary = sort_points.sort_point_pairs(bound_segment) - - # Ensure that the boundary is stored ccw - b_points = p_2d[:, boundary[:, 0]] - if not cg.is_ccw_polygon(b_points): - boundary = boundary[::-1] - - # The boundary is expressed in terms of the local point - # numbering. Convert to global numbering. - boundary_glob = p_ind_loc[boundary].astype('int') - poly_segments.append(boundary_glob) - poly_2_frac.append(fi) - # End of processing this sub-polygon - - # End of processing this fracture. - - # We have now split all fractures into non-intersecting polygons. The - # polygon boundary consists of fracture boundaries, fracture - # intersections and auxiliary edges. The latter must be identified, and - # added to the global edge list. - all_edges = np.empty((2, 0), dtype='int') - for p in poly_segments: - all_edges = np.hstack((all_edges, p)) - - # Find edges of the polygons that can be found in the list of existing - # edges - edge_exist, _ = setmembership.ismember_rows(all_edges, edges) - # Also check with the reverse ordering of edge points - edge_exist_reverse, _ = setmembership.ismember_rows(all_edges[::-1], - edges) - # The edge exists if it is found of the orderings - edge_exist = np.logical_or(edge_exist, edge_exist_reverse) - # Where in the list of all edges are the new ones located? - new_edge_ind = np.where([~i for i in edge_exist])[0] - new_edges = all_edges[:, new_edge_ind] - - # If we have found new edges (will most likely happen if there are - # intersecting fractures in the network), these should be added to the - # edge list, in a unique representation. - if new_edges.size > 0: - # Sort the new edges to root out the same edge being found twice - # from different directions - new_edges = np.sort(new_edges, axis=0) - # Find unique representation - unique_new_edges, _, _\ - = setmembership.unique_columns_tol(new_edges) - edges = np.hstack((edges, unique_new_edges)) - - logger.info('Fracture splitting done. Elapsed time %.5f', time.time() - - start_time) - - return all_p, edges, is_boundary_edge, poly_segments, poly_2_frac - def report_on_decomposition(self, do_print=True, verbose=None): """ From f6be8746f3d0389beddbae76f60e412946955dcb Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 10:46:29 +0200 Subject: [PATCH 034/118] fix bad function name in simplex fracture meshing. --- src/porepy/fracs/simplex.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/simplex.py b/src/porepy/fracs/simplex.py index 0d8aa33e5a..4825d55517 100644 --- a/src/porepy/fracs/simplex.py +++ b/src/porepy/fracs/simplex.py @@ -234,8 +234,8 @@ def triangle_grid(fracs, domain, tol=1e-4, **kwargs): # Gridding size if 'mesh_size' in kwargs.keys(): mesh_size, mesh_size_bound, pts_split, lines_split = \ - utils.determine_mesh_size_2d(pts_split, lines_split, - **kwargs['mesh_size']) + utils.determine_mesh_size(pts_split, lines_split, + **kwargs['mesh_size']) else: mesh_size = None mesh_size_bound = None From 0ef2922c6f5f627e6f2ae25a1a20bfbf81e1ff86 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 11:02:20 +0200 Subject: [PATCH 035/118] Simplex meshing of 3d fracture network tries to use new mesh size computation --- src/porepy/fracs/simplex.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/porepy/fracs/simplex.py b/src/porepy/fracs/simplex.py index 4825d55517..62c1bef73a 100644 --- a/src/porepy/fracs/simplex.py +++ b/src/porepy/fracs/simplex.py @@ -11,6 +11,7 @@ from porepy.fracs import fractures, utils import porepy.utils.comp_geom as cg from porepy.utils.setmembership import unique_columns_tol +import porepy def tetrahedral_grid(fracs=None, box=None, network=None, **kwargs): @@ -95,6 +96,14 @@ def tetrahedral_grid(fracs=None, box=None, network=None, **kwargs): else: print('Use existing intersections') + # If fields h_ideal and h_min are provided, try to estimate mesh sizes. + h_ideal = kwargs.get('h_ideal', None) + h_min = kwargs.get('h_min', None) + if h_ideal is not None and h_min is not None: + network.compute_distances(h_ideal, h_min) + # In this case we need to recompute intersection decomposition anyhow. + network.split_intersections() + if not hasattr(network, 'decomposition'): network.split_intersections() else: From 8f310c7628a1d715a18ed3961c7800bf298b94fa Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 11:02:30 +0200 Subject: [PATCH 036/118] Fix bugs related to automatic mesh size computation in FractureNetwork. Issue arised when no parameters were given. --- src/porepy/fracs/fractures.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 4a17b21fe9..2f681c7a37 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -941,6 +941,10 @@ def __init__(self, fractures=None, verbose=0, tol=1e-4): self.tol = tol self.verbose = verbose + # Initialize mesh size parameters as empty + self.h_min = None + self.h_ideal = None + def add(self, f): # Careful here, if we remove one fracture and then add, we'll have # identical indices. @@ -1812,12 +1816,14 @@ def _determine_mesh_size(self, **kwargs): mesh_size_bound = None return mesh_size, mesh_size_bound elif mode == 'distance': + if self.h_min is None or self.h_ideal is None: + return None, None p = self.decomposition['points'] num_pts = p.shape[1] dist = cg.dist_pointset(p, max_diag=True) mesh_size = np.min(dist, axis=1) - mesh_size = np.maximum(mesh_size, self.hmin * np.ones(num_pts)) + mesh_size = np.maximum(mesh_size, self.h_min * np.ones(num_pts)) mesh_size = np.minimum(mesh_size, self.h_ideal * np.ones(num_pts)) mesh_size_bound = self.h_ideal @@ -1825,10 +1831,10 @@ def _determine_mesh_size(self, **kwargs): else: raise ValueError('Unknown mesh size mode ' + mode) - def compute_distances(self, h_ideal, hmin): + def compute_distances(self, h_ideal=None, h_min=None): self.h_ideal = h_ideal - self.hmin = hmin + self.h_min = h_min isect_pt = np.zeros((3, 0), dtype=np.double) def dist_p(a, b): @@ -2011,6 +2017,8 @@ def to_gmsh(self, file_name, **kwargs): self.zero_d_pt = intersection_points if 'mesh_size' in kwargs.keys(): + # Legacy option, this should be removed. + print('Using old version of mesh size determination') mesh_size, mesh_size_bound = \ determine_mesh_size(self.decomposition['points'], **kwargs['mesh_size']) From 0d444845d8c87e5194d42840f66882ae1f818277 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 13:44:31 +0200 Subject: [PATCH 037/118] Fix of testing when adding point to a Fracture --- src/porepy/fracs/fractures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 2f681c7a37..ef747f3878 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -172,7 +172,7 @@ def add_points(self, p, check_convexity=True, tol=1e-4): if check_convexity: return self.check_convexity() and self.is_planar(tol) else: - return self.check_convexity() + return self.is_planar() def remove_points(self, ind, keep_orig=False): """ Remove points from the fracture definition From 8a0b5e8f5627d3c75bc4263c08cabe925994f44f Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 13:44:40 +0200 Subject: [PATCH 038/118] Skip test of convexity when adding points in fracture extrusion. --- src/porepy/fracs/extrusion.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index ea7e86b5ad..2274de328f 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -171,6 +171,7 @@ def _disc_radius_center(lengths, p0, p1): return radius, np.vstack((mid_point, depth)) +@profile def discs_from_exposure(pt, edges): """ Create fracture discs based on exposed lines in an outrcrop. @@ -213,7 +214,8 @@ def discs_from_exposure(pt, edges): # distribution of the points, but it is the only hard information we # have on the fracture f.add_points(np.vstack((np.hstack((p0[:, i], 0)), - np.hstack((p1[:, i], 0)))).T) + np.hstack((p1[:, i], 0)))).T, + check_convexity=False) fracs.append(Fracture(f.p, check_convexity=False)) return fracs From bd0079bbc42453e1a7041c5edc3b2e4a07b53d4c Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 14:14:21 +0200 Subject: [PATCH 039/118] Remove @profile from fracture extrusion module --- src/porepy/fracs/extrusion.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 2274de328f..0d31b7ff36 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -171,7 +171,6 @@ def _disc_radius_center(lengths, p0, p1): return radius, np.vstack((mid_point, depth)) -@profile def discs_from_exposure(pt, edges): """ Create fracture discs based on exposed lines in an outrcrop. From b039595a27bc7a7da6d84c7b8f71d5c53a674365 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 11 Oct 2017 16:10:05 +0200 Subject: [PATCH 040/118] Update install instructions to avoid newer version of gmsh. --- Install.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Install.md b/Install.md index 7ac08626ca..bbf0a37f7c 100644 --- a/Install.md +++ b/Install.md @@ -18,10 +18,12 @@ The expectation is that it should work out nicely. If you try this, please let u # Setting up GMSH PorePy currently depends on `GMSH` for meshing of fractured domains. To make this work, you need gmsh installed on your system, and PorePy needs to know where to look for it. -First, visit the [GMSH webpage](http://gmsh.info) and download a suitable version. +First, visit the [Gmsh webpage](http://gmsh.info) and download a suitable version. Extract, and move the binary (probably located in the subfolder gmsh-x.x.x-Linux/bin or similar) to whereever you prefer. -Note to Linux users: Although GMSH is available through the standard packaging tools, it tends to be hopelessly outdated, +NOTE: We have experienced that some fracture geometries are best handled by somewhat older versions of Gmsh (2.11 seems to be a good version) - newer versions are often result in error messages. Until this is resolved, we therefore recommend to use Gmsh 2.11 with PorePy. + +Note to Linux users: Although Gmsh is available through the standard packaging tools, it tends to be hopelessly outdated, and resulted in severe issues for the fracture meshing last time we checked. Use the GMSH web page instead. From b6186cb618f82ac7cea7c0de669d0fc8d387c982 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 10:14:48 +0200 Subject: [PATCH 041/118] Restructuring of fracture extrusion module. Expose treatment of individual fractures to user. --- src/porepy/fracs/extrusion.py | 84 ++++++++++++++++++++++++++--------- 1 file changed, 63 insertions(+), 21 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 0d31b7ff36..7a0f438494 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -203,22 +203,71 @@ def discs_from_exposure(pt, edges): fracs = [] for i in range(num_fracs): - # The simplest way of distributing points along the disc seems to be to - # create an elliptic fracture, and pick out the points. - f = EllipticFracture(center=center[:, i], major_axis=radius[i], - minor_axis=radius[i], dip_angle=np.pi/2, - strike_angle=strike_angle[i], major_axis_angle=0) - - # Add the points on the exposed surface. This creates an unequal - # distribution of the points, but it is the only hard information we - # have on the fracture - f.add_points(np.vstack((np.hstack((p0[:, i], 0)), - np.hstack((p1[:, i], 0)))).T, - check_convexity=False) - fracs.append(Fracture(f.p, check_convexity=False)) - + fracs.append(create_fracture(center[:, i], radius[i], np.pi/2, + strike_angle[i], p0[:, i], p1[:, i])) return fracs +def create_fracture(center, radius, dip, strike, exp_0, exp_1): + """ Create a single circular fracture consistent with a given exposure. + + The exposed points will be added to the fracture description. + + Parameters: + center (np.array-like, dim 3): Center of the fracture. + radius (double): Fracture radius. + dip (double): dip angle of the fracture. See EllipticFracture for + details. + strike (np.array-like, dim 3): Strike angle for rotation. See + EllipticFracture for details. + exp_0 (np.array-like, dim 2 or 3): First exposed point. Will be + contained in fracture polygon. + exp_1 (np.array-like, dim 2 or 3): Second exposed point. Will be + contained in fracture polygon. + + Returns: + Fracture: New fracture, according to the specifications. + + """ + + if exp_0.size == 2: + exp_0 = np.append(exp_0, 0) + + if exp_1.size == 2: + exp_1 = np.append(exp_1, 0) + + # The simplest way of distributing points along the disc seems to be to + # create an elliptic fracture, and pick out the points. + f = EllipticFracture(center=center, major_axis=radius, minor_axis=radius, + dip_angle=dip, strike_angle=strike, + major_axis_angle=0) + # Add the points on the exposed surface. This creates an unequal + # distribution of the points, but it is the only hard information we have + # on the fracture + f.add_points(np.vstack((exp_0, exp_1)).T, check_convexity=False) + # Not sure if f still shoudl be EllipticFracture here, or if we should + # create a new fracture with the same point distribution. + return f + + +def rotate_fracture(frac, vec, angle): + """ Rotate a fracture along a specified strike vector. + + Parameters: + frac (Fracture): To be rotated. Points are modified in-place. + vec (np.array-like): Rotation will be around this angle. + ang (double). Rotation angle. Measured in radians. + + """ + + if vec.size == 2: + vec = np.append(vec, 0) + + rot = cg.rot(angle, vec) + p = frac.p + center = np.mean(p, axis=1).reshape((-1, 1)) + frac.p = center + rot.dot(p - center) + frac.points_2_ccw() + def impose_inlcine(fracs, exposure, family, family_mean, family_std): """ Impose incline on the fractures from family-based parameters. @@ -240,13 +289,6 @@ def impose_inlcine(fracs, exposure, family, family_mean, family_std): each family. In radians. """ - def rotate_fracture(frac, vec, angle): - # Helper function to carry out rotation. - rot = cg.rot(angle, vec) - p = frac.p - center = np.mean(p, axis=1).reshape((-1, 1)) - frac.p = center + rot.dot(p - center) - frac.points_2_ccw() exposure = np.vstack((exposure, np.zeros(len(fracs)))) for fi, f in enumerate(fracs): From cd66c988ddaf9df49b0ef98390bb542fb88fd147 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 11:09:26 +0200 Subject: [PATCH 042/118] Allow prescribed angle to center in fracture extrusion. --- src/porepy/fracs/extrusion.py | 37 ++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 7a0f438494..a875338e5b 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -132,20 +132,23 @@ def fracture_length(pt, e): return np.sqrt(np.power(x1-x0, 2) + np.power(y1-y0, 2)) -def _disc_radius_center(lengths, p0, p1): +def _disc_radius_center(lengths, p0, p1, theta=None): """ Compute radius and center of a disc, based on the length of a chord through the disc, and assumptions on the location of the chord. - For the moment, it is assumed that the chord is struck at an arbitrary - hight of the circle. Also, we assume that the chord is horizontal. In the - context of an exposed fracture, this implies that the exposure is - horizontal (I believe), and that an arbitrary portion of the original - (disc-shaped) fracture has been eroded. + The relation between the exposure and the center of the fracture is + given by the theta, which gives the angle between a vertical line through + the disc center and the line through the disc center and any of the + exposure endpoints. If no values are given a random value is assigned, + corresponding to an arbitrary portion of the original (disc-shaped) + fracture has been eroded. Parameters: length (np.array, double, size: num_frac): Of the chords p0 (np.array, 2 x num_frac): One endpoint of fractures. p1 (np.array, 2 x num_frac): Second endpoint of fractures + angle (np.array, num_frac, optional): Angle determining disc center, + see description above. Defaults to random values. Returns: np.array, num_frac: Radius of discs @@ -155,12 +158,10 @@ def _disc_radius_center(lengths, p0, p1): num_frac = lengths.size - # Angle between a vertical line through the disc center and the line - # through the disc center and one of the fracture endpoints. Assumed to be - # uniform, for the lack of more information. # Restrict the angle in the interval (0.1, 0.9) * pi to avoid excessively # large fractures. - theta = np.pi * (0.1 + 0.8 * np.random.rand(num_frac)) + if theta is None: + theta = np.pi * (0.1 + 0.8 * np.random.rand(num_frac)) radius = 0.5 * lengths / np.sin(theta) @@ -171,18 +172,22 @@ def _disc_radius_center(lengths, p0, p1): return radius, np.vstack((mid_point, depth)) -def discs_from_exposure(pt, edges): +def discs_from_exposure(pt, edges, angle=None): """ Create fracture discs based on exposed lines in an outrcrop. - The outcrop is assumed to be in the xy-plane. + The outcrop is assumed to be in the xy-plane. The returned disc will be + vertical, and the points on the outcrop will be included in the polygon + representation. - The returned disc will be vertical, and the points on the outcrop will be - included in the polygon representation. The disc center is calculated using - disc_radius_center(), see that function for assumptions. + The location of the center is calculated from the angle, see + _disc_radius_center() for details. Parameters: pt (np.array, 2 x num_pts): Coordinates of exposed points. edges (np.array, 2 x num_fracs): Connections between fractures. + angle (np.array, num_fracs, optional): See above, and + _disc_radius_center() for description. If not provided, random + values will be drawn. Measured in radians. Returns: list of Fracture: One per fracture trace. @@ -198,7 +203,7 @@ def discs_from_exposure(pt, edges): v = p1 - p0 strike_angle = np.arctan2(v[1], v[0]) - radius, center = _disc_radius_center(lengths, p0, p1) + radius, center = _disc_radius_center(lengths, p0, p1, angle) fracs = [] From 955438e77d1a718191e3b9f798e24638b2566db8 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 11:28:44 +0200 Subject: [PATCH 043/118] Ensure extrusion does not produce unreasonable results. Avoid point contacts, and virtually infinite extent of fractures. --- src/porepy/fracs/extrusion.py | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index a875338e5b..ab93e318c1 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -161,7 +161,13 @@ def _disc_radius_center(lengths, p0, p1, theta=None): # Restrict the angle in the interval (0.1, 0.9) * pi to avoid excessively # large fractures. if theta is None: - theta = np.pi * (0.1 + 0.8 * np.random.rand(num_frac)) + rnd = np.random.rand(num_frac) + # Angles of pi/2 will read to point contacts that cannot be handled + # of the FractureNetwork. Point contacts also make little physical + # sense, so we vaoid them. + hit = np.abs(rnd - 0.5) < 0.01 + rnd[hit] += 0.01 + theta = np.pi * (0.1 + 0.8 * rnd) radius = 0.5 * lengths / np.sin(theta) @@ -186,8 +192,10 @@ def discs_from_exposure(pt, edges, angle=None): pt (np.array, 2 x num_pts): Coordinates of exposed points. edges (np.array, 2 x num_fracs): Connections between fractures. angle (np.array, num_fracs, optional): See above, and - _disc_radius_center() for description. If not provided, random - values will be drawn. Measured in radians. + _disc_radius_center() for description. Values very close to pi/2, 0 + and pi will be modified to avoid unphysical extruded fractures. If + not provided, random values will be drawn. Measured in radians. + Should be between 0 and pi. Returns: list of Fracture: One per fracture trace. @@ -203,6 +211,17 @@ def discs_from_exposure(pt, edges, angle=None): v = p1 - p0 strike_angle = np.arctan2(v[1], v[0]) + if angle is not None: + # Angles of pi/2 will give point contacts + hit = np.abs(angle - np.pi/2) < 0.01 + angle[hit] = angle[hit] + 0.01 + + # Angles of 0 and pi give infinite fractures. + hit = angle < 0.05 + angle[hit] = 0.05 + hit = np.pi - angle < 0.05 + angle[hit] = 0.05 + radius, center = _disc_radius_center(lengths, p0, p1, angle) fracs = [] From 1aca26c880ed36f3c7452a68b8df97718e9bb5f6 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 13:25:22 +0200 Subject: [PATCH 044/118] Minor improvement of EllipticFracture class --- src/porepy/fracs/fractures.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index ef747f3878..12a57226f3 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -856,6 +856,8 @@ def __init__(self, center, major_axis, minor_axis, major_axis_angle, np.pi/6) """ + if center.ndim == 1: + center = center.reshape((-1, 1)) self.center = center # First, populate polygon in the xy-plane @@ -880,7 +882,7 @@ def __init__(self, center, major_axis, minor_axis, major_axis_angle, dip_pts = dip_rot.dot(rot_ref_pts) # Set the points, and store them in a backup. - self.p = center[:, np.newaxis] + dip_pts + self.p = center + dip_pts self.orig_p = self.p.copy() # Compute normal vector From d0591a635a094d7bb7eb15c09bd3c9bc3d89fc26 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 13:26:13 +0200 Subject: [PATCH 045/118] Made function in fracture extrusion module public --- src/porepy/fracs/extrusion.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index ab93e318c1..60c5393aba 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -132,7 +132,7 @@ def fracture_length(pt, e): return np.sqrt(np.power(x1-x0, 2) + np.power(y1-y0, 2)) -def _disc_radius_center(lengths, p0, p1, theta=None): +def disc_radius_center(lengths, p0, p1, theta=None): """ Compute radius and center of a disc, based on the length of a chord through the disc, and assumptions on the location of the chord. @@ -186,7 +186,7 @@ def discs_from_exposure(pt, edges, angle=None): representation. The location of the center is calculated from the angle, see - _disc_radius_center() for details. + disc_radius_center() for details. Parameters: pt (np.array, 2 x num_pts): Coordinates of exposed points. @@ -222,7 +222,7 @@ def discs_from_exposure(pt, edges, angle=None): hit = np.pi - angle < 0.05 angle[hit] = 0.05 - radius, center = _disc_radius_center(lengths, p0, p1, angle) + radius, center = disc_radius_center(lengths, p0, p1, angle) fracs = [] From 92127680e9aa2cbfea3fd2903bc2386f1976ca13 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 13:27:17 +0200 Subject: [PATCH 046/118] Minor changes to fracture extrusion module. --- src/porepy/fracs/extrusion.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 60c5393aba..8ada058cdb 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -172,7 +172,8 @@ def disc_radius_center(lengths, p0, p1, theta=None): radius = 0.5 * lengths / np.sin(theta) # Midpoint in (x, y)-coordinate given as midpoint of exposed line - mid_point = 0.5 * (p0 + p1) + mid_point = 0.5 * (p0 + p1).reshape((2, -1)) + # z-coordinate from angle depth = radius * np.cos(theta) @@ -312,14 +313,20 @@ def impose_inlcine(fracs, exposure, family, family_mean, family_std): family_std (np.array, num_family): Standard deviation of incine for each family. In radians. + Returns: + np.array, size num_frac: Rotation angles. + """ exposure = np.vstack((exposure, np.zeros(len(fracs)))) + all_ang = np.zeros(len(fracs)) for fi, f in enumerate(fracs): fam = family[fi] ang = np.random.normal(loc=family_mean[fam], scale=family_std[fam]) rotate_fracture(f, exposure[:, fi], ang) + all_ang[fam] = ang + return all_ang def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): """ Cut a fracture by a plane, and confine it to one side of the plane. From 1bc8015e63f625e164a70fcd9d57d3984893e7e1 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 13:35:53 +0200 Subject: [PATCH 047/118] Improved cut of extruded fractures. In cases where the cut extends beyond the confining fracture, we calculate the missing distance to the intersection point. --- src/porepy/fracs/extrusion.py | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 8ada058cdb..6f24c38e98 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -343,6 +343,10 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): Returns: Fracture: A copy of the main fracture, cut by the other fracture. + double: In cases where one interseciton point extends beyond the other + fracture, this is the distance between the center of the other + fracture and the intersection point. If both intersections are + within the polygon, None will be returned. Raises: ValueError if the points in the other fracture is too close. This could @@ -431,4 +435,26 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): # plane are present in the final fracture. main_frac.add_points(isect_pt) - return main_frac + # If the main fracture is too large compared to the other, the cut line + # will extend beyond the confining plane. In these cases, compute the + # distance from the fracture center to the outside intersection point. This + # can be used to extend the other fracture so as to avoid such strange + # configurations. + other_center = other_frac.center.reshape((-1, 1)) + other_p = other_frac.p + rot = cg.project_plane_matrix(other_p - other_center) + + other_rot = rot.dot(other_p - other_center)[:2] + isect_rot = rot.dot(isect_pt - other_center)[:2] + + is_inside = cg.is_inside_polygon(other_rot, isect_rot) + # At one point (the exposed point) must be in the polygon of the other + # fracture. + assert is_inside.any() + + if not is_inside.all(): + hit = np.logical_not(is_inside) + r = np.sqrt(np.sum(isect_pt[:, hit]**2)) + return main_frac, r + else: + return main_frac, None From 732b7bfc60d988580731d937003eee3985f1c395 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 14:45:55 +0200 Subject: [PATCH 048/118] Changed keyword arguments in fracture extrusion functions. --- src/porepy/fracs/extrusion.py | 55 ++++++++++++++++++++++------------- 1 file changed, 34 insertions(+), 21 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 6f24c38e98..4b270366a8 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -177,9 +177,9 @@ def disc_radius_center(lengths, p0, p1, theta=None): # z-coordinate from angle depth = radius * np.cos(theta) - return radius, np.vstack((mid_point, depth)) + return radius, np.vstack((mid_point, depth)), theta -def discs_from_exposure(pt, edges, angle=None): +def discs_from_exposure(pt, edges, exposure_angle=None, **kwargs): """ Create fracture discs based on exposed lines in an outrcrop. The outcrop is assumed to be in the xy-plane. The returned disc will be @@ -192,8 +192,8 @@ def discs_from_exposure(pt, edges, angle=None): Parameters: pt (np.array, 2 x num_pts): Coordinates of exposed points. edges (np.array, 2 x num_fracs): Connections between fractures. - angle (np.array, num_fracs, optional): See above, and - _disc_radius_center() for description. Values very close to pi/2, 0 + exposure_angle (np.array, num_fracs, optional): See above, and + disc_radius_center() for description. Values very close to pi/2, 0 and pi will be modified to avoid unphysical extruded fractures. If not provided, random values will be drawn. Measured in radians. Should be between 0 and pi. @@ -212,25 +212,26 @@ def discs_from_exposure(pt, edges, angle=None): v = p1 - p0 strike_angle = np.arctan2(v[1], v[0]) - if angle is not None: + if exposure_angle is not None: # Angles of pi/2 will give point contacts - hit = np.abs(angle - np.pi/2) < 0.01 - angle[hit] = angle[hit] + 0.01 + hit = np.abs(exposure_angle - np.pi/2) < 0.01 + exposure_angle[hit] = exposure_angle[hit] + 0.01 # Angles of 0 and pi give infinite fractures. - hit = angle < 0.05 - angle[hit] = 0.05 - hit = np.pi - angle < 0.05 - angle[hit] = 0.05 + hit = exposure_angle < 0.05 + exposure_angle[hit] = 0.05 + hit = np.pi - exposure_angle < 0.05 + exposure_angle[hit] = 0.05 - radius, center = disc_radius_center(lengths, p0, p1, angle) + radius, center, ang = disc_radius_center(lengths, p0, p1, exposure_angle) fracs = [] for i in range(num_fracs): fracs.append(create_fracture(center[:, i], radius[i], np.pi/2, strike_angle[i], p0[:, i], p1[:, i])) - return fracs + return fracs, ang + def create_fracture(center, radius, dip, strike, exp_0, exp_1): """ Create a single circular fracture consistent with a given exposure. @@ -294,7 +295,8 @@ def rotate_fracture(frac, vec, angle): frac.points_2_ccw() -def impose_inlcine(fracs, exposure, family, family_mean, family_std): +def impose_inlcine(fracs, exposure, frac_family=None, family_mean_incline=None, + family_std_incline=None, **kwargs): """ Impose incline on the fractures from family-based parameters. The incline for each family is specified in terms of its mean and standard @@ -307,22 +309,33 @@ def impose_inlcine(fracs, exposure, family, family_mean, family_std): exposure (np.array, 3xnum_frac): Exposed line for each fracture (visible in outcrop). Rotation will be around this line. family (np.array, num_fracs): For each fracture, which family does it - belong to. - family_mean (np.array, num_family): Mean value of incline for each - family. In radians. - family_std (np.array, num_family): Standard deviation of incine for - each family. In radians. + belong to. If not provided, all fractures are considered to belong + to the same family. + family_mean_incline (np.array, num_family): Mean value of incline for each + family. In radians. Defaults to zero. + family_std_incline (np.array, num_family): Standard deviation of incine for + each family. In radians. Defaultst to zero. + + To set value for each fracture, set family = np.arange(len(fracs)), + family_mean_incline=prescribed_value, and family_std_incline=None. Returns: np.array, size num_frac: Rotation angles. """ + if frac_family is None: + frac_family = np.zeros(len(fracs)) + if family_mean_incline is None: + family_mean_incline = np.zeros(np.unique(frac_family).size) + if family_std_incline is None: + family_std_incline = np.zeros(np.unique(frac_family).size) exposure = np.vstack((exposure, np.zeros(len(fracs)))) all_ang = np.zeros(len(fracs)) for fi, f in enumerate(fracs): - fam = family[fi] - ang = np.random.normal(loc=family_mean[fam], scale=family_std[fam]) + fam = frac_family[fi] + ang = np.random.normal(loc=family_mean_incline[fam], + scale=family_std_incline[fam]) rotate_fracture(f, exposure[:, fi], ang) all_ang[fam] = ang From 57a5181b693e79161d5c2f679b31eb708385dac4 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 12 Oct 2017 14:46:31 +0200 Subject: [PATCH 049/118] Documentation and main function of fracture extrusion module. --- src/porepy/fracs/extrusion.py | 91 +++++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 4b270366a8..5b16b3e391 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -1,9 +1,100 @@ +""" Simple functionality for extrusion of fractures (2D->3D). The main entry +point is the function fractures_from_outcrop(), this wraps the other functions. + +The extrusion is carried out to be consistent with a set of 2D lines, e.g. an +outcrop. Fracture discs are constructed from the exposed lines. In case of +intersections in the exposed plane, the fractures will either cross +(if X-intersection), or the the abutting relation will be preserved (if +Y/T-intersection). If the extruded discs intersect outside the plane of +exposure, this will (almost surely, no checks are actually made) become a +X-intersection. + +The fractures can be assigned a dip from the vertical, again taken as +consistent with the exposed line. + +No attempts are made to create fractures that do not cross the confined plane. + +""" import numpy as np from porepy.utils import comp_geom as cg from porepy.fracs.fractures import EllipticFracture, Fracture +def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, **kwargs): + """ Create a set of fractures compatible with exposed lines in an outcrop. + + See module-level documentation for futher comments. + + Parameters: + pt (np.array, 2 x num_pts): Coordinates of start and endpoints of + extruded lines. + edges (np.array, n x num_pts): Connections between points. Should not + have their crossings removed before. Should have 2 or 3 rows - the + first two identifying start and endpoints, and the third + potentially identifying fracture family relations. + ensure_realistic_cut (boolean, defaults to True): If True, we ensure + that T-intersections do not have cut fractures that extend beyond + the confining fracture. May overrule user-supplied controls on + fracture sizes. + **kwargs: Potentially user defined options. Forwarded to + discs_from_exposure() and impose_inclines() + + Returns: + list of Fracture: Fracture planes. + + """ + # identify crossings + split_pt, split_edges = cg.remove_edge_crossings(pt, edges) + + # Find t-intersections + abutment_pts, prim_frac, sec_frac, other_pt = t_intersections(split_edges) + + # Calculate fracture lengths + lengths = fracture_length(pt, edges) + + # Extrude to fracture discs + fractures, extrude_ang = discs_from_exposure(pt, edges, **kwargs) + + p0 = pt[:, edges[0]] + p1 = pt[:, edges[1]] + exposure = p1 - p0 + + if edges.shape[0] > 2: + family = edges[2] + else: + family = None + + # Impose incline. + rot_ang = impose_inlcine(fractures, exposure, family, **kwargs) + + # Cut fractures + for prim, sec, p in zip(prim_frac, sec_frac, other_pt): + _, radius = cut_fracture_by_plane(fractures[sec], fractures[prim], + split_pt[:, p]) + # If specified, ensure that cuts in T-intersections appear realistic. + if ensure_realistic_cuts and radius is not None: + ang = np.arctan2(0.5*lengths[prim], radius) + + # Ensure the center of both fractures are on the same side of the + # exposed plane - if not, the cut will still be bad. + if extrude_ang[sec] > np.pi/2 and ang < np.pi/2: + ang = np.pi-ang + elif extrude_ang[sec] < np.pi/2 and ang > np.pi/2: + ang = np.pi-ang + + e0 = p0[:, prim] + e1 = p1[:, prim] + new_radius, center, _ = disc_radius_center(lengths[prim], e0, e1, + theta=ang) + strike = np.arctan2(e1[1] - e0[1], e1[0]- e0[0]) + f = create_fracture(center, new_radius, np.pi/2, strike, e0, e1) + rotate_fracture(f, e1-e0, rot_ang[prim]) + fractures[prim] = f + + return fractures + + def _intersection_by_num_node(edges, num): """ Find all edges involved in intersections with a certain number of intersecting lines. From 6a48324c62f4266e9555e91a5701f357ca085262 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 16 Oct 2017 10:55:59 +0200 Subject: [PATCH 050/118] Make bounding_box function in fracture importer public --- src/porepy/fracs/importer.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/importer.py b/src/porepy/fracs/importer.py index 437b476af7..c385d744e4 100644 --- a/src/porepy/fracs/importer.py +++ b/src/porepy/fracs/importer.py @@ -35,7 +35,7 @@ def from_csv(f_name, mesh_kwargs, domain=None, pause=False,\ # Define the domain as bounding-box if not defined if domain is None: overlap = kwargs.get('domain_overlap', 0) - domain = _bounding_box(pts, overlap) + domain = bounding_box(pts, overlap) if return_domain: return meshing.simplex_grid(f_set, domain, **mesh_kwargs), domain @@ -105,7 +105,7 @@ def fractures_from_csv(f_name, tagcols=None, **kwargs): #------------------------------------------------------------------ -def _bounding_box(pts, overlap=0): +def bounding_box(pts, overlap=0): """ Obtain a bounding box for a point cloud. Parameters: From a7d6d6c9673ef2fe300241486ae3d87fe3318061 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 16 Oct 2017 17:40:20 +0200 Subject: [PATCH 051/118] Bugfix in fracture extrusion method. --- src/porepy/fracs/extrusion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 5b16b3e391..53d3dc5774 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -154,7 +154,7 @@ def t_intersections(edges): for i, (pi, ei) in enumerate(zip(abutments, edges_of_abutments)): # Count number of occurences for each fracture associated with this # intersection. - fi_all = frac_num[edges_of_abutments] + fi_all = frac_num[ei] fi, count = np.unique(fi_all, return_counts=True) assert fi.size == 2 # Find the fracture number associated with main and abutting edge. @@ -165,7 +165,7 @@ def t_intersections(edges): primal_frac[i] = fi[0] sec_frac[i] = fi[1] # Also find the other point of the abutting edge - ind = np.where(fi_all == sec_frac[i])[1] + ind = np.where(fi_all == sec_frac[i]) ei_abut = ei[ind] assert ei_abut.size == 1 if edges[0, ei_abut] == pi: From 334b79edf10c952a297fef39ab554e8e05b152c4 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 17 Oct 2017 08:38:24 +0200 Subject: [PATCH 052/118] Bugfix in FractureNetwork computation of mesh size. Bad variable name. --- src/porepy/fracs/fractures.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 12a57226f3..d21b1d121a 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1877,7 +1877,7 @@ def dist_p(a, b): d_2 = dist_p(cp[pi, si], f.p[:, (si+1)%nfp]) # If the intersection point is not very close to any of # the points on the segment, we split the segment. - if d_1 > hmin and d_2 > hmin: + if d_1 > h_min and d_2 > h_min: np.insert(f.p, (si+1)%nfp, cp[pi, si], axis=1) # Take note of the intersecting fractures @@ -1912,7 +1912,7 @@ def dist_p(a, b): if d[mi] < h_ideal: d_1 = dist_p(cp_f[:, mi], f.p[:, si]) d_2 = dist_p(cp_f[:, mi], f.p[:, (si+1)%nfp]) - if d_1 > hmin and d_2 > hmin: + if d_1 > h_min and d_2 > h_min: np.insert(f.p, (si+1)%nfp, cp[:, pi], axis=1) # Finally, cover the case where the smallest distance is given From 2e8e8ab891a8f1d73bfee2bc1451afc4d92c7fe7 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 17 Oct 2017 08:42:51 +0200 Subject: [PATCH 053/118] Improved simplex meshing of 3d fracture networks. Deal with arguments being both np-arrays and Fractures. The previous version failed to deal with EllipticFractures. The current verison is not perfect either - better handling of isinstance is needed. --- src/porepy/fracs/simplex.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/porepy/fracs/simplex.py b/src/porepy/fracs/simplex.py index 62c1bef73a..9ce8e41e56 100644 --- a/src/porepy/fracs/simplex.py +++ b/src/porepy/fracs/simplex.py @@ -73,12 +73,19 @@ def tetrahedral_grid(fracs=None, box=None, network=None, **kwargs): frac_list = [] for f in fracs: - if isinstance(f, fractures.Fracture) \ - or isinstance(f, porepy.Fracture): + # Input can be either numpy arrays or predifined fractures. As a + # guide, we treat f as a fracture if it has an attribute p which is + # a numpy array. + # If f turns out not to be a fracture, strange errors will result + # as the further program tries to access non-existing methods. + # The correct treatment here would be several + # isinstance-statements, but that became less than elegant. To + # revisit. + if hasattr(f, 'p') and isinstance(f.p, np.ndarray): + # Convert the fractures from numpy representation to our 3D + # fracture data structure. frac_list.append(f) else: - # Convert the fractures from numpy representation to our 3D - # fracture data structure.. frac_list.append(fractures.Fracture(f)) # Combine the fractures into a network From eced078c6cd1d5699afe3d19c984f644878e91bc Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 17 Oct 2017 08:44:59 +0200 Subject: [PATCH 054/118] Allow meshing() of fractures to deal with FractureNetworks. Bit of an ugly hack, but it works. --- src/porepy/fracs/meshing.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/meshing.py b/src/porepy/fracs/meshing.py index b66833182e..5ebc9561fd 100644 --- a/src/porepy/fracs/meshing.py +++ b/src/porepy/fracs/meshing.py @@ -15,7 +15,7 @@ from porepy.utils import setmembership, mcolon -def simplex_grid(fracs, domain, **kwargs): +def simplex_grid(fracs=None, domain=None, network=None, **kwargs): """ Main function for grid generation. Creates a fractured simiplex grid in 2 or 3 dimensions. @@ -52,7 +52,9 @@ def simplex_grid(fracs, domain, **kwargs): gb = simplex_grid(fracs, domain) """ - if 'zmax' in domain: + if domain is None: + ndim = 3 + elif 'zmax' in domain: ndim = 3 elif 'ymax' in domain: ndim = 2 @@ -61,6 +63,8 @@ def simplex_grid(fracs, domain, **kwargs): # Call relevant method, depending on grid dimensions. if ndim == 2: + assert fracs is not None, '2d requires definition of fractures' + assert domain is not None, '2d requires definition of domain' # Convert the fracture to a fracture dictionary. if len(fracs) == 0: f_lines = np.zeros((2, 0)) @@ -71,7 +75,7 @@ def simplex_grid(fracs, domain, **kwargs): frac_dic = {'points': f_pts, 'edges': f_lines} grids = simplex.triangle_grid(frac_dic, domain, **kwargs) elif ndim == 3: - grids = simplex.tetrahedral_grid(fracs, domain, **kwargs) + grids = simplex.tetrahedral_grid(fracs, domain, network, **kwargs) else: raise ValueError('Only support for 2 and 3 dimensions') # Tag tip faces From 56c8474e5506a787701bf3fc177e8e90c532b256 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 17 Oct 2017 08:53:21 +0200 Subject: [PATCH 055/118] Always impose domain for fracture network before meshing. If no domain size is not specified, a box larger than the network is constructed. --- src/porepy/fracs/simplex.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/simplex.py b/src/porepy/fracs/simplex.py index 9ce8e41e56..b9a831a77c 100644 --- a/src/porepy/fracs/simplex.py +++ b/src/porepy/fracs/simplex.py @@ -92,9 +92,9 @@ def tetrahedral_grid(fracs=None, box=None, network=None, **kwargs): network = fractures.FractureNetwork(frac_list, verbose=verbose, tol=kwargs.get('tol', 1e-4)) - # Impose domain boundary. - if box is not None: - network.impose_external_boundary(box) + # Impose external boundary. If box is None, a domain size somewhat larger + # than the network will be assigned. + network.impose_external_boundary(box) # Find intersections and split them, preparing the way for dumping the # network to gmsh From b27fdb0413a9ae21f7db7955cda347ed0b96ab67 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 17 Oct 2017 14:53:30 +0200 Subject: [PATCH 056/118] Loosen tolerance in polygon segment intersection, comp geom. Strict tolerances gave problems for T-type intersections of fractures. --- src/porepy/utils/comp_geom.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 8c574b3bbc..a4e275c8d4 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1342,7 +1342,7 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): pt_2 = poly_2_rot[:, ind[i+1]] # Check if segment crosses z=0 in the rotated coordinates - if max(pt_1[2], pt_2[2]) < 0 or min(pt_1[2], pt_2[2]) > 0: + if max(pt_1[2], pt_2[2]) < -tol or min(pt_1[2], pt_2[2]) > tol: continue dx = pt_2[0] - pt_1[0] @@ -1357,7 +1357,7 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): # Sanity check. We have ruled out segments not crossing the # origin above. - assert t >= 0 and t <= 1 + assert t >= -tol and t <= 1+tol # x and y-coordinate for z=0 x0 = pt_1[0] + dx * t From b20d6042ed0be135cb89e79dbcfb29bbd26aa5c9 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 17 Oct 2017 14:54:39 +0200 Subject: [PATCH 057/118] Minor bugfix in FractureNetwork distance computation. Cut and paste errors --- src/porepy/fracs/fractures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index d21b1d121a..3d06991c3d 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1913,7 +1913,7 @@ def dist_p(a, b): d_1 = dist_p(cp_f[:, mi], f.p[:, si]) d_2 = dist_p(cp_f[:, mi], f.p[:, (si+1)%nfp]) if d_1 > h_min and d_2 > h_min: - np.insert(f.p, (si+1)%nfp, cp[:, pi], axis=1) + np.insert(f.p, (si+1)%nfp, cp_f[:, si], axis=1) # Finally, cover the case where the smallest distance is given # by a point. Points with false in_point should already be From 923487624b0cc1aa866101f2d3d39bfd5984c10e Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 17 Oct 2017 15:15:34 +0200 Subject: [PATCH 058/118] Improved test of point inside polygon. Removed tacit assumption that the points are sorted in a ccw order. --- src/porepy/utils/comp_geom.py | 8 ++++++++ test/unit/test_inside_polygon.py | 19 +++++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index a4e275c8d4..d207be7c29 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -811,6 +811,14 @@ def is_inside_polygon(poly, p, tol=0, default=False): else: pt = p + # The test uses is_ccw_polyline, and tacitly assumes that the polygon + # vertexes is sorted in a ccw fashion. If this is not the case, flip the + # order of the nodes on a copy, and use this for the testing. + # Note that if the nodes are not cw nor ccw (e.g. they are crossing), the + # test cannot be trusted anyhow. + if not is_ccw_polygon(poly): + poly = poly.copy()[:, ::-1] + poly_size = poly.shape[1] inside = np.ones(pt.shape[1], dtype=np.bool) diff --git a/test/unit/test_inside_polygon.py b/test/unit/test_inside_polygon.py index 2b33a701a0..8334a64074 100644 --- a/test/unit/test_inside_polygon.py +++ b/test/unit/test_inside_polygon.py @@ -52,6 +52,25 @@ def test_multiple_points(self): assert inside[0] assert not inside[1] + def test_large_polygon(self): + a = np.array([[ -2.04462568e-01, -1.88898782e-01, -1.65916617e-01, + -1.44576869e-01, -7.82444375e-02, -1.44018389e-16, + 7.82444375e-02, 1.03961083e-01, 1.44576869e-01, + 1.88898782e-01, 2.04462568e-01, 1.88898782e-01, + 1.44576869e-01, 7.82444375e-02, 1.44018389e-16, + -7.82444375e-02, -1.44576869e-01, -1.88898782e-01], + [ -1.10953147e-16, 7.82444375e-02, 1.19484803e-01, + 1.44576869e-01, 1.88898782e-01, 2.04462568e-01, + 1.88898782e-01, 1.76059749e-01, 1.44576869e-01, + 7.82444375e-02, 1.31179355e-16, -7.82444375e-02, + -1.44576869e-01, -1.88898782e-01, -2.04462568e-01, + -1.88898782e-01, -1.44576869e-01, -7.82444375e-02]]) + b = np.array([[ 0.1281648 , 0.04746067], + [-0.22076491, 0.16421546]]) + inside = cg.is_inside_polygon(a, b) + assert not inside[0] + assert inside[1] + if __name__ == '__main__': unittest.main() From bfbd72c5717c5ec8fd2a9c499bfb21e261d4e579 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 19 Oct 2017 13:14:12 +0200 Subject: [PATCH 059/118] unit tests for fracture extrusion. --- test/integration/test_fracture_extrusion.py | 39 +++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 test/integration/test_fracture_extrusion.py diff --git a/test/integration/test_fracture_extrusion.py b/test/integration/test_fracture_extrusion.py new file mode 100644 index 0000000000..19b3e2f441 --- /dev/null +++ b/test/integration/test_fracture_extrusion.py @@ -0,0 +1,39 @@ +import unittest +import numpy as np + +from porepy.fracs import extrusion + +class TestFractureExtrusion(unittest.TestCase): + + def test_cut_off(self): + # Two fractures, the second abutting. Check that the second fracture is + # cut correctly. + p = np.array([[0, 0], [2, 0], [1, 0], [1, 10]]).T + e = np.array([[0, 1], [2, 3]]).T + + fractures = extrusion.fractures_from_outcrop(p, e) + assert np.all(fractures[1].p[1] >= 0) + + def test_not_aligned_outcrop(self): + # The fractures are not aligned with the mesh + p = np.array([[0, 0], [2, 2], [1, 1], [1, 3]]).T + e = np.array([[0, 1], [2, 3]]).T + + tol = 1e-4 + + fractures = extrusion.fractures_from_outcrop(p, e, tol=tol) + # + assert np.all(fractures[1].p[1] >= fractures[1].p[0]-tol) + + + def test_no_family_given(self): + # Assign fracture family automatically + tol = 1e-2 + p = np.array([[0, 0], [2, 0], [1, 0], [1, 3]]).T + e = np.array([[0, 1], [2, 3]]).T + + fractures = extrusion.fractures_from_outcrop(p, e, tol=tol) + assert np.all(fractures[1].p[1] >= -tol) + + if __name__ == '__main__': + unittest.main() From 87d4c3c07f139665256aa26b1459090aa57330f0 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 19 Oct 2017 13:15:00 +0200 Subject: [PATCH 060/118] Fix of tolerance in comp geom test of ccw polyline. Try to allow for rounding errors in the computation. --- src/porepy/utils/comp_geom.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index d207be7c29..7be324714b 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -781,8 +781,7 @@ def is_ccw_polyline(p1, p2, p3, tol=0, default=False): if np.abs(cross_product) <= tol: return default - - return cross_product > 0 + return cross_product > -tol #----------------------------------------------------------------------------- From a35d5e8bc22b4d4c803c3ad9d1e40b5a6758f568 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 19 Oct 2017 13:16:48 +0200 Subject: [PATCH 061/118] Cleanup of fracture extrusion edge num vs fracture family. Concepts was confused before, better now. --- src/porepy/fracs/extrusion.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 53d3dc5774..664284eef0 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -21,7 +21,7 @@ from porepy.fracs.fractures import EllipticFracture, Fracture -def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, **kwargs): +def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, **kwargs): """ Create a set of fractures compatible with exposed lines in an outcrop. See module-level documentation for futher comments. @@ -29,10 +29,8 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, **kwargs): Parameters: pt (np.array, 2 x num_pts): Coordinates of start and endpoints of extruded lines. - edges (np.array, n x num_pts): Connections between points. Should not - have their crossings removed before. Should have 2 or 3 rows - the - first two identifying start and endpoints, and the third - potentially identifying fracture family relations. + edges (np.array, 2 x num_pts): Connections between points. Should not + have their crossings removed before. ensure_realistic_cut (boolean, defaults to True): If True, we ensure that T-intersections do not have cut fractures that extend beyond the confining fracture. May overrule user-supplied controls on @@ -44,6 +42,9 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, **kwargs): list of Fracture: Fracture planes. """ + assert edges.shape[0] == 2, 'Edges have two endpoints' + edges = np.vstack((edges, np.arange(edges.shape[1], dtype=np.int))) + # identify crossings split_pt, split_edges = cg.remove_edge_crossings(pt, edges) @@ -60,13 +61,8 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, **kwargs): p1 = pt[:, edges[1]] exposure = p1 - p0 - if edges.shape[0] > 2: - family = edges[2] - else: - family = None - # Impose incline. - rot_ang = impose_inlcine(fractures, exposure, family, **kwargs) + rot_ang = impose_inlcine(fractures, exposure, frac_family=family, **kwargs) # Cut fractures for prim, sec, p in zip(prim_frac, sec_frac, other_pt): From 4988d89d8433979d4ae3ea77e43af9ec64db3dc9 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 19 Oct 2017 13:17:56 +0200 Subject: [PATCH 062/118] Improved adding of outcrop points to fracture extrusion. If only exposed points are added, this often leads to point-connections. Add extra points to avoid this. --- src/porepy/fracs/extrusion.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 664284eef0..91dcf2bf75 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -84,7 +84,8 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, * new_radius, center, _ = disc_radius_center(lengths[prim], e0, e1, theta=ang) strike = np.arctan2(e1[1] - e0[1], e1[0]- e0[0]) - f = create_fracture(center, new_radius, np.pi/2, strike, e0, e1) + f = create_fracture(center, new_radius, np.pi/2, strike, + np.vstack((e0, e1)).T) rotate_fracture(f, e1-e0, rot_ang[prim]) fractures[prim] = f @@ -315,12 +316,17 @@ def discs_from_exposure(pt, edges, exposure_angle=None, **kwargs): fracs = [] for i in range(num_fracs): + z = 2 * center[2, i] + extra_point_depth = np.array([0, 0, z, z]) + extra_points = np.vstack((np.vstack((p0[:, i], p1[:, i], p0[:, i], + p1[:, i])).T, + extra_point_depth)) fracs.append(create_fracture(center[:, i], radius[i], np.pi/2, - strike_angle[i], p0[:, i], p1[:, i])) + strike_angle[i], extra_points)) return fracs, ang -def create_fracture(center, radius, dip, strike, exp_0, exp_1): +def create_fracture(center, radius, dip, strike, extra_points): """ Create a single circular fracture consistent with a given exposure. The exposed points will be added to the fracture description. @@ -332,21 +338,16 @@ def create_fracture(center, radius, dip, strike, exp_0, exp_1): details. strike (np.array-like, dim 3): Strike angle for rotation. See EllipticFracture for details. - exp_0 (np.array-like, dim 2 or 3): First exposed point. Will be - contained in fracture polygon. - exp_1 (np.array-like, dim 2 or 3): Second exposed point. Will be - contained in fracture polygon. + extra_points (np.array, 3xnpt): Extra points to be added to the + fracture. The points are assumed to lie on the ellipsis. Returns: Fracture: New fracture, according to the specifications. """ - - if exp_0.size == 2: - exp_0 = np.append(exp_0, 0) - - if exp_1.size == 2: - exp_1 = np.append(exp_1, 0) + if extra_points.shape[0] == 2: + extra_points = np.vstack((extra_points, + np.zeros(extra_points.shape[1]))) # The simplest way of distributing points along the disc seems to be to # create an elliptic fracture, and pick out the points. @@ -356,7 +357,7 @@ def create_fracture(center, radius, dip, strike, exp_0, exp_1): # Add the points on the exposed surface. This creates an unequal # distribution of the points, but it is the only hard information we have # on the fracture - f.add_points(np.vstack((exp_0, exp_1)).T, check_convexity=False) + f.add_points(extra_points, check_convexity=False) # Not sure if f still shoudl be EllipticFracture here, or if we should # create a new fracture with the same point distribution. return f From cdbc44dcb7b08a1658ec8bf7140213c5cb999668 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 19 Oct 2017 13:19:54 +0200 Subject: [PATCH 063/118] Minor bugfixes in fracture extrusion. --- src/porepy/fracs/extrusion.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 91dcf2bf75..5072ecb393 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -67,7 +67,7 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, * # Cut fractures for prim, sec, p in zip(prim_frac, sec_frac, other_pt): _, radius = cut_fracture_by_plane(fractures[sec], fractures[prim], - split_pt[:, p]) + split_pt[:, p], **kwargs) # If specified, ensure that cuts in T-intersections appear realistic. if ensure_realistic_cuts and radius is not None: ang = np.arctan2(0.5*lengths[prim], radius) @@ -412,7 +412,7 @@ def impose_inlcine(fracs, exposure, frac_family=None, family_mean_incline=None, """ if frac_family is None: - frac_family = np.zeros(len(fracs)) + frac_family = np.zeros(len(fracs), dtype=np.int) if family_mean_incline is None: family_mean_incline = np.zeros(np.unique(frac_family).size) if family_std_incline is None: @@ -548,7 +548,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): other_rot = rot.dot(other_p - other_center)[:2] isect_rot = rot.dot(isect_pt - other_center)[:2] - is_inside = cg.is_inside_polygon(other_rot, isect_rot) + is_inside = cg.is_inside_polygon(other_rot, isect_rot, tol, default=True) # At one point (the exposed point) must be in the polygon of the other # fracture. assert is_inside.any() From b55f9679682b0a72db54a18fe20b56893ce07f9f Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 19 Oct 2017 13:21:06 +0200 Subject: [PATCH 064/118] Fix of tolerance in unit test for fracture extrusion. --- test/integration/test_fracture_extrusion.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/integration/test_fracture_extrusion.py b/test/integration/test_fracture_extrusion.py index 19b3e2f441..ca4f5c8798 100644 --- a/test/integration/test_fracture_extrusion.py +++ b/test/integration/test_fracture_extrusion.py @@ -8,11 +8,12 @@ class TestFractureExtrusion(unittest.TestCase): def test_cut_off(self): # Two fractures, the second abutting. Check that the second fracture is # cut correctly. + tol = 1e-4 p = np.array([[0, 0], [2, 0], [1, 0], [1, 10]]).T e = np.array([[0, 1], [2, 3]]).T - fractures = extrusion.fractures_from_outcrop(p, e) - assert np.all(fractures[1].p[1] >= 0) + fractures = extrusion.fractures_from_outcrop(p, e, tol=tol) + assert np.all(fractures[1].p[1] > -tol) def test_not_aligned_outcrop(self): # The fractures are not aligned with the mesh From 66198d37235a761b9462d60e09d02a4551d513f7 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 20 Oct 2017 19:48:56 +0200 Subject: [PATCH 065/118] Input formating in comp geom file. Ensure the input in intersection finding between 3d segments always are 1d. --- src/porepy/utils/comp_geom.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 7be324714b..db41373466 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -997,10 +997,10 @@ def segments_intersect_3d(start_1, end_1, start_2, end_2, tol=1e-8): """ # Convert input to numpy if necessary - start_1 = np.asarray(start_1).astype(np.float) - end_1 = np.asarray(end_1).astype(np.float) - start_2 = np.asarray(start_2).astype(np.float) - end_2 = np.asarray(end_2).astype(np.float) + start_1 = np.asarray(start_1).astype(np.float).ravel() + end_1 = np.asarray(end_1).astype(np.float).ravel() + start_2 = np.asarray(start_2).astype(np.float).ravel() + end_2 = np.asarray(end_2).astype(np.float).ravel() # Short hand for component of start and end points, as well as vectors # along lines. From 8892a7415acfacc8ad1ddb5932de25b0b54d6544 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 20 Oct 2017 19:49:58 +0200 Subject: [PATCH 066/118] Tolerance tuning in comp geom polygon intersection. Try to be somewhat more forgiving on tolerances in intersection between polygon segments. Turned out to be important for fractures from extrusions, where T-intersection are found. --- src/porepy/utils/comp_geom.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index db41373466..6d75ee2fcd 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1303,9 +1303,9 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): # If the rotation of whole second point cloud lies on the same side of z=0, # there are no intersections. - if poly_2_rot[2].min() > 0: + if poly_2_rot[2].min() > tol: return None - elif poly_2_rot[2].max() < 0: + elif poly_2_rot[2].max() < -tol: return None # Check if the second plane is parallel to the first (same xy-plane) @@ -1364,7 +1364,8 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): # Sanity check. We have ruled out segments not crossing the # origin above. - assert t >= -tol and t <= 1+tol + if t < -tol or t > 1+tol: + continue # x and y-coordinate for z=0 x0 = pt_1[0] + dx * t From cd470e2ef27fbd2632082b6432f3b029d943de7b Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 20 Oct 2017 19:51:35 +0200 Subject: [PATCH 067/118] Improved polygon segment intersection when segment in polygon plane. Also unit tests. --- src/porepy/utils/comp_geom.py | 35 ++++++++++++++++- .../unit/test_polygon_segment_intersection.py | 39 +++++++++++++++++-- 2 files changed, 70 insertions(+), 4 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 6d75ee2fcd..82f3f0a671 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1381,9 +1381,42 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): # inverse rotation, 3) translate to original coordinate. isect = np.hstack((isect, irot.dot(_to3D(p_00)) + center_1)) + elif np.abs(pt_1[2]) < tol and np.abs(pt_2[2]) < tol: + # The segment lies completely within the polygon plane. + both_pts = np.vstack((pt_1, pt_2)).T + # Find points within tho polygon itself + inside = is_inside_polygon(poly_1_xy, both_pts[:2],tol=tol) + + if inside.all(): + # Both points are inside, add and go on + isect = np.hstack((isect, irot.dot(both_pts) + center_1)) + else: + # A single point is inside. Need to find the intersection between this line segment and the polygon + if inside.any(): + isect_loc = both_pts[:2, inside].reshape((2, -1)) + p1 = both_pts[:, inside] + p2 = both_pts[:, np.logical_not(inside)] + else: + isect_loc = np.empty((2, 0)) + p1 = both_pts[:, 0] + p2 = both_pts[:, 1] + poly_1_start = poly_1 + poly_1_end = np.roll(poly_1, 1, axis=1) + for j in range(poly_1.shape[1]): + ip = segments_intersect_3d(p1, p2, poly_1_start[:, j], + poly_1_end[:, j]) + if ip is not None: + isect_loc = np.hstack((isect_loc, ip[:2])) + + isect = np.hstack((isect, irot.dot(_to3D(isect_loc)) + center_1)) + if isect.shape[1] == 0: isect = None + # For points lying in the plane of poly_1, the same points may be found + # several times + if isect is not None: + isect, _, _ = setmembership.unique_columns_tol(isect, tol=tol) return isect @@ -2141,7 +2174,7 @@ def dist_segments_polygon(start, end, poly, tol=1e-5): start = orig_start end = orig_end - # Distance from endpoints to + # Distance from endpoints to d_start_poly, cp_s_p, s_in_poly = dist_points_polygon(start, poly) d_end_poly, cp_e_p, e_in_poly = dist_points_polygon(end, poly) diff --git a/test/unit/test_polygon_segment_intersection.py b/test/unit/test_polygon_segment_intersection.py index bd80a6cb37..6780917cd0 100644 --- a/test/unit/test_polygon_segment_intersection.py +++ b/test/unit/test_polygon_segment_intersection.py @@ -38,7 +38,7 @@ def test_mutual_intersection(self): p_i_known_2 = np.array([0, 0, 1.0]).reshape((-1, 1)) assert np.min(np.sum(np.abs(p_1_3 - p_i_known_1), axis=0)) < 1e-5 - # Then intersection of plane of 3 by edges of 1. + # Then intersection of plane of 3 by edges of 1. p_3_1 = cg.polygon_segment_intersect(p3, p1) p_i_known_2 = np.array([0, 0, 1.0]).reshape((-1, 1)) @@ -56,7 +56,7 @@ def test_mutual_intersection_not_at_origin(self): p_i_known_1 = np.array([0, 0, 0.5]).reshape((-1, 1)) + incr assert np.min(np.sum(np.abs(p_1_3 - p_i_known_1), axis=0)) < 1e-5 - # Then intersection of plane of 3 by edges of 1. + # Then intersection of plane of 3 by edges of 1. p_3_1 = cg.polygon_segment_intersect(p3, p1) p_i_known_2 = np.array([0, 0, 1.0]).reshape((-1, 1)) + incr @@ -95,7 +95,7 @@ def test_issue_16(self): frac2 = np.array([[2, 2, 2], [2, 4, 1], [1, 2, 4]]) - + isect_known = np.array([[2], [5/3], [2]]) isect = cg.polygon_segment_intersect(frac1, frac2) assert np.allclose(isect, isect_known) @@ -103,6 +103,39 @@ def test_issue_16(self): isect = cg.polygon_segment_intersect(frac1[:, [0, 2, 1]], frac2) assert np.allclose(isect, isect_known) + def test_segment_in_polygon_plane(self): + # One segment lies fully within a plane + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + p2 = np.array([[0.3, 0.5, 0], + [0.7, 0.5, 0], + [0.7, 0.5, 1], + [0.3, 0.5, 1]]).T + isect_known = np.array([[0.3, 0.5, 0], + [0.7, 0.5, 0]]).T + isect = cg.polygon_segment_intersect(p1, p2) + assert np.allclose(isect, isect_known) + + def test_segment_in_plane_but_not_in_polygon(self): + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + p2 = np.array([[1.3, 0.5, 0], + [1.7, 0.5, 0], + [1.7, 0.5, 1], + [1.3, 0.5, 1]]).T + isect = cg.polygon_segment_intersect(p1, p2) + assert isect is None + + def test_segment_partly_in_polygon(self): + # Segment in plane, one point inside, one point outside polygon + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + p2 = np.array([[0.3, 0.5, 0], + [1.7, 0.5, 0], + [0.7, 0.5, 1], + [0.3, 0.5, 1]]).T + isect_known = np.array([[0.3, 0.5, 0], + [1, 0.5, 0]]).T + isect = cg.polygon_segment_intersect(p1, p2) + assert np.allclose(isect, isect_known) + if __name__ == '__main__': unittest.main() From 44a44d1fc88690093ce969e1b4e33cc5e229a46f Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Sun, 22 Oct 2017 18:55:56 +0200 Subject: [PATCH 068/118] Minor updates to fracture extrusion model. --- src/porepy/fracs/extrusion.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 5072ecb393..8bededb005 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -46,7 +46,7 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, * edges = np.vstack((edges, np.arange(edges.shape[1], dtype=np.int))) # identify crossings - split_pt, split_edges = cg.remove_edge_crossings(pt, edges) + split_pt, split_edges = cg.remove_edge_crossings(pt, edges, **kwargs) # Find t-intersections abutment_pts, prim_frac, sec_frac, other_pt = t_intersections(split_edges) @@ -253,8 +253,8 @@ def disc_radius_center(lengths, p0, p1, theta=None): # Angles of pi/2 will read to point contacts that cannot be handled # of the FractureNetwork. Point contacts also make little physical # sense, so we vaoid them. - hit = np.abs(rnd - 0.5) < 0.01 - rnd[hit] += 0.01 + hit = np.abs(rnd - 0.5) < 0.05 + rnd[hit] += 0.05 theta = np.pi * (0.1 + 0.8 * rnd) radius = 0.5 * lengths / np.sin(theta) @@ -321,6 +321,7 @@ def discs_from_exposure(pt, edges, exposure_angle=None, **kwargs): extra_points = np.vstack((np.vstack((p0[:, i], p1[:, i], p0[:, i], p1[:, i])).T, extra_point_depth)) + fracs.append(create_fracture(center[:, i], radius[i], np.pi/2, strike_angle[i], extra_points)) return fracs, ang @@ -429,11 +430,13 @@ def impose_inlcine(fracs, exposure, frac_family=None, family_mean_incline=None, return all_ang -def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): +def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, + recompute_center=False): """ Cut a fracture by a plane, and confine it to one side of the plane. Intended use is to confine abutting fractures (T-intersections) to one side - of the fracture it hits. + of the fracture it hits. This is done by deleting points on the abutting + fracture. Parameters: main_frac (Fracture): The fracture to be cut. @@ -529,6 +532,9 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4): eliminate = np.where(sgn * right_sign < 0)[0] main_frac.remove_points(eliminate) + if recompute_center: + main_frac.compute_center() + # Add intersection points on the main fracture. One of these may already be # present, as the point of extrusion, but add_point will uniquify the point # cloud. From 58d95b694719844a8c2bc3d0c3c09eade8e4a4b4 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Sun, 22 Oct 2017 19:34:50 +0200 Subject: [PATCH 069/118] Opt-out of case in polygon segment intersection, comp geom. In some cases, e.g. fracture intersection, it is beneficial to disregard polygon segment intersection finds points on the boundary. Turned out to be important related to fracture extrusion. --- src/porepy/utils/comp_geom.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 82f3f0a671..5459c6ec8a 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1249,16 +1249,13 @@ def polygon_boundaries_intersect(poly_1, poly_2, tol=1e-8): #---------------------------------------------------------- -def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): +def polygon_segment_intersect(poly_1, poly_2, tol=1e-8, include_bound_pt=True): """ Find intersections between polygons embeded in 3D. The intersections are defined as between the interior of the first polygon and the boundary of the second. - TODO: - 1) Also cover case where the one polygon ends in the plane of the other. - Parameters: poly_1 (np.ndarray, 3xn1): Vertexes of polygon, assumed ordered as cw or ccw. @@ -1266,6 +1263,9 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): as cw or ccw. tol (double, optional): Tolerance for when two points are equal. Defaults to 1e-8. + include_bound_pt (boolean, optional): Include cases where a segment is + in the plane of the first ploygon, and the segment crosses the + polygon boundary. Defaults to True. Returns: np.ndarray, size 3 x num_isect, coordinates of intersection points; or @@ -1400,13 +1400,14 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8): isect_loc = np.empty((2, 0)) p1 = both_pts[:, 0] p2 = both_pts[:, 1] - poly_1_start = poly_1 - poly_1_end = np.roll(poly_1, 1, axis=1) - for j in range(poly_1.shape[1]): - ip = segments_intersect_3d(p1, p2, poly_1_start[:, j], - poly_1_end[:, j]) - if ip is not None: - isect_loc = np.hstack((isect_loc, ip[:2])) + if include_bound_pt: + poly_1_start = poly_1 + poly_1_end = np.roll(poly_1, 1, axis=1) + for j in range(poly_1.shape[1]): + ip = segments_intersect_3d(p1, p2, poly_1_start[:, j], + poly_1_end[:, j]) + if ip is not None: + isect_loc = np.hstack((isect_loc, ip[:2])) isect = np.hstack((isect, irot.dot(_to3D(isect_loc)) + center_1)) From e6d5d8b90431fc030cf8f6f41b8cc97cab4bf61c Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Sun, 22 Oct 2017 19:36:33 +0200 Subject: [PATCH 070/118] Fix fracture intersection finder for boundary points. When the new polygon segment intersection algorithm was used, points on the boundary were found by the method aimed at internal points. Missing information the segment is on the boundary, the further meshing algorithm failed. --- src/porepy/fracs/fractures.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 3d06991c3d..6729d0a782 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -368,8 +368,14 @@ def intersects(self, other, tol, check_point_contact=True): # segment of the other. # Compute intersections, with both polygons as first argument - isect_self_other = cg.polygon_segment_intersect(self.p, other.p, tol=tol) - isect_other_self = cg.polygon_segment_intersect(other.p, self.p, tol=tol) + # Do not include boundary points of the polygons here, these will be + # processed later + isect_self_other = cg.polygon_segment_intersect(self.p, other.p, + tol=tol, + include_bound_pt=False) + isect_other_self = cg.polygon_segment_intersect(other.p, self.p, + tol=tol, + include_bound_pt=False) # Process data if isect_self_other is not None: @@ -397,8 +403,8 @@ def intersects(self, other, tol, check_point_contact=True): # There at at the most two intersection points between the fractures # (assuming convexity). If two interior points are found, we can simply # cut it short here. - if int_points.shape[1] == 2: - return int_points, on_boundary_self, on_boundary_other +# if int_points.shape[1] == 2: +# return int_points, on_boundary_self, on_boundary_other #### # Next, check for intersections between the polygon boundaries @@ -1861,7 +1867,7 @@ def dist_p(a, b): else: isect_f.append(i.first.index) - # Assuming the fracture is convex, the closest point for + # Assuming the fracture is convex, the closest point for dist, cp = cg.dist_points_segments(i.coord, f.p, np.roll(f.p, 1, axis=1)) # Insert a (candidate) point only at the segment closest to the From d6aa00b6d0ab421cc06a344ee843c4250a68ab82 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 23 Oct 2017 11:58:26 +0200 Subject: [PATCH 071/118] Gitignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index 35052244b3..fe008156df 100644 --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,6 @@ *.coverage .cache/* bin/* +*.pvd +*.so +*.png From 243fbad51a44d57eb71c6ebb9cb10deca067292e Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 23 Oct 2017 15:47:25 +0200 Subject: [PATCH 072/118] Improvement of fracture intersection for segments within planes. Allow for more than two intersection points if the effective line is more or less colinear. --- src/porepy/fracs/fractures.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 6729d0a782..6e062f9b81 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -397,8 +397,21 @@ def intersects(self, other, tol, check_point_contact=True): if int_points.shape[1] > 1: int_points, _, _ \ = setmembership.unique_columns_tol(int_points, tol=tol) - # There should be at most two of these points - assert int_points.shape[1] <= 2 + + # There should be at most two of these points. + # In some cases, likely involving extrusion, several segments may lay + # essentially in the fracture plane, producing more than two segments. + # Thus, if more than two colinear points are found, pick out the first + # and last one. + if int_points.shape[1] > 2: + if cg.is_collinear(int_points, tol): + sort_ind = cg.argsort_point_on_line(int_points, tol) + int_points = int_points[:, [sort_ind[0], sort_ind[-1]]] + else: + # This is a bug + raise ValueError(''' Found more than two intersection between + fracture polygons. + ''') # There at at the most two intersection points between the fractures # (assuming convexity). If two interior points are found, we can simply From 76708c9d8e72cf4c9df37a4d173d04012fcaf4ab Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 23 Oct 2017 15:49:13 +0200 Subject: [PATCH 073/118] Remove check on fracture convexity in extrusion. Save computational time. --- src/porepy/fracs/extrusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 8bededb005..2d54bb30e1 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -540,7 +540,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, # cloud. # We add the points after elimination, to ensure that the points on the # plane are present in the final fracture. - main_frac.add_points(isect_pt) + main_frac.add_points(isect_pt, check_convexity=False) # If the main fracture is too large compared to the other, the cut line # will extend beyond the confining plane. In these cases, compute the From 1c245d0b5ab9589128bdd06e3f57cfca5af62307 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 08:48:10 +0200 Subject: [PATCH 074/118] Bugfix in polygon intersection algorithm. Logical error and poor treatment of case with segments embedded in plane of other polygon. --- src/porepy/utils/comp_geom.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 5459c6ec8a..9545378082 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1287,12 +1287,12 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8, include_bound_pt=True): # Obtain the rotation matrix that projects p1 to the xy-plane rot_p_1 = project_plane_matrix(poly_1) irot = rot_p_1.transpose() - poly_1_xy = rot_p_1.dot(poly_1) + poly_1_rot = rot_p_1.dot(poly_1) # Sanity check: The points should lay on a plane - assert np.all(np.abs(poly_1_xy[2]) < tol) + assert np.all(np.abs(poly_1_rot[2]) < tol) # Drop the z-coordinate - poly_1_xy = poly_1_xy[:2] + poly_1_xy = poly_1_rot[:2] # Make sure the xy-polygon is ccw. if not is_ccw_polygon(poly_1_xy): @@ -1400,11 +1400,14 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8, include_bound_pt=True): isect_loc = np.empty((2, 0)) p1 = both_pts[:, 0] p2 = both_pts[:, 1] - if include_bound_pt: - poly_1_start = poly_1 - poly_1_end = np.roll(poly_1, 1, axis=1) + + # If a single internal point is found + if isect_loc.shape[1] == 1 or include_bound_pt: + poly_1_start = poly_1_rot + poly_1_end = np.roll(poly_1_rot, 1, axis=1) for j in range(poly_1.shape[1]): - ip = segments_intersect_3d(p1, p2, poly_1_start[:, j], + ip = segments_intersect_3d(p1, p2, + poly_1_start[:, j], poly_1_end[:, j]) if ip is not None: isect_loc = np.hstack((isect_loc, ip[:2])) From 5e98dcca6e86954ca6637d9e4b70825ca2b7fc51 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 08:49:48 +0200 Subject: [PATCH 075/118] Tweak of fracture extrusion to limit point intersections. Try to avoid cases where the fracture center is too close to z=0, this will give a point intersection. --- src/porepy/fracs/extrusion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 2d54bb30e1..3e0179cd77 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -253,8 +253,8 @@ def disc_radius_center(lengths, p0, p1, theta=None): # Angles of pi/2 will read to point contacts that cannot be handled # of the FractureNetwork. Point contacts also make little physical # sense, so we vaoid them. - hit = np.abs(rnd - 0.5) < 0.05 - rnd[hit] += 0.05 + hit = np.abs(rnd - 0.5) < 0.1 + rnd[hit] += 0.1 theta = np.pi * (0.1 + 0.8 * rnd) radius = 0.5 * lengths / np.sin(theta) From 31e93e6cbf6be2057e5ed25f5ca0edc828200339 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 10:18:21 +0200 Subject: [PATCH 076/118] Update of functionality for comp goem polygon segment intersection. It turned out to be beneficial to include intersection points close to the polygon boundary. This commit has the side-effect that the method in some cases will effectively identify some boundary segment intersections as internal to one of the polygons. --- src/porepy/utils/comp_geom.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index 9545378082..c74b2b7bd4 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1254,7 +1254,10 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8, include_bound_pt=True): Find intersections between polygons embeded in 3D. The intersections are defined as between the interior of the first polygon - and the boundary of the second. + and the boundary of the second, although intersections on the boundary of + both polygons can also be picked up sometimes. If you need to distinguish + between the two, the safer option is to also call + polygon_boundary_intersect(), and compare the results. Parameters: poly_1 (np.ndarray, 3xn1): Vertexes of polygon, assumed ordered as cw or @@ -1376,7 +1379,11 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8, include_bound_pt=True): # Check if the first polygon encloses the point. If the # intersection is on the border, this will not be detected. - if is_inside_polygon(poly_1_xy, p_00, tol=tol): + # For fracture intersection, in particular for intrusion / + # Y-intersections, it turned out to be critical to have True + # as default here. This comes to the price of picking up + # some boundary points as well. + if is_inside_polygon(poly_1_xy, p_00, tol=tol, default=True): # Back to physical coordinates by 1) expand to 3D, 2) # inverse rotation, 3) translate to original coordinate. isect = np.hstack((isect, irot.dot(_to3D(p_00)) + From 7f51dbe34075411e0e2809f2fda1aaeec31f650f Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 10:22:09 +0200 Subject: [PATCH 077/118] Fracture intersection finder ignores single point intersections. Fingers crossed gmsh will agree with this. --- src/porepy/fracs/fractures.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 6e062f9b81..fb739652ab 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -433,9 +433,15 @@ def intersects(self, other, tol, check_point_contact=True): # points if len(bound_sect_self_other) == 0 and len(bound_sect_other_self) == 0: if check_point_contact and int_points.shape[1] == 1: - raise ValueError('Contact in a single point not allowed') + # Point contacts are not implemented. Give a warning, return no + # interseciton, and hope the meshing software is merciful + logger.warning("""Found a point contact between fracture %i + and %i at (%.5f, %.5f, %.5f)""", self.index, + other.index, *int_points) + return np.empty((3, 0)), on_boundary_self, on_boundary_other # None of the intersection points lay on the boundary - return int_points, on_boundary_self, on_boundary_other + else: + 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 = \ From 20b8015c75c88040e460d091e4c144659c81afe3 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 10:24:39 +0200 Subject: [PATCH 078/118] Fracture intersection more attentive to points close to polygon boundary. Allow for these, and treat possible duplicates further down. --- src/porepy/fracs/fractures.py | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index fb739652ab..0577b1ebe1 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -372,10 +372,10 @@ def intersects(self, other, tol, check_point_contact=True): # processed later isect_self_other = cg.polygon_segment_intersect(self.p, other.p, tol=tol, - include_bound_pt=False) + include_bound_pt=True) isect_other_self = cg.polygon_segment_intersect(other.p, self.p, tol=tol, - include_bound_pt=False) + include_bound_pt=True) # Process data if isect_self_other is not None: @@ -530,19 +530,22 @@ def intersects(self, other, tol, check_point_contact=True): # Special case if a segment is cut as one interior point (by its # vertex), and a boundary point on the same segment, this is a boundary # segment - if int_points.shape[1] == 1 and bound_pt_self.shape[1] == 1: - is_vert, vi = self.is_vertex(int_points) - seg_i = bound_sect_self_other[0][0] - nfp = self.p.shape[1] - if is_vert and (vi == seg_i or (seg_i+1) % nfp == vi): - on_boundary_self = True - if int_points.shape[1] == 1 and bound_pt_other.shape[1] == 1: - is_vert, vi = other.is_vertex(int_points) - seg_i = bound_sect_other_self[0][0] - nfp = other.p.shape[1] - if is_vert and (vi == seg_i or (seg_i+1) % nfp == vi): - on_boundary_other = True - + if bound_pt_self.shape[1] == 1: + for ind in range(int_points.shape[1]): + is_vert, vi = self.is_vertex(int_points[:, ind]) + seg_i = bound_sect_self_other[0][0] + nfp = self.p.shape[1] + if is_vert and (vi == seg_i or (seg_i+1) % nfp == vi): + on_boundary_self = True + if bound_pt_other.shape[1] == 1: + for ind in range(int_points.shape[1]): + is_vert, vi = other.is_vertex(int_points[:, ind]) + seg_i = bound_sect_other_self[0][0] + nfp = other.p.shape[1] + if is_vert and (vi == seg_i or (seg_i+1) % nfp == vi): + on_boundary_other = True + + bound_pt, _, _ = setmembership.unique_columns_tol(bound_pt, tol=tol) return bound_pt, on_boundary_self, on_boundary_other def _process_segment_isect(self, isect_bound, poly, tol): From 78649ec7d1476ba82bc0871909c718143f0e5b75 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 10:12:59 +0200 Subject: [PATCH 079/118] Update unit test for polygon segment intersection after code updates. --- .../unit/test_polygon_segment_intersection.py | 29 ++++++++++++++----- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/test/unit/test_polygon_segment_intersection.py b/test/unit/test_polygon_segment_intersection.py index 6780917cd0..db3293de21 100644 --- a/test/unit/test_polygon_segment_intersection.py +++ b/test/unit/test_polygon_segment_intersection.py @@ -76,18 +76,27 @@ def test_extension_would_intersect(self): assert isect is None def test_segments_intersect(self): - # Test where the planes intersect in a way where segments only touches - # segments, which does not qualify as intersection. + # Test where the planes intersect in a way where vertex only touches + # vertex. This is now updated to count as an intersection. p_1, _, _, p_4 = self.setup_polygons() + isect = cg.polygon_segment_intersect(p_1, p_4) - assert isect is None + + isect_known_1 = np.array([[0, 0, -1]]).T + isect_known_2 = np.array([[0, 0, 1]]).T + + assert np.min(np.sum(np.abs(isect - isect_known_1), axis=0)) < 1e-5 + assert np.min(np.sum(np.abs(isect - isect_known_2), axis=0)) < 1e-5 + # Also try the other way around isect = cg.polygon_segment_intersect(p_4, p_1) - assert isect is None + assert np.min(np.sum(np.abs(isect - isect_known_1), axis=0)) < 1e-5 + assert np.min(np.sum(np.abs(isect - isect_known_2), axis=0)) < 1e-5 def test_issue_16(self): # Test motivated from debuging Issue #16 (GitHub) - # Two polygons meet in a common vertex; return should be None + # After updates of the code, we should find both intersection at vertex, + # and internal to segments of both polygons frac1 = np.array([[1, 2, 4], [1, 4, 1], [2, 2, 2]]) @@ -96,12 +105,16 @@ def test_issue_16(self): [2, 4, 1], [1, 2, 4]]) - isect_known = np.array([[2], [5/3], [2]]) + # Segment + isect_known_1 = np.array([[2, 5/3, 2]]).T + isect_known_2 = np.array([[2, 4, 2]]).T isect = cg.polygon_segment_intersect(frac1, frac2) - assert np.allclose(isect, isect_known) + assert np.min(np.sum(np.abs(isect - isect_known_1), axis=0)) < 1e-5 + assert np.min(np.sum(np.abs(isect - isect_known_2), axis=0)) < 1e-5 isect = cg.polygon_segment_intersect(frac1[:, [0, 2, 1]], frac2) - assert np.allclose(isect, isect_known) + assert np.min(np.sum(np.abs(isect - isect_known_1), axis=0)) < 1e-5 + assert np.min(np.sum(np.abs(isect - isect_known_2), axis=0)) < 1e-5 def test_segment_in_polygon_plane(self): # One segment lies fully within a plane From de359b3a1946e002f38043c187be186fc0d627f9 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 10:25:57 +0200 Subject: [PATCH 080/118] Remove shortcuts and assertions in fracture intersection finder. These no longer are safe, as the distinction between internal and boundary intersections is somewhat blurred. --- src/porepy/fracs/fractures.py | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 0577b1ebe1..a0888bb42c 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -413,12 +413,6 @@ def intersects(self, other, tol, check_point_contact=True): fracture polygons. ''') - # There at at the most two intersection points between the fractures - # (assuming convexity). If two interior points are found, we can simply - # cut it short here. -# if int_points.shape[1] == 2: -# return int_points, on_boundary_self, on_boundary_other - #### # Next, check for intersections between the polygon boundaries bound_sect_self_other = cg.polygon_boundaries_intersect(self.p, @@ -455,9 +449,7 @@ def intersects(self, other, tol, check_point_contact=True): # Convex polygons can intersect each other in at most two points (and a # line inbetween) - if int_points.shape[1] > 1: - assert bound_pt_self.shape[1] == 0 and bound_pt_other.shape[1] == 0 - elif int_points.shape[1] == 1: + if int_points.shape[1] == 1: # There should be exactly one unique boundary point. bp = np.hstack((bound_pt_self, bound_pt_other)) if bp.shape[1] == 0: @@ -483,8 +475,6 @@ def intersects(self, other, tol, check_point_contact=True): # If a segment runs through the polygon, there should be no interior points. # This may be sensitive to tolerances, should be a useful test to gauge # that. - if self_cuts_through or other_cuts_through: - assert int_points.shape[1] == 0 # Storage for intersection points located on the boundary bound_pt = [] From 326abc773a4132ab51a868d69fe6c7ac764ea14d Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 10:21:21 +0200 Subject: [PATCH 081/118] Pass on kwargs between fracture extrusion methods --- src/porepy/fracs/extrusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 3e0179cd77..d1d0220df9 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -431,7 +431,7 @@ def impose_inlcine(fracs, exposure, frac_family=None, family_mean_incline=None, return all_ang def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, - recompute_center=False): + recompute_center=False, **kwargs): """ Cut a fracture by a plane, and confine it to one side of the plane. Intended use is to confine abutting fractures (T-intersections) to one side From 683a8a7305cf9c22449da373fe0b3cfea024c42b Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 10:43:23 +0200 Subject: [PATCH 082/118] Disband consistency test of cell-face relation in one fracture network test. What seems like a bug in older versions of Gmsh leads to faces on fracture surfaces that are not matched by 2d cells. Disband the test of consistency between faces and cells, and instead use the test to check that fracture intersection algorithm works well. --- examples/example3/test_simple_networks.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/examples/example3/test_simple_networks.py b/examples/example3/test_simple_networks.py index 75b38fefdf..2749f1a7b1 100644 --- a/examples/example3/test_simple_networks.py +++ b/examples/example3/test_simple_networks.py @@ -199,7 +199,12 @@ def test_T_intersection_within_plane(**kwargs): domain = {'xmin':-1, 'xmax': 2, 'ymin': -2, 'ymax': 2, 'zmin': -1, 'zmax':2} - grids = meshing.simplex_grid([p1, p2], domain) + # This test, when used with certain versions of gmsh (<2.15?) gives a + # a mismatch between 2d cells and 3d faces on fracture surfaces. The bug + # can be located in the .msh-file. To function as a test, we disband the + # test of cell-face relations. + grids = meshing.simplex_grid([p1, p2], domain, + ensure_matching_face_cell=False) def test_T_intersection_one_outside_plane(**kwargs): p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T From 4281f11be237b7cee8dda00df9bc542120c35e42 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 20:36:07 +0200 Subject: [PATCH 083/118] Remove unnecessary convexity check when fitting Fracture to domain boundary. We know that the planes are convex anyhow. --- src/porepy/fracs/fractures.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index a0888bb42c..8b0e64defc 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -688,20 +688,26 @@ def outside_box(p, bound_i): # by the box, and extend the two others to cover segments that run # through the extension of the plane only. west = Fracture(np.array([[x0_box, x0_box, x0_box, x0_box], - [y0, y1, y1, y0], [z0, z0, z1, z1]])) + [y0, y1, y1, y0], [z0, z0, z1, z1]]), + check_convexity=False) east = Fracture(np.array([[x1_box, x1_box, x1_box, x1_box], [y0, y1, y1, y0], - [z0, z0, z1, z1]])) + [z0, z0, z1, z1]]), + check_convexity=False) south = Fracture(np.array([[x0, x1, x1, x0], [y0_box, y0_box, y0_box, y0_box], - [z0, z0, z1, z1]])) + [z0, z0, z1, z1]]), + check_convexity=False) north = Fracture(np.array([[x0, x1, x1, x0], [y1_box, y1_box, y1_box, y1_box], - [z0, z0, z1, z1]])) + [z0, z0, z1, z1]]), + check_convexity=False) bottom = Fracture(np.array([[x0, x1, x1, x0], [y0, y0, y1, y1], - [z0_box, z0_box, z0_box, z0_box]])) + [z0_box, z0_box, z0_box, z0_box]]), + check_convexity=False) top = Fracture(np.array([[x0, x1, x1, x0], [y0, y0, y1, y1], - [z1_box, z1_box, z1_box, z1_box]])) + [z1_box, z1_box, z1_box, z1_box]]), + check_convexity=False) # Collect in a list to allow iteration bound_planes = [west, east, south, north, bottom, top] From ae931225f5fb5e4cc19b00e9cb07472c4ac81e8c Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 20:38:04 +0200 Subject: [PATCH 084/118] Bugfix in FractureNetwork distance computations. --- src/porepy/fracs/fractures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 8b0e64defc..e7a048f921 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1937,7 +1937,7 @@ def dist_p(a, b): d_1 = dist_p(cp_f[:, mi], f.p[:, si]) d_2 = dist_p(cp_f[:, mi], f.p[:, (si+1)%nfp]) if d_1 > h_min and d_2 > h_min: - np.insert(f.p, (si+1)%nfp, cp_f[:, si], axis=1) + np.insert(f.p, (si+1)%nfp, cp_f[:, mi], axis=1) # Finally, cover the case where the smallest distance is given # by a point. Points with false in_point should already be From 75acf0adb0c78fdaae27c5f49ffa47d94407a50c Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 25 Oct 2017 20:43:01 +0200 Subject: [PATCH 085/118] Bugfix in FractureNetwork intersection algorithm. When the original point definition was used, and the Fracture had been modified, this gave rise to intersections outside the (new) fracture polygon. --- src/porepy/fracs/fractures.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index e7a048f921..d5d2cbf412 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1209,8 +1209,8 @@ def _point_and_edge_lists(self): # will deal with coinciding points later. for fi, frac in enumerate(self._fractures): num_p = all_p.shape[1] - num_p_loc = frac.orig_p.shape[1] - all_p = np.hstack((all_p, frac.orig_p)) + num_p_loc = frac.p.shape[1] + all_p = np.hstack((all_p, frac.p)) loc_e = num_p + np.vstack((np.arange(num_p_loc), (np.arange(num_p_loc) + 1) % num_p_loc)) From 253b88bfaf7e74a03f069681869fd140dd823617 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 26 Oct 2017 09:07:48 +0200 Subject: [PATCH 086/118] Tuning of Fracture projection to 2d plane tolerance. --- src/porepy/fracs/fractures.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index d5d2cbf412..cd8974db2d 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -397,7 +397,7 @@ def intersects(self, other, tol, check_point_contact=True): if int_points.shape[1] > 1: int_points, _, _ \ = setmembership.unique_columns_tol(int_points, tol=tol) - + # There should be at most two of these points. # In some cases, likely involving extrusion, several segments may lay # essentially in the fracture plane, producing more than two segments. @@ -1613,8 +1613,8 @@ def _points_2_plane(self, p_loc, edges_loc, p_ind_loc): p_2d = rot.dot(p_loc) extent = p_2d.max(axis=1) - p_2d.min(axis=1) - - assert extent[2] < np.max(extent[:2]) * self.tol * 5 + lateral_extent = np.maximum(np.max(extent[2]), 1) + assert extent[2] < lateral_extent * self.tol * 10 # Dump third coordinate p_2d = p_2d[:2] From 6dbfaf1a3a53451ec26e5905f4e27d669302b43a Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 27 Oct 2017 10:25:16 +0200 Subject: [PATCH 087/118] 2d intersection removal optionally drops snapping to grid. --- src/porepy/utils/comp_geom.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index c74b2b7bd4..d3b35e79ea 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -384,7 +384,7 @@ def _split_edge(vertices, edges, edge_ind, new_pt, **kwargs): #------------------------------------------------------**kwargs------------------------# -def _add_point(vertices, pt, tol=1e-3, **kwargs): +def _add_point(vertices, pt, tol=1e-3, snap=True, **kwargs): """ Add a point to a point set, unless the point already exist in the set. @@ -414,8 +414,9 @@ def _add_point(vertices, pt, tol=1e-3, **kwargs): nd = vertices.shape[0] # Before comparing coordinates, snap both existing and new point to the # underlying grid - vertices = snap_to_grid(vertices, **kwargs) - pt = snap_to_grid(pt, **kwargs) + if snap: + vertices = snap_to_grid(vertices, **kwargs) + pt = snap_to_grid(pt, **kwargs) new_pt = np.empty((nd, 0)) ind = [] @@ -442,7 +443,8 @@ def _add_point(vertices, pt, tol=1e-3, **kwargs): #-----------------------------------------------------------------------------# -def remove_edge_crossings(vertices, edges, tol=1e-3, verbose=0, **kwargs): +def remove_edge_crossings(vertices, edges, tol=1e-3, verbose=0, snap=True, + **kwargs): """ Process a set of points and connections between them so that the result is an extended point set and new connections that do not intersect. @@ -487,8 +489,9 @@ def remove_edge_crossings(vertices, edges, tol=1e-3, verbose=0, **kwargs): # Add tolerance to kwargs, this is later passed to split_edges, and further # on. kwargs['tol'] = tol - - vertices = snap_to_grid(vertices, **kwargs) + kwargs['snap'] = snap + if snap: + vertices = snap_to_grid(vertices, **kwargs) # Field used for debugging of edge splits. To see the meaning of the values # of each split, look in the source code of split_edges. @@ -617,7 +620,8 @@ def __min_dist(p): if new_pt is None: logger.debug('No intersection found') else: - new_pt = snap_to_grid(new_pt, tol=tol) + if snap: + new_pt = snap_to_grid(new_pt, tol=tol) # The case of segment intersections need special treatment. if new_pt.shape[-1] == 1: logger.debug('Found intersection (%.5f, %.5f)', new_pt[0], From a990fba9679bbe65d38ea2cae137a71482b44955 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 27 Oct 2017 10:30:36 +0200 Subject: [PATCH 088/118] Avoid snapping to grid when searching for fracture intersections. This seems a more robust approach - searching for identical points are rather handled by unique_columns_tol() --- src/porepy/fracs/fractures.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index cd8974db2d..60e6ddf13e 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1262,7 +1262,7 @@ def _uniquify_points_and_edges(self, all_p, edges, edges_2_frac, logger.info("""Uniquify points and edges, starting with %i points, %i edges""", all_p.shape[1], edges.shape[1]) - all_p = cg.snap_to_grid(all_p, tol=self.tol) +# all_p = cg.snap_to_grid(all_p, tol=self.tol) # We now need to find points that occur in multiple places p_unique, unique_ind_p, \ @@ -1387,8 +1387,9 @@ def _remove_edge_intersections(self, all_p, edges, edges_2_frac, # obtain a more robust algorithm. Not sure about how to do this # consistent. p_new, edges_new = cg.remove_edge_crossings(p_2d, edges_2d, - tol=self.tol*5, - verbose=self.verbose) + tol=self.tol, + verbose=self.verbose, + snap=False) # Then, patch things up by converting new points to 3D, # From the design of the functions in cg, we know that new points From 081bee7378437fa4880bd38a3c7f8bc29123ed36 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 27 Oct 2017 11:46:36 +0200 Subject: [PATCH 089/118] Improved treatment of polygon segments intersecting near boundaries. In cases where the intersection point is just outside the boundary point, use the projection of the intersection point onto the polygon. --- src/porepy/utils/comp_geom.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index d3b35e79ea..f21ea01c8d 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1380,18 +1380,17 @@ def polygon_segment_intersect(poly_1, poly_2, tol=1e-8, include_bound_pt=True): # Representation as point p_00 = np.array([x0, y0]).reshape((-1, 1)) - # Check if the first polygon encloses the point. If the - # intersection is on the border, this will not be detected. - - # For fracture intersection, in particular for intrusion / - # Y-intersections, it turned out to be critical to have True - # as default here. This comes to the price of picking up - # some boundary points as well. - if is_inside_polygon(poly_1_xy, p_00, tol=tol, default=True): - # Back to physical coordinates by 1) expand to 3D, 2) - # inverse rotation, 3) translate to original coordinate. - isect = np.hstack((isect, irot.dot(_to3D(p_00)) + - center_1)) + # Check if the first polygon encloses the point. When applied + # to fracture intersections of T-type (segment embedded in the + # plane of another fracture), it turned out to be useful to be + # somewhat generous with the definition of the intersection. + # Therefore, allow for intersections that are slightly outside + # the polygon, and use the projection onto the polygon. + dist, cp, ins = dist_points_polygon(_to3D(p_00), + _to3D(poly_1_xy)) + if (dist[0] < tol and include_bound_pt) or dist[0] < 1e-12: + isect = np.hstack((isect, irot.dot(cp) + center_1)) + elif np.abs(pt_1[2]) < tol and np.abs(pt_2[2]) < tol: # The segment lies completely within the polygon plane. both_pts = np.vstack((pt_1, pt_2)).T From f3a2268a938cca10d61a29df59f9b64a9fe475b5 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 30 Oct 2017 12:18:10 +0100 Subject: [PATCH 090/118] Bugfix in Fracture intersection definition of boundary segments. --- src/porepy/fracs/fractures.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 60e6ddf13e..88a0b52c66 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -435,6 +435,25 @@ def intersects(self, other, tol, check_point_contact=True): return np.empty((3, 0)), on_boundary_self, on_boundary_other # None of the intersection points lay on the boundary else: + def point_on_segment(ip, poly): + # Check if a set of points are located on a single segment + if ip.shape[1] == 0: + return False + start = poly + end = np.roll(poly, 1, axis=1) + for si in range(start.shape[1]): + dist, cp = cg.dist_points_segments(ip, start[:, si], + end[:, si]) + if np.all(dist < tol): + return True + return False + + # The 'interior' points can still be on the boundary (naming + # of variables should be updated). The points form a boundary + # segment if they all lie on the a single segment of the + # fracture. + on_boundary_self = point_on_segment(int_points, self.p) + on_boundary_other = point_on_segment(int_points, other.p) return int_points, on_boundary_self, on_boundary_other # Else, we have boundary intersection, and need to process them From c37aa30ea989dd795bfcc51284defda067848f0c Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 30 Oct 2017 12:20:05 +0100 Subject: [PATCH 091/118] Minor fix of FractureNetwork method to get intersections of specified fracture. --- src/porepy/fracs/fractures.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 88a0b52c66..ac29fed571 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1017,12 +1017,15 @@ def intersections_of_fracture(self, frac): np.array (Intersection): Array of intersections """ - fi = frac.index + if isinstance(frac, int): + fi = frac + else: + fi = frac.index frac_arr = [] for i in self.intersections: if i.coord.size == 0: continue - if i.first == frac or i.second == frac: + if i.first.index == fi or i.second.index == fi: frac_arr.append(i) return frac_arr From e45fa2e4eaf1dd6e74bb0d310ad28e34c7986916 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 30 Oct 2017 12:20:59 +0100 Subject: [PATCH 092/118] Tuning of threshold in FractureNetwork test of plane surfaces. Not sure if this is a good option, but it works well for now. --- src/porepy/fracs/fractures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index ac29fed571..2c145cb7ab 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1637,7 +1637,7 @@ def _points_2_plane(self, p_loc, edges_loc, p_ind_loc): extent = p_2d.max(axis=1) - p_2d.min(axis=1) lateral_extent = np.maximum(np.max(extent[2]), 1) - assert extent[2] < lateral_extent * self.tol * 10 + assert extent[2] < lateral_extent * self.tol * 30 # Dump third coordinate p_2d = p_2d[:2] From 1e23d939d0ed58f30c1e6291d7a8ff41384ee5ca Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 30 Oct 2017 14:44:16 +0100 Subject: [PATCH 093/118] Fix of T-type Fracture intersection. New tests. Special case when the intersection points coincide with the vertexes of the abutting fracture in a T-configuration. --- examples/example3/test_simple_networks.py | 24 +++++++++++++++++++++++ src/porepy/fracs/fractures.py | 21 ++++++++++++++++---- 2 files changed, 41 insertions(+), 4 deletions(-) diff --git a/examples/example3/test_simple_networks.py b/examples/example3/test_simple_networks.py index 2749f1a7b1..d2def7aeb2 100644 --- a/examples/example3/test_simple_networks.py +++ b/examples/example3/test_simple_networks.py @@ -221,3 +221,27 @@ def test_T_intersection_both_outside_plane(**kwargs): domain = {'xmin':-1, 'xmax': 2, 'ymin': -2, 'ymax': 2, 'zmin': -1, 'zmax':2} grids = meshing.simplex_grid([p1, p2], domain) + +def test_T_intersection_both_on_boundary(**kwargs): + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + p2 = np.array([[0., 0.5, 0], [1, 0.5, 0], [0.5, 0.5, 1.]]).T + + domain = {'xmin':-1, 'xmax': 2, 'ymin': -2, 'ymax': 2, 'zmin': -1, + 'zmax':2} + grids = meshing.simplex_grid([p1, p2], domain) + +def test_T_intersection_one_boundary_one_outside(**kwargs): + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + p2 = np.array([[-0.2, 0.5, 0], [1, 0.5, 0], [0.5, 0.5, 1.]]).T + + domain = {'xmin':-1, 'xmax': 2, 'ymin': -2, 'ymax': 2, 'zmin': -1, + 'zmax':2} + grids = meshing.simplex_grid([p1, p2], domain) + +def test_T_intersection_one_boundary_one_inside(**kwargs): + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + p2 = np.array([[0.2, 0.5, 0], [1, 0.5, 0], [0.5, 0.5, 1.]]).T + + domain = {'xmin':-1, 'xmax': 2, 'ymin': -2, 'ymax': 2, 'zmin': -1, + 'zmax':2} + grids = meshing.simplex_grid([p1, p2], domain) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 2c145cb7ab..35e8d80585 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -458,11 +458,12 @@ def point_on_segment(ip, poly): # Else, we have boundary intersection, and need to process them bound_pt_self, self_segment, _, self_cuts_through = \ - self._process_segment_isect(bound_sect_self_other, self.p, tol) + self._process_segment_isect(bound_sect_self_other, self.p, tol, + other.p) bound_pt_other, other_segment, _, other_cuts_through = \ self._process_segment_isect(bound_sect_other_self, other.p, - tol) + tol, self.p) # Run some sanity checks @@ -557,7 +558,7 @@ def point_on_segment(ip, poly): bound_pt, _, _ = setmembership.unique_columns_tol(bound_pt, tol=tol) return bound_pt, on_boundary_self, on_boundary_other - def _process_segment_isect(self, isect_bound, poly, tol): + def _process_segment_isect(self, isect_bound, poly, tol, other_poly): """ Helper function to interpret result from polygon_boundaries_intersect @@ -582,11 +583,23 @@ def _process_segment_isect(self, isect_bound, poly, tol): # Counter for how many times each segment is found num_occ = np.zeros(poly.shape[1]) + # Number of vertexes on the other polygon - used below to identify + # whether the intersection coincides with these. + num_pt_other = other_poly.shape[1] + # Loop over all intersections, process information for bi in isect_bound: # Index of the intersecting segments - num_occ[bi[0]] += 1 + + # Special case: For T-type intersection where the vertexes of one + # polygon lies on the segment of the other, they do not count + # towards cutting through the segment. + expanded = np.hstack((other_poly, bi[2])) + other_unique, _, _ = setmembership.unique_columns_tol(expanded, + tol=tol) + if other_unique.shape[1] > num_pt_other: + num_occ[bi[0]] += 1 # Coordinates of the intersections ip = bi[2] From e765a57ea033ea403a1c19ced0d29f51ccba32e1 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 30 Oct 2017 15:42:13 +0100 Subject: [PATCH 094/118] Meshing of fractured domains reports progress --- src/porepy/fracs/meshing.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/meshing.py b/src/porepy/fracs/meshing.py index 5026750efb..4243ce633a 100644 --- a/src/porepy/fracs/meshing.py +++ b/src/porepy/fracs/meshing.py @@ -9,6 +9,7 @@ import numpy as np import scipy.sparse as sps import warnings +import time from porepy.fracs import structured, simplex, split_grid from porepy.grids.grid_bucket import GridBucket @@ -16,7 +17,7 @@ from porepy.utils import setmembership, mcolon -def simplex_grid(fracs=None, domain=None, network=None, **kwargs): +def simplex_grid(fracs=None, domain=None, network=None, verbose=0, **kwargs): """ Main function for grid generation. Creates a fractured simiplex grid in 2 or 3 dimensions. @@ -70,6 +71,10 @@ def simplex_grid(fracs=None, domain=None, network=None, **kwargs): else: raise ValueError('simplex_grid only supported for 2 or 3 dimensions') + if verbose > 0: + print('Construct mesh') + tm_msh = time.time() + tm_tot = time.time() # Call relevant method, depending on grid dimensions. if ndim == 2: assert fracs is not None, '2d requires definition of fractures' @@ -87,15 +92,38 @@ def simplex_grid(fracs=None, domain=None, network=None, **kwargs): grids = simplex.tetrahedral_grid(fracs, domain, network, **kwargs) else: raise ValueError('Only support for 2 and 3 dimensions') + + if verbose > 0: + print('Done. Elapsed time ' + str(time.time() - tm_msh)) + # Tag tip faces tag_faces(grids) # Assemble grids in a bucket + + if verbose > 0: + print('Assemble in bucket') + tm_bucket = time.time() gb = assemble_in_bucket(grids, **kwargs) + if verbose > 0: + print('Done. Elapsed time ' + str(time.time() - tm_bucket)) + print('Compute geometry') + tm_geom = time.time() + gb.compute_geometry() # Split the grids. + if verbose > 0: + print('Done. Elapsed time ' + str(time.time() - tm_geom)) + print('Split fractures') + tm_split = time.time() split_grid.split_fractures(gb, **kwargs) + if verbose > 0: + print('Done. Elapsed time ' + str(time.time() - tm_split)) gb.assign_node_ordering() + + if verbose > 0: + print('Mesh construction completed. Total time ' + str(time.time() - tm_tot)) + return gb #------------------------------------------------------------------------------# From 064846b03c1415c1c021720a1b8532f4b6bde2e8 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 30 Oct 2017 15:55:35 +0100 Subject: [PATCH 095/118] Minor fix in comp_geom. Scale determinant with length of fractures before declearing error. --- src/porepy/utils/comp_geom.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/utils/comp_geom.py b/src/porepy/utils/comp_geom.py index f21ea01c8d..6e6fc5ef96 100644 --- a/src/porepy/utils/comp_geom.py +++ b/src/porepy/utils/comp_geom.py @@ -1799,7 +1799,7 @@ def dist_two_segments(s1_start, s1_end, s2_start, s2_end): dot_2_starts = d2.dot(d_starts) discr = dot_1_1 * dot_2_2 - dot_1_2 ** 2 # Sanity check - assert discr >= -SMALL_TOLERANCE + assert discr >= -SMALL_TOLERANCE * dot_1_1 * dot_2_2 sc = sN = sD = discr tc = tN = tD = discr From 099e8f5fa77756143d51a743a4b802b34b166df9 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 31 Oct 2017 12:02:58 +0100 Subject: [PATCH 096/118] Tried to fix #54. Problems related to T-intersection, with interseciton through fracture vertexes. --- examples/example3/test_simple_networks.py | 18 ++++++++++ src/porepy/fracs/fractures.py | 43 ++++++++++++++++------- 2 files changed, 48 insertions(+), 13 deletions(-) diff --git a/examples/example3/test_simple_networks.py b/examples/example3/test_simple_networks.py index d2def7aeb2..acf1103393 100644 --- a/examples/example3/test_simple_networks.py +++ b/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) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 35e8d80585..7287164a39 100644 --- a/src/porepy/fracs/fractures.py +++ b/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: From 274e73d22d7a2cebddcf8d006d3164e5f4a90c61 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 31 Oct 2017 12:26:51 +0100 Subject: [PATCH 097/118] FractureNetwork distance computation outputs minimal distance between points. Can be used to guide mesh size parameters. --- src/porepy/fracs/fractures.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 7287164a39..a139c1ad37 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1902,6 +1902,7 @@ def _determine_mesh_size(self, **kwargs): num_pts = p.shape[1] dist = cg.dist_pointset(p, max_diag=True) mesh_size = np.min(dist, axis=1) + print('Minimal distance between points encountered is ' + str(np.min(dist))) mesh_size = np.maximum(mesh_size, self.h_min * np.ones(num_pts)) mesh_size = np.minimum(mesh_size, self.h_ideal * np.ones(num_pts)) From ea644c773d4780b31775cb5161a29f90b90d8db1 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 31 Oct 2017 12:52:34 +0100 Subject: [PATCH 098/118] Minor updates to FractureNetwork. Renamed function, comments. Also avoid that auxiliary points for mesh size control are added twice. --- src/porepy/fracs/fractures.py | 47 +++++++++++++++++++++++++---------- src/porepy/fracs/simplex.py | 2 +- 2 files changed, 35 insertions(+), 14 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index a139c1ad37..085b1c4bd0 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -1020,6 +1020,9 @@ def __init__(self, fractures=None, verbose=0, tol=1e-4): self.h_min = None self.h_ideal = None + # No auxiliary points have been added + self.auxiliary_points_added = False + def add(self, f): # Careful here, if we remove one fracture and then add, we'll have # identical indices. @@ -1896,6 +1899,7 @@ def _determine_mesh_size(self, **kwargs): return mesh_size, mesh_size_bound elif mode == 'distance': if self.h_min is None or self.h_ideal is None: + print('Found no information on mesh sizes. Returning') return None, None p = self.decomposition['points'] @@ -1911,11 +1915,39 @@ def _determine_mesh_size(self, **kwargs): else: raise ValueError('Unknown mesh size mode ' + mode) - def compute_distances(self, h_ideal=None, h_min=None): + def insert_auxiliary_points(self, h_ideal=None, h_min=None): + """ Insert auxiliary points on fracture edges. Used to guide gmsh mesh + size parameters. + + The function should only be called once to avoid insertion of multiple + layers of extra points, this will likely kill gmsh. + + The function is motivated by similar functionality for 2d domains, but + is considerably less mature. + + The function should be used in conjunction with _determine_mesh_size(), + called with mode='distance'. The ultimate goal is to set the mesh size + for geometrical points in Gmsh. To that end, this function inserts + additional points on the fracture boundaries. The mesh size is then + determined as the distance between all points in the fracture + description. + + Parameters: + h_ideal: Ideal mesh size. Will be added to all points that are + sufficiently far away from other points. + h_min: Minimal mesh size; we will make no attempts to enforce + even smaller mesh sizes upon Gmsh. + + """ + + + if self.auxiliary_points_added: + print('Auxiliary points already added. Returning.') + else: + self.auxiliary_points_added = True self.h_ideal = h_ideal self.h_min = h_min - isect_pt = np.zeros((3, 0), dtype=np.double) def dist_p(a, b): a = a.reshape((-1, 1)) @@ -1993,17 +2025,6 @@ def dist_p(a, b): if d_1 > h_min and d_2 > h_min: np.insert(f.p, (si+1)%nfp, cp_f[:, mi], axis=1) - # Finally, cover the case where the smallest distance is given - # by a point. Points with false in_point should already be - # covered by the above iteration over segments. - d, cp, in_poly = cg.dist_points_polygon(of.p, f.p) - for di, cpi, ip in zip(d, cp, in_poly): - # Closest point on segment is covered above - if not ip: - continue - if di < h_ideal: - pass - def distance_point_segment(self): pass diff --git a/src/porepy/fracs/simplex.py b/src/porepy/fracs/simplex.py index b9a831a77c..ceda87f3bb 100644 --- a/src/porepy/fracs/simplex.py +++ b/src/porepy/fracs/simplex.py @@ -107,7 +107,7 @@ def tetrahedral_grid(fracs=None, box=None, network=None, **kwargs): h_ideal = kwargs.get('h_ideal', None) h_min = kwargs.get('h_min', None) if h_ideal is not None and h_min is not None: - network.compute_distances(h_ideal, h_min) + network.insert_auxiliary_points(h_ideal, h_min) # In this case we need to recompute intersection decomposition anyhow. network.split_intersections() From c25dea28bf61428e14fa574527bbf1c810582643 Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 31 Oct 2017 14:31:49 +0100 Subject: [PATCH 099/118] Cut the running time of split grid with about 50% Biggest change was to call cell_nodes() only ones. I'm not 100% confident that this actually works. Also changed mcolon -> slice where the index only has size 1 --- src/porepy/fracs/meshing.py | 5 +++-- src/porepy/fracs/split_grid.py | 6 +++--- src/porepy/utils/sparse_mat.py | 7 ++++++- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/porepy/fracs/meshing.py b/src/porepy/fracs/meshing.py index 4243ce633a..f9281317b1 100644 --- a/src/porepy/fracs/meshing.py +++ b/src/porepy/fracs/meshing.py @@ -122,7 +122,8 @@ def simplex_grid(fracs=None, domain=None, network=None, verbose=0, **kwargs): gb.assign_node_ordering() if verbose > 0: - print('Mesh construction completed. Total time ' + str(time.time() - tm_tot)) + print('Mesh construction completed. Total time ' + + str(time.time() - tm_tot)) return gb @@ -396,7 +397,7 @@ def obtain_interdim_mappings(lg, fn, n_per_face, if not (np.all(is_mem) or np.all(~is_mem)): if ensure_matching_face_cell: raise ValueError( - '''Either all cells should have a corresponding face in a higher + '''Either all cells should have a corresponding face in a higher dim grid or no cells should have a corresponding face in a higher dim grid. This likely is related to gmsh behavior. ''') else: diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py index 8d10a590ea..c38d2c6b8f 100644 --- a/src/porepy/fracs/split_grid.py +++ b/src/porepy/fracs/split_grid.py @@ -387,16 +387,16 @@ def duplicate_nodes(g, nodes, offset): _, iv = sort_sub_list(g.face_nodes.indices, g.face_nodes.indptr) g.face_nodes = g.face_nodes.tocsr() - # Iterate over each internal node and split it according to the graph. # For each cell attached to the node, we check wich color the cell has. # All cells with the same color is then attached to a new copy of the # node. + cell_nodes = g.cell_nodes().tocsr() for node in nodes: # t_node takes into account the added nodes. t_node = node + node_count # Find cells connected to node - cells = sparse_mat.slice_indices(g.cell_nodes().tocsr(), t_node) + cells = sparse_mat.slice_indices(cell_nodes, node) # cell_nodes = g.cell_nodes().tocsr() # ind_ptr = cell_nodes.indptr # cells = cell_nodes.indices[ @@ -454,7 +454,7 @@ def duplicate_nodes(g, nodes, offset): def sort_sub_list(indices, indptr): ix = np.zeros(indices.size, dtype=int) for i in range(indptr.size - 1): - sub_ind = mcolon(indptr[i], indptr[i + 1]) + sub_ind = slice(indptr[i], indptr[i + 1]) loc_ix = np.argsort(indices[sub_ind]) ix[sub_ind] = loc_ix + indptr[i] indices = indices[ix] diff --git a/src/porepy/utils/sparse_mat.py b/src/porepy/utils/sparse_mat.py index 64835aef7e..c5b5bf9133 100644 --- a/src/porepy/utils/sparse_mat.py +++ b/src/porepy/utils/sparse_mat.py @@ -52,5 +52,10 @@ def slice_indices(A, slice_ind): rows = slice_indices(A, np.array([0,2,3])) """ assert A.getformat() == 'csc' or A.getformat() == 'csr' - indices = A.indices[mcolon(A.indptr[slice_ind], A.indptr[slice_ind + 1])] + if slice_ind.size == 1: + indices = A.indices[slice( + A.indptr[slice_ind], A.indptr[slice_ind + 1])] + else: + indices = A.indices[mcolon( + A.indptr[slice_ind], A.indptr[slice_ind + 1])] return indices From af4be97ea5d2d9746b082d353c91edae2a254502 Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 31 Oct 2017 14:55:45 +0100 Subject: [PATCH 100/118] Fixed unit-test that failed --- src/porepy/fracs/split_grid.py | 1 + src/porepy/utils/sparse_mat.py | 7 +++++-- test/unit/test_sparse_mat.py | 11 ++++++----- 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py index c38d2c6b8f..848bc6c505 100644 --- a/src/porepy/fracs/split_grid.py +++ b/src/porepy/fracs/split_grid.py @@ -396,6 +396,7 @@ def duplicate_nodes(g, nodes, offset): # t_node takes into account the added nodes. t_node = node + node_count # Find cells connected to node + cells = sparse_mat.slice_indices(cell_nodes, node) # cell_nodes = g.cell_nodes().tocsr() # ind_ptr = cell_nodes.indptr diff --git a/src/porepy/utils/sparse_mat.py b/src/porepy/utils/sparse_mat.py index c5b5bf9133..5bf988fec9 100644 --- a/src/porepy/utils/sparse_mat.py +++ b/src/porepy/utils/sparse_mat.py @@ -52,9 +52,12 @@ def slice_indices(A, slice_ind): rows = slice_indices(A, np.array([0,2,3])) """ assert A.getformat() == 'csc' or A.getformat() == 'csr' - if slice_ind.size == 1: + if isinstance(slice_ind, int): indices = A.indices[slice( - A.indptr[slice_ind], A.indptr[slice_ind + 1])] + A.indptr[int(slice_ind)], A.indptr[int(slice_ind + 1)])] + elif slice_ind.size == 1: + indices = A.indices[slice( + A.indptr[int(slice_ind)], A.indptr[int(slice_ind + 1)])] else: indices = A.indices[mcolon( A.indptr[slice_ind], A.indptr[slice_ind + 1])] diff --git a/test/unit/test_sparse_mat.py b/test/unit/test_sparse_mat.py index df43adebb0..3bcde967ca 100644 --- a/test/unit/test_sparse_mat.py +++ b/test/unit/test_sparse_mat.py @@ -28,8 +28,9 @@ def test_csr_slice(self): A = sps.csr_matrix(np.array([[0, 0, 0], [1, 0, 0], [0, 0, 3]])) + cols_0 = sparse_mat.slice_indices(A, np.array([0])) - cols_2 = sparse_mat.slice_indices(A, np.array([2])) + cols_2 = sparse_mat.slice_indices(A, 2) cols0_2 = sparse_mat.slice_indices(A, np.array([0, 1, 2])) assert cols_0.size == 0 @@ -41,8 +42,8 @@ def test_csc_slice(self): A = sps.csc_matrix(np.array([[0, 0, 0], [1, 0, 0], [0, 0, 3]])) - rows_0 = sparse_mat.slice_indices(A, np.array([0])) - rows_2 = sparse_mat.slice_indices(A, np.array([2])) + rows_0 = sparse_mat.slice_indices(A, np.array([0], dtype=int)) + rows_2 = sparse_mat.slice_indices(A, 2) rows0_2 = sparse_mat.slice_indices(A, np.array([0, 1, 2])) assert rows_0 == np.array([1]) @@ -67,8 +68,8 @@ def test_zero_columns(self): A0 = A.copy() A2 = A.copy() A0_2 = A.copy() - sparse_mat.zero_columns(A0, np.array([0])) - sparse_mat.zero_columns(A2, np.array([2])) + sparse_mat.zero_columns(A0, np.array([0], dtype=int)) + sparse_mat.zero_columns(A2, 2) sparse_mat.zero_columns(A0_2, np.array([0, 1, 2])) assert np.sum(A0 != A0_t) == 0 From 4f2737ec18857ffc6aa4acbb45225f4e5467f6f9 Mon Sep 17 00:00:00 2001 From: Runar Date: Tue, 31 Oct 2017 17:27:58 +0100 Subject: [PATCH 101/118] Improved sparse matrix slicing --- src/porepy/fracs/split_grid.py | 18 ++--------- src/porepy/utils/sparse_mat.py | 55 +++++++++++++++++++++++++++++++++- test/unit/test_sparse_mat.py | 55 ++++++++++++++++++++++++++++++++++ 3 files changed, 111 insertions(+), 17 deletions(-) diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py index 848bc6c505..8856cd706c 100644 --- a/src/porepy/fracs/split_grid.py +++ b/src/porepy/fracs/split_grid.py @@ -484,8 +484,8 @@ def find_cell_color(g, cells): """ c = np.sort(cells) # Local cell-face and face-node maps. - cf_sub, _ = __extract_submatrix(g.cell_faces, c) - child_cell_ind = np.array([-1] * g.num_cells, dtype=np.int) + cf_sub = sparse_mat.slice_mat(g.cell_faces, c) + child_cell_ind = -np.ones(g.num_cells, dtype=np.int) child_cell_ind[c] = np.arange(cf_sub.shape[1]) # Create a copy of the cell-face relation, so that we can modify it at @@ -537,17 +537,3 @@ def remove_nodes(g, rem): g.face_nodes = g.face_nodes[rows_to_keep, :] g.nodes = g.nodes[:, rows_to_keep] return g - - -def __extract_submatrix(mat, ind): - """ From a matrix, extract the column specified by ind. All zero columns - are stripped from the sub-matrix. Mappings from global to local row numbers - are also returned. - """ - sub_mat = mat[:, ind] - cols = sub_mat.indptr - rows = sub_mat.indices - data = sub_mat.data - unique_rows, rows_sub = np.unique(sub_mat.indices, - return_inverse=True) - return sps.csc_matrix((data, rows_sub, cols)), unique_rows diff --git a/src/porepy/utils/sparse_mat.py b/src/porepy/utils/sparse_mat.py index 5bf988fec9..028113e596 100644 --- a/src/porepy/utils/sparse_mat.py +++ b/src/porepy/utils/sparse_mat.py @@ -1,12 +1,15 @@ """ module for operations on sparse matrices """ +import numpy as np +import scipy.sparse as sps + from porepy.utils.mcolon import mcolon def zero_columns(A, cols): ''' - Function to zero out columns in matrix A. Note that this function does not + Function to zero out columns in matrix A. Note that this function does not change the sparcity structure of the matrix, it only changes the column values to 0 @@ -62,3 +65,53 @@ def slice_indices(A, slice_ind): indices = A.indices[mcolon( A.indptr[slice_ind], A.indptr[slice_ind + 1])] return indices + + +def slice_mat(A, ind): + """ + Function for slicing sparse matrix along rows or columns. + If A is a csc_matrix A will be sliced along columns, while if A is a + csr_matrix A will be sliced along the rows. + + Parameters + ---------- + A (scipy.sparse.csc/csr_matrix): A sparse matrix. + ind (np.array): Array containing indices to be sliced. + + Returns + ------- + A_sliced (scipy.sparse.csc/csr_matrix): The sliced matrix + if A is a csc_matrix A_sliced = A[:, ind] + if A is a csr_matrix A_slice = A[ind, :] + + Examples + -------- + A = sps.csc_matrix(np.eye(10)) + rows = slice_mat(A, np.array([0,2,3])) + """ + assert A.getformat() == 'csc' or A.getformat() == 'csr' + + if isinstance(ind, int): + N = 1 + indptr = np.zeros(2) + ind_slice = slice( + A.indptr[int(ind)], A.indptr[int(ind + 1)]) + elif ind.size == 1: + N = 1 + indptr = np.zeros(2) + ind_slice = slice( + A.indptr[int(ind)], A.indptr[int(ind + 1)]) + else: + N = ind.size + indptr = np.zeros(ind.size + 1) + ind_slice = mcolon( + A.indptr[ind], A.indptr[ind + 1]) + + indices = A.indices[ind_slice] + indptr[1:] = np.cumsum(A.indptr[ind + 1] - A.indptr[ind]) + data = A.data[ind_slice] + + if A.getformat() == 'csc': + return sps.csc_matrix((data, indices, indptr), shape=(A.shape[0], N)) + elif A.getformat() == 'csr': + return sps.csr_matrix((data, indices, indptr), shape=(N, A.shape[1])) diff --git a/test/unit/test_sparse_mat.py b/test/unit/test_sparse_mat.py index 3bcde967ca..acc0eecbae 100644 --- a/test/unit/test_sparse_mat.py +++ b/test/unit/test_sparse_mat.py @@ -76,5 +76,60 @@ def test_zero_columns(self): assert np.sum(A2 != A2_t) == 0 assert np.sum(A0_2 != A0_2_t) == 0 + #------------------ Test sliced_mat() ----------------------- + def test_sliced_mat_columns(self): + # Test slicing of csr_matrix + A = sps.csc_matrix(np.array([[0, 0, 0], + [1, 0, 0], + [0, 0, 3]])) + + A0_t = sps.csc_matrix(np.array([[0, 0], + [0, 0], + [0, 3]])) + A1_t = sps.csc_matrix(np.array([[0, 0], + [1, 0], + [0, 0]])) + A2_t = sps.csc_matrix(np.array([[0], + [0], + [3]])) + A3_t = sps.csc_matrix(np.array([[], + [], + []])) + + A0 = sparse_mat.slice_mat(A, np.array([1, 2], dtype=int)) + A1 = sparse_mat.slice_mat(A, np.array([0, 1])) + A2 = sparse_mat.slice_mat(A, 2) + A3 = sparse_mat.slice_mat(A, np.array([], dtype=np.int)) + import pdb + pdb.set_trace() + + assert np.sum(A0 != A0_t) == 0 + assert np.sum(A1 != A1_t) == 0 + assert np.sum(A2 != A2_t) == 0 + assert np.sum(A3 != A3_t) == 0 + + def test_sliced_mat_columns(self): + # Test slicing of csr_matrix + A = sps.csr_matrix(np.array([[0, 0, 0], + [1, 0, 0], + [0, 0, 3]])) + + A0_t = sps.csr_matrix(np.array([[1, 0, 0], + [0, 0, 3]])) + A1_t = sps.csr_matrix(np.array([[0, 0, 0], + [1, 0, 0]])) + A2_t = sps.csr_matrix(np.array([[0, 0, 3]])) + A3_t = sps.csr_matrix(np.atleast_2d(np.array([[], [], []])).T) + + A0 = sparse_mat.slice_mat(A, np.array([1, 2], dtype=int)) + A1 = sparse_mat.slice_mat(A, np.array([0, 1])) + A2 = sparse_mat.slice_mat(A, 2) + A3 = sparse_mat.slice_mat(A, np.array([], dtype=np.int)) + + assert np.sum(A0 != A0_t) == 0 + assert np.sum(A1 != A1_t) == 0 + assert np.sum(A2 != A2_t) == 0 + assert np.sum(A3 != A3_t) == 0 + if __name__ == '__main__': unittest.main() From e888b93ba3411b874e22088d2d9e0e1055f9db75 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 09:52:20 +0100 Subject: [PATCH 102/118] Improved T-intersection finder in Fracture extrusion. Case with different family names are now handled properly. --- src/porepy/fracs/extrusion.py | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index d1d0220df9..9301304cfa 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -116,7 +116,7 @@ def _intersection_by_num_node(edges, num): return crosses, edges_of_crosses -def t_intersections(edges): +def t_intersections(edges, remove_three_families=True): """ Find points involved in T-intersections. A t-intersection is defined as a point involved in three fracture segments, @@ -129,6 +129,9 @@ def t_intersections(edges): edges (np.array, 3 x n): Fractures. First two rows give indices of start and endpoints. Last row gives index of fracture that the segment belongs to. + remove_three_families (boolean, defaults to True): If True, + T-intersections with where fractures from three different families + meet will be removed. Returns: abutments (np.ndarray, int): indices of points that are @@ -144,16 +147,24 @@ def t_intersections(edges): frac_num = edges[-1] abutments, edges_of_abutments = _intersection_by_num_node(edges, 3) + # Initialize fields for abutments. num_abut = abutments.size primal_frac = np.zeros(num_abut, dtype=np.int) sec_frac = np.zeros(num_abut, dtype=np.int) other_point = np.zeros(num_abut, dtype=np.int) + + # If fractures meeting in a T-intersection all have different family names, + # these will be removed from the list. + remove = np.zeros(num_abut, dtype=np.bool) + for i, (pi, ei) in enumerate(zip(abutments, edges_of_abutments)): # Count number of occurences for each fracture associated with this # intersection. fi_all = frac_num[ei] fi, count = np.unique(fi_all, return_counts=True) - assert fi.size == 2 + if fi.size > 2: + remove[i] = 1 + continue # Find the fracture number associated with main and abutting edge. if count[0] == 1: primal_frac[i] = fi[1] @@ -170,6 +181,14 @@ def t_intersections(edges): else: other_point[i] = edges[0, ei_abut] + # Remove any T-intersections that did not belong to + if remove_three_families and remove.any(): + remove = np.where(remove)[0] + abutments = np.delete(abutments, remove) + primal_frac = np.delete(primal_frac, remove) + sec_frac = np.delete(sec_frac, remove) + other_point = np.delete(other_point, remove) + return abutments, primal_frac, sec_frac, other_point From 2ccb503cb6d81af887c3490b959cd865026979ba Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 1 Nov 2017 09:53:13 +0100 Subject: [PATCH 103/118] Finished optimizing split_nodes --- src/porepy/fracs/meshing.py | 4 ++-- src/porepy/fracs/split_grid.py | 25 ++++++++++++++++--------- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/src/porepy/fracs/meshing.py b/src/porepy/fracs/meshing.py index f9281317b1..bf88050014 100644 --- a/src/porepy/fracs/meshing.py +++ b/src/porepy/fracs/meshing.py @@ -352,7 +352,7 @@ def assemble_in_bucket(grids, **kwargs): cell_2_face, cell = obtain_interdim_mappings( lg, fn, n_per_face, **kwargs) face_cells = sps.csc_matrix( - (np.array([True] * cell.size), (cell, cell_2_face)), + (np.ones(cell.size, dtype=bool), (cell, cell_2_face)), (lg.num_cells, hg.num_faces)) # This if may be unnecessary, but better safe than sorry. @@ -402,6 +402,6 @@ def obtain_interdim_mappings(lg, fn, n_per_face, dim grid. This likely is related to gmsh behavior. ''') else: warnings.warn('''Found inconsistency between cells and higher - dimensional faces. Continuing, faces crossed''') + dimensional faces. Continuing, fingers crossed''') low_dim_cell = np.where(is_mem)[0] return cell_2_face, low_dim_cell diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py index 8856cd706c..2032d1fc93 100644 --- a/src/porepy/fracs/split_grid.py +++ b/src/porepy/fracs/split_grid.py @@ -7,7 +7,7 @@ import warnings from porepy.utils.half_space import half_space_int -from porepy.utils import sparse_mat +from porepy.utils import sparse_mat, setmembership from porepy.utils.graph import Graph from porepy.utils.mcolon import mcolon from porepy.grids.grid import Grid, FaceTag @@ -254,6 +254,7 @@ def update_face_cells(face_cells, face_id, i): # For the other lower-dim grids we just add zeros to conserve # the right matrix dimensions. for j, f_c in enumerate(face_cells): + assert f_c.getformat() == 'csc' if j == i: f_c = sps.hstack((f_c, f_c[:, face_id])) else: @@ -415,21 +416,27 @@ def duplicate_nodes(g, nodes, offset): new_nodes = np.repeat(g.nodes[:, t_node, None], colors.size, axis=1) faces = np.array([], dtype=int) face_pos = np.array([g.face_nodes.indptr[t_node]]) + assert g.cell_faces.getformat() == 'csc' + assert g.face_nodes.getformat() == 'csr' + faces_of_node_t = sparse_mat.slice_indices(g.face_nodes, t_node) for j in range(colors.size): # For each color we wish to add one node. First we find all faces that - # are connected to the fracture node, and have the correct cell + # are connected to the fracture node and have the correct cell # color - local_faces = (g.cell_faces[:, cells[ix == j]]).nonzero()[0] - local_faces = np.unique(local_faces) - con_to_node = np.ravel(g.face_nodes[t_node, local_faces].todense()) - faces = np.append(faces, local_faces[con_to_node]) - # These faces is then attached to new node number j. - face_pos = np.append(face_pos, face_pos[-1] + np.sum(con_to_node)) + colored_faces = sparse_mat.slice_indices( + g.cell_faces, cells[ix == j]) + colored_faces = np.unique(colored_faces) + is_colored = np.in1d( + faces_of_node_t, colored_faces, assume_unique=True) + + faces = np.append(faces, faces_of_node_t[is_colored]) + # These faces are then attached to new node number j. + face_pos = np.append(face_pos, face_pos[-1] + np.sum(is_colored)) # If an offset is given, we will change the position of the nodes. # We move the nodes a length of offset away from the fracture(s). if offset > 0 and colors.size > 1: new_nodes[:, j] -= avg_normal(g, - local_faces[con_to_node]) * offset + faces_of_node_t[is_colored]) * offset # The total number of faces should not have changed, only their # connection to nodes. We can therefore just update the indices and # indptr map. From 9e6e54668ae121a29b30e55727f4d3361aa26fad Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 12:10:34 +0100 Subject: [PATCH 104/118] Fix of Fracture rotation in extrusion module. Rotate around exposed point, rather than the fracture center. The new version will preserve the original outcrop of a fracture set. --- src/porepy/fracs/extrusion.py | 39 ++++++++++++++++++++++++----------- 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 9301304cfa..486aabe7c9 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -62,7 +62,8 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, * exposure = p1 - p0 # Impose incline. - rot_ang = impose_inlcine(fractures, exposure, frac_family=family, **kwargs) + rot_ang = impose_inlcine(fractures, exposure, p0, frac_family=family, + **kwargs) # Cut fractures for prim, sec, p in zip(prim_frac, sec_frac, other_pt): @@ -86,7 +87,7 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, * strike = np.arctan2(e1[1] - e0[1], e1[0]- e0[0]) f = create_fracture(center, new_radius, np.pi/2, strike, np.vstack((e0, e1)).T) - rotate_fracture(f, e1-e0, rot_ang[prim]) + rotate_fracture(f, e1-e0, rot_ang[prim], p0[:, prim]) fractures[prim] = f return fractures @@ -383,28 +384,40 @@ def create_fracture(center, radius, dip, strike, extra_points): return f -def rotate_fracture(frac, vec, angle): - """ Rotate a fracture along a specified strike vector. +def rotate_fracture(frac, vec, angle, exposure): + """ Rotate a fracture along a specified strike vector, and centered on a + given point on the fracture surface. + + Modification of the fracture coordinates is done in place. + + TODO: Move this to the fracture itself? Parameters: frac (Fracture): To be rotated. Points are modified in-place. - vec (np.array-like): Rotation will be around this angle. + vec (np.array-like): Rotation will be around this vector. ang (double). Rotation angle. Measured in radians. + exposure (np.array-like): Point on the strike vector, rotation will be + centered around the line running through this point. """ + vec = np.asarray(vec) + exposure = np.asarray(exposure) if vec.size == 2: vec = np.append(vec, 0) + if exposure.size == 2: + exposure = np.append(exposure, 0) + exposure = exposure.reshape((3, 1)) rot = cg.rot(angle, vec) p = frac.p - center = np.mean(p, axis=1).reshape((-1, 1)) - frac.p = center + rot.dot(p - center) + frac.p = exposure + rot.dot(p - exposure) frac.points_2_ccw() -def impose_inlcine(fracs, exposure, frac_family=None, family_mean_incline=None, - family_std_incline=None, **kwargs): +def impose_inlcine(fracs, exposure_line, exposure_point, frac_family=None, + family_mean_incline=None, family_std_incline=None, + **kwargs): """ Impose incline on the fractures from family-based parameters. The incline for each family is specified in terms of its mean and standard @@ -414,8 +427,10 @@ def impose_inlcine(fracs, exposure, frac_family=None, family_mean_incline=None, Parameters: fracs (list of Frature): Fractures to be inclined. - exposure (np.array, 3xnum_frac): Exposed line for each fracture + exposure_line (np.array, 3xnum_frac): Exposed line for each fracture (visible in outcrop). Rotation will be around this line. + exposure_point (np.array, 3xnum_frac): Point on exposure line. This + point will not be rotated, it's a fixture. family (np.array, num_fracs): For each fracture, which family does it belong to. If not provided, all fractures are considered to belong to the same family. @@ -438,13 +453,13 @@ def impose_inlcine(fracs, exposure, frac_family=None, family_mean_incline=None, if family_std_incline is None: family_std_incline = np.zeros(np.unique(frac_family).size) - exposure = np.vstack((exposure, np.zeros(len(fracs)))) + exposure_line = np.vstack((exposure_line, np.zeros(len(fracs)))) all_ang = np.zeros(len(fracs)) for fi, f in enumerate(fracs): fam = frac_family[fi] ang = np.random.normal(loc=family_mean_incline[fam], scale=family_std_incline[fam]) - rotate_fracture(f, exposure[:, fi], ang) + rotate_fracture(f, exposure_line[:, fi], ang, exposure_point[:, fi]) all_ang[fam] = ang return all_ang From 68387317723d87d903ff26d5bfcdc739e8785f72 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 12:16:39 +0100 Subject: [PATCH 105/118] Cut of fractures in extrusion handles case of non-intersecting fractures. --- src/porepy/fracs/extrusion.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 486aabe7c9..0801a0e9df 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -555,6 +555,10 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, isect_pt, _, _ = main_frac.intersects(aux_frac, tol) + if isect_pt.size == 0: + print('No intersection found in cutting of fractures') + return main_frac + # Next step is to eliminate points in the main fracture that are on the # wrong side of the other fracture. v = main_frac.p - other_frac.center.reshape((-1, 1)) From 5a99aaedd60b3c8f838e1b1700dd9790229a7553 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 15:52:16 +0100 Subject: [PATCH 106/118] Speed up unit tests involving fractures. No need to verify convexity in these cases. --- test/unit/test_fracture_center.py | 20 +++++++++---------- test/unit/test_fracture_intersect_boundary.py | 4 ++-- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/test/unit/test_fracture_center.py b/test/unit/test_fracture_center.py index 149afdf006..8cf0451325 100644 --- a/test/unit/test_fracture_center.py +++ b/test/unit/test_fracture_center.py @@ -9,31 +9,31 @@ class TestFractureCenters(unittest.TestCase): def test_frac_1(self): f_1 = Fracture(np.array([[0, 2, 2, 0], [0, 2, 2, 0], - [-1, -1, 1, 1]])) - c_known = np.array([1, 1, 0]) + [-1, -1, 1, 1]]), check_convexity=False) + c_known = np.array([1, 1, 0]).reshape((3, 1)) assert np.allclose(c_known, f_1.center) def test_frac_2(self): f_1 = Fracture(np.array([[0, 2, 2, 0], [0, 2, 2, 0], - [0, 0, 1, 1]])) - c_known = np.array([1, 1, 0.5]) + [0, 0, 1, 1]]), check_convexity=False) + c_known = np.array([1, 1, 0.5]).reshape((3, 1)) assert np.allclose(c_known, f_1.center) def test_frac_3(self): # Fracture plane defined by x + y + z = 1 f_1 = Fracture(np.array([[0, 1, 1, 0], [0, 0, 1, 1], - [1, 0, -1, 0]])) - c_known = np.array([0.5, 0.5, 0]) + [1, 0, -1, 0]]), check_convexity=False) + c_known = np.array([0.5, 0.5, 0]).reshape((3, 1)) assert np.allclose(c_known, f_1.center) def test_frac_4(self): # Fracture plane defined by x + y + z = 4 f_1 = Fracture(np.array([[0, 1, 1, 0], [0, 0, 1, 1], - [4, 3, 2, 3]])) - c_known = np.array([0.5, 0.5, 3]) + [4, 3, 2, 3]]), check_convexity=False) + c_known = np.array([0.5, 0.5, 3]).reshape((3, 1)) assert np.allclose(c_known, f_1.center) def test_frac_rand(self): @@ -41,8 +41,8 @@ def test_frac_rand(self): x = np.array([0, 1, 1, 0]) y = np.array([0, 0, 1, 1]) z = (r[0] - r[1] * x - r[2] * y) / r[3] - f = Fracture(np.vstack((x, y, z))) + f = Fracture(np.vstack((x, y, z)), check_convexity=False) z_cc = (r[0] - 0.5 * r[1] - 0.5 * r[2]) / r[3] - c_known = np.array([0.5, 0.5, z_cc]) + c_known = np.array([0.5, 0.5, z_cc]).reshape((3, 1)) assert np.allclose(c_known, f.center) diff --git a/test/unit/test_fracture_intersect_boundary.py b/test/unit/test_fracture_intersect_boundary.py index 4b7dc122ea..55267b6e5d 100644 --- a/test/unit/test_fracture_intersect_boundary.py +++ b/test/unit/test_fracture_intersect_boundary.py @@ -16,7 +16,7 @@ class TestFractureBoundaryIntersection(unittest.TestCase): def setup(self): self.f_1 = Fracture(np.array([[0, 1, 1, 0], [.5, .5, .5, .5], - [0, 0, 1, 1]])) + [0, 0, 1, 1]]), check_convexity=False) self.domain = {'xmin': 0, 'xmax': 1, 'ymin': 0, 'ymax': 1, 'zmin': 0, 'zmax': 1} @@ -96,7 +96,7 @@ def test_full_incline(self): p = np.array([[-0.5, 0.5, 0.5, -0.5], [0.5, 0.5, 1.5, 1.5], [-0.5, -0.5, 1, 1]]) - f = Fracture(p) + f = Fracture(p, check_convexity=False) network = FractureNetwork([f]) network.impose_external_boundary(self.domain) p_known = np.array([[0., 0.5, 0.5, 0], From 6fe5d248807f5197cb249dcc59bd73997abe0618 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 16:02:44 +0100 Subject: [PATCH 107/118] Unit tests for fracture extrusion. Not complete set of tests yet. --- test/unit/test_fracture_extrusion.py | 135 +++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 test/unit/test_fracture_extrusion.py diff --git a/test/unit/test_fracture_extrusion.py b/test/unit/test_fracture_extrusion.py new file mode 100644 index 0000000000..d06bc4ffe4 --- /dev/null +++ b/test/unit/test_fracture_extrusion.py @@ -0,0 +1,135 @@ +import unittest +import numpy as np + +from porepy.fracs.fractures import Fracture +from porepy.fracs import extrusion + +class TestFractureExtrusion(unittest.TestCase): + + def test_single_t_intersection(self): + # Three fractures meeting in a T, same family tag for two of them + edges = np.array([[0, 1, 0], [1, 2, 0], [1, 3, 1]]).T + abut, prim, sec, other = extrusion.t_intersections(edges) + + known_abut = [1] + known_prim = [0] + known_sec = [1] + known_other = [3] + + assert np.all(known_abut == abut) + assert np.all(known_prim == prim) + assert np.all(known_sec == sec) + assert np.all(known_other == other) + + def test_single_t_intersection_different_family(self): + # Three fractures meeting in a T. Different family on all, so no + # intersections found + edges = np.array([[0, 1, 0], [1, 2, 1], [1, 3, 2]]).T + abut, prim, sec, other = extrusion.t_intersections(edges) + + assert abut.size == 0 + assert prim.size == 0 + assert sec.size == 0 + assert other.size == 0 + + def test_single_t_intersection_disregard_family_tag(self): + # Three fractures meeting in a T. Different family on all, but we + # disregard family tags, so a T-intersection is found + edges = np.array([[0, 1, 0], [1, 2, 1], [1, 3, 2]]).T + abut, prim, sec, other = \ + extrusion.t_intersections(edges, remove_three_families=True) + known_abut = [1] + known_prim = [0] + known_sec = [1] + known_other = [3] + + assert np.all(known_abut == abut) + assert np.all(known_prim == prim) + assert np.all(known_sec == sec) + assert np.all(known_other == other) + + def test_disc_radius_center(self): + + p0 = np.array([[0], [0]]) + p1 = np.array([[0], [1]]) + length = np.array([1]) + theta = np.array([np.pi/2]) + + known_cc = np.array([[0], [0.5], [0]]) + + rad, cc, _ = extrusion.disc_radius_center(length, p0, p1, theta) + assert rad[0] == 0.5 + assert np.allclose(cc, known_cc) + + def compare_arrays(self, a, b): + assert a.shape == b.shape + + for i in range(a.shape[1]): + ai = a[:, i].reshape((3, 1)) + assert np.min(np.sum((b - ai)**2, axis=0)) < 1e-5 + bi = b[:, i].reshape((3, 1)) + assert np.min(np.sum((a - bi)**2, axis=0)) < 1e-5 + + def test_fracture_rotation_90_deg(self): + + p = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + f = Fracture(p, check_convexity=False) + + p_rot_1_known = np.array([[0, 0, 0], [1, 0, 0], [1, 0, 1], [0, 0, 1]]).T + exposure = np.array([[0], [0], [0]]) + extrusion.rotate_fracture(f, [1, 0, 0], np.pi/2, exposure) + self.compare_arrays(f.p, p_rot_1_known) + + p_rot_2_known = np.array([[0, 0.5, -0.5], [1, 0.5, -0.5], + [1, 0.5, 0.5], [0, 0.5, 0.5]]).T + exposure = np.array([[0], [0.5], [0]]) + f = Fracture(p, check_convexity=False) + extrusion.rotate_fracture(f, [1, 0, 0], np.pi/2, exposure) + self.compare_arrays(f.p, p_rot_2_known) + + def test_cut_fracture_no_intersection(self): + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]).T + f1 = Fracture(p1, check_convexity=False) + + p2 = np.array([[0.5, -0.5, 0.2], [0.5, 0.5, 0.2], [0.5, 0.5, 0.8], + [0.5, -0.5, 0.8]]).T + p_known = p2.copy() + f2 = Fracture(p2, check_convexity=False) + + ref_pt = np.array([[0], [1], [0]]) + extrusion.cut_fracture_by_plane(f2, f1, ref_pt) + + self.compare_arrays(p_known, f2.p) + + def test_cut_fracture_simple_intersection(self): + p1 = np.array([[0, 0, 0], [1, 0, 0], [1, 0, 1], [0, 0, 1]]).T + f1 = Fracture(p1, check_convexity=False) + + p2 = np.array([[0.5, -0.5, 0.2], [0.5, 0.5, 0.2], [0.5, 0.5, 0.8], + [0.5, -0.5, 0.8]]).T + p_known = np.array([[0.5, 0, 0.2], [0.5, 0.5, 0.2], [0.5, 0.5, 0.8], + [0.5, 0, 0.8]]).T + f2 = Fracture(p2, check_convexity=False) + + ref_pt = np.array([[0], [1], [0]]) + extrusion.cut_fracture_by_plane(f2, f1, ref_pt) + + self.compare_arrays(p_known, f2.p) + + def test_cut_fracture_one_inclined(self): + p1 = np.array([[0, 1, -0.5], [1, 1, -0.5], [1, -1, 1.5], [0, -1, 1.5]]).T + f1 = Fracture(p1, check_convexity=False) + + p2 = np.array([[0.5, -0.5, 0.2], [0.5, 0.5, 0.2], [0.5, 0.5, 0.8], + [0.5, -1.5, 0.8]]).T + p_known = np.array([[0.5, 0.3, 0.2], [0.5, 0.5, 0.2], [0.5, 0.5, 0.8], + [0.5, -0.3, 0.8]]).T + f2 = Fracture(p2, check_convexity=False) + + ref_pt = np.array([[0], [1], [0]]) + extrusion.cut_fracture_by_plane(f2, f1, ref_pt) + + self.compare_arrays(p_known, f2.p) + + if __name__ == '__main__': + unittest.main() From cf97b5c4cf57e5141ef33307000b49f790d967a6 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 16:03:59 +0100 Subject: [PATCH 108/118] Fix of fracture center computation. --- src/porepy/fracs/fractures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 085b1c4bd0..2feeb370af 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -262,7 +262,7 @@ def compute_centroid(self): center = np.sum(cc * area, axis=1) / np.sum(area) # Project back again. - self.center = rot.transpose().dot(np.append(center, z)) + self.center = rot.transpose().dot(np.append(center, z)).reshape((3, 1)) def as_sp_polygon(self, p=None): """ Represent polygon as a sympy object. From 36194c226a801d057d20186c551fbd593d860b66 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 16:05:47 +0100 Subject: [PATCH 109/118] Fix of __str__ of Fracture. No longer relient on sympy. --- src/porepy/fracs/fractures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 2feeb370af..4fa43bee8d 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -866,7 +866,7 @@ def outside_box(p, bound_i): break def __repr__(self): - return str(self.as_sp_polygon()) + return self.__str__() def __str__(self): s = 'Points: \n' From 60ce8804d3c771ccedd69378035998b319a0d2e1 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 16:06:26 +0100 Subject: [PATCH 110/118] Minor bugfixes in fracture extrusion. Returned wrong rotation angles. Formating. --- src/porepy/fracs/extrusion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index 0801a0e9df..cf6131ba4b 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -68,7 +68,7 @@ def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, * # Cut fractures for prim, sec, p in zip(prim_frac, sec_frac, other_pt): _, radius = cut_fracture_by_plane(fractures[sec], fractures[prim], - split_pt[:, p], **kwargs) + split_pt[:, p], **kwargs) # If specified, ensure that cuts in T-intersections appear realistic. if ensure_realistic_cuts and radius is not None: ang = np.arctan2(0.5*lengths[prim], radius) @@ -460,7 +460,7 @@ def impose_inlcine(fracs, exposure_line, exposure_point, frac_family=None, ang = np.random.normal(loc=family_mean_incline[fam], scale=family_std_incline[fam]) rotate_fracture(f, exposure_line[:, fi], ang, exposure_point[:, fi]) - all_ang[fam] = ang + all_ang[fi] = ang return all_ang From 303ab9812acbfe11325a3860ad29790ea2dbfaa9 Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 1 Nov 2017 16:59:29 +0100 Subject: [PATCH 111/118] Improved update_face_cells --- src/porepy/fracs/split_grid.py | 36 ++++++++-- src/porepy/utils/sparse_mat.py | 116 +++++++++++++++++++++++++++++++++ test/unit/test_sparse_mat.py | 76 ++++++++++++++++++++- 3 files changed, 219 insertions(+), 9 deletions(-) diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py index 2032d1fc93..db4264a4e2 100644 --- a/src/porepy/fracs/split_grid.py +++ b/src/porepy/fracs/split_grid.py @@ -172,7 +172,7 @@ def split_nodes(gh, gl, gh_2_gl_nodes, offset=0): node_count = duplicate_nodes(gh, nodes, offset) # We remove the old nodes. - #gh = remove_nodes(gh, nodes) + # gh = remove_nodes(gh, nodes) # Update the number of nodes gh.num_nodes = gh.num_nodes + node_count # - nodes.size @@ -224,9 +224,9 @@ def duplicate_faces(gh, face_cells): node_start = gh.face_nodes.indptr[frac_id] node_end = gh.face_nodes.indptr[frac_id + 1] - #frac_nodes = gh.face_nodes[:, frac_id] + # frac_nodes = gh.face_nodes[:, frac_id] - #gh.face_nodes = sps.hstack((gh.face_nodes, frac_nodes)) + # gh.face_nodes = sps.hstack((gh.face_nodes, frac_nodes)) # We also copy the attributes of the original faces. gh.num_faces += frac_id.size gh.face_normals = np.hstack( @@ -253,13 +253,29 @@ def update_face_cells(face_cells, face_id, i): # The duplications should also be associated with grid i. # For the other lower-dim grids we just add zeros to conserve # the right matrix dimensions. + if face_id.size == 0: + return face_cells + for j, f_c in enumerate(face_cells): assert f_c.getformat() == 'csc' if j == i: - f_c = sps.hstack((f_c, f_c[:, face_id])) + f_c_sliced = sparse_mat.slice_mat(f_c, face_id) + new_indptr = f_c_sliced.indptr + f_c.indptr[-1] + new_ind = f_c_sliced.indices + new_data = f_c_sliced.data + + f_c.indptr = np.append(f_c.indptr, new_indptr[1:]) + f_c.indices = np.append(f_c.indices, new_ind) + f_c.data = np.append(f_c.data, new_data) + f_c._shape = (f_c._shape[0], f_c._shape[1] + face_id.size) + #f_c = sps.hstack((f_c, f_c[:, face_id])) else: - empty = sps.csc_matrix((f_c.shape[0], face_id.size)) - f_c = sps.hstack((f_c, empty)) + new_indptr = f_c.indptr[-1] * \ + np.ones(face_id.size, dtype=f_c.indptr.dtype) + f_c.indptr = np.append(f_c.indptr, new_indptr) + f_c._shape = (f_c._shape[0], f_c._shape[1] + face_id.size) +# empty = sps.csc_matrix((f_c.shape[0], face_id.size)) +# f_c = sps.hstack((f_c, empty)) face_cells[j] = f_c return face_cells @@ -284,6 +300,7 @@ def update_cell_connectivity(g, face_id, normal, x0): """ # We find the cells attached to the tagged faces. + # g.cell_faces = g.cell_faces.tocsr() cell_frac = g.cell_faces[face_id, :] cell_face_id = np.argwhere(cell_frac) @@ -324,12 +341,16 @@ def update_cell_connectivity(g, face_id, normal, x0): data = np.ravel(g.cell_faces[np.ravel(face_id[row]), col]) cell_frac_right = sps.csc_matrix((data, (row, col)), (face_id.size, g.cell_faces.shape[1])) - g.cell_faces[face_id, :] = cell_frac_right + # assert g.cell_faces.getformat() == 'csr' + # sparse_mat.merge_matrices(g.cell_faces, cell_frac_right, face_id) + g.cell_faces[face_id, :] = cell_frac_right # And then we add the new left-faces to the cell_face map. We do not # change the sign of the matrix since we did not flip the normals. # This means that the normals of right and left cells point in the same # direction, but their cell_faces values have oposite signs. + # sparse_mat.stack_mat(g.cell_faces, cell_frac_left) + # g.cell_faces.tocsc() g.cell_faces = sps.vstack((g.cell_faces, cell_frac_left), format='csc') return 0 @@ -491,6 +512,7 @@ def find_cell_color(g, cells): """ c = np.sort(cells) # Local cell-face and face-node maps. + assert g.cell_faces.getformat() == 'csc' cf_sub = sparse_mat.slice_mat(g.cell_faces, c) child_cell_ind = -np.ones(g.num_cells, dtype=np.int) child_cell_ind[c] = np.arange(cf_sub.shape[1]) diff --git a/src/porepy/utils/sparse_mat.py b/src/porepy/utils/sparse_mat.py index 028113e596..a7a4e6e6e1 100644 --- a/src/porepy/utils/sparse_mat.py +++ b/src/porepy/utils/sparse_mat.py @@ -32,6 +32,122 @@ def zero_columns(A, cols): A.data[col_indptr] = 0 +def merge_matrices(A, B, lines): + ''' + Replace rows/coloms of matrix A with rows/cols of matrix B. + If A and B are csc matrices this function is equivalent with + A[:, lines] = B + If A and B are csr matrices this funciton is equivalent iwth + A[lines, :] = B + + Parameter + --------- + A (scipy.sparse.spmatrix): A sparce matrix + B (scipy.sparse.spmatrix): A sparce matrix + lines (ndarray): Lines of A to be replaced by B. + + Return + ------ + None + + + ''' + if A.getformat() != 'csc' and A.getformat() != 'csr': + raise ValueError('Need a csc or csr matrix') + elif A.getformat() != B.getformat(): + raise ValueError('A and B must be of same matrix type') + if A.getformat() == 'csc': + if A.shape[0] != B.shape[0]: + raise ValueError('A.shape[0] must equal B.shape[0]') + if A.getformat() == 'csr': + if A.shape[1] != B.shape[1]: + raise ValueError('A.shape[0] must equal B.shape[0]') + + if B.getformat() == 'csc': + if lines.size != B.shape[1]: + raise ValueError('B.shape[1] must equal size of lines') + if B.getformat() == 'csr': + if lines.size != B.shape[0]: + raise ValueError('B.shape[0] must equal size of lines') + + indptr = A.indptr + indices = A.indices + data = A.data + + ind_ix = mcolon(indptr[lines], indptr[lines + 1]) + + # First we remove the old data + num_rem = np.zeros(indptr.size, dtype=np.int32) + num_rem[lines + 1] = indptr[lines + 1] - indptr[lines] + num_rem = np.cumsum(num_rem, dtype=num_rem.dtype) + + indptr = indptr - num_rem + + keep = np.ones(A.data.size, dtype=np.bool) + keep[ind_ix] = False + indices = indices[keep] + data = data[keep] + + # Then we add the new + b_indptr = B.indptr + b_indices = B.indices + b_data = B.data + + num_added = np.zeros(indptr.size, dtype=np.int32) + num_added[lines + 1] = b_indptr[1:] - b_indptr[:-1] + num_added = np.cumsum(num_added, dtype=num_added.dtype) + + rep = np.diff(b_indptr) + indPos = np.repeat(indptr[lines], rep) + + A.indices = np.insert(indices, indPos, b_indices) + A.data = np.insert(data, indPos, b_data) + A.indptr = indptr + num_added + + +def stack_mat(A, B): + ''' + Stack matrix B at the end of matrix A. + If A and B are csc matrices this function is equivalent to + A = scipy.sparse.hstack((A, B)) + If A and B are csr matrices this function is equivalent to + A = scipy.sparse.vstack((A, B)) + + Parameters: + ----------- + A (scipy.sparse.spmatrix): A sparce matrix + B (scipy.sparse.spmatrix): A sparce matrix + + Return + ------ + None + + + ''' + if A.getformat() != 'csc' and A.getformat() != 'csr': + raise ValueError('Need a csc or csr matrix') + elif A.getformat() != B.getformat(): + raise ValueError('A and B must be of same matrix type') + if A.getformat() == 'csc': + if A.shape[0] != B.shape[0]: + raise ValueError('A.shape[0] must equal B.shape[0]') + if A.getformat() == 'csr': + if A.shape[1] != B.shape[1]: + raise ValueError('A.shape[0] must equal B.shape[0]') + + if B.indptr.size == 1: + return + + A.indptr = np.append(A.indptr, B.indptr[1:] + A.indptr[-1]) + A.indices = np.append(A.indices, B.indices) + A.data = np.append(A.data, B.data) + + if A.getformat() == 'csc': + A._shape = (A._shape[0], A._shape[1] + B._shape[1]) + if A.getformat() == 'csr': + A._shape = (A._shape[0] + B._shape[0], A._shape[1]) + + def slice_indices(A, slice_ind): """ Function for slicing sparse matrix along rows or columns. diff --git a/test/unit/test_sparse_mat.py b/test/unit/test_sparse_mat.py index acc0eecbae..8a868c1c4a 100644 --- a/test/unit/test_sparse_mat.py +++ b/test/unit/test_sparse_mat.py @@ -100,8 +100,6 @@ def test_sliced_mat_columns(self): A1 = sparse_mat.slice_mat(A, np.array([0, 1])) A2 = sparse_mat.slice_mat(A, 2) A3 = sparse_mat.slice_mat(A, np.array([], dtype=np.int)) - import pdb - pdb.set_trace() assert np.sum(A0 != A0_t) == 0 assert np.sum(A1 != A1_t) == 0 @@ -131,5 +129,79 @@ def test_sliced_mat_columns(self): assert np.sum(A2 != A2_t) == 0 assert np.sum(A3 != A3_t) == 0 + #------------------ Test stack_mat() ----------------------- + def test_stack_mat_columns(self): + # Test slicing of csr_matrix + A = sps.csc_matrix(np.array([[0, 0, 0], + [1, 0, 0], + [0, 0, 3]])) + + B = sps.csc_matrix(np.array([[0, 2], + [3, 1], + [1, 0]])) + + A_t = sps.csc_matrix(np.array([[0, 0, 0, 0, 2], + [1, 0, 0, 3, 1], + [0, 0, 3, 1, 0]])) + + sparse_mat.stack_mat(A, B) + + assert np.sum(A != A_t) == 0 + + def test_stack_mat_rows(self): + # Test slicing of csr_matrix + A = sps.csr_matrix(np.array([[0, 0, 0], + [1, 0, 0], + [0, 0, 3]])) + + B = sps.csr_matrix(np.array([[0, 2, 2], + [3, 1, 3]])) + + A_t = sps.csr_matrix(np.array([[0, 0, 0], + [1, 0, 0], + [0, 0, 3], + [0, 2, 2], + [3, 1, 3]])) + + sparse_mat.stack_mat(A, B) + + assert np.sum(A != A_t) == 0 + #------------------ Test merge_mat() ----------------------- + + def test_merge_mat_columns(self): + # Test slicing of csr_matrix + A = sps.csc_matrix(np.array([[0, 0, 0], + [1, 0, 0], + [0, 0, 3]])) + + B = sps.csc_matrix(np.array([[0, 2], + [3, 1], + [1, 0]])) + + A_t = sps.csc_matrix(np.array([[0, 0, 2], + [1, 3, 1], + [0, 1, 0]])) + + sparse_mat.merge_matrices(A, B, np.array([1, 2], dtype=np.int)) + + assert np.sum(A != A_t) == 0 + + def test_merge_mat_rows(self): + # Test slicing of csr_matrix + A = sps.csr_matrix(np.array([[0, 0, 0], + [1, 0, 0], + [0, 0, 3]])) + + B = sps.csr_matrix(np.array([[0, 2, 2], + [3, 1, 3]])) + + A_t = sps.csr_matrix(np.array([[0, 2, 2], + [3, 1, 3], + [0, 0, 3]])) + + sparse_mat.merge_matrices(A, B, np.array([0, 1], dtype=np.int)) + + assert np.sum(A != A_t) == 0 + if __name__ == '__main__': unittest.main() From d34fff4ef25bcbfc7281f226ff0a4e64284c478a Mon Sep 17 00:00:00 2001 From: Runar Date: Wed, 1 Nov 2017 17:01:21 +0100 Subject: [PATCH 112/118] Started improving update_cell_connectivity --- src/porepy/fracs/split_grid.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py index db4264a4e2..8ef786f047 100644 --- a/src/porepy/fracs/split_grid.py +++ b/src/porepy/fracs/split_grid.py @@ -300,7 +300,7 @@ def update_cell_connectivity(g, face_id, normal, x0): """ # We find the cells attached to the tagged faces. - # g.cell_faces = g.cell_faces.tocsr() + g.cell_faces = g.cell_faces.tocsr() cell_frac = g.cell_faces[face_id, :] cell_face_id = np.argwhere(cell_frac) @@ -339,18 +339,18 @@ def update_cell_connectivity(g, face_id, normal, x0): col = cell_face_id[~left_cell, 1] row = cell_face_id[~left_cell, 0] data = np.ravel(g.cell_faces[np.ravel(face_id[row]), col]) - cell_frac_right = sps.csc_matrix((data, (row, col)), + cell_frac_right = sps.csr_matrix((data, (row, col)), (face_id.size, g.cell_faces.shape[1])) - # assert g.cell_faces.getformat() == 'csr' - # sparse_mat.merge_matrices(g.cell_faces, cell_frac_right, face_id) - g.cell_faces[face_id, :] = cell_frac_right + assert g.cell_faces.getformat() == 'csr' + sparse_mat.merge_matrices(g.cell_faces, cell_frac_right, face_id) + # g.cell_faces[face_id, :] = cell_frac_right # And then we add the new left-faces to the cell_face map. We do not # change the sign of the matrix since we did not flip the normals. # This means that the normals of right and left cells point in the same # direction, but their cell_faces values have oposite signs. # sparse_mat.stack_mat(g.cell_faces, cell_frac_left) - # g.cell_faces.tocsc() + g.cell_faces.tocsc() g.cell_faces = sps.vstack((g.cell_faces, cell_frac_left), format='csc') return 0 From 47d7743f9404b39feac936144494c622ca43d751 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 19:54:55 +0100 Subject: [PATCH 113/118] Fix of log warning in Fracture intersection. Cannot print indices for fractures that has none. --- src/porepy/fracs/fractures.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 4fa43bee8d..62115850a8 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -427,11 +427,17 @@ def intersects(self, other, tol, check_point_contact=True): # points if len(bound_sect_self_other) == 0 and len(bound_sect_other_self) == 0: if check_point_contact and int_points.shape[1] == 1: - # Point contacts are not implemented. Give a warning, return no + # contacts are not implemented. Give a warning, return no # interseciton, and hope the meshing software is merciful - logger.warning("""Found a point contact between fracture %i - and %i at (%.5f, %.5f, %.5f)""", self.index, - other.index, *int_points) + if hasattr(self, 'index') and hasattr(other, 'index'): + logger.warning("""Found a point contact between fracture + %i + and %i at (%.5f, %.5f, %.5f)""", self.index, + other.index, *int_points) + else: + logger.warning("""Found point contact between fractures + with no index. Coordinate: (%.5f, %.5f, + %.5f)""", *int_points) return np.empty((3, 0)), on_boundary_self, on_boundary_other # None of the intersection points lay on the boundary else: From b9cdce2790bd6daa7d30233f862d8693b4b3ba27 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 19:56:00 +0100 Subject: [PATCH 114/118] Fixes of fracture extrusion model. Recompute cell centers as the fracture geometry changes. --- src/porepy/fracs/extrusion.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index cf6131ba4b..d63fa83d57 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -412,7 +412,9 @@ def rotate_fracture(frac, vec, angle, exposure): rot = cg.rot(angle, vec) p = frac.p frac.p = exposure + rot.dot(p - exposure) + frac.points_2_ccw() + frac.compute_centroid() def impose_inlcine(fracs, exposure_line, exposure_point, frac_family=None, @@ -465,7 +467,7 @@ def impose_inlcine(fracs, exposure_line, exposure_point, frac_family=None, return all_ang def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, - recompute_center=False, **kwargs): + recompute_center=True, **kwargs): """ Cut a fracture by a plane, and confine it to one side of the plane. Intended use is to confine abutting fractures (T-intersections) to one side @@ -557,7 +559,7 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, if isect_pt.size == 0: print('No intersection found in cutting of fractures') - return main_frac + return main_frac, None # Next step is to eliminate points in the main fracture that are on the # wrong side of the other fracture. @@ -570,8 +572,6 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, eliminate = np.where(sgn * right_sign < 0)[0] main_frac.remove_points(eliminate) - if recompute_center: - main_frac.compute_center() # Add intersection points on the main fracture. One of these may already be # present, as the point of extrusion, but add_point will uniquify the point @@ -580,6 +580,9 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, # plane are present in the final fracture. main_frac.add_points(isect_pt, check_convexity=False) + if recompute_center: + main_frac.compute_centroid() + # If the main fracture is too large compared to the other, the cut line # will extend beyond the confining plane. In these cases, compute the # distance from the fracture center to the outside intersection point. This From eac9dd91a47cdc471f9e71e8799e66f478f61e25 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Nov 2017 21:33:15 +0100 Subject: [PATCH 115/118] Fracture has method to compute normal vector, extrusion uses it. Bug when the fracture was rotated, but the normal vector not updated. --- src/porepy/fracs/extrusion.py | 1 + src/porepy/fracs/fractures.py | 7 ++++++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index d63fa83d57..f6b914afbc 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -415,6 +415,7 @@ def rotate_fracture(frac, vec, angle, exposure): frac.points_2_ccw() frac.compute_centroid() + frac.compute_normal() def impose_inlcine(fracs, exposure_line, exposure_point, frac_family=None, diff --git a/src/porepy/fracs/fractures.py b/src/porepy/fracs/fractures.py index 62115850a8..7d688f122a 100644 --- a/src/porepy/fracs/fractures.py +++ b/src/porepy/fracs/fractures.py @@ -52,7 +52,7 @@ def __init__(self, points, index=None, check_convexity=True): # Ensure the points are ccw self.points_2_ccw() self.compute_centroid() - self.normal = cg.compute_normal(points)[:, None] + self.compute_normal() self.orig_p = self.p.copy() @@ -264,6 +264,11 @@ def compute_centroid(self): # Project back again. self.center = rot.transpose().dot(np.append(center, z)).reshape((3, 1)) + def compute_normal(self): + """ Compute normal to the polygon. + """ + self.normal = cg.compute_normal(self.p)[:, None] + def as_sp_polygon(self, p=None): """ Represent polygon as a sympy object. From 00a2690886ee526b8d02dc249bd5743896ebdabf Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 2 Nov 2017 08:05:00 +0100 Subject: [PATCH 116/118] Commented on unresolved issue in fracture extrusion. --- src/porepy/fracs/extrusion.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/src/porepy/fracs/extrusion.py b/src/porepy/fracs/extrusion.py index f6b914afbc..17d2b5c5de 100644 --- a/src/porepy/fracs/extrusion.py +++ b/src/porepy/fracs/extrusion.py @@ -14,14 +14,23 @@ No attempts are made to create fractures that do not cross the confined plane. +KNOWN ISSUES: + * When the extrusion is applied with rotations relative to outcrop plane, + two fractures may meet in a point only. A warning is issued in this case. + The consequences in terms of meshing are hard to predict. To fix this, it + is likely necessary to constrain rotation angles to the extrusion angles + in T-intersections. + """ import numpy as np +import warnings from porepy.utils import comp_geom as cg from porepy.fracs.fractures import EllipticFracture, Fracture -def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, **kwargs): +def fractures_from_outcrop(pt, edges, ensure_realistic_cuts=True, family=None, + **kwargs): """ Create a set of fractures compatible with exposed lines in an outcrop. See module-level documentation for futher comments. @@ -558,8 +567,16 @@ def cut_fracture_by_plane(main_frac, other_frac, reference_point, tol=1e-4, isect_pt, _, _ = main_frac.intersects(aux_frac, tol) + # The extension of the abutting fracture will cross the other fracture + # with a certain angle to the vertical. If the other fracture is rotated + # with a similar angle, point contact results. if isect_pt.size == 0: - print('No intersection found in cutting of fractures') + warnings.warn("""No intersection found in cutting of fractures. This is + likely caused by an unfortunate combination of + extrusion and rotation angles, which created fractures + that only intersect in a single point (the outcrop + plane. Will try to continue, but this may cause + trouble for meshing etc.""") return main_frac, None # Next step is to eliminate points in the main fracture that are on the From bf7956dd0e359533b1d562cc89c600efdb69e28e Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 2 Nov 2017 08:46:45 +0100 Subject: [PATCH 117/118] Fracture extrusion unit test with rotation --- test/unit/test_fracture_extrusion.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/unit/test_fracture_extrusion.py b/test/unit/test_fracture_extrusion.py index d06bc4ffe4..735aca0864 100644 --- a/test/unit/test_fracture_extrusion.py +++ b/test/unit/test_fracture_extrusion.py @@ -131,5 +131,21 @@ def test_cut_fracture_one_inclined(self): self.compare_arrays(p_known, f2.p) + def test_rotation(self): + # Test a case with a fracture terminated in both ends (H-configuration) + # This turned out to be problematic for certain cases. + # The test is simply that it runs - should have something more + # insightful here + pt = np.array([[ 0.4587156 , 0.76605504, 0.60586278, 0.74934585, + 0.46570401, 0.55721055], + [ 0.28593274, 0.35321103, 0.31814406, 0.028759 , + 0.44178218, 0.30749384]]) + edges = np.array([[0, 1], [2, 3], [4, 5]]).T + + np.random.seed(7) + family = np.zeros(3, dtype=np.int) + extrusion.fractures_from_outcrop(pt, edges, family=family, + family_std_incline=[0.3]) + if __name__ == '__main__': unittest.main() From 305296326d87e054d1671cfa7a5fa7820e86c0a9 Mon Sep 17 00:00:00 2001 From: Runar Date: Thu, 2 Nov 2017 09:02:50 +0100 Subject: [PATCH 118/118] Finished improving splitting of faces --- src/porepy/fracs/split_grid.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/porepy/fracs/split_grid.py b/src/porepy/fracs/split_grid.py index 8ef786f047..7fdad2d9bc 100644 --- a/src/porepy/fracs/split_grid.py +++ b/src/porepy/fracs/split_grid.py @@ -329,7 +329,7 @@ def update_cell_connectivity(g, face_id, normal, x0): row = cell_face_id[left_cell, 0] data = np.ravel(g.cell_faces[np.ravel(face_id[row]), col]) assert data.size == face_id.size - cell_frac_left = sps.csc_matrix((data, (row, col)), + cell_frac_left = sps.csr_matrix((data, (row, col)), (face_id.size, g.cell_faces.shape[1])) # We now update the cell_faces map of the faces on the right side of @@ -343,15 +343,17 @@ def update_cell_connectivity(g, face_id, normal, x0): (face_id.size, g.cell_faces.shape[1])) assert g.cell_faces.getformat() == 'csr' + sparse_mat.merge_matrices(g.cell_faces, cell_frac_right, face_id) # g.cell_faces[face_id, :] = cell_frac_right + # And then we add the new left-faces to the cell_face map. We do not # change the sign of the matrix since we did not flip the normals. # This means that the normals of right and left cells point in the same # direction, but their cell_faces values have oposite signs. - # sparse_mat.stack_mat(g.cell_faces, cell_frac_left) - g.cell_faces.tocsc() - g.cell_faces = sps.vstack((g.cell_faces, cell_frac_left), format='csc') + sparse_mat.stack_mat(g.cell_faces, cell_frac_left) + g.cell_faces = g.cell_faces.tocsc() + #g.cell_faces = sps.vstack((g.cell_faces, cell_frac_left), format='csc') return 0 @@ -513,13 +515,15 @@ def find_cell_color(g, cells): c = np.sort(cells) # Local cell-face and face-node maps. assert g.cell_faces.getformat() == 'csc' - cf_sub = sparse_mat.slice_mat(g.cell_faces, c) + cell_faces = sparse_mat.slice_mat(g.cell_faces, c) child_cell_ind = -np.ones(g.num_cells, dtype=np.int) - child_cell_ind[c] = np.arange(cf_sub.shape[1]) + child_cell_ind[c] = np.arange(cell_faces.shape[1]) # Create a copy of the cell-face relation, so that we can modify it at # will - cell_faces = cf_sub.copy() + # com Runar: I don't think this is neccessary as slice_mat creates a copy +# cell_faces = cf_sub.copy() + # Direction of normal vector does not matter here, only 0s and 1s cell_faces.data = np.abs(cell_faces.data)