81_voronoi_packing.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_polybezier.hpp"
#include "shape_sphere.hpp"
#include "shape_trimesh.hpp"
#include "simulation.hpp"
#include "solver_gjk_pp.hpp"
#include <filesystem>
#include <iostream>
#include <omp.h>
#include <unordered_map>
 
using namespace netdem;
using namespace std;
 
STLModel GetKernelPolyhedron() {
  int num_verts{30 / 2 + 2}; // V = E-F+2
 
  UniformDistribution uniform_dist(0, Math::PI);
  auto theta_list = uniform_dist.Get(num_verts);
  auto phi_list = uniform_dist.Get(num_verts);
 
  VecXT<Vec3d> vertices;
  VecXT<Vec3i> facets;
  vertices.resize(num_verts);
  for (int i = 0; i < num_verts; i++) {
    vertices[i][0] = sin(theta_list[i]) * cos(phi_list[i] * 2);
    vertices[i][1] = sin(theta_list[i]) * sin(phi_list[i] * 2);
    vertices[i][2] = cos(theta_list[i]);
  }
 
  IGLWrapper::ConvexHull(vertices, &vertices, &facets);
  STLModel kernel_stl(vertices, facets);
 
  Vec3d origin{0, 0, 0};
  if (!kernel_stl.CheckEnclose(origin)) {
    return GetKernelPolyhedron();
  } else {
    return kernel_stl;
  }
}
 
void PrepareTemplates() {
  std::filesystem::create_directories(
      "local/polybezier_model/voronoi_packing/kernel_stl/");
 
  int num_p{0};
  for (int i = 0; i < 1000 * 1000; i++) {
    auto kernel_stl = GetKernelPolyhedron();
 
    Polybezier polybezier;
    polybezier.InitFromKernelSTL(kernel_stl);
 
    auto [bound_aabb_min, bound_aabb_max] = polybezier.GetBoundAABB();
    auto bound_len = bound_aabb_max - bound_aabb_min;
    auto len_max = max(bound_len[0], max(bound_len[1], bound_len[2]));
    auto len_min = min(bound_len[0], min(bound_len[1], bound_len[2]));
    auto len_mid =
        bound_len[0] + bound_len[1] + bound_len[2] - len_max - len_min;
    auto flatness = len_min / len_mid;
    auto elongation = len_mid / len_max;
 
    if (flatness > 0.7 && flatness < 0.9 && elongation > 0.6 &&
        elongation < 0.8) {
      num_p++;
      kernel_stl.SaveAsSTL(
          "local/polybezier_model/voronoi_packing/kernel_stl/" +
          to_string(num_p - 1) + ".stl");
    }
 
    if (num_p > 999) {
      break;
    }
  }
}
 
void PolybezierCase() {
  Simulation *sim = new Simulation();
 
  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(2.0e6, 1.0e6, 0.5, 0.0);
  sim->scene.InsertContactModel(&cnt_model);
  sim->scene.SetNumberOfMaterials(1);
  sim->scene.SetCollisionModel(0, 0, cnt_model.label);
 
  // voronoi cells
  STLModel stl_model;
  stl_model.InitFromSTL("data/box.stl");
  auto [vt_nodes, vt_cells, vt_seeds] =
      Voronoi::Solve(stl_model, 1000, 1000, 1.0e-3);
  Voronoi::SaveAsVTK("local/polybezier_model/voronoi_packing/voronoi.vtk",
                     vt_nodes, vt_cells, vt_seeds);
 
  for (int i = 0; i < 1000; i++) {
    STLModel kernel_stl;
    kernel_stl.InitFromSTL(
        "local/polybezier_model/voronoi_packing/kernel_stl/" + to_string(i) +
        ".stl");
 
    Polybezier polybezier;
    polybezier.InitFromKernelSTL(kernel_stl);
    polybezier.SetSize(0.08);
 
    auto shape_prt = sim->scene.InsertShape(&polybezier);
 
    Particle p;
    p.SetShape(shape_prt);
    p.SetPosition(vt_seeds[i][0], vt_seeds[i][1], vt_seeds[i][2]);
    sim->scene.InsertParticle(p);
  }
 
  // container
  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("local/polybezier_model/voronoi_packing/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);
 
  sim->Run(2.0);
 
  auto cnt_solver = sim->dem_solver.contact_solver_factory.GetCollisionSolver(
      sim->scene.particle_list[0], sim->scene.particle_list[1]);
  auto gjk_solver = static_cast<SolverGJKPP *>(cnt_solver);
  omp_set_num_threads(1);
  sim->dem_solver.Cycle(1);
 
  delete sim;
 
  MPI_Finalize();
}
 
void TrimeshCase() {
  Simulation *sim = new Simulation();
 
  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(2.0e6, 1.0e6, 0.5, 0.0);
  sim->scene.InsertContactModel(&cnt_model);
  sim->scene.SetNumberOfMaterials(1);
  sim->scene.SetCollisionModel(0, 0, cnt_model.label);
 
  // voronoi cells
  STLModel stl_model;
  stl_model.InitFromSTL("data/box.stl");
  auto [vt_nodes, vt_cells, vt_seeds] =
      Voronoi::Solve(stl_model, 1000, 1000, 1.0e-3);
  Voronoi::SaveAsVTK("local/polybezier_model/voronoi_packing/voronoi.vtk",
                     vt_nodes, vt_cells, vt_seeds);
 
  for (int i = 0; i < 1000; i++) {
    STLModel kernel_stl;
    kernel_stl.InitFromSTL(
        "local/polybezier_model/voronoi_packing/kernel_stl/" + to_string(i) +
        ".stl");
 
    Polybezier polybezier;
    polybezier.InitFromKernelSTL(kernel_stl);
    polybezier.SetSize(0.08);
 
    TriMesh trimesh;
    // trimesh.use_linked_patches = true;
    trimesh.InitFromSTL(polybezier.GetSTLModel());
    trimesh.MakeConvex();
    trimesh.AlignAxes();
    auto shape_prt = sim->scene.InsertShape(&trimesh);
 
    Particle p;
    p.SetShape(shape_prt);
    p.SetPosition(vt_seeds[i][0], vt_seeds[i][1], vt_seeds[i][2]);
    p.damp_numerical = 0.7;
    sim->scene.InsertParticle(p);
  }
 
  // container
  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("local/polybezier_model/voronoi_packing/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);
 
  sim->Run(2.0);
 
  for (auto shape_ptr : sim->scene.GetShapes()) {
    if (shape_ptr->shape_type == Shape::trimesh_convex) {
      auto trimesh_ptr = static_cast<TriMesh *>(shape_ptr);
      // trimesh_ptr->use_linked_patches = true;
      // trimesh_ptr->Init();
      trimesh_ptr->enable_debugging = true;
    }
  }
  auto cnt_solver = sim->dem_solver.contact_solver_factory.GetCollisionSolver(
      sim->scene.particle_list[0], sim->scene.particle_list[1]);
  auto gjk_solver = static_cast<SolverGJKPP *>(cnt_solver);
  omp_set_num_threads(1);
  sim->dem_solver.Cycle(1);
 
  delete sim;
 
  MPI_Finalize();
}
 
void MeshTemplates() {
  std::filesystem::create_directories(
      "local/polybezier_model/voronoi_packing/mesh_stl/");
 
  for (int i = 0; i < 1000; i++) {
    STLModel kernel_stl;
    kernel_stl.InitFromSTL(
        "local/polybezier_model/voronoi_packing/kernel_stl/" + to_string(i) +
        ".stl");
 
    Polybezier polybezier;
    polybezier.InitFromKernelSTL(kernel_stl);
    polybezier.SetSize(0.08);
 
    polybezier.GetSTLModel(2000).SaveAsSTL(
        "local/polybezier_model/voronoi_packing/mesh_stl/" + to_string(i) +
        ".stl");
  }
}
 
void VoronoiPacking(int argc, char **argv) {
  if (argc <= 2) {
    cout << "please specify the id of the task. \n"
         << "e.g.: netdem_example_polybezier_model 81 i" << endl;
    cout << ">> 1: PrepareTemplates" << endl;
    cout << ">> 2: PolybezierCase" << endl;
    cout << ">> 3: TrimeshCase" << endl;
    cout << ">> 4: MeshTemplates" << endl;
 
    return;
  }
 
  int id = atof(argv[2]);
  switch (id) {
  case 1:
    PrepareTemplates();
    break;
  case 2:
    PolybezierCase();
    break;
  case 3:
    TrimeshCase();
    break;
  case 4:
    MeshTemplates();
    break;
  default:
    break;
  }
 
  return;
}