"""
Class types:
- MomentumBalanceEquations defines subdomain and interface equations through the
terms entering. Momentum balance between opposing fracture interfaces is imposed.
Notes:
- The class MomentumBalanceEquations is a mixin class, and should be inherited by a
class that defines the variables and discretization.
- Refactoring needed for constitutive equations. Modularisation and moving to the
library.
"""
from __future__ import annotations
import logging
from functools import partial
from typing import Callable, Optional, Sequence
import numpy as np
import porepy as pp
from . import constitutive_laws
logger = logging.getLogger(__name__)
[docs]
class MomentumBalanceEquations(pp.BalanceEquation):
"""Class for momentum balance equations and fracture deformation equations."""
stress: Callable[[list[pp.Grid]], pp.ad.Operator]
"""Stress on the grid faces. Provided by a suitable mixin class that specifies the
physical laws governing the stress.
"""
mdg: pp.MixedDimensionalGrid
"""Mixed dimensional grid for the current model. Normally defined in a mixin
instance of :class:`~porepy.models.geometry.ModelGeometry`.
"""
interfaces_to_subdomains: Callable[[list[pp.MortarGrid]], list[pp.Grid]]
"""Map from interfaces to the adjacent subdomains. Normally defined in a mixin
instance of :class:`~porepy.models.geometry.ModelGeometry`.
"""
internal_boundary_normal_to_outwards: Callable[[list[pp.Grid], int], pp.ad.Matrix]
"""Switch interface normal vectors to point outwards from the subdomain. Normally
set by a mixin instance of :class:`porepy.models.geometry.ModelGeometry`.
"""
fracture_stress: Callable[[list[pp.MortarGrid]], pp.ad.Operator]
"""Stress on the fracture faces. Provided by a suitable mixin class that specifies
the physical laws governing the stress, see for instance
:class:`~porepy.models.constitutive_laws.LinearElasticMechanicalStress` or
:class:`~porepy.models.constitutive_laws.PressureStress`.
"""
basis: Callable[[Sequence[pp.GridLike], int], list[pp.ad.Matrix]]
"""Basis for the local coordinate system. Normally set by a mixin instance of
:class:`porepy.models.geometry.ModelGeometry`.
"""
normal_component: Callable[[list[pp.Grid]], pp.ad.Matrix]
"""Operator giving the normal component of vectors. Normally defined in a mixin
instance of :class:`~porepy.models.models.ModelGeometry`.
"""
tangential_component: Callable[[list[pp.Grid]], pp.ad.Matrix]
"""Operator giving the tangential component of vectors. Normally defined in a mixin
instance of :class:`~porepy.models.models.ModelGeometry`.
"""
displacement_jump: Callable[[list[pp.Grid]], pp.ad.Operator]
"""Operator giving the displacement jump on fracture grids. Normally defined in a
mixin instance of :class:`~porepy.models.models.ModelGeometry`.
"""
contact_traction: Callable[[list[pp.Grid]], pp.ad.MixedDimensionalVariable]
"""Contact traction variable. Normally defined in a mixin instance of
:class:`~porepy.models.momentum_balance.VariablesMomentumBalance`.
"""
gap: Callable[[list[pp.Grid]], pp.ad.Operator]
"""Gap of a fracture. Normally provided by a mixin instance of
:class:`~porepy.models.constitutive_laws.FracturedSolid`.
"""
friction_bound: Callable[[list[pp.Grid]], pp.ad.Operator]
"""Friction bound of a fracture. Normally provided by a mixin instance of
:class:`~porepy.models.constitutive_laws.FrictionBound`.
"""
contact_mechanics_numerical_constant: Callable[[list[pp.Grid]], pp.ad.Scalar]
"""Numerical constant for contact mechanics. Normally provided by a mixin instance
of :class:`~porepy.models.momuntum_balance.SolutionStrategyMomentumBalance`.
"""
[docs]
def set_equations(self):
"""Set equations for the subdomains and interfaces.
The following equations are set:
- Momentum balance in the matrix.
- Force balance between fracture interfaces.
- Deformation constraints for fractures, split into normal and tangential
part.
See individual equation methods for details.
"""
matrix_subdomains = self.mdg.subdomains(dim=self.nd)
fracture_subdomains = self.mdg.subdomains(dim=self.nd - 1)
interfaces = self.mdg.interfaces(dim=self.nd - 1)
matrix_eq = self.momentum_balance_equation(matrix_subdomains)
# We split the fracture deformation equations into two parts, for the normal and
# tangential components for convenience.
fracture_eq_normal = self.normal_fracture_deformation_equation(
fracture_subdomains
)
fracture_eq_tangential = self.tangential_fracture_deformation_equation(
fracture_subdomains
)
intf_eq = self.interface_force_balance_equation(interfaces)
self.equation_system.set_equation(
matrix_eq, matrix_subdomains, {"cells": self.nd}
)
self.equation_system.set_equation(
fracture_eq_normal, fracture_subdomains, {"cells": 1}
)
self.equation_system.set_equation(
fracture_eq_tangential, fracture_subdomains, {"cells": self.nd - 1}
)
self.equation_system.set_equation(intf_eq, interfaces, {"cells": self.nd})
[docs]
def momentum_balance_equation(self, subdomains: list[pp.Grid]):
"""Momentum balance equation in the matrix.
Inertial term is not included.
Parameters:
subdomains: List of subdomains where the force balance is defined. Only
known usage
is for the matrix domain(s).
Returns:
Operator for the force balance equation in the matrix.
"""
accumulation = self.inertia(subdomains)
stress = self.stress(subdomains)
body_force = self.body_force(subdomains)
equation = self.balance_equation(
subdomains, accumulation, stress, body_force, dim=self.nd
)
equation.set_name("momentum_balance_equation")
return equation
[docs]
def inertia(self, subdomains: list[pp.Grid]):
"""Inertial term [m^2/s].
Added here for completeness, but not used in the current implementation. Be
aware that the elasticity discretization has employed herein has, as far as we
know, never been used to solve a problem with inertia. Thus, if inertia is added
in a submodel, proceed with caution. In addition to overriding this method, it
would also be necessary to add the inertial term to the balance equation
:meth:`momentum_balance_equation`.
Parameters:
subdomains: List of subdomains where the inertial term is defined.
Returns:
Operator for the inertial term.
"""
return pp.ad.Scalar(0)
[docs]
def interface_force_balance_equation(
self,
interfaces: list[pp.MortarGrid],
) -> pp.ad.Operator:
"""Momentum balance equation at matrix-fracture interfaces.
Parameters:
interfaces: Fracture-matrix interfaces.
Returns:
Operator representing the force balance equation.
Raises:
ValueError: If an interface is not a fracture-matrix interface.
"""
# Check that the interface is a fracture-matrix interface.
for interface in interfaces:
if interface.dim != self.nd - 1:
raise ValueError("Interface must be a fracture-matrix interface.")
subdomains = self.interfaces_to_subdomains(interfaces)
# Split into matrix and fractures. Sort on dimension to allow for multiple
# matrix domains. Otherwise, we could have picked the first element.
matrix_subdomains = [sd for sd in subdomains if sd.dim == self.nd]
# Geometry related
mortar_projection = pp.ad.MortarProjections(
self.mdg, subdomains, interfaces, self.nd
)
proj = pp.ad.SubdomainProjections(subdomains, self.nd)
# Contact traction from primary grid and mortar displacements (via primary grid).
# Spelled out for clarity:
# 1) The sign of the stress on the matrix subdomain is corrected so that all
# stress components point outwards from the matrix (or inwards, EK is not
# completely sure, but the point is the consistency).
# 2) The stress is prolonged from the matrix subdomains to all subdomains seen
# by the mortar grid (that is, the matrix and the fracture).
# 3) The stress is projected to the mortar grid.
contact_from_primary_mortar = (
mortar_projection.primary_to_mortar_int
* proj.face_prolongation(matrix_subdomains)
* self.internal_boundary_normal_to_outwards(
matrix_subdomains, dim=self.nd # type: ignore[call-arg]
)
* self.stress(matrix_subdomains)
)
# Traction from the actual contact force.
traction_from_secondary = self.fracture_stress(interfaces)
# The force balance equation. Note that the force from the fracture is a
# traction, not a stress, and must be scaled with the area of the interface.
# This is not the case for the force from the matrix, which is a stress.
force_balance_eq: pp.ad.Operator = (
contact_from_primary_mortar
+ self.volume_integral(traction_from_secondary, interfaces, dim=self.nd)
)
force_balance_eq.set_name("interface_force_balance_equation")
return force_balance_eq
[docs]
def body_force(self, subdomains: list[pp.Grid]):
"""Body force integrated over the subdomain cells.
FIXME: See FluidMassBalanceEquations.fluid_source.
Parameters:
subdomains: List of subdomains where the body force is defined.
Returns:
Operator for the body force.
"""
num_cells = sum([sd.num_cells for sd in subdomains])
vals = np.zeros(num_cells * self.nd)
source = pp.ad.Array(vals, "body_force")
return source
[docs]
class ConstitutiveLawsMomentumBalance(
constitutive_laws.LinearElasticSolid,
constitutive_laws.FracturedSolid,
constitutive_laws.FrictionBound,
):
"""Class for constitutive equations for momentum balance equations."""
[docs]
def stress(self, subdomains: list[pp.Grid]):
"""Stress operator.
Parameters:
subdomains: List of subdomains where the stress is defined.
Returns:
Operator for the stress.
"""
# Method from constitutive library's LinearElasticRock.
return self.mechanical_stress(subdomains)
[docs]
class VariablesMomentumBalance:
"""
Variables for mixed-dimensional deformation:
Displacement in matrix and on fracture-matrix interfaces. Fracture contact
traction.
.. note::
Implementation postponed till Veljko's more convenient SystemManager is available.
"""
displacement_variable: str
"""Name of the primary variable representing the displacement in subdomains.
Normally defined in a mixin of instance
:class:`~porepy.models.momentum_balance.SolutionStrategyMomentumBalance`.
"""
interface_displacement_variable: str
"""Name of the primary variable representing the displacement on an interface.
Normally defined in a mixin of instance
:class:`~porepy.models.momentum_balance.SolutionStrategyMomentumBalance`.
"""
contact_traction_variable: str
"""Name of the primary variable representing the contact traction on a fracture
subdomain. Normally defined in a mixin of instance
:class:`~porepy.models.momentum_balance.SolutionStrategyMomentumBalance`.
"""
mdg: pp.MixedDimensionalGrid
"""Mixed dimensional grid for the current model. Normally defined in a mixin
instance of :class:`~porepy.models.geometry.ModelGeometry`.
"""
nd: int
"""Ambient dimension of the problem. Normally set by a mixin instance of
:class:`porepy.models.geometry.ModelGeometry`.
"""
subdomains_to_interfaces: Callable[[list[pp.Grid], list[int]], list[pp.MortarGrid]]
"""Map from subdomains to the adjacent interfaces. Normally defined in a mixin
instance of :class:`~porepy.models.geometry.ModelGeometry`.
"""
equation_system: pp.ad.EquationSystem
"""EquationSystem object for the current model. Normally defined in a mixin class
defining the solution strategy.
"""
local_coordinates: Callable[[list[pp.Grid]], pp.ad.Matrix]
"""Mapping to local coordinates. Normally defined in a mixin instance of
:class:`~porepy.models.geometry.ModelGeometry`.
"""
[docs]
def create_variables(self):
"""Set variables for the subdomains and interfaces.
The following variables are set:
- Displacement in the matrix.
- Displacement on fracture-matrix interfaces.
- Fracture contact traction.
See individual variable methods for details.
"""
self.equation_system.create_variables(
dof_info={"cells": self.nd},
name=self.displacement_variable,
subdomains=self.mdg.subdomains(dim=self.nd),
)
self.equation_system.create_variables(
dof_info={"cells": self.nd},
name=self.interface_displacement_variable,
interfaces=self.mdg.interfaces(dim=self.nd - 1),
)
self.equation_system.create_variables(
dof_info={"cells": self.nd},
name=self.contact_traction_variable,
subdomains=self.mdg.subdomains(dim=self.nd - 1),
)
[docs]
def displacement(self, subdomains: list[pp.Grid]):
"""Displacement in the matrix.
Parameters:
grids: List of subdomains or interface grids where the displacement is
defined. Should be the matrix subdomains.
Returns:
Variable for the displacement.
Raises:
ValueError: If the dimension of the subdomains is not equal to the ambient
dimension of the problem.
"""
if not all([sd.dim == self.nd for sd in subdomains]):
raise ValueError(
"Displacement is only defined in subdomains of dimension nd."
)
return self.equation_system.md_variable(self.displacement_variable, subdomains)
[docs]
def interface_displacement(self, interfaces: list[pp.MortarGrid]):
"""Displacement on fracture-matrix interfaces.
Parameters:
interfaces: List of interface grids where the displacement is defined.
Should be between the matrix and fractures.
Returns:
Variable for the displacement.
Raises:
ValueError: If the dimension of the interfaces is not equal to the ambient
dimension of the problem minus one.
"""
if not all([intf.dim == self.nd - 1 for intf in interfaces]):
raise ValueError(
"Interface displacement is only defined on interfaces of dimension "
"nd - 1."
)
return self.equation_system.md_variable(
self.interface_displacement_variable, interfaces
)
[docs]
def displacement_jump(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
"""Displacement jump on fracture-matrix interfaces.
Parameters:
subdomains: List of subdomains where the displacement jump is defined.
Should be a fracture subdomain.
Returns:
Operator for the displacement jump.
Raises:
AssertionError: If the subdomains are not fractures, i.e. have dimension
`nd - 1`.
"""
if not all([sd.dim == self.nd - 1 for sd in subdomains]):
raise ValueError("Displacement jump only defined on fractures")
interfaces = self.subdomains_to_interfaces(subdomains, [1])
# Only use matrix-fracture interfaces
interfaces = [intf for intf in interfaces if intf.dim == self.nd - 1]
mortar_projection = pp.ad.MortarProjections(
self.mdg, subdomains, interfaces, self.nd
)
# The displacement jmup is expressed in the local coordinates of the fracture.
# First use the sign of the mortar sides to get a difference, then map first
# from the interface to the fracture, and finally to the local coordinates.
rotated_jumps: pp.ad.Operator = (
self.local_coordinates(subdomains)
* mortar_projection.mortar_to_secondary_avg
* mortar_projection.sign_of_mortar_sides
* self.interface_displacement(interfaces)
)
rotated_jumps.set_name("Rotated_displacement_jump")
return rotated_jumps
[docs]
class SolutionStrategyMomentumBalance(pp.SolutionStrategy):
"""This is whatever is left of pp.ContactMechanics.
At some point, this will be refined to be a more sophisticated (modularised)
solution strategy class. More refactoring may be beneficial.
Parameters:
params: Parameters for the solution strategy.
"""
nd: int
"""Ambient dimension of the problem. Normally set by a mixin instance of
:class:`porepy.models.geometry.ModelGeometry`.
"""
solid: pp.SolidConstants
"""Solid constant object that takes care of scaling of solid-related quantities.
Normally, this is set by a mixin of instance
:class:`~porepy.models.solution_strategy.SolutionStrategy`.
"""
equation_system: pp.ad.EquationSystem
"""EquationSystem object for the current model. Normally defined in a mixin class
defining the solution strategy.
"""
stiffness_tensor: Callable[[pp.Grid], pp.FourthOrderTensor]
"""Function that returns the stiffness tensor of a subdomain. Normally provided by a
mixin of instance :class:`~porepy.models.constitutive_laws.LinearElasticSolid`.
"""
bc_type_mechanics: Callable[[pp.Grid], pp.BoundaryConditionVectorial]
"""Function that returns the boundary condition type for the momentum problem.
Normally provided by a mixin instance of
:class:`~porepy.models.momentum_balance.BoundaryConditionsMomentumBalance`.
"""
def __init__(self, params: Optional[dict] = None) -> None:
super().__init__(params)
# Variables
self.displacement_variable: str = "u"
"""Name of the displacement variable."""
self.interface_displacement_variable: str = "u_interface"
"""Name of the displacement variable on fracture-matrix interfaces."""
self.contact_traction_variable: str = "t"
"""Name of the contact traction variable."""
# Discretization
self.stress_keyword: str = "mechanics"
"""Keyword for stress term.
Used to access discretization parameters and store discretization matrices.
"""
[docs]
def initial_condition(self) -> None:
"""Set initial guess for the variables.
The displacement is set to zero in the Nd-domain, and at the fracture interfaces
The displacement jump is thereby also zero.
The contact pressure is set to zero in the tangential direction, and -1 (that
is, in contact) in the normal direction.
"""
# Zero for displacement and initial bc values for Biot
super().initial_condition()
# Contact as initial guess. Ensure traction is consistent with zero jump, which
# follows from the default zeros set for all variables, specifically interface
# displacement, by super method.
num_frac_cells = sum(
sd.num_cells for sd in self.mdg.subdomains(dim=self.nd - 1)
)
traction_vals = np.zeros((self.nd, num_frac_cells))
traction_vals[-1] = self.solid.convert_units(-1, "Pa")
self.equation_system.set_variable_values(
traction_vals.ravel("F"),
[self.contact_traction_variable],
to_state=True,
to_iterate=True,
)
[docs]
def set_discretization_parameters(self) -> None:
"""Set discretization parameters for the simulation."""
super().set_discretization_parameters()
for sd, data in self.mdg.subdomains(return_data=True):
if sd.dim == self.nd:
pp.initialize_data(
sd,
data,
self.stress_keyword,
{
"bc": self.bc_type_mechanics(sd),
"fourth_order_tensor": self.stiffness_tensor(sd),
},
)
def _is_nonlinear_problem(self) -> bool:
"""
If there is no fracture, the problem is usually linear. Overwrite this function
if e.g. parameter nonlinearities are included.
"""
return self.mdg.dim_min() < self.nd
[docs]
class BoundaryConditionsMomentumBalance:
"""Boundary conditions for the momentum balance."""
nd: int
"""Ambient dimension of the problem. Normally set by a mixin instance of
:class:`porepy.models.geometry.ModelGeometry`.
"""
domain_boundary_sides: Callable[[pp.Grid], pp.bounding_box.DomainSides]
[docs]
def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
"""Define type of boundary conditions.
Parameters:
sd: Subdomain grid.
Returns:
bc: Boundary condition representation. Dirichlet on all global boundaries,
Dirichlet also on fracture faces.
"""
# Define boundary faces.
boundary_faces = self.domain_boundary_sides(sd).all_bf
bc = pp.BoundaryConditionVectorial(sd, boundary_faces, "dir")
# Default internal BC is Neumann. We change to Dirichlet for the contact
# problem. I.e., the mortar variable represents the displacement on the fracture
# faces.
bc.internal_to_dirichlet(sd)
return bc
[docs]
def bc_values_mechanics(self, subdomains: list[pp.Grid]) -> pp.ad.Array:
"""Boundary values for the momentum balance.
Parameters:
subdomains: List of subdomains.
Returns:
bc_values: Array of boundary condition values, zero by default. If combined
with transient problems in e.g. Biot, this should be a
:class:`pp.ad.TimeDependentArray` (or a variable following BoundaryGrid
extension).
"""
num_faces = sum([sd.num_faces for sd in subdomains])
return pp.wrap_as_ad_array(0, num_faces * self.nd, "bc_vals_mechanics")
# Note that we ignore a mypy error here. There are some inconsistencies in the method
# definitions of the mixins, related to the enforcement of keyword-only arguments. The
# type Callable is poorly supported, except if protocols are used and we really do not
# want to go there. Specifically, method definitions that contains a *, for instance,
# def method(a: int, *, b: int) -> None: pass
# which should be types as Callable[[int, int], None], cannot be parsed by mypy.
# For this reason, we ignore the error here, and rely on the tests to catch any
# inconsistencies.
[docs]
class MomentumBalance( # type: ignore[misc]
MomentumBalanceEquations,
VariablesMomentumBalance,
ConstitutiveLawsMomentumBalance,
BoundaryConditionsMomentumBalance,
SolutionStrategyMomentumBalance,
pp.ModelGeometry,
pp.DataSavingMixin,
):
"""Class for mixed-dimensional momentum balance with contact mechanics."""
pass