Nonmatching grid tests

src/porepy/applications/md_grids/mdg_library.py

@@ -7,14 +7,15 @@
 from __future__ import annotations
 
 from pathlib import Path
-from typing import Literal, Optional, cast
+from typing import Literal, Optional, cast, Union
 
 import numpy as np
 
 import porepy as pp
 from porepy.fracs.fracture_network_2d import FractureNetwork2d
 from porepy.fracs.fracture_network_3d import FractureNetwork3d
 
+
 from . import domains, fracture_sets
 
 
@@ -24,6 +25,7 @@ def square_with_orthogonal_fractures(
     fracture_indices: list[int],
     fracture_endpoints: Optional[list[np.ndarray]] = None,
     size: pp.number = 1,
+    non_matching: bool = False,
     **meshing_kwargs,
 ) -> tuple[pp.MixedDimensionalGrid, FractureNetwork2d]:
     """Create a mixed-dimensional grid for a square domain with up to two orthogonal
@@ -39,9 +41,17 @@ def square_with_orthogonal_fractures(
             x coordinate of fracture two. Should have the same length as
             fracture_indices. If not provided, the endpoints will be set to [0, 1].
         size: Side length of square.
+        non_matching: If True, the fracture and 1d interface grids are refined to
+            generate a generally non-matching grid.
         **meshing_kwargs: Keyword arguments for meshing as used by
             :meth:`~porepy.grids.mdg_generation.create_mdg`.
 
+            If non_matching is True, the following additional arguments are also used:
+            - fracture_refinement_ratio: Refinement ratio for the fracture grids.
+                Defaults to 2.
+            - interface_refinement_ratio: Refinement ratio for the mortar grids.
+                Defaults to 2.
+
     Returns:
         Tuple containing a:
 
@@ -70,6 +80,89 @@ def square_with_orthogonal_fractures(
         FractureNetwork2d, pp.create_fracture_network(fractures, domain)
     )
     mdg = pp.create_mdg(grid_type, meshing_args, fracture_network, **meshing_kwargs)
+
+    # If a non-matching grid is requested, we refine the fracture and interface grids in
+    # the if-block below. However, if there are no fractures, we skip this step.
+    if non_matching and len(fracture_indices) > 0:
+        # Get the mesh refinement ratios for the fracture and interface grids..
+        fracture_refinement_ratio: pp.number = meshing_kwargs.get(
+            "fracture_refinement_ratio", 2
+        )
+        interface_refinement_ratio: pp.number = meshing_kwargs.get(
+            "interface_refinement_ratio", 2
+        )
+
+        # For the fracture grids, we can use the pp.refinement.refine_grid_1d function
+        # and impose the ratio directly. For the interface grids, the simplest way of
+        # refining is to create a new mixed-dimensional grid with the same geometry, but
+        # with refined rates, and then map the new grids to the old ones.
+
+        # Refine the mesh size arguments. Provided interface_refinement_ratio > 1, this
+        # will lead to a finer mesh than the one used to generate the old mdg.
+        if grid_type == "simplex":
+            # Loop over all the possible meshing arguments for the simplex grid, as
+            # specified by pp.create_mdg. Decrease the mesh size with the specified
+            # ratio.
+            for key in [
+                "cell_size",
+                "mesh_size_frac",
+                "mesh_size_min",
+                "mesh_size_bound",
+            ]:
+                if key in meshing_args:
+                    meshing_args[key] = meshing_args[key] / interface_refinement_ratio
+        elif grid_type == "cartesian":
+            # Loop over all the possible meshing arguments for the Cartesian grid.
+            # Refine.
+            for key in ["cell_size", "cell_size_x", "cell_size_y"]:
+                if key in meshing_args:
+                    meshing_args[key] = meshing_args[key] / interface_refinement_ratio
+        elif grid_type == "tensor_grid":
+            # This should also be possible to implement, but has not been prioritized.
+            raise NotImplementedError("Tensor grid not implemented yet.")
+
+        # Refine the fracture grids.
+        old_fracture_grids = mdg.subdomains(dim=1)
+        new_fracture_grids = [
+            pp.refinement.refine_grid_1d(g=old_grid, ratio=fracture_refinement_ratio)
+            for old_grid in old_fracture_grids
+        ]
+        # Create a mapping between the old and new fracture grids.
+        grid_map = dict(zip(old_fracture_grids, new_fracture_grids))
+
+        # Create a new mdg for the same domain, but with a modified mesh size.
+        mdg_new, _ = square_with_orthogonal_fractures(
+            grid_type,
+            meshing_args,
+            fracture_indices,
+            fracture_endpoints,
+            size,
+            non_matching=False,
+        )
+
+        intf_map = {}
+        # Loop over all the interfaces in the new grid. Find its corresponding interface
+        # in the old grid, and update the interface map. There is no point in updating
+        # the 0d mortar, since there is no room for refinement.
+        for intf in mdg_new.interfaces(dim=1):
+            # Then, for each interface, we fetch the secondary grid which belongs to it.
+            _, g_sec_new = mdg_new.interface_to_subdomain_pair(intf)
+
+            # We then loop through all the interfaces in the original grid (recipient).
+            for intf_old in mdg.interfaces(dim=1):
+                # Fetch the secondary grid of the interface in the old grid.
+                _, g_sec_old = mdg.interface_to_subdomain_pair(intf_old)
+
+                # Checking the fracture number of the secondary grid in the old mdg. If
+                # they are the same, i.e., if the fractures are the same ones, we update
+                # the interface map.
+                if g_sec_old.frac_num == g_sec_new.frac_num:
+                    intf_map.update({intf_old: intf})
+
+        # Finally replace the subdomains and interfaces in the original
+        # mixed-dimensional grid.
+        mdg.replace_subdomains_and_interfaces(sd_map=grid_map, intf_map=intf_map)
+
     return mdg, fracture_network
 
 

src/porepy/applications/md_grids/model_geometries.py

@@ -1,8 +1,12 @@
 from __future__ import annotations
 
+from typing import Union
+
 import numpy as np
 
 import porepy as pp
+from porepy.grids import mortar_grid
+from porepy.grids.mdg_generation import _preprocess_cartesian_args
 
 from . import domains, fracture_sets
 
@@ -155,3 +159,97 @@ def meshing_arguments(self) -> dict:
             "cell_size_min": 0.2 * ls,
         }
         return mesh_sizes
+
+
+class NonMatchingSquareDomainOrthogonalFractures(SquareDomainOrthogonalFractures):
+    """Create a non-matching mixed-dimensional grid of a square domain with up to two
+    orthogonal fractures.
+
+    The setup is similar to :class:`SquareDomainOrthogonalFractures`, but the
+    geometry allows for non-matching grids and different resolution for each grid.
+    """
+
+    def set_geometry(self) -> None:
+        """Define geometry and create a non-matching mixed-dimensional grid.
+
+        We here make a non-matching mixed-dimensional grid in the sense that neither of
+        the fracture grids, interface grids or the rock grid are matching. This is done
+        by refining the fracture and interface grids.
+
+        """
+
+        # The fracture grids are replaced by a refined version of the fractures which
+        # are already present in the grid. The interface grids are replaced by first
+        # creating a new mixed-dimensional grid with a higher resolution. The new
+        # mixed-dimensional grid will "donate" its interface grids to the original
+        # mixed-dimensional grid.
+
+        # First set the geometry and create a matching mixed-dimensional grid.
+        assert isinstance(self, pp.ModelGeometry)
+        super().set_geometry()  # type:ignore[safe-super]
+        self.mdg = pp.mdg_library.square_with_orthogonal_fractures(????, non_matching=True)
+
+        # Refine and replace fracture grids. First we fetch the old fracture grids:
+        old_fracture_grids = self.mdg.subdomains(dim=1)
+
+        # Ratios which we want to refine the fracture grids with. Default, if no ratios
+        # are provided, is to refine by a factor of 2.
+        ratios = self.params.get(
+            "fracture_refinement_ratios", np.ones(len(old_fracture_grids) * 2)
+        )
+
+        # The actual refinement of the fracture grids.
+        new_fracture_grids = [
+            pp.refinement.refine_grid_1d(g=old_grid, ratio=ratio)
+            for old_grid, ratio in zip(old_fracture_grids, ratios)
+        ]
+
+        # Create a mapping between the old and new fracture grids.
+        grid_map = dict(zip(old_fracture_grids, new_fracture_grids))
+
+        # Refine and replace interface grids.
+        # Directly refining the already existing interface grids is not straightforward.
+        # Instead, we create a new mixed-dimensional grid and donate the interface grids
+        # to the original mixed-dimensional grid.
+
+        # First we gather the number of cells in x- and y- direction, as well as the
+        # physical dimensions of the original mixed-dimensional grid:
+        domain = pp.Domain(self.domain.bounding_box)
+        meshing_args = self.meshing_arguments()
+        nx_cells, phys_dims, _ = _preprocess_cartesian_args(domain, meshing_args, {})
+
+        # Then we create a new mdg with the same fractures and physical dimensions as
+        # the old one, but with a higher resolution:
+        nx_cells_new = np.array(nx_cells) * self.params.get(
+            "interface_refinement_ratio", 2
+        )
+        fracs = [self.fractures[i].pts for i in range(len(self.fractures))]
+        mdg_new = pp.meshing.cart_grid(fracs, nx_cells_new, physdims=phys_dims)
+
+        # Loop through all the interfaces in the new mixed-dimensional grid (donor) such
+        # that we can fill intf_map with a map between old and new interface grids.
+        intf_map: dict[
+            pp.MortarGrid,
+            Union[pp.MortarGrid, dict[mortar_grid.MortarSides, pp.Grid]],
+        ] = {}
+        for intf in mdg_new.interfaces(dim=1):
+            # Then, for each interface, we fetch the secondary grid which belongs to it.
+            _, g_sec = mdg_new.interface_to_subdomain_pair(intf)
+
+            # We then loop through all the interfaces in the original grid (recipient).
+            for intf_coarse in self.mdg.interfaces(dim=1):
+                # Fetch the secondary grid of the interface in the original grid.
+                _, g_sec_coarse = self.mdg.interface_to_subdomain_pair(intf_coarse)
+
+                # Checking the fracture number of the secondary grid in the recipient
+                # mdg. If they are the same, i.e., if the fractures are the same ones,
+                # we update the interface map.
+                if g_sec_coarse.frac_num == g_sec.frac_num:
+                    intf_map.update({intf_coarse: intf})
+
+        # Finally replace the subdomains and interfaces in the original
+        # mixed-dimensional grid.
+        self.mdg.replace_subdomains_and_interfaces(sd_map=grid_map, intf_map=intf_map)
+
+        # Create projections between local and global coordinates for fracture grids.
+        pp.set_local_coordinate_projections(self.mdg)

src/porepy/grids/refinement.py

@@ -90,6 +90,9 @@ def refine_grid_1d(g: pp.Grid, ratio: int = 2) -> pp.Grid:
     cell_nodes = g.cell_nodes().tocsc()
     nodes, cells, _ = sparse_array_to_row_col_data(cell_nodes)
 
+    # Extract the fracture number.
+    frac_num = g.frac_num
+
     # Every cell will contribute (ratio - 1) new nodes.
     num_new_nodes = (ratio - 1) * g.num_cells + g.num_nodes
 
@@ -217,6 +220,9 @@ def refine_grid_1d(g: pp.Grid, ratio: int = 2) -> pp.Grid:
     g = Grid(1, x, face_nodes, cell_faces, "Refined 1d grid")
     g.compute_geometry()
 
+    # Keep the original fracture number.
+    g.frac_num = frac_num
+
     return g
 
 

tests/applications/md_grids/test_mdg_library.py

@@ -2,6 +2,7 @@
 Some of these tests are sensitive to meshing or node ordering. If this turns out to
 cause problems, we deactivate the corresponding asserts.
 """
+
 import pytest
 import numpy as np
 
@@ -45,6 +46,7 @@ def check_fracture_coordinates(self, cell_centers, face_centers):
             assert np.all(np.isclose(g.face_centers, face_centers[i]))
 
     def check_intersections(self, n_intersections):
+        """Check that the number of intersections is as expected."""
         assert len(self.mdg.subdomains(dim=0)) == n_intersections
 
     def check_domain(self, x_length, y_length, z_length=None):
@@ -113,6 +115,67 @@ def test_two_intersecting_custom_values(self):
 
         self.check_fracture_coordinates([cc0, cc1], [fc0, fc1])
 
+    @pytest.mark.parametrize("mesh_type", ["simplex", "cartesian"])
+    @pytest.mark.parametrize("num_fractures", [1, 2])
+    def test_two_intersecting_nonmatching(self, mesh_type, num_fractures):
+        """Test meshing of a 2d square domain to generate a non-matching grid.
+
+        This is not a very powerful test, but it does verify that the generated grids
+        have  different numbers of cells in the fractures and along the mortars, and
+        that the mortar grids have a different number of grids than the matrix faces
+        tagged as on a fracture.
+
+        Parameters:
+            mesh_type: The type of mesh to generate, either 'simplex' or 'cartesian'.
+            num_fractures: The number of fractures to include in the domain.
+
+        """
+        meshing_arguments = {"cell_size": 1 / 4}
+
+        # Define the fracture endpoints and indices according to the input.
+        fracture_endpoints = []
+        fracture_indices = []
+        if num_fractures > 0:
+            fracture_endpoints.append(np.array([1 / 4, 3 / 4]))
+            fracture_indices.append(0)
+        if num_fractures > 1:
+            fracture_endpoints.append(np.array([0, 1]))
+            fracture_indices.append(1)
+
+        # Generate
+        self.mdg, _ = pp.mdg_library.square_with_orthogonal_fractures(
+            mesh_type,
+            meshing_arguments,
+            fracture_indices=fracture_indices,
+            fracture_endpoints=fracture_endpoints,
+            non_matching=True,
+            **{"interface_refinement_ratio": 2, "fracture_refinement_ratio": 4},
+        )
+
+        # Number of faces in the 2d grid that are tagged as fracture faces.
+        num_fracture_faces_from_matrix = (
+            self.mdg.subdomains(dim=2)[0].tags["fracture_faces"].sum()
+        )
+
+        # Loop over the fractures, fetch the projections to the 2d grid, and count the
+        # number of non-zero entries in the projection matrix.
+        non_zero_projection_primary = 0
+        for mg in self.mdg.interfaces(dim=1):
+            proj_primary = mg.mortar_to_primary_avg()
+            non_zero_projection_primary += proj_primary.nnz
+        # A matching grid would have equality here. We have refined the mortar grid,
+        # hence there should be more items in the projection matrix than there are
+        # fracture faces in the matrix grid.
+        assert non_zero_projection_primary > num_fracture_faces_from_matrix
+
+        # For the mortar grid vs fracture grids, we can simply do a cell count. We do
+        # not know precisely how many cells there will be in each, but with the given
+        # refinement settings (above), they should at least not be equal.
+        for mg in self.mdg.interfaces(dim=1):
+            _, sd_secondary = self.mdg.interface_to_subdomain_pair(mg)
+            # Multiply by two since there are two sides of the mortar.
+            assert 2 * sd_secondary.num_cells != mg.num_cells
+
     def test_benchmark_regular(self):
         """Test the mdg generator for the regular case of the benchmark study.