81_golf_ball_impact.cpp

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

#include "data_dumper.hpp"
#include "deformable_particle.hpp"
#include "deformation_analysis.hpp"
#include "gen_pack.hpp"
#include "gen_wall_box_plane.hpp"
#include "gravity.hpp"
#include "model_hertz_mindlin.hpp"
#include "model_linear_spring.hpp"
#include "shape_ellipsoid.hpp"
#include "shape_sphere.hpp"
#include "simulation.hpp"
#include "solver_sdf_pw.hpp"
#include <filesystem>
#include <iostream>
#include <unordered_map>
 
using namespace netdem;
using namespace std;
 
void GolfBallImpact() {
  string result_dir = "tmp/examples/deformable_particle/golf_ball/";
  filesystem::create_directories(result_dir);
 
  // golf ball
  Sphere sphere(0.04);
  auto stl_model = sphere.GetSTLModel(500);
 
  TriMesh trimesh;
  trimesh.InitFromSTL(stl_model);
 
  DeformableParticle p;
 
  p.mesh_res = 20;
  p.SetShape(&trimesh);
  p.SetDensity(1150);
 
  // E = 107 MPa, nu = 0.45
  p.fem_simulator.neo_k = 3.57e8 / 5;
  p.fem_simulator.neo_mu = 3.69e7 / 5;
  p.fem_simulator.damp_coef = 0.3;
  p.fem_simulator.timestep = 1.0e-8;
 
  p.SetPosition(-0.05, 0, 0);
  p.SetVelocity(-67, 0, 0);
 
  // plane wall
  Plane plane(-0.1, 0, 0, 1, 0, 0);
  Wall w(&plane);
 
  // contact model and solver
  LinearSpring cnt_model = LinearSpring(1.0e8, 1.0e8, 0.0, 0.0);
  SolverSDFPW cnt_solver;
  cnt_solver.potential_type = 0;
  cnt_solver.solve_two_sides = true;
  double dem_timestep = 1.0e-7;
 
  // simulation
  p.SaveAsVTK(result_dir + "p_tetmesh_000.vtk");
  p.SaveSurfaceAsVTK(result_dir + "p_trimesh_000.vtk");
  for (int ti = 0; ti < 400; ti++) {
    for (int j = 0; j < 50; j++) {
      p.ClearForce();
 
      cnt_solver.Init(&p, &w);
      if (cnt_solver.Detect()) {
        ContactPW cnt = ContactPW(&p, &w);
        cnt.SetCollisionModel(&cnt_model);
 
        cnt_solver.ResolveInit(&cnt, dem_timestep);
        cnt.EvaluateForces(dem_timestep);
        p.ApplyContactForce(&cnt);
      }
      p.UpdateMotion(dem_timestep);
    }
 
    char filename[128];
    snprintf(filename, 128, "p_tetmesh_%03d.vtk", ti + 1);
    p.SaveAsVTK(result_dir + filename);
    snprintf(filename, 128, "p_trimesh_%03d.vtk", ti + 1);
    p.SaveSurfaceAsVTK(result_dir + filename);
  }
}