mychgpu.cpp

// =============================================================================
// PROJECT CHRONO - http://projectchrono.org
//
// Copyright (c) 2019 projectchrono.org
// All rights reserved.
//
// Use of this source code is governed by a BSD-style license that can be found
// in the LICENSE file at the top level of the distribution and at
// http://projectchrono.org/license-chrono.txt.
//
// =============================================================================
// Authors: Nic Olsen
// =============================================================================
// Chrono::Gpu demo using SMC method. A body whose geometry is described by an
// OBJ file is time-integrated in Chrono and interacts with a granular wave tank
// in Chrono::Gpu via the co-simulation framework.
// =============================================================================

#include <iostream>
#include <string>
#include <vector>

#include "chrono/assets/ChSphereShape.h"
#include "chrono/core/ChGlobal.h"
#include "chrono/physics/ChBody.h"
#include "chrono/physics/ChForce.h"
#include "chrono/physics/ChSystemSMC.h"
#include "chrono/timestepper/ChTimestepper.h"
#include "chrono/utils/ChUtilsCreators.h"
#include "chrono/utils/ChUtilsSamplers.h"

#include "chrono_gpu/ChGpuData.h"
#include "chrono_gpu/physics/ChSystemGpu.h"
#include "chrono_gpu/utils/ChGpuJsonParser.h"
#include "chrono_gpu/utils/ChGpuVisualization.h"

#include "chrono_thirdparty/filesystem/path.h"

using namespace chrono;
using namespace chrono::gpu;

// Output frequency
float out_fps = 50;

// Enable/disable run-time visualization (if Chrono::OpenGL is available)
bool render = true;
float render_fps = 2000;

void writeMeshFrames(std::ostringstream &outstream, ChBody &body,
                     std::string obj_name, float mesh_scaling) {
  outstream << obj_name << ",";

  // Get frame position
  ChFrame<> body_frame = body.GetFrame_REF_to_abs();
  ChQuaternion<> rot = body_frame.GetRot();
  ChVector<> pos = body_frame.GetPos();

  // Get basis vectors
  ChVector<> vx = rot.GetXaxis();
  ChVector<> vy = rot.GetYaxis();
  ChVector<> vz = rot.GetZaxis();

  // Output in order
  outstream << pos.x() << ",";
  outstream << pos.y() << ",";
  outstream << pos.z() << ",";
  outstream << vx.x() << ",";
  outstream << vx.y() << ",";
  outstream << vx.z() << ",";
  outstream << vy.x() << ",";
  outstream << vy.y() << ",";
  outstream << vy.z() << ",";
  outstream << vz.x() << ",";
  outstream << vz.y() << ",";
  outstream << vz.z() << ",";
  outstream << mesh_scaling << "," << mesh_scaling << "," << mesh_scaling;
  outstream << "\n";
}

int main(int argc, char *argv[]) {

  SetDataPath("../data/");

  ChGpuSimulationParameters params;
  if (argc != 2 || ParseJSON(gpu::GetDataFile(argv[1]), params) == false) {
    std::cout << "Usage:\n./mychgpu <json_file>" << std::endl;
    return 1;
  }

  float iteration_step = params.step_size;

  ChSystemGpuMesh gpu_sys(
      params.sphere_radius, params.sphere_density,
      make_float3(params.box_X, params.box_Y, params.box_Z));

  double fill_bottom = -params.box_Z / 2.0;
  double fill_top = params.box_Z / 4.0;

  chrono::utils::PDSampler<float> sampler(2.4f * params.sphere_radius);
  // chrono::utils::HCPSampler<float> sampler(2.05 * params.sphere_radius);

  // leave a 4cm margin at edges of sampling
  ChVector<> hdims(params.box_X / 2 - 4.0, params.box_Y / 2 - 4.0, 0);
  ChVector<> center(0, 0, fill_bottom + 2.0 * params.sphere_radius);
  std::vector<ChVector<float>> body_points;

  // Shift up for bottom of box
  center.z() += 3 * params.sphere_radius;
  while (center.z() < fill_top) {
    std::cout << "Create layer at " << center.z() << std::endl;
    auto points = sampler.SampleBox(center, hdims);
    body_points.insert(body_points.end(), points.begin(), points.end());
    center.z() += 2.05 * params.sphere_radius;
  }

  gpu_sys.SetParticlePositions(body_points);

  gpu_sys.SetBDFixed(true);
  std::function<double3(float)> pos_func_wave = [&params](float t) {
    double3 pos = {0, 0, 0};

    double t0 = 0.5;
    double freq = CH_C_PI / 4;

    if (t > t0) {
      pos.x = 0.1 * params.box_X * std::sin((t - t0) * freq);
    }
    return pos;
  };

  // gpu_sys.setBDWallsMotionFunction(pos_func_wave);

  gpu_sys.SetKn_SPH2SPH(params.normalStiffS2S);
  gpu_sys.SetKn_SPH2WALL(params.normalStiffS2W);
  gpu_sys.SetKn_SPH2MESH(params.normalStiffS2M);

  gpu_sys.SetGn_SPH2SPH(params.normalDampS2S);
  gpu_sys.SetGn_SPH2WALL(params.normalDampS2W);
  gpu_sys.SetGn_SPH2MESH(params.normalDampS2M);

  gpu_sys.SetKt_SPH2SPH(params.tangentStiffS2S);
  gpu_sys.SetKt_SPH2WALL(params.tangentStiffS2W);
  gpu_sys.SetKt_SPH2MESH(params.tangentStiffS2M);

  gpu_sys.SetGt_SPH2SPH(params.tangentDampS2S);
  gpu_sys.SetGt_SPH2WALL(params.tangentDampS2W);
  gpu_sys.SetGt_SPH2MESH(params.tangentDampS2M);

  gpu_sys.SetCohesionRatio(params.cohesion_ratio);
  gpu_sys.SetAdhesionRatio_SPH2MESH(params.adhesion_ratio_s2m);
  gpu_sys.SetAdhesionRatio_SPH2WALL(params.adhesion_ratio_s2w);

  gpu_sys.SetGravitationalAcceleration(
      ChVector<float>(params.grav_X, params.grav_Y, params.grav_Z));

  gpu_sys.SetFixedStepSize(params.step_size);
  gpu_sys.SetFrictionMode(CHGPU_FRICTION_MODE::MULTI_STEP);
  gpu_sys.SetTimeIntegrator(CHGPU_TIME_INTEGRATOR::CENTERED_DIFFERENCE);

  gpu_sys.SetStaticFrictionCoeff_SPH2SPH(params.static_friction_coeffS2S);
  gpu_sys.SetStaticFrictionCoeff_SPH2WALL(params.static_friction_coeffS2W);
  gpu_sys.SetStaticFrictionCoeff_SPH2MESH(params.static_friction_coeffS2M);

  // gpu_sys.SetRollingCoeff_SPH2SPH(params.rolling_friction_coeffS2S);
  // gpu_sys.SetRollingCoeff_SPH2WALL(params.rolling_friction_coeffS2W);
  // gpu_sys.SetRollingCoeff_SPH2MESH(params.rolling_friction_coeffS2M);

  std::string mesh_filename(GetChronoDataFile("sphere.obj"));
  std::vector<string> mesh_filenames(1, mesh_filename);

  std::vector<float3> mesh_translations(1, make_float3(0.f, 0.f, 0.f));

  float ball_radius = 20.f;
  std::vector<ChMatrix33<float>> mesh_rotscales(1,
                                                ChMatrix33<float>(ball_radius));

  float ball_density = params.sphere_density;
  float ball_mass = (float)(4.f * CH_C_PI * ball_radius * ball_radius *
                            ball_radius * ball_density / 3.f);
  std::vector<float> mesh_masses(1, ball_mass);

  gpu_sys.LoadMeshes(mesh_filenames, mesh_rotscales, mesh_translations,
                     mesh_masses);

  gpu_sys.SetOutputMode(params.write_mode);
  gpu_sys.SetVerbosity(params.verbose);

  std::string out_dir = "../";
  filesystem::create_directory(filesystem::path(out_dir));
  out_dir = out_dir + params.output_dir;
  filesystem::create_directory(filesystem::path(out_dir));

  std::cout << gpu_sys.GetNumMeshes() << " meshes" << std::endl;

  gpu_sys.Initialize();

  // Create rigid ball_body simulation
  ChSystemSMC sys_ball;
  sys_ball.SetContactForceModel(ChSystemSMC::ContactForceModel::Hooke);
  sys_ball.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
  sys_ball.Set_G_acc(ChVector<>(0, 0, -980));

  double inertia = 2.0 / 5.0 * ball_mass * ball_radius * ball_radius;
  ChVector<> ball_initial_pos(
      0, 0, fill_top + ball_radius + 2 * params.sphere_radius);

  std::shared_ptr<ChBody> ball_body(sys_ball.NewBody());
  ball_body->SetMass(ball_mass);
  ball_body->SetInertiaXX(ChVector<>(inertia, inertia, inertia));
  ball_body->SetPos(ball_initial_pos);
  auto sph = chrono_types::make_shared<ChSphereShape>();
  sph->GetSphereGeometry().rad = ball_radius;
  ball_body->AddAsset(sph);
  sys_ball.AddBody(ball_body);

  ChGpuVisualization gpu_vis(&gpu_sys, &sys_ball);
  if (render) {
    gpu_vis.SetTitle("Chrono::Gpu ball cosim demo");
    gpu_vis.SetCameraPosition(ChVector<>(0, -200, 100), ChVector<>(0, 0, 0));
    gpu_vis.SetCameraMoveScale(1.0f);
    gpu_vis.Initialize();
  }

  std::cout << "Output at    " << out_fps << " FPS" << std::endl;
  std::cout << "Rendering at " << render_fps << " FPS" << std::endl;
  unsigned int out_steps = (unsigned int)(1.0 / (out_fps * iteration_step));
  unsigned int render_steps =
      (unsigned int)(1.0 / (render_fps * iteration_step));
  unsigned int total_frames = (unsigned int)(params.time_end * out_fps);

  int currframe = 0;
  unsigned int curr_step = 0;

  clock_t start = std::clock();
  for (double t = 0; t < (double)params.time_end;
       t += iteration_step, curr_step++) {
    gpu_sys.ApplyMeshMotion(0, ball_body->GetPos(), ball_body->GetRot(),
                            ball_body->GetPos_dt(), ball_body->GetWvel_par());

    ChVector<> ball_force;
    ChVector<> ball_torque;
    gpu_sys.CollectMeshContactForces(0, ball_force, ball_torque);

    ball_body->Empty_forces_accumulators();
    ball_body->Accumulate_force(ball_force, ball_body->GetPos(), false);
    ball_body->Accumulate_torque(ball_torque, false);

    if (curr_step % out_steps == 0) {
      std::cout << "Output frame " << currframe + 1 << " of " << total_frames
                << std::endl;
      char filename[100];
      sprintf(filename, "%s/step%06d", out_dir.c_str(), currframe++);
      gpu_sys.WriteFile(std::string(filename));
      gpu_sys.WriteMeshes(std::string(filename));

      /*  // disable meshframes output, for it may be confusing for users
      dealing with Chrono::Gpu only std::string mesh_output =
      std::string(filename) + "_meshframes.csv"; std::ofstream
      meshfile(mesh_output); std::ostringstream outstream; outstream <<
      "mesh_name,dx,dy,dz,x1,x2,x3,y1,y2,y3,z1,z2,z3,sx,sy,sz\n";
      writeMeshFrames(outstream, *ball_body, mesh_filename, ball_radius);
      meshfile << outstream.str();
      */
    }

    if (render && curr_step % render_steps == 0) {
      if (gpu_vis.Render())
        break;
    }

    gpu_sys.AdvanceSimulation(iteration_step);
    sys_ball.DoStepDynamics(iteration_step);
  }

  clock_t end = std::clock();
  double total_time = ((double)(end - start)) / CLOCKS_PER_SEC;

  std::cout << "Time: " << total_time << " seconds" << std::endl;

  return 0;
}