thalfshear-darcy-forcechainmapper.py - 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
       ---
       thalfshear-darcy-forcechainmapper.py (2842B)
       ---
            1 #!/usr/bin/env python
            2 import matplotlib
            3 matplotlib.use('Agg')
            4 #matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
            5 matplotlib.rcParams.update({'font.size': 7, 'font.family': 'sans-serif'})
            6 matplotlib.rc('text', usetex=True)
            7 matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
            8 import shutil
            9 
           10 import os
           11 import sys
           12 import numpy
           13 import sphere
           14 import matplotlib.pyplot as plt
           15 
           16 sids =\
           17         ['halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2']
           18         #['halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.797e-06-ss=10000.0-A=70000.0-f=0.2']
           19 outformat = 'pdf'
           20 fluid = True
           21 #threshold = 100.0 # [N]
           22 
           23 
           24 
           25 for sid in sids:
           26 
           27     sim = sphere.sim(sid, fluid=fluid)
           28 
           29     #nsteps = 2
           30     #nsteps = 10
           31     nsteps = sim.status()
           32 
           33     t = numpy.empty(nsteps)
           34     n = numpy.empty(nsteps)
           35     coordinationnumber = numpy.empty(nsteps)
           36     nkept = numpy.empty(nsteps)
           37 
           38     for i in numpy.arange(nsteps):
           39         sim.readstep(i+1)
           40         t[i] = sim.currentTime()
           41         #n[i] = countLoadedContacts(sim, threshold)
           42 
           43         if i > 0:
           44             loaded_contacts_prev = numpy.copy(loaded_contacts)
           45             pairs_prev = numpy.copy(sim.pairs)
           46 
           47         loaded_contacts = sim.findLoadedContacts(
           48                 threshold = sim.currentNormalStress()*2.)
           49                 #sim.currentNormalStress()/1000.)
           50         n[i] = numpy.size(loaded_contacts)
           51 
           52         sim.findCoordinationNumber()
           53         coordinationnumber[i] = sim.findMeanCoordinationNumber()
           54 
           55         nfound = 0
           56         if i > 0:
           57             for a in loaded_contacts[0]:
           58                 for b in loaded_contacts_prev[0]:
           59                     if (sim.pairs[:,a] == pairs_prev[:,b]).all():
           60                         nfound += 1;
           61 
           62         nkept[i] = nfound
           63 
           64         print coordinationnumber[i]
           65         print nfound
           66 
           67     #fig = plt.figure(figsize=[8,8])
           68     fig = plt.figure(figsize=[3.5,7])
           69 
           70     ax1 = plt.subplot(3,1,1)
           71     #ax1.plot(t, n)
           72     ax1.semilogy(t, n, 'k')
           73     ax1.set_xlabel('Time $t$ [s]')
           74     #ax1.set_ylabel(\
           75             #'Heavily loaded contacts $||\\boldsymbol{{f}}_n|| \\geq {}$ N'.format(
           76                 #threshold))
           77     ax1.set_ylabel(\
           78             'Heavily loaded contacts $\sigma_c \\geq \\sigma_0\\times 4$ Pa')
           79             #'Heavily loaded contacts $||\\boldsymbol{{\sigma}}_c|| \\geq \\sigma_0\\times 4$ Pa')
           80 
           81     ax2 = plt.subplot(3,1,2)
           82     ax2.semilogy(t, n - nkept, 'k')
           83     ax2.set_xlabel('Time $t$ [s]')
           84     ax2.set_ylabel(\
           85             'New heavily loaded contacts')
           86             #'Heavily loaded contacts $||\\boldsymbol{{\sigma}}_c|| \\geq \\sigma_0\\times 4$ Pa')
           87 
           88     ax3 = plt.subplot(3,1,3)
           89     #ax3.semilogy(t, coordinationnumber)
           90     ax3.plot(t, coordinationnumber, 'k')
           91     ax3.set_xlabel('Time $t$ [s]')
           92     ax3.set_ylabel('Coordination number $z$')
           93 
           94 
           95 
           96     plt.savefig(sid + '-nloaded.' + outformat)