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)
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	makeDeck (const std::string &model, size_t dim, bool plane_strain, int influence, double L, double h, double horizon, double E, double nu, double Gc, const std::string &out)
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, bool deck_only=false)
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 41 of file main.cpp.

                                 {
  return geometry::breakBondsCrossingVerticalLines(
      *dem.d_fracture_p, dem.d_xRef, dem.d_neighPd, dem.d_ptId, 0,
      {x_left, x_right}, y_lo, y_hi);
}

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

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 57 of file main.cpp.

                                         {
  return geometry::breakBondsCrossingVerticalLines(
      *dem.d_fracture_p, dem.d_xRef, dem.d_neighPd, dem.d_ptId, 0,
      {-notch_half, notch_half}, y_tip, y_top);
}

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

Referenced by main().

Here is the call graph for this function:

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 48 of file main.cpp.

                                                 {
  return geometry::breakBondsInSlots(*dem.d_fracture_p, dem.d_xRef,
                                     dem.d_neighPd, dem.d_ptId, 0,
                                     {-notch_half, notch_half}, notch_w, y_lo,
                                     y_hi);
}

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

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 137 of file main.cpp.

                                                        {
  size_t n_broken = 0;
  for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
    const auto &xi = dem.d_xRef[i];
    for (size_t k = 0; k < dem.d_neighPd[i].size(); ++k) {
      const size_t j = dem.d_neighPd[i][k];
      const auto &xj = dem.d_xRef[j];
      if ((xi.d_x - x_crack) * (xj.d_x - x_crack) < 0.) {
        dem.d_fracture_p->setBondState(i, k, true);
        ++n_broken;
      }
    }
  }
  return n_broken;
}

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 153 of file main.cpp.

                                                                            {
  for (size_t i = 0; i < dem.d_xRef.size(); ++i) {
    const double ux = (dem.d_xRef[i].d_x >= x_crack) ? u_right : 0.;
    dem.d_u[i] = util::Point(ux, 0., 0.);
    dem.d_v[i] = util::Point();
    dem.d_x[i] = dem.d_xRef[i] + dem.d_u[i];
  }
}

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 241 of file main.cpp.

                                                                         {
  const double x0 = -0.5 * W, y0 = -0.5 * H, x1 = 0.5 * W, y1 = 0.5 * H;
  const double hw = 0.5 * notch_w;
  const double y_tip = y1 - notch_depth;
  const double y_top = y1 + 1.0e-6;
  geom::GeomData plate;
  plate.d_geomName = "complex";
  plate.d_geomComplexInfo = {
      std::vector<std::string>{"rectangle", "triangle", "triangle"},
      std::vector<std::string>{"plus", "minus", "minus"}};
  plate.d_geomParams = {
      // outer plate
      x0, y0, 0., x1, y1, 0.,
      // left V: top-left, top-right, tip
      -notch_half - hw, y_top, 0., -notch_half + hw, y_top, 0., -notch_half, y_tip, 0.,
      // right V
      notch_half - hw, y_top, 0., notch_half + hw, y_top, 0., notch_half, y_tip, 0.};
  return plate;
}

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 570 of file main.cpp.

                                                                              {
  auto model = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "velocity_verlet"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
  model["Self_Contact"] = "reference_gap"; // Bhat 2023 §4.4
  model["Bond_Break"] = "absolute_stretch";  // Bhat §3.1: |s| > |s0|
  model["Wall_Contact"] = "meshed";
 
  auto output = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage", "Damage_Bond",
                                "Damage_Z", "Particle_ID"})},
           {"Output_Interval", std::max<size_t>(1, num_steps / 10)},
           {"Debug", 1},
           {"Perform_FE_Out", false},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", num_steps},
           {"Tag_PP", ""},
           {"PVD_Collection", false}});
 
  json ic = {
      {"Sets", 1},
      {"Set_1",
       {{"Particle_List", std::vector<size_t>{1}},
        {"Constant_Velocity",
         {{"Velocity_Vector", std::vector<double>{0., -v_impact, 0.}}}}}}};
 
  // Paper: only the outer top ligaments are fixed; impactor is free.
  json outer = fixedTopOuterBC(W, H, notch_half);
  json bc_disp = {{"Sets", 2}, {"Set_1", outer["Set_1"]}, {"Set_2", outer["Set_2"]}};
 
  // Bhat specimen: OCC/Gmsh nonconvex plate with V-notch cutouts (§5, Fig. 4).
  geom::GeomData plate =
      bhatVNotchedPlateGeom(W, H, notch_half, notch_w, notch_depth);
  geom::GeomData impactor;
  impactor.d_geomName = "rectangle";
  impactor.d_geomParams = {-0.5 * Iw, -0.5 * Ih, 0., 0.5 * Iw, 0.5 * Ih, 0.};
  auto particle = inp::ParticleDeck::getParticleGeomExampleJson({plate, impactor});
 
  json mesh = {{"Sets", 2},
               {"Set_1", meshSetJson(mesh_plate, mesh_size)},
               {"Set_2", meshSetJson(mesh_impactor, mesh_size)}};
 
  // Bhat Table 1, constant micromodulus, ν = 1/3 (bond-based, 2D).
  const double nu_bhat = 1.0 / 3.0;
  const double c_bhat =
      6.0 * E / (M_PI * std::pow(horizon, 3.0) * (1.0 - nu_bhat));
  const double s0_bhat = std::sqrt(4.0 * M_PI * Gc / (9.0 * E * horizon));
 
  auto material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
  auto mat_bhat = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PMBBond"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"E", E},
           {"Compute_From_Classical", false},
           {"Influence_Function", json{{"Type", 0}}}});
  mat_bhat["Bond_Potential_Params"] = std::vector<double>{c_bhat, s0_bhat};
  mat_bhat["Influence_Function"] = {{"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
  material["Set_1"] = mat_bhat;
  // Free deformable striker (paper). Unbreakable elastic so the contact face
  // delivers bulk impulse; plate remains PMB with paper (c, s0).
  auto mat_strike = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDElasticBond"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", 0.},
           {"E", E},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  mat_strike["Influence_Function"] = {{"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
  material["Set_2"] = mat_strike;
 
  // Absolute Rc from mesh_size (not Factor×hMin).
  const double Kn_bhat = 18.0 * K / (M_PI * std::pow(horizon, 4.0));
  const double Rc_abs = Rc_factor * mesh_size;
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius", Rc_abs},
               {"Kn", Kn_bhat},
               {"Damping_On", false},
               {"Epsilon", 1.0},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 0.0},
               {"K", K}});
  contact_base["Kn"] = Kn_bhat;
  json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
  contact["Set_1_1"] = contact_base;
  contact["Set_1_2"] = contact_base;
  contact["Set_2_2"] = contact_base;
  contact["Damping_Law"] = "off";            // Bhat §6: β_d = 0
  contact["Friction_Law"] = "coulomb_simple"; // μ = 0 below
  contact["Set_1_1"]["Friction_Coefficient"] = 0.;
  contact["Set_1_2"]["Friction_Coefficient"] = 0.;
  contact["Set_2_2"]["Friction_Coefficient"] = 0.;
 
  const double cy_imp = 0.5 * H + 0.5 * Ih + gap;
  auto pgen = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pgen["Random_Rotation"] = false;
  pgen["Data"]["N"] = 2;
  pgen["Data"]["0"] = json{{"x", 0.},     {"y", 0.},     {"z", 0.},
                           {"theta", 0.}, {"s", 1.},     {"geom_id", 0},
                           {"mat_id", 0}, {"contact_id", 0}};
  pgen["Data"]["1"] = json{{"x", 0.},     {"y", cy_imp}, {"z", 0.},
                           {"theta", 0.}, {"s", 1.},     {"geom_id", 1},
                           {"mat_id", 1}, {"contact_id", 1}};
 
  return json{{"Model", model},
              {"Output", output},
              {"Displacement_BC", bc_disp},
              {"IC", ic},
              {"Particle", particle},
              {"Mesh", mesh},
              {"Material", material},
              {"Contact", contact},
              {"Neighbor", inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 5.0},
               {"Search_Interval", 1},
               {"Near_Bd_Nodes_Tol", 0.5}})},
              {"Particle_Generation", pgen}};
}

References bhatVNotchedPlateGeom(), geom::GeomData::d_geomName, geom::GeomData::d_geomParams, fixedTopOuterBC(), inp::ContactPairDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::OutputDeck::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 377 of file main.cpp.

                                                 {
  // Silling 2003 EMU KW: explicit central difference; brittle-microelastic PMB;
  // broken bond force = 0; plate boundaries load-free (no Displacement_BC);
  // rigid impactor of finite mass.
  // thickness > 0 runs his actual 3D case (200 x 100 x 9 grid, 9 mm plate).
  const bool dim3 = thickness > 0.;
  const size_t dim = dim3 ? 3 : 2;
  // Element-node connectivity is only used for strain/stress output and does not
  // support hexahedra, which is what the 3D structured grid produces.
  auto model = inp::ModelDeck::getExampleJson(
      json{{"Dimension", dim},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", !dim3},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
  model["Self_Contact"] = "none";
  model["Bond_Break"] = "tension";
  model["Wall_Contact"] = "meshed";
  // Silling Fig. 2: rigid cylinder of 1.57 kg. In 2D the analogue is the mass
  // per unit thickness of his 9 mm plate, M/t = 1.57/0.009 kg/m.
  model["Rigid_Particles"] = json::array(
      {json{{"Id", 1}, {"Mass", dim3 ? 1.57 : 1.57 / 0.009}}});
 
  auto output = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage", "Damage_Bond",
                                "Damage_Z", "Particle_ID"})},
           {"Output_Interval", std::max<size_t>(1, num_steps / 10)},
           {"Debug", 1},
           {"Perform_FE_Out", false},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", num_steps},
           {"Tag_PP", ""},
           {"PVD_Collection", false}});
 
  json ic = {
      {"Sets", 1},
      {"Set_1",
       {{"Particle_List", std::vector<size_t>{1}},
        {"Constant_Velocity",
         {{"Velocity_Vector", std::vector<double>{0., -v_impact, 0.}}}}}}};
 
  // The impactor is rigid but of finite mass (Model.Rigid_Particles above), so
  // it decelerates on contact instead of being driven through the plate.
  json bc_disp = {{"Sets", 1},
                  {"Set_1",
                   {{"Particle_List", std::vector<size_t>{1}},
                    {"Direction", std::vector<size_t>{1}},
                    {"Time_Function",
                     {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
                    {"Spatial_Function", {{"Type", "constant"}}},
                    {"Zero_Displacement", true}}}};
 
  // Plate is the full box on a structured grid with the two 1.5 mm notch slots
  // removed from the mesh, so they are real gaps (Fig. 2) rather than material
  // whose bonds have been cut.
  geom::GeomData plate;
  geom::GeomData impactor;
  if (dim3) {
    // Keep both bodies centred on z = 0 so particle generation places them
    // consistently (it positions objects by their centre).
    plate.d_geomName = "cuboid";
    plate.d_geomParams = {-0.5 * W, -0.5 * H, -0.5 * thickness,
                          0.5 * W,  0.5 * H,  0.5 * thickness};
    // Silling Fig. 4: the cylinder axis is along the impact direction and its
    // flat end face strikes the plate edge. Params: radius, centre of the
    // beginning cross-section, and the vector to the end cross-section.
    impactor.d_geomName = "cylinder";
    impactor.d_geomParams = {0.5 * Iw, 0., -0.5 * Ih, 0., 0., Ih, 0.};
  } else {
    plate.d_geomName = "rectangle";
    plate.d_geomParams = {-0.5 * W, -0.5 * H, 0., 0.5 * W, 0.5 * H, 0.};
    // In-plane section of that same cylinder: Iw wide, Ih long, flat face down.
    impactor.d_geomName = "rectangle";
    impactor.d_geomParams = {-0.5 * Iw, -0.5 * Ih, 0., 0.5 * Iw, 0.5 * Ih, 0.};
  }
  auto particle = inp::ParticleDeck::getParticleGeomExampleJson({plate, impactor});
 
  const auto voids =
      notchVoidBoxes(H, notch_half, notch_w, notch_depth,
                     dim3 ? -0.5 * thickness - 1.0e-9 : -1.0e-9,
                     dim3 ? 0.5 * thickness + 1.0e-9 : 1.0e-9);
  // Plate on Silling's equally spaced structured grid. In 2D the impactor is a
  // rectangle so it goes on the same grid; in 3D it is a cylinder, which the
  // uniform grid cannot represent, so it is meshed by Gmsh at the same spacing.
  json mesh = {{"Sets", 2},
               {"Set_1", uniformMeshSetJson(mesh_plate, mesh_size, voids)},
               {"Set_2", dim3 ? meshSetJson(mesh_impactor, mesh_size)
                              : uniformMeshSetJson(mesh_impactor, mesh_size)}};
 
  // Impactor is steel like the plate but unbreakable, and its nodal forces are
  // replaced by the rigid-body acceleration each step, so its bond stiffness
  // never enters the solution.
  auto material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
  material["Set_1"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PMBBond"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"E", E},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  material["Set_2"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDElasticBond"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", 0.},
           {"E", E},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  material["Set_1"]["Influence_Function"] = {
      {"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
  material["Set_2"]["Influence_Function"] = {
      {"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
 
  // Absolute Contact_Radius from the requested mesh size (not Factor×hMin).
  const double Rc_abs = Rc_factor * mesh_size;
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius", Rc_abs},
               {"Kn", Kn},
               {"Damping_On", false},
               {"Epsilon", 1.0},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 0.0},
               {"K", K}});
  contact_base["Kn"] = Kn;
  json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
  contact["Set_1_1"] = contact_base;
  contact["Set_1_2"] = contact_base;
  contact["Set_2_2"] = contact_base;
  contact["Damping_Law"] = "off";
  contact["Friction_Law"] = "coulomb_simple";
 
  const double cy_imp = 0.5 * H + 0.5 * Ih + gap; // flat face `gap` above the edge
  auto pgen = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pgen["Random_Rotation"] = false;
  pgen["Data"]["N"] = 2;
  pgen["Data"]["0"] = json{{"x", 0.},
                           {"y", 0.},
                           {"z", 0.},
                           {"theta", 0.},
                           {"s", 1.},
                           {"geom_id", 0},
                           {"mat_id", 0},
                           {"contact_id", 0}};
  pgen["Data"]["1"] = json{{"x", 0.},
                           {"y", cy_imp},
                           {"z", 0.},
                           {"theta", 0.},
                           {"s", 1.},
                           {"geom_id", 1},
                           {"mat_id", 1},
                           {"contact_id", 1}};
 
  // Plate has no Displacement_BC (Silling: load-free all around); the only
  // Displacement_BC drives the rigid impactor.
  return json{{"Model", model},
              {"Output", output},
              {"Displacement_BC", bc_disp},
              {"IC", ic},
              {"Particle", particle},
              {"Mesh", mesh},
              {"Material", material},
              {"Contact", contact},
              {"Neighbor", inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 5.0},
               {"Search_Interval", 1},
               {"Near_Bd_Nodes_Tol", 0.5}})},
              {"Particle_Generation", pgen}};
}

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::ContactPairDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::OutputDeck::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.

                                                          {
 
  const double s = 0.001;
  const double mesh_size = s / 5.0;
  const double horizon = 3.0 * mesh_size;
  const double particle_dist = 0.001;
 
  const double poisson1 = 0.25;
  const double rho1 = 1200.0;
  const double K1 = 2.16e+7;
  const double E1 = material::toE(K1, poisson1);
  const double G1 = material::toGE(E1, poisson1);
  const double Gc1 = 50.0;
 
  const double poisson2 = 0.25;
  const double rho2 = 1200.0;
  const double K2 = 2.16e+7;
  const double E2 = material::toE(K2, poisson2);
  const double G2 = material::toGE(E2, poisson2);
  const double Gc2 = 50.0;
 
  const double R_contact_factor = 0.95;
  const double Kn_11 = util::normalContactStiffness(K1, K1, horizon);
  const double Kn_22 = util::normalContactStiffness(K2, K2, horizon);
  const double Kn_12 = util::normalContactStiffness(K1, K2, horizon);
  const double beta_n_eps = 0.9;
  const double friction_coeff = 0.5;
  const double beta_n_factor = 100.0;
 
  const int model_dim = is3DGeometry(geomName) ? 3 : 2;
 
  /* Reference shapes at origin for both mesh groups (same as twop_circ_inbuilt). */
  const util::Point origin(0.0, 0.0, 0.0);
  const std::vector<double> gp_ref = geom::exampleGeomParams(geomName, origin, s);
 
  /* Same integration window and IC gap as twop_circ_inbuilt (stable dt). */
  const double final_time = 0.0001; // 0.012;
  const size_t num_steps = 1000; //36000;
  const size_t dt_out_n = num_steps / 4; // 10;
 
  auto modelDeckJson = inp::ModelDeck::getExampleJson(
      json{{"Dimension", static_cast<size_t>(model_dim)},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
 
  auto outputDeckJson = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path_for_deck},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage_Z", "Damage", "Particle_ID"})},
           {"Output_Interval", dt_out_n},
           {"Debug", 2},
           {"Perform_FE_Out", true},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", 1},
           {"Tag_PP", ""},
           {"PVD_Collection", true}});
 
  auto bcDeckJson = inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 1},
           {"IC_Sets", 1},
           {"Gravity", util::Point(0, -10, 0).toVec()}});
 
  bcDeckJson["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "Displacement_BC"},
           {"Particle_List", {0}},
           {"Direction", {1, 2}},
           {"Zero_Displacement", true}});
 
  const double free_fall_dist = particle_dist - horizon;
  const double free_fall_vel = -std::sqrt(2.0 * std::abs(-10.0) * std::max(1e-30, free_fall_dist));
  bcDeckJson["IC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "IC"},
           {"Particle_List", {1}},
           {"IC_Type", "Constant_Velocity"},
           {"IC_Vector", {0.0, free_fall_vel, 0.0}}});
 
  auto pDeckJson = json({});
 
  std::vector<geom::GeomData> pGeomVec(2);
  pGeomVec[0].d_geomName = geomName;
  pGeomVec[0].d_geomParams = gp_ref;
  pGeomVec[1].d_geomName = geomName;
  pGeomVec[1].d_geomParams = gp_ref;
 
  pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
 
  const std::string f1 = mesh_file_1.string();
  const std::string f2 = mesh_file_2.string();
  auto meshSet1 = json({{"File", f1},
                        {"CreateMesh",
                         {{"Flag", true},
                          {"Info", "gmsh_builtin_mesh"},
                          {"Mesh_Size", mesh_size},
                          {"Write_Mesh_File", true}}}});
  auto meshSet2 = json({{"File", f2},
                        {"CreateMesh",
                         {{"Flag", true},
                          {"Info", "gmsh_builtin_mesh"},
                          {"Mesh_Size", mesh_size},
                          {"Write_Mesh_File", true}}}});
  pDeckJson["Mesh"] = json({{"Sets", 2}, {"Set_1", meshSet1}, {"Set_2", meshSet2}});
 
  auto pMatJson = inp::ParticleDeck::getParticleMaterialExampleJson(2);
  pMatJson["Set_1"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho1},
           {"K", K1},
           {"G", G1},
           {"Gc", Gc1},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  pMatJson["Set_2"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho2},
           {"K", K2},
           {"G", G2},
           {"Gc", Gc2},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  pDeckJson["Material"] = pMatJson;
 
  auto pContactJson = inp::ParticleDeck::getParticleContactExampleJson(2);
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", R_contact_factor},
               {"Kn", Kn_11},
               {"Damping_On", false},
               {"Epsilon", beta_n_eps},
               {"Friction_On", false},
               {"Friction_Coeff", friction_coeff},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", beta_n_factor},
               {"K", 0.0}});
 
  pContactJson["Set_1_1"] = contact_base;
  pContactJson["Set_1_1"]["Kn"] = Kn_11;
  pContactJson["Set_1_2"] = contact_base;
  pContactJson["Set_1_2"]["Kn"] = Kn_12;
  pContactJson["Set_2_2"] = contact_base;
  pContactJson["Set_2_2"]["Kn"] = Kn_22;
  pDeckJson["Contact"] = pContactJson;
 
  pDeckJson["Neighbor"] = inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 10.0},
               {"Search_Interval", 40},
               {"Near_Bd_Nodes_Tol", 0.5}});
 
  auto pGenJson = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pGenJson["Data"]["N"] = 2;
 
  /*
   * Placement matches twop_circ: lower particle corner at (s,s[,s]); upper at
   * (s, 2*s + particle_dist + s, [s]) so vertical gap is particle_dist between
   * outer extents (same as R1=R2=s in twop_circ).
   */
  if (model_dim == 2) {
    pGenJson["Data"]["0"] = {
        {"x", s}, {"y", s}, {"z", 0.0},
        {"theta", 0.0}, {"s", 1.0},
        {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
        {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0}
    };
    pGenJson["Data"]["1"] = {
        {"x", s}, {"y", 2.0 * s + particle_dist + s}, {"z", 0.0},
        {"theta", M_PI}, {"s", 1.0},
        {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
        {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}
    };
  } else {
    pGenJson["Data"]["0"] = {
        {"x", s}, {"y", s}, {"z", s},
        {"theta", 0.0}, {"s", 1.0},
        {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
        {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0}
    };
    pGenJson["Data"]["1"] = {
        {"x", s}, {"y", 2.0 * s + particle_dist + s}, {"z", s},
        {"theta", M_PI}, {"s", 1.0},
        {"ax", 0.0}, {"ay", 0.0}, {"az", 1.0},
        {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}
    };
  }
  pDeckJson["Particle_Generation"] = pGenJson;
 
  return json({{"Model", modelDeckJson},
               {"Output", outputDeckJson},
               {"Force_BC", bcDeckJson["Force_BC"]},
               {"Displacement_BC", bcDeckJson["Displacement_BC"]},
               {"IC", bcDeckJson["IC"]},
               {"Particle", pDeckJson["Particle"]},
               {"Mesh", pDeckJson["Mesh"]},
               {"Material", pDeckJson["Material"]},
               {"Contact", pDeckJson["Contact"]},
               {"Neighbor", pDeckJson["Neighbor"]},
               {"Particle_Generation", pDeckJson["Particle_Generation"]}});
}

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

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.

                                      {
  const double R = 0.001;
  const double mesh_size = R / 4.0;
  const double horizon = 3.0 * mesh_size;
  const double h_est = 0.7 * mesh_size;
  const double Rc_est = 0.95 * h_est;
 
  // Start inside PP Rc and near the floor; soft Kn (below) keeps contact stable.
  const double gap_pp = 0.5 * Rc_est;
  const double cy = R + 0.35 * Rc_est;
  const double cx0 = -(R + 0.5 * gap_pp);
  const double cx1 = +(R + 0.5 * gap_pp);
 
  const double wall_half = 3.5 * R;
  const double wall_thick = std::max(2.0 * mesh_size, horizon);
 
  geom::GeomData grain;
  grain.d_geomName = "circle";
  grain.d_geomParams = {R, 0., 0., 0.};
  geom::GeomData floor;
  floor.d_geomName = "rectangle";
  floor.d_geomParams = {-wall_half, -wall_thick, 0., wall_half, 0., 0.};
 
  const double poisson = 0.25;
  const double rho = 1200.0;
  const double K = 2.16e7;
  const double E = material::toE(K, poisson);
  const double G = material::toGE(E, poisson);
  const double Gc = 50.0;
  const double Kn =
      util::normalContactStiffness(K, K, horizon);
 
  auto model = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
  model["MPI_Strategy"] = mpi_strategy;
  model["Wall_Contact"] = "meshed";
 
  auto output = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force",
                                "Particle_ID"})},
           {"Output_Interval", std::max<size_t>(1, num_steps)},
           {"Debug", 1},
           {"Perform_FE_Out", false},
           {"Compress_Type", "zlib"},
           {"Perform_Out", false},
           {"Test_Output_Interval", num_steps},
           {"Tag_PP", ""},
           {"PVD_Collection", false}});
 
  // Wall (particle index 2) fully fixed; downward IC on both grains.
  auto bc = inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 1},
           {"IC_Sets", 1},
           {"Gravity", util::Point(0, -10, 0).toVec()}});
  bc["Displacement_BC"]["Set_1"] = json{
      {"Particle_List", std::vector<size_t>{2}},
      {"Direction", std::vector<size_t>{1, 2}},
      {"Time_Function",
       {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
      {"Spatial_Function", {{"Type", "constant"}}},
      {"Zero_Displacement", true}};
  bc["IC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "IC"},
           {"Particle_List", {0, 1}},
           {"IC_Type", "Constant_Velocity"},
           {"IC_Vector", {0., -0.02, 0.}}});
 
  auto particle = inp::ParticleDeck::getParticleGeomExampleJson({grain, floor});
  auto mesh_set = [&](const std::filesystem::path &f) {
    return json{{"File", f.string()},
                {"CreateMesh",
                 {{"Flag", true},
                  {"Info", "gmsh_builtin_mesh"},
                  {"Mesh_Size", mesh_size},
                  {"Write_Mesh_File", true}}}};
  };
  json mesh = {{"Sets", 2},
               {"Set_1", mesh_set(mesh_cir)},
               {"Set_2", mesh_set(mesh_wall)}};
 
  json material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
  material["Set_1"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  material["Set_2"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
 
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", 0.95},
               {"Kn", Kn},
               {"Damping_On", false},
               {"Epsilon", 0.95},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 100.0},
               {"K", K}});
  // Soften contact spring so short MPI identity runs stay stable and bit-reproducible.
  contact_base["Kn_Factor"] = 1.0e-4;
  json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
  contact["Set_1_1"] = contact_base;
  contact["Set_1_2"] = contact_base;
  contact["Set_2_2"] = contact_base;
  contact["Damping_Law"] = "off";
  contact["Friction_Law"] = "coulomb_simple";
 
  auto pgen = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pgen["Random_Rotation"] = false;
  pgen["Data"]["N"] = 3;
  pgen["Data"]["0"] = json{{"x", cx0},
                           {"y", cy},
                           {"z", 0.},
                           {"theta", 0.},
                           {"s", 1.},
                           {"geom_id", 0},
                           {"mat_id", 0},
                           {"contact_id", 0}};
  pgen["Data"]["1"] = json{{"x", cx1},
                           {"y", cy},
                           {"z", 0.},
                           {"theta", 0.},
                           {"s", 1.},
                           {"geom_id", 0},
                           {"mat_id", 0},
                           {"contact_id", 0}};
  pgen["Data"]["2"] = json{{"x", 0.},
                           {"y", -0.5 * wall_thick},
                           {"z", 0.},
                           {"theta", 0.},
                           {"s", 1.},
                           {"geom_id", 1},
                           {"mat_id", 1},
                           {"contact_id", 1},
                           {"is_wall", true}};
 
  util::io::print(std::format(
      "mpi_identity_twop_wall: R={:.4g} mesh={:.4g} Rc~{:.4g} gap_pp={:.4g} "
      "cy={:.4g} T={:.3g} N={}\n",
      R, mesh_size, Rc_est, gap_pp, cy, final_time, num_steps));
 
  return json{{"Model", model},
              {"Output", output},
              {"Force_BC", bc["Force_BC"]},
              {"Displacement_BC", bc["Displacement_BC"]},
              {"IC", bc["IC"]},
              {"Particle", particle},
              {"Mesh", mesh},
              {"Material", material},
              {"Contact", contact},
              {"Neighbor",
               inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 5.0},
               {"Search_Interval", 1},
               {"Near_Bd_Nodes_Tol", 0.5}})},
              {"Particle_Generation", pgen}};
}

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

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.

                                               {
  auto model = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
  model["Wall_Contact"] = "analytical_plane";
  if (policy_combo) {
    model["Self_Contact"] = "reference_gap";
    model["Bond_Break"] = "absolute_stretch";
  }
 
  auto output = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path},
           {"Tags", std::vector<std::string>({"Displacement", "Force", "Particle_ID"})},
           {"Output_Interval", std::max<size_t>(1, num_steps / 5)},
           {"Debug", 1},
           {"Perform_FE_Out", false},
           {"Compress_Type", "zlib"},
           {"Perform_Out", false},
           {"Test_Output_Interval", num_steps},
           {"Tag_PP", ""},
           {"PVD_Collection", false}});
 
  auto bc = inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 1}});
  bc["Displacement_BC"]["Set_1"] = json{
      {"Particle_List", std::vector<size_t>{1}},
      {"Direction", std::vector<size_t>{1, 2}},
      {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
      {"Spatial_Function", {{"Type", "constant"}}},
      {"Zero_Displacement", true}};
 
  const double wall_half = 3.0 * R;
  const double wall_thick = mesh_size;
  geom::GeomData grain;
  grain.d_geomName = "circle";
  grain.d_geomParams = {R, 0., 0., 0.};
  geom::GeomData floor;
  floor.d_geomName = "rectangle";
  floor.d_geomParams = {-wall_half, -wall_thick, 0., wall_half, 0., 0.};
 
  auto particle = inp::ParticleDeck::getParticleGeomExampleJson({grain, floor});
 
  auto mesh_set = [&](const std::filesystem::path &f) {
    return json{{"File", f.string()},
                {"CreateMesh",
                 {{"Flag", true},
                  {"Info", "gmsh_builtin_mesh"},
                  {"Mesh_Size", mesh_size},
                  {"Write_Mesh_File", true}}}};
  };
  json mesh = {{"Sets", 2},
               {"Set_1", mesh_set(mesh_cir)},
               {"Set_2", mesh_set(mesh_wall)}};
 
  const double rho = 1200.0;
  const double K = 2.16e7;
  const double G = 1.296e7;
  const double Gc = 50.0;
  auto material = inp::ParticleDeck::getParticleMaterialExampleJson(2);
  material["Set_1"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  material["Set_2"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
 
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", Rc_factor},
               {"Kn", Kn},
               {"Damping_On", false},
               {"Epsilon", 1.0},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 100.0},
               {"K", K}});
  // Soften spring for a short stable settle (formula check uses this Kn too).
  contact_base["Kn"] = Kn;
  contact_base["Kn_Factor"] = 1.0;
  json contact = inp::ParticleDeck::getParticleContactExampleJson(2);
  contact["Set_1_1"] = contact_base;
  contact["Set_1_2"] = contact_base;
  contact["Set_2_2"] = contact_base;
  contact["Damping_Law"] = policy_combo ? "node" : "off";
  contact["Friction_Law"] = policy_combo ? "stick_slip" : "coulomb_simple";
 
  // Bottom of grain near the floor (y=0): slight overlap into Rc after setup.
  const double cy = R + 0.25 * mesh_size;
  auto pgen = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pgen["Random_Rotation"] = false;
  pgen["Data"]["N"] = 2;
  pgen["Data"]["0"] = json{{"x", 0.},
                           {"y", cy},
                           {"z", 0.},
                           {"theta", 0.},
                           {"s", 1.},
                           {"geom_id", 0},
                           {"mat_id", 0},
                           {"contact_id", 0}};
  pgen["Data"]["1"] = json{{"x", 0.},
                           {"y", -0.5 * wall_thick},
                           {"z", 0.},
                           {"theta", 0.},
                           {"s", 1.},
                           {"geom_id", 1},
                           {"mat_id", 1},
                           {"contact_id", 1},
                           {"is_wall", true}};
 
  return json{{"Model", model},
              {"Output", output},
              {"Displacement_BC", bc["Displacement_BC"]},
              {"Particle", particle},
              {"Mesh", mesh},
              {"Material", material},
              {"Contact", contact},
              {"Neighbor",
               inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 5.0},
               {"Search_Interval", 1},
               {"Near_Bd_Nodes_Tol", 0.5}})},
              {"Particle_Generation", pgen}};
}

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::BCDeck::getExampleJson(), inp::ContactPairDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::OutputDeck::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.

                                                                     {
  auto model = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Single_Particle"},
           {"Seed", 0}});
  model["Self_Contact"] = self_contact;
  model["Bond_Break"] = "tension";
 
  // Soft PD; large Gc so only the precrack is fractured.
  const double rho = 1200.0;
  const double K = 2.16e5;
  const double G = 1.296e5;
  const double Gc = 5.0e6;
 
  auto output = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path},
           {"Tags", std::vector<std::string>({"Displacement", "Force", "Damage_Z", "Particle_ID"})},
           {"Output_Interval", num_steps},
           {"Debug", 1},
           {"Perform_FE_Out", false},
           {"Compress_Type", "zlib"},
           {"Perform_Out", false},
           {"Test_Output_Interval", num_steps},
           {"Tag_PP", ""},
           {"PVD_Collection", false}});
 
  const double strip = 0.2 * L;
  json bc = json::object();
  bc["Displacement_BC"] = {
      {"Sets", 2},
      {"Set_1",
       {{"Region",
         {{"Geometry",
           {{"Type", "rectangle"},
            {"Parameters", std::vector<double>{-0.05 * L, -0.05 * L, 0., strip, L + 0.05 * L, 0.}}}}}},
        {"Direction", std::vector<size_t>{1, 2}},
        {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
        {"Spatial_Function", {{"Type", "constant"}}},
        {"Zero_Displacement", true}}},
      {"Set_2",
       {{"Region",
         {{"Geometry",
           {{"Type", "rectangle"},
            {"Parameters",
             std::vector<double>{L - strip, -0.05 * L, 0., L + 0.05 * L, L + 0.05 * L, 0.}}}}}},
        {"Direction", std::vector<size_t>{1}},
        {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
        {"Spatial_Function", {{"Type", "constant"}}},
        {"Zero_Displacement", false}}}};
 
  geom::GeomData rect;
  rect.d_geomName = "rectangle";
  rect.d_geomParams = {0., 0., 0., L, L, 0.};
  auto particle = inp::ParticleDeck::getParticleGeomExampleJson({rect});
 
  json mesh_set = {{"File", mesh_file.string()},
                   {"CreateMesh",
                    {{"Flag", true},
                     {"Info", "gmsh_builtin_mesh"},
                     {"Mesh_Size", mesh_size},
                     {"Write_Mesh_File", true}}}};
  json mesh = {{"Set_1", mesh_set}};
 
  auto material = inp::ParticleDeck::getParticleMaterialExampleJson(1);
  material["Set_1"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
 
  json root = {{"Model", model},
               {"Output", output},
               {"Particle", particle},
               {"Mesh", mesh},
               {"Material", material}};
  root.merge_patch(bc);
  return root;
}

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::MaterialDeck::getExampleJson(), inp::ModelDeck::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.

                                                      {
  // Short tip — point only, not a tall triangle. Narrower tip concentrates load.
  const double W = 0.0008;
  const double H = 0.0012;
 
  // Hollow ellipse (visibly elliptical ring).
  const double a_out = 0.0020;
  const double b_out = 0.0014;
  // Thin ring so a tip-seeded crack can bipartition.
  const double a_in = 0.00180;
  const double b_in = 0.00122;
  const double ell_theta = 0.0;
  const double tip_gap = 0.00008;
  const double tip_y = H;
  const double ell_cy = tip_y + tip_gap + b_out;
 
  const double rho_t = 1200.0;
  const double K_t = 2.16e7;
  const double nu_t = 0.25;
  const double E_t = material::toE(K_t, nu_t);
  const double G_t = material::toGE(E_t, nu_t);
  const double Gc_t = 200.0;
 
  const double rho_e = 1200.0;
  const double K_e = 2.16e7;
  const double nu_e = 0.25;
  const double E_e = material::toE(K_e, nu_e);
  const double G_e = material::toGE(E_e, nu_e);
  // Low Gc so tip damage runs into a diametral through-crack.
  const double Gc_e = 1.0;
 
  const double R_contact_factor = 0.95;
  const double Kn_tt = util::normalContactStiffness(K_t, K_t, horizon);
  const double Kn_ee = util::normalContactStiffness(K_e, K_e, horizon);
  const double Kn_te = util::normalContactStiffness(K_t, K_e, horizon);
 
  const double g = 10.0;
  const double ic_vy = -2.5;
 
  const size_t dt_out_n = std::max<size_t>(1, num_steps / 40);
 
  auto model = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
  model["Bond_Break"] = "tension";
  // No self-contact across broken bonds — otherwise the crack stays glued.
  model["Self_Contact"] = "none";
 
  auto output = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage_Z", "Damage",
                                "Particle_ID", "Fixity"})},
           {"Output_Interval", dt_out_n},
           {"Debug", 1},
           {"Perform_FE_Out", true},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", 1},
           {"Tag_PP", ""},
           {"PVD_Collection", true}});
 
  auto bc = inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 1},
           {"IC_Sets", 1},
           {"Gravity", util::Point(0., -g, 0.).toVec()}});
  bc["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "Displacement_BC"},
           {"Particle_List", {0}},
           {"Direction", {1, 2}},
           {"Zero_Displacement", true}});
 
  bc["IC"]["Set_1"] = {
      {"Particle_List", std::vector<size_t>{1}},
      {"Constant_Velocity", {{"Velocity_Vector", std::vector<double>{0., ic_vy, 0.}}}}};
 
  std::vector<geom::GeomData> pGeom(2);
  pGeom[0].d_geomName = "triangle";
  pGeom[0].d_geomParams = {-0.5 * W, 0., 0., 0.5 * W, 0., 0., 0., H, 0.};
 
  pGeom[1].d_geomName = "ellipse_minus_ellipse";
  pGeom[1].d_geomParams = {a_out, b_out, a_in, b_in, ell_theta, 0., 0., 0.};
 
  for (auto &g : pGeom)
    geom::createGeomObject(g);
 
  const util::Point c_tri = pGeom[0].d_geom_p->center();
 
  json pDeck = json::object();
  pDeck["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeom);
 
  pDeck["Mesh"] = {
      {"Sets", 2},
      {"Set_1",
       {{"File", mesh_tri.string()},
        {"CreateMesh",
         {{"Flag", true},
          {"Info", "gmsh_builtin_mesh"},
          {"Mesh_Size", mesh_size},
          {"Write_Mesh_File", true}}}}},
      {"Set_2",
       {{"File", mesh_ell.string()},
        {"CreateMesh",
         {{"Flag", true},
          {"Info", "gmsh_builtin_mesh"},
          {"Mesh_Size", mesh_size},
          {"Write_Mesh_File", true}}}}}};
 
  auto mat = inp::ParticleDeck::getParticleMaterialExampleJson(2);
  mat["Set_1"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho_t},
           {"K", K_t},
           {"G", G_t},
           {"Gc", Gc_t},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  mat["Set_2"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho_e},
           {"K", K_e},
           {"G", G_e},
           {"Gc", Gc_e},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  pDeck["Material"] = mat;
 
  auto contact = inp::ParticleDeck::getParticleContactExampleJson(2);
  json cpair_wall = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", R_contact_factor},
               {"Kn", Kn_tt},
               {"Damping_On", true},
               {"Epsilon", 0.95},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 100.0},
               {"K", K_t}});
  // Tip contact: enough Beta_n to drive a crack, damping on to limit spray.
  json cpair_tip = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", 0.90},
               {"Kn", Kn_te},
               {"Damping_On", true},
               {"Epsilon", 0.4},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 8.0},
               {"K", util::harmonicMean(K_t, K_e)}});
  json cpair_ell = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", R_contact_factor},
               {"Kn", Kn_ee},
               {"Damping_On", true},
               {"Epsilon", 0.95},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 100.0},
               {"K", K_e}});
  contact["Set_1_1"] = cpair_wall;
  contact["Set_1_1"]["Kn"] = Kn_tt;
  contact["Set_1_1"]["K"] = K_t;
  contact["Set_1_2"] = cpair_tip;
  contact["Set_1_2"]["Kn"] = Kn_te;
  contact["Set_1_2"]["K"] = util::harmonicMean(K_t, K_e);
  contact["Set_2_2"] = cpair_ell;
  contact["Set_2_2"]["Kn"] = Kn_ee;
  contact["Set_2_2"]["K"] = K_e;
  contact["Damping_Law"] = "com_and_node";
  contact["Friction_Law"] = "coulomb_simple";
  pDeck["Contact"] = contact;
 
  pDeck["Neighbor"] = inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 8.0},
               {"Search_Interval", 5},
               {"Near_Bd_Nodes_Tol", 0.5}});
 
  auto gen = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  gen["Random_Rotation"] = false;
  gen["Data"]["N"] = 2;
  gen["Data"]["0"] = {{"x", c_tri.d_x}, {"y", c_tri.d_y}, {"z", 0.0},
                      {"theta", 0.0},   {"s", 1.0},
                      {"geom_id", 0},   {"mat_id", 0},    {"contact_id", 0},
                      {"is_wall", true}};
  gen["Data"]["1"] = {{"x", 0.0}, {"y", ell_cy}, {"z", 0.0},
                      {"theta", 0.0}, {"s", 1.0},
                      {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}};
  pDeck["Particle_Generation"] = gen;
 
  json j = json::object();
  j["Model"] = model;
  j["Output"] = output;
  j["Force_BC"] = bc["Force_BC"];
  j["Displacement_BC"] = bc["Displacement_BC"];
  j["IC"] = bc["IC"];
  j["Particle"] = pDeck["Particle"];
  j["Mesh"] = pDeck["Mesh"];
  j["Material"] = pDeck["Material"];
  j["Contact"] = pDeck["Contact"];
  j["Neighbor"] = pDeck["Neighbor"];
  j["Particle_Generation"] = pDeck["Particle_Generation"];
  return j;
}

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

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.

                                            {
 
  util::io::InputParser input(argc, argv);
 
  const std::vector<double> center = {0.0, 0.0, 0.0};
  const double R = 0.001;
  const double mesh_size = R / 5.0;
  const double horizon = 2.0 * mesh_size;
 
  const double Lin = 0.05;
  const double Win = 0.04;
 
  constexpr double wall_top_inset = 1e-7;
  constexpr double annulus_inset = 1e-7;
  /* Clearance between void floor / plate and particle circum-bounds (~1–2 mesh). */
  constexpr double clearance_mesh = 1.5;
  /* Meshed triangles/hex/drum can extend slightly beyond circumcircle r in VTU. */
  const double geom_pad = 0.25 * mesh_size;
  const double plate_thickness = std::max(3.0 * mesh_size, 2.0 * mesh_size);
 
  const double w_drum2d = R * 0.2;
 
  const double final_time = 0.001;
  size_t num_steps = 400;
  if (input.cmdOptionExists("-numSteps"))
    num_steps = static_cast<size_t>(std::stoul(input.getCmdOption("-numSteps")));
 
  const size_t num_outputs = 10;
  const size_t dt_out_n = std::max<size_t>(1, num_steps / num_outputs);
  const size_t test_dt_out_n = std::max<size_t>(1, dt_out_n / 10);
 
  const double rho_wall = 600.;
  const double poisson_wall = 0.25;
  const double K_wall = 1.e+4;
  const double E_wall = material::toE(K_wall, poisson_wall);
  const double G_wall = material::toGE(E_wall, poisson_wall);
  const double KIc_wall = 5e+6;
  const double Gc_wall = material::toGc(KIc_wall, poisson_wall, E_wall);
 
  const double rho_p = 600.;
  const double poisson_p = poisson_wall;
  const double K_p = 5.e+3;
  const double E_p = material::toE(K_p, poisson_p);
  const double G_p = material::toGE(E_p, poisson_p);
  const double KIc_p = 5e+6;
  const double Gc_p = material::toGc(KIc_p, poisson_p, E_p);
 
  const double R_contact_factor = 0.95;
  const double padding = 1.1 * R_contact_factor * mesh_size;
  const int N_target = 500;
 
  const std::vector<double> in_rect = {center[0] - 0.5 * Lin, center[1] - 0.5 * Win, center[2],
                                       center[0] + 0.5 * Lin, center[1] + 0.5 * Win, center[2]};
 
  std::mt19937 gen(30);
  const double max_y = in_rect[4] - clearance_mesh * mesh_size;
  std::vector<PackedParticle> packed =
      generateParticleLocations(in_rect, max_y, mesh_size, R, N_target, padding, gen);
  if (packed.empty())
    throw std::runtime_error("compression_large_set_inbuilt_2mat2contact: particle pack is empty");
  for (auto &p : packed)
    p.z = 0.0;
 
  const double m = clearance_mesh * mesh_size;
  /* Seat bed: lowest circum-bottom = void_floor + m (void floor = seed inner bottom). */
  {
    const PackedBounds bb0 = packedAxisBounds(packed);
    const double delta_y = (in_rect[1] + m) - bb0.min_y;
    for (auto &p : packed)
      p.y += delta_y;
  }
 
  const PackedBounds bb = packedAxisBounds(packed);
 
  /* Inner void (Gmsh “remove” rectangle): circum-bounds ± m ± geom_pad in x; floor at seed; top = domain. */
  const double void_lox = bb.min_x - m - geom_pad;
  const double void_hix = bb.max_x + m + geom_pad;
  const double void_loy = in_rect[1];
  /* Moving plate: bottom = top of bed + clearance + pad so FE nodes stay under the plate. */
  const double plate_bottom_y = bb.max_y + m + geom_pad;
  double plate_top_y = plate_bottom_y + plate_thickness;
 
  /* Outer solid: wrap void with wall thickness; extend upward slightly past the plate. */
  const double side = std::max(1.5 * horizon, 2.0 * mesh_size);
  const double out_lox = void_lox - side;
  const double out_hix = void_hix + side;
  const double out_loy = void_loy - side;
  const double out_hiy = plate_top_y + wall_top_inset + 2.5 * mesh_size;
 
  std::vector<double> out_rect = {out_lox, out_loy, center[2], out_hix, out_hiy, center[2]};
 
  /* Open-top void: inner cut runs to outer top (inset for Gmsh strict-inside). */
  const double void_hi_y = out_rect[4] - annulus_inset;
  std::vector<double> remove_rect = {void_lox, void_loy, center[2], void_hix, void_hi_y, center[2]};
 
  util::io::print(std::format(
      "[compression_large_set_inbuilt_2mat2contact] bbox(center±r) x∈[{:.6f},{:.6f}] y∈[{:.6f},{:.6f}]; "
      "void x∈[{:.6f},{:.6f}] y_lo {:.6f}; plate y [{:.6f},{:.6f}]; void_hi {:.6f}; outer y∈[{:.6f},{:.6f}]\n",
      bb.min_x, bb.max_x, bb.min_y, bb.max_y, void_lox, void_hix, void_loy, plate_bottom_y, plate_top_y,
      void_hi_y, out_rect[1], out_rect[4]));
 
  std::vector<double> moving_rect = {void_lox, plate_bottom_y, center[2], void_hix, plate_top_y, center[2]};
 
  std::vector<double> fixed_container_params;
  fixed_container_params.insert(fixed_container_params.end(), remove_rect.begin(), remove_rect.end());
  fixed_container_params.insert(fixed_container_params.end(), out_rect.begin(), out_rect.end());
 
  const size_t n_pack = packed.size();
  const size_t n_wall_fixed = n_pack;
  const size_t n_wall_moving = n_pack + 1;
  const size_t n_total = n_pack + 2;
 
  /* Reference geometries at origin; packing uses circum-radius ~ R for all zones. */
  const double ell_a = 0.95 * R;
  const double ell_b = 0.58 * R;
  const double ell_theta = 0.35;
  const double sq_half = R / std::sqrt(2.0);
  /* Axis-aligned rectangle (~same circum extent as R). */
  const double rx = 0.82 * R;
  const double ry = 0.62 * R;
 
  std::vector<geom::GeomData> pGeomVec(10);
  pGeomVec[0].d_geomName = "circle";
  pGeomVec[0].d_geomParams = {R, center[0], center[1], center[2]};
  pGeomVec[1].d_geomName = "triangle";
  pGeomVec[1].d_geomParams = {R, center[0], center[1], center[2]};
  pGeomVec[2].d_geomName = "drum2d";
  pGeomVec[2].d_geomParams = {R, w_drum2d, center[0], center[1], center[2]};
  pGeomVec[3].d_geomName = "hexagon";
  pGeomVec[3].d_geomParams = {R, center[0], center[1], center[2]};
  pGeomVec[4].d_geomName = "ellipse";
  pGeomVec[4].d_geomParams = {ell_a, ell_b, ell_theta, center[0], center[1], center[2]};
  pGeomVec[5].d_geomName = "rectangle";
  pGeomVec[5].d_geomParams = {-rx, -ry, center[2], rx, ry, center[2]};
  pGeomVec[6].d_geomName = "square";
  pGeomVec[6].d_geomParams = {-sq_half, -sq_half, center[2], sq_half, sq_half, center[2]};
  pGeomVec[7].d_geomName = "circle_minus_circle";
  pGeomVec[7].d_geomParams = {center[0], center[1], center[2], R, 0.35 * R};
  pGeomVec[8].d_geomName = "rectangle_minus_rectangle";
  pGeomVec[8].d_geomParams = fixed_container_params;
  pGeomVec[9].d_geomName = "rectangle";
  pGeomVec[9].d_geomParams = moving_rect;
 
  for (auto &g : pGeomVec)
    geom::createGeomObject(g);
 
  /* Annulus composite centroid can pick up tiny numerical z; 2D setup keeps all sites on z = 0. */
  const util::Point cfix = pGeomVec[8].d_geom_p->center();
  const util::Point cmov = pGeomVec[9].d_geom_p->center();
  const util::Point site_wall_fixed(cfix.d_x, cfix.d_y, 0.0);
  const util::Point site_wall_moving(cmov.d_x, cmov.d_y, 0.0);
 
  auto modelDeckJson = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
 
  std::vector<std::string> out_tags = {"Displacement", "Velocity", "Force", "Force_Density", "Damage_Z",
                                         "Damage",      "Nodal_Volume", "Zone_ID", "Particle_ID", "Fixity",
                                         "Force_Fixity", "Contact_Nodes", "No_Fail_Node", "Boundary_Node_Flag"};
  /* Perform_FE_Out must be true so VTU includes element connectivity (appendMesh). If false, only
   * points are written (appendNodes) and ParaView’s default “Surface” view looks empty — use
   * Representation → Points, or enable this flag. Large packs often leave FE out. */
  auto outputDeckJson = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path_for_deck},
           {"Tags", out_tags},
           {"Output_Interval", dt_out_n},
           {"Debug", 2},
           {"Perform_FE_Out", true},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", test_dt_out_n},
           {"Tag_PP", "0"},
           {"PVD_Collection", true}});
 
  auto bcDeckJson = inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 2},
           {"Gravity", util::Point(0, -10, 0).toVec()}});
 
  bcDeckJson["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "Displacement_BC"},
           {"Particle_List", {n_wall_fixed}},
           {"Direction", {1, 2}},
           {"Zero_Displacement", true}});
 
  json set2;
  set2["Particle_List"] = std::vector<size_t>{n_wall_moving};
  set2["Direction"] = std::vector<size_t>{2};
  set2["Time_Function"] = json{{"Type", "linear"}, {"Parameters", std::vector<double>{-0.06}}};
  set2["Spatial_Function"] = json{{"Type", "constant"}};
  bcDeckJson["Displacement_BC"]["Set_2"] = set2;
 
  json pDeckJson = json::object();
 
  pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
 
  json meshRoot = json{{"Sets", 10}};
  const char *mesh_names[] = {"mesh_cir",      "mesh_tri",      "mesh_drum2d",   "mesh_hex",
                              "mesh_ellipse",  "mesh_rect",     "mesh_square",   "mesh_circirc",
                              "mesh_fixed_container", "mesh_moving_container"};
  for (int zi = 0; zi < 10; ++zi) {
    const std::string fname = (inp_dir / (std::string(mesh_names[zi]) + ".msh")).string();
    meshRoot["Set_" + std::to_string(zi + 1)] =
        json{{"File", fname},
             {"CreateMesh",
              json{{"Flag", true}, {"Info", "gmsh_builtin_mesh"}, {"Mesh_Size", mesh_size}, {"Write_Mesh_File", true}}}};
  }
  pDeckJson["Mesh"] = meshRoot;
 
  /* Two material laws: 0 = particle, 1 = wall (fixed + moving). geom_id still selects mesh shape (10 meshes). */
  json matRoot = json{{"Sets", 2}};
  matRoot["Set_1"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho_p},
           {"K", K_p},
           {"G", G_p},
           {"Gc", Gc_p},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  matRoot["Set_2"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho_wall},
           {"K", K_wall},
           {"G", G_wall},
           {"Gc", Gc_wall},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  pDeckJson["Material"] = matRoot;
 
  const double Kn_pp = KnFromBulk(K_p, K_p, horizon);
  const double Kn_pw = KnFromBulk(K_p, K_wall, horizon);
  const double Kn_ww = 0.0;
 
  const double beta_n_eps = 0.95;
  const double friction_coeff = 0.5;
  const bool damping_on = true;
  const bool friction_on = false;
  const double beta_n_factor = 100.;
  const double Kn_factor = 1.;
 
  /* Contact zones: 0 = granular, 1 = wall (fixed + moving). 2×2 matrix → three unique pairs (1-1, 1-2, 2-2). */
  const json j_contact_pp =
      contactPairJson(R_contact_factor, damping_on, friction_on, Kn_pp, beta_n_eps, friction_coeff, Kn_factor,
                      beta_n_factor);
  const json j_contact_pw =
      contactPairJson(R_contact_factor, damping_on, friction_on, Kn_pw, beta_n_eps, friction_coeff, Kn_factor,
                      beta_n_factor);
  const json j_contact_ww =
      contactPairJson(R_contact_factor, damping_on, friction_on, Kn_ww, beta_n_eps, friction_coeff, Kn_factor,
                      beta_n_factor);
 
  json contactRoot = inp::ContactDeck::getExampleJson(json{{"Sets", 2}});
  contactRoot["Set_1_1"] = j_contact_pp;
  contactRoot["Set_1_2"] = j_contact_pw;
  contactRoot["Set_2_2"] = j_contact_ww;
  pDeckJson["Contact"] = contactRoot;
 
  pDeckJson["Neighbor"] = inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 10.0},
               {"Search_Interval", 100},
               {"Near_Bd_Nodes_Tol", 0.5}});
 
  auto pGenJson = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pGenJson["Random_Rotation"] = false;
  pGenJson["Data"]["N"] = n_total;
 
  constexpr size_t k_mat_particle = 0;
  constexpr size_t k_mat_wall = 1;
  constexpr size_t k_contact_grains = 0;
  constexpr size_t k_contact_wall = 1;
 
  for (size_t pi = 0; pi < n_pack; ++pi) {
    const auto &p = packed[pi];
    pGenJson["Data"][std::to_string(pi)] = json{
        {"x", p.x},         {"y", p.y},         {"z", 0.0},         {"theta", p.theta},
        {"s", 1.0},         {"geom_id", static_cast<size_t>(p.zone)},
        {"mat_id", k_mat_particle},
        {"contact_id", k_contact_grains},
    };
  }
 
  pGenJson["Data"][std::to_string(n_pack)] = json{{"x", site_wall_fixed.d_x},
                                                  {"y", site_wall_fixed.d_y},
                                                  {"z", 0.0},
                                                  {"theta", 0.0},
                                                  {"s", 1.0},
                                                  {"geom_id", size_t(8)},
                                                  {"mat_id", k_mat_wall},
                                                  {"contact_id", k_contact_wall}};
 
  pGenJson["Data"][std::to_string(n_pack + 1)] = json{{"x", site_wall_moving.d_x},
                                                       {"y", site_wall_moving.d_y},
                                                       {"z", 0.0},
                                                       {"theta", 0.0},
                                                       {"s", 1.0},
                                                       {"geom_id", size_t(9)},
                                                       {"mat_id", k_mat_wall},
                                                       {"contact_id", k_contact_wall}};
 
  pDeckJson["Particle_Generation"] = pGenJson;
 
  /* Top-level comment only: ignored by inp::Input (only known sections are parsed). */
  return json{{"Comment",
               "compression_large_set_inbuilt"},
              {"Model", modelDeckJson},
              {"Output", outputDeckJson},
              {"Force_BC", bcDeckJson["Force_BC"]},
              {"Displacement_BC", bcDeckJson["Displacement_BC"]},
              {"Particle", pDeckJson["Particle"]},
              {"Mesh", pDeckJson["Mesh"]},
              {"Material", pDeckJson["Material"]},
              {"Contact", pDeckJson["Contact"]},
              {"Neighbor", pDeckJson["Neighbor"]},
              {"Particle_Generation", pDeckJson["Particle_Generation"]}};
}

References util::io::InputParser::cmdOptionExists(), contactPairJson(), geom::createGeomObject(), util::Point::d_x, util::Point::d_y, generateParticleLocations(), util::io::InputParser::getCmdOption(), inp::BCDeck::getExampleJson(), inp::ContactDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::OutputDeck::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(), material::toGE(), and util::Point::toVec().

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.

                                                            {
 
  if (ncols < 1 || nrows < 1)
    throw std::runtime_error("jha2021_comp: ncols and nrows must be >= 1");
  const int ngrains = ncols * nrows;
 
  const std::vector<double> center = {0.0, 0.0, 0.0};
  const double R = 0.001;
  const double mesh_size = R / 5.0;
  const double horizon = 3.0 * mesh_size;
  // Realized hmin on this disk is 0.7 lc. Start slightly outside Rc so contact
  // is not present at t=0; plate and gravity pull pairs into Rc.
  const double h_est = 0.7 * mesh_size;
  const double Rc_est = 0.95 * h_est;
  const double particle_padding = 1.15 * Rc_est;
  const double wpd = particle_padding;
  const double rwp = horizon + wpd;
  constexpr double wall_vy = -0.06;
 
  const double Lin = 2.0 * particle_padding + 2.0 * R +
                     static_cast<double>(ncols - 1) * (2.0 * R + particle_padding);
  const double Win = 2.0 * particle_padding + 2.0 * R +
                     static_cast<double>(nrows - 1) * (2.0 * R + particle_padding);
 
  const double wall_t = rwp - wpd;
  const std::vector<double> mw_rect = {center[0] - wpd, Win, center[2],
                                       Lin + wpd, Win + wall_t, center[2]};
  const std::vector<double> cup_channel = {center[0] - rwp, center[1] - rwp,
                                           Lin + rwp, Win + wall_t, wall_t,
                                           center[2]};
 
  std::vector<geom::GeomData> pGeomVec(3);
  pGeomVec[0].d_geomName = "circle";
  pGeomVec[0].d_geomParams = {R, center[0], center[1], center[2]};
  pGeomVec[1].d_geomName = "open_rect_channel_2d";
  pGeomVec[1].d_geomParams = cup_channel;
  pGeomVec[2].d_geomName = "rectangle";
  pGeomVec[2].d_geomParams = mw_rect;
 
  std::vector<PackedParticle> packed =
      generateCircularGrid(ncols, nrows, 0.0, 0.0, R, particle_padding);
  if (static_cast<int>(packed.size()) != ngrains)
    throw std::runtime_error("jha2021_comp: grid count mismatch");
 
  const double poisson = 0.25;
  const double rho = 1200.0;
  const double K = 2.16e+7;
  const double E = material::toE(K, poisson);
  const double G = material::toGE(E, poisson);
  const double Gc = 50.0;
  const double Kn = util::normalContactStiffness(K, K, horizon);
 
  const double dt = final_time / static_cast<double>(num_steps);
  const double c_wave = std::sqrt(E / rho);
  const double dt_cfl = h_est / c_wave;
  util::io::print(std::format(
      "jha2021_comp: N={} ({}x{}), gap={:.6e} (1.15 Rc_est), wpd={:.6e}, "
      "dt={:.6e} s, h/c={:.6e} s, dt/(h/c)={:.3f}\n",
      ngrains, ncols, nrows, particle_padding, wpd, dt, dt_cfl, dt / dt_cfl));
 
  const size_t num_outputs = 4;
  const size_t dt_out_n = std::max<size_t>(1, num_steps / num_outputs);
  const size_t test_dt_out_n = std::max<size_t>(1, dt_out_n / 10);
 
  for (auto &g : pGeomVec)
    geom::createGeomObject(g);
 
  const util::Point cfix = pGeomVec[1].d_geom_p->center();
  const util::Point cmov = pGeomVec[2].d_geom_p->center();
  const size_t n_pack = packed.size();
  const size_t n_wall_fixed = n_pack;
  const size_t n_wall_moving = n_pack + 1;
  const size_t n_total = n_pack + 2;
 
  auto modelDeckJson = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
  modelDeckJson["MPI_Strategy"] = mpi_strategy;
 
  std::vector<std::string> out_tags = {"Displacement", "Velocity", "Force", "Damage_Z",
                                       "Damage", "Particle_ID", "Contact_Nodes"};
  auto outputDeckJson = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path_for_deck},
           {"Tags", out_tags},
           {"Output_Interval", dt_out_n},
           {"Debug", 1},
           {"Perform_FE_Out", true},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", test_dt_out_n},
           {"Tag_PP", "0"},
           {"PVD_Collection", true}});
 
  auto bcDeckJson =
      inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 2},
           {"Gravity", util::Point(0, -10, 0).toVec()}});
  bcDeckJson["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "Displacement_BC"},
           {"Particle_List", {n_wall_fixed}},
           {"Direction", {1, 2}},
           {"Zero_Displacement", true}});
  json set2;
  set2["Particle_List"] = std::vector<size_t>{n_wall_moving};
  set2["Direction"] = std::vector<size_t>{2};
  set2["Time_Function"] = json{{"Type", "linear"}, {"Parameters", std::vector<double>{wall_vy}}};
  set2["Spatial_Function"] = json{{"Type", "constant"}};
  bcDeckJson["Displacement_BC"]["Set_2"] = set2;
 
  json pDeckJson = json::object();
  pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
 
  auto meshEntry = [&](const std::filesystem::path &f) -> json {
    if (file_mesh)
      return json{{"File", f.string()}};
    return json{{"File", f.string()},
                {"CreateMesh",
                 json{{"Flag", true},
                      {"Info", "gmsh_builtin_mesh"},
                      {"Mesh_Size", mesh_size},
                      {"Write_Mesh_File", write_meshes}}}};
  };
  pDeckJson["Mesh"] = json{{"Sets", 3},
                           {"Set_1", meshEntry(mesh_cir)},
                           {"Set_2", meshEntry(mesh_fixed)},
                           {"Set_3", meshEntry(mesh_moving)}};
 
  json matRoot = json{{"Sets", 2}};
  matRoot["Set_1"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  matRoot["Set_2"] =
      inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  pDeckJson["Material"] = matRoot;
 
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", 0.95},
               {"Kn", Kn},
               {"Damping_On", true},
               {"Epsilon", 0.95},
               {"Friction_On", false},
               {"Friction_Coeff", 0.0},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", 100.0},
               {"K", K}});
  json contactRoot = inp::ParticleDeck::getParticleContactExampleJson(2);
  contactRoot["Set_1_1"] = contact_base;
  contactRoot["Set_1_2"] = contact_base;
  contactRoot["Set_2_2"] = contact_base;
  pDeckJson["Contact"] = contactRoot;
  pDeckJson["Neighbor"] =
      inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 5.0},
               {"Search_Interval", search_interval},
               {"Near_Bd_Nodes_Tol", 0.5}});
 
  auto pGenJson = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pGenJson["Random_Rotation"] = false;
  pGenJson["Data"]["N"] = n_total;
  for (size_t pi = 0; pi < n_pack; ++pi) {
    const auto &p = packed[pi];
    pGenJson["Data"][std::to_string(pi)] = json{{"x", p.x},
                                                 {"y", p.y},
                                                 {"z", 0.0},
                                                 {"theta", 0.0},
                                                 {"s", 1.0},
                                                 {"geom_id", size_t(0)},
                                                 {"mat_id", size_t(0)},
                                                 {"contact_id", size_t(0)}};
  }
  pGenJson["Data"][std::to_string(n_pack)] = json{{"x", cfix.d_x},
                                                  {"y", cfix.d_y},
                                                  {"z", 0.0},
                                                  {"theta", 0.0},
                                                  {"s", 1.0},
                                                  {"geom_id", size_t(1)},
                                                  {"mat_id", size_t(1)},
                                                  {"contact_id", size_t(1)},
                                                  {"is_wall", true}};
  pGenJson["Data"][std::to_string(n_pack + 1)] = json{{"x", cmov.d_x},
                                                       {"y", cmov.d_y},
                                                       {"z", 0.0},
                                                       {"theta", 0.0},
                                                       {"s", 1.0},
                                                       {"geom_id", size_t(2)},
                                                       {"mat_id", size_t(1)},
                                                       {"contact_id", size_t(1)},
                                                       {"is_wall", true}};
  pDeckJson["Particle_Generation"] = pGenJson;
 
  util::io::print(std::format(
      "jha2021_comp: Lin={:.6f}, Win={:.6f}, horizon={:.6f}, plate v_y={}, Wall_Id={}\n",
      Lin, Win, horizon, wall_vy, n_wall_moving));
 
  return json{{"Comment", "jha2021_comp_contact"},
              {"Model", modelDeckJson},
              {"Output", outputDeckJson},
              {"Force_BC", bcDeckJson["Force_BC"]},
              {"Displacement_BC", bcDeckJson["Displacement_BC"]},
              {"Particle", pDeckJson["Particle"]},
              {"Mesh", pDeckJson["Mesh"]},
              {"Material", pDeckJson["Material"]},
              {"Contact", pDeckJson["Contact"]},
              {"Neighbor", pDeckJson["Neighbor"]},
              {"Particle_Generation", pDeckJson["Particle_Generation"]},
              {"Test", json{{"Test_Name", "compressive_test"},
                            {"Compressive_Test",
                             json{{"Wall_Id", n_wall_moving},
                                  {"Wall_Force_Direction", 2}}}}}};
}

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

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 95 of file main.cpp.

                                                            {
 
  const std::vector<double> center = {0.0, 0.0, 0.0};
  const double R1 = 0.001;
  const double R2 = 0.001;
  const double mesh_size =
      mesh_size_in > 0. ? mesh_size_in : std::min(R1, R2) / 5.0;
  const double horizon = horizon_in > 0. ? horizon_in : 3.0 * mesh_size;
  const double particle_dist = 0.001;
 
  const double poisson1 = 0.25;
  const double rho1 = 1200.0;
  const double K1 = 2.16e+7;
  const double E1 = material::toE(K1, poisson1);
  const double G1 = material::toGE(E1, poisson1);
  const double Gc1 = 50.0;
 
  const double poisson2 = 0.25;
  const double rho2 = 1200.0;
  const double K2 = 2.16e+7;
  const double E2 = material::toE(K2, poisson2);
  const double G2 = material::toGE(E2, poisson2);
  const double Gc2 = 50.0;
 
  const double R_contact_factor = 0.95;
  const double Kn_11 = util::normalContactStiffness(K1, K1, horizon);
  const double Kn_22 = util::normalContactStiffness(K2, K2, horizon);
  const double Kn_12 = util::normalContactStiffness(K1, K2, horizon);
  const double friction_coeff = 0.5;
 
  std::vector<double> p1_center = center;
  std::vector<double> p2_center = center;
 
  const double H0_drop = particle_dist; // paper drop height, 1 mm
  // Table 2 energy: remaining gap = horizon, rest of H0 already in IC velocity.
  const double top_gap = two_particle_test ? horizon : particle_dist;
 
  const size_t dt_out_n = num_steps / 10;
  auto modelDeckJson = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
  modelDeckJson["MPI_Strategy"] = mpi_strategy;
 
  auto outputDeckJson = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path_for_deck},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage_Z", "Damage", "Particle_ID"})},
           {"Output_Interval", dt_out_n},
           {"Debug", 2},
           {"Perform_FE_Out", true},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", 1},
           {"Tag_PP", ""},
           {"PVD_Collection", true}});
 
  auto bcDeckJson = inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 1},
           {"IC_Sets", 1},
           {"Gravity", util::Point(0, -10, 0).toVec()}});
 
  if (bottom_patch_bc) {
    // Particle 0 centered at (R1,R1); fix only a bottom strip (not the whole grain).
    geom::GeomData patch;
    patch.d_geomName = "rectangle";
    patch.d_geomParams = {R1 - 1.1 * R1, -0.1 * R1, 0., R1 + 1.1 * R1, 0.35 * R1,
                          0.};
    json j_geom;
    geom::writeGeometry(j_geom, patch);
    bcDeckJson["Displacement_BC"]["Set_1"] = json{
        {"Particle_List", json::array({0})},
        {"Region", {{"Geometry", j_geom}}},
        {"Direction", json::array({1, 2})},
        {"Zero_Displacement", true}};
  } else {
    bcDeckJson["Displacement_BC"]["Set_1"] =
        inp::BCBaseDeck::getExampleJson(
      json{{"Type", "Displacement_BC"},
           {"Particle_List", {0}},
           {"Direction", {1, 2}},
           {"Zero_Displacement", true}});
  }
 
  std::vector<double> ic_vel = {0.0, 0.0, 0.0};
  if (!zero_ic) {
    const double fallen =
        two_particle_test ? (H0_drop - top_gap) : (particle_dist - horizon);
    if (fallen > 0.)
      ic_vel[1] = -std::sqrt(2.0 * std::abs(-10.0) * fallen);
  }
  ic_vel[0] = ic_vx;
  bcDeckJson["IC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "IC"},
           {"Particle_List", {1}},
           {"IC_Type", "Constant_Velocity"},
           {"IC_Vector", ic_vel}});
 
  auto pDeckJson = json({});
 
  std::vector<geom::GeomData> pGeomVec(2);
  pGeomVec[0].d_geomName = "circle";
  pGeomVec[0].d_geomParams = {R1, p1_center[0], p1_center[1], p1_center[2]};
  pGeomVec[1].d_geomName = "circle";
  pGeomVec[1].d_geomParams = {R2, p2_center[0], p2_center[1], p2_center[2]};
 
  pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
 
  const std::string f1 = mesh_file_1.string();
  const std::string f2 = mesh_file_2.string();
  json meshSet1, meshSet2;
  if (file_mesh) {
    meshSet1 = json({{"File", f1}});
    meshSet2 = json({{"File", f2}});
  } else {
    meshSet1 = json({{"File", f1},
                     {"CreateMesh",
                      {{"Flag", true},
                       {"Info", "gmsh_builtin_mesh"},
                       {"Mesh_Size", mesh_size},
                       {"Write_Mesh_File", true}}}});
    meshSet2 = json({{"File", f2},
                     {"CreateMesh",
                      {{"Flag", true},
                       {"Info", "gmsh_builtin_mesh"},
                       {"Mesh_Size", mesh_size},
                       {"Write_Mesh_File", true}}}});
  }
  pDeckJson["Mesh"] = json({{"Sets", 2}, {"Set_1", meshSet1}, {"Set_2", meshSet2}});
 
  auto pMatJson = inp::ParticleDeck::getParticleMaterialExampleJson(2);
  pMatJson["Set_1"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho1},
           {"K", K1},
           {"G", G1},
           {"Gc", Gc1},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  pMatJson["Set_2"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho2},
           {"K", K2},
           {"G", G2},
           {"Gc", Gc2},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  pDeckJson["Material"] = pMatJson;
 
  auto pContactJson = inp::ParticleDeck::getParticleContactExampleJson(2);
  // v0.1.0 circ_damp: Friction_On: false leaves μ = 0 (coeff is not read).
  // Nonzero μ with the flag off still applies a tangential force in PairForce,
  // which gives the falling particle a lateral velocity (Table 2 test 4).
  const double mu =
      (friction_mu >= 0.) ? friction_mu
                          : (two_particle_test ? 0. : friction_coeff);
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", R_contact_factor},
               {"Kn", Kn_11},
               {"Damping_On", damping_on},
               {"Epsilon", eps_n},
               {"Friction_On", friction_on},
               {"Friction_Coeff", mu},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", beta_n_factor},
               {"K", K1}});
 
  pContactJson["Set_1_1"] = contact_base;
  pContactJson["Set_1_1"]["Kn"] = Kn_11;
  pContactJson["Set_1_1"]["K"] = K1;
  pContactJson["Set_1_2"] = contact_base;
  pContactJson["Set_1_2"]["Kn"] = Kn_12;
  pContactJson["Set_1_2"]["K"] = util::harmonicMean(K1, K2);
  pContactJson["Set_2_2"] = contact_base;
  pContactJson["Set_2_2"]["Kn"] = Kn_22;
  pContactJson["Set_2_2"]["K"] = K2;
  pContactJson["Damping_Law"] = damping_law;
  pContactJson["Friction_Law"] = friction_law;
  pDeckJson["Contact"] = pContactJson;
 
  pDeckJson["Neighbor"] = two_particle_test
      ? inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 5.0},
               {"Search_Interval", 1},
               {"Near_Bd_Nodes_Tol", 0.5}})
      : inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 10.0},
               {"Search_Interval", 40},
               {"Near_Bd_Nodes_Tol", 0.5}});
 
  auto pGenJson = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
  pGenJson["Random_Rotation"] = false;
  pGenJson["Data"]["N"] = 2;
  pGenJson["Data"]["0"] = {
      {"x", R1}, {"y", R1}, {"z", 0.0},
      {"theta", 0.0}, {"s", 1.0},
      {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0}
  };
  pGenJson["Data"]["1"] = {
      {"x", R1}, {"y", 2.0 * R1 + R2 + top_gap}, {"z", 0.0},
      {"theta", two_particle_test ? M_PI / 2.0 : M_PI}, {"s", 1.0},
      {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}
  };
  pDeckJson["Particle_Generation"] = pGenJson;
 
  auto j = json({{"Model", modelDeckJson},
               {"Output", outputDeckJson},
               {"Force_BC", bcDeckJson["Force_BC"]},
               {"Displacement_BC", bcDeckJson["Displacement_BC"]},
               {"IC", bcDeckJson["IC"]},
               {"Particle", pDeckJson["Particle"]},
               {"Mesh", pDeckJson["Mesh"]},
               {"Material", pDeckJson["Material"]},
               {"Contact", pDeckJson["Contact"]},
               {"Neighbor", pDeckJson["Neighbor"]},
               {"Particle_Generation", pDeckJson["Particle_Generation"]}});
  if (two_particle_test)
    j["Test"] = json{{"Test_Name", "two_particle"}};
  // bottom_patch_bc keeps PD force on particle 0 (deformable base).
  return j;
}

References geom::GeomData::d_geomName, geom::GeomData::d_geomParams, inp::BCDeck::getExampleJson(), inp::ContactPairDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::PGenDeck::getExampleJson(), inp::PNeighborDeck::getExampleJson(), inp::BCBaseDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), inp::ParticleDeck::getParticleGeomExampleJson(), inp::ParticleDeck::getParticleMaterialExampleJson(), util::harmonicMean(), util::normalContactStiffness(), material::toE(), material::toGE(), util::Point::toVec(), 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 292 of file main.cpp.

                                                                               {
  auto model = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Single_Particle"},
           {"Seed", 0}});
  // Trask §5: broken bond → weight 0 (no force). No self-contact.
  model["Self_Contact"] = "none";
  model["Bond_Break"] = "tension"; // literature PMB: break in tension only
 
  auto output = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_path},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage", "Damage_Z",
                                "Particle_ID", "Fixity"})},
           {"Output_Interval", std::max<size_t>(1, num_steps / 10)},
           {"Debug", 1},
           {"Perform_FE_Out", false},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", num_steps},
           {"Tag_PP", ""},
           {"PVD_Collection", false}});
 
  // Trask §6.2: top left/right of notches u=<0,0>; drive between notches
  // u=<0,-v t>; sides and bottom free (collar bond-break is not implemented).
  const double drive_h = std::max(2.0 * mesh_size, 0.002);
  const std::vector<double> drive_strip{-notch_half, 0.5 * H - drive_h, 0., notch_half,
                                        0.55 * H, 0.};
  json outer = fixedTopOuterBC(W, H, notch_half);
  json bc_disp = {
      {"Sets", 4},
      {"Set_1", outer["Set_1"]},
      {"Set_2", outer["Set_2"]},
      {"Set_3",
       {{"Particle_List", std::vector<size_t>{0}},
        {"Region", {{"Geometry", {{"Type", "rectangle"}, {"Parameters", drive_strip}}}}},
        {"Direction", std::vector<size_t>{1}},
        {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
        {"Spatial_Function", {{"Type", "constant"}}},
        {"Zero_Displacement", true}}},
      {"Set_4",
       {{"Particle_List", std::vector<size_t>{0}},
        {"Region", {{"Geometry", {{"Type", "rectangle"}, {"Parameters", drive_strip}}}}},
        {"Direction", std::vector<size_t>{2}},
        {"Time_Function",
         {{"Type", "linear"}, {"Parameters", std::vector<double>{-v_impact}}}},
        {"Spatial_Function", {{"Type", "constant"}}}}}};
 
  geom::GeomData plate =
      sillingNotchedPlateGeom(W, H, notch_half, notch_w, notch_depth);
  auto particle = inp::ParticleDeck::getParticleGeomExampleJson({plate});
  json mesh = {{"Sets", 1}, {"Set_1", meshSetJson(mesh_plate, mesh_size)}};
  auto material = inp::ParticleDeck::getParticleMaterialExampleJson(1);
  material["Set_1"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PMBBond"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho},
           {"K", K},
           {"G", G},
           {"Gc", Gc},
           {"E", E},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
  // Trask/Silling ω≡1 (ConstInfluence default a0=dim+1 would mis-scale c)
  material["Set_1"]["Influence_Function"] = {
      {"Type", 0}, {"Parameters", std::vector<double>{1.0}}};
 
  return json{{"Model", model},
              {"Output", output},
              {"Displacement_BC", bc_disp},
              {"Particle", particle},
              {"Mesh", mesh},
              {"Material", material}};
}

References fixedTopOuterBC(), inp::MaterialDeck::getExampleJson(), inp::ModelDeck::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 204 of file main.cpp.

                                                                        {
  Probe p;
  if (!dem.d_contact_p || !dem.d_contact_p->d_wallContact ||
      !dem.d_contact_p->d_pairForce)
    throw std::runtime_error("analytical wall: contact path not set up");
 
  const auto &contact = dem.d_particleDeck_p->d_contactDeck.getContact(0, 1);
  const double Kn = contact.d_Kn;
  const double Rc = contact.d_contactR;
 
  for (auto &f : dem.d_f)
    f = util::Point();
  dem.d_contact_p->d_wallContact->apply(dem, dem.d_contact_p->d_pairForce.get(),
                                        /*use_node_damping*/ false);
 
  for (size_t i = 0; i < dem.d_x.size(); ++i) {
    if (dem.d_ptId[i] != 0)
      continue;
    geom::WallContactHit hit;
    if (!wall_geom.wallContactQuery(dem.d_x[i], hit) || !hit.active)
      continue;
    if (!(hit.signed_gap < Rc))
      continue;
 
    ++p.n_contact_nodes;
    p.min_gap = std::min(p.min_gap, hit.signed_gap);
 
    auto scalar = Kn * (hit.signed_gap - Rc) * dem.d_vol[i];
    if (scalar > 0.)
      scalar = 0.;
    const util::Point f_ref = scalar * (-1. * hit.outward_n);
    const double err = (dem.d_f[i] - f_ref).length();
    p.max_force_err = std::max(p.max_force_err, err);
    p.max_force_ref = std::max(p.max_force_ref, f_ref.length());
  }
  if (p.n_contact_nodes == 0)
    p.min_gap = 0.;
  return p;
}

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.

                                                                                    {
  json j;
  j["Contact_Radius_Factor"] = R_contact_factor;
  if (Kn < 1e-20) {
    j["Kn"] = 0.0;
  } else {
    j["Kn"] = Kn;
  }
  j["Damping_On"] = damping_on;
  j["Epsilon"] = beta_n_eps;
  j["Friction_On"] = friction_on;
  j["Friction_Coeff"] = friction_coeff;
  j["Kn_Factor"] = Kn_factor;
  j["Beta_n_Factor"] = damping_on ? beta_n_factor : 0.0;
  return j;
}

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.

                                                                     {
  namespace fs = std::filesystem;
  fs::path n = fs::absolute(dir).lexically_normal();
  std::string s = n.string();
  if (!s.empty() && s.back() != '/' && s.back() != '\\')
    s += fs::path::preferred_separator;
  return s;
}

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.

                                                                          {
  const double p_rect_lox = p.x - p.r;
  const double p_rect_loy = p.y - p.r;
  const double p_rect_hix = p.x + p.r;
  const double p_rect_hiy = p.y + p.r;
  if (p_rect_lox < rect[0] + padding || p_rect_loy < rect[1] + padding || p_rect_hix > rect[3] - padding ||
      p_rect_hiy > rect[4] - padding)
    return true;
 
  for (const auto &q : existing) {
    const double dx = p.x - q.x;
    const double dy = p.y - q.y;
    const double dz = p.z - q.z;
    const double dist = std::sqrt(dx * dx + dy * dy + dz * dz);
    if (dist <= p.r + q.r + padding)
      return true;
  }
  return false;
}

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.

                                                      {
  const size_t n = data.d_u.size();
  std::vector<double> buf(6 * n, 0.);
  for (size_t i = 0; i < n; ++i) {
    if (!ownsNode(data, i))
      continue;
    buf[6 * i + 0] = data.d_u[i].d_x;
    buf[6 * i + 1] = data.d_u[i].d_y;
    buf[6 * i + 2] = data.d_u[i].d_z;
    buf[6 * i + 3] = data.d_v[i].d_x;
    buf[6 * i + 4] = data.d_v[i].d_y;
    buf[6 * i + 5] = data.d_v[i].d_z;
  }
  if (util::parallel::mpiSize() > 1)
    MPI_Allreduce(MPI_IN_PLACE, buf.data(), static_cast<int>(buf.size()),
                  MPI_DOUBLE, MPI_SUM, util::parallel::mpiComm());
  if (util::parallel::mpiRank() != 0)
    return;
 
  std::filesystem::create_directories(out_dir / "nodal");
  {
    std::ofstream xr(out_dir / "nodal" / "x_ref.bin", std::ios::binary);
    const uint32_t nn = static_cast<uint32_t>(n);
    xr.write(reinterpret_cast<const char *>(&nn), sizeof(nn));
    for (size_t i = 0; i < n; ++i) {
      const double p[3] = {data.d_xRef[i].d_x, data.d_xRef[i].d_y,
                           data.d_xRef[i].d_z};
      xr.write(reinterpret_cast<const char *>(p), sizeof(p));
    }
  }
  {
    std::ofstream os(out_dir / "nodal" / "final_uv.bin", std::ios::binary);
    const char magic[4] = {'P', 'D', 'U', 'V'};
    const uint32_t step32 = static_cast<uint32_t>(data.currentStep());
    const uint32_t nn = static_cast<uint32_t>(n);
    os.write(magic, 4);
    os.write(reinterpret_cast<const char *>(&step32), sizeof(step32));
    os.write(reinterpret_cast<const char *>(&nn), sizeof(nn));
    os.write(reinterpret_cast<const char *>(buf.data()),
             static_cast<std::streamsize>(buf.size() * sizeof(double)));
  }
}

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 110 of file main.cpp.

                                                                                      {
  CrackFit fit;
  std::vector<double> xs, ys;
 
  const int nbin = std::max(4, static_cast<int>(std::ceil((band_x - exclude_r) / bin_h)));
  for (int b = 0; b < nbin; ++b) {
    const double r0 = exclude_r + b * bin_h;
    const double r1 = std::min(band_x, r0 + bin_h);
    double best_phi = phi_cut;
    double best_dx = 0., best_dy = 0.;
    bool found = false;
    for (size_t i = 0; i < dem.d_xRef.size(); ++i) {
      if (dem.d_ptId[i] != 0)
        continue;
      const double phi = nodeDamageForFit(dem, i, use_bond_count);
      if (phi < best_phi)
        continue;
      const double dx_ref = outward_sign * (dem.d_xRef[i].d_x - tip_x);
      const double dy_ref = tip_y - dem.d_xRef[i].d_y;
      if (dx_ref < r0 || dx_ref >= r1)
        continue;
      if (dy_ref < -2.0 * bin_h || dy_ref > band_y)
        continue;
      const double dx = dem.d_xRef[i].d_x - tip_x;
      const double dy = tip_y - dem.d_xRef[i].d_y;
      best_phi = phi;
      best_dx = dx;
      best_dy = dy;
      found = true;
    }
    if (found) {
      xs.push_back(best_dx);
      ys.push_back(best_dy);
    }
  }
  fit.n_pts = xs.size();
  if (fit.n_pts < 4)
    return fit;
  double mx = 0., my = 0.;
  for (size_t i = 0; i < fit.n_pts; ++i) {
    mx += xs[i];
    my += ys[i];
  }
  mx /= static_cast<double>(fit.n_pts);
  my /= static_cast<double>(fit.n_pts);
  fit.mean_dx = mx;
  fit.mean_dy = my;
  double cxx = 0., cxy = 0., cyy = 0.;
  for (size_t i = 0; i < fit.n_pts; ++i) {
    const double x = xs[i] - mx;
    const double y = ys[i] - my;
    cxx += x * x;
    cxy += x * y;
    cyy += y * y;
  }
  const double trace = cxx + cyy;
  const double det = cxx * cyy - cxy * cxy;
  const double tmp = std::sqrt(std::max(0., 0.25 * trace * trace - det));
  const double l1 = 0.5 * trace + tmp;
  double vx = cxy;
  double vy = l1 - cxx;
  if (vx * vx + vy * vy < 1.e-30) {
    vx = l1 - cyy;
    vy = cxy;
  }
  const double angle_from_horiz =
      std::atan2(std::abs(vy), std::abs(vx)) * 180. / M_PI;
  fit.angle_to_notch_deg = 90. - angle_from_horiz;
  return fit;
}

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 262 of file main.cpp.

                                                            {
  const double fix_h = 0.08 * H;
  return json{
      {"Set_1",
       {{"Particle_List", std::vector<size_t>{0}},
        {"Region",
         {{"Geometry",
           {{"Type", "rectangle"},
            {"Parameters",
             std::vector<double>{-0.55 * W, 0.5 * H - fix_h, 0., -notch_half - 0.02 * W,
                                 0.55 * H, 0.}}}}}},
        {"Direction", std::vector<size_t>{1, 2}},
        {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
        {"Spatial_Function", {{"Type", "constant"}}},
        {"Zero_Displacement", true}}},
      {"Set_2",
       {{"Particle_List", std::vector<size_t>{0}},
        {"Region",
         {{"Geometry",
           {{"Type", "rectangle"},
            {"Parameters",
             std::vector<double>{notch_half + 0.02 * W, 0.5 * H - fix_h, 0., 0.55 * W, 0.55 * H,
                                 0.}}}}}},
        {"Direction", std::vector<size_t>{1, 2}},
        {"Time_Function", {{"Type", "constant"}, {"Parameters", std::vector<double>{0.}}}},
        {"Spatial_Function", {{"Type", "constant"}}},
        {"Zero_Displacement", true}}}};
}

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.

                                                                                    {
  std::vector<PackedParticle> particles;
  particles.reserve(static_cast<size_t>(ncols * nrows));
  const double step = 2.0 * R + gap;
  const double x0 = x_lo + gap + R;
  const double y0 = y_lo + gap + R;
  for (int j = 0; j < nrows; ++j) {
    for (int i = 0; i < ncols; ++i) {
      PackedParticle p;
      p.x = x0 + static_cast<double>(i) * step;
      p.y = y0 + static_cast<double>(j) * step;
      p.z = 0.0;
      p.r = R;
      p.theta = 0.0;
      particles.push_back(p);
    }
  }
  return particles;
}

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.

                                                                                       {
  std::vector<PackedParticle> particles;
  std::uniform_real_distribution<double> u_r(-0.1 * R, 0.1 * R);
  /* Do not use uniform_int_distribution: libstdc++ vs libc++ map mt19937
   * differently, so geom_id diverges while x,y,theta stay the same. */
 
  const double check_r = R;
  const int rows = static_cast<int>((max_y - in_rect[1]) / (2.0 * check_r));
  const double rect_L = in_rect[3] - in_rect[0];
  const int cols = static_cast<int>(rect_L / (2.0 * check_r));
 
  int counter = 0;
  double x_old = in_rect[0];
  double x_old_right = in_rect[3];
  double y_old = in_rect[1];
  const double cz = 0.0;
 
  std::vector<double> cy_accptd;
  cy_accptd.push_back(y_old);
 
  std::vector<double> row_rads_prev{R, R};
  for (int i = 0; i < rows; ++i) {
    if (i > 0)
      y_old = maxElem(cy_accptd) + maxElem(row_rads_prev);
 
    std::vector<double> row_rads{R};
 
    if (y_old + padding + maxElem(row_rads) >= max_y)
      break;
 
    int num_p_cols = 0;
    int j = 0;
    while (true) {
      if (num_p_cols > cols - 1 || j > 100 * N_target)
        break;
      if (counter >= N_target)
        break;
 
      if (j == 0) {
        x_old = in_rect[0];
        x_old_right = in_rect[3];
      }
 
      const int p_zone = static_cast<int>(gen() % 8);
      const double r0 = R;
      double r = r0 + u_r(gen);
 
      double cx = 0., cy = 0.;
      if (i % 2 == 0) {
        std::uniform_real_distribution<double> rph(-0.1 * r0, 0.05 * r0);
        std::uniform_real_distribution<double> rpv(-0.05 * r0, 0.05 * r0);
        const double cx0 = x_old_right - padding - r;
        cx = cx0 - rph(gen);
        cy = y_old + padding + r + rpv(gen);
      } else {
        std::uniform_real_distribution<double> rph(-0.05 * r0, 0.1 * r0);
        std::uniform_real_distribution<double> rpv(-0.05 * r0, 0.05 * r0);
        const double cx0 = x_old + padding + r;
        cx = cx0 + rph(gen);
        cy = y_old + padding + r + rpv(gen);
      }
 
      PackedParticle trial{p_zone, cx, cy, cz, r, 0.0};
      if (!doesParticleIntersect(trial, particles, in_rect, padding)) {
        std::uniform_real_distribution<double> orient(0.0, 2.0 * M_PI);
        trial.theta = orient(gen);
        particles.push_back(trial);
        row_rads.push_back(trial.r);
        cy_accptd.push_back(cy);
        if (i % 2 == 0)
          x_old_right = cx - trial.r;
        else
          x_old = cx + trial.r;
 
        ++counter;
        num_p_cols++;
      }
      ++j;
    }
    row_rads_prev = std::move(row_rads);
  }
  return particles;
}

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.

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

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.

                                                  {
  for (const std::string &n : geom::getAcceptableGeometries()) {
    if (n == g)
      return true;
  }
  return false;
}

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.

                                                        {
  return util::normalContactStiffness(Ka, Kb, horizon);
}

References util::normalContactStiffness().

Referenced by buildInputJson().

Here is the call graph for this function:

Here is the caller graph for this function:

makeDeck()

json anonymous_namespace{main.cpp}::makeDeck	(	const std::string &	model,
		size_t	dim,
		bool	plane_strain,
		int	influence,
		double	L,
		double	h,
		double	horizon,
		double	E,
		double	nu,
		double	Gc,
		const std::string &	out
	)

Definition at line 32 of file main.cpp.

                                                          {
  const double rho = 1000.;
  json j;
  j["Model"] = inp::ModelDeck::getExampleJson(
      json{{"Dimension", dim}, {"Final_Time", 1.0e-6}, {"Time_Steps", 1},
           {"Particle_Sim_Type", "Single_Particle"}});
  j["Output"] = inp::OutputDeck::getExampleJson(
      json{{"Path", out}, {"Perform_Out", false}, {"Debug", 0}});
  json particle;
  particle["Sets"] = 1;
  particle["Set_1"] =
      json{{"Type", dim == 2 ? "rectangle" : "cuboid"},
           {"Parameters", {-L / 2, -L / 2, dim == 2 ? 0. : -L / 2, L / 2, L / 2,
                           dim == 2 ? 0. : L / 2}}};
  j["Particle"] = particle;
  j["Mesh"] = json{{"Sets", 1},
                   {"Set_1", json{{"File", out + "mesh.vtu"},
                                  {"CreateMesh", json{{"Flag", true},
                                                      {"Info", "uniform"},
                                                      {"Mesh_Size", h}}}}}};
  json mat = inp::ParticleDeck::getParticleMaterialExampleJson(1);
  json m{{"Type", model}, {"Is_Plane_Strain", plane_strain},
         {"Horizon", horizon}, {"Density", rho}, {"E", E}, {"Gc", Gc},
         {"Compute_From_Classical", true},
         {"Influence_Function", json{{"Type", influence}}}};
  if (model == "PDState")
    m["G"] = E / (2. * (1. + nu));
  mat["Set_1"] = inp::MaterialDeck::getExampleJson(m);
  j["Material"] = mat;
  json contact = inp::ParticleDeck::getParticleContactExampleJson(1);
  contact["Set_1_1"] = inp::ContactPairDeck::getExampleJson(
      json{{"Contact_Radius", 0.9 * h}, {"Kn", 1.0}, {"Damping_On", false},
           {"Friction_On", false}, {"Friction_Coeff", 0.}, {"K", E}});
  j["Contact"] = contact;
  j["Particle_Generation"] = json{
      {"Method", "From_File"},
      {"Data", json{{"N", 1},
                    {"0", json{{"x", 0.}, {"y", 0.}, {"z", 0.}, {"theta", 0.},
                               {"s", 1.}, {"geom_id", 0}, {"mat_id", 0},
                               {"contact_id", 0}}}}}};
  return j;
}

References inp::ContactPairDeck::getExampleJson(), inp::MaterialDeck::getExampleJson(), inp::ModelDeck::getExampleJson(), inp::OutputDeck::getExampleJson(), inp::ParticleDeck::getParticleContactExampleJson(), and inp::ParticleDeck::getParticleMaterialExampleJson().

Referenced by main().

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 333 of file main.cpp.

                                                        {
  std::ifstream is(csv);
  if (!is)
    return 0.;
  std::string line;
  std::getline(is, line);
  double max_abs = 0.;
  while (std::getline(is, line)) {
    if (line.empty())
      continue;
    const auto c2 = line.rfind(',');
    if (c2 == std::string::npos)
      continue;
    max_abs = std::max(max_abs, std::abs(std::stod(line.substr(c2 + 1))));
  }
  return max_abs;
}

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.

                                           {
  double m = -std::numeric_limits<double>::infinity();
  for (double x : v)
    m = std::max(m, x);
  return m;
}

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 170 of file main.cpp.

                                                                     {
  CrackMetrics m;
  for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
    const auto &xi = dem.d_xRef[i];
    for (size_t k = 0; k < dem.d_neighPd[i].size(); ++k) {
      if (!dem.d_fracture_p->getBondState(i, k))
        continue;
      const size_t j = dem.d_neighPd[i][k];
      const auto &xj = dem.d_xRef[j];
      if ((xi.d_x - x_crack) * (xj.d_x - x_crack) >= 0.)
        continue;
      const double gap = (dem.d_x[j] - dem.d_x[i]).length();
      m.min_gap = std::min(m.min_gap, gap);
      m.max_gap = std::max(m.max_gap, gap);
      ++m.n_pairs;
    }
  }
  if (m.n_pairs == 0)
    m.min_gap = 0.;
 
  for (size_t i = 0; i < dem.d_xRef.size(); ++i) {
    const double x = dem.d_xRef[i].d_x;
    const double fx = dem.d_f[i].d_x * dem.d_vol[i];
    if (x < x_crack && x > x_crack - band)
      m.force_x_left += fx;
    if (x >= x_crack && x < x_crack + band)
      m.force_x_right += fx;
  }
  return m;
}

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 183 of file main.cpp.

                                                               {
  return json{{"File", f.string()},
              {"CreateMesh",
               {{"Flag", true},
                {"Info", "gmsh_builtin_mesh"},
                {"Mesh_Size", mesh_size},
                {"Write_Mesh_File", true}}}};
}

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 244 of file main.cpp.

                                                                     {
  double ming = 1.e300;
  size_t n = 0;
  for (size_t i = 0; i < dem.d_x.size(); ++i) {
    if (dem.d_ptId[i] != 0)
      continue;
    geom::WallContactHit hit;
    if (!wall_geom.wallContactQuery(dem.d_x[i], hit) || !hit.active)
      continue;
    ming = std::min(ming, hit.signed_gap);
    ++n;
  }
  return n ? ming : 0.;
}

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.

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

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 91 of file main.cpp.

                                                                                {
  if (use_bond_count) {
    if (!dem.d_phiBond.empty() && i < dem.d_phiBond.size())
      return static_cast<double>(dem.d_phiBond[i]);
    return nodePhi(dem, i);
  }
  if (!dem.d_phi.empty() && i < dem.d_phi.size())
    return static_cast<double>(dem.d_phi[i]);
  return nodePhi(dem, i);
}

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 79 of file main.cpp.

                                                  {
  const auto &nb = dem.d_neighPd[i];
  if (nb.empty())
    return 0.;
  size_t n_br = 0;
  for (size_t k = 0; k < nb.size(); ++k) {
    if (dem.d_fracture_p->getBondState(i, k))
      ++n_br;
  }
  return static_cast<double>(n_br) / static_cast<double>(nb.size());
}

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 207 of file main.cpp.

                                                                          {
  const double hw = 0.5 * notch_w;
  const double y_tip = 0.5 * H - notch_depth;
  const double y_hi = 0.5 * H + 1.0e-9;
  return {{-notch_half - hw, y_tip, z_lo, -notch_half + hw, y_hi, z_hi},
          {notch_half - hw, y_tip, z_lo, notch_half + hw, y_hi, z_hi}};
}

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.

                                                   {
  const int rank = util::parallel::mpiRank();
  if (data.d_pdDofMpi) {
    if (i >= data.d_pdNodePartition.size())
      return false;
    return static_cast<int>(data.d_pdNodePartition[i]) == rank;
  }
  for (const auto *p : data.d_particlesListTypeAll) {
    const size_t i0 = p->d_globStart;
    const size_t i1 = i0 + p->getNumNodes();
    if (i < i0 || i >= i1)
      continue;
    if (p->isWall())
      return rank == 0;
    return particle::isLocallyOwned(*p);
  }
  return rank == 0;
}

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.

                                                                       {
  PackedBounds b;
  b.min_x = std::numeric_limits<double>::infinity();
  b.max_x = -std::numeric_limits<double>::infinity();
  b.min_y = std::numeric_limits<double>::infinity();
  b.max_y = -std::numeric_limits<double>::infinity();
  for (const auto &p : packed) {
    b.min_x = std::min(b.min_x, p.x - p.r);
    b.max_x = std::max(b.max_x, p.x + p.r);
    b.min_y = std::min(b.min_y, p.y - p.r);
    b.max_y = std::max(b.max_y, p.y + p.r);
  }
  return b;
}

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,
		bool	deck_only = false
	)

Definition at line 207 of file main.cpp.

                                         {
  namespace fs = std::filesystem;
  const fs::path out_dir = run_dir / "out";
  const fs::path inp_dir = run_dir / "inp";
  fs::create_directories(out_dir);
  fs::create_directories(inp_dir);
 
  const fs::path mesh = inp_dir / "mesh_rect.msh";
  // Time fields unused for prescribed kinematics; keep Output_Interval valid.
  auto input_json =
      buildInputJson(directoryPathWithTrailingSep(out_dir), mesh, L, mesh_size, horizon, 1.0e-6,
                     2, self_contact);
  {
    std::ofstream os(inp_dir / "input.json");
    os << input_json.dump(2);
  }
 
  // -deckOnly stops after the deck is written. The deck comparison reads the
  // deck and does not need this driver to run the simulation.
  if (deck_only) {
    util::io::print("deck written; -deckOnly, not running\n");
    return LawResult{};
  }
 
  auto deck = std::make_shared<inp::Input>(input_json);
  PeriDEMModel dem(deck);
  dem.init();
 
  const double x_crack = 0.5 * L;
  const size_t n_broken = applyPrecrack(dem, x_crack);
  if (n_broken == 0)
    throw std::runtime_error(self_contact + ": precrack broke no bonds");
 
  const double Rc = dem.d_maxContactR > 0. ? dem.d_maxContactR : 0.95 * mesh_size;
  const double h = dem.d_hMax > 0. ? dem.d_hMax : mesh_size;
  util::io::print(std::format("{}: nodes={}, broken_bonds={}, h={}, Rc={}, horizon={}\n",
                              self_contact, dem.d_x.size(), n_broken, h, Rc, horizon));
 
  // Open: separate faces so current gap >> Rc (no self-contact yet).
  const double u_open = 2.0 * Rc;
  applyRigidOpenClose(dem, x_crack, u_open);
  dem.computeForces();
  const auto open_m = measure(dem, x_crack, 2.0 * horizon);
  if (!(open_m.min_gap > 1.5 * Rc))
    throw std::runtime_error(std::format("{}: open did not separate (min_gap={})", self_contact,
                                         open_m.min_gap));
 
  // Close into the window where reference_gap (natural R = r0) and
  // broken_bond_kn (natural R = Rc) disagree: Rc < R < typical r0.
  // Mid-horizon neighbor distances are ~O(horizon) >> Rc = 0.95 h.
  const double u_close = -0.5 * (Rc + 0.5 * horizon);
  applyRigidOpenClose(dem, x_crack, u_close);
  dem.computeForces();
  const auto close_m = measure(dem, x_crack, 2.0 * horizon);
  const double force_amp =
      0.5 * (std::abs(close_m.force_x_left) + std::abs(close_m.force_x_right));
 
  util::io::print(std::format("{}: open min_gap={}, close min_gap={}, force_amp={}\n", self_contact,
                              open_m.min_gap, close_m.min_gap, force_amp));
 
  dem.close();
 
  LawResult r;
  r.open_min_gap = open_m.min_gap;
  r.close_min_gap = close_m.min_gap;
  r.close_force_amp = force_amp;
  return r;
}

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:

sillingNotchedPlateGeom()

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

Definition at line 219 of file main.cpp.

                                                                           {
  const double x0 = -0.5 * W, y0 = -0.5 * H, x1 = 0.5 * W, y1 = 0.5 * H;
  const double hw = 0.5 * notch_w;
  const double y_tip = y1 - notch_depth;
  geom::GeomData plate;
  plate.d_geomName = "complex";
  plate.d_geomComplexInfo = {
      std::vector<std::string>{"rectangle", "rectangle", "rectangle"},
      std::vector<std::string>{"plus", "minus", "minus"}};
  plate.d_geomParams = {
      // outer plate
      x0, y0, 0., x1, y1, 0.,
      // left notch cutout (opens at top)
      -notch_half - hw, y_tip, 0., -notch_half + hw, y1 + 1.0e-6, 0.,
      // right notch cutout
      notch_half - hw, y_tip, 0., notch_half + hw, y1 + 1.0e-6, 0.};
  return plate;
}

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 195 of file main.cpp.

                                                                    {}) {
  json cm = {{"Flag", true},
             {"Info", "uniform"},
             {"Mesh_Size", mesh_size},
             {"Write_Mesh_File", true}};
  if (!voids.empty())
    cm["Void_Regions"] = voids;
  return json{{"File", f.string()}, {"CreateMesh", cm}};
}

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 325 of file main.cpp.

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

Referenced by main().

Here is the caller graph for this function: