tadd option to return either defined or effective normal stress - 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 ff9b1797587c11d741b216e15ea3ac180366b9b2
 (DIR) parent dfce171595e257cf25c6fc52744f673663865807
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Tue, 24 Feb 2015 10:48:27 +0100
       
       add option to return either defined or effective normal stress
       
       Diffstat:
         M python/halfshear-darcy-forcechainm… |       2 +-
         M python/sphere.py                    |      39 ++++++++++++++++++++++++-------
       
       2 files changed, 32 insertions(+), 9 deletions(-)
       ---
 (DIR) diff --git a/python/halfshear-darcy-forcechainmapper.py b/python/halfshear-darcy-forcechainmapper.py
       t@@ -82,7 +82,7 @@ for sid in sids:
            ax2.semilogy(t, n - nkept, 'k')
            ax2.set_xlabel('Time $t$ [s]')
            ax2.set_ylabel(\
       -            'New heavily loaded contacts $\sigma_c \\geq \\sigma_0\\times 4$ Pa')
       +            'New heavily loaded contacts')
                    #'Heavily loaded contacts $||\\boldsymbol{{\sigma}}_c|| \\geq \\sigma_0\\times 4$ Pa')
        
            ax3 = plt.subplot(3,1,3)
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -3736,16 +3736,24 @@ class sim:
                # Increment the number of bonds with one
                self.nb0 += 1
        
       -    def currentNormalStress(self):
       +    def currentNormalStress(self, type='defined'):
                '''
       -        Calculates the current magnitude of the top wall normal stress.
       +        Calculates the current magnitude of the defined or effective top wall
       +        normal stress.
       +
       +        :param type: Find the 'defined' (default) or 'effective' normal stress 
       +        :type type: str
        
                :returns: The current top wall normal stress in Pascal
                :return type: float
                '''
       -        return self.w_sigma0[0] \
       -                + self.w_sigma0_A[0] \
       -                *numpy.sin(2.0*numpy.pi*self.w_sigma0_f[0]*self.time_current[0])
       +        if type == 'defined':
       +            return self.w_sigma0[0] \
       +                    + self.w_sigma0_A[0] \
       +                    *numpy.sin(2.0*numpy.pi*self.w_sigma0_f[0]\
       +                    *self.time_current[0])
       +        elif type == 'effective':
       +            return self.w_force[0]/(self.L[0]*self.L[1])
        
            def volume(self, idx):
                '''
       t@@ -4378,6 +4386,23 @@ class sim:
                # 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
       +        is found by determining the total x-axis displacement of the upper,
       +        fixed particles.
       +
       +        :returns: The total shear displacement [m]
       +        :return type: float
       +
       +        See also: :func:`shearStrain()`
       +        '''
       +
       +        # Displacement of the upper, fixed particles in the shear direction
       +        #xdisp = self.time_current[0] * self.shearVel()
       +        fixvel = numpy.nonzero(self.fixvel > 0.0)
       +        return numpy.max(self.xyzsum[fixvel,0])
       +
            def shearStrain(self):
                '''
                Calculates and returns the current shear strain (gamma) value of the
       t@@ -4394,9 +4419,7 @@ class sim:
                w_x0 = self.w_x[0]
        
                # Displacement of the upper, fixed particles in the shear direction
       -        #xdisp = self.time_current[0] * self.shearVel()
       -        fixvel = numpy.nonzero(self.fixvel > 0.0)
       -        xdisp = numpy.max(self.xyzsum[fixvel,0])
       +        xdisp = self.shearDisplacement()
        
                # Return shear strain
                return xdisp/w_x0