fluid.py 4.03 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# Marblesbag - Copyright (C) <2010-2018>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
# 	
5
# List of the contributors to the development of Marblesbag: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
7
8
# Description and complete License: see LICENSE file.
# 	
# This program (Marblesbag) is free software: 
9
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
12
13
14
15
# 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
16
# GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
# 
18
19
# 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, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
22
23
24
25
26
from ctypes import *
import numpy
import signal
import os

27
dir_path = os.path.dirname(os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
28
lib = CDLL(os.path.join(dir_path,"libmbfluid2.so"))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
29
30
31
32

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


Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
34
class StrongBoundary(Structure):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
35
    _fields_ = [("tag",c_char_p),("field",c_int),("apply",BNDCB)]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
36
37
38

class fluid_problem :

Matthieu Constant's avatar
Matthieu Constant committed
39
    def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
40
        nsb = len(strong_boundaries)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
41
42
43
44
45
46
47
48
49
50
51
52
        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
53
        self.asb = asb
54
        lib.fluid_problem_new.restype = c_void_p
Matthieu Constant's avatar
Matthieu Constant committed
55
        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
56
57
58
59
60
61
62
        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
63
64
    def adapt_mesh(self, lcmax, lcmin, n_el) :
        lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
65
66
67
68

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

69
70
71
    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
72
73
74
    def import_vtk(self, filename) :
        lib.fluid_problem_import_vtk(self._fp, filename.encode())

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
75
    def compute_node_force(self, dt) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
76
77
        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
78
79
80
        return forces

    def implicit_euler(self, dt) :
81
        return lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
82

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
83
    def set_particles(self, mass, volume, position, velocity):
84
        def np2c(a) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
85
86
87
88
            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
89
        self.n_particles = mass.shape[0]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
90
        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
91

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
92

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
93
94
95
96
97
98
99
100
101
102
    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)