scontact.py 15.3 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
45
lib2 = np.ctypeslib.load_library("libscontact2",dir_path)
lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
46
47
48

is_64bits = sys.maxsize > 2**32

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

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

    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
64
        return np.ctypeslib.as_array(buf).reshape([-1,ncol])
65
66
67
68
69
70
71
    
    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
72
        return np.ctypeslib.as_array(buf)
73
74

    def __init__(self, dim) :
75
76
77
78
79
80
81
82
83
        """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
        """
84
85
86
87
88
89
90
91
92
93
94
95
96
97
        self._dim = dim
        if dim == 2 :
            self._lib = lib2
            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
98
99
100
101
102
    def __del__(self):
        """Delete the particle solver structure."""
        if(self._p is not None) :
            self._lib.particleProblemDelete(self._p)

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

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

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

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

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

Nathan Coppin's avatar
Nathan Coppin committed
126
127
128
129
130
131
132
    def omega(self):
        """Return the list of particle angular velocity."""
        return self._get_matrix("Omega", 1)

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

133
    def segments(self):
Nathan Coppin's avatar
Nathan Coppin committed
134
        return self._get_matrix("Segment",3*self._dim+1)
135
136
137
138
139

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

    def triangles(self):
Nathan Coppin's avatar
Nathan Coppin committed
140
141
142
143
144
145
146
147
148
149
        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)
150
151
152
153

    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
154
        disks = self._get_matrix("Disk",2*self._dim+1)
155
        diskTag = self._get_idx("DiskTag")
Nathan Coppin's avatar
Nathan Coppin committed
156
157
        disks[diskTag == tag, self._dim:2*self._dim] = v
        segments = self._get_matrix("Segment",3*self._dim)
158
        segmentTag = self._get_idx("SegmentTag")
Nathan Coppin's avatar
Nathan Coppin committed
159
        segments[segmentTag == tag, 2*self._dim:3*self._dim] = v
160
161
162
163
164
        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
165
166
    def set_tangent_boundary_velocity(self, tag, v):
        self._lib.particleProblemSetBoundaryTangentVelocity(self._p,tag.encode(),c_double(v))
167
168
169
    
    def set_friction_relaxation(self, relaxation):
        self._lib.particleProblemSetFrictionRelaxation(self._p,c_double(relaxation))
Nathan Coppin's avatar
Nathan Coppin committed
170

171
    def add_friction(self, static) :
Nathan Coppin's avatar
Nathan Coppin committed
172
        arg_type = c_double*len(static)
173
        self._lib.particleProblemAddFriction(self._p, arg_type(*static),len(static))
Nathan Coppin's avatar
Nathan Coppin committed
174

175
    def iterate(self, dt, forces,tol=1e-8) :
176
177
178
179
180
181
182
        """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
        """ 
183
        self._get_matrix("ExternalForces",self._dim)[:] = forces
Nathan Coppin's avatar
Nathan Coppin committed
184
        self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
185
186

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

Nathan Coppin's avatar
Nathan Coppin committed
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
        v = self.velocity()
        omega = self.omega()
        x = self.position()
        tag = self._get_idx("ParticleTag").reshape(-1,1)
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
            x = np.insert(x,self._dim,0,axis=1)
        data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Tag",tag)]
        VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp")

        xdisk = self._get_matrix("Disk",2*self._dim+2)[:,:self._dim]
        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)
        tags = np.array(np.hstack([self.disk_tag(),self.segment_tag(),self.triangle_tag()]),dtype=np.float64)
        VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1]))])

        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)
        
229
    def read_vtk(self,dirname,i) :
230
        """Read an existing output file to set the particle data."""
Nathan Coppin's avatar
Nathan Coppin committed
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
259
260
261
262
263
264
        x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
        print(d["Omega"])
        self._lib.particleProblemClearParticles(self._p)
        self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
                _np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(d["Tag"]))

        x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
        self._lib.particleProblemClearBoundaries(self._p)
        offsets = np.hstack([[0],cells["offsets"]])
        el = cells["connectivity"]
        tags = cdata["Tag"]
        for idx in np.where(cells["types"] == 1)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
            t = int(tags[idx,0])
            self._lib.particleProblemAddBoundaryDiskTagId(self._p, self._coord_type(*x0), c_double(0.), c_int(t))
        for idx in np.where(cells["types"] == 3)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
            x1 = x[el[offsets[idx]+1],:self._dim]
            t = int(tags[idx,0])
            self._lib.particleProblemAddBoundarySegmentTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(t))
        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]
            t = int(tags[idx,0])
            self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(t))
        _,_,_,_,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))
265
266
267
268
269
270
271
272
273
274

    def add_boundary_disk(self, x0, r, tag) :
        self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
        return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode())

    def add_boundary_segment(self, x0, x1, tag) :
        self._lib.particleProblemAddBoundarySegment.restype = c_size_t
        return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode())

    def add_boundary_triangle(self, x0, x1, x2, tag) :
275
        if self._dim != 3 :
276
277
278
279
            raise NameError("Triangle boundaries only available in 3D")
        self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t
        return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), tag.encode())

Nathan Coppin's avatar
Nathan Coppin committed
280
    def add_particle(self, x, r, m, tag):
281
282
283
284
285
286
        """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
287
        tag -- particle tag
288
        """
Nathan Coppin's avatar
Nathan Coppin committed
289
        self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), c_int(tag))
290
291

    def load_msh_boundaries(self, filename, tags, shift=None) :
292
293
294
295
296
297
298
        """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
        """
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
        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])
                        self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag)
                    v0 = el.vertices[1]
                    v1 = el.vertices[0]
                    self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag)
        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])
                            self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag)
                        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]),
                                stag)
                            
                    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]),
                          stag)