1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 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: Radu Serban, Hammad Mazhar, Arman Pazouki
13 // =============================================================================
14 //
15 // =============================================================================
16
17 #include "chrono/utils/ChUtilsGenerators.h"
18
19 namespace chrono {
20 namespace utils {
21
22 // =============================================================================
23 // Implementation of MixtureIngredient class
24 // =============================================================================
25
26 // Constructor: create a new mixture ingredient for the specified generator, of
27 // a given type and with a specified ratio in the final mixture.
28 // This constructor is private as a new mixture ingredient can only be created
29 // through a generator.
MixtureIngredient(Generator * generator,MixtureType type,double ratio)30 MixtureIngredient::MixtureIngredient(Generator* generator, MixtureType type, double ratio)
31 : m_generator(generator),
32 m_type(type),
33 m_ratio(ratio),
34 m_cumRatio(0),
35 m_defMaterialNSC(chrono_types::make_shared<ChMaterialSurfaceNSC>()),
36 m_defMaterialSMC(chrono_types::make_shared<ChMaterialSurfaceSMC>()),
37 m_defDensity(1),
38 m_defSize(ChVector<>(1, 1, 1)),
39 m_frictionDist(nullptr),
40 m_cohesionDist(nullptr),
41 m_youngDist(nullptr),
42 m_poissonDist(nullptr),
43 m_restitutionDist(nullptr),
44 m_densityDist(nullptr),
45 m_sizeDist(nullptr),
46 add_body_callback(nullptr) {}
47
48 // Destructor:: free the various distribution associated with this ingredient
~MixtureIngredient()49 MixtureIngredient::~MixtureIngredient() {
50 freeMaterialDist();
51 delete m_densityDist;
52 delete m_sizeDist;
53 }
54
55 // Set constant material properties for all objects based on this ingredient.
setDefaultMaterial(std::shared_ptr<ChMaterialSurface> mat)56 void MixtureIngredient::setDefaultMaterial(std::shared_ptr<ChMaterialSurface> mat) {
57 assert(mat->GetContactMethod() == m_generator->m_system->GetContactMethod());
58
59 if (mat->GetContactMethod() == ChContactMethod::NSC) {
60 m_defMaterialNSC = std::static_pointer_cast<ChMaterialSurfaceNSC>(mat);
61 } else {
62 m_defMaterialSMC = std::static_pointer_cast<ChMaterialSurfaceSMC>(mat);
63 }
64
65 freeMaterialDist();
66 }
67
setDefaultDensity(double density)68 void MixtureIngredient::setDefaultDensity(double density) {
69 m_defDensity = density;
70 delete m_densityDist;
71 m_densityDist = nullptr;
72 }
73
setDefaultSize(const ChVector<> & size)74 void MixtureIngredient::setDefaultSize(const ChVector<>& size) {
75 m_defSize = size;
76 delete m_sizeDist;
77 m_sizeDist = nullptr;
78 }
79
80 // Functions to set parameters of truncated normal distributions for the various
81 // attributes of an ingredient. Properties of objects created based on this
82 // ingredient will be obtained by drawing from the corresponding distribution.
83 // If no distribution is specified for a particular property, the default
84 // constant value is used.
setDistributionFriction(float friction_mean,float friction_stddev,float friction_min,float friction_max)85 void MixtureIngredient::setDistributionFriction(float friction_mean,
86 float friction_stddev,
87 float friction_min,
88 float friction_max) {
89 m_frictionDist = new std::normal_distribution<float>(friction_mean, friction_stddev);
90 m_minFriction = friction_min;
91 m_maxFriction = friction_max;
92 }
93
setDistributionCohesion(float cohesion_mean,float cohesion_stddev,float cohesion_min,float cohesion_max)94 void MixtureIngredient::setDistributionCohesion(float cohesion_mean,
95 float cohesion_stddev,
96 float cohesion_min,
97 float cohesion_max) {
98 m_cohesionDist = new std::normal_distribution<float>(cohesion_mean, cohesion_stddev);
99 m_minCohesion = cohesion_min;
100 m_maxCohesion = cohesion_max;
101 }
102
setDistributionYoung(float young_mean,float young_stddev,float young_min,float young_max)103 void MixtureIngredient::setDistributionYoung(float young_mean, float young_stddev, float young_min, float young_max) {
104 m_youngDist = new std::normal_distribution<float>(young_mean, young_stddev);
105 m_minYoung = young_min;
106 m_maxYoung = young_max;
107 }
108
setDistributionPoisson(float poisson_mean,float poisson_stddev,float poisson_min,float poisson_max)109 void MixtureIngredient::setDistributionPoisson(float poisson_mean,
110 float poisson_stddev,
111 float poisson_min,
112 float poisson_max) {
113 m_poissonDist = new std::normal_distribution<float>(poisson_mean, poisson_stddev);
114 m_minPoisson = poisson_min;
115 m_maxPoisson = poisson_max;
116 }
117
setDistributionRestitution(float restitution_mean,float restitution_stddev,float restitution_min,float restitution_max)118 void MixtureIngredient::setDistributionRestitution(float restitution_mean,
119 float restitution_stddev,
120 float restitution_min,
121 float restitution_max) {
122 m_restitutionDist = new std::normal_distribution<float>(restitution_mean, restitution_stddev);
123 m_minRestitution = restitution_min;
124 m_maxRestitution = restitution_max;
125 }
126
setDistributionDensity(double density_mean,double density_stddev,double density_min,double density_max)127 void MixtureIngredient::setDistributionDensity(double density_mean,
128 double density_stddev,
129 double density_min,
130 double density_max) {
131 m_densityDist = new std::normal_distribution<double>(density_mean, density_stddev);
132 m_minDensity = density_min;
133 m_maxDensity = density_max;
134 }
135
setDistributionSize(double size_mean,double size_stddev,const ChVector<> & size_min,const ChVector<> & size_max)136 void MixtureIngredient::setDistributionSize(double size_mean,
137 double size_stddev,
138 const ChVector<>& size_min,
139 const ChVector<>& size_max) {
140 m_sizeDist = new std::normal_distribution<double>(size_mean, size_stddev);
141 m_minSize = size_min;
142 m_maxSize = size_max;
143 }
144
145 // Utility function to delete all distributions associated with this ingredient.
freeMaterialDist()146 void MixtureIngredient::freeMaterialDist() {
147 delete m_frictionDist;
148 delete m_cohesionDist;
149 delete m_youngDist;
150 delete m_poissonDist;
151 delete m_restitutionDist;
152
153 m_frictionDist = nullptr;
154 m_cohesionDist = nullptr;
155 m_youngDist = nullptr;
156 m_poissonDist = nullptr;
157 m_restitutionDist = nullptr;
158 }
159
160 // Modify the specified NSC material surface based on attributes of this ingredient.
setMaterialProperties(std::shared_ptr<ChMaterialSurfaceNSC> mat)161 void MixtureIngredient::setMaterialProperties(std::shared_ptr<ChMaterialSurfaceNSC> mat) {
162 // Copy properties from the default material.
163 *mat = *m_defMaterialNSC;
164
165 // If using distributions for any of the supported properties, override those.
166 if (m_frictionDist)
167 mat->SetFriction(sampleTruncatedDist<float>(*m_frictionDist, m_minFriction, m_maxFriction));
168
169 if (m_cohesionDist)
170 mat->SetCohesion(sampleTruncatedDist<float>(*m_cohesionDist, m_minCohesion, m_maxCohesion));
171 }
172
173 // Modify the specified SMC material surface based on attributes of this ingredient.
setMaterialProperties(std::shared_ptr<ChMaterialSurfaceSMC> mat)174 void MixtureIngredient::setMaterialProperties(std::shared_ptr<ChMaterialSurfaceSMC> mat) {
175 // Copy properties from the default material.
176 *mat = *m_defMaterialSMC;
177
178 // If using distributions for any of the supported properties, override those.
179 if (m_youngDist)
180 mat->SetYoungModulus(sampleTruncatedDist<float>(*m_youngDist, m_minYoung, m_maxYoung));
181
182 if (m_poissonDist)
183 mat->SetPoissonRatio(sampleTruncatedDist<float>(*m_poissonDist, m_minPoisson, m_maxPoisson));
184
185 if (m_frictionDist)
186 mat->SetFriction(sampleTruncatedDist<float>(*m_frictionDist, m_minFriction, m_maxFriction));
187
188 if (m_restitutionDist)
189 mat->SetRestitution(sampleTruncatedDist<float>(*m_restitutionDist, m_minRestitution, m_maxRestitution));
190
191 if (m_cohesionDist)
192 mat->SetAdhesion(sampleTruncatedDist<float>(*m_cohesionDist, m_minCohesion, m_maxCohesion));
193 }
194
195 // Return a size for an object created based on attributes of this ingredient.
getSize()196 ChVector<> MixtureIngredient::getSize() {
197 if (m_sizeDist)
198 return ChVector<>(sampleTruncatedDist<double>(*m_sizeDist, m_minSize.x(), m_maxSize.x()),
199 sampleTruncatedDist<double>(*m_sizeDist, m_minSize.y(), m_maxSize.y()),
200 sampleTruncatedDist<double>(*m_sizeDist, m_minSize.z(), m_maxSize.z()));
201
202 return m_defSize;
203 }
204
205 // Return a density for an object created based on attributes of this ingredient.
getDensity()206 double MixtureIngredient::getDensity() {
207 if (m_densityDist)
208 return sampleTruncatedDist<double>(*m_densityDist, m_minDensity, m_maxDensity);
209
210 return m_defDensity;
211 }
212
213 // Calculate the volume and gyration tensor of an object of given size created
214 // from this ingredient type.
calcGeometricProps(const ChVector<> & size,double & volume,ChVector<> & gyration)215 void MixtureIngredient::calcGeometricProps(const ChVector<>& size, double& volume, ChVector<>& gyration) {
216 switch (m_type) {
217 case MixtureType::SPHERE:
218 volume = CalcSphereVolume(size.x());
219 gyration = CalcSphereGyration(size.x()).diagonal();
220 break;
221 case MixtureType::ELLIPSOID:
222 volume = CalcEllipsoidVolume(size);
223 gyration = CalcEllipsoidGyration(size).diagonal();
224 break;
225 case MixtureType::BOX:
226 volume = CalcBoxVolume(size);
227 gyration = CalcBoxGyration(size).diagonal();
228 break;
229 case MixtureType::CYLINDER:
230 volume = CalcCylinderVolume(size.x(), size.y());
231 gyration = CalcCylinderGyration(size.x(), size.y()).diagonal();
232 break;
233 case MixtureType::CONE:
234 volume = CalcConeVolume(size.x(), size.y());
235 gyration = CalcConeGyration(size.x(), size.y()).diagonal();
236 break;
237 case MixtureType::BISPHERE:
238 volume = CalcBiSphereVolume(size.x(), size.y());
239 gyration = CalcBiSphereGyration(size.x(), size.y()).diagonal();
240 break;
241 case MixtureType::CAPSULE:
242 volume = CalcCapsuleVolume(size.x(), size.y());
243 gyration = CalcCapsuleGyration(size.x(), size.y()).diagonal();
244 break;
245 case MixtureType::ROUNDEDCYLINDER:
246 volume = CalcRoundedCylinderVolume(size.x(), size.y(), size.z());
247 gyration = CalcRoundedCylinderGyration(size.x(), size.y(), size.z()).diagonal();
248 break;
249 }
250 }
251
252 // Calculate a necessary minimum separation based on the largest possible
253 // dimension of an object created based on attributes of this ingredient.
calcMinSeparation()254 double MixtureIngredient::calcMinSeparation() {
255 if (m_sizeDist)
256 return std::max(m_maxSize.x(), std::max(m_maxSize.y(), m_maxSize.z()));
257
258 return std::max(m_defSize.x(), std::max(m_defSize.y(), m_defSize.z()));
259 }
260
261 // =============================================================================
262 // Implementation of Generator class
263 // =============================================================================
264
265 // Constructor: create a generator for the specified system.
Generator(ChSystem * system)266 Generator::Generator(ChSystem* system)
267 : m_system(system), m_mixDist(0, 1), m_crtBodyId(0), m_totalNumBodies(0), m_totalMass(0), m_totalVolume(0) {
268 }
269
270 // Destructor
~Generator()271 Generator::~Generator() {
272 }
273
274 // Add a new ingredient to the current mixture by specifying its type
275 // and the ratio in the final mixture. A smart pointer to the new
276 // mixture ingredient is returned to allow modifying its properties.
AddMixtureIngredient(MixtureType type,double ratio)277 std::shared_ptr<MixtureIngredient> Generator::AddMixtureIngredient(MixtureType type, double ratio) {
278 m_mixture.push_back(chrono_types::make_shared<MixtureIngredient>(this, type, ratio));
279 return m_mixture.back();
280 }
281
282 // Create objects in a box domain using the given point sampler and separation distance and the current mixture settings.
283 // The types of objects created are selected randomly with probability
284 // proportional to the ratio of that ingredient in the mixture.
CreateObjectsBox(Sampler<double> & sampler,const ChVector<> & pos,const ChVector<> & hdims,const ChVector<> & vel)285 void Generator::CreateObjectsBox(Sampler<double>& sampler,
286 const ChVector<>& pos,
287 const ChVector<>& hdims,
288 const ChVector<>& vel) {
289 // Normalize the mixture ratios
290 normalizeMixture();
291
292 // Generate the object locations
293 double sep = sampler.GetSeparation();
294 if (m_system->GetContactMethod() == ChContactMethod::SMC)
295 sep = calcMinSeparation(sep);
296 sampler.SetSeparation(sep);
297
298 PointVector points = sampler.SampleBox(pos, hdims);
299 createObjects(points, vel);
300 }
301
302 // Create objects in a box domain sampled on a regular grid with
303 // specified separation distances and using the current mixture settings.
304 // The types of objects created are selected randomly with probability
305 // proportional to the ratio of that ingredient in the mixture.
CreateObjectsBox(const ChVector<> & dist,const ChVector<> & pos,const ChVector<> & hdims,const ChVector<> & vel)306 void Generator::CreateObjectsBox(const ChVector<>& dist,
307 const ChVector<>& pos,
308 const ChVector<>& hdims,
309 const ChVector<>& vel) {
310 // Normalize the mixture ratios
311 normalizeMixture();
312
313 // When using SMC, make sure there is no shape overlap.
314 ChVector<> distv;
315 if (m_system->GetContactMethod() == ChContactMethod::SMC)
316 distv = calcMinSeparation(dist);
317 else
318 distv = dist;
319
320 // Generate the object locations
321 GridSampler<> sampler(distv);
322 PointVector points;
323 points = sampler.SampleBox(pos, hdims);
324
325 createObjects(points, vel);
326 }
327
328 // Create objects in a cylindrical domain using the specified type of point
329 // sampler and separation distance and the current mixture settings.
330 // The types of objects created are selected randomly with probability
331 // proportional to the ratio of that ingredient in the mixture.
CreateObjectsCylinderX(Sampler<double> & sampler,const ChVector<> & pos,float radius,float halfHeight,const ChVector<> & vel)332 void Generator::CreateObjectsCylinderX(Sampler<double>& sampler,
333 const ChVector<>& pos,
334 float radius,
335 float halfHeight,
336 const ChVector<>& vel) {
337 // Normalize the mixture ratios
338 normalizeMixture();
339
340 // Generate the object locations
341 double sep = sampler.GetSeparation();
342 if (m_system->GetContactMethod() == ChContactMethod::SMC)
343 sep = calcMinSeparation(sep);
344 sampler.SetSeparation(sep);
345
346 PointVector points = sampler.SampleCylinderX(pos, radius, halfHeight);
347 createObjects(points, vel);
348 }
349
CreateObjectsCylinderY(Sampler<double> & sampler,const ChVector<> & pos,float radius,float halfHeight,const ChVector<> & vel)350 void Generator::CreateObjectsCylinderY(Sampler<double>& sampler,
351 const ChVector<>& pos,
352 float radius,
353 float halfHeight,
354 const ChVector<>& vel) {
355 // Normalize the mixture ratios
356 normalizeMixture();
357
358 // Generate the object locations
359 double sep = sampler.GetSeparation();
360 if (m_system->GetContactMethod() == ChContactMethod::SMC)
361 sep = calcMinSeparation(sep);
362 sampler.SetSeparation(sep);
363
364 PointVector points = sampler.SampleCylinderY(pos, radius, halfHeight);
365 createObjects(points, vel);
366 }
367
CreateObjectsCylinderZ(Sampler<double> & sampler,const ChVector<> & pos,float radius,float halfHeight,const ChVector<> & vel)368 void Generator::CreateObjectsCylinderZ(Sampler<double>& sampler,
369 const ChVector<>& pos,
370 float radius,
371 float halfHeight,
372 const ChVector<>& vel) {
373 // Normalize the mixture ratios
374 normalizeMixture();
375
376 // Generate the object locations
377 double sep = sampler.GetSeparation();
378 if (m_system->GetContactMethod() == ChContactMethod::SMC)
379 sep = calcMinSeparation(sep);
380 sampler.SetSeparation(sep);
381
382 PointVector points = sampler.SampleCylinderZ(pos, radius, halfHeight);
383 createObjects(points, vel);
384 }
385
CreateObjectsSphere(Sampler<double> & sampler,const ChVector<> & pos,float radius,const ChVector<> & vel)386 void Generator::CreateObjectsSphere(Sampler<double>& sampler,
387 const ChVector<>& pos,
388 float radius,
389 const ChVector<>& vel) {
390 // Normalize the mixture ratios
391 normalizeMixture();
392
393 // Generate the object locations
394 double sep = sampler.GetSeparation();
395 if (m_system->GetContactMethod() == ChContactMethod::SMC)
396 sep = calcMinSeparation(sep);
397 sampler.SetSeparation(sep);
398
399 PointVector points = sampler.SampleSphere(pos, radius);
400 createObjects(points, vel);
401 }
402
403 // Normalize the mixture ratios (so that their sum is 1) and calculate
404 // the exclusive scan of these ratios.
normalizeMixture()405 void Generator::normalizeMixture() {
406 if (m_mixture.empty()) {
407 AddMixtureIngredient(MixtureType::SPHERE, 1);
408 return;
409 }
410
411 double sum = 0;
412
413 for (int i = 0; i < m_mixture.size(); i++)
414 sum += m_mixture[i]->m_ratio;
415
416 for (int i = 0; i < m_mixture.size(); i++) {
417 m_mixture[i]->m_ratio /= sum;
418 if (i == 0)
419 m_mixture[i]->m_cumRatio = 0;
420 else
421 m_mixture[i]->m_cumRatio = m_mixture[i - 1]->m_cumRatio + m_mixture[i - 1]->m_ratio;
422 }
423 }
424
425 // Select one of the mixture types with probability proportional to that
426 // type's ratio in the mixture.
selectIngredient()427 int Generator::selectIngredient() {
428 double val = m_mixDist(rengine());
429
430 for (int i = (int)m_mixture.size() - 1; i >= 0; i--) {
431 if (val > m_mixture[i]->m_cumRatio)
432 return i;
433 }
434
435 return 0;
436 }
437
438 // Calculate the minimum distance to ensure contact shapes never overlap.
calcMinSeparation(double sep)439 double Generator::calcMinSeparation(double sep) {
440 for (int i = 0; i < m_mixture.size(); i++)
441 sep = std::max(sep, m_mixture[i]->calcMinSeparation());
442
443 return sep;
444 }
445
calcMinSeparation(const ChVector<> & sep)446 ChVector<> Generator::calcMinSeparation(const ChVector<>& sep) {
447 ChVector<> res;
448
449 for (int i = 0; i < m_mixture.size(); i++) {
450 double mix_sep = m_mixture[i]->calcMinSeparation();
451 res.x() = std::max(sep.x(), mix_sep);
452 res.y() = std::max(sep.y(), mix_sep);
453 res.z() = std::max(sep.z(), mix_sep);
454 }
455
456 return res;
457 }
458
459 // Create objects at the specified locations using the current mixture settings.
createObjects(const PointVector & points,const ChVector<> & vel)460 void Generator::createObjects(const PointVector& points, const ChVector<>& vel) {
461 bool check = false;
462 std::vector<bool> flags;
463 if (m_callback) {
464 flags.resize(points.size(), true);
465 m_callback->OnCreateObjects(points, flags);
466 check = true;
467 }
468
469 for (int i = 0; i < points.size(); i++) {
470 if (check && !flags[i])
471 continue;
472
473 // Select the type of object to be created.
474 int index = selectIngredient();
475
476 // Create a contact material consistent with the associated system and modify it based on attributes of the
477 // current ingredient.
478 std::shared_ptr<ChMaterialSurface> mat;
479 switch (m_system->GetContactMethod()) {
480 case ChContactMethod::NSC: {
481 auto matNSC = chrono_types::make_shared<ChMaterialSurfaceNSC>();
482 m_mixture[index]->setMaterialProperties(matNSC);
483 mat = matNSC;
484 break;
485 }
486 case ChContactMethod::SMC: {
487 auto matSMC = chrono_types::make_shared<ChMaterialSurfaceSMC>();
488 m_mixture[index]->setMaterialProperties(matSMC);
489 mat = matSMC;
490 break;
491 }
492 }
493
494 // Create the body (with appropriate collision model, consistent with the associated system)
495 ChBody* body = m_system->NewBody();
496
497 // Set identifier
498 body->SetIdentifier(m_crtBodyId++);
499
500 // Set position and orientation
501 body->SetPos(points[i]);
502 body->SetRot(ChQuaternion<>(1, 0, 0, 0));
503 body->SetPos_dt(vel);
504 body->SetBodyFixed(false);
505 body->SetCollide(true);
506
507 // Get size and density; calculate geometric properties
508 ChVector<> size = m_mixture[index]->getSize();
509 double density = m_mixture[index]->getDensity();
510 double volume;
511 ChVector<> gyration;
512 m_mixture[index]->calcGeometricProps(size, volume, gyration);
513 double mass = density * volume;
514
515 // Set mass properties
516 body->SetMass(mass);
517 body->SetInertiaXX(mass * gyration);
518 m_totalMass += mass;
519 m_totalVolume += volume;
520
521 // Add collision geometry
522 body->GetCollisionModel()->ClearModel();
523
524 switch (m_mixture[index]->m_type) {
525 case MixtureType::SPHERE:
526 AddSphereGeometry(body, mat, size.x());
527 break;
528 case MixtureType::ELLIPSOID:
529 AddEllipsoidGeometry(body, mat, size);
530 break;
531 case MixtureType::BOX:
532 AddBoxGeometry(body, mat, size);
533 break;
534 case MixtureType::CYLINDER:
535 AddCylinderGeometry(body, mat, size.x(), size.y());
536 break;
537 case MixtureType::CONE:
538 AddConeGeometry(body, mat, size.x(), size.y());
539 break;
540 case MixtureType::BISPHERE:
541 AddBiSphereGeometry(body, mat, size.x(), size.y());
542 break;
543 case MixtureType::CAPSULE:
544 AddCapsuleGeometry(body, mat, size.x(), size.y());
545 break;
546 case MixtureType::ROUNDEDCYLINDER:
547 AddRoundedCylinderGeometry(body, mat, size.x(), size.y(), size.z());
548 break;
549 }
550
551 body->GetCollisionModel()->BuildModel();
552
553 // Attach the body to the system and append to list of generated bodies.
554 std::shared_ptr<ChBody> bodyPtr(body);
555
556 m_system->AddBody(bodyPtr);
557
558 // If the callback pointer is set, call the function with the body pointer
559 if (m_mixture[index]->add_body_callback) {
560 m_mixture[index]->add_body_callback->OnAddBody(bodyPtr);
561 }
562
563 m_bodies.push_back(BodyInfo(m_mixture[index]->m_type, density, size, bodyPtr));
564 }
565
566 m_totalNumBodies += (unsigned int)points.size();
567 }
568
569 // Write body information to a CSV file
writeObjectInfo(const std::string & filename)570 void Generator::writeObjectInfo(const std::string& filename) {
571 CSV_writer csv;
572
573 for (int i = 0; i < m_bodies.size(); i++) {
574 csv << static_cast<int>(m_bodies[i].m_type);
575 csv << m_bodies[i].m_body->GetPos() << m_bodies[i].m_size;
576 csv << m_bodies[i].m_density << m_bodies[i].m_body->GetMass();
577
578 //// RADU: write collision shape information, including contact material properties?
579
580 csv << std::endl;
581 }
582
583 csv.write_to_file(filename);
584 }
585
586 } // namespace utils
587 } // namespace chrono
588