TP3.py 3.67 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
import os,sys

from pylmgc90.chipy import *

##### TP3 : apply pressure to particles  #####

# TODO 1 : modify the function 'particles_with_force'
#          in 'applyPressure' to really add a force
#          to the particles.
# input : number of layers, time step and confinment pressure
# output : a numpy array with 0 on particles without force
#          and 1 on particles with a force
#
# needed functions:
# - RBDY3_PutBodyVector on field 'Fext_'

def applyPressure(dt, nb_layers,pressure):

  has_p = np.zeros(RBDY3_GetNbRBDY3(),dtype=np.double)
 
  nbs = SPHER_GetNbSPHER()
  s2b = SPHER_GetPtrSPHER2RBDY3()
  x = np.empty([nbs,4],dtype=np.double)
  for i in xrange(nbs):
    x[i,:3] = RBDY3_GetBodyVector('Coorb',int(s2b[i,0]))[:3]
    x[i,3]  = SPHER_GetContactorRadius(i+1)

  r  = np.empty([x.shape[0]])
  for i in xrange(x.shape[0]):
    r[i]  = np.linalg.norm(x[i,:2])

  h0 = np.min(x[:,2]-x[:,3])
  h  = np.max(x[:,2]+x[:,3]) - h0
  dh = h/nb_layers

  layer = []
  for i in xrange(nb_layers):
    layer.append([])
  for i in xrange(s2b.shape[0]):
    layer[int((x[i,2]-h0)/dh)].append(i)

  # Add impulse on each sphere in the crown.
  # The impulse value is a force times a time.
  f_ext = np.zeros([6])

  for k in layer:
    k = np.array(k)
    r_out = np.max( r[k] + x[k,3] )
    dr    = 2.*np.max( x[k,3] )
    sic   = np.nonzero( r[k] > r_out - dr )
    # Compute the equivalent force applied by
    # the confinment pressure on the surface of
    # the cylinder of a layer
    # F = P * 2 * pi * Radius * H
    for i in sic[0]:
      has_p[k[i]] = 1.
      # compute the force to apply to each sphere
      # F_i = radius_i / sum(radius_i)
      # F_i_x = - F * X_x/norm(X)
      # F_i_y = - F * X_y/norm(X)
      # add F_i times time step to external forces

  return has_p

checkDirectories()

dt = 1.e-1
nb_steps = 1000

theta = 0.5

freq_detect = 1
tol = 0.166e-4
relax = 1.0
norm = 'Quad '
gs_it1 = 33
gs_it2 = 100
78
storage='Stored_Delassus_Loop          '
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
nlgs_SetWithQuickScramble()

freq_write   = 10
freq_display = 10
ref_radius   = 1.

confinent_p = 5.e6
nb_layers   = 10

SetDimension(3)

TimeEvolution_SetTimeStep(dt)
Integrator_InitTheta(theta)

ReadBehaviours()
ReadBodies()
LoadBehaviours()
ReadIniDof()
ReadDrivenDof()

LoadTactors()

ReadIniVlocRloc()

WriteBodies()
WriteBehaviours()
WriteDrivenDof()

OpenDisplayFiles()
OpenPostproFiles()

ComputeMass()

nbs = SPHER_GetNbSPHER()
s2b = SPHER_GetPtrSPHER2RBDY3()
x = np.empty([nbs,4],dtype=np.double)
for i in xrange(nbs):
  x[i,:3] = RBDY3_GetBodyVector('Coorb',int(s2b[i,0]))[:3]
  x[i,3]  = SPHER_GetContactorRadius(i+1)

h0 = np.min(x[:,2]-x[:,3])
h  = np.max(x[:,2]+x[:,3]) - h0
dh = h / nb_layers

slayerId = np.empty([nbs])
for i in xrange(nbs):
  slayerId[i] = int((x[i,2]-h0)/dh) + 1

layerID = {'SPHER':slayerId}
layerID['PLANx'] = np.zeros(PLANx_GetNbPLANx())

has_p = applyPressure(dt,nb_layers,confinent_p)

WriteDisplayFiles(freq_display,ref_radius,has_p=('rigid', has_p),layerID=('tactor',layerID))

for k in xrange(nb_steps):
  #
  IncrementStep()
  
  ComputeFext()
  # TODO 3:
  #
  # Call your function to get the returned array
  # and add the computed impulse
  
  ComputeBulk()
  ComputeFreeVelocity()
  
  SelectProxTactors(freq_detect)
  
  RecupRloc()
  ExSolver(storage,norm,tol,relax,gs_it1,gs_it2)
  StockRloc()
  
  ComputeDof()
  UpdateStep()
  
  WriteLastDof()
  WriteLastVlocRloc()

  WriteDisplayFiles(freq_display,ref_radius,has_p=('rigid', has_p),layerID=('tactor',layerID))

  # TODO 4:
  #
  # Break out of the time loop if
  # the maximum velocity of the sphere is
  # less than a chosen value
  #
  # needed functions:
  # - RBDY3_GetBodyVector for field 'V____'
  # - np.linalg.norm

CloseDisplayFiles()
ClosePostproFiles()