From bfe25d34f76f936e727013c4351f6e7c5c6057c6 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 2 Jul 2019 16:27:03 +0200 Subject: [PATCH] Deleted old implementation of contact mechanics. --- src/porepy/__init__.py | 3 +- src/porepy/numerics/fv/biot.py | 28 +- src/porepy/numerics/fv/mpsa.py | 26 +- .../interface_laws/elliptic_discretization.py | 293 -------- .../interface_laws/elliptic_interface_laws.py | 675 ------------------ .../test_elliptic_interface_laws.py | 505 ------------- test/unit/test_elliptic_interface_laws.py | 346 --------- 7 files changed, 46 insertions(+), 1830 deletions(-) delete mode 100644 test/integration/test_elliptic_interface_laws.py delete mode 100644 test/unit/test_elliptic_interface_laws.py diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py index 1cc6dea781..fb3c7e6f6f 100644 --- a/src/porepy/__init__.py +++ b/src/porepy/__init__.py @@ -52,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 diff --git a/src/porepy/numerics/fv/biot.py b/src/porepy/numerics/fv/biot.py index 799a006e80..faee14700f 100644 --- a/src/porepy/numerics/fv/biot.py +++ b/src/porepy/numerics/fv/biot.py @@ -792,11 +792,29 @@ def compute_stress(self, g, u, data): return stress -class GradP( - pp.numerics.interface_laws.elliptic_discretization.VectorEllipticDiscretization -): +class GradP(): """ Class for the pressure gradient term of the Biot equation. """ + + def __init__(self, keyword): + """ Set the discretization, with the keyword used for storing various + information associated with the discretization. + + Paramemeters: + keyword (str): Identifier of all information used for this + discretization. + """ + self.keyword = keyword + + def _key(self): + """ Get the keyword of this object, on a format friendly to access relevant + fields in the data dictionary + + Returns: + String, on the form self.keyword + '_'. + + """ + return self.keyword + "_" def ndof(self, g): """ Return the number of degrees of freedom associated to the method. @@ -974,9 +992,7 @@ def enforce_neumann_int_bound(self, *_): pass -class DivU( - pp.numerics.interface_laws.elliptic_discretization.VectorEllipticDiscretization -): +class DivU(): """ Class for the displacement divergence term of the Biot equation. """ diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index cabce8af2c..6bc69a1e63 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -19,9 +19,29 @@ logger = logging.getLogger(__name__) -class Mpsa( - pp.numerics.interface_laws.elliptic_discretization.VectorEllipticDiscretization -): +class Mpsa(): + + + def __init__(self, keyword): + """ Set the discretization, with the keyword used for storing various + information associated with the discretization. + + Paramemeters: + keyword (str): Identifier of all information used for this + discretization. + """ + self.keyword = keyword + + def _key(self): + """ Get the keyword of this object, on a format friendly to access relevant + fields in the data dictionary + + Returns: + String, on the form self.keyword + '_'. + + """ + return self.keyword + "_" + def ndof(self, g): """ Return the number of degrees of freedom associated to the method. diff --git a/src/porepy/numerics/interface_laws/elliptic_discretization.py b/src/porepy/numerics/interface_laws/elliptic_discretization.py index eabc74b9f8..d942091628 100644 --- a/src/porepy/numerics/interface_laws/elliptic_discretization.py +++ b/src/porepy/numerics/interface_laws/elliptic_discretization.py @@ -332,296 +332,3 @@ def enforce_neumann_int_bound( matrix (scipy.sparse.matrix): Discretization matrix to be modified. """ raise NotImplementedError("Method not implemented") - - -class VectorEllipticDiscretization: - """ Superclass for finite volume discretizations of the vector elliptic equation. - - Should not be used by itself, instead use a subclass that implements an - actual discretization method. Known subclass is Mpsa. - """ - - def __init__(self, keyword): - """ Set the discretization, with the keyword used for storing various - information associated with the discretization. - - Paramemeters: - keyword (str): Identifier of all information used for this - discretization. - """ - self.keyword = keyword - - def _key(self): - """ Get the keyword of this object, on a format friendly to access relevant - fields in the data dictionary - - Returns: - String, on the form self.keyword + '_'. - - """ - return self.keyword + "_" - - def ndof(self, g): - """ Abstract method. Return the number of degrees of freedom associated to the - method. - - Parameter: - g (grid): Computational grid - - Returns: - int: the number of degrees of freedom. - - """ - raise NotImplementedError("Method not implemented") - - def extract_displacement(self, g, solution_array, d): - """ Abstract method. Extract the displacement part of a solution. - - The implementation will depend on what the primary variables of the specific - implementation are. - - Parameters: - g (grid): To which the solution array belongs. - solution_array (np.array): Solution for this grid obtained from - either a mono-dimensional or a mixed-dimensional problem. - d (dictionary): Data dictionary associated with the grid. Not used, - but included for consistency reasons. - - Returns: - np.array (g.num_cells): Displacement solution vector. - - Raises: - NotImplementedError if the method - """ - raise NotImplementedError("Method not implemented") - - def extract_traction(self, g, solution_array, d): - """ Abstract method. Extract the traction part of a solution. - - The implementation will depend on what the primary variables of the specific - implementation are. - - Parameters: - g (grid): To which the solution array belongs. - solution_array (np.array): Solution for this grid obtained from - either a mono-dimensional or a mixed-dimensional problem. Will - correspond to the displacement solution. - d (dictionary): Data dictionary associated with the grid. - - Returns: - np.array (g.num_faces): Traction vector. - - """ - raise NotImplementedError() - - # ------------------------------------------------------------------------------# - - def assemble_matrix_rhs(self, g, data): - """ Return the matrix and right-hand side for a discretization of a second - order elliptic vector equation. - - - - Parameters: - g : grid, or a subclass, with geometry fields computed. - data: dictionary to store the data. For details on necessary keywords, - see method discretize() - discretize (boolean, optional): default True. Whether to discetize prior to - matrix assembly. If False, data should already contain discretization. - - Returns: - matrix: sparse csr (self.ndof, self.ndof) discretization matrix. - rhs: np.ndarray (self.ndof) Right-hand side which contains the boundary - conditions and the vector source term. - """ - raise NotImplementedError("Method not implemented") - - def assemble_matrix(self, g, data): - """Return the matrix for a discretization of a second order elliptic vector - equation. - - Parameters: - g (Grid): Computational grid, with geometry fields computed. - data (dictionary): With data stored. - - Returns: - scipy.sparse.csr_matrix: System matrix of this discretization. The - size of the matrix will depend on the specific discretization. - """ - raise NotImplementedError("Method not implemented") - - # ------------------------------------------------------------------------------# - - def assemble_rhs(self, g, data): - """ Return the right-hand side for a discretization of a second order elliptic - vector equation. - - Also discretize the necessary operators if the data dictionary does not - contain a discretization of the boundary equation. - - Parameters: - g (Grid): Computational grid, with geometry fields computed. - data (dictionary): With data stored. - - Returns: - np.ndarray: Right hand side vector with representation of boundary - conditions. The size of the vector will depend on the discretization. - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_stress( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """Assemble the contribution from an internal boundary, manifested as a - stress boundary condition. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a higher-dimensional domain. - - Implementations of this method will use an interplay between the grid - on the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_source( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """ Abstract method. Assemble the contribution from an internal - boundary, manifested as a body force term. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a lower-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_displacement_trace( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """ Abstract method. Assemble the contribution from an internal - boundary, manifested as a condition on the boundary displacement. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a higher-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_displacement_cell( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """ Abstract method. Assemble the contribution from an internal - boundary, manifested as a condition on the cell displacement. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a lower-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - """ - raise NotImplementedError("Method not implemented") - - def enforce_neumann_int_bound( - self, g_master, data_edge, matrix, swap_grid, self_ind - ): - """ Enforce Neumann boundary conditions on a given system matrix. - - The method is void for finite volume approaches, but is implemented - to be compatible with the general framework. - - Parameters: - g (Grid): On which the equation is discretized - data (dictionary): Of data related to the discretization. - matrix (scipy.sparse.matrix): Discretization matrix to be modified. - """ - raise NotImplementedError("Method not implemented") diff --git a/src/porepy/numerics/interface_laws/elliptic_interface_laws.py b/src/porepy/numerics/interface_laws/elliptic_interface_laws.py index c8059b7a7e..d9fde73923 100644 --- a/src/porepy/numerics/interface_laws/elliptic_interface_laws.py +++ b/src/porepy/numerics/interface_laws/elliptic_interface_laws.py @@ -395,678 +395,3 @@ def assemble_matrix_rhs( ) return matrix, rhs - - -# ------------------------------------------------------------------------------ - - -class RobinContact(object): - """ - Contact condition for elastic problem. This condition defines a Robin condition - for the stress and displacement jump between slave and master boundaries. g_slave - and g_master must have the same dimension. - - The contact condition is Newton's third law - \sigma \cdot n_slave = -\sigma \cdot n_master, - i.e., traction on the two sides must be equal and oposite, and a Robin-type condition - on the displacement jump - MW * \lambda + RW [u] = robin_rhs - where MW and RW are matrices of size (g_slave.dim, g.slave.dim), and - \labmda = \sigma \cdot \n_slave. - The jump operator [\cdot] is given by - [v] = v_slave - v_master, - and robin_rhs is a given rhs. - """ - - def __init__(self, keyword, discr_master, discr_slave=None): - self.keyword = keyword - if discr_slave is None: - discr_slave = discr_master - self.discr_master = discr_master - self.discr_slave = discr_slave - - def _key(self): - return self.keyword + "_" - - def _discretization_key(self): - return self._key() + pp.keywords.DISCRETIZATION - - def ndof(self, mg): - return (mg.dim + 1) * mg.num_cells - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ - Discretize the Mortar coupling. - We assume the following two sub-dictionaries to be present in the data_edge - dictionary: - parameter_dictionary, storing all parameters. - Stored in data[pp.PARAMETERS][self.keyword]. - matrix_dictionary, for storage of discretization matrices. - Stored in data[pp.DISCRETIZATION_MATRICES][self.keyword] - - parameter_dictionary contains the entries: - robin_weigth (list): a list of mortar_grid.num_cells np.ndarrays of - shape (mortar_grid.dim + 1, mortar_grid.dim + 1) giving the displacement - jump weight. - mortar_weigth (list): a list of mortar_grid.num_cells np.ndarrays of - shape (mortar_grid.dim + 1, mortar_grid.dim + 1) giving the mortar - - matrix_dictionary will be updated with the following entries: - mortar_weigth: sps.csc_matrix (mg.num_cells * mg.dim, mg.num_cells * mg.dim) - The weight matrix applied to the mortar variables. - robin_weigth: sps.csc_matrix (mg.num_cells * mg.dim, mg.num_cells * mg.dim) - The weight matrix applied to the displacement jump. - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - parameter_dictionary_edge = data_edge[pp.PARAMETERS][self.keyword] - - mortar_weight = sps.block_diag(parameter_dictionary_edge["mortar_weight"]) - robin_weight = sps.block_diag(parameter_dictionary_edge["robin_weight"]) - robin_rhs = parameter_dictionary_edge["robin_rhs"] - matrix_dictionary_edge["mortar_weight"] = mortar_weight - matrix_dictionary_edge["robin_weight"] = robin_weight - matrix_dictionary_edge["robin_rhs"] = robin_rhs - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - - self.discretize(g_master, g_slave, data_master, data_slave, data_edge) - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - dof_master = self.discr_master.ndof(g_master) - dof_slave = self.discr_slave.ndof(g_slave) - - if not dof_master == matrix[master_ind, master_ind].shape[1]: - raise ValueError( - """The number of dofs of the master discretization given - in RobinContact must match the number of dofs given by the matrix - """ - ) - elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: - raise ValueError( - """The number of dofs of the slave discretization given - in RobinContact must match the number of dofs given by the matrix - """ - ) - elif not self.ndof(mg) == matrix[master_ind, 2].shape[1]: - raise ValueError( - """The number of dofs of the edge discretization given - in RobinContact must match the number of dofs given by the matrix - """ - ) - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - cc_mortar = cc_master.copy() - - # EK: For some reason, the following lines were necessary to apease python - rhs_slave = np.empty(3, dtype=np.object) - rhs_slave[master_ind] = np.zeros(dof_master) - rhs_slave[slave_ind] = np.zeros(dof_slave) - rhs_slave[2] = np.zeros(self.ndof(mg)) - # I got some problems with pointers when doing rhs_master = rhs_slave.copy() - # so just reconstruct everything. - rhs_master = np.empty(3, dtype=np.object) - rhs_master[master_ind] = np.zeros(dof_master) - rhs_master[slave_ind] = np.zeros(dof_slave) - rhs_master[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # Obtain the displacement trace u_master - self.discr_master.assemble_int_bound_displacement_trace( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # set \sigma_master = -\lamba - self.discr_master.assemble_int_bound_stress( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # Obtain the displacement trace u_slave - self.discr_slave.assemble_int_bound_displacement_trace( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # set \sigma_slave = \lamba - self.discr_slave.assemble_int_bound_stress( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - # Then we flip the sign for the master displacement since the displacement - # jump is defined as u_slave - u_master - cc_master[2, master_ind] = -cc_master[2, master_ind] - rhs_master[2] = -rhs_master[2] - # Note that cc_master[2, 2] is fliped twice, first in Newton's third law, - # then for the displacement jump. - - # now, the matrix cc = cc_slave + cc_master expresses the stress and displacement - # continuities over the mortar grid. - # cc[0] -> stress_master = mortar_stress - # cc[1] -> stress_slave = -mortar_stress - # cc[2] -> mortar_weight * lambda + robin_weight * (u_slave - u_master) = robin_rhs - - # We don't want to enforce the displacement jump, but a Robin condition. - # We therefore add the mortar variable to the last equation. - cc_mortar[2, 2] = matrix_dictionary_edge["mortar_weight"] - - # The displacement jump is scaled by a matrix in the Robin condition: - robin_weight = matrix_dictionary_edge["robin_weight"] - cc_sm = cc_master + cc_slave - rhs = rhs_slave + rhs_master - rhs[2] = robin_weight * rhs[2] - for i in range(3): - cc_sm[2, i] = robin_weight * cc_sm[2, i] - - # Now define the complete Robin condition: - # mortar_weight * \lambda + "robin_weight" * [u] = robin_rhs - matrix += cc_sm + cc_mortar - rhs[2] += matrix_dictionary_edge["robin_rhs"] - - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs - - -class StressDisplacementContinuity(RobinContact): - """ - Contact condition for elastic problem. This condition defines continuity for - the stress and displacement jump between slave and master boundaries. g_slave - and g_master must have the same dimension. - - This contact condition is equivalent as if the slave and master domain was - a single connected domain (the discrete solution will be different as the - discretization will be slightly different). - """ - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ Nothing really to do here - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - pass - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - ---------- - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - - self.discretize(g_master, g_slave, data_master, data_slave, data_edge) - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - dof_master = self.discr_master.ndof(g_master) - dof_slave = self.discr_slave.ndof(g_slave) - - if not dof_master == matrix[master_ind, master_ind].shape[1]: - raise ValueError( - """The number of dofs of the master discretization given - in FluxPressureContinuity must match the number of dofs given by the matrix - """ - ) - elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: - raise ValueError( - """The number of dofs of the slave discretization given - in FluxPressureContinuity must match the number of dofs given by the matrix - """ - ) - elif not self.ndof(mg) == matrix[master_ind, 2].shape[1]: - raise ValueError( - """The number of dofs of the edge discretization given - in FluxPressureContinuity must match the number of dofs given by the matrix - """ - ) - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - # The rhs is just zeros - # EK: For some reason, the following lines were necessary to apease python - rhs_slave = np.empty(3, dtype=np.object) - rhs_slave[master_ind] = np.zeros(dof_master) - rhs_slave[slave_ind] = np.zeros(dof_slave) - rhs_slave[2] = np.zeros(self.ndof(mg)) - # I got some problems with pointers when doing rhs_master = rhs_slave.copy() - # so just reconstruct everything. - rhs_master = np.empty(3, dtype=np.object) - rhs_master[master_ind] = np.zeros(dof_master) - rhs_master[slave_ind] = np.zeros(dof_slave) - rhs_master[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # Obtain the displacement trace u_master - self.discr_master.assemble_int_bound_displacement_trace( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # set \sigma_master = -\lamba - self.discr_master.assemble_int_bound_stress( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # Obtain the displacement trace u_slave - self.discr_slave.assemble_int_bound_displacement_trace( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # set \sigma_slave = \lamba - self.discr_slave.assemble_int_bound_stress( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - # Then we flip the sign for the master displacement since the displacement - # jump is defined as u_slave - u_master - cc_master[2, master_ind] = -cc_master[2, master_ind] - rhs_master[2] = -rhs_master[2] - # Note that cc_master[2, 2] is fliped twice, first in Newton's third law, - # then for the displacement jump. - - # now, the matrix cc = cc_slave + cc_master expresses the stress and displacement - # continuities over the mortar grid. - # cc[0] -> stress_master = mortar_stress - # cc[1] -> stress_slave = -mortar_stress - # cc[2] -> u_slave - u_master = 0 - - matrix += cc_master + cc_slave - rhs = rhs_slave + rhs_master - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - - # Consider this terms only if the grids are of the same dimension - if g_master.dim == g_slave.dim: - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs - - -class RobinContactBiotPressure(RobinContact): - """ - This condition adds the fluid pressure contribution to the Robin contact condition. - The Robin condition says: - MW * lambda + RW * [u] = robin_rhs, - where MW (mortar_weight) and RW (robin_weight) are two matrices, and - [u] = u_slave - u_master is the displacement jump from the slave to the master. - In Biot the displacement on the contact boundary (u_slave and u_master) will be a - linear function of cell center displacement (u), mortar stress (lambda) and cell - centere fluid pressure (p): - A * u + B * p + C * lam = u_slave/u_master - This class adds the contribution B. - - To enforce the full continuity this interface law must be used in combination with - the RobinContact conditions which adds the contributions A and C - """ - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ Discretize the robin weight (RW) - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - parameter_dictionary_edge = data_edge[pp.PARAMETERS][self.keyword] - - robin_weight = sps.block_diag(parameter_dictionary_edge["robin_weight"]) - matrix_dictionary_edge["robin_weight"] = robin_weight - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - ---------- - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - - self.discretize(g_master, g_slave, data_master, data_slave, data_edge) - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - - # The rhs is just zeros - # EK: For some reason, the following lines were necessary to apease python - rhs = np.empty(3, dtype=np.object) - rhs[master_ind] = np.zeros(matrix[master_ind, master_ind].shape[1]) - rhs[slave_ind] = np.zeros(matrix[slave_ind, slave_ind].shape[1]) - rhs[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # Obtain the contribution of the cell centered pressure on the displacement - # trace u_master - self.discr_master.assemble_int_bound_displacement_trace( - g_master, data_master, data_edge, False, cc_master, matrix, rhs, master_ind - ) - # Equivalent for u_slave - self.discr_slave.assemble_int_bound_displacement_trace( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - # Then we flip the sign for the master displacement since the displacement - # jump is defined as u_slave - u_master - cc_master[2, master_ind] = -cc_master[2, master_ind] - - matrix += cc_master + cc_slave - - # The displacement jump is scaled by a matrix in the Robin condition: - robin_weight = matrix_dictionary_edge["robin_weight"] - - for i in range(3): - matrix[2, i] = robin_weight * matrix[2, i] - - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs - - -class DivU_StressMortar(RobinContactBiotPressure): - """ - This condition adds the stress mortar contribution to the div u term in the - fluid mass conservation equation of the Biot equations. When fractures are - present the divergence of u (div_u) will be a function of cell centere displacement, - boundary conditions and the stress mortar (lambda): - div_u = A * u + B * u_bc_val + C * lambda - The class adds the contribution C, while the DivU discretization adds A and B. - """ - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ Nothing really to do here - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - pass - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - ---------- - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - - # When we do time stepping in Biot the mortar variable from the previous - # time step will add a contribution to the rhs due to Backward Euler: - # \partial div_u / \partial_t = (\div_u^k - \div_u^{k-1})/dt. - rhs_slave = np.empty(3, dtype=np.object) - rhs_slave[master_ind] = np.zeros(matrix[master_ind, master_ind].shape[1]) - rhs_slave[slave_ind] = np.zeros(matrix[slave_ind, slave_ind].shape[1]) - rhs_slave[2] = np.zeros(self.ndof(mg)) - # I got some problems with pointers when doing rhs_master = rhs_slave.copy() - # so just reconstruct everything. - rhs_master = np.empty(3, dtype=np.object) - rhs_master[master_ind] = np.zeros(matrix[master_ind, master_ind].shape[1]) - rhs_master[slave_ind] = np.zeros(matrix[slave_ind, slave_ind].shape[1]) - rhs_master[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # lambda acts as a boundary condition on the div_u term. Assemble it for the master. - self.discr_master.assemble_int_bound_stress( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # Equivalent for the slave - self.discr_slave.assemble_int_bound_stress( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - rhs_master[master_ind] = -rhs_master[master_ind] - - matrix += cc_master + cc_slave - rhs = [s + m for s, m in zip(rhs_slave, rhs_master)] - - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs diff --git a/test/integration/test_elliptic_interface_laws.py b/test/integration/test_elliptic_interface_laws.py deleted file mode 100644 index d568d72bd6..0000000000 --- a/test/integration/test_elliptic_interface_laws.py +++ /dev/null @@ -1,505 +0,0 @@ -""" -Module for testing the vector elliptic couplings in the interface_laws. -""" - - -import numpy as np -import scipy.sparse as sps -import unittest - -import porepy as pp - - -class TestTwoGridCoupling(unittest.TestCase): - """ - In this test we set up a coupling between two grids of dimension 2: - g_slave g_master - |-----| | |------| - | | x | | - |-----| | |------| - g_mortar - There is one cell per grid and they are coupled together by a single mortar - variable. - - We define a random Robin and Mortar weights and test if we recover the - condition on the interface. - """ - - def test_robin_coupling(self): - """ - Test a Robin condition on the interface - """ - self.kw = "mech" - gb, _ = define_gb() - mortar_weight = np.random.rand(gb.dim_max()) - robin_weight = np.random.rand(gb.dim_max()) - rhs = np.random.rand(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb) - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - def test_continuity_coupling(self): - """ - Test a continuity condition on the interface. This is equivalent to - zero mortar weight and identity matrix for the robin weight. These - matrices are only used to check the solution. We check the solution - against a reference computed on a single grid including both g_s and g_m. - """ - self.kw = "mech" - gb, gb_full = define_gb() - # We assign weighs according to the condition. - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.zeros(gb.dim_max()) - # Assign data to coupling gb - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb, robin=False) - # Assign data to mono gb - self.assign_parameters(gb_full, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb_full, robin=False) - - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - assembler_full = pp.Assembler(gb_full) - matrix, rhs = assembler_full.assemble_matrix_rhs() - u_full = sps.linalg.spsolve(matrix, rhs) - # compare solutions - # We need to rearange the solutions because the ordering of the dofs are not the same - # Also, we don't have equality because the weak symmetry breaks when a cell has to many - # Neumann conditions (see comments in mpsa) - us = [] - ID = [] # to appease porpy 3.5 - for g, d in gb: - us.append(d[pp.STATE][self.kw]) - ID.append(g.grid_num - 1) - us = np.hstack([np.array(us)[ID].ravel()]) - IA = np.array([0, 1, 4, 5, 2, 3, 6, 7]) - sol = us[IA] - self.assertTrue(np.all(np.abs(sol - u_full) < 1e-4)) - - def check_solution(self, gb): - for e, d in gb.edges(): - mg = d["mortar_grid"] - gs, gm = gb.nodes_of_edge(e) - - ds = gb.node_props(gs) - dm = gb.node_props(gm) - - us = ds[pp.STATE][self.kw] - um = dm[pp.STATE][self.kw] - lam = d[pp.STATE][self.kw] - - bdcs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdcm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdfs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - bdfm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - - bc_val_s = ds[pp.PARAMETERS][self.kw]["bc_values"] - bc_val_m = dm[pp.PARAMETERS][self.kw]["bc_values"] - - RW = d[pp.PARAMETERS][self.kw]["robin_weight"] - MW = d[pp.PARAMETERS][self.kw]["mortar_weight"] - rhs = d[pp.PARAMETERS][self.kw]["robin_rhs"].reshape( - (gs.dim, -1), order="F" - ) - slv2mrt_nd = sps.kron(mg.slave_to_mortar_int(), sps.eye(gs.dim)).tocsr() - mstr2mrt_nd = sps.kron(mg.master_to_mortar_int(), sps.eye(gs.dim)).tocsr() - - hf2fs = pp.fvutils.map_hf_2_f(g=gs) / 2 - hf2fm = pp.fvutils.map_hf_2_f(g=gm) / 2 - jump_u = ( - slv2mrt_nd - * hf2fs - * (bdcs * us + bdfs * (bc_val_s + slv2mrt_nd.T * lam)) - - mstr2mrt_nd - * hf2fm - * (bdcm * um - bdfm * (bc_val_m + mstr2mrt_nd.T * lam)) - ).reshape((gs.dim, -1), order="F") - lam_nd = lam.reshape((gs.dim, -1), order="F") - for i in range(len(RW)): - rhs_robin = MW[i].dot(lam_nd[:, i]) + RW[i].dot(jump_u[:, i]) - self.assertTrue(np.allclose(rhs_robin, rhs[:, i])) - - def assign_discretization(self, gb, robin=True): - for _, d in gb: - d[pp.PRIMARY_VARIABLES] = {self.kw: {"cells": gb.dim_max()}} - d[pp.DISCRETIZATION] = {self.kw: {"mortar": pp.Mpsa(self.kw)}} - - if robin: - contact = pp.RobinContact(self.kw, pp.Mpsa(self.kw)) - else: - contact = pp.StressDisplacementContinuity(self.kw, pp.Mpsa(self.kw)) - for e, d in gb.edges(): - g1, g2 = gb.nodes_of_edge(e) - d[pp.PRIMARY_VARIABLES] = {self.kw: {"cells": gb.dim_max()}} - d[pp.COUPLING_DISCRETIZATION] = { - self.kw: { - g1: (self.kw, "mortar"), - g2: (self.kw, "mortar"), - e: (self.kw, contact), - } - } - - def assign_parameters(self, gb, mortar_weight, robin_weight, rhs): - for g, d in gb: - if g.grid_num == 1: - dir_faces = g.face_centers[0] < 1e-10 - elif g.grid_num == 2: - dir_faces = g.face_centers[0] > 2 - 1e-10 - elif g.grid_num == 3: - dir_faces = (g.face_centers[0] < 1e-10) + ( - g.face_centers[0] > 2 - 1e-10 - ) - - bc_val = np.zeros((g.dim, g.num_faces)) - bc_val[0, dir_faces] = 0.1 * g.face_centers[0, dir_faces] - bc = pp.BoundaryConditionVectorial(g, dir_faces, "dir") - C = pp.FourthOrderTensor( - gb.dim_max(), np.ones(g.num_cells), np.ones(g.num_cells) - ) - data = { - "bc": bc, - "bc_values": bc_val.ravel("F"), - "fourth_order_tensor": C, - "source": np.zeros(g.num_cells * g.dim), - "inverter": "python", - } - pp.initialize_data(g, d, self.kw, data) - - for _, d in gb.edges(): - mg = d["mortar_grid"] - MW = sps.diags(mortar_weight) - RW = sps.diags(robin_weight) - data = { - "mortar_weight": [MW] * mg.num_cells, - "robin_weight": [RW] * mg.num_cells, - "robin_rhs": np.tile(rhs, (mg.num_cells)), - } - pp.initialize_data(mg, d, self.kw, data) - - -class TestBiotTwoGridCoupling(unittest.TestCase): - """ - In this test we set up a coupling between two grids of dimension 2: - g_slave g_master - |-----| | |------| - | | x | | - |-----| | |------| - g_mortar - There is one cell per grid and they are coupled together by a single mortar - variable. - - We define a random Robin and Mortar weights and test if we recover the - condition on the interface. - """ - - def test_robin_coupling(self): - """ - Test a Robin condition on the interface - """ - self.kw = "mech" - self.kw_f = "flow" - gb, _ = define_gb() - mortar_weight = np.random.rand(gb.dim_max()) - robin_weight = np.random.rand(gb.dim_max()) - rhs = np.random.rand(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb) - - # Discretize - for g, d in gb: - pp.Biot( - self.kw, self.kw_f, vector_variable="u", scalar_variable="p" - ).discretize(g, d) - - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - def test_continuity_coupling(self): - """ - Test a continuity condition on the interface. This is equivalent to - zero mortar weight and identity matrix for the robin weight. These - matrices are only used to check the solution. - """ - self.kw = "mech" - self.kw_f = "flow" - gb, gb_full = define_gb() - # We assign weighs according to the condition. - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.zeros(gb.dim_max()) - # Assign data to coupling gb - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb) - # Assign data to mono gb - self.assign_parameters(gb_full, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb_full) - - # Discretize - for g, d in gb: - pp.Biot(self.kw, self.kw_f).discretize(g, d) - for g, d in gb_full: - pp.Biot(self.kw, self.kw_f).discretize(g, d) - - # Assemble and solve - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - assembler_full = pp.Assembler(gb_full) - matrix, rhs = assembler_full.assemble_matrix_rhs() - u_full = sps.linalg.spsolve(matrix, rhs) - assembler_full.distribute_variable(u_full) - # Compare solutions: - # We need to rearange the solutions because the ordering of the dofs are not the - # same for the two grids. - us = [] - ps = [] - ID = [] - for g, d in gb: - us.append(d[pp.STATE]["u"]) - ps.append(d[pp.STATE]["p"]) - ID.append(g.grid_num - 1) - us = np.hstack([np.array(us)[ID].ravel()]) - ps = np.hstack([np.array(ps)[ID].ravel()]) - IA = np.array([0, 1, 4, 5, 2, 3, 6, 7]) - IAp = np.array([0, 2, 1, 3]) - sol = np.hstack([us[IA], ps[IAp]]) - - us = [] - ps = [] - for g, d in gb_full: - us.append(d[pp.STATE]["u"]) - ps.append(d[pp.STATE]["p"]) - sol_full = np.hstack([np.array(us), np.array(ps)]) - # Note, we don't have equality because the weak symmetry breaks when a cell has to many - # Neumann conditions (see comments in mpsa) - self.assertTrue(np.all(np.abs(sol - sol_full) < 5e-4)) - - def check_solution(self, gb): - for e, d in gb.edges(): - mg = d["mortar_grid"] - gs, gm = gb.nodes_of_edge(e) - - ds = gb.node_props(gs) - dm = gb.node_props(gm) - - us = ds[pp.STATE]["u"] - um = dm[pp.STATE]["u"] - lam = d[pp.STATE]["lam_u"] - ps = ds[pp.STATE]["p"] - pm = dm[pp.STATE]["p"] - - bdcs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdcm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdfs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - bdfm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - bdps = ds[pp.DISCRETIZATION_MATRICES][self.kw][ - "bound_displacement_pressure" - ] - bdpm = dm[pp.DISCRETIZATION_MATRICES][self.kw][ - "bound_displacement_pressure" - ] - - bc_val_s = ds[pp.PARAMETERS][self.kw]["bc_values"] - bc_val_m = dm[pp.PARAMETERS][self.kw]["bc_values"] - - RW = d[pp.PARAMETERS][self.kw]["robin_weight"] - MW = d[pp.PARAMETERS][self.kw]["mortar_weight"] - rhs = d[pp.PARAMETERS][self.kw]["robin_rhs"].reshape( - (gs.dim, -1), order="F" - ) - slv2mrt_nd = sps.kron(mg.slave_to_mortar_int(), sps.eye(gs.dim)).tocsr() - mstr2mrt_nd = sps.kron(mg.master_to_mortar_int(), sps.eye(gs.dim)).tocsr() - - hf2fs = pp.fvutils.map_hf_2_f(g=gs) / 2 - hf2fm = pp.fvutils.map_hf_2_f(g=gm) / 2 - jump_u = ( - slv2mrt_nd - * hf2fs - * (bdcs * us + bdfs * (bc_val_s + slv2mrt_nd.T * lam) + bdps * ps) - - mstr2mrt_nd - * hf2fm - * (bdcm * um - bdfm * (bc_val_m + mstr2mrt_nd.T * lam) + bdpm * pm) - ).reshape((gs.dim, -1), order="F") - lam_nd = lam.reshape((gs.dim, -1), order="F") - for i in range(len(RW)): - rhs_robin = MW[i].dot(lam_nd[:, i]) + RW[i].dot(jump_u[:, i]) - self.assertTrue(np.allclose(rhs_robin, rhs[:, i])) - - def assign_discretization(self, gb): - gradP_disc = pp.GradP(self.kw) - divU_disc = pp.DivU(self.kw, variable="u", mortar_variable="lam_u") - - for g, d in gb: - d[pp.DISCRETIZATION] = { - "u": {"div_sigma": pp.Mpsa(self.kw)}, - "p": { - "flux": pp.Mpfa(self.kw_f), - "mass": pp.MassMatrix(self.kw_f), - "stab": pp.BiotStabilization(self.kw_f, variable="p"), - }, - "u_p": {"grad_p": gradP_disc}, - "p_u": {"div_u": divU_disc}, - } - - d[pp.PRIMARY_VARIABLES] = {"u": {"cells": g.dim}, "p": {"cells": 1}} - - gradP_disp = pp.numerics.interface_laws.elliptic_interface_laws.RobinContactBiotPressure( - self.kw, gradP_disc - ) - div_u_lam = pp.numerics.interface_laws.elliptic_interface_laws.DivU_StressMortar( - self.kw, divU_disc - ) - for e, d in gb.edges(): - g1, g2 = gb.nodes_of_edge(e) - d[pp.PRIMARY_VARIABLES] = {"lam_u": {"cells": 2}, "lam_p": {"cells": 1}} - d[pp.COUPLING_DISCRETIZATION] = { - "u_contribution": { - g1: ("u", "div_sigma"), - g2: ("u", "div_sigma"), - (g1, g2): ("lam_u", pp.RobinContact(self.kw, pp.Mpsa(self.kw))), - }, - "p_contribution_to_displacement": { - g1: ("p", "flux"), - g2: ("p", "flux"), - (g1, g2): ("lam_u", gradP_disp), - }, - "lam_u_contr_2_div_u": { - g1: ("p", "flux"), - g2: ("p", "flux"), - (g1, g2): ("lam_u", div_u_lam), - }, - "lam_p": { - g1: ("p", "flux"), - g2: ("p", "flux"), - (g1, g2): ( - "lam_p", - pp.FluxPressureContinuity(self.kw_f, pp.Mpfa(self.kw_f)), - ), - }, - } - - def assign_parameters(self, gb, mortar_weight, robin_weight, rhs): - for g, d in gb: - if g.grid_num == 1: - dir_faces = g.face_centers[0] < 1e-10 - elif g.grid_num == 2: - dir_faces = g.face_centers[0] > 2 - 1e-10 - elif g.grid_num == 3: - dir_faces = (g.face_centers[0] < 1e-10) + ( - g.face_centers[0] > 2 - 1e-10 - ) - bc_val = np.zeros((g.dim, g.num_faces)) - bc_val[0, dir_faces] = 0.1 * g.face_centers[0, dir_faces] - bc = pp.BoundaryConditionVectorial(g, dir_faces, "dir") - C = pp.FourthOrderTensor( - gb.dim_max(), np.ones(g.num_cells), np.ones(g.num_cells) - ) - alpha = 1 / np.pi - data = { - "bc": bc, - "bc_values": bc_val.ravel("F"), - "fourth_order_tensor": C, - "source": np.zeros(g.num_cells * g.dim), - "inverter": "python", - "biot_alpha": alpha, - } - data_f = { - "bc": pp.BoundaryCondition(g, g.get_boundary_faces(), "dir"), - "bc_values": np.zeros(g.num_faces), - "second_order_tensor": pp.SecondOrderTensor( - g.dim, np.ones(g.num_cells) - ), - "inverter": "python", - "aperture": np.ones(g.num_cells), - "biot_alpha": alpha, - "mass_weight": 1e-1, - } - pp.initialize_data(g, d, self.kw, data) - pp.initialize_data(g, d, self.kw_f, data_f) - state = { - "u": np.zeros(g.dim * g.num_cells), - "p": np.zeros(g.num_cells), - self.kw: {"bc_values": bc_val.ravel("F")}, - } - pp.set_state(d, state) - - for _, d in gb.edges(): - mg = d["mortar_grid"] - MW = sps.diags(mortar_weight) - RW = sps.diags(robin_weight) - data = { - "mortar_weight": [MW] * mg.num_cells, - "robin_weight": [RW] * mg.num_cells, - "robin_rhs": np.tile(rhs, (mg.num_cells)), - } - pp.initialize_data(mg, d, self.kw, data) - pp.set_state(d, {"lam_u": np.zeros((mg.dim + 1) * mg.num_cells)}) - - -def define_gb(): - """ - Construct grids - """ - g_s = pp.CartGrid([1, 2], [1, 2]) - g_m = pp.CartGrid([1, 2], [1, 2]) - g_full = pp.CartGrid([2, 2], [2, 2]) - g_m.nodes[0] += 1 - g_s.compute_geometry() - g_m.compute_geometry() - g_full.compute_geometry() - - g_s.grid_num = 1 - g_m.grid_num = 2 - g_full.grid_num = 3 - - gb = pp.GridBucket() - gb_full = pp.GridBucket() - gb.add_nodes([g_s, g_m]) - gb_full.add_nodes([g_full]) - - contact_s = np.where(g_s.face_centers[0] > 1 - 1e-10)[0] - contact_m = np.where(g_m.face_centers[0] < 1 + 1e-10)[0] - data = np.ones(contact_s.size, dtype=np.bool) - - shape = (g_s.num_faces, g_m.num_faces) - slave_master = sps.csc_matrix((data, (contact_m, contact_s)), shape=shape) - - mortar_grid, _, _ = pp.grids.partition.extract_subgrid(g_s, contact_s, faces=True) - - gb.add_edge([g_s, g_m], slave_master) - - gb.assign_node_ordering() - gb_full.assign_node_ordering() - - # Slave and master is defined by the node number. - # In python 3.5 the node-nombering does not follow the one given in gb.add_edge - # I guess also the face_face mapping given on the edge also should change, - # but this is not used - g_1, _ = gb.nodes_of_edge([g_s, g_m]) - if g_1.grid_num == 2: - g_m = g_s - g_s = g_1 - slave_master = slave_master.T - - mg = pp.BoundaryMortar(mortar_grid.dim, mortar_grid, slave_master.T) - gb.set_edge_prop([g_s, g_m], "mortar_grid", mg) - return gb, gb_full - - -if __name__ == "__main__": - unittest.main() diff --git a/test/unit/test_elliptic_interface_laws.py b/test/unit/test_elliptic_interface_laws.py deleted file mode 100644 index e2d3aac33d..0000000000 --- a/test/unit/test_elliptic_interface_laws.py +++ /dev/null @@ -1,346 +0,0 @@ -""" -Module for testing the vector elliptic couplings in the interface_laws. -""" - - -import numpy as np -import scipy.sparse as sps -import unittest - -import porepy as pp - - -class TestTwoGridCoupling(unittest.TestCase): - """ - In this test we set up a coupling between two grids of dimension 2: - g_slave g_master - |-----| | |------| - | | x | | - |-----| | |------| - g_mortar - There is one cell per grid and they are coupled together by a single mortar - variable. - """ - - def test_robin_assembly_master(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the master terms. - """ - self.kw = "mech" - gb = define_gb() - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.ones(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - varargs = get_variables(gb) - - robin_contact = pp.RobinContact(self.kw, MockId(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [2, 0, 2, 0, -2, 0], - [0, 2, 0, 2, 0, -2], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [-2, 0, 0, 0, 2, 0], - [0, -2, 0, 0, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 1, 1])) - ) - - def test_robin_assembly_slave(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the slave terms. - """ - self.kw = "mech" - gb = define_gb() - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.ones(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - varargs = get_variables(gb) - - robin_contact = pp.RobinContact(self.kw, MockZero(), MockId()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [2, 0, 2, 0, 2, 0], - [0, 2, 0, 2, 0, 2], - [0, 0, 2, 0, 2, 0], - [0, 0, 0, 2, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 1, 1])) - ) - - def test_robin_assembly_mortar(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the mortar terms. - """ - self.kw = "mech" - gb = define_gb() - mortar_weight = np.ones(gb.dim_max()) - robin_weight = np.zeros(gb.dim_max()) - rhs = np.ones(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - varargs = get_variables(gb) - - robin_contact = pp.RobinContact(self.kw, MockZero(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 1, 0], - [0, 0, 0, 0, 0, 1], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 1, 1])) - ) - - def test_continuity_assembly_master(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the master terms. - """ - self.kw = "mech" - gb = define_gb() - varargs = get_variables(gb) - robin_contact = pp.StressDisplacementContinuity(self.kw, MockId(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [2, 0, 2, 0, -2, 0], - [0, 2, 0, 2, 0, -2], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [-2, 0, 0, 0, 2, 0], - [0, -2, 0, 0, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 0, 0])) - ) - - def test_continuity_assembly_slave(self): - """ - We test the StressDisplacementContinuity interface law. This gives a continuity - of stress and displacement on the mortar grid. - - Test the assembly of the slave terms. - """ - self.kw = "mech" - gb = define_gb() - varargs = get_variables(gb) - - robin_contact = pp.StressDisplacementContinuity(self.kw, MockZero(), MockId()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [2, 0, 2, 0, 2, 0], - [0, 2, 0, 2, 0, 2], - [0, 0, 2, 0, 2, 0], - [0, 0, 0, 2, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 0, 0])) - ) - - def test_continuity_assembly_mortar(self): - """ - We test the StressDisplacementContinuity interface law. This gives a continuity - of stress and displacement on the mortar grid. - - Test the assembly of the mortar terms. - """ - self.kw = "mech" - gb = define_gb() - varargs = get_variables(gb) - - robin_contact = pp.StressDisplacementContinuity(self.kw, MockZero(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - ] - ) - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 0, 0])) - ) - - def assign_parameters(self, gb, mortar_weight, robin_weight, rhs): - for e, d in gb.edges(): - mg = d["mortar_grid"] - MW = sps.diags(mortar_weight) - RW = sps.diags(robin_weight) - data = { - "mortar_weight": [MW] * mg.num_cells, - "robin_weight": [RW] * mg.num_cells, - "robin_rhs": np.tile(rhs, (mg.num_cells)), - } - d[pp.PARAMETERS] = pp.Parameters(e, self.kw, data) - d[pp.DISCRETIZATION_MATRICES] = {self.kw: {}} - - -class MockId(object): - """ - returns an identity mapping - """ - - def ndof(self, g): - return g.dim * g.num_cells - - def assemble_int_bound_displacement_trace(self, *vargs): - identity_mapping(self, *vargs) - - def assemble_int_bound_stress(self, *vargs): - identity_mapping(self, *vargs) - - def enforce_neumann_int_bound(self, *varargs): - pass - - -class MockZero(object): - """ - Do nothing - """ - - def ndof(self, g): - return g.dim * g.num_cells - - def assemble_int_bound_displacement_trace(self, *vargs): - pass - - def assemble_int_bound_stress(self, *vargs): - pass - - def enforce_neumann_int_bound(self, *varargs): - pass - - -def identity_mapping(RC, g, data, data_edge, swap, cc, matrix, rhs, ind): - dof_master = g.dim * g.num_cells - dof_slave = g.dim * g.num_cells - dof_mortar = g.dim * data_edge["mortar_grid"].num_cells - - cc[ind, 0] += sps.diags(np.ones(dof_slave), shape=(dof_slave, dof_master)) - cc[ind, 1] += sps.diags(np.ones(dof_slave), shape=(dof_slave, dof_slave)) - cc[ind, 2] += sps.diags(np.ones(dof_slave), shape=(dof_slave, dof_mortar)) - - cc[2, ind] += sps.diags(np.ones(dof_mortar), shape=(dof_mortar, dof_master)) - cc[2, 2] += sps.diags(np.ones(dof_mortar), shape=(dof_mortar, dof_mortar)) - - -def get_variables(gb): - for e, data_edge in gb.edges(): - g_slave, g_master = gb.nodes_of_edge(e) - data_slave = gb.node_props(g_slave) - data_master = gb.node_props(g_master) - break - dof = np.array( - [ - g_master.dim * g_master.num_cells, - g_slave.dim * g_slave.num_cells, - g_master.dim * data_edge["mortar_grid"].num_cells, - ] - ) - matrix = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - matrix = matrix.reshape((3, 3)) - - rhs = np.empty(3, dtype=np.object) - rhs[0] = np.zeros(g_master.dim * g_master.num_cells) - rhs[1] = np.zeros(g_slave.dim * g_slave.num_cells) - rhs[2] = np.zeros(data_edge["mortar_grid"].num_cells * g_slave.dim) - - return g_master, g_slave, data_master, data_slave, data_edge, matrix - - -def define_gb(): - """ - Construct grids - """ - g1 = pp.CartGrid([1, 1]) - g2 = pp.CartGrid([1, 1]) - g1.compute_geometry() - g2.compute_geometry() - - g1.grid_num = 1 - g2.grid_num = 2 - - gb = pp.GridBucket() - gb.add_nodes([g1, g2]) - gb.add_edge([g1, g2], None) - mortar_grid = pp.Grid( - 1, - np.array([[0, 0, 0], [0, 1, 0]]).T, - sps.csc_matrix(([True, True], [0, 1], [0, 1, 2])), - sps.csc_matrix(([1, -1], [0, 1], [0, 2])), - "mortar_grid", - ) - mortar_grid.compute_geometry() - face_faces = sps.csc_matrix(([True], [0], [0, 0, 1, 1, 1]), shape=(4, 4)) - mg = pp.BoundaryMortar(1, mortar_grid, face_faces) - mg.num_cells = 1 - gb.set_edge_prop([g1, g2], "mortar_grid", mg) - return gb - - -if __name__ == "__main__": - unittest.main()