scontact.py 35.7 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
4
#
5
6
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
7
8
9
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
10
11
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
12
#
13
14
15
16
# 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.
17
#
18
# You should have received a copy of the GNU Lesser General Public License
19
# along with this program (see COPYING and COPYING.LESSER files).  If not,
20
21
# see <http://www.gnu.org/licenses/>.

22
23
24
25
26
"""Model for Immersed Granular Flow -- Particles user interface

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

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

31
from __future__ import division
Nathan Coppin's avatar
Nathan Coppin committed
32
from . import VTK
33
34
35
36
import shutil
import os
import sys
from ctypes import *
37
from ._tools import gmsh
38
39


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

43
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)
45
lib2fwr = np.ctypeslib.load_library("libscontact2fwr",dir_path)
Nathan Coppin's avatar
Nathan Coppin committed
46
lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
47
lib3fwr = np.ctypeslib.load_library("libscontact3fwr",dir_path)
48
49
50

is_64bits = sys.maxsize > 2**32

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

57
class ParticleProblem :
Nathan Coppin's avatar
Nathan Coppin committed
58
    """Creates the numerical structure containing all the physical particles that appear in the problem"""
59

Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
60
    def _get_array(self, fName, dtype) :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
61
        f = getattr(self._lib,"particle_problem_"+fName)
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
62
63
64
65
66
        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
67
        buf = (size*np.ctypeslib.as_ctypes_type(dtype)).from_address(ptr)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
68
        return np.ctypeslib.as_array(buf)
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
69

70
    def _get_matrix(self, fName, ncol) :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
71
        f = getattr(self._lib,"particle_problem_"+fName)
72
73
        f.restype = c_void_p
        ptr = f(self._p)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
74
75
        if ptr is None :
            return np.ndarray((0,ncol))
76
77
        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
78
        return np.ctypeslib.as_array(buf).reshape([-1,ncol])
79

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
80
    def __init__(self, dim) :
81
        """Builds the particles problem structure.
82
83

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
84
        dim -- Dimension of the problem
85
86

        Raises:
Nathan Coppin's avatar
Nathan Coppin committed
87
88
        ValueError -- If dimension differs from 2 or 3
        NameError -- If C builder cannot be found
89
        """
90
91
        self._dim = dim
        if dim == 2 :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
92
            self._lib = lib2
93
94
            self._coord_type = c_double*2
        elif dim == 3 :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
95
            self._lib = lib3
96
97
98
            self._coord_type = c_double*3
        else :
            raise ValueError("Dimension should be 2 or 3.")
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
99
100
        self._lib.particle_problem_new.restype = c_void_p
        self._p = c_void_p(self._lib.particle_problem_new())
101
102
        if self._p == None :
            raise NameError("Cannot create particle problem.")
103
        bndtype =[('material',np.int32),('tag',np.int32),('body',np.uint64)]
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
104
        self._particletype = np.dtype([('r',np.float64),('body',np.uint64),('x',np.float64,dim),('material',np.int32)], align=True)
105
106
        self._segmenttype = np.dtype(bndtype+[('x',np.float64,(2,dim))])
        self._triangletype = np.dtype(bndtype+[('x',np.float64,(3,dim))])
Michel Henry's avatar
Michel Henry committed
107
108
109
        self._periodicEntitytype = np.dtype([('etag', np.int32),('edim', np.int32),('periodic_transformation', np.float64, dim)])
        self._periodicSegmenttype = np.dtype([('entity_id', np.int64), ('p',np.float64,(2,dim))])
        self._periodicTriangletype = np.dtype([('entity_id', np.int64),('p', np.float64,(3,dim))])
110
        self._bodytype = np.dtype([('invertM', np.float64),('invertI', np.float64), ('theta', np.float64),('omega',np.float64),('position',np.float64,dim),('velocity',np.float64,dim)], align=True)
111
        self._contacttype = np.dtype([('o', np.uint64,(2,)), ('type', np.int32), ('space',np.float64),('r', np.float64,(2,self._dim)) , ('reaction', np.float64,(self._dim,))],align=True)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
112
113
114
115
        self.material2id = {}
        self.tag2id = {}
        self.id2tag = []
        self.id2material = []
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
116

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
117
    def __del__(self):
Nathan Coppin's avatar
Nathan Coppin committed
118
        """Deletes the particle solver structure."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
119
        if(self._p is not None) :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
120
            self._lib.particle_problem_delete(self._p)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121

122
    def volume(self):
Nathan Coppin's avatar
Nathan Coppin committed
123
        """Returns the list of particle volumes."""
124
        if self._dim == 2 :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
125
            return np.pi * self._get_array("particle", self._particletype)["r"].reshape(-1,1)**2
126
        else :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
127
            return 4./3.*np.pi * self._get_array("particle", self._particletype)["r"].reshape(-1,1)**3
128
129

    def r(self) :
Nathan Coppin's avatar
Nathan Coppin committed
130
        """Returns the list of particle radii."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
131
        return self._get_array("particle", self._particletype)["r"].reshape(-1,1)
132

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
133
    def body_invert_mass(self):
Nathan Coppin's avatar
Nathan Coppin committed
134
        """Returns the list of body masses."""
135
        return self._get_array("body", self._bodytype)["invertM"][:,None]
136

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
137
    def body_invert_inertia(self):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
138
        """Returns the list of body masses."""
139
        return self._get_array("body", self._bodytype)["invertI"][:,None]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
140

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
141
    def body_position(self):
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
142
        return self._get_array("body", self._bodytype)["position"]
143

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
144
    def position(self):
Nathan Coppin's avatar
Nathan Coppin committed
145
        ppos = np.zeros((self.n_particles(),self._dim))
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
146
147
        self._lib.particle_problem_get_particles_position(self._p, _np2c(ppos))
        ppos.flags.writeable = False
Nathan Coppin's avatar
Nathan Coppin committed
148
149
        return ppos

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
150
    def velocity(self):
Nathan Coppin's avatar
Nathan Coppin committed
151
        pvel = np.zeros((self.n_particles(),self._dim))
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
152
153
        self._lib.particle_problem_get_particles_velocity(self._p, _np2c(pvel))
        pvel.flags.writeable = False
Nathan Coppin's avatar
Nathan Coppin committed
154
155
        return pvel

Nathan Coppin's avatar
Nathan Coppin committed
156
157
158
159
160
    def delassus(self):
        delassus = np.zeros((self.n_particles(),self._dim,self._dim))
        self._lib.particle_problem_get_delassus(self._p,_np2c(delassus))
        return delassus

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
161
    def body_velocity(self):
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
162
        return self._get_array("body", self._bodytype)["velocity"]
Nathan Coppin's avatar
Nathan Coppin committed
163

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
164
    def body_omega(self) : 
165
        d = 1 if self._dim==2 else 3
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
166
167
168
169
        if self._dim == 2:
            return self._get_array("body", self._bodytype)["omega"][:,None]
        else:
            return self._get_array("body", self._bodytype)["omega"]
170

Nathan Coppin's avatar
Nathan Coppin committed
171
172
173
    def body_theta(self) :
        return self._get_array("body", self._bodytype)["theta"][:,None]

174
    def save_state(self) :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
175
176
        self._saved_velocity = np.copy(self.body_velocity())
        self._saved_position = np.copy(self.body_position())
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
177
178
179
        self._saved_segments = np.copy(self.segments())
        if self.dim() == 3 :
            self._saved_triangles = np.copy(self.triangles())
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
180
        self._saved_omega = np.copy(self.body_omega())
Nathan Coppin's avatar
Nathan Coppin committed
181
        self._saved_theta = np.copy(self.body_theta())
182
183

    def restore_state(self) :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
184
185
186
        self.body_velocity()[:] = self._saved_velocity
        self.body_position()[:] = self._saved_position
        self.body_omega()[:] = self._saved_omega
Nathan Coppin's avatar
Nathan Coppin committed
187
        self.body_theta()[:] = self._saved_theta
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
188
189
190
        self.segments()[:] = self._saved_segments
        if self.dim() == 3 :
            self.triangles()[:] = self._saved_triangles
191

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
192
    def material(self):
Nathan Coppin's avatar
Nathan Coppin committed
193
        """Returns the list of particle materials."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
194
        return self._get_array("particle", self._particletype)["material"].reshape(-1,1)
Nathan Coppin's avatar
Nathan Coppin committed
195

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
196
    def get_tag_id(self, tag="default"):
Nathan Coppin's avatar
Nathan Coppin committed
197
        """Returns the number associated to a string tag."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
        if tag not in self.tag2id.keys():
            n = len(self.tag2id)
            self.tag2id[tag] = n
            self.id2tag.append(tag)
            return n
        return self.tag2id[tag]

    def get_material_id(self, material):
        """Returns the number associated to a string tag."""
        if material not in self.material2id:
            n = len(self.material2id)
            self.material2id[material] = n
            self._lib.particle_problem_allocate_material(self._p, c_int(n+1))
            self.id2material.append(material)
            return n
        return self.material2id[material]
214
215

    def n_particles(self) :
Nathan Coppin's avatar
Nathan Coppin committed
216
        """Returns the number of particles."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
217
        return self._get_array("particle",self._particletype).shape[0]
218

Nathan Coppin's avatar
Nathan Coppin committed
219
220
    def n_bodies(self) :
        """Returns the number of bodies."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
221
        return self._get_array("body",self._bodytype).shape[0]
Nathan Coppin's avatar
Nathan Coppin committed
222

223
    def mu(self) :
Nathan Coppin's avatar
Nathan Coppin committed
224
        """Returns the matrix of the friction coefficients between materials."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
225
        return self._get_matrix("mu", len(self.material2id))
226

227
    def segments(self):
Nathan Coppin's avatar
Nathan Coppin committed
228
        """Returns the list of boundary segments."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
229
        return self._get_array("segment",self._segmenttype)
Michel Henry's avatar
Michel Henry committed
230

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
231
    def periodic_entity(self):
Michel Henry's avatar
Michel Henry committed
232
        """Returns the list of periodic entity."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
233
        return self._get_array("periodic_entity", self._periodicEntitytype)
Michel Henry's avatar
Michel Henry committed
234

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
235
    def periodic_segments(self):
Michel Henry's avatar
Michel Henry committed
236
        """Returns the list of periodic segments."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
237
        return self._get_array("periodic_segment", self._periodicSegmenttype)
Michel Henry's avatar
Michel Henry committed
238

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
239
    def periodic_triangles(self):
Michel Henry's avatar
Michel Henry committed
240
241
242
        """Returns the list of periodic triangles."""
        if self._dim == 2:
            return np.ndarray((0),self._periodicTriangletype)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
243
        return self._get_array("periodic_triangle", self._periodicTriangletype)
Michel Henry's avatar
Michel Henry committed
244

245
    def triangles(self):
Nathan Coppin's avatar
Nathan Coppin committed
246
        """Returns the list of boundary triangles."""
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
247
248
        if self._dim == 2:
            return np.ndarray((0),self._triangletype)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
249
        return self._get_array("triangle",self._triangletype)
250

251
252
    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
253
254
255
        Keyword arguments:
        tag -- Name of the boundary
        """
Nathan Coppin's avatar
Nathan Coppin committed
256
257
        self._lib.particleProblemGetTagId.restype = c_size_t
        tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
258
259
260
261
262
263
264
265
266
267
268
269
270
271
        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
272
        if self._dim == 3:
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
273
            f+= compute_fn_ft("particle_triangle",self.triangles(),tag)
274
        return -f
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
275

276
    def set_boundary_velocity(self, tag, v) :
Nathan Coppin's avatar
Nathan Coppin committed
277
278
279
280
281
        """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
        """
282
283
        self._lib.particleProblemGetTagId.restype = c_size_t
        tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
284
        disks = self.disks()
Nathan Coppin's avatar
Nathan Coppin committed
285
        disks["v"][disks["tag"]==tag] = v
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
286
        segments = self.segments()
Nathan Coppin's avatar
Nathan Coppin committed
287
        segments["v"][segments["tag"]==tag] = v
288
        if self._dim == 3 :
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
289
            triangles = self.triangles()
Nathan Coppin's avatar
Nathan Coppin committed
290
            triangles["v"][triangles["tag"]==tag] = v
291

Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
292
    def set_tangent_boundary_velocity(self, tag, vt):
293
        """ Sets the tangent velocity of the frictional boundary.
Nathan Coppin's avatar
Nathan Coppin committed
294
        Only in 2D
295
296
        Keyword arguments:
        tag -- Name of the boundary
Nathan Coppin's avatar
Nathan Coppin committed
297
        vt -- Velocity in the tangent direction
298
        """
Nathan Coppin's avatar
Nathan Coppin committed
299
        assert self._dim==2
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
300
301
302
303
        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
304
305
306
        disks["vt"][disks["tag"]==tag] = vt
        segments["vt"][segments["tag"]==tag] = vt

307
    def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
308
309
310
        """ Sets the friction coefficient between two materials.

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
311
312
313
        mu -- Value of the friction coefficient
        mat0 -- First material
        mat1 -- Second material
314
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
315
316
317
        imat0 = self.get_material_id(mat0)
        imat1 = self.get_material_id(mat1)
        self._lib.particle_problem_set_friction_coefficient(self._p, c_double(mu), c_int(imat0), c_int(imat1))
Nathan Coppin's avatar
Nathan Coppin committed
318

319
    def iterate(self, dt, forces,tol=1e-8,force_motion=0) :
Nathan Coppin's avatar
Nathan Coppin committed
320
        """Computes iteratively the set of velocities that respect the non-interpenetration constrain.
321
           Returns 1 if the computation converged, 0 otherwise.
322
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
323
324
325
326
        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
327
        """ 
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
328
329
        self._lib.particle_problem_iterate.restype = c_int
        return self._lib.particle_problem_iterate(self._p, c_double(np.max(self.r()) if self.r().size != 0 else tol), c_double(dt), c_double(tol), c_int(-1),c_int(force_motion),_np2c(forces))
330

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
331
    def contacts(self) :
332
        """Gives the contact forces.
Nathan Coppin's avatar
Nathan Coppin committed
333
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
334
        return self._get_array("contact", self._contacttype)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
335

336
    def computeStressTensor(self, nodes, radius):
Nathan Coppin's avatar
Nathan Coppin committed
337
338
339
340
341
        """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 
        """
342
        n_nodes = len(nodes[:,0])
343
        s = np.ndarray((n_nodes,self._dim**2))
344
        self._lib.particle_problem_compute_stress_tensor(self._p,_np2c(nodes[:,:self._dim]),c_double(radius),c_int(n_nodes),_np2c(s))
345
        return s
346

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
347
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
348
    def write_vtk(self, odir, i, t) :
Nathan Coppin's avatar
Nathan Coppin committed
349
350
351
352
353
        """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 
Nathan Coppin's avatar
Nathan Coppin committed
354
        """      
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
355
        bodies = self._get_array("body",self._bodytype)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
356
357
358
359
360
361
362
363
364
        v = self.body_velocity()
        x = self.body_position()
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
            x = np.insert(x,self._dim,0,axis=1)
        bodies = self._get_array("body", self._bodytype)
        data = [("Velocity",v),
                ("Omega",self.body_omega()),
                ("Theta",bodies["theta"][:,None]),
365
366
                ("iMass",bodies["invertM"][:,None]),
                ("iInertia",bodies["invertI"][:,None]),
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
367
                ]
368
369
        sint = np.sin(bodies["theta"])
        cost = np.cos(bodies["theta"])
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
370
371
372
373
374
        VTK.write(odir+"/bodies",i,t,None,x,data,vtktype="vtp")
        x = self.position()
        v = self.velocity()
        particles = self._get_array("particle", self._particletype)
        xrel = particles["x"]
Nathan Coppin's avatar
Nathan Coppin committed
375
376
377
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
            x = np.insert(x,self._dim,0,axis=1)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
378
379
380
381
382
383
384
            xrel = np.insert(xrel,self._dim,0,axis=1)
        data = [("Radius",particles["r"][:,None]),
                ("Body",particles["body"][:,None]),
                ("Material",particles["material"][:,None]),
                ("Velocity",v),
                ("RelativePosition",xrel)]
        tags = list (s.encode() for s in self.id2material)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
385
386
        mu = self.mu()
        VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags)),(b"FrictionCoefficients",mu.reshape([-1,1]))])
Michel Henry's avatar
Michel Henry committed
387

Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
388
389
        segments = self.segments()
        triangles = self.triangles()
390
391
392
393
394
395
396
397
398
399
        def absolute_coord(o):
            xb = bodies["position"][o["body"]][:,None,:]
            ct = cost[o["body"]][:,None]
            st = sint[o["body"]][:,None]
            xo = np.empty_like(o["x"])
            xo[:,:,0] = xb[:,:,0]+o["x"][:,:,0]*ct-o["x"][:,:,1]*st
            xo[:,:,1] = xb[:,:,1]+o["x"][:,:,0]*st+o["x"][:,:,1]*ct
            return xo
        xsegs = absolute_coord(segments).reshape([-1,self._dim])
        xtris = absolute_coord(triangles).reshape([-1,self._dim])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
400
401
        x = np.vstack([xsegs,xtris])
        xrel = np.vstack([segments["x"].reshape(-1,self._dim),triangles["x"].reshape(-1,self._dim)])
Nathan Coppin's avatar
Nathan Coppin committed
402
403
        if self._dim == 2 :
            x = np.insert(x,self._dim,0,axis=1)
404
            xrel = np.insert(xrel,self._dim,0,axis=1)
Nathan Coppin's avatar
Nathan Coppin committed
405
        connectivity = np.arange(0,x.shape[0],dtype=np.int32)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
406
407
408
409
        types = np.repeat(np.array([3,5],dtype=np.int32),[segments.shape[0],triangles.shape[0]])
        offsets = np.cumsum(np.repeat([2,3],[segments.shape[0],triangles.shape[0]]),dtype=np.int32)
        tags = np.array(np.hstack([segments["tag"],triangles["tag"]]))
        materials = np.array(np.hstack([segments["material"],triangles["material"]]))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
410
        bodies = np.array(np.hstack([segments["body"],triangles["body"]]))
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
411
        tagnames = list(s.encode() for s in self.id2tag)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
412
413
        cell_data = [("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1])),("Body",bodies.reshape([-1,1]))]
        VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,data=[("RelativePosition",xrel)],cell_data=cell_data,field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
414
415
416
        periodicEntity = self.periodic_entity()
        periodicSegments = self.periodic_segments()
        periodicTriangles = self.periodic_triangles()
Michel Henry's avatar
Michel Henry committed
417
        periodicX = np.vstack([periodicSegments["p"].reshape([-1,self._dim]), periodicTriangles["p"].reshape([-1,self._dim])])
Michel Henry's avatar
Michel Henry committed
418
        if self._dim == 2 :
Michel Henry's avatar
Michel Henry committed
419
            periodicX = np.insert(periodicX,self._dim, 0, axis=1)
Michel Henry's avatar
Michel Henry committed
420
        periodicConnectivity = np.arange(0,periodicX.shape[0], dtype = np.int32)
Michel Henry's avatar
Michel Henry committed
421
422
423
424
425
        periodicTypes = np.repeat(np.array([3,5],dtype=np.int32),[periodicSegments.shape[0], periodicTriangles.shape[0]])
        periodicOffsets = np.cumsum(np.repeat([2,3],[periodicSegments.shape[0], periodicTriangles.shape[0]]),dtype=np.int32)
        periodicEntityIds = np.array(np.hstack([periodicSegments["entity_id"], periodicTriangles["entity_id"]]))
        entityTransformation = np.array(periodicEntity["periodic_transformation"])
        
Michel Henry's avatar
Michel Henry committed
426
        VTK.write(odir+"/periodicBoundaries", i, t, (periodicTypes, periodicConnectivity, periodicOffsets),
Michel Henry's avatar
Michel Henry committed
427
                    periodicX, cell_data=[("Entity_id", periodicEntityIds.reshape([-1,1]))], field_data=[(b"Entity_transformation", entityTransformation)])
Michel Henry's avatar
Michel Henry committed
428

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
429
        ##Contacts
Nathan Coppin's avatar
Nathan Coppin committed
430
        fdata = []
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
431
        contacts = self.contacts()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
432
        for tid,name in enumerate(("particle_particle","particle_segment","particle_triangle")) :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
433
            o = contacts[contacts["type"]==tid]["o"]
434
            v = contacts[contacts["type"]==tid]["reaction"]
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
435
            nameb = name.encode()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
436
            fdata.append((nameb, v))
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
437
            fdata.append((nameb+b"_idx",o))
Nathan Coppin's avatar
Nathan Coppin committed
438
439
        x = np.array([[0.,0.,0.]])
        VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")
440
        VTK.write_multipart(odir+"/particle_problem",["particles_%05d.vtp"%i,"boundaries_%05d.vtu"%i,"periodicBoundaries_%05d.vtu"%i,"contacts_%05d.vtp"%i],i,t)
441

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
442
    def read_vtk(self,dirname,i) :
Nathan Coppin's avatar
Nathan Coppin committed
443
444
445
446
447
448
        """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
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
449
450
451
452
        x,_,data,_,_ = VTK.read(dirname+"/bodies_%05d.vtp"%i)
        self._lib.particle_problem_allocate_bodies(self._p, c_size_t(x.shape[0]))
        bodies = self._get_array("body",self._bodytype)
        bodies["theta"][:] = data["Theta"][:,0]
453
454
        bodies["invertM"][:] = data["iMass"][:,0]
        bodies["invertI"][:] = data["iInertia"][:,0]
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
455
456
457
458
459
460
        self.body_position()[:] = x[:,:self._dim]
        self.body_velocity()[:] = data["Velocity"][:,:self._dim]
        self.body_omega()[:] = data["Omega"]
        x,_,data,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
        self.id2material = list(s.decode() for s in VTK.string_array_decode(fdata["MaterialNames"]))
        self.material2id = dict((name,i) for i,name in enumerate(self.id2material))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
461
462
        self._lib.particle_problem_allocate_material(self._p, c_int(len(self.material2id)))
        self.mu()[:] = fdata["FrictionCoefficients"].reshape(self.mu().shape)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
463
464
465
466
467
468
        self._lib.particle_problem_allocate_particles(self._p, c_size_t(x.shape[0]))
        particles = self._get_array("particle",self._particletype)
        particles["r"][:] = data["Radius"][:,0]
        particles["body"][:] = data["Body"][:,0]
        particles["x"][:] = data["RelativePosition"][:,:self._dim]
        particles["material"][:] = data["Material"][:,0]
Nathan Coppin's avatar
Nathan Coppin committed
469
        x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
470
471
        self.id2tag = list(s.decode() for s in VTK.string_array_decode(fdata["TagNames"]))
        self.tag2id = dict((name,i) for i,name in enumerate(self.id2tag))
Nathan Coppin's avatar
Nathan Coppin committed
472
473
        offsets = np.hstack([[0],cells["offsets"]])
        el = cells["connectivity"]
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
474
475
476
477
        segidx = np.where(cells["types"] == 3)[0]
        self._lib.particle_problem_allocate_segments(self._p, c_size_t(segidx.shape[0]))
        segxidx = np.column_stack([el[offsets[segidx]], el[offsets[segidx]+1]])
        segments = self.segments()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
478
479
        segments["x"][:] = d["RelativePosition"][segxidx,:self._dim]
        segments["body"][:] = cdata["Body"][:,0]
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
480
481
482
483
484
485
486
487
488
        segments["material"] = cdata["Material"][segidx][:,0]
        segments["tag"] = cdata["Tag"][segidx][:,0]
        if self._dim == 3:
            trixidx = np.column_stack([el[offsets[segidx]], el[offsets[segidx]+1],el[offsets[segidx]+3]])
            triangles = self.triangles()
            triangles["p"] = x[el[trixidx],:self._dim]
            triangles["v"] = 0
            triangles["material"] = cdata["Material"][triidx][:,0]
            triangles["tag"] = cdata["Tag"][triidx][:,0]
Michel Henry's avatar
Michel Henry committed
489
490
491
        
        periodicX, periodicCells,_, periodicCellsData, periodicFieldsData = VTK.read(dirname+"/periodicBoundaries_%05d.vtu"%i)
        entityTransformation = periodicFieldsData["Entity_transformation"]
492
        entity_ids = periodicCellsData["Entity_id"].ravel()
Michel Henry's avatar
Michel Henry committed
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
        periodicOffsets = np.hstack([[0], periodicCells["offsets"]])
        periodicEl = periodicCells["connectivity"]
        # add entity
        for entity_id in range(np.shape(entityTransformation)[0]):
            trans = tuple(entityTransformation[entity_id])
            self.add_boundary_periodic_entity(entity_id, self._dim-1, trans)
        # add periodic segment
        for idx in np.where(periodicCells["types"] == 3)[0]:
            pidx = periodicEl[periodicOffsets[idx]]
            entity_id = entity_ids[pidx//2]
            x0 = periodicX[pidx,:self._dim]
            x1 = periodicX[pidx+1,:self._dim]
            self.add_boundary_periodic_segment(tuple(x0), tuple(x1),entity_id)
        # add periodic triangles
        for idx in np.where(periodicCells["types"] == 5)[0]:
            pidx = periodicEl[periodicOffsets[idx]]
            entity_id = entity_ids[pidx//3]
            x0 = periodicX[pidx,:self._dim]
            x1 = periodicX[pidx+1,:self._dim]
            x2 = periodicX[pidx+2,:self._dim]
            self.add_boundary_periodic_triangle(tuple(x0), tuple(x1), tuple(x2), entity_id)

Nathan Coppin's avatar
Nathan Coppin committed
515
        _,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
516
        ks = ("particle_particle","particle_segment","particle_triangle")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
517
        v = np.vstack(list([fdata[k] for k in ks]))
Nathan Coppin's avatar
Nathan Coppin committed
518
        o = np.vstack([fdata[k+"_idx"] for k in ks])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
519
        t = np.repeat([0,2,3],[fdata[k].shape[0] for k in ks])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
520
521
522
523
        self._lib.particle_problem_allocate_contacts(self._p, c_size_t(t.size))
        self.contacts()["reaction"] = v
        self.contacts()["o"] = o
        self.contacts()["type"] = t
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
524

525
    def add_boundary_periodic_entity(self, etag, edim, transform):
Michel Henry's avatar
Michel Henry committed
526
527
528
529
530
531
532
        """Adds a periodic entity.

        Keyword arguments:
        etag -- tag of the entity
        edim -- dimension of the entity
        transform -- tuple of the transformation to applied to the periodic entity
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
533
534
        self._lib.particle_problem_add_boundary_periodic_entity.restype = c_size_t
        return self._lib.particle_problem_add_boundary_periodic_entity(self._p, c_int(etag), c_int(edim), self._coord_type(*transform))
Michel Henry's avatar
Michel Henry committed
535

536
    def add_boundary_periodic_segment(self, x0, x1, etag) :
Michel Henry's avatar
Michel Henry committed
537
538
539
540
541
542
543
        """Adds a boundary periodic segment.

        Keyword arguments:
        x0 -- Tuple of the coordinates of the first endpoint
        x1 -- Tuple of the coordinates of the second endpoint
        tag -- entity tag
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
544
545
        self._lib.particle_problem_add_boundary_periodic_segment.restype = c_size_t
        return self._lib.particle_problem_add_boundary_periodic_segment(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(etag))
Michel Henry's avatar
Michel Henry committed
546

547
    def add_boundary_segment(self, x0, x1, tag, material="default") :
Nathan Coppin's avatar
Nathan Coppin committed
548
549
550
551
552
553
554
555
        """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
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
556
557
558
559
        imat = self.get_material_id(material)
        itag = self.get_tag_id(material)
        self._lib.particle_problem_add_boundary_segment.restype = c_size_t
        return self._lib.particle_problem_add_boundary_segment(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(itag), c_int(imat))
560

561
    def add_boundary_triangle(self, x0, x1, x2, tag, material="default") :
Nathan Coppin's avatar
Nathan Coppin committed
562
563
564
565
566
567
568
569
570
        """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
        """
571
        if self._dim != 3 :
572
            raise NameError("Triangle boundaries only available in 3D")
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
573
574
575
576
        imat = self.get_material_id(material)
        itag = self.get_tag_id(material)
        self._lib.particle_problem_add_boundary_triangle.restype = c_size_t
        return self._lib.particle_problem_add_boundary_triangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), c_int(itag), c_int(imat))
Michel Henry's avatar
Michel Henry committed
577

Michel Henry's avatar
Michel Henry committed
578
    def add_boundary_periodic_triangle(self, x0, x1, x2, etag):
Michel Henry's avatar
Michel Henry committed
579
580
581
582
583
584
585
586
        """Adds a boundary periodic 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 
        etag -- tag of the periodic entity 
        """
Michel Henry's avatar
Michel Henry committed
587
588
        if self._dim != 3 :
            raise NameError("Triangle boundaries only available in 3D")
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
589
590
        self._lib.particle_problem_add_boundary_periodic_triangle.restype = c_size_t
        return self._lib.particle_problem_add_boundary_periodic_triangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), c_int(etag))
Nathan Coppin's avatar
Nathan Coppin committed
591

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
592
593
594
    def add_body(self, x, invertM, invertI):
        self._lib.particle_problem_add_body.restype = c_size_t
        return self._lib.particle_problem_add_body(self._p, self._coord_type(*x), c_double(invertM), c_double(invertI))
Nathan Coppin's avatar
Nathan Coppin committed
595

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
596
    def add_particle_body(self, positions,radii,masses,material="default",forced=False):
Nathan Coppin's avatar
Nathan Coppin committed
597
        """Adds a body in the particle problem.
598
599

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
600
601
602
603
604
        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
605
        """
Nathan Coppin's avatar
Nathan Coppin committed
606
607
608
609
        if positions.shape[0] != radii.shape[0] or positions.shape[0]!=masses.shape[0] or masses.shape[0]!=radii.shape[0]: 
          raise NameError("Wrong lengths of particle data vectors.")

        x = np.sum(positions*masses.reshape((-1,1)),axis=0)/np.sum(masses)
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
610
        I = np.sum(masses*radii**2/2 + np.sum((positions-x)**2,axis=1)*masses)
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
611
612
        invertI = 0 if forced else 1/I
        invertM = 0 if forced else 1/np.sum(masses)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
613
        idd = self.add_body(x, invertM, invertI)
614
615
        if not isinstance(material,(list, np.ndarray)):
          material = [material for i in range(positions.shape[0])]
Nathan Coppin's avatar
Nathan Coppin committed
616
        for i in range(positions.shape[0]):
617
          self.add_particle_to_body(positions[i,:]-x,radii[i],idd, material[i])
Nathan Coppin's avatar
Nathan Coppin committed
618

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
619
    def add_particle_to_body(self, x, r, body, material="default"):
Nathan Coppin's avatar
Nathan Coppin committed
620
        """Adds a particle in the particle problem.
621
622

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
623
624
        x -- Tuple to set the centre particle position
        r -- Particle radius
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
625
        material -- Particle material
626
        """
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
627
        imat = self.get_material_id(material)
628
629
        if material.lower() == "ignore":
          imat = -1
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
630
        self._lib.particle_problem_add_particle(self._p, self._coord_type(*x), c_double(r), c_size_t(body), c_int(imat))
631

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
632
633
634
635
636
637
638
639
640
641
642
643
644
645
    def add_particle(self, x, r, m, tag="default",forced=False, inverse_inertia=None):
        """Add a particle in the particle container.

        Keyword arguments:
        x -- tuple to set the centre particle position
        r -- particle radius
        m -- particle mass
        tag -- particle material
        """
        if (not forced) and inverse_inertia is None:
            inverse_inertia = 2/(m*r**2) if self._dim == 2 else 5/(2*m*r**2)
        body = self.add_body(x, 0 if forced else 1/m, 0 if forced else inverse_inertia)
        self.add_particle_to_body([0]*self._dim, r, body, tag)

Nathan Coppin's avatar
Nathan Coppin committed
646
    
647
    def set_use_queue(self, use_queue=1):
Nathan Coppin's avatar
Nathan Coppin committed
648
      """Enables the use of the queue if 1 and disables it if 0."""
649
      self._lib.particle_problem_set_use_queue(self._p, c_int(use_queue))
650

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
651
    def _save_contacts(self):
Nathan Coppin's avatar
Nathan Coppin committed
652
      """Saves the contacts from the current time step."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
653
      self._saved_contacts = np.copy(self.contacts())
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
654
655

    def _restore_contacts(self):
Nathan Coppin's avatar
Nathan Coppin committed
656
      """Restores the saved contacts from the previous time step."""
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
657
658
659
      self._lib.particle_problem.allocate_contacts(self._p,c_size_t(self.saved_contact.shape[0]))
      self.contacts()[:] = self._saved_contacts
      del self._saved_contacts
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
660

Nathan Coppin's avatar
Nathan Coppin committed
661
    def remove_bodies_flag(self,flag) :
Nathan Coppin's avatar
Nathan Coppin committed
662
      """Removes particles based on given flag."""
Nathan Coppin's avatar
Nathan Coppin committed
663
      if flag.shape != (self.n_bodies(),) :
664
          raise NameError("size of flag array should be the number of particles")
Nathan Coppin's avatar
Nathan Coppin committed
665
      self._lib.particle_problem_remove_bodies(self._p, _np2c(flag,np.int32))
666

667
668

    def load_msh_boundaries(self, filename, tags=None, shift=None,material="default") :
Nathan Coppin's avatar
Nathan Coppin committed
669
        """Loads the numerical domain and set the physical boundaries the particles cannot go through.
670
671
        
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
672
673
674
675
        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
676
        """
677
678
679
        if not os.path.isfile(filename):
          print("Error : no such file as " + filename)
          exit(1)
680
        if shift is None :
681
            shift = [0]*self._dim
682
683
684
685
686
687
        shift = np.array(shift)
        gmsh.model.add("tmp")
        gmsh.open(filename)
        gmsh.model.mesh.renumberNodes()
        addv = set()
        adds = set()
688
        periodic_entities = set()
689
        _, x, _ = gmsh.model.mesh.getNodes()
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
690
        x = x.reshape([-1,3])[:,:self._dim]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
691
        self._lib.particle_problem_add_body.restype = c_size_t
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
692
        bndbody = self._lib.particle_problem_add_body(self._p, self._coord_type(*((0,)*self._dim)), c_double(0), c_double(0))
693
694
695
696
697
698
699
700
701
702

        def get_entity_name(edim, etag) :
            for tag in  gmsh.model.getPhysicalGroupsForEntity(edim, etag):
                name = gmsh.model.getPhysicalName(edim,tag)
                if name is not None :
                    return name
            return "none"

        def add_disk(t, stag) :
            if t in addv : return
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
703
            self.add_particle(x[t]+shift, 0, bndbody, material, forced=True)
704
705
706
707
708
709
710
711
712
713
714
715
716
            addv.add(t)

        def add_segment(t0, t1, stag) :
            key = (min(t0,t1),max(t0,t1))
            if key in adds : return
            adds.add(key)
            self.add_boundary_segment(x[t0]+shift, x[t1]+shift,
                                      stag, material)

        def add_triangle(t0, t1, t2, stag):
            self.add_boundary_triangle(x[t0]+shift, x[t1]+shift,x[t2]+shift,
                                        stag, material)

717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
        for dim, tag in gmsh.model.getEntities(self._dim-1) :
            ptag, cnodes, pnodes, _ = gmsh.model.mesh.getPeriodicNodes(dim, tag)
            if ptag == tag or len(cnodes) == 0: continue
            periodic_entities.add((dim,tag))
            periodic_entities.add((dim,ptag))
            trans = x[(pnodes-1)[0]] - x[(cnodes-1)[0]]
            self.add_boundary_periodic_entity(tag,dim,trans)
            self.add_boundary_periodic_entity(ptag,dim,-trans)
            pmap = dict((int(i)-1,int(j)-1) for i,j in zip (cnodes,pnodes))
            for etype, _, enodes in zip(*gmsh.model.mesh.getElements(dim,tag)):
                if dim == 1 :
                    if etype != 1 : continue
                    for l in (enodes-1).reshape([-1,2]):
                        self.add_boundary_periodic_segment(
                                x[l[0]]+shift, x[l[1]]+shift, tag)
                        self.add_boundary_periodic_segment(
                                x[pmap[l[0]]]+shift, x[pmap[l[1]]]+shift, ptag)
                else :
                    if etype != 2 : continue
                    for l in (enodes-1).reshape([-1,3]):
                        self.add_boundary_periodic_triangle(
                                x[l[0]]+shift, x[l[1]]+shift,
                                x[l[2]]+shift, tag)
                        self.add_boundary_periodic_triangle(
                                x[pmap[l[0]]]+shift, x[pmap[l[1]]]+shift,
                                x[pmap[l[2]]]+shift, ptag)
743
744
        for dim, tag in gmsh.model.getEntities(0) :
            stag = get_entity_name(dim,tag)
745
            if (dim,tag) in periodic_entities : continue
746
747
748
749
750
751
752
753
            if not (tags is None or stag in tags) : continue
            for etype, _, enodes in zip(*gmsh.model.mesh.getElements(dim,tag)):
                if etype != 15 : continue
                for t in enodes-1 :
                    add_disk(t, stag)

        for dim, tag in gmsh.model.getEntities(1) :
            stag = get_entity_name(dim,tag)
754
            if (dim,tag) in periodic_entities : continue
755
756
757
758
759
760
761
762
763
764
765
            if not (tags is None or stag in tags) : continue
            for etype, _, enodes in zip(*gmsh.model.mesh.getElements(dim,tag)):
                if etype != 1 : continue
                for t in enodes-1 :
                    add_disk(t, stag)
                for l in (enodes-1).reshape([-1,2]):
                    add_segment(l[0], l[1], stag)

        if self._dim == 3 :
            for dim, tag in gmsh.model.getEntities(2) :
                stag = get_entity_name(dim,tag)
766
                if (dim,tag) in periodic_entities : continue
767
768
769
770
771
772
773
774
775
776
777
778
779
                if not (tags is None or stag in tags) : continue
                for etype, _, enodes in zip(*gmsh.model.mesh.getElements(dim,tag)):
                    if etype != 2 : continue
                    for t in enodes-1 :
                        add_disk(t, stag)
                    for v0,v1,v2 in (enodes-1).reshape([-1,3]):
                        add_segment(v0, v1, stag)
                        add_segment(v1, v2, stag)
                        add_segment(v2, v0, stag)
                        add_triangle(v0, v1, v2, stag)

        gmsh.model.remove()
            
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
780
    def dim(self) :
781
        """Returns the dimension of the particle problem"""
782
        return self._dim