fluid3.py 3.43 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
2
3
4
5
from ctypes import *
import numpy
import signal
import os

Matthieu Constant's avatar
Matthieu Constant committed
6

7
dir_path = os.path.dirname(os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
lib = CDLL(os.path.join(dir_path,"libmbfluid3.so"))
Matthieu Constant's avatar
Matthieu Constant committed
9
10
11
12
13
14
15
16
17
18

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


BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
class StrongBoundary(Structure):
    _fields_ = [("tag",c_char_p),("field",c_int),("apply",BNDCB)]

class fluid_problem :

Matthieu Constant's avatar
Matthieu Constant committed
19
    def __init__(self, mesh_file_name, g, mu, rho, alpha, strong_boundaries, notComputeEpsilon):
Matthieu Constant's avatar
Matthieu Constant committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
        nsb = len(strong_boundaries)
        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])
        self.asb = asb
        lib.fluid_problem_new.restype = c_void_p
Matthieu Constant's avatar
Matthieu Constant committed
35
        self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(alpha), nsb, asb, c_int(notComputeEpsilon)))
Matthieu Constant's avatar
Matthieu Constant committed
36
37
38
39
40
41
42
        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
43
44
    def adapt_mesh(self, lcmax, lcmin, n_el, alpha, notComputeEpsilon) :
        lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(alpha), c_int(notComputeEpsilon))
Matthieu Constant's avatar
Matthieu Constant committed
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))

Matthieu Constant's avatar
visu 3d    
Matthieu Constant committed
49
50
    def export_vtk(self, output_dir, t, it) :
        lib.fluid_problem_export_vtk(self._fp, output_dir.encode(), c_double(t), c_int(it))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
52
53

    def import_vtk(self, filename) :
        lib.fluid_problem_import_vtk(self._fp, filename.encode())
Matthieu Constant's avatar
visu 3d    
Matthieu Constant committed
54
    
Matthieu Constant's avatar
Matthieu Constant committed
55
56
57
58
59
60
61
62
    def compute_node_force(self, dt) :
        def np2c(a) :
            return c_void_p(a.ctypes.data)
        forces = numpy.ndarray([self.n_particles,3],"d")
        lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), np2c(forces))
        return forces

    def implicit_euler(self, dt) :
63
        return lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
Matthieu Constant's avatar
Matthieu Constant committed
64

65
    def set_particles(self, mass, volume, position, velocity, elid=None):
Matthieu Constant's avatar
Matthieu Constant committed
66
67
        def np2c(a) :
            return c_void_p(numpy.ascontiguousarray(a).ctypes.data)
68
69
        def np2ci(a) :
            return c_void_p(numpy.require(a,"int","C").ctypes.data)
Matthieu Constant's avatar
Matthieu Constant committed
70
        self.n_particles = mass.shape[0]
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
        if elid is not None:
            lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),np2ci(elid))
        else :
            lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None)

    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)