00_sdf_demo.cpp

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

#include "dem_profiler.hpp"
#include "distribution_uniform.hpp"
#include "shape_level_set.hpp"
#include "shape_poly_super_ellipsoid.hpp"
#include "shape_poly_super_quadrics.hpp"
#include "shape_spherical_harmonics.hpp"
#include "shape_trimesh.hpp"
#include "stl_model.hpp"
#include "tetmesh.hpp"
#include "utils_math.hpp"
#include <filesystem>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
 
using namespace netdem;
using namespace std;
 
void SDFDemo(int argc, char **argv) {
  Shape *shape{nullptr};
 
  int shape_case;
  if (argc == 2) {
    shape_case = -1;
  } else {
    shape_case = atof(argv[2]);
  }
 
  switch (shape_case) {
  case 0: {
    // convex
    shape = new PolySuperEllipsoid(0.5, 1, 1.5, 0.5, 1, 1.5, 1.5, 1.5);
  } break;
  case 1: {
    // convex
    shape = new PolySuperQuadrics(0.5, 1.0, 2.5, 0.5, 1.7, 0.5, 0.5, 1.0, 1.5,
                                  1.0, 1.2, 0.8);
  } break;
  case 10: {
    // non-convex
    shape = new PolySuperEllipsoid(0.5, 1, 1.5, 0.5, 1, 1.5, 2.5, 2.5);
  } break;
  case 11: {
    // non-convex
    shape = new PolySuperQuadrics(0.5, 1.0, 2.5, 0.5, 1.7, 0.5, 0.5, 1.0, 1.5,
                                  2.0, 2.2, 1.8);
  } break;
  case 2: {
    auto sh = new SphericalHarmonics();
    sh->InitFromSTL("data/particle_template.stl");
    sh->SetSize(1.0);
    shape = sh;
  } break;
  case 3: {
    auto trimesh = new TriMesh();
    trimesh->InitFromSTL("data/copyleft/bolt_nut/bolt.stl");
    // trimesh->Decimate(8000);
    trimesh->AlignAxes();
    trimesh->SetSize(1.0);
    shape = trimesh;
  } break;
  case 4: {
    auto level_set = new LevelSet();
    level_set->InitFromSTL("data/particle_template.stl");
    level_set->AlignAxes();
    level_set->SetSize(1.0);
    shape = level_set;
  } break;
  default: {
    auto trimesh = new TriMesh();
    trimesh->InitFromSTL("data/copyleft/bolt_nut/bolt.stl");
    trimesh->AlignAxes();
    trimesh->SetSize(1.0);
    shape = trimesh;
  } break;
  }
 
  string root_dir = "tmp/examples/sdf_dem/sdf_demo/";
  filesystem::create_directories(root_dir);
 
  auto stl_model = shape->GetSTLModel(1000);
  stl_model.SaveAsSTL(root_dir + "surface.stl");
 
  double mesh_size =
      pow(stl_model.GetVolume() / Math::PI * 6.0, 1.0 / 3.0) / 20.0;
  TetMesh tetmesh(stl_model.vertices, stl_model.facets, mesh_size);
 
  stl_model.vertices = tetmesh.surface_nodes;
  stl_model.facets = tetmesh.surface_facets;
  stl_model.SaveAsSTL(root_dir + "surface_mesh.stl");
 
  // calculate the signed distance
  VecXT<double> dist_list;
  for (auto &node : tetmesh.nodes) {
    dist_list.push_back(shape->SignedDistance(node));
  }
 
  DEMProfiler dem_profiler;
  dem_profiler.Clear();
  dem_profiler.StartTimer(DEMProfiler::TimerType::custom);
 
  auto [bound_aabb_min, bound_aabb_max] = shape->GetBoundAABB();
  UniformDistribution uniform_dist_x(bound_aabb_min[0], bound_aabb_max[0]);
  UniformDistribution uniform_dist_y(bound_aabb_min[1], bound_aabb_max[1]);
  UniformDistribution uniform_dist_z(bound_aabb_min[2], bound_aabb_max[1]);
 
  int num_trials = 1e3;
  auto pos_x_list = uniform_dist_x.Get(num_trials);
  auto pos_y_list = uniform_dist_y.Get(num_trials);
  auto pos_z_list = uniform_dist_z.Get(num_trials);
 
  for (int i = 0; i < num_trials; i++) {
    for (auto &node : tetmesh.nodes) {
      Vec3d pos{pos_x_list[i], pos_y_list[i], pos_z_list[i]};
      shape->SignedDistance(pos);
    }
    if ((i + 1) % 100 == 0) {
      cout << i + 1 << endl;
    }
  }
  dem_profiler.EndTimer(DEMProfiler::TimerType::custom);
  cout << "no. nodes queried: " << tetmesh.nodes.size() << endl;
  printf("cost: %f microsecond \n",
         dem_profiler.timer_list[DEMProfiler::TimerType::custom] / 1.0 /
             num_trials / tetmesh.nodes.size());
 
  // save results
  tetmesh.SaveAsVTK(root_dir + "tetmesh.vtk");
 
  stringbuf buf;
  ostream os(&buf);
 
  os << "POINT_DATA " << dist_list.size() << endl;
 
  os << "SCALARS SDF float 1" << endl;
  os << "LOOKUP_TABLE default" << endl;
  for (int i = 0; i < dist_list.size(); i++) {
    os << dist_list[i] << endl;
  }
 
  ofstream outfile;
  outfile.open(root_dir + "tetmesh.vtk", ios::app);
  if (!outfile) {
    cout << "cannot open file: " << root_dir + "tetmesh.vtk" << endl;
  }
  outfile << buf.str();
  outfile.close();
 
  delete shape;
}