tshearVel as alias of shearVelocity - 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 3a0dffe3b2b4ac5ad6ad10473b1fb64f6741769d
 (DIR) parent 452b1fc5599e039d19dc96986a4c97cc79016887
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Mon,  2 Mar 2015 13:15:16 +0100
       
       shearVel as alias of shearVelocity
       
       Diffstat:
         M python/halfshear-darcy-combined.py  |      34 ++++++++++++++++++++++++++-----
         M python/halfshear-darcy-rate-depend… |      19 ++++++++++++++-----
         M python/sphere.py                    |      24 ++++++------------------
       
       3 files changed, 49 insertions(+), 28 deletions(-)
       ---
 (DIR) diff --git a/python/halfshear-darcy-combined.py b/python/halfshear-darcy-combined.py
       t@@ -12,6 +12,7 @@ import numpy
        import sphere
        import matplotlib.pyplot as plt
        import matplotlib.patches
       +import matplotlib.colors
        
        sid = 'halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2'
        outformat = 'pdf'
       t@@ -56,6 +57,7 @@ phi_bar   = numpy.empty_like(t)
        
        # mean horizontal porosity plot
        poros     = numpy.empty((sim.num[2], nsteps))
       +#xvel      = numpy.zeros((sim.num[2], nsteps))
        zpos_c    = numpy.empty(sim.num[2])
        dz = sim.L[2]/sim.num[2]
        for i in numpy.arange(sim.num[2]):
       t@@ -87,6 +89,21 @@ for i in numpy.arange(nsteps):
        
            poros[:,i]   = numpy.average(numpy.average(sim.phi, axis=0),axis=0)
        
       +    # calculate mean values of xvel
       +    '''
       +    dz = sim.L[2]/(sim.num[2])
       +    for iz in numpy.arange(sim.num[2]):
       +        z_bot = iz*dz
       +        z_top = (iz+1)*dz
       +        I = numpy.nonzero((sim.x[:,2] >= z_bot) & (sim.x[:,2] < z_top))
       +        if I[0].size > 0:
       +            #print I
       +            #xvel[iz,i] = numpy.mean(sim.vel[I,0])
       +            xvel[iz,i] = numpy.mean(numpy.abs(sim.vel[I,0]))
       +            #print numpy.mean(sim.vel[I,0])
       +            #xvel[iz,i] = numpy.abs(numpy.mean(sim.vel[I,0]))
       +            '''
       +
            if calculateforcechains:
                if i > 0:
                    loaded_contacts_prev = numpy.copy(loaded_contacts)
       t@@ -167,6 +184,7 @@ ax2.legend(lns, labs, loc='upper right', ncol=3,
                fancybox=True, framealpha=legend_alpha)
        ax1.set_ylim([-30, 200])
        #ax2.set_ylim(ax1.get_ylim())
       +ax2.set_ylim([-115,115])
        
        ax1.text(bbox_x, bbox_y, 'a',
                horizontalalignment=horizontalalignment,
       t@@ -239,18 +257,20 @@ ax7.text(bbox_x, bbox_y, 'd',
        
        ## ax9: porosity, ax10: unused
        ax9 = plt.subplot(5, 1, 5, sharex=ax1)
       -#poros_max = numpy.max(poros[0:sim.wall0iz(),:])
       -#poros_min = numpy.min(poros)
       -#poros_max = 0.44
        poros_max = 0.45
        poros_min = 0.37
        cmap = matplotlib.cm.get_cmap('Blues_r')
       +#cmap = matplotlib.cm.get_cmap('afmhot')
       +#im9 = ax9.pcolormesh(t, zpos_c, poros,
       +#zpos_c = zpos_c[:-1]
       +#xvel = xvel[:-1]
        im9 = ax9.pcolormesh(t, zpos_c, poros,
                cmap=cmap,
                #cmap=matplotlib.cm.get_cmap('bwr'),
                #cmap=matplotlib.cm.get_cmap('coolwarm'),
                #vmin=-p_ext, vmax=p_ext,
                vmin=poros_min, vmax=poros_max,
       +        #norm=matplotlib.colors.LogNorm(vmin=1.0e-10, vmax=xvel.max()),
                rasterized=True)
        ax9.set_ylim([zpos_c[0], sim.w_x[0]])
        ax9.set_ylabel('Vertical position [m]')
       t@@ -273,9 +293,13 @@ ax9.add_patch(matplotlib.patches.Rectangle(
            facecolor='white'))
        
        cb9 = plt.colorbar(im9, cax=cbaxes,
       -        ticks=[poros_min, poros_min + 0.5*(poros_max-poros_min), poros_max],
       +        ticks=[poros_min,
       +            poros_min + 0.5*(poros_max-poros_min),
       +            poros_max],
       +        #ticks=[xvel.min(), xvel.min() + 0.5*(xvel.max()-xvel.min()), xvel.max()],
                orientation='horizontal',
       -        extend='min', cmap=cmap)
       +        extend='min',
       +        cmap=cmap)
        cmap.set_under([8./255., 48./255., 107./255.])
        cb9.set_label('Mean horizontal porosity [-]')
        '''
 (DIR) diff --git a/python/halfshear-darcy-rate-dependence.py b/python/halfshear-darcy-rate-dependence.py
       t@@ -126,13 +126,23 @@ for sid in sids:
                    #c=shearstrain_nonzero[idxfit], linewidth=0.2,
                    #cmap=matplotlib.cm.get_cmap('afmhot'))
        
       -    ax1.plot(friction_fit, strainrate_fit)
       -    ax1.plot(friction_fit2, strainrate_fit2)
       +    ## plastic limit
       +    x = [0.3, 0.3]
       +    y = ax1.get_ylim()
       +    limitcolor = '#333333'
       +    ax1.plot(x, y, '--', linewidth=2, color=limitcolor)
       +    ax1.text(x[0]+0.03, 2.0e-4,
       +            'Yield strength', rotation=90, color=limitcolor,
       +            bbox={'fc':'#ffffff', 'pad':3, 'alpha':0.7})
       +
        
       +    ## Fit
       +    ax1.plot(friction_fit, strainrate_fit)
       +    #ax1.plot(friction_fit2, strainrate_fit2)
            ax1.annotate('$\\dot{\\gamma} = (\\tau/N)^{6.40}$',
                    xy = (friction_fit[40], strainrate_fit[40]),
       -            xytext = (0.32, 2.0e-7),
       -            arrowprops=dict(facecolor='blue', edgecolor='blue', shrink=0.05,
       +            xytext = (0.32+0.05, 2.0e-9),
       +            arrowprops=dict(facecolor='blue', edgecolor='blue', shrink=0.1,
                        width=1, headwidth=4, frac=0.2),)
                    #xytext = (friction_fit[50]+0.15, strainrate_fit[50]-1.0e-5))#,
                    #arrowprops=dict(facecolor='black', shrink=0.05),)
       t@@ -155,7 +165,6 @@ for sid in sids:
            #ax1.set_ylabel('Shear velocity [m/s]')
            ax1.set_ylabel('Shear strain rate $\\dot{\\gamma}$ [s$^{-1}$]')
        
       -
            plt.tight_layout()
            filename = sid + '-rate-dependence.' + outformat
            plt.savefig(filename)
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -4411,24 +4411,6 @@ class sim:
                video(self.sid, out_folder, video_format, graphics_folder, \
                      graphics_format, fps, qscale, bitrate, verbose)
        
       -    def shearVel(self):
       -        '''
       -        Calculates and returns the shear velocity (gamma_dot) of the
       -        experiment. The shear velocity is the x-axis velocity value of the
       -        upper particles.
       -
       -        :returns: The shear velocity applied by the upper, fixed particles [m/s]
       -        :return type: float
       -
       -        See also: :func:`shearStrainRate()` and :func:`shearStrain()`
       -        '''
       -
       -        # Find the fixed particles
       -        fixvel = numpy.nonzero(self.fixvel > 0.0)
       -
       -        # The shear velocity is the x-axis velocity value of the upper particles
       -        return self.vel[fixvel,0].max()
       -
            def shearDisplacement(self):
                '''
                Calculates and returns the current shear displacement. The displacement
       t@@ -4462,6 +4444,12 @@ class sim:
                fixvel = numpy.nonzero(self.fixvel > 0.0)
                return numpy.max(self.vel[fixvel,0])
        
       +    def shearVel(self):
       +        '''
       +        Alias of :func:`shearVelocity()`
       +        '''
       +        return self.shearVelocity()
       +
            def shearStrain(self):
                '''
                Calculates and returns the current shear strain (gamma) value of the