timplementing forcechains2d plot in sphere.py - 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 d56b993770731d00848516d7910fd41e9f9d7356
(DIR) parent df6ae0a224ab005821c2c6e95a35430975018268
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 23 Sep 2014 15:43:56 +0200
implementing forcechains2d plot in sphere.py
Diffstat:
M python/permeability-results.py | 2 ++
M python/shear-results-forces.py | 5 +++--
M python/shear-results.py | 12 ++++++++----
M python/shear-starter.py | 6 ++++--
M python/sphere.py | 54 +++++++++++++++++++++++++++++++
5 files changed, 71 insertions(+), 8 deletions(-)
---
(DIR) diff --git a/python/permeability-results.py b/python/permeability-results.py
t@@ -60,6 +60,8 @@ for c_grad_p in cvals:
sim.readlast(verbose=False)
Re[c][i] = numpy.mean(sim.ReynoldsNumber())
+ #sim.writeVTKall()
+
# find magnitude of fluid pressure force and total interaction force
'''
fp_magn = numpy.empty(sim.np)
(DIR) diff --git a/python/shear-results-forces.py b/python/shear-results-forces.py
t@@ -15,11 +15,12 @@ from matplotlib.ticker import MaxNLocator
#steps = [5, 10, 100]
#steps = [5, 10]
-steps = sys.argv[2:]
+steps = sys.argv[3:]
nsteps_avg = 5 # no. of steps to average over
sigma0 = float(sys.argv[1])
-c_grad_p = 1.0
+#c_grad_p = 1.0
+c_grad_p = float(sys.argv[2])
c_phi = 1.0
sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + \
(DIR) diff --git a/python/shear-results.py b/python/shear-results.py
t@@ -16,8 +16,8 @@ import matplotlib.pyplot as plt
#sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
#sigma0 = 10.0e3
sigma0 = float(sys.argv[1])
-#cvals = [1.0, 0.1]
-cvals = [1.0]
+cvals = [1.0, 0.1]
+#cvals = [1.0]
c_phi = 1.0
shear_strain = [[], [], []]
t@@ -34,7 +34,8 @@ sid = 'shear-sigma0=' + sys.argv[1] + '-hw'
sim = sphere.sim(sid)
sim.readlast(verbose=False)
sim.visualize('shear')
-shear_strain[0] = sim.shear_strain
+#shear_strain[0] = sim.shear_strain
+shear_strain[0] = numpy.arange(sim.status()+1)
friction[0] = sim.tau/sim.sigma_eff
dilation[0] = sim.dilation
t@@ -55,7 +56,8 @@ for c in numpy.arange(1,len(cvals)+1):
sim.readlast(verbose=False)
sim.visualize('shear')
- shear_strain[c] = sim.shear_strain
+ #shear_strain[c] = sim.shear_strain
+ shear_strain[c] = numpy.arange(sim.status()+1)
friction[c] = sim.tau/sim.sigma_eff
dilation[c] = sim.dilation
t@@ -116,6 +118,8 @@ ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]')
ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]')
ax3.set_ylabel('Fluid pressure $p_\\text{f}$ [kPa]')
+ax1.set_xlim([200,300])
+
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax2.get_xticklabels(), visible=False)
(DIR) diff --git a/python/shear-starter.py b/python/shear-starter.py
t@@ -27,7 +27,7 @@ sim.readlast()
if fluid:
sim.sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + str(c_phi) + \
- '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc-hw'
+ '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc-hw-noshear'
else:
sim.sid = 'shear-sigma0=' + str(sigma0) + '-hw'
t@@ -39,7 +39,8 @@ sim.cleanup()
sim.adjustUpperWall()
sim.zeroKinematics()
-sim.shear(1.0/20.0)
+#sim.shear(1.0/20.0)
+sim.shear(0.0)
if fluid:
#sim.num[2] *= 2
t@@ -51,6 +52,7 @@ sim.setFluidTopFixedPressure()
sim.setDEMstepsPerCFDstep(10)
sim.setMaxIterations(2e5)
sim.initTemporal(total = 20.0, file_dt = 0.01, epsilon=0.07)
+#sim.initTemporal(total = 20.0, file_dt = 0.01, epsilon=0.05)
sim.c_phi[0] = c_phi
sim.c_grad_p[0] = c_grad_p
sim.w_devs[0] = sigma0
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -3604,6 +3604,60 @@ class sim:
return self.shearStrainRate() * numpy.mean(self.radius) \
* numpy.sqrt(self.rho[0]/self.w_devs[0])
+ def findOverlaps(self):
+ '''
+ Find all particle-particle overlaps by a n^2 contact search. The
+ particle pair indexes and the distance of the overlaps is saved in the
+ object itself as the ``.pairs`` and ``.overlaps`` members.
+ '''
+
+ # Allocate big arrays instead of dynamically allocating during the
+ # iterative process. Blank positions will be erased afterwards
+ max_avg_contacts_per_particle = 16
+ self.pairs = numpy.empty(self.np*max_avg_contacts_per_particle, 2)
+ self.overlaps = numpy.empty_like(self.pairs)
+
+ p = 0
+ for i in numpy.arange(self.np):
+ for j in numpy.arange(self.np):
+
+ if i < j:
+ x_ij = self.x[i,:] - self.x[j,:]
+ x_ij_length = numpy.sqrt(x_ij.dot(x_ij))
+
+ delta_n = x_ij_length - (self.radius[i] + self.radius[j])
+
+ if delta_n < 0.0:
+ self.pairs[p,:] = [i,j]
+ self.overlaps[p] = delta_n
+ p += 1
+ self.pairs = self.pairs[:p+1,:]
+ self.overlaps = self.overlaps[:p+1]
+
+
+ def forcechains2d(self, lc=200.0, uc=650.0, outformat='png', axes=[0,2]):
+ '''
+ Visualizes the force chains in the system from the magnitude of the
+ normal contact forces, and produces an image of them. The force chains
+ are orthogonally projected into a 2d plane specified by the axes
+ parameter.
+
+ :param lc: Lower cutoff of contact forces. Contacts below are not
+ visualized
+ :type lc: float
+ :param uc: Upper cutoff of contact forces. Contacts above are
+ visualized with this value
+ :type uc: float
+ :param outformat: Format of output image. Possible values are
+ 'interactive', 'png', 'epslatex', 'epslatex-color'
+ :type outformat: str
+ :param axes: The coordinate system axes in the projection plane (0:x,
+ 1:y, 2:z), default 0,2.
+ :type axes: tuple
+ '''
+
+
+
def forcechains(self, lc=200.0, uc=650.0, outformat='png', disp='2d'):
'''
Visualizes the force chains in the system from the magnitude of the