fluid.py 2.88 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
3
4
5
6
from ctypes import *
import numpy
import signal
import os

dir_path = os.path.dirname(os.path.realpath(__file__))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
7
lib = CDLL(os.path.join(dir_path,"libmbfluid2.so"))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
9
10
11

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


Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
12
BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
13
class StrongBoundary(Structure):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
14
    _fields_ = [("tag",c_char_p),("field",c_int),("apply",BNDCB)]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
15
16
17

class fluid_problem :

18
    def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
19
        nsb = len(strong_boundaries)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
22
23
24
25
26
27
28
29
30
31
        class Bnd :
            def __init__(self, b) :
                self._b = b
            def apply(self, n, xp, vp) :
                x = numpy.frombuffer(cast(xp, POINTER(int(n)*2*c_double)).contents).reshape([n,2])
                v = numpy.frombuffer(cast(vp, POINTER(int(n)*c_double)).contents)
                if callable(self._b) :
                    v[:] = self._b(x)
                else :
                    v[:] = self._b

        asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(Bnd(i[2]).apply)) for i in strong_boundaries])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
32
        self.asb = asb
33
        lib.fluid_problem_new.restype = c_void_p
34
        self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(epsilon), nsb, asb))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
35
36
37
38
39
40
41
        if self._fp == None :
            raise NameError("cannot create fluid problem")

    def __del__(self):
        if(self._fp != None) :
            lib.fluid_problem_free(self._fp)

Matthieu Constant's avatar
Matthieu Constant committed
42
43
    def adapt_mesh(self, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el) :
        lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(gradPmin), c_double(gradPmax), c_double(lcmin), c_double(n_el))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
44
45
46
47
48

    def export(self, output_dir, t, it) :
        lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))

    def compute_node_force(self, dt) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
49
50
        forces = numpy.ndarray([self.n_particles,2],"d",order="C")
        lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), c_void_p(forces.ctypes.data))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
52
53
54
55
        return forces

    def implicit_euler(self, dt) :
        lib.fluid_problem_implicit_euler(self._fp, c_double(dt))

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
56
    def set_particles(self, mass, volume, position, velocity):
57
        def np2c(a) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
58
59
60
61
            tmp = numpy.require(a,"float","C")
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
62
        self.n_particles = mass.shape[0]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
63
        lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
64

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
65

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
66
67
68
69
70
71
72
73
74
75
    def particle_element_id(self) :
        f = getattr(lib,"fluid_problem_particle_element_id")
        f.restype = c_void_p
        fs = getattr(lib,"fluid_problem_n_particles")
        fs.restype = c_int
        ptr = f(self._fp)
        size = fs(self._fp)
        buf = (size*c_int).from_address(ptr)
        return numpy.ctypeslib.as_array(buf)