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];
1466 }
else if (rep_geom_p->d_name ==
"ellipse") {
1469 size_t num_params = 6;
1471 if (params.size() < num_params)
1472 params.resize(num_params);
1474 params[0] = particle_radius;
1475 params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0];
1476 for (
int dof = 0; dof < 3; dof++)
1477 params[dof + 2] = site[dof];
1479 params[5] = particle_orient;
1481 scale = params[0] / rep_geom_params[0];
1484 std::cerr << fmt::format(
"Error: PeriDEM supports following type "
1485 "of geometries for particles = {}\n",
1491 std::vector<std::string> vec_geom_type;
1492 std::vector<std::string> vec_geom_flag;
1494 vec_geom_flag, p_geom,
1495 d_modelDeck_p->d_dim,
false);
1502 for (
auto &p: d_particlesListTypeAll) {
1503 if (p->d_geom_p->d_name ==
"null" or
1506 auto bbox = p->getMeshP()->getBoundingBox();
1508 std::string geom_name =
"rectangle";
1509 if ( p->getMeshP()->getDimension() == 3)
1510 geom_name =
"cuboid";
1512 std::vector<double> geom_params(6, 0.);
1513 for (
size_t i=0; i<3; i++) {
1514 geom_params[i] = bbox.first[i];
1515 geom_params[i+3] = bbox.second[i];
1518 std::vector<std::string> vec_type;
1519 std::vector<std::string> vec_flag;
1526 p->getMeshP()->getDimension());
1535 for (
const auto *p : d_particlesListTypeAll) {
1537 auto h = p->getMeshSize();
1550 log(fmt::format(
"{}: Contact setup\n hmin = {:.6f}, hmax = {:.6f} \n",
1551 d_name, d_hMin, d_hMax), 1);
1555 std::vector<double> bulk_modulus;
1557 for (
size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
1559 double kappa_i = d_pDeck_p->d_particleZones[i].d_matDeck.d_matData.d_K;
1562 std::cerr <<
"Error: We need bulk modulus provided in input file.\n";
1563 std::cerr << d_pDeck_p->d_particleZones[i].printStr();
1567 bulk_modulus.push_back(kappa_i);
1570 for (
size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
1571 for (
size_t j = 0; j < d_cDeck_p->d_data.size(); j++) {
1578 if (d_maxContactR < deck->d_contactR)
1588 double log_e = std::log(deck->
d_eps);
1591 (-2. * log_e * std::sqrt(1. / (M_PI * M_PI + log_e * log_e)));
1593 log(fmt::format(
" contact_radius = {:.6f}, hmin = {:.6f}, Kn = {:5.3e}, "
1595 "betan = {:7.5f}, mu = {:.4f}, kappa = {:5.3e}\n",
1605 or d_modelDeck_p->d_populateElementNodeConnectivity) {
1608 for (
auto &p: d_referenceParticles) {
1609 auto &particle_mesh_p = p->getMeshP();
1610 if (!particle_mesh_p->d_encDataPopulated && particle_mesh_p->d_enc.empty()) {
1611 particle_mesh_p->readElementData(particle_mesh_p->d_filename);
1617 size_t totalQuadPoints = 0;
1618 for (
auto &p: d_particlesListTypeAll) {
1619 const auto &particle_mesh_p = p->getMeshP();
1626 elem =
new fe::TriElem(d_modelDeck_p->d_quadOrder);
1630 elem =
new fe::TetElem(d_modelDeck_p->d_quadOrder);
1632 std::cerr << fmt::format(
"Error: Can not compute strain/stress as the element "
1633 "type = {} is not yet supported in this routine.\n", particle_mesh_p->getElementType());
1637 p->d_globQuadStart = totalQuadPoints;
1638 totalQuadPoints += particle_mesh_p->getNumElements() *
1640 p->d_globQuadEnd = totalQuadPoints;
1642 std::cout << fmt::format(
"p->id() = {}, "
1643 "p->d_globQuadStart = {}, "
1644 "totalQuadPoints = {}, "
1645 "p->d_globQuadEnd = {}",
1646 p->getId(), p->d_globQuadStart,
1647 particle_mesh_p->getNumElements() *
1653 d_xQuadCur.resize(totalQuadPoints);
1654 d_strain.resize(totalQuadPoints);
1655 d_stress.resize(totalQuadPoints);
1661 d_neighPd.resize(d_x.size());
1663 auto t1 = steady_clock::now();
1666 tf::Taskflow taskflow;
1668 taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1, [
this](std::size_t i) {
1669 const auto &pi = this->d_ptId[i];
1670 double search_r = this->d_particlesListTypeAll[pi]->d_material_p->getHorizon();
1672 std::vector<size_t> neighs;
1673 std::vector<double> sqr_dist;
1674 if (this->d_nsearch_p->radiusSearchIncludeTag(this->d_x[i],
1679 this->d_ptId) > 0) {
1680 for (std::size_t j = 0; j < neighs.size(); ++j)
1681 if (neighs[j] != i && this->d_ptId[neighs[j]] == pi) {
1682 this->d_neighPd[i].push_back(size_t(neighs[j]));
1689 executor.run(taskflow).get();
1691 auto t2 = steady_clock::now();
1692 log(fmt::format(
"{}: Peridynamics neighbor update time = {}\n",
1698 auto update = updateContactNeighborSearchParameters();
1706 auto pt_cloud_update_time = d_nsearch_p->setInputCloud();
1707 setKeyData(
"pt_cloud_update_time", pt_cloud_update_time);
1708 appendKeyData(
"tree_compute_time", pt_cloud_update_time);
1709 appendKeyData(
"avg_tree_update_time", pt_cloud_update_time/d_infoN);
1711 if (d_neighC.size() != d_x.size())
1712 d_neighC.resize(d_x.size());
1715 tf::Taskflow taskflow;
1717 taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1,
1718 [
this](std::size_t i) {
1720 const auto &pi = this->d_ptId[i];
1721 const auto &pi_particle = this->d_particlesListTypeAll[pi];
1724 bool perform_search_based_on_particle = true;
1725 if (pi_particle->d_typeIndex == 1)
1726 perform_search_based_on_particle = false;
1728 if (pi_particle->d_allDofsConstrained or !pi_particle->d_computeForce)
1729 perform_search_based_on_particle = false;
1731 if (perform_search_based_on_particle) {
1733 std::vector<size_t> neighs;
1734 std::vector<double> sqr_dist;
1736 this->d_neighC[i].clear();
1738 auto n = this->d_nsearch_p->radiusSearchExcludeTag(
1740 this->d_contNeighSearchRadius,
1747 for (auto neigh: neighs) {
1749 this->d_neighC[i].push_back(neigh);
1756 executor.run(taskflow).get();
1760 d_neighWallNodes.resize(d_particlesListTypeAll.size());
1761 d_neighWallNodesDistance.resize(d_particlesListTypeAll.size());
1762 d_neighWallNodesCondensed.resize(d_particlesListTypeAll.size());
1764 for (
auto &pi : d_particlesListTypeParticle) {
1766 d_neighWallNodes[pi->getId()].resize(pi->getNumNodes());
1767 d_neighWallNodesDistance[pi->getId()].resize(pi->getNumNodes());
1772 tf::Taskflow taskflow;
1774 taskflow.for_each_index((std::size_t) 0,
1777 [
this, &pi](std::size_t i) {
1779 auto i_glob = pi->getNodeId(i);
1780 auto yi = this->d_x[i_glob];
1782 const std::vector<size_t> &neighs = this->d_neighC[i_glob];
1784 this->d_neighWallNodes[pi->getId()][i].clear();
1785 this->d_neighWallNodesDistance[pi->getId()][i].clear();
1787 for (const auto &j_id: neighs) {
1789 auto &ptIdj = this->d_ptId[j_id];
1790 auto &pj = this->getParticleFromAllList(
1794 if (pj->getTypeIndex() == 1) {
1795 this->d_neighWallNodes[pi->getId()][i].push_back(j_id);
1802 executor.run(taskflow).get();
1811 if (d_contNeighUpdateInterval == 0 and
1813 d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
1814 d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
1815 d_contNeighSearchRadius = d_maxContactR * d_pDeck_p->d_pNeighDeck.d_sFactor;
1822 appendKeyData(
"update_contact_neigh_search_params_init_call_count", 1);
1824 if (
int(getKeyData(
"update_contact_neigh_search_params_init_call_count")) == 1)
1827 if (
int(getKeyData(
"update_contact_neigh_search_params_init_call_count")) == 2) {
1828 d_contNeighTimestepCounter++;
1829 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
1834 if (d_modelDeck_p->d_isRestartActive and d_n == d_restartDeck_p->d_step) {
1836 d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
1839 if (d_contNeighUpdateInterval == 1) {
1841 d_contNeighSearchRadius = d_maxContactR;
1844 d_contNeighTimestepCounter++;
1845 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
1851 size_t update_param_interval =
1852 d_contNeighUpdateInterval > 5 ? size_t(
1853 0.2 * d_contNeighUpdateInterval) : 1;
1856 if (d_contNeighTimestepCounter > 0 and d_contNeighTimestepCounter % update_param_interval != 0) {
1858 d_contNeighTimestepCounter++;
1859 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
1863 for (
auto &pi : d_particlesListTypeAll) {
1865 pi->d_globStart, pi->d_globEnd);
1867 if (max_v_node > pi->d_globEnd or max_v_node < pi->d_globStart) {
1868 std::cerr << fmt::format(
"Error: max_v_node = {} for "
1869 "particle of id = {} is not in the limit.\n",
1870 max_v_node, pi->getId())
1871 <<
"Particle info = \n"
1873 <<
"\n\n Magnitude of velocity = "
1874 << d_vMag[max_v_node] <<
"\n";
1878 d_maxVelocityParticlesListTypeAll[pi->getId()]
1879 = d_vMag[max_v_node];
1886 auto up_interval_old = d_contNeighUpdateInterval;
1891 double safety_factor = d_pDeck_p->d_pNeighDeck.d_sFactor > 5 ? d_pDeck_p->d_pNeighDeck.d_sFactor : 10;
1892 auto max_search_r_from_contact_R = d_pDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
1893 auto max_search_r = d_maxVelocity * d_currentDt
1894 * d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval
1900 d_contNeighUpdateInterval = size_t(d_maxContactR/(d_maxVelocity * d_currentDt));
1901 if (up_interval_old > d_contNeighUpdateInterval) {
1903 log(fmt::format(
"Warning: Contact search radius based on velocity is greater than "
1904 "the max contact radius.\n"
1905 "Warning: Adjusting contact neighborlist update interval.\n"
1906 "{:>13} = {:4.6e}, time step = {}, "
1907 "velocity-based r = {:4.6e}, max contact r = {:4.6e}\n",
1908 "Time", d_time, d_n, max_search_r, max_search_r_from_contact_R),
1909 2, d_n % d_infoN == 0, 3);
1912 d_contNeighSearchRadius = max_search_r_from_contact_R;
1915 d_contNeighTimestepCounter = 0;
1917 if (d_contNeighUpdateInterval < 1) {
1918 d_contNeighUpdateInterval = 1;
1919 d_contNeighSearchRadius = d_maxContactR;
1924 d_contNeighSearchRadius = d_contNeighUpdateInterval < 2 ? d_maxContactR : max_search_r_from_contact_R;
1927 if (up_interval_old > d_contNeighUpdateInterval) {
1928 log(fmt::format(
" Contact neighbor parameters: \n"
1932 " {:48s} = {:4.6e}\n"
1933 " {:48s} = {:4.6e}\n"
1934 " {:48s} = {:4.6e}\n"
1935 " {:48s} = {:4.6e}\n"
1936 " {:48s} = {:4.6e}\n"
1937 " {:48s} = {:4.6e}\n",
1939 "contact neighbor update interval",
1940 d_contNeighUpdateInterval,
1941 "contact neighbor update time step counter",
1942 d_contNeighTimestepCounter,
1943 "search radius", d_contNeighSearchRadius,
1944 "max contact radius", d_maxContactR,
1945 "search radius factor", d_pDeck_p->d_pNeighDeck.d_sFactor,
1946 "max search r from velocity", max_search_r,
1947 "max search r from contact r", max_search_r_from_contact_R,
1948 "max velocity", d_maxVelocity),
1949 2, d_n % d_infoN == 0, 3);
1953 d_contNeighTimestepCounter++;
1954 return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
1966 float p = float(d_n) * 100. / d_modelDeck_p->d_Nt;
1967 int m = std::max(1,
int(d_modelDeck_p->d_Nt / 10));
1968 if (d_n % m == 0 &&
int(p) > 0)
1969 log(fmt::format(
"{}: Simulation {}% complete\n",
1974 log(fmt::format(
"{}: Output step = {}, time = {:.6f} \n",
1975 d_name, d_n, d_time),
1978 if (d_outputDeck_p->d_debug > 0 and getKeyData(
"debug_once") < 0) {
1980 setKeyData(
"debug_once", 1);
1984 std::ostringstream oss;
1985 oss << tabS <<
"*******************************************\n";
1986 oss << tabS <<
"Debug various input decks\n\n\n";
1987 oss << d_modelDeck_p->printStr(nt + 1);
1988 oss << d_pDeck_p->printStr(nt + 1);
1989 oss << d_cDeck_p->printStr(nt + 1);
1990 oss << tabS <<
"\n\n*******************************************\n";
1991 oss << tabS <<
"Debug particle data\n\n\n";
1992 oss << tabS <<
"Number of particles = " << d_particlesListTypeAll.size() << std::endl;
1993 oss << tabS <<
"Number of particle zones = " << d_zInfo.size() << std::endl;
1994 for (
auto zone : d_zInfo) {
2000 oss << tabS <<
"Number of walls = " << d_particlesListTypeWall.size() << std::endl;
2001 for (
auto &d_wall : d_particlesListTypeWall)
2002 oss << tabS <<
"Number of nodes in wall in zone " << d_wall->d_zoneId
2003 <<
" is " << d_wall->getNumNodes() << std::endl;
2005 oss << tabS <<
"h_min = " << d_hMin <<
", h_max = " << d_hMax << std::endl;
2010 size_t dt_out = d_outputDeck_p->d_dtOutCriteria;
2011 std::string out_filename = d_outputDeck_p->d_path +
"output_";
2012 if (d_outputDeck_p->d_tagPPFile.empty())
2013 out_filename = out_filename + std::to_string(d_n / dt_out);
2015 out_filename = out_filename + d_outputDeck_p->d_tagPPFile +
"_" + std::to_string(d_n / dt_out);
2018 if (d_outputDeck_p->d_performFEOut)
2019 writer.appendMesh(
this, d_outputDeck_p->d_outTags);
2021 writer.appendNodes(
this, d_outputDeck_p->d_outTags);
2023 writer.addTimeStep(d_time);
2031 if (d_particlesMatDataList.empty()) {
2032 for (
auto &p: d_particlesListTypeAll) {
2033 d_particlesMatDataList.push_back(p->getMaterial()->computeMaterialProperties(
2034 p->getMeshP()->getDimension()));
2038 for (
auto &p: d_particlesListTypeAll) {
2040 const auto particle_mesh_p = p->getMeshP();
2045 d_modelDeck_p->d_quadOrder);
2047 auto p_z_id = p->d_zoneId;
2048 auto isPlaneStrain = d_pDeck_p->d_particleZones[p_z_id].d_matDeck.d_isPlaneStrain;
2054 d_particlesMatDataList[p->getId()].d_nu,
2055 d_particlesMatDataList[p->getId()].d_lambda,
2056 d_particlesMatDataList[p->getId()].d_mu,
2058 d_modelDeck_p->d_quadOrder);
2062 out_filename = d_outputDeck_p->d_path +
"output_strain_";
2063 if (d_outputDeck_p->d_tagPPFile.empty())
2064 out_filename = out_filename + std::to_string(d_n / dt_out);
2066 out_filename = out_filename + d_outputDeck_p->d_tagPPFile +
"_" + std::to_string(d_n / dt_out);
2069 writer1.appendStrainStress(
this);
2070 writer1.addTimeStep(d_time);
2076 d_outputDeck_p->d_outTags)) {
2078 out_filename = d_outputDeck_p->d_path +
"particle_locations_";
2079 if (d_outputDeck_p->d_tagPPFile.empty())
2080 out_filename = out_filename + std::to_string(d_n / dt_out) +
".csv";
2082 out_filename = out_filename + d_outputDeck_p->d_tagPPFile
2083 +
"_" + std::to_string(d_n / dt_out) +
".csv";
2085 std::ofstream oss(out_filename);
2086 oss <<
"i, x, y, z, r\n";
2087 for (
const auto &p : d_particlesListTypeAll) {
2088 auto xc = p->getXCenter();
2089 oss << p->d_zoneId <<
", " << xc.d_x <<
", " << xc.d_y <<
", " << xc.d_z
2090 <<
", " << p->d_geom_p->boundingRadius() <<
"\n";
2098 bool continue_dt =
false;
2099 auto check_dt = d_outputDeck_p->d_dtTestOut;
2100 if ((d_n % check_dt == 0) && (d_n >= check_dt))
2107 const auto &p0 = this->d_particlesListTypeAll[0];
2108 const auto &p1 = this->d_particlesListTypeAll[1];
2111 const auto &xc0 = p0->getXCenter();
2112 const auto &xc1 = p1->getXCenter();
2113 const double &r = p0->d_geom_p->boundingRadius();
2115 const auto &contact = d_cDeck_p->getContact(p0->d_zoneId, p1->d_zoneId);
2116 double r_e = r + contact.d_contactR;
2118 double pen_dist = xc1.dist(xc0) - r_e - r;
2119 double contact_area_radius = 0.;
2121 contact_area_radius =
2122 std::sqrt(std::pow(r_e, 2.) - std::pow(r_e + pen_dist, 2.));
2125 contact_area_radius = 0.;
2129 double max_dist = xc1.d_y + p1->d_geom_p->boundingRadius();
2132 double max_y_loc = p1->getXLocal(0).d_y;
2134 for (
size_t i = 0; i < p1->getNumNodes(); i++)
2136 max_y_loc = p1->getXLocal(i).d_y;
2141 setKeyData(
"pen_dist", pen_dist);
2142 setKeyData(
"contact_area_radius", contact_area_radius);
2143 setKeyData(
"max_y", max_y);
2144 setKeyData(
"max_dist", max_dist);
2145 setKeyData(
"max_y_loc", max_y_loc);
2148 return fmt::format(
" Post-processing: max y = {:.6f} \n", max_y);
2153 if (d_outputDeck_p->d_outCriteria ==
"max_particle_dist" &&
2154 d_pDeck_p->d_testName ==
"two_particle") {
2158 const auto &xci = d_particlesListTypeAll[0]->getXCenter();
2159 const auto &xcj = d_particlesListTypeAll[1]->getXCenter();
2163 d_outputDeck_p->d_outCriteriaParams[0])) {
2165 if(d_ppFile.is_open())
2170 else if (d_outputDeck_p->d_outCriteria ==
"max_node_dist") {
2180 auto max_x = d_x[max_pt_and_index.second];
2184 d_outputDeck_p->d_outCriteriaParams[0])) {
2187 if(d_ppFile.is_open())
2190 log(fmt::format(
"{}: Terminating simulation as one of the failing"
2191 " criteria is met. Point ({:.6f}, {:.6f}, {:.6f}) is at "
2194 " allowed distance {:.6f}\n",
2195 d_name, max_x.d_x, max_x.d_y, max_x.d_z, max_x.length(),
2196 d_outputDeck_p->d_outCriteriaParams[0]));
2203 bool continue_dt =
false;
2204 auto check_dt = d_outputDeck_p->d_dtTestOut;
2205 if ((d_n % check_dt == 0) && (d_n >= check_dt))
2212 auto w_id = d_pDeck_p->d_particleIdCompressiveTest;
2213 auto f_dir = d_pDeck_p->d_particleForceDirectionCompressiveTest - 1;
2214 const auto &wall = d_particlesListTypeAll[w_id];
2217 auto dx = wall->getXLocal(0) - wall->getXRefLocal(0);
2218 double wall_penetration = dx[f_dir];
2221 double tot_reaction_force = 0.;
2222 for (
size_t i = 0; i < wall->getNumNodes(); i++) {
2223 tot_reaction_force += wall->getFLocal(i)[f_dir] * wall->getVolLocal(i);
2227 bool use_static_file =
true;
2228 if (use_static_file) {
2229 if (!d_ppFile.is_open()) {
2231 std::string tag_pp_file = d_outputDeck_p->d_tagPPFile.empty() ?
"0" : d_outputDeck_p->d_tagPPFile;
2232 std::string filename = d_outputDeck_p->d_path +
"pp_" +
2233 d_pDeck_p->d_testName +
"_" +
2234 tag_pp_file +
".csv";
2235 d_ppFile.open(filename.c_str(), std::ofstream::out | std::ofstream::app);
2237 d_ppFile <<
"t, delta, force \n";
2240 d_ppFile << fmt::format(
"%4.6e, %4.6e, %4.6e\n", d_time, wall_penetration,
2241 tot_reaction_force);
2244 setKeyData(
"wall_penetration", wall_penetration);
2245 setKeyData(
"tot_reaction_force", tot_reaction_force);
2247 return fmt::format(
" Post-processing: wall penetration = {:"
2249 "reaction force = {:5.3e} \n",
2250 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.