diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py index c5328f591..a05ebd925 100644 --- a/src/porepy/__init__.py +++ b/src/porepy/__init__.py @@ -242,7 +242,7 @@ # Modules from porepy.fracs import utils as frac_utils from porepy.fracs import meshing, fracture_importer -from porepy.grids import coarsening, partition, refinement +from porepy.grids import partition, refinement from porepy.numerics import displacement_correlation # Applications diff --git a/src/porepy/applications/md_grids/mdg_library.py b/src/porepy/applications/md_grids/mdg_library.py index 62164455f..0e94e3b96 100644 --- a/src/porepy/applications/md_grids/mdg_library.py +++ b/src/porepy/applications/md_grids/mdg_library.py @@ -249,17 +249,17 @@ def seven_fractures_one_L_intersection( def benchmark_regular_2d( - meshing_args: dict, is_coarse: bool = False, **meshing_kwargs + meshing_args: dict, **meshing_kwargs ) -> tuple[pp.MixedDimensionalGrid, FractureNetwork2d]: """ - Create a grid bucket for a domain containing the network introduced as example 2 of - Berre et al. 2018: Benchmarks for single-phase flow in fractured porous media. + Create a MixedDimensionalGrid for a domain containing the network introduced as + example 2 of Berre et al. 2018: Benchmarks for single-phase flow in fractured porous + media. Parameters: meshing_args: Dictionary containing at least "mesh_size_frac". If the optional values of "mesh_size_bound" and "mesh_size_min" are not provided, these are set by utils.set_mesh_sizes. - is_coarse: If True, coarsen the grid by volume. **meshing_kwargs: Keyword arguments for meshing as used by :meth:`~porepy.grids.mdg_generation.create_mdg`. @@ -281,8 +281,6 @@ def benchmark_regular_2d( ) mdg = pp.create_mdg("simplex", meshing_args, fracture_network, **meshing_kwargs) - if is_coarse: - pp.coarsening.coarsen(mdg, "by_volume") return mdg, fracture_network diff --git a/src/porepy/fracs/fracture_network_3d.py b/src/porepy/fracs/fracture_network_3d.py index 2d41f1d25..b67c71c9b 100644 --- a/src/porepy/fracs/fracture_network_3d.py +++ b/src/porepy/fracs/fracture_network_3d.py @@ -151,27 +151,6 @@ def __init__( """ - def add(self, fracture: pp.PlaneFracture) -> None: - """Add a fracture to the network. - - The fracture will be assigned a new index, higher than the maximum value - currently found in the network. - - Parameters: - fracture: Plane fracture to be added. - - """ - msg = "This functionality is deprecated and will be removed in a future version" - warnings.warn(msg, DeprecationWarning) - - ind = np.array([f.index for f in self.fractures]) - - if ind.size > 0: - fracture.set_index(np.max(ind) + 1) - else: - fracture.set_index(0) - self.fractures.append(fracture) - def copy(self) -> FractureNetwork3d: """Create a deep copy of the fracture network. @@ -1291,82 +1270,6 @@ def _remove_edge_intersections( all_p, edges, edges_2_frac, is_boundary_edge ) - def fractures_of_points(self, pts: np.ndarray) -> list[np.ndarray]: - """For a given point, find all fractures that refer to it. - - The point can be either a vertex or an internal point. - - Parameters: - pts: ``shape=(3, num_points)`` - - Coordinates of the points for which fractures should be found. - - Returns: - List of numpy arrays of integers. Each item, contains the indices of the - fractures associated to each point. - - """ - msg = "This functionality is deprecated and will be removed in a future version" - warnings.warn(msg, DeprecationWarning) - - fracs_of_points = [] - pts = np.atleast_1d(np.asarray(pts)) - for i in pts: - fracs_loc = [] - - # First identify edges that refer to the point - edge_ind = np.argwhere( - np.any(self.decomposition["edges"][:2] == i, axis=0) - ).ravel("F") - edges_loc = self.decomposition["edges"][:, edge_ind] - # Loop over all polygons. If their edges are found in edges_loc, - # store the corresponding fracture index - for poly_ind, poly in enumerate(self.decomposition["polygons"]): - ismem, _ = pp.array_operations.ismember_columns(edges_loc, poly) - if any(ismem): - fracs_loc.append(self.decomposition["polygon_frac"][poly_ind]) - fracs_of_points.append(list(np.unique(fracs_loc))) - - # Prepare for returning - fracs_of_points_out = [np.asarray(item) for item in fracs_of_points] - - return fracs_of_points_out - - def close_points(self, dist: float) -> list[tuple[int, int, float]]: - """Find pairs that are closer than the specified distance. - - In the set of points used to describe the fractures (after decomposition), - find pairs that are closer than a certain distance. - - This function is intended for debugging purposes. - - Parameters: - dist: Threshold distance. All points closer than this will be reported. - - Returns: - List of tuples, where each tuple contains indices of a set of close - points, and the distance between the points. - - The list is not symmetric. - If ``(a, b)`` is a member, ``(b, a)`` will not be. - - """ - msg = "This functionality is deprecated and will be removed in a future version" - warnings.warn(msg, DeprecationWarning) - - c_points = [] - - pt = self.decomposition["points"] - for pi in range(pt.shape[1]): - d = pp.distances.point_pointset(pt[:, pi], pt[:, pi + 1 :]) - ind = np.argwhere(d < dist).ravel("F") - for i in ind: - # Indices of close points, with an offset to compensate for slicing - # of the point cloud. - c_points.append((pi, i + pi + 1, d[i])) - - return c_points - def _verify_fractures_in_plane( self, p: np.ndarray, diff --git a/src/porepy/grids/coarsening.py b/src/porepy/grids/coarsening.py deleted file mode 100644 index 6ba66c47f..000000000 --- a/src/porepy/grids/coarsening.py +++ /dev/null @@ -1,806 +0,0 @@ -"""This module contains methods to coarsen grids. - -The main function is :func:`~porepy.grids.coarsening.coarsen` -(see there for more information). -""" - -from __future__ import annotations - -import warnings -from typing import Any, Optional, Union - -import numpy as np -import scipy.sparse as sps - -import porepy as pp -from porepy.grids import grid -from porepy.numerics.linalg.matrix_operations import sparse_array_to_row_col_data -from porepy.utils import grid_utils, tags - - -def coarsen( - grid: Union[pp.Grid, pp.MixedDimensionalGrid], method: str, **method_kwargs -) -> None: - """Create a coarse grid from a given grid. - - If a md-grid is passed, the procedure is applied to the grid of highest dimension. - - Notes: - - - The passed grid is modified. - - Do not call :meth:`~porepy.grids.grid.Grid.compute_geometry` afterwards. - - Parameters: - grid: The grid or mixed-dimensional grid to be coarsened. - method: A string defining the coarsening method. - - The available options are: - - - ``'by_volume'``: The coarsening is based on the cell volume. - - ``'by_tpfa'``: Uses the AMG method's coarse/fine-splittings based on - direct couplings - - method_kwargs: Arguments for each ``method``. - - For the ``'by_volume'``- method, see :func:`create_aggregations` for an - overview on admissible keyword arguments. - - For the ``'by_tpfa'``- method, see :func:`create_partition`. - Additionally, a keyword argument ``if_seeds`` of boolean type is supported. - If ``True``, :func:`generate_seeds` is called to create seeds. - - Raises: - ValueError: If ``method`` is not among the supported options. - - """ - msg = "This functionality is deprecated and will be removed in a future version" - warnings.warn(msg, DeprecationWarning) - - if method.lower() == "by_volume": - partition = create_aggregations(grid, **method_kwargs) - - elif method.lower() == "by_tpfa": - seeds = np.empty(0, dtype=int) - if method_kwargs.get("if_seeds", False): - seeds = generate_seeds(grid) - matrix = _tpfa_matrix(grid) - partition = create_partition(matrix, grid, seeds=seeds, **method_kwargs) # type: ignore - - else: - raise ValueError(f"Undefined method `{method}` for coarsening algorithm.") - - generate_coarse_grid(grid, partition) # type: ignore - - -def generate_coarse_grid( - grid: Union[pp.Grid, pp.MixedDimensionalGrid], - subdiv: Union[np.ndarray, dict[pp.Grid, tuple[Any, np.ndarray]]], -) -> None: - """Generates a coarse grid by clustering the cells according to the flags - given by ``subdiv``. - - ``subdiv`` must be as long as the number of cells in the original grid, it contains - integers (possibly not continuous) which represent the cell IDs in the final, - coarser mesh. I.e. it is a cell map from finer to coarser. - - If ``grid`` is a mixed-dimensional grid, the coarsening is applied to the higher - dimensional grid. - - Warning: - This method effectively overwrites every grid property computed by - :meth:`~porepy.grids.grid.Grid.compute_geometry`. Do not call that method after - calling this one. - - Note: - There is no check for disconnected cells in the final grid. - - Example: - - >>> grid = ... # some grid with 12 cells - >>> subdiv = np.array([0,0,1,1,2,2,3,3,4,4,5,5]) # coarser grid with 6 cells - >>> generate_coarse_grid(grid, subdiv) - - Parameters: - grid: A grid or mixed-dimensional grid. - subdiv: If ``grid`` is a single grid, a single array-like object as in above - example suffices. - - If ``grid`` is a mixed-dimensional grid, a dictionary containing per grid - (key) a 2-tuple, where the second entry is the partition map as seen above. - This special structure is passed by :func:`coarsen`. - - """ - msg = "This functionality is deprecated and will be removed in a future version" - warnings.warn(msg, DeprecationWarning) - - if isinstance(grid, pp.Grid): - if isinstance(subdiv, dict): - # If the subdiv is a dictionary with grid as a key (this can happen if we - # are forwarded here from coarsen), the input must be simplified. - subdiv = subdiv[grid][1] - _generate_coarse_grid_single(grid, subdiv, False) - - if isinstance(grid, pp.MixedDimensionalGrid): - _generate_coarse_grid_mdg(grid, subdiv) - - -def reorder_partition( - subdiv: Union[np.ndarray, dict[Any, tuple[Any, np.ndarray]]], -) -> Union[np.ndarray, dict[Any, tuple[Any, np.ndarray]]]: - """Re-order the partition IDs in order to obtain contiguous numbers. - - Parameters: - subdiv: A subdivision/partition as an array containing an ID for each cell, or - a dictionary containing the previous in a 2-tuple for any key. - - Return: - The subdivision stored in a contiguous way. - - """ - if isinstance(subdiv, dict): - for _, (_, partition) in subdiv.items(): - old_ids = np.unique(partition) - for new_id, old_id in enumerate(old_ids): - partition[partition == old_id] = new_id - else: - old_ids = np.unique(subdiv) - for new_id, old_id in enumerate(old_ids): - subdiv[subdiv == old_id] = new_id - - return subdiv - - -def generate_seeds(mdg: Union[pp.Grid, pp.MixedDimensionalGrid]) -> np.ndarray: - """Generates a priory seeds (cells) for coarsening a mixed-dimensional grid based - on topological information about the highest-dimensional grid. - - Parameters: - mdg: A grid or mixed-dimensional grid. - - Returns: - If ``mdg`` is a single grid, this function returns an empty array. - - If ``mdg`` is a mixed-dimensional grid, this function returns an initial seed - for the coarsening based on the mortar projections between the grid of highest - dimension and grids of co-dimension 1. - - """ - seeds = np.empty(0, dtype=int) - - if isinstance(mdg, grid.Grid): - return seeds - - # Extract the higher dimensional grid - g_h = mdg.subdomains(dim=mdg.dim_max())[0] - g_h_faces, g_h_cells, _ = sparse_array_to_row_col_data(g_h.cell_faces) - - # Extract the 1-codimensional grids - gs = mdg.subdomains(dim=mdg.dim_max() - 1) - - for g in gs: - tips = np.where(g.tags["tip_faces"])[0] - faces, cells, _ = sparse_array_to_row_col_data(g.cell_faces) - index = np.isin(faces, tips).nonzero()[0] - cells = np.unique(cells[index]) - - # recover the mapping between the secondary and the primary grid - mg = mdg.subdomain_pair_to_interface((g_h, g)) - primary_to_mortar = mg.primary_to_mortar_int() - secondary_to_mortar = mg.secondary_to_mortar_int() - # this is the old face_cells mapping - face_cells = secondary_to_mortar.T * primary_to_mortar - - interf_cells, interf_faces, _ = sparse_array_to_row_col_data(face_cells) - index = np.isin(interf_cells, cells).nonzero()[0] - index = np.isin(g_h_faces, interf_faces[index]).nonzero()[0] - seeds = np.concatenate((seeds, g_h_cells[index])) - - return seeds - - -def create_aggregations( - grid: Union[pp.Grid, pp.MixedDimensionalGrid], **kwargs -) -> dict[pp.Grid, tuple[pp.Grid, np.ndarray]]: - """Creates a cell partition based on their volumes. - - Parameters: - grid: A single grid or mixed-dimensional grid. - **kwargs: Following keyword arguments are supported: - - - ``'weight'``: A float serving as weight for the mean of the cell volumes. - Defaults to 1. - - Returns: - A dictionary containing a partition per grid. - - """ - # Extract the higher dimensional grids and store in a list - if isinstance(grid, pp.MixedDimensionalGrid): - grid_list: list[pp.Grid] = grid.subdomains(dim=grid.dim_max()) - elif isinstance(grid, pp.Grid): - grid_list = [grid] - else: - raise ValueError("Only subdomains and mixed-dimensional grids supported.") - - partition: dict[pp.Grid, tuple[pp.Grid, np.ndarray]] = {} - - for g in grid_list: - partition_local = -np.ones(g.num_cells, dtype=int) - - volumes = g.cell_volumes.copy() - volumes_checked = volumes.copy() - c2c = g.cell_connection_map() - - # Compute the inverse of the harmonic mean - weight = kwargs.get("weight", 1.0) - mean = weight * np.mean(volumes) - - new_id = 1 - while np.any(partition_local < 0): - # Consider the smallest element to be aggregated - cell_id = np.argmin(volumes_checked) - - # If the smallest fulfil the condition, stop the loop - if volumes[cell_id] > mean: - break - - do_it = True - old_cluster = np.array([cell_id]) - while do_it: - cluster = __get_neigh(old_cluster, c2c, partition_local) - volume = np.sum(volumes[cluster]) - if volume > mean or np.array_equal(old_cluster, cluster): - do_it = False - else: - old_cluster = cluster - - # If one of the element in the cluster has already a partition id, - # we uniform the ids - partition_cluster = partition_local[cluster] - has_coarse_id = partition_cluster > 0 - if np.any(has_coarse_id): - # For each coarse id in the cluster, rename the coarse ids in the - # partition_local - for partition_id in partition_cluster[has_coarse_id]: - which_partition_id = partition_local == partition_id - partition_local[which_partition_id] = new_id - volumes[which_partition_id] = volume - volumes_checked[which_partition_id] = volume - - # Update the data for the cluster - partition_local[cluster] = new_id - volumes[cluster] = volume - new_id += 1 - - volumes_checked[cluster] = np.inf - - volumes_checked = volumes.copy() - which_cell = volumes_checked < mean - volumes_checked[np.logical_not(which_cell)] = np.inf - - while np.any(which_cell): - cell_id = np.argmin(volumes_checked) - part_cell = partition_local[cell_id] - # Extract the neighbors of the current cell - loc = slice(c2c.indptr[cell_id], c2c.indptr[cell_id + 1]) - neighbors = np.setdiff1d(c2c.indices[loc], np.asarray(cell_id)) - part_neighbors = partition_local[neighbors] - neighbors = neighbors[part_neighbors != part_cell] - if neighbors.size == 0: - volumes_checked[:] = np.inf - which_cell = volumes_checked < mean - continue - smallest = np.argmin(volumes[neighbors]) - mask = partition_local == part_cell - partition_local[mask] = partition_local[neighbors[smallest]] - volumes[mask] = volumes[smallest] + volumes[cell_id] - volumes_checked[mask] = volumes[smallest] + volumes[cell_id] - which_cell = volumes_checked < mean - - # Fill up the cells which are left - has_not_coarse_id = partition_local < 0 - partition_local[has_not_coarse_id] = new_id + np.arange( - np.sum(has_not_coarse_id) - ) - - partition[g] = (g.copy(), partition_local) - - return partition - - -def create_partition( - A: sps.spmatrix, - g: Union[pp.Grid, pp.MixedDimensionalGrid], - seeds: Optional[np.ndarray] = None, - **kwargs, -) -> dict[pp.Grid, tuple[pp.Grid, np.ndarray]]: - """Create the partition based on an input matrix using the AMG - method's coarse/fine-splittings based on direct couplings. - - The standard values for ``cdepth`` and ``epsilon`` are taken from the reference - below. - - Example: - - >>> part = create_partition(tpfa_matrix(g)) - >>> g = generate_coarse_grid(g, part) - - References: - U. Trottenberg, C. W. Oosterlee, and A. Schuller (200): - Multigrid, Academic press. - - Parameters: - A: A sparse matrix used for the agglomeration. - g: A single grid or mixed-dimensional grid. - seeds: ``default=None`` - - A-priory defined cells of the coarser grid. - **kwargs: The following keyword arguments are supported: - - - ``'cdepth'``: A number to define the strength of the aggregation, i.e. a - a greater number results in lesser cells. Defaults to 2. - - ``'epsilon'``: A float representing the weight for the off-diagonal - entries to define the *strong negative coupling*. Defaults to 0.25. - - Returns: - A dictionary containing the a 2-tuple per grid. The 2-tuple contains the grid - with the highest dimension and the map from finer to coarser grid containing as - an array of agglomeration indices. - - If ``g`` is a single grid, the grid of highest dimension is ``g`` itself. - - """ - - cdepth = int(kwargs.get("cdepth", 2)) - epsilon = kwargs.get("epsilon", 0.25) - - # NOTE: Extract the higher dimensional grids, we suppose it is unique - if isinstance(g, pp.MixedDimensionalGrid): - g_high = g.subdomains(dim=g.dim_max())[0] - else: - g_high = g - - if A.size == 0: - return {g_high: (g_high.copy(), np.zeros(1))} - - Nc = A.shape[0] - - # For each node, which other nodes are strongly connected to it - ST = sps.lil_matrix((Nc, Nc), dtype=bool) - - # In the first instance, all cells are strongly connected to each other - At = A.T - - for i in np.arange(Nc): - loc = slice(At.indptr[i], At.indptr[i + 1]) - ci, vals = At.indices[loc], At.data[loc] - neg = vals < 0.0 - nvals = vals[neg] - nci = ci[neg] - minId = np.argmin(nvals) - ind = -nvals >= epsilon * np.abs(nvals[minId]) - ST[nci[ind], i] = True - - # Temporary field, will store connections of depth 1 - for _ in np.arange(2, cdepth + 1): - STold = ST.copy() - for j in np.arange(Nc): - rowj = np.array(STold.rows[j]) - if rowj.size == 0: - continue - row = np.hstack([STold.rows[r] for r in rowj]) - ST[j, np.concatenate((rowj, row))] = True - - ST.setdiag(False) - lmbda = np.array([len(s) for s in ST.rows]) - - # Define coarse nodes - candidate = np.ones(Nc, dtype=bool) - is_fine = np.zeros(Nc, dtype=bool) - is_coarse = np.zeros(Nc, dtype=bool) - - # cells that are not important for any other cells are on the fine scale. - for row_id, row in enumerate(ST.rows): - if not row: - is_fine[row_id] = True - candidate[row_id] = False - - ST = ST.tocsr() - it = 0 - while np.any(candidate): - i = np.argmax(lmbda) - is_coarse[i] = True - j = ST.indices[ST.indptr[i] : ST.indptr[i + 1]] - jf = j[candidate[j]] - is_fine[jf] = True - candidate[np.r_[i, jf]] = False - loop = ST.indices[ - pp.array_operations.expand_index_pointers(ST.indptr[jf], ST.indptr[jf + 1]) - ] - for row in np.unique(loop): - s = ST.indices[ST.indptr[row] : ST.indptr[row + 1]] - lmbda[row] = s[candidate[s]].size + 2 * s[is_fine[s]].size - lmbda[np.logical_not(candidate)] = -1 - it = it + 1 - - # Something went wrong during aggregation - assert it <= Nc - - del lmbda, ST - - if seeds is not None: - is_coarse[seeds] = True - is_fine[seeds] = False - - # If two neighbors are coarse, eliminate one of them without touching the - # seeds - c2c = np.abs(A) > 0 - c2c_rows, _, _ = sparse_array_to_row_col_data(c2c.transpose()) - - pairs = np.empty((0, 2), dtype=int) - for idx, it in enumerate(np.where(is_coarse)[0]): - loc = slice(c2c.indptr[it], c2c.indptr[it + 1]) - ind = np.setdiff1d(c2c_rows[loc], it) - cind = ind[is_coarse[ind]] - new_pair = np.stack((np.repeat(it, cind.size), cind), axis=-1) - pairs = np.append(pairs, new_pair, axis=0) - - # Remove one of the neighbors cells - if pairs.size: - pairs = np.unique( - np.sort(pairs, axis=1), axis=0, return_index=True, return_inverse=True - )[0] - for ij in pairs: - A_val = A[ij, ij].A.ravel() - ids = ij[np.argsort(A_val)] - if seeds is not None: - ids = np.setdiff1d(ids, seeds, assume_unique=True) - if ids.size: - is_coarse[ids[0]] = False - is_fine[ids[0]] = True - - coarse = np.where(is_coarse)[0] - - # Primal grid - NC = coarse.size - primal = sps.lil_matrix((NC, Nc), dtype=bool) - primal[np.arange(NC), coarse[np.arange(NC)]] = True - - connection = sps.lil_matrix((Nc, Nc), dtype=np.double) - for it in np.arange(Nc): - n = np.setdiff1d(c2c_rows[c2c.indptr[it] : c2c.indptr[it + 1]], it) - loc = slice(A.indptr[it], A.indptr[it + 1]) - A_idx, A_row = A.indices[loc], A.data[loc] - mask = A_idx != it - connection[it, n] = np.abs(A_row[mask] / A_row[np.logical_not(mask)]) - - connection = connection.tocsr() - - candidates_rep = np.ediff1d(connection.indptr) - candidates_idx = np.repeat(is_coarse, candidates_rep) - candidates = np.stack( - ( - connection.indices[candidates_idx], - np.repeat(np.arange(NC), candidates_rep[is_coarse]), - ), - axis=-1, - ) - - connection_idx = pp.array_operations.expand_index_pointers( - connection.indptr[coarse], connection.indptr[coarse + 1] - ) - - # The array "candidates" stores the nonzero rows and columns of the connection - # matrix. We utilize the COO sparse matrix, which has the same internal structure. - # Using the following constructor: coo_matrix((data, [row_idx, col_idx])) - # Then casting the format to CSR, which is conveniently used later. - vals = sps.coo_matrix( - (connection.data[connection_idx], candidates.T), shape=[Nc, NC] - ).tocsr() - - del candidates_rep, candidates_idx, connection_idx - - it = NC - not_found = np.logical_not(is_coarse) - # Process the strongest connection globally - while np.any(not_found): - np.argmax(vals.data) - vals.argmax(axis=0) - mcind = np.atleast_1d(np.squeeze(np.asarray(vals.argmax(axis=0)))) - mcval = -np.inf * np.ones(mcind.size) - for c, r in enumerate(mcind): - loc = slice(vals.indptr[r], vals.indptr[r + 1]) - vals_idx, vals_data = vals.indices[loc], vals.data[loc] - mask = vals_idx == c - if vals_idx.size == 0 or not np.any(mask): - continue - mcval[c] = vals_data[mask] - - mi = np.argmax(mcval) - nadd = mcind[mi] - - primal[mi, nadd] = True - it = it + 1 - if it > Nc + 5: - break - - not_found[nadd] = False - vals.data[vals.indptr[nadd] : vals.indptr[nadd + 1]] = 0 - - loc = slice(connection.indptr[nadd], connection.indptr[nadd + 1]) - nc = connection.indices[loc] - af = not_found[nc] - nc = nc[af] - nv = mcval[mi] * connection[nadd, :] - nv = nv.data[af] - if len(nc) > 0: - vals += sps.csr_matrix((nv, (nc, np.repeat(mi, len(nc)))), shape=(Nc, NC)) - - coarse, fine = primal.tocsr().nonzero() - partition = {g_high: (g_high.copy(), coarse[np.argsort(fine)])} - return partition - - -def _generate_coarse_grid_single( - g: pp.Grid, subdiv: np.ndarray, face_map: bool -) -> Union[np.ndarray, None]: - """Auxiliary function to create a coarsening for a given single grid. - - Parameters: - g: A single-dimensional grid. - subdiv: A coarsening map containing cell IDs of the coarser grid for each ID - in the finer grid - face_map: A bool indicating if the face map for the coarser grid should be - returned - - Returns: - If ``face_map`` is True, returns the cell-to-face map of the coarser grid. - - """ - - subdiv = np.asarray(subdiv) - assert subdiv.size == g.num_cells - - # declare the storage array to build the cell_faces map - cell_faces = np.empty(0, dtype=g.cell_faces.indptr.dtype) - cells = np.empty(0, dtype=cell_faces.dtype) - orient = np.empty(0, dtype=g.cell_faces.data.dtype) - - # declare the storage array to build the face_nodes map - face_nodes = np.empty(0, dtype=g.face_nodes.indptr.dtype) - nodes = np.empty(0, dtype=face_nodes.dtype) - visit = np.zeros(g.num_faces, dtype=bool) - - # compute the face_node indexes - num_nodes_per_face = g.face_nodes.indptr[1:] - g.face_nodes.indptr[:-1] - face_node_ind = pp.matrix_operations.rldecode( - np.arange(g.num_faces), num_nodes_per_face - ) - - cells_list = np.unique(subdiv) - cell_volumes = np.zeros(cells_list.size) - cell_centers = np.zeros((3, cells_list.size)) - - for cellId, cell in enumerate(cells_list): - # extract the cells of the original mesh associated to a specific label - cells_old = np.where(subdiv == cell)[0] - - # compute the volume - cell_volumes[cellId] = np.sum(g.cell_volumes[cells_old]) - cell_centers[:, cellId] = np.average(g.cell_centers[:, cells_old], axis=1) - - # reconstruct the cell_faces mapping - faces_old, _, orient_old = sparse_array_to_row_col_data( - g.cell_faces[:, cells_old] - ) - mask = np.ones(faces_old.size, dtype=bool) - mask[np.unique(faces_old, return_index=True)[1]] = False - # extract the indexes of the internal edges, to be discared - index = np.array( - [np.where(faces_old == f)[0] for f in faces_old[mask]], dtype=int - ).ravel() - faces_new = np.delete(faces_old, index) - cell_faces = np.r_[cell_faces, faces_new] - cells = np.r_[cells, np.repeat(cellId, faces_new.shape[0])] - orient = np.r_[orient, np.delete(orient_old, index)] - - # reconstruct the face_nodes mapping - # consider only the unvisited faces - not_visit = ~visit[faces_new] - if not_visit.size == 0 or np.all(~not_visit): - continue - # mask to consider only the external faces - mask = np.atleast_1d( - np.sum( - [face_node_ind == f for f in faces_new[not_visit]], - axis=0, - dtype=bool, - ) - ) - face_nodes = np.r_[face_nodes, face_node_ind[mask]] - - nodes_new = g.face_nodes.indices[mask] - nodes = np.r_[nodes, nodes_new] - visit[faces_new] = True - - # Rename the faces - cell_faces_unique = np.unique(cell_faces) - cell_faces_id = np.arange(cell_faces_unique.size, dtype=cell_faces.dtype) - cell_faces = np.array( - [cell_faces_id[np.where(cell_faces_unique == f)[0]] for f in cell_faces] - ).ravel() - shape = (cell_faces_unique.size, cells_list.size) - cell_faces = sps.csc_matrix((orient, (cell_faces, cells)), shape=shape) - - # Rename the nodes - face_nodes = np.array( - [cell_faces_id[np.where(cell_faces_unique == f)[0]] for f in face_nodes] - ).ravel() - nodes_list = np.unique(nodes) - nodes_id = np.arange(nodes_list.size, dtype=nodes.dtype) - nodes = np.array([nodes_id[np.where(nodes_list == n)[0]] for n in nodes]).ravel() - - # sort the nodes - nodes = nodes[np.argsort(face_nodes, kind="mergesort")] - data = np.ones(nodes.size, dtype=g.face_nodes.data.dtype) - indptr = np.r_[0, np.cumsum(np.bincount(face_nodes))] - face_nodes = sps.csc_matrix((data, nodes, indptr)) - - g.name = "Coarsened grid" - # store again the data in the same grid - g.history.append("coarse") - - g.nodes = g.nodes[:, nodes_list] - g.num_nodes = g.nodes.shape[1] - - g.face_nodes = face_nodes - g.num_faces = g.face_nodes.shape[1] - g.face_areas = g.face_areas[cell_faces_unique] - g.tags = tags.extract(g.tags, cell_faces_unique, tags.standard_face_tags()) - g.face_normals = g.face_normals[:, cell_faces_unique] - g.face_centers = g.face_centers[:, cell_faces_unique] - - g.cell_faces = cell_faces - g.num_cells = g.cell_faces.shape[1] - g.cell_volumes = cell_volumes - g.cell_centers = grid_utils.star_shape_cell_centers(g, as_nan=True) - is_nan = np.isnan(g.cell_centers[0, :]) - g.cell_centers[:, is_nan] = cell_centers[:, is_nan] - - if face_map: - return np.array([cell_faces_unique, cell_faces_id]) - else: - # Explicitly return None to make mypy happy. - return None - - -def _generate_coarse_grid_mdg( - mdg: pp.MixedDimensionalGrid, - subdiv: Union[np.ndarray, dict[pp.Grid, tuple[pp.Grid, np.ndarray]]], -): - """Auxiliary function to create a coarsening for a given mixed-dimensional grid. - - Parameters: - mdg: A mixed-dimensional grid. - subdiv: A subdivision containing the coarsening map for each grid in a 2-tuple. - - """ - - if not isinstance(subdiv, dict): - g = mdg.subdomains(dim=mdg.dim_max())[0] - subdiv = {g: subdiv} # type: ignore - - for g, (_, partition) in subdiv.items(): - # Construct the coarse grids - face_map = _generate_coarse_grid_single(g, partition, True) - assert face_map is not None # make mypy happy - - # Update all the primary_to_mortar_int for all the 'edges' connected to the grid - # We update also all the face_cells - for intf in mdg.subdomain_to_interfaces(g): - # The indices that need to be mapped to the new grid - - data = mdg.interface_data(intf) - - projections = [ - intf.primary_to_mortar_int().tocsr(), - intf.primary_to_mortar_avg().tocsr(), - ] - - for ind, mat in enumerate(projections): - indices = mat.indices - - # Map indices - mask = np.argsort(indices) - indices = np.isin(face_map[0, :], indices[mask]).nonzero()[0] - # Reverse the ordering - indices = indices[np.argsort(mask)] - - # Create the new matrix - shape = (mat.shape[0], g.num_faces) - projections[ind] = sps.csr_matrix( - (mat.data, indices, mat.indptr), shape=shape - ) - - # Update mortar projection - intf._primary_to_mortar_int = projections[0].tocsc() - intf._primary_to_mortar_avg = projections[1].tocsc() - - # Also update all other projections to primary - intf._set_projections(secondary=False) - - # update also the face_cells map - face_cells = data["face_cells"].tocsr() - indices = face_cells.indices - - # map indices - mask = np.argsort(indices) - indices = np.isin(face_map[0, :], indices[mask]).nonzero()[0] - face_cells.indices = indices[np.argsort(mask)] - - # update the map - data["face_cells"] = face_cells.tocsc() - - -def _tpfa_matrix( - g: Union[pp.Grid, pp.MixedDimensionalGrid], - perm: Optional[pp.SecondOrderTensor] = None, -) -> sps.spmatrix: - """Compute a two-point flux approximation for a given grid - - This is a helper method for :func:`create_partition`. - - Parameters: - g: A single grid or mixed-dimensional grid. - perm: ``default=None`` - - The permeability as a tensor. If not given, defaults to the unit tensor. - - Returns: - The TPFA matrix for given grid and permeability. - - """ - if isinstance(g, pp.MixedDimensionalGrid): - g = g.subdomains(dim=g.dim_max())[0] - - if perm is None: - perm = pp.SecondOrderTensor(np.ones(g.num_cells)) - - solver = pp.Tpfa("flow") - specified_parameters = { - "second_order_tensor": perm, - "bc": pp.BoundaryCondition(g, np.empty(0), ""), - } - data = pp.initialize_default_data(g, {}, "flow", specified_parameters) - solver.discretize(g, data) - flux, _ = solver.assemble_matrix_rhs(g, data) - return flux - - -def __get_neigh( - cells_id: np.ndarray, c2c: sps.spmatrix, partition: np.ndarray -) -> np.ndarray: - """An auxiliary function for :func:`create_aggregations` to get neighbouring cells. - - Parameters: - cells_id: An array containing cell IDs. - c2c: A sparse map between cells sharing a face. - partition: A partition map containing indices of the coarser grid for each cell - in the finer grid. - - Returns: - Neighbouring cell IDs in the ``partition``. - - """ - neighbors = np.empty(0, dtype=int) - - for cell_id in np.atleast_1d(cells_id): - # Extract the neighbors of the current cell - loc = slice(c2c.indptr[cell_id], c2c.indptr[cell_id + 1]) - neighbors = np.hstack((neighbors, c2c.indices[loc])) - - neighbors = np.unique(neighbors) - partition_neighbors = partition[neighbors] - - # Check if some neighbor has already a coarse id - return np.sort(neighbors[partition_neighbors < 0]) diff --git a/src/porepy/numerics/vem/dual_elliptic.py b/src/porepy/numerics/vem/dual_elliptic.py index a70e3a6a6..6e13d0372 100644 --- a/src/porepy/numerics/vem/dual_elliptic.py +++ b/src/porepy/numerics/vem/dual_elliptic.py @@ -378,305 +378,6 @@ def _velocity_dof( hat_E_int = sps.bmat([[U * hat_E_int], [sps.csr_matrix(shape)]]) return hat_E_int - def assemble_int_bound_flux( - self, - sd: pp.Grid, - data: dict, - intf: pp.MortarGrid, - data_edge: dict, - cc: np.ndarray, - matrix: np.ndarray, - rhs: np.ndarray, - self_ind: int, - use_secondary_proj: bool = False, - ) -> None: - """Abstract method. Assemble the contribution from an internal - boundary, manifested as a flux boundary condition. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose flux continuity on a higher-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g: Grid which the condition should be imposed on. - data: Data dictionary for the node in the mixed-dimensional grid. - data_edge: Data dictionary for the edge in the mixed-dimensional grid. - cc: (block matrix, 3x3) Block matrix for the coupling condition. - The first and second rows and columns are identified with the - primary and secondary side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix: (block matrix 3x3) Discretization matrix for the edge and the two - adjacent nodes. - rhs: (block_array 3x1) Right hand side contribution for the edge and the two - adjacent nodes. - self_ind: Index in cc and matrix associated with this node. Should be either - 1 or 2. - use_secondary_proj: If True, the secondary side projection operator is - used. Needed for periodic boundary conditions. - - """ - msg = """This function is deprecated and will be removed, most likely in the - second half of 2022. - - To assemble mixed-dimensional elliptic problems, the recommended solution is - either to use the models, or to use the automatic differentiation framework - directly. - """ - warn(msg, DeprecationWarning, stacklevel=2) - - # The matrix must be the VEM discretization matrix. - if use_secondary_proj: - proj = intf.mortar_to_secondary_int() - else: - proj = intf.mortar_to_primary_int() - - hat_E_int = self._velocity_dof(sd, intf, proj) - cc[self_ind, 2] += matrix[self_ind, self_ind] * hat_E_int - - def assemble_int_bound_source( - self, - sd: pp.Grid, - data: dict, - intf: pp.MortarGrid, - data_edge: dict, - cc: np.ndarray, - matrix: np.ndarray, - rhs: np.ndarray, - self_ind: int, - ) -> None: - """Abstract method. Assemble the contribution from an internal - boundary, manifested as a source term. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose flux continuity on a lower-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - primary and secondary side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - - """ - msg = """This function is deprecated and will be removed, most likely in the - second half of 2022. - - To assemble mixed-dimensional elliptic problems, the recommended solution is - either to use the models, or to use the automatic differentiation framework - directly. - """ - warn(msg, DeprecationWarning, stacklevel=2) - proj = intf.secondary_to_mortar_avg() - - A = proj.T - shape = (sd.num_faces, A.shape[1]) - cc[self_ind, 2] += sps.bmat([[sps.csr_matrix(shape)], [A]]) - - def assemble_int_bound_pressure_trace( - self, - sd: pp.Grid, - data: dict, - intf: pp.MortarGrid, - data_edge: dict, - cc: Optional[np.ndarray], - matrix: np.ndarray, - rhs: np.ndarray, - self_ind: int, - use_secondary_proj: bool = False, - assemble_matrix=True, - assemble_rhs=True, - ) -> None: - """Abstract method. Assemble the contribution from an internal - boundary, manifested as a condition on the boundary pressure. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose flux continuity on a higher-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g: Grid which the condition should be imposed on. - data: Data dictionary for the node in the mixed-dimensional grid. - data_edge: Data dictionary for the edge in the mixed-dimensional grid. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - primary and secondary side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind: Index in cc and matrix associated with this node. - Should be either 1 or 2. - use_secondary_proj: If True, the secondary side projection operator is - used. Needed for periodic boundary conditions. - - """ - msg = """This function is deprecated and will be removed, most likely in the - second half of 2022. - - To assemble mixed-dimensional elliptic problems, the recommended solution is - either to use the models, or to use the automatic differentiation framework - directly. - """ - warn(msg, DeprecationWarning, stacklevel=2) - - if use_secondary_proj: - proj = intf.mortar_to_secondary_int() - else: - proj = intf.mortar_to_primary_int() - - hat_E_int = self._velocity_dof(sd, intf, proj) - - assert cc is not None - - cc[2, self_ind] -= hat_E_int.T * matrix[self_ind, self_ind] - cc[2, 2] -= hat_E_int.T * matrix[self_ind, self_ind] * hat_E_int - - def assemble_int_bound_pressure_trace_rhs( - self, sd, data, data_edge, cc, rhs, self_ind, use_secondary_proj=False - ): - """Assemble the rhs contribution from an internal - boundary, manifested as a condition on the boundary pressure. - - For details, see self.assemble_int_bound_pressure_trace() - - Parameters: - g: Grid which the condition should be imposed on. - data: Data dictionary for the node in the mixed-dimensional grid. - data_edge: Data dictionary for the edge in the mixed-dimensional grid. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - primary and secondary side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind: Index in cc and matrix associated with this node. - Should be either 1 or 2. - use_secondary_proj: If True, the secondary side projection operator is - used. Needed for periodic boundary conditions. - - """ - # Nothing to do here. - pass - - def assemble_int_bound_pressure_trace_between_interfaces( - self, - g: pp.Grid, - data_grid: dict, - data_primary_edge, - data_secondary_edge, - cc: np.ndarray, - matrix: np.ndarray, - rhs: np.ndarray, - ) -> None: - """Assemble the contribution from an internal - boundary, manifested as a condition on the boundary pressure. - - No contribution for this method. - - Parameters: - g: Grid which the condition should be imposed on. - data_grid: Data dictionary for the node in the mixed-dimensional grid. - data_primary_edge: Data dictionary for the primary edge in the - mixed-dimensional grid. - data_secondary_edge: Data dictionary for the secondary edge in the - mixed-dimensional grid. - cc (block matrix, 3x3): Block matrix of size 3 x 3, whwere each block - represents coupling between variables on this interface. Index 0, 1 and - 2 represent the primary grid, the primary and secondary interface, - respectively. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Block matrix of size 3 x 1, representing the right - hand side of this coupling. Index 0, 1 and 2 represent the primary grid, - the primary and secondary interface, respectively. - - """ - pass - - def assemble_int_bound_pressure_cell( - self, - sd: pp.Grid, - data: dict, - intf: pp.MortarGrid, - data_edge: dict, - cc: np.ndarray, - matrix: np.ndarray, - rhs: np.ndarray, - self_ind: int, - ) -> None: - """Abstract method. Assemble the contribution from an internal - boundary, manifested as a condition on the cell pressure. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose flux continuity on a lower-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - sd (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - secondary side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - primary and secondary side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - - """ - msg = """This function is deprecated and will be removed, most likely in the - second half of 2022. - - To assemble mixed-dimensional elliptic problems, the recommended solution is - either to use the models, or to use the automatic differentiation framework - directly. - """ - warn(msg, DeprecationWarning, stacklevel=2) - - proj = intf.secondary_to_mortar_avg() - - A = proj.T - shape = (sd.num_faces, A.shape[1]) - - cc[2, self_ind] -= sps.bmat([[sps.csr_matrix(shape)], [A]]).T - def enforce_neumann_int_bound( self, sd: pp.Grid, diff --git a/tests/grids/test_coarsening.py b/tests/grids/test_coarsening.py deleted file mode 100644 index 70783e3fd..000000000 --- a/tests/grids/test_coarsening.py +++ /dev/null @@ -1,507 +0,0 @@ -import inspect -import sys - -import numpy as np -import pytest - -import porepy as pp -from porepy.grids import coarsening as co -from porepy.numerics.linalg.matrix_operations import sparse_array_to_row_col_data - - -class TestPartitioning: - def reference(self): - caller_method = inspect.stack()[1][3] - return pp.test_utils.reference_dense_arrays.test_coarsening[caller_method] - - # Skipped: The test fails after update to scipy 1.13 due to changes in the sparse - # matrix format. Since the underlying code is marked for deprecation, the test will - # not be updated. - @pytest.mark.xfail - def test_coarse_grid_2d(self): - g = pp.CartGrid([3, 2]) - g.compute_geometry() - co.generate_coarse_grid(g, [5, 2, 2, 5, 2, 2]) - - assert g.num_cells == 2 - assert g.num_faces == 12 - assert g.num_nodes == 11 - - pt = np.tile(np.array([2, 1, 0]), (g.nodes.shape[1], 1)).T - find = np.isclose(pt, g.nodes).all(axis=0) - assert not find.any() - - faces_cell0, _, orient_cell0 = sparse_array_to_row_col_data(g.cell_faces[:, 0]) - assert np.array_equal(faces_cell0, [1, 2, 4, 5, 7, 8, 10, 11]) - assert np.array_equal(orient_cell0, [-1, 1, -1, 1, -1, -1, 1, 1]) - - faces_cell1, _, orient_cell1 = sparse_array_to_row_col_data(g.cell_faces[:, 1]) - assert np.array_equal(faces_cell1, [0, 1, 3, 4, 6, 9]) - assert np.array_equal(orient_cell1, [-1, 1, -1, 1, -1, 1]) - - known = np.array( - [ - [0, 4], - [1, 5], - [3, 6], - [4, 7], - [5, 8], - [6, 10], - [1, 0], - [2, 1], - [3, 2], - [8, 7], - [9, 8], - [10, 9], - ] - ) - - for f in np.arange(g.num_faces): - assert np.array_equal( - sparse_array_to_row_col_data(g.face_nodes[:, f])[0], known[f, :] - ) - - # Skipped: The test fails after update to scipy 1.13 due to changes in the sparse - # matrix format. Since the underlying code is marked for deprecation, the test will - # not be updated. - @pytest.mark.xfail - def test_coarse_grid_3d(self): - g = pp.CartGrid([2, 2, 2]) - g.compute_geometry() - co.generate_coarse_grid(g, [0, 0, 0, 0, 1, 1, 2, 2]) - - assert g.num_cells == 3 - assert g.num_faces == 30 - assert g.num_nodes == 27 - - faces_cell0, _, orient_cell0 = sparse_array_to_row_col_data(g.cell_faces[:, 0]) - known = [0, 1, 2, 3, 8, 9, 10, 11, 18, 19, 20, 21, 22, 23, 24, 25] - assert np.array_equal(faces_cell0, known) - known = [-1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1] - assert np.array_equal(orient_cell0, known) - - faces_cell1, _, orient_cell1 = sparse_array_to_row_col_data(g.cell_faces[:, 1]) - known = [4, 5, 12, 13, 14, 15, 22, 23, 26, 27] - assert np.array_equal(faces_cell1, known) - known = [-1, 1, -1, -1, 1, 1, -1, -1, 1, 1] - assert np.array_equal(orient_cell1, known) - - faces_cell2, _, orient_cell2 = sparse_array_to_row_col_data(g.cell_faces[:, 2]) - known = [6, 7, 14, 15, 16, 17, 24, 25, 28, 29] - assert np.array_equal(faces_cell2, known) - known = [-1, 1, -1, -1, 1, 1, -1, -1, 1, 1] - assert np.array_equal(orient_cell2, known) - - reference = self.reference()["face_nodes"] - for f in np.arange(g.num_faces): - assert np.array_equal( - sparse_array_to_row_col_data(g.face_nodes[:, f])[0], reference[f, :] - ) - - def test_coarse_grid_2d_1d(self): - part = np.array([0, 0, 1, 1, 2, 0, 3, 1]) - f = np.array([[2, 2], [0, 2]]) - - mdg = pp.meshing.cart_grid([f], [4, 2]) - mdg.compute_geometry() - co.generate_coarse_grid(mdg, (None, part)) - - # Test - known = np.array([1, 5, 18, 19]) - - for intf in mdg.interfaces(): - faces = sparse_array_to_row_col_data(intf.primary_to_mortar_int())[1] - assert np.array_equal(faces, known) - - def test_coarse_grid_2d_1d_cross(self): - # NOTE: Since for python 2.7 and 3.5 the meshes in gridbucket may have - # non-fixed order, we need to exclude this test. - if sys.version_info >= (3, 6): - part = np.zeros(36) - part[[0, 1, 2, 6, 7]] = 1 - part[[8, 14, 13]] = 2 - part[[12, 18, 19]] = 3 - part[[24, 30, 31, 32]] = 4 - part[[21, 22, 23, 27, 28, 29, 33, 34, 35]] = 5 - part[[9]] = 6 - part[[15, 16, 17]] = 7 - part[[9, 10]] = 8 - part[[20, 26, 25]] = 9 - part[[3, 4, 5, 11]] = 10 - f1 = np.array([[3.0, 3.0], [1.0, 5.0]]) - f2 = np.array([[1.0, 5.0], [3.0, 3.0]]) - - mdg = pp.meshing.cart_grid([f1, f2], [6, 6]) - mdg.compute_geometry() - - cell_centers_1 = np.array( - [ - [3.00000000e00, 3.00000000e00, 3.00000000e00, 3.00000000e00], - [4.50000000e00, 3.50000000e00, 2.50000000e00, 1.50000000e00], - [-1.66533454e-16, -5.55111512e-17, 5.55111512e-17, 1.66533454e-16], - ] - ) - cell_centers_2 = np.array( - [ - [4.50000000e00, 3.50000000e00, 2.50000000e00, 1.50000000e00], - [3.00000000e00, 3.00000000e00, 3.00000000e00, 3.00000000e00], - [-1.66533454e-16, -5.55111512e-17, 5.55111512e-17, 1.66533454e-16], - ] - ) - - co.generate_coarse_grid(mdg, (None, part)) - - # Test - for intf in mdg.interfaces(): - faces = sparse_array_to_row_col_data(intf.primary_to_mortar_int())[1] - - sd_primary, sd_secondary = mdg.interface_to_subdomain_pair(intf) - - if sd_secondary.dim == 0 and sd_primary.dim == 1: - known = [2, 5] - - elif sd_secondary.dim == 1 and sd_primary.dim == 2: - if np.allclose(sd_secondary.cell_centers, cell_centers_1): - known = [5, 10, 14, 18, 52, 53, 54, 55] - elif np.allclose(sd_secondary.cell_centers, cell_centers_2): - known = [37, 38, 39, 40, 56, 57, 58, 59] - else: - raise ValueError("Grid not found") - - assert np.array_equal(faces, known) - - def test_coarse_grid_3d_2d(self): - f = np.array([[2.0, 2.0, 2.0, 2.0], [0.0, 2.0, 2.0, 0.0], [0.0, 0.0, 2.0, 2.0]]) - mdg = pp.meshing.cart_grid([f], [4, 2, 2]) - mdg.compute_geometry() - - g = mdg.subdomains(dim=mdg.dim_max())[0] - part = np.zeros(g.num_cells) - part[g.cell_centers[0, :] < 2.0] = 1 - co.generate_coarse_grid(mdg, (None, part)) - # Test - # Be carefull! If the indexing of any grids (including mg) change the hard-coded - # indexes may be wrong - known_indices = np.array([1, 3, 0, 2, 5, 7, 4, 6]) - known = np.array([1, 4, 7, 10, 44, 45, 46, 47]) - - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - - assert np.array_equal(indices, known_indices) - assert np.array_equal(faces, known) - - def test_coarse_grid_3d_2d_cross(self): - # NOTE: Since for python 2.7 and 3.5 the meshes in gridbucket may have - # non-fixed order, we need to exclude this test. - if sys.version_info >= (3, 6): - f1 = np.array( - [[3.0, 3.0, 3.0, 3.0], [1.0, 5.0, 5.0, 1.0], [1.0, 1.0, 5.0, 5.0]] - ) - f2 = np.array( - [[1.0, 5.0, 5.0, 1.0], [1.0, 1.0, 5.0, 5.0], [3.0, 3.0, 3.0, 3.0]] - ) - mdg = pp.meshing.cart_grid([f1, f2], [6, 6, 6]) - mdg.compute_geometry() - - g = mdg.subdomains(dim=mdg.dim_max())[0] - part = np.zeros(g.num_cells) - p1, p2 = g.cell_centers[0, :] < 3.0, g.cell_centers[2, :] < 3.0 - part[np.logical_and(p1, p2)] = 1 - part[np.logical_and(p1, np.logical_not(p2))] = 2 - part[np.logical_and(np.logical_not(p1), p2)] = 3 - part[np.logical_and(np.logical_not(p1), np.logical_not(p2))] = 4 - - co.generate_coarse_grid(mdg, (None, part)) - - cell_centers_1 = self.reference()["cell_centers_1"] - cell_centers_2 = self.reference()["cell_centers_2"] - cell_centers_3 = self.reference()["cell_centers_3"] - - # Test - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - sd_primary, sd_secondary = mdg.interface_to_subdomain_pair(intf) - - if sd_primary.dim == 2 and sd_secondary.dim == 1: - reference_indices = [0, 1, 2, 3, 4, 5, 6, 7] - reference_faces = [17, 12, 7, 2, 43, 42, 41, 40] - - if sd_primary.dim == 3 and sd_secondary.dim == 2: - if np.allclose(sd_secondary.cell_centers, cell_centers_1): - reference_faces = self.reference()["faces_1"] - reference_indices = self.reference()["indices_1"] - elif np.allclose(sd_secondary.cell_centers, cell_centers_2): - reference_faces = self.reference()["faces_2"] - reference_indices = self.reference()["indices_2"] - elif np.allclose(sd_secondary.cell_centers, cell_centers_3): - reference_faces = self.reference()["faces_3"] - reference_indices = self.reference()["indices_3"] - else: - raise ValueError("Grid not found") - - assert np.array_equal(indices, np.array(reference_indices)) - assert np.array_equal(faces, np.array(reference_faces)) - - def test_create_partition_2d_cart(self): - g = pp.CartGrid([5, 5]) - g.compute_geometry() - part = co.create_partition(co._tpfa_matrix(g), g) - known = np.array( - [0, 0, 0, 1, 1, 0, 0, 2, 1, 1, 3, 2, 2, 2, 1, 3, 3, 2, 4, 4, 3, 3, 4, 4, 4] - ) - part = next(iter(part.values()))[1] - assert np.array_equal(part, known) - - def test_create_partition_2d_tri(self): - g = pp.StructuredTriangleGrid([3, 2]) - g.compute_geometry() - part = co.create_partition(co._tpfa_matrix(g), g) - known = np.array([1, 1, 1, 0, 0, 1, 0, 2, 2, 0, 2, 2]) - known_map = np.array([4, 3, 7, 5, 11, 8, 1, 2, 10, 6, 12, 9]) - 1 - part = next(iter(part.values()))[1] - assert np.array_equal(part, known[known_map]) - - def test_create_partition_2d_cart_cdepth4(self): - g = pp.CartGrid([10, 10]) - g.compute_geometry() - part = co.create_partition(co._tpfa_matrix(g), g, cdepth=4) - - part = next(iter(part.values()))[1] - assert np.array_equal(part, self.reference()["partition"]) - - def test_create_partition_3d_cart(self): - g = pp.CartGrid([4, 4, 4]) - g.compute_geometry() - part = co.create_partition(co._tpfa_matrix(g), g) - part = next(iter(part.values()))[1] - assert np.array_equal(part, self.reference()["partition"]) - - def test_create_partition_2d_1d_test0(self): - mdg, _ = pp.mdg_library.square_with_orthogonal_fractures( - "cartesian", - meshing_args={"cell_size": 1 / 2}, - fracture_indices=[1], - ) - - part = co.create_partition(co._tpfa_matrix(mdg), mdg) - co.generate_coarse_grid(mdg, part) - - # Test - known_indices = np.array([1, 0, 3, 2]) - known = np.array([6, 7, 10, 11]) - - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - assert np.array_equal(faces, known) - assert np.array_equal(indices, known_indices) - - def test_create_partition_2d_1d_test1(self): - mdg, _ = pp.mdg_library.square_with_orthogonal_fractures( - "cartesian", - meshing_args={"cell_size": 1 / 2}, - fracture_indices=[1], - fracture_endpoints=[np.array([0.5, 0.0])], - ) - part = co.create_partition(co._tpfa_matrix(mdg), mdg) - co.generate_coarse_grid(mdg, part) - - # Test - known_indices = np.array([0, 1]) - known = np.array([6, 9]) - - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - assert np.array_equal(faces, known) - assert np.array_equal(indices, known_indices) - - def test_create_partition_2d_1d_test2(self): - mdg, _ = pp.mdg_library.square_with_orthogonal_fractures( - "cartesian", - meshing_args={"cell_size": 1 / 2}, - fracture_indices=[1], - fracture_endpoints=[np.array([0.0, 0.5])], - ) - - seeds = co.generate_seeds(mdg) - known_seeds = np.array([0, 2]) - assert np.array_equal(seeds, known_seeds) - - part = co.create_partition(co._tpfa_matrix(mdg), mdg, seeds=seeds) - co.generate_coarse_grid(mdg, part) - - # Test - known_indices = np.array([0, 1]) - known = np.array([6, 10]) - - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - assert np.array_equal(faces, known) - assert np.array_equal(indices, known_indices) - - def test_create_partition_2d_1d_test3(self): - mdg, _ = pp.mdg_library.square_with_orthogonal_fractures( - "cartesian", - meshing_args={"cell_size": 1 / 2}, - fracture_indices=[1], - fracture_endpoints=[np.array([0.5, 1.0])], - ) - - part = co.create_partition(co._tpfa_matrix(mdg), mdg) - co.generate_coarse_grid(mdg, part) - - # Test - known_indices = np.array([0, 1]) - known = np.array([6, 9]) - - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - assert np.array_equal(faces, known) - assert np.array_equal(indices, known_indices) - - def test_create_partition_2d_1d_test4(self): - mdg, _ = pp.mdg_library.square_with_orthogonal_fractures( - "cartesian", - meshing_args={"cell_size": 1 / 2}, - fracture_indices=[1], - fracture_endpoints=[np.array([0.5, 1.0])], - ) - - seeds = co.generate_seeds(mdg) - known_seeds = np.array([1, 3]) - assert np.array_equal(seeds, known_seeds) - - part = co.create_partition(co._tpfa_matrix(mdg), mdg, seeds=seeds) - co.generate_coarse_grid(mdg, part) - - # Test - known_indices = np.array([0, 1]) - known = np.array([7, 10]) - - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - assert np.array_equal(faces, known) - assert np.array_equal(indices, known_indices) - - def test_create_partition_2d_1d_cross_test5(self): - # NOTE: Since for python 2.7 and 3.5 the meshes in gridbucket may have - # non-fixed order, we need to exclude this test. - if sys.version_info >= (3, 6): - f1 = np.array([[3.0, 3.0], [1.0, 5.0]]) - f2 = np.array([[1.0, 5.0], [3.0, 3.0]]) - mdg = pp.meshing.cart_grid([f1, f2], [6, 6]) - mdg.compute_geometry() - - part = co.create_partition(co._tpfa_matrix(mdg), mdg, cdepth=3) - co.generate_coarse_grid(mdg, part) - - cell_centers_1 = np.array( - [ - [3.00000000e00, 3.00000000e00, 3.00000000e00, 3.00000000e00], - [4.50000000e00, 3.50000000e00, 2.50000000e00, 1.50000000e00], - [-1.66533454e-16, -5.55111512e-17, 5.55111512e-17, 1.66533454e-16], - ] - ) - cell_centers_2 = np.array( - [ - [4.50000000e00, 3.50000000e00, 2.50000000e00, 1.50000000e00], - [3.00000000e00, 3.00000000e00, 3.00000000e00, 3.00000000e00], - [-1.66533454e-16, -5.55111512e-17, 5.55111512e-17, 1.66533454e-16], - ] - ) - - # Test - for intf in mdg.interfaces(): - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - sd_primary, sd_secondary = mdg.interface_to_subdomain_pair(intf) - - if sd_primary.dim == 1 and sd_secondary.dim == 0: - known = [2, 5] - known_indices = [0, 1] - - elif sd_primary.dim == 2 and sd_secondary.dim == 1: - g = sd_secondary - - if np.allclose(g.cell_centers, cell_centers_1): - known = [4, 9, 12, 16, 44, 45, 46, 47] - known_indices = [3, 2, 1, 0, 7, 6, 5, 4] - elif np.allclose(g.cell_centers, cell_centers_2): - known = [31, 32, 33, 34, 48, 49, 50, 51] - known_indices = [3, 2, 1, 0, 7, 6, 5, 4] - else: - raise ValueError("Grid not found") - - assert np.array_equal(faces, np.array(known)) - assert np.array_equal(indices, np.array(known_indices)) - - def test_create_partition_2d_1d_cross_test6(self): - # NOTE: Since for python 2.7 and 3.5 the meshes in gridbucket may have - # non-fixed order, we need to exclude this test. - if sys.version_info >= (3, 6): - f1 = np.array([[3.0, 3.0], [1.0, 5.0]]) - f2 = np.array([[1.0, 5.0], [3.0, 3.0]]) - mdg = pp.meshing.cart_grid([f1, f2], [6, 6]) - mdg.compute_geometry() - - seeds = co.generate_seeds(mdg) - known_seeds = np.array([8, 9, 26, 27, 13, 16, 19, 22]) - assert np.array_equal(np.sort(seeds), np.sort(known_seeds)) - - part = co.create_partition(co._tpfa_matrix(mdg), mdg, cdepth=3, seeds=seeds) - co.generate_coarse_grid(mdg, part) - - cell_centers_1 = np.array( - [ - [3.00000000e00, 3.00000000e00, 3.00000000e00, 3.00000000e00], - [4.50000000e00, 3.50000000e00, 2.50000000e00, 1.50000000e00], - [-1.66533454e-16, -5.55111512e-17, 5.55111512e-17, 1.66533454e-16], - ] - ) - cell_centers_2 = np.array( - [ - [4.50000000e00, 3.50000000e00, 2.50000000e00, 1.50000000e00], - [3.00000000e00, 3.00000000e00, 3.00000000e00, 3.00000000e00], - [-1.66533454e-16, -5.55111512e-17, 5.55111512e-17, 1.66533454e-16], - ] - ) - - # Test - for intf in mdg.interfaces(): - sd_primary, sd_secondary = mdg.interface_to_subdomain_pair(intf) - - indices, faces, _ = sparse_array_to_row_col_data( - intf.primary_to_mortar_int() - ) - - if sd_secondary.dim == 0 and sd_primary.dim == 1: - known = [2, 5] - known_indices = [0, 1] - - if sd_secondary.dim == 1 and sd_primary.dim == 2: - if np.allclose(sd_secondary.cell_centers, cell_centers_1): - known = [5, 10, 14, 18, 52, 53, 54, 55] - known_indices = [3, 2, 1, 0, 7, 6, 5, 4] - elif np.allclose(sd_secondary.cell_centers, cell_centers_2): - known = [37, 38, 39, 40, 56, 57, 58, 59] - known_indices = [3, 2, 1, 0, 7, 6, 5, 4] - else: - raise ValueError("Grid not found") - - assert np.array_equal(faces, np.array(known)) - assert np.array_equal(indices, np.array(known_indices))