tAdded bondsRose, visualizing strike and dip of bonds - 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 cb39751b5def0109642e7bb7769cb07863792384
(DIR) parent e86480db55360ef61f8cb441c0160817fec8689a
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Tue, 12 Mar 2013 06:09:19 +0100
Added bondsRose, visualizing strike and dip of bonds
Diffstat:
M python/sphere.py | 53 +++++++++++++++++++++++++++++--
1 file changed, 51 insertions(+), 2 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -1234,7 +1234,7 @@ class Spherebin:
subprocess.call("cd .. && ./forcechains " + nd + "-f " + outformat + " -lc " + str(lc) + " -uc " + str(uc) + " input/" + self.sid + ".bin > python/tmp.gp", shell=True)
subprocess.call("gnuplot tmp.gp && rm tmp.bin && rm tmp.gp", shell=True)
- def forcechainsRose(self):
+ def forcechainsRose(self, lower_limit=0.25):
''' Visualize strike- and dip angles of the strongest force chains in a
rose plot.
'''
t@@ -1249,7 +1249,7 @@ class Spherebin:
f_n_max = numpy.amax(data[:,6])
# specify the lower limit of force chains to do statistics on
- f_n_lim = 0.25 * f_n_max * 0.6
+ f_n_lim = lower_limit * f_n_max * 0.6
# find the indexes of these contacts
I = numpy.nonzero(data[:,6] > f_n_lim)
t@@ -1298,6 +1298,55 @@ class Spherebin:
subprocess.call("rm fc-tmp.txt", shell=True)
+ def bondsRose(self):
+ ''' Visualize strike- and dip angles of the bond pairs in a rose plot.
+ '''
+ # loop through these contacts and find the strike and dip of the contacts
+ strikelist = [] # strike direction of the normal vector, [0:360[
+ diplist = [] # dip of the normal vector, [0:90]
+ for n in numpy.arange(self.nb0):
+
+ i = self.bonds[n,0]
+ j = self.bonds[n,1]
+
+ x1 = self.x[i,0]
+ y1 = self.x[i,1]
+ z1 = self.x[i,2]
+ x2 = self.x[j,0]
+ y2 = self.x[j,1]
+ z2 = self.x[j,2]
+
+ if (z1 < z2):
+ xlower = x1; ylower = y1; zlower = z1
+ xupper = x2; yupper = y2; zupper = z2
+ else :
+ xlower = x2; ylower = y2; zlower = z2
+ xupper = x1; yupper = y1; zupper = z1
+
+ # Vector pointing downwards
+ dx = xlower - xupper
+ dy = ylower - yupper
+ dz = zlower - zupper
+ dhoriz = numpy.sqrt(dx**2 + dy**2)
+
+ # Find dip angle
+ diplist.append(math.degrees(math.atan((zupper - zlower)/dhoriz)))
+
+ # Find strike angle
+ if (ylower >= yupper): # in first two quadrants
+ strikelist.append(math.acos(dx/dhoriz))
+ else :
+ strikelist.append(2.0*numpy.pi - math.acos(dx/dhoriz))
+
+
+ plt.figure(figsize=[4,4])
+ ax = plt.subplot(111, polar=True, axisbg="w")
+ ax.scatter(strikelist, diplist, c='k', marker='+')
+ ax.set_rmax(90)
+ ax.set_rticks([])
+ plt.savefig("bonds-" + self.sid + "-rose.pdf", transparent=True)
+
+
def thinsection_x1x3(self,
x2 = 'center',