PeriDEM 0.3.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
Data Structures | Functions
anonymous_namespace{main.cpp} Namespace Reference

Data Structures

class  ContactWitness
 
struct  CrackFit
 
struct  CrackMetrics
 
class  GrainContactProbe
 
struct  LawResult
 
struct  PackedBounds
 
struct  PackedParticle
 
struct  Probe
 

Functions

std::string directoryPathWithTrailingSep (const std::filesystem::path &dir)
 
json buildInputJson (const std::string &output_path, const std::filesystem::path &mesh_tri, const std::filesystem::path &mesh_ell, double final_time, size_t num_steps, double mesh_size, double horizon)
 
double maxElem (const std::vector< double > &v)
 
bool doesParticleIntersect (const PackedParticle &p, const std::vector< PackedParticle > &existing, const std::vector< double > &rect, double padding)
 
std::vector< PackedParticle > generateParticleLocations (const std::vector< double > &in_rect, double max_y, double mesh_size, double R, int N_target, double padding, std::mt19937 &gen)
 
PackedBounds packedAxisBounds (const std::vector< PackedParticle > &packed)
 
json contactPairJson (double R_contact_factor, bool damping_on, bool friction_on, double Kn, double beta_n_eps, double friction_coeff, double Kn_factor, double beta_n_factor)
 
double KnFromBulk (double Ka, double Kb, double horizon)
 
json buildInputJson (const std::string &output_path_for_deck, const std::filesystem::path &inp_dir, int argc, char *argv[])
 
std::vector< PackedParticle > generateCircularGrid (int ncols, int nrows, double x_lo, double y_lo, double R, double gap)
 
json buildInputJson (const std::string &output_path_for_deck, const std::filesystem::path &mesh_cir, const std::filesystem::path &mesh_fixed, const std::filesystem::path &mesh_moving, double final_time, size_t num_steps, bool file_mesh, bool write_meshes, size_t search_interval, int ncols, int nrows, const std::string &mpi_strategy="auto")
 
void writeLocations (const std::filesystem::path &csv, const std::vector< PackedParticle > &packed)
 
double maxAbsForceFromCsv (const std::filesystem::path &csv)
 
bool ownsNode (const data::ModelData &data, size_t i)
 
void dumpFinalNodal (const data::ModelData &data, const std::filesystem::path &out_dir)
 
json buildInputJson (const std::string &output_path, const std::filesystem::path &mesh_cir, const std::filesystem::path &mesh_wall, const std::string &mpi_strategy, double final_time, size_t num_steps)
 
bool segmentCrossesVertical (double x0, double y0, double x1, double y1, double x_line, double y_lo, double y_hi)
 
size_t applyNotches (PeriDEMModel &dem, double x_left, double x_right, double y_lo, double y_hi)
 
size_t applyNotchSlots (PeriDEMModel &dem, double notch_half, double notch_w, double y_lo, double y_hi)
 
size_t applyNotchMidplanes (PeriDEMModel &dem, double notch_half, double y_tip, double y_top)
 
double nodePhi (const PeriDEMModel &dem, size_t i)
 
double nodeDamageForFit (const PeriDEMModel &dem, size_t i, bool use_bond_count)
 
CrackFit fitCrackFromPhi (const PeriDEMModel &dem, double tip_x, double tip_y, double outward_sign, double phi_cut, double band_x, double band_y, double exclude_r, double bin_h, bool use_bond_count=false)
 
json meshSetJson (const std::filesystem::path &f, double mesh_size)
 
json uniformMeshSetJson (const std::filesystem::path &f, double mesh_size, const std::vector< std::vector< double > > &voids={})
 
std::vector< std::vector< double > > notchVoidBoxes (double H, double notch_half, double notch_w, double notch_depth, double z_lo, double z_hi)
 
geom::GeomData sillingNotchedPlateGeom (double W, double H, double notch_half, double notch_w, double notch_depth)
 
geom::GeomData bhatVNotchedPlateGeom (double W, double H, double notch_half, double notch_w, double notch_depth)
 
json fixedTopOuterBC (double W, double H, double notch_half)
 
json buildTraskInputJson (const std::string &output_path, const std::filesystem::path &mesh_plate, double W, double H, double notch_half, double notch_w, double notch_depth, double mesh_size, double horizon, double rho, double E, double K, double G, double Gc, double v_impact, double final_time, size_t num_steps)
 
json buildImpactInputJson (const std::string &output_path, const std::filesystem::path &mesh_plate, const std::filesystem::path &mesh_impactor, double W, double H, double notch_half, double notch_w, double notch_depth, double Iw, double Ih, double gap, double mesh_size, double horizon, double Rc_factor, double Kn, double rho, double E, double K, double G, double Gc, double v_impact, double final_time, size_t num_steps, double thickness=0.)
 
json buildBhatInputJson (const std::string &output_path, const std::filesystem::path &mesh_plate, const std::filesystem::path &mesh_impactor, double W, double H, double notch_half, double notch_w, double notch_depth, double Iw, double Ih, double gap, double mesh_size, double horizon, double Rc_factor, double rho, double E, double K, double G, double Gc, double v_impact, double final_time, size_t num_steps)
 
json buildInputJson (const std::string &output_path, const std::filesystem::path &mesh_file, double L, double mesh_size, double horizon, double final_time, size_t num_steps, const std::string &self_contact)
 
size_t applyPrecrack (PeriDEMModel &dem, double x_crack)
 
void applyRigidOpenClose (PeriDEMModel &dem, double x_crack, double u_right)
 
CrackMetrics measure (PeriDEMModel &dem, double x_crack, double band)
 
LawResult runLaw (const std::string &self_contact, const std::filesystem::path &run_dir, double L, double mesh_size, double horizon, unsigned n_threads)
 
json buildInputJson (const std::string &output_path_for_deck, const std::filesystem::path &mesh_file_1, const std::filesystem::path &mesh_file_2, double final_time, size_t num_steps, bool zero_ic, double mesh_size_in, double horizon_in, bool damping_on, double eps_n, bool two_particle_test, bool file_mesh, double beta_n_factor, const std::string &damping_law="com_and_node", const std::string &friction_law="coulomb_simple", bool friction_on=false, double friction_mu=-1., double ic_vx=0., bool bottom_patch_bc=false, const std::string &mpi_strategy="auto")
 
bool is3DGeometry (const std::string &g)
 
bool isAcceptableGeometryName (const std::string &g)
 
json buildInputJson (const std::string &geomName, const std::string &output_path_for_deck, const std::filesystem::path &mesh_file_1, const std::filesystem::path &mesh_file_2)
 
json buildInputJson (const std::string &output_path, const std::filesystem::path &mesh_cir, const std::filesystem::path &mesh_wall, double R, double mesh_size, double horizon, double Rc_factor, double Kn, double final_time, size_t num_steps, bool policy_combo=false)
 
Probe checkWallSpringOnly (PeriDEMModel &dem, geom::GeomObject &wall_geom)
 
double minGrainWallGap (PeriDEMModel &dem, geom::GeomObject &wall_geom)
 
bool near (double a, double b, double tol=1.e-12)
 

Function Documentation

◆ applyNotches()

size_t anonymous_namespace{main.cpp}::applyNotches ( PeriDEMModel &  dem,
double  x_left,
double  x_right,
double  y_lo,
double  y_hi 
)

Definition at line 59 of file main.cpp.

60 {
61 size_t n_broken = 0;
62 for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
63 if (dem.d_ptId[i] != 0)
64 continue;
65 const auto &xi = dem.d_xRef[i];
66 for (size_t k = 0; k < dem.d_neighPd[i].size(); ++k) {
67 const size_t j = dem.d_neighPd[i][k];
68 if (dem.d_ptId[j] != 0)
69 continue;
70 const auto &xj = dem.d_xRef[j];
71 if (segmentCrossesVertical(xi.d_x, xi.d_y, xj.d_x, xj.d_y, x_left, y_lo, y_hi) ||
72 segmentCrossesVertical(xi.d_x, xi.d_y, xj.d_x, xj.d_y, x_right, y_lo, y_hi)) {
73 dem.d_fracture_p->setBondState(i, k, true);
74 ++n_broken;
75 }
76 }
77 }
78 return n_broken;
79}
data::ModelData::d_xRef
std::vector< util::Point > d_xRef
reference positions of the nodes
Definition modelData.h:742
data::ModelData::d_fracture_p
std::unique_ptr< geometry::Fracture > d_fracture_p
Fracture state of bonds.
Definition modelData.h:736
data::ModelData::d_ptId
std::vector< size_t > d_ptId
Global node to particle id (walls are assigned id after last particle id)
Definition modelData.h:764
data::ModelData::d_neighPd
std::vector< std::vector< size_t > > d_neighPd
Neighbor data for peridynamic forces.
Definition modelData.h:770
anonymous_namespace{main.cpp}::segmentCrossesVertical
bool segmentCrossesVertical(double x0, double y0, double x1, double y1, double x_line, double y_lo, double y_hi)
Definition main.cpp:48

References data::ModelData::d_fracture_p, data::ModelData::d_neighPd, data::ModelData::d_ptId, data::ModelData::d_xRef, and segmentCrossesVertical().

Here is the call graph for this function:

◆ applyNotchMidplanes()

size_t anonymous_namespace{main.cpp}::applyNotchMidplanes ( PeriDEMModel &  dem,
double  notch_half,
double  y_tip,
double  y_top 
)

Definition at line 122 of file main.cpp.

123 {
124 auto crossesMid = [&](double x0, double y0, double x1, double y1, double xc) {
125 if ((x0 - xc) * (x1 - xc) >= 0.)
126 return false;
127 const double den = x1 - x0;
128 if (std::abs(den) < 1.0e-30)
129 return false;
130 const double t = (xc - x0) / den;
131 if (t <= 0. || t >= 1.)
132 return false;
133 const double y = y0 + t * (y1 - y0);
134 return y >= y_tip && y <= y_top;
135 };
136 size_t n_broken = 0;
137 for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
138 if (dem.d_ptId[i] != 0)
139 continue;
140 const auto &xi = dem.d_xRef[i];
141 for (size_t k = 0; k < dem.d_neighPd[i].size(); ++k) {
142 const size_t j = dem.d_neighPd[i][k];
143 if (dem.d_ptId[j] != 0)
144 continue;
145 const auto &xj = dem.d_xRef[j];
146 if (crossesMid(xi.d_x, xi.d_y, xj.d_x, xj.d_y, -notch_half) ||
147 crossesMid(xi.d_x, xi.d_y, xj.d_x, xj.d_y, notch_half)) {
148 dem.d_fracture_p->setBondState(i, k, true);
149 ++n_broken;
150 }
151 }
152 }
153 return n_broken;
154}

References data::ModelData::d_fracture_p, data::ModelData::d_neighPd, data::ModelData::d_ptId, and data::ModelData::d_xRef.

Referenced by main().

Here is the caller graph for this function:

◆ applyNotchSlots()

size_t anonymous_namespace{main.cpp}::applyNotchSlots ( PeriDEMModel &  dem,
double  notch_half,
double  notch_w,
double  y_lo,
double  y_hi 
)

Definition at line 82 of file main.cpp.

83 {
84 const double hw = 0.5 * notch_w;
85 auto crossesSlot = [&](double x0, double y0, double x1, double y1, double xc) {
86 const double xa = xc - hw, xb = xc + hw;
87 auto inSlot = [&](double x, double y) {
88 return x >= xa && x <= xb && y >= y_lo && y <= y_hi;
89 };
90 if (inSlot(x0, y0) || inSlot(x1, y1))
91 return true;
92 for (int s = 0; s <= 8; ++s) {
93 const double t = s / 8.0;
94 const double x = x0 + t * (x1 - x0);
95 const double y = y0 + t * (y1 - y0);
96 if (inSlot(x, y))
97 return true;
98 }
99 return segmentCrossesVertical(x0, y0, x1, y1, xc, y_lo, y_hi);
100 };
101 size_t n_broken = 0;
102 for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
103 if (dem.d_ptId[i] != 0)
104 continue;
105 const auto &xi = dem.d_xRef[i];
106 for (size_t k = 0; k < dem.d_neighPd[i].size(); ++k) {
107 const size_t j = dem.d_neighPd[i][k];
108 if (dem.d_ptId[j] != 0)
109 continue;
110 const auto &xj = dem.d_xRef[j];
111 if (crossesSlot(xi.d_x, xi.d_y, xj.d_x, xj.d_y, -notch_half) ||
112 crossesSlot(xi.d_x, xi.d_y, xj.d_x, xj.d_y, notch_half)) {
113 dem.d_fracture_p->setBondState(i, k, true);
114 ++n_broken;
115 }
116 }
117 }
118 return n_broken;
119}

References data::ModelData::d_fracture_p, data::ModelData::d_neighPd, data::ModelData::d_ptId, data::ModelData::d_xRef, and segmentCrossesVertical().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ applyPrecrack()

size_t anonymous_namespace{main.cpp}::applyPrecrack ( PeriDEMModel &  dem,
double  x_crack 
)

Definition at line 114 of file main.cpp.

114 {
115 size_t n_broken = 0;
116 for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
117 const auto &xi = dem.d_xRef[i];
118 for (size_t k = 0; k < dem.d_neighPd[i].size(); ++k) {
119 const size_t j = dem.d_neighPd[i][k];
120 const auto &xj = dem.d_xRef[j];
121 if ((xi.d_x - x_crack) * (xj.d_x - x_crack) < 0.) {
122 dem.d_fracture_p->setBondState(i, k, true);
123 ++n_broken;
124 }
125 }
126 }
127 return n_broken;
128}

References data::ModelData::d_fracture_p, data::ModelData::d_neighPd, and data::ModelData::d_xRef.

Referenced by runLaw().

Here is the caller graph for this function:

◆ applyRigidOpenClose()

void anonymous_namespace{main.cpp}::applyRigidOpenClose ( PeriDEMModel &  dem,
double  x_crack,
double  u_right 
)

Definition at line 130 of file main.cpp.

130 {
131 for (size_t i = 0; i < dem.d_xRef.size(); ++i) {
132 const double ux = (dem.d_xRef[i].d_x >= x_crack) ? u_right : 0.;
133 dem.d_u[i] = util::Point(ux, 0., 0.);
134 dem.d_v[i] = util::Point();
135 dem.d_x[i] = dem.d_xRef[i] + dem.d_u[i];
136 }
137}
data::ModelData::d_x
std::vector< util::Point > d_x
Current positions of the nodes.
Definition modelData.h:745
data::ModelData::d_u
std::vector< util::Point > d_u
Displacement of the nodes.
Definition modelData.h:748
data::ModelData::d_v
std::vector< util::Point > d_v
Velocity of the nodes.
Definition modelData.h:751
util::Point
A structure to represent 3d vectors.
Definition point.h:30

References data::ModelData::d_u, data::ModelData::d_v, data::ModelData::d_x, and data::ModelData::d_xRef.

Referenced by runLaw().

Here is the caller graph for this function:

◆ bhatVNotchedPlateGeom()

geom::GeomData anonymous_namespace{main.cpp}::bhatVNotchedPlateGeom ( double  W,
double  H,
double  notch_half,
double  notch_w,
double  notch_depth 
)

Definition at line 318 of file main.cpp.

319 {
320 const double x0 = -0.5 * W, y0 = -0.5 * H, x1 = 0.5 * W, y1 = 0.5 * H;
321 const double hw = 0.5 * notch_w;
322 const double y_tip = y1 - notch_depth;
323 const double y_top = y1 + 1.0e-6;
324 geom::GeomData plate;
325 plate.d_geomName = "complex";
326 plate.d_geomComplexInfo = {
327 std::vector<std::string>{"rectangle", "triangle", "triangle"},
328 std::vector<std::string>{"plus", "minus", "minus"}};
329 plate.d_geomParams = {
330 // outer plate
331 x0, y0, 0., x1, y1, 0.,
332 // left V: top-left, top-right, tip
333 -notch_half - hw, y_top, 0., -notch_half + hw, y_top, 0., -notch_half, y_tip, 0.,
334 // right V
335 notch_half - hw, y_top, 0., notch_half + hw, y_top, 0., notch_half, y_tip, 0.};
336 return plate;
337}
geom::GeomData
Input data for geometrical objects.
Definition geomObjectsUtil.h:22
geom::GeomData::d_geomParams
std::vector< double > d_geomParams
Zone parameters.
Definition geomObjectsUtil.h:28
geom::GeomData::d_geomName
std::string d_geomName
Zone type.
Definition geomObjectsUtil.h:25
geom::GeomData::d_geomComplexInfo
std::pair< std::vector< std::string >, std::vector< std::string > > d_geomComplexInfo
Zone geometry info if it is a complex type.
Definition geomObjectsUtil.h:34

References geom::GeomData::d_geomComplexInfo, geom::GeomData::d_geomName, and geom::GeomData::d_geomParams.

Referenced by buildBhatInputJson().

Here is the caller graph for this function:

◆ buildBhatInputJson()

json anonymous_namespace{main.cpp}::buildBhatInputJson ( const std::string &  output_path,
const std::filesystem::path &  mesh_plate,
const std::filesystem::path &  mesh_impactor,
double  W,
double  H,
double  notch_half,
double  notch_w,
double  notch_depth,
double  Iw,
double  Ih,
double  gap,
double  mesh_size,
double  horizon,
double  Rc_factor,
double  rho,
double  E,
double  K,
double  G,
double  Gc,
double  v_impact,
double  final_time,
size_t  num_steps 
)

Definition at line 581 of file main.cpp.

587 {
588 auto model = inp::ModelDeck::getExampleJson(2, final_time, num_steps, "finite_difference",
589 "velocity_verlet", true, 2, "Multi_Particle", 0);
590 model["Self_Contact"] = "reference_gap"; // Bhat 2023 §4.4
591 model["Bond_Break"] = "absolute_stretch"; // Bhat §3.1: |s| > |s0|
592 model["Wall_Contact"] = "meshed";
593
594 auto output = inp::OutputDeck::getExampleJson(
595 "vtu", output_path,
596 std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage", "Damage_Bond",
597 "Damage_Z", "Particle_ID"}),
598 std::max<size_t>(1, num_steps / 10), 1, false, "zlib", true, num_steps, "", false);
599
600 json ic = {
601 {"Sets", 1},
602 {"Set_1",
603 {{"Particle_List", std::vector<size_t>{1}},
604 {"Constant_Velocity",
605 {{"Velocity_Vector", std::vector<double>{0., -v_impact, 0.}}}}}}};
606
607 // Paper: only the outer top ligaments are fixed; impactor is free.
608 json outer = fixedTopOuterBC(W, H, notch_half);
609 json bc_disp = {{"Sets", 2}, {"Set_1", outer["Set_1"]}, {"Set_2", outer["Set_2"]}};
610
611 // Bhat specimen: OCC/Gmsh nonconvex plate with V-notch cutouts (§5, Fig. 4).
612 geom::GeomData plate =
613 bhatVNotchedPlateGeom(W, H, notch_half, notch_w, notch_depth);
614 geom::GeomData impactor;
615 impactor.d_geomName = "rectangle";
616 impactor.d_geomParams = {-0.5 * Iw, -0.5 * Ih, 0., 0.5 * Iw, 0.5 * Ih, 0.};
617 auto particle = inp::ParticleDeck::getParticleGeomExampleJson({plate, impactor});
618
619 json mesh = {{"Sets", 2},
620 {"Set_1", meshSetJson(mesh_plate, mesh_size)},
621 {"Set_2", meshSetJson(mesh_impactor, mesh_size)}};
622
623 // Bhat Table 1, constant micromodulus, ν = 1/3 (bond-based, 2D).
624 const double nu_bhat = 1.0 / 3.0;
625 const double c_bhat =
626 6.0 * E / (M_PI * std::pow(horizon, 3.0) * (1.0 - nu_bhat));
627 const double s0_bhat = std::sqrt(4.0 * M_PI * Gc / (9.0 * E * horizon));
628
629 auto material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
630 auto mat_bhat = inp::MaterialDeck::getExampleJson("PMBBond", false, horizon, 0, rho, K, G, Gc,
631 /*computeParamsFromElastic=*/false, 0, E);
632 mat_bhat["Bond_Potential_Params"] = std::vector<double>{c_bhat, s0_bhat};
633 mat_bhat["Influence_Function"] = {{"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
634 material["Set_1"] = mat_bhat;
635 // Free deformable striker (paper). Unbreakable elastic so the contact face
636 // delivers bulk impulse; plate remains PMB with paper (c, s0).
637 auto mat_strike = inp::MaterialDeck::getExampleJson("PDElasticBond", false, horizon, 0, rho, K, G,
638 0., true, 0, E);
639 mat_strike["Influence_Function"] = {{"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
640 material["Set_2"] = mat_strike;
641
642 // Absolute Rc from mesh_size (not Factor×hMin).
643 const double Kn_bhat = 18.0 * K / (M_PI * std::pow(horizon, 4.0));
644 const double Rc_abs = Rc_factor * mesh_size;
645 json contact_base = inp::ContactPairDeck::getExampleJson(
646 Rc_abs, /*computeContactR=*/false, /*damping*/ false, /*friction*/ false, Kn_bhat,
647 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, K);
648 contact_base["Kn"] = Kn_bhat;
649 json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
650 contact["Set_1_1"] = contact_base;
651 contact["Set_1_2"] = contact_base;
652 contact["Set_2_2"] = contact_base;
653 contact["Damping_Law"] = "off"; // Bhat §6: β_d = 0
654 contact["Friction_Law"] = "coulomb_simple"; // μ = 0 below
655 contact["Set_1_1"]["Friction_Coefficient"] = 0.;
656 contact["Set_1_2"]["Friction_Coefficient"] = 0.;
657 contact["Set_2_2"]["Friction_Coefficient"] = 0.;
658
659 const double cy_imp = 0.5 * H + 0.5 * Ih + gap;
660 auto pgen = inp::PGenDeck::getExampleJson("From_File");
661 pgen["Random_Rotation"] = false;
662 pgen["Data"]["N"] = 2;
663 pgen["Data"]["0"] = json{{"x", 0.}, {"y", 0.}, {"z", 0.},
664 {"theta", 0.}, {"s", 1.}, {"geom_id", 0},
665 {"mat_id", 0}, {"contact_id", 0}};
666 pgen["Data"]["1"] = json{{"x", 0.}, {"y", cy_imp}, {"z", 0.},
667 {"theta", 0.}, {"s", 1.}, {"geom_id", 1},
668 {"mat_id", 1}, {"contact_id", 1}};
669
670 return json{{"Model", model},
671 {"Output", output},
672 {"Displacement_BC", bc_disp},
673 {"IC", ic},
674 {"Particle", particle},
675 {"Mesh", mesh},
676 {"Material", material},
677 {"Contact", contact},
678 {"Neighbor", inp::PNeighborDeck::getExampleJson("simple_all", 5.0, 1, 0.5)},
679 {"Particle_Generation", pgen}};
680}
json
nlohmann::ordered_json json
Definition geomObjectsUtil.h:17
anonymous_namespace{main.cpp}::meshSetJson
json meshSetJson(const std::filesystem::path &f, double mesh_size)
Definition main.cpp:260
anonymous_namespace{main.cpp}::bhatVNotchedPlateGeom
geom::GeomData bhatVNotchedPlateGeom(double W, double H, double notch_half, double notch_w, double notch_depth)
Definition main.cpp:318
anonymous_namespace{main.cpp}::fixedTopOuterBC
json fixedTopOuterBC(double W, double H, double notch_half)
Definition main.cpp:339
mesh
Collection of methods and data related to finite element and mesh.
Definition mesh.cpp:28
particle
Collection of methods and data related to particle object.
Definition modelData.h:36
inp::ContactPairDeck::getExampleJson
static json getExampleJson(double contactR=0., bool computeContactR=true, bool dampingOn=true, bool frictionOn=true, double Kn=0., double eps=1., double mu=0., double KnFactor=1., double betanFactor=1., double deltaMax=1., double vMax=0., double K=0.)
Returns example JSON object for ModelDeck configuration.
Definition contactPairDeck.h:98
inp::MaterialDeck::getExampleJson
static json getExampleJson(std::string materialType="PDState", bool isPlainStrain=false, double horizon=-1., double horizonMeshRatio=-1., double density=1., double K=0., double G=0., double Gc=0., bool computeParamsFromElastic=true, size_t influenceFnType=0, double E=-1.)
Returns example JSON object for ModelDeck configuration.
Definition materialDeck.h:368
inp::ModelDeck::getExampleJson
static json getExampleJson(size_t dim=2, double tFinal=1.0, size_t Nt=10, std::string spatialDiscretization="finite_difference", std::string timeDiscretization="central_difference", bool populateElementNodeConnectivity=true, size_t quadOrder=2, std::string particleSimType="Multi_Particle", int seed=0)
Returns example JSON object for ModelDeck configuration.
Definition modelDeck.h:164
inp::OutputDeck::getExampleJson
static json getExampleJson(std::string outFormat="vtu", std::string path="./", std::vector< std::string > outTags={"Displacement"}, size_t outputInterval=1, size_t debug=2, bool performFEOut=true, std::string compressType="zlib", bool performOut=true, size_t dtTestOut=1, std::string tagPPFile="", bool pvdCollection=false)
Returns example JSON object for ModelDeck configuration.
Definition outputDeck.h:144
inp::PGenDeck::getExampleJson
static json getExampleJson(std::string genMethod="From_File")
Returns example JSON object for ModelDeck configuration.
Definition pGenDeck.h:65
inp::PNeighborDeck::getExampleJson
static json getExampleJson(std::string updateCriteria="simple_all", double sFactor=1., size_t neighUpdateInterval=1, double nearBdNodesTol=0.5)
Returns example JSON object for ModelDeck configuration.
Definition pNeighborDeck.h:62
inp::ParticleDeck::getParticleContactExampleJson
static json getParticleContactExampleJson(size_t nSets=0)
Returns example JSON object for ModelDeck configuration.
Definition particleDeck.h:259
inp::ParticleDeck::getParticleGeomExampleJson
static json getParticleGeomExampleJson(std::vector< geom::GeomData > pGeomVec=std::vector< geom::GeomData >())
Returns example JSON object for ModelDeck configuration.
Definition particleDeck.h:95
inp::ParticleDeck::getParticleMaterialExampleJson
static json getParticleMaterialExampleJson(size_t nSets=0)
Returns example JSON object for ModelDeck configuration.
Definition particleDeck.h:209

References bhatVNotchedPlateGeom(), geom::GeomData::d_geomName, geom::GeomData::d_geomParams, fixedTopOuterBC(), inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), and meshSetJson().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildImpactInputJson()

json anonymous_namespace{main.cpp}::buildImpactInputJson ( const std::string &  output_path,
const std::filesystem::path &  mesh_plate,
const std::filesystem::path &  mesh_impactor,
double  W,
double  H,
double  notch_half,
double  notch_w,
double  notch_depth,
double  Iw,
double  Ih,
double  gap,
double  mesh_size,
double  horizon,
double  Rc_factor,
double  Kn,
double  rho,
double  E,
double  K,
double  G,
double  Gc,
double  v_impact,
double  final_time,
size_t  num_steps,
double  thickness = 0. 
)

Definition at line 431 of file main.cpp.

438 {
439 // Silling 2003 EMU KW: explicit central difference; brittle-microelastic PMB;
440 // broken bond force = 0; plate boundaries load-free (no Displacement_BC);
441 // rigid impactor of finite mass.
442 // thickness > 0 runs his actual 3D case (200 x 100 x 9 grid, 9 mm plate).
443 const bool dim3 = thickness > 0.;
444 const size_t dim = dim3 ? 3 : 2;
445 // Element-node connectivity is only used for strain/stress output and does not
446 // support hexahedra, which is what the 3D structured grid produces.
447 auto model = inp::ModelDeck::getExampleJson(dim, final_time, num_steps, "finite_difference",
448 "central_difference", !dim3, 2, "Multi_Particle",
449 0);
450 model["Self_Contact"] = "none";
451 model["Bond_Break"] = "tension";
452 model["Wall_Contact"] = "meshed";
453 // Silling Fig. 2: rigid cylinder of 1.57 kg. In 2D the analogue is the mass
454 // per unit thickness of his 9 mm plate, M/t = 1.57/0.009 kg/m.
455 model["Rigid_Particles"] = json::array(
456 {json{{"Id", 1}, {"Mass", dim3 ? 1.57 : 1.57 / 0.009}}});
457
458 auto output = inp::OutputDeck::getExampleJson(
459 "vtu", output_path,
460 std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage", "Damage_Bond",
461 "Damage_Z", "Particle_ID"}),
462 std::max<size_t>(1, num_steps / 10), 1, false, "zlib", true, num_steps, "", false);
463
464 json ic = {
465 {"Sets", 1},
466 {"Set_1",
467 {{"Particle_List", std::vector<size_t>{1}},
468 {"Constant_Velocity",
469 {{"Velocity_Vector", std::vector<double>{0., -v_impact, 0.}}}}}}};
470
471 // The impactor is rigid but of finite mass (Model.Rigid_Particles above), so
472 // it decelerates on contact instead of being driven through the plate.
473 json bc_disp = {{"Sets", 1},
474 {"Set_1",
475 {{"Particle_List", std::vector<size_t>{1}},
476 {"Direction", std::vector<size_t>{1}},
477 {"Time_Function",
478 {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
479 {"Spatial_Function", {{"Type", "constant"}}},
480 {"Zero_Displacement", true}}}};
481
482 // Plate is the full box on a structured grid with the two 1.5 mm notch slots
483 // removed from the mesh, so they are real gaps (Fig. 2) rather than material
484 // whose bonds have merely been cut.
485 geom::GeomData plate;
486 geom::GeomData impactor;
487 if (dim3) {
488 // Keep both bodies centred on z = 0 so particle generation places them
489 // consistently (it positions objects by their centre).
490 plate.d_geomName = "cuboid";
491 plate.d_geomParams = {-0.5 * W, -0.5 * H, -0.5 * thickness,
492 0.5 * W, 0.5 * H, 0.5 * thickness};
493 // Silling Fig. 4: the cylinder axis is along the impact direction and its
494 // flat end face strikes the plate edge. Params: radius, centre of the
495 // beginning cross-section, and the vector to the end cross-section.
496 impactor.d_geomName = "cylinder";
497 impactor.d_geomParams = {0.5 * Iw, 0., -0.5 * Ih, 0., 0., Ih, 0.};
498 } else {
499 plate.d_geomName = "rectangle";
500 plate.d_geomParams = {-0.5 * W, -0.5 * H, 0., 0.5 * W, 0.5 * H, 0.};
501 // In-plane section of that same cylinder: Iw wide, Ih long, flat face down.
502 impactor.d_geomName = "rectangle";
503 impactor.d_geomParams = {-0.5 * Iw, -0.5 * Ih, 0., 0.5 * Iw, 0.5 * Ih, 0.};
504 }
505 auto particle = inp::ParticleDeck::getParticleGeomExampleJson({plate, impactor});
506
507 const auto voids =
508 notchVoidBoxes(H, notch_half, notch_w, notch_depth,
509 dim3 ? -0.5 * thickness - 1.0e-9 : -1.0e-9,
510 dim3 ? 0.5 * thickness + 1.0e-9 : 1.0e-9);
511 // Plate on Silling's equally spaced structured grid. In 2D the impactor is a
512 // rectangle so it goes on the same grid; in 3D it is a cylinder, which the
513 // uniform grid cannot represent, so it is meshed by Gmsh at the same spacing.
514 json mesh = {{"Sets", 2},
515 {"Set_1", uniformMeshSetJson(mesh_plate, mesh_size, voids)},
516 {"Set_2", dim3 ? meshSetJson(mesh_impactor, mesh_size)
517 : uniformMeshSetJson(mesh_impactor, mesh_size)}};
518
519 // Impactor is steel like the plate but unbreakable, and its nodal forces are
520 // replaced by the rigid-body acceleration each step, so its bond stiffness
521 // never enters the solution.
522 auto material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
523 material["Set_1"] = inp::MaterialDeck::getExampleJson(
524 "PMBBond", false, horizon, 0, rho, K, G, Gc, true, 0, E);
525 material["Set_2"] = inp::MaterialDeck::getExampleJson(
526 "PDElasticBond", false, horizon, 0, rho, K, G, 0., true, 0, E);
527 material["Set_1"]["Influence_Function"] = {
528 {"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
529 material["Set_2"]["Influence_Function"] = {
530 {"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
531
532 // Absolute Contact_Radius from the requested mesh size (not Factor×hMin).
533 const double Rc_abs = Rc_factor * mesh_size;
534 json contact_base = inp::ContactPairDeck::getExampleJson(
535 Rc_abs, /*computeContactR=*/false, /*damping*/ false, /*friction*/ false, Kn,
536 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, K);
537 contact_base["Kn"] = Kn;
538 json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
539 contact["Set_1_1"] = contact_base;
540 contact["Set_1_2"] = contact_base;
541 contact["Set_2_2"] = contact_base;
542 contact["Damping_Law"] = "off";
543 contact["Friction_Law"] = "coulomb_simple";
544
545 const double cy_imp = 0.5 * H + 0.5 * Ih + gap; // flat face `gap` above the edge
546 auto pgen = inp::PGenDeck::getExampleJson("From_File");
547 pgen["Random_Rotation"] = false;
548 pgen["Data"]["N"] = 2;
549 pgen["Data"]["0"] = json{{"x", 0.},
550 {"y", 0.},
551 {"z", 0.},
552 {"theta", 0.},
553 {"s", 1.},
554 {"geom_id", 0},
555 {"mat_id", 0},
556 {"contact_id", 0}};
557 pgen["Data"]["1"] = json{{"x", 0.},
558 {"y", cy_imp},
559 {"z", 0.},
560 {"theta", 0.},
561 {"s", 1.},
562 {"geom_id", 1},
563 {"mat_id", 1},
564 {"contact_id", 1}};
565
566 // Plate has no Displacement_BC (Silling: load-free all around); the only
567 // Displacement_BC drives the rigid impactor.
568 return json{{"Model", model},
569 {"Output", output},
570 {"Displacement_BC", bc_disp},
571 {"IC", ic},
572 {"Particle", particle},
573 {"Mesh", mesh},
574 {"Material", material},
575 {"Contact", contact},
576 {"Neighbor", inp::PNeighborDeck::getExampleJson("simple_all", 5.0, 1, 0.5)},
577 {"Particle_Generation", pgen}};
578}
anonymous_namespace{main.cpp}::notchVoidBoxes
std::vector< std::vector< double > > notchVoidBoxes(double H, double notch_half, double notch_w, double notch_depth, double z_lo, double z_hi)
Definition main.cpp:284
anonymous_namespace{main.cpp}::uniformMeshSetJson
json uniformMeshSetJson(const std::filesystem::path &f, double mesh_size, const std::vector< std::vector< double > > &voids={})
Definition main.cpp:272

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), meshSetJson(), notchVoidBoxes(), and uniformMeshSetJson().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [1/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  geomName,
const std::string &  output_path_for_deck,
const std::filesystem::path &  mesh_file_1,
const std::filesystem::path &  mesh_file_2 
)

Definition at line 65 of file main.cpp.

67 {
68
69 const double s = 0.001;
70 const double mesh_size = s / 5.0;
71 const double horizon = 3.0 * mesh_size;
72 const double particle_dist = 0.001;
73
74 const double poisson1 = 0.25;
75 const double rho1 = 1200.0;
76 const double K1 = 2.16e+7;
77 const double E1 = material::toE(K1, poisson1);
78 const double G1 = material::toGE(E1, poisson1);
79 const double Gc1 = 50.0;
80
81 const double poisson2 = 0.25;
82 const double rho2 = 1200.0;
83 const double K2 = 2.16e+7;
84 const double E2 = material::toE(K2, poisson2);
85 const double G2 = material::toGE(E2, poisson2);
86 const double Gc2 = 50.0;
87
88 const double R_contact_factor = 0.95;
89 const double Kn_11 = 18.0 * util::harmonicMean(K1, K1) / (M_PI * std::pow(horizon, 5));
90 const double Kn_22 = 18.0 * util::harmonicMean(K2, K2) / (M_PI * std::pow(horizon, 5));
91 const double Kn_12 = 18.0 * util::harmonicMean(K1, K2) / (M_PI * std::pow(horizon, 5));
92 const double beta_n_eps = 0.9;
93 const double friction_coeff = 0.5;
94 const double beta_n_factor = 100.0;
95
96 const int model_dim = is3DGeometry(geomName) ? 3 : 2;
97
98 /* Reference shapes at origin for both mesh groups (same as twop_circ_inbuilt). */
99 const util::Point origin(0.0, 0.0, 0.0);
100 const std::vector<double> gp_ref = geom::exampleGeomParams(geomName, origin, s);
101
102 /* Same integration window and IC gap as twop_circ_inbuilt (stable dt). */
103 const double final_time = 0.0001; // 0.012;
104 const size_t num_steps = 1000; //36000;
105 const size_t dt_out_n = num_steps / 4; // 10;
106
107 auto modelDeckJson = inp::ModelDeck::getExampleJson(static_cast<size_t>(model_dim), final_time, num_steps,
108 "finite_difference", "central_difference",
109 true, 2, "Multi_Particle", 0);
110
111 auto outputDeckJson = inp::OutputDeck::getExampleJson("vtu", output_path_for_deck,
112 std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage_Z", "Damage", "Particle_ID"}),
113 dt_out_n, 2, true, "zlib", true, 1, "", true);
114
115 auto bcDeckJson = inp::BCDeck::getExampleJson(0, 1, 1, true, util::Point(0, -10, 0));
116
117 bcDeckJson["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson("Displacement_BC", false, geom::GeomData(),
118 {0}, {}, "", {}, "", {},
119 {1, 2}, true, "", {});
120
121 const double free_fall_dist = particle_dist - horizon;
122 const double free_fall_vel = -std::sqrt(2.0 * std::abs(-10.0) * std::max(1e-30, free_fall_dist));
123 bcDeckJson["IC"]["Set_1"] = inp::BCBaseDeck::getExampleJson("IC", false, geom::GeomData(),
124 {1}, {}, "", {}, "", {},
125 {}, false, "Constant_Velocity", {0.0, free_fall_vel, 0.0});
126
127 auto pDeckJson = json({});
128
129 std::vector<geom::GeomData> pGeomVec(2);
130 pGeomVec[0].d_geomName = geomName;
131 pGeomVec[0].d_geomParams = gp_ref;
132 pGeomVec[1].d_geomName = geomName;
133 pGeomVec[1].d_geomParams = gp_ref;
134
135 pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
136
137 const std::string f1 = mesh_file_1.string();
138 const std::string f2 = mesh_file_2.string();
139 auto meshSet1 = json({{"File", f1},
140 {"CreateMesh",
141 {{"Flag", true},
142 {"Info", "gmsh_builtin_mesh"},
143 {"Mesh_Size", mesh_size},
144 {"Write_Mesh_File", true}}}});
145 auto meshSet2 = json({{"File", f2},
146 {"CreateMesh",
147 {{"Flag", true},
148 {"Info", "gmsh_builtin_mesh"},
149 {"Mesh_Size", mesh_size},
150 {"Write_Mesh_File", true}}}});
151 pDeckJson["Mesh"] = json({{"Sets", 2}, {"Set_1", meshSet1}, {"Set_2", meshSet2}});
152
153 auto pMatJson = inp::ParticleDeck::getParticleMaterialExampleJson(2);
154 pMatJson["Set_1"] = inp::MaterialDeck::getExampleJson("PDState", false, horizon,
155 0, rho1, K1, G1, Gc1, true, 1);
156 pMatJson["Set_2"] = inp::MaterialDeck::getExampleJson("PDState", false, horizon,
157 0, rho2, K2, G2, Gc2, true, 1);
158 pDeckJson["Material"] = pMatJson;
159
160 auto pContactJson = inp::ParticleDeck::getParticleContactExampleJson(2);
161 json contact_base = inp::ContactPairDeck::getExampleJson(R_contact_factor,
162 true, false, false,
163 Kn_11, beta_n_eps, friction_coeff, 1.0, beta_n_factor, 1.0, 0.0, 0.0);
164
165 pContactJson["Set_1_1"] = contact_base;
166 pContactJson["Set_1_1"]["Kn"] = Kn_11;
167 pContactJson["Set_1_2"] = contact_base;
168 pContactJson["Set_1_2"]["Kn"] = Kn_12;
169 pContactJson["Set_2_2"] = contact_base;
170 pContactJson["Set_2_2"]["Kn"] = Kn_22;
171 pDeckJson["Contact"] = pContactJson;
172
173 pDeckJson["Neighbor"] = inp::PNeighborDeck::getExampleJson("simple_all", 10.0, 40, 0.5);
174
175 auto pGenJson = inp::PGenDeck::getExampleJson("From_File");
176 pGenJson["Data"]["N"] = 2;
177
178 /*
179 * Placement matches twop_circ: lower particle corner at (s,s[,s]); upper at
180 * (s, 2*s + particle_dist + s, [s]) so vertical gap is particle_dist between
181 * outer extents (same as R1=R2=s in twop_circ).
182 */
183 if (model_dim == 2) {
184 pGenJson["Data"]["0"] = {
185 {"x", s}, {"y", s}, {"z", 0.0},
186 {"theta", 0.0}, {"s", 1.0},
187 {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
188 {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0}
189 };
190 pGenJson["Data"]["1"] = {
191 {"x", s}, {"y", 2.0 * s + particle_dist + s}, {"z", 0.0},
192 {"theta", M_PI}, {"s", 1.0},
193 {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
194 {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}
195 };
196 } else {
197 pGenJson["Data"]["0"] = {
198 {"x", s}, {"y", s}, {"z", s},
199 {"theta", 0.0}, {"s", 1.0},
200 {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
201 {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0}
202 };
203 pGenJson["Data"]["1"] = {
204 {"x", s}, {"y", 2.0 * s + particle_dist + s}, {"z", s},
205 {"theta", M_PI}, {"s", 1.0},
206 {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
207 {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}
208 };
209 }
210 pDeckJson["Particle_Generation"] = pGenJson;
211
212 return json({{"Model", modelDeckJson},
213 {"Output", outputDeckJson},
214 {"Force_BC", bcDeckJson["Force_BC"]},
215 {"Displacement_BC", bcDeckJson["Displacement_BC"]},
216 {"IC", bcDeckJson["IC"]},
217 {"Particle", pDeckJson["Particle"]},
218 {"Mesh", pDeckJson["Mesh"]},
219 {"Material", pDeckJson["Material"]},
220 {"Contact", pDeckJson["Contact"]},
221 {"Neighbor", pDeckJson["Neighbor"]},
222 {"Particle_Generation", pDeckJson["Particle_Generation"]}});
223}
anonymous_namespace{main.cpp}::is3DGeometry
bool is3DGeometry(const std::string &g)
Definition main.cpp:52
anonymous_namespace{testParallelCompLib.cpp}::f1
double f1(const double &x)
Definition testParallelCompLib.cpp:44
anonymous_namespace{testParallelCompLib.cpp}::f2
double f2(const double &x)
Definition testParallelCompLib.cpp:48
geom::exampleGeomParams
std::vector< double > exampleGeomParams(const std::string &geom_type, const util::Point &c, double s)
Canonical parameter vector for examples, mesh generation tests, and demos.
Definition geomObjectsUtil.cpp:1241
material::toGE
double toGE(double E, double nu)
Compute shear modulus from Young's modulus E and Poisson's ratio nu.
Definition materialUtil.h:87
material::toE
double toE(double K, double nu)
Compute Young's modulus E from Bulk modulus K and Poisson's ratio nu.
Definition materialUtil.h:51
util::harmonicMean
double harmonicMean(const double &m1, const double &m2)
Definition function.cpp:131
inp::BCBaseDeck::getExampleJson
static json getExampleJson(std::string type="Foce_BC", bool isRegionActive=false, geom::GeomData regionGeomData=geom::GeomData(), std::vector< size_t > pList=std::vector< size_t >(), std::vector< size_t > pNotList=std::vector< size_t >(), std::string timeFnType="", std::vector< double > timeFnParams=std::vector< double >(), std::string spatialFnType="", std::vector< double > spatialFnParams=std::vector< double >(), std::vector< size_t > direction=std::vector< size_t >(), bool isDisplacementZero=false, std::string icType="", std::vector< double > icVec=std::vector< double >())
Returns example JSON object for ModelDeck configuration.
Definition bcBaseDeck.h:157
inp::BCDeck::getExampleJson
static json getExampleJson(size_t nForceSets=0, size_t nDispSets=0, size_t nICSets=0, bool gravityActive=false, util::Point gravity=util::Point())
Returns example JSON object for ModelDeck configuration.
Definition bcDeck.h:69

References geom::exampleGeomParams(), inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::BCDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), util::harmonicMean(), is3DGeometry(), material::toE(), and material::toGE().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [2/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  output_path,
const std::filesystem::path &  mesh_cir,
const std::filesystem::path &  mesh_wall,
const std::string &  mpi_strategy,
double  final_time,
size_t  num_steps 
)

Definition at line 202 of file main.cpp.

206 {
207 const double R = 0.001;
208 const double mesh_size = R / 4.0;
209 const double horizon = 3.0 * mesh_size;
210 const double h_est = 0.7 * mesh_size;
211 const double Rc_est = 0.95 * h_est;
212
213 // Start inside PP Rc and near the floor; soft Kn (below) keeps contact stable.
214 const double gap_pp = 0.5 * Rc_est;
215 const double cy = R + 0.35 * Rc_est;
216 const double cx0 = -(R + 0.5 * gap_pp);
217 const double cx1 = +(R + 0.5 * gap_pp);
218
219 const double wall_half = 3.5 * R;
220 const double wall_thick = std::max(2.0 * mesh_size, horizon);
221
222 geom::GeomData grain;
223 grain.d_geomName = "circle";
224 grain.d_geomParams = {R, 0., 0., 0.};
225 geom::GeomData floor;
226 floor.d_geomName = "rectangle";
227 floor.d_geomParams = {-wall_half, -wall_thick, 0., wall_half, 0., 0.};
228
229 const double poisson = 0.25;
230 const double rho = 1200.0;
231 const double K = 2.16e7;
232 const double E = material::toE(K, poisson);
233 const double G = material::toGE(E, poisson);
234 const double Gc = 50.0;
235 const double Kn =
236 18.0 * util::harmonicMean(K, K) / (M_PI * std::pow(horizon, 5));
237
238 auto model = inp::ModelDeck::getExampleJson(
239 2, final_time, num_steps, "finite_difference", "central_difference", true,
240 2, "Multi_Particle", 0);
241 model["MPI_Strategy"] = mpi_strategy;
242 model["Wall_Contact"] = "meshed";
243
244 auto output = inp::OutputDeck::getExampleJson(
245 "vtu", output_path,
246 std::vector<std::string>({"Displacement", "Velocity", "Force",
247 "Particle_ID"}),
248 std::max<size_t>(1, num_steps), 1, false, "zlib", false, num_steps, "",
249 false);
250
251 // Wall (particle index 2) fully fixed; downward IC on both grains.
252 auto bc = inp::BCDeck::getExampleJson(0, 1, 1, true, util::Point(0, -10, 0));
253 bc["Displacement_BC"]["Set_1"] = json{
254 {"Particle_List", std::vector<size_t>{2}},
255 {"Direction", std::vector<size_t>{1, 2}},
256 {"Time_Function",
257 {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
258 {"Spatial_Function", {{"Type", "constant"}}},
259 {"Zero_Displacement", true}};
260 bc["IC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
261 "IC", false, geom::GeomData(), {0, 1}, {}, "", {}, "", {}, {}, false,
262 "Constant_Velocity", {0., -0.02, 0.});
263
264 auto particle = inp::ParticleDeck::getParticleGeomExampleJson({grain, floor});
265 auto mesh_set = [&](const std::filesystem::path &f) {
266 return json{{"File", f.string()},
267 {"CreateMesh",
268 {{"Flag", true},
269 {"Info", "gmsh_builtin_mesh"},
270 {"Mesh_Size", mesh_size},
271 {"Write_Mesh_File", true}}}};
272 };
273 json mesh = {{"Sets", 2},
274 {"Set_1", mesh_set(mesh_cir)},
275 {"Set_2", mesh_set(mesh_wall)}};
276
277 json material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
278 material["Set_1"] = inp::MaterialDeck::getExampleJson(
279 "PDState", false, horizon, 0, rho, K, G, Gc, true, 1);
280 material["Set_2"] = inp::MaterialDeck::getExampleJson(
281 "PDState", false, horizon, 0, rho, K, G, Gc, true, 1);
282
283 json contact_base = inp::ContactPairDeck::getExampleJson(
284 0.95, true, /*damping*/ false, /*friction*/ false, Kn, 0.95, 0.0, 1.0,
285 100.0, 1.0, 0.0, K);
286 // Soften contact spring so short MPI identity runs stay stable and bit-reproducible.
287 contact_base["Kn_Factor"] = 1.0e-4;
288 json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
289 contact["Set_1_1"] = contact_base;
290 contact["Set_1_2"] = contact_base;
291 contact["Set_2_2"] = contact_base;
292 contact["Damping_Law"] = "off";
293 contact["Friction_Law"] = "coulomb_simple";
294
295 auto pgen = inp::PGenDeck::getExampleJson("From_File");
296 pgen["Random_Rotation"] = false;
297 pgen["Data"]["N"] = 3;
298 pgen["Data"]["0"] = json{{"x", cx0},
299 {"y", cy},
300 {"z", 0.},
301 {"theta", 0.},
302 {"s", 1.},
303 {"geom_id", 0},
304 {"mat_id", 0},
305 {"contact_id", 0}};
306 pgen["Data"]["1"] = json{{"x", cx1},
307 {"y", cy},
308 {"z", 0.},
309 {"theta", 0.},
310 {"s", 1.},
311 {"geom_id", 0},
312 {"mat_id", 0},
313 {"contact_id", 0}};
314 pgen["Data"]["2"] = json{{"x", 0.},
315 {"y", -0.5 * wall_thick},
316 {"z", 0.},
317 {"theta", 0.},
318 {"s", 1.},
319 {"geom_id", 1},
320 {"mat_id", 1},
321 {"contact_id", 1},
322 {"is_wall", true}};
323
324 util::io::print(std::format(
325 "mpi_identity_twop_wall: R={:.4g} mesh={:.4g} Rc~{:.4g} gap_pp={:.4g} "
326 "cy={:.4g} T={:.3g} N={}\n",
327 R, mesh_size, Rc_est, gap_pp, cy, final_time, num_steps));
328
329 return json{{"Model", model},
330 {"Output", output},
331 {"Force_BC", bc["Force_BC"]},
332 {"Displacement_BC", bc["Displacement_BC"]},
333 {"IC", bc["IC"]},
334 {"Particle", particle},
335 {"Mesh", mesh},
336 {"Material", material},
337 {"Contact", contact},
338 {"Neighbor",
339 inp::PNeighborDeck::getExampleJson("simple_all", 5.0, 1, 0.5)},
340 {"Particle_Generation", pgen}};
341}
util::io::print
void print(const T &msg, int nt=print_default_tab, int printMpiRank=print_default_mpi_rank)
Prints formatted information.
Definition io.h:128

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::BCDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), util::harmonicMean(), util::io::print(), material::toE(), and material::toGE().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [3/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  output_path,
const std::filesystem::path &  mesh_cir,
const std::filesystem::path &  mesh_wall,
double  R,
double  mesh_size,
double  horizon,
double  Rc_factor,
double  Kn,
double  final_time,
size_t  num_steps,
bool  policy_combo = false 
)

Definition at line 44 of file main.cpp.

49 {
50 auto model = inp::ModelDeck::getExampleJson(2, final_time, num_steps,
51 "finite_difference",
52 "central_difference", true, 2,
53 "Multi_Particle", 0);
54 model["Wall_Contact"] = "analytical_plane";
55 if (policy_combo) {
56 model["Self_Contact"] = "reference_gap";
57 model["Bond_Break"] = "absolute_stretch";
58 }
59
60 auto output = inp::OutputDeck::getExampleJson(
61 "vtu", output_path,
62 std::vector<std::string>({"Displacement", "Force", "Particle_ID"}),
63 std::max<size_t>(1, num_steps / 5), 1, false, "zlib", false, num_steps, "",
64 false);
65
66 auto bc = inp::BCDeck::getExampleJson(0, 1, 0, false, util::Point());
67 bc["Displacement_BC"]["Set_1"] = json{
68 {"Particle_List", std::vector<size_t>{1}},
69 {"Direction", std::vector<size_t>{1, 2}},
70 {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
71 {"Spatial_Function", {{"Type", "constant"}}},
72 {"Zero_Displacement", true}};
73
74 const double wall_half = 3.0 * R;
75 const double wall_thick = mesh_size;
76 geom::GeomData grain;
77 grain.d_geomName = "circle";
78 grain.d_geomParams = {R, 0., 0., 0.};
79 geom::GeomData floor;
80 floor.d_geomName = "rectangle";
81 floor.d_geomParams = {-wall_half, -wall_thick, 0., wall_half, 0., 0.};
82
83 auto particle = inp::ParticleDeck::getParticleGeomExampleJson({grain, floor});
84
85 auto mesh_set = [&](const std::filesystem::path &f) {
86 return json{{"File", f.string()},
87 {"CreateMesh",
88 {{"Flag", true},
89 {"Info", "gmsh_builtin_mesh"},
90 {"Mesh_Size", mesh_size},
91 {"Write_Mesh_File", true}}}};
92 };
93 json mesh = {{"Sets", 2},
94 {"Set_1", mesh_set(mesh_cir)},
95 {"Set_2", mesh_set(mesh_wall)}};
96
97 const double rho = 1200.0;
98 const double K = 2.16e7;
99 const double G = 1.296e7;
100 const double Gc = 50.0;
101 auto material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
102 material["Set_1"] =
103 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho, K, G, Gc,
104 true, 1);
105 material["Set_2"] =
106 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho, K, G, Gc,
107 true, 1);
108
109 json contact_base = inp::ContactPairDeck::getExampleJson(
110 Rc_factor, true, /*damping*/ false, /*friction*/ false, Kn, 1.0, 0.0, 1.0,
111 100.0, 1.0, 0.0, K);
112 // Soften spring for a short stable settle (formula check uses this Kn too).
113 contact_base["Kn"] = Kn;
114 contact_base["Kn_Factor"] = 1.0;
115 json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
116 contact["Set_1_1"] = contact_base;
117 contact["Set_1_2"] = contact_base;
118 contact["Set_2_2"] = contact_base;
119 contact["Damping_Law"] = policy_combo ? "node" : "off";
120 contact["Friction_Law"] = policy_combo ? "stick_slip" : "coulomb_simple";
121
122 // Bottom of grain near the floor (y=0): slight overlap into Rc after setup.
123 const double cy = R + 0.25 * mesh_size;
124 auto pgen = inp::PGenDeck::getExampleJson("From_File");
125 pgen["Random_Rotation"] = false;
126 pgen["Data"]["N"] = 2;
127 pgen["Data"]["0"] = json{{"x", 0.},
128 {"y", cy},
129 {"z", 0.},
130 {"theta", 0.},
131 {"s", 1.},
132 {"geom_id", 0},
133 {"mat_id", 0},
134 {"contact_id", 0}};
135 pgen["Data"]["1"] = json{{"x", 0.},
136 {"y", -0.5 * wall_thick},
137 {"z", 0.},
138 {"theta", 0.},
139 {"s", 1.},
140 {"geom_id", 1},
141 {"mat_id", 1},
142 {"contact_id", 1},
143 {"is_wall", true}};
144
145 return json{{"Model", model},
146 {"Output", output},
147 {"Displacement_BC", bc["Displacement_BC"]},
148 {"Particle", particle},
149 {"Mesh", mesh},
150 {"Material", material},
151 {"Contact", contact},
152 {"Neighbor",
153 inp::PNeighborDeck::getExampleJson("simple_all", 5.0, 1, 0.5)},
154 {"Particle_Generation", pgen}};
155}

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::BCDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), and inp::ParticleDeck::getParticleMaterialExampleJson().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [4/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  output_path,
const std::filesystem::path &  mesh_file,
double  L,
double  mesh_size,
double  horizon,
double  final_time,
size_t  num_steps,
const std::string &  self_contact 
)

Definition at line 44 of file main.cpp.

46 {
47 auto model = inp::ModelDeck::getExampleJson(2, final_time, num_steps, "finite_difference",
48 "central_difference", true, 2, "Single_Particle",
49 0);
50 model["Self_Contact"] = self_contact;
51 model["Bond_Break"] = "tension";
52
53 // Soft PD; large Gc so only the precrack is fractured.
54 const double rho = 1200.0;
55 const double K = 2.16e5;
56 const double G = 1.296e5;
57 const double Gc = 5.0e6;
58
59 auto output = inp::OutputDeck::getExampleJson(
60 "vtu", output_path,
61 std::vector<std::string>({"Displacement", "Force", "Damage_Z", "Particle_ID"}),
62 num_steps, 1, false, "zlib", false, num_steps, "", false);
63
64 const double strip = 0.2 * L;
65 json bc = json::object();
66 bc["Displacement_BC"] = {
67 {"Sets", 2},
68 {"Set_1",
69 {{"Region",
70 {{"Geometry",
71 {{"Type", "rectangle"},
72 {"Parameters", std::vector<double>{-0.05 * L, -0.05 * L, 0., strip, L + 0.05 * L, 0.}}}}}},
73 {"Direction", std::vector<size_t>{1, 2}},
74 {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
75 {"Spatial_Function", {{"Type", "constant"}}},
76 {"Zero_Displacement", true}}},
77 {"Set_2",
78 {{"Region",
79 {{"Geometry",
80 {{"Type", "rectangle"},
81 {"Parameters",
82 std::vector<double>{L - strip, -0.05 * L, 0., L + 0.05 * L, L + 0.05 * L, 0.}}}}}},
83 {"Direction", std::vector<size_t>{1}},
84 {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
85 {"Spatial_Function", {{"Type", "constant"}}},
86 {"Zero_Displacement", false}}}};
87
88 geom::GeomData rect;
89 rect.d_geomName = "rectangle";
90 rect.d_geomParams = {0., 0., 0., L, L, 0.};
91 auto particle = inp::ParticleDeck::getParticleGeomExampleJson({rect});
92
93 json mesh_set = {{"File", mesh_file.string()},
94 {"CreateMesh",
95 {{"Flag", true},
96 {"Info", "gmsh_builtin_mesh"},
97 {"Mesh_Size", mesh_size},
98 {"Write_Mesh_File", true}}}};
99 json mesh = {{"Set_1", mesh_set}};
100
101 auto material = inp::ParticleDeck::getParticleMaterialExampleJson(1);
102 material["Set_1"] =
103 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho, K, G, Gc, true, 1);
104
105 json root = {{"Model", model},
106 {"Output", output},
107 {"Particle", particle},
108 {"Mesh", mesh},
109 {"Material", material}};
110 root.merge_patch(bc);
111 return root;
112}

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::ModelDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), and inp::ParticleDeck::getParticleMaterialExampleJson().

Referenced by runLaw().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [5/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  output_path,
const std::filesystem::path &  mesh_tri,
const std::filesystem::path &  mesh_ell,
double  final_time,
size_t  num_steps,
double  mesh_size,
double  horizon 
)

Definition at line 43 of file main.cpp.

45 {
46 // Short tip — point only, not a tall triangle. Narrower tip concentrates load.
47 const double W = 0.0008;
48 const double H = 0.0012;
49
50 // Hollow ellipse (visibly elliptical ring).
51 const double a_out = 0.0020;
52 const double b_out = 0.0014;
53 // Thin ring so a tip-seeded crack can bipartition.
54 const double a_in = 0.00180;
55 const double b_in = 0.00122;
56 const double ell_theta = 0.0;
57 const double tip_gap = 0.00008;
58 const double tip_y = H;
59 const double ell_cy = tip_y + tip_gap + b_out;
60
61 const double rho_t = 1200.0;
62 const double K_t = 2.16e7;
63 const double nu_t = 0.25;
64 const double E_t = material::toE(K_t, nu_t);
65 const double G_t = material::toGE(E_t, nu_t);
66 const double Gc_t = 200.0;
67
68 const double rho_e = 1200.0;
69 const double K_e = 2.16e7;
70 const double nu_e = 0.25;
71 const double E_e = material::toE(K_e, nu_e);
72 const double G_e = material::toGE(E_e, nu_e);
73 // Low Gc so tip damage runs into a diametral through-crack.
74 const double Gc_e = 1.0;
75
76 const double R_contact_factor = 0.95;
77 const double Kn_tt = 18.0 * util::harmonicMean(K_t, K_t) / (M_PI * std::pow(horizon, 5));
78 const double Kn_ee = 18.0 * util::harmonicMean(K_e, K_e) / (M_PI * std::pow(horizon, 5));
79 const double Kn_te = 18.0 * util::harmonicMean(K_t, K_e) / (M_PI * std::pow(horizon, 5));
80
81 const double g = 10.0;
82 const double ic_vy = -2.5;
83
84 const size_t dt_out_n = std::max<size_t>(1, num_steps / 40);
85
86 auto model = inp::ModelDeck::getExampleJson(2, final_time, num_steps, "finite_difference",
87 "central_difference", true, 2, "Multi_Particle", 0);
88 model["Bond_Break"] = "tension";
89 // No self-contact across broken bonds — otherwise the crack stays glued.
90 model["Self_Contact"] = "none";
91
92 auto output = inp::OutputDeck::getExampleJson(
93 "vtu", output_path,
94 std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage_Z", "Damage",
95 "Particle_ID", "Fixity"}),
96 dt_out_n, 1, true, "zlib", true, 1, "", true);
97
98 auto bc = inp::BCDeck::getExampleJson(0, 1, 1, true, util::Point(0., -g, 0.));
99 bc["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
100 "Displacement_BC", false, geom::GeomData(), {0}, {}, "", {}, "", {}, {1, 2}, true, "", {});
101
102 bc["IC"]["Set_1"] = {
103 {"Particle_List", std::vector<size_t>{1}},
104 {"Constant_Velocity", {{"Velocity_Vector", std::vector<double>{0., ic_vy, 0.}}}}};
105
106 std::vector<geom::GeomData> pGeom(2);
107 pGeom[0].d_geomName = "triangle";
108 pGeom[0].d_geomParams = {-0.5 * W, 0., 0., 0.5 * W, 0., 0., 0., H, 0.};
109
110 pGeom[1].d_geomName = "ellipse_minus_ellipse";
111 pGeom[1].d_geomParams = {a_out, b_out, a_in, b_in, ell_theta, 0., 0., 0.};
112
113 for (auto &g : pGeom)
114 geom::createGeomObject(g);
115
116 const util::Point c_tri = pGeom[0].d_geom_p->center();
117
118 json pDeck = json::object();
119 pDeck["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeom);
120
121 pDeck["Mesh"] = {
122 {"Sets", 2},
123 {"Set_1",
124 {{"File", mesh_tri.string()},
125 {"CreateMesh",
126 {{"Flag", true},
127 {"Info", "gmsh_builtin_mesh"},
128 {"Mesh_Size", mesh_size},
129 {"Write_Mesh_File", true}}}}},
130 {"Set_2",
131 {{"File", mesh_ell.string()},
132 {"CreateMesh",
133 {{"Flag", true},
134 {"Info", "gmsh_builtin_mesh"},
135 {"Mesh_Size", mesh_size},
136 {"Write_Mesh_File", true}}}}}};
137
138 auto mat = inp::ParticleDeck::getParticleMaterialExampleJson(2);
139 mat["Set_1"] =
140 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho_t, K_t, G_t, Gc_t, true, 1);
141 mat["Set_2"] =
142 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho_e, K_e, G_e, Gc_e, true, 1);
143 pDeck["Material"] = mat;
144
145 auto contact = inp::ParticleDeck::getParticleContactExampleJson(2);
146 json cpair_wall = inp::ContactPairDeck::getExampleJson(R_contact_factor, true, true, false, Kn_tt, 0.95,
147 0.0, 1.0, 100.0, 1.0, 0.0, K_t);
148 // Tip contact: enough Beta_n to drive a crack, damping on to limit spray.
149 json cpair_tip = inp::ContactPairDeck::getExampleJson(0.90, true, true, false, Kn_te, 0.4,
150 0.0, 1.0, 8.0, 1.0, 0.0,
151 util::harmonicMean(K_t, K_e));
152 json cpair_ell = inp::ContactPairDeck::getExampleJson(R_contact_factor, true, true, false, Kn_ee, 0.95,
153 0.0, 1.0, 100.0, 1.0, 0.0, K_e);
154 contact["Set_1_1"] = cpair_wall;
155 contact["Set_1_1"]["Kn"] = Kn_tt;
156 contact["Set_1_1"]["K"] = K_t;
157 contact["Set_1_2"] = cpair_tip;
158 contact["Set_1_2"]["Kn"] = Kn_te;
159 contact["Set_1_2"]["K"] = util::harmonicMean(K_t, K_e);
160 contact["Set_2_2"] = cpair_ell;
161 contact["Set_2_2"]["Kn"] = Kn_ee;
162 contact["Set_2_2"]["K"] = K_e;
163 contact["Damping_Law"] = "com_and_node";
164 contact["Friction_Law"] = "coulomb_simple";
165 pDeck["Contact"] = contact;
166
167 pDeck["Neighbor"] = inp::PNeighborDeck::getExampleJson("simple_all", 8.0, 5, 0.5);
168
169 auto gen = inp::PGenDeck::getExampleJson("From_File");
170 gen["Random_Rotation"] = false;
171 gen["Data"]["N"] = 2;
172 gen["Data"]["0"] = {{"x", c_tri.d_x}, {"y", c_tri.d_y}, {"z", 0.0},
173 {"theta", 0.0}, {"s", 1.0},
174 {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0},
175 {"is_wall", true}};
176 gen["Data"]["1"] = {{"x", 0.0}, {"y", ell_cy}, {"z", 0.0},
177 {"theta", 0.0}, {"s", 1.0},
178 {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}};
179 pDeck["Particle_Generation"] = gen;
180
181 json j = json::object();
182 j["Model"] = model;
183 j["Output"] = output;
184 j["Force_BC"] = bc["Force_BC"];
185 j["Displacement_BC"] = bc["Displacement_BC"];
186 j["IC"] = bc["IC"];
187 j["Particle"] = pDeck["Particle"];
188 j["Mesh"] = pDeck["Mesh"];
189 j["Material"] = pDeck["Material"];
190 j["Contact"] = pDeck["Contact"];
191 j["Neighbor"] = pDeck["Neighbor"];
192 j["Particle_Generation"] = pDeck["Particle_Generation"];
193 return j;
194}
geom::createGeomObject
void createGeomObject(const std::string &geom_type, const std::vector< double > &params, const std::vector< std::string > &vec_type, const std::vector< std::string > &vec_flag, std::shared_ptr< geom::GeomObject > &obj, bool perform_check)
Definition geomObjectsUtil.cpp:665
util::Point::d_y
double d_y
the y coordinate
Definition point.h:36
util::Point::d_x
double d_x
the x coordinate
Definition point.h:33

References geom::createGeomObject(), util::Point::d_x, util::Point::d_y, inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::BCDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), util::harmonicMean(), material::toE(), and material::toGE().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [6/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  output_path_for_deck,
const std::filesystem::path &  inp_dir,
int  argc,
char *  argv[] 
)

Definition at line 203 of file main.cpp.

204 {
205
206 util::io::InputParser input(argc, argv);
207
208 const std::vector<double> center = {0.0, 0.0, 0.0};
209 const double R = 0.001;
210 const double mesh_size = R / 5.0;
211 const double horizon = 2.0 * mesh_size;
212
213 const double Lin = 0.05;
214 const double Win = 0.04;
215
216 constexpr double wall_top_inset = 1e-7;
217 constexpr double annulus_inset = 1e-7;
218 /* Clearance between void floor / plate and particle circum-bounds (~1–2 mesh). */
219 constexpr double clearance_mesh = 1.5;
220 /* Meshed triangles/hex/drum can extend slightly beyond circumcircle r in VTU. */
221 const double geom_pad = 0.25 * mesh_size;
222 const double plate_thickness = std::max(3.0 * mesh_size, 2.0 * mesh_size);
223
224 const double w_drum2d = R * 0.2;
225
226 const double final_time = 0.001;
227 size_t num_steps = 400;
228 if (input.cmdOptionExists("-numSteps"))
229 num_steps = static_cast<size_t>(std::stoul(input.getCmdOption("-numSteps")));
230
231 const size_t num_outputs = 10;
232 const size_t dt_out_n = std::max<size_t>(1, num_steps / num_outputs);
233 const size_t test_dt_out_n = std::max<size_t>(1, dt_out_n / 10);
234
235 const double rho_wall = 600.;
236 const double poisson_wall = 0.25;
237 const double K_wall = 1.e+4;
238 const double E_wall = material::toE(K_wall, poisson_wall);
239 const double G_wall = material::toGE(E_wall, poisson_wall);
240 const double KIc_wall = 5e+6;
241 const double Gc_wall = material::toGc(KIc_wall, poisson_wall, E_wall);
242
243 const double rho_p = 600.;
244 const double poisson_p = poisson_wall;
245 const double K_p = 5.e+3;
246 const double E_p = material::toE(K_p, poisson_p);
247 const double G_p = material::toGE(E_p, poisson_p);
248 const double KIc_p = 5e+6;
249 const double Gc_p = material::toGc(KIc_p, poisson_p, E_p);
250
251 const double R_contact_factor = 0.95;
252 const double padding = 1.1 * R_contact_factor * mesh_size;
253 const int N_target = 500;
254
255 const std::vector<double> in_rect = {center[0] - 0.5 * Lin, center[1] - 0.5 * Win, center[2],
256 center[0] + 0.5 * Lin, center[1] + 0.5 * Win, center[2]};
257
258 std::mt19937 gen(30);
259 const double max_y = in_rect[4] - clearance_mesh * mesh_size;
260 std::vector<PackedParticle> packed =
261 generateParticleLocations(in_rect, max_y, mesh_size, R, N_target, padding, gen);
262 if (packed.empty())
263 throw std::runtime_error("compression_large_set_inbuilt_2mat2contact: particle pack is empty");
264 for (auto &p : packed)
265 p.z = 0.0;
266
267 const double m = clearance_mesh * mesh_size;
268 /* Seat bed: lowest circum-bottom = void_floor + m (void floor = seed inner bottom). */
269 {
270 const PackedBounds bb0 = packedAxisBounds(packed);
271 const double delta_y = (in_rect[1] + m) - bb0.min_y;
272 for (auto &p : packed)
273 p.y += delta_y;
274 }
275
276 const PackedBounds bb = packedAxisBounds(packed);
277
278 /* Inner void (Gmsh “remove” rectangle): circum-bounds ± m ± geom_pad in x; floor at seed; top = domain. */
279 const double void_lox = bb.min_x - m - geom_pad;
280 const double void_hix = bb.max_x + m + geom_pad;
281 const double void_loy = in_rect[1];
282 /* Moving plate: bottom = top of bed + clearance + pad so FE nodes stay under the plate. */
283 const double plate_bottom_y = bb.max_y + m + geom_pad;
284 double plate_top_y = plate_bottom_y + plate_thickness;
285
286 /* Outer solid: wrap void with wall thickness; extend upward just past the plate (no huge empty band). */
287 const double side = std::max(1.5 * horizon, 2.0 * mesh_size);
288 const double out_lox = void_lox - side;
289 const double out_hix = void_hix + side;
290 const double out_loy = void_loy - side;
291 const double out_hiy = plate_top_y + wall_top_inset + 2.5 * mesh_size;
292
293 std::vector<double> out_rect = {out_lox, out_loy, center[2], out_hix, out_hiy, center[2]};
294
295 /* Open-top void: inner cut runs to outer top (inset for Gmsh strict-inside). */
296 const double void_hi_y = out_rect[4] - annulus_inset;
297 std::vector<double> remove_rect = {void_lox, void_loy, center[2], void_hix, void_hi_y, center[2]};
298
299 util::io::print(std::format(
300 "[compression_large_set_inbuilt_2mat2contact] bbox(center±r) x∈[{:.6f},{:.6f}] y∈[{:.6f},{:.6f}]; "
301 "void x∈[{:.6f},{:.6f}] y_lo {:.6f}; plate y [{:.6f},{:.6f}]; void_hi {:.6f}; outer y∈[{:.6f},{:.6f}]\n",
302 bb.min_x, bb.max_x, bb.min_y, bb.max_y, void_lox, void_hix, void_loy, plate_bottom_y, plate_top_y,
303 void_hi_y, out_rect[1], out_rect[4]));
304
305 std::vector<double> moving_rect = {void_lox, plate_bottom_y, center[2], void_hix, plate_top_y, center[2]};
306
307 std::vector<double> fixed_container_params;
308 fixed_container_params.insert(fixed_container_params.end(), remove_rect.begin(), remove_rect.end());
309 fixed_container_params.insert(fixed_container_params.end(), out_rect.begin(), out_rect.end());
310
311 const size_t n_pack = packed.size();
312 const size_t n_wall_fixed = n_pack;
313 const size_t n_wall_moving = n_pack + 1;
314 const size_t n_total = n_pack + 2;
315
316 /* Reference geometries at origin; packing uses circum-radius ~ R for all zones. */
317 const double ell_a = 0.95 * R;
318 const double ell_b = 0.58 * R;
319 const double ell_theta = 0.35;
320 const double sq_half = R / std::sqrt(2.0);
321 /* Axis-aligned rectangle (~same circum extent as R). */
322 const double rx = 0.82 * R;
323 const double ry = 0.62 * R;
324
325 std::vector<geom::GeomData> pGeomVec(10);
326 pGeomVec[0].d_geomName = "circle";
327 pGeomVec[0].d_geomParams = {R, center[0], center[1], center[2]};
328 pGeomVec[1].d_geomName = "triangle";
329 pGeomVec[1].d_geomParams = {R, center[0], center[1], center[2]};
330 pGeomVec[2].d_geomName = "drum2d";
331 pGeomVec[2].d_geomParams = {R, w_drum2d, center[0], center[1], center[2]};
332 pGeomVec[3].d_geomName = "hexagon";
333 pGeomVec[3].d_geomParams = {R, center[0], center[1], center[2]};
334 pGeomVec[4].d_geomName = "ellipse";
335 pGeomVec[4].d_geomParams = {ell_a, ell_b, ell_theta, center[0], center[1], center[2]};
336 pGeomVec[5].d_geomName = "rectangle";
337 pGeomVec[5].d_geomParams = {-rx, -ry, center[2], rx, ry, center[2]};
338 pGeomVec[6].d_geomName = "square";
339 pGeomVec[6].d_geomParams = {-sq_half, -sq_half, center[2], sq_half, sq_half, center[2]};
340 pGeomVec[7].d_geomName = "circle_minus_circle";
341 pGeomVec[7].d_geomParams = {center[0], center[1], center[2], R, 0.35 * R};
342 pGeomVec[8].d_geomName = "rectangle_minus_rectangle";
343 pGeomVec[8].d_geomParams = fixed_container_params;
344 pGeomVec[9].d_geomName = "rectangle";
345 pGeomVec[9].d_geomParams = moving_rect;
346
347 for (auto &g : pGeomVec)
348 geom::createGeomObject(g);
349
350 /* Annulus composite centroid can pick up tiny numerical z; 2D setup keeps all sites on z = 0. */
351 const util::Point cfix = pGeomVec[8].d_geom_p->center();
352 const util::Point cmov = pGeomVec[9].d_geom_p->center();
353 const util::Point site_wall_fixed(cfix.d_x, cfix.d_y, 0.0);
354 const util::Point site_wall_moving(cmov.d_x, cmov.d_y, 0.0);
355
356 auto modelDeckJson = inp::ModelDeck::getExampleJson(2, final_time, num_steps, "finite_difference",
357 "central_difference", true, 2, "Multi_Particle", 0);
358
359 std::vector<std::string> out_tags = {"Displacement", "Velocity", "Force", "Force_Density", "Damage_Z",
360 "Damage", "Nodal_Volume", "Zone_ID", "Particle_ID", "Fixity",
361 "Force_Fixity", "Contact_Nodes", "No_Fail_Node", "Boundary_Node_Flag"};
362 /* Perform_FE_Out must be true so VTU includes element connectivity (appendMesh). If false, only
363 * points are written (appendNodes) and ParaView’s default “Surface” view looks empty — use
364 * Representation → Points, or enable this flag. Large packs often leave FE out. */
365 auto outputDeckJson = inp::OutputDeck::getExampleJson("vtu", output_path_for_deck, out_tags, dt_out_n, 2,
366 true, "zlib", true, test_dt_out_n, "0", true);
367
368 auto bcDeckJson = inp::BCDeck::getExampleJson(0, 2, 0, true, util::Point(0, -10, 0));
369
370 bcDeckJson["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson("Displacement_BC", false, geom::GeomData(),
371 {n_wall_fixed}, {}, "", {}, "", {},
372 {1, 2}, true, "", {});
373
374 json set2;
375 set2["Particle_List"] = std::vector<size_t>{n_wall_moving};
376 set2["Direction"] = std::vector<size_t>{2};
377 set2["Time_Function"] = json{{"Type", "linear"}, {"Parameters", std::vector<double>{-0.06}}};
378 set2["Spatial_Function"] = json{{"Type", "constant"}};
379 bcDeckJson["Displacement_BC"]["Set_2"] = set2;
380
381 json pDeckJson = json::object();
382
383 pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
384
385 json meshRoot = json{{"Sets", 10}};
386 const char *mesh_names[] = {"mesh_cir", "mesh_tri", "mesh_drum2d", "mesh_hex",
387 "mesh_ellipse", "mesh_rect", "mesh_square", "mesh_circirc",
388 "mesh_fixed_container", "mesh_moving_container"};
389 for (int zi = 0; zi < 10; ++zi) {
390 const std::string fname = (inp_dir / (std::string(mesh_names[zi]) + ".msh")).string();
391 meshRoot["Set_" + std::to_string(zi + 1)] =
392 json{{"File", fname},
393 {"CreateMesh",
394 json{{"Flag", true}, {"Info", "gmsh_builtin_mesh"}, {"Mesh_Size", mesh_size}, {"Write_Mesh_File", true}}}};
395 }
396 pDeckJson["Mesh"] = meshRoot;
397
398 /* Two material laws: 0 = particle, 1 = wall (fixed + moving). geom_id still selects mesh shape (10 meshes). */
399 json matRoot = json{{"Sets", 2}};
400 matRoot["Set_1"] = inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho_p, K_p, G_p, Gc_p, true, 1);
401 matRoot["Set_2"] =
402 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho_wall, K_wall, G_wall, Gc_wall, true, 1);
403 pDeckJson["Material"] = matRoot;
404
405 const double Kn_pp = KnFromBulk(K_p, K_p, horizon);
406 const double Kn_pw = KnFromBulk(K_p, K_wall, horizon);
407 const double Kn_ww = 0.0;
408
409 const double beta_n_eps = 0.95;
410 const double friction_coeff = 0.5;
411 const bool damping_on = true;
412 const bool friction_on = false;
413 const double beta_n_factor = 100.;
414 const double Kn_factor = 1.;
415
416 /* Contact zones: 0 = granular, 1 = wall (fixed + moving). 2×2 matrix → three unique pairs (1-1, 1-2, 2-2). */
417 const json j_contact_pp =
418 contactPairJson(R_contact_factor, damping_on, friction_on, Kn_pp, beta_n_eps, friction_coeff, Kn_factor,
419 beta_n_factor);
420 const json j_contact_pw =
421 contactPairJson(R_contact_factor, damping_on, friction_on, Kn_pw, beta_n_eps, friction_coeff, Kn_factor,
422 beta_n_factor);
423 const json j_contact_ww =
424 contactPairJson(R_contact_factor, damping_on, friction_on, Kn_ww, beta_n_eps, friction_coeff, Kn_factor,
425 beta_n_factor);
426
427 json contactRoot = inp::ContactDeck::getExampleJson(2);
428 contactRoot["Set_1_1"] = j_contact_pp;
429 contactRoot["Set_1_2"] = j_contact_pw;
430 contactRoot["Set_2_2"] = j_contact_ww;
431 pDeckJson["Contact"] = contactRoot;
432
433 pDeckJson["Neighbor"] = inp::PNeighborDeck::getExampleJson("simple_all", 10.0, 100, 0.5);
434
435 auto pGenJson = inp::PGenDeck::getExampleJson("From_File");
436 pGenJson["Random_Rotation"] = false;
437 pGenJson["Data"]["N"] = n_total;
438
439 constexpr size_t k_mat_particle = 0;
440 constexpr size_t k_mat_wall = 1;
441 constexpr size_t k_contact_grains = 0;
442 constexpr size_t k_contact_wall = 1;
443
444 for (size_t pi = 0; pi < n_pack; ++pi) {
445 const auto &p = packed[pi];
446 pGenJson["Data"][std::to_string(pi)] = json{
447 {"x", p.x}, {"y", p.y}, {"z", 0.0}, {"theta", p.theta},
448 {"s", 1.0}, {"geom_id", static_cast<size_t>(p.zone)},
449 {"mat_id", k_mat_particle},
450 {"contact_id", k_contact_grains},
451 };
452 }
453
454 pGenJson["Data"][std::to_string(n_pack)] = json{{"x", site_wall_fixed.d_x},
455 {"y", site_wall_fixed.d_y},
456 {"z", 0.0},
457 {"theta", 0.0},
458 {"s", 1.0},
459 {"geom_id", size_t(8)},
460 {"mat_id", k_mat_wall},
461 {"contact_id", k_contact_wall}};
462
463 pGenJson["Data"][std::to_string(n_pack + 1)] = json{{"x", site_wall_moving.d_x},
464 {"y", site_wall_moving.d_y},
465 {"z", 0.0},
466 {"theta", 0.0},
467 {"s", 1.0},
468 {"geom_id", size_t(9)},
469 {"mat_id", k_mat_wall},
470 {"contact_id", k_contact_wall}};
471
472 pDeckJson["Particle_Generation"] = pGenJson;
473
474 /* Top-level comment only: ignored by inp::Input (only known sections are parsed). */
475 return json{{"Comment",
476 "compression_large_set_inbuilt"},
477 {"Model", modelDeckJson},
478 {"Output", outputDeckJson},
479 {"Force_BC", bcDeckJson["Force_BC"]},
480 {"Displacement_BC", bcDeckJson["Displacement_BC"]},
481 {"Particle", pDeckJson["Particle"]},
482 {"Mesh", pDeckJson["Mesh"]},
483 {"Material", pDeckJson["Material"]},
484 {"Contact", pDeckJson["Contact"]},
485 {"Neighbor", pDeckJson["Neighbor"]},
486 {"Particle_Generation", pDeckJson["Particle_Generation"]}};
487}
util::io::InputParser
Input command line argument parser.
Definition inputParser.h:28
anonymous_namespace{main.cpp}::generateParticleLocations
std::vector< PackedParticle > generateParticleLocations(const std::vector< double > &in_rect, double max_y, double mesh_size, double R, int N_target, double padding, std::mt19937 &gen)
Definition main.cpp:75
anonymous_namespace{main.cpp}::packedAxisBounds
PackedBounds packedAxisBounds(const std::vector< PackedParticle > &packed)
Definition main.cpp:166
anonymous_namespace{main.cpp}::KnFromBulk
double KnFromBulk(double Ka, double Kb, double horizon)
Definition main.cpp:199
anonymous_namespace{main.cpp}::contactPairJson
json contactPairJson(double R_contact_factor, bool damping_on, bool friction_on, double Kn, double beta_n_eps, double friction_coeff, double Kn_factor, double beta_n_factor)
Definition main.cpp:181
material::toGc
double toGc(double KIc, double nu, double E)
Compute critical energy release rate Gc from critical stress-intensity factor KIc,...
Definition materialUtil.h:123
inp::ContactDeck::getExampleJson
static json getExampleJson(size_t nSets=0)
Returns example JSON object for ModelDeck configuration.
Definition contactDeck.h:77

References util::io::InputParser::cmdOptionExists(), contactPairJson(), geom::createGeomObject(), util::Point::d_x, util::Point::d_y, generateParticleLocations(), util::io::InputParser::getCmdOption(), inp::ModelDeck::getExampleJson(), inp::BCDeck::getExampleJson(), inp::ContactDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), KnFromBulk(), anonymous_namespace{main.cpp}::PackedBounds::max_x, anonymous_namespace{main.cpp}::PackedBounds::max_y, anonymous_namespace{main.cpp}::PackedBounds::min_x, anonymous_namespace{main.cpp}::PackedBounds::min_y, packedAxisBounds(), util::io::print(), material::toE(), material::toGc(), and material::toGE().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [7/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  output_path_for_deck,
const std::filesystem::path &  mesh_cir,
const std::filesystem::path &  mesh_fixed,
const std::filesystem::path &  mesh_moving,
double  final_time,
size_t  num_steps,
bool  file_mesh,
bool  write_meshes,
size_t  search_interval,
int  ncols,
int  nrows,
const std::string &  mpi_strategy = "auto" 
)

Definition at line 83 of file main.cpp.

89 {
90
91 if (ncols < 1 || nrows < 1)
92 throw std::runtime_error("jha2021_comp: ncols and nrows must be >= 1");
93 const int ngrains = ncols * nrows;
94
95 const std::vector<double> center = {0.0, 0.0, 0.0};
96 const double R = 0.001;
97 const double mesh_size = R / 5.0;
98 const double horizon = 3.0 * mesh_size;
99 // Realized hmin on this disk is ~0.7 lc. Start just outside Rc so contact
100 // is not present at t=0; plate + gravity pull a few pairs into Rc.
101 const double h_est = 0.7 * mesh_size;
102 const double Rc_est = 0.95 * h_est;
103 const double particle_padding = 1.15 * Rc_est;
104 const double wpd = particle_padding;
105 const double rwp = horizon + wpd;
106 constexpr double wall_vy = -0.06;
107
108 const double Lin = 2.0 * particle_padding + 2.0 * R +
109 static_cast<double>(ncols - 1) * (2.0 * R + particle_padding);
110 const double Win = 2.0 * particle_padding + 2.0 * R +
111 static_cast<double>(nrows - 1) * (2.0 * R + particle_padding);
112
113 const double wall_t = rwp - wpd;
114 const std::vector<double> mw_rect = {center[0] - wpd, Win, center[2],
115 Lin + wpd, Win + wall_t, center[2]};
116 const std::vector<double> cup_channel = {center[0] - rwp, center[1] - rwp,
117 Lin + rwp, Win + wall_t, wall_t,
118 center[2]};
119
120 std::vector<geom::GeomData> pGeomVec(3);
121 pGeomVec[0].d_geomName = "circle";
122 pGeomVec[0].d_geomParams = {R, center[0], center[1], center[2]};
123 pGeomVec[1].d_geomName = "open_rect_channel_2d";
124 pGeomVec[1].d_geomParams = cup_channel;
125 pGeomVec[2].d_geomName = "rectangle";
126 pGeomVec[2].d_geomParams = mw_rect;
127
128 std::vector<PackedParticle> packed =
129 generateCircularGrid(ncols, nrows, 0.0, 0.0, R, particle_padding);
130 if (static_cast<int>(packed.size()) != ngrains)
131 throw std::runtime_error("jha2021_comp: grid count mismatch");
132
133 const double poisson = 0.25;
134 const double rho = 1200.0;
135 const double K = 2.16e+7;
136 const double E = material::toE(K, poisson);
137 const double G = material::toGE(E, poisson);
138 const double Gc = 50.0;
139 const double Kn = 18.0 * util::harmonicMean(K, K) / (M_PI * std::pow(horizon, 5));
140
141 const double dt = final_time / static_cast<double>(num_steps);
142 const double c_wave = std::sqrt(E / rho);
143 const double dt_cfl = h_est / c_wave;
144 util::io::print(std::format(
145 "jha2021_comp: N={} ({}x{}), gap={:.6e} (1.15 Rc_est), wpd={:.6e}, "
146 "dt={:.6e} s, h/c={:.6e} s, dt/(h/c)={:.3f}\n",
147 ngrains, ncols, nrows, particle_padding, wpd, dt, dt_cfl, dt / dt_cfl));
148
149 const size_t num_outputs = 4;
150 const size_t dt_out_n = std::max<size_t>(1, num_steps / num_outputs);
151 const size_t test_dt_out_n = std::max<size_t>(1, dt_out_n / 10);
152
153 for (auto &g : pGeomVec)
154 geom::createGeomObject(g);
155
156 const util::Point cfix = pGeomVec[1].d_geom_p->center();
157 const util::Point cmov = pGeomVec[2].d_geom_p->center();
158 const size_t n_pack = packed.size();
159 const size_t n_wall_fixed = n_pack;
160 const size_t n_wall_moving = n_pack + 1;
161 const size_t n_total = n_pack + 2;
162
163 auto modelDeckJson = inp::ModelDeck::getExampleJson(
164 2, final_time, num_steps, "finite_difference", "central_difference", true, 2,
165 "Multi_Particle", 0);
166 modelDeckJson["MPI_Strategy"] = mpi_strategy;
167
168 std::vector<std::string> out_tags = {"Displacement", "Velocity", "Force", "Damage_Z",
169 "Damage", "Particle_ID", "Contact_Nodes"};
170 auto outputDeckJson = inp::OutputDeck::getExampleJson(
171 "vtu", output_path_for_deck, out_tags, dt_out_n, 1, true, "zlib", true,
172 test_dt_out_n, "0", true);
173
174 auto bcDeckJson =
175 inp::BCDeck::getExampleJson(0, 2, 0, true, util::Point(0, -10, 0));
176 bcDeckJson["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
177 "Displacement_BC", false, geom::GeomData(), {n_wall_fixed}, {}, "", {}, "", {},
178 {1, 2}, true, "", {});
179 json set2;
180 set2["Particle_List"] = std::vector<size_t>{n_wall_moving};
181 set2["Direction"] = std::vector<size_t>{2};
182 set2["Time_Function"] = json{{"Type", "linear"}, {"Parameters", std::vector<double>{wall_vy}}};
183 set2["Spatial_Function"] = json{{"Type", "constant"}};
184 bcDeckJson["Displacement_BC"]["Set_2"] = set2;
185
186 json pDeckJson = json::object();
187 pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
188
189 auto meshEntry = [&](const std::filesystem::path &f) -> json {
190 if (file_mesh)
191 return json{{"File", f.string()}};
192 return json{{"File", f.string()},
193 {"CreateMesh",
194 json{{"Flag", true},
195 {"Info", "gmsh_builtin_mesh"},
196 {"Mesh_Size", mesh_size},
197 {"Write_Mesh_File", write_meshes}}}};
198 };
199 pDeckJson["Mesh"] = json{{"Sets", 3},
200 {"Set_1", meshEntry(mesh_cir)},
201 {"Set_2", meshEntry(mesh_fixed)},
202 {"Set_3", meshEntry(mesh_moving)}};
203
204 json matRoot = json{{"Sets", 2}};
205 matRoot["Set_1"] =
206 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho, K, G, Gc, true, 1);
207 matRoot["Set_2"] =
208 inp::MaterialDeck::getExampleJson("PDState", false, horizon, 0, rho, K, G, Gc, true, 1);
209 pDeckJson["Material"] = matRoot;
210
211 json contact_base = inp::ContactPairDeck::getExampleJson(
212 0.95, true, true, false, Kn, 0.95, 0.0, 1.0, 100.0, 1.0, 0.0, K);
213 json contactRoot = inp::ParticleDeck::getParticleContactExampleJson(2);
214 contactRoot["Set_1_1"] = contact_base;
215 contactRoot["Set_1_2"] = contact_base;
216 contactRoot["Set_2_2"] = contact_base;
217 pDeckJson["Contact"] = contactRoot;
218 pDeckJson["Neighbor"] =
219 inp::PNeighborDeck::getExampleJson("simple_all", 5.0, search_interval, 0.5);
220
221 auto pGenJson = inp::PGenDeck::getExampleJson("From_File");
222 pGenJson["Random_Rotation"] = false;
223 pGenJson["Data"]["N"] = n_total;
224 for (size_t pi = 0; pi < n_pack; ++pi) {
225 const auto &p = packed[pi];
226 pGenJson["Data"][std::to_string(pi)] = json{{"x", p.x},
227 {"y", p.y},
228 {"z", 0.0},
229 {"theta", 0.0},
230 {"s", 1.0},
231 {"geom_id", size_t(0)},
232 {"mat_id", size_t(0)},
233 {"contact_id", size_t(0)}};
234 }
235 pGenJson["Data"][std::to_string(n_pack)] = json{{"x", cfix.d_x},
236 {"y", cfix.d_y},
237 {"z", 0.0},
238 {"theta", 0.0},
239 {"s", 1.0},
240 {"geom_id", size_t(1)},
241 {"mat_id", size_t(1)},
242 {"contact_id", size_t(1)},
243 {"is_wall", true}};
244 pGenJson["Data"][std::to_string(n_pack + 1)] = json{{"x", cmov.d_x},
245 {"y", cmov.d_y},
246 {"z", 0.0},
247 {"theta", 0.0},
248 {"s", 1.0},
249 {"geom_id", size_t(2)},
250 {"mat_id", size_t(1)},
251 {"contact_id", size_t(1)},
252 {"is_wall", true}};
253 pDeckJson["Particle_Generation"] = pGenJson;
254
255 util::io::print(std::format(
256 "jha2021_comp: Lin={:.6f}, Win={:.6f}, horizon={:.6f}, plate v_y={}, Wall_Id={}\n",
257 Lin, Win, horizon, wall_vy, n_wall_moving));
258
259 return json{{"Comment", "jha2021_comp_contact"},
260 {"Model", modelDeckJson},
261 {"Output", outputDeckJson},
262 {"Force_BC", bcDeckJson["Force_BC"]},
263 {"Displacement_BC", bcDeckJson["Displacement_BC"]},
264 {"Particle", pDeckJson["Particle"]},
265 {"Mesh", pDeckJson["Mesh"]},
266 {"Material", pDeckJson["Material"]},
267 {"Contact", pDeckJson["Contact"]},
268 {"Neighbor", pDeckJson["Neighbor"]},
269 {"Particle_Generation", pDeckJson["Particle_Generation"]},
270 {"Test", json{{"Test_Name", "compressive_test"},
271 {"Compressive_Test",
272 json{{"Wall_Id", n_wall_moving},
273 {"Wall_Force_Direction", 2}}}}}};
274}
anonymous_namespace{main.cpp}::generateCircularGrid
std::vector< PackedParticle > generateCircularGrid(int ncols, int nrows, double x_lo, double y_lo, double R, double gap)
Definition main.cpp:61

References geom::createGeomObject(), util::Point::d_x, util::Point::d_y, generateCircularGrid(), inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::BCDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), util::harmonicMean(), util::io::print(), material::toE(), and material::toGE().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildInputJson() [8/8]

json anonymous_namespace{main.cpp}::buildInputJson ( const std::string &  output_path_for_deck,
const std::filesystem::path &  mesh_file_1,
const std::filesystem::path &  mesh_file_2,
double  final_time,
size_t  num_steps,
bool  zero_ic,
double  mesh_size_in,
double  horizon_in,
bool  damping_on,
double  eps_n,
bool  two_particle_test,
bool  file_mesh,
double  beta_n_factor,
const std::string &  damping_law = "com_and_node",
const std::string &  friction_law = "coulomb_simple",
bool  friction_on = false,
double  friction_mu = -1.,
double  ic_vx = 0.,
bool  bottom_patch_bc = false,
const std::string &  mpi_strategy = "auto" 
)

Definition at line 91 of file main.cpp.

104 {
105
106 const std::vector<double> center = {0.0, 0.0, 0.0};
107 const double R1 = 0.001;
108 const double R2 = 0.001;
109 const double mesh_size =
110 mesh_size_in > 0. ? mesh_size_in : std::min(R1, R2) / 5.0;
111 const double horizon = horizon_in > 0. ? horizon_in : 3.0 * mesh_size;
112 const double particle_dist = 0.001;
113
114 const double poisson1 = 0.25;
115 const double rho1 = 1200.0;
116 const double K1 = 2.16e+7;
117 const double E1 = material::toE(K1, poisson1);
118 const double G1 = material::toGE(E1, poisson1);
119 const double Gc1 = 50.0;
120
121 const double poisson2 = 0.25;
122 const double rho2 = 1200.0;
123 const double K2 = 2.16e+7;
124 const double E2 = material::toE(K2, poisson2);
125 const double G2 = material::toGE(E2, poisson2);
126 const double Gc2 = 50.0;
127
128 const double R_contact_factor = 0.95;
129 const double Kn_11 = 18.0 * util::harmonicMean(K1, K1) / (M_PI * std::pow(horizon, 5));
130 const double Kn_22 = 18.0 * util::harmonicMean(K2, K2) / (M_PI * std::pow(horizon, 5));
131 const double Kn_12 = 18.0 * util::harmonicMean(K1, K2) / (M_PI * std::pow(horizon, 5));
132 const double friction_coeff = 0.5;
133
134 std::vector<double> p1_center = center;
135 std::vector<double> p2_center = center;
136
137 const double H0_drop = particle_dist; // paper drop height, 1 mm
138 // Table 2 energy: remaining gap = horizon, rest of H0 already in IC velocity.
139 const double top_gap = two_particle_test ? horizon : particle_dist;
140
141 const size_t dt_out_n = num_steps / 10;
142 auto modelDeckJson = inp::ModelDeck::getExampleJson(2, final_time, num_steps,
143 "finite_difference", "central_difference",
144 true, 2, "Multi_Particle", 0);
145 modelDeckJson["MPI_Strategy"] = mpi_strategy;
146
147 auto outputDeckJson = inp::OutputDeck::getExampleJson("vtu", output_path_for_deck,
148 std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage_Z", "Damage", "Particle_ID"}),
149 dt_out_n, 2, true, "zlib", true, 1, "", true);
150
151 auto bcDeckJson = inp::BCDeck::getExampleJson(0, 1, 1, true, util::Point(0, -10, 0));
152
153 if (bottom_patch_bc) {
154 // Particle 0 centered at (R1,R1); fix only a bottom strip (not the whole grain).
155 geom::GeomData patch;
156 patch.d_geomName = "rectangle";
157 patch.d_geomParams = {R1 - 1.1 * R1, -0.1 * R1, 0., R1 + 1.1 * R1, 0.35 * R1,
158 0.};
159 json j_geom;
160 geom::writeGeometry(j_geom, patch);
161 bcDeckJson["Displacement_BC"]["Set_1"] = json{
162 {"Particle_List", json::array({0})},
163 {"Region", {{"Geometry", j_geom}}},
164 {"Direction", json::array({1, 2})},
165 {"Zero_Displacement", true}};
166 } else {
167 bcDeckJson["Displacement_BC"]["Set_1"] =
168 inp::BCBaseDeck::getExampleJson("Displacement_BC", false, geom::GeomData(),
169 {0}, {}, "", {}, "", {}, {1, 2}, true, "",
170 {});
171 }
172
173 std::vector<double> ic_vel = {0.0, 0.0, 0.0};
174 if (!zero_ic) {
175 const double fallen =
176 two_particle_test ? (H0_drop - top_gap) : (particle_dist - horizon);
177 if (fallen > 0.)
178 ic_vel[1] = -std::sqrt(2.0 * std::abs(-10.0) * fallen);
179 }
180 ic_vel[0] = ic_vx;
181 bcDeckJson["IC"]["Set_1"] = inp::BCBaseDeck::getExampleJson("IC", false, geom::GeomData(),
182 {1}, {}, "", {}, "", {},
183 {}, false, "Constant_Velocity", ic_vel);
184
185 auto pDeckJson = json({});
186
187 std::vector<geom::GeomData> pGeomVec(2);
188 pGeomVec[0].d_geomName = "circle";
189 pGeomVec[0].d_geomParams = {R1, p1_center[0], p1_center[1], p1_center[2]};
190 pGeomVec[1].d_geomName = "circle";
191 pGeomVec[1].d_geomParams = {R2, p2_center[0], p2_center[1], p2_center[2]};
192
193 pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
194
195 const std::string f1 = mesh_file_1.string();
196 const std::string f2 = mesh_file_2.string();
197 json meshSet1, meshSet2;
198 if (file_mesh) {
199 meshSet1 = json({{"File", f1}});
200 meshSet2 = json({{"File", f2}});
201 } else {
202 meshSet1 = json({{"File", f1},
203 {"CreateMesh",
204 {{"Flag", true},
205 {"Info", "gmsh_builtin_mesh"},
206 {"Mesh_Size", mesh_size},
207 {"Write_Mesh_File", true}}}});
208 meshSet2 = json({{"File", f2},
209 {"CreateMesh",
210 {{"Flag", true},
211 {"Info", "gmsh_builtin_mesh"},
212 {"Mesh_Size", mesh_size},
213 {"Write_Mesh_File", true}}}});
214 }
215 pDeckJson["Mesh"] = json({{"Sets", 2}, {"Set_1", meshSet1}, {"Set_2", meshSet2}});
216
217 auto pMatJson = inp::ParticleDeck::getParticleMaterialExampleJson(2);
218 pMatJson["Set_1"] = inp::MaterialDeck::getExampleJson("PDState", false, horizon,
219 0, rho1, K1, G1, Gc1, true, 1);
220 pMatJson["Set_2"] = inp::MaterialDeck::getExampleJson("PDState", false, horizon,
221 0, rho2, K2, G2, Gc2, true, 1);
222 pDeckJson["Material"] = pMatJson;
223
224 auto pContactJson = inp::ParticleDeck::getParticleContactExampleJson(2);
225 // v0.1.0 circ_damp: Friction_On: false leaves μ = 0 (coeff is not read).
226 // Nonzero μ with the flag off still applies a tangential force in PairForce
227 // and walks the falling particle sideways (Table 2 test 4).
228 const double mu =
229 (friction_mu >= 0.) ? friction_mu
230 : (two_particle_test ? 0. : friction_coeff);
231 json contact_base = inp::ContactPairDeck::getExampleJson(
232 R_contact_factor, true, damping_on, friction_on, Kn_11, eps_n, mu,
233 1.0, beta_n_factor, 1.0, 0.0, K1);
234
235 pContactJson["Set_1_1"] = contact_base;
236 pContactJson["Set_1_1"]["Kn"] = Kn_11;
237 pContactJson["Set_1_1"]["K"] = K1;
238 pContactJson["Set_1_2"] = contact_base;
239 pContactJson["Set_1_2"]["Kn"] = Kn_12;
240 pContactJson["Set_1_2"]["K"] = util::harmonicMean(K1, K2);
241 pContactJson["Set_2_2"] = contact_base;
242 pContactJson["Set_2_2"]["Kn"] = Kn_22;
243 pContactJson["Set_2_2"]["K"] = K2;
244 pContactJson["Damping_Law"] = damping_law;
245 pContactJson["Friction_Law"] = friction_law;
246 pDeckJson["Contact"] = pContactJson;
247
248 pDeckJson["Neighbor"] = two_particle_test
249 ? inp::PNeighborDeck::getExampleJson("simple_all", 5.0, 1, 0.5)
250 : inp::PNeighborDeck::getExampleJson("simple_all", 10.0, 40, 0.5);
251
252 auto pGenJson = inp::PGenDeck::getExampleJson("From_File");
253 pGenJson["Random_Rotation"] = false;
254 pGenJson["Data"]["N"] = 2;
255 pGenJson["Data"]["0"] = {
256 {"x", R1}, {"y", R1}, {"z", 0.0},
257 {"theta", 0.0}, {"s", 1.0},
258 {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0}
259 };
260 pGenJson["Data"]["1"] = {
261 {"x", R1}, {"y", 2.0 * R1 + R2 + top_gap}, {"z", 0.0},
262 {"theta", two_particle_test ? M_PI / 2.0 : M_PI}, {"s", 1.0},
263 {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}
264 };
265 pDeckJson["Particle_Generation"] = pGenJson;
266
267 auto j = json({{"Model", modelDeckJson},
268 {"Output", outputDeckJson},
269 {"Force_BC", bcDeckJson["Force_BC"]},
270 {"Displacement_BC", bcDeckJson["Displacement_BC"]},
271 {"IC", bcDeckJson["IC"]},
272 {"Particle", pDeckJson["Particle"]},
273 {"Mesh", pDeckJson["Mesh"]},
274 {"Material", pDeckJson["Material"]},
275 {"Contact", pDeckJson["Contact"]},
276 {"Neighbor", pDeckJson["Neighbor"]},
277 {"Particle_Generation", pDeckJson["Particle_Generation"]}});
278 if (two_particle_test)
279 j["Test"] = json{{"Test_Name", "two_particle"}};
280 // bottom_patch_bc keeps PD force on particle 0 (deformable base).
281 return j;
282}
geom::writeGeometry
void writeGeometry(json &j, const geom::GeomData &geomData)
Definition geomObjectsUtil.cpp:1363
inp
Collection of methods and database related to input.
Definition pairForce.h:20

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::ContactPairDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::BCDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), util::harmonicMean(), material::toE(), material::toGE(), and geom::writeGeometry().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ buildTraskInputJson()

json anonymous_namespace{main.cpp}::buildTraskInputJson ( const std::string &  output_path,
const std::filesystem::path &  mesh_plate,
double  W,
double  H,
double  notch_half,
double  notch_w,
double  notch_depth,
double  mesh_size,
double  horizon,
double  rho,
double  E,
double  K,
double  G,
double  Gc,
double  v_impact,
double  final_time,
size_t  num_steps 
)

Definition at line 369 of file main.cpp.

373 {
374 auto model = inp::ModelDeck::getExampleJson(2, final_time, num_steps, "finite_difference",
375 "central_difference", true, 2, "Single_Particle",
376 0);
377 // Trask §5: broken bond → weight 0 (no force). No self-contact.
378 model["Self_Contact"] = "none";
379 model["Bond_Break"] = "tension"; // literature PMB: break in tension only
380
381 auto output = inp::OutputDeck::getExampleJson(
382 "vtu", output_path,
383 std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage", "Damage_Z",
384 "Particle_ID", "Fixity"}),
385 std::max<size_t>(1, num_steps / 10), 1, false, "zlib", true, num_steps, "", false);
386
387 // Trask §6.2: top left/right of notches u=<0,0>; drive between notches
388 // u=<0,-v t>; sides and bottom free (collar bond-break not yet implemented).
389 const double drive_h = std::max(2.0 * mesh_size, 0.002);
390 const std::vector<double> drive_strip{-notch_half, 0.5 * H - drive_h, 0., notch_half,
391 0.55 * H, 0.};
392 json outer = fixedTopOuterBC(W, H, notch_half);
393 json bc_disp = {
394 {"Sets", 4},
395 {"Set_1", outer["Set_1"]},
396 {"Set_2", outer["Set_2"]},
397 {"Set_3",
398 {{"Particle_List", std::vector<size_t>{0}},
399 {"Region", {{"Geometry", {{"Type", "rectangle"}, {"Parameters", drive_strip}}}}},
400 {"Direction", std::vector<size_t>{1}},
401 {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
402 {"Spatial_Function", {{"Type", "constant"}}},
403 {"Zero_Displacement", true}}},
404 {"Set_4",
405 {{"Particle_List", std::vector<size_t>{0}},
406 {"Region", {{"Geometry", {{"Type", "rectangle"}, {"Parameters", drive_strip}}}}},
407 {"Direction", std::vector<size_t>{2}},
408 {"Time_Function",
409 {{"Type", "linear"}, {"Parameters", std::vector<double>{-v_impact}}}},
410 {"Spatial_Function", {{"Type", "constant"}}}}}};
411
412 geom::GeomData plate =
413 sillingNotchedPlateGeom(W, H, notch_half, notch_w, notch_depth);
414 auto particle = inp::ParticleDeck::getParticleGeomExampleJson({plate});
415 json mesh = {{"Sets", 1}, {"Set_1", meshSetJson(mesh_plate, mesh_size)}};
416 auto material = inp::ParticleDeck::getParticleMaterialExampleJson(1);
417 material["Set_1"] = inp::MaterialDeck::getExampleJson(
418 "PMBBond", false, horizon, 0, rho, K, G, Gc, true, 0, E);
419 // Trask/Silling ω≡1 (ConstInfluence default a0=dim+1 would mis-scale c)
420 material["Set_1"]["Influence_Function"] = {
421 {"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
422
423 return json{{"Model", model},
424 {"Output", output},
425 {"Displacement_BC", bc_disp},
426 {"Particle", particle},
427 {"Mesh", mesh},
428 {"Material", material}};
429}
anonymous_namespace{main.cpp}::sillingNotchedPlateGeom
geom::GeomData sillingNotchedPlateGeom(double W, double H, double notch_half, double notch_w, double notch_depth)
Definition main.cpp:296

References fixedTopOuterBC(), inp::ModelDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), meshSetJson(), and sillingNotchedPlateGeom().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ checkWallSpringOnly()

Probe anonymous_namespace{main.cpp}::checkWallSpringOnly ( PeriDEMModel &  dem,
geom::GeomObject &  wall_geom 
)

Definition at line 164 of file main.cpp.

164 {
165 Probe p;
166 if (!dem.d_contact_p || !dem.d_contact_p->d_wallContact ||
167 !dem.d_contact_p->d_pairForce)
168 throw std::runtime_error("analytical wall: contact path not set up");
169
170 const auto &contact = dem.d_particleDeck_p->d_contactDeck.getContact(0, 1);
171 const double Kn = contact.d_Kn;
172 const double Rc = contact.d_contactR;
173
174 for (auto &f : dem.d_f)
175 f = util::Point();
176 dem.d_contact_p->d_wallContact->apply(dem, dem.d_contact_p->d_pairForce.get(),
177 /*use_node_damping*/ false);
178
179 for (size_t i = 0; i < dem.d_x.size(); ++i) {
180 if (dem.d_ptId[i] != 0)
181 continue;
182 geom::WallContactHit hit;
183 if (!wall_geom.wallContactQuery(dem.d_x[i], hit) || !hit.active)
184 continue;
185 if (!(hit.signed_gap < Rc))
186 continue;
187
188 ++p.n_contact_nodes;
189 p.min_gap = std::min(p.min_gap, hit.signed_gap);
190
191 auto scalar = Kn * (hit.signed_gap - Rc) * dem.d_vol[i];
192 if (scalar > 0.)
193 scalar = 0.;
194 const util::Point f_ref = scalar * (-1. * hit.outward_n);
195 const double err = (dem.d_f[i] - f_ref).length();
196 p.max_force_err = std::max(p.max_force_err, err);
197 p.max_force_ref = std::max(p.max_force_ref, f_ref.length());
198 }
199 if (p.n_contact_nodes == 0)
200 p.min_gap = 0.;
201 return p;
202}
PeriDEMModel::d_contact_p
std::unique_ptr< contact::Contact > d_contact_p
Definition periDEMModel.h:75
data::ModelData::d_f
std::vector< util::Point > d_f
Total force on the nodes.
Definition modelData.h:757
data::ModelData::d_particleDeck_p
std::shared_ptr< inp::ParticleDeck > d_particleDeck_p
Pointer to deck holding particle related data.
Definition modelData.h:637
data::ModelData::d_vol
std::vector< double > d_vol
Nodal volumes.
Definition modelData.h:760
geom::GeomObject::wallContactQuery
virtual bool wallContactQuery(const util::Point &x, WallContactHit &hit) const
Closest-point / signed-gap query for analytical wall contact. Default: unsupported (returns false,...
Definition geomObjects.h:345
util
Collection of methods useful in simulation.
Definition constants.h:14
anonymous_namespace{main.cpp}::Probe
Definition main.cpp:157
anonymous_namespace{main.cpp}::Probe::max_force_err
double max_force_err
Definition main.cpp:160
anonymous_namespace{main.cpp}::Probe::max_force_ref
double max_force_ref
Definition main.cpp:161
anonymous_namespace{main.cpp}::Probe::min_gap
double min_gap
Definition main.cpp:159
anonymous_namespace{main.cpp}::Probe::n_contact_nodes
size_t n_contact_nodes
Definition main.cpp:158
geom::WallContactHit
Result of a wall-contact query against a geom (analytical walls).
Definition geomObjects.h:112
geom::WallContactHit::outward_n
util::Point outward_n
Definition geomObjects.h:115
util::Point::length
double length() const
Computes the Euclidean length of the vector.
Definition point.h:124

References geom::WallContactHit::active, PeriDEMModel::d_contact_p, data::ModelData::d_f, data::ModelData::d_particleDeck_p, data::ModelData::d_ptId, data::ModelData::d_vol, data::ModelData::d_x, util::Point::length(), anonymous_namespace{main.cpp}::Probe::max_force_err, anonymous_namespace{main.cpp}::Probe::max_force_ref, anonymous_namespace{main.cpp}::Probe::min_gap, anonymous_namespace{main.cpp}::Probe::n_contact_nodes, geom::WallContactHit::outward_n, geom::WallContactHit::signed_gap, and geom::GeomObject::wallContactQuery().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ contactPairJson()

json anonymous_namespace{main.cpp}::contactPairJson ( double  R_contact_factor,
bool  damping_on,
bool  friction_on,
double  Kn,
double  beta_n_eps,
double  friction_coeff,
double  Kn_factor,
double  beta_n_factor 
)

Definition at line 181 of file main.cpp.

182 {
183 json j;
184 j["Contact_Radius_Factor"] = R_contact_factor;
185 if (Kn < 1e-20) {
186 j["Kn"] = 0.0;
187 } else {
188 j["Kn"] = Kn;
189 }
190 j["Damping_On"] = damping_on;
191 j["Epsilon"] = beta_n_eps;
192 j["Friction_On"] = friction_on;
193 j["Friction_Coeff"] = friction_coeff;
194 j["Kn_Factor"] = Kn_factor;
195 j["Beta_n_Factor"] = damping_on ? beta_n_factor : 0.0;
196 return j;
197}

Referenced by buildInputJson().

Here is the caller graph for this function:

◆ directoryPathWithTrailingSep()

std::string anonymous_namespace{main.cpp}::directoryPathWithTrailingSep ( const std::filesystem::path &  dir)

OutputDeck concatenates Path + "log.txt" without inserting a separator; Path must end with one.

Definition at line 34 of file main.cpp.

34 {
35 namespace fs = std::filesystem;
36 fs::path n = fs::absolute(dir).lexically_normal();
37 std::string s = n.string();
38 if (!s.empty() && s.back() != '/' && s.back() != '\\')
39 s += fs::path::preferred_separator;
40 return s;
41}

Referenced by main(), and runLaw().

Here is the caller graph for this function:

◆ doesParticleIntersect()

bool anonymous_namespace{main.cpp}::doesParticleIntersect ( const PackedParticle &  p,
const std::vector< PackedParticle > &  existing,
const std::vector< double > &  rect,
double  padding 
)

Definition at line 50 of file main.cpp.

51 {
52 const double p_rect_lox = p.x - p.r;
53 const double p_rect_loy = p.y - p.r;
54 const double p_rect_hix = p.x + p.r;
55 const double p_rect_hiy = p.y + p.r;
56 if (p_rect_lox < rect[0] + padding || p_rect_loy < rect[1] + padding || p_rect_hix > rect[3] - padding ||
57 p_rect_hiy > rect[4] - padding)
58 return true;
59
60 for (const auto &q : existing) {
61 const double dx = p.x - q.x;
62 const double dy = p.y - q.y;
63 const double dz = p.z - q.z;
64 const double dist = std::sqrt(dx * dx + dy * dy + dz * dz);
65 if (dist <= p.r + q.r + padding)
66 return true;
67 }
68 return false;
69}
anonymous_namespace{main.cpp}::PackedParticle::r
double r
Definition main.cpp:40
anonymous_namespace{main.cpp}::PackedParticle::y
double y
Definition main.cpp:40
anonymous_namespace{main.cpp}::PackedParticle::z
double z
Definition main.cpp:40
anonymous_namespace{main.cpp}::PackedParticle::x
double x
Definition main.cpp:40

References anonymous_namespace{main.cpp}::PackedParticle::r, anonymous_namespace{main.cpp}::PackedParticle::x, anonymous_namespace{main.cpp}::PackedParticle::y, and anonymous_namespace{main.cpp}::PackedParticle::z.

Referenced by generateParticleLocations().

Here is the caller graph for this function:

◆ dumpFinalNodal()

void anonymous_namespace{main.cpp}::dumpFinalNodal ( const data::ModelData &  data,
const std::filesystem::path &  out_dir 
)

Definition at line 74 of file main.cpp.

75 {
76 const size_t n = data.d_u.size();
77 std::vector<double> buf(6 * n, 0.);
78 for (size_t i = 0; i < n; ++i) {
79 if (!ownsNode(data, i))
80 continue;
81 buf[6 * i + 0] = data.d_u[i].d_x;
82 buf[6 * i + 1] = data.d_u[i].d_y;
83 buf[6 * i + 2] = data.d_u[i].d_z;
84 buf[6 * i + 3] = data.d_v[i].d_x;
85 buf[6 * i + 4] = data.d_v[i].d_y;
86 buf[6 * i + 5] = data.d_v[i].d_z;
87 }
88 if (util::parallel::mpiSize() > 1)
89 MPI_Allreduce(MPI_IN_PLACE, buf.data(), static_cast<int>(buf.size()),
90 MPI_DOUBLE, MPI_SUM, util::parallel::mpiComm());
91 if (util::parallel::mpiRank() != 0)
92 return;
93
94 std::filesystem::create_directories(out_dir / "nodal");
95 {
96 std::ofstream xr(out_dir / "nodal" / "x_ref.bin", std::ios::binary);
97 const uint32_t nn = static_cast<uint32_t>(n);
98 xr.write(reinterpret_cast<const char *>(&nn), sizeof(nn));
99 for (size_t i = 0; i < n; ++i) {
100 const double p[3] = {data.d_xRef[i].d_x, data.d_xRef[i].d_y,
101 data.d_xRef[i].d_z};
102 xr.write(reinterpret_cast<const char *>(p), sizeof(p));
103 }
104 }
105 {
106 std::ofstream os(out_dir / "nodal" / "final_uv.bin", std::ios::binary);
107 const char magic[4] = {'P', 'D', 'U', 'V'};
108 const uint32_t step32 = static_cast<uint32_t>(data.currentStep());
109 const uint32_t nn = static_cast<uint32_t>(n);
110 os.write(magic, 4);
111 os.write(reinterpret_cast<const char *>(&step32), sizeof(step32));
112 os.write(reinterpret_cast<const char *>(&nn), sizeof(nn));
113 os.write(reinterpret_cast<const char *>(buf.data()),
114 static_cast<std::streamsize>(buf.size() * sizeof(double)));
115 }
116}
anonymous_namespace{main.cpp}::ownsNode
bool ownsNode(const data::ModelData &data, size_t i)
Definition main.cpp:55
data
Definition contact.h:20
util::parallel::mpiSize
int mpiSize()
Get size (number) of processors.
Definition parallelUtil.cpp:88
util::parallel::mpiRank
int mpiRank()
get rank (id) of this processor
Definition parallelUtil.cpp:96
util::parallel::mpiComm
MPI_Comm mpiComm()
Get MPI comm.
Definition parallelUtil.cpp:104

References util::parallel::mpiComm(), util::parallel::mpiRank(), util::parallel::mpiSize(), and ownsNode().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ fitCrackFromPhi()

CrackFit anonymous_namespace{main.cpp}::fitCrackFromPhi ( const PeriDEMModel &  dem,
double  tip_x,
double  tip_y,
double  outward_sign,
double  phi_cut,
double  band_x,
double  band_y,
double  exclude_r,
double  bin_h,
bool  use_bond_count = false 
)

Definition at line 187 of file main.cpp.

189 {
190 CrackFit fit;
191 std::vector<double> xs, ys;
192
193 const int nbin = std::max(4, static_cast<int>(std::ceil((band_x - exclude_r) / bin_h)));
194 for (int b = 0; b < nbin; ++b) {
195 const double r0 = exclude_r + b * bin_h;
196 const double r1 = std::min(band_x, r0 + bin_h);
197 double best_phi = phi_cut;
198 double best_dx = 0., best_dy = 0.;
199 bool found = false;
200 for (size_t i = 0; i < dem.d_xRef.size(); ++i) {
201 if (dem.d_ptId[i] != 0)
202 continue;
203 const double phi = nodeDamageForFit(dem, i, use_bond_count);
204 if (phi < best_phi)
205 continue;
206 const double dx_ref = outward_sign * (dem.d_xRef[i].d_x - tip_x);
207 const double dy_ref = tip_y - dem.d_xRef[i].d_y;
208 if (dx_ref < r0 || dx_ref >= r1)
209 continue;
210 if (dy_ref < -2.0 * bin_h || dy_ref > band_y)
211 continue;
212 const double dx = dem.d_xRef[i].d_x - tip_x;
213 const double dy = tip_y - dem.d_xRef[i].d_y;
214 best_phi = phi;
215 best_dx = dx;
216 best_dy = dy;
217 found = true;
218 }
219 if (found) {
220 xs.push_back(best_dx);
221 ys.push_back(best_dy);
222 }
223 }
224 fit.n_pts = xs.size();
225 if (fit.n_pts < 4)
226 return fit;
227 double mx = 0., my = 0.;
228 for (size_t i = 0; i < fit.n_pts; ++i) {
229 mx += xs[i];
230 my += ys[i];
231 }
232 mx /= static_cast<double>(fit.n_pts);
233 my /= static_cast<double>(fit.n_pts);
234 fit.mean_dx = mx;
235 fit.mean_dy = my;
236 double cxx = 0., cxy = 0., cyy = 0.;
237 for (size_t i = 0; i < fit.n_pts; ++i) {
238 const double x = xs[i] - mx;
239 const double y = ys[i] - my;
240 cxx += x * x;
241 cxy += x * y;
242 cyy += y * y;
243 }
244 const double trace = cxx + cyy;
245 const double det = cxx * cyy - cxy * cxy;
246 const double tmp = std::sqrt(std::max(0., 0.25 * trace * trace - det));
247 const double l1 = 0.5 * trace + tmp;
248 double vx = cxy;
249 double vy = l1 - cxx;
250 if (vx * vx + vy * vy < 1.e-30) {
251 vx = l1 - cyy;
252 vy = cxy;
253 }
254 const double angle_from_horiz =
255 std::atan2(std::abs(vy), std::abs(vx)) * 180. / M_PI;
256 fit.angle_to_notch_deg = 90. - angle_from_horiz;
257 return fit;
258}
anonymous_namespace{main.cpp}::nodeDamageForFit
double nodeDamageForFit(const PeriDEMModel &dem, size_t i, bool use_bond_count)
Definition main.cpp:168
util::det
double det(const std::vector< std::vector< double > > &m)
Computes the determinant of matrix.
Definition matrix.cpp:75
anonymous_namespace{main.cpp}::CrackFit
Definition main.cpp:179
anonymous_namespace{main.cpp}::CrackFit::angle_to_notch_deg
double angle_to_notch_deg
Definition main.cpp:181
anonymous_namespace{main.cpp}::CrackFit::mean_dx
double mean_dx
Definition main.cpp:182
anonymous_namespace{main.cpp}::CrackFit::n_pts
size_t n_pts
Definition main.cpp:180
anonymous_namespace{main.cpp}::CrackFit::mean_dy
double mean_dy
Definition main.cpp:183

References anonymous_namespace{main.cpp}::CrackFit::angle_to_notch_deg, data::ModelData::d_ptId, data::ModelData::d_xRef, anonymous_namespace{main.cpp}::CrackFit::mean_dx, anonymous_namespace{main.cpp}::CrackFit::mean_dy, anonymous_namespace{main.cpp}::CrackFit::n_pts, and nodeDamageForFit().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ fixedTopOuterBC()

json anonymous_namespace{main.cpp}::fixedTopOuterBC ( double  W,
double  H,
double  notch_half 
)

Definition at line 339 of file main.cpp.

339 {
340 const double fix_h = 0.08 * H;
341 return json{
342 {"Set_1",
343 {{"Particle_List", std::vector<size_t>{0}},
344 {"Region",
345 {{"Geometry",
346 {{"Type", "rectangle"},
347 {"Parameters",
348 std::vector<double>{-0.55 * W, 0.5 * H - fix_h, 0., -notch_half - 0.02 * W,
349 0.55 * H, 0.}}}}}},
350 {"Direction", std::vector<size_t>{1, 2}},
351 {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
352 {"Spatial_Function", {{"Type", "constant"}}},
353 {"Zero_Displacement", true}}},
354 {"Set_2",
355 {{"Particle_List", std::vector<size_t>{0}},
356 {"Region",
357 {{"Geometry",
358 {{"Type", "rectangle"},
359 {"Parameters",
360 std::vector<double>{notch_half + 0.02 * W, 0.5 * H - fix_h, 0., 0.55 * W, 0.55 * H,
361 0.}}}}}},
362 {"Direction", std::vector<size_t>{1, 2}},
363 {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
364 {"Spatial_Function", {{"Type", "constant"}}},
365 {"Zero_Displacement", true}}}};
366}

Referenced by buildBhatInputJson(), and buildTraskInputJson().

Here is the caller graph for this function:

◆ generateCircularGrid()

std::vector< PackedParticle > anonymous_namespace{main.cpp}::generateCircularGrid ( int  ncols,
int  nrows,
double  x_lo,
double  y_lo,
double  R,
double  gap 
)

Definition at line 61 of file main.cpp.

62 {
63 std::vector<PackedParticle> particles;
64 particles.reserve(static_cast<size_t>(ncols * nrows));
65 const double step = 2.0 * R + gap;
66 const double x0 = x_lo + gap + R;
67 const double y0 = y_lo + gap + R;
68 for (int j = 0; j < nrows; ++j) {
69 for (int i = 0; i < ncols; ++i) {
70 PackedParticle p;
71 p.x = x0 + static_cast<double>(i) * step;
72 p.y = y0 + static_cast<double>(j) * step;
73 p.z = 0.0;
74 p.r = R;
75 p.theta = 0.0;
76 particles.push_back(p);
77 }
78 }
79 return particles;
80}
anonymous_namespace{main.cpp}::PackedParticle
Definition main.cpp:38
anonymous_namespace{main.cpp}::PackedParticle::theta
double theta
Definition main.cpp:40

Referenced by buildInputJson(), and main().

Here is the caller graph for this function:

◆ generateParticleLocations()

std::vector< PackedParticle > anonymous_namespace{main.cpp}::generateParticleLocations ( const std::vector< double > &  in_rect,
double  max_y,
double  mesh_size,
double  R,
int  N_target,
double  padding,
std::mt19937 &  gen 
)

Random row-wise packing (same spirit as problem_setup.py method 1) with std::mt19937(30). One characteristic radius R; zones 0–7 pick among eight geometry types.

Definition at line 75 of file main.cpp.

77 {
78 std::vector<PackedParticle> particles;
79 std::uniform_real_distribution<double> u_r(-0.1 * R, 0.1 * R);
80 /* Do not use uniform_int_distribution: libstdc++ vs libc++ map mt19937
81 * differently, so geom_id diverges while x,y,theta stay the same. */
82
83 const double check_r = R;
84 const int rows = static_cast<int>((max_y - in_rect[1]) / (2.0 * check_r));
85 const double rect_L = in_rect[3] - in_rect[0];
86 const int cols = static_cast<int>(rect_L / (2.0 * check_r));
87
88 int counter = 0;
89 double x_old = in_rect[0];
90 double x_old_right = in_rect[3];
91 double y_old = in_rect[1];
92 const double cz = 0.0;
93
94 std::vector<double> cy_accptd;
95 cy_accptd.push_back(y_old);
96
97 std::vector<double> row_rads_prev{R, R};
98 for (int i = 0; i < rows; ++i) {
99 if (i > 0)
100 y_old = maxElem(cy_accptd) + maxElem(row_rads_prev);
101
102 std::vector<double> row_rads{R};
103
104 if (y_old + padding + maxElem(row_rads) >= max_y)
105 break;
106
107 int num_p_cols = 0;
108 int j = 0;
109 while (true) {
110 if (num_p_cols > cols - 1 || j > 100 * N_target)
111 break;
112 if (counter >= N_target)
113 break;
114
115 if (j == 0) {
116 x_old = in_rect[0];
117 x_old_right = in_rect[3];
118 }
119
120 const int p_zone = static_cast<int>(gen() % 8);
121 const double r0 = R;
122 double r = r0 + u_r(gen);
123
124 double cx = 0., cy = 0.;
125 if (i % 2 == 0) {
126 std::uniform_real_distribution<double> rph(-0.1 * r0, 0.05 * r0);
127 std::uniform_real_distribution<double> rpv(-0.05 * r0, 0.05 * r0);
128 const double cx0 = x_old_right - padding - r;
129 cx = cx0 - rph(gen);
130 cy = y_old + padding + r + rpv(gen);
131 } else {
132 std::uniform_real_distribution<double> rph(-0.05 * r0, 0.1 * r0);
133 std::uniform_real_distribution<double> rpv(-0.05 * r0, 0.05 * r0);
134 const double cx0 = x_old + padding + r;
135 cx = cx0 + rph(gen);
136 cy = y_old + padding + r + rpv(gen);
137 }
138
139 PackedParticle trial{p_zone, cx, cy, cz, r, 0.0};
140 if (!doesParticleIntersect(trial, particles, in_rect, padding)) {
141 std::uniform_real_distribution<double> orient(0.0, 2.0 * M_PI);
142 trial.theta = orient(gen);
143 particles.push_back(trial);
144 row_rads.push_back(trial.r);
145 cy_accptd.push_back(cy);
146 if (i % 2 == 0)
147 x_old_right = cx - trial.r;
148 else
149 x_old = cx + trial.r;
150
151 ++counter;
152 num_p_cols++;
153 }
154 ++j;
155 }
156 row_rads_prev = std::move(row_rads);
157 }
158 return particles;
159}
anonymous_namespace{main.cpp}::doesParticleIntersect
bool doesParticleIntersect(const PackedParticle &p, const std::vector< PackedParticle > &existing, const std::vector< double > &rect, double padding)
Definition main.cpp:50
anonymous_namespace{main.cpp}::maxElem
double maxElem(const std::vector< double > &v)
Definition main.cpp:43

References doesParticleIntersect(), and maxElem().

Referenced by buildInputJson().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ is3DGeometry()

bool anonymous_namespace{main.cpp}::is3DGeometry ( const std::string &  g)

Definition at line 52 of file main.cpp.

52 {
53 return g == "sphere" || g == "ellipsoid" || g == "cube" || g == "cuboid" || g == "cylinder" ||
54 g == "sphere_minus_sphere" || g == "cuboid_minus_cuboid" || g == "open_cuboid_channel_3d";
55}

Referenced by buildInputJson().

Here is the caller graph for this function:

◆ isAcceptableGeometryName()

bool anonymous_namespace{main.cpp}::isAcceptableGeometryName ( const std::string &  g)

Definition at line 57 of file main.cpp.

57 {
58 for (const std::string &n : geom::getAcceptableGeometries()) {
59 if (n == g)
60 return true;
61 }
62 return false;
63}

References geom::getAcceptableGeometries().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ KnFromBulk()

double anonymous_namespace{main.cpp}::KnFromBulk ( double  Ka,
double  Kb,
double  horizon 
)

Definition at line 199 of file main.cpp.

199 {
200 return 18.0 * util::harmonicMean(Ka, Kb) / (M_PI * std::pow(horizon, 5));
201}

References util::harmonicMean().

Referenced by buildInputJson().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ maxAbsForceFromCsv()

double anonymous_namespace{main.cpp}::maxAbsForceFromCsv ( const std::filesystem::path &  csv)

Definition at line 284 of file main.cpp.

284 {
285 std::ifstream is(csv);
286 if (!is)
287 return 0.;
288 std::string line;
289 std::getline(is, line);
290 double max_abs = 0.;
291 while (std::getline(is, line)) {
292 if (line.empty())
293 continue;
294 const auto c2 = line.rfind(',');
295 if (c2 == std::string::npos)
296 continue;
297 max_abs = std::max(max_abs, std::abs(std::stod(line.substr(c2 + 1))));
298 }
299 return max_abs;
300}

Referenced by main().

Here is the caller graph for this function:

◆ maxElem()

double anonymous_namespace{main.cpp}::maxElem ( const std::vector< double > &  v)

Definition at line 43 of file main.cpp.

43 {
44 double m = -std::numeric_limits<double>::infinity();
45 for (double x : v)
46 m = std::max(m, x);
47 return m;
48}

Referenced by generateParticleLocations().

Here is the caller graph for this function:

◆ measure()

CrackMetrics anonymous_namespace{main.cpp}::measure ( PeriDEMModel &  dem,
double  x_crack,
double  band 
)

Definition at line 147 of file main.cpp.

147 {
148 CrackMetrics m;
149 for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
150 const auto &xi = dem.d_xRef[i];
151 for (size_t k = 0; k < dem.d_neighPd[i].size(); ++k) {
152 if (!dem.d_fracture_p->getBondState(i, k))
153 continue;
154 const size_t j = dem.d_neighPd[i][k];
155 const auto &xj = dem.d_xRef[j];
156 if ((xi.d_x - x_crack) * (xj.d_x - x_crack) >= 0.)
157 continue;
158 const double gap = (dem.d_x[j] - dem.d_x[i]).length();
159 m.min_gap = std::min(m.min_gap, gap);
160 m.max_gap = std::max(m.max_gap, gap);
161 ++m.n_pairs;
162 }
163 }
164 if (m.n_pairs == 0)
165 m.min_gap = 0.;
166
167 for (size_t i = 0; i < dem.d_xRef.size(); ++i) {
168 const double x = dem.d_xRef[i].d_x;
169 const double fx = dem.d_f[i].d_x * dem.d_vol[i];
170 if (x < x_crack && x > x_crack - band)
171 m.force_x_left += fx;
172 if (x >= x_crack && x < x_crack + band)
173 m.force_x_right += fx;
174 }
175 return m;
176}
anonymous_namespace{main.cpp}::CrackMetrics
Definition main.cpp:139
anonymous_namespace{main.cpp}::CrackMetrics::n_pairs
size_t n_pairs
Definition main.cpp:142
anonymous_namespace{main.cpp}::CrackMetrics::force_x_left
double force_x_left
Definition main.cpp:143
anonymous_namespace{main.cpp}::CrackMetrics::min_gap
double min_gap
Definition main.cpp:140
anonymous_namespace{main.cpp}::CrackMetrics::max_gap
double max_gap
Definition main.cpp:141
anonymous_namespace{main.cpp}::CrackMetrics::force_x_right
double force_x_right
Definition main.cpp:144

References data::ModelData::d_f, data::ModelData::d_fracture_p, data::ModelData::d_neighPd, data::ModelData::d_vol, data::ModelData::d_x, data::ModelData::d_xRef, anonymous_namespace{main.cpp}::CrackMetrics::force_x_left, anonymous_namespace{main.cpp}::CrackMetrics::force_x_right, anonymous_namespace{main.cpp}::CrackMetrics::max_gap, anonymous_namespace{main.cpp}::CrackMetrics::min_gap, and anonymous_namespace{main.cpp}::CrackMetrics::n_pairs.

Referenced by runLaw().

Here is the caller graph for this function:

◆ meshSetJson()

json anonymous_namespace{main.cpp}::meshSetJson ( const std::filesystem::path &  f,
double  mesh_size 
)

Definition at line 260 of file main.cpp.

260 {
261 return json{{"File", f.string()},
262 {"CreateMesh",
263 {{"Flag", true},
264 {"Info", "gmsh_builtin_mesh"},
265 {"Mesh_Size", mesh_size},
266 {"Write_Mesh_File", true}}}};
267}

Referenced by buildBhatInputJson(), buildImpactInputJson(), and buildTraskInputJson().

Here is the caller graph for this function:

◆ minGrainWallGap()

double anonymous_namespace{main.cpp}::minGrainWallGap ( PeriDEMModel &  dem,
geom::GeomObject &  wall_geom 
)

Definition at line 204 of file main.cpp.

204 {
205 double ming = 1.e300;
206 size_t n = 0;
207 for (size_t i = 0; i < dem.d_x.size(); ++i) {
208 if (dem.d_ptId[i] != 0)
209 continue;
210 geom::WallContactHit hit;
211 if (!wall_geom.wallContactQuery(dem.d_x[i], hit) || !hit.active)
212 continue;
213 ming = std::min(ming, hit.signed_gap);
214 ++n;
215 }
216 return n ? ming : 0.;
217}

References geom::WallContactHit::active, data::ModelData::d_ptId, data::ModelData::d_x, geom::WallContactHit::signed_gap, and geom::GeomObject::wallContactQuery().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ near()

bool anonymous_namespace{main.cpp}::near ( double  a,
double  b,
double  tol = 1.e-12 
)

Definition at line 25 of file main.cpp.

25 {
26 return std::abs(a - b) <= tol * (1. + std::abs(a) + std::abs(b));
27}

Referenced by main().

Here is the caller graph for this function:

◆ nodeDamageForFit()

double anonymous_namespace{main.cpp}::nodeDamageForFit ( const PeriDEMModel &  dem,
size_t  i,
bool  use_bond_count 
)

Definition at line 168 of file main.cpp.

168 {
169 if (use_bond_count) {
170 if (!dem.d_phiBond.empty() && i < dem.d_phiBond.size())
171 return static_cast<double>(dem.d_phiBond[i]);
172 return nodePhi(dem, i);
173 }
174 if (!dem.d_phi.empty() && i < dem.d_phi.size())
175 return static_cast<double>(dem.d_phi[i]);
176 return nodePhi(dem, i);
177}
data::ModelData::d_phi
std::vector< float > d_phi
Damage function at the nodes (volume-weighted, Silling 2000)
Definition modelData.h:834
data::ModelData::d_phiBond
std::vector< float > d_phiBond
Damage as broken-bond count fraction (Bhattacharya & Lipton 2023)
Definition modelData.h:837
anonymous_namespace{main.cpp}::nodePhi
double nodePhi(const PeriDEMModel &dem, size_t i)
Definition main.cpp:156

References data::ModelData::d_phi, data::ModelData::d_phiBond, and nodePhi().

Referenced by fitCrackFromPhi(), and main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ nodePhi()

double anonymous_namespace{main.cpp}::nodePhi ( const PeriDEMModel &  dem,
size_t  i 
)

Definition at line 156 of file main.cpp.

156 {
157 const auto &nb = dem.d_neighPd[i];
158 if (nb.empty())
159 return 0.;
160 size_t n_br = 0;
161 for (size_t k = 0; k < nb.size(); ++k) {
162 if (dem.d_fracture_p->getBondState(i, k))
163 ++n_br;
164 }
165 return static_cast<double>(n_br) / static_cast<double>(nb.size());
166}

References data::ModelData::d_fracture_p, and data::ModelData::d_neighPd.

Referenced by main(), and nodeDamageForFit().

Here is the caller graph for this function:

◆ notchVoidBoxes()

std::vector< std::vector< double > > anonymous_namespace{main.cpp}::notchVoidBoxes ( double  H,
double  notch_half,
double  notch_w,
double  notch_depth,
double  z_lo,
double  z_hi 
)

Definition at line 284 of file main.cpp.

286 {
287 const double hw = 0.5 * notch_w;
288 const double y_tip = 0.5 * H - notch_depth;
289 const double y_hi = 0.5 * H + 1.0e-9;
290 return {{-notch_half - hw, y_tip, z_lo, -notch_half + hw, y_hi, z_hi},
291 {notch_half - hw, y_tip, z_lo, notch_half + hw, y_hi, z_hi}};
292}

Referenced by buildImpactInputJson().

Here is the caller graph for this function:

◆ ownsNode()

bool anonymous_namespace{main.cpp}::ownsNode ( const data::ModelData &  data,
size_t  i 
)

Definition at line 55 of file main.cpp.

55 {
56 const int rank = util::parallel::mpiRank();
57 if (data.d_pdDofMpi) {
58 if (i >= data.d_pdNodePartition.size())
59 return false;
60 return static_cast<int>(data.d_pdNodePartition[i]) == rank;
61 }
62 for (const auto *p : data.d_particlesListTypeAll) {
63 const size_t i0 = p->d_globStart;
64 const size_t i1 = i0 + p->getNumNodes();
65 if (i < i0 || i >= i1)
66 continue;
67 if (p->isWall())
68 return rank == 0;
69 return particle::isLocallyOwned(*p);
70 }
71 return rank == 0;
72}
particle::isLocallyOwned
bool isLocallyOwned(const BaseParticle &p)
True if this rank updates / assembles forces for the particle. Walls are replicated on every rank....
Definition particleMpi.cpp:28

References particle::isLocallyOwned(), and util::parallel::mpiRank().

Referenced by dumpFinalNodal(), and anonymous_namespace{main.cpp}::GrainContactProbe::sampleMetric().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ packedAxisBounds()

PackedBounds anonymous_namespace{main.cpp}::packedAxisBounds ( const std::vector< PackedParticle > &  packed)

Definition at line 166 of file main.cpp.

166 {
167 PackedBounds b;
168 b.min_x = std::numeric_limits<double>::infinity();
169 b.max_x = -std::numeric_limits<double>::infinity();
170 b.min_y = std::numeric_limits<double>::infinity();
171 b.max_y = -std::numeric_limits<double>::infinity();
172 for (const auto &p : packed) {
173 b.min_x = std::min(b.min_x, p.x - p.r);
174 b.max_x = std::max(b.max_x, p.x + p.r);
175 b.min_y = std::min(b.min_y, p.y - p.r);
176 b.max_y = std::max(b.max_y, p.y + p.r);
177 }
178 return b;
179}
anonymous_namespace{main.cpp}::PackedBounds
Definition main.cpp:162
anonymous_namespace{main.cpp}::PackedBounds::max_y
double max_y
Definition main.cpp:163
anonymous_namespace{main.cpp}::PackedBounds::min_x
double min_x
Definition main.cpp:163
anonymous_namespace{main.cpp}::PackedBounds::min_y
double min_y
Definition main.cpp:163
anonymous_namespace{main.cpp}::PackedBounds::max_x
double max_x
Definition main.cpp:163

References anonymous_namespace{main.cpp}::PackedBounds::max_x, anonymous_namespace{main.cpp}::PackedBounds::max_y, anonymous_namespace{main.cpp}::PackedBounds::min_x, and anonymous_namespace{main.cpp}::PackedBounds::min_y.

Referenced by buildInputJson().

Here is the caller graph for this function:

◆ runLaw()

LawResult anonymous_namespace{main.cpp}::runLaw ( const std::string &  self_contact,
const std::filesystem::path &  run_dir,
double  L,
double  mesh_size,
double  horizon,
unsigned  n_threads 
)

Definition at line 184 of file main.cpp.

185 {
186 namespace fs = std::filesystem;
187 const fs::path out_dir = run_dir / "out";
188 const fs::path inp_dir = run_dir / "inp";
189 fs::create_directories(out_dir);
190 fs::create_directories(inp_dir);
191
192 const fs::path mesh = inp_dir / "mesh_rect.msh";
193 // Time fields unused for prescribed kinematics; keep Output_Interval valid.
194 auto input_json =
195 buildInputJson(directoryPathWithTrailingSep(out_dir), mesh, L, mesh_size, horizon, 1.0e-6,
196 2, self_contact);
197 {
198 std::ofstream os(inp_dir / "input.json");
199 os << input_json.dump(2);
200 }
201
202 auto deck = std::make_shared<inp::Input>(input_json);
203 PeriDEMModel dem(deck);
204 dem.init();
205
206 const double x_crack = 0.5 * L;
207 const size_t n_broken = applyPrecrack(dem, x_crack);
208 if (n_broken == 0)
209 throw std::runtime_error(self_contact + ": precrack broke no bonds");
210
211 const double Rc = dem.d_maxContactR > 0. ? dem.d_maxContactR : 0.95 * mesh_size;
212 const double h = dem.d_hMax > 0. ? dem.d_hMax : mesh_size;
213 util::io::print(std::format("{}: nodes={}, broken_bonds={}, h={}, Rc={}, horizon={}\n",
214 self_contact, dem.d_x.size(), n_broken, h, Rc, horizon));
215
216 // Open: separate faces so current gap >> Rc (no self-contact yet).
217 const double u_open = 2.0 * Rc;
218 applyRigidOpenClose(dem, x_crack, u_open);
219 dem.computeForces();
220 const auto open_m = measure(dem, x_crack, 2.0 * horizon);
221 if (!(open_m.min_gap > 1.5 * Rc))
222 throw std::runtime_error(std::format("{}: open did not separate (min_gap={})", self_contact,
223 open_m.min_gap));
224
225 // Close into the window where reference_gap (natural R = r0) and
226 // broken_bond_kn (natural R = Rc) disagree: Rc < R < typical r0.
227 // Mid-horizon neighbor distances are ~O(horizon) >> Rc = 0.95 h.
228 const double u_close = -0.5 * (Rc + 0.5 * horizon);
229 applyRigidOpenClose(dem, x_crack, u_close);
230 dem.computeForces();
231 const auto close_m = measure(dem, x_crack, 2.0 * horizon);
232 const double force_amp =
233 0.5 * (std::abs(close_m.force_x_left) + std::abs(close_m.force_x_right));
234
235 util::io::print(std::format("{}: open min_gap={}, close min_gap={}, force_amp={}\n", self_contact,
236 open_m.min_gap, close_m.min_gap, force_amp));
237
238 dem.close();
239
240 LawResult r;
241 r.open_min_gap = open_m.min_gap;
242 r.close_min_gap = close_m.min_gap;
243 r.close_force_amp = force_amp;
244 return r;
245}
anonymous_namespace{main.cpp}::directoryPathWithTrailingSep
std::string directoryPathWithTrailingSep(const std::filesystem::path &dir)
Definition main.cpp:34
anonymous_namespace{main.cpp}::buildInputJson
json buildInputJson(const std::string &output_path, const std::filesystem::path &mesh_file, double L, double mesh_size, double horizon, double final_time, size_t num_steps, const std::string &self_contact)
Definition main.cpp:44
anonymous_namespace{main.cpp}::measure
CrackMetrics measure(PeriDEMModel &dem, double x_crack, double band)
Definition main.cpp:147
anonymous_namespace{main.cpp}::applyRigidOpenClose
void applyRigidOpenClose(PeriDEMModel &dem, double x_crack, double u_right)
Definition main.cpp:130
anonymous_namespace{main.cpp}::applyPrecrack
size_t applyPrecrack(PeriDEMModel &dem, double x_crack)
Definition main.cpp:114

References applyPrecrack(), applyRigidOpenClose(), buildInputJson(), PeriDEMModel::close(), anonymous_namespace{main.cpp}::LawResult::close_force_amp, anonymous_namespace{main.cpp}::LawResult::close_min_gap, PeriDEMModel::computeForces(), data::ModelData::d_hMax, data::ModelData::d_maxContactR, data::ModelData::d_x, directoryPathWithTrailingSep(), PeriDEMModel::init(), measure(), anonymous_namespace{main.cpp}::LawResult::open_min_gap, and util::io::print().

Referenced by main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ segmentCrossesVertical()

bool anonymous_namespace{main.cpp}::segmentCrossesVertical ( double  x0,
double  y0,
double  x1,
double  y1,
double  x_line,
double  y_lo,
double  y_hi 
)

Definition at line 48 of file main.cpp.

49 {
50 if ((x0 - x_line) * (x1 - x_line) >= 0.)
51 return false;
52 const double t = (x_line - x0) / (x1 - x0);
53 if (t < 0. || t > 1.)
54 return false;
55 const double y = y0 + t * (y1 - y0);
56 return y >= y_lo && y <= y_hi;
57}

Referenced by applyNotches(), applyNotchSlots(), and main().

Here is the caller graph for this function:

◆ sillingNotchedPlateGeom()

geom::GeomData anonymous_namespace{main.cpp}::sillingNotchedPlateGeom ( double  W,
double  H,
double  notch_half,
double  notch_w,
double  notch_depth 
)

Definition at line 296 of file main.cpp.

297 {
298 const double x0 = -0.5 * W, y0 = -0.5 * H, x1 = 0.5 * W, y1 = 0.5 * H;
299 const double hw = 0.5 * notch_w;
300 const double y_tip = y1 - notch_depth;
301 geom::GeomData plate;
302 plate.d_geomName = "complex";
303 plate.d_geomComplexInfo = {
304 std::vector<std::string>{"rectangle", "rectangle", "rectangle"},
305 std::vector<std::string>{"plus", "minus", "minus"}};
306 plate.d_geomParams = {
307 // outer plate
308 x0, y0, 0., x1, y1, 0.,
309 // left notch cutout (opens at top)
310 -notch_half - hw, y_tip, 0., -notch_half + hw, y1 + 1.0e-6, 0.,
311 // right notch cutout
312 notch_half - hw, y_tip, 0., notch_half + hw, y1 + 1.0e-6, 0.};
313 return plate;
314}

References geom::GeomData::d_geomComplexInfo, geom::GeomData::d_geomName, and geom::GeomData::d_geomParams.

Referenced by buildTraskInputJson().

Here is the caller graph for this function:

◆ uniformMeshSetJson()

json anonymous_namespace{main.cpp}::uniformMeshSetJson ( const std::filesystem::path &  f,
double  mesh_size,
const std::vector< std::vector< double > > &  voids = {} 
)

Definition at line 272 of file main.cpp.

273 {}) {
274 json cm = {{"Flag", true},
275 {"Info", "uniform"},
276 {"Mesh_Size", mesh_size},
277 {"Write_Mesh_File", true}};
278 if (!voids.empty())
279 cm["Void_Regions"] = voids;
280 return json{{"File", f.string()}, {"CreateMesh", cm}};
281}

Referenced by buildImpactInputJson().

Here is the caller graph for this function:

◆ writeLocations()

void anonymous_namespace{main.cpp}::writeLocations ( const std::filesystem::path &  csv,
const std::vector< PackedParticle > &  packed 
)

Definition at line 276 of file main.cpp.

277 {
278 std::ofstream os(csv);
279 os << "i, x, y, z, r, o\n";
280 for (const auto &p : packed)
281 os << std::format("0, {:.6f}, {:.6f}, {:.6f}, {:.6f}, 0.000000\n", p.x, p.y, p.z, p.r);
282}

Referenced by main().

Here is the caller graph for this function: