diff --git a/README.md b/README.md index 609e89c..90ad730 100644 --- a/README.md +++ b/README.md @@ -8,4 +8,4 @@ If you want to preprocess your fracture network, make sure it is in the format: In order to perform meshing you need to have gmsh installed and it must be recognized as an environment variable (i.e., you must be able to call it from the command line with "gmsh"). -Only numpy and scipy (and gmsh for meshing) are required to run this preprocessing tool. +Only numpy, scipy, matplotlib, and igraph (and gmsh for meshing) are required to run this preprocessing tool. diff --git a/graph_code.py b/graph_code.py new file mode 100644 index 0000000..9c87b66 --- /dev/null +++ b/graph_code.py @@ -0,0 +1,912 @@ +""" +MIT License +Copyright (c) 2021 Stephan de Hoop S.dehoop-1@tudelft.nl + Denis Voskov D.V.Voskov@tudelft.nl + Delft University of Technology, the Netherlands +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +""" +import numpy as np +import matplotlib.pyplot as plt +import igraph +import copy +import matplotlib.colors as colors +import matplotlib.cm as cm +import os + + +class Graph: + """ + (weighted) Graph class containing a list of Vertices (x- and y-coordinates) and Edges (list of + (v_i, v_j) pairs), where v_i and v_j are vertex i and j respectively) + """ + def __init__(self, max_vertices=int(1e5), max_edges=int(1e5), edge_data=None, matrix_perm=1000, fracture_aperture=1e-4): + self.vertices = np.ones((max_vertices, 2), dtype=float) * np.inf + self.edge_to_vertex = np.ones((max_edges, 2), dtype=int) * -1 + self.apertures = np.ones((max_edges,), dtype=float) * fracture_aperture + self.heat_transfer_mult = np.ones((max_edges,), dtype=float) + self.volumes = np.zeros((max_edges,), dtype=float) + self.active_edges = np.zeros((max_edges,), dtype=bool) + self.active_vertices = np.zeros((max_vertices,), dtype=bool) + self.vertex_to_edge = [] + self.num_edges = 0 + self.matrix_perm = matrix_perm * 1e-15 + self.matrix_eff_aperture = np.sqrt(12 * self.matrix_perm) + + if edge_data is not None: + # Construct graph based on edge_data: + self.add_multiple_edges(edge_data) + + def add_edge(self, x1, y1, x2, y2): + """ + Add edge to self.edge_to_vertex matrix [vertex_id_i, vertex_id_j], and adds new edge to self.vertex_to_edge for + vertex_id_i and vertex_id_j + :param x1: x-coordinate of vertex_id_i + :param y1: y-coordinate of vertex_id_i + :param x2: x-coordinate of vertex_id_j + :param y2: y-coordinate of vertex_id_j + :return: + """ + # Find or insert both vertices from edge: + vertex_id_1 = self.find_or_insert_vertex(x1, y1) + vertex_id_2 = self.find_or_insert_vertex(x2, y2) + + if vertex_id_1 == vertex_id_2: + return 0 + + edge_id = self.num_edges + self.edge_to_vertex[edge_id, 0] = vertex_id_1 + self.edge_to_vertex[edge_id, 1] = vertex_id_2 + self.num_edges += 1 + self.active_edges[edge_id] = True + self.volumes[edge_id] = self.apertures[edge_id] * np.linalg.norm(self.vertices[vertex_id_1].flatten() - + self.vertices[vertex_id_2].flatten()) + + self.vertex_to_edge[vertex_id_1].append(edge_id) + self.vertex_to_edge[vertex_id_2].append(edge_id) + return 0 + + def add_multiple_edges(self, edge_data): + """ + Loop over edgedata to construct graph + :param edge_data: N x 4 array that contains all the edges to be added to the domain, for should be: + [[x1, y1, x2, y2], ..., []] + :return: + """ + for edge in edge_data: + self.add_edge(edge[0], edge[1], edge[2], edge[3]) + return 0 + + def get_num_vertices(self): + return len(self.vertex_to_edge) + + def get_num_edges(self): + return self.num_edges + + def insert_vertex(self, x: float, y: float): + """ + Insert new vertex into self.vertices and update number of existing vertices + NOTE: This method currently assumes that the vertex with coordinates (x,y) is not yet in self.vertices!!! + :param x: x-coordinate of vertex + :param y: y-coordinate of vertex + :return: the new id of vertex + """ + vertex_id = self.get_num_vertices() + self.vertices[vertex_id, 0] = x + self.vertices[vertex_id, 1] = y + self.vertex_to_edge.append([]) + self.active_vertices[vertex_id] = True + return vertex_id + + def get_vertex_id(self, x, y): + """ + Get vertex id if already in the domain, otherwise return 0, also catches if there exists duplicated vertices in + the domain (i.e., when results.size > 1) + :param x: x-coordinate of vertex + :param y: y-coordinate of vertex + :return: + """ + # linear search is a potential problem + # need to introduce additional, sorted by coordinate array, coord_to_vertex_id + result = np.where((self.vertices[0:self.get_num_vertices()] == (x, y)).all(axis=1))[0] + # assert (result.size == 1) + if result.size == 1: + return result[0] + elif result.size == 0: + return None + else: + # Duplicate vertex... + print('Duplicate vertex found in self.get_vertex_id with coordinates [x,y] = [{:}, {:}]'.format(x,y)) + return np.NaN + + def find_or_insert_vertex(self, x, y): + """ + Flexible method which finds the id of an existing vertex or adds a potentially new vertex to the graph + :param x: x-coordinate of vertex + :param y: y-coordinate of vertex + :return: + """ + vertex_id = self.get_vertex_id(x, y) + if vertex_id is None: + # New vertex: + vertex_id = self.insert_vertex(x, y) + return vertex_id + + def remove_edges(self, ghost_edge_to_vertices, ghosted_edges): + """ + Removes edges that were flagged as ghosts (edge_to_vertex remains unchanged in size, to simplify bookkeeping, + edges are removed from vertex_to_edge to keep node degree accurate + :param ghost_edge_to_vertices: list of vertices for each edge that is removed after merging vertex + :param ghosted_edges: list of edges that are removed after merging vertex + :return: + """ + for ii in range(len(ghosted_edges)): + edge = ghosted_edges[ii] + for jj in ghost_edge_to_vertices[ii]: + self.vertex_to_edge[jj].remove(edge) + + self.edge_to_vertex[edge] = -1 + self.active_edges[edge] = False + + def find_conflicting_edges(self, vertex_id_from, vertex_id_to): + """ + For each edge leaving vertex_from, check what is the status after merging and flag for ghosting if edge is + collapsed or overlaps after merging + :param vertex_id_from: vertex that is removed after merging + :param vertex_id_to: vertex that is merged into + :return: + """ + ghosted_edges = [] + ghost_edge_to_vertices = [] + status_after_merge = [None] * len(self.vertex_to_edge[vertex_id_from]) + edge_id_after_merge = [None] * len(self.vertex_to_edge[vertex_id_from]) + edge_to_vertex_after_merge = np.ones((len(self.vertex_to_edge[vertex_id_from]), 2), dtype=int) * -1 + edges = self.vertex_to_edge[vertex_id_from] + for ii, edge in enumerate(edges): + status_after_merge[ii], edge_to_vertex_after_merge[ii], edge_id_after_merge[ii] = \ + self.status_edge_after_merge(vertex_id_from, vertex_id_to, edge) + if status_after_merge[ii] == 'collapsed' or status_after_merge[ii] == 'overlap': + # Edge will be ghosted: + ghosted_edges.append(edge) + ghost_edge_to_vertices.append(list(self.edge_to_vertex[edge])) + + # self.visualize_sub_graph(self.edge_to_vertex[list(set(self.vertex_to_edge[vertex_id_from] + self.vertex_to_edge[vertex_id_to]))]) + return ghost_edge_to_vertices, ghosted_edges, status_after_merge, edge_to_vertex_after_merge, edge_id_after_merge + + def merge_vertices(self, vertex_id_from, vertex_id_to, char_len, correct_aperture): + """ + Merge two vertices in domain while preserving connections and ghosting the merged vertex + :param vertex_id_from: vertex that is merged in vertex_id_to (is ghosted) + :param vertex_id_to: target vertex which vertex_id_from is merged into (stays in domain) + :return: + """ + # ghost any edge that "collapses" (i.e., the edge is shared by both vertices and has zero length after + # vertex-merging) or if the edge will overlap after merging. + ghost_edge_to_vertices, ghosted_edges, status_after_merging, edge_to_vertex_after_merging, edge_id_after_merge = \ + self.find_conflicting_edges(vertex_id_from, vertex_id_to) + if correct_aperture: + self.update_volumes_and_apertures(vertex_id_from, vertex_id_to, status_after_merging, edge_to_vertex_after_merging, edge_id_after_merge, char_len) + self.remove_edges(ghost_edge_to_vertices, ghosted_edges) + + # 1. preserve connections to ghost-vertex + [self.vertex_to_edge[vertex_id_to].append(x) for x in self.vertex_to_edge[vertex_id_from] if + x not in self.vertex_to_edge[vertex_id_to]] + # self.vertex_to_edge[vertex_id_to] = set(self.vertex_to_edge[vertex_id_to] + self.vertex_to_edge[vertex_id_from]) + + # 2. apply ghosting to vertex: + self.vertices[vertex_id_from] = np.inf + self.active_vertices[vertex_id_from] = False + + # 3. update all instances of ghost-vertex in edge_to_vertex to new vertex_id + slice_edge_to_vertex = self.edge_to_vertex[0:self.get_num_edges()] + slice_edge_to_vertex[slice_edge_to_vertex == vertex_id_from] = vertex_id_to + return 0 + + def check_distance_constraint(self, new_vertex, existing_vertices, char_len, merge_threshold): + """ + Computes the distance between all vertices already in domain and to be added vertex + :param new_vertex: to be added vertex + :param existing_vertices: vertices already approved + :param char_len: radius for algebraic constraint + :param merge_threshold: h-factor which is recommended between [0.5, 0.86] + :return: + """ + assert 0.5 <= merge_threshold <= 0.86, "Choose threshold on closed interval [0.5, 0.86]" + dist_vec = np.linalg.norm(self.vertices[new_vertex] - self.vertices[existing_vertices], axis=1) + argmin_id = np.argmin(dist_vec) + fixed_vertex = existing_vertices[argmin_id] + + if dist_vec[argmin_id] < char_len * merge_threshold: + return fixed_vertex + else: + return -1 + + def closest_point_method(self, order_discr, char_len, merge_threshold, correct_aperture): + """ + Main sequential algorithm which performs cleaning of the fracture network based on algebraic constraint + :param order_discr: order in which vertices are sequentially added to the domain and checked for constraint + :param char_len: radius within which vertices are violating constraint + :param merge_threshold: h-factor which is recommended between [0.5, 0.86] + :param correct_aperture: boolean for applying aperture correction or not + :return: + """ + assert 0.5 <= merge_threshold <= 0.86, "Choose threshold on closed interval [0.5, 0.86]" + count = 0 + for new_vertex in order_discr[1:]: + count += 1 + if not self.active_vertices[new_vertex]: + continue + fixed_vertex = self.check_distance_constraint(new_vertex, order_discr[:count], char_len, merge_threshold) + if fixed_vertex >= 0: + self.merge_vertices(new_vertex, fixed_vertex, char_len, correct_aperture) + return 0 + + def calc_angles_vectors_2D(self, vec_m, vec_p, angle_360=False): + """ + Computes the angle between any two vectors in 2D + :param vec_m: first vector + :param vec_p: second vector + :param angle_360: boolean for returning angle in 360 or 180 180 degrees spectrum + :return: + """ + vec_m = vec_m / np.linalg.norm(vec_m) + vec_p = vec_p / np.linalg.norm(vec_p) + dot_product = min(1, max(-1, np.dot(vec_m, vec_p))) + + if angle_360: + det_product = vec_m[0] * vec_p[1] - vec_m[1] * vec_p[0] + angle_full = np.arctan2(-det_product, dot_product) * 180 / np.pi + if angle_full < 0: + return 360 + angle_full + else: + return angle_full + else: + return np.arccos(dot_product) * 180 / np.pi + + def calc_angle_edges(self, edge_1, edge_2, angle_360=False): + """ + Computes the angle between any two edges + :param edge_1: first edge (id, not (x,y)-coordinates! + :param edge_2: second edge (id, not (x,y)-coordinates! + :param angle_360: boolean for returning angle in 360 or 180 180 degrees spectrum + :return: + """ + common_vertex = np.intersect1d(self.edge_to_vertex[edge_1], self.edge_to_vertex[edge_2]) + other_vertex_1 = self.edge_to_vertex[edge_1, self.edge_to_vertex[edge_1] != common_vertex] + other_vertex_2 = self.edge_to_vertex[edge_2, self.edge_to_vertex[edge_2] != common_vertex] + vec_m = self.vertices[common_vertex].flatten() - self.vertices[other_vertex_1].flatten() + vec_p = self.vertices[common_vertex].flatten() - self.vertices[other_vertex_2].flatten() + return self.calc_angles_vectors_2D(vec_m, vec_p, angle_360=angle_360) + + def check_angles_leaving_vertices(self, vertex_id): + """ + Calculates all the angles between all vertices leaving vertex_id + :param vertex_id: integer with id for vertex + :return: + """ + edges = np.array(self.vertex_to_edge[vertex_id]) + num_edges = len(edges) + if num_edges == 2: + return self.calc_angle_edges(edges[0], edges[1]), np.array([[edges[0], edges[1]]]) + elif num_edges == 1: + return 360, None + + angles = np.zeros((num_edges,)) + + # Determine order, starting with 0 angle at positive y-axis [0, 1]: + angles_y_axis = np.zeros((num_edges,)) + vec_p = np.array([1, 0]) + for ii, edge in enumerate(edges): + other_vertex = self.edge_to_vertex[edge][self.edge_to_vertex[edge] != vertex_id] + vec_m = self.vertices[other_vertex].flatten() - self.vertices[vertex_id].flatten() + angles_y_axis[ii] = self.calc_angles_vectors_2D(vec_m, vec_p, angle_360=True) + + sort_edges_ids = edges[np.argsort(angles_y_axis)] + edge_pair = np.zeros((num_edges, 2), dtype=int) + for ii in range(num_edges - 1): + edge_pair[ii] = [sort_edges_ids[ii + 1], sort_edges_ids[ii]] + angles[ii] = self.calc_angle_edges(sort_edges_ids[ii + 1], sort_edges_ids[ii], angle_360=True) + edge_pair[-1] = [sort_edges_ids[0], sort_edges_ids[-1]] + angles[-1] = self.calc_angle_edges(sort_edges_ids[0], sort_edges_ids[-1], angle_360=True) + + if abs(np.sum(angles) - 360) > 1e-4: + print('Mmmhhh...') + return angles, edge_pair + + def straighten_edges(self, tolerance_angle, char_len, correct_aperture): + """ + Method which straightens fractures if within certain tolerance from a straight line + :param tolerance_angle: deviation from straight line + :param char_len: radius of algebraic constraint + :param correct_aperture: boolean if applying aperture correction + :return: + """ + for ii in range(self.get_num_vertices()): + if not self.active_vertices[ii]: + continue + edges = self.vertex_to_edge[ii] + if len(edges) == 2: + angle = self.calc_angle_edges(edges[0], edges[1]) + + if np.abs(angle - 180) <= tolerance_angle: + # Edge straight enough to be merged into closest vertex: + other_vertices = self.edge_to_vertex[edges].flatten()[np.where(self.edge_to_vertex[edges].flatten() != ii)[0]] + dist_vertices = np.linalg.norm(self.vertices[ii] - self.vertices[other_vertices], axis=1) + closests_id = np.argmin(dist_vertices) + self.merge_vertices(vertex_id_from=ii, vertex_id_to=other_vertices[closests_id], char_len=char_len, correct_aperture=correct_aperture) + return 0 + + def remove_small_angles(self, tolerange_small_angle, char_len, correct_aperture): + """ + Method which removes small angles which might result in "skinny" triangles in meshed results + :param tolerange_small_angle: max. allowable angle between two fracture segments + :param char_len: radius of algebraic constraint + :param correct_aperture: boolean if applying aperture correction + :return: + """ + active_ids = np.where(self.active_vertices)[0] + + for vertex_id in active_ids: + if not self.active_vertices[vertex_id]: + continue + + angles, edge_ids = self.check_angles_leaving_vertices(vertex_id) + while np.any(angles < tolerange_small_angle): + # Determine which angle is the smallest and merge two edges into one, logic: keep largest edge and merge + # non-shared vertex into closest vertex, then update angles and continue + edge_pair = edge_ids[np.argmin(angles)] + other_vertices = self.edge_to_vertex[edge_pair].flatten()[self.edge_to_vertex[edge_pair].flatten() != vertex_id] + length_vec = np.linalg.norm(self.vertices[vertex_id] - self.vertices[other_vertices], axis=1) + vertex_from = other_vertices[np.argmin(length_vec)] + large_edge = edge_pair[np.argmax(length_vec)] + + dist_vec = np.linalg.norm(self.vertices[vertex_from] - self.vertices[self.edge_to_vertex[large_edge]], axis=1) + vertex_to = self.edge_to_vertex[large_edge][np.argmin(dist_vec)] + self.merge_vertices(vertex_from, vertex_to, char_len, correct_aperture) + angles, edge_ids = self.check_angles_leaving_vertices(vertex_id) + return 0 + + def simplify_graph(self, order_discr, char_len, merge_treshold=0.66, + tolerange_small_angle=20, small_angle_iter=2, tolerange_straight_angle=7.5, correct_aperture=True, straighten_edges=False): + """ + Method which performs actual simplification of the graph + :param order_discr: order in which vertices are sequentially checked on algebraic constraint + :param char_len: radius of algebraic constraint + :param merge_treshold: h-factor which is recommended between [0.5, 0.86] + :param tolerange_small_angle: max. allowable angle between two fracture segments + :param small_angle_iter: number of times to run small angle correction, merging vertices changes angle and might + introduce new accute intersection + :param tolerange_straight_angle: allowable deviation from straight line + :param correct_aperture: boolean to applying aperture correction or not + :param straighten_edges: boolean to straighten edges (calling method which does this) after preprocessing to + speed up succesive gridding + :return: + """ + assert 0.5 <= merge_treshold <= 0.86, "Choose threshold on closed interval [0.5, 0.86]" + self.closest_point_method(order_discr, char_len, merge_treshold, correct_aperture) + for ii in range(small_angle_iter): + self.remove_small_angles(tolerange_small_angle, char_len, correct_aperture) + if straighten_edges: + self.straighten_edges(tolerange_straight_angle, char_len, correct_aperture) + return 0 + + def visualize_edges(self, vertex_id): + """ + Method which plots fracture network around vertex_id + :param vertex_id: integer id for current vertex + :return: + """ + edges = self.vertex_to_edge[vertex_id] + leaving_vertices = self.edge_to_vertex[edges].flatten()[np.where(self.edge_to_vertex[edges].flatten() != vertex_id)[0]] + plt.figure() + plt.plot( + np.vstack((self.vertices[self.edge_to_vertex[edges, 0], 0].T, + self.vertices[self.edge_to_vertex[edges, 1], 0].T)), + np.vstack((self.vertices[self.edge_to_vertex[edges, 0], 1].T, + self.vertices[self.edge_to_vertex[edges, 1], 1].T)), + color='black' + ) + plt.axis('equal') + plt.plot(self.vertices[leaving_vertices, 0], self.vertices[leaving_vertices, 1], '.', color='red') + plt.plot(self.vertices[vertex_id, 0], self.vertices[vertex_id, 1], '.', color='blue') + plt.show() + + def visualize_graph(self): + """ + Method which visualize current graph state + :return: + """ + num_edges = self.get_num_edges() + num_vertices = self.get_num_vertices() + plt.figure() + plt.plot( + np.vstack((self.vertices[self.edge_to_vertex[:num_edges, 0], 0].T, + self.vertices[self.edge_to_vertex[:num_edges, 1], 0].T)), + np.vstack((self.vertices[self.edge_to_vertex[:num_edges, 0], 1].T, + self.vertices[self.edge_to_vertex[:num_edges, 1], 1].T)), + color='black' + ) + plt.axis('equal') + plt.plot(self.vertices[:num_vertices, 0], self.vertices[:num_vertices, 1], '.', color='red') + plt.show() + return 0 + + def visualize_sub_graph(self, edges): + """ + Visualize subgraph based on list of edges + :param edges: subset of self.edge_to_vertex containing vertex pairs of sub-graph + :return: + """ + vertices = list(set(edges.flatten())) + epsilon = 0.1 + plt.figure() + for edge in edges: + plt.plot( + np.vstack((self.vertices[edge[0], 0].T, + self.vertices[edge[1], 0].T)), + np.vstack((self.vertices[edge[0], 1].T, + self.vertices[edge[1], 1].T)), + color='black' + ) + plt.text((self.vertices[edge[0], 0] + self.vertices[edge[1], 0]) / 2 + epsilon, + (self.vertices[edge[0], 1] + self.vertices[edge[1], 1]) / 2 + epsilon, + str(np.intersect1d(self.vertex_to_edge[edge[0]], self.vertex_to_edge[edge[1]])[0]), fontsize=10) + plt.axis('equal') + for ii in vertices: + plt.plot(self.vertices[ii, 0], self.vertices[ii, 1], '.', color='red') + plt.text(self.vertices[ii, 0] + epsilon, self.vertices[ii, 1] + epsilon, + str(ii), fontsize=10) + plt.show() + return 0 + + def check_collapsed_edges(self): + """ + Safety method to check graph consistency on collapsed edges that might occur in domain + :return: + """ + check_result = False + + active_edge_ids = np.where(self.active_edges)[0] + for edge in active_edge_ids: + vertices = self.edge_to_vertex[edge] + if vertices[0] == vertices[1]: + check_result = True + break + return check_result + + def check_duplicate_edges(self): + """ + Safety method to check graph consistency on duplicate edges that might occur in domain + :return: + """ + check_result = False + slice_edge_to_vertex = self.edge_to_vertex[0:self.get_num_edges()] + unique_edges = np.unique(slice_edge_to_vertex, axis=0) + if unique_edges.shape[0] != slice_edge_to_vertex.shape[0]: + check_result = True + return check_result + + def check_graph_consistency(self): + """ + Check consistency of the graph + :return: + """ + self.check_collapsed_edges() + self.check_duplicate_edges() + return 0 + + def connectivity_merging_vertices(self, vertex_from, vertex_to, radius): + """ + Compute connectivity between vertex_from and vertex_to using a sub-graph and Dijkstra's shortest path + :param vertex_from: id for vertex which is merged + :param vertex_to: id for vertex which is fixed + :param radius: radius around which vertices and edges are extract for the subgraph, choosing this parameter too + large can result in unexpected results or long preprocessing time + :return: resistance (L / aperture^2), is infinity if no direct connection exists + """ + # Check if nodes are first-degree connected + common_edge = np.intersect1d(self.vertex_to_edge[vertex_from], self.vertex_to_edge[vertex_to]) + if common_edge.size != 0: + return np.linalg.norm(self.vertices[vertex_from].flatten() - self.vertices[vertex_to].flatten()) / (self.apertures[common_edge[0]] ** 2) + else: + # Extract subpgrah of vertices and edges and their weights (within a radius of char_len * X around vertices) + vertices_in_radius = np.where(np.linalg.norm(self.vertices[:self.get_num_vertices()] - + self.vertices[vertex_from], axis=1) < radius)[0] + + edges_in_radius = [] + for ii in vertices_in_radius: + edges_in_radius += self.vertex_to_edge[ii] + edges_in_radius = list(set(edges_in_radius)) + vertices_in_radius = list(set(self.edge_to_vertex[edges_in_radius].flatten())) + if vertex_to not in vertices_in_radius: + return np.inf + + edges = self.edge_to_vertex[edges_in_radius] + edge_weights = np.linalg.norm(self.vertices[edges[:, 0]] - self.vertices[edges[:, 1]], axis=1) / \ + self.apertures[edges_in_radius] ** 2 # m to km and m for aperture to cm (note: ** 2) in order to avoid extremely large values + + # Construct subgraph with weights in igraph and calculate shortest_path between vertices + g = igraph.Graph(edges=edges, edge_attrs={'weight': edge_weights}) + dist_from_to = g.shortest_paths(source=vertex_from, target=vertex_to, weights=edge_weights) + return dist_from_to[0][0] + + def status_edge_after_merge(self, vertex_from, vertex_to, edge): + """ + Determines status of edge leaving vertex_from when applying merge into vertex_to + :param vertex_from: vertex id which is merged + :param vertex_to: vertex id which remains in domain + :param edge: id of edge which is leaving vertex_from + :return: status of edge (str), list of two vertices making up new edge, id of new edge (or old, if simply + extended) + """ + if np.all(np.sort(self.edge_to_vertex[edge]) == np.sort(np.array([vertex_from, vertex_to]))): + # Edge will collapse: + return 'collapsed', np.array([vertex_to, vertex_to]), -1 + else: + # Loop over all edges in vertex_to and check of the new edge will overlap + new_edge = np.sort(np.array([vertex_to, self.edge_to_vertex[edge, self.edge_to_vertex[edge] != vertex_from][0]], dtype=int)) + overlap = False + new_edge_id = -1 + for curr_edge in self.vertex_to_edge[vertex_to]: + if np.all(np.sort(self.edge_to_vertex[curr_edge]) == new_edge): + overlap = True + new_edge_id = curr_edge + if overlap: + return 'overlap', new_edge, new_edge_id + else: + return 'extension', new_edge, edge + + def calc_effective_aperture_and_heat_transfer_sequential(self, new_edge, collapsed_edge_id, extended_edge_id, resistance, vertex_from, vertex_to): + """ + Calculate effect aperture and heat transfer for Type 2 corrections (similar to sequential resistors in circuit) + :param new_edge: list of pair of vertices making up new edge + :param collapsed_edge_id: list with edge id if edge of other edge is it collapsed (is empty if no edge connects + to this vertex_from + :param extended_edge_id: id of the edge that is currently evaluated ( + :param resistance: L / aperture^2 , computed using Dijkstra's shortest path (== inf if it doesn't exist) + :param vertex_from: id vertex which is merged + :param vertex_to: id vertex which stays in domain + :return: effective_aperture, effective_heattransfer + """ + len_new_edge = np.linalg.norm(self.vertices[new_edge[0]].flatten() - self.vertices[new_edge[1]].flatten()) + len_old_edge = np.linalg.norm(self.vertices[self.edge_to_vertex[extended_edge_id, 0]].flatten() - + self.vertices[self.edge_to_vertex[extended_edge_id, 1]].flatten()) + + dist_gap = np.linalg.norm(self.vertices[vertex_from].flatten() - self.vertices[vertex_to].flatten()) + weights = np.array([len_old_edge, dist_gap]) / (len_old_edge + dist_gap) + n_pow = 9.5608 * weights[1] + 1.18024 # fitted from tested different gaps and different lc!!! + if resistance == np.inf: + if self.matrix_eff_aperture < self.apertures[extended_edge_id]: + eff_aperture = 1 / ((1 - weights[1] ** n_pow) / (self.apertures[extended_edge_id]) + (weights[1] ** n_pow) / (self.matrix_eff_aperture)) + else: + weights = np.array([len_old_edge, dist_gap]) / (len_old_edge + dist_gap) + eff_aperture = (self.apertures[extended_edge_id] * weights[0]) + (self.matrix_eff_aperture * weights[1]) + eff_heat_transfer = 1 + + else: + dist_gap = np.max((len_new_edge - len_old_edge, 0.1)) + weights = np.array([len_old_edge, dist_gap]) / (len_old_edge + dist_gap) + eff_gap_aper = np.sqrt(dist_gap / resistance) + if eff_gap_aper > np.max(self.apertures): + eff_gap_aper = np.max(self.apertures) + eff_aperture = 1 / (weights[0] / (self.apertures[extended_edge_id]) + weights[1] / eff_gap_aper) + if collapsed_edge_id.size == 0: + len_other_edge = 0 + else: + len_other_edge = np.linalg.norm(self.vertices[self.edge_to_vertex[collapsed_edge_id[0], 0]].flatten() - + self.vertices[self.edge_to_vertex[collapsed_edge_id[0], 1]].flatten()) + eff_heat_transfer = (len_other_edge + len_old_edge) / len_new_edge + + return eff_aperture, eff_heat_transfer + + def calc_effective_aperture_and_heat_transfer_parallel(self, old_edge_id, new_edge_id): + """ + Simple arithmetic mean weighted by length of the two edges which are overlapping after merge (similar to + parallel resistors in circuit) + :param old_edge_id: id of the old edge after merging + :param new_edge_id: id of the new edge after merging + :return: effective_aperture, effective_heattransfer + """ + len_new_edge = np.linalg.norm(self.vertices[self.edge_to_vertex[new_edge_id, 0]].flatten() - + self.vertices[self.edge_to_vertex[new_edge_id, 1]].flatten()) + len_old_edge = np.linalg.norm(self.vertices[self.edge_to_vertex[old_edge_id, 0]].flatten() - + self.vertices[self.edge_to_vertex[old_edge_id, 1]].flatten()) + eff_aperture = np.sqrt(len_new_edge * (self.apertures[new_edge_id] ** 2 / len_new_edge + + self.apertures[old_edge_id] ** 2 / len_old_edge)) + eff_heat_transfer = (len_new_edge + len_old_edge) / len_new_edge + return eff_aperture, eff_heat_transfer + + def update_volume_collapsed_edge(self, vertex_from, vertex_to, edge_id): + """ + Updates volumes after an edge is collapsed, basically distributes the volumes to all connecting edges to + vertex_from + :param vertex_from: id vertex which is merged + :param vertex_to: id vertex which remains in domain + :param edge_id: id of edge currently under consideration + :return: + """ + # Find if neighbouring vertices have edges, if both have then distribute 50% to both vertices and edges respectively + edges_vertex_from = copy.deepcopy(self.vertex_to_edge[vertex_from]) + edges_vertex_to = copy.deepcopy(self.vertex_to_edge[vertex_to]) + volume_collapsed_edge = copy.deepcopy(self.volumes[edge_id]) + if edge_id in edges_vertex_from and edge_id in edges_vertex_to: + other_edges_vertex_from = edges_vertex_from + other_edges_vertex_from.remove(edge_id) + for edge in other_edges_vertex_from: + self.volumes[edge] += volume_collapsed_edge / len(other_edges_vertex_from) / 2 + + other_edges_vertex_to = edges_vertex_to + other_edges_vertex_to.remove(edge_id) + for edge in other_edges_vertex_to: + self.volumes[edge] += volume_collapsed_edge / len(other_edges_vertex_to) / 2 + + elif edge_id in edges_vertex_from: + other_edges_vertex_from = edges_vertex_from + other_edges_vertex_from.remove(edge_id) + for edge in other_edges_vertex_from: + self.volumes[edge] += volume_collapsed_edge / len(other_edges_vertex_from) + + elif edge_id in edges_vertex_to: + other_edges_vertex_to = edges_vertex_to + other_edges_vertex_to.remove(edge_id) + for edge in other_edges_vertex_to: + self.volumes[edge] += volume_collapsed_edge / len(other_edges_vertex_to) + self.volumes[edge_id] = 0 + return 0 + + def update_volume_overlap_edge(self, old_edge_id, new_edge_id): + """ + Update volume of overlapping and removed edge + :param old_edge_id: id of edge which is removed after merging + :param new_edge_id: id of edge which remains in domain after merging + :return: + """ + # Add volume old edge to new edge + self.volumes[new_edge_id] += copy.deepcopy(self.volumes[old_edge_id]) + self.volumes[old_edge_id] = 0 + return 0 + + def update_volumes_and_apertures(self, vertex_from, vertex_to, status_after_merging, + edge_to_vertex_after_merging, edge_id_after_merge, char_len): + """ + Main function that updates volumes and apertures after merging + :param vertex_from: id vertex which gets merged + :param vertex_to: id vertex which remains in domain (target vertex) + :param status_after_merging: status of edges connecting to vertex_from after merging + :param edge_to_vertex_after_merging: pair of vertices of edges leaving vertex_from after merging + :param edge_id_after_merge: ids of edges leaving vertex_from after merging + :param char_len: radius of algebraic constraint + :return: + """ + # First check if any of the edges leaving vertex_from is collapsing, distribute volume first, then loop + # over the remaining edges to update properties accordingly + resistance = self.connectivity_merging_vertices(vertex_from, vertex_to, char_len * 2.5) + if 'collapsed' in status_after_merging: + self.update_volume_collapsed_edge(vertex_from, vertex_to, + self.vertex_to_edge[vertex_from][status_after_merging.index('collapsed')]) + + for ii, status_edge in enumerate(status_after_merging): + if status_edge == 'collapsed': + continue + elif status_edge == 'overlap': + # Parallel: + old_edge_id = self.vertex_to_edge[vertex_from][ii] + new_edge_id = edge_id_after_merge[ii] + eff_aperture, eff_heat_transfer = self.calc_effective_aperture_and_heat_transfer_parallel(old_edge_id, new_edge_id) + self.apertures[new_edge_id] = eff_aperture + self.heat_transfer_mult[new_edge_id] = eff_heat_transfer + self.update_volume_overlap_edge(old_edge_id, new_edge_id) + elif status_edge == 'extension': + # Sequential: + new_edge = edge_to_vertex_after_merging[ii] + collapsed_edge_id = np.intersect1d(self.vertex_to_edge[vertex_from], self.vertex_to_edge[vertex_to]) + extended_edge_id = self.vertex_to_edge[vertex_from][ii] + eff_aperture, eff_heat_transfer = self.calc_effective_aperture_and_heat_transfer_sequential( + new_edge, collapsed_edge_id, extended_edge_id, resistance, vertex_from, vertex_to) + self.apertures[extended_edge_id] = eff_aperture + self.heat_transfer_mult[extended_edge_id] = eff_heat_transfer + # self.visualize_sub_graph(self.edge_to_vertex[list(set(self.vertex_to_edge[vertex_from] + self.vertex_to_edge[vertex_to]))]) + return 0 + + def visualize_graph_with_volume_weights(self, min_val=None, max_val=None): + """ + Visualze edges in graph by volume + :param min_val: minimum volume for scaling colorbar and data + :param max_val: maximum volume for scaling colorbar and data + :return: + """ + fracs = copy.deepcopy(self.volumes[self.active_edges]) + if min_val is None: + min_val = np.min(fracs) + + if max_val is None: + max_val = np.min(fracs) + + fracs[fracs < min_val] = min_val + fracs[fracs > max_val] = max_val + norm = colors.Normalize(min_val, max_val) + colors_aper = cm.viridis(norm(fracs)) + + plt.figure() + for jj in range(self.get_num_edges()): + if not self.active_edges[jj]: + continue + plt.plot( + np.vstack((self.vertices[self.edge_to_vertex[jj, 0], 0].T, + self.vertices[self.edge_to_vertex[jj, 1], 0].T)), + np.vstack((self.vertices[self.edge_to_vertex[jj, 0], 1].T, + self.vertices[self.edge_to_vertex[jj, 1], 1].T)), + color=colors_aper[jj, :-1] + ) + plt.axis('equal') + plt.title('Volume Weights') + plt.show() + return + + def visualize_graph_with_aperture_weights(self, min_val=None, max_val=None): + """ + Visualze edges in graph by aperture + :param min_val: minimum aperture for scaling colorbar and data + :param max_val: maximum aperture for scaling colorbar and data + :return: + """ + fracs = copy.deepcopy(self.apertures[self.active_edges]) + if min_val is None: + min_val = np.min(fracs) + + if max_val is None: + max_val = np.min(fracs) + + fracs[fracs < min_val] = min_val + fracs[fracs > max_val] = max_val + norm = colors.Normalize(min_val, max_val) + colors_aper = cm.viridis(norm(fracs)) + + plt.figure() + for jj in range(self.get_num_edges()): + if not self.active_edges[jj]: + continue + plt.plot( + np.vstack((self.vertices[self.edge_to_vertex[jj, 0], 0].T, + self.vertices[self.edge_to_vertex[jj, 1], 0].T)), + np.vstack((self.vertices[self.edge_to_vertex[jj, 0], 1].T, + self.vertices[self.edge_to_vertex[jj, 1], 1].T)), + color=colors_aper[jj, :-1] + ) + plt.axis('equal') + plt.title('Aperture Weights') + plt.show() + return + + +def create_geo_file(act_frac_sys, filename, decimals, + height_res, char_len, box_data, char_len_boundary, export_frac=True): + """ + Creates geo file which serves as input to gmsh + :param act_frac_sys: list of fractures in domain in format [[x1, y1, x2, y2], [...], ..., [...]] + :param filename: name of the resulting geo-file + :param decimals: data is rounded off to this number of decimals + :param height_res: height of the resulting 1-layer 3D reservoir + :param char_len: characteristic length of the resulting mesh + :param box_data: coordinates of the box-data around the fracture network + :param char_len_boundary: characteristic length of mesh elements at the boundary + :param export_frac: boolean which exports fractures into the meshed file + :return: + """ + act_frac_sys = np.round(act_frac_sys * 10 ** decimals) * 10 ** (-decimals) + num_segm_tot = act_frac_sys.shape[0] + unique_nodes = np.unique(np.vstack((act_frac_sys[:, :2], act_frac_sys[:, 2:])), axis=0) + num_nodes_tot = unique_nodes.shape[0] + f = open(filename, "w+") + + # Note: always end statement in GMSH with ";" + # Note: comments in GMSH are as in C(++) "//" + f.write('// Geo file which meshes the input mesh from act_frac_sys.\n') + f.write('// Change mesh-elements size by varying "lc" below.\n\n') + + # Can specify the type of meshing algorithm for 2D meshing here: + # f.write('// Specify meshing algorithm:\n') + # f.write('-algo meshadapt;\n\n') + + # Set some parameters in the model: + f.write('lc = {:1.3f};\n'.format(char_len)) + f.write('lc_box = {:1.3f};\n'.format(char_len_boundary)) + f.write('height_res = {:4.3f};\n\n'.format(height_res)) + + # Allocate memory for points_created array and counters: + points_created = np.zeros((num_nodes_tot,), dtype=bool) + line_count = 0 + point_count = 0 + + for ii in act_frac_sys: + # Take two points per segment and write to .geo file: + # e.g. point: Point(1) = {.1, 0, 0, lc}; + nodes = np.zeros((2,), dtype=int) + nodes[0] = np.where(np.logical_and(ii[0] == unique_nodes[:, 0], ii[1] == unique_nodes[:, 1]))[0] + nodes[1] = np.where(np.logical_and(ii[2] == unique_nodes[:, 0], ii[3] == unique_nodes[:, 1]))[0] + + # Check if first point is already created, if not, add it: + if not points_created[nodes[0]]: + points_created[nodes[0]] = True + point_count += 1 + f.write('Point({:d}) = {{{:8.5f}, {:8.5f}, 0, lc}};\n'.format(nodes[0] + 1, + unique_nodes[nodes[0], 0], + unique_nodes[nodes[0], 1])) + + if not points_created[nodes[1]]: + points_created[nodes[1]] = True + point_count += 1 + f.write('Point({:d}) = {{{:8.5f}, {:8.5f}, 0, lc}};\n'.format(nodes[1] + 1, + unique_nodes[nodes[1], 0], + unique_nodes[nodes[1], 1])) + + line_count += 1 + f.write('Line({:d}) = {{{:d}, {:d}}};\n\n'.format(line_count, nodes[0] + 1, nodes[1] + 1)) + + # Store some internal variables for gmsh (used later after extrude): + f.write('num_points_frac = newp - 1;\n') + f.write('num_lines_frac = newl - 1;\n\n') + + # Write the box_data (box around fracture network in which we embed the fractures) + f.write('// Extra points for boundary of domain:\n') + for ii in range(4): + # For every corner of the box: + point_count += 1 + f.write('Point({:d}) = {{{:8.5f}, {:8.5f}, 0, lc_box}};\n'.format(point_count, + box_data[ii, 0], box_data[ii, 1])) + + # Add four lines for each side of the box: + f.write('\n// Extra lines for boundary of domain:\n') + line_count += 1 + f.write('Line({:d}) = {{{:d}, {:d}}};\n'.format(line_count, point_count - 3, point_count - 2)) + + line_count += 1 + f.write('Line({:d}) = {{{:d}, {:d}}};\n'.format(line_count, point_count - 2, point_count - 1)) + + line_count += 1 + f.write('Line({:d}) = {{{:d}, {:d}}};\n'.format(line_count, point_count - 1, point_count - 0)) + + line_count += 1 + f.write('Line({:d}) = {{{:d}, {:d}}};\n'.format(line_count, point_count - 0, point_count - 3)) + + # Make Curve loop for the boundary: + f.write('\n// Create line loop for boundary surface:\n') + f.write('Curve Loop(1) = {{{:d}, {:d}, {:d}, {:d}}};\n'.format(line_count - 3, line_count - 2, + line_count - 1, line_count)) + f.write('Plane Surface(1) = {1};\n\n') + f.write('Curve{1:num_lines_frac} In Surface{1};\n') + + # Extrude model to pseuo-3D: + f.write('\n// Extrude surface with embedded features\n') + f.write('Extrude {0, 0, height_res}{ Surface {1}; Layers{1}; Recombine;}\n') + f.write('Physical Volume("matrix", 9991) = {1};\n') + + f.write('num_surfaces_before = news;\n') + # f.write('Extrude {0, 0, height_res}{ Line {1:num_lines_frac}; Layers{1}; Recombine;}\n') + f.write('num_surfaces_after = news - 1;\n') + f.write('num_surfaces_fracs = num_surfaces_after - num_surfaces_before;\n\n') + for ii in range(act_frac_sys.shape[0]): + # f.write('Physical Surface({:d}) = {{num_surfaces_before + {:d}}};\n'.format(90000 + ii, ii)) + f.write('Extrude {{0, 0, height_res}}{{ Line {{{:d}}}; Layers{{1}}; Recombine;}}\n'.format(ii + 1)) + if export_frac: + f.write('Physical Surface({:d}) = {{news - 1}};\n'.format(90000 + ii)) + + # Create mesh and perform coherency check: + f.write('Mesh 3; // Generalte 3D mesh\n') + f.write('Coherence Mesh; // Remove duplicate entities\n') + f.write('Mesh.MshFileVersion = 2.1;\n') + f.close() + return 0 diff --git a/main.py b/main.py index db8a535..873b716 100644 --- a/main.py +++ b/main.py @@ -1,50 +1,33 @@ import numpy as np from multiprocessing import freeze_support -from preprocessing import frac_preprocessing +from preprocessing_code import frac_preprocessing +import os def main(): # File names and directories: - output_dir = 'FracData_test_with_inter' - filename_base = 'FracData_test_with_inter' - frac_data_raw = np.genfromtxt(filename_base + '.txt') + output_dir = 'output_Whitby' + filename_base = 'Whitby' + frac_data_raw = np.genfromtxt(os.path.join('Datasets', 'Whitby', 'Whitby_with_intersections.txt')) + # frac_data_raw = np.genfromtxt(os.path.join('Datasets', 'Whitby', 'Whitby_raw.txt')) # Input parameters for cleaning procedure - char_len_vec = np.array([10, 20, 40]) # vector containing the desired accuracy at which to process the network [m] - # NOTE: if you do not want to process in a hierarchical sense, scaler (int) input will work too and will process - # to a single characteristic length - - angle_tol_straighten = 5 # tolerance for straightening fracture segments [degrees] + char_len = 16 # characteristic length for cleaning and mesh generation [m] + angle_tol_straighten = 7.5 # tolerance for straightening fracture segments [degrees] merge_threshold = 0.86 # tolerance for merging nodes in algebraic constraint, values on interval [0.5, 0.86] [-] angle_tol_remove_segm = np.arctan(0.35) * 180 / np.pi # tolerance for removing accute intersections, values on interval [15, 25] [degrees] decimals = 7 # in order to remove duplicates we need to have fixed number of decimals - reuse_clean = True # if True will re-use cleaned network as starting point for hierarchical cleaning - straighten_first = True # if True will straighten fractures first before removing accute angles - mesh_clean = True # need gmsh installed and callable from command line in order to mesh!!! - mesh_raw = True # need gmsh installed and callable from command line in order to mesh!!! - num_partition_x = 1 # number of partitions for parallel implementation of intersection finding algorithm - num_partition_y = 1 # " ... " - - # Reservoir properties: - height_res = 50 - origin_x = 0 - origin_y = 0 - margin = 25 - length_x = np.max(frac_data_raw[:, [0, 2]]) - length_y = np.max(frac_data_raw[:, [1, 3]]) - - # Boxdata is only used in meshing and generation of geo file - box_data = np.array([[origin_x - margin, origin_y - margin], - [length_x + margin, origin_y - margin], - [length_x + margin, length_y + margin], - [origin_x - margin, length_y + margin]]) - - frac_preprocessing(frac_data_raw=frac_data_raw, char_len_vec=char_len_vec, output_dir=output_dir, - filename_base=filename_base, height_res=height_res, box_data=box_data, - angle_tol_straighten=angle_tol_straighten, merge_threshold=merge_threshold, - angle_tol_remove_segm=angle_tol_remove_segm, decimals=decimals, reuse_clean=reuse_clean, - straighten_first=straighten_first, mesh_clean=mesh_clean, mesh_raw=mesh_raw, - num_partition_x=num_partition_x, num_partition_y=num_partition_y) + mesh_clean = False # need gmsh installed and callable from command line in order to mesh!!! + mesh_raw = False # need gmsh installed and callable from command line in order to mesh!!! + num_partition_x = 4 # number of partitions for parallel implementation of intersection finding algorithm + num_partition_y = 4 # " ... " + + frac_preprocessing(frac_data_raw, char_len, output_dir=output_dir, filename_base=filename_base, merge_threshold=merge_threshold, + height_res=50, angle_tol_small_intersect=angle_tol_remove_segm, apertures_raw=None, box_data=None, margin=25, + mesh_clean=mesh_clean, mesh_raw=mesh_raw, angle_tol_straighten=angle_tol_straighten, straighten_after_cln=True, decimals=decimals, + tolerance_zero=1e-10, tolerance_intersect=1e-10, calc_intersections_before=False, calc_intersections_after=False, + num_partition_x=num_partition_x, num_partition_y=num_partition_y, partition_fractures_in_segms=True, matrix_perm=1, correct_aperture=False, + small_angle_iter=2, char_len_mult=1, char_len_boundary=None, main_algo_iters=1) if __name__ == "__main__": diff --git a/preprocessing_code.py b/preprocessing_code.py new file mode 100644 index 0000000..d57fa4f --- /dev/null +++ b/preprocessing_code.py @@ -0,0 +1,316 @@ +from graph_code import Graph, create_geo_file +from multiprocessing import Process, freeze_support +import numpy as np +from calc_intersections_segm_parallel import calc_intersections_segm_parallel +import os +import sys + + +def frac_preprocessing(frac_data_raw, char_len, output_dir='', filename_base='output', merge_threshold=0.66, + height_res=50, angle_tol_small_intersect=20, apertures_raw=None, box_data=None, margin=25, mesh_clean=False, + mesh_raw=False, angle_tol_straighten=7, straighten_after_cln=False, decimals=5, tolerance_zero=1e-10, + tolerance_intersect=1e-10, calc_intersections_before=False, calc_intersections_after=True, num_partition_x=1, + num_partition_y=1, partition_fractures_in_segms=True, matrix_perm=1, correct_aperture=False, + small_angle_iter=0, char_len_mult=1, char_len_boundary=None, main_algo_iters=2): + """ + Main fracture preprocessing code, most arguments are optional, but can be tweaked. Please see: + doi.org/10.1002/essoar.10507519.1 for more explanation on the theory behind the code and some results. + :param frac_data_raw: input fracture network data in the form [[x1, y1, x2, y2], ..., [...]] [m] + :param char_len: minimum allowable distance between any two vertices in the domain (characteristic length) [m] + :param output_dir: directory of the output (cleaned) fracture network and potential .geo and .msh results + :param filename_base: base name used for all the output files + :param merge_threshold: coefficient chosen on the closed interval from [0.5, 0.86] + :param height_res: height of the resulting 1-layer 3D reservoir [m] + :param angle_tol_small_intersect: minimum angle in degrees between two intersecting fractures + :param apertures_raw: list of apertures for each fracture segment [m] + :param box_data: coordinates of the bounding box (in order, from bottom left --> bottom right --> top right --> top left) [m] + :param margin: margin around the fracture network, used in case no bounding box is given + :param mesh_clean: boolean, if True then will call GMSH to mesh the cleaned fracture network in a triangular mesh + :param mesh_raw: boolean, if True then will call GMSH to mesh the raw fracture network in a triangular mesh + :param angle_tol_straighten: allowable deviation from a straight line for straightening algorithm + :param straighten_after_cln: boolean, if True then straighten the fractures after cleaning + :param decimals: number of decimals used to round-off input data + :param tolerance_zero: anything below this threshold is considered absolute zero! + :param calc_intersections_before: boolean, if True calculates intersections between fractures segments before cleaning + :param calc_intersections_after: boolean, if True calculates intersections between fractures segments after cleaning + :param num_partition_x: number of partitions used in fracture intersection calculation in x-direction (parallel computing required) + :param num_partition_y: number of partitions used in fracture intersection calculation in y-direction (parallel computing required) + :param partition_fractures_in_segms: boolean, if True partitions the fracture network into smaller subsegments of length char_len + :param matrix_perm: matrix permeability [mD] + :param correct_aperture: boolean, if True apply aperture correction method + :param small_angle_iter: number of iterations with small-angle correction algorithm + :param char_len_mult: multiplier for mesh characteristic length + :param char_len_boundary: characteristic mesh length on the boundary (before multiplier) + :param main_algo_iters: number of times the main cleaning algorithm is run + :return: + """ + assert 0.5 <= merge_threshold <= 0.86, "Choose threshold on closed interval [0.5, 0.86]" + if apertures_raw is None: + apertures_raw = np.ones((frac_data_raw.shape[0],)) * 1e-4 + + if char_len_boundary is None: + char_len_boundary = char_len + + if isinstance(frac_data_raw, str): + frac_data_raw = np.genfromtxt(frac_data_raw) + assert frac_data_raw.shape[1] == 4, "Data in wrong format, need N rows x 4 columns" + + print('--------------------------------------') + print('START preprocessing fracture network') + tot_partitions = num_partition_x * num_partition_y + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + frac_data_raw = np.round(frac_data_raw * 10 ** decimals) * 10 ** (-decimals) + print('Remove segments of zero length and duplicate segments') + frac_data_raw, apertures_raw = extract_unique_segms(frac_data_raw, apertures_raw) + len_raw_sys = np.sqrt((frac_data_raw[:, 0] - frac_data_raw[:, 2]) ** 2 + + (frac_data_raw[:, 1] - frac_data_raw[:, 3]) ** 2) + + print('Number of fracture segments: {:}'.format(frac_data_raw.shape[0])) + print('Min fracture segment length: {:}'.format(np.min(len_raw_sys))) + print('Max fracture segment length: {:}'.format(np.max(len_raw_sys))) + print('Mean fracture segment length: {:}'.format(np.mean(len_raw_sys))) + print('Cleaning length(s): {:}\n'.format(char_len)) + # -------------------------------------------------------------------------- + + act_frac_sys = frac_data_raw + apertures = apertures_raw + if calc_intersections_before: + print('START calculating initial intersections raw input fracture network') + print('\tNOTE: unoptimized!, can take long for very large networks') + # First find all intersections: + system_out_par, frac_order_vec_par, partition_lines = \ + calc_intersections_segm_parallel(frac_data_raw, apertures_raw, tolerance_intersect, + tolerance_zero, num_partition_x, num_partition_y) + + # Stack output from all domain together: + act_frac_sys = system_out_par[0] + apertures = frac_order_vec_par[0] + for ii in range(1, tot_partitions): + act_frac_sys = np.vstack((act_frac_sys, system_out_par[ii])) + apertures = np.hstack((apertures, frac_order_vec_par[ii])) + + print('DONE calculating initial intersections raw input fracture network') + if act_frac_sys.shape != frac_data_raw.shape: + num_intersections = int((act_frac_sys.shape[0] - frac_data_raw.shape[0]) / 2) + print('\tFound {:} intersections in raw input fracture network\n'.format(num_intersections)) + else: + print('\tNo intersections found in raw input fracture network\n') + + print('Remove duplicated segments\n') + act_frac_sys, apertures = extract_unique_segms(act_frac_sys, apertures) + act_frac_sys_cln = act_frac_sys + apertures_cln = apertures + act_frac_sys_raw = np.array(act_frac_sys, copy=True) + apertures_raw = np.array(apertures_cln, copy=True) + + if np.unique(apertures_cln).size == 1: + order_cleaning_segms = np.argsort(-np.sqrt((act_frac_sys[:, 0] - act_frac_sys[:, 2])**2 + (act_frac_sys[:, 1] - act_frac_sys[:, 3])**2)) + else: + order_cleaning_segms = np.argsort(-apertures_cln) + + if partition_fractures_in_segms: + act_frac_sys_cln, order_cleaning_segms, apertures_cln = segment_fractures(act_frac_sys=act_frac_sys_cln, + char_len=char_len, + order_discr_segms=order_cleaning_segms, + decimals=decimals, + aperture=apertures_cln) + + # -------------------------------------------------------------------------- + print('START constructing graph') + my_graph = Graph(matrix_perm=matrix_perm) + my_graph.add_multiple_edges(act_frac_sys_cln) + my_graph.apertures[np.where(my_graph.active_edges[:my_graph.get_num_edges()] == True)[0]] = apertures_cln + print('DONE constructing graph\n') + + print('START main cleaning loop for l_f={:}'.format(char_len)) + # print('\tNOTE: unoptimized!, can take long for very large networks or very small l_f') + for ii in range(main_algo_iters): + my_graph.simplify_graph(order_discr=order_cleaning_segms, char_len=char_len, merge_treshold=merge_threshold, + tolerange_small_angle=angle_tol_small_intersect, small_angle_iter=small_angle_iter, + tolerange_straight_angle=angle_tol_straighten, + correct_aperture=correct_aperture, straighten_edges=straighten_after_cln) + print('DONE main cleaning loop for l_f={:}\n'.format(char_len)) + + active_edges = np.where(my_graph.active_edges[:my_graph.get_num_edges()] == True)[0] + num_act_frac = len(active_edges) + act_frac_sys_cln = np.zeros((num_act_frac, 4)) + act_frac_sys_cln[:, 0] = my_graph.vertices[my_graph.edge_to_vertex[my_graph.active_edges, 0], 0] + act_frac_sys_cln[:, 1] = my_graph.vertices[my_graph.edge_to_vertex[my_graph.active_edges, 0], 1] + act_frac_sys_cln[:, 2] = my_graph.vertices[my_graph.edge_to_vertex[my_graph.active_edges, 1], 0] + act_frac_sys_cln[:, 3] = my_graph.vertices[my_graph.edge_to_vertex[my_graph.active_edges, 1], 1] + order_segms_after_cleaning = order_cleaning_segms[active_edges] + apertures_cln = my_graph.apertures[active_edges] + + # -------------------------------------------------------------------------- + if calc_intersections_after: + print('START calculating intersections clean fracture network for coherent mesh') + print('\tNOTE: unoptimized!, can take long for very large networks') + # First find all intersections: + num_part_x = 1 + num_part_y = 1 + num_part_tot = num_part_x * num_part_y + system_out_par, frac_order_vec_par, partition_lines = \ + calc_intersections_segm_parallel(act_frac_sys_cln, apertures_cln, tolerance_intersect, + tolerance_zero, number_partitions_x=num_part_x, number_partitions_y=num_part_y) + + # Stack output from all domain together: + act_frac_sys_cln = system_out_par[0] + apertures_cln = frac_order_vec_par[0] + for ii in range(1, num_part_tot): + act_frac_sys_cln = np.vstack((act_frac_sys_cln, system_out_par[ii])) + apertures_cln = np.hstack((apertures_cln, frac_order_vec_par[ii])) + print('DONE calculating intersections clean fracture network for coherent mesh\n') + + # -------------------------------------------------------------------------- + if box_data is None: + x_min = np.min(np.min(frac_data_raw[:, [0, 2]])) - margin + y_min = np.min(np.min(frac_data_raw[:, [1, 3]])) - margin + x_max = np.max(np.max(frac_data_raw[:, [0, 2]])) + margin + y_max = np.max(np.max(frac_data_raw[:, [1, 3]])) + margin + box_data = np.array([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max]]) + + # Filenames for meshing: + filename_geo_cln = os.path.join(output_dir, filename_base + '_mergefac_' + str(merge_threshold) + '_clean_lc_' + str( + char_len) + '.geo') + filename_out_cln = os.path.join(output_dir, filename_base + '_mergefac_' + str(merge_threshold) + '_clean_lc_' + str( + char_len) + '.msh') + filename_geo_raw = os.path.join(output_dir, filename_base + '_raw_lc_' + str(char_len) + '.geo') + filename_out_raw = os.path.join(output_dir, filename_base + '_raw_lc_' + str(char_len) + '.msh') + + # -------------------------------------------------------------------------- + print('START writing clean fracture system to file') + filename_clean = os.path.join(output_dir, filename_base + '_mergefac_' + str(merge_threshold) + '_clean_lc_' + str( + char_len) + '_fracsys.txt') + f = open(filename_clean, "w+") + for frac in act_frac_sys_cln: + f.write('{:9.5f} {:9.5f} {:9.5f} {:9.5f}\n'.format(frac[0], frac[1], frac[2], frac[3])) + f.close() + + filename_aper_clean = os.path.join(output_dir, filename_base + '_mergefac_' + str(merge_threshold) + '_clean_lc_' + str( + char_len) + '_aperture.txt') + f = open(filename_aper_clean, "w+") + for aper in apertures_cln: + f.write('{:16.15f} \n'.format(aper)) + f.close() + print('DONE writing clean fracture system to file\n') + + print('START creating geo-file for cleaned network (input for gmsh)') + create_geo_file(act_frac_sys=act_frac_sys_cln, filename=filename_geo_cln, decimals=decimals, + height_res=height_res, char_len=char_len * char_len_mult, box_data=box_data, + char_len_boundary=char_len_boundary * char_len_mult) + print('DONE creating geo-file for cleaned network (input for gmsh)\n') + + if mesh_clean: + print('START meshing cleaned network') + print('\tNOTE: In gmsh you need to have under Options -> Geometry -> General -> uncheck "Remove duplicate ..." otherwise meshing will crash/take too long') + os.system("gmsh {:s} -o {:s} -save".format(filename_geo_cln, filename_out_cln)) + print('DONE meshing cleaned network\n') + + # -------------------------------------------------------------------------- + print('START writing raw fracture system to file') + filename_raw = os.path.join(output_dir, filename_base + '_raw_lc_' + str(char_len) + '_fracsys.txt') + f = open(filename_raw, "w+") + for frac in act_frac_sys_raw: + f.write('{:9.5f} {:9.5f} {:9.5f} {:9.5f}\n'.format(frac[0], frac[1], frac[2], frac[3])) + f.close() + + filename_aper_raw = os.path.join(output_dir, filename_base + '_raw_lc_' + str(char_len) + '_aperture.txt') + f = open(filename_aper_clean, "w+") + for aper in apertures_raw: + f.write('{:16.15f} \n'.format(aper)) + f.close() + print('DONE writing raw fracture system to file\n') + + print('START creating geo-file for raw network (input for gmsh)') + create_geo_file(act_frac_sys=act_frac_sys_raw, filename=filename_geo_raw, decimals=decimals, + height_res=height_res, char_len=char_len * char_len_mult, box_data=box_data, + char_len_boundary=char_len_boundary * char_len_mult) + print('DONE creating geo-file for cleaned network (input for gmsh)\n') + + if mesh_raw: + print('START meshing raw network') + print( + '\tNOTE: In gmsh you need to have under Options -> Geometry -> General -> uncheck "Remove duplicate ..." otherwise meshing will crash/take too long') + os.system("gmsh {:s} -o {:s} -save".format(filename_geo_raw, filename_out_raw)) + print('DONE meshing raw network\n') + + print('Preprocessing succesfully finished') + print('-----------------------------------') + return 0 + + +def extract_unique_segms(act_frac_sys, apertures, remove_small_segms=True, tolerance_zero=1e-10): + """ + Extracts unique fracture segments and apertures + :param act_frac_sys: input fracture network data in the form [[x1, y1, x2, y2], ..., [...]] [m] + :param apertures: list of apertures for each fracture segment [m] + :param remove_small_segms: boolean, if True also removes zero-length fracture segments + :param tolerance_zero: any value below this will be interpreted as zero! + :return: unique_segms, unique_apertures + """ + # Extract unique indices: + dummy_sys = np.array(act_frac_sys, copy=True) + indices = np.where(act_frac_sys[:, 0] > act_frac_sys[:, 2])[0] + dummy_sys[indices, 0:2] = act_frac_sys[indices, 2:] + dummy_sys[indices, 2:] = act_frac_sys[indices, :2] + dummy_sys, unique_indices = np.unique(dummy_sys, axis=0, return_index=True) + dummy_aper = apertures[unique_indices] + + if remove_small_segms: + len_segm = np.sqrt((dummy_sys[:, 0] - dummy_sys[:, 2])**2 + (dummy_sys[:, 1] - dummy_sys[:, 3])**2) + dummy_sys = dummy_sys[len_segm > tolerance_zero, :] + dummy_aper = dummy_aper[len_segm > tolerance_zero] + return dummy_sys, dummy_aper + + +def segment_fractures(act_frac_sys, char_len, order_discr_segms, decimals=5, aperture=None): + """ + Perform partitioning into smaller subsegments (around lc size) + :param act_frac_sys: input fracture network data in the form [[x1, y1, x2, y2], ..., [...]] [m] + :param char_len: length of segments after partitioning [m] + :param order_discr_segms: order in which fracture segments are partitioned and stored in final array + :param decimals: number of decimals used in rounding off + :param aperture: array of apertures for each fracture segment + :return: segmented_fracture_system, new_order_array, new_aperture_array + """ + length_segms = np.sqrt((act_frac_sys[:, 0] - act_frac_sys[:, 2]) ** 2 + + (act_frac_sys[:, 1] - act_frac_sys[:, 3]) ** 2) + num_new_segms = int(np.sum(np.round(length_segms / char_len)) + np.sum(np.round(length_segms / char_len) == 0)) + act_frac_sys_new = np.zeros((num_new_segms, 4)) + aperture_segm = np.zeros((num_new_segms,)) + order_discr_segms_new = np.zeros((num_new_segms,), dtype=int) + ith_segm = 0 + + for ii in order_discr_segms: + size_segm = int(max(1, np.round(length_segms[ii] / char_len))) + id_vec = np.arange(0, size_segm) + act_frac_sys_new[ith_segm:(ith_segm + size_segm), 0] = act_frac_sys[ii, 0] + id_vec / size_segm * ( + act_frac_sys[ii, 2] - act_frac_sys[ii, 0]) + act_frac_sys_new[ith_segm:(ith_segm + size_segm), 1] = act_frac_sys[ii, 1] + id_vec / size_segm * ( + act_frac_sys[ii, 3] - act_frac_sys[ii, 1]) + act_frac_sys_new[ith_segm:(ith_segm + size_segm), 2] = act_frac_sys[ii, 0] + (id_vec + 1) / size_segm * ( + act_frac_sys[ii, 2] - act_frac_sys[ii, 0]) + act_frac_sys_new[ith_segm:(ith_segm + size_segm), 3] = act_frac_sys[ii, 1] + (id_vec + 1) / size_segm * ( + act_frac_sys[ii, 3] - act_frac_sys[ii, 1]) + + if aperture is not None: + aperture_segm[ith_segm:(ith_segm + size_segm)] = aperture[ii] + + # Update order of fractures: + order_discr_segms_new[ith_segm:(ith_segm + size_segm)] = np.arange(0, size_segm) + ith_segm + + ith_segm += size_segm + + act_frac_sys_new = np.round(act_frac_sys_new * 10 ** decimals) * 10 ** (-decimals) + if aperture is not None: + return act_frac_sys_new, order_discr_segms_new, aperture_segm + else: + return act_frac_sys_new, order_discr_segms_new + + +if __name__ == "__main__": + freeze_support() + frac_preprocessing(sys.argv[1], sys.argv[2]) diff --git a/requirements.txt b/requirements.txt index 6bad103..c4b83f7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,4 @@ numpy scipy +igraph +matplotlib