tadd particle indexes to forcechain output - 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 4032f0705e116464014078e8e0eb5779e4de7042
(DIR) parent 66d59542ce89a0a6c7ed15ab7a1918e017a16a90
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 21 May 2015 16:25:53 +0200
add particle indexes to forcechain output
Diffstat:
M python/halfshear-darcy-creep-dynam… | 40 +++++++++++++++++++++++++++++--
M python/sphere.py | 2 +-
M src/sphere.cpp | 49 +++++++++++++++----------------
3 files changed, 63 insertions(+), 28 deletions(-)
---
(DIR) diff --git a/python/halfshear-darcy-creep-dynamics.py b/python/halfshear-darcy-creep-dynamics.py
t@@ -30,10 +30,12 @@ matplotlib.rc('grid', linestyle=':', linewidth=0.2)
outformat='pdf'
scatter=False
-plotContacts=True
+plotContacts=False
#plotContacts=False
plotForceChains=True
#plotForceChains=False
+plotHistograms=False
+plotCombinedHistogram=True
#sids = ['halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4000.0-f=0.2']
sids = ['halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4375.0-f=0.2']
t@@ -232,7 +234,7 @@ for sid in sids:
axsc1.set_rticks([])
#axsc1.grid(False)
- else:
+ if plotHistograms:
axhistload1 = fig.add_axes([Lx[sc], Ly[sc], dx, dy*.5])
axhistload1.hist(datalists[sc][:,6], alpha=0.75, facecolor='gray', log=True)
#plt.yscale('log', nonposy='clip')
t@@ -307,3 +309,37 @@ for sid in sids:
plt.close()
shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
print(filename)
+
+if plotCombinedHistogram:
+ fig = plt.figure(figsize=[3.5, 3.5])
+
+ ax1 = plt.subplot(1, 1, 1)
+
+ f_min = 0.0
+ f_max = 0.0
+ for sc in range(len(Lx)):
+ if f_max < numpy.max(datalists[sc][:,6]):
+ f_max = numpy.max(datalists[sc][:,6])
+ f_max = numpy.ceil(f_max/10.)*10.
+ print(f_max)
+
+ for sc in range(len(Lx)):
+ #axhistload1.hist(datalists[sc][:,6], alpha=0.75, facecolor='gray', log=True)
+ hist, bin_edges = numpy.histogram(datalists[sc][:,6],
+ #numpy.arange(8),
+ range=(f_min, f_max))
+
+ ax1.semilogy((bin_edges[1:] - bin_edges[:-1])/2 + bin_edges[:-1],
+ hist, 'o-',
+ label='$N$ = {:3.1f} kPa'.format(Ns[sc]/1000.))
+
+ #plt.yscale('log', nonposy='clip')
+ ax1.set_xlabel('Contact load [N]')
+ ax1.set_ylabel('Number of contacts')
+ ax1.legend(loc='upper right', fancybox=True, framealpha=1.0)
+
+ filename = sid + '-creep-dynamics-hist.' + outformat
+ plt.savefig(filename)
+ plt.close()
+ shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
+ print(filename)
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -5393,7 +5393,7 @@ class sim:
plt.axvline(90. - theta_sigma1, color='k', linestyle='dashed',
linewidth=1)
plt.xlim([0, 90.])
- plt.ylim([0, self.np[0]/100])
+ plt.ylim([0, self.np[0]/10])
#plt.xlabel('$\\boldsymbol{f}_\text{n}$ [N]')
plt.xlabel('Contact angle [deg]')
plt.ylabel('Count $N$')
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -13,7 +13,7 @@
// Constructor: Reads an input binary, and optionally checks
// and reports the values
-DEM::DEM(const std::string inputbin,
+DEM::DEM(const std::string inputbin,
const int verbosity,
const int checkVals,
const int dry,
t@@ -26,7 +26,7 @@ DEM::DEM(const std::string inputbin,
using std::cout;
using std::cerr;
- // Extract sid from input binary filename
+ // Extract sid from input binary filename
size_t dotpos = inputbin.rfind('.');
size_t slashpos = inputbin.rfind('/');
if (slashpos - dotpos < 1) {
t@@ -359,21 +359,21 @@ void DEM::reportValues()
if (verbose == 1)
cout << "double";
} else {
- cerr << "Error! Chosen precision not available. Check datatypes.h"
+ cerr << "Error! Chosen precision not available. Check datatypes.h"
<< endl;
exit(1);
}
cout << " precision\n";
- cout << " - Timestep length: time.dt = "
+ cout << " - Timestep length: time.dt = "
<< time.dt << " s\n"
- << " - Start at time: time.current = "
+ << " - Start at time: time.current = "
<< time.current << " s\n"
- << " - Total sim. time: time.total = "
+ << " - Total sim. time: time.total = "
<< time.total << " s\n"
- << " - File output interval: time.file_dt = "
+ << " - File output interval: time.file_dt = "
<< time.file_dt << " s\n"
- << " - Start at step count: time.step_count = "
+ << " - Start at step count: time.step_count = "
<< time.step_count << endl;
if (params.contactmodel == 1)
t@@ -420,8 +420,8 @@ void DEM::reportValues()
cout << grid.num[0];
else if (nd == 2)
cout << grid.num[0] << " * " << grid.num[1];
- else
- cout << grid.num[0] << " * "
+ else
+ cout << grid.num[0] << " * "
<< grid.num[1] << " * "
<< grid.num[2];
cout << " cells\n";
t@@ -473,9 +473,9 @@ Float DEM::r_max()
// Calculate the porosity with depth, and write to file in output directory
void DEM::porosity(const int z_slices)
{
- // The porosity value is higher at the boundaries due
+ // The porosity value is higher at the boundaries due
// to the no-flux BCs.
-
+
Float top;
if (walls.nw > 0)
top = walls.nx->w;
t@@ -521,7 +521,7 @@ void DEM::porosity(const int z_slices)
// 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;
t@@ -548,7 +548,7 @@ void DEM::porosity(const int z_slices)
// and the centre is above the boundary
else if (z_slice_low < z_sphere_centre && z_slice_low > z_sphere_low) {
- // Subtract the volume of the sphere,
+ // Subtract the volume of the sphere,
// then add the volume of the spherical cap below
V_void -= V_sphere + sphericalCap(z_slice_low - z_sphere_low, radius);
}
t@@ -564,15 +564,12 @@ void DEM::porosity(const int z_slices)
// and the centre is below the boundary
else if (z_slice_high > z_sphere_centre && z_slice_high < z_sphere_high) {
- // Subtract the volume of the sphere,
+ // Subtract the volume of the sphere,
// then add the volume of the spherical cap above
V_void -= V_sphere + sphericalCap(z_sphere_high - z_slice_high, radius);
}
-
-
}
-
}
// Save the mid z-point
t@@ -589,7 +586,7 @@ void DEM::porosity(const int z_slices)
// Report values to stdout
//std::cout << "z-pos" << '\t' << "porosity" << '\n';
for (int i = 0; i<z_slices; ++i) {
- std::cout << z_pos[i] << '\t' << porosity[i] << '\n';
+ std::cout << z_pos[i] << '\t' << porosity[i] << '\n';
}
}
t@@ -661,7 +658,7 @@ Float3 DEM::maxPos()
}
-// Finds all overlaps between particles,
+// Finds all overlaps between particles,
// returns the indexes as a 2-row vector and saves
// the overlap size
void DEM::findOverlaps(
t@@ -715,7 +712,7 @@ void DEM::findOverlaps(
}
// Calculate force chains and generate visualization script
-void DEM::forcechains(const std::string format, const int threedim,
+void DEM::forcechains(const std::string format, const int threedim,
const double lower_cutoff,
const double upper_cutoff)
{
t@@ -750,7 +747,8 @@ void DEM::forcechains(const std::string format, const int threedim,
if (threedim == 1)
cout << "y_2, [m]\t";
cout << "z_2, [m]\t";
- cout << "||f_n||, [N]" << endl;
+ cout << "||f_n||, [N]"
+ << "i\tj" << endl;
} else {
t@@ -760,7 +758,7 @@ void DEM::forcechains(const std::string format, const int threedim,
std::replace(s.begin(), s.end(), '.', '-');
// Write Gnuplot header
- cout << "#!/usr/bin/env gnuplot\n"
+ 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"
t@@ -843,7 +841,8 @@ void DEM::forcechains(const std::string format, const int threedim,
if (threedim == 1)
cout << k.x[j].y << '\t';
cout << k.x[j].z << '\t';
- cout << f_n << endl;
+ cout << f_n
+ << i << '\t' << j << endl;
} else {
// Gnuplot output
t@@ -864,7 +863,7 @@ void DEM::forcechains(const std::string format, const int threedim,
cout << k.x[j].z;
cout << " nohead "
<< "lw " << ratio * thickness_scaling
- << " lc palette cb " << f_n
+ << " lc palette cb " << f_n
<< endl;
}
}