00_contact_test_sphere.cpp

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

#include "contact_solver_factory.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) {
  stringbuf buf;
  ostream os(&buf);
  int os_width = 24;
  os.setf(ios::scientific);
  os.precision(15);
  for (auto &data_row : data) {
    for (int j = 0; j < data_row.size() - 1; j++) {
      os.width(os_width);
      os << data_row[j] << ", ";
    }
    os.width(os_width);
    os << data_row.back() << endl;
  }
 
  ofstream outfile;
  outfile.open(filename);
  if (!outfile) {
    cout << "cannot open or create file: " << filename << endl;
    abort();
  }
  outfile << buf.str();
  outfile.close();
  cout << "data saved to: " << filename << endl;
}
 
void ContactTestSphere(int potential_case, int surface_node_num) {
  // load particle
  SphericalHarmonics sh;
  sh.a_nm[0] = sqrt(Math::PI);
  sh.UpdateShapeProperties();
 
  sh.EnableSurfaceNodes();
  sh.SetSurfaceNodeNum(surface_node_num);
  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.4, 0, 1, 0);
  obj_p2.SetRodrigues(0.6, 0, 1, 0);
 
  // solver settings
  SolverSDFPP cnt_solver;
  cnt_solver.potential_type = potential_case;
  cnt_solver.use_equivalent_stiffness = false;
  double dt = 1.0e-4;
 
  LinearSpring cnt_model(1, 1, 0, 0);
  ContactPP *cnt = nullptr;
 
  // gradually push particle 2 onto particle 1
  VecXT<VecXT<double>> cnt_data;
  for (int i = 0; i < 61; i++) {
    // clear the particle force and moment
    obj_p1.ClearForce();
    obj_p1.ClearMoment();
 
    obj_p2.ClearForce();
    obj_p2.ClearMoment();
 
    // push down particle 2
    obj_p2.SetPosition(0, 0, 1.01 - i * 0.001);
 
    // contact detection and resolution procedure
    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, dt);
      } else {
        cnt_solver.ResolveUpdate(cnt, dt);
      }
 
      cnt->EvaluateForces(dt);
      cnt->ApplyToParticle();
 
      double len_n = 1.0 - obj_p2.pos[2];
 
      cout << len_n << ", " << cnt->collision_entries.size() << ",    "
           << obj_p2.force << ",    " << obj_p2.moment << 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]});
    }
  }
 
  // gradually move particle 2 away particle 1
  for (int i = 0; i < 61; i++) {
    obj_p1.ClearForce();
    obj_p1.ClearMoment();
 
    obj_p2.ClearForce();
    obj_p2.ClearMoment();
 
    obj_p2.SetPosition(0, 0, 1.01 - 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, dt);
      } else {
        cnt_solver.ResolveUpdate(cnt, dt);
      }
 
      cnt->EvaluateForces(dt);
      cnt->ApplyToParticle();
 
      double len_n = 1.0 - obj_p2.pos[2];
 
      cout << len_n << ", " << cnt->collision_entries.size() << ",    "
           << obj_p2.force << ",    " << obj_p2.moment << 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]});
    }
  }
 
  string root_dir = "local/potential_models/contact_test/";
  filesystem::create_directories(root_dir);
  switch (potential_case) {
  case SolverSDFPP::PotentialType::linear:
    SaveDataset(root_dir + "sphere_linear.txt", cnt_data);
    break;
  case SolverSDFPP::PotentialType::hertz:
    SaveDataset(root_dir + "sphere_hertz.txt", cnt_data);
    break;
  default:
    SaveDataset(root_dir + "sphere_hertz.txt", cnt_data);
    break;
  }
 
  if (cnt != nullptr) {
    delete cnt;
    cnt = nullptr;
  }
}