material.h

Go to the documentation of this file.

/*
 * -------------------------------------------
 * Copyright (c) 2021 - 2026 Prashant K. Jha
 * -------------------------------------------
 * PeriDEM https://github.com/prashjha/PeriDEM
 *
 * Distributed under the Boost Software License, Version 1.0. (See accompanying
 * file LICENSE)
 */
 
#ifndef MATERIAL_PARTILCE_MATERIAL_H
#define MATERIAL_PARTILCE_MATERIAL_H
 
#include "influenceFn.h"
#include "util/io.h"
#include <stdexcept>
#include "inp/materialDeck.h"
#include "util/function.h"
#include "util/point.h" // definition of Point
#include <limits>
#include <memory>
#include <string>
#include <vector>
 
namespace {
 
size_t dimension = 0;
 
bool is_plane_strain = false;
 
std::shared_ptr<material::BaseInfluenceFn> influence_fn;
 
double getGlobalInfFn(const double &r) {
    return influence_fn->getInfFn(r);
}
 
double getGlobalMoment(const size_t &i) {return influence_fn->getMoment(i);}
 
double getModelBulkModulus(const size_t &dim, const bool &plane_strain,
                           const double &K, const double &G) {
  if (dim != 2)
    return K;
  return plane_strain ? K + G / 3. : 9. * K * G / (3. * K + 4. * G);
}
}
 
namespace material {
 
class Material {
 
public:
  explicit Material(std::string name = "") : d_name(name) {}
 
  virtual ~Material() {}
 
  std::string name() { return d_name; }
 
  size_t getDimension() const { return dimension; }
 
  bool isPlaneStrain() const { return is_plane_strain; }
 
  virtual bool isStateActive() const = 0;
 
  virtual std::pair<double, double>
  getBondEF(const double &r, const double &s, bool &fs,
            const bool &break_bonds) const = 0;
 
  virtual std::pair<double, double>
  getBondEF(const double &r, const double &s, bool &fs, const double
  &mx, const double &thetax) const = 0;
 
  virtual util::Point getBondForceDirection(const util::Point &dx,
                                             const util::Point &du) const = 0;
 
  virtual double getS(const util::Point &dx, const util::Point &du) const = 0;
 
  virtual double getSc(const double &r) const = 0;
 
  virtual double getBreakSc(const double &r) const { return getSc(r); }
 
  virtual double getDilatationEnergyDensity(const double &thetax) const {
    return 0.;
  }
 
  virtual double getDensity() const = 0;
 
  virtual double getInfFn(const double &r) const = 0;
 
  virtual double getMoment(const size_t &i) const = 0;
 
  virtual double getHorizon() const = 0;
 
  virtual inp::MatData computeMaterialProperties(const size_t &dim) const = 0;
 
  virtual std::string printStr(int nt, int lvl) const {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- particle::Material --------" << std::endl
        << std::endl;
    oss << tabS << "Abstract class of peridynamic materials" << std::endl;
    oss << tabS << "See RnpMaterial and PmbMaterial for implementation"
        << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  virtual void print(int nt, int lvl) const { std::cout << printStr(nt, lvl); }
 
  virtual void print() const { print(0, 0); }
 
private:
  std::string d_name;
};
 
class RnpMaterial : public Material {
 
public:
  RnpMaterial(inp::MaterialDeck &deck, const size_t &dim, const double &horizon)
      : Material("RNPBond"), d_horizon(horizon), d_density(deck.d_density),
        d_C(0.), d_beta(0.), d_rbar(0.), d_invFactor(0.),
        d_factorSc(deck.d_checkScFactor),
        d_irrevBondBreak(deck.d_irreversibleBondBreak) {
 
    // set global fields
    if (dimension != dim)
      dimension = dim;
 
    if (is_plane_strain != deck.d_isPlaneStrain)
      is_plane_strain = deck.d_isPlaneStrain;
 
    // create influence function
    if (deck.d_influenceFnType == 0) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::ConstInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 1) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::LinearInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 2) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::GaussianInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Influence function type = "
          << deck.d_influenceFnType
          << " is invalid.\n");
    }
 
    if (dim == 1)
      d_invFactor = std::pow(horizon, 2) * 2.;
    else if (dim == 2)
      d_invFactor = std::pow(horizon, 3) * M_PI;
    else if (dim == 3)
      d_invFactor = std::pow(horizon, 4) * 4. * M_PI / 3.;
 
    // check if we need to compute the material parameters
    if (deck.d_computeParamsFromElastic)
      computeParameters(deck, dim);
    else {
      d_C = deck.d_bondPotentialParams[0];
      d_beta = deck.d_bondPotentialParams[1];
      d_rbar = std::sqrt(0.5 / d_beta);
    }
  };
 
  bool isStateActive() const override { return false; };
 
  std::pair<double, double> getBondEF(const double &r, const double &s,
                                      bool &fs,
                                      const bool &break_bonds) const override {
 
    if (break_bonds && !fs && util::isGreater(std::abs(s), getBreakSc(r)))
      fs = true;
 
    // intact bonds follow the nonlinear potential; broken bonds keep the
    // saturated energy and carry no force
    if (!fs)
      return std::make_pair(
          getInfFn(r) * d_C *
              (1. - std::exp(-d_beta * r * s * s)) / d_invFactor,
          getInfFn(r) * 4. * s * d_C * d_beta *
              std::exp(-d_beta * r * s * s) / d_invFactor);
    else
      return std::make_pair(getInfFn(r) * d_C / d_invFactor, 0.);
  };
 
  std::pair<double, double>
  getBondEF(const double &r, const double &s, bool &fs, const double
  &mx, const double &thetax) const override {
 
    return this->getBondEF(r, s, fs, true);
  };
 
  util::Point getBondForceDirection(const util::Point &dx,
                                     const util::Point &du) const override {
    return dx / dx.length();
  };
 
  double getS(const util::Point &dx, const util::Point &du) const override {
    return dx.dot(du) / dx.dot(dx);
  };
 
  double getSc(const double &r) const override {
    return d_rbar / std::sqrt(r);
  };
 
  double getBreakSc(const double &r) const override {
    if (!d_irrevBondBreak)
      return std::numeric_limits<double>::max();
    return d_factorSc * getSc(r);
  };
 
  double getDensity() const override { return d_density; };
 
  double getInfFn(const double &r) const override {
    return getGlobalInfFn(r / d_horizon);
  };
 
  double getMoment(const size_t &i) const override {
    return getGlobalMoment(i);
  };
 
  double getHorizon() const override { return d_horizon; };
 
  inp::MatData computeMaterialProperties(const size_t &dim) const override {
 
    auto data = inp::MatData();
 
    // set Poisson's ratio to 1/4
    data.d_nu = 0.25;
 
    // get moment of influence function
    double M = getMoment(dim);
 
    // inverse of computeParameters
    if (dim == 2) {
      data.d_Gc = 4. * M * d_C / M_PI;
      data.d_lambda = d_C * M * d_beta / 2.;
    } else if (dim == 3) {
      data.d_Gc = 3. * M * d_C / 2.;
      data.d_lambda = 2. * d_C * M * d_beta / 5.;
    }
    data.d_mu = data.d_lambda;
    data.d_G = data.d_lambda;
    data.d_E = data.toELambda(data.d_lambda,
                              data.d_nu);
    data.d_K =
        data.toK(data.d_E, data.d_nu);
    data.d_KIc = data.toKIc(
        data.d_Gc, data.d_nu, data.d_E);
 
    return data;
  };
 
  std::string printStr(int nt, int lvl) const override {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- particle::RnpMaterial --------" << std::endl
        << std::endl;
    oss << tabS << "State active = " << 0 << std::endl;
    oss << tabS << "Horizon = " << d_horizon << std::endl;
    oss << tabS << "Influence fn address = " << influence_fn.get() << std::endl;
    oss << tabS << "Influence fn info: " << std::endl;
    oss << influence_fn->printStr(nt + 1, lvl);
    oss << tabS << "Peridynamic parameters: " << std::endl;
    oss << tabS << "  C = " << d_C << std::endl;
    oss << tabS << "  beta = " << d_beta << std::endl;
    oss << tabS << "  r_bar = " << d_rbar << std::endl;
    oss << tabS << "  inv_factor = " << d_invFactor << std::endl;
    oss << tabS << "  factor_Sc = " << d_factorSc << std::endl;
    oss << tabS << "  irrev_bond_breaking = " << d_irrevBondBreak << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  void print(int nt, int lvl) const override {
    std::cout << printStr(nt, lvl);
  }
 
  void print() const override { print(0, 0); }
 
private:
  void computeParameters(inp::MaterialDeck &deck, const size_t &dim) {
    //
    // Need following elastic and fracture properties
    // 1. E or K
    // 2. Gc or KIc
    // For bond-based, Poisson's ratio is fixed to 1/4, so 2D must be plane
    // strain (Lipton 2016; Jha 2025, sec. 7.2).
    //
    if (dim == 2 && !is_plane_strain) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: RNP calibration needs nu = 1/4, which in 2D means plane "
             "strain. Set Is_Plane_Strain = true.\n");
    }
    if (util::isLess(deck.d_matData.d_E, 0.) &&
        util::isLess(deck.d_matData.d_K, 0.)) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Require either Young's modulus E or Bulk modulus K"
          " to compute the RNP bond-based peridynamic parameters.\n");
    }
    if (util::isGreater(deck.d_matData.d_E, 0.) &&
        util::isGreater(deck.d_matData.d_K, 0.)) {
      std::cout << "Warning: Both Young's modulus E and Bulk modulus K are "
                   "provided.\n";
      std::cout << "Warning: To compute the RNP bond-based peridynamic "
                   "parameters, we only require one of those.\n";
      std::cout
          << "Warning: Selecting Young's modulus to compute parameters.\n";
    }
 
    if (util::isLess(deck.d_matData.d_Gc, 0.) &&
        util::isLess(deck.d_matData.d_KIc, 0.)) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Require either critical energy release rate Gc or "
          "critical stress intensity factor KIc to compute the RNP "
          "bond-based peridynamic parameters.\n");
    } else if (util::isGreater(deck.d_matData.d_Gc, 0.) &&
               util::isGreater(deck.d_matData.d_KIc, 0.)) {
      std::cout << "Warning: Both critical energy release rate Gc and critical "
                   "stress intensity factor KIc are provided.\n";
      std::cout << "Warning: To compute the RNP bond-based peridynamic "
                   "parameters, we only require one of those.\n";
      std::cout << "Warning: Selecting critical energy release rate Gc to "
                   "compute parameters.\n";
    }
 
    // set Poisson's ratio to 1/4
    if (util::isGreater(deck.d_matData.d_nu, 0.) &&
        std::abs(deck.d_matData.d_nu - 0.25) > 1.0e-12)
      std::cout << "Warning: RNP bond-based model fixes nu = 1/4; ignoring "
                   "nu = " << deck.d_matData.d_nu << ".\n";
    deck.d_matData.d_nu = 0.25;
 
    // compute E if not provided or K if not provided
    if (deck.d_matData.d_E > 0.)
      deck.d_matData.d_K =
          deck.d_matData.toK(deck.d_matData.d_E, deck.d_matData.d_nu);
 
    if (deck.d_matData.d_K > 0. && deck.d_matData.d_E < 0.)
      deck.d_matData.d_E =
          deck.d_matData.toE(deck.d_matData.d_K, deck.d_matData.d_nu);
 
    if (deck.d_matData.d_Gc > 0.)
      deck.d_matData.d_KIc = deck.d_matData.toKIc(
          deck.d_matData.d_Gc, deck.d_matData.d_nu, deck.d_matData.d_E);
 
    if (deck.d_matData.d_KIc > 0. && deck.d_matData.d_Gc < 0.)
      deck.d_matData.d_Gc = deck.d_matData.toGc(
          deck.d_matData.d_KIc, deck.d_matData.d_nu, deck.d_matData.d_E);
 
    // compute lame parameter
    deck.d_matData.d_lambda =
        deck.d_matData.toLambdaE(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_G =
        deck.d_matData.toGE(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_mu = deck.d_matData.d_G;
 
    // get moment of influence function
    double M = getMoment(dim);
 
    // Small strain gives lambda = mu = C beta M / 2 (2D) and 2 C beta M / 5
    // (3D); Gc = 4 C M / pi (2D) and 3 C M / 2 (3D). Jha 2025 eq. 83 takes
    // Lipton's lambda, which is half the physical one (Lipton 2014, eq. 3.10).
    if (dim == 2) {
      d_C = M_PI * deck.d_matData.d_Gc / (4. * M);
      d_beta = 2. * deck.d_matData.d_lambda / (d_C * M);
    } else if (dim == 3) {
      d_C = 2. * deck.d_matData.d_Gc / (3. * M);
      d_beta = 5. * deck.d_matData.d_lambda / (2. * d_C * M);
    }
 
    d_rbar = std::sqrt(0.5 / d_beta);
  };
 
private:
 
  double d_horizon;
 
  double d_density;
 
  double d_C;
 
  double d_beta;
 
  double d_rbar;
 
  double d_invFactor;
 
  double d_factorSc;
 
  bool d_irrevBondBreak;
};
 
class PmbMaterial : public Material {
 
public:
  PmbMaterial(inp::MaterialDeck &deck, const size_t &dim, const double &horizon)
      : Material("PMBBond"), d_horizon(horizon), d_density(deck.d_density),
        d_c(0.), d_s0(0.), d_Jscale(1.) {
 
    // set global fields
    if (dimension != dim)
      dimension = dim;
 
    if (is_plane_strain != deck.d_isPlaneStrain)
      is_plane_strain = deck.d_isPlaneStrain;
 
    // create influence function
    if (deck.d_influenceFnType == 0) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::ConstInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 1) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::LinearInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 2) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::GaussianInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Influence function type = "
          << deck.d_influenceFnType
          << " is invalid.\n");
    }
 
    // check if we need to compute the material parameters
    if (deck.d_computeParamsFromElastic)
      computeParameters(deck, dim);
    else {
      d_c = deck.d_bondPotentialParams[0];
      d_s0 = deck.d_bondPotentialParams[1];
    }
  };
 
  bool isStateActive() const override { return false; };
 
  std::pair<double, double> getBondEF(const double &r, const double &s,
                                      bool &fs,
                                      const bool &break_bonds) const override {
 
    if (!break_bonds)
      return std::make_pair(getInfFn(r) * 0.25 * d_c * s *
                                s * r,
                            getInfFn(r) * d_c * s);
 
    // Fracture: literature PMB is tension-only (s > s0). Absolute stretch is
    // applied in pdForce when Model.Bond_Break = absolute_stretch.
    if (!fs && util::isGreater(s, d_s0 + 1.0e-10))
      fs = true;
 
    // if bond is not fractured, return energy and force from nonlinear
    // potential otherwise return energy of fractured bond, and zero force
    if (!fs)
      return std::make_pair(getInfFn(r) * 0.25 * d_c * s *
                                s * r,
                            getInfFn(r) * d_c * s);
    else
      return std::make_pair(
          getInfFn(r) * 0.25 * d_c * d_s0 * d_s0 * r, 0.);
  };
 
  std::pair<double, double>
  getBondEF(const double &r, const double &s, bool &fs, const double
  &mx, const double &thetax) const override {
 
    return this->getBondEF(r, s, fs, true);
  };
 
  util::Point getBondForceDirection(const util::Point &dx,
                                     const util::Point &du) const override {
    return (dx + du) / (dx + du).length();
  };
 
  double getS(const util::Point &dx, const util::Point &du) const override {
    return ((dx + du).length() - dx.length()) / dx.length();
  };
 
  double getSc(const double &r) const override { return d_s0; };
 
  double getDensity() const override { return d_density; };
 
  double getInfFn(const double &r) const override {
    return getGlobalInfFn(r / d_horizon);
  };
 
  double getMoment(const size_t &i) const override {
    return getGlobalMoment(i);
  };
 
  double getHorizon() const override { return d_horizon; };
 
  inp::MatData computeMaterialProperties(const size_t &dim) const override {
 
    auto data = inp::MatData();
    data.d_nu = (dim == 2 && !is_plane_strain) ? 1. / 3. : 0.25;
 
    // inverse of computeParameters, with c for J = 1
    const double h = d_horizon;
    const double c = d_c * d_Jscale;
    if (dim == 2) {
      data.d_E = is_plane_strain ? 5. * M_PI * std::pow(h, 3) * c / 48.
                                 : M_PI * std::pow(h, 3) * c / 9.;
      data.d_K = data.toK(data.d_E, data.d_nu);
      data.d_Gc = c * d_s0 * d_s0 * std::pow(h, 4) / 4.;
    } else if (dim == 3) {
      data.d_K = M_PI * std::pow(h, 4) * c / 18.;
      data.d_E = data.toE(data.d_K, data.d_nu);
      data.d_Gc = M_PI * c * d_s0 * d_s0 * std::pow(h, 5) / 10.;
    }
    data.d_G = data.toGE(data.d_E, data.d_nu);
    data.d_lambda = data.toLambdaE(data.d_E, data.d_nu);
    data.d_mu = data.d_G;
    data.d_KIc = data.toKIc(data.d_Gc, data.d_nu, data.d_E);
    return data;
  };
 
  std::string printStr(int nt, int lvl) const override {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- particle::PmbMaterial --------" << std::endl
        << std::endl;
    oss << tabS << "State active = " << 0 << std::endl;
    oss << tabS << "Horizon = " << d_horizon << std::endl;
    oss << tabS << "Influence fn address = " << influence_fn.get() << std::endl;
    oss << tabS << "Influence fn info: " << std::endl;
    oss << influence_fn->printStr(nt + 1, lvl);
    oss << tabS << "Peridynamic parameters: " << std::endl;
    oss << tabS << "  c = " << d_c << std::endl;
    oss << tabS << "  s0 = " << d_s0 << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  };
 
  void print(int nt, int lvl) const override {
    std::cout << printStr(nt, lvl);
  };
 
  void print() const override { print(0, 0); };
 
private:
  void computeParameters(inp::MaterialDeck &deck, const size_t &dim) {
    //
    // Need following elastic and fracture properties
    // 1. E or K
    // 2. Gc or KIc
    // Bond-based Poisson ratio is locked by dimension / plane-strain flag
    // (Silling & Askari 2005; Madenci & Oterkus 2014; Ha & Bobaru 2010).
    //
    if (util::isLess(deck.d_matData.d_E, 0.) &&
        util::isLess(deck.d_matData.d_K, 0.)) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Require either Young's modulus E or Bulk modulus K"
          " to compute the PMB bond-based peridynamic parameters.\n");
    }
    if (util::isGreater(deck.d_matData.d_E, 0.) &&
        util::isGreater(deck.d_matData.d_K, 0.)) {
      std::cout << "Warning: Both Young's modulus E and Bulk modulus K are "
                   "provided.\n";
      std::cout << "Warning: Selecting Young's modulus E for PMB calibration.\n";
    }
 
    if (util::isLess(deck.d_matData.d_Gc, 0.) &&
        util::isLess(deck.d_matData.d_KIc, 0.)) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Require either critical energy release rate Gc or "
          "critical stress intensity factor KIc to compute the PMB "
          "bond-based peridynamic parameters.\n");
    } else if (util::isGreater(deck.d_matData.d_Gc, 0.) &&
               util::isGreater(deck.d_matData.d_KIc, 0.)) {
      std::cout << "Warning: Both Gc and KIc provided; selecting Gc.\n";
    }
 
    // Bond-based PMB fixes nu: 1/4 in 3D and plane strain, 1/3 in plane
    // stress. E is matched.
    const double nu_pmb = (dim == 2 && !is_plane_strain) ? 1. / 3. : 0.25;
    if (util::isGreater(deck.d_matData.d_nu, 0.) &&
        std::abs(deck.d_matData.d_nu - nu_pmb) > 1.0e-12)
      std::cout << "Warning: PMB bond-based model fixes nu = " << nu_pmb
                << "; ignoring nu = " << deck.d_matData.d_nu << ".\n";
    deck.d_matData.d_nu = nu_pmb;
 
    // Prefer E; else recover E from K with the fixed nu.
    if (deck.d_matData.d_E < 0. && deck.d_matData.d_K > 0.)
      deck.d_matData.d_E =
          deck.d_matData.toE(deck.d_matData.d_K, deck.d_matData.d_nu);
    if (deck.d_matData.d_E > 0.)
      deck.d_matData.d_K =
          deck.d_matData.toK(deck.d_matData.d_E, deck.d_matData.d_nu);
 
    if (deck.d_matData.d_Gc > 0.)
      deck.d_matData.d_KIc = deck.d_matData.toKIc(
          deck.d_matData.d_Gc, deck.d_matData.d_nu, deck.d_matData.d_E);
 
    if (deck.d_matData.d_KIc > 0. && deck.d_matData.d_Gc < 0.)
      deck.d_matData.d_Gc = deck.d_matData.toGc(
          deck.d_matData.d_KIc, deck.d_matData.d_nu, deck.d_matData.d_E);
 
    deck.d_matData.d_lambda =
        deck.d_matData.toLambdaE(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_G =
        deck.d_matData.toGE(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_mu = deck.d_matData.d_G;
 
    // The closed forms for c and s0 hold for a constant influence function.
    if (deck.d_influenceFnType != 0) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: PMB parameters from E and Gc need a constant influence "
             "function. Set Influence_Function Type = 0.\n");
    }
    d_Jscale = deck.d_influenceFnParams.empty() ? double(dim + 1)
                                                : deck.d_influenceFnParams[0];
    if (!(d_Jscale > 0.)) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: PMB influence scale must be > 0.\n");
    }
 
    // c and s0 for J = 1: Silling and Askari 2005 (3D), Ha and Bobaru 2010
    // (2D plane stress); plane strain uses the same energy match with
    // nu = 1/4. Then J = a0 is absorbed into c.
    const double Gc = deck.d_matData.d_Gc;
    const double E = deck.d_matData.d_E;
    const double h = d_horizon;
    if (dim == 2) {
      d_c = is_plane_strain ? 48. * E / (5. * M_PI * std::pow(h, 3))
                            : 9. * E / (M_PI * std::pow(h, 3));
      d_s0 = std::sqrt(4. * Gc / (d_c * std::pow(h, 4)));
    } else if (dim == 3) {
      d_c = 18. * deck.d_matData.d_K / (M_PI * std::pow(h, 4));
      d_s0 = std::sqrt(10. * Gc / (M_PI * d_c * std::pow(h, 5)));
    } else {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: PMB computeParameters: unsupported dim=" << dim
          << "\n");
    }
    d_c /= d_Jscale;
  };
 
private:
  double d_horizon;
 
  double d_density;
 
  double d_c;
 
  double d_s0;
 
  double d_Jscale;
 
};
 
class PdElastic : public Material {
 
public:
  PdElastic(inp::MaterialDeck &deck, const size_t &dim, const double &horizon)
      : Material("PDElasticBond"), d_horizon(horizon),
        d_density(deck.d_density), d_c(0.) {
 
    // set global fields
    if (dimension != dim)
      dimension = dim;
 
    if (is_plane_strain != deck.d_isPlaneStrain)
      is_plane_strain = deck.d_isPlaneStrain;
 
    // create influence function
    if (deck.d_influenceFnType == 0) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::ConstInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 1) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::LinearInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 2) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::GaussianInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Influence function type = "
          << deck.d_influenceFnType
          << " is invalid.\n");
    }
 
    // check if we need to compute the material parameters
    if (deck.d_computeParamsFromElastic)
      computeParameters(deck, dim);
    else
      d_c = deck.d_bondPotentialParams[0];
  };
 
  bool isStateActive() const override { return false; };
 
  std::pair<double, double> getBondEF(const double &r, const double &s,
                                      bool &fs,
                                      const bool &break_bonds) const override {
 
    return std::make_pair(getInfFn(r) * 0.5 * d_c * s *
                            s * r,
                          getInfFn(r) * d_c * s);
  };
 
  std::pair<double, double>
  getBondEF(const double &r, const double &s, bool &fs, const double
  &mx, const double &thetax) const override {
 
    return this->getBondEF(r, s, fs, true);
  };
 
  util::Point getBondForceDirection(const util::Point &dx,
                                     const util::Point &du) const override {
    return (dx + du) / (dx + du).length();
  };
 
  double getS(const util::Point &dx, const util::Point &du) const override {
    return ((dx + du).length() - dx.length()) / dx.length();
  };
 
  double getSc(const double &r) const override { return std::numeric_limits<double>::max(); };
 
  double getDensity() const override { return d_density; };
 
  double getInfFn(const double &r) const override {
    return getGlobalInfFn(r / d_horizon);
  };
 
  double getMoment(const size_t &i) const override {
    return getGlobalMoment(i);
  };
 
  double getHorizon() const override { return d_horizon; };
 
  inp::MatData computeMaterialProperties(const size_t &dim) const override {
 
    auto data = inp::MatData();
 
    if (dim == 2 && !is_plane_strain)
      data.d_nu = 1.0 / 3.0;
    else
      data.d_nu = 0.25;
 
    const double h = d_horizon;
    if (dim == 2 && !is_plane_strain) {
      data.d_E = d_c * (M_PI * std::pow(h, 3.0)) / 9.0;
    } else if (dim == 2 && is_plane_strain) {
      data.d_E = d_c * (5.0 * M_PI * std::pow(h, 3.0)) / 48.0;
    } else if (dim == 3) {
      data.d_K = d_c * (M_PI * std::pow(h, 4.0)) / 18.0;
      data.d_E = data.toE(data.d_K, data.d_nu);
    }
    if (dim == 2)
      data.d_K = data.toK(data.d_E, data.d_nu);
    data.d_lambda = data.toLambdaE(data.d_E, data.d_nu);
    data.d_mu = data.toGE(data.d_E, data.d_nu);
    data.d_G = data.d_mu;
 
    return data;
  };
 
  std::string printStr(int nt, int lvl) const override {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- particle::PdElastic --------" << std::endl
        << std::endl;
    oss << tabS << "State active = " << 0 << std::endl;
    oss << tabS << "Horizon = " << d_horizon << std::endl;
    oss << tabS << "Influence fn address = " << influence_fn.get() << std::endl;
    oss << tabS << "Influence fn info: " << std::endl;
    oss << influence_fn->printStr(nt + 1, lvl);
    oss << tabS << "Peridynamic parameters: " << std::endl;
    oss << tabS << "  c = " << d_c << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  };
 
  void print(int nt, int lvl) const override {
    std::cout << printStr(nt, lvl);
  };
 
  void print() const override { print(0, 0); };
 
private:
  void computeParameters(inp::MaterialDeck &deck, const size_t &dim) {
    // Elastic-only bond-based micromodulus (same δ powers as PMB).
    if (util::isLess(deck.d_matData.d_E, 0.) &&
        util::isLess(deck.d_matData.d_K, 0.)) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Require either Young's modulus E or Bulk modulus K"
          " to compute the PdElastic parameters.\n");
    }
    if (util::isGreater(deck.d_matData.d_E, 0.) &&
        util::isGreater(deck.d_matData.d_K, 0.)) {
      std::cout << "Warning: Both E and K provided; selecting E for PdElastic.\n";
    }
 
    // Emmrich / Trask micromodulus (J≡1). Scale if ConstInfluence a0≠1.
    deck.d_matData.d_nu = 0.25;
    if (deck.d_matData.d_E < 0. && deck.d_matData.d_K > 0.)
      deck.d_matData.d_E =
          deck.d_matData.toE(deck.d_matData.d_K, deck.d_matData.d_nu);
    if (deck.d_matData.d_E > 0.)
      deck.d_matData.d_K =
          deck.d_matData.toK(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_lambda =
        deck.d_matData.toLambdaE(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_G =
        deck.d_matData.toGE(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_mu = deck.d_matData.d_G;
 
    double J_scale = 1.0;
    if (deck.d_influenceFnType == 0)
      J_scale = deck.d_influenceFnParams.empty() ? double(dim + 1)
                                                 : deck.d_influenceFnParams[0];
    else if (deck.d_influenceFnType == 1 || deck.d_influenceFnType == 2)
      J_scale = std::max(getMoment(0), 1.0e-30);
 
    const double h = d_horizon;
    const double kappa = deck.d_matData.d_K;
    if (dim == 2)
      d_c = 72.0 * kappa / (5.0 * M_PI * std::pow(h, 3.0));
    else if (dim == 3)
      d_c = 18.0 * kappa / (M_PI * std::pow(h, 4.0));
    else {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: PdElastic computeParameters: unsupported dim=" << dim
          << "\n");
    }
    d_c /= J_scale;
  };
 
private:
  double d_horizon;
 
  double d_density;
 
  double d_c;
 
};
 
class PdState : public Material {
 
public:
  PdState(inp::MaterialDeck &deck, const size_t &dim, const double &horizon)
      : Material("PDState"), d_horizon(horizon), d_density(deck.d_density),
        d_K(0.), d_G(0.), d_kappa(0.), d_s0(0.), d_dim(double(dim)),
        d_alphaFactor(dim == 2 ? 8. : 15.) {
 
    // set global fields
    if (dimension != dim)
      dimension = dim;
 
    if (is_plane_strain != deck.d_isPlaneStrain)
      is_plane_strain = deck.d_isPlaneStrain;
 
    // create influence function
    if (deck.d_influenceFnType == 0) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::ConstInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 1) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::LinearInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else if (deck.d_influenceFnType == 2) {
      if (influence_fn == nullptr)
        influence_fn = std::make_shared<material::GaussianInfluenceFn>(
            deck.d_influenceFnParams, dim);
    }
    else {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Influence function type = "
          << deck.d_influenceFnType
          << " is invalid.\n");
    }
 
    // check if we need to compute the material parameters
    if (deck.d_computeParamsFromElastic)
      computeParameters(deck, dim);
    else {
      d_K = deck.d_bondPotentialParams[0];
      d_G = deck.d_bondPotentialParams[1];
      d_s0 = deck.d_bondPotentialParams[2];
      d_kappa = getModelBulkModulus(dim, is_plane_strain, d_K, d_G);
    }
  };
 
  bool isStateActive() const override { return true; };
 
  std::pair<double, double> getBondEF(const double &r, const double &s,
                                      bool &fs,
                                      const bool &break_bonds) const override {
 
    return {0., 0.};
  };
 
  std::pair<double, double>
  getBondEF(const double &r, const double &s, bool &fs, const double
  &mx, const double &thetax) const override {
 
    if (fs)
      return {0., 0.};
 
    double J = getInfFn(r);
    double change_length = s * r;
 
    double alpha = d_alphaFactor * d_G / mx;
    double factor = (d_dim * d_kappa / mx) - alpha / d_dim;
    double e_dev = change_length - thetax * r / d_dim;
 
    return {0.5 * alpha * J * e_dev * e_dev,
            J * (r * thetax * factor + change_length * alpha)};
  };
 
  double getDilatationEnergyDensity(const double &thetax) const override {
    return 0.5 * d_kappa * thetax * thetax;
  };
 
  util::Point getBondForceDirection(const util::Point &dx,
                                     const util::Point &du) const override {
    return (dx + du) / (dx + du).length();
  };
 
  double getS(const util::Point &dx, const util::Point &du) const override {
    return ((dx + du).length() - dx.length()) / dx.length();
  };
 
  double getSc(const double &r) const override { return d_s0; };
 
  double getDensity() const override { return d_density; };
 
  double getInfFn(const double &r) const override {
    return getGlobalInfFn(r / d_horizon);
  };
 
  double getMoment(const size_t &i) const override {
    return getGlobalMoment(i);
  };
 
  double getHorizon() const override { return d_horizon; };
 
  inp::MatData computeMaterialProperties(const size_t &dim) const override {
 
    auto data = inp::MatData();
 
    // we already have G and K
    data.d_G = d_G;
    data.d_K = d_K;
 
    // get Poisson ratio and Young's modulus
    data.d_nu = (3. * d_K - 2. * d_G) / (2. * (3. * d_K + d_G));
    data.d_E = data.toE(d_K, data.d_nu);
 
    // get lame parameters
    data.d_lambda = data.toLambdaE(data.d_E, data.d_nu);
    data.d_mu = d_G;
 
    // get Gc from s0 (inverse of computeParameters)
    data.d_Gc = d_s0 * d_s0 * getCriticalStretchDensity(dim);
 
    // KIc
    data.d_KIc = data.toKIc(
        data.d_Gc, data.d_nu, data.d_E);
 
    return data;
  };
 
  std::string printStr(int nt, int lvl) const override {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- particle::PdState --------" << std::endl
        << std::endl;
    oss << tabS << "State active = " << 1 << std::endl;
    oss << tabS << "Horizon = " << d_horizon << std::endl;
    oss << tabS << "Influence fn address = " << influence_fn.get() << std::endl;
    oss << tabS << "Influence fn info: " << std::endl;
    oss << influence_fn->printStr(nt + 1, lvl);
    oss << tabS << "Peridynamic parameters: " << std::endl;
    oss << tabS << "  K = " << d_K << std::endl;
    oss << tabS << "  G = " << d_G << std::endl;
    oss << tabS << "  s0 = " << d_s0 << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  };
 
  void print(int nt, int lvl) const override {
    std::cout << printStr(nt, lvl);
  };
 
  void print() const override { print(0, 0); };
 
private:
  void computeParameters(inp::MaterialDeck &deck, const size_t &dim) {
    //
    // Need following elastic and fracture properties
    // 1. E or K
    // 2. Poisson ratio or shear modulus
    // 3. Gc or KIc
    //
    bool found_E = false;
    bool found_K = false;
    bool found_G = false;
    bool found_nu = false;
    size_t num_props = 0;
 
    found_E = util::isGreater(deck.d_matData.d_E, 0.);
    if (found_E)
      num_props++;
 
    found_K = util::isGreater(deck.d_matData.d_K, 0.);
    if (found_K)
      num_props++;
 
    found_G = util::isGreater(deck.d_matData.d_G, 0.);
    if (found_G)
      num_props++;
 
    found_nu = util::isGreater(deck.d_matData.d_nu, 0.);
    if (found_nu)
      num_props++;
 
    if (num_props != 2) {
      std::ostringstream oss;
      oss << "Error: Require two different elastic properties for the "
                   "PdState material. Pairs supported are (E, K), (E, G), "
                   "(E, nu), (K, G).\n";
      oss << deck.printStr(0, 0) << "\n";
      throw std::runtime_error(
          util::io::Msg()
          << oss.str());
    }
 
    if (util::isLess(deck.d_matData.d_Gc, 0.) &&
        util::isLess(deck.d_matData.d_KIc, 0.)) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: Require either critical energy release rate Gc or "
          "critical stress intensity factor KIc to compute the RNP "
          "bond-based peridynamic parameters.\n");
    } else if (util::isGreater(deck.d_matData.d_Gc, 0.) &&
               util::isGreater(deck.d_matData.d_KIc, 0.)) {
      std::cout << "Warning: Both critical energy release rate Gc and critical "
                   "stress intensity factor KIc are provided.\n";
      std::cout << "Warning: To compute the RNP bond-based peridynamic "
                   "parameters, we only require one of those.\n";
      std::cout << "Warning: Selecting critical energy release rate Gc to "
                   "compute parameters.\n";
    }
 
    // compute nu if not provided
    if (!found_nu) {
      if (found_E and found_G)
        deck.d_matData.d_nu =
            (0.5 * deck.d_matData.d_E / deck.d_matData.d_G) - 1.;
 
      if (found_E and found_K)
        deck.d_matData.d_nu =
            (3. * deck.d_matData.d_K - deck.d_matData.d_E) /
            (6. * deck.d_matData.d_K);
 
      if (found_G and found_K)
        deck.d_matData.d_nu =
            (3. * deck.d_matData.d_K - 2. * deck.d_matData.d_G) /
            (2. * (3. * deck.d_matData.d_K + deck.d_matData.d_G));
    }
    found_nu = true;
 
    // compute E if not provided
    if (!found_E) {
      if (found_K)
        deck.d_matData.d_E =
            deck.d_matData.toE(deck.d_matData.d_K, deck.d_matData.d_nu);
 
      if (found_G)
        deck.d_matData.d_E =
            2. * deck.d_matData.d_G * (1. + deck.d_matData.d_nu);
    }
    found_E = true;
 
    // compute K if not provided
    if (!found_K)
      deck.d_matData.d_K =
          deck.d_matData.toK(deck.d_matData.d_E, deck.d_matData.d_nu);
 
    found_K = true;
 
 
    // compute G if not provided
    if (!found_G)
      deck.d_matData.d_G =
          deck.d_matData.toGE(deck.d_matData.d_E, deck.d_matData.d_nu);
 
    found_G = true;
 
    // compute Gc (if not provided) and KIc (if not provided)
    if (deck.d_matData.d_Gc > 0.)
      deck.d_matData.d_KIc = deck.d_matData.toKIc(
          deck.d_matData.d_Gc, deck.d_matData.d_nu, deck.d_matData.d_E);
 
    if (deck.d_matData.d_KIc > 0. && deck.d_matData.d_Gc < 0.)
      deck.d_matData.d_Gc = deck.d_matData.toGc(
          deck.d_matData.d_KIc, deck.d_matData.d_nu, deck.d_matData.d_E);
 
    // compute lame parameter
    deck.d_matData.d_lambda =
        deck.d_matData.toLambdaE(deck.d_matData.d_E, deck.d_matData.d_nu);
    deck.d_matData.d_mu = deck.d_matData.d_G;
 
    // compute peridynamic parameters
    d_K = deck.d_matData.d_K;
    d_G = deck.d_matData.d_G;
    d_kappa = getModelBulkModulus(dim, is_plane_strain, d_K, d_G);
 
    // The closed form for s0 holds for a constant influence function.
    if (deck.d_influenceFnType != 0) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: PDState critical stretch from Gc needs a constant "
             "influence function. Set Influence_Function Type = 0.\n");
    }
    if (dim != 2 && dim != 3) {
      throw std::runtime_error(
          util::io::Msg()
          << "Error: PdState computeParameters: unsupported dim=" << dim
          << "\n");
    }
    d_s0 = std::sqrt(deck.d_matData.d_Gc / getCriticalStretchDensity(dim));
  };
 
  double getCriticalStretchDensity(const size_t &dim) const {
    if (dim == 2)
      return (6.0 * d_G / M_PI +
              16.0 * (d_kappa - 2.0 * d_G) / (9.0 * M_PI * M_PI)) *
             d_horizon;
    return (3. * d_G + std::pow(3. / 4., 4) * (d_K - 5. * d_G / 3.)) *
           d_horizon;
  };
 
private:
  double d_horizon;
 
  double d_density;
 
  double d_K;
 
  double d_G;
 
  double d_kappa;
 
  double d_s0;
 
  double d_dim;
 
  double d_alphaFactor;
 
};
 
} // namespace material
 
#endif // MATERIAL_PD_MATERIAL_H