Compliance to BCs and new face tags for the files merged from develop…

… in previous commit

src/porepy/fracs/non_conforming.py

@@ -610,7 +610,7 @@ def update_face_tags(g, delete_faces, new_faces):
 
     Delete tags for old faces, and add new tags for their replacements.
 
-    If the grid has no face_tags, no change is done
+    If the grid has no face tags, no change is done
 
     Parameters:
         g (grid): To be modified

src/porepy/numerics/fv/mpsa.py

@@ -14,7 +14,6 @@
 from porepy.numerics.fv import fvutils
 from porepy.utils import matrix_compression, mcolon, sparse_mat
 from porepy.grids import structured, partition
-from porepy.grids.grid import FaceTag
 from porepy.params import tensor, bc
 from porepy.numerics.mixed_dim.solver import Solver
 
@@ -129,11 +128,13 @@ def rhs(self, g, bound_stress, bc_val, f):
 
 #------------------------------------------------------------------------------#
 
+
 class FracturedMpsa(Mpsa):
     """
     Subclass of MPSA for discretizing a fractured domain. Adds DOFs on each
     fracture face which describe the fracture deformation.
     """
+
     def __init__(self, **kwargs):
         Mpsa.__init__(self, **kwargs)
         assert hasattr(self, 'physics'), 'Mpsa must assign physics'
@@ -152,10 +153,9 @@ def ndof(self, g):
         dof: the number of degrees of freedom.
 
         """
-        num_fracs = np.sum(g.has_face_tag(FaceTag.FRACTURE))
+        num_fracs = np.sum(g.tags['fracture_faces'])
         return g.dim * (g.num_cells + num_fracs)
 
-
     def matrix_rhs(self, g, data, discretize=True):
         """
         Return the matrix and right-hand side for a discretization of a second
@@ -191,14 +191,14 @@ def matrix_rhs(self, g, data, discretize=True):
         L, b_l = self.given_slip_distance(g, stress, bound_stress)
 
         bc_val = data['param'].get_bc_val(self)
-        
-        frac_faces = np.matlib.repmat(g.has_face_tag(FaceTag.FRACTURE), 3, 1)
+
+        frac_faces = np.matlib.repmat(g.tags['fracture_faces'], 3, 1)
         assert np.all(bc_val[frac_faces.ravel('F')] == 0), \
             '''Fracture should have zero boundary condition. Set slip by
                Parameters.set_slip_distance'''
 
-        slip_distance = data['param'].get_slip_distance()        
-        
+        slip_distance = data['param'].get_slip_distance()
+
         A = sps.vstack((A_e, L), format='csr')
         rhs = np.hstack((b_e * bc_val, b_l * slip_distance))
 
@@ -219,13 +219,14 @@ def traction(self, g, data, sol):
         T : array (g.dim * g.num_faces)
             traction on each face
 
-        """        
+        """
         bc_val = data['param'].get_bc_val(self.physics).copy()
         frac_disp = self.extract_frac_u(g, sol)
         cell_disp = self.extract_u(g, sol)
 
         frac_faces = (g.frac_pairs).ravel('C')
-        frac_ind = mcolon.mcolon(g.dim * frac_faces, g.dim * frac_faces + g.dim)
+        frac_ind = mcolon.mcolon(
+            g.dim * frac_faces, g.dim * frac_faces + g.dim)
 
         bc_val[frac_ind] = frac_disp
 
@@ -295,20 +296,18 @@ def discretize_fractures(self, g, data, faces=None, **kwargs):
         data: dictionary to store the data.
         """
 
-
-
         #    dir_bound = g.get_boundary_faces()
         #    bound = bc.BoundaryCondition(g, dir_bound, ['dir'] * dir_bound.size)
 
-        frac_faces = g.has_face_tag(FaceTag.FRACTURE)
+        frac_faces = g.tags['fracture_faces']
 
         bound = data['param'].get_bc(self)
-        is_dir  = bound.is_dir
+        is_dir = bound.is_dir
         if not np.all(is_dir[frac_faces]):
             is_dir[frac_faces] = True
 
         bound = bc.BoundaryCondition(g, is_dir, 'dir')
-        
+
         assert np.all(bound.is_dir[frac_faces]), \
             'fractures must be set as dirichlet boundary faces'
 
@@ -361,14 +360,13 @@ def discretize_fractures(self, g, data, faces=None, **kwargs):
         d_b = -b_external
         # sps.csr_matrix((int_b_left.size, g.num_faces * g.dim))
         d_t = -sgn_left * bound_stress_external[int_b_left] \
-             - sgn_right * bound_stress_external[int_b_right]
+            - sgn_right * bound_stress_external[int_b_right]
 
         b_matrix = sps.vstack((d_b, d_t), format='csr')
-                
+
         data['b_e'] = b_matrix
         data['A_e'] = A
 
-
     def given_traction(self, g, stress, bound_stress, faces=None, **kwargs):
         # we find the matrix indices of the fracture
         if faces is None:
@@ -533,8 +531,9 @@ def mpsa(g, constit, bound, eta=None, inverter=None, max_memory=None,
         # For the moment nothing to do here, just call main mpfa method for the
         # 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)
+        # this seems excessive
+        stress, bound_stress = _mpsa_local(
+            g, constit, bound, eta=eta, inverter=inverter)
     else:
         # Estimate number of partitions necessary based on prescribed memory
         # usage
@@ -549,7 +548,7 @@ def mpsa(g, constit, bound, eta=None, inverter=None, max_memory=None,
         # Empty fields for stress and bound_stress. Will be expanded as we go.
         # Implementation note: It should be relatively straightforward to
         # estimate the memory need of stress (face_nodes -> node_cells ->
-        # unique). 
+        # unique).
         stress = sps.csr_matrix((g.num_faces * g.dim, g.num_cells * g.dim))
         bound_stress = sps.csr_matrix((g.num_faces * g.dim,
                                        g.num_faces * g.dim))
@@ -570,7 +569,7 @@ def mpsa(g, constit, bound, eta=None, inverter=None, max_memory=None,
 
             # Perform local discretization.
             loc_stress, loc_bound_stress, loc_faces \
-                = mpsa_partial(g, constit, bound, eta=eta, inverter=inverter, 
+                = mpsa_partial(g, constit, bound, eta=eta, inverter=inverter,
                                nodes=active_nodes)
 
             # Eliminate contribution from faces already covered
@@ -586,7 +585,7 @@ def mpsa(g, constit, bound, eta=None, inverter=None, max_memory=None,
     return stress, bound_stress
 
 
-def mpsa_partial(g, constit, bound, eta=0, inverter='numba', cells=None, 
+def mpsa_partial(g, constit, bound, eta=0, inverter='numba', cells=None,
                  faces=None, nodes=None):
     """
     Run an MPFA discretization on subgrid, and return discretization in terms
@@ -669,7 +668,7 @@ def mpsa_partial(g, constit, bound, eta=0, inverter='numba', cells=None,
 
     # Discretization of sub-problem
     stress_loc, bound_stress_loc = _mpsa_local(sub_g, loc_c, loc_bnd,
-                                           eta=eta, inverter=inverter)
+                                               eta=eta, inverter=inverter)
 
     face_map, cell_map = fvutils.map_subgrid_to_grid(g, l2g_faces, l2g_cells,
                                                      is_vector=True)
@@ -688,6 +687,7 @@ def mpsa_partial(g, constit, bound, eta=0, inverter='numba', cells=None,
 
     return stress_glob, bound_stress_glob, active_faces
 
+
 def _mpsa_local(g, constit, bound, eta=0, inverter='numba'):
     """
     Actual implementation of the MPSA W-method. To calculate the MPSA
@@ -769,7 +769,7 @@ def _mpsa_local(g, constit, bound, eta=0, inverter='numba'):
 
     # Output should be on face-level (not sub-face)
     hf2f = fvutils.map_hf_2_f(subcell_topology.fno_unique,
-                   	      subcell_topology.subfno_unique, nd)
+                              subcell_topology.subfno_unique, nd)
 
     # Stress discretization
     stress = hf2f * hook_igrad * rhs_cells
@@ -784,7 +784,7 @@ def _mpsa_local(g, constit, bound, eta=0, inverter='numba'):
 
 
 def mpsa_elasticity(g, constit, subcell_topology, bound_exclusion, eta,
-                      inverter):
+                    inverter):
     """
     This is the function where the real discretization takes place. It contains
     the parts that are common for elasticity and poro-elasticity, and was thus
@@ -884,7 +884,7 @@ def _estimate_peak_memory_mpsa(g):
     nd = g.dim
     num_cell_nodes = g.cell_nodes().sum(axis=1).A
 
-    # Number of unknowns around a vertex: nd^2 per cell that share the vertex 
+    # Number of unknowns around a vertex: nd^2 per cell that share the vertex
     # for pressure gradients, and one per cell (cell center pressure)
     num_grad_unknowns = nd**2 * num_cell_nodes
 
@@ -963,7 +963,7 @@ def _split_stiffness_matrix(constit):
     # necessary
     csym = 0 * constit.copy().c
     casym = constit.copy().c
-    
+
     # The copy constructor for the stiffness matrix will represent all
     # dimensions as 3d. If dim==2, delete the redundant rows and columns
     if dim == 2 and csym.shape[0] == 9:
@@ -1502,4 +1502,3 @@ def _sign_matrix(g, faces):
     sgn = sps.diags(sgn_d, 0)
 
     return sgn
-