03_friction_test_sphere.cpp

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

#include "netdem.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);
 
void FrictionTestSphere(int potential_case) {
  // load particle
  SphericalHarmonics sh;
  sh.a_nm[0] = sqrt(Math::PI);
  sh.UpdateShapeProperties();
 
  sh.EnableSurfaceNodes();
  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);
 
  obj_p2.SetPosition(0, 0, 0.95);
  obj_p2.SetVelocity(0.5, 0, 0);
 
  // solver settings
  SolverSDFPP cnt_solver;
  cnt_solver.potential_type = potential_case;
  cnt_solver.use_equivalent_stiffness = false;
 
  LinearSpring cnt_model(1, 1, 0.3, 0);
  ContactPP *cnt = nullptr;
 
  double timestep = 1.0e-4;
 
  VecXT<VecXT<double>> cnt_data;
  for (int i = 0; i < 1001; i++) {
    // clear the particle force and moment
    obj_p1.ClearForce();
    obj_p1.ClearMoment();
 
    obj_p2.ClearForce();
    obj_p2.ClearMoment();
 
    // translational moving
    obj_p2.UpdateMotion(timestep);
 
    // 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, timestep);
      } else {
        cnt_solver.ResolveUpdate(cnt, timestep);
      }
 
      cnt->EvaluateForces(timestep);
      cnt->ApplyToParticle();
 
      double pos_x = obj_p2.pos[0];
 
      Vec3d cnt_force{0, 0, 0}, cnt_moment{0, 0, 0};
      Vec3d dir_n{0, 0, 0}, dir_s{0, 0, 0};
      for (auto &cnt_entry : cnt->collision_entries) {
        cnt_force[0] += cnt_entry.cnt_forces.force[0];
        cnt_force[1] += cnt_entry.cnt_forces.force[1];
        cnt_force[2] += cnt_entry.cnt_forces.force[2];
 
        dir_n[0] += cnt_entry.cnt_forces.force_n[0];
        dir_n[1] += cnt_entry.cnt_forces.force_n[1];
        dir_n[2] += cnt_entry.cnt_forces.force_n[2];
 
        dir_s[0] += cnt_entry.cnt_forces.force_t[0];
        dir_s[1] += cnt_entry.cnt_forces.force_t[1];
        dir_s[2] += cnt_entry.cnt_forces.force_t[2];
      }
      Math::Normalize(&dir_n);
      Math::Normalize(&dir_s);
 
      cout << pos_x << ", " << cnt->collision_entries.size() << ",    "
           << cnt_force << ",    " << dir_n << ",    " << dir_s << endl;
 
      cnt_data.push_back({pos_x, double(cnt->collision_entries.size()),
                          cnt_force[0], cnt_force[1], cnt_force[2], dir_n[0],
                          dir_n[1], dir_n[2], dir_s[0], dir_s[1], dir_s[2]});
    } else {
      double pos_x = obj_p2.pos[0];
      cnt_data.push_back({pos_x, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0});
    }
  }
 
  string root_dir = "local/potential_models/friction_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;
  }
}