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: Conlain Kelly, Nic Olsen, Ruochun Zhang, Dan Negrut, Radu Serban
13 // =============================================================================
14
15 #include <vector>
16 #include <algorithm>
17 #include <fstream>
18 #include <string>
19 #include <cmath>
20
21 #include "chrono/core/ChGlobal.h"
22 #include "chrono/core/ChVector.h"
23 #include "chrono/core/ChQuaternion.h"
24 #include "chrono/core/ChMatrix33.h"
25
26 #include "chrono_gpu/physics/ChSystemGpuMesh_impl.h"
27 #include "chrono_gpu/utils/ChGpuUtilities.h"
28
29 namespace chrono {
30 namespace gpu {
31
ChSystemGpuMesh_impl(float sphere_rad,float density,float3 boxDims,float3 O)32 ChSystemGpuMesh_impl::ChSystemGpuMesh_impl(float sphere_rad, float density, float3 boxDims, float3 O)
33 : ChSystemGpu_impl(sphere_rad, density, boxDims, O),
34 K_n_s2m_UU(0),
35 K_t_s2m_UU(0),
36 Gamma_n_s2m_UU(0),
37 Gamma_t_s2m_UU(0),
38 YoungsModulus_mesh_UU(0.0),
39 COR_mesh_UU(0.0),
40 use_mat_based(false),
41 PoissonRatio_mesh_UU(0.0),
42 rolling_coeff_s2m_UU(0),
43 spinning_coeff_s2m_UU(0),
44 adhesion_s2m_over_gravity(0) {
45 // Allocate triangle collision parameters
46 gpuErrchk(cudaMallocManaged(&tri_params, sizeof(MeshParams), cudaMemAttachGlobal));
47
48 // Allocate the device soup storage
49 gpuErrchk(cudaMallocManaged(&meshSoup, sizeof(TriangleSoup), cudaMemAttachGlobal));
50 // start with no triangles
51 meshSoup->nTrianglesInSoup = 0;
52 meshSoup->numTriangleFamilies = 0;
53
54 tri_params->static_friction_coeff_s2m = 0;
55 }
56
~ChSystemGpuMesh_impl()57 ChSystemGpuMesh_impl::~ChSystemGpuMesh_impl() {
58 // work to do here
59 cleanupTriMesh();
60 }
get_max_K() const61 double ChSystemGpuMesh_impl::get_max_K() const {
62 double maxK;
63 if (tri_params->use_mat_based == true) {
64 // material pparameter sigma for different surface
65 // see reference eq 2.13, 2.14 in Book Contact Force Model, Flores and Lankarani
66 double sigma_sphere = (1 - std::pow(PoissonRatio_sphere_UU, 2)) / YoungsModulus_sphere_UU;
67 double sigma_wall = (1 - std::pow(PoissonRatio_wall_UU, 2)) / YoungsModulus_wall_UU;
68 double sigma_mesh = (1 - std::pow(PoissonRatio_mesh_UU, 2)) / YoungsModulus_mesh_UU;
69
70 maxK = 4.0 / (3.0 * (sigma_sphere + std::min(std::min(sigma_sphere, sigma_wall), sigma_mesh))) *
71 std::sqrt(sphere_radius_UU);
72 INFO_PRINTF("Use material based contact force model, maximum effective stiffnes is %e\n", maxK);
73 return maxK;
74
75 } else {
76 maxK = std::max(std::max(K_n_s2s_UU, K_n_s2w_UU), K_n_s2m_UU);
77 INFO_PRINTF("Use user defined contact force model, maximum effective stiffnes is %e\n", maxK);
78 return maxK;
79 }
80 }
81
combineMaterialSurface()82 void ChSystemGpuMesh_impl::combineMaterialSurface() {
83 // effective youngs modulus and shear modulus in user units
84 double E_eff_s2m_uu, G_eff_s2m_uu;
85
86 materialPropertyCombine(YoungsModulus_sphere_UU, YoungsModulus_mesh_UU, PoissonRatio_sphere_UU,
87 PoissonRatio_mesh_UU, E_eff_s2m_uu, G_eff_s2m_uu);
88
89 // unit conversion of youngs modulus
90 double E_SU2UU = this->MASS_SU2UU / (this->LENGTH_SU2UU * this->TIME_SU2UU * this->TIME_SU2UU);
91
92 // convert effective youngs modulus and shear modulus to simulation units
93 tri_params->E_eff_s2m_SU = (float)E_eff_s2m_uu / E_SU2UU;
94 tri_params->G_eff_s2m_SU = (float)G_eff_s2m_uu / E_SU2UU;
95
96 // assign coefficient of restitution in simulation units
97 tri_params->COR_s2m_SU = std::min(COR_sphere_UU, COR_mesh_UU);
98 }
99
initializeTriangles()100 void ChSystemGpuMesh_impl::initializeTriangles() {
101 tri_params->use_mat_based = use_mat_based;
102
103 if (use_mat_based == false) {
104 double K_SU2UU = MASS_SU2UU / (TIME_SU2UU * TIME_SU2UU);
105 double GAMMA_SU2UU = 1. / TIME_SU2UU;
106
107 tri_params->K_n_s2m_SU = (float)(K_n_s2m_UU / K_SU2UU);
108 tri_params->K_t_s2m_SU = (float)(K_t_s2m_UU / K_SU2UU);
109
110 tri_params->Gamma_n_s2m_SU = (float)(Gamma_n_s2m_UU / GAMMA_SU2UU);
111 tri_params->Gamma_t_s2m_SU = (float)(Gamma_t_s2m_UU / GAMMA_SU2UU);
112
113 } else {
114 combineMaterialSurface();
115 }
116
117 double magGravAcc = sqrt(X_accGrav * X_accGrav + Y_accGrav * Y_accGrav + Z_accGrav * Z_accGrav);
118 tri_params->adhesionAcc_s2m = (float)(adhesion_s2m_over_gravity * magGravAcc);
119
120 for (unsigned int fam = 0; fam < meshSoup->numTriangleFamilies; fam++) {
121 meshSoup->familyMass_SU[fam] = (float)(meshSoup->familyMass_SU[fam] / MASS_SU2UU);
122 }
123
124 tri_params->rolling_coeff_s2m_SU = (float)rolling_coeff_s2m_UU;
125
126 const double meshRot[4] = {1., 0., 0., 0.};
127 for (unsigned int fam = 0; fam < meshSoup->numTriangleFamilies; fam++) {
128 generate_rot_matrix<float>(meshRot, tri_params->fam_frame_broad[fam].rot_mat);
129 tri_params->fam_frame_broad[fam].pos[0] = (float)0.0;
130 tri_params->fam_frame_broad[fam].pos[1] = (float)0.0;
131 tri_params->fam_frame_broad[fam].pos[2] = (float)0.0;
132
133 generate_rot_matrix<double>(meshRot, tri_params->fam_frame_narrow[fam].rot_mat);
134 tri_params->fam_frame_narrow[fam].pos[0] = (double)0.0;
135 tri_params->fam_frame_narrow[fam].pos[1] = (double)0.0;
136 tri_params->fam_frame_narrow[fam].pos[2] = (double)0.0;
137 }
138
139 TRACK_VECTOR_RESIZE(SD_numTrianglesTouching, nSDs, "SD_numTrianglesTouching", 0);
140 TRACK_VECTOR_RESIZE(SD_TrianglesCompositeOffsets, nSDs, "SD_TrianglesCompositeOffsets", 0);
141
142 TRACK_VECTOR_RESIZE(Triangle_NumSDsTouching, meshSoup->nTrianglesInSoup, "Triangle_NumSDsTouching", 0);
143 TRACK_VECTOR_RESIZE(Triangle_SDsCompositeOffsets, meshSoup->nTrianglesInSoup, "Triangle_SDsCompositeOffsets", 0);
144
145 // TODO do we have a good heuristic???
146 // this gets resized on-the-fly every timestep
147 TRACK_VECTOR_RESIZE(SD_trianglesInEachSD_composite, 0, "SD_trianglesInEachSD_composite", 0);
148 }
149
150 // p = pos + rot_mat * p
ApplyFrameTransform(float3 & p,float * pos,float * rot_mat)151 void ChSystemGpuMesh_impl::ApplyFrameTransform(float3& p, float* pos, float* rot_mat) {
152 float3 result;
153
154 // Apply rotation matrix to point
155 result.x = pos[0] + rot_mat[0] * p.x + rot_mat[1] * p.y + rot_mat[2] * p.z;
156 result.y = pos[1] + rot_mat[3] * p.x + rot_mat[4] * p.y + rot_mat[5] * p.z;
157 result.z = pos[2] + rot_mat[6] * p.x + rot_mat[7] * p.y + rot_mat[8] * p.z;
158
159 // overwrite p only at the end
160 p = result;
161 }
162
cleanupTriMesh()163 void ChSystemGpuMesh_impl::cleanupTriMesh() {
164 cudaFree(meshSoup->triangleFamily_ID);
165 cudaFree(meshSoup->familyMass_SU);
166
167 cudaFree(meshSoup->node1);
168 cudaFree(meshSoup->node2);
169 cudaFree(meshSoup->node3);
170
171 cudaFree(meshSoup->vel);
172 cudaFree(meshSoup->omega);
173
174 cudaFree(meshSoup->generalizedForcesPerFamily);
175 cudaFree(tri_params->fam_frame_broad);
176 cudaFree(tri_params->fam_frame_narrow);
177 cudaFree(meshSoup);
178 cudaFree(tri_params);
179 }
180
ApplyMeshMotion(unsigned int mesh_id,const double * pos,const double * rot,const double * lin_vel,const double * ang_vel)181 void ChSystemGpuMesh_impl::ApplyMeshMotion(unsigned int mesh_id,
182 const double* pos,
183 const double* rot,
184 const double* lin_vel,
185 const double* ang_vel) {
186 // Set position and orientation
187 tri_params->fam_frame_broad[mesh_id].pos[0] = (float)pos[0];
188 tri_params->fam_frame_broad[mesh_id].pos[1] = (float)pos[1];
189 tri_params->fam_frame_broad[mesh_id].pos[2] = (float)pos[2];
190 generate_rot_matrix<float>(rot, tri_params->fam_frame_broad[mesh_id].rot_mat);
191
192 tri_params->fam_frame_narrow[mesh_id].pos[0] = pos[0];
193 tri_params->fam_frame_narrow[mesh_id].pos[1] = pos[1];
194 tri_params->fam_frame_narrow[mesh_id].pos[2] = pos[2];
195 generate_rot_matrix<double>(rot, tri_params->fam_frame_narrow[mesh_id].rot_mat);
196
197 // Set linear and angular velocity
198 const float C_V = (float)(TIME_SU2UU / LENGTH_SU2UU);
199 meshSoup->vel[mesh_id] = make_float3(C_V * (float)lin_vel[0], C_V * (float)lin_vel[1], C_V * (float)lin_vel[2]);
200
201 const float C_O = (float)TIME_SU2UU;
202 meshSoup->omega[mesh_id] = make_float3(C_O * (float)ang_vel[0], C_O * (float)ang_vel[1], C_O * (float)ang_vel[2]);
203 }
204
205 template <typename T>
generate_rot_matrix(const double * ep,T * rot_mat)206 void ChSystemGpuMesh_impl::generate_rot_matrix(const double* ep, T* rot_mat) {
207 rot_mat[0] = (T)(2 * (ep[0] * ep[0] + ep[1] * ep[1] - 0.5));
208 rot_mat[1] = (T)(2 * (ep[1] * ep[2] - ep[0] * ep[3]));
209 rot_mat[2] = (T)(2 * (ep[1] * ep[3] + ep[0] * ep[2]));
210 rot_mat[3] = (T)(2 * (ep[1] * ep[2] + ep[0] * ep[3]));
211 rot_mat[4] = (T)(2 * (ep[0] * ep[0] + ep[2] * ep[2] - 0.5));
212 rot_mat[5] = (T)(2 * (ep[2] * ep[3] - ep[0] * ep[1]));
213 rot_mat[6] = (T)(2 * (ep[1] * ep[3] - ep[0] * ep[2]));
214 rot_mat[7] = (T)(2 * (ep[2] * ep[3] + ep[0] * ep[1]));
215 rot_mat[8] = (T)(2 * (ep[0] * ep[0] + ep[3] * ep[3] - 0.5));
216 }
217
218 } // namespace gpu
219 } // namespace chrono
220