model::DEMModel Class Reference

A class for discrete element particle simulation with peridynamic model More...

#include <demModel.h>

Inheritance diagram for model::DEMModel:

Collaboration diagram for model::DEMModel:

Public Member Functions
	DEMModel (inp::Input *deck, std::string modelName="DEMModel")
	Constructor.
void	log (std::ostringstream &oss, int priority=0, bool check_condition=true, int override_priority=-1, bool screen_out=false)
	Prints message if any of these two conditions are true.
void	log (const std::string &str, int priority=0, bool check_condition=true, int override_priority=-1, bool screen_out=false)
	Prints message if any of these two conditions are true.
virtual void	run (inp::Input *deck)
	Main driver to simulate.
Methods to initialize and run model
virtual void	restart (inp::Input *deck)
	Restarts the simulation from previous state.
virtual void	init ()
	Initialize remaining data members.
virtual void	close ()
	Closure operations.
Methods to implement explicit time integration
virtual void	integrate ()
	Perform time integration.
virtual void	integrateStep ()
	Performs one time step.
virtual void	integrateCD ()
	Perform time integration using central-difference scheme.
virtual void	integrateVerlet ()
	Perform time integration using velocity verlet scheme.
Compute methods
virtual void	computeForces ()
	Computes peridynamic forces and contact forces.
virtual void	computePeridynamicForces ()
	Computes peridynamic forces.
virtual void	computeExternalForces ()
	Computes external/boundary condition forces.
virtual void	computeExternalDisplacementBC ()
	Applies displacement boundary conditions.
virtual void	computeContactForces ()
	Computes contact forces.
virtual void	applyInitialCondition ()
	Applies initial condition.
Setup methods
virtual void	createParticles ()
	Creates particles in a given container.
virtual void	createParticlesFromFile (size_t z, std::shared_ptr< particle::RefParticle > ref_p)
	Creates particles in a Hexagonal arrangement.
virtual void	createParticleUsingParticleZoneGeomObject (size_t z, std::shared_ptr< particle::RefParticle > ref_p)
	Creates particles in a given container.
virtual double	createGeometryAtSite (const double &particle_radius, const double &particle_orient, const util::Point &site, const std::vector< double > &rep_geom_params, const std::shared_ptr< util::geometry::GeomObject > &rep_geom_p, std::shared_ptr< util::geometry::GeomObject > &p_geom)
	Creates geometrical object for a particle given particle radius, orientation, and site location.
virtual void	updateGeometryObjectsPostInit ()
	Update varioud geometry objects associated with container, particles, and reference particles.
virtual void	updateContactNeighborlist ()
	Update neighborlist for contact.
virtual void	updatePeridynamicNeighborlist ()
	Update neighborlist for peridynamics force.
virtual void	updateNeighborlistCombine ()
	Update neighborlist for contact and peridynamics force.
virtual void	setupContact ()
	Creates particles in a given container.
virtual void	setupQuadratureData ()
	Sets up quadrature data.
virtual bool	updateContactNeighborSearchParameters ()
	Update contact neighbor search parameters.
Methods to handle output and debug
virtual void	output ()
	Output the snapshot of data at current time step.
std::string	ppTwoParticleTest ()
	Function that handles post-processing for two particle collision test and returns maximum vertical displacement of particle.
std::string	ppCompressiveTest ()
	Function that handles post-processing for compressive test of particulate media by rigid wall and returns wall penetration and reaction force on wall.
virtual void	checkStop ()
	Checks if simulation should be stopped due to abnormal state of system.
Public Member Functions inherited from model::ModelData
	ModelData (inp::Input *deck)
	Constructor.
const particle::BaseParticle *	getParticleFromAllList (size_t i) const
	Get pointer to base particle.
particle::BaseParticle *&	getParticleFromAllList (size_t i)
	Get pointer to base particle.
const particle::BaseParticle *	getParticleFromParticleList (size_t i) const
	Get pointer to particle (excluding wall)
particle::BaseParticle *&	getParticleFromParticleList (size_t i)
const particle::BaseParticle *	getParticleFromWallList (size_t i) const
	Get pointer to wall.
particle::BaseParticle *&	getParticleFromWallList (size_t i)
double	getDensity (size_t i)
	Get density of particle.
double	getHorizon (size_t i)
	Get horizon of particle.
size_t &	getPtId (size_t i)
	Get particle id given the location in particle list.
const size_t &	getPtId (size_t i) const
	Get particle id given the location in particle list.
void	setPtId (size_t i, const size_t &id)
	Set particle id given the location in particle list.
double	getKeyData (std::string key, bool issue_err=false)
	Get data for a key.
void	appendKeyData (std::string key, double data, bool issue_err=false)
	Append value to data associated with key.
void	setKeyData (std::string key, double data, bool issue_err=false)
	Set value to data associated with key.
util::Point &	getXRef (size_t i)
	Get reference coordinate of the node.
const util::Point &	getXRef (size_t i) const
	Get reference coordinate of the node.
void	setXRef (size_t i, const util::Point &x)
	Set reference coordinate of the node.
void	addXRef (size_t i, const util::Point &x)
	Add reference coordinate of the node.
void	setXRef (size_t i, int dof, double x)
	Set specific reference coordinate of the node.
void	addXRef (size_t i, int dof, double x)
	Add specific reference coordinate of the node.
util::Point &	getX (size_t i)
	Get current coordinate of the node.
const util::Point &	getX (size_t i) const
	Get current coordinate of the node.
void	setX (size_t i, const util::Point &x)
	Set current coordinate of the node.
void	addX (size_t i, const util::Point &x)
	Add current coordinate of the node.
void	setX (size_t i, int dof, double x)
	Set specific current coordinate of the node.
void	addX (size_t i, int dof, double x)
	Add to specific current coordinate of the node.
util::Point &	getU (size_t i)
	Get displacement of the node.
const util::Point &	getU (size_t i) const
	Get displacement of the node.
void	setU (size_t i, const util::Point &u)
	Set displacement of the node.
void	addU (size_t i, const util::Point &u)
	Add to displacement of the node.
void	setU (size_t i, int dof, double u)
	Set displacement of the node.
void	addU (size_t i, int dof, double u)
	Add to displacement of the node.
util::Point &	getV (size_t i)
	Get velocity of the node.
const util::Point &	getV (size_t i) const
	Get velocity of the node.
void	setV (size_t i, const util::Point &v)
	Set velocity of the node.
void	addV (size_t i, const util::Point &v)
	Add to velocity of the node.
void	setV (size_t i, int dof, double v)
	Set velocity of the node.
void	addV (size_t i, int dof, double v)
	Add to velocity of the node.
util::Point &	getF (size_t i)
	Get force of the node.
const util::Point &	getF (size_t i) const
	Get force of the node.
void	setF (size_t i, const util::Point &f)
	Set force of the node.
void	addF (size_t i, const util::Point &f)
	Add to force of the node.
void	setF (size_t i, int dof, double f)
	Set force of the node.
void	addF (size_t i, int dof, double f)
	Add to force of the node.
double &	getVol (size_t i)
	Get volume of the node.
const double &	getVol (size_t i) const
	Get volume of the node.
void	setVol (size_t i, const double &vol)
	Set volume of the node.
void	addVol (size_t i, const double &vol)
	Add to volume of the node.
uint8_t &	getFix (size_t i)
	Get fixity of the node.
const uint8_t &	getFix (size_t i) const
	Get fixity of the node.
void	setFix (size_t i, const unsigned int &dof, const bool &flag)
	Set fixity of the node.
double &	getMx (size_t i)
	Get weighted-volume (mx) of the node.
const double &	getMx (size_t i) const
	Get weighted-volume (mx) of the node.
void	setMx (size_t i, const double &mx)
	Set weighted-volume (mx) of the node.
void	addMx (size_t i, const double &mx)
	Add to weighted-volume (mx) of the node.
double &	getThetax (size_t i)
	Get volumetric deformation (thetax) of the node.
const double &	getThetax (size_t i) const
	Get volumetric deformation (thetax) of the node.
void	setThetax (size_t i, const double &thetax)
	Set volumetric deformation (thetax) of the node.
void	addThetax (size_t i, const double &thetax)
	Add to volumetric deformation (thetax) of the node.

Data Fields
std::string	d_name
	Name of the model for logging purposes (useful if other classes are built on top of this class)
Data Fields inherited from model::ModelData
size_t	d_n
	Current time step.
double	d_time
	Current time.
double	d_currentDt
	Current timestep.
size_t	d_infoN
	Print log step interval.
std::map< std::string, double >	d_dbgData
	Debug data.
std::ofstream	d_ppFile
	File stream to output information.
inp::Input *	d_input_p
	Pointer to Input object.
std::shared_ptr< inp::ModelDeck >	d_modelDeck_p
	Model deck.
std::shared_ptr< inp::RestartDeck >	d_restartDeck_p
	Restart deck.
std::shared_ptr< inp::OutputDeck >	d_outputDeck_p
	Output deck.
std::shared_ptr< inp::ParticleDeck >	d_pDeck_p
	Particle deck.
std::shared_ptr< inp::ContactDeck >	d_cDeck_p
	Contact deck.
bool	d_stop
	flag to stop the simulation midway
double	d_hMax
	Minimum mesh over all particles and walls.
double	d_hMin
	Maximum mesh over all particles and walls.
double	d_maxContactR
	Maximum contact radius between over pairs of particles and walls.
size_t	d_contNeighUpdateInterval
	Neighborlist update interval.
size_t	d_contNeighTimestepCounter
	Contact neighborlist time step counter.
double	d_contNeighSearchRadius
	Neighborlist contact search radius (multiple of d_maxContactR). This variable will be updated during simulation based on maximum velocity.
std::vector< std::shared_ptr< particle::RefParticle > >	d_referenceParticles
	Pointer to reference particle.
std::vector< particle::BaseParticle * >	d_particlesListTypeAll
	List of particles + walls.
std::vector< particle::BaseParticle * >	d_particlesListTypeParticle
	List of particles.
std::vector< particle::BaseParticle * >	d_particlesListTypeWall
	List of walls.
std::vector< inp::MatData >	d_particlesMatDataList
	List of particle material data. Only populated if needed for calculation of stress or other quantities.
std::vector< double >	d_maxVelocityParticlesListTypeAll
	Maximum velocity among all nodes in the particle for each particle.
double	d_maxVelocity
	Maximum velocity among all nodes.
std::vector< std::vector< size_t > >	d_zInfo
	Zone information of particles. For zone 0, d_zInfo[0] = [n1, n2] where n1 is the index at which the particle in this zone starts in d_particlesListTypeAll and n2 is the index + 1, where index is at which the particle in this zone ends.
std::unique_ptr< loading::ParticleULoading >	d_uLoading_p
	Pointer to displacement Loading object.
std::unique_ptr< loading::ParticleFLoading >	d_fLoading_p
	Pointer to force Loading object.
std::unique_ptr< geometry::Fracture >	d_fracture_p
	Fracture state of bonds.
std::unique_ptr< NSearch >	d_nsearch_p
	Pointer to nsearch.
std::vector< util::Point >	d_xRef
	reference positions of the nodes
std::vector< util::Point >	d_x
	Current positions of the nodes.
std::vector< util::Point >	d_u
	Displacement of the nodes.
std::vector< util::Point >	d_v
	Velocity of the nodes.
std::vector< double >	d_vMag
	Magnitude of velocity of the nodes.
std::vector< util::Point >	d_f
	Total force on the nodes.
std::vector< double >	d_vol
	Nodal volumes.
std::vector< size_t >	d_ptId
	Global node to particle id (walls are assigned id after last particle id)
std::vector< std::vector< size_t > >	d_neighC
	Neighbor data for contact forces.
std::vector< std::vector< size_t > >	d_neighPd
	Neighbor data for peridynamic forces.
std::vector< std::vector< float > >	d_neighPdSqdDist
	Square distance neighbor data for peridynamic forces.
std::vector< std::vector< std::vector< size_t > > >	d_neighWallNodes
	Neighbor data for contact between particle and walls.
std::vector< std::vector< std::vector< double > > >	d_neighWallNodesDistance
	Neighbor data (distance) for contact between particle and walls.
std::vector< std::vector< size_t > >	d_neighWallNodesCondensed
	Neighbor data for contact between particle and walls condensed into single vector for each particle.
std::vector< uint8_t >	d_fix
	Vector of fixity mask of each node.
std::vector< uint8_t >	d_forceFixity
	Vector of fixity mask of each node for force.
std::vector< double >	d_thetaX
	Dilation.
std::vector< double >	d_mX
	Weighted volume.
std::vector< size_t >	d_fPdCompNodes
	List of global nodes on which force (peridynamic/internal) is to be computed.
std::vector< size_t >	d_fContCompNodes
	List of global nodes on which force (contact) is to be computed.
std::vector< float >	d_Z
	Damage at nodes.
std::vector< float >	d_e
	Energy of the nodes.
std::vector< float >	d_w
	Work done on each of the nodes.
std::vector< float >	d_phi
	Damage function \( \phi \) at the nodes.
std::vector< float >	d_eF
	Fracture energy of the nodes.
std::vector< float >	d_eFB
	Bond-based fracture energy of the nodes.
std::vector< util::Point >	d_xQuadCur
	Current position of quadrature points.
std::vector< util::SymMatrix3 >	d_strain
	Strain in elements (values at quadrature points)
std::vector< util::SymMatrix3 >	d_stress
	Stress in elements (values at quadrature points)
float	d_te
	Total internal energy.
float	d_tw
	Total work done.
float	d_tk
	Total kinetic energy.
float	d_teF
	Total fracture energy.
float	d_teFB
	Total bond-based fracture energy.

Detailed Description

A class for discrete element particle simulation with peridynamic model

We consider central difference for time integration.

This class acts as a holder of lower rank classes, such as set of particles, neighborlist, etc., and uses the methods and data of the lower rank classes to perform calculation.

Definition at line 32 of file demModel.h.

Constructor & Destructor Documentation

DEMModel()

model::DEMModel::DEMModel ( inp::Input *  deck, std::string  modelName = "DEMModel"  )

explicit

Constructor.

Parameters

deck	The input deck

Definition at line 44 of file demModel.cpp.

  : ModelData(deck),
    d_name(modelName) {
 
  // initialize logger
  util::io::initLogger(d_outputDeck_p->d_debug,
                       d_outputDeck_p->d_path + "log.txt");
}

References model::ModelData::d_outputDeck_p, and util::io::initLogger().

Here is the call graph for this function:

Member Function Documentation

applyInitialCondition()

void model::DEMModel::applyInitialCondition ( )

virtual

Applies initial condition.

Definition at line 1068 of file demModel.cpp.

                                          {
 
  log("Applying initial condition \n", 3);
 
  if (!d_pDeck_p->d_icDeck.d_icActive)
    return;
 
  const auto ic_v = d_pDeck_p->d_icDeck.d_icVec;
  const auto ic_p_list = d_pDeck_p->d_icDeck.d_pList;
 
  // add specified velocity to particle
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index((std::size_t) 0,
                          ic_p_list.size(),
                          (std::size_t) 1,
                          [this, ic_v, ic_p_list](std::size_t i) {
      auto &p = this->d_particlesListTypeAll[ic_p_list[i]];
 
      // velocity
      for (size_t j = 0; j < p->getNumNodes(); j++)
        p->setVLocal(j, ic_v);
    } // loop over particles
  ); // for_each
 
  executor.run(taskflow).get();
}

References util::parallel::getNThreads().

Referenced by twoparticle_demo::Model::integrate().

Here is the call graph for this function:

Here is the caller graph for this function:

checkStop()

void model::DEMModel::checkStop ( )

virtual

Checks if simulation should be stopped due to abnormal state of system.

Definition at line 2151 of file demModel.cpp.

                              {
 
  if (d_outputDeck_p->d_outCriteria == "max_particle_dist" &&
      d_pDeck_p->d_testName == "two_particle") {
 
    // compute max distance between two particles
    // current center position
    const auto &xci = d_particlesListTypeAll[0]->getXCenter();
    const auto &xcj = d_particlesListTypeAll[1]->getXCenter();
 
    // check
    if (util::isGreater(xci.dist(xcj),
                        d_outputDeck_p->d_outCriteriaParams[0])) {
 
      if(d_ppFile.is_open())
        d_ppFile.close();
      exit(1);
    }
  }
  else if (d_outputDeck_p->d_outCriteria == "max_node_dist") {
 
    //    static int msg_printed = 0;
    //    if (msg_printed == 0) {
    //      std::cout << "Check = " << d_outputDeck_p->d_outCriteria
    //              << " is no longer supported. In future, this test will be implemented when function util::methods::maxLength() is defined." << std::endl;
    //      msg_printed = 1;
    //    }
    //exit(EXIT_FAILURE);
    auto max_pt_and_index = util::methods::maxLengthAndMaxLengthIndex(d_x);
    auto max_x = d_x[max_pt_and_index.second];
 
    // check
    if (util::isGreater(max_pt_and_index.first,
                        d_outputDeck_p->d_outCriteriaParams[0])) {
 
      // close open file
      if(d_ppFile.is_open())
        d_ppFile.close();
 
      log(fmt::format("{}: Terminating simulation as one of the failing"
                      " criteria is met. Point ({:.6f}, {:.6f}, {:.6f}) is at "
                      "distance {:.6f} "
                      "more than"
                      " allowed distance {:.6f}\n",
                      d_name, max_x.d_x, max_x.d_y, max_x.d_z, max_x.length(),
                      d_outputDeck_p->d_outCriteriaParams[0]));
      exit(1);
    }
  }
}

References util::isGreater(), and util::methods::maxLengthAndMaxLengthIndex().

Referenced by twoparticle_demo::Model::integrate().

Here is the call graph for this function:

Here is the caller graph for this function:

close()

void model::DEMModel::close ( )

virtual

Closure operations.

Definition at line 104 of file demModel.cpp.

                          {
  if (d_ppFile.is_open())
    d_ppFile.close();
}

Referenced by twoparticle_demo::Model::run().

Here is the caller graph for this function:

computeContactForces()

void model::DEMModel::computeContactForces ( )

virtual

Computes contact forces.

Definition at line 820 of file demModel.cpp.

                                         {
 
  log("    Computing normal contact force \n", 3);
 
  // Description:
  // 1. Normal contact is applied between nodes of particles and walls
  // 2. Normal damping is applied between particle centers
  // 3. Normal damping is applied between nodes of particle and wall pairs
 
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index((std::size_t) 0,
                          d_fContCompNodes.size(),
                          (std::size_t) 1,
                          [this](std::size_t II) {
 
                              auto i = this->d_fContCompNodes[II];
 
                              // local variable to hold force
                              util::Point force_i = util::Point();
                              double scalar_f = 0.;
 
                              const auto &ptIdi = this->getPtId(i);
                              auto &pi = this->getParticleFromAllList(ptIdi);
                              double horizon = pi->d_material_p->getHorizon();
                              double search_r = this->d_maxContactR;
 
                              // particle data
                              double rhoi = pi->getDensity();
 
                              const auto &yi = this->d_x[i]; // current coordinates
                              const auto &ui = this->d_u[i];
                              const auto &vi = this->d_v[i];
                              const auto &voli = this->d_vol[i];
 
                              const std::vector<size_t> &neighs = this->d_neighC[i];
 
                              if (neighs.size() > 0) {
 
                                for (const auto &j_id: neighs) {
 
                                  //auto &j_id = neighs[j];
                                  const auto &yj = this->d_x[j_id]; // current coordinates
                                  double Rji = (yj - yi).length();
                                  auto &ptIdj = this->d_ptId[j_id];
                                  auto &pj = this->getParticleFromAllList(ptIdj);
                                  double rhoj = pj->getDensity();
 
                                  bool both_walls =
                                          (pi->getTypeIndex() == 1 and pj->getTypeIndex() == 1);
 
                                  if (j_id != i) {
                                    if (ptIdj != ptIdi && !both_walls) {
 
                                      // apply particle-particle or particle-wall contact here
                                      const auto &contact =
                                              d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
 
                                      if (util::isLess(Rji, contact.d_contactR)) {
 
                                        auto yji = this->d_x[j_id] - yi;
                                        auto volj = this->d_vol[j_id];
                                        auto vji = this->d_v[j_id] - vi;
 
                                        // resolve velocity vector in normal and tangential components
                                        auto en = yji / Rji;
                                        auto vn_mag = (vji * en);
                                        auto et = vji - vn_mag * en;
                                        if (util::isGreater(et.length(), 0.))
                                          et = et / et.length();
                                        else
                                          et = util::Point();
 
                                        // Formula using bulk modulus and horizon
                                        scalar_f = contact.d_Kn * (Rji - contact.d_contactR) *
                                                   volj; // divided by voli
                                        if (scalar_f > 0.)
                                          scalar_f = 0.;
                                        force_i += scalar_f * en;
 
                                        // compute friction force (since f < 0, |f| = -f)
                                        force_i += contact.d_mu * scalar_f * et;
 
                                        // if particle-wall pair, apply damping contact here <--
                                        // doesnt seem to work
                                        bool node_lvl_damp = false;
                                        // if (pi->getTypeIndex() == 0 and pj->getTypeIndex() == 1)
                                        //   node_lvl_damp = true;
 
                                        if (node_lvl_damp) {
                                          // apply damping at the node level
                                          auto meq = util::equivalentMass(rhoi * voli, rhoj * volj);
                                          auto beta_n =
                                                  contact.d_betan *
                                                  std::sqrt(contact.d_kappa * contact.d_contactR * meq);
 
                                          auto &pii = this->d_particlesListTypeAll[pi->getId()];
                                          vji = this->d_v[j_id] - pii->getVCenter();
                                          vn_mag = (vji * en);
                                          if (vn_mag > 0.)
                                            vn_mag = 0.;
                                          force_i += beta_n * vn_mag * en / voli;
                                        }
                                      } // within contact radius
                                    }   // particle-particle contact
                                  }     // if j_id is not i
                                }       // loop over neighbors
                              }         // contact neighbor
 
                              this->d_f[i] += force_i;
                          }
  ); // for_each
 
  executor.run(taskflow).get();
 
 
  // damping force
  log("    Computing normal damping force \n", 3);
  for (auto &pi : d_particlesListTypeParticle) {
 
    auto pi_id = pi->getId();
 
    double Ri = pi->d_geom_p->boundingRadius();
    double vol_pi = M_PI * Ri * Ri;
    auto pi_xc = pi->getXCenter();
    auto pi_vc = pi->getVCenter();
    auto rhoi = pi->getDensity();
    util::Point force_i = util::Point();
 
    // particle-particle
    for (auto &pj : this->d_particlesListTypeParticle) {
      if (pj->getId() != pi->getId()) {
        auto Rj = pj->d_geom_p->boundingRadius();
        auto xc_ji = pj->getXCenter() - pi_xc;
        auto dist_xcji = xc_ji.length();
 
        const auto &contact = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
 
        if (util::isLess(dist_xcji, Rj + Ri + 1.01 * contact.d_contactR)) {
 
          auto vol_pj = M_PI * Rj * Rj;
          auto rhoj = pj->getDensity();
          // equivalent mass
          auto meq = util::equivalentMass(rhoi * vol_pi, rhoj * vol_pj);
 
          // beta_n
          auto beta_n = contact.d_betan *
                        std::sqrt(contact.d_kappa * contact.d_contactR * meq);
 
          // center-center vector
          auto hat_xc_ji = util::Point();
          if (util::isGreater(dist_xcji, 0.))
            hat_xc_ji = xc_ji / dist_xcji;
          else
            hat_xc_ji = util::Point();
 
          // center-center velocity
          auto vc_ji = pj->getVCenter() - pi_vc;
          auto vc_mag = vc_ji * hat_xc_ji;
          if (vc_mag > 0.)
            vc_mag = 0.;
 
          // force at node of pi
          force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
        } // if within contact distance
      }   // if not same particles
    }     // other particles
 
    // particle-wall
    // Step 1: Create list of wall nodes that are within the Rc distance
    // of at least one of the particle
    // This is done already in updateContactNeighborList()
 
    // step 2 - condensed wall nodes into one vector (has to be done serially
    d_neighWallNodesCondensed[pi->getId()].clear();
    {
      for (size_t j=0; j<d_neighWallNodes[pi_id].size(); j++) {
 
        const auto &j_id = pi->getNodeId(j);
        const auto &yj = this->d_x[j_id];
 
        for (size_t k=0; k<d_neighWallNodes[pi_id][j].size(); k++) {
 
          const auto &k_id = d_neighWallNodes[pi_id][j][k];
          const auto &pk = d_particlesListTypeAll[d_ptId[k_id]];
 
          double Rjk = (this->d_x[k_id] - yj).length();
 
          const auto &contact =
                  d_cDeck_p->getContact(pi->d_zoneId, pk->d_zoneId);
 
          if (util::isLess(Rjk, contact.d_contactR))
            util::methods::addToList(k_id, d_neighWallNodesCondensed[pi_id]);
 
        } // loop over k
      } // loop over j
    } // step 2
 
    // now loop over wall nodes and add force to center of particle
    for (auto &j : d_neighWallNodesCondensed[pi_id]) {
 
      auto &ptIdj = this->d_ptId[j];
      auto &pj = this->d_particlesListTypeAll[ptIdj];
      auto rhoj = pj->getDensity();
      auto volj = this->d_vol[j];
      auto meq = rhoi * vol_pi;
      //auto meq = util::equivalentMass(rhoi * vol_pi, rhoj * volj);
 
      const auto &contact
              = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
 
      // beta_n
      auto beta_n = contact.d_betan *
                    std::sqrt(contact.d_kappa * contact.d_contactR * meq);
 
      // center-node vector
      auto xc_ji = this->d_x[j] - pi_xc;
      auto hat_xc_ji = util::Point();
      if (util::isGreater(xc_ji.length(), 0.))
        hat_xc_ji = xc_ji / xc_ji.length();
 
      // center-node velocity
      auto vc_ji = this->d_v[j] - pi_vc;
      auto vc_mag = vc_ji * hat_xc_ji;
      if (vc_mag > 0.)
        vc_mag = 0.;
 
      // force at node of pi
      force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
    }
 
    // distribute force_i to all nodes of particle pi
    {
      tf::Executor executor(util::parallel::getNThreads());
      tf::Taskflow taskflow;
 
      taskflow.for_each_index((std::size_t) 0, pi->getNumNodes(), (std::size_t) 1,
                              [this, pi, force_i](std::size_t i) {
                                  this->d_f[pi->getNodeId(i)] += force_i;
                              }
      ); // for_each
 
      executor.run(taskflow).get();
    }
  } // loop over particle for damping
}

References util::parallel::getNThreads().

Here is the call graph for this function:

computeExternalDisplacementBC()

void model::DEMModel::computeExternalDisplacementBC ( )

virtual

Applies displacement boundary conditions.

Definition at line 814 of file demModel.cpp.

                                                  {
  log("    Computing external displacement bc \n", 3);
  for (auto &p : d_particlesListTypeAll)
    d_uLoading_p->apply(d_time, p); // applied in parallel
}

Referenced by twoparticle_demo::Model::integrate().

Here is the caller graph for this function:

computeExternalForces()

void model::DEMModel::computeExternalForces ( )

virtual

Computes external/boundary condition forces.

Definition at line 791 of file demModel.cpp.

                                          {
 
  log("    Computing external force \n", 3);
 
  auto gravity = d_pDeck_p->d_gravity;
 
  if (gravity.length() > 1.0E-8) {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;
 
    taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t)1, [this, gravity](std::size_t i) {
          this->d_f[i] += this->getDensity(i) * gravity;
      } // loop over particles
    ); // for_each
 
    executor.run(taskflow).get();
  }
 
  //
  for (auto &p : d_particlesListTypeAll)
    d_fLoading_p->apply(d_time, p); // applied in parallel
}

References util::parallel::getNThreads().

Referenced by peridynamics::Model::computeForces(), and twoparticle_demo::Model::integrate().

Here is the call graph for this function:

Here is the caller graph for this function:

computeForces()

void model::DEMModel::computeForces ( )

virtual

Computes peridynamic forces and contact forces.

Reimplemented in peridynamics::Model.

Definition at line 462 of file demModel.cpp.

                                  {
 
  bool dbg_condition = d_n % d_infoN == 0;
 
  log("  Compute forces \n", 2, dbg_condition, 3);
 
  // reset force
  auto t1 = steady_clock::now();
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index(
    (std::size_t) 0, d_x.size(), (std::size_t) 1,
      [this](std::size_t i) { this->d_f[i] = util::Point(); }
  ); // for_each
 
  executor.run(taskflow).get();
  auto force_reset_time = util::methods::timeDiff(t1, steady_clock::now());
 
  // compute peridynamic forces
  t1 = steady_clock::now();
  computePeridynamicForces();
  auto pd_time = util::methods::timeDiff(t1, steady_clock::now());
  appendKeyData("pd_compute_time", pd_time);
  appendKeyData("avg_peridynamics_force_time", pd_time/d_infoN);
 
  float current_contact_neigh_update_time = 0;
  float contact_time = 0;
  if (d_input_p->isMultiParticle()) {
    // update contact neighborlist
    t1 = steady_clock::now();
    updateContactNeighborlist();
    auto current_contact_neigh_update_time = util::methods::timeDiff(t1,
                                                                     steady_clock::now());
    appendKeyData("contact_neigh_update_time",
                  current_contact_neigh_update_time);
    appendKeyData("avg_contact_neigh_update_time",
                  current_contact_neigh_update_time / d_infoN);
 
    // compute contact forces between particles
    t1 = steady_clock::now();
    computeContactForces();
    auto contact_time = util::methods::timeDiff(t1, steady_clock::now());
    appendKeyData("contact_compute_time", contact_time);
    appendKeyData("avg_contact_force_time", contact_time / d_infoN);
  }
 
  // Compute external forces
  t1 = steady_clock::now();
  computeExternalForces();
  auto extf_time = util::methods::timeDiff(t1, steady_clock::now());
  appendKeyData("extf_compute_time", extf_time);
  appendKeyData("avg_extf_compute_time", extf_time/d_infoN);
 
  // output avg time info
  if (dbg_condition) {
    if (d_input_p->isMultiParticle()) {
      log(fmt::format("    Avg time (ms): \n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n",
                      "tree update", size_t(getKeyData("avg_tree_update_time")),
                      "contact neigh update",
                      size_t(getKeyData("avg_contact_neigh_update_time")),
                      "contact force",
                      size_t(getKeyData("avg_contact_force_time")),
                      "total contact", size_t(getKeyData("avg_tree_update_time")
                                              + getKeyData(
                          "avg_contact_neigh_update_time")
                                              + getKeyData(
                          "avg_contact_force_time")),
                      "peridynamics force",
                      size_t(getKeyData("avg_peridynamics_force_time")),
                      "external force",
                      size_t(getKeyData("avg_extf_compute_time") / d_infoN)),
          2, dbg_condition, 3);
 
      appendKeyData("avg_tree_update_time", 0.);
      appendKeyData("avg_contact_neigh_update_time", 0.);
      appendKeyData("avg_contact_force_time", 0.);
      appendKeyData("avg_peridynamics_force_time", 0.);
      appendKeyData("avg_extf_compute_time", 0.);
    }
    else {
      log(fmt::format("    Avg time (ms): \n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n",
                      "peridynamics force", size_t(getKeyData("avg_peridynamics_force_time")),
                      "external force", size_t(getKeyData("avg_extf_compute_time")/d_infoN)),
          2, dbg_condition, 3);
 
      appendKeyData("avg_peridynamics_force_time", 0.);
      appendKeyData("avg_extf_compute_time", 0.);
    }
  }
 
  log(fmt::format("    {:50s} = {:8d} \n",
                  "Force reset time (ms)",
                  size_t(force_reset_time)
      ),
      2, dbg_condition, 3);
 
  log(fmt::format("    {:50s} = {:8d} \n",
                  "External force time (ms)",
                  size_t(extf_time)
      ),
      2, dbg_condition, 3);
 
  log(fmt::format("    {:50s} = {:8d} \n",
                  "Peridynamics force time (ms)",
                  size_t(pd_time)
      ),
      2, dbg_condition, 3);
 
  if (d_input_p->isMultiParticle()) {
 
    log(fmt::format("    {:50s} = {:8d} \n",
                    "Point cloud update time (ms)",
                    size_t(getKeyData("pt_cloud_update_time"))
        ),
        2, dbg_condition, 3);
 
    log(fmt::format("    {:50s} = {:8d} \n",
                    "Contact neighborlist update time (ms)",
                    size_t(current_contact_neigh_update_time)
        ),
        2, dbg_condition, 3);
 
    log(fmt::format("    {:50s} = {:8d} \n",
                    "Contact force time (ms)",
                    size_t(contact_time)
        ),
        2, dbg_condition, 3);
  }
 
}

References util::parallel::getNThreads(), and util::methods::timeDiff().

Here is the call graph for this function:

computePeridynamicForces()

void model::DEMModel::computePeridynamicForces ( )

virtual

Computes peridynamic forces.

Definition at line 602 of file demModel.cpp.

                                             {
 
  log("    Computing peridynamic force \n", 3);
 
  const auto dim = d_modelDeck_p->d_dim;
  const bool is_state = d_particlesListTypeAll[0]->getMaterial()->isStateActive();
 
  // 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, d_fPdCompNodes.size(), (std::size_t) 1, [this](std::size_t II) {
        auto i = this->d_fPdCompNodes[II];
 
        const auto rho = this->getDensity(i);
        const auto &fix = this->d_fix[i];
        const auto &ptId = this->getPtId(i);
        auto &pi = this->getParticleFromAllList(ptId);
 
        if (pi->d_material_p->isStateActive()) {
 
          const double horizon = pi->getHorizon();
          const double mesh_size = pi->getMeshSize();
          const auto &xi = this->d_xRef[i];
          const auto &ui = this->d_u[i];
 
          // update bond state and compute thetax
          const auto &m = this->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 : this->d_neighPd[i]) {
 
            const auto &xj = this->d_xRef[j];
            const auto &uj = this->d_u[j];
            double rji = (xj - xi).length();
            // double rji = std::sqrt(this->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 = this->d_fracture_p->getBondState(i, k);
            if (!fs && util::isGreater(std::abs(s), sc + 1.0e-10))
              fs = true;
            this->d_fracture_p->setBondState(i, k, fs);
 
            if (!fs) {
 
              // get corrected volume of node j
              auto volj = this->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
 
          this->d_thetaX[i] = 3. * theta / m;
        } // if it is state-based
      } // loop over nodes
    ); // for_each
 
    executor.run(taskflow).get();
  }
 
  // compute the internal forces
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index(
    (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1, [this](std::size_t II) {
      auto i = this->d_fPdCompNodes[II];
 
      // local variable to hold force
      util::Point force_i = util::Point();
      double scalar_f = 0.;
 
      // for damage
      float Zi = 0.;
 
      const auto rhoi = this->getDensity(i);
      const auto &ptIdi = this->getPtId(i);
      auto &pi = this->getParticleFromAllList(ptIdi);
 
      const double horizon = pi->getHorizon();
      const double mesh_size = pi->getMeshSize();
      const auto &xi = this->d_xRef[i];
      const auto &ui = this->d_u[i];
      const auto &mi = this->d_mX[i];
      const auto &thetai = this->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 : this->d_neighPd[i]) {
          auto fs = this->d_fracture_p->getBondState(i, k);
          const auto &xj = this->d_xRef[j];
          const auto &uj = this->d_u[j];
          auto volj = this->d_vol[j];
          double rji = (xj - xi).length();
          double Sji = pi->d_material_p->getS(xj - xi, uj - ui);
 
          if (!fs) {
            const auto &mj = this->d_mX[j];
            const auto &thetaj = this->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;
 
              force_i += scalar_f * pi->d_material_p->getBondForceDirection(
                                        xj - xi, uj - ui);
            } // if state-based
            else {
 
              // Debug
              bool break_bonds = true;
 
              auto ef =
                  pi->d_material_p->getBondEF(rji, Sji, fs, break_bonds);
              this->d_fracture_p->setBondState(i, k, fs);
 
              // compute the contribution of bond force to force at i
              scalar_f = ef.second * volj;
 
              force_i += scalar_f * pi->d_material_p->getBondForceDirection(
                                        xj - xi, uj - ui);
            } // if bond-based
          }   // if bond not broken
          else {
            // add normal contact force
            auto yji = xj + uj - (xi + ui);
            auto Rji = yji.length();
            scalar_f = pi->d_Kn * volj * (Rji - pi->d_Rc) / Rji;
            if (scalar_f > 0.)
              scalar_f = 0.;
            force_i += scalar_f * yji;
          } // 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);
 
          k++;
        } // loop over neighbors
 
      } // peridynamic force
 
      // update force (we remove any force from
      // previous steps and add peridynamics force)
      this->d_f[i] = force_i;
 
      this->d_Z[i] = Zi;
    }
  ); // for_each
 
  executor.run(taskflow).get();
}

References util::parallel::getNThreads().

Referenced by peridynamics::Model::computeForces().

Here is the call graph for this function:

Here is the caller graph for this function:

createGeometryAtSite()

double model::DEMModel::createGeometryAtSite ( const double &  particle_radius, const double &  particle_orient, const util::Point &  site, const std::vector< double > &  rep_geom_params, const std::shared_ptr< util::geometry::GeomObject > &  rep_geom_p, std::shared_ptr< util::geometry::GeomObject > &  p_geom  )

virtual

Creates geometrical object for a particle given particle radius, orientation, and site location.

Definition at line 1377 of file demModel.cpp.

                                                                                          {
  std::vector<double> params;
  for (auto x : rep_geom_params)
    params.push_back(x);
 
  // scale to be used in transforming the geometry of representative particle to get this particle
  double scale = 1.;
 
  if (util::methods::isTagInList(rep_geom_p->d_name,
                                 util::geometry::acceptable_geometries)) {
 
    if (util::methods::isTagInList(rep_geom_p->d_name,
                                   {"circle", "sphere", "hexagon",
                                    "triangle", "square", "cube"})) {
 
      // case - objects requiring four parameters
      // here 'triangle' is a uniform triangle (see constructor of Triangle)
      // 'hexagon' is a hexagon with axis (1, 0, 0)
      size_t num_params = 4;
 
      if (params.size() < num_params)
        params.resize(num_params);
      params[0] = particle_radius;
      for (int dof = 0; dof < 3; dof++)
        params[dof + 1] = site[dof];
 
      scale = particle_radius / rep_geom_params[0];
    }
    else if (rep_geom_p->d_name == "drum2d") {
 
      // case - objects requiring five parameters
      size_t num_params = 5;
 
      if (params.size() < num_params)
        params.resize(num_params);
 
      params[0] = particle_radius; // bigger length along x-direction
      params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0]; // neck length along x-direction
      for (int dof = 0; dof < 3; dof++)
        params[dof + 2] = site[dof];
 
      scale = params[0] / rep_geom_params[0];
    }
    else if (rep_geom_p->d_name == "rectangle") {
 
      // case - objects requiring five parameters
      size_t num_params = 5;
 
      if (params.size() < num_params)
        params.resize(num_params);
 
      params[0] = particle_radius; // length along x-direction
      params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0]; // length along y-direction
      for (int dof = 0; dof < 3; dof++)
        params[dof + 2] = site[dof];
 
      scale = params[0] / rep_geom_params[0];
    }
    else if (rep_geom_p->d_name == "cuboid") {
 
      // case - objects requiring six parameters
 
      if (params.size() < 6)
        params.resize(6);
 
 
      // reference geom object length
      const auto ref_Lx = rep_geom_params[3] - rep_geom_params[0];
      const auto ref_Ly = rep_geom_params[2] - rep_geom_params[1];
      const auto ref_Lz = rep_geom_params[4] - rep_geom_params[2];
 
      const auto Lx = particle_radius; // length is x-direction
      const auto Ly = particle_radius * ref_Ly / ref_Lx; // length in y-direction
      const auto Lz = particle_radius * ref_Lz / ref_Lx; // length in z-direction
 
      params[0] = site[0] - 0.5*Lx;
      params[1] = site[1] - 0.5*Ly;
      params[2] = site[2] - 0.5*Lz;
      params[3] = site[0] + 0.5*Lx;
      params[4] = site[1] + 0.5*Ly;
      params[5] = site[2] + 0.5*Lz;
 
      scale = params[0] / rep_geom_params[0];
    }
  } else if (rep_geom_p->d_name == "ellipse") {
 
    // case - objects requiring six parameters
    size_t num_params = 6;
 
    if (params.size() < num_params)
      params.resize(num_params);
 
    params[0] = particle_radius;
    params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0];
    for (int dof = 0; dof < 3; dof++)
      params[dof + 2] = site[dof];
 
    params[5] = particle_orient;
 
    scale = params[0] / rep_geom_params[0];
  }
  else {
    std::cerr << fmt::format("Error: PeriDEM supports following type "
                             "of geometries for particles = {}\n",
                             util::io::printStr(util::geometry::acceptable_geometries));
    exit(EXIT_FAILURE);
  }
 
  // create geometry now
  std::vector<std::string> vec_geom_type;
  std::vector<std::string> vec_geom_flag;
  util::geometry::createGeomObject(rep_geom_p->d_name, params, vec_geom_type,
                                   vec_geom_flag, p_geom,
                                   d_modelDeck_p->d_dim, false);
 
  return scale;
}

References util::geometry::acceptable_geometries, util::geometry::createGeomObject(), util::methods::isTagInList(), and util::io::printStr().

Here is the call graph for this function:

createParticles()

void model::DEMModel::createParticles ( )

virtual

Creates particles in a given container.

Definition at line 1097 of file demModel.cpp.

                                    {
 
  d_particlesListTypeParticle.resize(0);
  d_particlesListTypeAll.resize(0);
  d_particlesListTypeWall.resize(0);
  d_referenceParticles.clear();
 
  // loop over all particle zones
  for (size_t z = 0; z < d_pDeck_p->d_particleZones.size(); z++) {
 
    // does this particle zone correspond to particle or wall
    bool is_wall = d_pDeck_p->d_particleZones[z].d_isWall;
    std::string particle_type = d_pDeck_p->d_zoneToParticleORWallDeck[z].first;
    if (is_wall and particle_type != "wall") {
      std::cerr << fmt::format("Error: String d_zoneToParticleORWallDeck[z].first for zone z = {} "
                               "should be 'wall'.\n", z);
      exit(EXIT_FAILURE);
    }
    if (!is_wall and particle_type != "particle") {
      std::cerr << fmt::format("Error: String d_zoneToParticleORWallDeck[z].first for zone z = {} "
                               "should be 'particle'.\n", z);
      exit(EXIT_FAILURE);
    }
 
    // get current size of particles data
    auto psize = d_particlesListTypeAll.size();
 
    // get particle zone
    auto &pz = d_pDeck_p->d_particleZones[z];
 
    // get zone id
    auto z_id = pz.d_zone.d_zoneId;
    if (z_id != z) {
      std::cerr << fmt::format("Error: d_zoneId = {} in ParticleZone for "
                               "z = {} should be equal to z = {}.\n",
                               z_id, z, z);
      exit(EXIT_FAILURE);
    }
 
    // read mesh data
    log(d_name + ": Creating mesh for reference particle in zone = " +
        std::to_string(z_id) + "\n");
    std::shared_ptr<fe::Mesh> mesh;
    if (!pz.d_meshDeck.d_createMesh) {
      mesh = std::make_shared<fe::Mesh>(&pz.d_meshDeck);
    }
    else {
      const auto &geomData = pz.d_meshDeck.d_createMeshGeomData;
      if (pz.d_meshDeck.d_createMeshInfo == "uniform"
          and geomData.d_geomName == "rectangle") {
 
        // get the geometrical details
        std::pair<std::vector<double>, std::vector<double>> box;
        std::vector<size_t> nGrid(3, 0);
 
        for (size_t i=0; i<3; i++) {
          box.first.push_back(geomData.d_geomParams[i]);
          box.second.push_back(geomData.d_geomParams[i+3]);
 
          nGrid[i] = size_t((geomData.d_geomParams[i+3] - geomData.d_geomParams[i])/pz.d_meshDeck.d_h);
 
          std::cout << fmt::format("box.first[i] = {}, "
                                   "box.second[i] = {}, "
                                   "nGrid[i] = {}\n",
                                   box.first[i],
                                   box.second[i],
                                   nGrid[i]);
        }
        fe::Mesh temp_mesh;
        fe::createUniformMesh(&temp_mesh,
                              d_modelDeck_p->d_dim,
                              box,
                              nGrid);
        mesh = std::make_shared<fe::Mesh>(temp_mesh);
      }
      else {
        std::cerr << "Error: Currently, we can only support in-built uniform mesh for rectangles.\n";
        exit(EXIT_FAILURE);
      }
    }
 
    // create the reference particle
    log(d_name + ": Creating reference particle in zone = " +
        std::to_string(z_id) + "\n");
 
    // get representative particle for this zone
    auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
    auto rep_geom_params = pz.d_particleGeomData.d_geomParams;
 
    auto ref_p = std::make_shared<particle::RefParticle>(
            d_referenceParticles.size(),
            static_cast<std::shared_ptr<ModelData>>(this),
            rep_geom_p,
            mesh);
 
    d_referenceParticles.emplace_back(ref_p);
 
    // check the particle generation method
    log(d_name + ": Creating particles in zone = " +
                  std::to_string(z_id) + "\n");
 
    if (pz.d_genMethod == "From_File") {
      createParticlesFromFile(z, ref_p);
    }
    else {
      if (pz.d_createParticleUsingParticleZoneGeomObject or !d_input_p->isMultiParticle()) {
        createParticleUsingParticleZoneGeomObject(z, ref_p);
      }
      else {
        std::cerr << "Error: Particle generation method = " << pz.d_genMethod <<
                  " not recognized.\n";
        exit(1);
      }
    }
 
    // get new size of data
    auto psize_new = d_particlesListTypeAll.size();
 
    // store this in zone-info
    d_zInfo.emplace_back(std::vector<size_t>{psize, psize_new, z_id});
  }
}

References fe::createUniformMesh(), and fe::Mesh::d_dim.

Here is the call graph for this function:

createParticlesFromFile()

void model::DEMModel::createParticlesFromFile ( size_t  z, std::shared_ptr< particle::RefParticle >  ref_p  )

virtual

Creates particles in a Hexagonal arrangement.

How the generation works:

1. We translate the reference particle to location given in the file.
2. We scale the reference particle such that the bounding radius of scaled particle is same as radius in the data. So scaling factor is radius/radius_ref_particle
3. We rotate the particle by amount "orient". In case of "loc_rad", we apply random orientation.

Parameters

z	Zone id
ref_p	Shared pointer to reference particle from which we need to create new particles

Definition at line 1266 of file demModel.cpp.

                                                        {
 
  log(d_name + ": Creating particle from file\n", 1);
 
  // get particle zone
  auto &pz = d_pDeck_p->d_particleZones[z];
 
  // get zone id
  auto z_id = pz.d_zone.d_zoneId;
 
  // read file which contains location of centers of particle, zone id, and
  // radius of particle
  std::vector<util::Point> centers;
  std::vector<double> rads;
  std::vector<double> orients;
  if (pz.d_particleFileDataType == "loc_rad") {
    rw::reader::readParticleCsvFile(pz.d_particleFile, d_modelDeck_p->d_dim,
                                    &centers, &rads, z_id);
 
    util::DistributionSample<UniformDistribution> uniform_dist(
        0., 1., d_modelDeck_p->d_seed);
 
    if (d_pDeck_p->d_testName == "two_particle") {
      for (size_t i = 0; i < rads.size(); i++)
        orients.push_back((double(i)) * M_PI);
    } else {
      for (size_t i = 0; i < rads.size(); i++)
        orients.push_back(
            util::transform_to_uniform_dist(0., 2. * M_PI, uniform_dist()));
    }
  }
  else if (pz.d_particleFileDataType == "loc_rad_orient") {
    rw::reader::readParticleWithOrientCsvFile(pz.d_particleFile,
                                              d_modelDeck_p->d_dim, &centers,
                                              &rads, &orients, z_id);
  }
 
  log(fmt::format("zone_id: {}, rads: {}, orients: {}, centers: {} \n", z_id,
                    util::io::printStr(rads), util::io::printStr(orients),
                    util::io::printStr(centers)), 2);
 
  // get representative particle for this zone
  const auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
  auto rep_geom_params = pz.d_particleGeomData.d_geomParams;
 
  // get zone bounding box
  std::pair<util::Point, util::Point> box = rep_geom_p->box();
 
  size_t p_counter = 0;
  size_t p_old_size = d_particlesListTypeAll.size();
  for (const auto &site : centers) {
 
    double particle_radius = rads[p_counter];
    double particle_orient = orients[p_counter];
 
    // create geometrical object
    std::shared_ptr<util::geometry::GeomObject> p_geom;
    auto scale = createGeometryAtSite(particle_radius,
                         particle_orient,
                         site,
                         rep_geom_params,
                         rep_geom_p,
                         p_geom);
 
    // create transform
    auto p_transform = particle::ParticleTransform(
            site, util::Point(0., 0., 1.), particle_orient,
            scale); //particle_radius / ref_p->getParticleRadius());
 
    if (p_transform.d_scale < 1.E-8) {
      std::cerr << "Error: check scale in transform. "
                << " Scale: " << particle_radius / ref_p->getParticleRadius()
                << " p rad: " << particle_radius
                << " ref p rad: " << ref_p->getParticleRadius()
                << p_transform.printStr();
      exit(1);
    }
 
    // finally create dem particle at this site
    //auto particle_id = p_counter + p_old_size;
    auto p = new particle::BaseParticle(
            pz.d_isWall ? "wall" : "particle",
            d_particlesListTypeAll.size(),
            pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
            z_id,
            ref_p->getDimension(),
            pz.d_particleDescription,
            pz.d_isWall,
            pz.d_allDofsConstrained,
            ref_p->getNumNodes(),
            0.,
            static_cast<std::shared_ptr<ModelData>>(this),
            ref_p,
            p_geom,
            p_transform,
            ref_p->getMeshP(),
            pz.d_matDeck,
            true);
 
    // push p to list
    if (pz.d_isWall)
      d_particlesListTypeWall.push_back(p);
    else
      d_particlesListTypeParticle.push_back(p);
 
    d_particlesListTypeAll.push_back(p);
    p_counter++;
  }
}

References util::io::printStr(), rw::reader::readParticleCsvFile(), rw::reader::readParticleWithOrientCsvFile(), and util::transform_to_uniform_dist().

Here is the call graph for this function:

createParticleUsingParticleZoneGeomObject()

void model::DEMModel::createParticleUsingParticleZoneGeomObject ( size_t  z, std::shared_ptr< particle::RefParticle >  ref_p  )

virtual

Creates particles in a given container.

Definition at line 1220 of file demModel.cpp.

                                                  {
 
  log(d_name + ": Creating particle using Particle Zone Geometry Object\n", 1);
 
  // get particle zone
  auto &pz = d_pDeck_p->d_particleZones[z];
 
  // get zone id
  auto z_id = pz.d_zone.d_zoneId;
 
  // ref_p has geometry and mesh which will be used in creating this particle
  // we need to create identity transform
  auto p_transform = particle::ParticleTransform();
 
  // create particle
  auto p = new particle::BaseParticle(
          pz.d_isWall ? "wall" : "particle",
          d_particlesListTypeAll.size(),
          pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
          z_id,
          ref_p->getDimension(),
          pz.d_particleDescription,
          pz.d_isWall,
          pz.d_allDofsConstrained,
          ref_p->getNumNodes(),
          0.,
          static_cast<std::shared_ptr<ModelData>>(this),
          ref_p,
          ref_p->getGeomP(),
          p_transform,
          ref_p->getMeshP(),
          pz.d_matDeck,
          true);
 
  // push p to list
  if (pz.d_isWall)
    d_particlesListTypeWall.push_back(p);
  else
    d_particlesListTypeParticle.push_back(p);
 
  d_particlesListTypeAll.push_back(p);
 
}

init()

void model::DEMModel::init ( )

virtual

Initialize remaining data members.

Definition at line 109 of file demModel.cpp.

                         {
 
  // init time step
  d_n = 0;
  d_time = 0.;
  if (d_outputDeck_p->d_dtTestOut == 0)
    d_outputDeck_p->d_dtTestOut = d_outputDeck_p->d_dtOut / 10;
  d_infoN = d_outputDeck_p->d_dtOut;
 
  // debug/information variables
  {
    appendKeyData("debug_once", -1);
    appendKeyData("update_contact_neigh_search_params_init_call_count", 0);
    appendKeyData("tree_compute_time", 0);
    appendKeyData("contact_compute_time", 0);
    appendKeyData("contact_neigh_update_time", 0);
    appendKeyData("peridynamics_neigh_update_time", 0);
    appendKeyData("pd_compute_time", 0);
    appendKeyData("extf_compute_time", 0);
    appendKeyData("integrate_compute_time", 0);
    appendKeyData("pt_cloud_update_time", 0);
    appendKeyData("avg_tree_update_time", 0);
    appendKeyData("avg_contact_neigh_update_time", 0);
    appendKeyData("avg_contact_force_time", 0);
    appendKeyData("avg_peridynamics_force_time", 0);
    appendKeyData("avg_extf_compute_time", 0);
    appendKeyData("pen_dist", 0);
    appendKeyData("max_y", 0);
    appendKeyData("contact_area_radius", 0);
  }
 
 
  auto t1 = steady_clock::now();
  auto t2 = steady_clock::now();
  log(d_name + ": Initializing objects.\n");
 
  // create particles
  log(d_name + ": Creating particles.\n");
  createParticles();
 
  log(d_name + ": Creating maximum velocity data for particles.\n");
  d_maxVelocityParticlesListTypeAll
          = std::vector<double>(d_particlesListTypeAll.size(), 0.);
  d_maxVelocity = util::methods::max(d_maxVelocityParticlesListTypeAll);
 
  // setup contact
  if (d_input_p->isMultiParticle()) {
    log(d_name + ": Setting up contact.\n");
    setupContact();
  }
 
  // setup element-node connectivity data if needed
  log(d_name + ": Setting up element-node connectivity data for strain/stress.\n");
  setupQuadratureData();
 
  // create search object
  log(d_name + ": Creating neighbor search tree.\n");
 
  // create tree object
  d_nsearch_p = std::make_unique<NSearch>(d_x, d_outputDeck_p->d_debug);
 
  // setup tree
  double set_tree_time = d_nsearch_p->setInputCloud();
  log(fmt::format("{}: Tree setup time (ms) = {}. \n", d_name, set_tree_time));
 
  // create neighborlists
  log(d_name + ": Creating neighborlist for peridynamics.\n");
  t1 = steady_clock::now();
  updatePeridynamicNeighborlist();
  t2 = steady_clock::now();
  appendKeyData("peridynamics_neigh_update_time", util::methods::timeDiff(t1, t2));
 
  if (d_input_p->isMultiParticle()) {
    log(d_name + ": Creating neighborlist for contact.\n");
    d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
    d_contNeighSearchRadius = d_pDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
    t1 = steady_clock::now();
    updateContactNeighborlist();
    t2 = steady_clock::now();
    appendKeyData("contact_neigh_update_time", util::methods::timeDiff(t1, t2));
  }
 
  // create peridynamic bonds
  log(d_name + ": Creating peridynamics bonds.\n");
  d_fracture_p = std::make_unique<geometry::Fracture>(&d_x, &d_neighPd);
 
  // compute quantities in state-based simulations
  log(d_name + ": Compute state-based peridynamic quantities.\n");
  material::computeStateMx(this, true);
 
  // initialize loading class
  log(d_name + ": Initializing displacement loading object.\n");
  d_uLoading_p =
      std::make_unique<loading::ParticleULoading>(d_pDeck_p->d_dispDeck);
  for (auto &p : d_particlesListTypeAll)
    d_uLoading_p->setFixity(p);
 
  log(d_name + ": Initializing force loading object.\n");
  d_fLoading_p =
      std::make_unique<loading::ParticleFLoading>(d_pDeck_p->d_forceDeck);
 
  // if all dofs of particle is fixed, then mark it so that we do not
  // compute force
  // MAYBE NOT as we may be interested in reaction forces
  //  for (auto &p : d_particlesListTypeAll)
  //    p->checkFixityForForce(); // TODO implement
 
  // if this is a two-particle test, we set the force calculation off in
  // first particle
  if (d_pDeck_p->d_testName == "two_particle") {
    d_particlesListTypeAll[0]->d_computeForce = false;
  }
 
  log(fmt::format("{}: Total particles = {}. \n",
                  d_name, d_particlesListTypeAll.size()));
 
  for (const auto &p : d_particlesListTypeAll)
    if (!p->d_computeForce)
      log(fmt::format("{}: Force OFF in Particle i = {}. \n", d_name, p->getId()));
 
  log(d_name + ": Creating list of nodes on which force is to be computed.\n");
  // TODO for now we simply look at particle/wall and check if we compute
  //  force on any of its node. Later, one can have control on individual
  //  nodes of particle/wall and remove from d_fCompNodes if no force is to
  //  be computed on them
  for (size_t i = 0; i < d_x.size(); i++) {
    const auto &ptId = d_ptId[i];
    const auto &pi = getParticleFromAllList(ptId);
    if (pi->d_computeForce) {
      d_fContCompNodes.push_back(i);
      d_fPdCompNodes.push_back(i);
    }
  }
 
  // initialize remaining fields (if any)
  d_Z = std::vector<float>(d_x.size(), 0.);
 
  t2 = steady_clock::now();
  log(fmt::format("{}: Total setup time (ms) = {}. \n",
                  d_name, util::methods::timeDiff(t1, t2)));
 
  // compute complexity information
  size_t free_dofs = 0;
  for (const auto &f : d_fix) {
    for (size_t dof = 0; dof < 3; dof++)
      if (util::methods::isFree(f, dof))
        free_dofs++;
  }
  log(fmt::format("{}: Computational complexity information \n"
                  "  Total number of particles = {}, number of "
                  "particles = {}, number of walls = {}, \n"
                  "  number of dofs = {}, number of free dofs = {}. \n",
                  d_name, d_particlesListTypeAll.size(),
                  d_particlesListTypeParticle.size(),
                  d_particlesListTypeWall.size(),
                  3 * d_x.size(),
                  free_dofs));
}

References material::computeStateMx(), util::methods::isFree(), util::methods::max(), and util::methods::timeDiff().

Referenced by twoparticle_demo::Model::run().

Here is the call graph for this function:

Here is the caller graph for this function:

integrate()

void model::DEMModel::integrate ( )

virtual

Perform time integration.

Reimplemented in twoparticle_demo::Model.

Definition at line 268 of file demModel.cpp.

                              {
 
  // perform output at the beginning
  if (d_n == 0 && d_outputDeck_p->d_performOut) {
    log(fmt::format("{}: Output step = {}, time = {:.6f} \n", d_name, d_n, d_time),
        2);
    output();
  }
 
  // apply initial condition
  if (d_n == 0)
    applyInitialCondition();
 
  // apply loading
  computeExternalDisplacementBC();
  computeForces();
 
  for (size_t i = d_n; i < d_modelDeck_p->d_Nt; i++) {
 
    log(fmt::format("{}: Time step: {}, time: {:8.6f}, steps completed = {}%\n",
                    d_name,
                    i,
                    d_time,
                    float(i) * 100. / d_modelDeck_p->d_Nt),
        2, d_n % d_infoN == 0, 3);
    
    auto t1 = steady_clock::now();
    log("Integrating\n", false, 0, 3);
    integrateStep();
    double integrate_time =
        util::methods::timeDiff(t1, steady_clock::now());
 
    appendKeyData("integrate_compute_time", integrate_time, true);
 
    log(fmt::format("  Integration time (ms) = {}\n", integrate_time), 2, d_n % d_infoN == 0, 3);
 
    if (d_pDeck_p->d_testName == "two_particle") {
 
      // compute location of maximum shear stress and also compute
      // penetration length
      auto msg = ppTwoParticleTest();
      log(msg, 2, d_n % d_infoN == 0, 3);
    } else if (d_pDeck_p->d_testName == "compressive_test") {
      auto msg = ppCompressiveTest();
      log(msg, 2, d_n % d_infoN == 0, 3);
    }
 
    // handle general output
    if ((d_n % d_outputDeck_p->d_dtOut == 0) &&
        (d_n >= d_outputDeck_p->d_dtOut) && d_outputDeck_p->d_performOut) {
      output();
    }
 
    // check for stop
    checkStop();
 
  } // loop over time steps
 
  log(fmt::format(
          "{}: Total compute time information (s) \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n",
          d_name,
          "Time integration", getKeyData("integrate_compute_time") * 1.e-6,
          "Peridynamics force", getKeyData("pd_compute_time") * 1.e-6,
          "Contact force", getKeyData("contact_compute_time") * 1.e-6,
          "Search tree update", getKeyData("tree_compute_time") * 1.e-6,
          "External force", getKeyData("extf_compute_time") * 1.e-6)
          );
}

References util::methods::timeDiff().

Here is the call graph for this function:

integrateCD()

void model::DEMModel::integrateCD ( )

virtual

Perform time integration using central-difference scheme.

Central difference scheme

\[ u_{new} = \Delta t^2 (f_{int} + f_{ext}) / \rho + \Delta t v_{old} + u_{old} \]

\[ v_{new} = \frac{u_{new} - u_{old}}{\Delta t}. \]

Definition at line 349 of file demModel.cpp.

                                {
 
  // update velocity and displacement
  d_currentDt = d_modelDeck_p->d_dt;
  const auto dim = d_modelDeck_p->d_dim;
 
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  // update current position, displacement, and velocity of nodes
  taskflow.for_each_index(
    (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
      [this, dim](std::size_t II) {
        auto i = this->d_fPdCompNodes[II];
 
        const auto rho = this->getDensity(i);
        const auto &fix = this->d_fix[i];
 
        for (int dof = 0; dof < dim; dof++) {
          if (util::methods::isFree(fix, dof)) {
            this->d_v[i][dof] += (this->d_currentDt / rho) * this->d_f[i][dof];
            this->d_u[i][dof] += this->d_currentDt * this->d_v[i][dof];
            this->d_x[i][dof] += this->d_currentDt * this->d_v[i][dof];
          }
        }
 
        this->d_vMag[i] = this->d_v[i].length();
      } // loop over nodes
  ); // for_each
 
  executor.run(taskflow).get();
 
  // advance time
  d_n++;
  d_time += d_currentDt;
 
  // update displacement bc
  computeExternalDisplacementBC();
 
  // compute force
  computeForces();
}

References util::parallel::getNThreads().

Here is the call graph for this function:

integrateStep()

void model::DEMModel::integrateStep ( )

virtual

Performs one time step.

Depending on the time-stepping scheme specified in the input file, this will either call integrateCD or integrateVerlet.

Definition at line 342 of file demModel.cpp.

                                  {
  if (d_modelDeck_p->d_timeDiscretization == "central_difference")
    integrateCD();
  else if (d_modelDeck_p->d_timeDiscretization == "velocity_verlet")
    integrateVerlet();
}

Referenced by twoparticle_demo::Model::integrate().

Here is the caller graph for this function:

integrateVerlet()

void model::DEMModel::integrateVerlet ( )

virtual

Perform time integration using velocity verlet scheme.

Definition at line 392 of file demModel.cpp.

                                    {
 
  // update velocity and displacement
  d_currentDt = d_modelDeck_p->d_dt;
  const auto dim = d_modelDeck_p->d_dim;
 
  // update current position, displacement, and velocity of nodes
  {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;
 
    taskflow.for_each_index(
      (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
        [this, dim](std::size_t II) {
          auto i = this->d_fPdCompNodes[II];
 
          const auto rho = this->getDensity(i);
          const auto &fix = this->d_fix[i];
 
          for (int dof = 0; dof < dim; dof++) {
            if (util::methods::isFree(fix, dof)) {
              this->d_v[i][dof] += 0.5 * (this->d_currentDt / rho) * this->d_f[i][dof];
              this->d_u[i][dof] += this->d_currentDt * this->d_v[i][dof];
              this->d_x[i][dof] += this->d_currentDt * this->d_v[i][dof];
            }
 
            this->d_vMag[i] = this->d_v[i].length();
          }
        } // loop over nodes
    ); // for_each
 
    executor.run(taskflow).get();
  }
 
  // advance time
  d_n++;
  d_time += d_currentDt;
 
  // update displacement bc
  computeExternalDisplacementBC();
 
  // compute force
  computeForces();
 
  // update velocity of nodes
  {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;
 
    taskflow.for_each_index(
      (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
      [this, dim](std::size_t II) {
        auto i = this->d_fPdCompNodes[II];
 
        const auto rho = this->getDensity(i);
        const auto &fix = this->d_fix[i];
        for (int dof = 0; dof < dim; dof++) {
          if (util::methods::isFree(fix, dof)) {
            this->d_v[i][dof] += 0.5 * (this->d_currentDt / rho) * this->d_f[i][dof];
          }
 
          this->d_vMag[i] = this->d_v[i].length();
        }
      } // loop over nodes
    ); // for_each
 
    executor.run(taskflow).get();
  }
}

References util::parallel::getNThreads().

Here is the call graph for this function:

log() [1/2]

void model::DEMModel::log	(	const std::string &	str,
		int	priority = 0,
		bool	check_condition = true,
		int	override_priority = -1,
		bool	screen_out = false
	)

Prints message if any of these two conditions are true.

1. if check_condition == true and dbg_lvl > priority OR
2. dbg_lvl > override_priority

Parameters

str	Message string
priority	Priority of message
check_condition	Specify whether condition for logging/printing message is to be checked
override_priority	Specify new priority
screen_out	Specify whether to print the message also to screen (using std::cout)

Definition at line 61 of file demModel.cpp.

                                           {
  int op = override_priority == -1 ? priority : override_priority;  
  if ((check_condition and d_outputDeck_p->d_debug > priority) or d_outputDeck_p->d_debug > op)
    util::io::log(str, screen_out);
}

References util::io::log().

Here is the call graph for this function:

log() [2/2]

void model::DEMModel::log	(	std::ostringstream &	oss,
		int	priority = 0,
		bool	check_condition = true,
		int	override_priority = -1,
		bool	screen_out = false
	)

Prints message if any of these two conditions are true.

1. if check_condition == true and dbg_lvl > priority OR
2. dbg_lvl > override_priority

Parameters

oss	Stream
priority	Priority of message
check_condition	Specify whether condition for logging/printing message is to be checked
override_priority	Specify new priority
screen_out	Specify whether to print the message also to screen (using std::cout)

Definition at line 53 of file demModel.cpp.

                                           {
  int op = override_priority == -1 ? priority : override_priority;
  //if (d_outputDeck_p->d_debug > priority)
  if ((check_condition and d_outputDeck_p->d_debug > priority) or d_outputDeck_p->d_debug > op)
    util::io::log(oss, screen_out);
}

References util::io::log().

Referenced by peridynamics::Model::computeForces(), twoparticle_demo::Model::integrate(), twoparticle_demo::Model::run(), and twoparticle_demo::Model::twoParticleTestPenetrationDist().

Here is the call graph for this function:

Here is the caller graph for this function:

output()

void model::DEMModel::output ( )

virtual

Output the snapshot of data at current time step.

Definition at line 1962 of file demModel.cpp.

                           {
 
  // write out % completion of simulation at 10% interval
  {
    float p = float(d_n) * 100. / d_modelDeck_p->d_Nt;
    int m = std::max(1, int(d_modelDeck_p->d_Nt / 10));
    if (d_n % m == 0 && int(p) > 0)
      log(fmt::format("{}: Simulation {}% complete\n",
                      d_name, int(p)));
    ;
  }
 
  log(fmt::format("{}: Output step = {}, time = {:.6f} \n",
                  d_name, d_n, d_time),
      2);
 
  if (d_outputDeck_p->d_debug > 0 and getKeyData("debug_once") < 0) {
 
    setKeyData("debug_once", 1);
 
    size_t nt = 1;
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "*******************************************\n";
    oss << tabS << "Debug various input decks\n\n\n";
    oss << d_modelDeck_p->printStr(nt + 1);
    oss << d_pDeck_p->printStr(nt + 1);
    oss << d_cDeck_p->printStr(nt + 1);
    oss << tabS << "\n\n*******************************************\n";
    oss << tabS << "Debug particle data\n\n\n";
    oss << tabS << "Number of particles = " << d_particlesListTypeAll.size() << std::endl;
    oss << tabS << "Number of particle zones = " << d_zInfo.size() << std::endl;
    for (auto zone : d_zInfo) {
      oss << tabS << "zone of d_zInfo: " << util::io::printStr(zone)
          << std::endl;
    }
 
    // wall info
    oss << tabS << "Number of walls = " << d_particlesListTypeWall.size() << std::endl;
    for (auto &d_wall : d_particlesListTypeWall)
      oss << tabS << "Number of nodes in wall in zone " << d_wall->d_zoneId
          << " is " << d_wall->getNumNodes() << std::endl;
 
    oss << tabS << "h_min = " << d_hMin << ", h_max = " << d_hMax << std::endl;
 
    log(oss, 2);
  } // end of debug
 
  size_t dt_out = d_outputDeck_p->d_dtOutCriteria;
  std::string out_filename = d_outputDeck_p->d_path + "output_";
  if (d_outputDeck_p->d_tagPPFile.empty())
    out_filename = out_filename + std::to_string(d_n / dt_out);
  else
    out_filename = out_filename + d_outputDeck_p->d_tagPPFile + "_" + std::to_string(d_n / dt_out);
 
  auto writer = rw::writer::VtkParticleWriter(out_filename);
  if (d_outputDeck_p->d_performFEOut)
    writer.appendMesh(this, d_outputDeck_p->d_outTags);
  else
    writer.appendNodes(this, d_outputDeck_p->d_outTags);
 
  writer.addTimeStep(d_time);
  writer.close();
 
  if (util::methods::isTagInList("Strain_Stress", d_outputDeck_p->d_outTags)) {
 
    // compute current position of quadrature points and strain/stress data
    {
      // if particle mat data is not computed, compute them
      if (d_particlesMatDataList.empty()) {
        for (auto &p: d_particlesListTypeAll) {
          d_particlesMatDataList.push_back(p->getMaterial()->computeMaterialProperties(
                  p->getMeshP()->getDimension()));
        }
      }
 
      for (auto &p: d_particlesListTypeAll) {
 
        const auto particle_mesh_p = p->getMeshP();
 
        fe::getCurrentQuadPoints(particle_mesh_p.get(), d_xRef, d_u, d_xQuadCur,
                                 p->d_globStart,
                                 p->d_globQuadStart,
                                 d_modelDeck_p->d_quadOrder);
 
        auto p_z_id = p->d_zoneId;
        auto isPlaneStrain = d_pDeck_p->d_particleZones[p_z_id].d_matDeck.d_isPlaneStrain;
        fe::getStrainStress(particle_mesh_p.get(), d_xRef, d_u,
                            isPlaneStrain,
                            d_strain, d_stress,
                            p->d_globStart,
                            p->d_globQuadStart,
                            d_particlesMatDataList[p->getId()].d_nu,
                            d_particlesMatDataList[p->getId()].d_lambda,
                            d_particlesMatDataList[p->getId()].d_mu,
                            true,
                            d_modelDeck_p->d_quadOrder);
      } // for loop over particles
    } // compute strain/stress block
 
    out_filename = d_outputDeck_p->d_path + "output_strain_";
    if (d_outputDeck_p->d_tagPPFile.empty())
      out_filename = out_filename + std::to_string(d_n / dt_out);
    else
      out_filename = out_filename + d_outputDeck_p->d_tagPPFile + "_" + std::to_string(d_n / dt_out);
 
    auto writer1 = rw::writer::VtkParticleWriter(out_filename);
    writer1.appendStrainStress(this);
    writer1.addTimeStep(d_time);
    writer1.close();
  }
 
  // output particle locations to csv file
  if (util::methods::isTagInList("Particle_Locations",
                                 d_outputDeck_p->d_outTags)) {
 
    out_filename = d_outputDeck_p->d_path + "particle_locations_";
    if (d_outputDeck_p->d_tagPPFile.empty())
      out_filename = out_filename + std::to_string(d_n / dt_out) + ".csv";
    else
      out_filename = out_filename + d_outputDeck_p->d_tagPPFile
                      + "_" + std::to_string(d_n / dt_out) + ".csv";
 
    std::ofstream oss(out_filename);
    oss << "i, x, y, z, r\n";
    for (const auto &p : d_particlesListTypeAll) {
      auto xc = p->getXCenter();
      oss << p->d_zoneId << ", " << xc.d_x << ", " << xc.d_y << ", " << xc.d_z
          << ", " << p->d_geom_p->boundingRadius() << "\n";
    }
    oss.close();
  }
}

References fe::getCurrentQuadPoints(), fe::getStrainStress(), util::io::getTabS(), util::methods::isTagInList(), and util::io::printStr().

Referenced by twoparticle_demo::Model::integrate().

Here is the call graph for this function:

Here is the caller graph for this function:

ppCompressiveTest()

std::string model::DEMModel::ppCompressiveTest

(

)

Function that handles post-processing for compressive test of particulate media by rigid wall and returns wall penetration and reaction force on wall.

Returns

string Reports wall penetration and reaction force

Definition at line 2202 of file demModel.cpp.

                                           {
  bool continue_dt = false;
  auto check_dt = d_outputDeck_p->d_dtTestOut;
  if ((d_n % check_dt == 0) && (d_n >= check_dt))
    continue_dt = true;
 
  if (!continue_dt)
    return "";
 
  // get wall
  auto w_id = d_pDeck_p->d_particleIdCompressiveTest;
  auto f_dir = d_pDeck_p->d_particleForceDirectionCompressiveTest - 1;
  const auto &wall = d_particlesListTypeAll[w_id];
 
  // find the penetration of the wall from it's original location
  auto dx = wall->getXLocal(0) - wall->getXRefLocal(0);
  double wall_penetration = dx[f_dir];
 
  // get the total reaction force on wall along the direction of loading
  double tot_reaction_force = 0.;
  for (size_t i = 0; i < wall->getNumNodes(); i++) {
    tot_reaction_force += wall->getFLocal(i)[f_dir] * wall->getVolLocal(i);
  }
 
  // open file and write
  bool use_static_file = true;
  if (use_static_file) {
    if (!d_ppFile.is_open()) {
 
      std::string tag_pp_file = d_outputDeck_p->d_tagPPFile.empty() ? "0" : d_outputDeck_p->d_tagPPFile;
      std::string filename = d_outputDeck_p->d_path + "pp_" +
                             d_pDeck_p->d_testName + "_" +
                             tag_pp_file + ".csv";
      d_ppFile.open(filename.c_str(), std::ofstream::out | std::ofstream::app);
 
      d_ppFile << "t, delta, force \n";
    }
 
    d_ppFile << fmt::format("%4.6e, %4.6e, %4.6e\n", d_time, wall_penetration,
            tot_reaction_force);
  }
 
  setKeyData("wall_penetration", wall_penetration);
  setKeyData("tot_reaction_force", tot_reaction_force);
 
  return fmt::format("  Post-processing: wall penetration = {:"
                     ".6f}, "
                     "reaction force = {:5.3e} \n",
                     wall_penetration, tot_reaction_force);
}

ppTwoParticleTest()

std::string model::DEMModel::ppTwoParticleTest

(

)

Function that handles post-processing for two particle collision test and returns maximum vertical displacement of particle.

Returns

string Reports maximum vertical displacement

Definition at line 2096 of file demModel.cpp.

                                           {
 
  bool continue_dt = false;
  auto check_dt = d_outputDeck_p->d_dtTestOut;
  if ((d_n % check_dt == 0) && (d_n >= check_dt))
    continue_dt = true;
 
  if (!continue_dt)
    return "";
 
  // get alias for particles
  const auto &p0 = this->d_particlesListTypeAll[0];
  const auto &p1 = this->d_particlesListTypeAll[1];
 
  // get penetration distance
  const auto &xc0 = p0->getXCenter();
  const auto &xc1 = p1->getXCenter();
  const double &r = p0->d_geom_p->boundingRadius();
 
  const auto &contact = d_cDeck_p->getContact(p0->d_zoneId, p1->d_zoneId);
  double r_e = r + contact.d_contactR;
 
  double pen_dist = xc1.dist(xc0) - r_e - r;
  double contact_area_radius = 0.;
  if (util::isLess(pen_dist, 0.))
    contact_area_radius =
        std::sqrt(std::pow(r_e, 2.) - std::pow(r_e + pen_dist, 2.));
  else if (util::isGreater(pen_dist, 0.)) {
    pen_dist = 0.;
    contact_area_radius = 0.;
  }
 
  // get max distance of second particle (i.e. the y-coord of center + radius)
  double max_dist = xc1.d_y + p1->d_geom_p->boundingRadius();
 
  // compute maximum y coordinate of particle 2
  double max_y_loc = p1->getXLocal(0).d_y;
  double max_y = 0.;
  for (size_t i = 0; i < p1->getNumNodes(); i++)
    if (util::isLess(max_y_loc, p1->getXLocal(i).d_y))
      max_y_loc = p1->getXLocal(i).d_y;
 
  if (util::isLess(max_y, max_y_loc))
    max_y = max_y_loc;
 
  setKeyData("pen_dist", pen_dist);
  setKeyData("contact_area_radius", contact_area_radius);
  setKeyData("max_y", max_y);
  setKeyData("max_dist", max_dist);
  setKeyData("max_y_loc", max_y_loc);
 
 
  return fmt::format("  Post-processing: max y = {:.6f} \n", max_y);
}

References util::isGreater(), and util::isLess().

Here is the call graph for this function:

restart()

void model::DEMModel::restart ( inp::Input *  deck )

virtual

Restarts the simulation from previous state.

Parameters

deck	Input deck

Definition at line 84 of file demModel.cpp.

                                          {
 
  log(d_name + ": Restarting the simulation\n");
 
  // set time step to step specified in restart deck
  d_n = d_restartDeck_p->d_step;
  d_time = double(d_n) * d_modelDeck_p->d_dt;
  log(fmt::format("  Restart step = {}, time = {:.6f} \n", d_n, d_time));
 
  // get backup of reference configuration
  std::vector<util::Point> x_ref(d_x.size(), util::Point());
  for (auto &x : d_x)
    x_ref.push_back(x);
 
  // read displacement and velocity from restart file
  log("  Reading data from restart file = " + d_restartDeck_p->d_file + " \n");
  auto reader = rw::reader::VtkParticleReader(d_restartDeck_p->d_file);
  reader.readNodes(this);
}

run()

void model::DEMModel::run ( inp::Input *  deck )

virtual

Main driver to simulate.

Parameters

deck	Input deck

Reimplemented in twoparticle_demo::Model.

Definition at line 68 of file demModel.cpp.

                                      {
 
  // initialize data
  init();
 
  // check for restart
  if (d_modelDeck_p->d_isRestartActive)
    restart(deck);
 
  // integrate in time
  integrate();
 
  // close
  close();
}

Referenced by main(), and test::testPeriDEM().

Here is the caller graph for this function:

setupContact()

void model::DEMModel::setupContact ( )

virtual

Creates particles in a given container.

Definition at line 1531 of file demModel.cpp.

                                 {
 
  // loop over all particle zones and get minimum value of mesh size
  size_t c = 0;
  for (const auto *p : d_particlesListTypeAll) {
 
    auto h = p->getMeshSize();
    if (c == 0) {
      d_hMin = h;
      d_hMax = h;
      c++;
    }
 
    if (util::isGreater(d_hMin, h))
      d_hMin = h;
    if (util::isGreater(h, d_hMax))
      d_hMax = h;
  }
 
  log(fmt::format("{}: Contact setup\n  hmin = {:.6f}, hmax = {:.6f} \n",
                  d_name, d_hMin, d_hMax), 1);
 
  d_maxContactR = 0.;
  // precompute bulk modulus of all zones
  std::vector<double> bulk_modulus;
  // NOTE - d_data.size() and d_zoneVec.size() are equal
  for (size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
 
    double kappa_i = d_pDeck_p->d_particleZones[i].d_matDeck.d_matData.d_K;
 
    if (kappa_i < 0.) {
      std::cerr << "Error: We need bulk modulus provided in input file.\n";
      std::cerr << d_pDeck_p->d_particleZones[i].printStr();
      exit(1);
    }
 
    bulk_modulus.push_back(kappa_i);
  }
 
  for (size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
    for (size_t j = 0; j < d_cDeck_p->d_data.size(); j++) {
 
      inp::ContactPairDeck *deck = &(d_cDeck_p->d_data[i][j]);
 
      if (deck->d_computeContactR)
        deck->d_contactR *= d_hMin;
 
      if (d_maxContactR < deck->d_contactR)
        d_maxContactR = deck->d_contactR;
 
      // get effective bulk modulus for pair of zones and store it
      deck->d_kappa = util::equivalentMass(bulk_modulus[i], bulk_modulus[j]);
 
      // Kn
      deck->d_Kn *= deck->d_KnFactor;
 
      // Beta n
      double log_e = std::log(deck->d_eps);
      deck->d_betan =
          deck->d_betanFactor *
          (-2. * log_e * std::sqrt(1. / (M_PI * M_PI + log_e * log_e)));
 
      log(fmt::format("  contact_radius = {:.6f}, hmin = {:.6f}, Kn = {:5.3e}, "
                      "Vmax = {:5.3e}, "
                      "betan = {:7.5f}, mu = {:.4f}, kappa = {:5.3e}\n",
                      deck->d_contactR, d_hMin, deck->d_Kn, deck->d_vMax,
                      deck->d_betan, deck->d_mu, deck->d_kappa), 2);
    }
  }
}

References inp::ContactPairDeck::d_betan, inp::ContactPairDeck::d_betanFactor, inp::ContactPairDeck::d_computeContactR, inp::ContactPairDeck::d_contactR, inp::ContactPairDeck::d_eps, inp::ContactPairDeck::d_kappa, inp::ContactPairDeck::d_Kn, inp::ContactPairDeck::d_KnFactor, inp::ContactPairDeck::d_mu, inp::ContactPairDeck::d_vMax, util::equivalentMass(), and util::isGreater().

Here is the call graph for this function:

setupQuadratureData()

void model::DEMModel::setupQuadratureData ( )

virtual

Sets up quadrature data.

Definition at line 1602 of file demModel.cpp.

                                        {
 
  if (util::methods::isTagInList("Strain_Stress", d_outputDeck_p->d_outTags)
      or d_modelDeck_p->d_populateElementNodeConnectivity) {
 
    // read element-node connectivity data if not done
    for (auto &p: d_referenceParticles) {
      auto &particle_mesh_p = p->getMeshP();
      if (!particle_mesh_p->d_encDataPopulated && particle_mesh_p->d_enc.empty()) {
        particle_mesh_p->readElementData(particle_mesh_p->d_filename);
      }
    }
 
    // setup quadrature point and strain/stress data
    // we need to know size of the data
    size_t totalQuadPoints = 0;
    for (auto &p: d_particlesListTypeAll) {
      const auto &particle_mesh_p = p->getMeshP();
 
      // get Quadrature
      fe::BaseElem *elem;
      if (particle_mesh_p->getElementType() == util::vtk_type_line)
        elem = new fe::LineElem(d_modelDeck_p->d_quadOrder);
      else if (particle_mesh_p->getElementType() == util::vtk_type_triangle)
        elem = new fe::TriElem(d_modelDeck_p->d_quadOrder);
      else if (particle_mesh_p->getElementType() == util::vtk_type_quad)
        elem = new fe::QuadElem(d_modelDeck_p->d_quadOrder);
      else if (particle_mesh_p->getElementType() == util::vtk_type_tetra)
        elem = new fe::TetElem(d_modelDeck_p->d_quadOrder);
      else {
        std::cerr << fmt::format("Error: Can not compute strain/stress as the element "
                                 "type = {} is not yet supported in this routine.\n", particle_mesh_p->getElementType());
        exit(EXIT_FAILURE);
      }
 
      p->d_globQuadStart = totalQuadPoints;
      totalQuadPoints += particle_mesh_p->getNumElements() *
                         elem->getNumQuadPoints();
      p->d_globQuadEnd = totalQuadPoints;
 
      std::cout << fmt::format("p->id() = {}, "
                               "p->d_globQuadStart = {}, "
                               "totalQuadPoints = {}, "
                               "p->d_globQuadEnd = {}",
                               p->getId(), p->d_globQuadStart,
                               particle_mesh_p->getNumElements() *
                               elem->getNumQuadPoints(), p->d_globQuadEnd)
                << std::endl;
    }
 
    // resize data
    d_xQuadCur.resize(totalQuadPoints);
    d_strain.resize(totalQuadPoints);
    d_stress.resize(totalQuadPoints);
  } // setting up quadrature data
}

References fe::BaseElem::getNumQuadPoints(), util::methods::isTagInList(), util::vtk_type_line, util::vtk_type_quad, util::vtk_type_tetra, and util::vtk_type_triangle.

Here is the call graph for this function:

updateContactNeighborlist()

void model::DEMModel::updateContactNeighborlist ( )

virtual

Update neighborlist for contact.

Definition at line 1696 of file demModel.cpp.

                                              {
 
  auto update = updateContactNeighborSearchParameters();
 
  if (!update)
    return;
 
  // update contact neighborlist
 
  // update the point cloud (make sure that d_x is updated along with displacement)
  auto pt_cloud_update_time = d_nsearch_p->setInputCloud();
  setKeyData("pt_cloud_update_time", pt_cloud_update_time);
  appendKeyData("tree_compute_time", pt_cloud_update_time);
  appendKeyData("avg_tree_update_time", pt_cloud_update_time/d_infoN);
 
  if (d_neighC.size() != d_x.size())
    d_neighC.resize(d_x.size());
 
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1,
                          [this](std::size_t i) {
 
    const auto &pi = this->d_ptId[i];
    const auto &pi_particle = this->d_particlesListTypeAll[pi];
 
    // search?
    bool perform_search_based_on_particle = true;
    if (pi_particle->d_typeIndex == 1) // wall
      perform_search_based_on_particle = false;
 
    if (pi_particle->d_allDofsConstrained or !pi_particle->d_computeForce)
      perform_search_based_on_particle = false;
 
    if (perform_search_based_on_particle) {
 
      std::vector<size_t> neighs;
      std::vector<double> sqr_dist;
 
      this->d_neighC[i].clear();
 
      auto n = this->d_nsearch_p->radiusSearchExcludeTag(
              this->d_x[i],
              this->d_contNeighSearchRadius,
              neighs,
              sqr_dist,
              this->d_ptId[i],
              this->d_ptId);
 
      if (n > 0) {
        for (auto neigh: neighs) {
          if (neigh != i)
            this->d_neighC[i].push_back(neigh);
        }
      }
    }
}
  ); // for_each
 
  executor.run(taskflow).get();
 
 
  // handle particle-wall neighborlist (based on the d_neighC that we already computed)
  d_neighWallNodes.resize(d_particlesListTypeAll.size());
  d_neighWallNodesDistance.resize(d_particlesListTypeAll.size());
  d_neighWallNodesCondensed.resize(d_particlesListTypeAll.size());
 
  for (auto &pi : d_particlesListTypeParticle) {
 
    d_neighWallNodes[pi->getId()].resize(pi->getNumNodes());
    d_neighWallNodesDistance[pi->getId()].resize(pi->getNumNodes());
 
    // get all wall nodes that are within contact distance to the nodes of this particle
    {
      tf::Executor executor(util::parallel::getNThreads());
      tf::Taskflow taskflow;
 
      taskflow.for_each_index((std::size_t) 0,
                              pi->getNumNodes(),
                              (std::size_t) 1,
                              [this, &pi](std::size_t i) {
 
            auto i_glob = pi->getNodeId(i);
            auto yi = this->d_x[i_glob];
 
            const std::vector<size_t> &neighs = this->d_neighC[i_glob];
 
            this->d_neighWallNodes[pi->getId()][i].clear();
            this->d_neighWallNodesDistance[pi->getId()][i].clear();
 
            for (const auto &j_id: neighs) {
 
              auto &ptIdj = this->d_ptId[j_id];
              auto &pj = this->getParticleFromAllList(
                      ptIdj);
 
              // we are only interested in nodes from wall
              if (pj->getTypeIndex() == 1) {
                  this->d_neighWallNodes[pi->getId()][i].push_back(j_id);
                  //this->d_neighWallNodesDistance[pi->getId()][i].push_back(Rji);
              }
            }
        }
      ); // for_each
 
      executor.run(taskflow).get();
    }
  } // loop over particles
 
}

References util::parallel::getNThreads().

Here is the call graph for this function:

updateContactNeighborSearchParameters()

bool model::DEMModel::updateContactNeighborSearchParameters ( )

virtual

Update contact neighbor search parameters.

Definition at line 1808 of file demModel.cpp.

                                                          {
 
  // initialize parameters
  if (d_contNeighUpdateInterval == 0 and
      util::isLess(d_contNeighSearchRadius, 1.e-16)) {
    d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
    d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
    d_contNeighSearchRadius = d_maxContactR * d_pDeck_p->d_pNeighDeck.d_sFactor;
  }
 
  // at d_n = 0, this function will be called twice because updateContactNeighborlist() will be
  // called twice: one inside init() and second inside computeForces()
  // so to match d_n and d_contNeighTimestepCounter in the initial stage of simulation, we need to handle the special case
  if (d_n == 0) {
    appendKeyData("update_contact_neigh_search_params_init_call_count", 1);
 
    if (int(getKeyData("update_contact_neigh_search_params_init_call_count")) == 1)
      return true;
 
    if (int(getKeyData("update_contact_neigh_search_params_init_call_count")) == 2) {
      d_contNeighTimestepCounter++;
      return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
    }
  }
 
  // handle case of restart
  if (d_modelDeck_p->d_isRestartActive and d_n == d_restartDeck_p->d_step) {
    // assign correct value for restart step
    d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
  }
 
  if (d_contNeighUpdateInterval == 1) {
    // further optimization of parameters is not possible
    d_contNeighSearchRadius = d_maxContactR;
 
    // update counter and return condition for contact search
    d_contNeighTimestepCounter++;
    return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
  }
 
  // check if we should proceed with parameter update
  // param update is done at smaller interval than the search itself to avoid
  // scenarios where particles suddenly move with a high velocity
  size_t update_param_interval =
          d_contNeighUpdateInterval > 5 ? size_t(
                  0.2 * d_contNeighUpdateInterval) : 1;
 
  // check if we ought to update search parameters; if not, return
  if (d_contNeighTimestepCounter > 0 and d_contNeighTimestepCounter % update_param_interval != 0) {
    // update counter and return condition for contact search
    d_contNeighTimestepCounter++;
    return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
  }
 
  // first update the maximum velocity in all particles
  for (auto &pi : d_particlesListTypeAll) {
    auto max_v_node = util::methods::maxIndex(d_vMag,
                                              pi->d_globStart, pi->d_globEnd);
 
    if (max_v_node > pi->d_globEnd or max_v_node < pi->d_globStart) {
      std::cerr << fmt::format("Error: max_v_node = {} for "
                               "particle of id = {} is not in the limit.\n",
                               max_v_node, pi->getId())
                << "Particle info = \n"
                << pi->printStr()
                << "\n\n Magnitude of velocity = "
                << d_vMag[max_v_node] << "\n";
      exit(EXIT_FAILURE);
    }
 
    d_maxVelocityParticlesListTypeAll[pi->getId()]
            = d_vMag[max_v_node];
  }
 
  // find max velocity among all particles
  d_maxVelocity = util::methods::max(d_maxVelocityParticlesListTypeAll);
 
  // now we find the best parameters for contact search
  auto up_interval_old = d_contNeighUpdateInterval;
 
  // TO ensure that in d_neighUpdateInterval time steps, the search radius is above the
  // distance traveled by object with velocity d_maxVelocity
  // also multiply by a safety factor
  double safety_factor = d_pDeck_p->d_pNeighDeck.d_sFactor > 5 ? d_pDeck_p->d_pNeighDeck.d_sFactor : 10;
  auto max_search_r_from_contact_R = d_pDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
  auto max_search_r = d_maxVelocity * d_currentDt
                      * d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval
                      * safety_factor;
 
 
  if (util::isGreater(max_search_r, max_search_r_from_contact_R )) {
 
    d_contNeighUpdateInterval = size_t(d_maxContactR/(d_maxVelocity * d_currentDt));
    if (up_interval_old > d_contNeighUpdateInterval) {
      // issue warning
      log(fmt::format("Warning: Contact search radius based on velocity is greater than "
                      "the max contact radius.\n"
                      "Warning: Adjusting contact neighborlist update interval.\n"
                      "{:>13} = {:4.6e}, time step = {}, "
                      "velocity-based r = {:4.6e}, max contact r = {:4.6e}\n",
                      "Time", d_time, d_n, max_search_r, max_search_r_from_contact_R),
          2, d_n % d_infoN == 0, 3);
    }
 
    d_contNeighSearchRadius = max_search_r_from_contact_R;
    // reset time step counter for contact so that the contact list is updated in the current time step
    // and the update cycle starts from the current time step
    d_contNeighTimestepCounter = 0;
 
    if (d_contNeighUpdateInterval < 1) {
      d_contNeighUpdateInterval = 1;
      d_contNeighSearchRadius = d_maxContactR;
    }
  }
  else {
    // update search radius
    d_contNeighSearchRadius = d_contNeighUpdateInterval < 2 ? d_maxContactR : max_search_r_from_contact_R;
  }
 
  if (up_interval_old > d_contNeighUpdateInterval) {
    log(fmt::format("    Contact neighbor parameters: \n"
                    "      {:48s} = {:d}\n"
                    "      {:48s} = {:d}\n"
                    "      {:48s} = {:d}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n",
                    "time step", d_n,
                    "contact neighbor update interval",
                    d_contNeighUpdateInterval,
                    "contact neighbor update time step counter",
                    d_contNeighTimestepCounter,
                    "search radius", d_contNeighSearchRadius,
                    "max contact radius", d_maxContactR,
                    "search radius factor", d_pDeck_p->d_pNeighDeck.d_sFactor,
                    "max search r from velocity", max_search_r,
                    "max search r from contact r", max_search_r_from_contact_R,
                    "max velocity", d_maxVelocity),
        2, d_n % d_infoN == 0, 3);
  }
 
  // update counter and return condition for contact search
  d_contNeighTimestepCounter++;
  return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
}

References util::isGreater(), util::isLess(), util::methods::max(), and util::methods::maxIndex().

Here is the call graph for this function:

updateGeometryObjectsPostInit()

void model::DEMModel::updateGeometryObjectsPostInit ( )

virtual

Update varioud geometry objects associated with container, particles, and reference particles.

Definition at line 1500 of file demModel.cpp.

                                                  {
 
  for (auto &p: d_particlesListTypeAll) {
    if (p->d_geom_p->d_name == "null" or
        util::methods::isTagInList("copy_from_container", p->d_geom_p->d_tags)) {
      // update geometry of particle based on bounding box
      auto bbox = p->getMeshP()->getBoundingBox();
 
      std::string geom_name = "rectangle";
      if ( p->getMeshP()->getDimension() == 3)
        geom_name = "cuboid";
 
      std::vector<double> geom_params(6, 0.);
      for (size_t i=0; i<3; i++) {
        geom_params[i] = bbox.first[i];
        geom_params[i+3] = bbox.second[i];
      }
 
      std::vector<std::string> vec_type;
      std::vector<std::string> vec_flag;
 
      util::geometry::createGeomObject(geom_name,
                                       geom_params,
                                       vec_type,
                                       vec_flag,
                                       p->d_geom_p,
                                       p->getMeshP()->getDimension());
    }
  }
}

References util::geometry::createGeomObject(), and util::methods::isTagInList().

Here is the call graph for this function:

updateNeighborlistCombine()

void model::DEMModel::updateNeighborlistCombine ( )

virtual

Update neighborlist for contact and peridynamics force.

Definition at line 1957 of file demModel.cpp.

                                              {
  // Not used
  return;
}

updatePeridynamicNeighborlist()

void model::DEMModel::updatePeridynamicNeighborlist ( )

virtual

Update neighborlist for peridynamics force.

Definition at line 1659 of file demModel.cpp.

                                                  {
 
  d_neighPd.resize(d_x.size());
  // d_neighPdSqdDist.resize(d_x.size());
  auto t1 = steady_clock::now();
 
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1, [this](std::size_t i) {
      const auto &pi = this->d_ptId[i];
      double search_r = this->d_particlesListTypeAll[pi]->d_material_p->getHorizon();
 
      std::vector<size_t> neighs;
      std::vector<double> sqr_dist;
      if (this->d_nsearch_p->radiusSearchIncludeTag(this->d_x[i],
                                                    search_r,
                                                    neighs,
                                                    sqr_dist,
                                                    this->d_ptId[i],
                                                    this->d_ptId) > 0) {
        for (std::size_t j = 0; j < neighs.size(); ++j)
          if (neighs[j] != i && this->d_ptId[neighs[j]] == pi) {
            this->d_neighPd[i].push_back(size_t(neighs[j]));
            // this->d_neighPdSqdDist[i].push_back(sqr_dist[j]);
          }
      }
    }
  ); // for_each
 
  executor.run(taskflow).get();
 
  auto t2 = steady_clock::now();
  log(fmt::format("{}: Peridynamics neighbor update time = {}\n",
                  d_name, util::methods::timeDiff(t1, t2)), 2);
}

References util::parallel::getNThreads().

Here is the call graph for this function:

Field Documentation

d_name

std::string model::DEMModel::d_name

Name of the model for logging purposes (useful if other classes are built on top of this class)

Definition at line 42 of file demModel.h.

Referenced by twoparticle_demo::Model::integrate(), and twoparticle_demo::Model::run().

---

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

• /home/user/project/src/model/dem/demModel.h
• /home/user/project/src/model/dem/demModel.cpp