scontact.py 28.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
35
36
37
38
39
40
41
import shutil
import os
import sys
from ctypes import *
import signal

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42
43
44
import numpy as np
assert(np.lib.NumpyVersion(np.__version__) >= "1.17.0")

45
dir_path = os.path.dirname(os.path.realpath(__file__))
Nathan Coppin's avatar
Nathan Coppin committed
46
lib2 = np.ctypeslib.load_library("libscontact2",dir_path)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
47
lib2f = np.ctypeslib.load_library("libscontact2f",dir_path)
48
lib2fwr = np.ctypeslib.load_library("libscontact2fwr",dir_path)
Nathan Coppin's avatar
Nathan Coppin committed
49
lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
Nathan Coppin's avatar
Nathan Coppin committed
50
lib3f = np.ctypeslib.load_library("libscontact3f",dir_path)
51
lib3fwr = np.ctypeslib.load_library("libscontact3fwr",dir_path)
52
53
54

is_64bits = sys.maxsize > 2**32

Nathan Coppin's avatar
Nathan Coppin committed
55
56
57
58
59
60
def _np2c(a,dtype=np.float64) :
    tmp = np.require(a,dtype,"C")
    r = c_void_p(tmp.ctypes.data)
    r.tmp = tmp
    return r

Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
61

62
class ParticleProblem :
Nathan Coppin's avatar
Nathan Coppin committed
63
    """Creates the numerical structure containing all the physical particles that appear in the problem"""
64

Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
65
66
67
68
69
70
71
    def _get_array(self, fName, dtype) :
        f = getattr(self._lib,"particleProblem"+fName)
        f.restype = c_void_p
        ptr = f(self._p)
        if ptr is None :
            return np.ndarray((0),dtype)
        size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//dtype.itemsize
72
        buf = (size*np.ctypeslib.as_ctypes_type(dtype)).from_address(ptr)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
73
        return np.ctypeslib.as_array(buf)
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
74

75
76
77
78
    def _get_matrix(self, fName, ncol) :
        f = getattr(self._lib,"particleProblem"+fName)
        f.restype = c_void_p
        ptr = f(self._p)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
79
80
        if ptr is None :
            return np.ndarray((0,ncol))
81
82
        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
83
        return np.ctypeslib.as_array(buf).reshape([-1,ncol])
84
85
86
87
88
    
    def _get_idx(self, fName) :
        f = getattr(self._lib,"particleProblem"+fName)
        f.restype = c_void_p
        ptr = f(self._p)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
89
90
        if ptr is None :
            return np.ndarray((0,),dtype=np.int32)
91
92
        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
93
        return np.ctypeslib.as_array(buf)
94

95
    def __init__(self, dim, friction_enabled=False, rotation_enabled=True) :
96
        """Builds the particles problem structure.
97
98

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
99
        dim -- Dimension of the problem
100
101

        Raises:
Nathan Coppin's avatar
Nathan Coppin committed
102
103
        ValueError -- If dimension differs from 2 or 3
        NameError -- If C builder cannot be found
104
        """
105
        self._dim = dim
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
106
        self._friction_enabled = friction_enabled
107
        self._rotation_enabled = False if not friction_enabled else rotation_enabled
108
        if dim == 2 :
Nathan Coppin's avatar
Nathan Coppin committed
109
            if friction_enabled:
110
              self._lib = lib2f if rotation_enabled else lib2fwr
Nathan Coppin's avatar
Nathan Coppin committed
111
112
113
            else:
              self._lib = lib2
            
114
115
            self._coord_type = c_double*2
        elif dim == 3 :
116
            self._lib = (lib3f if rotation_enabled else lib3fwr) if (friction_enabled) else lib3
117
118
119
120
121
122
123
124
            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
clean    
Jonathan Lambrechts committed
125
126
127
128
129
130
131
132
133
134
        bndtype =[('v',np.float64,dim)]
        if friction_enabled :
            bndtype.append(('vt',np.float64))
            if dim == 3 :
                bndtype.append(('vs',np.float64))
        bndtype.extend([('material',np.int32),('tag',np.int32)])
        self._disktype = np.dtype(bndtype+[('x',np.float64,dim),('r',np.float64)])
        self._segmenttype = np.dtype(bndtype+[('p',np.float64,(2,dim))])
        self._triangletype = np.dtype(bndtype+[('p',np.float64,(3,dim))])

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
135
    def __del__(self):
Nathan Coppin's avatar
Nathan Coppin committed
136
        """Deletes the particle solver structure."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
137
138
139
        if(self._p is not None) :
            self._lib.particleProblemDelete(self._p)

140
    def volume(self):
Nathan Coppin's avatar
Nathan Coppin committed
141
        """Returns the list of particle volumes."""
142
        if self._dim == 2 :
Nathan Coppin's avatar
Nathan Coppin committed
143
            return np.pi * self._get_matrix("Particle", 2)[:, [0]]**2
144
        else :
Nathan Coppin's avatar
Nathan Coppin committed
145
            return 4./3.*np.pi * self._get_matrix("Particle", 2)[:, [0]]**3
146
147

    def r(self) :
Nathan Coppin's avatar
Nathan Coppin committed
148
        """Returns the list of particle radii."""
Nathan Coppin's avatar
Nathan Coppin committed
149
        return self._get_matrix("Particle", 2)[:, 0][:,None]
150
151

    def mass(self):
Nathan Coppin's avatar
Nathan Coppin committed
152
        """Returns the list of particle masses."""
Nathan Coppin's avatar
Nathan Coppin committed
153
        return self._get_matrix("Particle", 2)[:, 1][:,None]
154

155
156
157
158
159
160
161
162
    def position(self):
        return self._get_matrix("Position", self._dim)

    def velocity(self):
        return self._get_matrix("Velocity", self._dim)

    def omega(self) : 
        d = 1 if self._dim==2 else 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
163
        if self._friction_enabled and self._rotation_enabled :
164
            return(self._get_matrix("Omega", d))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
165
        else :
166
            return np.zeros((self.n_particles(),d))
167

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
168
    def set_state(self, state) :
169
        """Sets the positions, velocities and angular velocities of the particle problem"""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
170
171
172
        self._get_matrix("Position", self._dim)[:] = state.x
        self._get_matrix("Velocity", self._dim)[:] = state.v
        if self._friction_enabled and self._rotation_enabled :
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
173
            self._get_matrix("Omega", 1 if self._dim==2 else 3)[:] = state.omega
174

175
176
177
178
179
180
181
182
183
184
185
186
    def save_state(self) :
        self._saved_velocity = np.copy(self.velocity())
        self._saved_position = np.copy(self.position())
        if self._friction_enabled and self._rotation_enabled :
            self._saved_omega = np.copy(self.omega())

    def restore_state(self) :
        self.velocity()[:] = self._saved_velocity
        self.position()[:] = self._saved_position
        if self._friction_enabled and self._rotation_enabled :
            self.omega()[:] = self._saved_omega

Nathan Coppin's avatar
Nathan Coppin committed
187
    def particle_material(self):
Nathan Coppin's avatar
Nathan Coppin committed
188
        """Returns the list of particle materials."""
Nathan Coppin's avatar
Nathan Coppin committed
189
        return self._get_idx("ParticleMaterial")
Nathan Coppin's avatar
Nathan Coppin committed
190
191

    def disk_tag(self):
Nathan Coppin's avatar
Nathan Coppin committed
192
        """Returns the list of boundary disk tag numbers."""
Nathan Coppin's avatar
Nathan Coppin committed
193
194
        return self._get_idx("DiskTag")

Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
195
    def disks(self) :
Nathan Coppin's avatar
Nathan Coppin committed
196
        """Returns the list of boundary disks."""
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
197
        return self._get_array("Disk",self._disktype)
198

199
    def forced_flag(self):
Nathan Coppin's avatar
Nathan Coppin committed
200
        """Returns the list of particle flags indicating whether they are fixed or not."""
201
202
203
        return self._get_idx("ForcedFlag")

    def get_tag_id(self, tag = "default"):
Nathan Coppin's avatar
Nathan Coppin committed
204
        """Returns the number associated to a string tag."""
205
206
207
        return self._lib.particleProblemGetMaterialTagId(self._p, tag.encode())

    def n_particles(self) :
Nathan Coppin's avatar
Nathan Coppin committed
208
        """Returns the number of particles."""
209
210
211
212
        self._lib.particleProblemNParticle.restype = c_size_t
        return self._lib.particleProblemNParticle(self._p)

    def mu(self) :
Nathan Coppin's avatar
Nathan Coppin committed
213
214
        """Returns the matrix of the friction coefficients between materials."""
        return self._get_matrix("Mu", self._lib.particleProblemNMaterial(self._p))
215

216
    def segments(self):
Nathan Coppin's avatar
Nathan Coppin committed
217
        """Returns the list of boundary segments."""
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
218
        return self._get_array("Segment",self._segmenttype)
219

220
    def triangles(self):
Nathan Coppin's avatar
Nathan Coppin committed
221
        """Returns the list of boundary triangles."""
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
222
223
224
        if self._dim == 2:
            return np.ndarray((0),self._triangletype)
        return self._get_array("Triangle",self._triangletype)
225
    
226
    def contact_forces(self):
Nathan Coppin's avatar
Nathan Coppin committed
227
        """Returns the contact forces on each grain."""
Nathan Coppin's avatar
Nathan Coppin committed
228
        return (self._get_matrix("ContactForces",self._dim))
229

230
231
    def get_boundary_forces(self,tag="default") :
        """Returns the net normal and tangential forces acting on a boundary because of the contacts.
Nathan Coppin's avatar
Nathan Coppin committed
232
233
234
        Keyword arguments:
        tag -- Name of the boundary
        """
Nathan Coppin's avatar
Nathan Coppin committed
235
236
        self._lib.particleProblemGetTagId.restype = c_size_t
        tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
        def compute_fn_ft(contact_name,objects,tag) :
            o,v,b =self.get_contacts(contact_name)
            keep = objects["tag"][o[:,0]] == tag
            m = self.mass()[o[keep,1]]
            if m.shape[0] == 0 :
                return np.array([[0.]*self._dim]*2)
            fn = np.sum(m*b[keep,0,:]*v[keep,:][:,[0]],axis=0)
            ft = np.sum(m*b[keep,1,:]*v[keep,:][:,[1]],axis=0)
            if self._dim == 3:
                ft += np.sum(m*b[keep,2,:]*v[keep,:][:,[2]],axis=0)
            return np.array([fn,ft])
        f = compute_fn_ft("particle_disk",self.disks(),tag)
        f2 = compute_fn_ft("particle_segment",self.segments(),tag)
        f += f2
Nathan Coppin's avatar
Nathan Coppin committed
251
        if self._dim == 3:
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
252
            f+= compute_fn_ft("particle_triangle",self.triangles(),tag)
253
        return f
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
254

255
    def set_boundary_velocity(self, tag, v) :
Nathan Coppin's avatar
Nathan Coppin committed
256
257
258
259
260
        """Sets the velocity of a boundary to a given value.
        Keyword arguments:
        tag -- Name of the boundary
        v -- Velocity vector to be given to the boundary
        """
261
262
        self._lib.particleProblemGetTagId.restype = c_size_t
        tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
263
        disks = self.disks()
Nathan Coppin's avatar
Nathan Coppin committed
264
        disks["v"][disks["tag"]==tag] = v
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
265
        segments = self.segments()
Nathan Coppin's avatar
Nathan Coppin committed
266
        segments["v"][segments["tag"]==tag] = v
267
        if self._dim == 3 :
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
268
            triangles = self.triangles()
Nathan Coppin's avatar
Nathan Coppin committed
269
            triangles["v"][triangles["tag"]==tag] = v
270

Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
271
    def set_tangent_boundary_velocity(self, tag, vt):
272
        """ Sets the tangent velocity of the frictional boundary.
Nathan Coppin's avatar
Nathan Coppin committed
273
        Only in 2D
274
275
        Keyword arguments:
        tag -- Name of the boundary
Nathan Coppin's avatar
Nathan Coppin committed
276
        vt -- Velocity in the tangent direction
277
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
278
        assert self._friction_enabled
Nathan Coppin's avatar
Nathan Coppin committed
279
        assert self._dim==2
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
280
281
282
283
        self._lib.particleProblemGetTagId.restype = c_size_t
        tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
        disks = self.disks()
        segments = self.segments()
Nathan Coppin's avatar
Nathan Coppin committed
284
285
286
        disks["vt"][disks["tag"]==tag] = vt
        segments["vt"][segments["tag"]==tag] = vt

287
    def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
288
289
290
        """ Sets the friction coefficient between two materials.

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
291
292
293
        mu -- Value of the friction coefficient
        mat0 -- First material
        mat1 -- Second material
294
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
295
        assert self._friction_enabled
296
        self._lib.particleProblemSetFrictionCoefficient(self._p, c_double(mu), mat0.encode(),mat1.encode())
Nathan Coppin's avatar
Nathan Coppin committed
297

298
    def iterate(self, dt, forces,tol=1e-8,force_motion=0) :
Nathan Coppin's avatar
Nathan Coppin committed
299
        """Computes iteratively the set of velocities that respect the non-interpenetration constrain.
300
           Returns 1 if the computation converged, 0 otherwise.
301
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
302
303
304
305
        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
        force_motion -- Optional argument allowing to move the grains if convergence has not been reached when set to 1
306
        """ 
307
        self._lib.particleProblemIterate.restype = c_int
Nathan Coppin's avatar
Nathan Coppin committed
308
        return self._lib.particleProblemIterate(self._p, c_double(np.min(self.r()) if self.r().size != 0 else tol), c_double(dt), c_double(tol), c_int(-1),c_int(force_motion),_np2c(forces))
309

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
310
    def get_contacts(self,contact_type) :
311
        """Gives the contact forces. Warning : during the resolution of the contacts, 
Nathan Coppin's avatar
Nathan Coppin committed
312
313
           the considered quantities are impulses, while forces are written and read.
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
314
315
        ctype = {"particle_particle":0,"particle_disk":1,"particle_segment":2,"particle_triangle":3}[contact_type]
        n = self._lib.particleProblemContactN(self._p,c_int(ctype))
316
317
        basis = np.ndarray((n,self._dim*self._dim),dtype=np.float64,order="C")
        v = np.ndarray((n,self._dim),dtype=np.float64,order="C")
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
318
319
        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))
320
321
322
        self._lib.particleProblemContactBasis(self._p,c_int(ctype),c_void_p(basis.ctypes.data))
        basis = basis.reshape(n,self._dim, self._dim)
        return o,v,basis
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
323

324
    def computeStressTensor(self, nodes, radius):
Nathan Coppin's avatar
Nathan Coppin committed
325
326
327
328
329
        """Computes the stress tensor of the granular material
        Keyword arguments:
        nodes -- Array of nodes at which to compute the stress tensor
        r -- Radius within which the contact forces will be averaged 
        """
330
        n_nodes = len(nodes[:,0])
331
        s = np.ndarray((n_nodes,self._dim**2))
332
333
        self._lib.particleProblemComputeStressTensor(self._p,_np2c(nodes[:,:self._dim]),c_double(radius),c_int(n_nodes),_np2c(s))
        return s
334

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
335
    def write_vtk(self, odir, i, t) :
Nathan Coppin's avatar
Nathan Coppin committed
336
337
338
339
340
341
        """Writes output files for post-visualization.
        Keyword arguments:
        odir -- Directory in which to write the file
        i -- Number of the fiel to write
        t -- Time at which the simulation is 
        """
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
342
343
        state = self.state()
        v = state.v
344
        omega = state.omega
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
345
        x = state.x
346
        material = self._get_idx("ParticleMaterial").reshape(-1,1)
347
        forced = self._get_idx("ForcedFlag").reshape(-1,1)
Nathan Coppin's avatar
Nathan Coppin committed
348
349
350
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
            x = np.insert(x,self._dim,0,axis=1)
351
        data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Material",material),("ForcedFlag",forced),("Contacts",self.contact_forces())]
352
353
354
355
        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))])
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
356
357
358
359
        disks = self.disks()
        segments = self.segments()
        triangles = self.triangles()
        x = np.vstack([disks["x"],segments["p"].reshape([-1,self._dim]),triangles["p"].reshape([-1,self._dim])])
Nathan Coppin's avatar
Nathan Coppin committed
360
361
362
        if self._dim == 2 :
            x = np.insert(x,self._dim,0,axis=1)
        connectivity = np.arange(0,x.shape[0],dtype=np.int32)
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
363
364
365
366
        types = np.repeat(np.array([1,3,5],dtype=np.int32),[disks.shape[0],segments.shape[0],triangles.shape[0]])
        offsets = np.cumsum(np.repeat([1,2,3],[disks.shape[0],segments.shape[0],triangles.shape[0]]),dtype=np.int32)
        tags = np.array(np.hstack([disks["tag"],segments["tag"],triangles["tag"]]))
        materials = np.array(np.hstack([disks["material"],segments["material"],triangles["material"]]))
367
368
369
370
        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
371
372
        self._lib.particleProblemContactN.restype = c_size_t
        fdata = []
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
373
        for name in ("particle_particle","particle_disk","particle_segment","particle_triangle") :
374
            o,v,basis = self.get_contacts(name)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
375
376
377
            nameb = name.encode()
            fdata.append((nameb,v[:,[0]]))
            fdata.append((nameb+b"_t",v[:,[1]]))
378
            fdata.append((nameb+b"_dir_n",basis[:,0,:]))
379
380
            if self._dim == 3:
                fdata.append((nameb+b"_s",v[:,[2]]))
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
381
            fdata.append((nameb+b"_idx",o))
Nathan Coppin's avatar
Nathan Coppin committed
382
383
384
385
        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)
        
386
    def read_vtk(self,dirname,i,contact_factor=1) :
Nathan Coppin's avatar
Nathan Coppin committed
387
388
389
390
391
392
        """Reads an existing output file to set the particle data.
        Keyword arguments:
        dirname -- Path to the directory of the file to read
        i -- Number of the file to read
        contact_factor -- Factor that determines how to take the read contacts into account
        """
393
        self._lib.particleProblemClearBoundaries(self._p)
Nathan Coppin's avatar
Nathan Coppin committed
394
        x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
395
396
        matnames =  VTK.string_array_decode(fdata["MaterialNames"])
        matmap = {}
397
        forced = d["ForcedFlag"] if ("ForcedFlag" in d) else np.zeros_like(d["Mass"]).astype(int)
398
399
400
        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
401
402
        self._lib.particleProblemClearParticles(self._p)
        self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
Nathan Coppin's avatar
Nathan Coppin committed
403
                _np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(partmat,dtype=np.int32),_np2c(forced,dtype=np.int32),_np2c(d["Contacts"][:,:self._dim]*contact_factor))
Nathan Coppin's avatar
Nathan Coppin committed
404
        x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
405
406
407
408
        tagnames =  VTK.string_array_decode(fdata["TagNames"])
        tagmap = {}
        for j,n in enumerate(tagnames) :
            tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
409
        
Nathan Coppin's avatar
Nathan Coppin committed
410
411
        offsets = np.hstack([[0],cells["offsets"]])
        el = cells["connectivity"]
412
413
        tags = np.vectorize(tagmap.get)(cdata["Tag"])
        materials = cdata["Material"]
Nathan Coppin's avatar
Nathan Coppin committed
414
415
        for idx in np.where(cells["types"] == 1)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
416
            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
417
418
419
        for idx in np.where(cells["types"] == 3)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
            x1 = x[el[offsets[idx]+1],:self._dim]
420
            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
421
422
423
424
        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]
425
            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
426
427
        _,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
        ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
428
        v = []
429
        basis = []
430
431
        v.append(np.vstack([fdata[k] for k in ks]))
        v.append(np.vstack([fdata[k+"_t"] for k in ks]))
432
433
434
435
436
437
438
439
440
        basis.append(np.vstack([fdata[k+"_dir_n"] for k in ks]))
        if self._dim == 2:
            basis.append(np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[1,0]]*np.array([-1, 1]))
        else:
            dirt = np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[2,0,1]]*np.array([-1, 1, -1]) - np.vstack([fdata[k+"_dir_n"] for k in ks])*np.einsum("ij,ij->i",np.vstack([fdata[k+"_dir_n"] for k in ks]),np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[2,0,1]]*np.array([-1, 1, -1]))[:,np.newaxis]
            nnorm = np.linalg.norm(dirt,axis=1)
            nnorm[nnorm==0] = 1
            basis.append(dirt/nnorm[:,np.newaxis])
            basis.append(np.cross(np.vstack([fdata[k+"_dir_n"] for k in ks]),dirt/nnorm[:,np.newaxis]))
441
            v.append(np.vstack([fdata[k+"_s"] for k in ks]))
442
        basis = np.column_stack(basis)
443
        v = np.column_stack(v)*contact_factor
Nathan Coppin's avatar
Nathan Coppin committed
444
445
        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])
446
        self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(basis),_np2c(t,np.int32))
447

448
    def add_boundary_disk(self, x0, r, tag, material="default") :
Nathan Coppin's avatar
Nathan Coppin committed
449
450
451
452
453
454
455
456
        """Adds a boundary disk.

        Keyword arguments:
        x0 -- Tuple of the coordinates of the centre
        r -- Disk radius
        tag -- Disk tag
        material -- Disk material
        """
457
        self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
458
        return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode(),material.encode())
459

460
    def add_boundary_segment(self, x0, x1, tag, material="default") :
Nathan Coppin's avatar
Nathan Coppin committed
461
462
463
464
465
466
467
468
        """Adds a boundary segment.

        Keyword arguments:
        x0 -- Tuple of the coordinates of the first endpoint
        x1 -- Tuple of the coordinates of the second endpoint
        tag -- Segment tag
        material -- Segment material
        """
469
        self._lib.particleProblemAddBoundarySegment.restype = c_size_t
470
        return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode(),material.encode())
471

472
    def add_boundary_triangle(self, x0, x1, x2, tag, material="default") :
Nathan Coppin's avatar
Nathan Coppin committed
473
474
475
476
477
478
479
480
481
        """Adds a boundary triangle.

        Keyword arguments:
        x0 -- Tuple of the coordinates of the first summit 
        x1 -- Tuple of the coordinates of the second summit 
        x2 -- Tuple of the coordinates of the third summit 
        tag -- Triangle tag
        material -- Triangle material
        """
482
        if self._dim != 3 :
483
484
            raise NameError("Triangle boundaries only available in 3D")
        self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t
485
        return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), tag.encode(),material.encode())
486

487
    def add_particle(self, x, r, m, tag="default",forced = 0):
488
        """Adds a particle in the particle problem.
489
490

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
491
492
493
494
495
        x -- Tuple to set the centre particle position
        r -- Particle radius
        m -- Particle mass
        tag -- Particle material
        forced -- Flag indicating whether the particle is forced or not
496
        """
497
498
        self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), tag.encode(), c_int(forced))

499
    def add_particles(self, x, r, m, v=None, omega=None, tag="default", forced=None, contact_forces=None):
500
        """Adds particles in the particle problem.
501
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
502
503
504
505
506
507
508
509
        x -- Array with the particles's centers position
        r -- Array with the particles's radii
        m -- Array with the particles's masses
        v -- Array with the particles's velocities
        omega -- Array with the particles's angular velocities
        tag -- Particle material
        forced -- Array of flags indicating whether the particles are forced or not
        contact_forces -- Array of contact forces between the particles
510
        """
511
512
513
514
515
516
517
518
519
520
521
        n_particles = len(m)
        tags = np.ones(n_particles) * self.get_tag_id(tag)
        if omega is None :
         omega = np.zeros((1 if self._dim == 2 else 3)*n_particles)
        if v is None :
         v = np.zeros(self._dim*n_particles)
        if forced is None :
         forced = np.zeros(n_particles)
        if contact_forces is None :
         contact_forces = np.zeros(self._dim*n_particles)
        self._lib.particleProblemAddParticles(self._p, c_int(n_particles), _np2c(x), _np2c(m), _np2c(r),_np2c(v), _np2c(omega), _np2c(tags,dtype=np.int32), _np2c(forced,dtype=np.int32), _np2c(contact_forces))
522

523
    def set_use_queue(self, use_queue=1):
Nathan Coppin's avatar
Nathan Coppin committed
524
      """Enables the use of the queue if 1 and disables it if 0."""
525
526
      self._lib.particleProblemSetUseQueue(self._p, c_int(use_queue))

Nathan Coppin's avatar
Nathan Coppin committed
527
    def get_tagnames(self):
Nathan Coppin's avatar
Nathan Coppin committed
528
      """Returns the names of the boundaries """
529
530
531
532
      self._lib.particleProblemGetTagName.restype = c_char_p
      ntagname = self._lib.particleProblemNTag(self._p)
      tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)]) 
      return tagnames
Nathan Coppin's avatar
Nathan Coppin committed
533
534

    def clear_boundaries(self):
Nathan Coppin's avatar
Nathan Coppin committed
535
      """Removes the boundaries."""
536
      self._lib.particleProblemClearBoundaries(self._p)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
537
    def _save_contacts(self):
Nathan Coppin's avatar
Nathan Coppin committed
538
      """Saves the contacts from the current time step."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
539
540
541
      self._lib.particleProblemSaveContacts(self._p)

    def _restore_contacts(self):
Nathan Coppin's avatar
Nathan Coppin committed
542
      """Restores the saved contacts from the previous time step."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
543
544
      self._lib.particleProblemRestoreContacts(self._p)

545
    def remove_particles_flag(self,flag) :
Nathan Coppin's avatar
Nathan Coppin committed
546
      """Removes particles based on given flag."""
547
548
549
      if flag.shape != (self.n_particles(),) :
          raise NameError("size of flag array should be the number of particles")
      self._lib.particleProblemRemoveParticles(self._p, _np2c(flag,np.int32))
550

551
    def load_msh_boundaries(self, filename, tags, shift=None,material="default") :
Nathan Coppin's avatar
Nathan Coppin committed
552
        """Loads the numerical domain and set the physical boundaries the particles cannot go through.
553
554
        
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
555
556
557
558
        filename -- Name of the msh file
        tags -- Tags of physical boundary defined in the msh file
        shift -- Optional argument to shift the numerical domain
        material -- Material of the boundary
559
        """
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
        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])
576
                        self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag,material)
577
578
                    v0 = el.vertices[1]
                    v1 = el.vertices[0]
579
                    self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag,material)
580
581
582
583
584
585
586
587
588
589
590
591
592
        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])
593
                            self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag,material)
594
595
596
597
598
599
                        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]),
600
                                stag,material)
601
602
603
604
605
606
607
608
                            
                    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]),
609
                          stag,material)
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
610
611
    
    def dim(self) :
612
        """Returns the dimension of the particle problem"""
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
613
        return self._dim