tvisualize function moved to sim class - 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 11b2e20f1da8ed07bdff711f0a65d009073af6a3
(DIR) parent 46da6614cd60cb329d7e0551e117db6fa605f2a0
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 27 Mar 2014 09:07:38 +0100
visualize function moved to sim class
Diffstat:
M python/sphere.py | 728 +++++++++++++++----------------
1 file changed, 362 insertions(+), 366 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -3852,6 +3852,368 @@ class sim:
'''
self.maxiter = numpy.asarray(maxiter)
+ def visualize(self, method = 'energy', savefig = True, outformat = 'png'):
+ '''
+ Visualize output from the simulation, where the temporal progress is
+ of interest. The output will be saved in the current folder with a name
+ combining the simulation id of the simulation, and the visualization
+ method.
+
+ :param method: The type of plot to render. Possible values are 'energy',
+ 'walls', 'triaxial' and 'shear'
+ :type method: str
+ :param savefig: Save the image instead of showing it on screen
+ :type savefig: bool
+ :param outformat: The output format of the plot data. This can be an
+ image format, or in text ('txt').
+ '''
+
+ lastfile = self.status()
+ sb = sim(sid = self.sid, np = self.np, nw = self.nw, fluid = self.fluid)
+
+ ### Plotting
+ if (outformat != 'txt'):
+ fig = plt.figure(figsize=(15,10),dpi=300)
+
+ if method == 'energy':
+
+ # Allocate arrays
+ Epot = numpy.zeros(lastfile+1)
+ Ekin = numpy.zeros(lastfile+1)
+ Erot = numpy.zeros(lastfile+1)
+ Es = numpy.zeros(lastfile+1)
+ Ev = numpy.zeros(lastfile+1)
+ Es_dot = numpy.zeros(lastfile+1)
+ Ev_dot = numpy.zeros(lastfile+1)
+ Ebondpot = numpy.zeros(lastfile+1)
+ Esum = numpy.zeros(lastfile+1)
+
+ # Read energy values from simulation binaries
+ for i in range(lastfile+1):
+ sb.readstep(i, verbose = False)
+
+ Epot[i] = sb.energy("pot")
+ Ekin[i] = sb.energy("kin")
+ Erot[i] = sb.energy("rot")
+ Es[i] = sb.energy("shear")
+ Ev[i] = sb.energy("visc_n")
+ Es_dot[i] = sb.energy("shearrate")
+ Ev_dot[i] = sb.energy("visc_n_rate")
+ Ebondpot[i] = sb.energy("bondpot")
+ Esum[i] = Epot[i] + Ekin[i] + Erot[i] + Es[i] + Ev[i] +\
+ Ebondpot[i]
+
+ t = numpy.linspace(0.0, sb.time_current, lastfile+1)
+
+ if (outformat != 'txt'):
+ # Potential energy
+ ax1 = plt.subplot2grid((2,5),(0,0))
+ ax1.set_xlabel('Time [s]')
+ ax1.set_ylabel('Total potential energy [J]')
+ ax1.plot(t, Epot, '+-')
+ ax1.grid()
+
+ # Kinetic energy
+ ax2 = plt.subplot2grid((2,5),(0,1))
+ ax2.set_xlabel('Time [s]')
+ ax2.set_ylabel('Total kinetic energy [J]')
+ ax2.plot(t, Ekin, '+-')
+ ax2.grid()
+
+ # Rotational energy
+ ax3 = plt.subplot2grid((2,5),(0,2))
+ ax3.set_xlabel('Time [s]')
+ ax3.set_ylabel('Total rotational energy [J]')
+ ax3.plot(t, Erot, '+-')
+ ax3.grid()
+
+ # Bond energy
+ ax4 = plt.subplot2grid((2,5),(0,3))
+ ax4.set_xlabel('Time [s]')
+ ax4.set_ylabel('Bond energy [J]')
+ ax4.plot(t, Ebondpot, '+-')
+ ax4.grid()
+
+ # Total energy
+ ax5 = plt.subplot2grid((2,5),(0,4))
+ ax5.set_xlabel('Time [s]')
+ ax5.set_ylabel('Total energy [J]')
+ ax5.plot(t, Esum, '+-')
+ ax5.grid()
+
+ # Shear energy rate
+ ax6 = plt.subplot2grid((2,5),(1,0))
+ ax6.set_xlabel('Time [s]')
+ ax6.set_ylabel('Frictional dissipation rate [W]')
+ ax6.plot(t, Es_dot, '+-')
+ ax6.grid()
+
+ # Shear energy
+ ax7 = plt.subplot2grid((2,5),(1,1))
+ ax7.set_xlabel('Time [s]')
+ ax7.set_ylabel('Total frictional dissipation [J]')
+ ax7.plot(t, Es, '+-')
+ ax7.grid()
+
+ # Visc_n energy rate
+ ax8 = plt.subplot2grid((2,5),(1,2))
+ ax8.set_xlabel('Time [s]')
+ ax8.set_ylabel('Viscous dissipation rate [W]')
+ ax8.plot(t, Ev_dot, '+-')
+ ax8.grid()
+
+ # Visc_n energy
+ ax9 = plt.subplot2grid((2,5),(1,3))
+ ax9.set_xlabel('Time [s]')
+ ax9.set_ylabel('Total viscous dissipation [J]')
+ ax9.plot(t, Ev, '+-')
+ ax9.grid()
+
+ # Combined view
+ ax10 = plt.subplot2grid((2,5),(1,4))
+ ax10.set_xlabel('Time [s]')
+ ax10.set_ylabel('Energy [J]')
+ ax10.plot(t, Epot, '+-g')
+ ax10.plot(t, Ekin, '+-b')
+ ax10.plot(t, Erot, '+-r')
+ ax10.legend(('$\sum E_{pot}$','$\sum E_{kin}$',\
+ '$\sum E_{rot}$'), 'upper right', shadow=True)
+ ax10.grid()
+
+ fig.tight_layout()
+
+ elif method == 'walls':
+
+ # Read energy values from simulation binaries
+ for i in range(lastfile+1):
+ sb.readstep(i, verbose=False)
+
+ # Allocate arrays on first run
+ if (i == 0):
+ wforce = numpy.zeros((lastfile+1)*sb.nw[0],\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
+ wvel = numpy.zeros((lastfile+1)*sb.nw[0],\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
+ wpos = numpy.zeros((lastfile+1)*sb.nw[0],\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
+ wdevs = numpy.zeros((lastfile+1)*sb.nw[0],\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
+ maxpos = numpy.zeros((lastfile+1), dtype=numpy.float64)
+ logstress = numpy.zeros((lastfile+1), dtype=numpy.float64)
+ voidratio = numpy.zeros((lastfile+1), dtype=numpy.float64)
+
+ wforce[i] = sb.w_force[0]
+ wvel[i] = sb.w_vel[0]
+ wpos[i] = sb.w_x[0]
+ wdevs[i] = sb.w_devs[0]
+ maxpos[i] = numpy.max(sb.x[:,2]+sb.radius)
+ logstress[i] =\
+ numpy.log((sb.w_force[0]/(sb.L[0]*sb.L[1]))/1000.0)
+ voidratio[i] = sb.voidRatio()
+
+ t = numpy.linspace(0.0, sb.time_current, lastfile+1)
+
+ # Plotting
+ if (outformat != 'txt'):
+ # linear plot of time vs. wall position
+ ax1 = plt.subplot2grid((2,2),(0,0))
+ ax1.set_xlabel('Time [s]')
+ ax1.set_ylabel('Position [m]')
+ ax1.plot(t, wpos, '+-', label="upper wall")
+ ax1.plot(t, maxpos, '+-', label="heighest particle")
+ ax1.legend()
+ ax1.grid()
+
+ #ax2 = plt.subplot2grid((2,2),(1,0))
+ #ax2.set_xlabel('Time [s]')
+ #ax2.set_ylabel('Force [N]')
+ #ax2.plot(t, wforce, '+-')
+
+ # semilog plot of log stress vs. void ratio
+ ax2 = plt.subplot2grid((2,2),(1,0))
+ ax2.set_xlabel('log deviatoric stress [kPa]')
+ ax2.set_ylabel('Void ratio [-]')
+ ax2.plot(logstress, voidratio, '+-')
+ ax2.grid()
+
+ # linear plot of time vs. wall velocity
+ ax3 = plt.subplot2grid((2,2),(0,1))
+ ax3.set_xlabel('Time [s]')
+ ax3.set_ylabel('Velocity [m/s]')
+ ax3.plot(t, wvel, '+-')
+ ax3.grid()
+
+ # linear plot of time vs. deviatoric stress
+ ax4 = plt.subplot2grid((2,2),(1,1))
+ ax4.set_xlabel('Time [s]')
+ ax4.set_ylabel('Deviatoric stress [Pa]')
+ ax4.plot(t, wdevs, '+-', label="$\sigma_0$")
+ ax4.plot(t, wforce/(sb.L[0]*sb.L[1]), '+-', label="$\sigma'$")
+ ax4.legend(loc=4)
+ ax4.grid()
+
+ elif method == 'triaxial':
+
+ # Read energy values from simulation binaries
+ for i in range(lastfile+1):
+ sb.readstep(i, verbose = False)
+
+ vol = (sb.w_x[0]-sb.origo[2]) * (sb.w_x[1]-sb.w_x[2]) \
+ * (sb.w_x[3] - sb.w_x[4])
+
+ # Allocate arrays on first run
+ if (i == 0):
+ axial_strain = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ deviatoric_stress =\
+ numpy.zeros(lastfile+1, dtype=numpy.float64)
+ volumetric_strain =\
+ numpy.zeros(lastfile+1, dtype=numpy.float64)
+
+ w0pos0 = sb.w_x[0]
+ vol0 = vol
+
+ sigma1 = sb.w_force[0]/\
+ ((sb.w_x[1]-sb.w_x[2])*(sb.w_x[3]-sb.w_x[4]))
+
+ axial_strain[i] = (w0pos0 - sb.w_x[0])/w0pos0
+ volumetric_strain[i] = (vol0-vol)/vol0
+ deviatoric_stress[i] = sigma1 / sb.w_devs[1]
+
+ #print(lastfile)
+ #print(axial_strain)
+ #print(deviatoric_stress)
+ #print(volumetric_strain)
+
+ # Plotting
+ if (outformat != 'txt'):
+
+ # linear plot of deviatoric stress
+ ax1 = plt.subplot2grid((2,1),(0,0))
+ ax1.set_xlabel('Axial strain, $\gamma_1$, [-]')
+ ax1.set_ylabel('Deviatoric stress, $\sigma_1 - \sigma_3$, [Pa]')
+ ax1.plot(axial_strain, deviatoric_stress, '+-')
+ #ax1.legend()
+ ax1.grid()
+
+ #ax2 = plt.subplot2grid((2,2),(1,0))
+ #ax2.set_xlabel('Time [s]')
+ #ax2.set_ylabel('Force [N]')
+ #ax2.plot(t, wforce, '+-')
+
+ # semilog plot of log stress vs. void ratio
+ ax2 = plt.subplot2grid((2,1),(1,0))
+ ax2.set_xlabel('Axial strain, $\gamma_1$ [-]')
+ ax2.set_ylabel('Volumetric strain, $\gamma_v$, [-]')
+ ax2.plot(axial_strain, volumetric_strain, '+-')
+ ax2.grid()
+
+
+ elif method == 'shear':
+
+ # Read stress values from simulation binaries
+ for i in range(lastfile+1):
+ sb.readstep(i, verbose = False)
+
+ # First iteration: Allocate arrays and find constant values
+ if (i == 0):
+ # Shear displacement
+ xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
+
+ # Normal stress
+ sigma_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
+
+ # Normal stress
+ sigma_def = numpy.zeros(lastfile+1, dtype=numpy.float64)
+
+ # Shear stress
+ tau = numpy.zeros(lastfile+1, dtype=numpy.float64)
+
+ # Upper wall position
+ dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
+
+ # Upper wall position
+ tau_u = 0.0 # Peak shear stress
+ # Shear strain value of peak sh. stress
+ tau_u_shearstrain = 0.0
+
+ fixvel = numpy.nonzero(sb.fixvel > 0.0)
+ #fixvel_upper = numpy.nonzero(sb.vel[fixvel,0] > 0.0)
+ shearvel = sb.vel[fixvel,0].max()
+ w_x0 = sb.w_x[0] # Original height
+ A = sb.L[0] * sb.L[1] # Upper surface area
+
+ # Summation of shear stress contributions
+ for j in fixvel[0]:
+ if (sb.vel[j,0] > 0.0):
+ tau[i] += -sb.force[j,0]
+
+ if (i > 0):
+ xdisp[i] = xdisp[i-1] + sb.time_file_dt[0] * shearvel
+ sigma_eff[i] = sb.w_force[0] / A
+ sigma_def[i] = sb.w_devs[0]
+ dilation[i] = sb.w_x[0] - w_x0 # dilation in meters
+ #dilation[i] = (sb.w_x[0] - w_x0)/w_x0 * 100.0 # dilation in percent
+
+ # Test if this was the max. shear stress
+ if (tau[i] > tau_u):
+ tau_u = tau[i]
+ tau_u_shearstrain = xdisp[i]/w_x0
+
+
+ # Plot stresses
+ if (outformat != 'txt'):
+ shearinfo = "$\\tau_u$ = {:.3} Pa at $\gamma$ = {:.3}".format(\
+ tau_u, tau_u_shearstrain)
+ fig.text(0.5, 0.03, shearinfo, horizontalalignment='center',
+ fontproperties=FontProperties(size=14))
+ ax1 = plt.subplot2grid((1, 2), (0, 0))
+ ax1.set_xlabel('Shear strain [-]')
+ ax1.set_ylabel('Stress [Pa]')
+ ax1.plot(xdisp / w_x0, sigma_eff, '+-g', label="$\sigma'$")
+ ax1.plot(xdisp / w_x0, sigma_def, '+-b', label="$\sigma_0$")
+ ax1.plot(xdisp / w_x0, tau, '+-r', label="$\\tau$")
+ ax1.legend(loc=4)
+ ax1.grid()
+
+ # Plot dilation
+ ax2 = plt.subplot2grid((1,2),(0,1))
+ ax2.set_xlabel('Shear strain [-]')
+ ax2.set_ylabel('Dilation [m]')
+ ax2.plot(xdisp/w_x0, dilation, '+-')
+ ax2.grid()
+
+ else :
+ # Write values to textfile
+ filename = "shear-stresses-{0}.txt".format(self.sid)
+ #print("Writing stress data to " + filename)
+ fh = None
+ try :
+ fh = open(filename, "w")
+ L = sb.L[2] - sb.origo[2] # Initial height
+ for i in range(lastfile+1):
+ # format: shear distance [mm], sigma [kPa], tau [kPa],
+ # Dilation [%]
+ fh.write("{0}\t{1}\t{2}\t{3}\n".format(xdisp[i],
+ sigma_eff[i]/1000.0,
+ tau[i]/1000.0,
+ dilation[i]))
+ finally :
+ if fh is not None:
+ fh.close()
+
+
+ else :
+ print("Visualization type '" + method + "' not understood")
+
+
+ # Optional save of figure
+ if (outformat != 'txt'):
+ if (savefig == True):
+ fig.savefig("{0}-{1}.{2}".format(self.sid, method, outformat))
+ fig.clf()
+ else :
+ plt.show()
+
def convert(graphics_format = 'png', folder = '../img_out'):
'''
t@@ -4020,372 +4382,6 @@ def thinsectionVideo(project,
+ out_folder + "/" + project + "-ts-x1x3." + video_format, \
shell=True)
-def visualize(project, method = 'energy', savefig = True, outformat = 'png'):
- '''
- Visualize output from the target project, where the temporal progress is of
- interest. The output will be saved in the current folder with a name
- combining the simulation id of the project, and the visualization method.
-
- :param project: The simulation id of the project to render
- :type project: str
- :param method: The type of plot to render. Possible values are 'energy',
- 'walls', 'triaxial' and 'shear'
- :type method: str
- :param savefig: Save the image instead of showing it on screen
- :type savefig: bool
- :param outformat: The output format of the plot data. This can be an image
- format, or in text ('txt').
- '''
-
- lastfile = status(project)
-
- ### Plotting
- if (outformat != 'txt'):
- fig = plt.figure(figsize=(15,10),dpi=300)
-
- if method == 'energy':
-
- # Allocate arrays
- Epot = numpy.zeros(lastfile+1)
- Ekin = numpy.zeros(lastfile+1)
- Erot = numpy.zeros(lastfile+1)
- Es = numpy.zeros(lastfile+1)
- Ev = numpy.zeros(lastfile+1)
- Es_dot = numpy.zeros(lastfile+1)
- Ev_dot = numpy.zeros(lastfile+1)
- Ebondpot = numpy.zeros(lastfile+1)
- Esum = numpy.zeros(lastfile+1)
-
- # Read energy values from project binaries
- sb = sim(fluid = self.fluid)
- for i in range(lastfile+1):
- fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
- sb.readbin(fn, verbose = False)
-
- Epot[i] = sb.energy("pot")
- Ekin[i] = sb.energy("kin")
- Erot[i] = sb.energy("rot")
- Es[i] = sb.energy("shear")
- Ev[i] = sb.energy("visc_n")
- Es_dot[i] = sb.energy("shearrate")
- Ev_dot[i] = sb.energy("visc_n_rate")
- Ebondpot[i] = sb.energy("bondpot")
- Esum[i] = Epot[i] + Ekin[i] + Erot[i] + Es[i] + Ev[i] + Ebondpot[i]
-
- t = numpy.linspace(0.0, sb.time_current, lastfile+1)
-
- if (outformat != 'txt'):
- # Potential energy
- ax1 = plt.subplot2grid((2,5),(0,0))
- ax1.set_xlabel('Time [s]')
- ax1.set_ylabel('Total potential energy [J]')
- ax1.plot(t, Epot, '+-')
- ax1.grid()
-
- # Kinetic energy
- ax2 = plt.subplot2grid((2,5),(0,1))
- ax2.set_xlabel('Time [s]')
- ax2.set_ylabel('Total kinetic energy [J]')
- ax2.plot(t, Ekin, '+-')
- ax2.grid()
-
- # Rotational energy
- ax3 = plt.subplot2grid((2,5),(0,2))
- ax3.set_xlabel('Time [s]')
- ax3.set_ylabel('Total rotational energy [J]')
- ax3.plot(t, Erot, '+-')
- ax3.grid()
-
- # Bond energy
- ax4 = plt.subplot2grid((2,5),(0,3))
- ax4.set_xlabel('Time [s]')
- ax4.set_ylabel('Bond energy [J]')
- ax4.plot(t, Ebondpot, '+-')
- ax4.grid()
-
- # Total energy
- ax5 = plt.subplot2grid((2,5),(0,4))
- ax5.set_xlabel('Time [s]')
- ax5.set_ylabel('Total energy [J]')
- ax5.plot(t, Esum, '+-')
- ax5.grid()
-
- # Shear energy rate
- ax6 = plt.subplot2grid((2,5),(1,0))
- ax6.set_xlabel('Time [s]')
- ax6.set_ylabel('Frictional dissipation rate [W]')
- ax6.plot(t, Es_dot, '+-')
- ax6.grid()
-
- # Shear energy
- ax7 = plt.subplot2grid((2,5),(1,1))
- ax7.set_xlabel('Time [s]')
- ax7.set_ylabel('Total frictional dissipation [J]')
- ax7.plot(t, Es, '+-')
- ax7.grid()
-
- # Visc_n energy rate
- ax8 = plt.subplot2grid((2,5),(1,2))
- ax8.set_xlabel('Time [s]')
- ax8.set_ylabel('Viscous dissipation rate [W]')
- ax8.plot(t, Ev_dot, '+-')
- ax8.grid()
-
- # Visc_n energy
- ax9 = plt.subplot2grid((2,5),(1,3))
- ax9.set_xlabel('Time [s]')
- ax9.set_ylabel('Total viscous dissipation [J]')
- ax9.plot(t, Ev, '+-')
- ax9.grid()
-
-
- # Combined view
- ax10 = plt.subplot2grid((2,5),(1,4))
- ax10.set_xlabel('Time [s]')
- ax10.set_ylabel('Energy [J]')
- ax10.plot(t, Epot, '+-g')
- ax10.plot(t, Ekin, '+-b')
- ax10.plot(t, Erot, '+-r')
- ax10.legend(('$\sum E_{pot}$','$\sum E_{kin}$','$\sum E_{rot}$'),\
- 'upper right', shadow=True)
- ax10.grid()
-
- fig.tight_layout()
-
- elif method == 'walls':
-
- # Read energy values from project binaries
- sb = sim(fluid = self.fluid)
- for i in range(lastfile+1):
- fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
- sb.readbin(fn, verbose = False)
-
- # Allocate arrays on first run
- if (i == 0):
- wforce = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
- wvel = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
- wpos = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
- wdevs = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
- maxpos = numpy.zeros((lastfile+1), dtype=numpy.float64)
- logstress = numpy.zeros((lastfile+1), dtype=numpy.float64)
- voidratio = numpy.zeros((lastfile+1), dtype=numpy.float64)
-
- wforce[i] = sb.w_force[0]
- wvel[i] = sb.w_vel[0]
- wpos[i] = sb.w_x[0]
- wdevs[i] = sb.w_devs[0]
- maxpos[i] = numpy.max(sb.x[:,2]+sb.radius)
- logstress[i] = numpy.log((sb.w_force[0]/(sb.L[0]*sb.L[1]))/1000.0)
- voidratio[i] = sb.voidRatio()
-
-
- t = numpy.linspace(0.0, sb.time_current, lastfile+1)
-
- # Plotting
- if (outformat != 'txt'):
- # linear plot of time vs. wall position
- ax1 = plt.subplot2grid((2,2),(0,0))
- ax1.set_xlabel('Time [s]')
- ax1.set_ylabel('Position [m]')
- ax1.plot(t, wpos, '+-', label="upper wall")
- ax1.plot(t, maxpos, '+-', label="heighest particle")
- ax1.legend()
- ax1.grid()
-
- #ax2 = plt.subplot2grid((2,2),(1,0))
- #ax2.set_xlabel('Time [s]')
- #ax2.set_ylabel('Force [N]')
- #ax2.plot(t, wforce, '+-')
-
- # semilog plot of log stress vs. void ratio
- ax2 = plt.subplot2grid((2,2),(1,0))
- ax2.set_xlabel('log deviatoric stress [kPa]')
- ax2.set_ylabel('Void ratio [-]')
- ax2.plot(logstress, voidratio, '+-')
- ax2.grid()
-
- # linear plot of time vs. wall velocity
- ax3 = plt.subplot2grid((2,2),(0,1))
- ax3.set_xlabel('Time [s]')
- ax3.set_ylabel('Velocity [m/s]')
- ax3.plot(t, wvel, '+-')
- ax3.grid()
-
- # linear plot of time vs. deviatoric stress
- ax4 = plt.subplot2grid((2,2),(1,1))
- ax4.set_xlabel('Time [s]')
- ax4.set_ylabel('Deviatoric stress [Pa]')
- ax4.plot(t, wdevs, '+-', label="$\sigma_0$")
- ax4.plot(t, wforce/(sb.L[0]*sb.L[1]), '+-', label="$\sigma'$")
- ax4.legend(loc=4)
- ax4.grid()
-
- elif method == 'triaxial':
-
- # Read energy values from project binaries
- sb = sim(fluid = self.fluid)
- for i in range(lastfile+1):
- fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
- sb.readbin(fn, verbose = False)
-
- vol = (sb.w_x[0]-sb.origo[2]) * (sb.w_x[1]-sb.w_x[2]) \
- * (sb.w_x[3] - sb.w_x[4])
-
- # Allocate arrays on first run
- if (i == 0):
- axial_strain = numpy.zeros(lastfile+1, dtype=numpy.float64)
- deviatoric_stress = numpy.zeros(lastfile+1, dtype=numpy.float64)
- volumetric_strain = numpy.zeros(lastfile+1, dtype=numpy.float64)
-
- w0pos0 = sb.w_x[0]
- vol0 = vol
-
- sigma1 = sb.w_force[0]/((sb.w_x[1]-sb.w_x[2])*(sb.w_x[3]-sb.w_x[4]))
-
- axial_strain[i] = (w0pos0 - sb.w_x[0])/w0pos0
- volumetric_strain[i] = (vol0-vol)/vol0
- deviatoric_stress[i] = sigma1 / sb.w_devs[1]
-
- #print(lastfile)
- #print(axial_strain)
- #print(deviatoric_stress)
- #print(volumetric_strain)
-
- # Plotting
- if (outformat != 'txt'):
-
- # linear plot of deviatoric stress
- ax1 = plt.subplot2grid((2,1),(0,0))
- ax1.set_xlabel('Axial strain, $\gamma_1$, [-]')
- ax1.set_ylabel('Deviatoric stress, $\sigma_1 - \sigma_3$, [Pa]')
- ax1.plot(axial_strain, deviatoric_stress, '+-')
- #ax1.legend()
- ax1.grid()
-
- #ax2 = plt.subplot2grid((2,2),(1,0))
- #ax2.set_xlabel('Time [s]')
- #ax2.set_ylabel('Force [N]')
- #ax2.plot(t, wforce, '+-')
-
- # semilog plot of log stress vs. void ratio
- ax2 = plt.subplot2grid((2,1),(1,0))
- ax2.set_xlabel('Axial strain, $\gamma_1$ [-]')
- ax2.set_ylabel('Volumetric strain, $\gamma_v$, [-]')
- ax2.plot(axial_strain, volumetric_strain, '+-')
- ax2.grid()
-
-
- elif method == 'shear':
-
- sb = sim(fluid = self.fluid)
- # Read stress values from project binaries
- for i in range(lastfile+1):
- fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
- sb.readbin(fn, verbose = False)
-
- # First iteration: Allocate arrays and find constant values
- if (i == 0):
- # Shear displacement
- xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
-
- # Normal stress
- sigma_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
-
- # Normal stress
- sigma_def = numpy.zeros(lastfile+1, dtype=numpy.float64)
-
- # Shear stress
- tau = numpy.zeros(lastfile+1, dtype=numpy.float64)
-
- # Upper wall position
- dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
-
- # Upper wall position
- tau_u = 0.0 # Peak shear stress
- tau_u_shearstrain = 0.0 # Shear strain value of peak sh. stress
-
- fixvel = numpy.nonzero(sb.fixvel > 0.0)
- #fixvel_upper = numpy.nonzero(sb.vel[fixvel,0] > 0.0)
- shearvel = sb.vel[fixvel,0].max()
- w_x0 = sb.w_x[0] # Original height
- A = sb.L[0] * sb.L[1] # Upper surface area
-
- # Summation of shear stress contributions
- for j in fixvel[0]:
- if (sb.vel[j,0] > 0.0):
- tau[i] += -sb.force[j,0]
-
- if (i > 0):
- xdisp[i] = xdisp[i-1] + sb.time_file_dt[0] * shearvel
- sigma_eff[i] = sb.w_force[0] / A
- sigma_def[i] = sb.w_devs[0]
- dilation[i] = sb.w_x[0] - w_x0 # dilation in meters
- #dilation[i] = (sb.w_x[0] - w_x0)/w_x0 * 100.0 # dilation in percent
-
- # Test if this was the max. shear stress
- if (tau[i] > tau_u):
- tau_u = tau[i]
- tau_u_shearstrain = xdisp[i]/w_x0
-
-
- # Plot stresses
- if (outformat != 'txt'):
- shearinfo = "$\\tau_u$ = {:.3} Pa at $\gamma$ = {:.3}".format(\
- tau_u, tau_u_shearstrain)
- fig.text(0.5, 0.03, shearinfo, horizontalalignment='center',
- fontproperties=FontProperties(size=14))
- ax1 = plt.subplot2grid((1, 2), (0, 0))
- ax1.set_xlabel('Shear strain [-]')
- ax1.set_ylabel('Stress [Pa]')
- ax1.plot(xdisp / w_x0, sigma_eff, '+-g', label="$\sigma'$")
- ax1.plot(xdisp / w_x0, sigma_def, '+-b', label="$\sigma_0$")
- ax1.plot(xdisp / w_x0, tau, '+-r', label="$\\tau$")
- ax1.legend(loc=4)
- ax1.grid()
-
- # Plot dilation
- ax2 = plt.subplot2grid((1,2),(0,1))
- ax2.set_xlabel('Shear strain [-]')
- ax2.set_ylabel('Dilation [m]')
- ax2.plot(xdisp/w_x0, dilation, '+-')
- ax2.grid()
-
- else :
- # Write values to textfile
- filename = "shear-stresses-{0}.txt".format(project)
- #print("Writing stress data to " + filename)
- fh = None
- try :
- fh = open(filename, "w")
- L = sb.L[2] - sb.origo[2] # Initial height
- for i in range(lastfile+1):
- # format: shear distance [mm], sigma [kPa], tau [kPa],
- # Dilation [%]
- fh.write("{0}\t{1}\t{2}\t{3}\n".format(xdisp[i],
- sigma_eff[i]/1000.0,
- tau[i]/1000.0,
- dilation[i]))
- finally :
- if fh is not None:
- fh.close()
-
-
- else :
- print("Visualization type '" + method + "' not understood")
-
-
- # Optional save of figure
- if (outformat != 'txt'):
- if (savefig == True):
- fig.savefig("{0}-{1}.{2}".format(project, method, outformat))
- fig.clf()
- else :
- plt.show()
-
def run(binary, verbose=True, hideinputfile=False):
'''
Execute ``sphere`` with target binary file as input.