tForce chain visualization implemented - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit af8c54c428fc06f48214ddbd8df52189fefd967f
(DIR) parent 21d1bf93a237295c7fa054f9cf3bf3e073150796
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Tue, 15 Jan 2013 11:48:24 +0100
Force chain visualization implemented
Diffstat:
M src/forcechains.cpp | 27 ++++++++++++++++++---------
M src/raytracer.cuh | 37 +++----------------------------
M src/sphere.cpp | 236 +++++++++++++++++++++++--------
M src/sphere.h | 15 ++++++++++-----
4 files changed, 209 insertions(+), 106 deletions(-)
---
(DIR) diff --git a/src/forcechains.cpp b/src/forcechains.cpp
t@@ -29,7 +29,8 @@ int main(const int argc, const char *argv[])
// Default values
int verbose = 0;
int nfiles = 0; // number of input files
- int slices = 10; // number of vertical slices
+ int threedim = 1; // 0 if 2d, 1 if 3d
+ std::string format = "interactive"; // gnuplot terminal type
// Process input parameters
int i;
t@@ -39,19 +40,21 @@ int main(const int argc, const char *argv[])
// Display help if requested
if (argvi == "-h" || argvi == "--help") {
- std::cout << argv[0] << ": sphere porosity calculator\n"
- << "Usage: " << argv[0] << " [OPTION[S]]... [FILE1 ...]\nOptions:\n"
- << "-h, --help\t\tprint help\n"
- << "-V, --version\t\tprint version information and exit\n"
- << "-v, --verbose\t\tdisplay in-/output file names\n"
- << "The force chain asymptote script(s) are stored in the output/ folder"
+ std::cout << argv[0] << ": sphere force chain visualizer\n"
+ << "Usage: " << argv[0] << " [OPTION[S]]... [FILE1 ...] > outputfile\nOptions:\n"
+ << "-h, --help\t\tPrint help\n"
+ << "-V, --version\t\tPrint version information and exit\n"
+ << "-v, --verbose\t\tDisplay in-/output file names\n"
+ << "-f, --format\t\tOutput format to stdout, interactive default. Possible values:\n"
+ << "\t\t\tinteractive, png, epslatex, epslatex-color"
+ << "-2d\t\t\twrite output as 2d coordinates (3d default)"
<< std::endl;
return 0; // Exit with success
}
// Display version with fancy ASCII art
else if (argvi == "-V" || argvi == "--version") {
- std::cout << "Force chai calculator, sphere version " << VERS
+ std::cout << "Force chain calculator, sphere version " << VERS
<< std::endl;
return 0;
}
t@@ -59,6 +62,12 @@ int main(const int argc, const char *argv[])
else if (argvi == "-v" || argvi == "--verbose")
verbose = 1;
+ else if (argvi == "-f" || argvi == "--format")
+ format = argv[++i];
+
+ else if (argvi == "-2d")
+ threedim = 0;
+
// The rest of the values must be input binary files
else {
nfiles++;
t@@ -70,7 +79,7 @@ int main(const int argc, const char *argv[])
DEM dem(argvi, verbose, 0, 0, 0);
// Calculate porosity and save as file
- dem.forcechains();
+ dem.forcechains(format, threedim);
}
}
(DIR) diff --git a/src/raytracer.cuh b/src/raytracer.cuh
t@@ -361,39 +361,6 @@ __host__ void DEM::rt_transferFromGlobalDeviceMemory(void)
checkForCudaErrors("During rt_transferFromGlobalDeviceMemory()");
}
-// Find the max. spatial positions of the particles, and return these as a vector
-__host__ float3 DEM::maxPos()
-{
- int i;
- float3 shared_max = make_float3(0.0f, 0.0f, 0.0f);
-
-#pragma omp parallel if(np > 100)
- {
- // Max. val per thread
- float3 max = make_float3(0.0f, 0.0f, 0.0f);
-
-#pragma omp for nowait
- // Find max val. per thread
- for (i = 0; i<np; ++i) {
- max.x = std::max(max.x, (float)k.x[i].x);
- max.y = std::max(max.y, (float)k.x[i].y);
- max.z = std::max(max.z, (float)k.x[i].z);
- }
-
- // Find total max, by comparing one thread with the
- // shared result, one at a time
-#pragma omp critical
- {
- shared_max.x = std::max(shared_max.x, max.x);
- shared_max.y = std::max(shared_max.y, max.y);
- shared_max.z = std::max(shared_max.z, max.z);
- }
- }
-
- // Return final result
- return shared_max;
-}
-
// Wrapper for the rt kernel
__host__ void DEM::render(
const int method,
t@@ -438,7 +405,9 @@ __host__ void DEM::render(
// Initialize camera values and transfer to constant memory
float imgw = grid.L[0]*1.35f; // Screen width in world coordinates
- float3 lookat = maxPos() / 2.0f; // Look at the centre of the mean positions
+ Float3 maxpos = maxPos();
+ // Look at the centre of the mean positions
+ float3 lookat = make_float3(maxpos.x, maxpos.y, maxpos.z) / 2.0f;
float3 eye = make_float3(
grid.L[0] * 2.3f,
grid.L[1] * -5.0f,
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -4,6 +4,7 @@
#include <cstdlib>
#include <cmath>
#include <vector>
+#include <algorithm>
#include "typedefs.h"
#include "datatypes.h"
t@@ -444,6 +445,73 @@ void DEM::porosity(const int z_slices)
}
+// Find the min. spatial positions of the particles, and return these as a vector
+Float3 DEM::minPos()
+{
+ unsigned int i;
+ Float3 shared_min = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
+
+#pragma omp parallel if(np > 100)
+ {
+ // Max. val per thread
+ Float3 min = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
+
+#pragma omp for nowait
+ // Find min val. per thread
+ for (i = 0; i<np; ++i) {
+ min.x = std::min(min.x, k.x[i].x);
+ min.y = std::min(min.y, k.x[i].y);
+ min.z = std::min(min.z, k.x[i].z);
+ }
+
+ // Find total min, by comparing one thread with the
+ // shared result, one at a time
+#pragma omp critical
+ {
+ shared_min.x = std::min(shared_min.x, min.x);
+ shared_min.y = std::min(shared_min.y, min.y);
+ shared_min.z = std::min(shared_min.z, min.z);
+ }
+ }
+
+ // Return final result
+ return shared_min;
+}
+
+// Find the max. spatial positions of the particles, and return these as a vector
+Float3 DEM::maxPos()
+{
+ unsigned int i;
+ Float3 shared_max = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
+
+#pragma omp parallel if(np > 100)
+ {
+ // Max. val per thread
+ Float3 max = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
+
+#pragma omp for nowait
+ // Find max val. per thread
+ for (i = 0; i<np; ++i) {
+ max.x = std::max(max.x, k.x[i].x);
+ max.y = std::max(max.y, k.x[i].y);
+ max.z = std::max(max.z, k.x[i].z);
+ }
+
+ // Find total max, by comparing one thread with the
+ // shared result, one at a time
+#pragma omp critical
+ {
+ shared_max.x = std::max(shared_max.x, max.x);
+ shared_max.y = std::max(shared_max.y, max.y);
+ shared_max.z = std::max(shared_max.z, max.z);
+ }
+ }
+
+ // Return final result
+ return shared_max;
+}
+
+
// Finds all overlaps between particles,
// returns the indexes as a 2-row vector and saves
// the overlap size
t@@ -497,7 +565,7 @@ void DEM::findOverlaps(
}
// Calculate force chains and generate visualization script
-void DEM::forcechains()
+void DEM::forcechains(const std::string format, const int threedim)
{
using std::cout;
using std::endl;
t@@ -507,76 +575,128 @@ void DEM::forcechains()
std::vector< Float > delta_n_ij;
findOverlaps(ij, delta_n_ij);
- // Write Asymptote header
- //cout << "import three;\nsize(600);" << endl;
-
- // Write Gnuplot header
- //cout << "#!/usr/bin/env gnuplot" << endl;
-
- // Write Matlab header
- //cout << "plot(";
+ // Find minimum position
+ Float3 x_min = minPos();
+ Float3 x_max = maxPos();
+
+ // Find largest overlap, used for scaling the line thicknesses
+ Float delta_n_min = *std::min_element(delta_n_ij.begin(), delta_n_ij.end());
+ Float f_n_max = -params.k_n * delta_n_min;
+
+ // Define limits of visualization [0;1]
+ //Float lim_low = 0.1;
+ Float lim_low = 0.05;
+ Float lim_high = 0.25;
+
+ if (format == "txt") {
+ // Write text header
+ cout << "x_1, [m]\t";
+ if (threedim == 1)
+ cout << "y_1, [m]\t";
+ cout << "z_1, [m]\t";
+ cout << "x_2, [m]\t";
+ if (threedim == 1)
+ cout << "y_2, [m]\t";
+ cout << "z_2, [m]\t";
+ cout << "||f_n||, [N]" << endl;
+
+
+ } else {
+ // Write Gnuplot header
+ cout << "#!/usr/bin/env gnuplot\n"
+ << "# This Gnuplot script is automatically generated using\n"
+ << "# the forcechain utility in sphere. For more information,\n"
+ << "# see https://github.com/anders-dc/sphere\n"
+ << "set size ratio -1\n";
+ if (format == "png")
+ cout << "set term pngcairo size 50 cm,40 cm\n";
+ else if (format == "epslatex")
+ cout << "set term epslatex size 8.6 cm, 8.6 cm\n";
+ else if (format == "epslatex-color")
+ cout << "set term epslatex color size 8.6 cm, 8.6 cm\n";
+ cout << "set xlabel '$x^1$, [m]'\n"
+ << "set ylabel '$x^2$, [m]'\n"
+ << "set zlabel '$x^3$, [m]'\n"
+ << "set cblabel '$||f_n||$, [Pa]'\n"
+ << "set xyplane at " << x_min.z << '\n'
+ << "set pm3d\n"
+ << "set view 90.0,0.0\n"
+ //<< "set palette defined (0 'gray', 0.5 'blue', 1 'red')\n"
+ //<< "set palette defined (0 'white', 0.5 'gray', 1 'red')\n"
+ << "set palette defined ( 1 '#000fff', 2 '#0090ff', 3 '#0fffee', 4 '#90ff70', 5 '#ffee00', 6 '#ff7000', 7 '#ee0000', 8 '#7f0000')\n"
+ << "set cbrange [" << f_n_max*lim_low << ':' << f_n_max*lim_high << "]\n"
+ << endl;
+ }
// Loop over found contacts, report to stdout
unsigned int n, i, j;
- Float delta_n;
+ Float delta_n, f_n, ratio;
+ std::string color;
for (n=0; n<ij.size(); ++n) {
// Get contact particle indexes
i = ij[n][0];
j = ij[n][1];
- //cout << "(i,j) = " << i << ", " << j << endl;
+ // Overlap size
delta_n = delta_n_ij[n];
- /*cout << "Contact n = " << n
- << ": (i,j) = ("
- << i << ","
- << j << "), delta_n = " << delta_n
- << endl;*/
-
- // Asymptote output
- /*cout << "path3 g=("
- << k.x[i].x << ','
- << k.x[i].y << ','
- << k.x[i].z << ")..("
- << k.x[j].x << ','
- << k.x[j].y << ','
- << k.x[j].z << ");\ndraw(g);" << endl;*/
-
- // Gnuplot output
- /*cout << "set arrow " << n+1 << " from "
- << k.x[i].x << ','
- << k.x[i].y << ','
- << k.x[i].z << " to "
- << k.x[j].x << ','
- << k.x[j].y << ','
- << k.x[j].z << " nohead "
- << endl;*/
-
- // Matlab output
- /*cout << "plot::Line3d("
- << '[' << k.x[i].x << ','
- << k.x[i].y << ','
- << k.x[i].z << "], "
- << '[' << k.x[j].x << ','
- << k.x[j].y << ','
- << k.x[j].z << "]),";*/
-
- cout << k.x[i].x << '\t'
- << k.x[i].y << '\t'
- << k.x[i].z << '\t'
- << k.x[j].x << '\t'
- << k.x[j].y << '\t'
- << k.x[j].z << '\t'
- << -delta_n*params.k_n << endl;
+ // Normal force on contact
+ f_n = -params.k_n * delta_n;
+
+ // Line weight
+ ratio = f_n/f_n_max;
+
+ if (format == "txt") {
+
+ // Text output
+ cout << k.x[i].x << '\t';
+ if (threedim == 1)
+ cout << k.x[i].y << '\t';
+ cout << k.x[i].z << '\t';
+ cout << k.x[j].x << '\t';
+ if (threedim == 1)
+ cout << k.x[j].y << '\t';
+ cout << k.x[j].z << '\t';
+ cout << -delta_n*params.k_n << endl;
+ } else {
+
+ // Gnuplot output
+ // Save contact pairs if they are above the lower limit
+ // and not fixed at their horizontal velocity
+ if (ratio > lim_low && (k.vel[i].w + k.vel[j].w) == 0.0) {
+
+ // Plot contact as arrow without tip
+ cout << "set arrow " << n+1 << " from "
+ << k.x[i].x << ',';
+ if (threedim == 1)
+ cout << k.x[i].y << ',';
+ cout << k.x[i].z;
+ cout << " to " << k.x[j].x << ',';
+ if (threedim == 1)
+ cout << k.x[j].y, ',';
+ cout << k.x[j].z;
+ cout << " nohead "
+ << "lw " << ratio * 12.0
+ << " lc palette cb " << f_n
+ << endl;
+ }
+ }
+
}
- // Write Gnuplot footer
- /*cout << "splot "
- << '[' << grid.origo[0] << ':' << grid.L[0] << "] "
- << '[' << grid.origo[1] << ':' << grid.L[1] << "] "
- << '[' << grid.origo[2] << ':' << grid.L[2] << "] "
- << "NaN notitle" << endl;*/
+ if (format != "txt") {
+ // Write Gnuplot footer
+ if (threedim == 1)
+ cout << "splot ";
+ else
+ cout << "plot ";
+
+ cout << '[' << x_min.x << ':' << x_max.x << "] ";
+ if (threedim == 1)
+ cout << '[' << x_min.y << ':' << x_max.y << "] ";
+ cout << '[' << x_min.z << ':' << x_max.z << "] " << "NaN notitle" << endl;
+ }
}
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -121,9 +121,6 @@ class DEM {
void transferFromGlobalDeviceMemory(void);
void rt_transferFromGlobalDeviceMemory(void);
- // Find and return the max. position of any particle in each dimension
- float3 maxPos(void);
-
// Find and return the max. radius
Float r_max(void);
t@@ -179,14 +176,22 @@ class DEM {
// Calculate porosity with depth and save as text file
void porosity(const int z_slices = 10);
+ // find and return the min. position of any particle in each dimension
+ Float3 minPos(void);
+
+ // find and return the max. position of any particle in each dimension
+ Float3 maxPos(void);
+
// Find particle-particle intersections, saves the indexes
// and the overlap sizes
void findOverlaps(
std::vector< std::vector<unsigned int> > &ij,
std::vector< Float > &delta_n_ij);
- // Calculate force chains and save as asymptote script
- void forcechains(void);
+ // Calculate force chains and save as Gnuplot script
+ void forcechains(
+ const std::string format = "interactive",
+ const int threedim = 1);
};