tporosity function seems to work - 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 2e13e8b98b7b64cc9915b87bf4785d0a45b26295
(DIR) parent 177aefc917e9e7c3111de2c77b60f63e01c3b959
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Mon, 17 Dec 2012 15:17:39 +0100
porosity function seems to work
Diffstat:
M src/file_io.cpp | 29 +++++++++++++++++++++++++++++
M src/porosity.cpp | 15 +++++++++++++--
M src/sphere.cpp | 70 ++++++++++++++++++++-----------
M src/sphere.h | 9 ++++++++-
4 files changed, 95 insertions(+), 28 deletions(-)
---
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -372,4 +372,33 @@ void DEM::writePPM(const char *target)
ofs.close();
}
+// Write write depth vs. porosity values to file
+void DEM::writePorosities(
+ const char *target,
+ const int z_slices,
+ const Float *z_pos,
+ const Float *porosity)
+{
+ // Open output file
+ std::ofstream ofs(target);
+ if (!ofs) {
+ std::cerr << "Could not create output porosity file '"
+ << target << std::endl;
+ exit(1); // Return unsuccessful exit status
+ }
+
+ if (verbose == 1)
+ std::cout << " Saving porosities: " << target << std::endl;
+
+ // Write pixel array to ppm image file
+ for (int i=0; i<z_slices; ++i)
+ ofs << z_pos[i] << '\t' << porosity[i] << '\n';
+
+ // Close file if it is still open
+ if (ofs.is_open())
+ ofs.close();
+}
+
+
+
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/porosity.cpp b/src/porosity.cpp
t@@ -29,6 +29,7 @@ int main(const int argc, const char *argv[])
// Default values
int verbose = 1;
int nfiles = 0; // number of input files
+ int slices = 10; // number of vertical slices
// Process input parameters
int i;
t@@ -42,7 +43,8 @@ int main(const int argc, const char *argv[])
<< "Usage: " << argv[0] << " [OPTION[S]]... [FILE1 ...]\nOptions:\n"
<< "-h, --help\t\tprint help\n"
<< "-V, --version\t\tprint version information and exit\n"
- << "-q, --quiet\t\tDo not display in-/output file names\n"
+ << "-q, --quiet\t\tdo not display in-/output file names\n"
+ << "-s. --slices\t\tnumber of vertical slices to find porosity within\n"
<< "The porosity values are stored in the output/ folder"
<< std::endl;
return 0; // Exit with success
t@@ -58,6 +60,15 @@ int main(const int argc, const char *argv[])
else if (argvi == "-q" || argvi == "--quiet")
verbose = 0;
+ else if (argvi == "-s" || argvi == "--slices") {
+ slices = atoi(argv[++i]);
+ if (slices < 1) {
+ std::cerr << "Error: The number of slices must be a positive, real number (was "
+ << slices << ")" << std::endl;
+ return 1;
+ }
+ }
+
// The rest of the values must be input binary files
else {
nfiles++;
t@@ -69,7 +80,7 @@ int main(const int argc, const char *argv[])
DEM dem(argvi, verbose, 0, 0, 0);
// Calculate porosity and save as file
- dem.porosity();
+ dem.porosity(slices);
}
}
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -309,18 +309,26 @@ Float sphericalCap(const Float h, const Float r)
// Calculate the porosity with depth, and write to file in output directory
void DEM::porosity(const int z_slices)
{
+ Float top;
+ if (walls.nw > 0)
+ top = walls.nx->w;
+ else
+ top = grid.L[2];
// Calculate depth slice thickness
- Float h_slice = (walls.nx->w - grid.origo[2]) / (Float)z_slices;
+ Float h_slice = (top - grid.origo[2]) / (Float)z_slices;
// Calculate slice volume
Float V_slice = h_slice * grid.L[0] * grid.L[1];
+ // Array of depth values
+ Float z_pos[z_slices];
+
// Array of porosity values
Float porosity[z_slices];
// Loop over vertical slices
-//#pragma omp parallel for if(np > 100)
+#pragma omp parallel for if(np > 100)
for (int iz = 0; iz<z_slices; ++iz) {
// The void volume equals the slice volume, with the
t@@ -328,75 +336,87 @@ void DEM::porosity(const int z_slices)
Float V_void = V_slice;
// Bottom and top position of depth slice
- Float z_low = iz * h_slice;
- Float z_high = z_low + h_slice;
+ Float z_slice_low = iz * h_slice;
+ Float z_slice_high = z_slice_low + h_slice;
// Loop over particles to see whether they are inside of the slice
for (unsigned int i = 0; i<np; ++i) {
- Float z_pos = k.x[i].z;
+ // Read particle values
+ Float z_sphere_centre = k.x[i].z;
Float radius = k.x[i].w;
+
+ // Save vertical positions of particle boundaries
+ Float z_sphere_low = z_sphere_centre - radius;
+ Float z_sphere_high = z_sphere_centre + radius;
// Sphere volume
Float V_sphere = 4.0/3.0 * M_PI * radius * radius * radius;
-
// If the sphere is inside the slice and not intersecting the
// boundaries, subtract the entire sphere volume
- if (z_low < (z_pos - radius) && (z_pos + radius) < z_high) {
+ if (z_slice_low < z_sphere_low && z_sphere_high < z_slice_high) {
V_void -= V_sphere;
- std::cout << "sphere " << i << " entirely inside volume" << std::endl;
+
} else {
// If the sphere intersects with the lower boundary,
// and the centre is below the boundary
- if (z_low > z_pos && z_low < (z_pos + radius)) {
+ if (z_slice_low > z_sphere_centre && z_slice_low < z_sphere_high) {
// Subtract the volume of a spherical cap
- V_void -= sphericalCap(z_pos + radius - z_low, radius);
-
- std::cout << "sphere " << i << " intersects lower boundary on upper hemisphere" << std::endl;
+ V_void -= sphericalCap(z_sphere_high - z_slice_low, radius);
}
// If the sphere intersects with the lower boundary,
// and the centre is above the boundary
- else if (z_low < z_pos && z_low > (z_pos - radius)) {
+ else if (z_slice_low < z_sphere_centre && z_slice_low > z_sphere_low) {
// Subtract the volume of the sphere,
// then add the volume of the spherical cap below
- V_void -= V_sphere + sphericalCap(z_low - z_pos - radius, radius);
- std::cout << "sphere " << i << " intersects lower boundary on lower hemisphere" << std::endl;
+ V_void -= V_sphere + sphericalCap(z_slice_low - z_sphere_low, radius);
}
// If the sphere intersects with the upper boundary,
// and the centre is below the boundary
- if (z_high > z_pos && z_high < (z_pos + radius)) {
+ if (z_slice_high > z_sphere_centre && z_slice_high < z_sphere_high) {
// Subtract the volume of the sphere,
// then add the volume of the spherical cap above
- V_void -= V_sphere + sphericalCap(z_pos + radius - z_high, radius);
- std::cout << "sphere " << i << " intersects upper boundary on upper hemisphere" << std::endl;
+ V_void -= V_sphere + sphericalCap(z_sphere_high - z_slice_high, radius);
}
-
+
// If the sphere intersects with the upper boundary,
// and the centre is above the boundary
- else if (z_high < z_pos && z_high > (z_pos - radius)) {
+ else if (z_slice_high < z_sphere_centre && z_slice_high > z_sphere_low) {
// Subtract the volume of the spherical cap below
- V_void -= sphericalCap(z_high - z_pos - radius, radius);
- std::cout << "sphere " << i << " intersects upper boundary on lower hemisphere" << std::endl;
+ V_void -= sphericalCap(z_slice_high - z_sphere_low, radius);
}
+
}
}
+ // Save the mid z-point
+ z_pos[iz] = z_slice_low + 0.5*h_slice;
+
// Save the porosity
porosity[iz] = V_void / V_slice;
+
+ // Report values to stdout
+ /*
+ std::cout << iz << ": V_void = " << V_void
+ << "\tV_slice = " << V_slice
+ << "\tporosity = " << V_void/V_slice
+ << '\n' << std::endl;
+ */
}
- // Display results from the top down
- for (int i=z_slices-1; i>=0; --i)
- std::cout << porosity[i] << std::endl;
+ // Save results to text file
+ writePorosities(("output/" + sid + "-porosity.txt").c_str(), z_slices, z_pos, porosity);
+ //for (int i=z_slices-1; i>=0; --i)
+ //std::cout << porosity[i] << std::endl;
}
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -122,6 +122,13 @@ class DEM {
// Find and return the max. position of any particle in each dimension
float3 maxPos(void);
+ // Write porosities found in porosity() to text file
+ void writePorosities(
+ const char *target,
+ const int z_slices,
+ const Float *z_pos,
+ const Float *porosity);
+
// Values and functions accessible from the outside
public:
t@@ -153,7 +160,7 @@ class DEM {
void startTime(void);
// Render particles using raytracing
- void render(//const char *target,
+ void render(
const int method = 1,
const float maxval = 1.0e3,
const float focalLength = 1.0,