lmgc90Interface.py 16.5 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
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not, 
# see <http://www.gnu.org/licenses/>.

from pylmgc90 import chipy
import numpy as np
import os
import math
Frédéric Dubois's avatar
Frédéric Dubois committed
26
import sys
27

Frédéric Dubois's avatar
Frédéric Dubois committed
28
class ParticleProblem(object):
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
  """
  A class definition holding data of lmgc90
  in a compatible form to use coupling with
  gmsh
  
  Attributs are :

  - number of objects (disks|sphere) 
  - volume of objects (disks|sphere)
  - reference map on objects to models (diskx2rbdy2 | spher2rbdy3)

  The methods are :
  
  - iterate()  : do one time iteration
  - position() : give back position of all spher
Matthieu Constant's avatar
Matthieu Constant committed
44
  - velocity() : give back velocity of all spher
45
46
47
48
49
50
51
52
  - externalForces() : update external forces on particles
  
  - getMeanRadius() : return the mean radius of particles
  - getVolume() : give back the volume of all spher 

  - writeHeader() : write header file for plot
  """

53
  def __init__(self,dim,period=0,xp=0,yp=0,restart=0) :
54
      """
Frédéric Dubois's avatar
Frédéric Dubois committed
55
      Initialize lmgc90
56
57
58
59
60
61
62
63
64
      """

      self._dim=dim
      
      chipy.Initialize()
      
      chipy.checkDirectories()
      
      chipy.SetDimension(dim)
65
      chipy.TimeEvolution_SetTimeStep(0.)      
66
67
68
69
70
71
72
73
74
75
76
77
      chipy.Integrator_InitTheta(0.5)

      chipy.utilities_logMes('READ BEHAVIOURS')
      chipy.ReadBehaviours()
      
      chipy.utilities_logMes('READ BODIES')
      chipy.ReadBodies()
      
      chipy.utilities_logMes('LOAD BEHAVIOURS')
      chipy.LoadBehaviours()
      
      chipy.utilities_logMes('READ INI DOF')
78
      chipy.ReadIniDof(restart)
79
80
81
82
83
84
85
86
      
      chipy.utilities_logMes('READ DRIVEN DOF')
      chipy.ReadDrivenDof()
      
      chipy.utilities_logMes('LOAD TACTORS')
      chipy.LoadTactors()
      
      chipy.utilities_logMes('READ INI Vloc Rloc')
87
      chipy.ReadIniVlocRloc(restart)
88

Frédéric Dubois's avatar
Frédéric Dubois committed
89
90
91
      if (period):
        chipy.SetPeriodicCondition(xperiod=xp,yperiod=yp)
      
92
93
94
      chipy.PRPRx_UseCpCundallDetection(30)

      chipy.OpenPostproFiles()
95
      chipy.OpenDisplayFiles(restart+1)
96
      
97
      chipy.WriteDisplayFiles()
98
99
100
101
102
103
104

      chipy.ComputeMass()

      if dim == 2 :

        self._nbDisk   = chipy.DISKx_GetNbDISKx()
        self._d2bMap   = chipy.DISKx_GetPtrDISKx2RBDY2()
105
        self._nbPoly = chipy.POLYG_GetNbPOLYG()
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
        self._p2bMap = chipy.POLYG_GetPtrPOLYG2RBDY2()
        self._position = np.zeros([self._nbDisk,3], 'd')
        self._velocity = np.zeros([self._nbDisk,3], 'd')
        self._externalF= np.zeros([self._nbDisk,3], 'd')
        self._volume   = np.zeros([self._nbDisk,1], 'd')
        self._mass   = np.zeros([self._nbDisk,1], 'd')
        self._r   = np.zeros([self._nbDisk,1], 'd')
        for i in range(self._nbDisk):
          self._r[i] = chipy.DISKx_GetContactorRadius(i+1)
          self._mass[i] = chipy.RBDY2_GetBodyMass(int(self._d2bMap[i,0]))
        self._volume = np.pi * self._r**2 
        self._tag2bnd = {}
        self._tag2id = {}
        for i in range(chipy.RBDY2_GetNbRBDY2()):
            for t in range(chipy.RBDY2_GetNbContactor(i+1)):
                tag = chipy.RBDY2_GetContactorColor(i+1,1).rstrip('_')
                self._tag2id.setdefault(tag, []).append(i)
              
124
        self._refRadius=np.amin(self._r)
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
      
      elif dim == 3:  
      
        self._nbSpher   = chipy.SPHER_GetNbSPHER()
        self._d2bMap   = chipy.SPHER_GetPtrSPHER2BDYTY()
        self._position = np.zeros([self._nbSpher,3], 'd')
        self._velocity = np.zeros([self._nbSpher,6], 'd')
        self._externalF= np.zeros([self._nbSpher,6], 'd')
        self._volume   = np.zeros([self._nbSpher,1], 'd')
        self._mass   = np.zeros([self._nbSpher,1], 'd')
        self._r   = np.zeros([self._nbSpher,1], 'd')      
        for i in range(self._nbSpher):
          self._r[i] = chipy.SPHER_GetContactorRadius(i+1)
          self._volume[i] = np.pi * chipy.SPHER_GetContactorRadius(i+1)**3 * 4./3
          self._mass[i] = chipy.RBDY3_GetMass(int(self._d2bMap[i,0]))

        self._refRadius=np.amin(self._r)
        
        self._tag2bnd = {}
        self._tag2id = {}
        for i in range(chipy.RBDY3_GetNbRBDY3()):
            for t in range(chipy.RBDY3_GetNbContactor(i+1)):
                tag = chipy.RBDY3_GetContactorColor(i+1,1).rstrip('_')
                self._tag2id.setdefault(tag, []).append(i)

Frédéric Dubois's avatar
Frédéric Dubois committed
150
151
152
      else :
        chipy.utilities_logMes('dim must be 2 or 3')
        sys.exit()
153
                
154
155
156
157
158
159
160
161
162
      self._freq_detect   = 1
      self._solver_params = { 'type'  :'Stored_Delassus_Loops         ',
                              'norm'  : 'Quad ',
                              'conv'  : 1e-5,
                              'relax' : 1.,
                              'gsit1' : 200,
                              'gsit2' : 20
                            }
   
163
  def setBoundaryVelocity(self, tag, v) :
164
165
      self._tag2bnd[tag] = v

166
167
168
  def write_vtk(self, filemame, t, i) :
    chipy.WriteOutDof(1)
    chipy.WriteOutVlocRloc(1)
169
170
171
    chipy.WritePostproFiles()
    chipy.WriteDisplayFiles(1, self._refRadius)

172
  def iterate(self, dt, forces, tol=None, gsit1=None, gsit2=None):
173
      """
Frédéric Dubois's avatar
Frédéric Dubois committed
174
      Do one step of a lmgc90 computation.
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
      """

      if tol is not None : self._solver_params['conv'] = tol          
      if gsit1 is not None : self._solver_params['gsit1'] = gsit1
      if gsit2 is not None : self._solver_params['gsit2'] = gsit2          

      
      chipy.TimeEvolution_SetTimeStep(dt)

      self._externalF[:,:self._dim] = forces

      chipy.IncrementStep()
      chipy.ComputeFext()
      if (self._dim == 3):
        for tag, v in self._tag2bnd.items() :
            for i in self._tag2id[tag] :
                for j, iv in enumerate(v):
                  chipy.RBDY3_SetVlocyDrivenDof(i+1,j+1,iv)
  

        for i in range(self._nbSpher):
          chipy.RBDY3_PutBodyVector('Fext_', int(self._d2bMap[i,0]), self._externalF[i,:])

      elif (self._dim == 2):
        for tag, v in self._tag2bnd.items() :
           for i in self._tag2id[tag] :
              for j, iv in enumerate(v):
                 chipy.RBDY2_SetVlocyDrivenDof(i+1,j+1,iv)
  

        for i in range(self._nbDisk):
          chipy.RBDY2_PutBodyVector('Fext_', int(self._d2bMap[i,0]), self._externalF[i,:])        

      chipy.ComputeBulk()
        
      chipy.ComputeFreeVelocity()
  
      chipy.SelectProxTactors(self._freq_detect)
  
      chipy.RecupRloc()
      chipy.ExSolver(**self._solver_params)
      chipy.StockRloc()
  
      chipy.ComputeDof()
  
      chipy.UpdateStep()

  def volume(self):
      """ Return the volume of the  particles """
      return self._volume

  def mass(self):
      """ Return the mass of the  particles """
      return self._mass

  def r(self):
      """ Return the radius of the  particles """
      return self._r

  def position(self):
      """ Get current position of contactors and return it """

      if self._dim == 3:
        for i in range(self._nbSpher):
          self._position[i,:] = chipy.SPHER_GetContactorCoor(i+1)[:3]
        return self._position
      
      elif  self._dim == 2:
        for i in range(self._nbDisk):
          self._position[i,:] = chipy.RBDY2_GetBodyVector('Coor_',int(self._d2bMap[i,0]))
        self._position[i,2] = 0.
        return self._position[:,:2]
      
  def velocity(self):
      """ Get current velocity of body of contactor and return it
      Beware : it should not work very well with clusters !
      """
      if self._dim == 3:
        for i in range(self._nbSpher):
          self._velocity[i,:] = chipy.RBDY3_GetBodyVector('V____',int(self._d2bMap[i,0]))
        return self._velocity[:,:3]

      elif self._dim==2:
        for i in range(self._nbDisk):
Matthieu Constant's avatar
Matthieu Constant committed
259
260
261
          self._velocity[i,:] = chipy.RBDY2_GetBodyVector('V____',int(self._d2bMap[i,0]))
        self._velocity[i,2] = 0.
        return self._velocity[:,:2]
262
263
264
265
266
267
268
269
270
271
272
273
274
275
        
      
  def externalForces(self):
      """ Get an external forces array
      """
      return self._externalF[:,:self._dim]

  def __del__(self):
      """
      """
      chipy.ClosePostproFiles()
      chipy.CloseDisplayFiles()
      chipy.Finalize()

Frédéric Dubois's avatar
Frédéric Dubois committed
276
def scontactTolmgc90(dirname, dim, it=0, fric = 0., assume_box=False):
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
    from migflow import scontact
    from pylmgc90 import pre
    datbox_path = 'DATBOX'
    if not os.path.isdir(datbox_path):
      os.mkdir(datbox_path)

    sc = scontact.ParticleProblem(dim)
    sc.read_vtk(dirname,it)

    x = sc.position()
    v = sc.velocity()
    v[:] = 0.

    r = sc.r()
    vol = sc.volume()

    rho_particle = np.mean(sc.mass()/vol)

    avs = pre.avatars()
    mat = pre.materials()
    mod = pre.models()
    svs = pre.see_tables()
    tac = pre.tact_behavs()


    if dim == 3 :
    
      mater = pre.material(name='STONE',materialType='RIGID',density=rho_particle)
      model = pre.model(name='rigid', physics='MECAx', element='Rxx3D',dimension=3)
      mat.addMaterial(mater)
      mod.addModel(model)

      for i in range(r.size) :
         P = pre.rigidSphere( r=r[i], center=x[i], model=model, material=mater, color='INxxx')
         P.imposeInitValue(component=[1,2,3],value=list(v[i,:]))
         avs.addAvatar(P)

      clb_fric = pre.tact_behav(name='iqsc0',law='IQS_CLB',fric=fric)
      tac     += clb_fric

      svs += pre.see_table(CorpsCandidat   ="RBDY3", candidat   ="SPHER", colorCandidat   ='INxxx',
                           CorpsAntagoniste="RBDY3", antagoniste="SPHER", colorAntagoniste='INxxx',
                           behav=clb_fric, alert=r.min())
       
      ep=np.amin(r)   
      tri = sc.triangles()

      if (assume_box) :
        xmin=min(np.amin(tri[:,0]),np.amin(tri[:,3]), np.amin(tri[:,6]))
        xmax=min(np.amax(tri[:,0]),np.amax(tri[:,3]), np.amax(tri[:,6]))       
        ymin=min(np.amin(tri[:,1]),np.amin(tri[:,4]), np.amin(tri[:,7]))
        ymax=min(np.amax(tri[:,1]),np.amax(tri[:,4]), np.amax(tri[:,7]))       
        zmin=min(np.amin(tri[:,2]),np.amin(tri[:,5]), np.amin(tri[:,8]))
        zmax=min(np.amax(tri[:,3]),np.amax(tri[:,5]), np.amax(tri[:,8]))

        lx = xmax - xmin
        ly = ymax - ymin
        lz = zmax - zmin

        pxmin = pre.rigidPlan(axe1=lz*0.505, axe2=ly*0.505, axe3=ep, \
                              center=[xmin-ep,ymin+0.5*ly,zmin+0.5*lz],
                              model=model, material=mater, color='WALLx')
        pxmin.rotate(description='axis',center=pxmin.nodes[1].coor,axis=[0.,1.,0.],alpha=math.pi*0.5)
        pxmin.imposeDrivenDof(component=[1,2,3,4,5,6],dofty='vlocy')
        avs+=pxmin
        pxmax = pre.rigidPlan(axe1=lz*0.505, axe2=ly*0.505, axe3=ep, \
                              center=[xmax+ep,ymin+0.5*ly,zmin+0.5*lz],
                              model=model, material=mater, color='WALLx')
        pxmax.rotate(description='axis',center=pxmax.nodes[1].coor,axis=[0.,1.,0.],alpha=math.pi*0.5)
        pxmax.imposeDrivenDof(component=[1,2,3,4,5,6],dofty='vlocy')      
        avs+=pxmax

        pymin = pre.rigidPlan(axe1=lx*0.505, axe2=lz*0.505, axe3=ep, \
                              center=[xmin+0.5*lx,ymin-ep,zmin+0.5*lz],
                              model=model, material=mater, color='WALLx')
        pymin.rotate(description='axis',center=pymin.nodes[1].coor,axis=[1.,0.,0.],alpha=math.pi*0.5)
        pymin.imposeDrivenDof(component=[1,2,3,4,5,6],dofty='vlocy')      
        avs+=pymin
        
        pymax = pre.rigidPlan(axe1=lx*0.505, axe2=lz*0.505, axe3=ep, \
                              center=[xmin+0.5*lx,ymax+ep,zmin+0.5*lz],
                              model=model, material=mater, color='WALLx')
        pymax.rotate(description='axis',center=pymax.nodes[1].coor,axis=[1.,0.,0.],alpha=math.pi*0.5)
        pymax.imposeDrivenDof(component=[1,2,3,4,5,6],dofty='vlocy')      
        avs+=pymax
      
363
        pzmin = pre.rigidPlan(axe1=lx*0.505, axe2=ly*0.505, axe3=ep, \
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
                             center=[xmin+0.5*lx,ymin+0.5*ly,zmin-ep],
                             model=model, material=mater, color='WALLx')
        pzmin.imposeDrivenDof(component=[1,2,3,4,5,6],dofty='vlocy')      
        avs+=pzmin
        
        pzmax = pre.rigidPlan(axe1=lx*0.505, axe2=ly*0.505, axe3=ep, \
                              center=[xmin+0.5*lx,ymin+0.5*ly,zmax+ep],
                              model=model, material=mater, color='WALLx')
        pzmax.imposeDrivenDof(component=[1,2,3,4,5,6],dofty='vlocy')      
        avs+=pzmax
      
      
        svs += pre.see_table(CorpsCandidat   ="RBDY3", candidat   ="SPHER", colorCandidat   ='INxxx',
                             CorpsAntagoniste="RBDY3", antagoniste="PLANx", colorAntagoniste='WALLx',
                             behav=clb_fric, alert=r.min())
      else :
        bndtags = set()
        for i in range(tri.shape[0]) :
          x0 = np.array(tri[i, 0:3])
          x1 = np.array(tri[i, 3:6])
          x2 = np.array(tri[i, 6:9])
          tag = 'xxxxx'
          tag += "_" * (5 - len(tag))
          bndtags.add(tag)
          t0 = (x0 - x1)/np.linalg.norm(x0 - x1)
          t1 = (x0 - x2)/np.linalg.norm(x0 - x2)
          n = np.cross(t0, t1)
          n *= ep / np.linalg.norm(n)

          vs = np.array([x0-n, x1 - n, x2 - n, x0 + n, x1 + n, x2 + n])
          fs = np.array([[0, 2, 1], [3, 4, 5], [0, 1, 3], [3, 1, 4], [1, 5, 4], [1, 2, 5], [0, 3, 2], [3, 5, 2]]) + 1
          av = pre.rigidPolyhedron(model,mater,np.zeros([3]),color=tag,generation_type="full",vertices=vs, faces=fs)
          av.imposeDrivenDof(component=[1,2,3,4,5,6],dofty='vlocy')
          avs.addAvatar(av)
  
        for tag in bndtags :
          svs += pre.see_table(CorpsCandidat   ="RBDY3", candidat   ="SPHER", colorCandidat   ='INxxx',
                               CorpsAntagoniste="RBDY3", antagoniste="POLYR", colorAntagoniste=tag,
                               behav=clb_fric, alert=r.min())

    elif dim==2 :
      
      mater = pre.material(name='STONE',materialType='RIGID',density=rho_particle)
      model = pre.model(name='rigid', physics='MECAx', element='Rxx2D',dimension=2)
      mat.addMaterial(mater)
      mod.addModel(model)

      for i in range(r.size) :
         P = pre.rigidDisk( r=r[i], center=x[i], model=model, material=mater, color='INxxx')
         P.imposeInitValue(component=[1,2],value=list(v[i,:]))
         avs.addAvatar(P)

416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
      seg = sc.segments()
      bndtags = set()
      for i in range(seg.shape[0]) :
        x0 = np.array(seg[i, 0:2])
        x1 = np.array(seg[i, 2:4])
        tag = 'xxxxx'
        tag += "_" * (5 - len(tag))
        bndtags.add(tag)
        t = (x0 - x1)/np.linalg.norm(x0 - x1)
        n = np.array([-t[1], t[0]]) * 1e-4
        vs = np.array([x0-2*n, x0, x1, x1 -2*n])
        av = pre.rigidPolygon(model,mater,np.zeros([2]),color=tag,generation_type="full",vertices=vs)
        av.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
        avs.addAvatar(av)
        
431
432
433
434
435
436
437
438
439
440
441
      clb_fric = pre.tact_behav(name='iqsc0',law='IQS_CLB',fric=fric)
      tac     += clb_fric

      svs += pre.see_table(CorpsCandidat   ="RBDY2", candidat   ="DISKx", colorCandidat   ='INxxx',
                           CorpsAntagoniste="RBDY2", antagoniste="DISKx", colorAntagoniste='INxxx',
                           behav=clb_fric, alert=r.min())

      for tag in bndtags :
        svs += pre.see_table(CorpsCandidat   ="RBDY2", candidat   ="DISKx", colorCandidat   ='INxxx',
                             CorpsAntagoniste="RBDY2", antagoniste="POLYG", colorAntagoniste=tag,
                             behav=clb_fric, alert=r.min())
442

443
444
445
446
447
448
449
450
451
452
453
454
455
          
    # file writting
    pre.writeBodies(avs,chemin=datbox_path)
    pre.writeDrvDof(avs,chemin=datbox_path)
    pre.writeDofIni(avs,chemin=datbox_path)
    pre.writeModels(mod,chemin=datbox_path)
    pre.writeBulkBehav(mat,chemin=datbox_path, gravy=[0.,0.,0.])
    pre.writeTactBehav(tac,svs,chemin=datbox_path)
    pre.writeVlocRlocIni(chemin=datbox_path)
    
    post = pre.postpro_commands()
    post.addCommand(pre.postpro_command(name='SOLVER INFORMATIONS', step=1))
    pre.writePostpro(commands=post, parts=avs, path=datbox_path)