fluid.py

# Marblesbag - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (Marblesbag) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# 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
# GNU Lesser General Public License for more details.
# 
# 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, 
# see <http://www.gnu.org/licenses/>.

from ctypes import *
import numpy
import signal
import os

dir_path = os.path.dirname(os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
lib = CDLL(os.path.join(dir_path,"libmbfluid2.so"))

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


BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))

class fluid_problem :

    def __init__(self, g, mu, rho):
        self.strong_cb_ref = []
        def np2c(a) :
            tmp = numpy.require(a,"float","C")
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
        lib.fluid_problem_new.restype = c_void_p
        n_fluids = numpy.require(rho,"float","C").reshape([-1]).shape[0]
        self._fp = c_void_p(lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho)))
        if self._fp == None :
            raise NameError("cannot create fluid problem")

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

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

    def set_weak_boundary(self, tag, bndtype) :
        lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bndtype.encode())

    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)

    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))

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

    def export_vtk(self, output_dir, t, it) :
        lib.fluid_problem_export_vtk(self._fp, output_dir.encode(), c_double(t), c_int(it))

    def import_vtk(self, filename) :
        lib.fluid_problem_import_vtk(self._fp, filename.encode())

    def compute_node_force(self, dt) :
        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))
        return forces

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

    def set_particles(self, mass, volume, position, velocity, reload = False):
        def np2c(a) :
            tmp = numpy.require(a,"float","C")
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
        self.n_particles = mass.shape[0]
        lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None,c_int(reload))

    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);

    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)