00_hertz_contact.cpp

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

#include "deformable_particle.hpp"
#include "gen_wall_box_plane.hpp"
#include "model_linear_spring.hpp"
#include "shape_sphere.hpp"
#include "solver_sdf_pw.hpp"
#include <filesystem>
#include <string>
 
using namespace netdem;
using namespace std;
 
void HertzContact() {
  TriMesh trimesh;
  trimesh.InitFromSTL("data/sphere.stl");
  trimesh.SetSize(1.0);
 
  DeformableParticle p;
  p.SetShape(&trimesh);
 
  Plane plane_1(0, 0, -0.45, 0, 0, 1);
  Wall w_1 = Wall(&plane_1);
 
  Plane plane_2(0, 0, 0.45, 0, 0, -1);
  Wall w_2 = Wall(&plane_2);
 
  LinearSpring cnt_model = LinearSpring(2.0e6, 0.0, 0.0, 0.0);
  SolverSDFPW cnt_solver(&p, &w_1);
  cnt_solver.potential_type = 0;
 
  double dt = 1.0e-4;
  p.fem_simulator.gravity_coef[2] = 0.0;
  p.fem_simulator.timestep = 1.0e-5;
 
  string root_dir = "tmp/examples/deformable_particle/hertz_contact/";
  filesystem::create_directories(root_dir);
  p.SaveAsVTK(root_dir + "particle_000.vtk");
 
  double fz = 0;
  for (int ti = 0; ti < 100; ti++) {
    for (int tj = 0; tj < 50; tj++) {
      p.ClearForce();
 
      cnt_solver.Init(&p, &w_1);
      if (cnt_solver.Detect()) {
        ContactPW cnt = ContactPW(&p, &w_1);
        cnt.SetCollisionModel(&cnt_model);
 
        cnt_solver.ResolveInit(&cnt, dt);
        cnt.EvaluateForces(dt);
        for (int ii = 0; ii < cnt.collision_entries.size(); ii++) {
          p.AddForce(cnt.collision_entries[ii].cnt_geoms.node_id,
                     cnt.collision_entries[ii].cnt_forces.force);
        }
      }
      cnt_solver.Init(&p, &w_2);
      fz = 0;
      if (cnt_solver.Detect()) {
        ContactPW cnt = ContactPW(&p, &w_2);
        cnt.SetCollisionModel(&cnt_model);
 
        cnt_solver.ResolveInit(&cnt, dt);
        cnt.EvaluateForces(dt);
        for (int ii = 0; ii < cnt.collision_entries.size(); ii++) {
          p.AddForce(cnt.collision_entries[ii].cnt_geoms.node_id,
                     cnt.collision_entries[ii].cnt_forces.force);
          fz += cnt.collision_entries[ii].cnt_forces.force[2];
        }
      }
 
      p.UpdateMotion(dt);
    }
 
    cout << "contact force z:" << abs(fz) << endl;
    char filename[128];
    snprintf(filename, 128, "particle_%03d.vtk", ti + 1);
    p.SaveAsVTK(root_dir + filename);
  }
}