fluid.py 5.13 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
35
36

class fluid_problem :

37
38
    def __init__(self, g, mu, rho):
        self.strong_cb_ref = []
Matthieu Constant's avatar
Matthieu Constant committed
39
40
41
42
43
        def np2c(a) :
            tmp = numpy.require(a,"float","C")
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
44
        lib.fluid_problem_new.restype = c_void_p
45
        n_fluids = numpy.require(rho,"float","C").reshape([-1]).shape[0]
46
        self._fp = c_void_p(lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho)))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
47
48
49
50
51
52
53
        if self._fp == None :
            raise NameError("cannot create fluid problem")

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

54
55
56
    def load_msh(self, mesh_file_name) :
        lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
57
58
59
    def set_weak_boundary(self, tag, bndtype) :
        lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bndtype.encode())

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
    def set_strong_boundary(self, tag, field,callback_or_value, row=None) :
        if row is None :
            row = field
        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
        bndcb = BNDCB(Bnd(callback_or_value).apply)
        self.strong_cb_ref.append(bndcb)
        lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), field, row, bndcb)

Matthieu Constant's avatar
Matthieu Constant committed
77
78
    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
79
80
81
82

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

83
84
85
    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
86
87
88
    def import_vtk(self, filename) :
        lib.fluid_problem_import_vtk(self._fp, filename.encode())

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
89
    def compute_node_force(self, dt) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
90
91
        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
92
93
94
        return forces

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

Matthieu Constant's avatar
Matthieu Constant committed
97
    def set_particles(self, mass, volume, position, velocity, reload = False):
98
        def np2c(a) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
99
100
101
102
            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
103
        self.n_particles = mass.shape[0]
Matthieu Constant's avatar
Matthieu Constant committed
104
        lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None,c_int(reload))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
105

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    def _get_matrix(self, f_name, nrow, ncol) :
        f = getattr(lib,"fluid_problem_"+f_name)
        f.restype = POINTER(c_double)
        return numpy.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))

    def solution(self) :
        return self._get_matrix("solution", self.n_nodes(), self.n_fields())

    def coordinates(self) :
        return self._get_matrix("coordinates",self.n_nodes(), 3)
    
    def n_fields(self):
        lib.fluid_problem_n_fields.restype = c_int
        return lib.fluid_problem_n_fields(self._fp);

    def n_nodes(self):
        lib.fluid_problem_n_nodes.restype = c_int
        return lib.fluid_problem_n_nodes(self._fp);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
124

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
125
126
127
128
129
130
131
132
133
134
    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)