main.cpp File Reference

#include "fracture/prenotch.h"
#include "inp/deckIncludes.h"
#include "material/materialUtil.h"
#include "periDEMModel.h"
#include "time_int/integrator.h"
#include "util/function.h"
#include "util/io.h"
#include "util/parallelUtil.h"
#include "util/point.h"
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <filesystem>
#include <format>
#include <fstream>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>

Include dependency graph for main.cpp:

Go to the source code of this file.

Data Structures
struct	anonymous_namespace{main.cpp}::CrackFit

Namespaces
namespace	anonymous_namespace{main.cpp}

Functions
size_t	anonymous_namespace{main.cpp}::applyNotches (PeriDEMModel &dem, double x_left, double x_right, double y_lo, double y_hi)
size_t	anonymous_namespace{main.cpp}::applyNotchSlots (PeriDEMModel &dem, double notch_half, double notch_w, double y_lo, double y_hi)
size_t	anonymous_namespace{main.cpp}::applyNotchMidplanes (PeriDEMModel &dem, double notch_half, double y_tip, double y_top)
std::string	anonymous_namespace{main.cpp}::directoryPathWithTrailingSep (const std::filesystem::path &dir)
double	anonymous_namespace{main.cpp}::nodePhi (const PeriDEMModel &dem, size_t i)
double	anonymous_namespace{main.cpp}::nodeDamageForFit (const PeriDEMModel &dem, size_t i, bool use_bond_count)
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)
json	anonymous_namespace{main.cpp}::meshSetJson (const std::filesystem::path &f, double mesh_size)
json	anonymous_namespace{main.cpp}::uniformMeshSetJson (const std::filesystem::path &f, double mesh_size, const std::vector< std::vector< double > > &voids={})
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)
geom::GeomData	anonymous_namespace{main.cpp}::sillingNotchedPlateGeom (double W, double H, double notch_half, double notch_w, double notch_depth)
geom::GeomData	anonymous_namespace{main.cpp}::bhatVNotchedPlateGeom (double W, double H, double notch_half, double notch_w, double notch_depth)
json	anonymous_namespace{main.cpp}::fixedTopOuterBC (double W, double H, double notch_half)
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)
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.)
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)
int	main (int argc, char *argv[])

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 717 of file main.cpp.

                                 {
  util::parallel::initMpi(argc, argv);
  util::io::InputParser input(argc, argv);
 
  unsigned n_threads = 1;
  if (input.cmdOptionExists("-nThreads"))
    n_threads = static_cast<unsigned>(std::stoi(input.getCmdOption("-nThreads")));
  util::parallel::initNThreads(n_threads);
 
  const bool quick = input.cmdOptionExists("-quick");
  const bool trask = input.cmdOptionExists("-traskDisp");
  const bool bhat = input.cmdOptionExists("-bhatKW");
  // Silling's actual case is 3D: a 9 mm plate on a 200 x 100 x 9 grid.
  const bool dim3 = input.cmdOptionExists("-dim3");
 
  namespace fs = std::filesystem;
  // Always write under the binary cwd (build/linux/.../notched_impact_inbuilt/) unless overridden.
  std::string run_tag =
      quick ? "quick" : (trask ? "trask" : (bhat ? "bhat" : (dim3 ? "silling3d" : "lit")));
  // New campaign outputs go under runs_new/ (do not mix with archived runs/).
  fs::path base = fs::current_path() / "runs_new" / run_tag;
  if (input.cmdOptionExists("-outputDir"))
    base = input.getCmdOption("-outputDir");
  const fs::path out_dir = base / "out";
  const fs::path inp_dir = base / "inp";
  fs::create_directories(out_dir);
  fs::create_directories(inp_dir);
  std::cout << std::format("notched_impact: outputDir={}\n", base.string());
 
  // Fig. 2 Silling 2003: 200×100 mm plate, two 1.5 mm notches depth 50 mm, tip-to-tip 50 mm.
  double W = 0.200;
  double H = 0.100;
  double notch_depth = 0.050;
  double notch_half = 0.025;   // notch centerline ±25 mm
  double notch_w = 0.0015;     // Silling 1.5 mm open gap
  double plate_thickness = 0.009; // Silling Fig. 2: 9 mm (only used with -dim3)
  // Impactor: Silling Fig. 2/4 give a cylinder of 1.57 kg striking edge-on, as
  // wide as the 50 mm ligament between the notches. At ρ = 8000 that mass fixes
  // its length: 1.57/(8000·π·0.025²) = 0.100 m. Bhat Fig. 4(b) draws exactly
  // this 50 × 100 mm section, giving the same width and length.
  double Iw = 0.050;
  double Ih = 0.100;
  // Silling 2003 EMU grid 200×100×9 on the mm plate → h ≈ 1 mm; Trask KW: δ = 3h.
  double mesh_size = 0.001;
  double horizon = 3.0 * mesh_size; // 3 mm
  double Rc_factor = 0.95;
  // The impactor must start outside the contact radius, which the solver sets to
  // Rc_factor·h_min: any closer and the first step applies a saturated repulsive
  // force that fires the struck nodes off at kilometres per second. Nothing acts
  // on the impactor while it closes this gap, so its only effect is to delay
  // contact by gap/v, which the crack-speed fit measures from damage arrival.
  double gap = 1.5 * Rc_factor * mesh_size;
 
  // Table 2 M1 ≡ Silling maraging steel (E=191 GPa, ν≈0.3).
  // Gc from KIc≈90 MPa√m: plane-stress KIc²/E = 42408 J/m² (Trask/Silling handbook).
  double rho = 8000.0;
  double E = 191.0e9;
  double Kbulk = 159.2e9;
  double nu = 0.5 * (1.0 - E / (3.0 * Kbulk));
  double G = E / (2.0 * (1.0 + nu));
  double Gc = 42408.0;
  double v_impact = 32.0;
  double dt = 2.5e-9; // Bhat §6.1
  // Silling's cracks run all the way to the free edges (75 mm) at ~900 m/s, so
  // ~85 µs of propagation on top of the free flight and the initiation delay.
  double final_time = 1.7e-4; // 170 µs
  if (input.cmdOptionExists("-Gc"))
    Gc = std::stod(input.getCmdOption("-Gc"));
  if (input.cmdOptionExists("-finalTime"))
    final_time = std::stod(input.getCmdOption("-finalTime"));
  if (input.cmdOptionExists("-vImpact"))
    v_impact = std::stod(input.getCmdOption("-vImpact"));
  if (input.cmdOptionExists("-meshSize")) {
    mesh_size = std::stod(input.getCmdOption("-meshSize"));
    horizon = 3.0 * mesh_size;
  }
  if (input.cmdOptionExists("-horizon"))
    horizon = std::stod(input.getCmdOption("-horizon"));
  if (input.cmdOptionExists("-horizonFactor"))
    horizon = std::stod(input.getCmdOption("-horizonFactor")) * mesh_size;
 
  if (quick) {
    W = 0.040;
    H = 0.020;
    notch_depth = 0.5 * H;
    notch_half = 0.125 * W;
    mesh_size = H / 16.0;
    notch_w = std::max(mesh_size, 0.1 * notch_half);
    Iw = 0.008;
    Ih = 0.004;
    gap = 1.5 * Rc_factor * mesh_size;
    horizon = 3.0 * mesh_size;
    rho = 1200.0;
    E = 1.23e9;
    Kbulk = 2.0e9;
    nu = 0.5 * (1.0 - E / (3.0 * Kbulk));
    G = E / (2.0 * (1.0 + nu));
    Gc = 424.0;
    dt = 1.0e-8;
    final_time = 1.0e-4;
  }
 
  // Contact stiffness, Bhat Eq. (4.1). The solver's contact force density is
  // Kn·V_j·overlap, so Kn carries one power of δ per spatial dimension of the
  // node weight: 18k/(πδ⁴) with 2D weights h²·1 m, 18k/(πδ⁵) with 3D weights h³.
  // Using the 3D form in 2D makes contact 1/δ ≈ 333× too stiff.
  const double Kn =
      util::normalContactStiffness(Kbulk, Kbulk, horizon, dim3 ? 5 : 4);
  const size_t num_steps = static_cast<size_t>(std::llround(final_time / dt));
 
  json input_json;
  if (trask) {
    input_json = buildTraskInputJson(directoryPathWithTrailingSep(out_dir),
                                     inp_dir / "mesh_plate.msh", W, H, notch_half, notch_w,
                                     notch_depth, mesh_size, horizon, rho, E, Kbulk, G, Gc,
                                     v_impact, final_time, num_steps);
  } else if (bhat) {
    input_json = buildBhatInputJson(
        directoryPathWithTrailingSep(out_dir), inp_dir / "mesh_plate.msh",
        inp_dir / "mesh_impactor.msh", W, H, notch_half, notch_w, notch_depth, Iw, Ih, gap,
        mesh_size, horizon, Rc_factor, rho, E, Kbulk, G, Gc, v_impact, final_time, num_steps);
  } else {
    input_json = buildImpactInputJson(
        directoryPathWithTrailingSep(out_dir), inp_dir / "mesh_plate.msh",
        inp_dir / "mesh_impactor.msh", W, H, notch_half, notch_w, notch_depth, Iw, Ih, gap,
        mesh_size, horizon, Rc_factor, Kn, rho, E, Kbulk, G, Gc, v_impact, final_time,
        num_steps, dim3 ? plate_thickness : 0.);
  }
  // Optional override (e.g. Bhat paper default is reference_gap; -selfContact none
  // isolates prenotch+void setups from restorative self-contact across the notch).
  if (input.cmdOptionExists("-selfContact")) {
    const auto sc = input.getCmdOption("-selfContact");
    input_json["Model"]["Self_Contact"] = sc;
    std::cout << std::format("notched_impact: Self_Contact override -> {}\n", sc);
  }
  if (input.cmdOptionExists("-bondBreak")) {
    const auto bb = input.getCmdOption("-bondBreak");
    input_json["Model"]["Bond_Break"] = bb;
    std::cout << std::format("notched_impact: Bond_Break override -> {}\n", bb);
  }
  // The deck fixes Output_Interval at num_steps/10, which is right for a
  // check run and far too coarse for a frame sequence. Overriding it costs
  // nothing and leaves the physics untouched.
  if (input.cmdOptionExists("-outInterval")) {
    const size_t oi =
        static_cast<size_t>(std::stoull(input.getCmdOption("-outInterval")));
    input_json["Output"]["Output_Interval"] = std::max<size_t>(1, oi);
    std::cout << std::format("notched_impact: Output_Interval override -> {}\n", oi);
  }
  if (input.cmdOptionExists("-knScale")) {
    const double s = std::stod(input.getCmdOption("-knScale"));
    for (auto &kv : input_json["Contact"].items()) {
      if (!kv.value().is_object() || !kv.value().contains("Kn"))
        continue;
      kv.value()["Kn"] = kv.value()["Kn"].get<double>() * s;
    }
    std::cout << std::format("notched_impact: Contact Kn scaled by {}\n", s);
  }
  {
    std::ofstream os(inp_dir / "input.json");
    os << input_json.dump(2);
  }
 
  // -deckOnly stops after the deck is written. The Python comparison reads
  // the deck and does not need this driver to run the simulation.
  if (input.cmdOptionExists("-deckOnly")) {
    util::io::print("deck written; -deckOnly, not running\n");
    util::parallel::finalizeMpi();
    return 0;
  }
 
  auto deck = std::make_shared<inp::Input>(input_json);
  PeriDEMModel dem(deck);
  dem.init();
 
  const double y_top = 0.5 * H;
  const double y_tip = y_top - notch_depth;
  // Prenotch bonds: Bhat §5 removes bonds whose segment leaves the nonconvex
  // V-domain (midplane cuts through each V tip); Silling/Trask: slot spans.
  const size_t n_pre =
      bhat ? applyNotchMidplanes(dem, notch_half, y_tip, y_top + 0.01 * H)
           : applyNotchSlots(dem, notch_half, notch_w, y_tip, y_top + 0.01 * H);
  if (n_pre < 10) {
    std::cerr << "notched_impact: expected prenotch bonds, got " << n_pre << "\n";
    return 1;
  }
 
  const char *mode_str =
      quick ? "quick"
            : (trask ? "trask_disp"
                     : (bhat ? "bhat_kw" : (dim3 ? "silling_3d" : "silling_impact")));
  std::cout << std::format(
      "notched_impact: mode={} nodes={} prenotch={} h={:.4e} ε={:.4e} notch_w={:.4e} "
      "Gc={:.0f} E={:.3e} v={:.1f} Nt={} T={:.3e}\n",
      mode_str, dem.d_x.size(), n_pre, mesh_size, horizon, notch_w, Gc, E, v_impact, num_steps,
      final_time);
 
  // Drive the loop ourselves (same sequence as time_int::Integrator::integrate)
  // so we can record when each node first becomes damaged. That arrival-time
  // field gives the crack speed, which Silling reports as ~900 m/s.
  std::vector<float> arrival(dem.d_x.size(), -1.f);
  const size_t sample_every = std::max<size_t>(1, num_steps / 400);
  auto sampleArrivals = [&]() {
    for (size_t i = 0; i < dem.d_xRef.size(); ++i) {
      if (arrival[i] >= 0.f || dem.d_ptId[i] != 0)
        continue;
      if (nodePhi(dem, i) >= 0.30)
        arrival[i] = static_cast<float>(dem.d_time);
    }
  };
 
  dem.applyInitialCondition();
  if (dem.performOutput())
    dem.output();
  dem.setCurrentDt(dem.timeStep());
  dem.applyDisplacementBC();
  dem.computeForces();
  time_int::applyRigidBodyConstraint(dem);
  while (dem.currentStep() < dem.numTimeSteps()) {
    dem.integrateStep();
    if (dem.shouldOutput())
      dem.output();
    if (dem.currentStep() % sample_every == 0)
      sampleArrivals();
    dem.checkStop();
  }
  sampleArrivals();
 
  double vmax = 0., phi_max = 0.;
  size_t n_phi30 = 0;
  for (size_t i = 0; i < dem.d_v.size(); ++i) {
    vmax = std::max(vmax, dem.d_v[i].length());
    if (dem.d_ptId[i] != 0)
      continue;
    const double phi = nodePhi(dem, i);
    phi_max = std::max(phi_max, phi);
    if (phi >= 0.30)
      ++n_phi30;
  }
  const float zmax =
      dem.d_Z.empty() ? 0.f : *std::max_element(dem.d_Z.begin(), dem.d_Z.end());
  std::cout << std::format(
      "notched_impact diag: max|v|={:.3e} maxφ={:.3f} maxZ={:.3f} n(φ≥0.3)={}\n", vmax,
      phi_max, zmax, n_phi30);
 
  if (!(vmax < 8.0e3) || !std::isfinite(vmax)) {
    std::cerr << "notched_impact: blow-up, max |v| = " << vmax << "\n";
    return 1;
  }
 
  size_t n_extra = 0;
  for (size_t i = 0; i < dem.d_neighPd.size(); ++i) {
    if (dem.d_ptId[i] != 0)
      continue;
    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];
      if (dem.d_ptId[j] != 0)
        continue;
      const auto &xi = dem.d_xRef[i];
      const auto &xj = dem.d_xRef[j];
      const bool on_notch =
          geometry::segmentCrossesVerticalLine(xi.d_x, xi.d_y, xj.d_x, xj.d_y,
                                               -notch_half, y_tip,
                                               y_top + 0.01 * H) ||
          geometry::segmentCrossesVerticalLine(xi.d_x, xi.d_y, xj.d_x, xj.d_y,
                                               notch_half, y_tip,
                                               y_top + 0.01 * H);
      if (!on_notch)
        ++n_extra;
    }
  }
  if (n_extra < 20) {
    std::cerr << "notched_impact: expected new broken bonds, got " << n_extra
              << " (prenotch=" << n_pre << ")\n";
    return 1;
  }
 
  // Dump the damage field so the crack-path metric can be re-fitted offline
  // without repeating the simulation.
  {
    std::ofstream csv(base / "damage.csv");
    csv << "x_ref,y_ref,z_ref,x_cur,y_cur,z_cur,phi,phi_bond,arrival\n";
    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=*/false);
      const double phi_b = nodeDamageForFit(dem, i, /*use_bond_count=*/true);
      if (std::max(phi, phi_b) < 0.05)
        continue;
      csv << std::format("{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.4f},{:.4f},{:.6e}\n",
                         dem.d_xRef[i].d_x, dem.d_xRef[i].d_y, dem.d_xRef[i].d_z,
                         dem.d_x[i].d_x, dem.d_x[i].d_y, dem.d_x[i].d_z, phi, phi_b,
                         double(arrival[i]));
    }
  }
 
  const double band_x = 0.45 * W;
  const double band_y = 0.55 * H;
  const double phi_cut = quick ? 0.20 : 0.35;
  const double exclude_r = 2.0 * mesh_size;
  const bool use_bond_dmg = bhat;
  auto left = fitCrackFromPhi(dem, -notch_half, y_tip, -1.0, phi_cut, band_x, band_y, exclude_r,
                              mesh_size, use_bond_dmg);
  auto right = fitCrackFromPhi(dem, notch_half, y_tip, +1.0, phi_cut, band_x, band_y, exclude_r,
                               mesh_size, use_bond_dmg);
  if (left.n_pts < 4 || right.n_pts < 4) {
    std::cerr << "notched_impact: insufficient φ-ridge for path fit (L=" << left.n_pts
              << " R=" << right.n_pts << ")\n";
    return 1;
  }
  if (!(left.mean_dx < 0. && left.mean_dy > 0. && right.mean_dx > 0. && right.mean_dy > 0.)) {
    std::cerr << std::format(
        "notched_impact: tip damage not outward/down "
        "(L dx={:.3e} dy={:.3e} R dx={:.3e} dy={:.3e})\n",
        left.mean_dx, left.mean_dy, right.mean_dx, right.mean_dy);
    return 1;
  }
 
  const double ang = 0.5 * (left.angle_to_notch_deg + right.angle_to_notch_deg);
  std::cout << std::format(
      "notched_impact PASS: prenotch={} new_broken={} ang_to_notch={:.1f}° "
      "(L={:.1f}/{} R={:.1f}/{}) max|v|={:.3e}\n",
      n_pre, n_extra, ang, left.angle_to_notch_deg, left.n_pts, right.angle_to_notch_deg,
      right.n_pts, vmax);
 
  dem.close();
  return 0;
}

References PeriDEMModel::applyDisplacementBC(), PeriDEMModel::applyInitialCondition(), anonymous_namespace{main.cpp}::applyNotchMidplanes(), anonymous_namespace{main.cpp}::applyNotchSlots(), time_int::applyRigidBodyConstraint(), anonymous_namespace{main.cpp}::buildBhatInputJson(), anonymous_namespace{main.cpp}::buildImpactInputJson(), anonymous_namespace{main.cpp}::buildTraskInputJson(), PeriDEMModel::checkStop(), PeriDEMModel::close(), util::io::InputParser::cmdOptionExists(), PeriDEMModel::computeForces(), data::ModelData::currentStep(), data::ModelData::d_fracture_p, data::ModelData::d_neighPd, data::ModelData::d_ptId, data::ModelData::d_time, data::ModelData::d_v, data::ModelData::d_x, data::ModelData::d_xRef, data::ModelData::d_Z, util::parallel::finalizeMpi(), anonymous_namespace{main.cpp}::fitCrackFromPhi(), util::io::InputParser::getCmdOption(), PeriDEMModel::init(), util::parallel::initMpi(), util::parallel::initNThreads(), PeriDEMModel::integrateStep(), main(), anonymous_namespace{main.cpp}::nodeDamageForFit(), anonymous_namespace{main.cpp}::nodePhi(), util::normalContactStiffness(), data::ModelData::numTimeSteps(), PeriDEMModel::output(), data::ModelData::performOutput(), util::io::print(), geometry::segmentCrossesVerticalLine(), data::ModelData::setCurrentDt(), data::ModelData::shouldOutput(), and data::ModelData::timeStep().

Here is the call graph for this function: