Coulomb.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <Eigen/Core>
#include <utility>
 
#include "FractureModelBase.h"
#include "NumLib/NewtonRaphson.h"
#include "ParameterLib/Parameter.h"
 
namespace MaterialLib
{
namespace Fracture
{
namespace Coulomb
{
template <int DisplacementDim>
struct StateVariables
    : public FractureModelBase<DisplacementDim>::MaterialStateVariables
{
    void setInitialConditions() { w_p = w_p_prev; }
 
    void pushBackState() override { w_p_prev = w_p; }

    Eigen::Matrix<double, DisplacementDim, 1> w_p =
        Eigen::Matrix<double, DisplacementDim, 1>::Zero();
 
    // Initial values from previous timestep
    Eigen::Matrix<double, DisplacementDim, 1> w_p_prev;  
 
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
 
template <int DisplacementDim>
class Coulomb final : public FractureModelBase<DisplacementDim>
{
public:
    std::unique_ptr<
        typename FractureModelBase<DisplacementDim>::MaterialStateVariables>
    createMaterialStateVariables() override
    {
        return std::make_unique<StateVariables<DisplacementDim>>();
    }
 
public:
    struct MaterialProperties
    {
        using P = ParameterLib::Parameter<double>;
        using X = ParameterLib::SpatialPosition;
 
        MaterialProperties(
                P const& normal_stiffness_, P const& shear_stiffness_,
                P const& friction_angle_, P const& dilatancy_angle_,
                P const& cohesion_)
            : normal_stiffness(normal_stiffness_), shear_stiffness(shear_stiffness_),
              friction_angle(friction_angle_), dilatancy_angle(dilatancy_angle_),
              cohesion(cohesion_)
        {
        }

        P const& normal_stiffness;
        P const& shear_stiffness;
        P const& friction_angle;
        P const& dilatancy_angle;
        P const& cohesion;
    };
 
 
public:
    explicit Coulomb(
        NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
        double const penalty_aperture_cutoff,
        bool const tension_cutoff,
        MaterialProperties material_properties)
        : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
          _penalty_aperture_cutoff(penalty_aperture_cutoff),
          _tension_cutoff(tension_cutoff),
          _mp(std::move(material_properties))
    {
    }

    void computeConstitutiveRelation(
        double const t,
        ParameterLib::SpatialPosition const& x,
        double const aperture0,
        Eigen::Ref<Eigen::VectorXd const>
            sigma0,
        Eigen::Ref<Eigen::VectorXd const>
            w_prev,
        Eigen::Ref<Eigen::VectorXd const>
            w,
        Eigen::Ref<Eigen::VectorXd const>
            sigma_prev,
        Eigen::Ref<Eigen::VectorXd>
            sigma,
        Eigen::Ref<Eigen::MatrixXd>
            Kep,
        typename FractureModelBase<DisplacementDim>::MaterialStateVariables&
            material_state_variables) override;
 
private:
    NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;

    double const _penalty_aperture_cutoff;

    bool const _tension_cutoff;
 
    MaterialProperties _mp;
};
 
}  // namespace Coulomb
}  // namespace Fracture
}  // namespace MaterialLib
 
namespace MaterialLib
{
namespace Fracture
{
namespace Coulomb
{
extern template class Coulomb<2>;
extern template class Coulomb<3>;
}  // namespace Coulomb
}  // namespace Fracture
}  // namespace MaterialLib