37#include <fmt/format.h>
40#include <taskflow/taskflow/taskflow.hpp>
41#include <taskflow/taskflow/algorithm/for_each.hpp>
53void model::DEMModel::log(std::ostringstream &oss,
int priority,
bool check_condition,
int override_priority,
55 int op = override_priority == -1 ? priority : override_priority;
57 if ((check_condition and d_outputDeck_p->d_debug > priority) or d_outputDeck_p->d_debug > op)
61void model::DEMModel::log(
const std::string &str,
int priority,
bool check_condition,
int override_priority,
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)
74 if (d_modelDeck_p->d_isRestartActive)
86 log(d_name +
": Restarting the simulation\n");
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));
94 std::vector<util::Point> x_ref(d_x.size(),
util::Point());
99 log(
" Reading data from restart file = " + d_restartDeck_p->d_file +
" \n");
101 reader.readNodes(
this);
105 if (d_ppFile.is_open())
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;
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);
141 auto t1 = steady_clock::now();
142 auto t2 = steady_clock::now();
143 log(d_name +
": Initializing objects.\n");
146 log(d_name +
": Creating particles.\n");
149 log(d_name +
": Creating maximum velocity data for particles.\n");
150 d_maxVelocityParticlesListTypeAll
151 = std::vector<double>(d_particlesListTypeAll.size(), 0.);
155 if (d_input_p->isMultiParticle()) {
156 log(d_name +
": Setting up contact.\n");
161 log(d_name +
": Setting up element-node connectivity data for strain/stress.\n");
162 setupQuadratureData();
165 log(d_name +
": Creating neighbor search tree.\n");
168 d_nsearch_p = std::make_unique<NSearch>(d_x, d_outputDeck_p->d_debug);
171 double set_tree_time = d_nsearch_p->setInputCloud();
172 log(fmt::format(
"{}: Tree setup time (ms) = {}. \n", d_name, set_tree_time));
175 log(d_name +
": Creating neighborlist for peridynamics.\n");
176 t1 = steady_clock::now();
177 updatePeridynamicNeighborlist();
178 t2 = steady_clock::now();
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();
186 updateContactNeighborlist();
187 t2 = steady_clock::now();
192 log(d_name +
": Creating peridynamics bonds.\n");
193 d_fracture_p = std::make_unique<geometry::Fracture>(&d_x, &d_neighPd);
196 log(d_name +
": Compute state-based peridynamic quantities.\n");
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);
206 log(d_name +
": Initializing force loading object.\n");
208 std::make_unique<loading::ParticleFLoading>(d_pDeck_p->d_forceDeck);
218 if (d_pDeck_p->d_testName ==
"two_particle") {
219 d_particlesListTypeAll[0]->d_computeForce =
false;
222 log(fmt::format(
"{}: Total particles = {}. \n",
223 d_name, d_particlesListTypeAll.size()));
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()));
229 log(d_name +
": Creating list of nodes on which force is to be computed.\n");
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);
244 d_Z = std::vector<float>(d_x.size(), 0.);
246 t2 = steady_clock::now();
247 log(fmt::format(
"{}: Total setup time (ms) = {}. \n",
251 size_t free_dofs = 0;
252 for (
const auto &f : d_fix) {
253 for (
size_t dof = 0; dof < 3; dof++)
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",
261 d_name, d_particlesListTypeAll.size(),
262 d_particlesListTypeParticle.size(),
263 d_particlesListTypeWall.size(),
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),
279 applyInitialCondition();
282 computeExternalDisplacementBC();
285 for (
size_t i = d_n; i < d_modelDeck_p->d_Nt; i++) {
287 log(fmt::format(
"{}: Time step: {}, time: {:8.6f}, steps completed = {}%\n",
291 float(i) * 100. / d_modelDeck_p->d_Nt),
292 2, d_n % d_infoN == 0, 3);
294 auto t1 = steady_clock::now();
295 log(
"Integrating\n",
false, 0, 3);
297 double integrate_time =
300 appendKeyData(
"integrate_compute_time", integrate_time,
true);
302 log(fmt::format(
" Integration time (ms) = {}\n", integrate_time), 2, d_n % d_infoN == 0, 3);
304 if (d_pDeck_p->d_testName ==
"two_particle") {
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);
316 if ((d_n % d_outputDeck_p->d_dtOut == 0) &&
317 (d_n >= d_outputDeck_p->d_dtOut) && d_outputDeck_p->d_performOut) {
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",
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)
343 if (d_modelDeck_p->d_timeDiscretization ==
"central_difference")
345 else if (d_modelDeck_p->d_timeDiscretization ==
"velocity_verlet")
352 d_currentDt = d_modelDeck_p->d_dt;
353 const auto dim = d_modelDeck_p->d_dim;
356 tf::Taskflow taskflow;
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];
364 const auto rho = this->getDensity(i);
365 const auto &fix = this->d_fix[i];
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];
375 this->d_vMag[i] = this->d_v[i].length();
379 executor.run(taskflow).get();
383 d_time += d_currentDt;
386 computeExternalDisplacementBC();
395 d_currentDt = d_modelDeck_p->d_dt;
396 const auto dim = d_modelDeck_p->d_dim;
401 tf::Taskflow taskflow;
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];
408 const auto rho = this->getDensity(i);
409 const auto &fix = this->d_fix[i];
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];
418 this->d_vMag[i] = this->d_v[i].length();
423 executor.run(taskflow).get();
428 d_time += d_currentDt;
431 computeExternalDisplacementBC();
439 tf::Taskflow taskflow;
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];
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];
453 this->d_vMag[i] = this->d_v[i].length();
458 executor.run(taskflow).get();
464 bool dbg_condition = d_n % d_infoN == 0;
466 log(
" Compute forces \n", 2, dbg_condition, 3);
469 auto t1 = steady_clock::now();
471 tf::Taskflow taskflow;
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(); }
478 executor.run(taskflow).get();
482 t1 = steady_clock::now();
483 computePeridynamicForces();
485 appendKeyData(
"pd_compute_time", pd_time);
486 appendKeyData(
"avg_peridynamics_force_time", pd_time/d_infoN);
488 float current_contact_neigh_update_time = 0;
489 float contact_time = 0;
490 if (d_input_p->isMultiParticle()) {
492 t1 = steady_clock::now();
493 updateContactNeighborlist();
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);
502 t1 = steady_clock::now();
503 computeContactForces();
505 appendKeyData(
"contact_compute_time", contact_time);
506 appendKeyData(
"avg_contact_force_time", contact_time / d_infoN);
510 t1 = steady_clock::now();
511 computeExternalForces();
513 appendKeyData(
"extf_compute_time", extf_time);
514 appendKeyData(
"avg_extf_compute_time", extf_time/d_infoN);
518 if (d_input_p->isMultiParticle()) {
519 log(fmt::format(
" Avg time (ms): \n"
526 "tree update",
size_t(getKeyData(
"avg_tree_update_time")),
527 "contact neigh update",
528 size_t(getKeyData(
"avg_contact_neigh_update_time")),
530 size_t(getKeyData(
"avg_contact_force_time")),
531 "total contact",
size_t(getKeyData(
"avg_tree_update_time")
533 "avg_contact_neigh_update_time")
535 "avg_contact_force_time")),
536 "peridynamics force",
537 size_t(getKeyData(
"avg_peridynamics_force_time")),
539 size_t(getKeyData(
"avg_extf_compute_time") / d_infoN)),
540 2, dbg_condition, 3);
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.);
549 log(fmt::format(
" Avg time (ms): \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);
556 appendKeyData(
"avg_peridynamics_force_time", 0.);
557 appendKeyData(
"avg_extf_compute_time", 0.);
561 log(fmt::format(
" {:50s} = {:8d} \n",
562 "Force reset time (ms)",
563 size_t(force_reset_time)
565 2, dbg_condition, 3);
567 log(fmt::format(
" {:50s} = {:8d} \n",
568 "External force time (ms)",
571 2, dbg_condition, 3);
573 log(fmt::format(
" {:50s} = {:8d} \n",
574 "Peridynamics force time (ms)",
577 2, dbg_condition, 3);
579 if (d_input_p->isMultiParticle()) {
581 log(fmt::format(
" {:50s} = {:8d} \n",
582 "Point cloud update time (ms)",
583 size_t(getKeyData(
"pt_cloud_update_time"))
585 2, dbg_condition, 3);
587 log(fmt::format(
" {:50s} = {:8d} \n",
588 "Contact neighborlist update time (ms)",
589 size_t(current_contact_neigh_update_time)
591 2, dbg_condition, 3);
593 log(fmt::format(
" {:50s} = {:8d} \n",
594 "Contact force time (ms)",
597 2, dbg_condition, 3);
604 log(
" Computing peridynamic force \n", 3);
606 const auto dim = d_modelDeck_p->d_dim;
607 const bool is_state = d_particlesListTypeAll[0]->getMaterial()->isStateActive();
613 tf::Taskflow taskflow;
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];
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);
624 if (pi->d_material_p->isStateActive()) {
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];
632 const auto &m = this->d_mX[i];
636 auto check_up = horizon + 0.5 * mesh_size;
637 auto check_low = horizon - 0.5 * mesh_size;
640 for (size_t j : this->d_neighPd[i]) {
642 const auto &xj = this->d_xRef[j];
643 const auto &uj = this->d_u[j];
644 double rji = (xj - xi).length();
646 double change_length = (xj - xi + uj - ui).length() - rji;
649 double s = change_length / rji;
650 double sc = pi->d_material_p->getSc(rji);
653 auto fs = this->d_fracture_p->getBondState(i, k);
654 if (!fs && util::isGreater(std::abs(s), sc + 1.0e-10))
656 this->d_fracture_p->setBondState(i, k, fs);
661 auto volj = this->d_vol[j];
663 if (util::isGreater(rji, check_low))
664 volj *= (check_up - rji) / mesh_size;
666 theta += rji * change_length * pi->d_material_p->getInfFn(rji) *
673 this->d_thetaX[i] = 3. * theta / m;
678 executor.run(taskflow).get();
683 tf::Taskflow taskflow;
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];
690 util::Point force_i = util::Point();
691 double scalar_f = 0.;
696 const auto rhoi = this->getDensity(i);
697 const auto &ptIdi = this->getPtId(i);
698 auto &pi = this->getParticleFromAllList(ptIdi);
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];
708 auto check_up = horizon + 0.5 * mesh_size;
709 auto check_low = horizon - 0.5 * mesh_size;
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);
723 const auto &mj = this->d_mX[j];
724 const auto &thetaj = this->d_thetaX[j];
727 if (util::isGreater(rji, check_low))
728 volj *= (check_up - rji) / mesh_size;
731 if (pi->d_material_p->isStateActive()) {
734 pi->d_material_p->getBondEF(rji, Sji, fs, mi, thetai);
736 pi->d_material_p->getBondEF(rji, Sji, fs, mj, thetaj);
739 scalar_f = (ef_i.second + ef_j.second) * volj;
741 force_i += scalar_f * pi->d_material_p->getBondForceDirection(
747 bool break_bonds = true;
750 pi->d_material_p->getBondEF(rji, Sji, fs, break_bonds);
751 this->d_fracture_p->setBondState(i, k, fs);
754 scalar_f = ef.second * volj;
756 force_i += scalar_f * pi->d_material_p->getBondForceDirection(
762 auto yji = xj + uj - (xi + ui);
763 auto Rji = yji.length();
764 scalar_f = pi->d_Kn * volj * (Rji - pi->d_Rc) / Rji;
767 force_i += scalar_f * yji;
771 auto Sc = pi->d_material_p->getSc(rji);
772 if (util::isGreater(std::abs(Sji / Sc), Zi))
773 Zi = std::abs(Sji / Sc);
782 this->d_f[i] = force_i;
788 executor.run(taskflow).get();
793 log(
" Computing external force \n", 3);
795 auto gravity = d_pDeck_p->d_gravity;
797 if (gravity.length() > 1.0E-8) {
799 tf::Taskflow taskflow;
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;
806 executor.run(taskflow).get();
810 for (
auto &p : d_particlesListTypeAll)
811 d_fLoading_p->apply(d_time, p);
815 log(
" Computing external displacement bc \n", 3);
816 for (
auto &p : d_particlesListTypeAll)
817 d_uLoading_p->apply(d_time, p);
822 log(
" Computing normal contact force \n", 3);
830 tf::Taskflow taskflow;
832 taskflow.for_each_index((std::size_t) 0,
833 d_fContCompNodes.size(),
835 [
this](std::size_t II) {
837 auto i = this->d_fContCompNodes[II];
840 util::Point force_i = util::Point();
841 double scalar_f = 0.;
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;
849 double rhoi = pi->getDensity();
851 const auto &yi = this->d_x[i];
852 const auto &ui = this->d_u[i];
853 const auto &vi = this->d_v[i];
854 const auto &voli = this->d_vol[i];
856 const std::vector<size_t> &neighs = this->d_neighC[i];
858 if (neighs.size() > 0) {
860 for (const auto &j_id: neighs) {
863 const auto &yj = this->d_x[j_id];
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();
870 (pi->getTypeIndex() == 1 and pj->getTypeIndex() == 1);
873 if (ptIdj != ptIdi && !both_walls) {
876 const auto &contact =
877 d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
879 if (util::isLess(Rji, contact.d_contactR)) {
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;
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();
895 scalar_f = contact.d_Kn * (Rji - contact.d_contactR) *
899 force_i += scalar_f * en;
902 force_i += contact.d_mu * scalar_f * et;
906 bool node_lvl_damp = false;
912 auto meq = util::equivalentMass(rhoi * voli, rhoj * volj);
915 std::sqrt(contact.d_kappa * contact.d_contactR * meq);
917 auto &pii = this->d_particlesListTypeAll[pi->getId()];
918 vji = this->d_v[j_id] - pii->getVCenter();
922 force_i += beta_n * vn_mag * en / voli;
930 this->d_f[i] += force_i;
934 executor.run(taskflow).get();
938 log(
" Computing normal damping force \n", 3);
939 for (
auto &pi : d_particlesListTypeParticle) {
941 auto pi_id = pi->getId();
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();
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();
957 const auto &contact = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
959 if (
util::isLess(dist_xcji, Rj + Ri + 1.01 * contact.d_contactR)) {
961 auto vol_pj = M_PI * Rj * Rj;
962 auto rhoj = pj->getDensity();
967 auto beta_n = contact.d_betan *
968 std::sqrt(contact.d_kappa * contact.d_contactR * meq);
973 hat_xc_ji = xc_ji / dist_xcji;
978 auto vc_ji = pj->getVCenter() - pi_vc;
979 auto vc_mag = vc_ji * hat_xc_ji;
984 force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
995 d_neighWallNodesCondensed[pi->getId()].clear();
997 for (
size_t j=0; j<d_neighWallNodes[pi_id].size(); j++) {
999 const auto &j_id = pi->getNodeId(j);
1000 const auto &yj = this->d_x[j_id];
1002 for (
size_t k=0; k<d_neighWallNodes[pi_id][j].size(); k++) {
1004 const auto &k_id = d_neighWallNodes[pi_id][j][k];
1005 const auto &pk = d_particlesListTypeAll[d_ptId[k_id]];
1007 double Rjk = (this->d_x[k_id] - yj).length();
1009 const auto &contact =
1010 d_cDeck_p->getContact(pi->d_zoneId, pk->d_zoneId);
1020 for (
auto &j : d_neighWallNodesCondensed[pi_id]) {
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;
1030 = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);
1033 auto beta_n = contact.d_betan *
1034 std::sqrt(contact.d_kappa * contact.d_contactR * meq);
1037 auto xc_ji = this->d_x[j] - pi_xc;
1040 hat_xc_ji = xc_ji / xc_ji.length();
1043 auto vc_ji = this->d_v[j] - pi_vc;
1044 auto vc_mag = vc_ji * hat_xc_ji;
1049 force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
1055 tf::Taskflow taskflow;
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;
1063 executor.run(taskflow).get();
1070 log(
"Applying initial condition \n", 3);
1072 if (!d_pDeck_p->d_icDeck.d_icActive)
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;
1080 tf::Taskflow taskflow;
1082 taskflow.for_each_index((std::size_t) 0,
1085 [
this, ic_v, ic_p_list](std::size_t i) {
1086 auto &p = this->d_particlesListTypeAll[ic_p_list[i]];
1089 for (size_t j = 0; j < p->getNumNodes(); j++)
1090 p->setVLocal(j, ic_v);
1094 executor.run(taskflow).get();
1099 d_particlesListTypeParticle.resize(0);
1100 d_particlesListTypeAll.resize(0);
1101 d_particlesListTypeWall.resize(0);
1102 d_referenceParticles.clear();
1105 for (
size_t z = 0; z < d_pDeck_p->d_particleZones.size(); z++) {
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);
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);
1122 auto psize = d_particlesListTypeAll.size();
1125 auto &pz = d_pDeck_p->d_particleZones[z];
1128 auto z_id = pz.d_zone.d_zoneId;
1130 std::cerr << fmt::format(
"Error: d_zoneId = {} in ParticleZone for "
1131 "z = {} should be equal to z = {}.\n",
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);
1144 const auto &geomData = pz.d_meshDeck.d_createMeshGeomData;
1145 if (pz.d_meshDeck.d_createMeshInfo ==
"uniform"
1146 and geomData.d_geomName ==
"rectangle") {
1149 std::pair<std::vector<double>, std::vector<double>> box;
1150 std::vector<size_t> nGrid(3, 0);
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]);
1156 nGrid[i] = size_t((geomData.d_geomParams[i+3] - geomData.d_geomParams[i])/pz.d_meshDeck.d_h);
1158 std::cout << fmt::format(
"box.first[i] = {}, "
1159 "box.second[i] = {}, "
1167 d_modelDeck_p->
d_dim,
1170 mesh = std::make_shared<fe::Mesh>(temp_mesh);
1173 std::cerr <<
"Error: Currently, we can only support in-built uniform mesh for rectangles.\n";
1179 log(d_name +
": Creating reference particle in zone = " +
1180 std::to_string(z_id) +
"\n");
1183 auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
1184 auto rep_geom_params = pz.d_particleGeomData.d_geomParams;
1186 auto ref_p = std::make_shared<particle::RefParticle>(
1187 d_referenceParticles.size(),
1188 static_cast<std::shared_ptr<ModelData>
>(
this),
1192 d_referenceParticles.emplace_back(ref_p);
1195 log(d_name +
": Creating particles in zone = " +
1196 std::to_string(z_id) +
"\n");
1198 if (pz.d_genMethod ==
"From_File") {
1199 createParticlesFromFile(z, ref_p);
1202 if (pz.d_createParticleUsingParticleZoneGeomObject or !d_input_p->isMultiParticle()) {
1203 createParticleUsingParticleZoneGeomObject(z, ref_p);
1206 std::cerr <<
"Error: Particle generation method = " << pz.d_genMethod <<
1207 " not recognized.\n";
1213 auto psize_new = d_particlesListTypeAll.size();
1216 d_zInfo.emplace_back(std::vector<size_t>{psize, psize_new, z_id});
1222 std::shared_ptr<particle::RefParticle> ref_p) {
1224 log(d_name +
": Creating particle using Particle Zone Geometry Object\n", 1);
1227 auto &pz = d_pDeck_p->d_particleZones[z];
1230 auto z_id = pz.d_zone.d_zoneId;
1238 pz.d_isWall ?
"wall" :
"particle",
1239 d_particlesListTypeAll.size(),
1240 pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
1242 ref_p->getDimension(),
1243 pz.d_particleDescription,
1245 pz.d_allDofsConstrained,
1246 ref_p->getNumNodes(),
1248 static_cast<std::shared_ptr<ModelData>
>(
this),
1258 d_particlesListTypeWall.push_back(p);
1260 d_particlesListTypeParticle.push_back(p);
1262 d_particlesListTypeAll.push_back(p);
1267 size_t z, std::shared_ptr<particle::RefParticle> ref_p) {
1269 log(d_name +
": Creating particle from file\n", 1);
1272 auto &pz = d_pDeck_p->d_particleZones[z];
1275 auto z_id = pz.d_zone.d_zoneId;
1279 std::vector<util::Point> centers;
1280 std::vector<double> rads;
1281 std::vector<double> orients;
1282 if (pz.d_particleFileDataType ==
"loc_rad") {
1284 ¢ers, &rads, z_id);
1287 0., 1., d_modelDeck_p->d_seed);
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);
1293 for (
size_t i = 0; i < rads.size(); i++)
1298 else if (pz.d_particleFileDataType ==
"loc_rad_orient") {
1300 d_modelDeck_p->d_dim, ¢ers,
1301 &rads, &orients, z_id);
1304 log(fmt::format(
"zone_id: {}, rads: {}, orients: {}, centers: {} \n", z_id,
1309 const auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
1310 auto rep_geom_params = pz.d_particleGeomData.d_geomParams;
1313 std::pair<util::Point, util::Point> box = rep_geom_p->box();
1315 size_t p_counter = 0;
1316 size_t p_old_size = d_particlesListTypeAll.size();
1317 for (
const auto &site : centers) {
1319 double particle_radius = rads[p_counter];
1320 double particle_orient = orients[p_counter];
1323 std::shared_ptr<util::geometry::GeomObject> p_geom;
1324 auto scale = createGeometryAtSite(particle_radius,
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();
1348 pz.d_isWall ?
"wall" :
"particle",
1349 d_particlesListTypeAll.size(),
1350 pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
1352 ref_p->getDimension(),
1353 pz.d_particleDescription,
1355 pz.d_allDofsConstrained,
1356 ref_p->getNumNodes(),
1358 static_cast<std::shared_ptr<ModelData>
>(
this),
1368 d_particlesListTypeWall.push_back(p);
1370 d_particlesListTypeParticle.push_back(p);
1372 d_particlesListTypeAll.push_back(p);
1378 const double &particle_orient,
1380 const std::vector<double> &rep_geom_params,
1381 const std::shared_ptr<util::geometry::GeomObject> &rep_geom_p,
1382 std::shared_ptr<util::geometry::GeomObject> &p_geom) {
1383 std::vector<double> params;
1384 for (
auto x : rep_geom_params)
1385 params.push_back(x);
1394 {
"circle",
"sphere",
"hexagon",
1395 "triangle",
"square",
"cube"})) {
1400 size_t num_params = 4;
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];
1408 scale = particle_radius / rep_geom_params[0];
1410 else if (rep_geom_p->d_name ==
"drum2d") {
1413 size_t num_params = 5;
1415 if (params.size() < num_params)
1416 params.resize(num_params);
1418 params[0] = particle_radius;
1419 params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0];
1420 for (
int dof = 0; dof < 3; dof++)
1421 params[dof + 2] = site[dof];
1423 scale = params[0] / rep_geom_params[0];
1425 else if (rep_geom_p->d_name ==
"rectangle") {
1428 size_t num_params = 5;
1430 if (params.size() < num_params)
1431 params.resize(num_params);
1433 params[0] = particle_radius;
1434 params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0];
1435 for (
int dof = 0; dof < 3; dof++)
1436 params[dof + 2] = site[dof];
1438 scale = params[0] / rep_geom_params[0];
1440 else if (rep_geom_p->d_name ==
"cuboid") {
1444 if (params.size() < 6)
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];
1453 const auto Lx = particle_radius;
1454 const auto Ly = particle_radius * ref_Ly / ref_Lx;
1455 const auto Lz = particle_radius * ref_Lz / ref_Lx;
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;
1464 scale = params[0] / rep_geom_params[0];
1467 std::cerr << fmt::format(
"Error: PeriDEM supports following type "
1468 "of geometries for particles = {}\n",
1474 std::vector<std::string> vec_geom_type;
1475 std::vector<std::string> vec_geom_flag;
1477 vec_geom_flag, p_geom,
1478 d_modelDeck_p->d_dim,
false);
1485 for (
auto &p: d_particlesListTypeAll) {
1486 if (p->d_geom_p->d_name ==
"null" or
1489 auto bbox = p->getMeshP()->getBoundingBox();
1491 std::string geom_name =
"rectangle";
1492 if ( p->getMeshP()->getDimension() == 3)
1493 geom_name =
"cuboid";
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];
1501 std::vector<std::string> vec_type;
1502 std::vector<std::string> vec_flag;
1509 p->getMeshP()->getDimension());
1518 for (
const auto *p : d_particlesListTypeAll) {
1520 auto h = p->getMeshSize();
1533 log(fmt::format(
"{}: Contact setup\n hmin = {:.6f}, hmax = {:.6f} \n",
1534 d_name, d_hMin, d_hMax), 1);
1538 std::vector<double> bulk_modulus;
1540 for (
size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
1542 double kappa_i = d_pDeck_p->d_particleZones[i].d_matDeck.d_matData.d_K;
1545 std::cerr <<
"Error: We need bulk modulus provided in input file.\n";
1546 std::cerr << d_pDeck_p->d_particleZones[i].printStr();
1550 bulk_modulus.push_back(kappa_i);
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++) {
1561 if (d_maxContactR < deck->d_contactR)
1571 double log_e = std::log(deck->
d_eps);
1574 (-2. * log_e * std::sqrt(1. / (M_PI * M_PI + log_e * log_e)));
1576 log(fmt::format(
" contact_radius = {:.6f}, hmin = {:.6f}, Kn = {:5.3e}, "
1578 "betan = {:7.5f}, mu = {:.4f}, kappa = {:5.3e}\n",
1588 or d_modelDeck_p->d_populateElementNodeConnectivity) {
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);
1600 size_t totalQuadPoints = 0;
1601 for (
auto &p: d_particlesListTypeAll) {
1602 const auto &particle_mesh_p = p->getMeshP();
1609 elem =
new fe::TriElem(d_modelDeck_p->d_quadOrder);
1613 elem =
new fe::TetElem(d_modelDeck_p->d_quadOrder);
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());
1620 p->d_globQuadStart = totalQuadPoints;
1621 totalQuadPoints += particle_mesh_p->getNumElements() *
1623 p->d_globQuadEnd = totalQuadPoints;
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() *
1636 d_xQuadCur.resize(totalQuadPoints);
1637 d_strain.resize(totalQuadPoints);
1638 d_stress.resize(totalQuadPoints);
1644 d_neighPd.resize(d_x.size());
1646 auto t1 = steady_clock::now();
1649 tf::Taskflow taskflow;
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();
1655 std::vector<size_t> neighs;
1656 std::vector<double> sqr_dist;
1657 if (this->d_nsearch_p->radiusSearchIncludeTag(this->d_x[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]));
1672 executor.run(taskflow).get();
1674 auto t2 = steady_clock::now();
1675 log(fmt::format(
"{}: Peridynamics neighbor update time = {}\n",
1681 auto update = updateContactNeighborSearchParameters();
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);
1694 if (d_neighC.size() != d_x.size())
1695 d_neighC.resize(d_x.size());
1698 tf::Taskflow taskflow;
1700 taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1,
1701 [
this](std::size_t i) {
1703 const auto &pi = this->d_ptId[i];
1704 const auto &pi_particle = this->d_particlesListTypeAll[pi];
1707 bool perform_search_based_on_particle = true;
1708 if (pi_particle->d_typeIndex == 1)
1709 perform_search_based_on_particle = false;
1711 if (pi_particle->d_allDofsConstrained or !pi_particle->d_computeForce)
1712 perform_search_based_on_particle = false;
1714 if (perform_search_based_on_particle) {
1716 std::vector<size_t> neighs;
1717 std::vector<double> sqr_dist;
1719 this->d_neighC[i].clear();
1721 auto n = this->d_nsearch_p->radiusSearchExcludeTag(
1723 this->d_contNeighSearchRadius,
1730 for (auto neigh: neighs) {
1732 this->d_neighC[i].push_back(neigh);
1739 executor.run(taskflow).get();
1743 d_neighWallNodes.resize(d_particlesListTypeAll.size());
1744 d_neighWallNodesDistance.resize(d_particlesListTypeAll.size());
1745 d_neighWallNodesCondensed.resize(d_particlesListTypeAll.size());
1747 for (
auto &pi : d_particlesListTypeParticle) {
1749 d_neighWallNodes[pi->getId()].resize(pi->getNumNodes());
1750 d_neighWallNodesDistance[pi->getId()].resize(pi->getNumNodes());
1755 tf::Taskflow taskflow;
1757 taskflow.for_each_index((std::size_t) 0,
1760 [
this, &pi](std::size_t i) {
1762 auto i_glob = pi->getNodeId(i);
1763 auto yi = this->d_x[i_glob];
1765 const std::vector<size_t> &neighs = this->d_neighC[i_glob];
1767 this->d_neighWallNodes[pi->getId()][i].clear();
1768 this->d_neighWallNodesDistance[pi->getId()][i].clear();
1770 for (const auto &j_id: neighs) {
1772 auto &ptIdj = this->d_ptId[j_id];
1773 auto &pj = this->getParticleFromAllList(
1777 if (pj->getTypeIndex() == 1) {
1778 this->d_neighWallNodes[pi->getId()][i].push_back(j_id);
1785 executor.run(taskflow).get();
1794 if (d_contNeighUpdateInterval == 0 and
1796 d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
1797 d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
1798 d_contNeighSearchRadius = d_maxContactR * d_pDeck_p->d_pNeighDeck.d_sFactor;
1805 appendKeyData(
"update_contact_neigh_search_params_init_call_count", 1);
1807 if (
int(getKeyData(
"update_contact_neigh_search_params_init_call_count")) == 1)
1810 if (
int(getKeyData(
"update_contact_neigh_search_params_init_call_count")) == 2) {
1811 d_contNeighTimestepCounter++;
1812 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
1817 if (d_modelDeck_p->d_isRestartActive and d_n == d_restartDeck_p->d_step) {
1819 d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
1822 if (d_contNeighUpdateInterval == 1) {
1824 d_contNeighSearchRadius = d_maxContactR;
1827 d_contNeighTimestepCounter++;
1828 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
1834 size_t update_param_interval =
1835 d_contNeighUpdateInterval > 5 ? size_t(
1836 0.2 * d_contNeighUpdateInterval) : 1;
1839 if (d_contNeighTimestepCounter > 0 and d_contNeighTimestepCounter % update_param_interval != 0) {
1841 d_contNeighTimestepCounter++;
1842 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
1846 for (
auto &pi : d_particlesListTypeAll) {
1848 pi->d_globStart, pi->d_globEnd);
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"
1856 <<
"\n\n Magnitude of velocity = "
1857 << d_vMag[max_v_node] <<
"\n";
1861 d_maxVelocityParticlesListTypeAll[pi->getId()]
1862 = d_vMag[max_v_node];
1869 auto up_interval_old = d_contNeighUpdateInterval;
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
1883 d_contNeighUpdateInterval = size_t(d_maxContactR/(d_maxVelocity * d_currentDt));
1884 if (up_interval_old > d_contNeighUpdateInterval) {
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);
1895 d_contNeighSearchRadius = max_search_r_from_contact_R;
1898 d_contNeighTimestepCounter = 0;
1900 if (d_contNeighUpdateInterval < 1) {
1901 d_contNeighUpdateInterval = 1;
1902 d_contNeighSearchRadius = d_maxContactR;
1907 d_contNeighSearchRadius = d_contNeighUpdateInterval < 2 ? d_maxContactR : max_search_r_from_contact_R;
1910 if (up_interval_old > d_contNeighUpdateInterval) {
1911 log(fmt::format(
" Contact neighbor parameters: \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",
1922 "contact neighbor update interval",
1923 d_contNeighUpdateInterval,
1924 "contact neighbor update time step counter",
1925 d_contNeighTimestepCounter,
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);
1936 d_contNeighTimestepCounter++;
1937 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
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",
1957 log(fmt::format(
"{}: Output step = {}, time = {:.6f} \n",
1958 d_name, d_n, d_time),
1961 if (d_outputDeck_p->d_debug > 0 and getKeyData(
"debug_once") < 0) {
1963 setKeyData(
"debug_once", 1);
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) {
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;
1988 oss << tabS <<
"h_min = " << d_hMin <<
", h_max = " << d_hMax << std::endl;
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);
1998 out_filename = out_filename + d_outputDeck_p->d_tagPPFile +
"_" + std::to_string(d_n / dt_out);
2001 if (d_outputDeck_p->d_performFEOut)
2002 writer.appendMesh(
this, d_outputDeck_p->d_outTags);
2004 writer.appendNodes(
this, d_outputDeck_p->d_outTags);
2006 writer.addTimeStep(d_time);
2014 if (d_particlesMatDataList.empty()) {
2015 for (
auto &p: d_particlesListTypeAll) {
2016 d_particlesMatDataList.push_back(p->getMaterial()->computeMaterialProperties(
2017 p->getMeshP()->getDimension()));
2021 for (
auto &p: d_particlesListTypeAll) {
2023 const auto particle_mesh_p = p->getMeshP();
2028 d_modelDeck_p->d_quadOrder);
2030 auto p_z_id = p->d_zoneId;
2031 auto isPlaneStrain = d_pDeck_p->d_particleZones[p_z_id].d_matDeck.d_isPlaneStrain;
2037 d_particlesMatDataList[p->getId()].d_nu,
2038 d_particlesMatDataList[p->getId()].d_lambda,
2039 d_particlesMatDataList[p->getId()].d_mu,
2041 d_modelDeck_p->d_quadOrder);
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);
2049 out_filename = out_filename + d_outputDeck_p->d_tagPPFile +
"_" + std::to_string(d_n / dt_out);
2052 writer1.appendStrainStress(
this);
2053 writer1.addTimeStep(d_time);
2059 d_outputDeck_p->d_outTags)) {
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";
2065 out_filename = out_filename + d_outputDeck_p->d_tagPPFile
2066 +
"_" + std::to_string(d_n / dt_out) +
".csv";
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";
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))
2090 const auto &p0 = this->d_particlesListTypeAll[0];
2091 const auto &p1 = this->d_particlesListTypeAll[1];
2094 const auto &xc0 = p0->getXCenter();
2095 const auto &xc1 = p1->getXCenter();
2096 const double &r = p0->d_geom_p->boundingRadius();
2098 const auto &contact = d_cDeck_p->getContact(p0->d_zoneId, p1->d_zoneId);
2099 double r_e = r + contact.d_contactR;
2101 double pen_dist = xc1.dist(xc0) - r_e - r;
2102 double contact_area_radius = 0.;
2104 contact_area_radius =
2105 std::sqrt(std::pow(r_e, 2.) - std::pow(r_e + pen_dist, 2.));
2108 contact_area_radius = 0.;
2112 double max_dist = xc1.d_y + p1->d_geom_p->boundingRadius();
2115 double max_y_loc = p1->getXLocal(0).d_y;
2117 for (
size_t i = 0; i < p1->getNumNodes(); i++)
2119 max_y_loc = p1->getXLocal(i).d_y;
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);
2131 return fmt::format(
" Post-processing: max y = {:.6f} \n", max_y);
2136 if (d_outputDeck_p->d_outCriteria ==
"max_particle_dist" &&
2137 d_pDeck_p->d_testName ==
"two_particle") {
2141 const auto &xci = d_particlesListTypeAll[0]->getXCenter();
2142 const auto &xcj = d_particlesListTypeAll[1]->getXCenter();
2146 d_outputDeck_p->d_outCriteriaParams[0])) {
2148 if(d_ppFile.is_open())
2153 else if (d_outputDeck_p->d_outCriteria ==
"max_node_dist") {
2163 auto max_x = d_x[max_pt_and_index.second];
2167 d_outputDeck_p->d_outCriteriaParams[0])) {
2170 if(d_ppFile.is_open())
2173 log(fmt::format(
"{}: Terminating simulation as one of the failing"
2174 " criteria is met. Point ({:.6f}, {:.6f}, {:.6f}) is at "
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]));
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))
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];
2200 auto dx = wall->getXLocal(0) - wall->getXRefLocal(0);
2201 double wall_penetration = dx[f_dir];
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);
2210 bool use_static_file =
true;
2211 if (use_static_file) {
2212 if (!d_ppFile.is_open()) {
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);
2220 d_ppFile <<
"t, delta, force \n";
2223 d_ppFile << fmt::format(
"%4.6e, %4.6e, %4.6e\n", d_time, wall_penetration,
2224 tot_reaction_force);
2227 setKeyData(
"wall_penetration", wall_penetration);
2228 setKeyData(
"tot_reaction_force", tot_reaction_force);
2230 return fmt::format(
" Post-processing: wall penetration = {:"
2232 "reaction force = {:5.3e} \n",
2233 wall_penetration, tot_reaction_force);
A base class which provides methods to map points to/from reference element and to compute quadrature...
size_t getNumQuadPoints()
Get number of quadrature points in the data.
A class for mapping and quadrature related operations for linear 2-node line element.
size_t d_dim
Dimension of the mesh.
A class for mapping and quadrature related operations for bi-linear quadrangle element.
A class for mapping and quadrature related operations for linear tetrahedron element.
A class for mapping and quadrature related operations for linear triangle element.
virtual void computeExternalForces()
Computes external/boundary condition forces.
DEMModel(inp::Input *deck, std::string modelName="DEMModel")
Constructor.
virtual void integrate()
Perform time integration.
virtual void setupContact()
Creates particles in a given container.
std::string ppCompressiveTest()
Function that handles post-processing for compressive test of particulate media by rigid wall and ret...
virtual void updatePeridynamicNeighborlist()
Update neighborlist for peridynamics force.
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.
virtual void integrateStep()
Performs one time step.
virtual void computePeridynamicForces()
Computes peridynamic forces.
virtual void updateContactNeighborlist()
Update neighborlist for contact.
virtual void computeContactForces()
Computes contact forces.
virtual void close()
Closure operations.
virtual void updateGeometryObjectsPostInit()
Update varioud geometry objects associated with container, particles, and reference particles.
virtual void updateNeighborlistCombine()
Update neighborlist for contact and peridynamics force.
std::string ppTwoParticleTest()
Function that handles post-processing for two particle collision test and returns maximum vertical di...
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.
virtual void integrateVerlet()
Perform time integration using velocity verlet scheme.
virtual void createParticlesFromFile(size_t z, std::shared_ptr< particle::RefParticle > ref_p)
Creates particles in a Hexagonal arrangement.
virtual void init()
Initialize remaining data members.
virtual void computeForces()
Computes peridynamic forces and contact forces.
virtual void checkStop()
Checks if simulation should be stopped due to abnormal state of system.
virtual void restart(inp::Input *deck)
Restarts the simulation from previous state.
virtual void integrateCD()
Perform time integration using central-difference scheme.
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 bool updateContactNeighborSearchParameters()
Update contact neighbor search parameters.
virtual void createParticles()
Creates particles in a given container.
virtual void run(inp::Input *deck)
Main driver to simulate.
virtual void createParticleUsingParticleZoneGeomObject(size_t z, std::shared_ptr< particle::RefParticle > ref_p)
Creates particles in a given container.
virtual void setupQuadratureData()
Sets up quadrature data.
A class to store model data.
std::shared_ptr< inp::OutputDeck > d_outputDeck_p
Output deck.
A class to store particle geometry, nodal discretization, and methods.
A vtk writer for simple point data and complex fem mesh data.
A vtk writer for simple point data and complex fem mesh data.
Templated probability distribution.
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.
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.
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.
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...
void computeStateMx(model::ModelData *model, bool compute_in_parallel=false)
Computes the moment term in state-based peridynamic formulation.
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...
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.
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 > ¶ms, 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 getTabS(int nt)
Returns tab spaces of given size.
std::string printStr(const T &msg, int nt=print_default_tab)
Returns formatted string for output.
void initLogger(int debug_level=logger_default_debug_lvl, std::string filename="")
Initializes the logger.
void log(std::ostringstream &oss, bool screen_out=false, int printMpiRank=print_default_mpi_rank)
Global method to log the message.
std::pair< double, size_t > maxLengthAndMaxLengthIndex(const std::vector< util::Point > &data)
Returns the maximum length of point and index from list of points.
T max(const std::vector< T > &data)
Returns the maximum from list of data.
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.
bool isTagInList(const std::string &tag, const std::vector< std::string > &tags)
Returns true if tag is found in the list of tags.
bool isFree(const int &i, const unsigned int &dof)
Returns true if degree of freedom is free.
void addToList(const T &i, std::vector< T > &list)
Add element to the list.
size_t maxIndex(const std::vector< T > &data)
Returns the index corresponding to maximum from list of data.
unsigned int getNThreads()
Get number of threads to be used by taskflow.
bool isGreater(const double &a, const double &b)
Returns true if a > b.
double equivalentMass(const double &m1, const double &m2)
Compute harmonic mean of m1 and m2.
double transform_to_uniform_dist(double min, double max, double sample)
Transform sample from U(0,1) to U(a,b)
bool isLess(const double &a, const double &b)
Returns true if a < b.
A struct that stores transformation parameters and provides method to transform the particle....
A structure to represent 3d vectors.