tsphere.cpp - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tsphere.cpp (28489B)
---
1 #include <iostream>
2 #include <string>
3 #include <cstdio>
4 #include <cstdlib>
5 #include <cmath>
6 #include <vector>
7 #include <algorithm>
8
9 #include "typedefs.h"
10 #include "datatypes.h"
11 #include "constants.h"
12 #include "sphere.h"
13
14 // Constructor: Reads an input binary, and optionally checks
15 // and reports the values
16 DEM::DEM(const std::string inputbin,
17 const int verbosity,
18 const int checkVals,
19 const int dry,
20 const int initCuda,
21 const int transferConstMem,
22 const int fluidFlow,
23 const int device)
24 : verbose(verbosity), device(device), fluid(fluidFlow)
25 {
26 using std::cout;
27 using std::cerr;
28
29 // Extract sid from input binary filename
30 size_t dotpos = inputbin.rfind('.');
31 size_t slashpos = inputbin.rfind('/');
32 if (slashpos - dotpos < 1) {
33 std::cerr << "Error! Unable to extract simulation id "
34 << "from input file name.\n";
35 }
36 sid = inputbin.substr(slashpos+1, dotpos-slashpos-1);
37
38 // Read target input binary
39 readbin(inputbin.c_str());
40
41 // Check numeric values of chosen parameters
42 if (checkVals == 1)
43 checkValues();
44
45 // Report data values
46 if (dry == 1)
47 reportValues();
48
49 // If this is a dry run, exit
50 if (dry == 1)
51 exit(0);
52
53 if (fluid == 1) {
54 if (cfd_solver == 0)
55 initNS();
56 else if (cfd_solver == 1)
57 initDarcy();
58 else {
59 std::cerr << "DEM::DEM Error: Value of cfd_solver not understood ("
60 << cfd_solver << ")" << std::endl;
61 exit(1);
62 }
63 }
64
65 if (initCuda == 1) {
66
67 // Initialize CUDA
68 initializeGPU();
69
70 if (transferConstMem == 1) {
71 // Copy constant data to constant device memory
72 transferToConstantDeviceMemory();
73 }
74
75 if (fluid == 1) {
76 if (cfd_solver == 0)
77 initNSmemDev();
78 else if (cfd_solver == 1)
79 initDarcyMemDev();
80 else {
81 std::cerr
82 << "DEM::DEM Error: Value of cfd_solver not understood ("
83 << cfd_solver << ")" << std::endl;
84 exit(1);
85 }
86 }
87
88 // Allocate device memory for particle variables,
89 // tied to previously declared pointers in structures
90 allocateGlobalDeviceMemory();
91
92 // Transfer data from host to gpu device memory
93 transferToGlobalDeviceMemory();
94 }
95 }
96
97 // Destructor: Liberates dynamically allocated host memory
98 DEM::~DEM(void)
99 {
100 if (verbose == 1)
101 std::cout << "Freeing host memory: ";
102 delete[] k.x;
103 delete[] k.xyzsum;
104 delete[] k.vel;
105 delete[] k.force;
106 delete[] k.angpos;
107 delete[] k.angvel;
108 delete[] k.torque;
109 delete[] k.color;
110 delete[] e.es_dot;
111 delete[] e.es;
112 delete[] e.ev_dot;
113 delete[] e.ev;
114 delete[] e.p;
115 delete[] walls.nx;
116 delete[] walls.mvfd;
117 if (verbose == 1)
118 std::cout << "Done" << std::endl;
119 }
120
121 void checkIfNaN(Float3 vec, std::string description, unsigned int idx)
122 {
123 if (vec.x != vec.x || vec.y != vec.y || vec.z != vec.z) {
124 std::cerr << "Error: Particle " << idx << " has a "
125 << description << " with one or more NaN values: "
126 << vec.x << ", " << vec.y << ", " << vec.z << std::endl;
127 exit(1);
128 }
129 }
130
131 void checkIfNaN(Float4 vec, std::string description, unsigned int idx)
132 {
133 if (vec.x != vec.x || vec.y != vec.y ||
134 vec.z != vec.z || vec.w != vec.w) {
135 std::cerr << "Error: Particle " << idx << " has a "
136 << description << " with one or more NaN values: "
137 << vec.x << ", " << vec.y << ", " << vec.z << ", " << vec.w
138 << std::endl;
139 exit(1);
140 }
141 }
142
143
144 // Check numeric values of selected parameters
145 void DEM::checkValues(void)
146 {
147 using std::cerr;
148 using std::cout;
149 using std::endl;
150
151 unsigned int i;
152
153 // Check the number of dimensions
154 if (nd != ND) {
155 cerr << "Error: nd = " << nd << ", ND = " << ND << endl;
156 exit(1);
157 }
158
159 // Check the number of possible contacts
160 if (NC < 1) {
161 cerr << "Error: NC = " << NC << '\n';
162 exit(1);
163 } else if (NC < 8) {
164 cerr << "Warning: NC has a low value (" << NC << "). "
165 << "Consider increasing it in 'constants.h'" << endl;
166 }
167
168 // Check that we have a positive number of particles
169 if (np < 1 && verbose == 1) {
170 cout << "Info: No particles are being simulated (np = " << np
171 << ")" << endl;
172 }
173
174 // Check that the current time
175 if (time.current > time.total || time.current < 0.0) {
176 cerr << "Error: time.current = " << time.current
177 << " s, time.total = " << time.total << " s" << endl;
178 exit(1);
179 }
180
181 // Check origo placement
182 if (grid.origo[0] < 0.0 || grid.origo[1] < 0.0 || grid.origo[2] < 0.0) {
183 cerr << "Error: Negative values in grid.origo are known to cause "
184 << "problems. \n"
185 << "grid.origo[0] = " << grid.origo[0] << " m, "
186 << "grid.origo[1] = " << grid.origo[1] << " m, "
187 << "grid.origo[2] = " << grid.origo[2] << " m." << endl;
188 exit(1);
189 }
190
191 // Check world size
192 if (grid.L[0] <= 0.0 || grid.L[1] <= 0.0 || grid.L[2] <= 0.0) {
193 cerr << "Error: grid.L[0] = " << grid.L[0] << " m, "
194 << "grid.L[1] = " << grid.L[1] << " m, "
195 << "grid.L[2] = " << grid.L[2] << " m." << endl;
196 exit(1);
197 }
198
199 // Check grid size
200 if (grid.num[0] <= 0 || grid.num[1] <= 0 || grid.num[2] <= 0) {
201 cerr << "Error: grid.num[0] = " << grid.num[0] << ", "
202 << "grid.num[1] = " << grid.num[1] << ", "
203 << "grid.num[2] = " << grid.num[2] << "." << endl;
204 exit(1);
205 }
206
207 // Check grid size again
208 if (grid.periodic == 2 && grid.num[0] < 3) {
209 cerr << "Error: When 1st dimension boundaries are periodic, "
210 << "there must be at least 3 cells in that dimension." << endl;
211 exit(1);
212 }
213
214 if (grid.periodic == 1 && (grid.num[0] < 3 || grid.num[1] < 3)) {
215 cerr << "Error: When 1st and 2nd dimension boundaries are periodic, "
216 << "there must be at least 3 cells in each of those dimensions."
217 << endl;
218 exit(1);
219 }
220
221
222 // Per-particle checks
223 Float4 x, vel, acc, angpos, angvel, angacc;
224 for (i = 0; i < np; ++i) {
225
226 // Read values. Accelerations can't be checked by default, since these
227 // aren't initialized when this function is run after reading the input
228 // file from disk
229 x = k.x[i];
230 vel = k.vel[i];
231 //acc = k.acc[i];
232 angpos = k.angpos[i];
233 angvel = k.angvel[i];
234 //angacc = k.angacc[i];
235
236 // Check that radii are positive values
237 if (x.w <= 0.0) {
238 cerr << "Error: Particle " << i << " has a radius of "
239 << k.x[i].w << " m." << endl;
240 exit(1);
241 }
242
243 checkIfNaN(x, "position", i);
244 checkIfNaN(vel, "velocity", i);
245 //checkIfNaN(angpos, "angular position", i);
246 //checkIfNaN(angvel, "angular velocity", i);
247
248 // Check that all particles are inside of the grid
249 if (x.x < grid.origo[0] ||
250 x.y < grid.origo[1] ||
251 x.z < grid.origo[2] ||
252 x.x > grid.L[0] ||
253 x.y > grid.L[1] ||
254 x.z > grid.L[2]) {
255 cerr << "Error: Particle " << i << " is outside of "
256 << "the computational grid\n"
257 << "k.x[i] = ["
258 << x.x << ", "
259 << x.y << ", "
260 << x.z << "]\n"
261 << "grid.origo = ["
262 << grid.origo[0] << ", "
263 << grid.origo[1] << ", "
264 << grid.origo[2] << "], "
265 << "grid.L = ["
266 << grid.L[0] << ", "
267 << grid.L[1] << ", "
268 << grid.L[2] << "]."
269 << endl;
270 exit(1);
271
272 }
273 }
274
275 // If present, check that the upper wall is above all particles
276 if (walls.nw > 0) {
277 Float z_max = 0.0;
278 for (i = 0; i < np; ++i) {
279 if (k.x[i].z > z_max)
280 z_max = k.x[i].z;
281 }
282
283 if (walls.nx[0].w < z_max) {
284 cerr << "Error: One or more particles have centres above "
285 << "the upper, dynamic wall" << endl;
286 exit(1);
287 }
288 }
289
290
291
292 // Check constant, global parameters
293 if (params.k_n <= 0.0) {
294 cerr << "Error: k_n = " << params.k_n << " N/m" << endl;
295 exit(1);
296 }
297
298 if (params.E < 0.0) {
299 cerr << "Error: E = " << params.E << " N/m^3" << endl;
300 exit(1);
301 }
302
303 if (params.rho <= 0.0) {
304 cerr << "Error: rho = " << params.rho << " kg/m3" << endl;
305 exit(1);
306 }
307
308 if (fluid == 1) {
309
310 // Navier-Stokes tests
311 if (cfd_solver == 0) {
312 if (ns.rho_f <= 0.0) {
313 cerr << "Error: rho = " << params.rho << " kg/m3" << endl;
314 exit(1);
315 }
316 }
317
318 // Darcy tests
319 else if (cfd_solver == 1) {
320 if (darcy.rho_f <= 0.0) {
321 cerr << "Error: rho_f = " << darcy.rho_f << " kg/m3" << endl;
322 exit(1);
323 }
324
325 if (darcy.mu <= 0.0) {
326 cerr << "Error: mu = " << darcy.mu << " Pa s" << endl;
327 exit(1);
328 }
329
330 if (darcy.beta_f <= 0.0) {
331 cerr << "Error: beta_f = " << darcy.beta_f << " 1/Pa" << endl;
332 exit(1);
333 }
334
335 if (darcy.k_c <= 0.0) {
336 cerr << "Error: k_c = " << darcy.k_c << " m*m" << endl;
337 exit(1);
338 }
339 } else {
340 cerr << "Solver type value not understood (cfd_solver = "
341 << cfd_solver << endl;
342 exit(1);
343 }
344 }
345 }
346
347
348 // Report key parameter values to stdout
349 void DEM::reportValues()
350 {
351 using std::cout;
352 using std::cerr;
353 using std::endl;
354
355 cout << " - Number of dimensions: nd = " << nd << '\n'
356 << " - Number of particles: np = " << np << endl;
357
358 // Check precision choice
359 cout << " - Compiled for ";
360 if (sizeof(Float) == sizeof(float)) {
361 if (verbose == 1)
362 cout << "single";
363 } else if (sizeof(Float) == sizeof(double)) {
364 if (verbose == 1)
365 cout << "double";
366 } else {
367 cerr << "Error! Chosen precision not available. Check datatypes.h"
368 << endl;
369 exit(1);
370 }
371 cout << " precision\n";
372
373 cout << " - Timestep length: time.dt = "
374 << time.dt << " s\n"
375 << " - Start at time: time.current = "
376 << time.current << " s\n"
377 << " - Total sim. time: time.total = "
378 << time.total << " s\n"
379 << " - File output interval: time.file_dt = "
380 << time.file_dt << " s\n"
381 << " - Start at step count: time.step_count = "
382 << time.step_count << endl;
383
384 if (params.contactmodel == 1)
385 cout << " - Contact model: Linear-elastic-viscous (n), visco-frictional (t)\n";
386 else if (params.contactmodel == 2)
387 cout << " - Contact model: Linear-elastic-visco-frictional\n";
388 else if (params.contactmodel == 3)
389 cout << " - Contact model: Nonlinear-elastic-visco-frictional\n";
390 else {
391 cerr << "Error: Contact model value (" << params.contactmodel << ") not understood.\n";
392 exit(1);
393 }
394
395 if (params.E > 0.01)
396 cout << " - Using Young's modulus for contact stiffness\n";
397 else
398 cout << " - Using global value for contact stiffness\n";
399
400 if (params.kappa > 0.0 && params.V_b > 0.0)
401 cout << " - Capillary cohesion enabled\n";
402
403 cout << " - Number of dynamic walls: " << walls.nw << '\n';
404
405 if (grid.periodic == 1)
406 cout << " - 1st and 2nd dim. boundaries: Periodic\n";
407 else if (grid.periodic == 2)
408 cout << " - 1st dim. boundaries: Visco-frictional walls\n";
409 else
410 cout << " - 1st and 2nd dim. boundaries: Visco-frictional walls\n";
411
412 if (walls.nw > 0) {
413 cout << " - Top BC: ";
414 if (walls.wmode[0] == 0)
415 cout << "Fixed\n";
416 else if (walls.wmode[0] == 1)
417 cout << "Deviatoric stress, "
418 << walls.mvfd[0].w << " Pa\n";
419 else if (walls.wmode[0] == 2)
420 cout << "Velocity, "
421 << walls.mvfd[0].y << " m/s\n";
422 else if (walls.wmode[0] == 3)
423 cout << "Shear stress, "
424 << walls.tau_x[0] << " Pa\n";
425 else {
426 cerr << "Top BC not recognized!\n";
427 exit(1);
428 }
429 }
430
431 cout << " - Grid: ";
432 if (nd == 1)
433 cout << grid.num[0];
434 else if (nd == 2)
435 cout << grid.num[0] << " * " << grid.num[1];
436 else
437 cout << grid.num[0] << " * "
438 << grid.num[1] << " * "
439 << grid.num[2];
440 cout << " cells\n";
441 cout << " - Domain size: ";
442 if (nd == 1)
443 cout << grid.L[0];
444 else if (nd == 2)
445 cout << grid.L[0] << " * " << grid.L[1];
446 else
447 cout << grid.L[0] << " * "
448 << grid.L[1] << " * "
449 << grid.L[2];
450 cout << " m\n";
451
452 cout << " - No. of particle bonds: " << params.nb0 << endl;
453
454 if (fluid == 1) {
455 cout << " - Fluid solver: ";
456 if (cfd_solver == 0) {
457 cout << "Navier-Stokes";
458 } else if (cfd_solver == 1) {
459 cout << "Darcy";
460 } else {
461 cout << "Unknown";
462 }
463 cout << endl;
464 }
465 }
466
467 // Returns the volume of a spherical cap
468 Float sphericalCap(const Float h, const Float r)
469 {
470 return M_PI * h * h / 3.0 * (3.0 * r - h);
471 }
472
473 // Returns the max. radius of any particle
474 Float DEM::r_max()
475 {
476 Float r_max = 0.0;
477 Float r;
478 for (unsigned int i=0; i<np; ++i) {
479 r = k.x[i].w;
480 if (r > r_max)
481 r_max = r;
482 }
483 return r;
484 }
485
486 // Calculate the porosity with depth, and write to file in output directory
487 void DEM::porosity(const int z_slices)
488 {
489 // The porosity value is higher at the boundaries due
490 // to the no-flux BCs.
491
492 Float top;
493 if (walls.nw > 0)
494 top = walls.nx->w;
495 else
496 top = grid.L[2];
497
498 // Calculate depth slice thickness
499 Float h_slice = (top - grid.origo[2]) / (Float)z_slices;
500
501 // Check that the vertical slice height does not exceed the
502 // max particle diameter, since this function doesn't work
503 // if the sphere intersects more than 1 boundary
504 if (h_slice <= r_max()*2.0) {
505 std::cerr << "Error! The number of z-slices is too high."
506 << std::endl;
507 exit(1);
508 }
509
510 // Calculate slice volume
511 Float V_slice = h_slice * grid.L[0] * grid.L[1];
512
513 // Array of depth values
514 Float z_pos[z_slices];
515
516 // Array of porosity values
517 Float porosity[z_slices];
518
519 // Loop over vertical slices
520 #pragma omp parallel for if(np > 100)
521 for (int iz = 0; iz<z_slices; ++iz) {
522
523 // The void volume equals the slice volume, with the
524 // grain volumes subtracted
525 Float V_void = V_slice;
526
527 // Bottom and top position of depth slice
528 Float z_slice_low = iz * h_slice;
529 Float z_slice_high = z_slice_low + h_slice;
530
531 // Loop over particles to see whether they are inside of the slice
532 for (unsigned int i = 0; i<np; ++i) {
533
534 // Read particle values
535 Float z_sphere_centre = k.x[i].z;
536 Float radius = k.x[i].w;
537
538 // Save vertical positions of particle boundaries
539 Float z_sphere_low = z_sphere_centre - radius;
540 Float z_sphere_high = z_sphere_centre + radius;
541
542 // Sphere volume
543 Float V_sphere = 4.0/3.0 * M_PI * radius * radius * radius;
544
545 // If the sphere is inside the slice and not intersecting the
546 // boundaries, subtract the entire sphere volume
547 if (z_slice_low <= z_sphere_low && z_sphere_high <= z_slice_high) {
548 V_void -= V_sphere;
549
550 } else {
551
552 // If the sphere intersects with the lower boundary,
553 // and the centre is below the boundary
554 if (z_slice_low >= z_sphere_centre && z_slice_low <= z_sphere_high) {
555
556 // Subtract the volume of a spherical cap
557 V_void -= sphericalCap(z_sphere_high - z_slice_low, radius);
558 }
559
560 // If the sphere intersects with the lower boundary,
561 // and the centre is above the boundary
562 else if (z_slice_low < z_sphere_centre && z_slice_low > z_sphere_low) {
563
564 // Subtract the volume of the sphere,
565 // then add the volume of the spherical cap below
566 V_void -= V_sphere + sphericalCap(z_slice_low - z_sphere_low, radius);
567 }
568
569 // If the sphere intersects with the upper boundary,
570 // and the centre is above the boundary
571 if (z_slice_high <= z_sphere_centre && z_slice_high >= z_sphere_low) {
572
573 // Subtract the volume of the spherical cap below
574 V_void -= sphericalCap(z_slice_high - z_sphere_low, radius);
575 }
576 // If the sphere intersects with the upper boundary,
577 // and the centre is below the boundary
578 else if (z_slice_high > z_sphere_centre && z_slice_high < z_sphere_high) {
579
580 // Subtract the volume of the sphere,
581 // then add the volume of the spherical cap above
582 V_void -= V_sphere + sphericalCap(z_sphere_high - z_slice_high, radius);
583 }
584
585 }
586 }
587
588 // Save the mid z-point
589 z_pos[iz] = z_slice_low + 0.5*h_slice;
590
591 // Save the porosity
592 porosity[iz] = V_void / V_slice;
593
594 }
595
596 // Save results to text file
597 //writePorosities(("output/" + sid + "-porosity.txt").c_str(), z_slices, z_pos, porosity);
598
599 // Report values to stdout
600 //std::cout << "z-pos" << '\t' << "porosity" << '\n';
601 for (int i = 0; i<z_slices; ++i) {
602 std::cout << z_pos[i] << '\t' << porosity[i] << '\n';
603 }
604
605 }
606
607 // Find the min. spatial positions of the particles, and return these as a vector
608 Float3 DEM::minPos()
609 {
610 unsigned int i;
611 Float3 shared_min = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
612
613 #pragma omp parallel if(np > 100)
614 {
615 // Max. val per thread
616 Float3 min = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
617
618 #pragma omp for nowait
619 // Find min val. per thread
620 for (i = 0; i<np; ++i) {
621 min.x = std::min(min.x, k.x[i].x);
622 min.y = std::min(min.y, k.x[i].y);
623 min.z = std::min(min.z, k.x[i].z);
624 }
625
626 // Find total min, by comparing one thread with the
627 // shared result, one at a time
628 #pragma omp critical
629 {
630 shared_min.x = std::min(shared_min.x, min.x);
631 shared_min.y = std::min(shared_min.y, min.y);
632 shared_min.z = std::min(shared_min.z, min.z);
633 }
634 }
635
636 // Return final result
637 return shared_min;
638 }
639
640 // Find the max. spatial positions of the particles, and return these as a vector
641 Float3 DEM::maxPos()
642 {
643 unsigned int i;
644 Float3 shared_max = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
645
646 #pragma omp parallel if(np > 100)
647 {
648 // Max. val per thread
649 Float3 max = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
650
651 #pragma omp for nowait
652 // Find max val. per thread
653 for (i = 0; i<np; ++i) {
654 max.x = std::max(max.x, k.x[i].x);
655 max.y = std::max(max.y, k.x[i].y);
656 max.z = std::max(max.z, k.x[i].z);
657 }
658
659 // Find total max, by comparing one thread with the
660 // shared result, one at a time
661 #pragma omp critical
662 {
663 shared_max.x = std::max(shared_max.x, max.x);
664 shared_max.y = std::max(shared_max.y, max.y);
665 shared_max.z = std::max(shared_max.z, max.z);
666 }
667 }
668
669 // Return final result
670 return shared_max;
671 }
672
673
674 // Finds all overlaps between particles,
675 // returns the indexes as a 2-row vector and saves
676 // the overlap size
677 void DEM::findOverlaps(
678 std::vector< std::vector<unsigned int> > &ij,
679 std::vector< Float > &delta_n_ij)
680 {
681 unsigned int i, j;
682 Float4 x_i, x_j; // radius is in .w component of struct
683 Float3 x_ij;
684 Float x_ij_length, delta_n;
685
686 // Loop over particles, find intersections
687 for (i=0; i<np; ++i) {
688
689 for (j=0; j<np; ++j) {
690
691 // Only check once par particle pair
692 if (i < j) {
693
694 x_i = k.x[i];
695 x_j = k.x[j];
696
697 x_ij = MAKE_FLOAT3(
698 x_j.x - x_i.x,
699 x_j.y - x_i.y,
700 x_j.z - x_i.z);
701
702 x_ij_length = sqrt(
703 x_ij.x * x_ij.x +
704 x_ij.y * x_ij.y +
705 x_ij.z * x_ij.z);
706
707 // Distance between spheres
708 delta_n = x_ij_length - (x_i.w + x_j.w);
709
710 // Is there overlap?
711 // Do not include if both particles are fixed
712 if (delta_n < 0.0 && (k.vel[i].w == 0 && k.vel[j].w == 0)) {
713
714 // Store particle indexes and delta_n
715 std::vector<unsigned int> ij_pair(2);
716 ij_pair[0] = i;
717 ij_pair[1] = j;
718 ij.push_back(ij_pair);
719 delta_n_ij.push_back(delta_n);
720
721 }
722 }
723 }
724 }
725 }
726
727 // Calculate force chains and generate visualization script
728 void DEM::forcechains(const std::string format, const int threedim,
729 const double lower_cutoff,
730 const double upper_cutoff)
731 {
732 using std::cout;
733 using std::endl;
734
735 // Loop over all particles, find intersections
736 std::vector< std::vector<unsigned int> > ij;
737 std::vector< Float > delta_n_ij;
738 findOverlaps(ij, delta_n_ij);
739
740 // Find minimum position
741 Float3 x_min = minPos();
742 Float3 x_max = maxPos();
743
744 // Find largest overlap, used for scaling the line thicknesses
745 Float delta_n_min = *std::min_element(delta_n_ij.begin(), delta_n_ij.end());
746 Float f_n_max = -params.k_n * delta_n_min;
747
748 // Define limits of visualization [0;1]
749 //Float lim_low = 0.1;
750 //Float lim_low = 0.15;
751 //Float lim_high = 0.25;
752
753 if (format == "txt") {
754 // Write text header
755 cout << "x_1, [m]\t";
756 if (threedim == 1)
757 cout << "y_1, [m]\t";
758 cout << "z_1, [m]\t";
759 cout << "x_2, [m]\t";
760 if (threedim == 1)
761 cout << "y_2, [m]\t";
762 cout << "z_2, [m]\t";
763 cout << "||f_n||, [N]"
764 << "i\tj" << endl;
765
766
767 } else {
768
769 // Format sid so LaTeX won't encounter problems with the extension
770 std::string s = sid;
771 std::replace(s.begin(), s.end(), '.', '-');
772
773 // Write Gnuplot header
774 cout << "#!/usr/bin/env gnuplot\n"
775 << "# This Gnuplot script is automatically generated using\n"
776 << "# the forcechain utility in sphere. For more information,\n"
777 << "# see https://github.com/anders-dc/sphere\n"
778 << "set size ratio -1\n";
779 if (format == "png") {
780 cout << "set term pngcairo size 30 cm,20 cm\n";
781 cout << "set out '" << s << "-fc.png'\n";
782 } else if (format == "epslatex") {
783 cout << "set term epslatex size 8.6 cm, 5.6 cm\n";
784 //cout << "set out 'plots/" << s << "-fc.tex'\n";
785 cout << "set out '" << s << "-fc.tex'\n";
786 } else if (format == "epslatex-color") {
787 //cout << "set term epslatex color size 12 cm, 8.6 cm\n";
788 cout << "set term epslatex color size 8.6 cm, 5.6 cm\n";
789 cout << "set out 'plots/" << s << "-fc.tex'\n";
790 //cout << "set out '" << s << "-fc.tex'\n";
791 }
792 cout << "set xlabel '\\sffamily $x_1$, [m]'\n";
793 if (threedim == 1) {
794 cout << "set ylabel '\\sffamily $x_2$, [m]'\n"
795 << "set zlabel '\\sffamily $x_3$, [m]' offset 2\n";
796 //<< "set zlabel '\\sffamily $x_3$, [m]' offset 0\n";
797 } else
798 //cout << "set ylabel '\\sffamily $x_3$, [m]' offset 0\n";
799 cout << "set ylabel '\\sffamily $x_3$, [m]' offset 2\n";
800
801 cout << "set cblabel '\\sffamily $||\\boldsymbol{f}_n||$, [N]'\n"
802 << "set xyplane at " << x_min.z << '\n'
803 << "set pm3d\n"
804 << "set view 90.0,0.0\n"
805 //<< "set palette defined (0 'gray', 0.5 'blue', 1 'red')\n"
806 //<< "set palette defined (0 'white', 0.5 'gray', 1 'red')\n"
807
808 // MATLAB style jet colorscale
809 //<< "set palette defined ( 1 '#000fff', 2 '#0090ff', 3 '#0fffee', 4 '#90ff70', 5 '#ffee00', 6 '#ff7000', 7 '#ee0000', 8 '#7f0000')\n"
810
811 // Light gray to black
812 << "set palette defined ( 1 '#d3d3d3', 2 '#000000')\n"
813
814 // White to black (useful when using lc = 0)
815 //<< "set palette defined ( 1 '#ffffff', 2 '#000000')\n"
816
817 //<< "set cbrange [" << f_n_max*lim_low << ':' << f_n_max*lim_high << "]\n"
818 << "set cbrange [" << lower_cutoff << ':' << upper_cutoff << "]\n"
819 << endl;
820 }
821
822 // Loop over found contacts, report to stdout
823 unsigned int n, i, j;
824 Float delta_n, f_n, ratio;
825 std::string color;
826 //const double thickness_scaling = 8.0;
827 const double thickness_scaling = 16.0;
828 for (n=0; n<ij.size(); ++n) {
829
830 // Get contact particle indexes
831 i = ij[n][0];
832 j = ij[n][1];
833
834 // Overlap size
835 delta_n = delta_n_ij[n];
836
837 // Normal force on contact
838 f_n = -params.k_n * delta_n;
839
840 // TODO: Use Young's modulus if it is greater than 0
841
842 if (f_n < lower_cutoff)
843 continue; // skip the rest of this iteration
844
845 // Line weight
846 ratio = f_n/f_n_max;
847
848 if (format == "txt") {
849
850 // Text output
851 cout << k.x[i].x << '\t';
852 if (threedim == 1)
853 cout << k.x[i].y << '\t';
854 cout << k.x[i].z << '\t';
855 cout << k.x[j].x << '\t';
856 if (threedim == 1)
857 cout << k.x[j].y << '\t';
858 cout << k.x[j].z << '\t';
859 cout << f_n
860 << i << '\t' << j << endl;
861 } else {
862
863 // Gnuplot output
864 // Save contact pairs if they are above the lower limit
865 // and not fixed at their horizontal velocity
866 //if (ratio > lim_low && (k.vel[i].w + k.vel[j].w) == 0.0) {
867 if (f_n > lower_cutoff && (k.vel[i].w + k.vel[j].w) == 0.0) {
868
869 // Plot contact as arrow without tip
870 cout << "set arrow " << n+1 << " from "
871 << k.x[i].x << ',';
872 if (threedim == 1)
873 cout << k.x[i].y << ',';
874 cout << k.x[i].z;
875 cout << " to " << k.x[j].x << ',';
876 if (threedim == 1)
877 cout << k.x[j].y << ',';
878 cout << k.x[j].z;
879 cout << " nohead "
880 << "lw " << ratio * thickness_scaling
881 << " lc palette cb " << f_n
882 << endl;
883 }
884 }
885
886 }
887
888 if (format != "txt") {
889 // Write Gnuplot footer
890 if (threedim == 1)
891 cout << "splot ";
892 else
893 cout << "plot ";
894
895 cout << '[' << x_min.x << ':' << x_max.x << "] ";
896 if (threedim == 1)
897 cout << '[' << x_min.y << ':' << x_max.y << "] ";
898 cout << '[' << x_min.z << ':' << x_max.z << "] " << "NaN notitle" << endl;
899 }
900 }
901
902 // Loop over all particles, find intersections. Print particle indexes,
903 // overlap distance and normal force
904 void DEM::printContacts()
905 {
906 std::vector< std::vector<unsigned int> > ij;
907 std::vector< Float > delta_n_ij;
908 findOverlaps(ij, delta_n_ij);
909
910 for (unsigned int n=0; n<ij.size(); ++n)
911 std::cout << ij[n][0] << '\t' << ij[n][1] << '\t' << delta_n_ij[n]
912 << std::endl;
913 }
914
915 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4