02_energy_test.cpp

This is an example of how to use the netdem library.

#include "contact_solver_factory.hpp"
#include "data_dumper.hpp"
#include "gen_pack.hpp"
#include "gen_wall_box_plane.hpp"
#include "gravity.hpp"
#include "model_linear_spring.hpp"
#include "particle_energy_evaluator.hpp"
#include "shape_spherical_harmonics.hpp"
#include "simulation.hpp"
#include "solver_gjk_pp.hpp"
#include "solver_gjk_pw.hpp"
#include "solver_sdf_pp.hpp"
#include "solver_sdf_pw.hpp"
#include <filesystem>
#include <iostream>
#include <unordered_map>
using namespace netdem;
using namespace std;
 
void SaveDataset(string filename, VecXT<VecXT<double>> &data);
 
void EnergyTest(int potential_case, int shape_case, double timestep) {
  Simulation *sim = new Simulation();
  sim->enable_logging = false;
 
  sim->domain_manager.SetBound(-0.6, -0.6, -0.6, 0.6, 0.6, 0.6);
  sim->domain_manager.SetCellSpacing(0.3, 0.3, 0.3);
 
  LinearSpring cnt_model = LinearSpring(1.0e7, 1.0e6, 0.3, 0.0);
  sim->scene.InsertContactModel(&cnt_model);
  sim->scene.SetNumberOfMaterials(1);
  sim->scene.SetCollisionModel(0, 0, cnt_model.label);
 
  SphericalHarmonics sh;
  switch (shape_case) {
  case 0:
    // sphere case
    sh.a_nm[0] = 0.1 * sqrt(Math::PI);
    sh.UpdateShapeProperties();
    break;
  case 1:
    // trimesh case
    sh.InitFromSTL("data/particle_template.stl");
    sh.SetSize(0.1);
    break;
  default:
    // sphere case
    cout << "please specify 0 for sphere, 1 for trimesh" << endl;
    abort();
    break;
  }
  sim->scene.InsertShape(&sh);
 
  VecXT<Particle> particle_list = PackGenerator::GetGridPack(
      1, 1, 1, 0, 0, 0, 7, 7, 7, sim->scene.GetShapes());
  sim->scene.InsertParticle(particle_list);
 
  WallBoxPlane wall_box(1, 1, 1, 0, 0, 0);
  wall_box.ImportToScene(&(sim->scene));
 
  // manipulate the solver settings
  sim->dem_solver.contact_solver_factory.settings.sdf_potential_type =
      potential_case;
  sim->dem_solver.contact_solver_factory.settings.sdf_use_equivalent_stiffness =
      false;
 
  // modifiers
  Gravity grav;
  grav.Init(sim);
  sim->modifier_manager.Insert(&grav);
  sim->modifier_manager.Enable(grav.label);
 
  DataDumper data_dumper;
  data_dumper.Init(sim);
  data_dumper.SetRootPath("local/potential_models/energy_test/out/");
  data_dumper.SetSaveByCycles(100);
  data_dumper.dump_mesh = true;
  data_dumper.dump_wall_info = true;
  data_dumper.dump_contact_info = true;
  sim->modifier_manager.Insert(&data_dumper);
  sim->modifier_manager.Enable(data_dumper.label);
  data_dumper.SaveShapeInfoAsSTL();
 
  ParticleEnergyEvaluator pe_evaluator;
  pe_evaluator.Init(sim);
  pe_evaluator.SetParticleFromScene();
  auto ev = sim->modifier_manager.Insert(&pe_evaluator);
  auto pe_evaluator_ptr = static_cast<ParticleEnergyEvaluator *>(ev);
  sim->modifier_manager.Enable(pe_evaluator.label);
 
  sim->dem_solver.timestep = timestep;
  VecXT<VecXT<double>> energy_data;
  for (int i = 0; i < 300; i++) {
    sim->Run(0.01);
 
    // two particle case
    // energy_data.push_back(
    //     {pe_evaluator_ptr->particle_energy_list[0].total,
    //      pe_evaluator_ptr->particle_energy_list[0].translational,
    //      pe_evaluator_ptr->particle_energy_list[0].rotational,
    //      pe_evaluator_ptr->particle_energy_list[0].gravitational,
    //      pe_evaluator_ptr->particle_energy_list[1].total,
    //      pe_evaluator_ptr->particle_energy_list[1].translational,
    //      pe_evaluator_ptr->particle_energy_list[1].rotational,
    //      pe_evaluator_ptr->particle_energy_list[1].gravitational});
 
    // multiple particle case
    auto tmp_energy = pe_evaluator_ptr->GetEnergy();
    energy_data.push_back({tmp_energy.total, tmp_energy.translational,
                           tmp_energy.rotational, tmp_energy.gravitational});
  }
 
  string root_dir = "local/potential_models/energy_test/";
  filesystem::create_directories(root_dir);
  switch (shape_case) {
  case 0:
    switch (potential_case) {
    case SolverSDFPP::PotentialType::linear:
      SaveDataset(root_dir + "sphere_linear.txt", energy_data);
      break;
    case SolverSDFPP::PotentialType::hertz:
      SaveDataset(root_dir + "sphere_hertz.txt", energy_data);
      break;
    default:
      SaveDataset(root_dir + "sphere_hertz.txt", energy_data);
      break;
    }
    break;
  case 1:
    switch (potential_case) {
    case SolverSDFPP::PotentialType::linear:
      SaveDataset(root_dir + "trimesh_linear.txt", energy_data);
      break;
    case SolverSDFPP::PotentialType::hertz:
      SaveDataset(root_dir + "trimesh_hertz.txt", energy_data);
      break;
    default:
      SaveDataset(root_dir + "trimesh_hertz.txt", energy_data);
      break;
    }
    break;
  default:
    break;
  }
 
  sh.SaveAsSTL(root_dir + "shape.stl");
 
  delete sim;
}