timproved plots with means and porosities - 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 ce6e548cf6900e585d97dbe8b5482a428c96600e
 (DIR) parent 41f960f604e07ec961ae9fa30468295627f7eb4e
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Tue,  9 Sep 2014 15:40:37 +0200
       
       improved plots with means and porosities
       
       Diffstat:
         M python/shear-results-forces.py      |      56 ++++++++++++++++++++++++++-----
       
       1 file changed, 48 insertions(+), 8 deletions(-)
       ---
 (DIR) diff --git a/python/shear-results-forces.py b/python/shear-results-forces.py
       t@@ -44,6 +44,15 @@ f_pf  = numpy.zeros_like(xdisp)
        # pressure - hydrostatic pressure
        dev_p = numpy.zeros((len(steps), sim.num[2]))
        
       +# mean porosity
       +phi_bar = numpy.zeros((len(steps), sim.num[2]))
       +
       +# mean per-particle values
       +xdisp_mean = numpy.zeros((len(steps), sim.num[2]))
       +f_pf_mean = numpy.zeros((len(steps), sim.num[2]))
       +
       +shear_strain = numpy.zeros(len(steps))
       +
        s = 0
        for step in steps:
        
       t@@ -73,6 +82,21 @@ for step in steps:
                            numpy.average(numpy.average(sim.p_f, axis=0), axis=0)\
                            /nsteps_avg
        
       +            phi_bar[s,:] += \
       +                    numpy.average(numpy.average(sim.phi, axis=0), axis=0)\
       +                    /nsteps_avg
       +
       +            shear_strain[s] += sim.shearStrain()/nsteps_avg
       +
       +        # calculate mean values of xdisp and f_pf
       +        for iz in numpy.arange(sim.num[2]):
       +            z_bot = iz*dz
       +            z_top = (iz+1)*dz
       +            I = numpy.nonzero((zpos_p[s,:] >= z_bot) & (zpos_p[s,:] < z_top))
       +            if len(I) > 0:
       +                xdisp_mean[s,iz] = numpy.mean(xdisp[s,I])
       +                f_pf_mean[s,iz] = numpy.mean(f_pf[s,I])
       +
            else:
                print(sid + ' not found')
            s += 1
       t@@ -83,10 +107,20 @@ for s in numpy.arange(len(steps)):
            ax1 = plt.subplot((s+1)*100 + 31)
            ax2 = plt.subplot((s+1)*100 + 32, sharey=ax1)
            ax3 = plt.subplot((s+1)*100 + 33, sharey=ax1)
       +    ax4 = ax3.twiny()
       +
       +    ax1.plot(xdisp[s], zpos_p[s], '+', color = '#888888')
       +    ax1.plot(xdisp_mean[s], zpos_c[s], color = 'k')
       +
       +    ax2.plot(f_pf[s],  zpos_p[s], '+', color = '#888888')
       +    ax2.plot(f_pf_mean[s], zpos_c[s], color = 'k')
        
       -    ax1.plot(xdisp[s], zpos_p[s], '+')
       -    ax2.plot(f_pf[s],  zpos_p[s], '+')
       -    ax3.plot(dev_p[s]/1000.0, zpos_c[s])
       +    ax3.plot(dev_p[s]/1000.0, zpos_c[s], 'k')
       +
       +    phicolor = '#666666'
       +    ax4.plot(phi_bar[s], zpos_c[s], '--', color = phicolor)
       +    for tl in ax4.get_xticklabels():
       +        tl.set_color(phicolor)
        
            max_z = numpy.max(zpos_p)
            ax1.set_ylim([0, max_z])
       t@@ -98,26 +132,32 @@ for s in numpy.arange(len(steps)):
            #plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
        
            ax1.set_ylabel('Vertical position $z$ [m]')
       -    ax1.set_xlabel('$x^3_p$ [m]')
       +    ax1.set_xlabel('$x^3_\\text{p}$ [m]')
            ax2.set_xlabel('$\\boldsymbol{f}_\\text{pf}$ [N]')
            ax3.set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
       +    ax4.set_xlabel('$\\bar{\\phi}$ [-]', color=phicolor)
            plt.setp(ax2.get_yticklabels(), visible=False)
            plt.setp(ax3.get_yticklabels(), visible=False)
        
       -    ax1.get_xaxis().set_major_locator(MaxNLocator(nbins=4))
       -    ax2.get_xaxis().set_major_locator(MaxNLocator(nbins=4))
       -    ax3.get_xaxis().set_major_locator(MaxNLocator(nbins=4))
       +    ax1.get_xaxis().set_major_locator(MaxNLocator(nbins=5))
       +    ax2.get_xaxis().set_major_locator(MaxNLocator(nbins=5))
       +    ax3.get_xaxis().set_major_locator(MaxNLocator(nbins=5))
        
            plt.setp(ax1.xaxis.get_majorticklabels(), rotation=90)
            plt.setp(ax2.xaxis.get_majorticklabels(), rotation=90)
            plt.setp(ax3.xaxis.get_majorticklabels(), rotation=90)
       +    plt.setp(ax4.xaxis.get_majorticklabels(), rotation=90)
       +
       +    fig.text(0.1, 0.9,
       +            'Shear strain $\\gamma = %.3f$' % (shear_strain[s]),
       +            horizontalalignment='left', fontsize=22)
            #ax1.grid()
            #ax2.grid()
            #ax1.legend(loc='lower right', prop={'size':18})
            #ax2.legend(loc='lower right', prop={'size':18})
        
        plt.tight_layout()
       -plt.subplots_adjust(wspace = .001)
       +plt.subplots_adjust(wspace = .05)
        plt.MaxNLocator(nbins=4)
        filename = 'shear-10kPa-forces.pdf'
        plt.savefig(filename)