Nonmatching grid tests

src/porepy/applications/md_grids/model_geometries.py

@@ -155,3 +155,96 @@ def meshing_arguments(self) -> dict:
             "cell_size_min": 0.2 * ls,
         }
         return mesh_sizes
+
+
+class NonMatchingSquareDomainOrthogonalFractures(SquareDomainOrthogonalFractures):
+    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 grids and interfaces.
+
+        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.
+
+        There are a few things to be aware of when creating a non-matching
+        mixed-dimensional grid:
+            * If _both_ the interface and fracture grids are to be refined/coarsened,
+              you must be aware of:
+                * If you first refine (and replace) the fracture grids, you must also
+                  update the fracture numbering such that the new fracture grids have
+                  the same fracture number as the old ones. This is because the fracture
+                  numbering is used to identify the correct interface grids later.
+                  TODO: Make sure that the refine-grid-functions also copy frac_num
+                  attribute such that the abovementioned is no longer a problem.
+            * Ensure that the "donor" mixed-dimensional grid is physically the same as
+              the "recipient" mixed-dimensional grid. Fractures must be located at the
+              same place and the physical dimension of the grids must be the same.
+
+        """
+        super().set_geometry()
+
+        # Refine and replace fracture grids:
+        old_fracture_grids = self.mdg.subdomains(dim=1)
+
+        # Ratios which we want to refine the fracture grids with.
+        ratios = []
+        for i in range(len(old_fracture_grids)):
+            ratios.append(i + 2)
+
+        new_fracture_grids = [
+            pp.refinement.refine_grid_1d(g=old_grid, ratio=ratio)
+            for old_grid, ratio in zip(old_fracture_grids, ratios)
+        ]
+
+        grid_map = dict(zip(old_fracture_grids, new_fracture_grids))
+
+        # Refine and replace interface grids:
+        # We first create a new and more refined mixed-dimensional grid.
+        def mdg_func(self, nx=2, ny=2) -> pp.MixedDimensionalGrid:
+            """Generate a refined version of an already existing mixed-dimensional grid.
+
+            Parameters:
+                nx: Number of cells in x-direction.
+                ny: Number of cells in y-direction.
+
+            """
+            fracs = [self.fractures[i].pts for i in range(len(self.fractures))]
+            domain = self.domain.bounding_box
+            md_grid = pp.meshing.cart_grid(
+                fracs, np.array([nx, ny]), physdims=[domain["xmax"], domain["ymax"]]
+            )
+            return md_grid
+
+        mdg_new = mdg_func(self, nx=6, ny=6)
+        g_new = mdg_new.subdomains(dim=2)[0]
+        g_new.compute_geometry()
+
+        intf_map = {}
+
+        # First we loop through all the interfaces in the new mixed-dimensional grid
+        # (donor).
+        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)
+
+                # Checkinc the fracture number of the secondary grid in the recipient
+                # mdg. If they are the same, i.e., 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. Both can be done at the same time:
+        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,10 @@ 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 +221,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/models/test_fluid_mass_balance.py

@@ -27,6 +27,7 @@
 from porepy.applications.md_grids.model_geometries import (
     CubeDomainOrthogonalFractures,
     SquareDomainOrthogonalFractures,
+    NonMatchingSquareDomainOrthogonalFractures,
 )
 from porepy.applications.test_utils import models, well_models
 from porepy.models.fluid_mass_balance import SinglePhaseFlow
@@ -37,6 +38,7 @@
     extended_water_values_for_testing as water_values,
 )
 
+
 @pytest.fixture(scope="function")
 def model_setup():
     """Minimal compressible single-phase flow setup with two intersecting fractures.
@@ -462,16 +464,24 @@ def test_ad_operator_methods_single_phase_flow(
         {"m": 2, "kg": 3, "s": 1, "K": 1},
     ],
 )
-def test_unit_conversion(units):
+@pytest.mark.parametrize(
+    "grid",
+    [
+        SquareDomainOrthogonalFractures,
+        NonMatchingSquareDomainOrthogonalFractures,
+    ],
+)
+def test_unit_conversion(units, grid):
     """Test that solution is independent of units.
 
     Parameters:
         units (dict): Dictionary with keys as those in
             :class:`~pp.compositional.materials.Constants`.
+        grid: Mixed dimensional grid class which the test is performed on.
 
     """
 
-    class Model(SquareDomainOrthogonalFractures, SinglePhaseFlow):
+    class Model(grid, SinglePhaseFlow):
         """Single phase flow model in a domain with two intersecting fractures."""
 
         def bc_values_pressure(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray: