demModel.cpp

/*
 * -------------------------------------------
 * Copyright (c) 2021 - 2024 Prashant K. Jha
 * -------------------------------------------
 * PeriDEM https://github.com/prashjha/PeriDEM
 *
 * Distributed under the Boost Software License, Version 1.0. (See accompanying
 * file LICENSE)
 */

#include "demModel.h"

// utils
#include "particle/baseParticle.h"
#include "material/materialUtil.h"
#include "util/function.h"
#include "util/geomObjectsUtil.h"
#include "util/matrix.h"
#include "util/methods.h"
#include "util/point.h"
#include "inp/pdecks/contactDeck.h"
#include "rw/reader.h"
#include "util/function.h"
#include "util/geom.h"
#include "util/methods.h"
#include "util/randomDist.h"
#include "util/parallelUtil.h"
#include "inp/decks/materialDeck.h"
#include "inp/decks/modelDeck.h"
#include "inp/decks/outputDeck.h"
#include "inp/decks/restartDeck.h"
#include "rw/vtkParticleWriter.h"
#include "rw/vtkParticleReader.h"
#include "fe/elemIncludes.h"
#include "fe/meshUtil.h"

#include <fmt/format.h>
#include <random>

#include <taskflow/taskflow/taskflow.hpp>
#include <taskflow/taskflow/algorithm/for_each.hpp>


model::DEMModel::DEMModel(inp::Input *deck, std::string modelName)
  : ModelData(deck),
    d_name(modelName) {

  // initialize logger
  util::io::initLogger(d_outputDeck_p->d_debug,
                       d_outputDeck_p->d_path + "log.txt");
}

void model::DEMModel::log(std::ostringstream &oss, int priority, bool check_condition, int override_priority,
                          bool screen_out) {
  int op = override_priority == -1 ? priority : override_priority;
  //if (d_outputDeck_p->d_debug > priority)
  if ((check_condition and d_outputDeck_p->d_debug > priority) or d_outputDeck_p->d_debug > op)
    util::io::log(oss, screen_out);
}

void model::DEMModel::log(const std::string &str, int priority, bool check_condition, int override_priority,
                          bool screen_out) {
  int op = override_priority == -1 ? priority : override_priority;  
  if ((check_condition and d_outputDeck_p->d_debug > priority) or d_outputDeck_p->d_debug > op)
    util::io::log(str, screen_out);
}

void model::DEMModel::run(inp::Input *deck) {

  // initialize data
  init();

  // check for restart
  if (d_modelDeck_p->d_isRestartActive)
    restart(deck);

  // integrate in time
  integrate();

  // close
  close();
}

void model::DEMModel::restart(inp::Input *deck) {

  log(d_name + ": Restarting the simulation\n");

  // set time step to step specified in restart deck
  d_n = d_restartDeck_p->d_step;
  d_time = double(d_n) * d_modelDeck_p->d_dt;
  log(fmt::format("  Restart step = {}, time = {:.6f} \n", d_n, d_time));

  // get backup of reference configuration
  std::vector<util::Point> x_ref(d_x.size(), util::Point());
  for (auto &x : d_x)
    x_ref.push_back(x);

  // read displacement and velocity from restart file
  log("  Reading data from restart file = " + d_restartDeck_p->d_file + " \n");
  auto reader = rw::reader::VtkParticleReader(d_restartDeck_p->d_file);
  reader.readNodes(this);
}

void model::DEMModel::close() {
  if (d_ppFile.is_open())
    d_ppFile.close();
}

void model::DEMModel::init() {

  // init time step
  d_n = 0;
  d_time = 0.;
  if (d_outputDeck_p->d_dtTestOut == 0)
    d_outputDeck_p->d_dtTestOut = d_outputDeck_p->d_dtOut / 10;
  d_infoN = d_outputDeck_p->d_dtOut;

  // debug/information variables
  {
    appendKeyData("debug_once", -1);
    appendKeyData("update_contact_neigh_search_params_init_call_count", 0);
    appendKeyData("tree_compute_time", 0);
    appendKeyData("contact_compute_time", 0);
    appendKeyData("contact_neigh_update_time", 0);
    appendKeyData("peridynamics_neigh_update_time", 0);
    appendKeyData("pd_compute_time", 0);
    appendKeyData("extf_compute_time", 0);
    appendKeyData("integrate_compute_time", 0);
    appendKeyData("pt_cloud_update_time", 0);
    appendKeyData("avg_tree_update_time", 0);
    appendKeyData("avg_contact_neigh_update_time", 0);
    appendKeyData("avg_contact_force_time", 0);
    appendKeyData("avg_peridynamics_force_time", 0);
    appendKeyData("avg_extf_compute_time", 0);
    appendKeyData("pen_dist", 0);
    appendKeyData("max_y", 0);
    appendKeyData("contact_area_radius", 0);
  }


  auto t1 = steady_clock::now();
  auto t2 = steady_clock::now();
  log(d_name + ": Initializing objects.\n");

  // create particles
  log(d_name + ": Creating particles.\n");
  createParticles();

  log(d_name + ": Creating maximum velocity data for particles.\n");
  d_maxVelocityParticlesListTypeAll
          = std::vector<double>(d_particlesListTypeAll.size(), 0.);
  d_maxVelocity = util::methods::max(d_maxVelocityParticlesListTypeAll);

  // setup contact
  if (d_input_p->isMultiParticle()) {
    log(d_name + ": Setting up contact.\n");
    setupContact();
  }

  // setup element-node connectivity data if needed
  log(d_name + ": Setting up element-node connectivity data for strain/stress.\n");
  setupQuadratureData();

  // create search object
  log(d_name + ": Creating neighbor search tree.\n");

  // create tree object
  d_nsearch_p = std::make_unique<NSearch>(d_x, d_outputDeck_p->d_debug);

  // setup tree
  double set_tree_time = d_nsearch_p->setInputCloud();
  log(fmt::format("{}: Tree setup time (ms) = {}. \n", d_name, set_tree_time));

  // create neighborlists
  log(d_name + ": Creating neighborlist for peridynamics.\n");
  t1 = steady_clock::now();
  updatePeridynamicNeighborlist();
  t2 = steady_clock::now();
  appendKeyData("peridynamics_neigh_update_time", util::methods::timeDiff(t1, t2));

  if (d_input_p->isMultiParticle()) {
    log(d_name + ": Creating neighborlist for contact.\n");
    d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
    d_contNeighSearchRadius = d_pDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
    t1 = steady_clock::now();
    updateContactNeighborlist();
    t2 = steady_clock::now();
    appendKeyData("contact_neigh_update_time", util::methods::timeDiff(t1, t2));
  }

  // create peridynamic bonds
  log(d_name + ": Creating peridynamics bonds.\n");
  d_fracture_p = std::make_unique<geometry::Fracture>(&d_x, &d_neighPd);

  // compute quantities in state-based simulations
  log(d_name + ": Compute state-based peridynamic quantities.\n");
  material::computeStateMx(this, true);

  // initialize loading class
  log(d_name + ": Initializing displacement loading object.\n");
  d_uLoading_p =
      std::make_unique<loading::ParticleULoading>(d_pDeck_p->d_dispDeck);
  for (auto &p : d_particlesListTypeAll)
    d_uLoading_p->setFixity(p);

  log(d_name + ": Initializing force loading object.\n");
  d_fLoading_p =
      std::make_unique<loading::ParticleFLoading>(d_pDeck_p->d_forceDeck);

  // if all dofs of particle is fixed, then mark it so that we do not
  // compute force
  // MAYBE NOT as we may be interested in reaction forces
  //  for (auto &p : d_particlesListTypeAll)
  //    p->checkFixityForForce(); // TODO implement

  // if this is a two-particle test, we set the force calculation off in
  // first particle
  if (d_pDeck_p->d_testName == "two_particle") {
    d_particlesListTypeAll[0]->d_computeForce = false;
  }

  log(fmt::format("{}: Total particles = {}. \n",
                  d_name, d_particlesListTypeAll.size()));

  for (const auto &p : d_particlesListTypeAll)
    if (!p->d_computeForce)
      log(fmt::format("{}: Force OFF in Particle i = {}. \n", d_name, p->getId()));

  log(d_name + ": Creating list of nodes on which force is to be computed.\n");
  // TODO for now we simply look at particle/wall and check if we compute
  //  force on any of its node. Later, one can have control on individual
  //  nodes of particle/wall and remove from d_fCompNodes if no force is to
  //  be computed on them
  for (size_t i = 0; i < d_x.size(); i++) {
    const auto &ptId = d_ptId[i];
    const auto &pi = getParticleFromAllList(ptId);
    if (pi->d_computeForce) {
      d_fContCompNodes.push_back(i);
      d_fPdCompNodes.push_back(i);
    }
  }

  // initialize remaining fields (if any)
  d_Z = std::vector<float>(d_x.size(), 0.);

  t2 = steady_clock::now();
  log(fmt::format("{}: Total setup time (ms) = {}. \n",
                  d_name, util::methods::timeDiff(t1, t2)));

  // compute complexity information
  size_t free_dofs = 0;
  for (const auto &f : d_fix) {
    for (size_t dof = 0; dof < 3; dof++)
      if (util::methods::isFree(f, dof))
        free_dofs++;
  }
  log(fmt::format("{}: Computational complexity information \n"
                  "  Total number of particles = {}, number of "
                  "particles = {}, number of walls = {}, \n"
                  "  number of dofs = {}, number of free dofs = {}. \n",
                  d_name, d_particlesListTypeAll.size(),
                  d_particlesListTypeParticle.size(),
                  d_particlesListTypeWall.size(),
                  3 * d_x.size(),
                  free_dofs));
}

void model::DEMModel::integrate() {

  // perform output at the beginning
  if (d_n == 0 && d_outputDeck_p->d_performOut) {
    log(fmt::format("{}: Output step = {}, time = {:.6f} \n", d_name, d_n, d_time),
        2);
    output();
  }

  // apply initial condition
  if (d_n == 0)
    applyInitialCondition();

  // apply loading
  computeExternalDisplacementBC();
  computeForces();

  for (size_t i = d_n; i < d_modelDeck_p->d_Nt; i++) {

    log(fmt::format("{}: Time step: {}, time: {:8.6f}, steps completed = {}%\n",
                    d_name,
                    i,
                    d_time,
                    float(i) * 100. / d_modelDeck_p->d_Nt),
        2, d_n % d_infoN == 0, 3);
    
    auto t1 = steady_clock::now();
    log("Integrating\n", false, 0, 3);
    integrateStep();
    double integrate_time =
        util::methods::timeDiff(t1, steady_clock::now());

    appendKeyData("integrate_compute_time", integrate_time, true);

    log(fmt::format("  Integration time (ms) = {}\n", integrate_time), 2, d_n % d_infoN == 0, 3);

    if (d_pDeck_p->d_testName == "two_particle") {

      // compute location of maximum shear stress and also compute
      // penetration length
      auto msg = ppTwoParticleTest();
      log(msg, 2, d_n % d_infoN == 0, 3);
    } else if (d_pDeck_p->d_testName == "compressive_test") {
      auto msg = ppCompressiveTest();
      log(msg, 2, d_n % d_infoN == 0, 3);
    }

    // handle general output
    if ((d_n % d_outputDeck_p->d_dtOut == 0) &&
        (d_n >= d_outputDeck_p->d_dtOut) && d_outputDeck_p->d_performOut) {
      output();
    }

    // check for stop
    checkStop();

  } // loop over time steps

  log(fmt::format(
          "{}: Total compute time information (s) \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n"
          "  {:22s} = {:8.2f} \n",
          d_name,
          "Time integration", getKeyData("integrate_compute_time") * 1.e-6,
          "Peridynamics force", getKeyData("pd_compute_time") * 1.e-6,
          "Contact force", getKeyData("contact_compute_time") * 1.e-6,
          "Search tree update", getKeyData("tree_compute_time") * 1.e-6,
          "External force", getKeyData("extf_compute_time") * 1.e-6)
          );
}

void model::DEMModel::integrateStep() {
  if (d_modelDeck_p->d_timeDiscretization == "central_difference")
    integrateCD();
  else if (d_modelDeck_p->d_timeDiscretization == "velocity_verlet")
    integrateVerlet();
}

void model::DEMModel::integrateCD() {

  // update velocity and displacement
  d_currentDt = d_modelDeck_p->d_dt;
  const auto dim = d_modelDeck_p->d_dim;

  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;

  // update current position, displacement, and velocity of nodes
  taskflow.for_each_index(
    (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
      [this, dim](std::size_t II) {
        auto i = this->d_fPdCompNodes[II];

        const auto rho = this->getDensity(i);
        const auto &fix = this->d_fix[i];

        for (int dof = 0; dof < dim; dof++) {
          if (util::methods::isFree(fix, dof)) {
            this->d_v[i][dof] += (this->d_currentDt / rho) * this->d_f[i][dof];
            this->d_u[i][dof] += this->d_currentDt * this->d_v[i][dof];
            this->d_x[i][dof] += this->d_currentDt * this->d_v[i][dof];
          }
        }

        this->d_vMag[i] = this->d_v[i].length();
      } // loop over nodes
  ); // for_each

  executor.run(taskflow).get();

  // advance time
  d_n++;
  d_time += d_currentDt;

  // update displacement bc
  computeExternalDisplacementBC();

  // compute force
  computeForces();
}

void model::DEMModel::integrateVerlet() {

  // update velocity and displacement
  d_currentDt = d_modelDeck_p->d_dt;
  const auto dim = d_modelDeck_p->d_dim;

  // update current position, displacement, and velocity of nodes
  {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;

    taskflow.for_each_index(
      (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
        [this, dim](std::size_t II) {
          auto i = this->d_fPdCompNodes[II];

          const auto rho = this->getDensity(i);
          const auto &fix = this->d_fix[i];

          for (int dof = 0; dof < dim; dof++) {
            if (util::methods::isFree(fix, dof)) {
              this->d_v[i][dof] += 0.5 * (this->d_currentDt / rho) * this->d_f[i][dof];
              this->d_u[i][dof] += this->d_currentDt * this->d_v[i][dof];
              this->d_x[i][dof] += this->d_currentDt * this->d_v[i][dof];
            }

            this->d_vMag[i] = this->d_v[i].length();
          }
        } // loop over nodes
    ); // for_each

    executor.run(taskflow).get();
  }

  // advance time
  d_n++;
  d_time += d_currentDt;

  // update displacement bc
  computeExternalDisplacementBC();

  // compute force
  computeForces();

  // update velocity of nodes
  {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;

    taskflow.for_each_index(
      (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
      [this, dim](std::size_t II) {
        auto i = this->d_fPdCompNodes[II];

        const auto rho = this->getDensity(i);
        const auto &fix = this->d_fix[i];
        for (int dof = 0; dof < dim; dof++) {
          if (util::methods::isFree(fix, dof)) {
            this->d_v[i][dof] += 0.5 * (this->d_currentDt / rho) * this->d_f[i][dof];
          }

          this->d_vMag[i] = this->d_v[i].length();
        }
      } // loop over nodes
    ); // for_each

    executor.run(taskflow).get();
  }
}

void model::DEMModel::computeForces() {

  bool dbg_condition = d_n % d_infoN == 0;

  log("  Compute forces \n", 2, dbg_condition, 3);

  // reset force
  auto t1 = steady_clock::now();
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;

  taskflow.for_each_index(
    (std::size_t) 0, d_x.size(), (std::size_t) 1,
      [this](std::size_t i) { this->d_f[i] = util::Point(); }
  ); // for_each

  executor.run(taskflow).get();
  auto force_reset_time = util::methods::timeDiff(t1, steady_clock::now());

  // compute peridynamic forces
  t1 = steady_clock::now();
  computePeridynamicForces();
  auto pd_time = util::methods::timeDiff(t1, steady_clock::now());
  appendKeyData("pd_compute_time", pd_time);
  appendKeyData("avg_peridynamics_force_time", pd_time/d_infoN);

  float current_contact_neigh_update_time = 0;
  float contact_time = 0;
  if (d_input_p->isMultiParticle()) {
    // update contact neighborlist
    t1 = steady_clock::now();
    updateContactNeighborlist();
    auto current_contact_neigh_update_time = util::methods::timeDiff(t1,
                                                                     steady_clock::now());
    appendKeyData("contact_neigh_update_time",
                  current_contact_neigh_update_time);
    appendKeyData("avg_contact_neigh_update_time",
                  current_contact_neigh_update_time / d_infoN);

    // compute contact forces between particles
    t1 = steady_clock::now();
    computeContactForces();
    auto contact_time = util::methods::timeDiff(t1, steady_clock::now());
    appendKeyData("contact_compute_time", contact_time);
    appendKeyData("avg_contact_force_time", contact_time / d_infoN);
  }

  // Compute external forces
  t1 = steady_clock::now();
  computeExternalForces();
  auto extf_time = util::methods::timeDiff(t1, steady_clock::now());
  appendKeyData("extf_compute_time", extf_time);
  appendKeyData("avg_extf_compute_time", extf_time/d_infoN);

  // output avg time info
  if (dbg_condition) {
    if (d_input_p->isMultiParticle()) {
      log(fmt::format("    Avg time (ms): \n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n",
                      "tree update", size_t(getKeyData("avg_tree_update_time")),
                      "contact neigh update",
                      size_t(getKeyData("avg_contact_neigh_update_time")),
                      "contact force",
                      size_t(getKeyData("avg_contact_force_time")),
                      "total contact", size_t(getKeyData("avg_tree_update_time")
                                              + getKeyData(
                          "avg_contact_neigh_update_time")
                                              + getKeyData(
                          "avg_contact_force_time")),
                      "peridynamics force",
                      size_t(getKeyData("avg_peridynamics_force_time")),
                      "external force",
                      size_t(getKeyData("avg_extf_compute_time") / d_infoN)),
          2, dbg_condition, 3);

      appendKeyData("avg_tree_update_time", 0.);
      appendKeyData("avg_contact_neigh_update_time", 0.);
      appendKeyData("avg_contact_force_time", 0.);
      appendKeyData("avg_peridynamics_force_time", 0.);
      appendKeyData("avg_extf_compute_time", 0.);
    }
    else {
      log(fmt::format("    Avg time (ms): \n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n",
                      "peridynamics force", size_t(getKeyData("avg_peridynamics_force_time")),
                      "external force", size_t(getKeyData("avg_extf_compute_time")/d_infoN)),
          2, dbg_condition, 3);

      appendKeyData("avg_peridynamics_force_time", 0.);
      appendKeyData("avg_extf_compute_time", 0.);
    }
  }

  log(fmt::format("    {:50s} = {:8d} \n",
                  "Force reset time (ms)",
                  size_t(force_reset_time)
      ),
      2, dbg_condition, 3);

  log(fmt::format("    {:50s} = {:8d} \n",
                  "External force time (ms)",
                  size_t(extf_time)
      ),
      2, dbg_condition, 3);

  log(fmt::format("    {:50s} = {:8d} \n",
                  "Peridynamics force time (ms)",
                  size_t(pd_time)
      ),
      2, dbg_condition, 3);

  if (d_input_p->isMultiParticle()) {

    log(fmt::format("    {:50s} = {:8d} \n",
                    "Point cloud update time (ms)",
                    size_t(getKeyData("pt_cloud_update_time"))
        ),
        2, dbg_condition, 3);

    log(fmt::format("    {:50s} = {:8d} \n",
                    "Contact neighborlist update time (ms)",
                    size_t(current_contact_neigh_update_time)
        ),
        2, dbg_condition, 3);

    log(fmt::format("    {:50s} = {:8d} \n",
                    "Contact force time (ms)",
                    size_t(contact_time)
        ),
        2, dbg_condition, 3);
  }

}

void model::DEMModel::computePeridynamicForces() {

  log("    Computing peridynamic force \n", 3);

  const auto dim = d_modelDeck_p->d_dim;
  const bool is_state = d_particlesListTypeAll[0]->getMaterial()->isStateActive();

  // compute state-based helper quantities
  if (is_state) {

    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;

    taskflow.for_each_index(
      (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1, [this](std::size_t II) {
        auto i = this->d_fPdCompNodes[II];

        const auto rho = this->getDensity(i);
        const auto &fix = this->d_fix[i];
        const auto &ptId = this->getPtId(i);
        auto &pi = this->getParticleFromAllList(ptId);

        if (pi->d_material_p->isStateActive()) {

          const double horizon = pi->getHorizon();
          const double mesh_size = pi->getMeshSize();
          const auto &xi = this->d_xRef[i];
          const auto &ui = this->d_u[i];

          // update bond state and compute thetax
          const auto &m = this->d_mX[i];
          double theta = 0.;

          // upper and lower bound for volume correction
          auto check_up = horizon + 0.5 * mesh_size;
          auto check_low = horizon - 0.5 * mesh_size;

          size_t k = 0;
          for (size_t j : this->d_neighPd[i]) {

            const auto &xj = this->d_xRef[j];
            const auto &uj = this->d_u[j];
            double rji = (xj - xi).length();
            // double rji = std::sqrt(this->d_neighPdSqdDist[i][k]);
            double change_length = (xj - xi + uj - ui).length() - rji;

            // step 1: update the bond state
            double s = change_length / rji;
            double sc = pi->d_material_p->getSc(rji);

            // get fracture state, modify, and set
            auto fs = this->d_fracture_p->getBondState(i, k);
            if (!fs && util::isGreater(std::abs(s), sc + 1.0e-10))
              fs = true;
            this->d_fracture_p->setBondState(i, k, fs);

            if (!fs) {

              // get corrected volume of node j
              auto volj = this->d_vol[j];

              if (util::isGreater(rji, check_low))
                volj *= (check_up - rji) / mesh_size;

              theta += rji * change_length * pi->d_material_p->getInfFn(rji) *
                        volj;
            } // if bond is not broken

            k += 1;
          } // loop over neighbors

          this->d_thetaX[i] = 3. * theta / m;
        } // if it is state-based
      } // loop over nodes
    ); // for_each

    executor.run(taskflow).get();
  }

  // compute the internal forces
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;

  taskflow.for_each_index(
    (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1, [this](std::size_t II) {
      auto i = this->d_fPdCompNodes[II];

      // local variable to hold force
      util::Point force_i = util::Point();
      double scalar_f = 0.;

      // for damage
      float Zi = 0.;

      const auto rhoi = this->getDensity(i);
      const auto &ptIdi = this->getPtId(i);
      auto &pi = this->getParticleFromAllList(ptIdi);

      const double horizon = pi->getHorizon();
      const double mesh_size = pi->getMeshSize();
      const auto &xi = this->d_xRef[i];
      const auto &ui = this->d_u[i];
      const auto &mi = this->d_mX[i];
      const auto &thetai = this->d_thetaX[i];

      // upper and lower bound for volume correction
      auto check_up = horizon + 0.5 * mesh_size;
      auto check_low = horizon - 0.5 * mesh_size;

      // loop over neighbors
      {
        size_t k = 0;
        for (size_t j : this->d_neighPd[i]) {
          auto fs = this->d_fracture_p->getBondState(i, k);
          const auto &xj = this->d_xRef[j];
          const auto &uj = this->d_u[j];
          auto volj = this->d_vol[j];
          double rji = (xj - xi).length();
          double Sji = pi->d_material_p->getS(xj - xi, uj - ui);

          if (!fs) {
            const auto &mj = this->d_mX[j];
            const auto &thetaj = this->d_thetaX[j];

            // get corrected volume of node j
            if (util::isGreater(rji, check_low))
              volj *= (check_up - rji) / mesh_size;

            // handle two cases differently
            if (pi->d_material_p->isStateActive()) {

              auto ef_i =
                  pi->d_material_p->getBondEF(rji, Sji, fs, mi, thetai);
              auto ef_j =
                  pi->d_material_p->getBondEF(rji, Sji, fs, mj, thetaj);

              // compute the contribution of bond force to force at i
              scalar_f = (ef_i.second + ef_j.second) * volj;

              force_i += scalar_f * pi->d_material_p->getBondForceDirection(
                                        xj - xi, uj - ui);
            } // if state-based
            else {

              // Debug
              bool break_bonds = true;

              auto ef =
                  pi->d_material_p->getBondEF(rji, Sji, fs, break_bonds);
              this->d_fracture_p->setBondState(i, k, fs);

              // compute the contribution of bond force to force at i
              scalar_f = ef.second * volj;

              force_i += scalar_f * pi->d_material_p->getBondForceDirection(
                                        xj - xi, uj - ui);
            } // if bond-based
          }   // if bond not broken
          else {
            // add normal contact force
            auto yji = xj + uj - (xi + ui);
            auto Rji = yji.length();
            scalar_f = pi->d_Kn * volj * (Rji - pi->d_Rc) / Rji;
            if (scalar_f > 0.)
              scalar_f = 0.;
            force_i += scalar_f * yji;
          } // if bond is broken

          // calculate damage
          auto Sc = pi->d_material_p->getSc(rji);
          if (util::isGreater(std::abs(Sji / Sc), Zi))
            Zi = std::abs(Sji / Sc);

          k++;
        } // loop over neighbors

      } // peridynamic force

      // update force (we remove any force from
      // previous steps and add peridynamics force)
      this->d_f[i] = force_i;

      this->d_Z[i] = Zi;
    }
  ); // for_each

  executor.run(taskflow).get();
}

void model::DEMModel::computeExternalForces() {

  log("    Computing external force \n", 3);

  auto gravity = d_pDeck_p->d_gravity;

  if (gravity.length() > 1.0E-8) {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;

    taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t)1, [this, gravity](std::size_t i) {
          this->d_f[i] += this->getDensity(i) * gravity;
      } // loop over particles
    ); // for_each

    executor.run(taskflow).get();
  }

  //
  for (auto &p : d_particlesListTypeAll)
    d_fLoading_p->apply(d_time, p); // applied in parallel
}

void model::DEMModel::computeExternalDisplacementBC() {
  log("    Computing external displacement bc \n", 3);
  for (auto &p : d_particlesListTypeAll)
    d_uLoading_p->apply(d_time, p); // applied in parallel
}

void model::DEMModel::computeContactForces() {

  log("    Computing normal contact force \n", 3);

  // Description:
  // 1. Normal contact is applied between nodes of particles and walls
  // 2. Normal damping is applied between particle centers
  // 3. Normal damping is applied between nodes of particle and wall pairs

  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;

  taskflow.for_each_index((std::size_t) 0,
                          d_fContCompNodes.size(),
                          (std::size_t) 1,
                          [this](std::size_t II) {

                              auto i = this->d_fContCompNodes[II];

                              // local variable to hold force
                              util::Point force_i = util::Point();
                              double scalar_f = 0.;

                              const auto &ptIdi = this->getPtId(i);
                              auto &pi = this->getParticleFromAllList(ptIdi);
                              double horizon = pi->d_material_p->getHorizon();
                              double search_r = this->d_maxContactR;

                              // particle data
                              double rhoi = pi->getDensity();

                              const auto &yi = this->d_x[i]; // current coordinates
                              const auto &ui = this->d_u[i];
                              const auto &vi = this->d_v[i];
                              const auto &voli = this->d_vol[i];

                              const std::vector<size_t> &neighs = this->d_neighC[i];

                              if (neighs.size() > 0) {

                                for (const auto &j_id: neighs) {

                                  //auto &j_id = neighs[j];
                                  const auto &yj = this->d_x[j_id]; // current coordinates
                                  double Rji = (yj - yi).length();
                                  auto &ptIdj = this->d_ptId[j_id];
                                  auto &pj = this->getParticleFromAllList(ptIdj);
                                  double rhoj = pj->getDensity();

                                  bool both_walls =
                                          (pi->getTypeIndex() == 1 and pj->getTypeIndex() == 1);

                                  if (j_id != i) {
                                    if (ptIdj != ptIdi && !both_walls) {

                                      // apply particle-particle or particle-wall contact here
                                      const auto &contact =
                                              d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);

                                      if (util::isLess(Rji, contact.d_contactR)) {

                                        auto yji = this->d_x[j_id] - yi;
                                        auto volj = this->d_vol[j_id];
                                        auto vji = this->d_v[j_id] - vi;

                                        // resolve velocity vector in normal and tangential components
                                        auto en = yji / Rji;
                                        auto vn_mag = (vji * en);
                                        auto et = vji - vn_mag * en;
                                        if (util::isGreater(et.length(), 0.))
                                          et = et / et.length();
                                        else
                                          et = util::Point();

                                        // Formula using bulk modulus and horizon
                                        scalar_f = contact.d_Kn * (Rji - contact.d_contactR) *
                                                   volj; // divided by voli
                                        if (scalar_f > 0.)
                                          scalar_f = 0.;
                                        force_i += scalar_f * en;

                                        // compute friction force (since f < 0, |f| = -f)
                                        force_i += contact.d_mu * scalar_f * et;

                                        // if particle-wall pair, apply damping contact here <--
                                        // doesnt seem to work
                                        bool node_lvl_damp = false;
                                        // if (pi->getTypeIndex() == 0 and pj->getTypeIndex() == 1)
                                        //   node_lvl_damp = true;

                                        if (node_lvl_damp) {
                                          // apply damping at the node level
                                          auto meq = util::equivalentMass(rhoi * voli, rhoj * volj);
                                          auto beta_n =
                                                  contact.d_betan *
                                                  std::sqrt(contact.d_kappa * contact.d_contactR * meq);

                                          auto &pii = this->d_particlesListTypeAll[pi->getId()];
                                          vji = this->d_v[j_id] - pii->getVCenter();
                                          vn_mag = (vji * en);
                                          if (vn_mag > 0.)
                                            vn_mag = 0.;
                                          force_i += beta_n * vn_mag * en / voli;
                                        }
                                      } // within contact radius
                                    }   // particle-particle contact
                                  }     // if j_id is not i
                                }       // loop over neighbors
                              }         // contact neighbor

                              this->d_f[i] += force_i;
                          }
  ); // for_each

  executor.run(taskflow).get();


  // damping force
  log("    Computing normal damping force \n", 3);
  for (auto &pi : d_particlesListTypeParticle) {

    auto pi_id = pi->getId();

    double Ri = pi->d_geom_p->boundingRadius();
    double vol_pi = M_PI * Ri * Ri;
    auto pi_xc = pi->getXCenter();
    auto pi_vc = pi->getVCenter();
    auto rhoi = pi->getDensity();
    util::Point force_i = util::Point();

    // particle-particle
    for (auto &pj : this->d_particlesListTypeParticle) {
      if (pj->getId() != pi->getId()) {
        auto Rj = pj->d_geom_p->boundingRadius();
        auto xc_ji = pj->getXCenter() - pi_xc;
        auto dist_xcji = xc_ji.length();

        const auto &contact = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);

        if (util::isLess(dist_xcji, Rj + Ri + 1.01 * contact.d_contactR)) {

          auto vol_pj = M_PI * Rj * Rj;
          auto rhoj = pj->getDensity();
          // equivalent mass
          auto meq = util::equivalentMass(rhoi * vol_pi, rhoj * vol_pj);

          // beta_n
          auto beta_n = contact.d_betan *
                        std::sqrt(contact.d_kappa * contact.d_contactR * meq);

          // center-center vector
          auto hat_xc_ji = util::Point();
          if (util::isGreater(dist_xcji, 0.))
            hat_xc_ji = xc_ji / dist_xcji;
          else
            hat_xc_ji = util::Point();

          // center-center velocity
          auto vc_ji = pj->getVCenter() - pi_vc;
          auto vc_mag = vc_ji * hat_xc_ji;
          if (vc_mag > 0.)
            vc_mag = 0.;

          // force at node of pi
          force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
        } // if within contact distance
      }   // if not same particles
    }     // other particles

    // particle-wall
    // Step 1: Create list of wall nodes that are within the Rc distance
    // of at least one of the particle
    // This is done already in updateContactNeighborList()

    // step 2 - condensed wall nodes into one vector (has to be done serially
    d_neighWallNodesCondensed[pi->getId()].clear();
    {
      for (size_t j=0; j<d_neighWallNodes[pi_id].size(); j++) {

        const auto &j_id = pi->getNodeId(j);
        const auto &yj = this->d_x[j_id];

        for (size_t k=0; k<d_neighWallNodes[pi_id][j].size(); k++) {

          const auto &k_id = d_neighWallNodes[pi_id][j][k];
          const auto &pk = d_particlesListTypeAll[d_ptId[k_id]];

          double Rjk = (this->d_x[k_id] - yj).length();

          const auto &contact =
                  d_cDeck_p->getContact(pi->d_zoneId, pk->d_zoneId);

          if (util::isLess(Rjk, contact.d_contactR))
            util::methods::addToList(k_id, d_neighWallNodesCondensed[pi_id]);

        } // loop over k
      } // loop over j
    } // step 2

    // now loop over wall nodes and add force to center of particle
    for (auto &j : d_neighWallNodesCondensed[pi_id]) {

      auto &ptIdj = this->d_ptId[j];
      auto &pj = this->d_particlesListTypeAll[ptIdj];
      auto rhoj = pj->getDensity();
      auto volj = this->d_vol[j];
      auto meq = rhoi * vol_pi;
      //auto meq = util::equivalentMass(rhoi * vol_pi, rhoj * volj);

      const auto &contact
              = d_cDeck_p->getContact(pi->d_zoneId, pj->d_zoneId);

      // beta_n
      auto beta_n = contact.d_betan *
                    std::sqrt(contact.d_kappa * contact.d_contactR * meq);

      // center-node vector
      auto xc_ji = this->d_x[j] - pi_xc;
      auto hat_xc_ji = util::Point();
      if (util::isGreater(xc_ji.length(), 0.))
        hat_xc_ji = xc_ji / xc_ji.length();

      // center-node velocity
      auto vc_ji = this->d_v[j] - pi_vc;
      auto vc_mag = vc_ji * hat_xc_ji;
      if (vc_mag > 0.)
        vc_mag = 0.;

      // force at node of pi
      force_i += beta_n * vc_mag * hat_xc_ji / vol_pi;
    }

    // distribute force_i to all nodes of particle pi
    {
      tf::Executor executor(util::parallel::getNThreads());
      tf::Taskflow taskflow;

      taskflow.for_each_index((std::size_t) 0, pi->getNumNodes(), (std::size_t) 1,
                              [this, pi, force_i](std::size_t i) {
                                  this->d_f[pi->getNodeId(i)] += force_i;
                              }
      ); // for_each

      executor.run(taskflow).get();
    }
  } // loop over particle for damping
}

void model::DEMModel::applyInitialCondition() {

  log("Applying initial condition \n", 3);

  if (!d_pDeck_p->d_icDeck.d_icActive)
    return;

  const auto ic_v = d_pDeck_p->d_icDeck.d_icVec;
  const auto ic_p_list = d_pDeck_p->d_icDeck.d_pList;

  // add specified velocity to particle
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;

  taskflow.for_each_index((std::size_t) 0,
                          ic_p_list.size(),
                          (std::size_t) 1,
                          [this, ic_v, ic_p_list](std::size_t i) {
      auto &p = this->d_particlesListTypeAll[ic_p_list[i]];

      // velocity
      for (size_t j = 0; j < p->getNumNodes(); j++)
        p->setVLocal(j, ic_v);
    } // loop over particles
  ); // for_each

  executor.run(taskflow).get();
}

void model::DEMModel::createParticles() {

  d_particlesListTypeParticle.resize(0);
  d_particlesListTypeAll.resize(0);
  d_particlesListTypeWall.resize(0);
  d_referenceParticles.clear();

  // loop over all particle zones
  for (size_t z = 0; z < d_pDeck_p->d_particleZones.size(); z++) {

    // does this particle zone correspond to particle or wall
    bool is_wall = d_pDeck_p->d_particleZones[z].d_isWall;
    std::string particle_type = d_pDeck_p->d_zoneToParticleORWallDeck[z].first;
    if (is_wall and particle_type != "wall") {
      std::cerr << fmt::format("Error: String d_zoneToParticleORWallDeck[z].first for zone z = {} "
                               "should be 'wall'.\n", z);
      exit(EXIT_FAILURE);
    }
    if (!is_wall and particle_type != "particle") {
      std::cerr << fmt::format("Error: String d_zoneToParticleORWallDeck[z].first for zone z = {} "
                               "should be 'particle'.\n", z);
      exit(EXIT_FAILURE);
    }

    // get current size of particles data
    auto psize = d_particlesListTypeAll.size();

    // get particle zone
    auto &pz = d_pDeck_p->d_particleZones[z];

    // get zone id
    auto z_id = pz.d_zone.d_zoneId;
    if (z_id != z) {
      std::cerr << fmt::format("Error: d_zoneId = {} in ParticleZone for "
                               "z = {} should be equal to z = {}.\n",
                               z_id, z, z);
      exit(EXIT_FAILURE);
    }

    // read mesh data
    log(d_name + ": Creating mesh for reference particle in zone = " +
        std::to_string(z_id) + "\n");
    std::shared_ptr<fe::Mesh> mesh;
    if (!pz.d_meshDeck.d_createMesh) {
      mesh = std::make_shared<fe::Mesh>(&pz.d_meshDeck);
    }
    else {
      const auto &geomData = pz.d_meshDeck.d_createMeshGeomData;
      if (pz.d_meshDeck.d_createMeshInfo == "uniform"
          and geomData.d_geomName == "rectangle") {

        // get the geometrical details
        std::pair<std::vector<double>, std::vector<double>> box;
        std::vector<size_t> nGrid(3, 0);

        for (size_t i=0; i<3; i++) {
          box.first.push_back(geomData.d_geomParams[i]);
          box.second.push_back(geomData.d_geomParams[i+3]);

          nGrid[i] = size_t((geomData.d_geomParams[i+3] - geomData.d_geomParams[i])/pz.d_meshDeck.d_h);

          std::cout << fmt::format("box.first[i] = {}, "
                                   "box.second[i] = {}, "
                                   "nGrid[i] = {}\n",
                                   box.first[i],
                                   box.second[i],
                                   nGrid[i]);
        }
        fe::Mesh temp_mesh;
        fe::createUniformMesh(&temp_mesh,
                              d_modelDeck_p->d_dim,
                              box,
                              nGrid);
        mesh = std::make_shared<fe::Mesh>(temp_mesh);
      }
      else {
        std::cerr << "Error: Currently, we can only support in-built uniform mesh for rectangles.\n";
        exit(EXIT_FAILURE);
      }
    }

    // create the reference particle
    log(d_name + ": Creating reference particle in zone = " +
        std::to_string(z_id) + "\n");

    // get representative particle for this zone
    auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
    auto rep_geom_params = pz.d_particleGeomData.d_geomParams;

    auto ref_p = std::make_shared<particle::RefParticle>(
            d_referenceParticles.size(),
            static_cast<std::shared_ptr<ModelData>>(this),
            rep_geom_p,
            mesh);

    d_referenceParticles.emplace_back(ref_p);

    // check the particle generation method
    log(d_name + ": Creating particles in zone = " +
                  std::to_string(z_id) + "\n");

    if (pz.d_genMethod == "From_File") {
      createParticlesFromFile(z, ref_p);
    }
    else {
      if (pz.d_createParticleUsingParticleZoneGeomObject or !d_input_p->isMultiParticle()) {
        createParticleUsingParticleZoneGeomObject(z, ref_p);
      }
      else {
        std::cerr << "Error: Particle generation method = " << pz.d_genMethod <<
                  " not recognized.\n";
        exit(1);
      }
    }

    // get new size of data
    auto psize_new = d_particlesListTypeAll.size();

    // store this in zone-info
    d_zInfo.emplace_back(std::vector<size_t>{psize, psize_new, z_id});
  }
}

void model::DEMModel::createParticleUsingParticleZoneGeomObject(
        size_t z,
        std::shared_ptr<particle::RefParticle> ref_p) {

  log(d_name + ": Creating particle using Particle Zone Geometry Object\n", 1);

  // get particle zone
  auto &pz = d_pDeck_p->d_particleZones[z];

  // get zone id
  auto z_id = pz.d_zone.d_zoneId;

  // ref_p has geometry and mesh which will be used in creating this particle
  // we need to create identity transform
  auto p_transform = particle::ParticleTransform();

  // create particle
  auto p = new particle::BaseParticle(
          pz.d_isWall ? "wall" : "particle",
          d_particlesListTypeAll.size(),
          pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
          z_id,
          ref_p->getDimension(),
          pz.d_particleDescription,
          pz.d_isWall,
          pz.d_allDofsConstrained,
          ref_p->getNumNodes(),
          0.,
          static_cast<std::shared_ptr<ModelData>>(this),
          ref_p,
          ref_p->getGeomP(),
          p_transform,
          ref_p->getMeshP(),
          pz.d_matDeck,
          true);

  // push p to list
  if (pz.d_isWall)
    d_particlesListTypeWall.push_back(p);
  else
    d_particlesListTypeParticle.push_back(p);

  d_particlesListTypeAll.push_back(p);

}

void model::DEMModel::createParticlesFromFile(
    size_t z, std::shared_ptr<particle::RefParticle> ref_p) {

  log(d_name + ": Creating particle from file\n", 1);

  // get particle zone
  auto &pz = d_pDeck_p->d_particleZones[z];

  // get zone id
  auto z_id = pz.d_zone.d_zoneId;

  // read file which contains location of centers of particle, zone id, and
  // radius of particle
  std::vector<util::Point> centers;
  std::vector<double> rads;
  std::vector<double> orients;
  if (pz.d_particleFileDataType == "loc_rad") {
    rw::reader::readParticleCsvFile(pz.d_particleFile, d_modelDeck_p->d_dim,
                                    &centers, &rads, z_id);

    util::DistributionSample<UniformDistribution> uniform_dist(
        0., 1., d_modelDeck_p->d_seed);

    if (d_pDeck_p->d_testName == "two_particle") {
      for (size_t i = 0; i < rads.size(); i++)
        orients.push_back((double(i)) * M_PI);
    } else {
      for (size_t i = 0; i < rads.size(); i++)
        orients.push_back(
            util::transform_to_uniform_dist(0., 2. * M_PI, uniform_dist()));
    }
  }
  else if (pz.d_particleFileDataType == "loc_rad_orient") {
    rw::reader::readParticleWithOrientCsvFile(pz.d_particleFile,
                                              d_modelDeck_p->d_dim, &centers,
                                              &rads, &orients, z_id);
  }

  log(fmt::format("zone_id: {}, rads: {}, orients: {}, centers: {} \n", z_id,
                    util::io::printStr(rads), util::io::printStr(orients),
                    util::io::printStr(centers)), 2);

  // get representative particle for this zone
  const auto &rep_geom_p = pz.d_particleGeomData.d_geom_p;
  auto rep_geom_params = pz.d_particleGeomData.d_geomParams;

  // get zone bounding box
  std::pair<util::Point, util::Point> box = rep_geom_p->box();

  size_t p_counter = 0;
  size_t p_old_size = d_particlesListTypeAll.size();
  for (const auto &site : centers) {

    double particle_radius = rads[p_counter];
    double particle_orient = orients[p_counter];

    // create geometrical object
    std::shared_ptr<util::geometry::GeomObject> p_geom;
    createGeometryAtSite(particle_radius,
                         particle_orient,
                         site,
                         rep_geom_params,
                         rep_geom_p,
                         p_geom);

    // create transform
    auto p_transform = particle::ParticleTransform(
            site, util::Point(0., 0., 1.), particle_orient,
            particle_radius / ref_p->getParticleRadius());

    if (p_transform.d_scale < 1.E-8) {
      std::cerr << "Error: check scale in transform. "
                << " Scale: " << particle_radius / ref_p->getParticleRadius()
                << " p rad: " << particle_radius
                << " ref p rad: " << ref_p->getParticleRadius()
                << p_transform.printStr();
      exit(1);
    }

    // finally create dem particle at this site
    //auto particle_id = p_counter + p_old_size;
    auto p = new particle::BaseParticle(
            pz.d_isWall ? "wall" : "particle",
            d_particlesListTypeAll.size(),
            pz.d_isWall ? d_particlesListTypeWall.size() : d_particlesListTypeParticle.size(),
            z_id,
            ref_p->getDimension(),
            pz.d_particleDescription,
            pz.d_isWall,
            pz.d_allDofsConstrained,
            ref_p->getNumNodes(),
            0.,
            static_cast<std::shared_ptr<ModelData>>(this),
            ref_p,
            p_geom,
            p_transform,
            ref_p->getMeshP(),
            pz.d_matDeck,
            true);

    // push p to list
    if (pz.d_isWall)
      d_particlesListTypeWall.push_back(p);
    else
      d_particlesListTypeParticle.push_back(p);

    d_particlesListTypeAll.push_back(p);
    p_counter++;
  }
}

void model::DEMModel::createGeometryAtSite(const double &particle_radius,
                                           const double &particle_orient,
                                           const util::Point &site,
                                           const std::vector<double> &rep_geom_params,
                                           const std::shared_ptr<util::geometry::GeomObject> &rep_geom_p,
                                           std::shared_ptr<util::geometry::GeomObject> &p_geom) {
  std::vector<double> params;
  for (auto x : rep_geom_params)
    params.push_back(x);

  if (util::methods::isTagInList(rep_geom_p->d_name,
                                 util::geometry::acceptable_geometries)) {

    if (util::methods::isTagInList(rep_geom_p->d_name,
                                   {"circle", "sphere", "hexagon",
                                    "triangle", "square", "cube"})) {

      // case - objects requiring four parameters
      // here 'triangle' is a uniform triangle (see constructor of Triangle)
      // 'hexagon' is a hexagon with axis (1, 0, 0)
      size_t num_params = 4;

      if (params.size() < num_params)
        params.resize(num_params);
      params[0] = particle_radius;
      for (int dof = 0; dof < 3; dof++)
        params[dof + 1] = site[dof];
    }
    else if (rep_geom_p->d_name == "drum2d") {

      // case - objects requiring five parameters
      size_t num_params = 5;

      if (params.size() < num_params)
        params.resize(num_params);

      params[0] = particle_radius; // biggger length along x-direction
      params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0]; // neck length along x-direction
      for (int dof = 0; dof < 3; dof++)
        params[dof + 2] = site[dof];
    }
    else if (rep_geom_p->d_name == "rectangle") {

      // case - objects requiring five parameters
      size_t num_params = 5;

      if (params.size() < num_params)
        params.resize(num_params);

      params[0] = particle_radius; // length along x-direction
      params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0]; // length along y-direction
      for (int dof = 0; dof < 3; dof++)
        params[dof + 2] = site[dof];
    }
    else if (rep_geom_p->d_name == "cuboid") {

      // case - objects requiring six parameters

      if (params.size() < 6)
        params.resize(6);

      params[0] = particle_radius; // length is x-direction
      params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0]; // length in y-direction
      params[2] = particle_radius * rep_geom_params[2] / rep_geom_params[0]; // length in z-direction
      for (int dof = 0; dof < 3; dof++)
        params[dof + 2] = site[dof];
    }
  } else {
    std::cerr << fmt::format("Error: PeriDEM supports following type "
                             "of geometries for particles = {}\n",
                             util::io::printStr(util::geometry::acceptable_geometries));
    exit(EXIT_FAILURE);
  }

  // create geometry now
  std::vector<std::string> vec_geom_type;
  std::vector<std::string> vec_geom_flag;
  util::geometry::createGeomObject(rep_geom_p->d_name, params, vec_geom_type,
                                   vec_geom_flag, p_geom,
                                   d_modelDeck_p->d_dim, false);
}

void model::DEMModel::updateGeometryObjectsPostInit() {

  for (auto &p: d_particlesListTypeAll) {
    if (p->d_geom_p->d_name == "null" or
        util::methods::isTagInList("copy_from_container", p->d_geom_p->d_tags)) {
      // update geometry of particle based on bounding box
      auto bbox = p->getMeshP()->getBoundingBox();

      std::string geom_name = "rectangle";
      if ( p->getMeshP()->getDimension() == 3)
        geom_name = "cuboid";

      std::vector<double> geom_params(6, 0.);
      for (size_t i=0; i<3; i++) {
        geom_params[i] = bbox.first[i];
        geom_params[i+3] = bbox.second[i];
      }

      std::vector<std::string> vec_type;
      std::vector<std::string> vec_flag;

      util::geometry::createGeomObject(geom_name,
                                       geom_params,
                                       vec_type,
                                       vec_flag,
                                       p->d_geom_p,
                                       p->getMeshP()->getDimension());
    }
  }
}

void model::DEMModel::setupContact() {

  // loop over all particle zones and get minimum value of mesh size
  size_t c = 0;
  for (const auto *p : d_particlesListTypeAll) {

    auto h = p->getMeshSize();
    if (c == 0) {
      d_hMin = h;
      d_hMax = h;
      c++;
    }

    if (util::isGreater(d_hMin, h))
      d_hMin = h;
    if (util::isGreater(h, d_hMax))
      d_hMax = h;
  }

  log(fmt::format("{}: Contact setup\n  hmin = {:.6f}, hmax = {:.6f} \n",
                  d_name, d_hMin, d_hMax), 1);

  d_maxContactR = 0.;
  // precompute bulk modulus of all zones
  std::vector<double> bulk_modulus;
  // NOTE - d_data.size() and d_zoneVec.size() are equal
  for (size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {

    double kappa_i = d_pDeck_p->d_particleZones[i].d_matDeck.d_matData.d_K;

    if (kappa_i < 0.) {
      std::cerr << "Error: We need bulk modulus provided in input file.\n";
      std::cerr << d_pDeck_p->d_particleZones[i].printStr();
      exit(1);
    }

    bulk_modulus.push_back(kappa_i);
  }

  for (size_t i = 0; i < d_cDeck_p->d_data.size(); i++) {
    for (size_t j = 0; j < d_cDeck_p->d_data.size(); j++) {

      inp::ContactPairDeck *deck = &(d_cDeck_p->d_data[i][j]);

      if (deck->d_computeContactR)
        deck->d_contactR *= d_hMin;

      if (d_maxContactR < deck->d_contactR)
        d_maxContactR = deck->d_contactR;

      // get effective bulk modulus for pair of zones and store it
      deck->d_kappa = util::equivalentMass(bulk_modulus[i], bulk_modulus[j]);

      // Kn
      deck->d_Kn *= deck->d_KnFactor;

      // Beta n
      double log_e = std::log(deck->d_eps);
      deck->d_betan =
          deck->d_betanFactor *
          (-2. * log_e * std::sqrt(1. / (M_PI * M_PI + log_e * log_e)));

      log(fmt::format("  contact_radius = {:.6f}, hmin = {:.6f}, Kn = {:5.3e}, "
                      "Vmax = {:5.3e}, "
                      "betan = {:7.5f}, mu = {:.4f}, kappa = {:5.3e}\n",
                      deck->d_contactR, d_hMin, deck->d_Kn, deck->d_vMax,
                      deck->d_betan, deck->d_mu, deck->d_kappa), 2);
    }
  }
}

void model::DEMModel::setupQuadratureData() {

  if (util::methods::isTagInList("Strain_Stress", d_outputDeck_p->d_outTags)
      or d_modelDeck_p->d_populateElementNodeConnectivity) {

    // read element-node connectivity data if not done
    for (auto &p: d_referenceParticles) {
      auto &particle_mesh_p = p->getMeshP();
      if (!particle_mesh_p->d_encDataPopulated && particle_mesh_p->d_enc.empty()) {
        particle_mesh_p->readElementData(particle_mesh_p->d_filename);
      }
    }

    // setup quadrature point and strain/stress data
    // we need to know size of the data
    size_t totalQuadPoints = 0;
    for (auto &p: d_particlesListTypeAll) {
      const auto &particle_mesh_p = p->getMeshP();

      // get Quadrature
      fe::BaseElem *elem;
      if (particle_mesh_p->getElementType() == util::vtk_type_line)
        elem = new fe::LineElem(d_modelDeck_p->d_quadOrder);
      else if (particle_mesh_p->getElementType() == util::vtk_type_triangle)
        elem = new fe::TriElem(d_modelDeck_p->d_quadOrder);
      else if (particle_mesh_p->getElementType() == util::vtk_type_quad)
        elem = new fe::QuadElem(d_modelDeck_p->d_quadOrder);
      else if (particle_mesh_p->getElementType() == util::vtk_type_tetra)
        elem = new fe::TetElem(d_modelDeck_p->d_quadOrder);
      else {
        std::cerr << fmt::format("Error: Can not compute strain/stress as the element "
                                 "type = {} is not yet supported in this routine.\n", particle_mesh_p->getElementType());
        exit(EXIT_FAILURE);
      }

      p->d_globQuadStart = totalQuadPoints;
      totalQuadPoints += particle_mesh_p->getNumElements() *
                         elem->getNumQuadPoints();
      p->d_globQuadEnd = totalQuadPoints;

      std::cout << fmt::format("p->id() = {}, "
                               "p->d_globQuadStart = {}, "
                               "totalQuadPoints = {}, "
                               "p->d_globQuadEnd = {}",
                               p->getId(), p->d_globQuadStart,
                               particle_mesh_p->getNumElements() *
                               elem->getNumQuadPoints(), p->d_globQuadEnd)
                << std::endl;
    }

    // resize data
    d_xQuadCur.resize(totalQuadPoints);
    d_strain.resize(totalQuadPoints);
    d_stress.resize(totalQuadPoints);
  } // setting up quadrature data
}

void model::DEMModel::updatePeridynamicNeighborlist() {

  d_neighPd.resize(d_x.size());
  // d_neighPdSqdDist.resize(d_x.size());
  auto t1 = steady_clock::now();

  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;

  taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1, [this](std::size_t i) {
      const auto &pi = this->d_ptId[i];
      double search_r = this->d_particlesListTypeAll[pi]->d_material_p->getHorizon();

      std::vector<size_t> neighs;
      std::vector<double> sqr_dist;
      if (this->d_nsearch_p->radiusSearchIncludeTag(this->d_x[i],
                                                    search_r,
                                                    neighs,
                                                    sqr_dist,
                                                    this->d_ptId[i],
                                                    this->d_ptId) > 0) {
        for (std::size_t j = 0; j < neighs.size(); ++j)
          if (neighs[j] != i && this->d_ptId[neighs[j]] == pi) {
            this->d_neighPd[i].push_back(size_t(neighs[j]));
            // this->d_neighPdSqdDist[i].push_back(sqr_dist[j]);
          }
      }
    }
  ); // for_each

  executor.run(taskflow).get();

  auto t2 = steady_clock::now();
  log(fmt::format("{}: Peridynamics neighbor update time = {}\n",
                  d_name, util::methods::timeDiff(t1, t2)), 2);
}

void model::DEMModel::updateContactNeighborlist() {

  auto update = updateContactNeighborSearchParameters();

  if (!update)
    return;

  // update contact neighborlist

  // update the point cloud (make sure that d_x is updated along with displacement)
  auto pt_cloud_update_time = d_nsearch_p->setInputCloud();
  setKeyData("pt_cloud_update_time", pt_cloud_update_time);
  appendKeyData("tree_compute_time", pt_cloud_update_time);
  appendKeyData("avg_tree_update_time", pt_cloud_update_time/d_infoN);

  if (d_neighC.size() != d_x.size())
    d_neighC.resize(d_x.size());

  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;

  taskflow.for_each_index((std::size_t) 0, d_x.size(), (std::size_t) 1,
                          [this](std::size_t i) {

    const auto &pi = this->d_ptId[i];
    const auto &pi_particle = this->d_particlesListTypeAll[pi];

    // search?
    bool perform_search_based_on_particle = true;
    if (pi_particle->d_typeIndex == 1) // wall
      perform_search_based_on_particle = false;

    if (pi_particle->d_allDofsConstrained or !pi_particle->d_computeForce)
      perform_search_based_on_particle = false;

    if (perform_search_based_on_particle) {

      std::vector<size_t> neighs;
      std::vector<double> sqr_dist;

      this->d_neighC[i].clear();

      auto n = this->d_nsearch_p->radiusSearchExcludeTag(
              this->d_x[i],
              this->d_contNeighSearchRadius,
              neighs,
              sqr_dist,
              this->d_ptId[i],
              this->d_ptId);

      if (n > 0) {
        for (auto neigh: neighs) {
          if (neigh != i)
            this->d_neighC[i].push_back(neigh);
        }
      }
    }
}
  ); // for_each

  executor.run(taskflow).get();


  // handle particle-wall neighborlist (based on the d_neighC that we already computed)
  d_neighWallNodes.resize(d_particlesListTypeAll.size());
  d_neighWallNodesDistance.resize(d_particlesListTypeAll.size());
  d_neighWallNodesCondensed.resize(d_particlesListTypeAll.size());

  for (auto &pi : d_particlesListTypeParticle) {

    d_neighWallNodes[pi->getId()].resize(pi->getNumNodes());
    d_neighWallNodesDistance[pi->getId()].resize(pi->getNumNodes());

    // get all wall nodes that are within contact distance to the nodes of this particle
    {
      tf::Executor executor(util::parallel::getNThreads());
      tf::Taskflow taskflow;

      taskflow.for_each_index((std::size_t) 0,
                              pi->getNumNodes(),
                              (std::size_t) 1,
                              [this, &pi](std::size_t i) {

            auto i_glob = pi->getNodeId(i);
            auto yi = this->d_x[i_glob];

            const std::vector<size_t> &neighs = this->d_neighC[i_glob];

            this->d_neighWallNodes[pi->getId()][i].clear();
            this->d_neighWallNodesDistance[pi->getId()][i].clear();

            for (const auto &j_id: neighs) {

              auto &ptIdj = this->d_ptId[j_id];
              auto &pj = this->getParticleFromAllList(
                      ptIdj);

              // we are only interested in nodes from wall
              if (pj->getTypeIndex() == 1) {
                  this->d_neighWallNodes[pi->getId()][i].push_back(j_id);
                  //this->d_neighWallNodesDistance[pi->getId()][i].push_back(Rji);
              }
            }
        }
      ); // for_each

      executor.run(taskflow).get();
    }
  } // loop over particles

}

bool model::DEMModel::updateContactNeighborSearchParameters() {

  // initialize parameters
  if (d_contNeighUpdateInterval == 0 and
      util::isLess(d_contNeighSearchRadius, 1.e-16)) {
    d_contNeighUpdateInterval = d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval;
    d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
    d_contNeighSearchRadius = d_maxContactR * d_pDeck_p->d_pNeighDeck.d_sFactor;
  }

  // at d_n = 0, this function will be called twice because updateContactNeighborlist() will be
  // called twice: one inside init() and second inside computeForces()
  // so to match d_n and d_contNeighTimestepCounter in the initial stage of simulation, we need to handle the special case
  if (d_n == 0) {
    appendKeyData("update_contact_neigh_search_params_init_call_count", 1);

    if (int(getKeyData("update_contact_neigh_search_params_init_call_count")) == 1)
      return true;

    if (int(getKeyData("update_contact_neigh_search_params_init_call_count")) == 2) {
      d_contNeighTimestepCounter++;
      return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
    }
  }

  // handle case of restart
  if (d_modelDeck_p->d_isRestartActive and d_n == d_restartDeck_p->d_step) {
    // assign correct value for restart step
    d_contNeighTimestepCounter = d_n % d_contNeighUpdateInterval;
  }

  if (d_contNeighUpdateInterval == 1) {
    // further optimization of parameters is not possible
    d_contNeighSearchRadius = d_maxContactR;

    // update counter and return condition for contact search
    d_contNeighTimestepCounter++;
    return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
  }

  // check if we should proceed with parameter update
  // param update is done at smaller interval than the search itself to avoid
  // scenarios where particles suddenly move with a high velocity
  size_t update_param_interval =
          d_contNeighUpdateInterval > 5 ? size_t(
                  0.2 * d_contNeighUpdateInterval) : 1;

  // check if we ought to update search parameters; if not, return
  if (d_contNeighTimestepCounter > 0 and d_contNeighTimestepCounter % update_param_interval != 0) {
    // update counter and return condition for contact search
    d_contNeighTimestepCounter++;
    return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
  }

  // first update the maximum velocity in all particles
  for (auto &pi : d_particlesListTypeAll) {
    auto max_v_node = util::methods::maxIndex(d_vMag,
                                              pi->d_globStart, pi->d_globEnd);

    if (max_v_node > pi->d_globEnd or max_v_node < pi->d_globStart) {
      std::cerr << fmt::format("Error: max_v_node = {} for "
                               "particle of id = {} is not in the limit.\n",
                               max_v_node, pi->getId())
                << "Particle info = \n"
                << pi->printStr()
                << "\n\n Magnitude of velocity = "
                << d_vMag[max_v_node] << "\n";
      exit(EXIT_FAILURE);
    }

    d_maxVelocityParticlesListTypeAll[pi->getId()]
            = d_vMag[max_v_node];
  }

  // find max velocity among all particles
  d_maxVelocity = util::methods::max(d_maxVelocityParticlesListTypeAll);

  // now we find the best parameters for contact search
  auto up_interval_old = d_contNeighUpdateInterval;

  // TO ensure that in d_neighUpdateInterval time steps, the search radius is above the
  // distance traveled by object with velocity d_maxVelocity
  // also multiply by a safety factor
  double safety_factor = d_pDeck_p->d_pNeighDeck.d_sFactor > 5 ? d_pDeck_p->d_pNeighDeck.d_sFactor : 10;
  auto max_search_r_from_contact_R = d_pDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
  auto max_search_r = d_maxVelocity * d_currentDt
                      * d_pDeck_p->d_pNeighDeck.d_neighUpdateInterval
                      * safety_factor;


  if (util::isGreater(max_search_r, max_search_r_from_contact_R )) {

    d_contNeighUpdateInterval = size_t(d_maxContactR/(d_maxVelocity * d_currentDt));
    if (up_interval_old > d_contNeighUpdateInterval) {
      // issue warning
      log(fmt::format("Warning: Contact search radius based on velocity is greater than "
                      "the max contact radius.\n"
                      "Warning: Adjusting contact neighborlist update interval.\n"
                      "{:>13} = {:4.6e}, time step = {}, "
                      "velocity-based r = {:4.6e}, max contact r = {:4.6e}\n",
                      "Time", d_time, d_n, max_search_r, max_search_r_from_contact_R),
          2, d_n % d_infoN == 0, 3);
    }

    d_contNeighSearchRadius = max_search_r_from_contact_R;
    // reset time step counter for contact so that the contact list is updated in the current time step
    // and the update cycle starts from the current time step
    d_contNeighTimestepCounter = 0;

    if (d_contNeighUpdateInterval < 1) {
      d_contNeighUpdateInterval = 1;
      d_contNeighSearchRadius = d_maxContactR;
    }
  }
  else {
    // update search radius
    d_contNeighSearchRadius = d_contNeighUpdateInterval < 2 ? d_maxContactR : max_search_r_from_contact_R;
  }

  if (up_interval_old > d_contNeighUpdateInterval) {
    log(fmt::format("    Contact neighbor parameters: \n"
                    "      {:48s} = {:d}\n"
                    "      {:48s} = {:d}\n"
                    "      {:48s} = {:d}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n"
                    "      {:48s} = {:4.6e}\n",
                    "time step", d_n,
                    "contact neighbor update interval",
                    d_contNeighUpdateInterval,
                    "contact neighbor update time step counter",
                    d_contNeighTimestepCounter,
                    "search radius", d_contNeighSearchRadius,
                    "max contact radius", d_maxContactR,
                    "search radius factor", d_pDeck_p->d_pNeighDeck.d_sFactor,
                    "max search r from velocity", max_search_r,
                    "max search r from contact r", max_search_r_from_contact_R,
                    "max velocity", d_maxVelocity),
        2, d_n % d_infoN == 0, 3);
  }

  // update counter and return condition for contact search
  d_contNeighTimestepCounter++;
  return (d_contNeighTimestepCounter - 1) % d_contNeighUpdateInterval == 0;
}

void model::DEMModel::updateNeighborlistCombine() {
  // Not used
  return;
}

void model::DEMModel::output() {

  // write out % completion of simulation at 10% interval
  {
    float p = float(d_n) * 100. / d_modelDeck_p->d_Nt;
    int m = std::max(1, int(d_modelDeck_p->d_Nt / 10));
    if (d_n % m == 0 && int(p) > 0)
      log(fmt::format("{}: Simulation {}% complete\n",
                      d_name, int(p)));
    ;
  }

  log(fmt::format("{}: Output step = {}, time = {:.6f} \n",
                  d_name, d_n, d_time),
      2);

  if (d_outputDeck_p->d_debug > 0 and getKeyData("debug_once") < 0) {

    setKeyData("debug_once", 1);

    size_t nt = 1;
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "*******************************************\n";
    oss << tabS << "Debug various input decks\n\n\n";
    oss << d_modelDeck_p->printStr(nt + 1);
    oss << d_pDeck_p->printStr(nt + 1);
    oss << d_cDeck_p->printStr(nt + 1);
    oss << tabS << "\n\n*******************************************\n";
    oss << tabS << "Debug particle data\n\n\n";
    oss << tabS << "Number of particles = " << d_particlesListTypeAll.size() << std::endl;
    oss << tabS << "Number of particle zones = " << d_zInfo.size() << std::endl;
    for (auto zone : d_zInfo) {
      oss << tabS << "zone of d_zInfo: " << util::io::printStr(zone)
          << std::endl;
    }

    // wall info
    oss << tabS << "Number of walls = " << d_particlesListTypeWall.size() << std::endl;
    for (auto &d_wall : d_particlesListTypeWall)
      oss << tabS << "Number of nodes in wall in zone " << d_wall->d_zoneId
          << " is " << d_wall->getNumNodes() << std::endl;

    oss << tabS << "h_min = " << d_hMin << ", h_max = " << d_hMax << std::endl;

    log(oss, 2);
  } // end of debug

  size_t dt_out = d_outputDeck_p->d_dtOutCriteria;
  std::string out_filename = d_outputDeck_p->d_path + "output_";
  if (d_outputDeck_p->d_tagPPFile.empty())
    out_filename = out_filename + std::to_string(d_n / dt_out);
  else
    out_filename = out_filename + d_outputDeck_p->d_tagPPFile + "_" + std::to_string(d_n / dt_out);

  auto writer = rw::writer::VtkParticleWriter(out_filename);
  if (d_outputDeck_p->d_performFEOut)
    writer.appendMesh(this, d_outputDeck_p->d_outTags);
  else
    writer.appendNodes(this, d_outputDeck_p->d_outTags);

  writer.addTimeStep(d_time);
  writer.close();

  if (util::methods::isTagInList("Strain_Stress", d_outputDeck_p->d_outTags)) {

    // compute current position of quadrature points and strain/stress data
    {
      // if particle mat data is not computed, compute them
      if (d_particlesMatDataList.empty()) {
        for (auto &p: d_particlesListTypeAll) {
          d_particlesMatDataList.push_back(p->getMaterial()->computeMaterialProperties(
                  p->getMeshP()->getDimension()));
        }
      }

      for (auto &p: d_particlesListTypeAll) {

        const auto particle_mesh_p = p->getMeshP();

        fe::getCurrentQuadPoints(particle_mesh_p.get(), d_xRef, d_u, d_xQuadCur,
                                 p->d_globStart,
                                 p->d_globQuadStart,
                                 d_modelDeck_p->d_quadOrder);

        auto p_z_id = p->d_zoneId;
        auto isPlaneStrain = d_pDeck_p->d_particleZones[p_z_id].d_matDeck.d_isPlaneStrain;
        fe::getStrainStress(particle_mesh_p.get(), d_xRef, d_u,
                            isPlaneStrain,
                            d_strain, d_stress,
                            p->d_globStart,
                            p->d_globQuadStart,
                            d_particlesMatDataList[p->getId()].d_nu,
                            d_particlesMatDataList[p->getId()].d_lambda,
                            d_particlesMatDataList[p->getId()].d_mu,
                            true,
                            d_modelDeck_p->d_quadOrder);
      } // for loop over particles
    } // compute strain/stress block

    out_filename = d_outputDeck_p->d_path + "output_strain_";
    if (d_outputDeck_p->d_tagPPFile.empty())
      out_filename = out_filename + std::to_string(d_n / dt_out);
    else
      out_filename = out_filename + d_outputDeck_p->d_tagPPFile + "_" + std::to_string(d_n / dt_out);

    auto writer1 = rw::writer::VtkParticleWriter(out_filename);
    writer1.appendStrainStress(this);
    writer1.addTimeStep(d_time);
    writer1.close();
  }

  // output particle locations to csv file
  if (util::methods::isTagInList("Particle_Locations",
                                 d_outputDeck_p->d_outTags)) {

    out_filename = d_outputDeck_p->d_path + "particle_locations_";
    if (d_outputDeck_p->d_tagPPFile.empty())
      out_filename = out_filename + std::to_string(d_n / dt_out) + ".csv";
    else
      out_filename = out_filename + d_outputDeck_p->d_tagPPFile
                      + "_" + std::to_string(d_n / dt_out) + ".csv";

    std::ofstream oss(out_filename);
    oss << "i, x, y, z, r\n";
    for (const auto &p : d_particlesListTypeAll) {
      auto xc = p->getXCenter();
      oss << p->d_zoneId << ", " << xc.d_x << ", " << xc.d_y << ", " << xc.d_z
          << ", " << p->d_geom_p->boundingRadius() << "\n";
    }
    oss.close();
  }
}

std::string model::DEMModel::ppTwoParticleTest() {

  bool continue_dt = false;
  auto check_dt = d_outputDeck_p->d_dtTestOut;
  if ((d_n % check_dt == 0) && (d_n >= check_dt))
    continue_dt = true;

  if (!continue_dt)
    return "";

  // get alias for particles
  const auto &p0 = this->d_particlesListTypeAll[0];
  const auto &p1 = this->d_particlesListTypeAll[1];

  // get penetration distance
  const auto &xc0 = p0->getXCenter();
  const auto &xc1 = p1->getXCenter();
  const double &r = p0->d_geom_p->boundingRadius();

  const auto &contact = d_cDeck_p->getContact(p0->d_zoneId, p1->d_zoneId);
  double r_e = r + contact.d_contactR;

  double pen_dist = xc1.dist(xc0) - r_e - r;
  double contact_area_radius = 0.;
  if (util::isLess(pen_dist, 0.))
    contact_area_radius =
        std::sqrt(std::pow(r_e, 2.) - std::pow(r_e + pen_dist, 2.));
  else if (util::isGreater(pen_dist, 0.)) {
    pen_dist = 0.;
    contact_area_radius = 0.;
  }

  // get max distance of second particle (i.e. the y-coord of center + radius)
  double max_dist = xc1.d_y + p1->d_geom_p->boundingRadius();

  // compute maximum y coordinate of particle 2
  double max_y_loc = p1->getXLocal(0).d_y;
  double max_y = 0.;
  for (size_t i = 0; i < p1->getNumNodes(); i++)
    if (util::isLess(max_y_loc, p1->getXLocal(i).d_y))
      max_y_loc = p1->getXLocal(i).d_y;

  if (util::isLess(max_y, max_y_loc))
    max_y = max_y_loc;

  setKeyData("pen_dist", pen_dist);
  setKeyData("contact_area_radius", contact_area_radius);
  setKeyData("max_y", max_y);
  setKeyData("max_dist", max_dist);
  setKeyData("max_y_loc", max_y_loc);


  return fmt::format("  Post-processing: max y = {:.6f} \n", max_y);
}

void model::DEMModel::checkStop() {

  if (d_outputDeck_p->d_outCriteria == "max_particle_dist" &&
      d_pDeck_p->d_testName == "two_particle") {

    // compute max distance between two particles
    // current center position
    const auto &xci = d_particlesListTypeAll[0]->getXCenter();
    const auto &xcj = d_particlesListTypeAll[1]->getXCenter();

    // check
    if (util::isGreater(xci.dist(xcj),
                        d_outputDeck_p->d_outCriteriaParams[0])) {

      if(d_ppFile.is_open())
        d_ppFile.close();
      exit(1);
    }
  }
  else if (d_outputDeck_p->d_outCriteria == "max_node_dist") {

    //    static int msg_printed = 0;
    //    if (msg_printed == 0) {
    //      std::cout << "Check = " << d_outputDeck_p->d_outCriteria
    //              << " is no longer supported. In future, this test will be implemented when function util::methods::maxLength() is defined." << std::endl;
    //      msg_printed = 1;
    //    }
    //exit(EXIT_FAILURE);
    auto max_pt_and_index = util::methods::maxLengthAndMaxLengthIndex(d_x);
    auto max_x = d_x[max_pt_and_index.second];

    // check
    if (util::isGreater(max_pt_and_index.first,
                        d_outputDeck_p->d_outCriteriaParams[0])) {

      // close open file
      if(d_ppFile.is_open())
        d_ppFile.close();

      log(fmt::format("{}: Terminating simulation as one of the failing"
                      " criteria is met. Point ({:.6f}, {:.6f}, {:.6f}) is at "
                      "distance {:.6f} "
                      "more than"
                      " allowed distance {:.6f}\n",
                      d_name, max_x.d_x, max_x.d_y, max_x.d_z, max_x.length(),
                      d_outputDeck_p->d_outCriteriaParams[0]));
      exit(1);
    }
  }
}

std::string model::DEMModel::ppCompressiveTest() {
  bool continue_dt = false;
  auto check_dt = d_outputDeck_p->d_dtTestOut;
  if ((d_n % check_dt == 0) && (d_n >= check_dt))
    continue_dt = true;

  if (!continue_dt)
    return "";

  // get wall
  auto w_id = d_pDeck_p->d_particleIdCompressiveTest;
  auto f_dir = d_pDeck_p->d_particleForceDirectionCompressiveTest - 1;
  const auto &wall = d_particlesListTypeAll[w_id];

  // find the penetration of the wall from it's original location
  auto dx = wall->getXLocal(0) - wall->getXRefLocal(0);
  double wall_penetration = dx[f_dir];

  // get the total reaction force on wall along the direction of loading
  double tot_reaction_force = 0.;
  for (size_t i = 0; i < wall->getNumNodes(); i++) {
    tot_reaction_force += wall->getFLocal(i)[f_dir] * wall->getVolLocal(i);
  }

  // open file and write
  bool use_static_file = true;
  if (use_static_file) {
    if (!d_ppFile.is_open()) {

      std::string tag_pp_file = d_outputDeck_p->d_tagPPFile.empty() ? "0" : d_outputDeck_p->d_tagPPFile;
      std::string filename = d_outputDeck_p->d_path + "pp_" +
                             d_pDeck_p->d_testName + "_" +
                             tag_pp_file + ".csv";
      d_ppFile.open(filename.c_str(), std::ofstream::out | std::ofstream::app);

      d_ppFile << "t, delta, force \n";
    }

    d_ppFile << fmt::format("%4.6e, %4.6e, %4.6e\n", d_time, wall_penetration,
            tot_reaction_force);
  }

  setKeyData("wall_penetration", wall_penetration);
  setKeyData("tot_reaction_force", tot_reaction_force);

  return fmt::format("  Post-processing: wall penetration = {:"
                     ".6f}, "
                     "reaction force = {:5.3e} \n",
                     wall_penetration, tot_reaction_force);
}