91_triaixal_comp.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 "gen_wall_box_plate.hpp"
#include "gravity.hpp"
#include "model_volume_based.hpp"
#include "shape_trimesh.hpp"
#include "simulation.hpp"
#include "solver_ann_pp.hpp"
#include "solver_ann_pplane.hpp"
#include "wall_motion_control.hpp"
#include "wall_servo_control.hpp"
#include <iostream>
#include <unordered_map>
using namespace netdem;
using namespace std;
 
void TriaxialComp() {
  // new simulation
  Simulation *sim = new Simulation();
 
  // domain settings
  sim->domain_manager.SetBound(-0.03, -0.03, -0.05, 0.03, 0.03, 0.13);
  auto domain = sim->domain_manager.GetSelfSubDomain();
  domain->cell_manager.SetSpacing(0.02, 0.02, 0.02);
 
  // define contact model
  VolumeBased cnt_model = VolumeBased(1.0e16, 1.0e5, 0.5, 0.0);
  sim->scene.InsertContactModel(&cnt_model);
  sim->scene.SetNumberOfMaterials(1);
  sim->scene.SetCollisionModel(0, 0, cnt_model.label);
 
  // shape template
  TriMesh tri_mesh;
  tri_mesh.InitFromSTL("data/particle_template.stl");
  tri_mesh.Decimate(200);
  tri_mesh.AlignAxes();
  tri_mesh.SetSize(0.003);
  sim->scene.InsertShape(&tri_mesh);
 
  // get a pack of particles
  VecXT<Particle> particle_list = PackGenerator::GetGridPack(
      0.04, 0.04, 0.04, 0, 0, 0.10, 8, 8, 8, sim->scene.GetShapes());
  for (auto &p : particle_list) {
    p.SetDensity(265000);
    p.damp_numerical = 0.7;
  }
 
  // get a box of wall with plane shapes, and insert the wall and shapes into
  // the scene
  WallBoxPlane wall_box(0.04, 0.04, 0.16, 0, 0, 0.04);
  wall_box.ImportToScene(&(sim->scene));
 
  // manipulate the solver settings
  string root_dir = "local/examples/netdem/ann_models/";
  SolverANNPP cnt_solver_pp;
  cnt_solver_pp.Init(root_dir + "trimesh_trimesh/ann_classifier.xml",
                     root_dir + "trimesh_trimesh/ann_regressor.xml");
  SolverANNPPlane cnt_solver_pw;
  cnt_solver_pw.Init(root_dir + "trimesh_plane/ann_classifier.xml",
                     root_dir + "trimesh_plane/ann_regressor.xml");
 
  int solver_id_pp =
      sim->dem_solver.contact_solver_factory.InsertSolver(&cnt_solver_pp);
  int solver_id_pw =
      sim->dem_solver.contact_solver_factory.InsertSolver(&cnt_solver_pw);
 
  auto shape_list = sim->scene.GetShapes();
  for (auto shape_ptr_1 : shape_list) {
    for (auto shape_ptr_2 : shape_list) {
      if (shape_ptr_1->shape_type == Shape::Type::trimesh &&
          shape_ptr_2->shape_type == Shape::Type::trimesh) {
        shape_ptr_1->use_customized_solver = true;
        shape_ptr_2->use_customized_solver = true;
        sim->dem_solver.contact_solver_factory.CustomizeSolverPP(
            shape_ptr_1->id, shape_ptr_2->id, solver_id_pp);
      } else if (shape_ptr_1->shape_type == Shape::Type::trimesh &&
                 shape_ptr_2->shape_type == Shape::Type::plane) {
        shape_ptr_1->use_customized_solver = true;
        shape_ptr_2->use_customized_solver = true;
        sim->dem_solver.contact_solver_factory.CustomizeSolverPW(
            shape_ptr_1->id, shape_ptr_2->id, solver_id_pw);
      }
    }
  }
 
  // add a gravity modifier
  Gravity grav;
  grav.Init(sim);
  sim->modifier_manager.Insert(&grav);
  sim->modifier_manager.Enable(grav.label);
 
  // add a data dumper modifier
  DataDumper data_dumper;
  data_dumper.Init(sim);
  data_dumper.SetRootPath("tmp/netdem/triaxial_comp/");
  data_dumper.SetSaveByCycles(1000);
  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);
 
  // random packing
  sim->dem_solver.timestep = 1e-4;
  for (int i = 0; i < 10; i++) {
    sim->scene.InsertParticle(particle_list);
    sim->Run(0.1);
  }
  sim->Run(0.5);
 
  // servo control for top wall, move the top wall right above the packing
  sim->scene.wall_list[5]->SetPosition(0, 0, -0.08);
  WallServoControl servo_top(1.0e5, 1.0);
  servo_top.label = "servo_top";
  servo_top.enable_auto_area = true;
  servo_top.Init(sim);
  servo_top.target_pressure = 1e3;
  servo_top.study_rate = 0.1;
  servo_top.tol = 0.95;
  servo_top.vel_max = 1.0;
  servo_top.SetWall({5});
  auto tmp = sim->modifier_manager.Insert(&servo_top);
  sim->modifier_manager.Enable(servo_top.label);
  auto servo_top_ptr = static_cast<WallServoControl *>(tmp);
  do {
    sim->Run(0.1);
  } while (!servo_top_ptr->achieved);
 
  // servo control for lateral walls, isotropic compression
  servo_top_ptr->target_pressure = 1e5;
  servo_top.study_rate = 0.5;
  servo_top.tol = 0.05;
 
  WallServoControl servo_lateral(1.0e5, 1.0);
  servo_lateral.label = "servo_lateral";
  servo_lateral.enable_auto_area = true;
  servo_lateral.Init(sim);
  servo_lateral.target_pressure = 1e5;
  servo_lateral.vel_max = 1.0;
  servo_lateral.SetWall({0, 1, 2, 3});
  tmp = sim->modifier_manager.Insert(&servo_lateral);
  sim->modifier_manager.Enable(servo_lateral.label);
  auto servo_lateral_ptr = static_cast<WallServoControl *>(tmp);
 
  // keep loading until target pressure is achieved
  do {
    sim->Run(0.1);
  } while (!servo_top_ptr->achieved || !servo_lateral.achieved);
  sim->Run(0.1);
 
  // disable servo control of the top wall
  sim->modifier_manager.Disable(servo_top.label);
  servo_lateral_ptr->enable_warning = true;
 
  // add deformation drived loading for top wall
  WallMotionControl disp_control;
  disp_control.Init(sim);
  disp_control.SetWall({5});
  disp_control.SetVelocity(0, 0, -0.001);
  sim->modifier_manager.Insert(&disp_control);
  sim->modifier_manager.Enable(disp_control.label);
 
  sim->Run(20);
 
  delete sim;
}