scontact.py 17.6 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
129
    def omega(self):
        """Return the list of particle angular velocity."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
130
        assert self._friction_enabled
Nathan Coppin's avatar
Nathan Coppin committed
131
132
133
134
135
        return self._get_matrix("Omega", 1)

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

136
137
138
139
    def disk_material(self):
        return self._get_idx("DiskMaterial")


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

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

146
147
148
    def segment_material(self):
        return self._get_idx("SegmentMaterial")

149
    def triangles(self):
Nathan Coppin's avatar
Nathan Coppin committed
150
151
152
153
154
155
156
157
158
159
        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)
160

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

167
168
169
    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
170
        disks = self._get_matrix("Disk",2*self._dim+1)
171
        diskTag = self._get_idx("DiskTag")
Nathan Coppin's avatar
Nathan Coppin committed
172
173
        disks[diskTag == tag, self._dim:2*self._dim] = v
        segments = self._get_matrix("Segment",3*self._dim)
174
        segmentTag = self._get_idx("SegmentTag")
Nathan Coppin's avatar
Nathan Coppin committed
175
        segments[segmentTag == tag, 2*self._dim:3*self._dim] = v
176
177
178
179
180
        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
181
    def set_tangent_boundary_velocity(self, tag, v):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
182
        assert self._friction_enabled
Nathan Coppin's avatar
Nathan Coppin committed
183
        self._lib.particleProblemSetBoundaryTangentVelocity(self._p,tag.encode(),c_double(v))
184
185
    
    def set_friction_relaxation(self, relaxation):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
186
        assert self._friction_enabled
187
        self._lib.particleProblemSetFrictionRelaxation(self._p,c_double(relaxation))
Nathan Coppin's avatar
Nathan Coppin committed
188

189
    def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
190
        assert self._friction_enabled
191
        self._lib.particleProblemSetFrictionCoefficient(self._p, c_double(mu), mat0.encode(),mat1.encode())
Nathan Coppin's avatar
Nathan Coppin committed
192

193
    def iterate(self, dt, forces,tol=1e-8) :
194
195
196
197
198
199
200
        """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
        """ 
201
        self._get_matrix("ExternalForces",self._dim)[:] = forces
Nathan Coppin's avatar
Nathan Coppin committed
202
        self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
203
204

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

Nathan Coppin's avatar
Nathan Coppin committed
207
        v = self.velocity()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
208
209
210
211
        if self._friction_enabled :
            omega = self.omega()
        else :
            omega = np.zeros((v.shape[0],1))
Nathan Coppin's avatar
Nathan Coppin committed
212
        x = self.position()
213
        material = self._get_idx("ParticleMaterial").reshape(-1,1)
Nathan Coppin's avatar
Nathan Coppin committed
214
215
216
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
            x = np.insert(x,self._dim,0,axis=1)
217
218
219
220
221
        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
222

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
223
        xdisk = self._get_matrix("Disk",2*self._dim+(2 if self._friction_enabled else 1))[:,:self._dim]
Nathan Coppin's avatar
Nathan Coppin committed
224
225
226
227
228
229
230
231
232
233
234
        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)
235
236
237
238
239
240
        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
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256

        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)
        
257
    def read_vtk(self,dirname,i) :
258
        """Read an existing output file to set the particle data."""
Nathan Coppin's avatar
Nathan Coppin committed
259
        x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
260
261
262
263
264
        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
265
        self._lib.particleProblemClearParticles(self._p)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
266

Nathan Coppin's avatar
Nathan Coppin committed
267
        self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
268
                _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
269
270

        x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
271
272
273
274
        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
275
276
277
        self._lib.particleProblemClearBoundaries(self._p)
        offsets = np.hstack([[0],cells["offsets"]])
        el = cells["connectivity"]
278
279
        tags = np.vectorize(tagmap.get)(cdata["Tag"])
        materials = cdata["Material"]
Nathan Coppin's avatar
Nathan Coppin committed
280
281
        for idx in np.where(cells["types"] == 1)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
282
            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
283
284
285
        for idx in np.where(cells["types"] == 3)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
            x1 = x[el[offsets[idx]+1],:self._dim]
286
            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
287
288
289
290
        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]
291
            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
292
293
294
295
296
297
298
299
        _,_,_,_,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))
300

301
    def add_boundary_disk(self, x0, r, tag, material="default") :
302
        self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
303
        return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode(),material.encode())
304

305
    def add_boundary_segment(self, x0, x1, tag, material="default") :
306
        self._lib.particleProblemAddBoundarySegment.restype = c_size_t
307
        return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode(),material.encode())
308

309
    def add_boundary_triangle(self, x0, x1, x2, tag, material="default") :
310
        if self._dim != 3 :
311
312
            raise NameError("Triangle boundaries only available in 3D")
        self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t
313
        return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), tag.encode(),material.encode())
314

Nathan Coppin's avatar
Nathan Coppin committed
315
    def add_particle(self, x, r, m, tag):
316
317
318
319
320
321
        """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
322
        tag -- particle tag
323
        """
324
        self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), tag.encode())
325

326
    def load_msh_boundaries(self, filename, tags, shift=None,material="default") :
327
328
329
330
331
332
333
        """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
        """
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
        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])
350
                        self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag,material)
351
352
                    v0 = el.vertices[1]
                    v1 = el.vertices[0]
353
                    self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag,material)
354
355
356
357
358
359
360
361
362
363
364
365
366
        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])
367
                            self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag,material)
368
369
370
371
372
373
                        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]),
374
                                stag,material)
375
376
377
378
379
380
381
382
                            
                    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]),
383
                          stag,material)