text: /**
   * \file
   *
   * \copyright
   * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
   *            Distributed under a Modified BSD License.
   *              See accompanying file LICENSE.txt or
   *              http://www.opengeosys.org/project/license
   */
  
  #pragma once
  
  #include <spdlog/spdlog.h>
  
  #include <Eigen/Eigenvalues>
  
  #include "MaterialLib/SolidModels/Ehlers.h"
  
  namespace ProcessLib
  {
  namespace SmallDeformationNonlocal
  {
  /// Computes the damage internal material variable explicitly based on the
  /// results obtained from the local stress return algorithm.
  inline double calculateDamage(double const kappa_d, double const alpha_d,
                                double const beta_d)
  {
      double const damage = (1 - beta_d) * (1 - std::exp(-kappa_d / alpha_d));
  
      if (damage < 0. || damage >= 1.)
      {
          OGS_FATAL("Damage value {:g} outside of [0,1) interval.", damage);
      }
  
      return damage;
  }
  
  template <int DisplacementDim, typename KelvinVectorType>
  double calculateDamageKappaD(
      double const eps_p_eff_diff,
      KelvinVectorType const& sigma,
      double const kappa_d_prev,
      double const h_d,
      MaterialLib::Solids::Ehlers::MaterialProperties const& mp)
  {
      // Default case of the rate problem. Updated below if volumetric plastic
      // strain rate is positive (dilatancy).
  
      // non-const for Eigen solver.
      auto sigma_tensor = MathLib::KelvinVector::kelvinVectorToTensor(sigma);
  
      Eigen::EigenSolver<decltype(sigma_tensor)> eigen_solver(sigma_tensor);
      auto const principal_stress = real(eigen_solver.eigenvalues().array());
      double const prod_stress = std::sqrt(principal_stress.square().sum());
  
      // Brittleness decrease with confinement for the nonlinear flow rule.
      // ATTENTION: For linear flow rule -> constant brittleness.
      double const tensile_strength =
          std::sqrt(3.0) * mp.kappa / (1 + std::sqrt(3.0) * mp.beta);
      double const r_s = prod_stress / tensile_strength - 1.;
  
      // Compute normalizing strain.
      double const x_s = [](double const h_d, double const r_s) {
          if (r_s < 0)
          {
              return 1.;
          }
          if (r_s <= 1)
          {
              return 1. + h_d * r_s * r_s;
          }
          return 1. - 3 * h_d + 4 * h_d * std::sqrt(r_s);
      }(h_d, r_s);
  
      return kappa_d_prev + eps_p_eff_diff / x_s;
  }
  }  // namespace SmallDeformationNonlocal
  }  // namespace ProcessLib
  