PeriDEM 0.2.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
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::BaseParticlegetParticleFromAllList (size_t i) const
 Get pointer to base particle.
 
particle::BaseParticle *& getParticleFromAllList (size_t i)
 Get pointer to base particle.
 
const particle::BaseParticlegetParticleFromParticleList (size_t i) const
 Get pointer to particle (excluding wall)
 
particle::BaseParticle *& getParticleFromParticleList (size_t i)
 
const particle::BaseParticlegetParticleFromWallList (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::PointgetXRef (size_t i)
 Get reference coordinate of the node.
 
const util::PointgetXRef (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::PointgetX (size_t i)
 Get current coordinate of the node.
 
const util::PointgetX (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::PointgetU (size_t i)
 Get displacement of the node.
 
const util::PointgetU (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::PointgetV (size_t i)
 Get velocity of the node.
 
const util::PointgetV (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::PointgetF (size_t i)
 Get force of the node.
 
const util::PointgetF (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::Inputd_input_p
 Pointer to Input object.
 
std::shared_ptr< inp::ModelDeckd_modelDeck_p
 Model deck.
 
std::shared_ptr< inp::RestartDeckd_restartDeck_p
 Restart deck.
 
std::shared_ptr< inp::OutputDeckd_outputDeck_p
 Output deck.
 
std::shared_ptr< inp::ParticleDeckd_pDeck_p
 Particle deck.
 
std::shared_ptr< inp::ContactDeckd_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::MatDatad_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::ParticleULoadingd_uLoading_p
 Pointer to displacement Loading object.
 
std::unique_ptr< loading::ParticleFLoadingd_fLoading_p
 Pointer to force Loading object.
 
std::unique_ptr< geometry::Fractured_fracture_p
 Fracture state of bonds.
 
std::unique_ptr< NSearchd_nsearch_p
 Pointer to nsearch.
 
std::vector< util::Pointd_xRef
 reference positions of the nodes
 
std::vector< util::Pointd_x
 Current positions of the nodes.
 
std::vector< util::Pointd_u
 Displacement of the nodes.
 
std::vector< util::Pointd_v
 Velocity of the nodes.
 
std::vector< double > d_vMag
 Magnitude of velocity of the nodes.
 
std::vector< util::Pointd_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::Pointd_xQuadCur
 Current position of quadrature points.
 
std::vector< util::SymMatrix3d_strain
 Strain in elements (values at quadrature points)
 
std::vector< util::SymMatrix3d_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
deckThe input deck

Definition at line 44 of file demModel.cpp.

45 : ModelData(deck),
46 d_name(modelName) {
47
48 // initialize logger
50 d_outputDeck_p->d_path + "log.txt");
51}
std::string d_name
Name of the model for logging purposes (useful if other classes are built on top of this class)
Definition demModel.h:42
ModelData(inp::Input *deck)
Constructor.
Definition modelData.h:54
std::shared_ptr< inp::OutputDeck > d_outputDeck_p
Output deck.
Definition modelData.h:558
void initLogger(int debug_level=logger_default_debug_lvl, std::string filename="")
Initializes the logger.
Definition io.cpp:15

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.

1068 {
1069
1070 log("Applying initial condition \n", 3);
1071
1072 if (!d_pDeck_p->d_icDeck.d_icActive)
1073 return;
1074
1075 const auto ic_v = d_pDeck_p->d_icDeck.d_icVec;
1076 const auto ic_p_list = d_pDeck_p->d_icDeck.d_pList;
1077
1078 // add specified velocity to particle
1079 tf::Executor executor(util::parallel::getNThreads());
1080 tf::Taskflow taskflow;
1081
1082 taskflow.for_each_index((std::size_t) 0,
1083 ic_p_list.size(),
1084 (std::size_t) 1,
1085 [this, ic_v, ic_p_list](std::size_t i) {
1086 auto &p = this->d_particlesListTypeAll[ic_p_list[i]];
1087
1088 // velocity
1089 for (size_t j = 0; j < p->getNumNodes(); j++)
1090 p->setVLocal(j, ic_v);
1091 } // loop over particles
1092 ); // for_each
1093
1094 executor.run(taskflow).get();
1095}
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.
Definition demModel.cpp:53
std::shared_ptr< inp::ParticleDeck > d_pDeck_p
Particle deck.
Definition modelData.h:561
unsigned int getNThreads()
Get number of threads to be used by taskflow.

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 2134 of file demModel.cpp.

2134 {
2135
2136 if (d_outputDeck_p->d_outCriteria == "max_particle_dist" &&
2137 d_pDeck_p->d_testName == "two_particle") {
2138
2139 // compute max distance between two particles
2140 // current center position
2141 const auto &xci = d_particlesListTypeAll[0]->getXCenter();
2142 const auto &xcj = d_particlesListTypeAll[1]->getXCenter();
2143
2144 // check
2145 if (util::isGreater(xci.dist(xcj),
2146 d_outputDeck_p->d_outCriteriaParams[0])) {
2147
2148 if(d_ppFile.is_open())
2149 d_ppFile.close();
2150 exit(1);
2151 }
2152 }
2153 else if (d_outputDeck_p->d_outCriteria == "max_node_dist") {
2154
2155 // static int msg_printed = 0;
2156 // if (msg_printed == 0) {
2157 // std::cout << "Check = " << d_outputDeck_p->d_outCriteria
2158 // << " is no longer supported. In future, this test will be implemented when function util::methods::maxLength() is defined." << std::endl;
2159 // msg_printed = 1;
2160 // }
2161 //exit(EXIT_FAILURE);
2162 auto max_pt_and_index = util::methods::maxLengthAndMaxLengthIndex(d_x);
2163 auto max_x = d_x[max_pt_and_index.second];
2164
2165 // check
2166 if (util::isGreater(max_pt_and_index.first,
2167 d_outputDeck_p->d_outCriteriaParams[0])) {
2168
2169 // close open file
2170 if(d_ppFile.is_open())
2171 d_ppFile.close();
2172
2173 log(fmt::format("{}: Terminating simulation as one of the failing"
2174 " criteria is met. Point ({:.6f}, {:.6f}, {:.6f}) is at "
2175 "distance {:.6f} "
2176 "more than"
2177 " allowed distance {:.6f}\n",
2178 d_name, max_x.d_x, max_x.d_y, max_x.d_z, max_x.length(),
2179 d_outputDeck_p->d_outCriteriaParams[0]));
2180 exit(1);
2181 }
2182 }
2183}
std::vector< particle::BaseParticle * > d_particlesListTypeAll
List of particles + walls.
Definition modelData.h:592
std::vector< util::Point > d_x
Current positions of the nodes.
Definition modelData.h:633
std::ofstream d_ppFile
File stream to output information.
Definition modelData.h:546
std::pair< double, size_t > maxLengthAndMaxLengthIndex(const std::vector< util::Point > &data)
Returns the maximum length of point and index from list of points.
Definition methods.h:208
bool isGreater(const double &a, const double &b)
Returns true if a > b.
Definition function.cpp:15

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.

104 {
105 if (d_ppFile.is_open())
106 d_ppFile.close();
107}

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.

820 {
821
822 log(" Computing normal contact force \n", 3);
823
824 // Description:
825 // 1. Normal contact is applied between nodes of particles and walls
826 // 2. Normal damping is applied between particle centers
827 // 3. Normal damping is applied between nodes of particle and wall pairs
828
829 tf::Executor executor(util::parallel::getNThreads());
830 tf::Taskflow taskflow;
831
832 taskflow.for_each_index((std::size_t) 0,
833 d_fContCompNodes.size(),
834 (std::size_t) 1,
835 [this](std::size_t II) {
836
837 auto i = this->d_fContCompNodes[II];
838
839 // local variable to hold force
840 util::Point force_i = util::Point();
841 double scalar_f = 0.;
842
843 const auto &ptIdi = this->getPtId(i);
844 auto &pi = this->getParticleFromAllList(ptIdi);
845 double horizon = pi->d_material_p->getHorizon();
846 double search_r = this->d_maxContactR;
847
848 // particle data
849 double rhoi = pi->getDensity();
850
851 const auto &yi = this->d_x[i]; // current coordinates
852 const auto &ui = this->d_u[i];
853 const auto &vi = this->d_v[i];
854 const auto &voli = this->d_vol[i];
855
856 const std::vector<size_t> &neighs = this->d_neighC[i];
857
858 if (neighs.size() > 0) {
859
860 for (const auto &j_id: neighs) {
861
862 //auto &j_id = neighs[j];
863 const auto &yj = this->d_x[j_id]; // current coordinates
864 double Rji = (yj - yi).length();
865 auto &ptIdj = this->d_ptId[j_id];
866 auto &pj = this->getParticleFromAllList(ptIdj);
867 double rhoj = pj->getDensity();
868
869 bool both_walls =
870 (pi->getTypeIndex() == 1 and pj->getTypeIndex() == 1);
871
872 if (j_id != i) {
873 if (ptIdj != ptIdi && !both_walls) {
874
875 // apply particle-particle or particle-wall contact here
876 const auto &contact =
877 d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
878
879 if (util::isLess(Rji, contact.d_contactR)) {
880
881 auto yji = this->d_x[j_id] - yi;
882 auto volj = this->d_vol[j_id];
883 auto vji = this->d_v[j_id] - vi;
884
885 // resolve velocity vector in normal and tangential components
886 auto en = yji / Rji;
887 auto vn_mag = (vji * en);
888 auto et = vji - vn_mag * en;
889 if (util::isGreater(et.length(), 0.))
890 et = et / et.length();
891 else
892 et = util::Point();
893
894 // Formula using bulk modulus and horizon
895 scalar_f = contact.d_Kn * (Rji - contact.d_contactR) *
896 volj; // divided by voli
897 if (scalar_f > 0.)
898 scalar_f = 0.;
899 force_i += scalar_f * en;
900
901 // compute friction force (since f < 0, |f| = -f)
902 force_i += contact.d_mu * scalar_f * et;
903
904 // if particle-wall pair, apply damping contact here <--
905 // doesnt seem to work
906 bool node_lvl_damp = false;
907 // if (pi->getTypeIndex() == 0 and pj->getTypeIndex() == 1)
908 // node_lvl_damp = true;
909
910 if (node_lvl_damp) {
911 // apply damping at the node level
912 auto meq = util::equivalentMass(rhoi * voli, rhoj * volj);
913 auto beta_n =
914 contact.d_betan *
915 std::sqrt(contact.d_kappa * contact.d_contactR * meq);
916
917 auto &pii = this->d_particlesListTypeAll[pi->getId()];
918 vji = this->d_v[j_id] - pii->getVCenter();
919 vn_mag = (vji * en);
920 if (vn_mag > 0.)
921 vn_mag = 0.;
922 force_i += beta_n * vn_mag * en / voli;
923 }
924 } // within contact radius
925 } // particle-particle contact
926 } // if j_id is not i
927 } // loop over neighbors
928 } // contact neighbor
929
930 this->d_f[i] += force_i;
931 }
932 ); // for_each
933
934 executor.run(taskflow).get();
935
936
937 // damping force
938 log(" Computing normal damping force \n", 3);
939 for (auto &pi : d_particlesListTypeParticle) {
940
941 auto pi_id = pi->getId();
942
943 double Ri = pi->d_geom_p->boundingRadius();
944 double vol_pi = M_PI * Ri * Ri;
945 auto pi_xc = pi->getXCenter();
946 auto pi_vc = pi->getVCenter();
947 auto rhoi = pi->getDensity();
948 util::Point force_i = util::Point();
949
950 // particle-particle
951 for (auto &pj : this->d_particlesListTypeParticle) {
952 if (pj->getId() != pi->getId()) {
953 auto Rj = pj->d_geom_p->boundingRadius();
954 auto xc_ji = pj->getXCenter() - pi_xc;
955 auto dist_xcji = xc_ji.length();
956
957 const auto &contact = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
958
959 if (util::isLess(dist_xcji, Rj + Ri + 1.01 * contact.d_contactR)) {
960
961 auto vol_pj = M_PI * Rj * Rj;
962 auto rhoj = pj->getDensity();
963 // equivalent mass
964 auto meq = util::equivalentMass(rhoi * vol_pi, rhoj * vol_pj);
965
966 // beta_n
967 auto beta_n = contact.d_betan *
968 std::sqrt(contact.d_kappa * contact.d_contactR * meq);
969
970 // center-center vector
971 auto hat_xc_ji = util::Point();
972 if (util::isGreater(dist_xcji, 0.))
973 hat_xc_ji = xc_ji / dist_xcji;
974 else
975 hat_xc_ji = util::Point();
976
977 // center-center velocity
978 auto vc_ji = pj->getVCenter() - pi_vc;
979 auto vc_mag = vc_ji * hat_xc_ji;
980 if (vc_mag > 0.)
981 vc_mag = 0.;
982
983 // force at node of pi
984 force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
985 } // if within contact distance
986 } // if not same particles
987 } // other particles
988
989 // particle-wall
990 // Step 1: Create list of wall nodes that are within the Rc distance
991 // of at least one of the particle
992 // This is done already in updateContactNeighborList()
993
994 // step 2 - condensed wall nodes into one vector (has to be done serially
995 d_neighWallNodesCondensed[pi->getId()].clear();
996 {
997 for (size_t j=0; j<d_neighWallNodes[pi_id].size(); j++) {
998
999 const auto &j_id = pi->getNodeId(j);
1000 const auto &yj = this->d_x[j_id];
1001
1002 for (size_t k=0; k<d_neighWallNodes[pi_id][j].size(); k++) {
1003
1004 const auto &k_id = d_neighWallNodes[pi_id][j][k];
1005 const auto &pk = d_particlesListTypeAll[d_ptId[k_id]];
1006
1007 double Rjk = (this->d_x[k_id] - yj).length();
1008
1009 const auto &contact =
1010 d_cDeck_p->getContact(pi->d_zoneId, pk->d_zoneId);
1011
1012 if (util::isLess(Rjk, contact.d_contactR))
1014
1015 } // loop over k
1016 } // loop over j
1017 } // step 2
1018
1019 // now loop over wall nodes and add force to center of particle
1020 for (auto &j : d_neighWallNodesCondensed[pi_id]) {
1021
1022 auto &ptIdj = this->d_ptId[j];
1023 auto &pj = this->d_particlesListTypeAll[ptIdj];
1024 auto rhoj = pj->getDensity();
1025 auto volj = this->d_vol[j];
1026 auto meq = rhoi * vol_pi;
1027 //auto meq = util::equivalentMass(rhoi * vol_pi, rhoj * volj);
1028
1029 const auto &contact
1030 = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
1031
1032 // beta_n
1033 auto beta_n = contact.d_betan *
1034 std::sqrt(contact.d_kappa * contact.d_contactR * meq);
1035
1036 // center-node vector
1037 auto xc_ji = this->d_x[j] - pi_xc;
1038 auto hat_xc_ji = util::Point();
1039 if (util::isGreater(xc_ji.length(), 0.))
1040 hat_xc_ji = xc_ji / xc_ji.length();
1041
1042 // center-node velocity
1043 auto vc_ji = this->d_v[j] - pi_vc;
1044 auto vc_mag = vc_ji * hat_xc_ji;
1045 if (vc_mag > 0.)
1046 vc_mag = 0.;
1047
1048 // force at node of pi
1049 force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
1050 }
1051
1052 // distribute force_i to all nodes of particle pi
1053 {
1054 tf::Executor executor(util::parallel::getNThreads());
1055 tf::Taskflow taskflow;
1056
1057 taskflow.for_each_index((std::size_t) 0, pi->getNumNodes(), (std::size_t) 1,
1058 [this, pi, force_i](std::size_t i) {
1059 this->d_f[pi->getNodeId(i)] += force_i;
1060 }
1061 ); // for_each
1062
1063 executor.run(taskflow).get();
1064 }
1065 } // loop over particle for damping
1066}
std::vector< size_t > d_fContCompNodes
List of global nodes on which force (contact) is to be computed.
Definition modelData.h:703
std::shared_ptr< inp::ContactDeck > d_cDeck_p
Contact deck.
Definition modelData.h:564
std::vector< util::Point > d_v
Velocity of the nodes.
Definition modelData.h:639
std::vector< std::vector< std::vector< size_t > > > d_neighWallNodes
Neighbor data for contact between particle and walls.
Definition modelData.h:664
std::vector< std::vector< size_t > > d_neighWallNodesCondensed
Neighbor data for contact between particle and walls condensed into single vector for each particle.
Definition modelData.h:670
std::vector< double > d_vol
Nodal volumes.
Definition modelData.h:648
std::vector< size_t > d_ptId
Global node to particle id (walls are assigned id after last particle id)
Definition modelData.h:652
std::vector< particle::BaseParticle * > d_particlesListTypeParticle
List of particles.
Definition modelData.h:595
void addToList(const T &i, std::vector< T > &list)
Add element to the list.
Definition methods.h:289
double equivalentMass(const double &m1, const double &m2)
Compute harmonic mean of m1 and m2.
Definition function.cpp:127
bool isLess(const double &a, const double &b)
Returns true if a < b.
Definition function.cpp:20
A structure to represent 3d vectors.
Definition point.h:30

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.

814 {
815 log(" Computing external displacement bc \n", 3);
816 for (auto &p : d_particlesListTypeAll)
817 d_uLoading_p->apply(d_time, p); // applied in parallel
818}
std::unique_ptr< loading::ParticleULoading > d_uLoading_p
Pointer to displacement Loading object.
Definition modelData.h:618
double d_time
Current time.
Definition modelData.h:534

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.

791 {
792
793 log(" Computing external force \n", 3);
794
795 auto gravity = d_pDeck_p->d_gravity;
796
797 if (gravity.length() > 1.0E-8) {
798 tf::Executor executor(util::parallel::getNThreads());
799 tf::Taskflow taskflow;
800
801 taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t)1, [this, gravity](std::size_t i) {
802 this->d_f[i] += this->getDensity(i) * gravity;
803 } // loop over particles
804 ); // for_each
805
806 executor.run(taskflow).get();
807 }
808
809 //
810 for (auto &p : d_particlesListTypeAll)
811 d_fLoading_p->apply(d_time, p); // applied in parallel
812}
std::unique_ptr< loading::ParticleFLoading > d_fLoading_p
Pointer to force Loading object.
Definition modelData.h:621

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.

462 {
463
464 bool dbg_condition = d_n % d_infoN == 0;
465
466 log(" Compute forces \n", 2, dbg_condition, 3);
467
468 // reset force
469 auto t1 = steady_clock::now();
470 tf::Executor executor(util::parallel::getNThreads());
471 tf::Taskflow taskflow;
472
473 taskflow.for_each_index(
474 (std::size_t) 0, d_x.size(), (std::size_t) 1,
475 [this](std::size_t i) { this->d_f[i] = util::Point(); }
476 ); // for_each
477
478 executor.run(taskflow).get();
479 auto force_reset_time = util::methods::timeDiff(t1, steady_clock::now());
480
481 // compute peridynamic forces
482 t1 = steady_clock::now();
484 auto pd_time = util::methods::timeDiff(t1, steady_clock::now());
485 appendKeyData("pd_compute_time", pd_time);
486 appendKeyData("avg_peridynamics_force_time", pd_time/d_infoN);
487
488 float current_contact_neigh_update_time = 0;
489 float contact_time = 0;
490 if (d_input_p->isMultiParticle()) {
491 // update contact neighborlist
492 t1 = steady_clock::now();
494 auto current_contact_neigh_update_time = util::methods::timeDiff(t1,
495 steady_clock::now());
496 appendKeyData("contact_neigh_update_time",
497 current_contact_neigh_update_time);
498 appendKeyData("avg_contact_neigh_update_time",
499 current_contact_neigh_update_time / d_infoN);
500
501 // compute contact forces between particles
502 t1 = steady_clock::now();
504 auto contact_time = util::methods::timeDiff(t1, steady_clock::now());
505 appendKeyData("contact_compute_time", contact_time);
506 appendKeyData("avg_contact_force_time", contact_time / d_infoN);
507 }
508
509 // Compute external forces
510 t1 = steady_clock::now();
512 auto extf_time = util::methods::timeDiff(t1, steady_clock::now());
513 appendKeyData("extf_compute_time", extf_time);
514 appendKeyData("avg_extf_compute_time", extf_time/d_infoN);
515
516 // output avg time info
517 if (dbg_condition) {
518 if (d_input_p->isMultiParticle()) {
519 log(fmt::format(" Avg time (ms): \n"
520 " {:48s} = {:8d}\n"
521 " {:48s} = {:8d}\n"
522 " {:48s} = {:8d}\n"
523 " {:48s} = {:8d}\n"
524 " {:48s} = {:8d}\n"
525 " {:48s} = {:8d}\n",
526 "tree update", size_t(getKeyData("avg_tree_update_time")),
527 "contact neigh update",
528 size_t(getKeyData("avg_contact_neigh_update_time")),
529 "contact force",
530 size_t(getKeyData("avg_contact_force_time")),
531 "total contact", size_t(getKeyData("avg_tree_update_time")
532 + getKeyData(
533 "avg_contact_neigh_update_time")
534 + getKeyData(
535 "avg_contact_force_time")),
536 "peridynamics force",
537 size_t(getKeyData("avg_peridynamics_force_time")),
538 "external force",
539 size_t(getKeyData("avg_extf_compute_time") / d_infoN)),
540 2, dbg_condition, 3);
541
542 appendKeyData("avg_tree_update_time", 0.);
543 appendKeyData("avg_contact_neigh_update_time", 0.);
544 appendKeyData("avg_contact_force_time", 0.);
545 appendKeyData("avg_peridynamics_force_time", 0.);
546 appendKeyData("avg_extf_compute_time", 0.);
547 }
548 else {
549 log(fmt::format(" Avg time (ms): \n"
550 " {:48s} = {:8d}\n"
551 " {:48s} = {:8d}\n",
552 "peridynamics force", size_t(getKeyData("avg_peridynamics_force_time")),
553 "external force", size_t(getKeyData("avg_extf_compute_time")/d_infoN)),
554 2, dbg_condition, 3);
555
556 appendKeyData("avg_peridynamics_force_time", 0.);
557 appendKeyData("avg_extf_compute_time", 0.);
558 }
559 }
560
561 log(fmt::format(" {:50s} = {:8d} \n",
562 "Force reset time (ms)",
563 size_t(force_reset_time)
564 ),
565 2, dbg_condition, 3);
566
567 log(fmt::format(" {:50s} = {:8d} \n",
568 "External force time (ms)",
569 size_t(extf_time)
570 ),
571 2, dbg_condition, 3);
572
573 log(fmt::format(" {:50s} = {:8d} \n",
574 "Peridynamics force time (ms)",
575 size_t(pd_time)
576 ),
577 2, dbg_condition, 3);
578
579 if (d_input_p->isMultiParticle()) {
580
581 log(fmt::format(" {:50s} = {:8d} \n",
582 "Point cloud update time (ms)",
583 size_t(getKeyData("pt_cloud_update_time"))
584 ),
585 2, dbg_condition, 3);
586
587 log(fmt::format(" {:50s} = {:8d} \n",
588 "Contact neighborlist update time (ms)",
589 size_t(current_contact_neigh_update_time)
590 ),
591 2, dbg_condition, 3);
592
593 log(fmt::format(" {:50s} = {:8d} \n",
594 "Contact force time (ms)",
595 size_t(contact_time)
596 ),
597 2, dbg_condition, 3);
598 }
599
600}
bool isMultiParticle()
Get particle simulation type.
Definition input.cpp:154
virtual void computeExternalForces()
Computes external/boundary condition forces.
Definition demModel.cpp:791
virtual void computePeridynamicForces()
Computes peridynamic forces.
Definition demModel.cpp:602
virtual void updateContactNeighborlist()
Update neighborlist for contact.
virtual void computeContactForces()
Computes contact forces.
Definition demModel.cpp:820
double getKeyData(std::string key, bool issue_err=false)
Get data for a key.
Definition modelData.h:145
size_t d_infoN
Print log step interval.
Definition modelData.h:540
size_t d_n
Current time step.
Definition modelData.h:531
inp::Input * d_input_p
Pointer to Input object.
Definition modelData.h:549
void appendKeyData(std::string key, double data, bool issue_err=false)
Append value to data associated with key.
Definition modelData.h:154
float timeDiff(std::chrono::steady_clock::time_point begin, std::chrono::steady_clock::time_point end, std::string unit="microseconds")
Returns difference between two times.
Definition methods.h:304

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.

602 {
603
604 log(" Computing peridynamic force \n", 3);
605
606 const auto dim = d_modelDeck_p->d_dim;
607 const bool is_state = d_particlesListTypeAll[0]->getMaterial()->isStateActive();
608
609 // compute state-based helper quantities
610 if (is_state) {
611
612 tf::Executor executor(util::parallel::getNThreads());
613 tf::Taskflow taskflow;
614
615 taskflow.for_each_index(
616 (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1, [this](std::size_t II) {
617 auto i = this->d_fPdCompNodes[II];
618
619 const auto rho = this->getDensity(i);
620 const auto &fix = this->d_fix[i];
621 const auto &ptId = this->getPtId(i);
622 auto &pi = this->getParticleFromAllList(ptId);
623
624 if (pi->d_material_p->isStateActive()) {
625
626 const double horizon = pi->getHorizon();
627 const double mesh_size = pi->getMeshSize();
628 const auto &xi = this->d_xRef[i];
629 const auto &ui = this->d_u[i];
630
631 // update bond state and compute thetax
632 const auto &m = this->d_mX[i];
633 double theta = 0.;
634
635 // upper and lower bound for volume correction
636 auto check_up = horizon + 0.5 * mesh_size;
637 auto check_low = horizon - 0.5 * mesh_size;
638
639 size_t k = 0;
640 for (size_t j : this->d_neighPd[i]) {
641
642 const auto &xj = this->d_xRef[j];
643 const auto &uj = this->d_u[j];
644 double rji = (xj - xi).length();
645 // double rji = std::sqrt(this->d_neighPdSqdDist[i][k]);
646 double change_length = (xj - xi + uj - ui).length() - rji;
647
648 // step 1: update the bond state
649 double s = change_length / rji;
650 double sc = pi->d_material_p->getSc(rji);
651
652 // get fracture state, modify, and set
653 auto fs = this->d_fracture_p->getBondState(i, k);
654 if (!fs && util::isGreater(std::abs(s), sc + 1.0e-10))
655 fs = true;
656 this->d_fracture_p->setBondState(i, k, fs);
657
658 if (!fs) {
659
660 // get corrected volume of node j
661 auto volj = this->d_vol[j];
662
663 if (util::isGreater(rji, check_low))
664 volj *= (check_up - rji) / mesh_size;
665
666 theta += rji * change_length * pi->d_material_p->getInfFn(rji) *
667 volj;
668 } // if bond is not broken
669
670 k += 1;
671 } // loop over neighbors
672
673 this->d_thetaX[i] = 3. * theta / m;
674 } // if it is state-based
675 } // loop over nodes
676 ); // for_each
677
678 executor.run(taskflow).get();
679 }
680
681 // compute the internal forces
682 tf::Executor executor(util::parallel::getNThreads());
683 tf::Taskflow taskflow;
684
685 taskflow.for_each_index(
686 (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1, [this](std::size_t II) {
687 auto i = this->d_fPdCompNodes[II];
688
689 // local variable to hold force
690 util::Point force_i = util::Point();
691 double scalar_f = 0.;
692
693 // for damage
694 float Zi = 0.;
695
696 const auto rhoi = this->getDensity(i);
697 const auto &ptIdi = this->getPtId(i);
698 auto &pi = this->getParticleFromAllList(ptIdi);
699
700 const double horizon = pi->getHorizon();
701 const double mesh_size = pi->getMeshSize();
702 const auto &xi = this->d_xRef[i];
703 const auto &ui = this->d_u[i];
704 const auto &mi = this->d_mX[i];
705 const auto &thetai = this->d_thetaX[i];
706
707 // upper and lower bound for volume correction
708 auto check_up = horizon + 0.5 * mesh_size;
709 auto check_low = horizon - 0.5 * mesh_size;
710
711 // loop over neighbors
712 {
713 size_t k = 0;
714 for (size_t j : this->d_neighPd[i]) {
715 auto fs = this->d_fracture_p->getBondState(i, k);
716 const auto &xj = this->d_xRef[j];
717 const auto &uj = this->d_u[j];
718 auto volj = this->d_vol[j];
719 double rji = (xj - xi).length();
720 double Sji = pi->d_material_p->getS(xj - xi, uj - ui);
721
722 if (!fs) {
723 const auto &mj = this->d_mX[j];
724 const auto &thetaj = this->d_thetaX[j];
725
726 // get corrected volume of node j
727 if (util::isGreater(rji, check_low))
728 volj *= (check_up - rji) / mesh_size;
729
730 // handle two cases differently
731 if (pi->d_material_p->isStateActive()) {
732
733 auto ef_i =
734 pi->d_material_p->getBondEF(rji, Sji, fs, mi, thetai);
735 auto ef_j =
736 pi->d_material_p->getBondEF(rji, Sji, fs, mj, thetaj);
737
738 // compute the contribution of bond force to force at i
739 scalar_f = (ef_i.second + ef_j.second) * volj;
740
741 force_i += scalar_f * pi->d_material_p->getBondForceDirection(
742 xj - xi, uj - ui);
743 } // if state-based
744 else {
745
746 // Debug
747 bool break_bonds = true;
748
749 auto ef =
750 pi->d_material_p->getBondEF(rji, Sji, fs, break_bonds);
751 this->d_fracture_p->setBondState(i, k, fs);
752
753 // compute the contribution of bond force to force at i
754 scalar_f = ef.second * volj;
755
756 force_i += scalar_f * pi->d_material_p->getBondForceDirection(
757 xj - xi, uj - ui);
758 } // if bond-based
759 } // if bond not broken
760 else {
761 // add normal contact force
762 auto yji = xj + uj - (xi + ui);
763 auto Rji = yji.length();
764 scalar_f = pi->d_Kn * volj * (Rji - pi->d_Rc) / Rji;
765 if (scalar_f > 0.)
766 scalar_f = 0.;
767 force_i += scalar_f * yji;
768 } // if bond is broken
769
770 // calculate damage
771 auto Sc = pi->d_material_p->getSc(rji);
772 if (util::isGreater(std::abs(Sji / Sc), Zi))
773 Zi = std::abs(Sji / Sc);
774
775 k++;
776 } // loop over neighbors
777
778 } // peridynamic force
779
780 // update force (we remove any force from
781 // previous steps and add peridynamics force)
782 this->d_f[i] = force_i;
783
784 this->d_Z[i] = Zi;
785 }
786 ); // for_each
787
788 executor.run(taskflow).get();
789}
std::vector< float > d_Z
Damage at nodes.
Definition modelData.h:706
std::shared_ptr< inp::ModelDeck > d_modelDeck_p
Model deck.
Definition modelData.h:552
std::vector< size_t > d_fPdCompNodes
List of global nodes on which force (peridynamic/internal) is to be computed.
Definition modelData.h:700

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.

1382 {
1383 std::vector<double> params;
1384 for (auto x : rep_geom_params)
1385 params.push_back(x);
1386
1387 // scale to be used in transforming the geometry of representative particle to get this particle
1388 double scale = 1.;
1389
1390 if (util::methods::isTagInList(rep_geom_p->d_name,
1392
1393 if (util::methods::isTagInList(rep_geom_p->d_name,
1394 {"circle", "sphere", "hexagon",
1395 "triangle", "square", "cube"})) {
1396
1397 // case - objects requiring four parameters
1398 // here 'triangle' is a uniform triangle (see constructor of Triangle)
1399 // 'hexagon' is a hexagon with axis (1, 0, 0)
1400 size_t num_params = 4;
1401
1402 if (params.size() < num_params)
1403 params.resize(num_params);
1404 params[0] = particle_radius;
1405 for (int dof = 0; dof < 3; dof++)
1406 params[dof + 1] = site[dof];
1407
1408 scale = particle_radius / rep_geom_params[0];
1409 }
1410 else if (rep_geom_p->d_name == "drum2d") {
1411
1412 // case - objects requiring five parameters
1413 size_t num_params = 5;
1414
1415 if (params.size() < num_params)
1416 params.resize(num_params);
1417
1418 params[0] = particle_radius; // bigger length along x-direction
1419 params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0]; // neck length along x-direction
1420 for (int dof = 0; dof < 3; dof++)
1421 params[dof + 2] = site[dof];
1422
1423 scale = params[0] / rep_geom_params[0];
1424 }
1425 else if (rep_geom_p->d_name == "rectangle") {
1426
1427 // case - objects requiring five parameters
1428 size_t num_params = 5;
1429
1430 if (params.size() < num_params)
1431 params.resize(num_params);
1432
1433 params[0] = particle_radius; // length along x-direction
1434 params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0]; // length along y-direction
1435 for (int dof = 0; dof < 3; dof++)
1436 params[dof + 2] = site[dof];
1437
1438 scale = params[0] / rep_geom_params[0];
1439 }
1440 else if (rep_geom_p->d_name == "cuboid") {
1441
1442 // case - objects requiring six parameters
1443
1444 if (params.size() < 6)
1445 params.resize(6);
1446
1447
1448 // reference geom object length
1449 const auto ref_Lx = rep_geom_params[3] - rep_geom_params[0];
1450 const auto ref_Ly = rep_geom_params[2] - rep_geom_params[1];
1451 const auto ref_Lz = rep_geom_params[4] - rep_geom_params[2];
1452
1453 const auto Lx = particle_radius; // length is x-direction
1454 const auto Ly = particle_radius * ref_Ly / ref_Lx; // length in y-direction
1455 const auto Lz = particle_radius * ref_Lz / ref_Lx; // length in z-direction
1456
1457 params[0] = site[0] - 0.5*Lx;
1458 params[1] = site[1] - 0.5*Ly;
1459 params[2] = site[2] - 0.5*Lz;
1460 params[3] = site[0] + 0.5*Lx;
1461 params[4] = site[1] + 0.5*Ly;
1462 params[5] = site[2] + 0.5*Lz;
1463
1464 scale = params[0] / rep_geom_params[0];
1465 }
1466 } else {
1467 std::cerr << fmt::format("Error: PeriDEM supports following type "
1468 "of geometries for particles = {}\n",
1470 exit(EXIT_FAILURE);
1471 }
1472
1473 // create geometry now
1474 std::vector<std::string> vec_geom_type;
1475 std::vector<std::string> vec_geom_flag;
1476 util::geometry::createGeomObject(rep_geom_p->d_name, params, vec_geom_type,
1477 vec_geom_flag, p_geom,
1478 d_modelDeck_p->d_dim, false);
1479
1480 return scale;
1481}
const std::vector< std::string > acceptable_geometries
List of acceptable geometries for particles in PeriDEM.
void createGeomObject(const std::string &geom_type, const std::vector< double > &params, const std::vector< std::string > &vec_type, const std::vector< std::string > &vec_flag, std::shared_ptr< util::geometry::GeomObject > &obj, const size_t &dim, bool perform_check=true)
std::string printStr(const T &msg, int nt=print_default_tab)
Returns formatted string for output.
Definition io.h:54
bool isTagInList(const std::string &tag, const std::vector< std::string > &tags)
Returns true if tag is found in the list of tags.
Definition methods.h:279

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.

1097 {
1098
1100 d_particlesListTypeAll.resize(0);
1101 d_particlesListTypeWall.resize(0);
1102 d_referenceParticles.clear();
1103
1104 // loop over all particle zones
1105 for (size_t z = 0; z < d_pDeck_p->d_particleZones.size(); z++) {
1106
1107 // does this particle zone correspond to particle or wall
1108 bool is_wall = d_pDeck_p->d_particleZones[z].d_isWall;
1109 std::string particle_type = d_pDeck_p->d_zoneToParticleORWallDeck[z].first;
1110 if (is_wall and particle_type != "wall") {
1111 std::cerr << fmt::format("Error: String d_zoneToParticleORWallDeck[z].first for zone z = {} "
1112 "should be 'wall'.\n", z);
1113 exit(EXIT_FAILURE);
1114 }
1115 if (!is_wall and particle_type != "particle") {
1116 std::cerr << fmt::format("Error: String d_zoneToParticleORWallDeck[z].first for zone z = {} "
1117 "should be 'particle'.\n", z);
1118 exit(EXIT_FAILURE);
1119 }
1120
1121 // get current size of particles data
1122 auto psize = d_particlesListTypeAll.size();
1123
1124 // get particle zone
1125 auto &pz = d_pDeck_p->d_particleZones[z];
1126
1127 // get zone id
1128 auto z_id = pz.d_zone.d_zoneId;
1129 if (z_id != z) {
1130 std::cerr << fmt::format("Error: d_zoneId = {} in ParticleZone for "
1131 "z = {} should be equal to z = {}.\n",
1132 z_id, z, z);
1133 exit(EXIT_FAILURE);
1134 }
1135
1136 // read mesh data
1137 log(d_name + ": Creating mesh for reference particle in zone = " +
1138 std::to_string(z_id) + "\n");
1139 std::shared_ptr<fe::Mesh> mesh;
1140 if (!pz.d_meshDeck.d_createMesh) {
1141 mesh = std::make_shared<fe::Mesh>(&pz.d_meshDeck);
1142 }
1143 else {
1144 const auto &geomData = pz.d_meshDeck.d_createMeshGeomData;
1145 if (pz.d_meshDeck.d_createMeshInfo == "uniform"
1146 and geomData.d_geomName == "rectangle") {
1147
1148 // get the geometrical details
1149 std::pair<std::vector<double>, std::vector<double>> box;
1150 std::vector<size_t> nGrid(3, 0);
1151
1152 for (size_t i=0; i<3; i++) {
1153 box.first.push_back(geomData.d_geomParams[i]);
1154 box.second.push_back(geomData.d_geomParams[i+3]);
1155
1156 nGrid[i] = size_t((geomData.d_geomParams[i+3] - geomData.d_geomParams[i])/pz.d_meshDeck.d_h);
1157
1158 std::cout << fmt::format("box.first[i] = {}, "
1159 "box.second[i] = {}, "
1160 "nGrid[i] = {}\n",
1161 box.first[i],
1162 box.second[i],
1163 nGrid[i]);
1164 }
1165 fe::Mesh temp_mesh;
1166 fe::createUniformMesh(&temp_mesh,
1167 d_modelDeck_p->d_dim,
1168 box,
1169 nGrid);
1170 mesh = std::make_shared<fe::Mesh>(temp_mesh);
1171 }
1172 else {
1173 std::cerr << "Error: Currently, we can only support in-built uniform mesh for rectangles.\n";
1174 exit(EXIT_FAILURE);
1175 }
1176 }
1177
1178 // create the reference particle
1179 log(d_name + ": Creating reference particle in zone = " +
1180 std::to_string(z_id) + "\n");
1181
1182 // get representative particle for this zone
1183 auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
1184 auto rep_geom_params = pz.d_particleGeomData.d_geomParams;
1185
1186 auto ref_p = std::make_shared<particle::RefParticle>(
1187 d_referenceParticles.size(),
1188 static_cast<std::shared_ptr<ModelData>>(this),
1189 rep_geom_p,
1190 mesh);
1191
1192 d_referenceParticles.emplace_back(ref_p);
1193
1194 // check the particle generation method
1195 log(d_name + ": Creating particles in zone = " +
1196 std::to_string(z_id) + "\n");
1197
1198 if (pz.d_genMethod == "From_File") {
1199 createParticlesFromFile(z, ref_p);
1200 }
1201 else {
1202 if (pz.d_createParticleUsingParticleZoneGeomObject or !d_input_p->isMultiParticle()) {
1204 }
1205 else {
1206 std::cerr << "Error: Particle generation method = " << pz.d_genMethod <<
1207 " not recognized.\n";
1208 exit(1);
1209 }
1210 }
1211
1212 // get new size of data
1213 auto psize_new = d_particlesListTypeAll.size();
1214
1215 // store this in zone-info
1216 d_zInfo.emplace_back(std::vector<size_t>{psize, psize_new, z_id});
1217 }
1218}
A class for mesh data.
Definition mesh.h:51
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.
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 p...
Definition modelData.h:615
std::vector< std::shared_ptr< particle::RefParticle > > d_referenceParticles
Pointer to reference particle.
Definition modelData.h:589
std::vector< particle::BaseParticle * > d_particlesListTypeWall
List of walls.
Definition modelData.h:598
void createUniformMesh(fe::Mesh *mesh_p, size_t dim, std::pair< std::vector< double >, std::vector< double > > box, std::vector< size_t > nGrid)
Creates uniform mesh for rectangle/cuboid domain.
Definition meshUtil.cpp:23

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
zZone id
ref_pShared pointer to reference particle from which we need to create new particles

Definition at line 1266 of file demModel.cpp.

1267 {
1268
1269 log(d_name + ": Creating particle from file\n", 1);
1270
1271 // get particle zone
1272 auto &pz = d_pDeck_p->d_particleZones[z];
1273
1274 // get zone id
1275 auto z_id = pz.d_zone.d_zoneId;
1276
1277 // read file which contains location of centers of particle, zone id, and
1278 // radius of particle
1279 std::vector<util::Point> centers;
1280 std::vector<double> rads;
1281 std::vector<double> orients;
1282 if (pz.d_particleFileDataType == "loc_rad") {
1283 rw::reader::readParticleCsvFile(pz.d_particleFile, d_modelDeck_p->d_dim,
1284 &centers, &rads, z_id);
1285
1287 0., 1., d_modelDeck_p->d_seed);
1288
1289 if (d_pDeck_p->d_testName == "two_particle") {
1290 for (size_t i = 0; i < rads.size(); i++)
1291 orients.push_back((double(i)) * M_PI);
1292 } else {
1293 for (size_t i = 0; i < rads.size(); i++)
1294 orients.push_back(
1295 util::transform_to_uniform_dist(0., 2. * M_PI, uniform_dist()));
1296 }
1297 }
1298 else if (pz.d_particleFileDataType == "loc_rad_orient") {
1300 d_modelDeck_p->d_dim, &centers,
1301 &rads, &orients, z_id);
1302 }
1303
1304 log(fmt::format("zone_id: {}, rads: {}, orients: {}, centers: {} \n", z_id,
1305 util::io::printStr(rads), util::io::printStr(orients),
1306 util::io::printStr(centers)), 2);
1307
1308 // get representative particle for this zone
1309 const auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
1310 auto rep_geom_params = pz.d_particleGeomData.d_geomParams;
1311
1312 // get zone bounding box
1313 std::pair<util::Point, util::Point> box = rep_geom_p->box();
1314
1315 size_t p_counter = 0;
1316 size_t p_old_size = d_particlesListTypeAll.size();
1317 for (const auto &site : centers) {
1318
1319 double particle_radius = rads[p_counter];
1320 double particle_orient = orients[p_counter];
1321
1322 // create geometrical object
1323 std::shared_ptr<util::geometry::GeomObject> p_geom;
1324 auto scale = createGeometryAtSite(particle_radius,
1325 particle_orient,
1326 site,
1327 rep_geom_params,
1328 rep_geom_p,
1329 p_geom);
1330
1331 // create transform
1332 auto p_transform = particle::ParticleTransform(
1333 site, util::Point(0., 0., 1.), particle_orient,
1334 scale); //particle_radius / ref_p->getParticleRadius());
1335
1336 if (p_transform.d_scale < 1.E-8) {
1337 std::cerr << "Error: check scale in transform. "
1338 << " Scale: " << particle_radius / ref_p->getParticleRadius()
1339 << " p rad: " << particle_radius
1340 << " ref p rad: " << ref_p->getParticleRadius()
1341 << p_transform.printStr();
1342 exit(1);
1343 }
1344
1345 // finally create dem particle at this site
1346 //auto particle_id = p_counter + p_old_size;
1347 auto p = new particle::BaseParticle(
1348 pz.d_isWall ? "wall" : "particle",
1350 pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
1351 z_id,
1352 ref_p->getDimension(),
1353 pz.d_particleDescription,
1354 pz.d_isWall,
1355 pz.d_allDofsConstrained,
1356 ref_p->getNumNodes(),
1357 0.,
1358 static_cast<std::shared_ptr<ModelData>>(this),
1359 ref_p,
1360 p_geom,
1361 p_transform,
1362 ref_p->getMeshP(),
1363 pz.d_matDeck,
1364 true);
1365
1366 // push p to list
1367 if (pz.d_isWall)
1368 d_particlesListTypeWall.push_back(p);
1369 else
1370 d_particlesListTypeParticle.push_back(p);
1371
1372 d_particlesListTypeAll.push_back(p);
1373 p_counter++;
1374 }
1375}
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.
A class to store particle geometry, nodal discretization, and methods.
Templated probability distribution.
Definition randomDist.h:90
list p_geom
zone info we have two zones
void readParticleWithOrientCsvFile(const std::string &filename, size_t dim, std::vector< util::Point > *nodes, std::vector< double > *rads, std::vector< double > *orients, const size_t &zone)
Reads particles center location, radius, and zone id. In this case, file also provides initial orient...
Definition reader.cpp:101
void readParticleCsvFile(const std::string &filename, size_t dim, std::vector< util::Point > *nodes, std::vector< double > *rads, std::vector< size_t > *zones)
Reads particles center location, radius, and zone id.
Definition reader.cpp:59
double transform_to_uniform_dist(double min, double max, double sample)
Transform sample from U(0,1) to U(a,b)
Definition randomDist.h:80
A struct that stores transformation parameters and provides method to transform the particle....

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.

1222 {
1223
1224 log(d_name + ": Creating particle using Particle Zone Geometry Object\n", 1);
1225
1226 // get particle zone
1227 auto &pz = d_pDeck_p->d_particleZones[z];
1228
1229 // get zone id
1230 auto z_id = pz.d_zone.d_zoneId;
1231
1232 // ref_p has geometry and mesh which will be used in creating this particle
1233 // we need to create identity transform
1234 auto p_transform = particle::ParticleTransform();
1235
1236 // create particle
1237 auto p = new particle::BaseParticle(
1238 pz.d_isWall ? "wall" : "particle",
1240 pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
1241 z_id,
1242 ref_p->getDimension(),
1243 pz.d_particleDescription,
1244 pz.d_isWall,
1245 pz.d_allDofsConstrained,
1246 ref_p->getNumNodes(),
1247 0.,
1248 static_cast<std::shared_ptr<ModelData>>(this),
1249 ref_p,
1250 ref_p->getGeomP(),
1251 p_transform,
1252 ref_p->getMeshP(),
1253 pz.d_matDeck,
1254 true);
1255
1256 // push p to list
1257 if (pz.d_isWall)
1258 d_particlesListTypeWall.push_back(p);
1259 else
1260 d_particlesListTypeParticle.push_back(p);
1261
1262 d_particlesListTypeAll.push_back(p);
1263
1264}

◆ init()

void model::DEMModel::init ( )
virtual

Initialize remaining data members.

Definition at line 109 of file demModel.cpp.

109 {
110
111 // init time step
112 d_n = 0;
113 d_time = 0.;
114 if (d_outputDeck_p->d_dtTestOut == 0)
115 d_outputDeck_p->d_dtTestOut = d_outputDeck_p->d_dtOut / 10;
116 d_infoN = d_outputDeck_p->d_dtOut;
117
118 // debug/information variables
119 {
120 appendKeyData("debug_once", -1);
121 appendKeyData("update_contact_neigh_search_params_init_call_count", 0);
122 appendKeyData("tree_compute_time", 0);
123 appendKeyData("contact_compute_time", 0);
124 appendKeyData("contact_neigh_update_time", 0);
125 appendKeyData("peridynamics_neigh_update_time", 0);
126 appendKeyData("pd_compute_time", 0);
127 appendKeyData("extf_compute_time", 0);
128 appendKeyData("integrate_compute_time", 0);
129 appendKeyData("pt_cloud_update_time", 0);
130 appendKeyData("avg_tree_update_time", 0);
131 appendKeyData("avg_contact_neigh_update_time", 0);
132 appendKeyData("avg_contact_force_time", 0);
133 appendKeyData("avg_peridynamics_force_time", 0);
134 appendKeyData("avg_extf_compute_time", 0);
135 appendKeyData("pen_dist", 0);
136 appendKeyData("max_y", 0);
137 appendKeyData("contact_area_radius", 0);
138 }
139
140
141 auto t1 = steady_clock::now();
142 auto t2 = steady_clock::now();
143 log(d_name + ": Initializing objects.\n");
144
145 // create particles
146 log(d_name + ": Creating particles.\n");
148
149 log(d_name + ": Creating maximum velocity data for particles.\n");
151 = std::vector<double>(d_particlesListTypeAll.size(), 0.);
153
154 // setup contact
155 if (d_input_p->isMultiParticle()) {
156 log(d_name + ": Setting up contact.\n");
157 setupContact();
158 }
159
160 // setup element-node connectivity data if needed
161 log(d_name + ": Setting up element-node connectivity data for strain/stress.\n");
163
164 // create search object
165 log(d_name + ": Creating neighbor search tree.\n");
166
167 // create tree object
168 d_nsearch_p = std::make_unique<NSearch>(d_x, d_outputDeck_p->d_debug);
169
170 // setup tree
171 double set_tree_time = d_nsearch_p->setInputCloud();
172 log(fmt::format("{}: Tree setup time (ms) = {}. \n", d_name, set_tree_time));
173
174 // create neighborlists
175 log(d_name + ": Creating neighborlist for peridynamics.\n");
176 t1 = steady_clock::now();
178 t2 = steady_clock::now();
179 appendKeyData("peridynamics_neigh_update_time", util::methods::timeDiff(t1, t2));
180
181 if (d_input_p->isMultiParticle()) {
182 log(d_name + ": Creating neighborlist for contact.\n");
183 d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
184 d_contNeighSearchRadius = d_pDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
185 t1 = steady_clock::now();
187 t2 = steady_clock::now();
188 appendKeyData("contact_neigh_update_time", util::methods::timeDiff(t1, t2));
189 }
190
191 // create peridynamic bonds
192 log(d_name + ": Creating peridynamics bonds.\n");
193 d_fracture_p = std::make_unique<geometry::Fracture>(&d_x, &d_neighPd);
194
195 // compute quantities in state-based simulations
196 log(d_name + ": Compute state-based peridynamic quantities.\n");
197 material::computeStateMx(this, true);
198
199 // initialize loading class
200 log(d_name + ": Initializing displacement loading object.\n");
202 std::make_unique<loading::ParticleULoading>(d_pDeck_p->d_dispDeck);
203 for (auto &p : d_particlesListTypeAll)
204 d_uLoading_p->setFixity(p);
205
206 log(d_name + ": Initializing force loading object.\n");
208 std::make_unique<loading::ParticleFLoading>(d_pDeck_p->d_forceDeck);
209
210 // if all dofs of particle is fixed, then mark it so that we do not
211 // compute force
212 // MAYBE NOT as we may be interested in reaction forces
213 // for (auto &p : d_particlesListTypeAll)
214 // p->checkFixityForForce(); // TODO implement
215
216 // if this is a two-particle test, we set the force calculation off in
217 // first particle
218 if (d_pDeck_p->d_testName == "two_particle") {
219 d_particlesListTypeAll[0]->d_computeForce = false;
220 }
221
222 log(fmt::format("{}: Total particles = {}. \n",
224
225 for (const auto &p : d_particlesListTypeAll)
226 if (!p->d_computeForce)
227 log(fmt::format("{}: Force OFF in Particle i = {}. \n", d_name, p->getId()));
228
229 log(d_name + ": Creating list of nodes on which force is to be computed.\n");
230 // TODO for now we simply look at particle/wall and check if we compute
231 // force on any of its node. Later, one can have control on individual
232 // nodes of particle/wall and remove from d_fCompNodes if no force is to
233 // be computed on them
234 for (size_t i = 0; i < d_x.size(); i++) {
235 const auto &ptId = d_ptId[i];
236 const auto &pi = getParticleFromAllList(ptId);
237 if (pi->d_computeForce) {
238 d_fContCompNodes.push_back(i);
239 d_fPdCompNodes.push_back(i);
240 }
241 }
242
243 // initialize remaining fields (if any)
244 d_Z = std::vector<float>(d_x.size(), 0.);
245
246 t2 = steady_clock::now();
247 log(fmt::format("{}: Total setup time (ms) = {}. \n",
249
250 // compute complexity information
251 size_t free_dofs = 0;
252 for (const auto &f : d_fix) {
253 for (size_t dof = 0; dof < 3; dof++)
254 if (util::methods::isFree(f, dof))
255 free_dofs++;
256 }
257 log(fmt::format("{}: Computational complexity information \n"
258 " Total number of particles = {}, number of "
259 "particles = {}, number of walls = {}, \n"
260 " number of dofs = {}, number of free dofs = {}. \n",
264 3 * d_x.size(),
265 free_dofs));
266}
virtual void setupContact()
Creates particles in a given container.
virtual void updatePeridynamicNeighborlist()
Update neighborlist for peridynamics force.
virtual void createParticles()
Creates particles in a given container.
virtual void setupQuadratureData()
Sets up quadrature data.
std::vector< uint8_t > d_fix
Vector of fixity mask of each node.
Definition modelData.h:680
std::vector< double > d_maxVelocityParticlesListTypeAll
Maximum velocity among all nodes in the particle for each particle.
Definition modelData.h:607
std::unique_ptr< NSearch > d_nsearch_p
Pointer to nsearch.
Definition modelData.h:627
double d_contNeighSearchRadius
Neighborlist contact search radius (multiple of d_maxContactR). This variable will be updated during ...
Definition modelData.h:586
size_t d_contNeighUpdateInterval
Neighborlist update interval.
Definition modelData.h:580
double d_maxContactR
Maximum contact radius between over pairs of particles and walls.
Definition modelData.h:577
double d_maxVelocity
Maximum velocity among all nodes.
Definition modelData.h:610
const particle::BaseParticle * getParticleFromAllList(size_t i) const
Get pointer to base particle.
Definition modelData.h:76
std::unique_ptr< geometry::Fracture > d_fracture_p
Fracture state of bonds.
Definition modelData.h:624
std::vector< std::vector< size_t > > d_neighPd
Neighbor data for peridynamic forces.
Definition modelData.h:658
void computeStateMx(model::ModelData *model, bool compute_in_parallel=false)
Computes the moment term in state-based peridynamic formulation.
T max(const std::vector< T > &data)
Returns the maximum from list of data.
Definition methods.h:74
bool isFree(const int &i, const unsigned int &dof)
Returns true if degree of freedom is free.
Definition methods.h:249

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.

268 {
269
270 // perform output at the beginning
271 if (d_n == 0 && d_outputDeck_p->d_performOut) {
272 log(fmt::format("{}: Output step = {}, time = {:.6f} \n", d_name, d_n, d_time),
273 2);
274 output();
275 }
276
277 // apply initial condition
278 if (d_n == 0)
280
281 // apply loading
284
285 for (size_t i = d_n; i < d_modelDeck_p->d_Nt; i++) {
286
287 log(fmt::format("{}: Time step: {}, time: {:8.6f}, steps completed = {}%\n",
288 d_name,
289 i,
290 d_time,
291 float(i) * 100. / d_modelDeck_p->d_Nt),
292 2, d_n % d_infoN == 0, 3);
293
294 auto t1 = steady_clock::now();
295 log("Integrating\n", false, 0, 3);
297 double integrate_time =
298 util::methods::timeDiff(t1, steady_clock::now());
299
300 appendKeyData("integrate_compute_time", integrate_time, true);
301
302 log(fmt::format(" Integration time (ms) = {}\n", integrate_time), 2, d_n % d_infoN == 0, 3);
303
304 if (d_pDeck_p->d_testName == "two_particle") {
305
306 // compute location of maximum shear stress and also compute
307 // penetration length
308 auto msg = ppTwoParticleTest();
309 log(msg, 2, d_n % d_infoN == 0, 3);
310 } else if (d_pDeck_p->d_testName == "compressive_test") {
311 auto msg = ppCompressiveTest();
312 log(msg, 2, d_n % d_infoN == 0, 3);
313 }
314
315 // handle general output
316 if ((d_n % d_outputDeck_p->d_dtOut == 0) &&
317 (d_n >= d_outputDeck_p->d_dtOut) && d_outputDeck_p->d_performOut) {
318 output();
319 }
320
321 // check for stop
322 checkStop();
323
324 } // loop over time steps
325
326 log(fmt::format(
327 "{}: Total compute time information (s) \n"
328 " {:22s} = {:8.2f} \n"
329 " {:22s} = {:8.2f} \n"
330 " {:22s} = {:8.2f} \n"
331 " {:22s} = {:8.2f} \n"
332 " {:22s} = {:8.2f} \n",
333 d_name,
334 "Time integration", getKeyData("integrate_compute_time") * 1.e-6,
335 "Peridynamics force", getKeyData("pd_compute_time") * 1.e-6,
336 "Contact force", getKeyData("contact_compute_time") * 1.e-6,
337 "Search tree update", getKeyData("tree_compute_time") * 1.e-6,
338 "External force", getKeyData("extf_compute_time") * 1.e-6)
339 );
340}
std::string ppCompressiveTest()
Function that handles post-processing for compressive test of particulate media by rigid wall and ret...
virtual void output()
Output the snapshot of data at current time step.
virtual void applyInitialCondition()
Applies initial condition.
virtual void computeExternalDisplacementBC()
Applies displacement boundary conditions.
Definition demModel.cpp:814
virtual void integrateStep()
Performs one time step.
Definition demModel.cpp:342
std::string ppTwoParticleTest()
Function that handles post-processing for two particle collision test and returns maximum vertical di...
virtual void computeForces()
Computes peridynamic forces and contact forces.
Definition demModel.cpp:462
virtual void checkStop()
Checks if simulation should be stopped due to abnormal state of system.

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.

349 {
350
351 // update velocity and displacement
353 const auto dim = d_modelDeck_p->d_dim;
354
355 tf::Executor executor(util::parallel::getNThreads());
356 tf::Taskflow taskflow;
357
358 // update current position, displacement, and velocity of nodes
359 taskflow.for_each_index(
360 (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
361 [this, dim](std::size_t II) {
362 auto i = this->d_fPdCompNodes[II];
363
364 const auto rho = this->getDensity(i);
365 const auto &fix = this->d_fix[i];
366
367 for (int dof = 0; dof < dim; dof++) {
368 if (util::methods::isFree(fix, dof)) {
369 this->d_v[i][dof] += (this->d_currentDt / rho) * this->d_f[i][dof];
370 this->d_u[i][dof] += this->d_currentDt * this->d_v[i][dof];
371 this->d_x[i][dof] += this->d_currentDt * this->d_v[i][dof];
372 }
373 }
374
375 this->d_vMag[i] = this->d_v[i].length();
376 } // loop over nodes
377 ); // for_each
378
379 executor.run(taskflow).get();
380
381 // advance time
382 d_n++;
384
385 // update displacement bc
387
388 // compute force
390}
double d_currentDt
Current timestep.
Definition modelData.h:537

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.

342 {
343 if (d_modelDeck_p->d_timeDiscretization == "central_difference")
344 integrateCD();
345 else if (d_modelDeck_p->d_timeDiscretization == "velocity_verlet")
347}
virtual void integrateVerlet()
Perform time integration using velocity verlet scheme.
Definition demModel.cpp:392
virtual void integrateCD()
Perform time integration using central-difference scheme.
Definition demModel.cpp:349

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.

392 {
393
394 // update velocity and displacement
396 const auto dim = d_modelDeck_p->d_dim;
397
398 // update current position, displacement, and velocity of nodes
399 {
400 tf::Executor executor(util::parallel::getNThreads());
401 tf::Taskflow taskflow;
402
403 taskflow.for_each_index(
404 (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
405 [this, dim](std::size_t II) {
406 auto i = this->d_fPdCompNodes[II];
407
408 const auto rho = this->getDensity(i);
409 const auto &fix = this->d_fix[i];
410
411 for (int dof = 0; dof < dim; dof++) {
412 if (util::methods::isFree(fix, dof)) {
413 this->d_v[i][dof] += 0.5 * (this->d_currentDt / rho) * this->d_f[i][dof];
414 this->d_u[i][dof] += this->d_currentDt * this->d_v[i][dof];
415 this->d_x[i][dof] += this->d_currentDt * this->d_v[i][dof];
416 }
417
418 this->d_vMag[i] = this->d_v[i].length();
419 }
420 } // loop over nodes
421 ); // for_each
422
423 executor.run(taskflow).get();
424 }
425
426 // advance time
427 d_n++;
429
430 // update displacement bc
432
433 // compute force
435
436 // update velocity of nodes
437 {
438 tf::Executor executor(util::parallel::getNThreads());
439 tf::Taskflow taskflow;
440
441 taskflow.for_each_index(
442 (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
443 [this, dim](std::size_t II) {
444 auto i = this->d_fPdCompNodes[II];
445
446 const auto rho = this->getDensity(i);
447 const auto &fix = this->d_fix[i];
448 for (int dof = 0; dof < dim; dof++) {
449 if (util::methods::isFree(fix, dof)) {
450 this->d_v[i][dof] += 0.5 * (this->d_currentDt / rho) * this->d_f[i][dof];
451 }
452
453 this->d_vMag[i] = this->d_v[i].length();
454 }
455 } // loop over nodes
456 ); // for_each
457
458 executor.run(taskflow).get();
459 }
460}

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
strMessage string
priorityPriority of message
check_conditionSpecify whether condition for logging/printing message is to be checked
override_prioritySpecify new priority
screen_outSpecify whether to print the message also to screen (using std::cout)

Definition at line 61 of file demModel.cpp.

62 {
63 int op = override_priority == -1 ? priority : override_priority;
64 if ((check_condition and d_outputDeck_p->d_debug > priority) or d_outputDeck_p->d_debug > op)
65 util::io::log(str, screen_out);
66}
void log(std::ostringstream &oss, bool screen_out=false, int printMpiRank=print_default_mpi_rank)
Global method to log the message.
Definition io.cpp:38

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
ossStream
priorityPriority of message
check_conditionSpecify whether condition for logging/printing message is to be checked
override_prioritySpecify new priority
screen_outSpecify whether to print the message also to screen (using std::cout)

Definition at line 53 of file demModel.cpp.

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

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 1945 of file demModel.cpp.

1945 {
1946
1947 // write out % completion of simulation at 10% interval
1948 {
1949 float p = float(d_n) * 100. / d_modelDeck_p->d_Nt;
1950 int m = std::max(1, int(d_modelDeck_p->d_Nt / 10));
1951 if (d_n % m == 0 && int(p) > 0)
1952 log(fmt::format("{}: Simulation {}% complete\n",
1953 d_name, int(p)));
1954 ;
1955 }
1956
1957 log(fmt::format("{}: Output step = {}, time = {:.6f} \n",
1958 d_name, d_n, d_time),
1959 2);
1960
1961 if (d_outputDeck_p->d_debug > 0 and getKeyData("debug_once") < 0) {
1962
1963 setKeyData("debug_once", 1);
1964
1965 size_t nt = 1;
1966 auto tabS = util::io::getTabS(nt);
1967 std::ostringstream oss;
1968 oss << tabS << "*******************************************\n";
1969 oss << tabS << "Debug various input decks\n\n\n";
1970 oss << d_modelDeck_p->printStr(nt + 1);
1971 oss << d_pDeck_p->printStr(nt + 1);
1972 oss << d_cDeck_p->printStr(nt + 1);
1973 oss << tabS << "\n\n*******************************************\n";
1974 oss << tabS << "Debug particle data\n\n\n";
1975 oss << tabS << "Number of particles = " << d_particlesListTypeAll.size() << std::endl;
1976 oss << tabS << "Number of particle zones = " << d_zInfo.size() << std::endl;
1977 for (auto zone : d_zInfo) {
1978 oss << tabS << "zone of d_zInfo: " << util::io::printStr(zone)
1979 << std::endl;
1980 }
1981
1982 // wall info
1983 oss << tabS << "Number of walls = " << d_particlesListTypeWall.size() << std::endl;
1984 for (auto &d_wall : d_particlesListTypeWall)
1985 oss << tabS << "Number of nodes in wall in zone " << d_wall->d_zoneId
1986 << " is " << d_wall->getNumNodes() << std::endl;
1987
1988 oss << tabS << "h_min = " << d_hMin << ", h_max = " << d_hMax << std::endl;
1989
1990 log(oss, 2);
1991 } // end of debug
1992
1993 size_t dt_out = d_outputDeck_p->d_dtOutCriteria;
1994 std::string out_filename = d_outputDeck_p->d_path + "output_";
1995 if (d_outputDeck_p->d_tagPPFile.empty())
1996 out_filename = out_filename + std::to_string(d_n / dt_out);
1997 else
1998 out_filename = out_filename + d_outputDeck_p->d_tagPPFile + "_" + std::to_string(d_n / dt_out);
1999
2000 auto writer = rw::writer::VtkParticleWriter(out_filename);
2001 if (d_outputDeck_p->d_performFEOut)
2002 writer.appendMesh(this, d_outputDeck_p->d_outTags);
2003 else
2004 writer.appendNodes(this, d_outputDeck_p->d_outTags);
2005
2006 writer.addTimeStep(d_time);
2007 writer.close();
2008
2009 if (util::methods::isTagInList("Strain_Stress", d_outputDeck_p->d_outTags)) {
2010
2011 // compute current position of quadrature points and strain/stress data
2012 {
2013 // if particle mat data is not computed, compute them
2014 if (d_particlesMatDataList.empty()) {
2015 for (auto &p: d_particlesListTypeAll) {
2016 d_particlesMatDataList.push_back(p->getMaterial()->computeMaterialProperties(
2017 p->getMeshP()->getDimension()));
2018 }
2019 }
2020
2021 for (auto &p: d_particlesListTypeAll) {
2022
2023 const auto particle_mesh_p = p->getMeshP();
2024
2025 fe::getCurrentQuadPoints(particle_mesh_p.get(), d_xRef, d_u, d_xQuadCur,
2026 p->d_globStart,
2027 p->d_globQuadStart,
2028 d_modelDeck_p->d_quadOrder);
2029
2030 auto p_z_id = p->d_zoneId;
2031 auto isPlaneStrain = d_pDeck_p->d_particleZones[p_z_id].d_matDeck.d_isPlaneStrain;
2032 fe::getStrainStress(particle_mesh_p.get(), d_xRef, d_u,
2033 isPlaneStrain,
2035 p->d_globStart,
2036 p->d_globQuadStart,
2037 d_particlesMatDataList[p->getId()].d_nu,
2038 d_particlesMatDataList[p->getId()].d_lambda,
2039 d_particlesMatDataList[p->getId()].d_mu,
2040 true,
2041 d_modelDeck_p->d_quadOrder);
2042 } // for loop over particles
2043 } // compute strain/stress block
2044
2045 out_filename = d_outputDeck_p->d_path + "output_strain_";
2046 if (d_outputDeck_p->d_tagPPFile.empty())
2047 out_filename = out_filename + std::to_string(d_n / dt_out);
2048 else
2049 out_filename = out_filename + d_outputDeck_p->d_tagPPFile + "_" + std::to_string(d_n / dt_out);
2050
2051 auto writer1 = rw::writer::VtkParticleWriter(out_filename);
2052 writer1.appendStrainStress(this);
2053 writer1.addTimeStep(d_time);
2054 writer1.close();
2055 }
2056
2057 // output particle locations to csv file
2058 if (util::methods::isTagInList("Particle_Locations",
2059 d_outputDeck_p->d_outTags)) {
2060
2061 out_filename = d_outputDeck_p->d_path + "particle_locations_";
2062 if (d_outputDeck_p->d_tagPPFile.empty())
2063 out_filename = out_filename + std::to_string(d_n / dt_out) + ".csv";
2064 else
2065 out_filename = out_filename + d_outputDeck_p->d_tagPPFile
2066 + "_" + std::to_string(d_n / dt_out) + ".csv";
2067
2068 std::ofstream oss(out_filename);
2069 oss << "i, x, y, z, r\n";
2070 for (const auto &p : d_particlesListTypeAll) {
2071 auto xc = p->getXCenter();
2072 oss << p->d_zoneId << ", " << xc.d_x << ", " << xc.d_y << ", " << xc.d_z
2073 << ", " << p->d_geom_p->boundingRadius() << "\n";
2074 }
2075 oss.close();
2076 }
2077}
double d_hMax
Minimum mesh over all particles and walls.
Definition modelData.h:571
std::vector< util::Point > d_xRef
reference positions of the nodes
Definition modelData.h:630
std::vector< util::Point > d_u
Displacement of the nodes.
Definition modelData.h:636
std::vector< util::SymMatrix3 > d_stress
Stress in elements (values at quadrature points)
Definition modelData.h:737
std::vector< util::Point > d_xQuadCur
Current position of quadrature points.
Definition modelData.h:731
void setKeyData(std::string key, double data, bool issue_err=false)
Set value to data associated with key.
Definition modelData.h:163
std::vector< util::SymMatrix3 > d_strain
Strain in elements (values at quadrature points)
Definition modelData.h:734
double d_hMin
Maximum mesh over all particles and walls.
Definition modelData.h:574
std::vector< inp::MatData > d_particlesMatDataList
List of particle material data. Only populated if needed for calculation of stress or other quantitie...
Definition modelData.h:604
A vtk writer for simple point data and complex fem mesh data.
void getStrainStress(const fe::Mesh *mesh_p, const std::vector< util::Point > &xRef, const std::vector< util::Point > &u, bool isPlaneStrain, std::vector< util::SymMatrix3 > &strain, std::vector< util::SymMatrix3 > &stress, size_t iNodeStart=0, size_t iStrainStart=0, double nu=0., double lambda=0., double mu=0., bool computeStress=false, size_t quadOrder=1)
Strain and stress at quadrature points in the mesh.
Definition meshUtil.cpp:280
void getCurrentQuadPoints(const fe::Mesh *mesh_p, const std::vector< util::Point > &xRef, const std::vector< util::Point > &u, std::vector< util::Point > &xQuadCur, size_t iNodeStart=0, size_t iQuadStart=0, size_t quadOrder=1)
Get current location of quadrature points of elements in the mesh. This function expects mesh has ele...
Definition meshUtil.cpp:176
std::string getTabS(int nt)
Returns tab spaces of given size.
Definition io.h:40

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 2185 of file demModel.cpp.

2185 {
2186 bool continue_dt = false;
2187 auto check_dt = d_outputDeck_p->d_dtTestOut;
2188 if ((d_n % check_dt == 0) && (d_n >= check_dt))
2189 continue_dt = true;
2190
2191 if (!continue_dt)
2192 return "";
2193
2194 // get wall
2195 auto w_id = d_pDeck_p->d_particleIdCompressiveTest;
2196 auto f_dir = d_pDeck_p->d_particleForceDirectionCompressiveTest - 1;
2197 const auto &wall = d_particlesListTypeAll[w_id];
2198
2199 // find the penetration of the wall from it's original location
2200 auto dx = wall->getXLocal(0) - wall->getXRefLocal(0);
2201 double wall_penetration = dx[f_dir];
2202
2203 // get the total reaction force on wall along the direction of loading
2204 double tot_reaction_force = 0.;
2205 for (size_t i = 0; i < wall->getNumNodes(); i++) {
2206 tot_reaction_force += wall->getFLocal(i)[f_dir] * wall->getVolLocal(i);
2207 }
2208
2209 // open file and write
2210 bool use_static_file = true;
2211 if (use_static_file) {
2212 if (!d_ppFile.is_open()) {
2213
2214 std::string tag_pp_file = d_outputDeck_p->d_tagPPFile.empty() ? "0" : d_outputDeck_p->d_tagPPFile;
2215 std::string filename = d_outputDeck_p->d_path + "pp_" +
2216 d_pDeck_p->d_testName + "_" +
2217 tag_pp_file + ".csv";
2218 d_ppFile.open(filename.c_str(), std::ofstream::out | std::ofstream::app);
2219
2220 d_ppFile << "t, delta, force \n";
2221 }
2222
2223 d_ppFile << fmt::format("%4.6e, %4.6e, %4.6e\n", d_time, wall_penetration,
2224 tot_reaction_force);
2225 }
2226
2227 setKeyData("wall_penetration", wall_penetration);
2228 setKeyData("tot_reaction_force", tot_reaction_force);
2229
2230 return fmt::format(" Post-processing: wall penetration = {:"
2231 ".6f}, "
2232 "reaction force = {:5.3e} \n",
2233 wall_penetration, tot_reaction_force);
2234}

◆ 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 2079 of file demModel.cpp.

2079 {
2080
2081 bool continue_dt = false;
2082 auto check_dt = d_outputDeck_p->d_dtTestOut;
2083 if ((d_n % check_dt == 0) && (d_n >= check_dt))
2084 continue_dt = true;
2085
2086 if (!continue_dt)
2087 return "";
2088
2089 // get alias for particles
2090 const auto &p0 = this->d_particlesListTypeAll[0];
2091 const auto &p1 = this->d_particlesListTypeAll[1];
2092
2093 // get penetration distance
2094 const auto &xc0 = p0->getXCenter();
2095 const auto &xc1 = p1->getXCenter();
2096 const double &r = p0->d_geom_p->boundingRadius();
2097
2098 const auto &contact = d_cDeck_p->getContact(p0->d_zoneId, p1->d_zoneId);
2099 double r_e = r + contact.d_contactR;
2100
2101 double pen_dist = xc1.dist(xc0) - r_e - r;
2102 double contact_area_radius = 0.;
2103 if (util::isLess(pen_dist, 0.))
2104 contact_area_radius =
2105 std::sqrt(std::pow(r_e, 2.) - std::pow(r_e + pen_dist, 2.));
2106 else if (util::isGreater(pen_dist, 0.)) {
2107 pen_dist = 0.;
2108 contact_area_radius = 0.;
2109 }
2110
2111 // get max distance of second particle (i.e. the y-coord of center + radius)
2112 double max_dist = xc1.d_y + p1->d_geom_p->boundingRadius();
2113
2114 // compute maximum y coordinate of particle 2
2115 double max_y_loc = p1->getXLocal(0).d_y;
2116 double max_y = 0.;
2117 for (size_t i = 0; i < p1->getNumNodes(); i++)
2118 if (util::isLess(max_y_loc, p1->getXLocal(i).d_y))
2119 max_y_loc = p1->getXLocal(i).d_y;
2120
2121 if (util::isLess(max_y, max_y_loc))
2122 max_y = max_y_loc;
2123
2124 setKeyData("pen_dist", pen_dist);
2125 setKeyData("contact_area_radius", contact_area_radius);
2126 setKeyData("max_y", max_y);
2127 setKeyData("max_dist", max_dist);
2128 setKeyData("max_y_loc", max_y_loc);
2129
2130
2131 return fmt::format(" Post-processing: max y = {:.6f} \n", max_y);
2132}

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
deckInput deck

Definition at line 84 of file demModel.cpp.

84 {
85
86 log(d_name + ": Restarting the simulation\n");
87
88 // set time step to step specified in restart deck
89 d_n = d_restartDeck_p->d_step;
90 d_time = double(d_n) * d_modelDeck_p->d_dt;
91 log(fmt::format(" Restart step = {}, time = {:.6f} \n", d_n, d_time));
92
93 // get backup of reference configuration
94 std::vector<util::Point> x_ref(d_x.size(), util::Point());
95 for (auto &x : d_x)
96 x_ref.push_back(x);
97
98 // read displacement and velocity from restart file
99 log(" Reading data from restart file = " + d_restartDeck_p->d_file + " \n");
100 auto reader = rw::reader::VtkParticleReader(d_restartDeck_p->d_file);
101 reader.readNodes(this);
102}
std::shared_ptr< inp::RestartDeck > d_restartDeck_p
Restart deck.
Definition modelData.h:555
A vtk writer for simple point data and complex fem mesh data.

◆ run()

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

Main driver to simulate.

Parameters
deckInput deck

Reimplemented in twoparticle_demo::Model.

Definition at line 68 of file demModel.cpp.

68 {
69
70 // initialize data
71 init();
72
73 // check for restart
74 if (d_modelDeck_p->d_isRestartActive)
75 restart(deck);
76
77 // integrate in time
78 integrate();
79
80 // close
81 close();
82}
virtual void integrate()
Perform time integration.
Definition demModel.cpp:268
virtual void close()
Closure operations.
Definition demModel.cpp:104
virtual void init()
Initialize remaining data members.
Definition demModel.cpp:109
virtual void restart(inp::Input *deck)
Restarts the simulation from previous state.
Definition demModel.cpp:84

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 1514 of file demModel.cpp.

1514 {
1515
1516 // loop over all particle zones and get minimum value of mesh size
1517 size_t c = 0;
1518 for (const auto *p : d_particlesListTypeAll) {
1519
1520 auto h = p->getMeshSize();
1521 if (c == 0) {
1522 d_hMin = h;
1523 d_hMax = h;
1524 c++;
1525 }
1526
1527 if (util::isGreater(d_hMin, h))
1528 d_hMin = h;
1529 if (util::isGreater(h, d_hMax))
1530 d_hMax = h;
1531 }
1532
1533 log(fmt::format("{}: Contact setup\n hmin = {:.6f}, hmax = {:.6f} \n",
1534 d_name, d_hMin, d_hMax), 1);
1535
1536 d_maxContactR = 0.;
1537 // precompute bulk modulus of all zones
1538 std::vector<double> bulk_modulus;
1539 // NOTE - d_data.size() and d_zoneVec.size() are equal
1540 for (size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
1541
1542 double kappa_i = d_pDeck_p->d_particleZones[i].d_matDeck.d_matData.d_K;
1543
1544 if (kappa_i < 0.) {
1545 std::cerr << "Error: We need bulk modulus provided in input file.\n";
1546 std::cerr << d_pDeck_p->d_particleZones[i].printStr();
1547 exit(1);
1548 }
1549
1550 bulk_modulus.push_back(kappa_i);
1551 }
1552
1553 for (size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
1554 for (size_t j = 0; j < d_cDeck_p->d_data.size(); j++) {
1555
1556 inp::ContactPairDeck *deck = &(d_cDeck_p->d_data[i][j]);
1557
1558 if (deck->d_computeContactR)
1559 deck->d_contactR *= d_hMin;
1560
1561 if (d_maxContactR < deck->d_contactR)
1562 d_maxContactR = deck->d_contactR;
1563
1564 // get effective bulk modulus for pair of zones and store it
1565 deck->d_kappa = util::equivalentMass(bulk_modulus[i], bulk_modulus[j]);
1566
1567 // Kn
1568 deck->d_Kn *= deck->d_KnFactor;
1569
1570 // Beta n
1571 double log_e = std::log(deck->d_eps);
1572 deck->d_betan =
1573 deck->d_betanFactor *
1574 (-2. * log_e * std::sqrt(1. / (M_PI * M_PI + log_e * log_e)));
1575
1576 log(fmt::format(" contact_radius = {:.6f}, hmin = {:.6f}, Kn = {:5.3e}, "
1577 "Vmax = {:5.3e}, "
1578 "betan = {:7.5f}, mu = {:.4f}, kappa = {:5.3e}\n",
1579 deck->d_contactR, d_hMin, deck->d_Kn, deck->d_vMax,
1580 deck->d_betan, deck->d_mu, deck->d_kappa), 2);
1581 }
1582 }
1583}
float h
mesh size (use smaller radius to decide)
Structure to read and store particle-particle contact related input data.
Definition contactDeck.h:23
double d_contactR
contact radius
Definition contactDeck.h:26
double d_eps
parameters for normal damping force
Definition contactDeck.h:40
bool d_computeContactR
Flag that indicates whether contact radius is to be computed.
Definition contactDeck.h:29
double d_kappa
parameters for frictional force
Definition contactDeck.h:51
double d_betan
parameters for normal damping force
Definition contactDeck.h:41
double d_betanFactor
parameters for frictional force
Definition contactDeck.h:50
double d_vMax
parameters for normal force
Definition contactDeck.h:33
double d_KnFactor
parameters for frictional force
Definition contactDeck.h:49
double d_Kn
parameters for normal force
Definition contactDeck.h:35
double d_mu
parameters for frictional force
Definition contactDeck.h:46

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 1585 of file demModel.cpp.

1585 {
1586
1587 if (util::methods::isTagInList("Strain_Stress", d_outputDeck_p->d_outTags)
1588 or d_modelDeck_p->d_populateElementNodeConnectivity) {
1589
1590 // read element-node connectivity data if not done
1591 for (auto &p: d_referenceParticles) {
1592 auto &particle_mesh_p = p->getMeshP();
1593 if (!particle_mesh_p->d_encDataPopulated && particle_mesh_p->d_enc.empty()) {
1594 particle_mesh_p->readElementData(particle_mesh_p->d_filename);
1595 }
1596 }
1597
1598 // setup quadrature point and strain/stress data
1599 // we need to know size of the data
1600 size_t totalQuadPoints = 0;
1601 for (auto &p: d_particlesListTypeAll) {
1602 const auto &particle_mesh_p = p->getMeshP();
1603
1604 // get Quadrature
1605 fe::BaseElem *elem;
1606 if (particle_mesh_p->getElementType() == util::vtk_type_line)
1607 elem = new fe::LineElem(d_modelDeck_p->d_quadOrder);
1608 else if (particle_mesh_p->getElementType() == util::vtk_type_triangle)
1609 elem = new fe::TriElem(d_modelDeck_p->d_quadOrder);
1610 else if (particle_mesh_p->getElementType() == util::vtk_type_quad)
1611 elem = new fe::QuadElem(d_modelDeck_p->d_quadOrder);
1612 else if (particle_mesh_p->getElementType() == util::vtk_type_tetra)
1613 elem = new fe::TetElem(d_modelDeck_p->d_quadOrder);
1614 else {
1615 std::cerr << fmt::format("Error: Can not compute strain/stress as the element "
1616 "type = {} is not yet supported in this routine.\n", particle_mesh_p->getElementType());
1617 exit(EXIT_FAILURE);
1618 }
1619
1620 p->d_globQuadStart = totalQuadPoints;
1621 totalQuadPoints += particle_mesh_p->getNumElements() *
1622 elem->getNumQuadPoints();
1623 p->d_globQuadEnd = totalQuadPoints;
1624
1625 std::cout << fmt::format("p->id() = {}, "
1626 "p->d_globQuadStart = {}, "
1627 "totalQuadPoints = {}, "
1628 "p->d_globQuadEnd = {}",
1629 p->getId(), p->d_globQuadStart,
1630 particle_mesh_p->getNumElements() *
1631 elem->getNumQuadPoints(), p->d_globQuadEnd)
1632 << std::endl;
1633 }
1634
1635 // resize data
1636 d_xQuadCur.resize(totalQuadPoints);
1637 d_strain.resize(totalQuadPoints);
1638 d_stress.resize(totalQuadPoints);
1639 } // setting up quadrature data
1640}
A base class which provides methods to map points to/from reference element and to compute quadrature...
Definition baseElem.h:84
size_t getNumQuadPoints()
Get number of quadrature points in the data.
Definition baseElem.h:111
A class for mapping and quadrature related operations for linear 2-node line element.
Definition lineElem.h:49
A class for mapping and quadrature related operations for bi-linear quadrangle element.
Definition quadElem.h:64
A class for mapping and quadrature related operations for linear tetrahedron element.
Definition tetElem.h:141
A class for mapping and quadrature related operations for linear triangle element.
Definition triElem.h:91
static const int vtk_type_triangle
Integer flag for triangle element.
static const int vtk_type_quad
Integer flag for quad element.
static const int vtk_type_tetra
Integer flag for tetrahedron element.
static const int vtk_type_line
Integer flag for line element.

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 1679 of file demModel.cpp.

1679 {
1680
1682
1683 if (!update)
1684 return;
1685
1686 // update contact neighborlist
1687
1688 // update the point cloud (make sure that d_x is updated along with displacement)
1689 auto pt_cloud_update_time = d_nsearch_p->setInputCloud();
1690 setKeyData("pt_cloud_update_time", pt_cloud_update_time);
1691 appendKeyData("tree_compute_time", pt_cloud_update_time);
1692 appendKeyData("avg_tree_update_time", pt_cloud_update_time/d_infoN);
1693
1694 if (d_neighC.size() != d_x.size())
1695 d_neighC.resize(d_x.size());
1696
1697 tf::Executor executor(util::parallel::getNThreads());
1698 tf::Taskflow taskflow;
1699
1700 taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1,
1701 [this](std::size_t i) {
1702
1703 const auto &pi = this->d_ptId[i];
1704 const auto &pi_particle = this->d_particlesListTypeAll[pi];
1705
1706 // search?
1707 bool perform_search_based_on_particle = true;
1708 if (pi_particle->d_typeIndex == 1) // wall
1709 perform_search_based_on_particle = false;
1710
1711 if (pi_particle->d_allDofsConstrained or !pi_particle->d_computeForce)
1712 perform_search_based_on_particle = false;
1713
1714 if (perform_search_based_on_particle) {
1715
1716 std::vector<size_t> neighs;
1717 std::vector<double> sqr_dist;
1718
1719 this->d_neighC[i].clear();
1720
1721 auto n = this->d_nsearch_p->radiusSearchExcludeTag(
1722 this->d_x[i],
1723 this->d_contNeighSearchRadius,
1724 neighs,
1725 sqr_dist,
1726 this->d_ptId[i],
1727 this->d_ptId);
1728
1729 if (n > 0) {
1730 for (auto neigh: neighs) {
1731 if (neigh != i)
1732 this->d_neighC[i].push_back(neigh);
1733 }
1734 }
1735 }
1736}
1737 ); // for_each
1738
1739 executor.run(taskflow).get();
1740
1741
1742 // handle particle-wall neighborlist (based on the d_neighC that we already computed)
1746
1747 for (auto &pi : d_particlesListTypeParticle) {
1748
1749 d_neighWallNodes[pi->getId()].resize(pi->getNumNodes());
1750 d_neighWallNodesDistance[pi->getId()].resize(pi->getNumNodes());
1751
1752 // get all wall nodes that are within contact distance to the nodes of this particle
1753 {
1754 tf::Executor executor(util::parallel::getNThreads());
1755 tf::Taskflow taskflow;
1756
1757 taskflow.for_each_index((std::size_t) 0,
1758 pi->getNumNodes(),
1759 (std::size_t) 1,
1760 [this, &pi](std::size_t i) {
1761
1762 auto i_glob = pi->getNodeId(i);
1763 auto yi = this->d_x[i_glob];
1764
1765 const std::vector<size_t> &neighs = this->d_neighC[i_glob];
1766
1767 this->d_neighWallNodes[pi->getId()][i].clear();
1768 this->d_neighWallNodesDistance[pi->getId()][i].clear();
1769
1770 for (const auto &j_id: neighs) {
1771
1772 auto &ptIdj = this->d_ptId[j_id];
1773 auto &pj = this->getParticleFromAllList(
1774 ptIdj);
1775
1776 // we are only interested in nodes from wall
1777 if (pj->getTypeIndex() == 1) {
1778 this->d_neighWallNodes[pi->getId()][i].push_back(j_id);
1779 //this->d_neighWallNodesDistance[pi->getId()][i].push_back(Rji);
1780 }
1781 }
1782 }
1783 ); // for_each
1784
1785 executor.run(taskflow).get();
1786 }
1787 } // loop over particles
1788
1789}
virtual bool updateContactNeighborSearchParameters()
Update contact neighbor search parameters.
std::vector< std::vector< size_t > > d_neighC
Neighbor data for contact forces.
Definition modelData.h:655
std::vector< std::vector< std::vector< double > > > d_neighWallNodesDistance
Neighbor data (distance) for contact between particle and walls.
Definition modelData.h:667

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 1791 of file demModel.cpp.

1791 {
1792
1793 // initialize parameters
1794 if (d_contNeighUpdateInterval == 0 and
1796 d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
1798 d_contNeighSearchRadius = d_maxContactR * d_pDeck_p->d_pNeighDeck.d_sFactor;
1799 }
1800
1801 // at d_n = 0, this function will be called twice because updateContactNeighborlist() will be
1802 // called twice: one inside init() and second inside computeForces()
1803 // so to match d_n and d_contNeighTimestepCounter in the initial stage of simulation, we need to handle the special case
1804 if (d_n == 0) {
1805 appendKeyData("update_contact_neigh_search_params_init_call_count", 1);
1806
1807 if (int(getKeyData("update_contact_neigh_search_params_init_call_count")) == 1)
1808 return true;
1809
1810 if (int(getKeyData("update_contact_neigh_search_params_init_call_count")) == 2) {
1813 }
1814 }
1815
1816 // handle case of restart
1817 if (d_modelDeck_p->d_isRestartActive and d_n == d_restartDeck_p->d_step) {
1818 // assign correct value for restart step
1820 }
1821
1822 if (d_contNeighUpdateInterval == 1) {
1823 // further optimization of parameters is not possible
1825
1826 // update counter and return condition for contact search
1829 }
1830
1831 // check if we should proceed with parameter update
1832 // param update is done at smaller interval than the search itself to avoid
1833 // scenarios where particles suddenly move with a high velocity
1834 size_t update_param_interval =
1835 d_contNeighUpdateInterval > 5 ? size_t(
1836 0.2 * d_contNeighUpdateInterval) : 1;
1837
1838 // check if we ought to update search parameters; if not, return
1839 if (d_contNeighTimestepCounter > 0 and d_contNeighTimestepCounter % update_param_interval != 0) {
1840 // update counter and return condition for contact search
1843 }
1844
1845 // first update the maximum velocity in all particles
1846 for (auto &pi : d_particlesListTypeAll) {
1847 auto max_v_node = util::methods::maxIndex(d_vMag,
1848 pi->d_globStart, pi->d_globEnd);
1849
1850 if (max_v_node > pi->d_globEnd or max_v_node < pi->d_globStart) {
1851 std::cerr << fmt::format("Error: max_v_node = {} for "
1852 "particle of id = {} is not in the limit.\n",
1853 max_v_node, pi->getId())
1854 << "Particle info = \n"
1855 << pi->printStr()
1856 << "\n\n Magnitude of velocity = "
1857 << d_vMag[max_v_node] << "\n";
1858 exit(EXIT_FAILURE);
1859 }
1860
1862 = d_vMag[max_v_node];
1863 }
1864
1865 // find max velocity among all particles
1867
1868 // now we find the best parameters for contact search
1869 auto up_interval_old = d_contNeighUpdateInterval;
1870
1871 // TO ensure that in d_neighUpdateInterval time steps, the search radius is above the
1872 // distance traveled by object with velocity d_maxVelocity
1873 // also multiply by a safety factor
1874 double safety_factor = d_pDeck_p->d_pNeighDeck.d_sFactor > 5 ? d_pDeck_p->d_pNeighDeck.d_sFactor : 10;
1875 auto max_search_r_from_contact_R = d_pDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
1876 auto max_search_r = d_maxVelocity * d_currentDt
1877 * d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval
1878 * safety_factor;
1879
1880
1881 if (util::isGreater(max_search_r, max_search_r_from_contact_R )) {
1882
1884 if (up_interval_old > d_contNeighUpdateInterval) {
1885 // issue warning
1886 log(fmt::format("Warning: Contact search radius based on velocity is greater than "
1887 "the max contact radius.\n"
1888 "Warning: Adjusting contact neighborlist update interval.\n"
1889 "{:>13} = {:4.6e}, time step = {}, "
1890 "velocity-based r = {:4.6e}, max contact r = {:4.6e}\n",
1891 "Time", d_time, d_n, max_search_r, max_search_r_from_contact_R),
1892 2, d_n % d_infoN == 0, 3);
1893 }
1894
1895 d_contNeighSearchRadius = max_search_r_from_contact_R;
1896 // reset time step counter for contact so that the contact list is updated in the current time step
1897 // and the update cycle starts from the current time step
1899
1900 if (d_contNeighUpdateInterval < 1) {
1903 }
1904 }
1905 else {
1906 // update search radius
1907 d_contNeighSearchRadius = d_contNeighUpdateInterval < 2 ? d_maxContactR : max_search_r_from_contact_R;
1908 }
1909
1910 if (up_interval_old > d_contNeighUpdateInterval) {
1911 log(fmt::format(" Contact neighbor parameters: \n"
1912 " {:48s} = {:d}\n"
1913 " {:48s} = {:d}\n"
1914 " {:48s} = {:d}\n"
1915 " {:48s} = {:4.6e}\n"
1916 " {:48s} = {:4.6e}\n"
1917 " {:48s} = {:4.6e}\n"
1918 " {:48s} = {:4.6e}\n"
1919 " {:48s} = {:4.6e}\n"
1920 " {:48s} = {:4.6e}\n",
1921 "time step", d_n,
1922 "contact neighbor update interval",
1924 "contact neighbor update time step counter",
1926 "search radius", d_contNeighSearchRadius,
1927 "max contact radius", d_maxContactR,
1928 "search radius factor", d_pDeck_p->d_pNeighDeck.d_sFactor,
1929 "max search r from velocity", max_search_r,
1930 "max search r from contact r", max_search_r_from_contact_R,
1931 "max velocity", d_maxVelocity),
1932 2, d_n % d_infoN == 0, 3);
1933 }
1934
1935 // update counter and return condition for contact search
1938}
std::vector< double > d_vMag
Magnitude of velocity of the nodes.
Definition modelData.h:642
size_t d_contNeighTimestepCounter
Contact neighborlist time step counter.
Definition modelData.h:583
size_t maxIndex(const std::vector< T > &data)
Returns the index corresponding to maximum from list of data.
Definition methods.h:38

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 1483 of file demModel.cpp.

1483 {
1484
1485 for (auto &p: d_particlesListTypeAll) {
1486 if (p->d_geom_p->d_name == "null" or
1487 util::methods::isTagInList("copy_from_container", p->d_geom_p->d_tags)) {
1488 // update geometry of particle based on bounding box
1489 auto bbox = p->getMeshP()->getBoundingBox();
1490
1491 std::string geom_name = "rectangle";
1492 if ( p->getMeshP()->getDimension() == 3)
1493 geom_name = "cuboid";
1494
1495 std::vector<double> geom_params(6, 0.);
1496 for (size_t i=0; i<3; i++) {
1497 geom_params[i] = bbox.first[i];
1498 geom_params[i+3] = bbox.second[i];
1499 }
1500
1501 std::vector<std::string> vec_type;
1502 std::vector<std::string> vec_flag;
1503
1505 geom_params,
1506 vec_type,
1507 vec_flag,
1508 p->d_geom_p,
1509 p->getMeshP()->getDimension());
1510 }
1511 }
1512}

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 1940 of file demModel.cpp.

1940 {
1941 // Not used
1942 return;
1943}

◆ updatePeridynamicNeighborlist()

void model::DEMModel::updatePeridynamicNeighborlist ( )
virtual

Update neighborlist for peridynamics force.

Definition at line 1642 of file demModel.cpp.

1642 {
1643
1644 d_neighPd.resize(d_x.size());
1645 // d_neighPdSqdDist.resize(d_x.size());
1646 auto t1 = steady_clock::now();
1647
1648 tf::Executor executor(util::parallel::getNThreads());
1649 tf::Taskflow taskflow;
1650
1651 taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1, [this](std::size_t i) {
1652 const auto &pi = this->d_ptId[i];
1653 double search_r = this->d_particlesListTypeAll[pi]->d_material_p->getHorizon();
1654
1655 std::vector<size_t> neighs;
1656 std::vector<double> sqr_dist;
1657 if (this->d_nsearch_p->radiusSearchIncludeTag(this->d_x[i],
1658 search_r,
1659 neighs,
1660 sqr_dist,
1661 this->d_ptId[i],
1662 this->d_ptId) > 0) {
1663 for (std::size_t j = 0; j < neighs.size(); ++j)
1664 if (neighs[j] != i && this->d_ptId[neighs[j]] == pi) {
1665 this->d_neighPd[i].push_back(size_t(neighs[j]));
1666 // this->d_neighPdSqdDist[i].push_back(sqr_dist[j]);
1667 }
1668 }
1669 }
1670 ); // for_each
1671
1672 executor.run(taskflow).get();
1673
1674 auto t2 = steady_clock::now();
1675 log(fmt::format("{}: Peridynamics neighbor update time = {}\n",
1676 d_name, util::methods::timeDiff(t1, t2)), 2);
1677}

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: