01_sdf_error.cpp

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

#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 SDFError(int argc, char **argv) {
  string root_dir = "tmp/examples/sdf_dem/sdf_error/";
  filesystem::create_directories(root_dir);
 
  PolySuperEllipsoid ellipsoid(0.5, 1, 1.5, 0.5, 1, 1.5, 1.5, 1.5);
 
  auto stl_model = ellipsoid.GetSTLModel(1000);
  stl_model.SaveAsVTK(root_dir + "surface.vtk");
 
  TriMesh trimesh;
  trimesh.InitFromSTL(stl_model);
 
  double mesh_size =
      pow(stl_model.GetVolume() / Math::PI * 6.0, 1.0 / 3.0) / 30.0;
  TetMesh tetmesh(stl_model.vertices, stl_model.facets, mesh_size);
 
  // calculate the signed distance
  VecXT<Vec3d> dist_list;
  for (auto &node : tetmesh.nodes) {
    auto dist_ellipsoid = ellipsoid.SignedDistance(node);
    auto dist_trimesh = trimesh.SignedDistance(node);
    dist_list.push_back(
        {dist_ellipsoid, dist_trimesh, abs(dist_ellipsoid - dist_trimesh)});
  }
 
  // save the results
  tetmesh.SaveAsVTK(root_dir + "tetmesh.vtk");
 
  stringbuf buf;
  ostream os(&buf);
 
  os << "POINT_DATA " << dist_list.size() << endl;
 
  os << "SCALARS SDF float 3" << endl;
  os << "LOOKUP_TABLE default" << endl;
  for (int i = 0; i < dist_list.size(); i++) {
    os << dist_list[i][0] << " " << dist_list[i][1] << " " << dist_list[i][2]
       << 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();
}