1_gjk_solver.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 "shape_cylinder.hpp"
#include "shape_ellipsoid.hpp"
#include "shape_poly_super_ellipsoid.hpp"
#include "shape_poly_super_quadrics.hpp"
#include "shape_sphere.hpp"
#include "shape_trimesh.hpp"
#include "simulation.hpp"
#include <iostream>
#include <unordered_map>
using namespace netdem;
using namespace std;
 
int GJKSolver(int argc, char **argv) {
  int self_rank, total_rank;
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &self_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &total_rank);
 
  Simulation *sim = new Simulation();
 
  sim->domain_manager.SetBound(-0.6, -0.6, -0.6, 0.6, 0.6, 0.6);
  sim->domain_manager.SetDecomposition(total_rank, 1, 1);
  sim->domain_manager.SetCellSpacing(0.3, 0.3, 0.3);
 
  LinearSpring cnt_model = LinearSpring(2.0e6, 1.0e6, 0.5, 0.01);
  sim->scene.InsertContactModel(&cnt_model);
  sim->scene.SetNumberOfMaterials(1);
  sim->scene.SetCollisionModel(0, 0, cnt_model.label);
 
  int shape_case;
  if (argc == 2) {
    shape_case = -1;
  } else {
    shape_case = atof(argv[2]);
  }
 
  switch (shape_case) {
  case 0: {
    Ellipsoid ellipsoid(1.0, 1.0, 0.5);
    ellipsoid.SetSize(0.1);
    sim->scene.InsertShape(&ellipsoid);
 
    ellipsoid = Ellipsoid(1.0, 1.0, 2.0);
    ellipsoid.SetSize(0.1);
    sim->scene.InsertShape(&ellipsoid);
  } break;
  case 1: {
    Cylinder cylinder(0.2, 1.0);
    cylinder.SetSize(0.1);
    sim->scene.InsertShape(&cylinder);
 
    cylinder = Cylinder(1.0, 0.5);
    cylinder.SetSize(0.1);
    sim->scene.InsertShape(&cylinder);
  } break;
  case 2: {
    PolySuperEllipsoid ellipsoid(0.5, 1.0, 2.5, 0.5, 1.7, 0.5, 1.0, 1.5);
    ellipsoid.SetSize(0.1);
    sim->scene.InsertShape(&ellipsoid);
  } break;
  case 3: {
    PolySuperQuadrics quadrics(0.5, 1.0, 2.5, 0.5, 1.7, 0.5, 0.5, 1.0, 1.5, 1.0,
                               1.2, 0.8);
    quadrics.SetSize(0.1);
    sim->scene.InsertShape(&quadrics);
  } break;
  case 4: {
    TriMesh tri_mesh;
    tri_mesh.InitFromSTL("data/particle_template.stl");
    tri_mesh.Decimate(200);
    tri_mesh.MakeConvex();
    tri_mesh.AlignAxes();
    tri_mesh.SetSize(0.1);
    sim->scene.InsertShape(&tri_mesh);
  } break;
  default: {
    TriMesh tri_mesh;
    tri_mesh.InitFromSTL("data/particle_template.stl");
    tri_mesh.Decimate(200);
    tri_mesh.MakeConvex();
    tri_mesh.AlignAxes();
    tri_mesh.SetSize(0.1);
    sim->scene.InsertShape(&tri_mesh);
  } break;
  }
 
  VecXT<Particle> particle_list = PackGenerator::GetGridPack(
      1, 1, 0.2, 0, 0, 0.4, 5, 5, 1, sim->scene.GetShapes());
  for (auto &p : particle_list) {
    p.SetVelocity(0, 0, -2.0);
    p.damp_numerical = 0.0;
    p.damp_viscous = 0.0;
  }
 
  WallBoxPlane wall_box(1, 1, 1, 0, 0, 0);
  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("tmp/out/");
  data_dumper.SetSaveByCycles(100);
  data_dumper.SaveShapeInfoAsSTL();
  data_dumper.dump_mesh = true;
  sim->modifier_manager.Insert(&data_dumper);
  sim->modifier_manager.Enable(data_dumper.label);
 
  for (int i = 0; i < 20; i++) {
    sim->scene.InsertParticle(particle_list);
    sim->Run(0.1);
  }
  sim->Run(2.0);
 
  delete sim;
 
  MPI_Finalize();
  return 0;
}