scontact2.py 4.08 KB
Newer Older
1
from __future__ import division
2
3
from . import mesh
import shutil
4
5
6
import numpy
import os
import sys
7
8
9
10
from ctypes import *
import signal

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

12
dir_path = os.path.dirname(os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
13
14
15
16
17
18
19
lib = CDLL(os.path.join(dir_path,"libscontact2.so"))

coord_type = c_double*2
is_64bits = sys.maxsize > 2**32

class ParticleProblem :

20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
    def _get_matrix(self, fName, ncol) :
        f = getattr(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)
        return numpy.ctypeslib.as_array(buf).reshape([-1,ncol])
    
    def _get_idx(self, fName) :
        f = getattr(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)
        return numpy.ctypeslib.as_array(buf)

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
36
    def __init__(self) :
37
        lib.particleProblemNew.restype = c_void_p
38
39
        self._p = c_void_p(lib.particleProblemNew())
        if self._p == None :
40
41
            raise NameError("cannot create particle problem")

42
43
    def volume(self):
        return numpy.pi * self._get_matrix("Particle", 2)[:, [0]]**2
44

45
46
    def r(self) :
        return self._get_matrix("Particle", 2)[:, 0]
47

48
49
    def mass(self):
        return self._get_matrix("Particle", 2)[:, [1]]
50

51
52
    def position(self):
        return self._get_matrix("Position", 2)
53

54
55
56
    def velocity(self):
        return self._get_matrix("Velocity", 2)

Matthieu Constant's avatar
Matthieu Constant committed
57
58
59
60
61
62
    def segments(self):
        return self._get_matrix("Segment",6)

    def segmentsTag(self):
        return self._get_idx("SegmentTag")
    
63
64
65
66
67
68
69
70
71
72
73
74
    def set_boundary_velocity(self, tag, v) :
        lib.particleProblemGetTagId.restype = c_size_t
        tag = lib.particleProblemGetTagId(self._p, tag.encode())
        disks = self._get_matrix("Disk",5)
        diskTag = self._get_idx("DiskTag")
        disks[diskTag == tag, 2:4] = v
        segments = self._get_matrix("Segment",6)
        segmentTag = self._get_idx("SegmentTag")
        segments[segmentTag == tag, 4:6] = v

    def iterate(self, dt, forces) :
        self._get_matrix("ExternalForces",2)[:] = forces
Matthieu Constant's avatar
Matthieu Constant committed
75
        lib.particleProblemIterate(self._p, c_double(numpy.min(self.r())), c_double(dt), c_double(1e-8), c_int(-1))
76
77

    def write_vtk(self, odir, i, t) :
78
        lib.particleProblemWriteVtk(self._p, odir.encode(), i)
79

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
81
82
    def read_vtk(self,dirname,i) :
        lib.particleProblemReadVtk(self._p, dirname.encode(),i)

83
    def add_boundary_disk(self, x0, r, tag) :
84
        lib.particleProblemAddBoundaryDisk.restype = c_size_t
85
        return lib.particleProblemAddBoundaryDisk(self._p, coord_type(*x0), c_double(r), tag.encode())
86

87
    def add_boundary_segment(self, x0, x1, tag) :
88
        lib.particleProblemAddBoundarySegment.restype = c_size_t
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
        return lib.particleProblemAddBoundarySegment(self._p, coord_type(*x0), coord_type(*x1), tag.encode())

    def add_particle(self, x, r, m):
        lib.particleProblemAddParticle(self._p, coord_type(*x), c_double(r), c_double(m))

    def load_msh_boundaries(self, filename, tags, shift=[0,0]) :
        m = mesh.Mesh(filename)
        self._addv = set()
        pnum = [m.getPhysicalNumber(1, i) for i in tags]
        for ent in m.entities :
            if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
                continue
            tag = 0 if not ent.physicals else ent.physicals[0]
            stag = m.getPhysicalName(ent.dimension, tag) or str(tag)
            for i, el in enumerate(ent.elements) :
                for v in el.vertices :
                    if v[3] in self._addv :
                        continue
                    self._addv.add(v[3])
                    self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag)
                v0 = el.vertices[1]
                v1 = el.vertices[0]
                self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag)
112
113