Merge pull request pmgbergen#285 from pmgbergen/primal_contact

Contact mechanics

examples/example5/test_upwind_examples.py

@@ -38,6 +38,8 @@ def test_upwind_example0(self, if_export=False):
         time_step = advect.cfl(g, data)
         data[pp.PARAMETERS]["transport"]["time_step"] = time_step
 
+        advect.discretize(g, data)
+
         U, rhs = advect.assemble_matrix_rhs(g, data)
         rhs = time_step * rhs
         U = time_step * U
@@ -92,6 +94,7 @@ def test_upwind_example1(self, if_export=False):
         time_step = advect.cfl(g, data)
         data[pp.PARAMETERS]["transport"]["time_step"] = time_step
 
+        advect.discretize(g, data)
         U, rhs = advect.assemble_matrix_rhs(g, data)
         rhs = time_step * rhs
         U = time_step * U
@@ -168,6 +171,7 @@ def funp_ex(pt):
         }
         data = pp.initialize_default_data(g, {}, "flow", specified_parameters)
         solver = pp.MVEM("flow")
+        solver.discretize(g, data)
         D_flow, b_flow = solver.assemble_matrix_rhs(g, data)
 
         solver_source = pp.DualScalarSource("flow")
@@ -194,6 +198,7 @@ def funp_ex(pt):
 
         # Advect solver
         advect = pp.Upwind("transport")
+        advect.discretize(g, data)
 
         U, rhs = advect.assemble_matrix_rhs(g, data)
         time_step = advect.cfl(g, data)

src/porepy/__init__.py

@@ -26,12 +26,14 @@
 __all__ = []
 
 # Numerics
+from porepy.numerics.discretization import VoidDiscretization
+
 # Control volume, elliptic
-from porepy.numerics.fv.mpsa import Mpsa, FracturedMpsa
+from porepy.numerics.fv.mpsa import Mpsa
 from porepy.numerics.fv.fv_elliptic import FVElliptic
 from porepy.numerics.fv.tpfa import Tpfa
 from porepy.numerics.fv.mpfa import Mpfa
-from porepy.numerics.fv.biot import Biot, GradP, DivD, BiotStabilization
+from porepy.numerics.fv.biot import Biot, GradP, DivU, BiotStabilization
 from porepy.numerics.fv.source import ScalarSource
 
 # Virtual elements, elliptic
@@ -50,9 +52,8 @@
 from porepy.numerics.interface_laws.elliptic_interface_laws import (
     RobinCoupling,
     FluxPressureContinuity,
-    RobinContact,
-    StressDisplacementContinuity,
 )
+
 from porepy.numerics.interface_laws.cell_dof_face_dof_map import CellDofFaceDofMap
 from porepy.numerics.mixed_dim.assembler import Assembler
 
@@ -64,6 +65,16 @@
 from porepy.numerics.fv.mass_matrix import MassMatrix
 from porepy.numerics.fv.mass_matrix import InvMassMatrix
 
+# Contact mechanics
+from porepy.numerics.interface_laws.contact_mechanics_interface_laws import (
+    PrimalContactCoupling,
+    DivUCoupling,
+    MatrixScalarToForceBalance,
+    FractureScalarToForceBalance,
+)
+from porepy.numerics.contact_mechanics.contact_conditions import ColoumbContact
+from porepy.numerics.contact_mechanics import contact_conditions
+
 # Grids
 from porepy.grids.grid import Grid
 from porepy.grids.fv_sub_grid import FvSubGrid
@@ -114,7 +125,8 @@
 from porepy.fracs import meshing, fracture_importer, mortars
 from porepy.grids import structured, simplex, coarsening, partition, refinement
 from porepy.numerics.fv import fvutils
-from porepy.utils import error
+from porepy.utils import error, grid_utils
+from porepy.utils.tangential_normal_projection import TangentialNormalProjection
 
 # Constants, units and keywords
 from porepy.utils.common_constants import *

src/porepy/geometry/map_geometry.py

@@ -31,13 +31,18 @@ def force_point_collinearity(pts):
     return pts[:, 0, np.newaxis] * (1 - dist) + pts[:, end, np.newaxis] * dist
 
 
-def map_grid(g, tol=1e-5):
+def map_grid(g, tol=1e-5, R=None):
     """ If a 2d or a 1d grid is passed, the function return the cell_centers,
     face_normals, and face_centers using local coordinates. If a 3d grid is
     passed nothing is applied. The return vectors have a reduced number of rows.
 
     Parameters:
     g (grid): the grid.
+    tol (double, optional): Tolerance used to check that the grid is linear or planar.
+        Defaults to 1e-5.
+    R (np.array size 3x3, optional ): Rotation matrix. The first dim rows should map
+        vectors onto the tangential space of the grid. If not provided, a rotation
+        matrix will be computed.
 
     Returns:
     cell_centers: (g.dim x g.num_cells) the mapped centers of the cells.
@@ -52,9 +57,11 @@ def map_grid(g, tol=1e-5):
     face_normals = g.face_normals
     face_centers = g.face_centers
     nodes = g.nodes
-    R = np.eye(3)
 
     if g.dim == 0 or g.dim == 3:
+        if R is None:
+            R = np.eye(3)
+
         return (
             cell_centers,
             face_normals,
@@ -64,12 +71,12 @@ def map_grid(g, tol=1e-5):
             nodes,
         )
 
-    if g.dim == 1 or g.dim == 2:
-
-        if g.dim == 2:
-            R = project_plane_matrix(g.nodes, tol=tol)
-        else:
-            R = project_line_matrix(g.nodes, tol=tol)
+    else:  # g.dim == 1 or g.dim == 2:
+        if R is None:
+            if g.dim == 2:
+                R = project_plane_matrix(g.nodes, tol=tol)
+            else:
+                R = project_line_matrix(g.nodes, tol=tol)
 
         face_centers = np.dot(R, face_centers)
 

src/porepy/grids/grid.py

@@ -743,6 +743,34 @@ def cell_connection_map(self):
 
         return c2c
 
+    def sign_of_faces(self, faces):
+        """ Get the direction of the normal vector (inward or outwards from a cell)
+        of faces. Only boundary faces are permissible.
+
+        Parameters:
+            faces: (ndarray) indices of faces that you want to know the sign for. The
+                faces must be boundary faces.
+
+        Returns:
+            (ndarray) the sign of the faces
+
+        Raises:
+            ValueError if a target face is internal.
+
+        """
+
+        IA = np.argsort(faces)
+        IC = np.argsort(IA)
+
+        fi, _, sgn = sps.find(self.cell_faces[faces[IA], :])
+        if fi.size != faces.size:
+            raise ValueError("sign of internal faces does not make sense")
+
+        I = np.argsort(fi)
+        sgn = sgn[I]
+        sgn = sgn[IC]
+        return sgn
+
     def bounding_box(self):
         """
         Return the bounding box of the grid.