pd Namespace Reference

Data Structures
class	BrokenBondKnSelfContact
	Default: Kn * volj * capped_gap / R along yji with gap = R - Rc. More...
class	NoneSelfContact
	No force on broken bonds (Trask: broken → weight 0). More...
class	ReferenceGapSelfContact
class	SelfContact
	Intra-body force between nodes when the PD bond is inactive. More...

Functions
void	computeForces (data::ModelData &data)
	Assemble nodal peridynamic force (constitutive response stays in material).
bool	dofMpiEnabled (const data::ModelData &data)
	True when nodal DOF MPI is active (mpiSize > 1, MPI_Strategy=dof).
void	setupDofPartition (data::ModelData &data)
	Metis-partition nodes on the PD neighbor graph and build ghost plans. Call after d_neighPd is ready and before (or while) building d_fPdCompNodes. Works for Single_Particle and Multi_Particle when MPI_Strategy is dof. No-op for single-rank runs or other strategies.
void	exchangeGhostDisplacement (data::ModelData &data)
	Halo-exchange nodal displacements (and current x) for PD ghosts.
void	exchangeGhostTheta (data::ModelData &data)
	Halo-exchange state-based dilatation d_thetaX for PD ghosts.
std::unique_ptr< SelfContact >	makeSelfContact (const std::string &name)
	Build self-contact law named by Model.Self_Contact.

Function Documentation

computeForces()

void pd::computeForces

(

data::ModelData &

data

)

Assemble nodal peridynamic force (constitutive response stays in material).

Definition at line 28 of file pdForce.cpp.

                                          {
 
 
  util::io::log(3, "    Computing peridynamic force \n");
 
  const auto dim = data.d_modelDeck_p->d_dim;
  const bool is_state = data.d_particlesListTypeAll[0]->getMaterial()->isStateActive();
  const auto selfContact =
      pd::makeSelfContact(data.d_modelDeck_p->d_selfContact);
  const bool abs_stretch_break =
      (data.d_modelDeck_p->d_bondBreak == "absolute_stretch");
 
  // compute state-based helper quantities
  if (is_state) {
 
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;
 
    taskflow.for_each_index(
      (std::size_t) 0, data.d_fPdCompNodes.size(), (std::size_t) 1,
      [&data, abs_stretch_break, dim](std::size_t II) {
        auto i = data.d_fPdCompNodes[II];
 
        const auto rho = data.getDensity(i);
        const auto &fix = data.d_fix[i];
        const auto &ptId = data.getPtId(i);
        auto &pi = data.getParticleFromAllList(ptId);
 
        if (pi->d_material_p->isStateActive()) {
 
          const double horizon = pi->getHorizon();
          const double mesh_size = pi->getMeshSize();
          const auto &xi = data.d_xRef[i];
          const auto &ui = data.d_u[i];
 
          // update bond state and compute thetax
          const auto &m = data.d_mX[i];
          double theta = 0.;
 
          // upper and lower bound for volume correction
          auto check_up = horizon + 0.5 * mesh_size;
          auto check_low = horizon - 0.5 * mesh_size;
 
          size_t k = 0;
          for (size_t j : data.d_neighPd[i]) {
 
            const auto &xj = data.d_xRef[j];
            const auto &uj = data.d_u[j];
            double rji = (xj - xi).length();
            if (!(rji > 0.)) {
              k += 1;
              continue;
            }
            // double rji = std::sqrt(data.d_neighPdSqdDist[i][k]);
            double change_length = (xj - xi + uj - ui).length() - rji;
 
            // step 1: update the bond state
            double s = change_length / rji;
            double sc = pi->d_material_p->getSc(rji);
 
            // get fracture state, modify, and set
            auto fs = data.d_fracture_p->getBondState(i, k);
            // Model.Bond_Break: tension (s > sc) or absolute_stretch (|s| > sc).
            if (!fs) {
              const bool broke = abs_stretch_break
                                     ? util::isGreater(std::abs(s), sc + 1.0e-10)
                                     : util::isGreater(s, sc + 1.0e-10);
              if (broke)
                fs = true;
            }
            data.d_fracture_p->setBondState(i, k, fs);
 
            if (!fs) {
 
              // get corrected volume of node j
              auto volj = data.d_vol[j];
 
              if (util::isGreater(rji, check_low))
                volj *= (check_up - rji) / mesh_size;
 
              theta += rji * change_length * pi->d_material_p->getInfFn(rji) *
                        volj;
            } // if bond is not broken
 
            k += 1;
          } // loop over neighbors
 
          // dilatation: 3/m in 3D, 2/m in 2D (Yang et al. 2024, eqs. 8-9)
          data.d_thetaX[i] = double(dim) * theta / m;
        } // if it is state-based
      } // loop over nodes
    ); // for_each
 
    executor.run(taskflow).get();
 
    // Ghost nodes need owner dilatations before the force pass.
    pd::exchangeGhostTheta(data);
  }
 
  // strain energy density per node
  if (data.d_e.size() != data.d_x.size())
    data.d_e.resize(data.d_x.size(), 0.f);
 
  // compute the internal forces
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index(
    (std::size_t) 0, data.d_fPdCompNodes.size(), (std::size_t) 1,
    [&data, selfContact = selfContact.get(), abs_stretch_break](std::size_t II) {
      auto i = data.d_fPdCompNodes[II];
 
      // local variable to hold force
      util::Point force_i = util::Point();
      double scalar_f = 0.;
      double energy_i = 0.;
 
      // for damage
      float Zi = 0.;
      // Silling (2000) damage: phi = 1 - ∫ mu dV' / ∫ dV' over the horizon
      // (reported by Silling 2003 and Trask). Bhattacharya & Lipton instead
      // define damage as the broken-bond count fraction, kept separately.
      double vol_all = 0.;
      double vol_intact = 0.;
      size_t n_bonds = 0;
      size_t n_broken = 0;
 
      const auto rhoi = data.getDensity(i);
      const auto &ptIdi = data.getPtId(i);
      auto &pi = data.getParticleFromAllList(ptIdi);
 
      const double horizon = pi->getHorizon();
      const double mesh_size = pi->getMeshSize();
      const auto &xi = data.d_xRef[i];
      const auto &ui = data.d_u[i];
      const auto &mi = data.d_mX[i];
      const auto &thetai = data.d_thetaX[i];
 
      // upper and lower bound for volume correction
      auto check_up = horizon + 0.5 * mesh_size;
      auto check_low = horizon - 0.5 * mesh_size;
 
      // loop over neighbors
      {
        size_t k = 0;
        for (size_t j : data.d_neighPd[i]) {
          auto fs = data.d_fracture_p->getBondState(i, k);
          const auto &xj = data.d_xRef[j];
          const auto &uj = data.d_u[j];
          auto volj = data.d_vol[j];
          double rji = (xj - xi).length();
          if (!(rji > 0.)) {
            k++;
            continue;
          }
          double Sji = pi->d_material_p->getS(xj - xi, uj - ui);
 
          if (!fs) {
            const auto &mj = data.d_mX[j];
            const auto &thetaj = data.d_thetaX[j];
 
            // get corrected volume of node j
            if (util::isGreater(rji, check_low))
              volj *= (check_up - rji) / mesh_size;
 
            // handle two cases differently
            if (pi->d_material_p->isStateActive()) {
 
              auto ef_i =
                  pi->d_material_p->getBondEF(rji, Sji, fs, mi, thetai);
              auto ef_j =
                  pi->d_material_p->getBondEF(rji, Sji, fs, mj, thetaj);
 
              // compute the contribution of bond force to force at i
              scalar_f = (ef_i.second + ef_j.second) * volj;
              energy_i += ef_i.first * volj;
 
              force_i += scalar_f * pi->d_material_p->getBondForceDirection(
                                        xj - xi, uj - ui);
            } // if state-based
            else {
 
              // Bond-based materials (e.g. PMB) historically broke inside
              // getBondEF using |s|>Sc. Honor Model.Bond_Break here instead:
              // tension (default, literature PMB) vs absolute_stretch.
              const double Sc = pi->d_material_p->getBreakSc(rji);
              const bool broke = abs_stretch_break
                                     ? util::isGreater(std::abs(Sji), Sc + 1.0e-10)
                                     : util::isGreater(Sji, Sc + 1.0e-10);
              if (broke)
                fs = true;
              data.d_fracture_p->setBondState(i, k, fs);
 
              if (!fs) {
                auto ef =
                    pi->d_material_p->getBondEF(rji, Sji, fs, /*break_bonds=*/false);
                scalar_f = ef.second * volj;
                energy_i += ef.first * volj;
                force_i += scalar_f * pi->d_material_p->getBondForceDirection(
                                          xj - xi, uj - ui);
              } else {
                const auto yji = xj + uj - (xi + ui);
                force_i +=
                    selfContact->force(yji, volj, pi->d_Kn, pi->d_Rc, rji);
              }
            } // if bond-based
          }   // if bond not broken
          else {
            const auto yji = xj + uj - (xi + ui);
            force_i +=
                selfContact->force(yji, volj, pi->d_Kn, pi->d_Rc, rji);
          } // if bond is broken
 
          // calculate damage
          auto Sc = pi->d_material_p->getSc(rji);
          if (util::isGreater(std::abs(Sji / Sc), Zi))
            Zi = std::abs(Sji / Sc);
 
          // damage function phi (bond-state weighted by neighbor volume)
          vol_all += data.d_vol[j];
          n_bonds++;
          if (!data.d_fracture_p->getBondState(i, k))
            vol_intact += data.d_vol[j];
          else
            n_broken++;
 
          k++;
        } // loop over neighbors
 
      } // peridynamic force
 
      // update force (we remove any force from
      // previous steps and add peridynamics force)
      data.d_f[i] = force_i;
      if (pi->d_material_p->isStateActive())
        energy_i += pi->d_material_p->getDilatationEnergyDensity(thetai);
      data.d_e[i] = static_cast<float>(energy_i);
 
      data.d_Z[i] = Zi;
      if (!data.d_phi.empty())
        data.d_phi[i] =
            (vol_all > 0.) ? static_cast<float>(1. - vol_intact / vol_all) : 0.f;
      if (!data.d_phiBond.empty())
        data.d_phiBond[i] =
            (n_bonds > 0) ? static_cast<float>(double(n_broken) / double(n_bonds))
                          : 0.f;
    }
  ); // for_each
 
  executor.run(taskflow).get();
 
  // total strain energy over the nodes computed here
  double te = 0.;
  for (auto i : data.d_fPdCompNodes)
    te += double(data.d_e[i]) * data.d_vol[i];
  data.d_te = static_cast<float>(te);
}

References exchangeGhostTheta(), util::parallel::getNThreads(), util::io::log(), and makeSelfContact().

Referenced by PeriDEMModel::computeForces().

Here is the call graph for this function:

Here is the caller graph for this function:

dofMpiEnabled()

bool pd::dofMpiEnabled

(

const data::ModelData &

data

)

True when nodal DOF MPI is active (mpiSize > 1, MPI_Strategy=dof).

Definition at line 29 of file pdMpi.cpp.

                                                {
  return data.d_pdDofMpi;
}

exchangeGhostDisplacement()

void pd::exchangeGhostDisplacement

(

data::ModelData &

data

)

Halo-exchange nodal displacements (and current x) for PD ghosts.

Definition at line 286 of file pdMpi.cpp.

                                                      {
  if (!data.d_pdDofMpi)
    return;
  using clock = std::chrono::steady_clock;
  const auto t0 = clock::now();
 
  // DOF-MPI: each rank owns a subset of nodes. Before contact (and for PD
  // neighbors outside the local halo), every rank needs a consistent full
  // copy of u and v on Multi_Particle. Single_Particle has no inter-body
  // contact — PD halo exchange is enough.
  // Always refresh d_vMag after rewriting d_v (contact search uses vMag).
  if (data.d_input_p && data.d_input_p->isMultiParticle()) {
    syncAllOwnedPoints(data, data.d_u);
    syncAllOwnedPoints(data, data.d_v);
    for (size_t i = 0; i < data.d_x.size(); ++i)
      data.d_x[i] = data.d_xRef[i] + data.d_u[i];
    refreshVMagFromV(data);
  } else {
    exchangePoints(data, data.d_u);
    exchangePoints(data, data.d_v);
    for (int r = 0; r < util::parallel::mpiSize(); ++r) {
      for (int id : data.d_pdGhostNeedFrom[static_cast<size_t>(r)]) {
        const size_t i = static_cast<size_t>(id);
        data.d_x[i] = data.d_xRef[i] + data.d_u[i];
      }
    }
  }
  data.appendKeyData("pd_mpi_disp_exchange_time",
                     util::methods::timeDiff(t0, clock::now()));
}

References anonymous_namespace{pdMpi.cpp}::exchangePoints(), util::parallel::mpiSize(), refreshVMagFromV(), anonymous_namespace{pdMpi.cpp}::syncAllOwnedPoints(), and util::methods::timeDiff().

Referenced by PeriDEMModel::computeForces().

Here is the call graph for this function:

Here is the caller graph for this function:

exchangeGhostTheta()

void pd::exchangeGhostTheta

(

data::ModelData &

data

)

Halo-exchange state-based dilatation d_thetaX for PD ghosts.

Definition at line 317 of file pdMpi.cpp.

                                               {
  if (!data.d_pdDofMpi)
    return;
  if (data.d_thetaX.empty())
    return;
  using clock = std::chrono::steady_clock;
  const auto t0 = clock::now();
  if (data.d_input_p && data.d_input_p->isMultiParticle()) {
    const int rank = util::parallel::mpiRank();
    const size_t n = data.d_thetaX.size();
    std::vector<double> buf(n, 0.);
    for (size_t i = 0; i < n; ++i) {
      if (static_cast<int>(data.d_pdNodePartition[i]) == rank)
        buf[i] = data.d_thetaX[i];
    }
    MPI_Allreduce(MPI_IN_PLACE, buf.data(), static_cast<int>(n), MPI_DOUBLE,
                  MPI_SUM, util::parallel::mpiComm());
    data.d_thetaX.swap(buf);
  } else {
    exchangeDoubles(data, data.d_thetaX);
  }
  data.appendKeyData("pd_mpi_theta_exchange_time",
                     util::methods::timeDiff(t0, clock::now()));
}

References anonymous_namespace{pdMpi.cpp}::exchangeDoubles(), util::parallel::mpiComm(), util::parallel::mpiRank(), and util::methods::timeDiff().

Referenced by computeForces().

Here is the call graph for this function:

Here is the caller graph for this function:

makeSelfContact()

std::unique_ptr< SelfContact > pd::makeSelfContact

(

const std::string &

name

)

Build self-contact law named by Model.Self_Contact.

Definition at line 68 of file selfContact.cpp.

                                                                  {
  if (name == "broken_bond_kn")
    return std::make_unique<BrokenBondKnSelfContact>();
  if (name == "reference_gap")
    return std::make_unique<ReferenceGapSelfContact>();
  if (name == "none")
    return std::make_unique<NoneSelfContact>();
 
  throw std::runtime_error(
      "Unknown Model.Self_Contact '" + name +
      "'. Supported: broken_bond_kn, reference_gap, none.");
}

Referenced by computeForces(), and main().

Here is the caller graph for this function:

setupDofPartition()

void pd::setupDofPartition

(

data::ModelData &

data

)

Metis-partition nodes on the PD neighbor graph and build ghost plans. Call after d_neighPd is ready and before (or while) building d_fPdCompNodes. Works for Single_Particle and Multi_Particle when MPI_Strategy is dof. No-op for single-rank runs or other strategies.

Definition at line 233 of file pdMpi.cpp.

                                              {
  data.d_pdDofMpi = false;
  data.d_pdGrainAligned = false;
  data.d_pdNodePartition.clear();
  data.d_pdGhostNeedFrom.clear();
  data.d_pdGhostServeTo.clear();
 
  if (!util::parallel::isMpiEnabled())
    return;
  const std::string strategy = particle::resolvedMpiStrategy(data);
  if (strategy != "dof")
    return;
  const int size = util::parallel::mpiSize();
  if (size <= 1)
    return;
  if (data.d_neighPd.empty() || data.d_neighPd.size() != data.d_x.size())
    return;
 
  data.d_pdDofMpi = true;
  data.d_pdGrainAligned = false; // DOF-MPI = node owners, never whole-grain brick
  const int rank = util::parallel::mpiRank();
  const size_t n_nodes = data.d_x.size();
 
  // Partition nodes across ranks (graph partition). Walls / Multi_Particle do
  // not change the mode: DOF-MPI always distributes nodes.
  if (rank == 0) {
    mesh::metisGraphPartition("metis_kway", data.d_neighPd,
                              data.d_pdNodePartition,
                              static_cast<size_t>(size));
  } else {
    data.d_pdNodePartition.assign(n_nodes, 0);
  }
  MPI_Bcast(data.d_pdNodePartition.data(), static_cast<int>(n_nodes),
            MPI_UNSIGNED_LONG, 0, util::parallel::mpiComm());
 
  size_t n_owned = 0;
  for (size_t i = 0; i < n_nodes; ++i)
    if (static_cast<int>(data.d_pdNodePartition[i]) == rank)
      ++n_owned;
 
  buildGhostPlan(data);
 
  util::io::print(std::format(
      "DOF-MPI: rank {}/{} owns {}/{} nodes\n", rank, size, n_owned, n_nodes));
}

References anonymous_namespace{pdMpi.cpp}::buildGhostPlan(), util::parallel::isMpiEnabled(), mesh::metisGraphPartition(), util::parallel::mpiComm(), util::parallel::mpiRank(), util::parallel::mpiSize(), util::io::print(), and particle::resolvedMpiStrategy().

Referenced by PeriDEMModel::init().

Here is the call graph for this function:

Here is the caller graph for this function: