1_membrane_case.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 "gravity.hpp"
#include "membrane_wall.hpp"
#include "model_linear_spring.hpp"
#include "shape_sphere.hpp"
#include "simulation.hpp"
#include "wall_motion_control.hpp"
#include <iostream>
#include <unordered_map>
 
using namespace netdem;
using namespace std;
 
int MembraneCase(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.01, 0.01, 0.01);
 
  // // solver
  // sim->dem_solver.contact_solver_factory.settings.solver_type =
  //     ContactSolverSettings::SolverType::sdf;
  // sim->dem_solver.contact_solver_factory.settings.sdf_potential_type = 0;
 
  // material model
  auto cnt_model_0 = LinearSpring(2.0e5, 1.0e5, 0.3, 0.0);
  cnt_model_0.label = "particle_&_wall";
  auto cnt_model_0_prt = sim->scene.InsertContactModel(&cnt_model_0);
 
  auto cnt_model_1 = LinearSpring(2.0e5, 1.0e5, 0.5, 0.0);
  cnt_model_1.label = "membrane";
  auto cnt_model_1_prt = sim->scene.InsertContactModel(&cnt_model_1);
 
  sim->scene.SetNumberOfMaterials(2);
  sim->scene.SetCollisionModel(0, 0, cnt_model_0_prt);
  sim->scene.SetCollisionModel(0, 1, cnt_model_1_prt);
  sim->scene.SetCollisionModel(1, 1, cnt_model_1_prt);
 
  // particle shape
  auto sphere = Sphere(0.004);
  sim->scene.InsertShape(&sphere);
 
  // pack of particles to be inserted
  auto particle_list = PackGenerator::GetGridPack(
      0.035, 0.035, 0.06, 0, 0, 0.12, 5, 5, 8, sim->scene.GetShapes());
  for (auto &p : particle_list) {
    p.SetDensity(2650e3);
    p.SetVelocity(0, 0, -0.6);
    p.damp_numerical = 0.7;
  }
 
  // bottom wall
  auto plane_bot = Plane(0, 0, -0.05, 0, 0, 1);
  plane_bot.SetExtent(0.1); // make sure the plane extent is inside the domain
  auto plane_bot_ptr = sim->scene.InsertShape(&plane_bot);
  auto wall_bot = Wall(plane_bot_ptr);
  auto wall_bot_ptr = sim->scene.InsertWall(wall_bot);
 
  // membrane
  auto membrane_wall = MembraneWall(0.025, 0.2, 0.0025);
  membrane_wall.neo_k = 6.94e6;
  membrane_wall.neo_mu = 5.21e6;
  membrane_wall.density = 500e3;
  membrane_wall.damp_coef = 0.7;
  membrane_wall.timestep = 1.0e-5;
  membrane_wall.enable_deformation = false; // no deformation during packing
  membrane_wall.dump_info = true;
  membrane_wall.Init(sim);
  membrane_wall.SetRootPath("tmp/out_tri_comp/");
  membrane_wall.SetSaveByCycles(100);
  auto modifier_ptr = sim->modifier_manager.Insert(&membrane_wall);
  auto membrane_wall_ptr = static_cast<MembraneWall *>(modifier_ptr);
  sim->modifier_manager.Enable(membrane_wall.label);
 
  for (auto &w : membrane_wall_ptr->wall_list) {
    w->material_type = 1;
  }
 
  // gravity
  sim->scene.gravity_coef[2] = -9.81;
  auto grav = Gravity();
  grav.Init(sim);
  sim->modifier_manager.Insert(&grav);
  sim->modifier_manager.Enable(grav.label);
 
  // out particle, wall and contact data
  auto data_dumper = DataDumper();
  data_dumper.Init(sim);
  data_dumper.SetRootPath("tmp/out_tri_comp/");
  data_dumper.SetSaveByCycles(100);
  data_dumper.dump_contact_info = true;
  data_dumper.dump_wall_info = true;
  sim->modifier_manager.Insert(&data_dumper);
  sim->modifier_manager.Enable(data_dumper.label);
 
  // insert 20 packs of particles
  sim->dem_solver.timestep = 1.0e-4;
  for (int i = 0; i < 20; i++) {
    sim->scene.InsertParticle(particle_list);
    sim->Run(0.1);
  }
 
  // rest for 1.0 s
  sim->Run(1.0);
 
  // remove extra particles
  sim->scene.gravity_coef[2] = 0.0;
  for (auto p_ptr : sim->scene.particle_list) {
    if (p_ptr->pos[2] > 0.05 - 0.5 * p_ptr->shape->GetSize()) {
      sim->scene.RemoveParticle(p_ptr);
    }
  }
 
  // add top wall
  auto plane_top = Plane(0, 0, 0.05, 0, 0, -1);
  plane_top.SetExtent(0.1);
  auto plane_top_ptr = sim->scene.InsertShape(&plane_top);
  auto wall_top = Wall(plane_top_ptr);
  auto wall_top_ptr = sim->scene.InsertWall(wall_top);
 
  // reset the membrane dimensions
  membrane_wall_ptr->SetDimensions(0.025, 0.10);
 
  // fix the bottom and top of the membrane wall
  for (int i = 0; i < membrane_wall_ptr->nodes.size(); i++) {
    if (abs(membrane_wall_ptr->nodes[i][2] + 0.5 * membrane_wall_ptr->height) <
            1.0e-4 ||
        abs(membrane_wall_ptr->nodes[i][2] - 0.5 * membrane_wall_ptr->height) <
            1.0e-4) {
      membrane_wall_ptr->SetBCNodalVelocity(i, 0, 0, 0, 1, 1, 1);
    }
  }
 
  // gradually increase confining pressure
  membrane_wall_ptr->enable_deformation = true;
  for (int i = 0; i < 10; i++) {
    membrane_wall_ptr->SetPressure(-1.0e4 * (i + 1));
    sim->Run(0.1);
  }
 
  // add deformation drived loading for top wall
  auto disp_control = WallMotionControl();
  disp_control.Init(sim);
  disp_control.SetWall({wall_top_ptr->id});
  disp_control.SetVelocity(0, 0, -0.001);
  sim->modifier_manager.Insert(&disp_control);
  sim->modifier_manager.Enable(disp_control.label);
 
  // press the top of membrane to move with the compression top wall
  for (int i = 0; i < membrane_wall_ptr->nodes.size(); i++) {
    if (abs(membrane_wall_ptr->nodes[i][2] - 0.5 * membrane_wall_ptr->height) <
        1.0e-4) {
      membrane_wall_ptr->SetBCNodalVelocity(i, 0, 0, -0.001, 1, 1, 1);
    }
  }
 
  // sim->Run(2.0);
  sim->Run(20.0);
 
  delete sim;
 
  MPI_Finalize();
  return 0;
}