11_contact_test_trimesh.cpp

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

#include "contact_solver_factory.hpp"
#include "cork_wrapper.hpp"
#include "distribution_uniform.hpp"
#include "igl_wrapper.hpp"
#include "model_linear_spring.hpp"
#include "particle.hpp"
#include "shape_spherical_harmonics.hpp"
#include "solver_sdf_pp.hpp"
#include "spherical_voronoi.hpp"
#include "utils_math.hpp"
#include <filesystem>
#include <fstream>
#include <iostream>
#include <random>
#include <sstream>
#include <string>
using namespace netdem;
using namespace std;
 
void SaveDataset(string filename, VecXT<VecXT<double>> &data);
 
double GetContactVolume(STLModel &stl_a, STLModel &stl_b, double z_coord) {
  auto va = stl_a.vertices;
  auto fa = stl_a.facets;
 
  auto vb = stl_b.vertices;
  auto fb = stl_b.facets;
  for (auto &vert : vb) {
    vert[2] += z_coord;
  }
 
  VecXT<Vec3d> vab;
  VecXT<Vec3i> fab;
  CorkWrapper::MeshIntersect(va, fa, vb, fb, &vab, &fab);
 
  return STLModel::GetVolume(vab, fab);
}
 
void ContactTestTrimesh(int potential_case, int surface_node_num) {
  // load particle
  SphericalHarmonics sh;
  sh.InitFromSTL("data/particle_template.stl");
  sh.SetSize(1.0);
 
  sh.EnableSurfaceNodes();
  auto surf_stl = sh.GetSTLModel(surface_node_num);
  sh.SetSurfaceNodes(surf_stl);
  cout << "shape size: " << sh.GetSize() << endl;
 
  Particle obj_p1 = Particle(&sh);
  Particle obj_p2 = Particle(&sh);
  cout << "particle created ... " << endl;
 
  // try another particle configuration
  obj_p1.SetRodrigues(0.3, 0, 1, 0);
  obj_p2.SetRodrigues(1.2, 0, 1, 0);
 
  // obtain a mesh representation of the nodes
  auto facets = surf_stl.facets;
 
  auto nodes = sh.GetSurfaceNodes();
  for (auto &vert : nodes) {
    vert = Math::Rotate(vert, obj_p1.quaternion);
  }
  STLModel stl_model_1(nodes, facets);
 
  nodes = sh.GetSurfaceNodes();
  for (auto &vert : nodes) {
    vert = Math::Rotate(vert, obj_p2.quaternion);
  }
  STLModel stl_model_2(nodes, facets);
 
  // solver settings
  SolverSDFPP cnt_solver;
  cnt_solver.potential_type = potential_case;
 
  LinearSpring cnt_model(1, 1, 0, 0);
  ContactPP *cnt = nullptr;
 
  // gradually push particle 2 onto particle 1
  double pos_z0 = 1.015;
  VecXT<VecXT<double>> cnt_data;
  for (int i = 0; i < 61; i++) {
    obj_p1.ClearForce();
    obj_p1.ClearMoment();
 
    obj_p2.ClearForce();
    obj_p2.ClearMoment();
 
    obj_p2.SetPosition(0, 0, pos_z0 - i * 0.001);
 
    cnt_solver.Init(&obj_p2, &obj_p1);
    if (cnt_solver.Detect()) {
      if (cnt == nullptr) {
        cnt = new ContactPP(&obj_p2, &obj_p1);
        cnt->SetCollisionModel(&cnt_model);
        cnt_solver.ResolveInit(cnt, 1.0e-4);
      } else {
        cnt_solver.ResolveUpdate(cnt, 1.0e-4);
      }
 
      cnt->EvaluateForces(1.0e-4);
      cnt->ApplyToParticle();
 
      double len_n = 1.0 - obj_p2.pos[2];
      double vol = GetContactVolume(stl_model_1, stl_model_2, obj_p2.pos[2]);
 
      cout << len_n << ", " << cnt->collision_entries.size() << ", "
           << obj_p2.force[0] << ", " << obj_p2.force[1] << ", "
           << obj_p2.force[2] << ", " << obj_p2.moment[0] << ", "
           << obj_p2.moment[1] << ", " << obj_p2.moment[2] << ", " << vol
           << endl;
 
      cnt_data.push_back({len_n, double(cnt->collision_entries.size()),
                          obj_p2.force[0], obj_p2.force[1], obj_p2.force[2],
                          obj_p2.moment[0], obj_p2.moment[1], obj_p2.moment[2],
                          vol});
    }
  }
 
  for (int i = 0; i < 61; i++) {
    obj_p1.ClearForce();
    obj_p1.ClearMoment();
 
    obj_p2.ClearForce();
    obj_p2.ClearMoment();
 
    obj_p2.SetPosition(0, 0, pos_z0 - 60 * 0.001 + i * 0.001);
 
    cnt_solver.Init(&obj_p2, &obj_p1);
 
    if (cnt_solver.Detect()) {
      if (cnt == nullptr) {
        cnt = new ContactPP(&obj_p2, &obj_p1);
        cnt->SetCollisionModel(&cnt_model);
        cnt_solver.ResolveInit(cnt, 1.0e-4);
      } else {
        cnt_solver.ResolveUpdate(cnt, 1.0e-4);
      }
 
      cnt->EvaluateForces(1.0e-4);
      cnt->ApplyToParticle();
 
      double len_n = 1.0 - obj_p2.pos[2];
      double vol = GetContactVolume(stl_model_1, stl_model_2, obj_p2.pos[2]);
 
      cout << len_n << ", " << cnt->collision_entries.size() << ", "
           << obj_p2.force[0] << ", " << obj_p2.force[1] << ", "
           << obj_p2.force[2] << ", " << obj_p2.moment[0] << ", "
           << obj_p2.moment[1] << ", " << obj_p2.moment[2] << ", " << vol
           << endl;
 
      cnt_data.push_back({len_n, double(cnt->collision_entries.size()),
                          obj_p2.force[0], obj_p2.force[1], obj_p2.force[2],
                          obj_p2.moment[0], obj_p2.moment[1], obj_p2.moment[2],
                          vol});
    }
  }
 
  string root_dir = "tmp/examples/sdf_dem/contact_test/";
  filesystem::create_directories(root_dir);
  switch (potential_case) {
  case SolverSDFPP::PotentialType::linear:
    SaveDataset(root_dir + "trimesh_linear.txt", cnt_data);
    break;
  case SolverSDFPP::PotentialType::hertz:
    SaveDataset(root_dir + "trimesh_hertz.txt", cnt_data);
    break;
  default:
    SaveDataset(root_dir + "trimesh_hertz.txt", cnt_data);
    break;
  }
 
  stl_model_1.SaveAsSTL(root_dir + "p1.stl");
  for (auto &vert : stl_model_2.vertices) {
    vert[2] += pos_z0 - 60 * 0.001;
  }
  stl_model_2.SaveAsSTL(root_dir + "p2.stl");
 
  if (cnt != nullptr) {
    delete cnt;
    cnt = nullptr;
  }
}