material::RnpMaterial Class Reference

A class providing methods to compute energy density and force of peridynamic material. More...

#include <material.h>

Inheritance diagram for material::RnpMaterial:

Collaboration diagram for material::RnpMaterial:

Public Member Functions
	RnpMaterial (inp::MaterialDeck &deck, const size_t &dim, const double &horizon)
	Constructor.
bool	isStateActive () const override
	Returns true if state-based potential is active.
std::pair< double, double >	getBondEF (const double &r, const double &s, bool &fs, const bool &break_bonds) const override
	Returns energy and force between bond due to pairwise interaction.
std::pair< double, double >	getBondEF (const double &r, const double &s, bool &fs, const double &mx, const double &thetax) const override
	Returns energy and force between bond due to state-based model.
util::Point	getBondForceDirection (const util::Point &dx, const util::Point &du) const override
	Returns the unit vector along which bond-force acts.
double	getS (const util::Point &dx, const util::Point &du) const override
	Returns the bond strain.
double	getSc (const double &r) const override
	Returns critical bond strain.
double	getBreakSc (const double &r) const override
	Returns bond strain beyond which the bond is marked broken.
double	getDensity () const override
	Returns the density of the material.
double	getInfFn (const double &r) const override
	Returns the value of influence function.
double	getMoment (const size_t &i) const override
	Returns the moment of influence function.
double	getHorizon () const override
	Returns horizon.
inp::MatData	computeMaterialProperties (const size_t &dim) const override
	Computes elastic and fracture material properties and returns the data.
std::string	printStr (int nt, int lvl) const override
	Returns the string containing printable information about the object.
void	print (int nt, int lvl) const override
	Prints the information about the object.
void	print () const override
	Prints the information about the object.
Public Member Functions inherited from material::Material
	Material (std::string name="")
	Constructor.
virtual	~Material ()
	Destructor.
std::string	name ()
	Returns name of the material.
size_t	getDimension () const
	Returns dimension of the problem.
bool	isPlaneStrain () const
	Returns plane-strain condition.
virtual double	getDilatationEnergyDensity (const double &thetax) const
	Returns the dilatational part of the strain energy density.

Private Member Functions
void	computeParameters (inp::MaterialDeck &deck, const size_t &dim)
	Compute material model parameters.

Private Attributes
double	d_horizon
	Horizon.
double	d_density
	Density.
double	d_rbar
	Inflection point of nonlinear function = \( 1/\sqrt{2\beta}\).
double	d_invFactor
	Inverse of factor = \( \epsilon |B_\epsilon(0)|\).
double	d_factorSc
	Factor to multiply to critical strain to check if bond is fractured.
bool	d_irrevBondBreak
	Flag which indicates if the breaking of bond is irreversible.
Material parameters
double	d_C
	Parameter C.
double	d_beta
	Parameter \( \beta \).

Detailed Description

A class providing methods to compute energy density and force of peridynamic material.

Definition at line 293 of file material.h.

Constructor & Destructor Documentation

RnpMaterial()

material::RnpMaterial::RnpMaterial ( inp::MaterialDeck &  deck, const size_t &  dim, const double &  horizon  )

inline

Constructor.

Parameters

deck	Input deck which contains user-specified information
dim	Dimension
horizon	Horizon

Definition at line 302 of file material.h.

      : 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);
    }
  };

References computeParameters(), d_beta, inp::MaterialDeck::d_bondPotentialParams, d_C, inp::MaterialDeck::d_computeParamsFromElastic, inp::MaterialDeck::d_influenceFnParams, inp::MaterialDeck::d_influenceFnType, d_invFactor, inp::MaterialDeck::d_isPlaneStrain, and d_rbar.

Here is the call graph for this function:

Member Function Documentation

computeMaterialProperties()

inp::MatData material::RnpMaterial::computeMaterialProperties ( const size_t &  dim ) const

inlineoverridevirtual

Computes elastic and fracture material properties and returns the data.

Parameters

dim	Dimension of the problem

Returns

Data Material data

Implements material::Material.

Definition at line 512 of file material.h.

                                                                         {
 
    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;
  };

References d_beta, d_C, and getMoment().

Here is the call graph for this function:

computeParameters()

void material::RnpMaterial::computeParameters ( inp::MaterialDeck &  deck, const size_t &  dim  )

inlineprivate

Compute material model parameters.

Parameters

deck	MaterialDeck
dim	Dimension of the domain

Definition at line 592 of file material.h.

                                                                   {
    //
    // 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);
  };

References d_beta, d_C, inp::MatData::d_E, inp::MatData::d_G, inp::MatData::d_Gc, inp::MatData::d_K, inp::MatData::d_KIc, inp::MatData::d_lambda, inp::MaterialDeck::d_matData, inp::MatData::d_mu, inp::MatData::d_nu, d_rbar, getMoment(), util::isGreater(), util::isLess(), inp::MatData::toE(), inp::MatData::toGc(), inp::MatData::toGE(), inp::MatData::toK(), inp::MatData::toKIc(), and inp::MatData::toLambdaE().

Referenced by RnpMaterial().

Here is the call graph for this function:

Here is the caller graph for this function:

getBondEF() [1/2]

std::pair< double, double > material::RnpMaterial::getBondEF ( const double &  r, const double &  s, bool &  fs, const bool &  break_bonds  ) const

inlineoverridevirtual

Returns energy and force between bond due to pairwise interaction.

Peridynamic energy at point \( x \) is

\[ e(x) = \frac{1}{|B_\epsilon(0)|} \int_{B_\epsilon(x)} \frac{J^\epsilon(|y-x|)}{\epsilon} \psi(|y-x|S^2) dy \]

and force at point x is

\[ f(x) = \frac{4}{|B_\epsilon(0)|} \int_{B_\epsilon(x)} \frac{J^\epsilon(|y-x|)}{\epsilon} \psi'(|y-x|S^2) S \frac{y-x}{|y-x|} dy, \]

where \( \psi(r) = C(1-\exp(-\beta r))\) (Lipton 2016; Jha 2025, eqs. 8-9). The total energy is \( \int e(x) dx \), so the pair force is twice the derivative of the bond energy.

For given initial bond length \( r \) and bond strain \( s\), this function returns pair of

\[ \hat{e} = \frac{J^\epsilon(r)}{\epsilon |B_\epsilon(0)|} \psi(r s^2) \]

and

\[ \hat{f} = \frac{4 J^\epsilon(r) s}{\epsilon |B_\epsilon(0)|} \psi'(r s^2). \]

Parameters

r	Reference (initial) bond length
s	Bond strain
fs	Bond fracture state
break_bonds	If true, also update the fracture state of the bond

Returns

value Pair of energy and force

Implements material::Material.

Definition at line 387 of file material.h.

                                                                              {
 
    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.);
  };

References d_beta, d_C, d_invFactor, getBreakSc(), getInfFn(), and util::isGreater().

Referenced by getBondEF().

Here is the call graph for this function:

Here is the caller graph for this function:

getBondEF() [2/2]

std::pair< double, double > material::RnpMaterial::getBondEF ( const double &  r, const double &  s, bool &  fs, const double &  mx, const double &  thetax  ) const

inlineoverridevirtual

Returns energy and force between bond due to state-based model.

Parameters

r	Reference (initial) bond length
s	Bond strain
fs	Bond fracture state
mx	Weighted volume at node
thetax	Dilation

Returns

value Pair of energy and force

Implements material::Material.

Definition at line 417 of file material.h.

                                            {
 
    return this->getBondEF(r, s, fs, true);
  };

References getBondEF().

Here is the call graph for this function:

getBondForceDirection()

util::Point material::RnpMaterial::getBondForceDirection ( const util::Point &  dx, const util::Point &  du  ) const

inlineoverridevirtual

Returns the unit vector along which bond-force acts.

Parameters

dx	Reference bond vector
du	Difference of displacement

Returns

vector Unit vector

Implements material::Material.

Definition at line 430 of file material.h.

                                                                         {
    return dx / dx.length();
  };

References util::Point::length().

Here is the call graph for this function:

getBreakSc()

double material::RnpMaterial::getBreakSc ( const double &  r ) const

inlineoverridevirtual

Returns bond strain beyond which the bond is marked broken.

The force softens past the critical strain; the bond is marked broken only past factor_Sc times it, and never if breaking is reversible.

Parameters

r	Reference length of bond

Returns

strain Breaking strain

Reimplemented from material::Material.

Definition at line 464 of file material.h.

                                                    {
    if (!d_irrevBondBreak)
      return std::numeric_limits<double>::max();
    return d_factorSc * getSc(r);
  };

References d_factorSc, d_irrevBondBreak, and getSc().

Referenced by getBondEF().

Here is the call graph for this function:

Here is the caller graph for this function:

getDensity()

double material::RnpMaterial::getDensity ( ) const

inlineoverridevirtual

Returns the density of the material.

Returns

density Density of the material

Implements material::Material.

Definition at line 474 of file material.h.

{ return d_density; };

References d_density.

getHorizon()

double material::RnpMaterial::getHorizon ( ) const

inlineoverridevirtual

Returns horizon.

Returns

horizon Horizon

Implements material::Material.

Definition at line 503 of file material.h.

{ return d_horizon; };

References d_horizon.

getInfFn()

double material::RnpMaterial::getInfFn ( const double &  r ) const

inlineoverridevirtual

Returns the value of influence function.

Parameters

r	Reference (initial) bond length

Returns

value Influence function at r

Implements material::Material.

Definition at line 482 of file material.h.

                                                  {
    return getGlobalInfFn(r / d_horizon);
  };

References d_horizon.

Referenced by getBondEF().

Here is the caller graph for this function:

getMoment()

double material::RnpMaterial::getMoment ( const size_t &  i ) const

inlineoverridevirtual

Returns the moment of influence function.

If \( J(r) \) is the influence function for \( r\in [0,1)\) then \( i^{th}\) moment is given by

\[ M_i = \int_0^1 J(r) r^i dr. \]

Parameters

i	ith moment

Returns

value Moment

Implements material::Material.

Definition at line 495 of file material.h.

                                                   {
    return getGlobalMoment(i);
  };

Referenced by computeMaterialProperties(), and computeParameters().

Here is the caller graph for this function:

getS()

double material::RnpMaterial::getS ( const util::Point &  dx, const util::Point &  du  ) const

inlineoverridevirtual

Returns the bond strain.

Parameters

dx	Reference bond vector
du	Difference of displacement

Returns

strain Bond strain \( S = \frac{du \cdot dx}{|dx|^2} \)

Implements material::Material.

Definition at line 441 of file material.h.

                                                                       {
    return dx.dot(du) / dx.dot(dx);
  };

References util::Point::dot().

Here is the call graph for this function:

getSc()

double material::RnpMaterial::getSc ( const double &  r ) const

inlineoverridevirtual

Returns critical bond strain.

Parameters

r	Reference length of bond

Returns

strain Critical strain

Implements material::Material.

Definition at line 451 of file material.h.

                                               {
    return d_rbar / std::sqrt(r);
  };

References d_rbar.

Referenced by getBreakSc().

Here is the caller graph for this function:

isStateActive()

bool material::RnpMaterial::isStateActive ( ) const

inlineoverridevirtual

Returns true if state-based potential is active.

Returns

bool True/false

Implements material::Material.

Definition at line 357 of file material.h.

{ return false; };

print() [1/2]

void material::RnpMaterial::print ( ) const

inlineoverridevirtual

Prints the information about the object.

Reimplemented from material::Material.

Definition at line 583 of file material.h.

{ print(0, 0); }

References print().

Referenced by print().

Here is the call graph for this function:

Here is the caller graph for this function:

print() [2/2]

void material::RnpMaterial::print ( int  nt, int  lvl  ) const

inlineoverridevirtual

Prints the information about the object.

Parameters

nt	Number of tabs to append before printing
lvl	Information level (higher means more information)

Reimplemented from material::Material.

Definition at line 578 of file material.h.

                                             {
    std::cout << printStr(nt, lvl);
  }

References printStr().

Here is the call graph for this function:

printStr()

std::string material::RnpMaterial::printStr ( int  nt, int  lvl  ) const

inlineoverridevirtual

Returns the string containing printable information about the object.

Parameters

nt	Number of tabs to append before printing
lvl	Information level (higher means more information)

Returns

string String containing printable information about the object

Reimplemented from material::Material.

Definition at line 549 of file material.h.

                                                     {
 
    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();
  }

References d_beta, d_C, d_factorSc, d_horizon, d_invFactor, d_irrevBondBreak, d_rbar, and util::io::getTabS().

Referenced by print().

Here is the call graph for this function:

Here is the caller graph for this function:

Field Documentation

d_beta

double material::RnpMaterial::d_beta

private

Parameter \( \beta \).

Definition at line 705 of file material.h.

Referenced by computeMaterialProperties(), computeParameters(), getBondEF(), printStr(), and RnpMaterial().

d_C

double material::RnpMaterial::d_C

private

Parameter C.

Definition at line 702 of file material.h.

Referenced by computeMaterialProperties(), computeParameters(), getBondEF(), printStr(), and RnpMaterial().

d_density

double material::RnpMaterial::d_density

private

Density.

Definition at line 694 of file material.h.

Referenced by getDensity().

d_factorSc

double material::RnpMaterial::d_factorSc

private

Factor to multiply to critical strain to check if bond is fractured.

For nonlinear model, we consider bond is broken when it exceeds 10 times of critical strain. Typical value of factor is 10.

Definition at line 721 of file material.h.

Referenced by getBreakSc(), and printStr().

d_horizon

double material::RnpMaterial::d_horizon

private

Horizon.

Definition at line 691 of file material.h.

Referenced by getHorizon(), getInfFn(), and printStr().

d_invFactor

double material::RnpMaterial::d_invFactor

private

Inverse of factor = \( \epsilon |B_\epsilon(0)|\).

Definition at line 713 of file material.h.

Referenced by getBondEF(), printStr(), and RnpMaterial().

d_irrevBondBreak

bool material::RnpMaterial::d_irrevBondBreak

private

Flag which indicates if the breaking of bond is irreversible.

Definition at line 724 of file material.h.

Referenced by getBreakSc(), and printStr().

d_rbar

double material::RnpMaterial::d_rbar

private

Inflection point of nonlinear function = \( 1/\sqrt{2\beta}\).

Definition at line 710 of file material.h.

Referenced by computeParameters(), getSc(), printStr(), and RnpMaterial().

---

The documentation for this class was generated from the following file:

• /home/user/project/src/material/mparticle/material.h