1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2019 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Nic Olsen, Ruochun Zhang, Dan Negrut, Radu Serban
13 // =============================================================================
14
15 #include <string>
16 #include <cmath>
17
18 #include "chrono_gpu/physics/ChSystemGpu.h"
19 #include "chrono_gpu/physics/ChSystemGpu_impl.h"
20 #include "chrono_gpu/physics/ChSystemGpuMesh_impl.h"
21
22 #include "chrono_gpu/utils/ChGpuUtilities.h"
23
24 namespace chrono {
25 namespace gpu {
26
27 // -----------------------------------------------------------------------------
28
ChSystemGpu(float sphere_rad,float density,const ChVector<float> & boxDims,ChVector<float> O)29 ChSystemGpu::ChSystemGpu(float sphere_rad, float density, const ChVector<float>& boxDims, ChVector<float> O) {
30 m_sys = new ChSystemGpu_impl(sphere_rad, density, make_float3(boxDims.x(), boxDims.y(), boxDims.z()),
31 make_float3(O.x(), O.y(), O.z()));
32 }
33
ChSystemGpu(const std::string & checkpoint)34 ChSystemGpu::ChSystemGpu(const std::string& checkpoint) {
35 m_sys = new ChSystemGpu_impl(1.f, 1.f, make_float3(100, 100, 100), make_float3(0, 0, 0));
36 ReadCheckpointFile(checkpoint, true);
37 }
38
ChSystemGpuMesh(float sphere_rad,float density,const ChVector<float> & boxDims,ChVector<float> O)39 ChSystemGpuMesh::ChSystemGpuMesh(float sphere_rad, float density, const ChVector<float>& boxDims, ChVector<float> O)
40 : mesh_verbosity(CHGPU_MESH_VERBOSITY::QUIET) {
41 m_sys = new ChSystemGpuMesh_impl(sphere_rad, density, make_float3(boxDims.x(), boxDims.y(), boxDims.z()),
42 make_float3(O.x(), O.y(), O.z()));
43 }
44
ChSystemGpuMesh(const std::string & checkpoint)45 ChSystemGpuMesh::ChSystemGpuMesh(const std::string& checkpoint) : mesh_verbosity(CHGPU_MESH_VERBOSITY::QUIET) {
46 m_sys = new ChSystemGpuMesh_impl(1.f, 1.f, make_float3(100, 100, 100), make_float3(0, 0, 0));
47 ReadCheckpointFile(checkpoint, true);
48 }
49
~ChSystemGpu()50 ChSystemGpu::~ChSystemGpu() {
51 delete m_sys;
52 }
53
~ChSystemGpuMesh()54 ChSystemGpuMesh::~ChSystemGpuMesh() {}
55
56 // -----------------------------------------------------------------------------
57
SetGravitationalAcceleration(const ChVector<float> & g)58 void ChSystemGpu::SetGravitationalAcceleration(const ChVector<float>& g) {
59 m_sys->X_accGrav = (float)g.x();
60 m_sys->Y_accGrav = (float)g.y();
61 m_sys->Z_accGrav = (float)g.z();
62 }
63
SetGravitationalAcceleration(const float3 g)64 void ChSystemGpu::SetGravitationalAcceleration(const float3 g) {
65 m_sys->X_accGrav = g.x;
66 m_sys->Y_accGrav = g.y;
67 m_sys->Z_accGrav = g.z;
68 }
69
SetBDFixed(bool fixed)70 void ChSystemGpu::SetBDFixed(bool fixed) {
71 m_sys->BD_is_fixed = fixed;
72 }
73
SetBDCenter(const ChVector<float> & O)74 void ChSystemGpu::SetBDCenter(const ChVector<float>& O) {
75 m_sys->user_coord_O_X = O.x();
76 m_sys->user_coord_O_Y = O.y();
77 m_sys->user_coord_O_Z = O.z();
78 }
79
SetParticleFixed(const std::vector<bool> & fixed)80 void ChSystemGpu::SetParticleFixed(const std::vector<bool>& fixed) {
81 m_sys->user_sphere_fixed = fixed;
82 }
83
84 // Set particle output file format
SetParticleOutputMode(CHGPU_OUTPUT_MODE mode)85 void ChSystemGpu::SetParticleOutputMode(CHGPU_OUTPUT_MODE mode) {
86 m_sys->file_write_mode = mode;
87 }
88
89 // Set particle output file content
SetParticleOutputFlags(unsigned int flags)90 void ChSystemGpu::SetParticleOutputFlags(unsigned int flags) {
91 m_sys->output_flags = flags;
92 }
93
SetFixedStepSize(float size_UU)94 void ChSystemGpu::SetFixedStepSize(float size_UU) {
95 m_sys->stepSize_UU = size_UU;
96 }
97
EnableMinLength(bool useMinLen)98 void ChSystemGpu::EnableMinLength(bool useMinLen) {
99 m_sys->use_min_length_unit = useMinLen;
100 }
101
SetTimeIntegrator(CHGPU_TIME_INTEGRATOR new_integrator)102 void ChSystemGpu::SetTimeIntegrator(CHGPU_TIME_INTEGRATOR new_integrator) {
103 m_sys->gran_params->time_integrator = new_integrator;
104 m_sys->time_integrator = new_integrator;
105 }
106
SetFrictionMode(CHGPU_FRICTION_MODE new_mode)107 void ChSystemGpu::SetFrictionMode(CHGPU_FRICTION_MODE new_mode) {
108 m_sys->gran_params->friction_mode = new_mode;
109 }
110
SetRollingMode(CHGPU_ROLLING_MODE new_mode)111 void ChSystemGpu::SetRollingMode(CHGPU_ROLLING_MODE new_mode) {
112 m_sys->gran_params->rolling_mode = new_mode;
113 }
114
SetDefragmentOnInitialize(bool defragment)115 void ChSystemGpu::SetDefragmentOnInitialize(bool defragment) {
116 m_sys->defragment_on_start = defragment;
117 }
118
SetRecordingContactInfo(bool record)119 void ChSystemGpu::SetRecordingContactInfo(bool record) {
120 m_sys->gran_params->recording_contactInfo = record;
121 }
122
SetMaxSafeVelocity_SU(float max_vel)123 void ChSystemGpu::SetMaxSafeVelocity_SU(float max_vel) {
124 m_sys->gran_params->max_safe_vel = max_vel;
125 }
126
SetPsiFactors(unsigned int psi_T,unsigned int psi_L,float psi_R)127 void ChSystemGpu::SetPsiFactors(unsigned int psi_T, unsigned int psi_L, float psi_R) {
128 m_sys->psi_T = psi_T;
129 m_sys->psi_L = psi_L;
130 m_sys->psi_R = psi_R;
131 }
132
SetPsiT(unsigned int psi_T)133 void ChSystemGpu::SetPsiT(unsigned int psi_T) {
134 m_sys->psi_T = psi_T;
135 }
136
SetPsiL(unsigned int psi_L)137 void ChSystemGpu::SetPsiL(unsigned int psi_L) {
138 m_sys->psi_L = psi_L;
139 }
140
SetPsiR(float psi_R)141 void ChSystemGpu::SetPsiR(float psi_R) {
142 m_sys->psi_R = psi_R;
143 }
144
SetSimTime(float time)145 void ChSystemGpu::SetSimTime(float time) {
146 m_sys->elapsedSimTime = time;
147 }
148
149 // -----------------------------------------------------------------------------
150
SetStaticFrictionCoeff_SPH2SPH(float mu)151 void ChSystemGpu::SetStaticFrictionCoeff_SPH2SPH(float mu) {
152 m_sys->gran_params->static_friction_coeff_s2s = mu;
153 }
154
SetStaticFrictionCoeff_SPH2WALL(float mu)155 void ChSystemGpu::SetStaticFrictionCoeff_SPH2WALL(float mu) {
156 m_sys->gran_params->static_friction_coeff_s2w = mu;
157 }
158
SetRollingCoeff_SPH2SPH(float mu)159 void ChSystemGpu::SetRollingCoeff_SPH2SPH(float mu) {
160 m_sys->rolling_coeff_s2s_UU = mu;
161 }
162
SetRollingCoeff_SPH2WALL(float mu)163 void ChSystemGpu::SetRollingCoeff_SPH2WALL(float mu) {
164 m_sys->rolling_coeff_s2w_UU = mu;
165 }
166
SetSpinningCoeff_SPH2SPH(float mu)167 void ChSystemGpu::SetSpinningCoeff_SPH2SPH(float mu) {
168 m_sys->spinning_coeff_s2s_UU = mu;
169 }
170
SetSpinningCoeff_SPH2WALL(float mu)171 void ChSystemGpu::SetSpinningCoeff_SPH2WALL(float mu) {
172 m_sys->spinning_coeff_s2w_UU = mu;
173 }
174
SetKn_SPH2SPH(double someValue)175 void ChSystemGpu::SetKn_SPH2SPH(double someValue) {
176 m_sys->K_n_s2s_UU = someValue;
177 }
178
SetKn_SPH2WALL(double someValue)179 void ChSystemGpu::SetKn_SPH2WALL(double someValue) {
180 m_sys->K_n_s2w_UU = someValue;
181 }
182
SetGn_SPH2SPH(double someValue)183 void ChSystemGpu::SetGn_SPH2SPH(double someValue) {
184 m_sys->Gamma_n_s2s_UU = someValue;
185 }
186
SetGn_SPH2WALL(double someValue)187 void ChSystemGpu::SetGn_SPH2WALL(double someValue) {
188 m_sys->Gamma_n_s2w_UU = someValue;
189 }
190
SetKt_SPH2SPH(double someValue)191 void ChSystemGpu::SetKt_SPH2SPH(double someValue) {
192 m_sys->K_t_s2s_UU = someValue;
193 }
194
SetGt_SPH2SPH(double someValue)195 void ChSystemGpu::SetGt_SPH2SPH(double someValue) {
196 m_sys->Gamma_t_s2s_UU = someValue;
197 }
198
SetKt_SPH2WALL(double someValue)199 void ChSystemGpu::SetKt_SPH2WALL(double someValue) {
200 m_sys->K_t_s2w_UU = someValue;
201 }
202
SetGt_SPH2WALL(double someValue)203 void ChSystemGpu::SetGt_SPH2WALL(double someValue) {
204 m_sys->Gamma_t_s2w_UU = someValue;
205 }
206
SetCohesionRatio(float someValue)207 void ChSystemGpu::SetCohesionRatio(float someValue) {
208 m_sys->cohesion_over_gravity = someValue;
209 }
210
SetAdhesionRatio_SPH2WALL(float someValue)211 void ChSystemGpu::SetAdhesionRatio_SPH2WALL(float someValue) {
212 m_sys->adhesion_s2w_over_gravity = someValue;
213 }
214
UseMaterialBasedModel(bool val)215 void ChSystemGpu::UseMaterialBasedModel(bool val) {
216 m_sys->use_mat_based = val;
217 m_sys->gran_params->use_mat_based = val;
218 }
219
SetYoungModulus_SPH(double someValue)220 void ChSystemGpu::SetYoungModulus_SPH(double someValue) {
221 m_sys->YoungsModulus_sphere_UU = someValue;
222 }
223
SetYoungModulus_WALL(double someValue)224 void ChSystemGpu::SetYoungModulus_WALL(double someValue) {
225 m_sys->YoungsModulus_wall_UU = someValue;
226 }
227
SetPoissonRatio_SPH(double someValue)228 void ChSystemGpu::SetPoissonRatio_SPH(double someValue) {
229 m_sys->PoissonRatio_sphere_UU = someValue;
230 }
231
SetPoissonRatio_WALL(double someValue)232 void ChSystemGpu::SetPoissonRatio_WALL(double someValue) {
233 m_sys->PoissonRatio_wall_UU = someValue;
234 }
235
SetRestitution_SPH(double someValue)236 void ChSystemGpu::SetRestitution_SPH(double someValue) {
237 m_sys->COR_sphere_UU = someValue;
238 }
239
SetRestitution_WALL(double someValue)240 void ChSystemGpu::SetRestitution_WALL(double someValue) {
241 m_sys->COR_wall_UU = someValue;
242 }
243
244 // -----------------------------------------------------------------------------
245
UseMaterialBasedModel(bool val)246 void ChSystemGpuMesh::UseMaterialBasedModel(bool val) {
247 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
248 sys_trimesh->tri_params->use_mat_based = val;
249 sys_trimesh->gran_params->use_mat_based = val;
250 sys_trimesh->use_mat_based = val;
251 }
252
SetStaticFrictionCoeff_SPH2MESH(float mu)253 void ChSystemGpuMesh::SetStaticFrictionCoeff_SPH2MESH(float mu) {
254 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
255 sys_trimesh->tri_params->static_friction_coeff_s2m = mu;
256 }
257
SetRollingCoeff_SPH2MESH(float mu)258 void ChSystemGpuMesh::SetRollingCoeff_SPH2MESH(float mu) {
259 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
260 sys_trimesh->rolling_coeff_s2m_UU = mu;
261 }
262
SetSpinningCoeff_SPH2MESH(float mu)263 void ChSystemGpuMesh::SetSpinningCoeff_SPH2MESH(float mu) {
264 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
265 sys_trimesh->spinning_coeff_s2m_UU = mu;
266 }
267
SetKn_SPH2MESH(double someValue)268 void ChSystemGpuMesh::SetKn_SPH2MESH(double someValue) {
269 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
270 sys_trimesh->K_n_s2m_UU = someValue;
271 }
272
SetGn_SPH2MESH(double someValue)273 void ChSystemGpuMesh::SetGn_SPH2MESH(double someValue) {
274 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
275 sys_trimesh->Gamma_n_s2m_UU = someValue;
276 }
277
SetKt_SPH2MESH(double someValue)278 void ChSystemGpuMesh::SetKt_SPH2MESH(double someValue) {
279 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
280 sys_trimesh->K_t_s2m_UU = someValue;
281 }
282
SetGt_SPH2MESH(double someValue)283 void ChSystemGpuMesh::SetGt_SPH2MESH(double someValue) {
284 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
285 sys_trimesh->Gamma_t_s2m_UU = someValue;
286 }
287
SetAdhesionRatio_SPH2MESH(float someValue)288 void ChSystemGpuMesh::SetAdhesionRatio_SPH2MESH(float someValue) {
289 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
290 sys_trimesh->adhesion_s2m_over_gravity = someValue;
291 }
292
SetYoungModulus_MESH(double someValue)293 void ChSystemGpuMesh::SetYoungModulus_MESH(double someValue) {
294 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
295 sys_trimesh->YoungsModulus_mesh_UU = someValue;
296 }
297
SetPoissonRatio_MESH(double someValue)298 void ChSystemGpuMesh::SetPoissonRatio_MESH(double someValue) {
299 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
300 sys_trimesh->PoissonRatio_mesh_UU = someValue;
301 }
302
SetRestitution_MESH(double someValue)303 void ChSystemGpuMesh::SetRestitution_MESH(double someValue) {
304 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
305 sys_trimesh->COR_mesh_UU = someValue;
306 }
307
308 // -----------------------------------------------------------------------------
309
SetVerbosity(CHGPU_VERBOSITY level)310 void ChSystemGpu::SetVerbosity(CHGPU_VERBOSITY level) {
311 m_sys->verbosity = level;
312 }
313
SetMeshVerbosity(CHGPU_MESH_VERBOSITY level)314 void ChSystemGpuMesh::SetMeshVerbosity(CHGPU_MESH_VERBOSITY level) {
315 mesh_verbosity = level;
316 }
317
318 // -----------------------------------------------------------------------------
319
CreateBCSphere(const ChVector<float> & center,float radius,bool outward_normal,bool track_forces,float sphere_mass)320 size_t ChSystemGpu::CreateBCSphere(const ChVector<float>& center,
321 float radius,
322 bool outward_normal,
323 bool track_forces,
324 float sphere_mass) {
325 float sph_center[3] = {center.x(), center.y(), center.z()};
326 return m_sys->CreateBCSphere(sph_center, radius, outward_normal, track_forces, sphere_mass);
327 }
328
SetBCSpherePosition(size_t sphere_bc_id,const ChVector<float> & pos)329 void ChSystemGpu::SetBCSpherePosition(size_t sphere_bc_id, const ChVector<float>& pos) {
330 m_sys->SetBCSpherePosition(sphere_bc_id, make_float3(pos.x(), pos.y(), pos.z()));
331 }
332
GetBCSpherePosition(size_t sphere_bc_id) const333 ChVector<float> ChSystemGpu::GetBCSpherePosition(size_t sphere_bc_id) const {
334 float3 pos = m_sys->GetBCSpherePosition(sphere_bc_id);
335 return ChVector<float>(pos.x, pos.y, pos.z);
336 }
337
SetBCSphereVelocity(size_t sphere_bc_id,const ChVector<float> & velo)338 void ChSystemGpu::SetBCSphereVelocity(size_t sphere_bc_id, const ChVector<float>& velo) {
339 m_sys->SetBCSphereVelocity(sphere_bc_id, make_float3(velo.x(), velo.y(), velo.z()));
340 }
341
GetBCSphereVelocity(size_t sphere_bc_id) const342 ChVector<float> ChSystemGpu::GetBCSphereVelocity(size_t sphere_bc_id) const {
343 float3 velo = m_sys->GetBCSphereVelocity(sphere_bc_id);
344 return ChVector<float>(velo.x, velo.y, velo.z);
345 }
346
CreateBCConeZ(const ChVector<float> & tip,float slope,float hmax,float hmin,bool outward_normal,bool track_forces)347 size_t ChSystemGpu::CreateBCConeZ(const ChVector<float>& tip,
348 float slope,
349 float hmax,
350 float hmin,
351 bool outward_normal,
352 bool track_forces) {
353 float cone_tip[3] = {tip.x(), tip.y(), tip.z()};
354 return m_sys->CreateBCConeZ(cone_tip, slope, hmax, hmin, outward_normal, track_forces);
355 }
356
CreateBCPlane(const ChVector<float> & pos,const ChVector<float> & normal,bool track_forces)357 size_t ChSystemGpu::CreateBCPlane(const ChVector<float>& pos, const ChVector<float>& normal, bool track_forces) {
358 float plane_pos[3] = {pos.x(), pos.y(), pos.z()};
359 float plane_nrm[3] = {normal.x(), normal.y(), normal.z()};
360
361 return m_sys->CreateBCPlane(plane_pos, plane_nrm, track_forces);
362 }
363
364 // customized plate for angle of repose test, plate has a limited width in y direction, and can move in y direction
CreateCustomizedPlate(const ChVector<float> & pos_center,const ChVector<float> & normal,float hdim_y)365 size_t ChSystemGpu::CreateCustomizedPlate(const ChVector<float>& pos_center,
366 const ChVector<float>& normal,
367 float hdim_y) {
368 float plate_pos[3] = {pos_center.x(), pos_center.y(), pos_center.z()};
369 float plate_nrm[3] = {normal.x(), normal.y(), normal.z()};
370
371 return m_sys->CreateCustomizedPlate(plate_pos, plate_nrm, hdim_y);
372 }
373
CreateBCCylinderZ(const ChVector<float> & center,float radius,bool outward_normal,bool track_forces)374 size_t ChSystemGpu::CreateBCCylinderZ(const ChVector<float>& center,
375 float radius,
376 bool outward_normal,
377 bool track_forces) {
378 float cyl_center[3] = {center.x(), center.y(), center.z()};
379 return m_sys->CreateBCCylinderZ(cyl_center, radius, outward_normal, track_forces);
380 }
381
DisableBCbyID(size_t BC_id)382 bool ChSystemGpu::DisableBCbyID(size_t BC_id) {
383 return m_sys->DisableBCbyID(BC_id);
384 }
385
EnableBCbyID(size_t BC_id)386 bool ChSystemGpu::EnableBCbyID(size_t BC_id) {
387 return m_sys->EnableBCbyID(BC_id);
388 }
389
SetBCOffsetFunction(size_t BC_id,const GranPositionFunction & offset_function)390 bool ChSystemGpu::SetBCOffsetFunction(size_t BC_id, const GranPositionFunction& offset_function) {
391 return m_sys->SetBCOffsetFunction(BC_id, offset_function);
392 }
393
setBDWallsMotionFunction(const GranPositionFunction & pos_fn)394 void ChSystemGpu::setBDWallsMotionFunction(const GranPositionFunction& pos_fn) {
395 m_sys->setBDWallsMotionFunction(pos_fn);
396 }
397
398 // -----------------------------------------------------------------------------
399
GetNumSDs() const400 unsigned int ChSystemGpu::GetNumSDs() const {
401 return m_sys->nSDs;
402 }
403
GetBCPlanePosition(size_t plane_id) const404 ChVector<float> ChSystemGpu::GetBCPlanePosition(size_t plane_id) const {
405 // todo: throw an error if BC not a plane type
406 float3 pos = m_sys->GetBCPlanePosition(plane_id);
407 return ChVector<float>(pos.x, pos.y, pos.z);
408 }
409
SetBCPlaneRotation(size_t plane_id,ChVector<double> center,ChVector<double> omega)410 void ChSystemGpu::SetBCPlaneRotation(size_t plane_id, ChVector<double> center, ChVector<double> omega) {
411 double3 rotation_center = make_double3(center.x(), center.y(), center.z());
412 double3 rotation_omega = make_double3(omega.x(), omega.y(), omega.z());
413 m_sys->SetBCPlaneRotation(plane_id, rotation_center, rotation_omega);
414 }
415
GetBCReactionForces(size_t BC_id,ChVector<float> & force) const416 bool ChSystemGpu::GetBCReactionForces(size_t BC_id, ChVector<float>& force) const {
417 float3 frc;
418 bool ret = m_sys->GetBCReactionForces(BC_id, frc);
419 force = ChVector<float>(frc.x, frc.y, frc.z);
420 return ret;
421 }
422
GetNumContacts() const423 int ChSystemGpu::GetNumContacts() const {
424 return m_sys->GetNumContacts();
425 }
426
GetSimTime() const427 float ChSystemGpu::GetSimTime() const {
428 return m_sys->elapsedSimTime;
429 }
430
GetNumParticles() const431 size_t ChSystemGpu::GetNumParticles() const {
432 return m_sys->nSpheres;
433 }
434
GetParticleRadius() const435 float ChSystemGpu::GetParticleRadius() const {
436 return m_sys->sphere_radius_UU;
437 }
438
GetParticlePosition(int nSphere) const439 ChVector<float> ChSystemGpu::GetParticlePosition(int nSphere) const {
440 float3 pos = m_sys->GetParticlePosition(nSphere);
441 return ChVector<float>(pos.x, pos.y, pos.z);
442 }
443
SetParticlePosition(int nSphere,const ChVector<double> pos)444 void ChSystemGpu::SetParticlePosition(int nSphere, const ChVector<double> pos) {
445 double3 position = make_double3(pos.x(), pos.y(), pos.z());
446 m_sys->SetParticlePosition(nSphere, position);
447 }
448
GetParticleVelocity(int nSphere) const449 ChVector<float> ChSystemGpu::GetParticleVelocity(int nSphere) const {
450 float3 vel = m_sys->GetParticleLinVelocity(nSphere);
451 return ChVector<float>(vel.x, vel.y, vel.z);
452 }
453
SetParticleVelocity(int nSphere,const ChVector<double> velo)454 void ChSystemGpu::SetParticleVelocity(int nSphere, const ChVector<double> velo) {
455 double3 velocity = make_double3(velo.x(), velo.y(), velo.z());
456 m_sys->SetParticleVelocity(nSphere, velocity);
457 }
458
GetParticleAngVelocity(int nSphere) const459 ChVector<float> ChSystemGpu::GetParticleAngVelocity(int nSphere) const {
460 if (m_sys->gran_params->friction_mode == CHGPU_FRICTION_MODE::FRICTIONLESS)
461 return ChVector<float>(0);
462
463 float3 omega = m_sys->GetParticleAngVelocity(nSphere);
464 return ChVector<float>(omega.x, omega.y, omega.z);
465 }
466
GetParticleLinAcc(int nSphere) const467 ChVector<float> ChSystemGpu::GetParticleLinAcc(int nSphere) const {
468 float3 acc = m_sys->GetParticleLinAcc(nSphere);
469 return ChVector<float>(acc.x, acc.y, acc.z);
470 }
471
IsFixed(int nSphere) const472 bool ChSystemGpu::IsFixed(int nSphere) const {
473 return m_sys->IsFixed(nSphere);
474 }
475
GetParticlesKineticEnergy() const476 float ChSystemGpu::GetParticlesKineticEnergy() const {
477 float KE = (float)(m_sys->ComputeTotalKE());
478 return KE;
479 }
480
GetMaxParticleZ() const481 double ChSystemGpu::GetMaxParticleZ() const {
482 return m_sys->GetMaxParticleZ(true);
483 }
484
GetMinParticleZ() const485 double ChSystemGpu::GetMinParticleZ() const {
486 // Under the hood, GetMaxParticleZ(false) returns the lowest Z.
487 return m_sys->GetMaxParticleZ(false);
488 }
489
EstimateMemUsage() const490 size_t ChSystemGpu::EstimateMemUsage() const {
491 return m_sys->EstimateMemUsage();
492 }
493
494 // -----------------------------------------------------------------------------
495
AddMesh(std::shared_ptr<geometry::ChTriangleMeshConnected> mesh,float mass)496 unsigned int ChSystemGpuMesh::AddMesh(std::shared_ptr<geometry::ChTriangleMeshConnected> mesh, float mass) {
497 unsigned int id = static_cast<unsigned int>(m_meshes.size());
498 m_meshes.push_back(mesh);
499 m_mesh_masses.push_back(mass);
500
501 return id;
502 }
503
AddMesh(const std::string & filename,const ChVector<float> & translation,const ChMatrix33<float> & rotscale,float mass)504 unsigned int ChSystemGpuMesh::AddMesh(const std::string& filename,
505 const ChVector<float>& translation,
506 const ChMatrix33<float>& rotscale,
507 float mass) {
508 auto mesh = chrono_types::make_shared<geometry::ChTriangleMeshConnected>();
509 bool flag = mesh->LoadWavefrontMesh(filename, true, false);
510 if (!flag)
511 CHGPU_ERROR("ERROR! Mesh %s failed to load in!\n", filename.c_str());
512 if (mesh->getNumTriangles() == 0)
513 printf("WARNING: Mesh %s has no triangles!\n", filename.c_str());
514 mesh->Transform(translation, rotscale.cast<double>());
515
516 unsigned int id = static_cast<unsigned int>(m_meshes.size());
517 m_meshes.push_back(mesh);
518 m_mesh_masses.push_back(mass);
519
520 return id;
521 }
522
AddMeshes(const std::vector<std::string> & objfilenames,const std::vector<ChVector<float>> & translations,const std::vector<ChMatrix33<float>> & rotscales,const std::vector<float> & masses)523 std::vector<unsigned int> ChSystemGpuMesh::AddMeshes(const std::vector<std::string>& objfilenames,
524 const std::vector<ChVector<float>>& translations,
525 const std::vector<ChMatrix33<float>>& rotscales,
526 const std::vector<float>& masses) {
527 unsigned int size = (unsigned int)objfilenames.size();
528 if (size != rotscales.size() || size != translations.size() || size != masses.size())
529 CHGPU_ERROR("ERROR! Mesh loading vectors must all have same size!\n");
530 if (size == 0)
531 printf("WARNING: No meshes provided!\n");
532
533 std::vector<unsigned int> ids(size);
534 for (unsigned int i = 0; i < size; i++) {
535 ids[i] = AddMesh(objfilenames[i], translations[i], rotscales[i], masses[i]);
536 }
537
538 return ids;
539 }
540
SetMeshes()541 void ChSystemGpuMesh::SetMeshes() {
542 int nTriangles = 0;
543 for (const auto& mesh : m_meshes)
544 nTriangles += mesh->getNumTriangles();
545
546 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
547 ChSystemGpuMesh_impl::TriangleSoup* pMeshSoup = sys_trimesh->getMeshSoup();
548 pMeshSoup->nTrianglesInSoup = nTriangles;
549 if (nTriangles != 0) {
550 // Allocate all of the requisite pointers
551 gpuErrchk(
552 cudaMallocManaged(&pMeshSoup->triangleFamily_ID, nTriangles * sizeof(unsigned int), cudaMemAttachGlobal));
553
554 gpuErrchk(cudaMallocManaged(&pMeshSoup->node1, nTriangles * sizeof(float3), cudaMemAttachGlobal));
555 gpuErrchk(cudaMallocManaged(&pMeshSoup->node2, nTriangles * sizeof(float3), cudaMemAttachGlobal));
556 gpuErrchk(cudaMallocManaged(&pMeshSoup->node3, nTriangles * sizeof(float3), cudaMemAttachGlobal));
557 }
558
559 MESH_INFO_PRINTF("Done allocating nodes for %d triangles\n", nTriangles);
560
561 // Setup the clean copy of the mesh soup from the obj file data
562 unsigned int family = 0;
563 unsigned int tri_i = 0;
564 // for each obj file data set
565 for (const auto& mesh : m_meshes) {
566 for (int i = 0; i < mesh->getNumTriangles(); i++) {
567 geometry::ChTriangle tri = mesh->getTriangle(i);
568
569 pMeshSoup->node1[tri_i] = make_float3((float)tri.p1.x(), (float)tri.p1.y(), (float)tri.p1.z());
570 pMeshSoup->node2[tri_i] = make_float3((float)tri.p2.x(), (float)tri.p2.y(), (float)tri.p2.z());
571 pMeshSoup->node3[tri_i] = make_float3((float)tri.p3.x(), (float)tri.p3.y(), (float)tri.p3.z());
572
573 pMeshSoup->triangleFamily_ID[tri_i] = family;
574
575 // If we wish to correct surface orientation based on given vertex normals, rather than using RHR...
576 if (use_mesh_normals) {
577 int normal_i = mesh->m_face_n_indices.at(i).x(); // normals at each vertex of this triangle
578 ChVector<double> normal = mesh->m_normals.at(normal_i);
579
580 // Generate normal using RHR from nodes 1, 2, and 3
581 ChVector<double> AB = tri.p2 - tri.p1;
582 ChVector<double> AC = tri.p3 - tri.p1;
583 ChVector<double> cross;
584 cross.Cross(AB, AC);
585
586 // If the normal created by a RHR traversal is not correct, switch two vertices
587 if (cross.Dot(normal) < 0) {
588 std::swap(pMeshSoup->node2[tri_i], pMeshSoup->node3[tri_i]);
589 }
590 }
591
592 tri_i++;
593 }
594 family++;
595 MESH_INFO_PRINTF("Done writing family %d\n", family);
596 }
597
598 pMeshSoup->numTriangleFamilies = family;
599
600 if (pMeshSoup->nTrianglesInSoup != 0) {
601 gpuErrchk(cudaMallocManaged(&pMeshSoup->familyMass_SU, family * sizeof(float), cudaMemAttachGlobal));
602
603 for (unsigned int i = 0; i < family; i++) {
604 // NOTE The SU conversion is done in initialize after the scaling is determined
605 pMeshSoup->familyMass_SU[i] = m_mesh_masses[i];
606 }
607
608 gpuErrchk(cudaMallocManaged(&pMeshSoup->generalizedForcesPerFamily,
609 6 * pMeshSoup->numTriangleFamilies * sizeof(float), cudaMemAttachGlobal));
610 // Allocate memory for the float and double frames
611 gpuErrchk(cudaMallocManaged(&sys_trimesh->getTriParams()->fam_frame_broad,
612 pMeshSoup->numTriangleFamilies * sizeof(ChSystemGpuMesh_impl::MeshFrame<float>),
613 cudaMemAttachGlobal));
614 gpuErrchk(cudaMallocManaged(&sys_trimesh->getTriParams()->fam_frame_narrow,
615 pMeshSoup->numTriangleFamilies * sizeof(ChSystemGpuMesh_impl::MeshFrame<double>),
616 cudaMemAttachGlobal));
617
618 // Allocate memory for linear and angular velocity
619 gpuErrchk(
620 cudaMallocManaged(&pMeshSoup->vel, pMeshSoup->numTriangleFamilies * sizeof(float3), cudaMemAttachGlobal));
621 gpuErrchk(
622 cudaMallocManaged(&pMeshSoup->omega, pMeshSoup->numTriangleFamilies * sizeof(float3), cudaMemAttachGlobal));
623
624 for (unsigned int i = 0; i < family; i++) {
625 pMeshSoup->vel[i] = make_float3(0, 0, 0);
626 pMeshSoup->omega[i] = make_float3(0, 0, 0);
627 }
628 }
629 }
630
ApplyMeshMotion(unsigned int mesh_id,const ChVector<> & pos,const ChQuaternion<> & rot,const ChVector<> & lin_vel,const ChVector<> & ang_vel)631 void ChSystemGpuMesh::ApplyMeshMotion(unsigned int mesh_id,
632 const ChVector<>& pos,
633 const ChQuaternion<>& rot,
634 const ChVector<>& lin_vel,
635 const ChVector<>& ang_vel) {
636 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
637 sys_trimesh->ApplyMeshMotion(mesh_id, pos.data(), rot.data(), lin_vel.data(), ang_vel.data());
638 }
639
640 // -----------------------------------------------------------------------------
641
GetNumMeshes() const642 unsigned int ChSystemGpuMesh::GetNumMeshes() const {
643 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
644 return sys_trimesh->meshSoup->numTriangleFamilies;
645 }
646
EnableMeshCollision(bool val)647 void ChSystemGpuMesh::EnableMeshCollision(bool val) {
648 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
649 sys_trimesh->mesh_collision_enabled = val;
650 }
651
652 // -----------------------------------------------------------------------------
653
convertChVector2Float3Vec(const std::vector<ChVector<float>> & points,std::vector<float3> & pointsFloat3)654 static void convertChVector2Float3Vec(const std::vector<ChVector<float>>& points, std::vector<float3>& pointsFloat3) {
655 size_t nPoints = points.size();
656 pointsFloat3.resize(nPoints);
657 for (size_t index = 0; index < nPoints; index++) {
658 pointsFloat3.at(index).x = points.at(index)[0];
659 pointsFloat3.at(index).y = points.at(index)[1];
660 pointsFloat3.at(index).z = points.at(index)[2];
661 }
662 }
663
664 // Initialize particle positions, velocity and angular velocity in user units
SetParticles(const std::vector<ChVector<float>> & points,const std::vector<ChVector<float>> & vels,const std::vector<ChVector<float>> & ang_vel)665 void ChSystemGpu::SetParticles(const std::vector<ChVector<float>>& points,
666 const std::vector<ChVector<float>>& vels,
667 const std::vector<ChVector<float>>& ang_vel) {
668 std::vector<float3> pointsFloat3;
669 std::vector<float3> velsFloat3;
670 std::vector<float3> angVelsFloat3;
671 convertChVector2Float3Vec(points, pointsFloat3);
672 convertChVector2Float3Vec(vels, velsFloat3);
673 convertChVector2Float3Vec(ang_vel, angVelsFloat3);
674 m_sys->SetParticles(pointsFloat3, velsFloat3, angVelsFloat3);
675 }
676
677 // Use the CSV header to determine its content for data parsing format
quarryCsvFormat(const std::string & line)678 inline unsigned int quarryCsvFormat(const std::string& line) {
679 unsigned int formatMode = 0;
680 // Use the header to determine the information being loaded in.
681 // Note that the order they appear is guaranteed to be like below
682 if (line.find("x,y,z") == std::string::npos)
683 CHGPU_ERROR("ERROR! Checkpoint file (kinematics) does not contain xyz information!\n");
684 if (line.find("vx,vy,vz") != std::string::npos)
685 formatMode += VEL_COMPONENTS;
686 if (line.find("absv") != std::string::npos)
687 formatMode += ABSV;
688 if (line.find("fixed") != std::string::npos)
689 formatMode += FIXITY;
690 if (line.find("wx,wy,wz") != std::string::npos)
691 formatMode += ANG_VEL_COMPONENTS;
692 return formatMode;
693 }
694
695 // User the history header to determine the friction/contact info we'll load
quarryHistoryFormat(const std::string & line,unsigned int max_partner)696 inline unsigned int quarryHistoryFormat(const std::string& line, unsigned int max_partner) {
697 unsigned int formatMode = 0;
698 std::istringstream iss1(line);
699 if (line.find("partners") != std::string::npos) {
700 formatMode += 1;
701 iss1 >> max_partner; // this one is a string... not the number we're after
702 iss1 >> max_partner;
703 } else
704 printf(
705 "WARNING! The friction history either has no header, or its header indicates there is no history to "
706 "load "
707 "in!\n");
708 if (line.find("history") != std::string::npos)
709 formatMode += 2;
710 return formatMode;
711 }
712
713 // Read in particle positions, velocity and angular velocity through a csv-formatted ifstream
ReadCsvParticles(std::ifstream & ifile,unsigned int totRow)714 void ChSystemGpu::ReadCsvParticles(std::ifstream& ifile, unsigned int totRow) {
715 // Read the header to know the input format
716 std::string line;
717 std::getline(ifile, line);
718 while (line.find_first_not_of(' ') == std::string::npos)
719 std::getline(ifile, line);
720 unsigned int formatMode = quarryCsvFormat(line);
721
722 // Parse in every data line
723 std::string number;
724 unsigned int numRow = 0;
725 std::vector<float3> pointsFloat3, velsFloat3, angVelsFloat3;
726 std::vector<bool> fixity;
727 // don't always assume we know the num of rows, and we expand buffer every time batchSize is hit
728 const unsigned int batchSize = 1e6;
729 while (std::getline(ifile, line)) {
730 // if empty line, just keep going
731 if (line.find_first_not_of(' ') == std::string::npos)
732 continue;
733
734 std::istringstream iss1(line);
735 if (numRow % batchSize == 0) {
736 pointsFloat3.resize(pointsFloat3.size() + batchSize);
737 velsFloat3.resize(velsFloat3.size() + batchSize, make_float3(0, 0, 0));
738 fixity.resize(fixity.size() + batchSize, false);
739 angVelsFloat3.resize(angVelsFloat3.size() + batchSize, make_float3(0, 0, 0));
740 }
741 // Read xyz
742 getline(iss1, number, ',');
743 pointsFloat3.at(numRow).x = std::stof(number);
744 getline(iss1, number, ',');
745 pointsFloat3.at(numRow).y = std::stof(number);
746 getline(iss1, number, ',');
747 pointsFloat3.at(numRow).z = std::stof(number);
748 // Read velocity
749 if (formatMode & VEL_COMPONENTS) {
750 getline(iss1, number, ',');
751 velsFloat3.at(numRow).x = std::stof(number);
752 getline(iss1, number, ',');
753 velsFloat3.at(numRow).y = std::stof(number);
754 getline(iss1, number, ',');
755 velsFloat3.at(numRow).z = std::stof(number);
756 }
757 // Read absv, but we don't need it
758 if (formatMode & ABSV)
759 getline(iss1, number, ',');
760 // Read fixity
761 if (formatMode & FIXITY) {
762 getline(iss1, number, ',');
763 fixity.at(numRow) = (bool)std::stoi(number);
764 }
765 // Read angular velocity
766 if (formatMode & ANG_VEL_COMPONENTS) {
767 getline(iss1, number, ',');
768 angVelsFloat3.at(numRow).x = std::stof(number);
769 getline(iss1, number, ',');
770 angVelsFloat3.at(numRow).y = std::stof(number);
771 getline(iss1, number, ',');
772 angVelsFloat3.at(numRow).z = std::stof(number);
773 }
774 numRow++;
775 // If all spheres loaded, break. Don't keep going as it's not necessarily the EoF.
776 if (numRow >= totRow)
777 break;
778 }
779
780 // Final resize
781 pointsFloat3.resize(numRow);
782 velsFloat3.resize(numRow);
783 angVelsFloat3.resize(numRow);
784 fixity.resize(numRow);
785
786 // Feed data to the system
787 // We are directly using low-level m_sys functions here, since we're already working with float3
788 // But we can choose to call user panel SetParticles and SetParticleFixed as well
789 m_sys->SetParticles(pointsFloat3, velsFloat3, angVelsFloat3);
790 m_sys->user_sphere_fixed = fixity;
791
792 // Last thing: if this function runs successfully, then we automatically disable the defragment process on
793 // simulation initialization. A re-started simulation would better not have its particle order re-arranged. But
794 // if the user does not like it, they can still force the defragment by calling a method.
795 if (numRow > 0)
796 m_sys->defragment_on_start = false;
797 }
798
799 // Read in particle friction contact history through a hst-formatted ifstream
ReadHstHistory(std::ifstream & ifile,unsigned int totItem)800 void ChSystemGpu::ReadHstHistory(std::ifstream& ifile, unsigned int totItem) {
801 // Find the header line
802 std::string line;
803 std::getline(ifile, line);
804 while (line.find_first_not_of(' ') == std::string::npos)
805 std::getline(ifile, line);
806
807 // The header will tell the max number of partners (use as offset),
808 // but MAX_SPHERES_TOUCHED_BY_SPHERE is the monodispersity standard
809 unsigned int max_partner = MAX_SPHERES_TOUCHED_BY_SPHERE;
810 unsigned int formatMode = quarryHistoryFormat(line, max_partner);
811
812 // Parse in every data line
813 float3 history_UU;
814 unsigned int numItem = 0;
815 std::vector<unsigned int> partner;
816 std::vector<float3> history;
817 // don't always assume we know the num of items (spheres), and we expand buffer every time batchSize is hit
818 const unsigned int batchSize = 1e6;
819 while (std::getline(ifile, line)) {
820 // if empty line, just keep going
821 if (line.find_first_not_of(' ') == std::string::npos)
822 continue;
823
824 std::istringstream iss1(line);
825 if (numItem % batchSize == 0) {
826 partner.resize(partner.size() + max_partner * batchSize);
827 history.resize(history.size() + max_partner * batchSize, make_float3(0, 0, 0));
828 }
829
830 // Read partner map
831 if (formatMode & 1) {
832 for (unsigned int i = 0; i < max_partner; i++)
833 iss1 >> partner[max_partner * numItem + i];
834 }
835
836 // Read contact history
837 if (formatMode & 2) {
838 for (unsigned int i = 0; i < max_partner; i++) {
839 iss1 >> history_UU.x;
840 iss1 >> history_UU.y;
841 iss1 >> history_UU.z;
842 history[max_partner * numItem + i] = history_UU;
843 }
844 }
845
846 numItem++;
847 // If all spheres loaded, break. Don't keep going as it's not necessarily the EoF.
848 if (numItem >= totItem)
849 break;
850 }
851
852 // Final resize
853 partner.resize(max_partner * numItem);
854 history.resize(max_partner * numItem);
855
856 // Feed the contact/history info to the system. We are directly referencing low-level m_sys structs here.
857 // Because there are no user panel funtions (that loads contact/history) which can be used in a simulation
858 // script. Those utils can be added, but I doubt their usefulness.
859 m_sys->user_partner_map = partner;
860 m_sys->user_friction_history = history;
861
862 // Last thing: if this function runs successfully, then we automatically disable the defragment process on
863 // simulation initialization. A re-started simulation would better not have its particle order re-arranged. But
864 // if the user does not like it, they can still force the defragment by calling a method.
865 if (numItem > 0)
866 m_sys->defragment_on_start = false;
867 }
868
869 // Read in particle friction contact history through a file.
870 // Now it works with a custom format only. Potentially more formats can be added in the future.
ReadContactHistoryFile(const std::string & infilename)871 void ChSystemGpu::ReadContactHistoryFile(const std::string& infilename) {
872 std::ifstream ifile(infilename.c_str());
873 if (!ifile.is_open()) {
874 CHGPU_ERROR("ERROR! Checkpoint file (contact history) did not open successfully!\n");
875 }
876
877 ReadHstHistory(ifile, UINT_MAX);
878 }
879
880 // Read in particle positions, velocity and angular velocity through a file.
881 // Now it works with csv format only. Potentially more formats can be added in the future.
ReadParticleFile(const std::string & infilename)882 void ChSystemGpu::ReadParticleFile(const std::string& infilename) {
883 std::ifstream ifile(infilename.c_str());
884 if (!ifile.is_open()) {
885 CHGPU_ERROR("ERROR! Checkpoint file (kinematics) did not open successfully!\n");
886 }
887
888 ReadCsvParticles(ifile, UINT_MAX);
889 }
890
891 // A smaller hasher that helps determine the indentifier type.
892 // It is powerful enough for our purpose and good-looking. We did not use built-in functions (such as string_view?)
893 // because that could require C++17, also I feel it would make the code look longer. In any case, if this small
894 // hasher is not sufficient anymore in future updates, we can spot that during compilation.
hash_charr(const char * s,int off=0)895 constexpr unsigned int hash_charr(const char* s, int off = 0) {
896 return !s[off] ? 7001 : (hash_charr(s, off + 1) * 33) ^ s[off];
897 }
operator ""_(const char * s,size_t)898 constexpr inline unsigned int operator"" _(const char* s, size_t) {
899 return hash_charr(s);
900 }
diff(float a,float b)901 bool diff(float a, float b) {
902 return std::abs(a - b) > 1e-6f;
903 }
diff(float3 a,float3 b)904 bool diff(float3 a, float3 b) {
905 return std::abs(a.x - b.x) > 1e-6f || std::abs(a.y - b.y) > 1e-6f || std::abs(a.z - b.z) > 1e-6f;
906 }
907
908 // Use hash to find matching indentifier and load parameters. Return 1 if found no matching paramter to set, return
909 // 0 if status normal
SetParamsFromIdentifier(const std::string & identifier,std::istringstream & iss1,bool overwrite)910 bool ChSystemGpu::SetParamsFromIdentifier(const std::string& identifier, std::istringstream& iss1, bool overwrite) {
911 unsigned int i; // integer holder
912 float f; // float holder
913 float3 f3; // float3 holder
914 double d; // double holder
915 bool b; // bool holder
916 bool anomaly = false; // flag unknown identifier
917 bool incst = false; // flag parameter changes compared to current system
918
919 switch (hash_charr(identifier.c_str())) {
920 case ("density"_):
921 iss1 >> f;
922 incst = diff(m_sys->sphere_density_UU, f);
923 m_sys->sphere_density_UU = f;
924 break;
925 case ("radius"_):
926 iss1 >> f;
927 incst = diff(m_sys->sphere_radius_UU, f);
928 m_sys->sphere_radius_UU = f;
929 break;
930 case ("boxSize"_):
931 iss1 >> f3.x;
932 iss1 >> f3.y;
933 iss1 >> f3.z;
934 incst = diff(make_float3(m_sys->box_size_X, m_sys->box_size_Y, m_sys->box_size_Z), f3);
935 m_sys->box_size_X = f3.x;
936 m_sys->box_size_Y = f3.y;
937 m_sys->box_size_Z = f3.z;
938 break;
939 case ("BDFixed"_):
940 iss1 >> b;
941 SetBDFixed(b);
942 break;
943 case ("BDCenter"_):
944 iss1 >> f3.x;
945 iss1 >> f3.y;
946 iss1 >> f3.z;
947 incst = diff(make_float3(m_sys->user_coord_O_X, m_sys->user_coord_O_Y, m_sys->user_coord_O_Z), f3);
948 m_sys->user_coord_O_X = f3.x;
949 m_sys->user_coord_O_Y = f3.y;
950 m_sys->user_coord_O_Z = f3.z;
951 break;
952 case ("verbosity"_):
953 iss1 >> i;
954 SetVerbosity(static_cast<CHGPU_VERBOSITY>(i));
955 break;
956 case ("useMinLengthUnit"_):
957 iss1 >> b;
958 EnableMinLength(b);
959 break;
960 case ("recordContactInfo"_):
961 iss1 >> b;
962 SetRecordingContactInfo(b);
963 break;
964 case ("particleFileMode"_):
965 iss1 >> i;
966 SetParticleOutputMode(static_cast<CHGPU_OUTPUT_MODE>(i));
967 break;
968 case ("particleFileFlags"_):
969 iss1 >> i;
970 // ParticleOutputFlags is already int, instead of a enum class
971 SetParticleOutputFlags(i);
972 break;
973 case ("fixedStepSize"_):
974 iss1 >> f;
975 SetFixedStepSize(f);
976 break;
977 case ("cohesionOverG"_):
978 iss1 >> f;
979 SetCohesionRatio(f);
980 break;
981 case ("adhesionOverG_s2w"_):
982 iss1 >> f;
983 SetAdhesionRatio_SPH2WALL(f);
984 break;
985 case ("G"_):
986 iss1 >> f3.x;
987 iss1 >> f3.y;
988 iss1 >> f3.z;
989 SetGravitationalAcceleration(f3);
990 break;
991 case ("elapsedTime"_):
992 iss1 >> f;
993 SetSimTime(f);
994 break;
995 case ("K_n_s2s"_):
996 iss1 >> d;
997 SetKn_SPH2SPH(d);
998 break;
999 case ("K_n_s2w"_):
1000 iss1 >> d;
1001 SetKn_SPH2WALL(d);
1002 break;
1003 case ("K_t_s2s"_):
1004 iss1 >> d;
1005 SetKt_SPH2SPH(d);
1006 break;
1007 case ("K_t_s2w"_):
1008 iss1 >> d;
1009 SetKt_SPH2WALL(d);
1010 break;
1011 case ("G_n_s2s"_):
1012 iss1 >> d;
1013 SetGn_SPH2SPH(d);
1014 break;
1015 case ("G_n_s2w"_):
1016 iss1 >> d;
1017 SetGn_SPH2WALL(d);
1018 break;
1019 case ("G_t_s2s"_):
1020 iss1 >> d;
1021 SetGt_SPH2SPH(d);
1022 break;
1023 case ("G_t_s2w"_):
1024 iss1 >> d;
1025 SetGt_SPH2WALL(d);
1026 break;
1027 case ("RollingCoeff_s2s"_):
1028 iss1 >> d;
1029 SetRollingCoeff_SPH2SPH(d);
1030 break;
1031 case ("RollingCoeff_s2w"_):
1032 iss1 >> d;
1033 SetRollingCoeff_SPH2WALL(d);
1034 break;
1035 case ("SpinningCoeff_s2s"_):
1036 iss1 >> d;
1037 SetSpinningCoeff_SPH2SPH(d);
1038 break;
1039 case ("SpinningCoeff_s2w"_):
1040 iss1 >> d;
1041 SetSpinningCoeff_SPH2WALL(d);
1042 break;
1043 case ("StaticFrictionCoeff_s2s"_):
1044 iss1 >> d;
1045 SetStaticFrictionCoeff_SPH2SPH(d);
1046 break;
1047 case ("StaticFrictionCoeff_s2w"_):
1048 iss1 >> d;
1049 SetStaticFrictionCoeff_SPH2WALL(d);
1050 break;
1051 case ("PsiT"_):
1052 iss1 >> i;
1053 SetPsiT(i);
1054 break;
1055 case ("PsiL"_):
1056 iss1 >> i;
1057 SetPsiL(i);
1058 break;
1059 case ("PsiR"_):
1060 iss1 >> f;
1061 SetPsiR(f);
1062 break;
1063 case ("frictionMode"_):
1064 iss1 >> i;
1065 SetFrictionMode(static_cast<CHGPU_FRICTION_MODE>(i));
1066 break;
1067 case ("rollingMode"_):
1068 iss1 >> i;
1069 SetRollingMode(static_cast<CHGPU_ROLLING_MODE>(i));
1070 break;
1071 case ("timeIntegrator"_):
1072 iss1 >> i;
1073 SetTimeIntegrator(static_cast<CHGPU_TIME_INTEGRATOR>(i));
1074 break;
1075 case ("maxSafeVelSU"_):
1076 iss1 >> f;
1077 SetMaxSafeVelocity_SU(f);
1078 break;
1079 default:
1080 anomaly = true;
1081 }
1082
1083 if (incst && !overwrite)
1084 CHGPU_ERROR(
1085 "ERROR! Parameter \"%s\" is inconsistent with the current simulation system.\n"
1086 "If you wish to construct a simulation systen from scratch using this checkpoint file, then you "
1087 "should supply this file as the constructor parameter.\nExiting...\n",
1088 identifier.c_str());
1089
1090 return anomaly;
1091 }
1092
1093 // Read in simulation parameters. Returns the total number of particles. If instructed to overwrite, then overwrite
1094 // cuurent simulation parameters with the values in the checkpoint file; else, when an inconsistency is found, throw
1095 // an error.
ReadDatParams(std::ifstream & ifile,bool overwrite)1096 unsigned int ChSystemGpu::ReadDatParams(std::ifstream& ifile, bool overwrite) {
1097 std::string line;
1098 unsigned int nSpheres = 0;
1099 std::getline(ifile, line);
1100 // Get rid of possible empty lines, before the real party starts
1101 while (line.find_first_not_of(' ') == std::string::npos)
1102 std::getline(ifile, line);
1103
1104 std::string identifier;
1105 while (line != "ParamsEnd") {
1106 std::istringstream iss1(line);
1107 // Grab the param identifier
1108 getline(iss1, identifier, ':');
1109 // Call corresponding "Set" methods based on identifier
1110 if (identifier == "nSpheres") {
1111 iss1 >> nSpheres;
1112 if (!overwrite && (nSpheres - m_sys->nSpheres != 0))
1113 CHGPU_ERROR(
1114 "ERROR! Number of particles in checkpoint file is inconsistent with the current system.\n"
1115 "If you wish to construct a simulation systen from scratch using this checkpoint file, then "
1116 "you "
1117 "should supply this file as the constructor parameter.\nExiting...\n");
1118 } else {
1119 bool anomaly = SetParamsFromIdentifier(identifier, iss1, overwrite);
1120 if (anomaly)
1121 printf("WARNING! %s is an unknown parameter, skipped.\n", identifier.c_str());
1122 }
1123
1124 // Load next line
1125 std::getline(ifile, line);
1126 }
1127
1128 if (nSpheres == 0)
1129 printf(
1130 "WARNING! The checkpoint file indicates there are 0 particles to be loaded. If this is intended, the "
1131 "user "
1132 "must initialize the particles by themselves.\n");
1133 return nSpheres;
1134 }
1135
1136 // Read in a (Chrono::Gpu generated) checkpoint file to restart a simulation. It calls ReadHstHistory,
1137 // ReadCsvParticles and ReadDatParams to parse in the entire checkpoint file.
ReadCheckpointFile(const std::string & infilename,bool overwrite)1138 void ChSystemGpu::ReadCheckpointFile(const std::string& infilename, bool overwrite) {
1139 // Open the file
1140 std::ifstream ifile(infilename.c_str());
1141 if (!ifile.is_open()) {
1142 CHGPU_ERROR("ERROR! Checkpoint file did not open successfully!\n");
1143 }
1144
1145 // Let the user know we are loading...
1146 printf("Reading checkpoint data from file \"%s\"...\n", infilename.c_str());
1147
1148 // Process the header
1149 std::string line;
1150 std::getline(ifile, line);
1151 while (line.find_first_not_of(' ') == std::string::npos)
1152 std::getline(ifile, line);
1153 if (line != "ChSystemGpu")
1154 printf(
1155 "WARNING! Header of checkpoint file indicates this is not for ChSystemGpu (and it seems to be for %s)! "
1156 "There may still be parameters defaulted after loading this checkpoint file, and you may want to set "
1157 "them "
1158 "manually.\n",
1159 line.c_str());
1160
1161 // Load and set all simulation parameters. Keep tab of nSpheres, we'll need that soon.
1162 unsigned int nSpheres;
1163 nSpheres = ReadDatParams(ifile, overwrite);
1164
1165 // Now load the rest, particle kinematics info, contact pair/friction history, everything
1166 while (std::getline(ifile, line)) {
1167 if (line.find_first_not_of(' ') == std::string::npos)
1168 continue;
1169
1170 // If not an empty line, then check the header to call a subroutine
1171 if (line == "CsvParticles")
1172 ReadCsvParticles(ifile, nSpheres);
1173 else if (line == "HstHistory") {
1174 if (m_sys->gran_params->friction_mode != CHGPU_FRICTION_MODE::FRICTIONLESS) {
1175 ReadHstHistory(ifile, nSpheres);
1176 }
1177 }
1178 }
1179 }
1180
1181 // -----------------------------------------------------------------------------
1182
1183 // GpuMesh veriosn of reading checkpointed params using hashed identifier.
SetParamsFromIdentifier(const std::string & identifier,std::istringstream & iss1,bool overwrite)1184 bool ChSystemGpuMesh::SetParamsFromIdentifier(const std::string& identifier, std::istringstream& iss1, bool overwrite) {
1185 float f; // float holder
1186 double d; // double holder
1187 bool b; // bool holder
1188 bool anomaly = false; // flag unknown identifier
1189
1190 // Is the parameter one of the ChSystemGpu's? If not, we then search it in ChSystemGpuMesh.
1191 // The parent method returns true if it did not find a match.
1192 if (ChSystemGpu::SetParamsFromIdentifier(identifier, iss1, overwrite)) {
1193 switch (hash_charr(identifier.c_str())) {
1194 case ("K_n_s2m"_):
1195 iss1 >> d;
1196 SetKn_SPH2MESH(d);
1197 break;
1198 case ("K_t_s2m"_):
1199 iss1 >> d;
1200 SetKt_SPH2MESH(d);
1201 break;
1202 case ("G_n_s2m"_):
1203 iss1 >> d;
1204 SetGn_SPH2MESH(d);
1205 break;
1206 case ("G_t_s2m"_):
1207 iss1 >> d;
1208 SetGt_SPH2MESH(d);
1209 break;
1210 case ("RollingCoeff_s2m"_):
1211 iss1 >> f;
1212 SetRollingCoeff_SPH2MESH(f);
1213 break;
1214 case ("SpinningCoeff_s2m"_):
1215 iss1 >> f;
1216 SetSpinningCoeff_SPH2MESH(f);
1217 break;
1218 case ("StaticFrictionCoeff_s2m"_):
1219 iss1 >> f;
1220 SetStaticFrictionCoeff_SPH2MESH(f);
1221 break;
1222 case ("MeshCollisionEnabled"_):
1223 iss1 >> b;
1224 EnableMeshCollision(b);
1225 break;
1226 case ("adhesionOverG_s2m"_):
1227 iss1 >> f;
1228 SetAdhesionRatio_SPH2MESH(f);
1229 break;
1230 default:
1231 anomaly = true;
1232 // printf("WARNING! %s is an unknown parameter, skipped.\n", identifier.c_str());
1233 }
1234 }
1235 return anomaly;
1236 }
1237
1238 // GpuMesh version of checkpointing loading subroutine.
ReadCheckpointFile(const std::string & infilename,bool overwrite)1239 void ChSystemGpuMesh::ReadCheckpointFile(const std::string& infilename, bool overwrite) {
1240 // Open the file
1241 std::ifstream ifile(infilename.c_str());
1242 if (!ifile.is_open()) {
1243 CHGPU_ERROR("ERROR! Checkpoint file did not open successfully!\n");
1244 }
1245
1246 // Let the user know we are loading...
1247 printf("Reading checkpoint data from file \"%s\"...\n", infilename.c_str());
1248
1249 // Process the header
1250 std::string line;
1251 std::getline(ifile, line);
1252 while (line.find_first_not_of(' ') == std::string::npos)
1253 std::getline(ifile, line);
1254 if (line != "ChSystemGpuMesh")
1255 printf(
1256 "WARNING! Header of checkpoint file indicates this is not for ChSystemGpuMesh (and it seems to be for "
1257 "%s)! "
1258 "There may still be parameters defaulted after loading this checkpoint file, and you may want to set "
1259 "them "
1260 "manually.\n",
1261 line.c_str());
1262
1263 // Load and set all simulation parameters. Keep tab of nSpheres, we'll need that soon.
1264 unsigned int nSpheres;
1265 nSpheres = ReadDatParams(ifile, overwrite);
1266
1267 // Now load the rest, particle kinematics info, contact pair/friction history, everything
1268 while (std::getline(ifile, line)) {
1269 if (line.find_first_not_of(' ') == std::string::npos)
1270 continue;
1271
1272 // If not an empty line, then check the header to call a subroutine
1273 if (line == "CsvParticles")
1274 ReadCsvParticles(ifile, nSpheres);
1275 else if (line == "HstHistory") {
1276 if (m_sys->gran_params->friction_mode != CHGPU_FRICTION_MODE::FRICTIONLESS) {
1277 ReadHstHistory(ifile, nSpheres);
1278 }
1279 }
1280 // else, TODO it may need to read mesh-related info. This may be coming in the future.
1281 }
1282 }
1283
Initialize()1284 void ChSystemGpuMesh::Initialize() {
1285 if (m_meshes.size() > 0)
1286 SetMeshes();
1287 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1288 sys_trimesh->initializeSpheres();
1289 sys_trimesh->initializeTriangles();
1290 }
1291
InitializeMeshes()1292 void ChSystemGpuMesh::InitializeMeshes() {
1293 if (m_meshes.size() > 0)
1294 SetMeshes();
1295 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1296 sys_trimesh->initializeTriangles();
1297 }
1298
Initialize()1299 void ChSystemGpu::Initialize() {
1300 m_sys->initializeSpheres();
1301 if (m_sys->verbosity == CHGPU_VERBOSITY::INFO || m_sys->verbosity == CHGPU_VERBOSITY::METRICS) {
1302 printf("Approx mem usage is %s\n", pretty_format_bytes(EstimateMemUsage()).c_str());
1303 }
1304 }
1305
1306 // -----------------------------------------------------------------------------
1307
AdvanceSimulation(float duration)1308 double ChSystemGpuMesh::AdvanceSimulation(float duration) {
1309 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1310 return sys_trimesh->AdvanceSimulation(duration);
1311 }
1312
AdvanceSimulation(float duration)1313 double ChSystemGpu::AdvanceSimulation(float duration) {
1314 return m_sys->AdvanceSimulation(duration);
1315 }
1316
1317 // -----------------------------------------------------------------------------
1318
1319 template <typename Enumeration>
as_uint(Enumeration const value)1320 auto as_uint(Enumeration const value) -> typename std::underlying_type<Enumeration>::type {
1321 return static_cast<typename std::underlying_type<Enumeration>::type>(value);
1322 }
1323
WriteCheckpointParams(std::ofstream & cpFile) const1324 void ChSystemGpu::WriteCheckpointParams(std::ofstream& cpFile) const {
1325 std::ostringstream paramStream;
1326
1327 paramStream << "nSpheres: " << GetNumParticles() << "\n";
1328 paramStream << "density: " << m_sys->sphere_density_UU << "\n";
1329 paramStream << "radius: " << m_sys->sphere_radius_UU << "\n";
1330 paramStream << "boxSize: " << m_sys->box_size_X << " " << m_sys->box_size_Y << " " << m_sys->box_size_Z << "\n";
1331 paramStream << "BDFixed: " << (int)(m_sys->BD_is_fixed) << "\n";
1332 paramStream << "BDCenter: " << m_sys->user_coord_O_X << " " << m_sys->user_coord_O_Y << " " << m_sys->user_coord_O_Z
1333 << "\n";
1334 paramStream << "verbosity: " << as_uint(m_sys->verbosity) << "\n";
1335 paramStream << "useMinLengthUnit: " << (int)(m_sys->use_min_length_unit) << "\n";
1336 paramStream << "recordContactInfo: " << (int)(m_sys->gran_params->recording_contactInfo) << "\n";
1337 paramStream << "particleFileMode: " << as_uint(m_sys->file_write_mode) << "\n";
1338 // returned OutputFlags is already an int
1339 paramStream << "particleFileFlags: " << m_sys->output_flags << "\n";
1340 paramStream << "fixedStepSize: " << m_sys->stepSize_UU << "\n";
1341 // It is cohesion-over-gravity and adhesion-over-gravity
1342 paramStream << "cohesionOverG: " << m_sys->cohesion_over_gravity << "\n";
1343 paramStream << "adhesionOverG_s2w: " << m_sys->adhesion_s2w_over_gravity << "\n";
1344 paramStream << "G: " << m_sys->X_accGrav << " " << m_sys->Y_accGrav << " " << m_sys->Z_accGrav << "\n";
1345 paramStream << "elapsedTime: " << GetSimTime() << "\n";
1346
1347 paramStream << "K_n_s2s: " << m_sys->K_n_s2s_UU << "\n";
1348 paramStream << "K_n_s2w: " << m_sys->K_n_s2w_UU << "\n";
1349 paramStream << "K_t_s2s: " << m_sys->K_t_s2s_UU << "\n";
1350 paramStream << "K_t_s2w: " << m_sys->K_t_s2w_UU << "\n";
1351 paramStream << "G_n_s2s: " << m_sys->Gamma_n_s2s_UU << "\n";
1352 paramStream << "G_n_s2w: " << m_sys->Gamma_n_s2w_UU << "\n";
1353 paramStream << "G_t_s2s: " << m_sys->Gamma_t_s2s_UU << "\n";
1354 paramStream << "G_t_s2w: " << m_sys->Gamma_t_s2w_UU << "\n";
1355 paramStream << "RollingCoeff_s2s: " << m_sys->rolling_coeff_s2s_UU << "\n";
1356 paramStream << "RollingCoeff_s2w: " << m_sys->rolling_coeff_s2w_UU << "\n";
1357 paramStream << "SpinningCoeff_s2s: " << m_sys->spinning_coeff_s2s_UU << "\n";
1358 paramStream << "SpinningCoeff_s2w: " << m_sys->spinning_coeff_s2w_UU << "\n";
1359 paramStream << "StaticFrictionCoeff_s2s: " << m_sys->gran_params->static_friction_coeff_s2s << "\n";
1360 paramStream << "StaticFrictionCoeff_s2w: " << m_sys->gran_params->static_friction_coeff_s2w << "\n";
1361
1362 paramStream << "PsiT: " << m_sys->psi_T << "\n";
1363 paramStream << "PsiL: " << m_sys->psi_L << "\n";
1364 paramStream << "PsiR: " << m_sys->psi_R << "\n";
1365 paramStream << "frictionMode: " << as_uint(m_sys->gran_params->friction_mode) << "\n";
1366 paramStream << "rollingMode: " << as_uint(m_sys->gran_params->rolling_mode) << "\n";
1367
1368 // Notify the user that the support for CHUNG is limited. It still runs but loses some information in the
1369 // checkpointing process.
1370 paramStream << "timeIntegrator: " << as_uint(m_sys->time_integrator) << "\n";
1371 if (m_sys->time_integrator == CHGPU_TIME_INTEGRATOR::CHUNG)
1372 printf(
1373 "WARNING! CHUNG integrator is partially supported in this version. The acceleration and angular "
1374 "acceleration info is not stored in the checkpoint file, leading to a potential change in physics. You "
1375 "can "
1376 "consider using other integrators if checkpointing is needed, or wait for a fix in the next "
1377 "version.\n");
1378
1379 paramStream << "maxSafeVelSU: " << m_sys->gran_params->max_safe_vel << "\n";
1380
1381 // In the near future, TODO cache all extra boundaries here as well
1382
1383 cpFile << paramStream.str();
1384 }
1385
WriteCheckpointFile(const std::string & outfilename)1386 void ChSystemGpu::WriteCheckpointFile(const std::string& outfilename) {
1387 printf("Writing checkpoint data to file \"%s\"\n", outfilename.c_str());
1388 std::ofstream cpFile(outfilename, std::ios::out);
1389
1390 // Header, indicating the system (so meshed system will have something different)
1391 cpFile << std::string("ChSystemGpu\n");
1392
1393 // Simulation params go right after the system name
1394 WriteCheckpointParams(cpFile);
1395 cpFile << std::string("ParamsEnd\n");
1396 cpFile << std::string("\n");
1397
1398 // Then, the particle kinematics info
1399 cpFile << std::string("CsvParticles\n");
1400 // In checkpointing, we want all particle info written, instead of selected output that the user can enforce via
1401 // setting OutputFlags. Therefore, we temporarily set OutputFlags to maximum, then revert after this writting is
1402 // done.
1403 unsigned int outFlags = m_sys->output_flags;
1404 SetParticleOutputFlags(VEL_COMPONENTS | FIXITY | ANG_VEL_COMPONENTS);
1405 WriteCsvParticles(cpFile);
1406 SetParticleOutputFlags(outFlags);
1407 cpFile << std::string("\n");
1408
1409 // Then, write contact pair/friction history
1410 if (m_sys->gran_params->friction_mode != CHGPU_FRICTION_MODE::FRICTIONLESS) {
1411 cpFile << std::string("HstHistory\n");
1412 WriteHstHistory(cpFile);
1413 cpFile << std::string("\n");
1414 }
1415 }
1416
WriteRawParticles(std::ofstream & ptFile) const1417 void ChSystemGpu::WriteRawParticles(std::ofstream& ptFile) const {
1418 m_sys->WriteRawParticles(ptFile);
1419 }
1420
WriteCsvParticles(std::ofstream & ptFile) const1421 void ChSystemGpu::WriteCsvParticles(std::ofstream& ptFile) const {
1422 m_sys->WriteCsvParticles(ptFile);
1423 }
1424
WriteChPFParticles(std::ofstream & ptFile) const1425 void ChSystemGpu::WriteChPFParticles(std::ofstream& ptFile) const {
1426 m_sys->WriteChPFParticles(ptFile);
1427 }
1428
1429 #ifdef USE_HDF5
WriteH5Particles(H5::H5File ptFile) const1430 void ChSystemGpu::WriteH5Particles(H5::H5File ptFile) const {
1431 m_sys->WriteH5Particles(ptFile);
1432 }
1433 #endif
1434
WriteParticleFile(const std::string & outfilename) const1435 void ChSystemGpu::WriteParticleFile(const std::string& outfilename) const {
1436 // The file writes are a pretty big slowdown in CSV mode
1437 if (m_sys->file_write_mode == CHGPU_OUTPUT_MODE::BINARY) {
1438 std::ofstream ptFile(outfilename, std::ios::out | std::ios::binary);
1439 WriteRawParticles(ptFile);
1440 } else if (m_sys->file_write_mode == CHGPU_OUTPUT_MODE::CSV) {
1441 std::ofstream ptFile(outfilename, std::ios::out);
1442 WriteCsvParticles(ptFile);
1443 } else if (m_sys->file_write_mode == CHGPU_OUTPUT_MODE::CHPF) {
1444 std::ofstream ptFile(outfilename, std::ios::out | std::ios::binary);
1445 WriteChPFParticles(ptFile);
1446 } else if (m_sys->file_write_mode == CHGPU_OUTPUT_MODE::HDF5) {
1447 #ifdef USE_HDF5
1448 H5::H5File ptFile(outfilename.c_str(), H5F_ACC_TRUNC);
1449 WriteH5Particles(ptFile);
1450 #else
1451 CHGPU_ERROR("ERROR! HDF5 Installation not found. Recompile with HDF5.\n");
1452 #endif
1453 }
1454 }
1455
WriteContactInfoFile(const std::string & outfilename) const1456 void ChSystemGpu::WriteContactInfoFile(const std::string& outfilename) const {
1457 m_sys->WriteContactInfoFile(outfilename);
1458 }
1459
WriteHstHistory(std::ofstream & histFile) const1460 void ChSystemGpu::WriteHstHistory(std::ofstream& histFile) const {
1461 // Dump to a stream, write to file only at end
1462 std::ostringstream outstrstream;
1463
1464 // Figure out the write format first
1465 // 1: write contact_partners_map
1466 // 2: write contact_history_map
1467 unsigned int formatMode = 0;
1468 switch (m_sys->gran_params->friction_mode) {
1469 case (CHGPU_FRICTION_MODE::FRICTIONLESS):
1470 printf("WARNING! Currently using FRICTIONLESS model. There is no contact history to write!\n");
1471 return;
1472 case (CHGPU_FRICTION_MODE::SINGLE_STEP):
1473 formatMode += 1;
1474 break;
1475 case (CHGPU_FRICTION_MODE::MULTI_STEP):
1476 formatMode += (1 + 2);
1477 break;
1478 default:
1479 CHGPU_ERROR("ERROR! Unknown friction model, failed to write contact history file!\n");
1480 }
1481
1482 if (formatMode & 1)
1483 outstrstream << "partners " << MAX_SPHERES_TOUCHED_BY_SPHERE;
1484 if (formatMode & 2)
1485 outstrstream << " history " << MAX_SPHERES_TOUCHED_BY_SPHERE;
1486 outstrstream << "\n";
1487
1488 // We'll use space-separated formatting, as I found it more convenient to parse in and looks better.
1489 // Forget about CSV conventions, history info is not meant to be used by third-party tools anyway.
1490 float3 history_UU;
1491 for (unsigned int n = 0; n < m_sys->nSpheres; n++) {
1492 // Write contact_partners_map
1493 if (formatMode & 1) {
1494 for (unsigned int i = 0; i < MAX_SPHERES_TOUCHED_BY_SPHERE; i++)
1495 outstrstream << m_sys->contact_partners_map[MAX_SPHERES_TOUCHED_BY_SPHERE * n + i] << " ";
1496 }
1497 // Write write contact_history_map
1498 if (formatMode & 2) {
1499 for (unsigned int i = 0; i < MAX_SPHERES_TOUCHED_BY_SPHERE; i++) {
1500 history_UU.x =
1501 m_sys->contact_history_map[MAX_SPHERES_TOUCHED_BY_SPHERE * n + i].x * m_sys->LENGTH_SU2UU;
1502 history_UU.y =
1503 m_sys->contact_history_map[MAX_SPHERES_TOUCHED_BY_SPHERE * n + i].y * m_sys->LENGTH_SU2UU;
1504 history_UU.z =
1505 m_sys->contact_history_map[MAX_SPHERES_TOUCHED_BY_SPHERE * n + i].z * m_sys->LENGTH_SU2UU;
1506 outstrstream << history_UU.x << " " << history_UU.y << " " << history_UU.z << " ";
1507 }
1508 }
1509 outstrstream << "\n";
1510 }
1511
1512 histFile << outstrstream.str();
1513 }
1514
WriteContactHistoryFile(const std::string & outfilename) const1515 void ChSystemGpu::WriteContactHistoryFile(const std::string& outfilename) const {
1516 printf("Writing contact pair/history data to file \"%s\"\n", outfilename.c_str());
1517 std::ofstream histFile(outfilename, std::ios::out);
1518
1519 // Call a subroutine to write. This subroutine can be used in checkpointing as well.
1520 WriteHstHistory(histFile);
1521 }
1522
1523 // -----------------------------------------------------------------------------
1524
1525 // GpuMesh version of checkpoint params writing. Including extra params that only a meshed system has.
WriteCheckpointMeshParams(std::ofstream & cpFile) const1526 void ChSystemGpuMesh::WriteCheckpointMeshParams(std::ofstream& cpFile) const {
1527 std::ostringstream paramStream;
1528 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1529
1530 paramStream << "adhesionOverG_s2m: " << sys_trimesh->adhesion_s2m_over_gravity << "\n";
1531 paramStream << "K_n_s2m: " << sys_trimesh->K_n_s2m_UU << "\n";
1532 paramStream << "K_t_s2m: " << sys_trimesh->K_t_s2m_UU << "\n";
1533 paramStream << "G_n_s2m: " << sys_trimesh->Gamma_n_s2m_UU << "\n";
1534 paramStream << "G_t_s2m: " << sys_trimesh->Gamma_t_s2m_UU << "\n";
1535 paramStream << "RollingCoeff_s2m: " << sys_trimesh->rolling_coeff_s2m_UU << "\n";
1536 paramStream << "SpinningCoeff_s2m: " << sys_trimesh->spinning_coeff_s2m_UU << "\n";
1537 paramStream << "StaticFrictionCoeff_s2m: " << sys_trimesh->tri_params->static_friction_coeff_s2m << "\n";
1538 paramStream << "MeshCollisionEnabled: " << sys_trimesh->mesh_collision_enabled << "\n";
1539
1540 cpFile << paramStream.str();
1541 }
1542 // GpuMesh version of checkpoint writing
WriteCheckpointFile(const std::string & outfilename)1543 void ChSystemGpuMesh::WriteCheckpointFile(const std::string& outfilename) {
1544 printf("Writing checkpoint data to file \"%s\"\n", outfilename.c_str());
1545 std::ofstream cpFile(outfilename, std::ios::out);
1546
1547 // Header, indicating the system
1548 cpFile << std::string("ChSystemGpuMesh\n");
1549
1550 // Simulation params go right after the system name
1551 WriteCheckpointParams(cpFile);
1552 WriteCheckpointMeshParams(cpFile); // Meshed system has extra params
1553 cpFile << std::string("ParamsEnd\n");
1554 cpFile << std::string("\n");
1555
1556 // Then, the particle kinematics info
1557 cpFile << std::string("CsvParticles\n");
1558 // In checkpointing, we want all particle info written, instead of selected output that the user can enforce via
1559 // setting OutputFlags. Therefore, we temporarily set OutputFlags to maximum, then revert after this writting is
1560 // done.
1561 unsigned int outFlags = m_sys->output_flags;
1562 SetParticleOutputFlags(VEL_COMPONENTS | FIXITY | ANG_VEL_COMPONENTS);
1563 WriteCsvParticles(cpFile);
1564 SetParticleOutputFlags(outFlags);
1565 cpFile << std::string("\n");
1566
1567 // Then, write contact pair/friction history
1568 if (m_sys->gran_params->friction_mode != CHGPU_FRICTION_MODE::FRICTIONLESS) {
1569 cpFile << std::string("HstHistory\n");
1570 WriteHstHistory(cpFile);
1571 cpFile << std::string("\n");
1572 }
1573
1574 // In the future, TODO checkpointing mesh and mesh translation/rotation goes here.
1575 }
1576
WriteMesh(const std::string & outfilename,unsigned int i) const1577 void ChSystemGpuMesh::WriteMesh(const std::string& outfilename, unsigned int i) const {
1578 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1579 if (sys_trimesh->file_write_mode == CHGPU_OUTPUT_MODE::NONE) {
1580 return;
1581 }
1582 if (i >= m_meshes.size()) {
1583 printf("WARNING: attempted to write mesh %u, yet only %zu meshes present. No mesh file generated.\n", i,
1584 m_meshes.size());
1585 return;
1586 }
1587
1588 std::string ofile;
1589 if (outfilename.substr(outfilename.length() - std::min(outfilename.length(), (size_t)4)) != ".vtk" &&
1590 outfilename.substr(outfilename.length() - std::min(outfilename.length(), (size_t)4)) != ".VTK")
1591 ofile = outfilename + ".vtk";
1592 else
1593 ofile = outfilename;
1594 std::ofstream outfile(ofile, std::ios::out);
1595 std::ostringstream ostream;
1596 ostream << "# vtk DataFile Version 2.0\n";
1597 ostream << "VTK from simulation\n";
1598 ostream << "ASCII\n";
1599 ostream << "\n\n";
1600
1601 ostream << "DATASET UNSTRUCTURED_GRID\n";
1602
1603 const auto& mmesh = m_meshes.at(i);
1604
1605 // Writing vertices
1606 ostream << "POINTS " << mmesh->getCoordsVertices().size() << " float" << std::endl;
1607 for (auto& v : mmesh->getCoordsVertices()) {
1608 float3 point = make_float3(v.x(), v.y(), v.z());
1609 sys_trimesh->ApplyFrameTransform(point, sys_trimesh->tri_params->fam_frame_broad[i].pos,
1610 sys_trimesh->tri_params->fam_frame_broad[i].rot_mat);
1611 ostream << point.x << " " << point.y << " " << point.z << std::endl;
1612 }
1613
1614 // Writing faces
1615 ostream << "\n\n";
1616 ostream << "CELLS " << mmesh->getIndicesVertexes().size() << " " << 4 * mmesh->getIndicesVertexes().size()
1617 << std::endl;
1618 for (auto& f : mmesh->getIndicesVertexes())
1619 ostream << "3 " << f.x() << " " << f.y() << " " << f.z() << std::endl;
1620
1621 // Writing face types. Type 5 is generally triangles
1622 ostream << "\n\n";
1623 ostream << "CELL_TYPES " << mmesh->getIndicesVertexes().size() << std::endl;
1624 auto nfaces = mmesh->getIndicesVertexes().size();
1625 for (size_t j = 0; j < nfaces; j++)
1626 ostream << "5 " << std::endl;
1627
1628 outfile << ostream.str();
1629 }
1630
1631 /// get index list of neighbors
getNeighbors(int ID,std::vector<int> & neighborList)1632 void ChSystemGpu::getNeighbors(int ID, std::vector<int>& neighborList) {
1633 m_sys->getNeighbors(ID, neighborList);
1634 }
1635
1636 /// Get rolling friction torque between body i and j, return 0 if not in contact
getRollingFrictionTorque(int i,int j)1637 ChVector<float> ChSystemGpu::getRollingFrictionTorque(int i, int j) {
1638 float3 m_roll = m_sys->getRollingFrictionTorque(i, j);
1639 return ChVector<float>(m_roll.x, m_roll.y, m_roll.z);
1640 }
1641
1642 /// Get v_rot for rolling friction
getRollingVrot(int i,int j)1643 ChVector<float> ChSystemGpu::getRollingVrot(int i, int j) {
1644 float3 vrot = m_sys->getRollingVrot(i, j);
1645 return ChVector<float>(vrot.x, vrot.y, vrot.z);
1646 }
1647
1648 /// get contact char time
getRollingCharContactTime(int i,int j)1649 float ChSystemGpu::getRollingCharContactTime(int i, int j) {
1650 return m_sys->getRollingCharContactTime(i, j);
1651 }
1652
1653 /// Get tangential friction force between body i and j, return 0 if not in contact
getSlidingFrictionForce(int i,int j)1654 ChVector<float> ChSystemGpu::getSlidingFrictionForce(int i, int j) {
1655 float3 fr = m_sys->getSlidingFrictionForce(i, j);
1656 return ChVector<float>(fr.x, fr.y, fr.z);
1657 }
1658
1659 /// Get normal friction force between body i and j, return 0 if not in contact
getNormalForce(int i,int j)1660 ChVector<float> ChSystemGpu::getNormalForce(int i, int j) {
1661 float3 N = m_sys->getNormalForce(i, j);
1662 return ChVector<float>(N.x, N.y, N.z);
1663 }
1664
WriteMeshes(const std::string & outfilename) const1665 void ChSystemGpuMesh::WriteMeshes(const std::string& outfilename) const {
1666 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1667 if (sys_trimesh->file_write_mode == CHGPU_OUTPUT_MODE::NONE) {
1668 return;
1669 }
1670 if (m_meshes.size() == 0) {
1671 printf(
1672 "WARNING: attempted to write meshes to file yet no mesh found in system cache. No mesh file "
1673 "generated.\n");
1674 return;
1675 }
1676
1677 std::vector<unsigned int> vertexOffset(m_meshes.size() + 1, 0);
1678 size_t total_f = 0;
1679 size_t total_v = 0;
1680
1681 // printf("Writing %zu mesh(es)...\n", m_meshes.size());
1682 std::string ofile;
1683 if (outfilename.substr(outfilename.length() - std::min(outfilename.length(), (size_t)4)) != ".vtk" &&
1684 outfilename.substr(outfilename.length() - std::min(outfilename.length(), (size_t)4)) != ".VTK")
1685 ofile = outfilename + ".vtk";
1686 else
1687 ofile = outfilename;
1688 std::ofstream outfile(ofile, std::ios::out);
1689 std::ostringstream ostream;
1690 ostream << "# vtk DataFile Version 2.0\n";
1691 ostream << "VTK from simulation\n";
1692 ostream << "ASCII\n";
1693 ostream << "\n\n";
1694
1695 ostream << "DATASET UNSTRUCTURED_GRID\n";
1696
1697 // Prescan the V and F: to write all meshes to one file, we need vertex number offset info
1698 unsigned int mesh_num = 0;
1699 for (const auto& mmesh : m_meshes) {
1700 vertexOffset[mesh_num + 1] = (unsigned int)mmesh->getCoordsVertices().size();
1701 total_v += mmesh->getCoordsVertices().size();
1702 total_f += mmesh->getIndicesVertexes().size();
1703 mesh_num++;
1704 }
1705 for (unsigned int i = 1; i < m_meshes.size(); i++)
1706 vertexOffset[i] = vertexOffset[i] + vertexOffset[i - 1];
1707
1708 // Writing vertices
1709 ostream << "POINTS " << total_v << " float" << std::endl;
1710 mesh_num = 0;
1711 for (const auto& mmesh : m_meshes) {
1712 for (auto& v : mmesh->getCoordsVertices()) {
1713 float3 point = make_float3(v.x(), v.y(), v.z());
1714 sys_trimesh->ApplyFrameTransform(point, sys_trimesh->tri_params->fam_frame_broad[mesh_num].pos,
1715 sys_trimesh->tri_params->fam_frame_broad[mesh_num].rot_mat);
1716 ostream << point.x << " " << point.y << " " << point.z << std::endl;
1717 }
1718 mesh_num++;
1719 }
1720
1721 // Writing faces
1722 ostream << "\n\n";
1723 ostream << "CELLS " << total_f << " " << 4 * total_f << std::endl;
1724 mesh_num = 0;
1725 for (const auto& mmesh : m_meshes) {
1726 for (auto& f : mmesh->getIndicesVertexes()) {
1727 ostream << "3 " << f.x() + vertexOffset[mesh_num] << " " << f.y() + vertexOffset[mesh_num] << " "
1728 << f.z() + vertexOffset[mesh_num] << std::endl;
1729 }
1730 mesh_num++;
1731 }
1732
1733 // Writing face types. Type 5 is generally triangles
1734 ostream << "\n\n";
1735 ostream << "CELL_TYPES " << total_f << std::endl;
1736 for (const auto& mmesh : m_meshes) {
1737 auto nfaces = mmesh->getIndicesVertexes().size();
1738 for (size_t j = 0; j < nfaces; j++)
1739 ostream << "5 " << std::endl;
1740 }
1741
1742 outfile << ostream.str();
1743 }
1744
1745 // -----------------------------------------------------------------------------
1746
CollectMeshContactForces(std::vector<ChVector<>> & forces,std::vector<ChVector<>> & torques)1747 void ChSystemGpuMesh::CollectMeshContactForces(std::vector<ChVector<>>& forces, std::vector<ChVector<>>& torques) {
1748 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1749 unsigned int nmeshes = sys_trimesh->meshSoup->numTriangleFamilies;
1750 double force_factor = sys_trimesh->FORCE_SU2UU;
1751 double torque_factor = sys_trimesh->TORQUE_SU2UU;
1752
1753 forces.resize(nmeshes);
1754 torques.resize(nmeshes);
1755
1756 // Pull directly from unified memory
1757 for (unsigned int i = 0; i < nmeshes; i++) {
1758 double fx = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * i + 0];
1759 double fy = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * i + 1];
1760 double fz = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * i + 2];
1761 forces[i] = ChVector<>(fx, fy, fz) * force_factor; // Divide by C_F to go from SU to UU
1762
1763 double tx = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * i + 3];
1764 double ty = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * i + 4];
1765 double tz = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * i + 5];
1766 torques[i] = ChVector<>(tx, ty, tz) * torque_factor; // Divide by C_TAU to go from SU to UU
1767 }
1768 }
1769
CollectMeshContactForces(int mesh,ChVector<> & force,ChVector<> & torque)1770 void ChSystemGpuMesh::CollectMeshContactForces(int mesh, ChVector<>& force, ChVector<>& torque) {
1771 ChSystemGpuMesh_impl* sys_trimesh = static_cast<ChSystemGpuMesh_impl*>(m_sys);
1772 double force_factor = sys_trimesh->FORCE_SU2UU;
1773 double torque_factor = sys_trimesh->TORQUE_SU2UU;
1774
1775 double fx = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * mesh + 0];
1776 double fy = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * mesh + 1];
1777 double fz = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * mesh + 2];
1778 force = ChVector<>(fx, fy, fz) * force_factor; // Divide by C_F to go from SU to UU
1779
1780 double tx = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * mesh + 3];
1781 double ty = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * mesh + 4];
1782 double tz = sys_trimesh->meshSoup->generalizedForcesPerFamily[6 * mesh + 5];
1783 torque = ChVector<>(tx, ty, tz) * torque_factor; // Divide by C_TAU to go from SU to UU
1784 }
1785
1786 // -----------------------------------------------------------------------------
1787
1788 } // namespace gpu
1789 } // namespace chrono
1790