tImplemented simple force chain visualization with asymptote - 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 e0ae4137b4d832533ca332402af6cec28d43b8fa
(DIR) parent 220baf91b8b2819a8b4708244e4ef7a2ce762105
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Mon, 14 Jan 2013 12:47:51 +0100
Implemented simple force chain visualization with asymptote
Diffstat:
M src/CMakeLists.txt | 6 ++++++
A src/forcechains.cpp | 89 +++++++++++++++++++++++++++++++
M src/sphere.cpp | 93 +++++++++++++++++++++++++++++++
M src/sphere.h | 11 +++++++++++
4 files changed, 199 insertions(+), 0 deletions(-)
---
(DIR) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
t@@ -16,4 +16,10 @@ SET(CUDA_NVCC_FLAGS "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,c
# Rule to build executable program
CUDA_ADD_EXECUTABLE(../sphere main.cpp file_io.cpp sphere.cpp device.cu utility.cu)
CUDA_ADD_EXECUTABLE(../porosity porosity.cpp file_io.cpp sphere.cpp device.cu utility.cu)
+CUDA_ADD_EXECUTABLE(../forcechains forcechains.cpp file_io.cpp sphere.cpp device.cu utility.cu)
+#ADD_EXECUTABLE(unittests boost-unit-tests.cpp sphere.cpp)
+#TARGET_LINK_LIBRARIES(unittests
+# ${Boost_FILESYSTEM_LIBRARY}
+# ${Boost_SYSTEM_LIBRARY}
+# ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
(DIR) diff --git a/src/forcechains.cpp b/src/forcechains.cpp
t@@ -0,0 +1,89 @@
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+/* SPHERE source code by Anders Damsgaard Christensen, 2010-12, */
+/* a 3D Discrete Element Method algorithm with CUDA GPU acceleration. */
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+
+// Licence: GNU Public License (GPL) v. 3. See license.txt.
+// See doc/sphere-doc.pdf for full documentation.
+// Compile with GNU make by typing 'make' in the src/ directory.
+// SPHERE is called from the command line with './sphere_<architecture> projectname'
+
+
+// Including library files
+#include <iostream>
+#include <string>
+#include <cstdlib>
+
+// Including user files
+#include "constants.h"
+#include "datatypes.h"
+#include "sphere.h"
+
+//////////////////
+// MAIN ROUTINE //
+//////////////////
+// The main loop returns the value 0 to the shell, if the program terminated
+// successfully, and 1 if an error occured which caused the program to crash.
+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
+
+ // Process input parameters
+ int i;
+ for (i=1; i<argc; ++i) { // skip argv[0]
+
+ std::string argvi = std::string(argv[i]);
+
+ // 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::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::endl;
+ return 0;
+ }
+
+ else if (argvi == "-v" || argvi == "--verbose")
+ verbose = 1;
+
+ // The rest of the values must be input binary files
+ else {
+ nfiles++;
+
+ if (verbose == 1)
+ std::cout << argv[0] << ": processing input file: " << argvi << std::endl;
+
+ // Create DEM class, read data from input binary, check values
+ DEM dem(argvi, verbose, 0, 0, 0);
+
+ // Calculate porosity and save as file
+ dem.forcechains();
+
+ }
+ }
+
+ // Check whether there are input files specified
+ if (!argv[0] || argc == 1 || nfiles == 0) {
+ std::cerr << argv[0] << ": missing input binary file\n"
+ << "See `" << argv[0] << " --help` for more information"
+ << std::endl;
+ return 1; // Return unsuccessful exit status
+ }
+
+ return 0; // Return successfull exit status
+}
+// END OF FILE
+// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -3,6 +3,7 @@
#include <cstdio>
#include <cstdlib>
#include <cmath>
+#include <vector>
#include "typedefs.h"
#include "datatypes.h"
t@@ -443,4 +444,96 @@ void DEM::porosity(const int z_slices)
}
+// Finds all overlaps between particles,
+// returns the indexes as a 2-row vector and saves
+// the overlap size
+void DEM::findOverlaps(
+ std::vector< std::vector<unsigned int> > &ij,
+ std::vector< Float > &delta_n_ij)
+{
+ unsigned int i, j;
+ Float4 x_i, x_j; // radius is in .w component of struct
+ Float3 x_ij;
+ Float x_ij_length, delta_n;
+
+ // Loop over particles, find intersections
+ for (i=0; i<np; ++i) {
+
+ for (j=0; j<np; ++j) {
+
+ // Only check once par particle pair
+ if (i < j) {
+
+ x_i = k.x[i];
+ x_j = k.x[j];
+
+ x_ij = MAKE_FLOAT3(
+ x_j.x - x_i.x,
+ x_j.y - x_i.y,
+ x_j.z - x_i.z);
+
+ x_ij_length = sqrt(
+ x_ij.x * x_ij.x +
+ x_ij.y * x_ij.y +
+ x_ij.z * x_ij.z);
+
+ // Distance between spheres
+ delta_n = x_ij_length - (x_i.w + x_j.w);
+
+ // Is there overlap?
+ if (delta_n < 0.0) {
+
+ // Store particle indexes and delta_n
+ ij.push_back(std::vector<unsigned int> (i, j));
+ delta_n_ij.push_back(delta_n);
+
+ }
+ }
+ }
+ }
+}
+
+// Calculate force chains and generate visualization script
+void DEM::forcechains()
+{
+ using std::cout;
+ using std::endl;
+
+ // Loop over all particles, find intersections
+ std::vector< std::vector<unsigned int> > ij;
+ std::vector< Float > delta_n_ij;
+ findOverlaps(ij, delta_n_ij);
+
+ // Write Asymptote header
+ cout << "import three; \n size(600);" << endl;
+
+ // Loop over found contacts, report to stdout
+ unsigned int n, i, j;
+ std::vector<unsigned int> ijs(1);
+ Float delta_n;
+ for (n=0; n<ij.size(); ++n) {
+
+ ijs = ij[n];
+ i = ijs[0];
+ j = ijs[1];
+
+ delta_n = delta_n_ij[n];
+
+ /*cout << "Contact n = " << n
+ << ": (i,j) = ("
+ << i << ","
+ << j << "), delta_n = " << delta_n
+ << endl;*/
+
+ 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 << "); \n draw(g);" << endl;
+ }
+
+}
+
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -2,6 +2,8 @@
#ifndef SPHERE_H_
#define SPHERE_H_
+#include <vector>
+
#include "datatypes.h"
// DEM class
t@@ -177,6 +179,15 @@ class DEM {
// Calculate porosity with depth and save as text file
void porosity(const int z_slices = 10);
+ // 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);
+
};
#endif