scontact.py 18.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 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/>.

22
23
24
25
26
27
28
29
30
"""Model for Immersed Granular Flow -- Particles user interface

    Contact: jonathan.lambrechts@uclouvain.be
    Webpage: www.migflow.be

    MigFlow computes immersed granular flows using an unresolved FEM-DEM model. 
    The discrete phase is computed by updating iteratively the particle velocities until a set of velocities respecting the non-interpenetration contact law for the next time step is found
"""

31
32
from __future__ import division
from . import mesh
Nathan Coppin's avatar
Nathan Coppin committed
33
from . import VTK
34
import shutil
Nathan Coppin's avatar
Nathan Coppin committed
35
import numpy as np
36
37
38
39
40
41
42
43
import os
import sys
from ctypes import *
import signal

signal.signal(signal.SIGINT, signal.SIG_DFL)

dir_path = os.path.dirname(os.path.realpath(__file__))
Nathan Coppin's avatar
Nathan Coppin committed
44
lib2 = np.ctypeslib.load_library("libscontact2",dir_path)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
45
lib2f = np.ctypeslib.load_library("libscontact2f",dir_path)
Nathan Coppin's avatar
Nathan Coppin committed
46
lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
47
48
49

is_64bits = sys.maxsize > 2**32

Nathan Coppin's avatar
Nathan Coppin committed
50
51
52
53
54
55
def _np2c(a,dtype=np.float64) :
    tmp = np.require(a,dtype,"C")
    r = c_void_p(tmp.ctypes.data)
    r.tmp = tmp
    return r

56
class ParticleProblem :
57
    """Create the numerical structure containing all the physical particles that appear in the problem"""
58
59
60
61
62
63
64

    def _get_matrix(self, fName, ncol) :
        f = getattr(self._lib,"particleProblem"+fName)
        f.restype = c_void_p
        ptr = f(self._p)
        size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//8
        buf = (size*c_double).from_address(ptr)
Nathan Coppin's avatar
Nathan Coppin committed
65
        return np.ctypeslib.as_array(buf).reshape([-1,ncol])
66
67
68
69
70
71
72
    
    def _get_idx(self, fName) :
        f = getattr(self._lib,"particleProblem"+fName)
        f.restype = c_void_p
        ptr = f(self._p)
        size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//4
        buf = (size*c_int).from_address(ptr)
Nathan Coppin's avatar
Nathan Coppin committed
73
        return np.ctypeslib.as_array(buf)
74

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
75
    def __init__(self, dim, friction_enabled=False) :
76
77
78
79
80
81
82
83
84
        """Build the particles contanier structure.

        Keyword arguments:
        dim -- dimension of the problem

        Raises:
        ValueError -- if dimension differs from 2 or 3
        NameError -- if C builder cannot be found
        """
85
        self._dim = dim
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
86
        self._friction_enabled = friction_enabled
87
        if dim == 2 :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
88
            self._lib = lib2f if friction_enabled else lib2
89
90
91
92
93
94
95
96
97
98
99
            self._coord_type = c_double*2
        elif dim == 3 :
            self._lib = lib3
            self._coord_type = c_double*3
        else :
            raise ValueError("Dimension should be 2 or 3.")
        self._lib.particleProblemNew.restype = c_void_p
        self._p = c_void_p(self._lib.particleProblemNew())
        if self._p == None :
            raise NameError("Cannot create particle problem.")

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
100
101
102
103
104
    def __del__(self):
        """Delete the particle solver structure."""
        if(self._p is not None) :
            self._lib.particleProblemDelete(self._p)

105
    def volume(self):
106
        """Return the list of particle volumes."""
107
        if self._dim == 2 :
Nathan Coppin's avatar
Nathan Coppin committed
108
            return np.pi * self._get_matrix("Particle", 2)[:, [0]]**2
109
        else :
Nathan Coppin's avatar
Nathan Coppin committed
110
            return 4./3.*np.pi * self._get_matrix("Particle", 2)[:, [0]]**3
111
112

    def r(self) :
113
        """Return the list of particle radii."""
Nathan Coppin's avatar
Nathan Coppin committed
114
        return self._get_matrix("Particle", 2)[:, 0][:,None]
115
116

    def mass(self):
117
        """Return the list of particle masses."""
Nathan Coppin's avatar
Nathan Coppin committed
118
        return self._get_matrix("Particle", 2)[:, 1][:,None]
119
120

    def position(self):
121
        """Return the list of particle centre position vectors."""
122
123
124
        return self._get_matrix("Position", self._dim)

    def velocity(self):
125
        """Return the list of particle velocity vectors."""
126
127
        return self._get_matrix("Velocity", self._dim)

Nathan Coppin's avatar
Nathan Coppin committed
128
    def omega(self):
129
130
131
        """Return the list of particle angular velocity.
           Only in 2D when friction is enabled.
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
132
        assert self._friction_enabled
Nathan Coppin's avatar
Nathan Coppin committed
133
134
135
136
137
        return self._get_matrix("Omega", 1)

    def disk_tag(self):
        return self._get_idx("DiskTag")

138
139
140
141
    def disk_material(self):
        return self._get_idx("DiskMaterial")


142
    def segments(self):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
143
        return self._get_matrix("Segment",3*self._dim+(1 if self._friction_enabled else 0))
144
145
146
147

    def segment_tag(self):
        return self._get_idx("SegmentTag")

148
149
150
    def segment_material(self):
        return self._get_idx("SegmentMaterial")

151
    def triangles(self):
Nathan Coppin's avatar
Nathan Coppin committed
152
153
154
155
156
157
158
159
160
161
        if self._dim == 3 :
            return self._get_matrix("Triangle",12)
        else :
            return np.ndarray([0,12])

    def triangle_tag(self):
        if self._dim == 3 :
            return self._get_idx("TriangleTag")
        else :
            return np.array([],np.int32)
162

163
164
165
166
167
168
    def triangle_material(self):
        if self._dim == 3 :
            return self._get_idx("TriangleMaterial")
        else :
            return np.array([],np.int32)

169
170
171
    def set_boundary_velocity(self, tag, v) :
        self._lib.particleProblemGetTagId.restype = c_size_t
        tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
Nathan Coppin's avatar
Nathan Coppin committed
172
        disks = self._get_matrix("Disk",2*self._dim+1)
173
        diskTag = self._get_idx("DiskTag")
Nathan Coppin's avatar
Nathan Coppin committed
174
175
        disks[diskTag == tag, self._dim:2*self._dim] = v
        segments = self._get_matrix("Segment",3*self._dim)
176
        segmentTag = self._get_idx("SegmentTag")
Nathan Coppin's avatar
Nathan Coppin committed
177
        segments[segmentTag == tag, 2*self._dim:3*self._dim] = v
178
179
180
181
182
        if self._dim == 3 :
            triangles = self._get_matrix("Triangle",12)
            triangleTag = self._get_idx("TriangleTag")
            triangles[triangleTag == tag, 9:12] = v

Nathan Coppin's avatar
Nathan Coppin committed
183
    def set_tangent_boundary_velocity(self, tag, v):
184
185
186
187
188
189
190
        """ Sets the tangent velocity of the frictional boundary.
            Only in 2D when friction is enabled.

        Keyword arguments:
        tag -- Name of the boundary
        v -- velocity
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
191
        assert self._friction_enabled
Nathan Coppin's avatar
Nathan Coppin committed
192
        self._lib.particleProblemSetBoundaryTangentVelocity(self._p,tag.encode(),c_double(v))
193
    
194
195
196
197
198
199
    def set_use_queue(self, use_queue = True):
       """ Sets the use of the queue in the resolution of the contacts."""
       if use_queue:
         self._lib.particleProblemSetUseQueue(self._p,c_int(1))
       else:
         self._lib.particleProblemSetUseQueue(self._p,c_int(0))
Nathan Coppin's avatar
Nathan Coppin committed
200

201
    def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
202
203
204
205
206
207
208
209
        """ Sets the friction coefficient between two materials.
            Only in 2D when friction is enabled.

        Keyword arguments:
        mu -- value of the friction coefficient
        mat0 -- first material
        mat1 -- second material
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
210
        assert self._friction_enabled
211
        self._lib.particleProblemSetFrictionCoefficient(self._p, c_double(mu), mat0.encode(),mat1.encode())
Nathan Coppin's avatar
Nathan Coppin committed
212

213
    def iterate(self, dt, forces,tol=1e-8) :
214
215
216
217
218
219
220
        """Compute iteratively the set of velocities that respect the non-interpenetration constrain.

        Keyword arguments:
        dt -- numerical time step
        forces -- list of vectors containing the forces applied on the particles
        tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
        """ 
221
        self._get_matrix("ExternalForces",self._dim)[:] = forces
222
        self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
223
224

    def write_vtk(self, odir, i, t) :
225
        """Write output files for post-visualization."""
226

Nathan Coppin's avatar
Nathan Coppin committed
227
        v = self.velocity()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
228
229
230
231
        if self._friction_enabled :
            omega = self.omega()
        else :
            omega = np.zeros((v.shape[0],1))
Nathan Coppin's avatar
Nathan Coppin committed
232
        x = self.position()
233
        material = self._get_idx("ParticleMaterial").reshape(-1,1)
Nathan Coppin's avatar
Nathan Coppin committed
234
235
236
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
            x = np.insert(x,self._dim,0,axis=1)
237
238
239
240
241
        data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Material",material)]
        nmat = self._lib.particleProblemNMaterial(self._p)
        self._lib.particleProblemGetMaterialTagName.restype = c_char_p
        tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
        VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags))])
Nathan Coppin's avatar
Nathan Coppin committed
242

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
243
        xdisk = self._get_matrix("Disk",2*self._dim+(2 if self._friction_enabled else 1))[:,:self._dim]
Nathan Coppin's avatar
Nathan Coppin committed
244
245
246
247
248
249
250
251
252
253
254
        nd = xdisk.shape[0]
        xsegment = self.segments()[:,:self._dim*2].reshape([-1,self._dim])
        ns = xsegment.shape[0]/2
        xtriangles = self.triangles()[:,:self._dim*3].reshape([-1,self._dim])
        nt = xtriangles.shape[0]/3
        x = np.vstack([xdisk,xsegment,xtriangles])
        if self._dim == 2 :
            x = np.insert(x,self._dim,0,axis=1)
        connectivity = np.arange(0,x.shape[0],dtype=np.int32)
        types = np.repeat(np.array([1,3,5],dtype=np.int32),[nd,ns,nt])
        offsets = np.cumsum(np.repeat([1,2,3],[nd,ns,nt]),dtype=np.int32)
255
256
257
258
259
260
        tags = np.array(np.hstack([self.disk_tag(),self.segment_tag(),self.triangle_tag()]))
        materials = np.array(np.hstack([self.disk_material(),self.segment_material(),self.triangle_material()]))
        ntagname = self._lib.particleProblemNTag(self._p)
        self._lib.particleProblemGetTagName.restype = c_char_p
        tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)])
        VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1]))],field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
Nathan Coppin's avatar
Nathan Coppin committed
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276

        self._lib.particleProblemContactN.restype = c_size_t
        fdata = []
        for (ctype,name) in ((0,b"particle_particle"),(1,b"particle_disk"),(2,b"particle_segment"),(3,b"particle_triangle")) :
            n = self._lib.particleProblemContactN(self._p,c_int(ctype))
            v = np.ndarray((n,2),dtype=np.float64,order="C")
            o = np.ndarray((n,2),dtype=np.uint64,order="C")
            self._lib.particleProblemContact(self._p,c_int(ctype),c_void_p(o.ctypes.data),c_void_p(v.ctypes.data))
            fdata.append((name,v[:,[0]]))
            fdata.append((name+b"_t",v[:,[1]]))
            fdata.append((name+b"_idx",o))
        x = np.array([[0.,0.,0.]])
        VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")

        VTK.write_multipart(odir+"/particle_problem",["particles_%05d.vtp"%i,"boundaries_%05d.vtu"%i,"contacts_%05d.vtp"%i],i,t)
        
277
    def read_vtk(self,dirname,i) :
278
        """Read an existing output file to set the particle data."""
Nathan Coppin's avatar
Nathan Coppin committed
279
        x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
280
281
282
283
284
        matnames =  VTK.string_array_decode(fdata["MaterialNames"])
        matmap = {}
        for j,n in enumerate(matnames) :
            matmap[j] = self._lib.particleProblemGetMaterialTagId(self._p,n)
        partmat =np.vectorize(matmap.get)(d["Material"])
Nathan Coppin's avatar
Nathan Coppin committed
285
        self._lib.particleProblemClearParticles(self._p)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
286

Nathan Coppin's avatar
Nathan Coppin committed
287
        self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
288
                _np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(partmat,dtype=np.int32))
Nathan Coppin's avatar
Nathan Coppin committed
289
290

        x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
291
292
293
294
        tagnames =  VTK.string_array_decode(fdata["TagNames"])
        tagmap = {}
        for j,n in enumerate(tagnames) :
            tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
Nathan Coppin's avatar
Nathan Coppin committed
295
296
297
        self._lib.particleProblemClearBoundaries(self._p)
        offsets = np.hstack([[0],cells["offsets"]])
        el = cells["connectivity"]
298
299
        tags = np.vectorize(tagmap.get)(cdata["Tag"])
        materials = cdata["Material"]
Nathan Coppin's avatar
Nathan Coppin committed
300
301
        for idx in np.where(cells["types"] == 1)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
302
            self._lib.particleProblemAddBoundaryDiskTagId(self._p, self._coord_type(*x0), c_double(0.), c_int(tags[idx,0]), c_int(materials[idx,0]))
Nathan Coppin's avatar
Nathan Coppin committed
303
304
305
        for idx in np.where(cells["types"] == 3)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
            x1 = x[el[offsets[idx]+1],:self._dim]
306
            self._lib.particleProblemAddBoundarySegmentTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(tags[idx,0]), c_int(materials[idx,0]))
Nathan Coppin's avatar
Nathan Coppin committed
307
308
309
310
        for idx in np.where(cells["types"] == 5)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
            x1 = x[el[offsets[idx]+1],:self._dim]
            x2 = x[el[offsets[idx]+2],:self._dim]
311
            self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(tags[idx,0]), c_int(materials[idx,0]))
Nathan Coppin's avatar
Nathan Coppin committed
312
313
314
315
316
317
318
319
        _,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
        ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
        v = np.vstack([fdata[k] for k in ks])
        vt = np.vstack([fdata[k+"_t"] for k in ks])
        v = np.hstack([v,vt])
        o = np.vstack([fdata[k+"_idx"] for k in ks])
        t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
        self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(t,np.int32))
320

321
    def add_boundary_disk(self, x0, r, tag, material="default") :
322
        self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
323
        return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode(),material.encode())
324

325
    def add_boundary_segment(self, x0, x1, tag, material="default") :
326
        self._lib.particleProblemAddBoundarySegment.restype = c_size_t
327
        return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode(),material.encode())
328

329
    def add_boundary_triangle(self, x0, x1, x2, tag, material="default") :
330
        if self._dim != 3 :
331
332
            raise NameError("Triangle boundaries only available in 3D")
        self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t
333
        return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), tag.encode(),material.encode())
334

335
    def add_particle(self, x, r, m, tag="default"):
336
337
338
339
340
341
        """Add a particle in the particle container.

        Keyword arguments:
        x -- tuple to set the centre particle position
        r -- particle radius
        m -- particle mass
Nathan Coppin's avatar
Nathan Coppin committed
342
        tag -- particle tag
343
        """
344
        self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), tag.encode())
345

346
    def load_msh_boundaries(self, filename, tags, shift=None,material="default") :
347
348
349
350
351
352
        """Load the numerical domain and set the physical boundaries the particles cannot go through.
        
        Keyword arguments:
        filename -- name of the msh file
        tags -- tags of physical boundary defined in the msh file
        shift -- optional argument to shift the numerical domain
353
        material -- material of the boundary
354
        """
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
        if shift == None :
            shift = [0]*self._dim
        m = mesh.Mesh(filename)
        self._addv = set()
        if self._dim == 2 :
            pnum = [m.getPhysicalNumber(1, i) for i in tags]
            for ent in m.entities :
                if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
                    continue
                tag = 0 if not ent.physicals else ent.physicals[0]
                stag = m.getPhysicalName(ent.dimension, tag) or str(tag)
                for i, el in enumerate(ent.elements) :
                    for v in el.vertices :
                        if v[3] in self._addv :
                            continue
                        self._addv.add(v[3])
371
                        self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag,material)
372
373
                    v0 = el.vertices[1]
                    v1 = el.vertices[0]
374
                    self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag,material)
375
376
377
378
379
380
381
382
383
384
385
386
387
        else :
            self._adde = set()
            pnum = [m.getPhysicalNumber(2, i) for i in tags]
            for ent in m.entities :
                if ent.dimension !=2 or not [i for i in pnum if i in ent.physicals]:
                    continue
                tag = 0 if not ent.physicals else ent.physicals[0]
                stag = m.getPhysicalName(ent.dimension, tag) or str(tag)
                for i, el in enumerate(ent.elements) :
                    for j in range(3):
                        v0 = el.vertices[j]
                        if not v0[3] in self._addv :
                            self._addv.add(v0[3])
388
                            self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag,material)
389
390
391
392
393
394
                        v1 = el.vertices[(j+1)%3]
                        if (not (v0[3],v1[3]) in self._adde) and (not (v1[3],v0[3]) in self._adde) :
                            self._adde.add((v0[3],v1[3]))
                            self.add_boundary_segment(
                                (v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),
                                (v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
395
                                stag,material)
396
397
398
399
400
401
402
403
                            
                    v0 = el.vertices[1]
                    v1 = el.vertices[0]
                    v2 = el.vertices[2]
                    self.add_boundary_triangle(
                        (v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),
                        (v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
                        (v2[0]+shift[0],v2[1]+shift[1],v2[2]+shift[2]),
404
                          stag,material)