landslide_main.cpp

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

#include "cgal_wrapper.hpp"
#include "cork_wrapper.hpp"
#include "data_dumper.hpp"
#include "gen_pack.hpp"
#include "gen_wall_box_plane.hpp"
#include "gravity.hpp"
#include "model_linear_spring.hpp"
#include "shape_sphere.hpp"
#include "shape_triangle.hpp"
#include "simulation.hpp"
#include "stl_model.hpp"
#include "utils_io.hpp"
#include <iostream>
#include <unordered_map>
 
using namespace netdem;
using namespace std;
 
// reconstruct stl mesh based DEM (digital elevation model) data
STLModel GetSTLFromDEM(string const &dem_file, double spacing) {
  auto dem_data = IO::ImportData(dem_file, 6);
 
  int num_rows = dem_data.size();
  int num_cols = dem_data[0].size();
 
  double length = (num_cols - 1) * spacing;
  double width = (num_rows - 1) * spacing;
 
  STLModel stl_model;
 
  stl_model.vertices.reserve(num_rows * num_cols);
  for (int i = 0; i < num_rows; i++) {
    for (int j = 0; j < num_cols; j++) {
      stl_model.vertices.push_back({-0.5 * length + j * spacing,
                                    0.5 * width - i * spacing, dem_data[i][j]});
    }
  }
 
  for (int i = 0; i < num_rows - 1; i++) {
    for (int j = 0; j < num_cols - 1; j++) {
      auto vid_00 = i * num_cols + j;
      auto vid_01 = i * num_cols + j + 1;
      auto vid_10 = (i + 1) * num_cols + j;
      auto vid_11 = (i + 1) * num_cols + j + 1;
 
      // the elevation -9999 indicate points out of computational domain
      if (stl_model.vertices[vid_00][2] > -1000 &&
          stl_model.vertices[vid_01][2] > -1000 &&
          stl_model.vertices[vid_10][2] > -1000 &&
          stl_model.vertices[vid_11][2] > -1000) {
        stl_model.facets.push_back({vid_00, vid_10, vid_01});
        stl_model.facets.push_back({vid_01, vid_10, vid_11});
      }
    }
  }
 
  stl_model.RemoveUnreferencedVertices();
 
  return stl_model;
}
 
// reconstruct closed stl mesh based DEM (digital elevation model) data
STLModel GetSolidSTLFromDEM(string const &dem_file, double spacing) {
  auto dem_data = IO::ImportData(dem_file, 6);
 
  int num_rows = dem_data.size();
  int num_cols = dem_data[0].size();
 
  double length = (num_cols - 1) * spacing;
  double width = (num_rows - 1) * spacing;
 
  STLModel stl_model;
  stl_model.vertices.reserve(num_rows * num_cols * 2);
 
  double min_dem{Math::Infinity};
  for (int i = 0; i < num_rows; i++) {
    for (int j = 0; j < num_cols; j++) {
      if (dem_data[i][j] > -1000) {
        min_dem = min(min_dem, dem_data[i][j]);
      }
      stl_model.vertices.push_back({-0.5 * length + j * spacing,
                                    0.5 * width - i * spacing, dem_data[i][j]});
    }
  }
  for (int i = 0; i < num_rows; i++) {
    for (int j = 0; j < num_cols; j++) {
      if (dem_data[i][j] < -1000) {
        stl_model.vertices[i * num_cols + j][2] = min_dem - 0.5;
      }
      stl_model.vertices.push_back({-0.5 * length + j * spacing,
                                    0.5 * width - i * spacing, min_dem - 1.0});
    }
  }
 
  int num_dem_points = num_rows * num_cols;
  for (int i = 0; i < num_rows - 1; i++) {
    for (int j = 0; j < num_cols - 1; j++) {
      auto vid_00 = i * num_cols + j;
      auto vid_01 = i * num_cols + j + 1;
      auto vid_10 = (i + 1) * num_cols + j;
      auto vid_11 = (i + 1) * num_cols + j + 1;
 
      stl_model.facets.push_back({vid_00, vid_10, vid_01});
      stl_model.facets.push_back({vid_01, vid_10, vid_11});
 
      stl_model.facets.push_back({vid_00 + num_dem_points,
                                  vid_01 + num_dem_points,
                                  vid_10 + num_dem_points});
      stl_model.facets.push_back({vid_01 + num_dem_points,
                                  vid_11 + num_dem_points,
                                  vid_10 + num_dem_points});
    }
  }
 
  for (int j = 0; j < num_cols - 1; j++) {
    auto vid_00 = j;
    auto vid_01 = j + 1;
    auto vid_10 = vid_00 + num_dem_points;
    auto vid_11 = vid_01 + num_dem_points;
 
    stl_model.facets.push_back({vid_00, vid_01, vid_10});
    stl_model.facets.push_back({vid_10, vid_01, vid_11});
 
    stl_model.facets.push_back({vid_00 + (num_dem_points - num_cols),
                                vid_10 + (num_dem_points - num_cols),
                                vid_01 + (num_dem_points - num_cols)});
    stl_model.facets.push_back({vid_10 + (num_dem_points - num_cols),
                                vid_11 + (num_dem_points - num_cols),
                                vid_01 + (num_dem_points - num_cols)});
  }
 
  for (int i = 0; i < num_rows - 1; i++) {
    auto vid_00 = i * num_cols;
    auto vid_01 = (i + 1) * num_cols;
    auto vid_10 = vid_00 + num_dem_points;
    auto vid_11 = vid_01 + num_dem_points;
 
    stl_model.facets.push_back({vid_00, vid_10, vid_01});
    stl_model.facets.push_back({vid_10, vid_11, vid_01});
 
    stl_model.facets.push_back(
        {vid_00 + num_cols - 1, vid_01 + num_cols - 1, vid_10 + num_cols - 1});
    stl_model.facets.push_back(
        {vid_10 + num_cols - 1, vid_01 + num_cols - 1, vid_11 + num_cols - 1});
  }
 
  return stl_model;
}
 
// use pre difference post to get the slide region
STLModel GetSlideRegion(string const &dem_file_pre,
                        string const &dem_file_slide, double spacing) {
  auto stl_pre = GetSolidSTLFromDEM(dem_file_pre, 1.0);
  auto stl_slide = GetSolidSTLFromDEM(dem_file_slide, 1.0);
 
  for (int i = 0; i < stl_pre.vertices.size() / 2; i++) {
    if (abs(stl_pre.vertices[i][2] - stl_slide.vertices[i][2]) < spacing) {
      stl_slide.vertices[i][2] = stl_pre.vertices[i][2] + 0.5;
    }
  }
 
  for (int i = stl_pre.vertices.size() / 2; i < stl_pre.vertices.size(); i++) {
    stl_slide.vertices[i][2] += 150;
  }
  for (auto &facet : stl_slide.facets) {
    std::swap(facet[1], facet[2]);
  }
 
  STLModel stl_slide_region;
  IGLWrapper::MeshIntersect(
      stl_pre.vertices, stl_pre.facets, stl_slide.vertices, stl_slide.facets,
      &(stl_slide_region.vertices), &(stl_slide_region.facets));
  return stl_slide_region;
}
 
// discretize slide region.
VecXT<Vec3d> DiscretizeSlideRegion(string const &dem_file_pre,
                                   string const &dem_file_slide,
                                   double spacing) {
  auto dem_data_pre = IO::ImportData(dem_file_pre, 6);
  auto dem_data_slide = IO::ImportData(dem_file_slide, 6);
 
  int num_rows = dem_data_pre.size();
  int num_cols = dem_data_pre[0].size();
 
  double length = (num_cols - 1) * spacing;
  double width = (num_rows - 1) * spacing;
 
  VecXT<Vec3d> pos_list;
  for (int i = 1; i < num_rows - 1; i++) {
    for (int j = 1; j < num_cols - 1; j++) {
      double z_min = -Math::Infinity;
      double z_max = Math::Infinity;
      for (int ii = -1; ii <= 1; ii++) {
        for (int jj = -1; jj <= 1; jj++) {
          z_min = max(z_min, dem_data_slide[i + ii][j + jj]);
          z_max = min(z_max, dem_data_pre[i + ii][j + jj]);
        }
      }
 
      for (int ii = 0; ii < (z_max - z_min) / spacing - 1; ii++) {
        pos_list.push_back({-0.5 * length + j * spacing,
                            0.5 * width - i * spacing,
                            z_min + (ii + 0.5) * spacing});
      }
    }
  }
  return pos_list;
}
 
int main(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(-150, -150, 1550, 150, 150, 1750);
  sim->domain_manager.SetDecomposition(total_rank, 1, 1);
  sim->domain_manager.SetCellSpacing(5.0, 5.0, 5.0);
 
  LinearSpring cnt_model = LinearSpring(2.0e5, 1.0e5, 0.5, 0.0);
  cnt_model.label = "cnt_model_pp";
  sim->scene.InsertContactModel(&cnt_model);
 
  cnt_model = LinearSpring(2.0e7, 1.0e7, 0.5, 0.0);
  cnt_model.label = "cnt_model_pw";
  sim->scene.InsertContactModel(&cnt_model);
 
  sim->scene.SetNumberOfMaterials(2);
  sim->scene.SetCollisionModel(0, 0, "cnt_model_pp");
  sim->scene.SetCollisionModel(1, 0, "cnt_model_pw");
  sim->scene.SetCollisionModel(1, 1, "cnt_model_pw");
 
  // insert terrain triangles and walls to the scene
  filesystem::create_directories("tmp/examples/landslide/");
  auto terrain_stl =
      GetSTLFromDEM("data/copyleft/landslide/dc7_dem_1411slide.txt", 1.0);
  terrain_stl.SaveAsSTL("tmp/examples/landslide/terrain_stl.stl");
 
  auto terrain_stl_post =
      GetSTLFromDEM("data/copyleft/landslide/dc7_dem_1411post.txt", 1.0);
  terrain_stl_post.SaveAsSTL("tmp/examples/landslide/terrain_stl_post.stl");
 
  for (auto &facet : terrain_stl.facets) {
    Triangle triangle(terrain_stl.vertices[facet[0]],
                      terrain_stl.vertices[facet[1]],
                      terrain_stl.vertices[facet[2]]);
    auto shape_ptr = sim->scene.InsertShape(&triangle);
    Wall wall(shape_ptr);
    wall.material_type = 1;
    sim->scene.InsertWall(wall);
  }
 
  // // for visualize the slide region
  // auto slide_stl =
  //     GetSlideRegion("data/copyleft/landslide/dc7_dem_1411pre.txt",
  //                    "data/copyleft/landslide/dc7_dem_1411slide.txt", 1.0);
  // slide_stl.SaveAsSTL("tmp/examples/landslide/slide_stl.stl");
 
  // discretize and insert particles
  auto pos_list = DiscretizeSlideRegion(
      "data/copyleft/landslide/dc7_dem_1411pre.txt",
      "data/copyleft/landslide/dc7_dem_1411slide.txt", 1.0);
  cout << "num particles: " << pos_list.size() << endl;
 
  /* write particle position to txt
  stringbuf buf;
  ostream os(&buf);
 
  for (auto &pos : pos_list) {
    os.width(12);
    os << pos[0] << " ";
    os.width(12);
    os << pos[1] << " ";
    os.width(12);
    os << pos[2] << endl;
  }
 
  ofstream outfile;
  outfile.open("tmp/examples/landslide/pos_list.txt");
  outfile << buf.str();
  outfile.close();
  // */
 
  Sphere sphere(1.0);
  sphere.SetSkinFactor(0.5);
  auto sphere_ptr = sim->scene.InsertShape(&sphere);
  for (int i = 0; i < pos_list.size(); i++) {
    auto &pos = pos_list[i];
    Particle p(sphere_ptr);
    p.SetPosition(pos[0], pos[1], pos[2]);
    p.damp_numerical = 0.7;
    sim->scene.InsertParticle(p);
  }
 
  // gravity
  Gravity grav;
  grav.Init(sim);
  sim->modifier_manager.Insert(&grav);
  sim->modifier_manager.Enable(grav.label);
 
  // save results
  DataDumper data_dumper;
  data_dumper.Init(sim);
  data_dumper.SetRootPath("tmp/examples/landslide/");
  data_dumper.SetSaveByCycles(100);
  // data_dumper.dump_wall_info = true;
  // data_dumper.dump_shape_info = true;
  sim->modifier_manager.Insert(&data_dumper);
  sim->modifier_manager.Enable(data_dumper.label);
 
  // run
  sim->dem_solver.timestep = 1.0e-3;
  sim->Run(40.0);
 
  // finilize
  delete sim;
 
  MPI_Finalize();
  return 0;
}