Merge 4e28c49 into c020429

src/porepy/numerics/fv/fvutils.py

@@ -159,7 +159,9 @@ def pair_over_subfaces_nd(self, other):
 # ------------------------ End of class SubcellTopology ----------------------
 
 
-def compute_dist_face_cell(g, subcell_topology, eta):
+def compute_dist_face_cell(
+    g, subcell_topology, eta, eta_at_bnd=False, return_paired=True
+):
     """
     Compute vectors from cell centers continuity points on each sub-face.
 
@@ -190,16 +192,20 @@ def compute_dist_face_cell(g, subcell_topology, eta):
     cols += matrix_compression.rldecode(np.cumsum(blocksz) - blocksz[0], blocksz)
     eta_vec = eta * np.ones(subcell_topology.fno.size)
     # Set eta values to zero at the boundary
-    bnd = np.in1d(subcell_topology.fno, g.get_all_boundary_faces())
-    eta_vec[bnd] = 0
+    if not eta_at_bnd:
+        bnd = np.in1d(subcell_topology.fno, g.get_all_boundary_faces())
+        eta_vec[bnd] = 0
     cp = g.face_centers[:, subcell_topology.fno] + eta_vec * (
         g.nodes[:, subcell_topology.nno] - g.face_centers[:, subcell_topology.fno]
     )
     dist = cp - g.cell_centers[:, subcell_topology.cno]
 
     ind_ptr = np.hstack((np.arange(0, cols.size, dims), cols.size))
     mat = sps.csr_matrix((dist.ravel("F"), cols.ravel("F"), ind_ptr))
-    return subcell_topology.pair_over_subfaces(mat)
+    if return_paired:
+        return subcell_topology.pair_over_subfaces(mat)
+    else:
+        return mat
 
 
 def determine_eta(g):

src/porepy/numerics/fv/mpsa.py

@@ -533,6 +533,8 @@ def mpsa(
     eta=None,
     inverter=None,
     max_memory=None,
+    hf_disp=False,
+    hf_eta=None,
     **kwargs
 ):
     """
@@ -575,7 +577,13 @@ def mpsa(
             If the **estimated** memory need is larger than the provided
             threshold, the discretization will be split into an appropriate
             number of sub-calculations, using mpsa_partial().
-
+        hf_disp (bool) False: If true two matrices hf_cell, hf_bound is also returned such
+            that hf_cell * U + hf_bound * u_bound gives the reconstructed displacement
+            at the point on the face hf_eta. U is the cell centered displacement and
+            u_bound the boundary conditions
+        hf_eta (float) None: The point of displacment on the half-faces. hf_eta=0 gives the
+            displacement at the face centers while hf_eta=1 gives the displacements at
+            the nodes. If None is given, the continuity points eta will be used.
     Returns:
         scipy.sparse.csr_matrix (shape num_faces, num_cells): stress
             discretization, in the form of mapping from cell displacement to
@@ -591,6 +599,18 @@ def mpsa(
             Incorporation as a right hand side in linear system by
             multiplication with divergence operator.
             NOTE: The stresses are ordered facewise (first s_x_1, s_y_1 etc)
+        If hf_disp is True the following will also be returned
+        scipy.sparse.csr_matrix (g.dim*shape num_hfaces, g.dim*num_cells):
+            displacement reconstruction for the displacement at the half faces. This is
+            the contribution from the cell-center displacements.
+            NOTE: The half-face displacements are ordered cell wise
+            (U_x_0, U_x_1, ..., U_x_n, U_y0, U_y1, ...)
+        scipy.sparse.csr_matrix (g.dim*shape num_hfaces, g.dim*num_faces):
+            displacement reconstruction for the displacement at the half faces.
+            This is the contribution from the boundary conditions.
+            NOTE: The half-face displacements are ordered cell wise
+            (U_x_0, U_x_1, ..., U_x_n, U_y0, U_y1, ...)
+
 
     Example:
         # Set up a Cartesian grid
@@ -632,10 +652,26 @@ def mpsa(
         # entire grid.
         # TODO: We may want to estimate the memory need, and give a warning if
         # this seems excessive
-        stress, bound_stress = _mpsa_local(
-            g, constit, bound, eta=eta, inverter=inverter, robin_weight=robin_weight
-        )
+        if hf_disp:
+            stress, bound_stress, hf_cell, hf_bound = _mpsa_local(
+                g,
+                constit,
+                bound,
+                eta=eta,
+                inverter=inverter,
+                robin_weight=robin_weight,
+                hf_disp=hf_disp,
+                hf_eta=hf_eta,
+            )
+        else:
+            stress, bound_stress = _mpsa_local(
+                g, constit, bound, eta=eta, inverter=inverter, robin_weight=robin_weight
+            )
     else:
+        if hf_disp:
+            raise NotImplementedError(
+                "No half face reconstruction for partial discretization yet"
+            )
         # Estimate number of partitions necessary based on prescribed memory
         # usage
         peak_mem = _estimate_peak_memory_mpsa(g)
@@ -688,7 +724,10 @@ def mpsa(
             stress += loc_stress
             bound_stress += loc_bound_stress
 
-    return stress, bound_stress
+    if hf_disp:
+        return stress, bound_stress, hf_cell, hf_bound
+    else:
+        return stress, bound_stress
 
 
 def mpsa_partial(
@@ -814,7 +853,16 @@ def mpsa_partial(
     return stress_glob, bound_stress_glob, active_faces
 
 
-def _mpsa_local(g, constit, bound, eta=0, robin_weight=None, inverter="numba"):
+def _mpsa_local(
+    g,
+    constit,
+    bound,
+    eta=0,
+    robin_weight=None,
+    inverter="numba",
+    hf_disp=False,
+    hf_eta=None,
+):
     """
     Actual implementation of the MPSA W-method. To calculate the MPSA
     discretization on a grid, either call this method, or, to respect the
@@ -897,7 +945,9 @@ def _mpsa_local(g, constit, bound, eta=0, robin_weight=None, inverter="numba"):
 
     hook_igrad = hook * igrad
     # NOTE: This is the point where we expect to reach peak memory need.
-    del hook, igrad
+    del hook
+    if not hf_disp:
+        del igrad
 
     # Output should be on face-level (not sub-face)
     hf2f = fvutils.map_hf_2_f(
@@ -917,7 +967,19 @@ def _mpsa_local(g, constit, bound, eta=0, robin_weight=None, inverter="numba"):
     bound_stress = hf2f * hook_igrad * rhs_bound
     stress, bound_stress = _zero_neu_rows(g, stress, bound_stress, bound)
 
-    return stress, bound_stress
+    if hf_disp:
+        eta_at_bnd = True
+        if hf_eta is None:
+            hf_eta = eta
+            eta_at_bnd = False
+        dist_grad, cell_centers = reconstruct_displacement(
+            g, subcell_topology, hf_eta, eta_at_bnd=eta_at_bnd
+        )
+        hf_cell = dist_grad * igrad * rhs_cells + cell_centers
+        hf_bound = dist_grad * igrad * rhs_bound
+        return stress, bound_stress, hf_cell, hf_bound
+    else:
+        return stress, bound_stress
 
 
 def mpsa_elasticity(
@@ -1048,6 +1110,48 @@ def mpsa_elasticity(
     return hook, igrad, rhs_cells, cell_node_blocks, hook_normal
 
 
+def reconstruct_displacement(g, subcell_topology, eta=None, eta_at_bnd=False):
+    """
+    Function for reconstructing the displacement at the half faces given the
+    local gradients. reconstruct_displacement(...) construct a matrix that
+    includes the distances from the cell centers to the evaluation point eta.
+    """
+    if eta is None:
+        eta = fvutils.determine_eta(g)
+
+    # Calculate the distance from the cell centers to continuity points
+    D_g = fvutils.compute_dist_face_cell(
+        g, subcell_topology, eta, eta_at_bnd=eta_at_bnd, return_paired=False
+    )
+    # expand indices to x-y-z
+    D_g = sps.kron(sps.eye(g.dim), D_g)
+    D_g = D_g.tocsr()
+
+    # Get a mapping from cell centers to half-faces
+    D_c = sps.coo_matrix(
+        (
+            np.ones(subcell_topology.subhfno.size),
+            (subcell_topology.subhfno, subcell_topology.cno),
+        )
+    ).tocsr()
+    # Expand indices to x-y-z
+    D_c = sps.kron(sps.eye(g.dim), D_c)
+    D_c = D_c.tocsc()
+    # book keeping
+    cell_node_blocks, _ = matrix_compression.rlencode(
+        np.vstack((subcell_topology.cno, subcell_topology.nno))
+    )
+    num_sub_cells = cell_node_blocks[0].size
+    # The column ordering of the displacement equilibrium equations are
+    # formed as a Kronecker product of scalar equations. Bring them to the
+    # same form as that applied in the force balance equations
+    dist_grad, cell_centers = __rearange_columns_displacement_eqs(
+        D_g, D_c, num_sub_cells, g.dim
+    )
+
+    return dist_grad, cell_centers
+
+
 # -----------------------------------------------------------------------------
 #
 # Below here are helper functions, which tend to be less than well documented.

test/integration/test_mpsa_basis_transformation.py

@@ -9,6 +9,7 @@
 import porepy as pp
 
 
+
 class TestMpsaBoundRhs(unittest.TestCase):
     def test_neu(self):
         g = pp.StructuredTriangleGrid([1, 1])

test/integration/test_mpsa_gradient_reconstruction.py

@@ -0,0 +1,155 @@
+import numpy as np
+import scipy.sparse as sps
+import unittest
+
+import porepy as pp
+
+
+class MpsaReconstructDisplacement(unittest.TestCase):
+    def test_cart_2d(self):
+        """
+        Test that mpsa gives out the correct matrices for
+        reconstruction of the displacement at the faces
+        """
+        nx = 1
+        ny = 1
+        g = pp.CartGrid([nx, ny], physdims=[2, 2])
+        g.compute_geometry()
+
+        lam = np.array([1])
+        mu = np.array([2])
+        k = pp.FourthOrderTensor(g.dim, mu, lam)
+
+        bc = pp.BoundaryCondition(g)
+        _, _, grad_cell, grad_bound = pp.numerics.fv.mpsa.mpsa(
+            g, k, bc, hf_disp=True, inverter="python"
+        )
+
+        grad_bound_known = np.array(
+            [
+                [0.10416667, 0., 0., 0., 0., -0.02083333, 0., 0.],
+                [0., 0., 0., 0., 0.25, 0., 0., 0.],
+                [0., 0., 0.10416667, 0., 0., 0.02083333, 0., 0.],
+                [0., 0., 0., 0., 0.25, 0., 0., 0.],
+                [0.10416667, 0., 0., 0., 0., 0., 0., 0.02083333],
+                [0., 0., 0., 0., 0., 0., 0.25, 0.],
+                [0., 0., 0.10416667, 0., 0., 0., 0., -0.02083333],
+                [0., 0., 0., 0., 0., 0., 0.25, 0.],
+                [0., 0.25, 0., 0., 0., 0., 0., 0.],
+                [-0.02083333, 0., 0., 0., 0., 0.10416667, 0., 0.],
+                [0., 0., 0., 0.25, 0., 0., 0., 0.],
+                [0., 0., 0.02083333, 0., 0., 0.10416667, 0., 0.],
+                [0., 0.25, 0., 0., 0., 0., 0., 0.],
+                [0.02083333, 0., 0., 0., 0., 0., 0., 0.10416667],
+                [0., 0., 0., 0.25, 0., 0., 0., 0.],
+                [0., 0., -0.02083333, 0., 0., 0., 0., 0.10416667],
+            ]
+        )
+
+        grad_cell_known = np.array(
+            [
+                [1., 0.],
+                [1., 0.],
+                [1., 0.],
+                [1., 0.],
+                [1., 0.],
+                [1., 0.],
+                [1., 0.],
+                [1., 0.],
+                [0., 1.],
+                [0., 1.],
+                [0., 1.],
+                [0., 1.],
+                [0., 1.],
+                [0., 1.],
+                [0., 1.],
+                [0., 1.],
+            ]
+        )
+
+        self.assertTrue(np.all(np.abs(grad_bound - grad_bound_known) < 1e-7))
+        self.assertTrue(np.all(np.abs(grad_cell - grad_cell_known) < 1e-12))
+
+    def test_simplex_3d_dirichlet(self):
+        """
+        Test that we retrieve a linear solution exactly
+        """
+        nx = 2
+        ny = 2
+        nz = 2
+        g = pp.StructuredTetrahedralGrid([nx, ny, nz], physdims=[1, 1, 1])
+        g.compute_geometry()
+
+        np.random.seed(2)
+
+        lam = np.ones(g.num_cells)
+        mu = np.ones(g.num_cells)
+        k = pp.FourthOrderTensor(g.dim, mu, lam)
+
+        s_t = pp.fvutils.SubcellTopology(g)
+
+        bc = pp.BoundaryConditionVectorial(g)
+        bc.is_dir[:, g.get_all_boundary_faces()] = True
+        bc.is_neu[bc.is_dir] = False
+
+        x0 = np.array([[1, 2, 3]]).T
+        u_b = g.face_centers + x0
+
+        stress, bound_stress, grad_cell, grad_bound = pp.numerics.fv.mpsa.mpsa(
+            g, k, bc, eta=0, hf_disp=True, inverter="python"
+        )
+
+        div = pp.fvutils.vector_divergence(g)
+
+        U = sps.linalg.spsolve(div * stress, -div * bound_stress * u_b.ravel("F"))
+
+        U_hf = (grad_cell * U + grad_bound * u_b.ravel("F")).reshape((g.dim, -1))
+
+        _, IA = np.unique(s_t.fno, True)
+        U_f = U_hf[:, IA]
+
+        U = U.reshape((g.dim, -1), order="F")
+
+        self.assertTrue(np.all(np.abs(U - g.cell_centers - x0) < 1e-10))
+        self.assertTrue(np.all(np.abs(U_f - g.face_centers - x0) < 1e-10))
+
+    def test_simplex_3d_boundary(self):
+        """
+        Even if we do not get exact solution at interiour we should be able to
+        retrieve the boundary conditions
+        """
+        nx = 2
+        ny = 2
+        nz = 2
+        g = pp.StructuredTetrahedralGrid([nx, ny, nz], physdims=[1, 1, 1])
+        g.compute_geometry()
+
+        np.random.seed(2)
+
+        lam = 10 * np.random.rand(g.num_cells)
+        mu = 10 * np.random.rand(g.num_cells)
+        k = pp.FourthOrderTensor(g.dim, mu, lam)
+
+        s_t = pp.fvutils.SubcellTopology(g)
+
+        bc = pp.BoundaryConditionVectorial(g)
+        dir_ind = g.get_all_boundary_faces()[[0, 2, 5, 8, 10, 13, 15, 21]]
+        bc.is_dir[:, dir_ind] = True
+        bc.is_neu[bc.is_dir] = False
+
+        u_b = np.random.randn(g.face_centers.shape[0], g.face_centers.shape[1])
+
+        stress, bound_stress, grad_cell, grad_bound = pp.numerics.fv.mpsa.mpsa(
+            g, k, bc, eta=0, hf_disp=True, inverter="python"
+        )
+
+        div = pp.fvutils.vector_divergence(g)
+
+        U = sps.linalg.spsolve(div * stress, -div * bound_stress * u_b.ravel("F"))
+
+        U_hf = (grad_cell * U + grad_bound * u_b.ravel("F")).reshape((g.dim, -1))
+
+        _, IA = np.unique(s_t.fno, True)
+        U_f = U_hf[:, IA]
+
+        self.assertTrue(np.all(np.abs(U_f[:, dir_ind] - u_b[:, dir_ind]) < 1e-10))