Merge pull request #1346 from keileg/matrix_slicer

Introduce matrix slicer objects

Author: keileg

src/porepy/applications/test_utils/arrays.py

@@ -3,6 +3,8 @@
 import numpy as np
 import scipy.sparse as sps
 
+import porepy as pp
+
 
 def compare_arrays(
     a: np.ndarray, b: np.ndarray, tol: float = 1e-4, sort: bool = True
@@ -68,3 +70,37 @@ def compare_matrices(m1: sps.spmatrix, m2: sps.spmatrix, tol: float = 1e-10) ->
         if np.max(np.abs(d.data)) > tol:
             return False
     return True
+
+
+def projection_matrix_from_array_slicers(
+    slicers: pp.matrix_operations.ArraySlicer | list[pp.matrix_operations.ArraySlicer],
+    dim: int,
+) -> sps.coo_matrix:
+    """Recover a projection matrix from a one or multiple matrix slicers.
+
+    For a set of slicers {P_1, P_2, ..., P_n}, obtain the matrix P that corresponds to
+
+        P = P_1 + P_2 + ... + P_n.
+
+    Parameters:
+        slicers: One or more ArraySlicers.
+        dim: Dimension of the domain space of the slicer.
+
+    Returns:
+        Matrix representation of the basis vector.
+
+    """
+    # Always deal with a list of slicers.
+    if isinstance(slicers, pp.matrix_operations.ArraySlicer):
+        slicers = [slicers]
+
+    # Initialize result.
+    result = None
+
+    for slicer in slicers:
+        if result is None:
+            result = slicer @ np.eye(dim)
+        else:
+            result += slicer @ np.eye(dim)
+
+    return sps.coo_matrix(result)

src/porepy/fracs/split_grid.py

@@ -445,7 +445,7 @@ def _update_face_cells(
         if j == i:
             # We hit the target neighbor.
             # Pick out the part of f_c to be used with this neighbor.
-            f_c_sliced = pp.matrix_operations.slice_mat(f_c, face_id)
+            f_c_sliced = pp.matrix_operations.slice_sparse_matrix(f_c, face_id)
             # The new face-cell relations are added to the end of the matrix (since
             # the faces were added to the end of the face arrays in the
             # higher-dimensional grid). Columns (face-indices in the higher
@@ -1075,14 +1075,10 @@ def _find_cell_color(sd: pp.Grid, cells: np.ndarray) -> np.ndarray:
     c = np.sort(cells)
     # Local cell-face and face-node maps.
     assert sd.cell_faces.getformat() == "csc"
-    cell_faces = pp.matrix_operations.slice_mat(sd.cell_faces, c)
+    cell_faces = pp.matrix_operations.slice_sparse_matrix(sd.cell_faces, c)
     child_cell_ind = -np.ones(sd.num_cells, dtype=int)
     child_cell_ind[c] = np.arange(cell_faces.shape[1])
 
-    # Create a copy of the cell-face relation, so that we can modify it at will.
-    # RB: I don't think this is necessary as slice_mat creates a copy cell_faces =
-    # cf_sub.copy()
-
     # Direction of normal vector does not matter here, only 0s and 1s
     cell_faces.data = np.abs(cell_faces.data)
 

src/porepy/grids/partition.py

@@ -525,7 +525,7 @@ def _extract_submatrix(
     """
     if mat.format != "csc":
         raise ValueError("To extract columns from a matrix, it must be csc")
-    sub_mat = pp.matrix_operations.slice_mat(mat, ind)
+    sub_mat = pp.matrix_operations.slice_sparse_matrix(mat, ind)
     cols = sub_mat.indptr
     data = sub_mat.data
     unique_rows, rows_sub = np.unique(sub_mat.indices, return_inverse=True)

src/porepy/models/constitutive_laws.py

@@ -96,7 +96,7 @@ def elastic_displacement_jump(self, subdomains: list[pp.Grid]) -> pp.ad.Operator
         # respectively, to nd dimensions.
         basis = self.basis(subdomains, dim=self.nd)
         local_basis = self.basis(subdomains, dim=self.nd - 1)
-        tangential_to_nd = pp.ad.sum_operator_list(
+        tangential_to_nd = pp.ad.sum_projection_list(
             [e_nd @ e_f.T for e_nd, e_f in zip(basis[:-1], local_basis)]
         )
         normal_to_nd = basis[-1]
@@ -1117,7 +1117,7 @@ def interface_vector_source_darcy_flux(
         # composed by summation).
         normals_times_source = normals * vector_source
         # Then sum over the nd dimensions. The result will in effect be a matrix.
-        nd_to_scalar_sum = pp.ad.sum_operator_list(
+        nd_to_scalar_sum = pp.ad.sum_projection_list(
             [e.T for e in self.basis(interfaces, dim=self.nd)]
         )
         # Finally, the dot product between normal vectors and the vector source. This
@@ -3104,9 +3104,7 @@ def fracture_pressure_stress(
 
         # Expands from cell-wise scalar to vector. Equivalent to the :math:`\mathbf{I}p`
         # operation.
-        scalar_to_nd = pp.ad.sum_operator_list(
-            [e_i for e_i in self.basis(interfaces, dim=self.nd)]
-        )
+        scalar_to_nd = pp.ad.sum_projection_list(self.basis(interfaces, dim=self.nd))
         # Spelled out, from the right: Project the pressure from the fracture to the
         # mortar, expand to an nd-vector, and multiply with the outwards normal vector.
         stress = outwards_normal * (

src/porepy/models/contact_mechanics.py

@@ -181,14 +181,12 @@ def tangential_fracture_deformation_equation(
         tangential_basis = self.basis(subdomains, dim=self.nd - 1)
 
         # To map a scalar to the tangential plane, we need to sum the basis vectors. The
-        # individual basis functions have shape (Nc * (self.nd - 1), Nc), where Nc is
-        # the total number of cells in the subdomain. The sum will have the same shape,
-        # but the row corresponding to each cell will be non-zero in all rows
-        # corresponding to the tangential basis vectors of this cell. EK: mypy insists
-        # that the argument to sum should be a list of booleans. Ignore this error.
-        scalar_to_tangential = pp.ad.sum_operator_list(
-            [e_i for e_i in tangential_basis]
-        )
+        # individual basis vectors can be represented as projection matrices of shape
+        # (Nc * (self.nd - 1), Nc), where Nc is the total number of cells in the
+        # subdomain. The matrix representation of the sum has the same shape, but the
+        # row corresponding to each cell will be non-zero in all rows corresponding to
+        # the tangential basis vectors of this cell.
+        scalar_to_tangential = pp.ad.sum_projection_list(tangential_basis)
 
         # Variables: The tangential component of the contact traction and the plastic
         # displacement jump.
@@ -207,39 +205,28 @@ def tangential_fracture_deformation_equation(
         f_norm = pp.ad.Function(partial(pp.ad.l2_norm, self.nd - 1), "norm_function")
 
         # The numerical constant is used to loosen the sensitivity in the transition
-        # between sticking and sliding.
-        # Expanding using only left multiplication to with scalar_to_tangential does not
-        # work for an array, unlike the operators below. Arrays need right
-        # multiplication as well.
+        # between sticking and sliding. Expanding using only left multiplication to with
+        # scalar_to_tangential does not work for an array, unlike the operators below.
+        # Arrays need right multiplication as well.
         c_num_as_scalar = self.contact_mechanics_numerical_constant(subdomains)
 
-        # The numerical parameter is a cell-wise scalar which must be extended to a
-        # vector quantity to be used in the equation (multiplied from the right).
-        # Spelled out, from the right: Restrict the vector quantity to one dimension in
-        # the tangential plane (e_i.T), multiply with the numerical parameter, prolong
-        # to the full vector quantity (e_i), and sum over all all directions in the
-        # tangential plane.
-        c_num = pp.ad.sum_operator_list(
-            [e_i * c_num_as_scalar * e_i.T for e_i in tangential_basis]
-        )
-
-        # Combine the above into expressions that enter the equation. c_num will
-        # effectively be a sum of SparseArrays, thus we use a matrix-vector product @
-        tangential_sum = t_t + c_num @ u_t_increment
+        # The numerical parameter is a cell-wise scalar, or a single scalar common for
+        # all cells. In both cases, it must be extended to a vector quantity to be used
+        # in the equation (multiplied from the right). Do this by multiplying with the
+        # sum of the tangential basis vectors. Then take a Hadamard product with the
+        # tangential displacement jump and add to the tangential component of the
+        # contact traction to arrive at the expression that enters the equation.
+        basis_sum = pp.ad.sum_projection_list(tangential_basis)
+        tangential_sum = t_t + (basis_sum @ c_num_as_scalar) * u_t_increment
 
         norm_tangential_sum = f_norm(tangential_sum)
         norm_tangential_sum.set_name("norm_tangential")
 
         b_p = f_max(self.friction_bound(subdomains), zeros_frac)
         b_p.set_name("bp")
 
-        # Remove parentheses to make the equation more readable if possible. The product
-        # between (the SparseArray) scalar_to_tangential and b_p is of matrix-vector
-        # type (thus @), and the result is then multiplied elementwise with
-        # tangential_sum.
         bp_tang = (scalar_to_tangential @ b_p) * tangential_sum
 
-        # For the use of @, see previous comment.
         maxbp_abs = scalar_to_tangential @ f_max(b_p, norm_tangential_sum)
 
         # The characteristic function below reads "1 if (abs(b_p) < tol) else 0".
@@ -538,13 +525,12 @@ def contact_mechanics_open_state_characteristic(
         tangential_basis = self.basis(subdomains, dim=self.nd - 1)
 
         # To map a scalar to the tangential plane, we need to sum the basis vectors. The
-        # individual basis functions have shape (Nc * (self.nd - 1), Nc), where Nc is
-        # the total number of cells in the subdomain. The sum will have the same shape,
-        # but the row corresponding to each cell will be non-zero in all rows
-        # corresponding to the tangential basis vectors of this cell.
-        scalar_to_tangential = pp.ad.sum_operator_list(
-            [e_i for e_i in tangential_basis]
-        )
+        # individual basis functions can be represented as projection matrices of size
+        # (Nc * (self.nd - 1), Nc), where Nc is # the total number of cells in the
+        # subdomain. The sum of basis vectors can likewise be represented as a matrix of
+        # the same shape, but the row corresponding to each cell will be non-zero in all
+        # rows corresponding to the tangential basis vectors of this cell.
+        scalar_to_tangential = pp.ad.sum_projection_list(tangential_basis)
 
         # With the active set method, the performance of the Newton solver is sensitive
         # to changes in state between sticking and sliding. To reduce the sensitivity to

src/porepy/models/geometry.py

@@ -278,12 +278,13 @@ def wrap_grid_attribute(
         array.set_name(f"Array wrapping attribute {attr} on {len(grids)} grids.")
         return array
 
-    def basis(self, grids: Sequence[pp.GridLike], dim: int) -> list[pp.ad.SparseArray]:
+    def basis(self, grids: Sequence[pp.GridLike], dim: int) -> list[pp.ad.Projection]:
         """Return a cell-wise basis for all subdomains.
 
-        The basis is represented as a list of matrices, each of which represents a
-        basis function. The individual matrices have shape ``Nc * dim, Nc`` where
-        ``Nc`` is the total number of cells in the subdomains.
+        The basis is represented as a list of projections, each of which represents a
+        basis function. The individiual basis functions can be represented as a
+        projection matrix, of shape ``Nc * dim, Nc`` where ``Nc`` is the total number of
+        cells in the subdomains.
 
         Examples:
             To extend a cell-wise scalar to a vector field, use
@@ -307,32 +308,33 @@ def basis(self, grids: Sequence[pp.GridLike], dim: int) -> list[pp.ad.SparseArra
             function.
 
         """
-        # Collect the basis functions for each dimension
-        basis: list[pp.ad.SparseArray] = []
+        # Collect the basis functions for each dimension.
+        basis: list[pp.ad.Projection] = []
         for i in range(dim):
             basis.append(self.e_i(grids, i=i, dim=dim))
-        # Stack the basis functions horizontally
+        # Stack the basis functions horizontally.
         return basis
 
     def e_i(
         self, grids: Sequence[pp.GridLike], *, i: int, dim: int
-    ) -> pp.ad.SparseArray:
+    ) -> pp.ad.Projection:
         """Return a cell-wise basis function in a specified dimension.
 
         It is assumed that the grids are embedded in a space of dimension dim and
         aligned with the coordinate axes, that is, the reference space of the grid.
         Moreover, the grid is assumed to be planar.
 
         Example:
-            For a grid with two cells, and with `i=1` and `dim=3`, the returned
-            basis will be (after conversion to a numpy array)
+            For a grid with two cells, and with `i=1` and `dim=3`, the returned basis
+            will be a Projection that is equivalent to applying the following projection
+            matrix:
             .. code-block:: python
                 array([[0., 0.],
-                    [1., 0.],
-                    [0., 0.],
-                    [0., 0.],
-                    [0., 1.],
-                    [0., 0.]])
+                       [1., 0.],
+                       [0., 0.],
+                       [0., 0.],
+                       [0., 1.],
+                       [0., 0.]])
 
         See also:
             :meth:`basis` for the construction of a full basis.
@@ -343,16 +345,12 @@ def e_i(
             dim: Dimension of the functions.
 
         Returns:
-            pp.ad.SparseArray: Ad representation of a matrix with the basis
-            functions as columns.
+            Ad projection that represents a basis function.
 
         Raises:
             ValueError: If i is larger than dim - 1.
 
         """
-        # TODO: Should we expand this to grids not aligned with the coordinate axes, and
-        # possibly unify with ``porepy.utils.projections.TangentialNormalProjection``?
-        # This is not a priority for the moment, though.
 
         if dim is None:
             dim = self.nd
@@ -361,15 +359,18 @@ def e_i(
         if i >= dim:
             raise ValueError("Basis function index out of range")
 
-        # Construct a single vector, and later stack it to a matrix
-        # Collect the basis functions for each dimension
-        e_i = np.zeros((dim, 1))
-        e_i[i] = 1
         # Expand to cell-wise column vectors.
         num_cells = sum([g.num_cells for g in grids])
-        # Expand to a matrix.
-        mat = sps.kron(sps.eye(num_cells), e_i)
-        return pp.ad.SparseArray(mat)
+        range_ind = np.arange(i, dim * num_cells, dim)
+
+        slicer = pp.ad.Projection(
+            domain_indices=np.arange(num_cells),
+            range_indices=range_ind,
+            range_size=num_cells * dim,
+            domain_size=num_cells,
+        )
+
+        return slicer
 
     # Local basis related methods
     def tangential_component(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
@@ -394,7 +395,7 @@ def tangential_component(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         # it in the tangential basis. The two operations are combined in a single
         # operator composed right to left: v will be hit by first e_i.T (row vector) and
         # secondly t_i (column vector).
-        op: pp.ad.Operator = pp.ad.sum_operator_list(
+        op: pp.ad.Operator = pp.ad.sum_projection_list(
             [
                 self.e_i(subdomains, i=i, dim=self.nd - 1)
                 @ self.e_i(subdomains, i=i, dim=self.nd).T
@@ -404,7 +405,7 @@ def tangential_component(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         op.set_name("tangential_component")
         return op
 
-    def normal_component(self, subdomains: list[pp.Grid]) -> pp.ad.SparseArray:
+    def normal_component(self, subdomains: list[pp.Grid]) -> pp.ad.Projection:
         """Compute the normal component of a vector field.
 
         The normal space is defined according to the local coordinates of the
@@ -421,8 +422,8 @@ def normal_component(self, subdomains: list[pp.Grid]) -> pp.ad.SparseArray:
             subdomains: List of grids on which the vector field is defined.
 
         Returns:
-            Matrix extracting normal component of the vector field and expressing it
-            in normal basis. The size of the matrix is `(Nc, Nc * self.nd)`, where
+            Projection extracting normal component of the vector field and expressing it
+            in normal basis. The size of the projection is `(Nc, Nc * self.nd)`, where
             `Nc` is the total number of cells in the subdomains.
 
         """

src/porepy/models/protocol.py

@@ -198,7 +198,7 @@ def wrap_grid_attribute(
 
         def basis(
             self, grids: Sequence[pp.GridLike], dim: int
-        ) -> list[pp.ad.SparseArray]:
+        ) -> list[pp.ad.Projection]:
             """Return a cell-wise basis for all subdomains.
 
             The basis is represented as a list of matrices, each of which represents a
@@ -230,7 +230,7 @@ def basis(
 
         def e_i(
             self, grids: Sequence[pp.GridLike], *, i: int, dim: int
-        ) -> pp.ad.SparseArray:
+        ) -> pp.ad.Projection:
             """Return a cell-wise basis function in a specified dimension.
 
             It is assumed that the grids are embedded in a space of dimension dim and
@@ -257,8 +257,7 @@ def e_i(
                 dim: Dimension of the functions.
 
             Returns:
-                pp.ad.SparseArray: Ad representation of a matrix with the basis
-                functions as columns.
+                Ad projection that represents a basis function.
 
             Raises:
                 ValueError: If i is larger than dim - 1.
@@ -283,7 +282,7 @@ def tangential_component(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
 
             """
 
-        def normal_component(self, subdomains: list[pp.Grid]) -> pp.ad.SparseArray:
+        def normal_component(self, subdomains: list[pp.Grid]) -> pp.ad.Projection:
             """Compute the normal component of a vector field.
 
             The normal space is defined according to the local coordinates of the
@@ -300,9 +299,9 @@ def normal_component(self, subdomains: list[pp.Grid]) -> pp.ad.SparseArray:
                 subdomains: List of grids on which the vector field is defined.
 
             Returns:
-                Matrix extracting normal component of the vector field and expressing it
-                in normal basis. The size of the matrix is `(Nc, Nc * self.nd)`, where
-                `Nc` is the total number of cells in the subdomains.
+                Projection extracting normal component of the vector field and
+                expressing it in normal basis. The size of the projection is `(Nc, Nc *
+                self.nd)`, where `Nc` is the total number of cells in the subdomains.
 
             """