scontact.py 20 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not, 
# see <http://www.gnu.org/licenses/>.

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

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

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

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

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

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

is_64bits = sys.maxsize > 2**32

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

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

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

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

        Keyword arguments:
        dim -- dimension of the problem

        Raises:
        ValueError -- if dimension differs from 2 or 3
        NameError -- if C builder cannot be found
        """
85
        self._dim = dim
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
86
        self._friction_enabled = friction_enabled
87
        if dim == 2 :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
88
            self._lib = lib2f if friction_enabled else lib2
89
90
91
92
93
94
95
96
            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())
97
        self._fc = None
98
99
100
        if self._p == None :
            raise NameError("Cannot create particle problem.")

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

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

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

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

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

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

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

136
137
138
139
    def stressTensor(self):
        """Return the stress tensor for post visualization."""
        return self._get_matrix("StressTensor",self._dim*self._dim)

Nathan Coppin's avatar
Nathan Coppin committed
140
141
142
    def disk_tag(self):
        return self._get_idx("DiskTag")

143
144
145
146
    def disk_material(self):
        return self._get_idx("DiskMaterial")


147
    def segments(self):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
        return self._get_matrix("Segment",3*self._dim+(1 if self._friction_enabled else 0))
149
150
151
152

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

153
154
155
    def segment_material(self):
        return self._get_idx("SegmentMaterial")

156
    def triangles(self):
Nathan Coppin's avatar
Nathan Coppin committed
157
158
159
160
161
162
163
164
165
166
        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)
167

168
169
170
171
172
    def triangle_material(self):
        if self._dim == 3 :
            return self._get_idx("TriangleMaterial")
        else :
            return np.array([],np.int32)
173
    
174
    def contact_forces(self):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
175
176
        if self._fc is None:
            self._fc = np.zeros_like(self.velocity())
177
        return(self._fc)
178

179
180
181
    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
182
        disks = self._get_matrix("Disk",2*self._dim+1)
183
        diskTag = self._get_idx("DiskTag")
Nathan Coppin's avatar
Nathan Coppin committed
184
185
        disks[diskTag == tag, self._dim:2*self._dim] = v
        segments = self._get_matrix("Segment",3*self._dim)
186
        segmentTag = self._get_idx("SegmentTag")
Nathan Coppin's avatar
Nathan Coppin committed
187
        segments[segmentTag == tag, 2*self._dim:3*self._dim] = v
188
189
190
191
192
        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
193
    def set_tangent_boundary_velocity(self, tag, v):
194
195
196
197
198
199
200
        """ Sets the tangent velocity of the frictional boundary.
            Only in 2D when friction is enabled.

        Keyword arguments:
        tag -- Name of the boundary
        v -- velocity
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
201
        assert self._friction_enabled
Nathan Coppin's avatar
Nathan Coppin committed
202
        self._lib.particleProblemSetBoundaryTangentVelocity(self._p,tag.encode(),c_double(v))
203
204
    
    def set_friction_relaxation(self, relaxation):
205
206
207
       """ Sets the relaxation parameter for the resolution of the frictional problem.
            Only in 2D when friction is enabled.

208
209
210
211
212
       Keyword arguments:
       relaxation -- relaxation parameter, between 0 and 1
       """
       assert self._friction_enabled
       self._lib.particleProblemSetFrictionRelaxation(self._p,c_double(relaxation))
Nathan Coppin's avatar
Nathan Coppin committed
213

214
    def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
215
216
217
218
219
220
221
222
        """ Sets the friction coefficient between two materials.
            Only in 2D when friction is enabled.

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

226
    def iterate(self, dt, forces,tol=1e-8) :
227
228
229
230
231
232
233
        """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
        """ 
234
        vn = self.velocity().copy()
235
        self._get_matrix("ExternalForces",self._dim)[:] = forces
Nathan Coppin's avatar
Nathan Coppin committed
236
        self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
237
        self._fc = (self.velocity()-vn)*self.mass()/dt - forces
238

239
240
241
242
243
244
    def computeStressTensor(self, nodes, radius):
        """Compute the stres tensor of the granular material"""
        n_nodes = len(nodes[:,0])
        self._lib.particleProblemComputeStressTensor(self._p,_np2c(nodes[:,:self._dim]),c_double(radius),c_int(n_nodes))

    def write_vtk(self, odir, i, t, nodes=None, elements=None) :
245
        """Write output files for post-visualization."""
246

Nathan Coppin's avatar
Nathan Coppin committed
247
        v = self.velocity()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
248
249
250
251
        if self._friction_enabled :
            omega = self.omega()
        else :
            omega = np.zeros((v.shape[0],1))
Nathan Coppin's avatar
Nathan Coppin committed
252
        x = self.position()
253
        material = self._get_idx("ParticleMaterial").reshape(-1,1)
Nathan Coppin's avatar
Nathan Coppin committed
254
255
256
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
            x = np.insert(x,self._dim,0,axis=1)
257
        data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Material",material),("Contacts",self.contact_forces())]
258
259
260
261
        nmat = self._lib.particleProblemNMaterial(self._p)
        self._lib.particleProblemGetMaterialTagName.restype = c_char_p
        tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
        VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags))])
Nathan Coppin's avatar
Nathan Coppin committed
262

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
263
        xdisk = self._get_matrix("Disk",2*self._dim+(2 if self._friction_enabled else 1))[:,:self._dim]
Nathan Coppin's avatar
Nathan Coppin committed
264
265
266
267
268
269
270
271
272
273
274
        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)
275
276
277
278
279
280
        tags = np.array(np.hstack([self.disk_tag(),self.segment_tag(),self.triangle_tag()]))
        materials = np.array(np.hstack([self.disk_material(),self.segment_material(),self.triangle_material()]))
        ntagname = self._lib.particleProblemNTag(self._p)
        self._lib.particleProblemGetTagName.restype = c_char_p
        tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)])
        VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1]))],field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
Nathan Coppin's avatar
Nathan Coppin committed
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295

        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)
296
297
298
299
300
301
302
303
304
305
306

        if nodes is not None:
            stressTensor = self.stressTensor()
            data = [("stressTensor",stressTensor)]
            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)])
            connectivity = elements
            types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
            offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0]))
            VTK.write(odir+"/post",i,t,(types,connectivity,offsets),nodes,data,vtktype="vtu")
Nathan Coppin's avatar
Nathan Coppin committed
307
        
308
    def read_vtk(self,dirname,i,contact_factor=1) :
309
        """Read an existing output file to set the particle data."""
Nathan Coppin's avatar
Nathan Coppin committed
310
        x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
311
312
313
314
315
        matnames =  VTK.string_array_decode(fdata["MaterialNames"])
        matmap = {}
        for j,n in enumerate(matnames) :
            matmap[j] = self._lib.particleProblemGetMaterialTagId(self._p,n)
        partmat =np.vectorize(matmap.get)(d["Material"])
Nathan Coppin's avatar
Nathan Coppin committed
316
        self._lib.particleProblemClearParticles(self._p)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
317

Nathan Coppin's avatar
Nathan Coppin committed
318
        self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
319
                _np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(partmat,dtype=np.int32))
320
        self._fc = d["Contacts"][:,:self._dim]*contact_factor
Nathan Coppin's avatar
Nathan Coppin committed
321
322

        x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
323
324
325
326
        tagnames =  VTK.string_array_decode(fdata["TagNames"])
        tagmap = {}
        for j,n in enumerate(tagnames) :
            tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
Nathan Coppin's avatar
Nathan Coppin committed
327
328
329
        self._lib.particleProblemClearBoundaries(self._p)
        offsets = np.hstack([[0],cells["offsets"]])
        el = cells["connectivity"]
330
331
        tags = np.vectorize(tagmap.get)(cdata["Tag"])
        materials = cdata["Material"]
Nathan Coppin's avatar
Nathan Coppin committed
332
333
        for idx in np.where(cells["types"] == 1)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
334
            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
335
336
337
        for idx in np.where(cells["types"] == 3)[0] :
            x0 = x[el[offsets[idx]],:self._dim]
            x1 = x[el[offsets[idx]+1],:self._dim]
338
            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
339
340
341
342
        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]
343
            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
344
345
346
347
        _,_,_,_,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])
348
        v = np.hstack([v,vt])*contact_factor
Nathan Coppin's avatar
Nathan Coppin committed
349
350
351
        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))
352

353
    def add_boundary_disk(self, x0, r, tag, material="default") :
354
        self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
355
        return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode(),material.encode())
356

357
    def add_boundary_segment(self, x0, x1, tag, material="default") :
358
        self._lib.particleProblemAddBoundarySegment.restype = c_size_t
359
        return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode(),material.encode())
360

361
    def add_boundary_triangle(self, x0, x1, x2, tag, material="default") :
362
        if self._dim != 3 :
363
364
            raise NameError("Triangle boundaries only available in 3D")
        self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t
365
        return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), tag.encode(),material.encode())
366

367
    def add_particle(self, x, r, m, tag="default"):
368
369
370
371
372
373
        """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
374
        tag -- particle tag
375
        """
376
        self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), tag.encode())
377

378
379
380
381
382
    def set_use_queue(self, use_queue=1):
      """Enables the use of the queue if 1 and disables it if 0."""
      
      self._lib.particleProblemSetUseQueue(self._p, c_int(use_queue))

383
    def load_msh_boundaries(self, filename, tags, shift=None,material="default") :
384
385
386
387
388
389
        """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
390
        material -- material of the boundary
391
        """
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
        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])
408
                        self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag,material)
409
410
                    v0 = el.vertices[1]
                    v1 = el.vertices[0]
411
                    self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag,material)
412
413
414
415
416
417
418
419
420
421
422
423
424
        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])
425
                            self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag,material)
426
427
428
429
430
431
                        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]),
432
                                stag,material)
433
434
435
436
437
438
439
440
                            
                    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]),
441
                          stag,material)