tapply normal stress, do not create visualization image if method isn't understood - 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 8e72814af7eb96dd3b6673e10ce7567023829e5d
 (DIR) parent ebbddb828cd9e760f1a9ba991298b15fc1cc8d34
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Thu,  2 Oct 2014 13:52:53 +0200
       
       apply normal stress, do not create visualization image if method isn't understood
       
       Diffstat:
         M python/halfshear-starter.py         |      31 ++++++++++++++++++-------------
         M python/shear-results-pressures.py   |      52 ++++++++++++++++++++-----------
         M python/sphere.py                    |       2 +-
       
       3 files changed, 53 insertions(+), 32 deletions(-)
       ---
 (DIR) diff --git a/python/halfshear-starter.py b/python/halfshear-starter.py
       t@@ -21,13 +21,17 @@ else:
            
        #sim = sphere.sim('diffusivity-sigma0=' + str(sigma0) +'-c_phi=1.0-c_grad_p=1.0',
        #        fluid=True)
       -sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-hw-relax', fluid=False)
       +sim = sphere.sim('halfshear-sigma0=' + str(sigma0), fluid=False)
       +print('Input: ' + sim.sid)
        sim.readlast()
        
       -print(sim.sid)
        sim.fluid = fluid
       +if fluid:
       +    sim.id('halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p) + '-shear')
       +else:
       +    sim.id('halfshear-sigma0=' + str(sigma0) + '-shear')
        
       -sim.checkerboardColors(nx=6,ny=6,nz=6)
       +sim.checkerboardColors(nx=6,ny=3,nz=6)
        sim.cleanup()
        sim.adjustUpperWall()
        sim.zeroKinematics()
       t@@ -39,24 +43,25 @@ if fluid:
            #sim.L[2] *= 2.0
            sim.initFluid(mu = 1.787e-6, p = 600.0e3, hydrostatic = True)
            #sim.initFluid(mu = 17.87e-4, p = 1.0e5, hydrostatic = True)
       -sim.setFluidBottomNoFlow()
       -sim.setFluidTopFixedPressure()
       -sim.setDEMstepsPerCFDstep(10)
       -sim.setMaxIterations(2e5)
       +    sim.setFluidBottomNoFlow()
       +    sim.setFluidTopFixedPressure()
       +    sim.setDEMstepsPerCFDstep(100)
       +    sim.setMaxIterations(2e5)
       +    sim.c_phi[0] = c_phi
       +    sim.c_grad_p[0] = c_grad_p
       +
        sim.initTemporal(total = 20.0, file_dt = 0.01, epsilon=0.07)
       -sim.c_phi[0] = c_phi
       -sim.c_grad_p[0] = c_grad_p
        sim.w_devs[0] = sigma0
       -#sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
       +sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
        sim.mu_s[0] = 0.5
        sim.mu_d[0] = 0.5
        sim.setDampingNormal(0.0)
        sim.setDampingTangential(0.0)
        
        # Fix lowermost particles
       -dz = sim.L[2]/sim.num[2]
       -I = numpy.nonzero(sim.x[:,2] < 1.5*dz)
       -sim.fixvel[I] = 1
       +#dz = sim.L[2]/sim.num[2]
       +#I = numpy.nonzero(sim.x[:,2] < 1.5*dz)
       +#sim.fixvel[I] = 1
        
        sim.run(dry=True)
        sim.run(device=device)
 (DIR) diff --git a/python/shear-results-pressures.py b/python/shear-results-pressures.py
       t@@ -34,7 +34,7 @@ for i in numpy.arange(sim.num[2]):
        shear_strain = numpy.zeros(sim.status())
        
        dev_pres = numpy.zeros((sim.num[2], sim.status()))
       -pres_static = numpy.zeros_like(dev_pres)
       +pres_static = numpy.ones_like(dev_pres)*100.0e3
        pres = numpy.zeros_like(dev_pres)
        
        for i in numpy.arange(sim.status()):
       t@@ -43,47 +43,63 @@ for i in numpy.arange(sim.status()):
        
            pres[:,i] = numpy.average(numpy.average(sim.p_f, axis=0), axis=0)
        
       -    for z in numpy.arange(0, sim.w_x[0]+1):
       +    wall0_iz = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))
       +    for z in numpy.arange(0, wall0_iz+1):
                pres_static[z,i] = \
                        (sim.w_x[0] - zpos_c[z])*sim.rho_f*numpy.abs(sim.g[2])\
                        + sim.p_f[0,0,-1]
       +        #pres_static[z,i] = zpos_c[z]
       +        #pres_static[z,i] = z
        
            shear_strain[i] = sim.shearStrain()
        
        dev_pres = pres - pres_static
        
       -#fig = plt.figure(figsize=(8,6))
       -fig = plt.figure(figsize=(8,15))
       +fig = plt.figure(figsize=(8,6))
       +#fig = plt.figure(figsize=(8,15))
        
       -ax1 = plt.subplot(311)
       -ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, rasterized=True)
       +min_p = numpy.min(dev_pres)/1000.0
       +#max_p = numpy.min(dev_pres)
       +max_p = numpy.abs(min_p)
       +
       +cmap = matplotlib.colors.ListedColormap(['b', 'w', 'r'])
       +bounds = [min_p, 0, max_p]
       +norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
       +
       +#ax1 = plt.subplot(311)
       +ax1 = plt.subplot(111)
       +im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, rasterized=True,
       +        cmap=cmap, norm=norm)
        ax1.set_xlim([0, shear_strain[-1]])
        ax1.set_ylim([zpos_c[0], sim.w_x[0]])
        ax1.set_xlabel('Shear strain $\\gamma$ [-]')
        ax1.set_ylabel('Vertical position $z$ [m]')
       -cb = plt.colorbar()
       -cb.set_label('Deviatoric pressure $p_\\text{f}$ [kPa]')
       -cb.solids.set_rasterized(True)
       +#cb1 = plt.colorbar(im1, boundaries=[min_p, numpy.abs(min_p)])
       +cb1 = plt.colorbar(cmap=cmap, norm=norm)
       +cb1.set_label('Deviatoric pressure $p_\\text{f}$ [kPa]')
       +cb1.solids.set_rasterized(True)
        
       +'''
        ax2 = plt.subplot(312)
       -ax2.pcolormesh(shear_strain, zpos_c, pres/1000.0, rasterized=True)
       +im2 = ax2.pcolormesh(shear_strain, zpos_c, pres/1000.0, rasterized=True)
        ax2.set_xlim([0, shear_strain[-1]])
        ax2.set_ylim([zpos_c[0], sim.w_x[0]])
        ax2.set_xlabel('Shear strain $\\gamma$ [-]')
        ax2.set_ylabel('Vertical position $z$ [m]')
       -cb = plt.colorbar()
       -cb.set_label('Pressure $p_\\text{f}$ [kPa]')
       -cb.solids.set_rasterized(True)
       +cb2 = plt.colorbar(im2)
       +cb2.set_label('Pressure $p_\\text{f}$ [kPa]')
       +cb2.solids.set_rasterized(True)
        
       -ax3 = plt.subplot(312)
       -ax3.pcolormesh(shear_strain, zpos_c, pres_static/1000.0, rasterized=True)
       +ax3 = plt.subplot(313)
       +im3 = ax3.pcolormesh(shear_strain, zpos_c, pres_static/1000.0, rasterized=True)
        ax3.set_xlim([0, shear_strain[-1]])
        ax3.set_ylim([zpos_c[0], sim.w_x[0]])
        ax3.set_xlabel('Shear strain $\\gamma$ [-]')
        ax3.set_ylabel('Vertical position $z$ [m]')
       -cb = plt.colorbar()
       -cb.set_label('Pressure $p_\\text{f}$ [kPa]')
       -cb.solids.set_rasterized(True)
       +cb3 = plt.colorbar(im3)
       +cb3.set_label('Static Pressure $p_\\text{f}$ [kPa]')
       +cb3.solids.set_rasterized(True)
       +'''
        
        
        #plt.MaxNLocator(nbins=4)
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -5185,7 +5185,7 @@ class sim:
        
                else :
                    print("Visualization type '" + method + "' not understood")
       -
       +            return
        
                # Optional save of figure
                if (outformat != 'txt'):