Main model class to handle demo two-particle test implementation. Goal is to show that it is easy to specialize the DEMModel class for different scenarios. More...

Inheritance diagram for twoparticle_demo::Model:

Collaboration diagram for twoparticle_demo::Model:

Public Member Functions
	Model (inp::Input *deck)
	Constructor.

void	run (inp::Input *deck) override
	Runs simulation by executing steps such as init(), integrate(), and close().

void	integrate () override
	Perform time step integration.

void	twoParticleTest ()
	Handles postprocessing for two-particle test.

void	twoParticleTestPenetrationDist ()
	Computes penetration distance of top particle into bottom particle.

void	twoParticleTestMaxShearStress ()
	Computes maximum shear stress and its location in particle.

Public Member Functions inherited from model::DEMModel
	DEMModel (inp::Input *deck, std::string modelName="DEMModel")
	Constructor.

void	log (std::ostringstream &oss, int priority=0, bool check_condition=true, int override_priority=-1, bool screen_out=false)
	Prints message if any of these two conditions are true.

void	log (const std::string &str, int priority=0, bool check_condition=true, int override_priority=-1, bool screen_out=false)
	Prints message if any of these two conditions are true.

virtual void	restart (inp::Input *deck)
	Restarts the simulation from previous state.

virtual void	init ()
	Initialize remaining data members.

virtual void	close ()
	Closure operations.

virtual void	integrateStep ()
	Performs one time step.

virtual void	integrateCD ()
	Perform time integration using central-difference scheme.

virtual void	integrateVerlet ()
	Perform time integration using velocity verlet scheme.

virtual void	computeForces ()
	Computes peridynamic forces and contact forces.

virtual void	computePeridynamicForces ()
	Computes peridynamic forces.

virtual void	computeExternalForces ()
	Computes external/boundary condition forces.

virtual void	computeExternalDisplacementBC ()
	Applies displacement boundary conditions.

virtual void	computeContactForces ()
	Computes contact forces.

virtual void	applyInitialCondition ()
	Applies initial condition.

virtual void	createParticles ()
	Creates particles in a given container.

virtual void	createParticlesFromFile (size_t z, std::shared_ptr< particle::RefParticle > ref_p)
	Creates particles in a Hexagonal arrangement.

virtual void	createParticleUsingParticleZoneGeomObject (size_t z, std::shared_ptr< particle::RefParticle > ref_p)
	Creates particles in a given container.

virtual double	createGeometryAtSite (const double &particle_radius, const double &particle_orient, const util::Point &site, const std::vector< double > &rep_geom_params, const std::shared_ptr< util::geometry::GeomObject > &rep_geom_p, std::shared_ptr< util::geometry::GeomObject > &p_geom)
	Creates geometrical object for a particle given particle radius, orientation, and site location.

virtual void	updateGeometryObjectsPostInit ()
	Update varioud geometry objects associated with container, particles, and reference particles.

virtual void	updateContactNeighborlist ()
	Update neighborlist for contact.

virtual void	updatePeridynamicNeighborlist ()
	Update neighborlist for peridynamics force.

virtual void	updateNeighborlistCombine ()
	Update neighborlist for contact and peridynamics force.

virtual void	setupContact ()
	Creates particles in a given container.

virtual void	setupQuadratureData ()
	Sets up quadrature data.

virtual bool	updateContactNeighborSearchParameters ()
	Update contact neighbor search parameters.

virtual void	output ()
	Output the snapshot of data at current time step.

std::string	ppTwoParticleTest ()
	Function that handles post-processing for two particle collision test and returns maximum vertical displacement of particle.

std::string	ppCompressiveTest ()
	Function that handles post-processing for compressive test of particulate media by rigid wall and returns wall penetration and reaction force on wall.

virtual void	checkStop ()
	Checks if simulation should be stopped due to abnormal state of system.

Public Member Functions inherited from model::ModelData
	ModelData (inp::Input *deck)
	Constructor.

const particle::BaseParticle *	getParticleFromAllList (size_t i) const
	Get pointer to base particle.

particle::BaseParticle *&	getParticleFromAllList (size_t i)
	Get pointer to base particle.

const particle::BaseParticle *	getParticleFromParticleList (size_t i) const
	Get pointer to particle (excluding wall)

particle::BaseParticle *&	getParticleFromParticleList (size_t i)

const particle::BaseParticle *	getParticleFromWallList (size_t i) const
	Get pointer to wall.

particle::BaseParticle *&	getParticleFromWallList (size_t i)

double	getDensity (size_t i)
	Get density of particle.

double	getHorizon (size_t i)
	Get horizon of particle.

size_t &	getPtId (size_t i)
	Get particle id given the location in particle list.

const size_t &	getPtId (size_t i) const
	Get particle id given the location in particle list.

void	setPtId (size_t i, const size_t &id)
	Set particle id given the location in particle list.

double	getKeyData (std::string key, bool issue_err=false)
	Get data for a key.

void	appendKeyData (std::string key, double data, bool issue_err=false)
	Append value to data associated with key.

void	setKeyData (std::string key, double data, bool issue_err=false)
	Set value to data associated with key.

util::Point &	getXRef (size_t i)
	Get reference coordinate of the node.

const util::Point &	getXRef (size_t i) const
	Get reference coordinate of the node.

void	setXRef (size_t i, const util::Point &x)
	Set reference coordinate of the node.

void	addXRef (size_t i, const util::Point &x)
	Add reference coordinate of the node.

void	setXRef (size_t i, int dof, double x)
	Set specific reference coordinate of the node.

void	addXRef (size_t i, int dof, double x)
	Add specific reference coordinate of the node.

util::Point &	getX (size_t i)
	Get current coordinate of the node.

const util::Point &	getX (size_t i) const
	Get current coordinate of the node.

void	setX (size_t i, const util::Point &x)
	Set current coordinate of the node.

void	addX (size_t i, const util::Point &x)
	Add current coordinate of the node.

void	setX (size_t i, int dof, double x)
	Set specific current coordinate of the node.

void	addX (size_t i, int dof, double x)
	Add to specific current coordinate of the node.

util::Point &	getU (size_t i)
	Get displacement of the node.

const util::Point &	getU (size_t i) const
	Get displacement of the node.

void	setU (size_t i, const util::Point &u)
	Set displacement of the node.

void	addU (size_t i, const util::Point &u)
	Add to displacement of the node.

void	setU (size_t i, int dof, double u)
	Set displacement of the node.

void	addU (size_t i, int dof, double u)
	Add to displacement of the node.

util::Point &	getV (size_t i)
	Get velocity of the node.

const util::Point &	getV (size_t i) const
	Get velocity of the node.

void	setV (size_t i, const util::Point &v)
	Set velocity of the node.

void	addV (size_t i, const util::Point &v)
	Add to velocity of the node.

void	setV (size_t i, int dof, double v)
	Set velocity of the node.

void	addV (size_t i, int dof, double v)
	Add to velocity of the node.

util::Point &	getF (size_t i)
	Get force of the node.

const util::Point &	getF (size_t i) const
	Get force of the node.

void	setF (size_t i, const util::Point &f)
	Set force of the node.

void	addF (size_t i, const util::Point &f)
	Add to force of the node.

void	setF (size_t i, int dof, double f)
	Set force of the node.

void	addF (size_t i, int dof, double f)
	Add to force of the node.

double &	getVol (size_t i)
	Get volume of the node.

const double &	getVol (size_t i) const
	Get volume of the node.

void	setVol (size_t i, const double &vol)
	Set volume of the node.

void	addVol (size_t i, const double &vol)
	Add to volume of the node.

uint8_t &	getFix (size_t i)
	Get fixity of the node.

const uint8_t &	getFix (size_t i) const
	Get fixity of the node.

void	setFix (size_t i, const unsigned int &dof, const bool &flag)
	Set fixity of the node.

double &	getMx (size_t i)
	Get weighted-volume (mx) of the node.

const double &	getMx (size_t i) const
	Get weighted-volume (mx) of the node.

void	setMx (size_t i, const double &mx)
	Set weighted-volume (mx) of the node.

void	addMx (size_t i, const double &mx)
	Add to weighted-volume (mx) of the node.

double &	getThetax (size_t i)
	Get volumetric deformation (thetax) of the node.

const double &	getThetax (size_t i) const
	Get volumetric deformation (thetax) of the node.

void	setThetax (size_t i, const double &thetax)
	Set volumetric deformation (thetax) of the node.

void	addThetax (size_t i, const double &thetax)
	Add to volumetric deformation (thetax) of the node.

Data Fields
double	d_penDist = 0.
	Penetration distance of top particle into bottom particle.

double	d_contactAreaRadius = 0.
	Contact area radius.

double	d_maxDist = 0.
	Maximum vertical distance of top particle in initial configuration.

double	d_maxStress = 0.
	Maximum stress.

double	d_maxStressLocRef = 0.
	Maximum stress location in reference configuration.

double	d_maxStressLocCur = 0.
	Maximum stress location in current configuration.

double	d_maxY = 0.
	Current maximum vertical distance of top particle.

double	d_contactAreaRadiusIdeal = 0.
	Ideal contact area radius from Hertz theory.

double	d_maxStressIdeal = 0.
	Ideal maximum stress.

double	d_maxStressLocRefIdeal = 0.
	Ideal maximum stress location location in reference configuration.

Data Fields inherited from model::DEMModel
std::string	d_name
	Name of the model for logging purposes (useful if other classes are built on top of this class)

Data Fields inherited from model::ModelData
size_t	d_n
	Current time step.

double	d_time
	Current time.

double	d_currentDt
	Current timestep.

size_t	d_infoN
	Print log step interval.

std::map< std::string, double >	d_dbgData
	Debug data.

std::ofstream	d_ppFile
	File stream to output information.

inp::Input *	d_input_p
	Pointer to Input object.

std::shared_ptr< inp::ModelDeck >	d_modelDeck_p
	Model deck.

std::shared_ptr< inp::RestartDeck >	d_restartDeck_p
	Restart deck.

std::shared_ptr< inp::OutputDeck >	d_outputDeck_p
	Output deck.

std::shared_ptr< inp::ParticleDeck >	d_pDeck_p
	Particle deck.

std::shared_ptr< inp::ContactDeck >	d_cDeck_p
	Contact deck.

bool	d_stop
	flag to stop the simulation midway

double	d_hMax
	Minimum mesh over all particles and walls.

double	d_hMin
	Maximum mesh over all particles and walls.

double	d_maxContactR
	Maximum contact radius between over pairs of particles and walls.

size_t	d_contNeighUpdateInterval
	Neighborlist update interval.

size_t	d_contNeighTimestepCounter
	Contact neighborlist time step counter.

double	d_contNeighSearchRadius
	Neighborlist contact search radius (multiple of d_maxContactR). This variable will be updated during simulation based on maximum velocity.

std::vector< std::shared_ptr< particle::RefParticle > >	d_referenceParticles
	Pointer to reference particle.

std::vector< particle::BaseParticle * >	d_particlesListTypeAll
	List of particles + walls.

std::vector< particle::BaseParticle * >	d_particlesListTypeParticle
	List of particles.

std::vector< particle::BaseParticle * >	d_particlesListTypeWall
	List of walls.

std::vector< inp::MatData >	d_particlesMatDataList
	List of particle material data. Only populated if needed for calculation of stress or other quantities.

std::vector< double >	d_maxVelocityParticlesListTypeAll
	Maximum velocity among all nodes in the particle for each particle.

double	d_maxVelocity
	Maximum velocity among all nodes.

std::vector< std::vector< size_t > >	d_zInfo
	Zone information of particles. For zone 0, d_zInfo[0] = [n1, n2] where n1 is the index at which the particle in this zone starts in d_particlesListTypeAll and n2 is the index + 1, where index is at which the particle in this zone ends.

std::unique_ptr< loading::ParticleULoading >	d_uLoading_p
	Pointer to displacement Loading object.

std::unique_ptr< loading::ParticleFLoading >	d_fLoading_p
	Pointer to force Loading object.

std::unique_ptr< geometry::Fracture >	d_fracture_p
	Fracture state of bonds.

std::unique_ptr< NSearch >	d_nsearch_p
	Pointer to nsearch.

std::vector< util::Point >	d_xRef
	reference positions of the nodes

std::vector< util::Point >	d_x
	Current positions of the nodes.

std::vector< util::Point >	d_u
	Displacement of the nodes.

std::vector< util::Point >	d_v
	Velocity of the nodes.

std::vector< double >	d_vMag
	Magnitude of velocity of the nodes.

std::vector< util::Point >	d_f
	Total force on the nodes.

std::vector< double >	d_vol
	Nodal volumes.

std::vector< size_t >	d_ptId
	Global node to particle id (walls are assigned id after last particle id)

std::vector< std::vector< size_t > >	d_neighC
	Neighbor data for contact forces.

std::vector< std::vector< size_t > >	d_neighPd
	Neighbor data for peridynamic forces.

std::vector< std::vector< float > >	d_neighPdSqdDist
	Square distance neighbor data for peridynamic forces.

std::vector< std::vector< std::vector< size_t > > >	d_neighWallNodes
	Neighbor data for contact between particle and walls.

std::vector< std::vector< std::vector< double > > >	d_neighWallNodesDistance
	Neighbor data (distance) for contact between particle and walls.

std::vector< std::vector< size_t > >	d_neighWallNodesCondensed
	Neighbor data for contact between particle and walls condensed into single vector for each particle.

std::vector< uint8_t >	d_fix
	Vector of fixity mask of each node.

std::vector< uint8_t >	d_forceFixity
	Vector of fixity mask of each node for force.

std::vector< double >	d_thetaX
	Dilation.

std::vector< double >	d_mX
	Weighted volume.

std::vector< size_t >	d_fPdCompNodes
	List of global nodes on which force (peridynamic/internal) is to be computed.

std::vector< size_t >	d_fContCompNodes
	List of global nodes on which force (contact) is to be computed.

std::vector< float >	d_Z
	Damage at nodes.

std::vector< float >	d_e
	Energy of the nodes.

std::vector< float >	d_w
	Work done on each of the nodes.

std::vector< float >	d_phi
	Damage function \( \phi \) at the nodes.

std::vector< float >	d_eF
	Fracture energy of the nodes.

std::vector< float >	d_eFB
	Bond-based fracture energy of the nodes.

std::vector< util::Point >	d_xQuadCur
	Current position of quadrature points.

std::vector< util::SymMatrix3 >	d_strain
	Strain in elements (values at quadrature points)

std::vector< util::SymMatrix3 >	d_stress
	Stress in elements (values at quadrature points)

float	d_te
	Total internal energy.

float	d_tw
	Total work done.

float	d_tk
	Total kinetic energy.

float	d_teF
	Total fracture energy.

float	d_teFB
	Total bond-based fracture energy.

Detailed Description

Main model class to handle demo two-particle test implementation. Goal is to show that it is easy to specialize the DEMModel class for different scenarios.

Definition at line 55 of file main.cpp.

Constructor & Destructor Documentation

◆ Model()

twoparticle_demo::Model::Model ( inp::Input * deck )

inlineexplicit

Constructor.

Parameters

deck	Input deck

Definition at line 64 of file main.cpp.

                                 : model::DEMModel(deck, "twoparticle_demo::Model") {
 
    if (d_ppFile.is_open())
      d_ppFile.close();
 
    std::string filename =
        d_outputDeck_p->d_path + "pp_" + d_outputDeck_p->d_tagPPFile + ".csv";
    d_ppFile = std::ofstream(filename, std::ios_base::out);
    d_ppFile << "t, delta, cont_area_r, s_loc, s_val, max_dist, "
                "cont_area_r_ideal, s_loc_ideal, s_val_ideal" << std::flush;
  }

References model::ModelData::d_outputDeck_p, and model::ModelData::d_ppFile.

Member Function Documentation

◆ integrate()

void twoparticle_demo::Model::integrate ( )

inlineoverridevirtual

Perform time step integration.

Reimplemented from model::DEMModel.

Definition at line 132 of file main.cpp.

                            {
    // 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();
    computeExternalForces();
 
    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);
 
      // NOTE: If there is need for different time-stepping scheme, one can
      // define another function similar to these two below
      auto t1 = steady_clock::now();
      integrateStep();
      double integrate_time =
              util::methods::timeDiff(t1, steady_clock::now());
 
      log(fmt::format("  Integration time (ms) = {}\n", integrate_time),
          2, d_n % d_infoN == 0, 3);
 
      if (d_pDeck_p->d_testName == "two_particle") {
        // NOTE: The purpose of this app 'twop' is to show that if we have
        // specific post-processing requirement from the outcome of
        // simulation -- eg in two-particle test, we may be interested in the
        // maximum of y-coordinate of top particle to measure the damping
        // effect or the maximum shear stress -- we wrap the base DEMModel
        // class as we did here and add a specific post-processing function
        // like 'ppTwoParticleTest()'.
        twoParticleTest();
      }
 
      // 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 (we may want to terminate the simulation early if the
      // results are garbage, or if some other criteria is met)
      // NOTE: this too can be specific to particular application
      checkStop();
    } // loop over time steps
  }

References model::DEMModel::applyInitialCondition(), model::DEMModel::checkStop(), model::DEMModel::computeExternalDisplacementBC(), model::DEMModel::computeExternalForces(), model::ModelData::d_infoN, model::ModelData::d_modelDeck_p, model::ModelData::d_n, model::DEMModel::d_name, model::ModelData::d_outputDeck_p, model::ModelData::d_pDeck_p, model::ModelData::d_time, model::DEMModel::integrateStep(), model::DEMModel::log(), model::DEMModel::output(), util::methods::timeDiff(), and twoParticleTest().

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ run()

void twoparticle_demo::Model::run ( inp::Input * deck )

inlineoverridevirtual

Runs simulation by executing steps such as init(), integrate(), and close().

Parameters

deck	Input deck

Reimplemented from model::DEMModel.

Definition at line 81 of file main.cpp.

                                    {
 
    log(d_name + ": Running TwoParticle_Demo app \n");
 
    // initialize data
    init();
 
    // check if init() successfully created quadrature data which we need for postprocessing
    {
      size_t totalQuadPoints = 0;
      for (auto &p: d_particlesListTypeAll) {
        const auto &particle_mesh_p = p->d_rp_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);
        }
 
        totalQuadPoints += particle_mesh_p->getNumElements() *
                           elem->getNumQuadPoints();
      }
 
      if (d_xQuadCur.size() != totalQuadPoints
            or d_strain.size() != totalQuadPoints
            or d_stress.size() != totalQuadPoints) {
        std::cerr << "Error: DEMModel::init() did not initialize quadrature data.\n";
        exit(EXIT_FAILURE);
      }
    }
 
    // integrate in time
    integrate();
 
    // close
    close();
  }

References model::DEMModel::close(), model::ModelData::d_modelDeck_p, model::DEMModel::d_name, model::ModelData::d_particlesListTypeAll, model::ModelData::d_strain, model::ModelData::d_stress, model::ModelData::d_xQuadCur, fe::BaseElem::getNumQuadPoints(), model::DEMModel::init(), integrate(), model::DEMModel::log(), util::vtk_type_line, util::vtk_type_quad, util::vtk_type_tetra, and util::vtk_type_triangle.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ twoParticleTest()

void twoparticle_demo::Model::twoParticleTest ( )

inline

Handles postprocessing for two-particle test.

Definition at line 193 of file main.cpp.

                         {
 
    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;
 
    // compute QoIs
    twoParticleTestPenetrationDist();
    twoParticleTestMaxShearStress();
 
    // log
    d_ppFile << d_time << ", " << -d_penDist << ", " << d_contactAreaRadius
             << ", " << d_maxStressLocRef << ", " << d_maxStress << ", "
             << d_maxDist << ", " << d_contactAreaRadiusIdeal << ", "
             << d_maxStressLocRefIdeal << ", " << d_maxStressIdeal
             << std::endl;
  }

References d_contactAreaRadius, d_contactAreaRadiusIdeal, d_maxDist, d_maxStress, d_maxStressIdeal, d_maxStressLocRef, d_maxStressLocRefIdeal, model::ModelData::d_n, model::ModelData::d_outputDeck_p, d_penDist, model::ModelData::d_ppFile, model::ModelData::d_time, twoParticleTestMaxShearStress(), and twoParticleTestPenetrationDist().

Referenced by integrate().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ twoParticleTestMaxShearStress()

void twoparticle_demo::Model::twoParticleTestMaxShearStress ( )

inline

Computes maximum shear stress and its location in particle.

Definition at line 281 of file main.cpp.

                                       {
 
    // compute maximum shear stress and where it occurs
    double max_stress_t = 0.;
    auto max_stress_loc_cur_t = util::Point();
    auto max_stress_loc_ref_t = util::Point();
 
    // 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()));
      }
    }
 
    // compute stress and strain
    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);
 
      double p_max_stress = 0.;
      auto p_max_stress_loc_cur = util::Point();
      auto p_max_stress_loc_ref = util::Point();
      fe::getMaxShearStressAndLoc(p->getMeshP().get(), d_xRef, d_u, d_stress,
                                  p_max_stress,
                                  p_max_stress_loc_ref,
                                  p_max_stress_loc_cur,
                                  p->d_globStart, p->d_globQuadStart, d_modelDeck_p->d_quadOrder);
 
      if (util::isGreater(p_max_stress, max_stress_t)) {
        max_stress_t = p_max_stress;
        auto p_center_node_id = p->d_globStart + p->d_rp_p->getCenterNodeId();
        max_stress_loc_ref_t = p_max_stress_loc_ref - d_xRef[p_center_node_id];
        max_stress_loc_cur_t = p_max_stress_loc_cur - d_x[p_center_node_id];
      }
    }
 
    d_maxStress = max_stress_t;
    d_maxStressLocRef = max_stress_loc_ref_t.length();
    d_maxStressLocCur = max_stress_loc_cur_t.length();
  }

References d_maxStress, d_maxStressLocCur, d_maxStressLocRef, model::ModelData::d_modelDeck_p, model::ModelData::d_particlesListTypeAll, model::ModelData::d_particlesMatDataList, model::ModelData::d_pDeck_p, model::ModelData::d_strain, model::ModelData::d_stress, model::ModelData::d_u, model::ModelData::d_x, model::ModelData::d_xQuadCur, model::ModelData::d_xRef, fe::getCurrentQuadPoints(), fe::getMaxShearStressAndLoc(), fe::getStrainStress(), and util::isGreater().

Referenced by twoParticleTest().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ twoParticleTestPenetrationDist()

void twoparticle_demo::Model::twoParticleTestPenetrationDist ( )

inline

Computes penetration distance of top particle into bottom particle.

Definition at line 218 of file main.cpp.

                                        {
 
    // 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;
 
    d_penDist = xc1.dist(xc0) - r_e - r;
    if (util::isLess(d_penDist, 0.))
      d_contactAreaRadius =
              std::sqrt(std::pow(r_e, 2.) - std::pow(r_e + d_penDist, 2.));
    else if (util::isGreater(d_penDist, 0.)) {
      d_penDist = 0.;
      d_contactAreaRadius = 0.;
    }
 
    // get max distance of second particle (i.e. the y-coord of center + radius)
    d_maxDist = xc1.d_y + p1->d_geom_p->boundingRadius();
 
    // compute maximum y coordinate of particle 2
    double max_y_loc = p1->getXLocal(0).d_y;
    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(d_maxY, max_y_loc))
      d_maxY = max_y_loc;
 
    log(fmt::format("max y: {} \n", d_maxY), 2, d_n % d_infoN == 0, 3);
 
    // compute ideal values
    static int contact_pp_ideal = -1;
    if (contact_pp_ideal == -1) {
      double mass = p1->getDensity() * M_PI * std::pow(r, 2.);
 
      const auto &mat_data =
              p1->getMaterial()->computeMaterialProperties(d_modelDeck_p->d_dim);
 
      d_contactAreaRadiusIdeal = 3. * mass * std::abs(d_pDeck_p->d_gravity[1]) *
                                 2. * r * (1. - std::pow(mat_data.d_nu, 2.)) /
                                 (4. * mat_data.d_E);
      d_contactAreaRadiusIdeal = std::pow(d_contactAreaRadiusIdeal, 1. / 3.);
 
      d_maxStressLocRefIdeal = r - 0.48 * d_contactAreaRadiusIdeal;
 
      d_maxStressIdeal =
              0.93 * mass * std::abs(d_pDeck_p->d_gravity[1]) /
              (2. * M_PI * d_contactAreaRadiusIdeal * d_contactAreaRadiusIdeal);
 
      contact_pp_ideal = 0;
    }
  }