91_packing_ellipsoid.cpp

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

#include "data_dumper.hpp"
#include "gen_pack.hpp"
#include "gen_wall_box_plane.hpp"
#include "gen_wall_box_plate.hpp"
#include "gravity.hpp"
#include "model_linear_spring.hpp"
#include "model_volume_based.hpp"
#include "shape_cylinder.hpp"
#include "shape_ellipsoid.hpp"
#include "shape_level_set.hpp"
#include "shape_poly_super_ellipsoid.hpp"
#include "shape_poly_super_quadrics.hpp"
#include "shape_sphere.hpp"
#include "shape_spherical_harmonics.hpp"
#include "shape_trimesh.hpp"
#include "simulation.hpp"
#include <filesystem>
#include <iostream>
#include <string>
#include <unordered_map>
using namespace netdem;
using namespace std;
 
void PackingEllipsoid_Sub(double eta, double xi, string root_dir) {
  Simulation *sim = new Simulation();
 
  sim->domain_manager.SetBound(-0.6, -0.6, -0.6, 1.6, 1.6, 1.6);
  sim->domain_manager.SetCellSpacing(0.3, 0.3, 0.3);
 
  sim->dem_solver.contact_solver_factory.settings.solver_type =
      ContactSolverSettings::SolverType::sdf;
  sim->dem_solver.contact_solver_factory.settings.sdf_potential_type = 1;
  LinearSpring contact_model = LinearSpring(2.0e3, 1.0e4, 0.0, 0.0);
  sim->scene.InsertContactModel(&contact_model);
  sim->scene.SetNumberOfMaterials(1);
  sim->scene.SetCollisionModel(0, 0, contact_model.label);
 
  // note eta is acctually 1.0/eta. Zhao etal 2011 vs. Zhao 2017
  PolySuperEllipsoid ellipsoid(eta, eta, eta, eta, 1.0, 1.0, xi, xi);
  ellipsoid.SetSize(0.05);
  sim->scene.InsertShape(&ellipsoid);
 
  VecXT<Particle> particle_list = PackGenerator::GetGridPack(
      1, 1, 0.2, 0, 0, 1.4, 10, 10, 2, sim->scene.GetShapes());
  for (auto &p : particle_list) {
    p.SetVelocity(0, 0, -2.0);
    p.damp_numerical = 0.7;
  }
 
  WallBoxPlane wall_box(1, 1, 2, 0, 0, 0.5);
  wall_box.ImportToScene(&(sim->scene));
 
  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(root_dir);
  data_dumper.SetSaveByCycles(100);
  data_dumper.SaveShapeInfoAsSTL();
  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);
 
  for (int i = 0; i < 40; i++) {
    sim->scene.InsertParticle(particle_list);
    sim->Run(0.1);
  }
  sim->Run(1.0);
 
  delete sim;
}
 
int PackingEllipsoid(int argc, char **argv) {
  string root_dir;
 
  for (int i = 0; i < 5; i++) {
    double eta = 0.4 + 0.1 * i;
    double xi = 1.0;
    root_dir = "tmp/examples/sdf_dem/packing/eta_v1_" + to_string(i) + "/";
    filesystem::create_directories(root_dir);
    PackingEllipsoid_Sub(eta, xi, root_dir);
  }
 
  for (int i = 0; i < 5; i++) {
    double eta = 1.0 / (0.4 + 0.1 * i);
    double xi = 1.0;
    root_dir = "tmp/examples/sdf_dem/packing/eta_v2_" + to_string(i) + "/";
    filesystem::create_directories(root_dir);
    PackingEllipsoid_Sub(eta, xi, root_dir);
  }
 
  for (int i = 0; i < 10; i++) {
    double eta = 1.0;
    double xi = 0.4 + 0.1 * i;
    root_dir = "tmp/examples/sdf_dem/packing/zeta_" + to_string(i) + "/";
    filesystem::create_directories(root_dir);
    PackingEllipsoid_Sub(eta, xi, root_dir);
  }
 
  return 0;
}