fluid.pyx 3.42 KB
Newer Older
1
import numpy
2
import signal
3
4
from cython.view cimport array as cvarray

5
6
7
8
9
10
11
12
signal.signal(signal.SIGINT, signal.SIG_DFL)

cdef extern from "linear_system.h":
    ctypedef struct StrongBoundary :
        const char *tag;
        int field;
        double v;

13
14
15
16
17
cdef extern from "fluid_problem.h":
    ctypedef struct FluidProblem :
        double mu
        double epsilon 
        double rho
18
19
        int n_particles
        double *particle_force
20
21
22
        pass
    int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double it);
    void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
23
    FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_boundary, StrongBoundary *boundaries);
24
    void fluid_problem_free(FluidProblem *problem);
25
    void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
26
27

cdef extern from "linear_system.h":
28
    void initialize_linear_solver(int argc, char **argv);
29
30

include "openmpihack.pxi"
31
32
33
34
35
36
37
38
import sys
from libc.stdlib cimport malloc, free
cdef char **argv = <char**>malloc(sizeof(char*)*(len(sys.argv)));
argvbyte = [c.encode() for c in sys.argv]
for i,c in enumerate(argvbyte) :
    argv[i] = argvbyte[i]
initialize_linear_solver(len(sys.argv), argv)
free(argv)
39
40
41

cdef class fluid_problem :
    cdef FluidProblem *_fp
42
43
44
45
46
47
48
49
50
    def __cinit__(self, mesh_file_name, mu, rho, epsilon, strong_boundaries):
        cdef StrongBoundary* bnds = <StrongBoundary*> malloc(sizeof(StrongBoundary)*len(strong_boundaries));
        tags = [bnd[0].encode() for bnd in strong_boundaries]
        for i, bnd in enumerate(strong_boundaries) :
            bnds[i].tag = <char*>tags[i]
            bnds[i].field = bnd[1]
            bnds[i].v = bnd[2]
        self._fp = fluid_problem_new(mesh_file_name.encode(), mu, rho, epsilon, len(strong_boundaries),bnds)
        free(bnds)
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
        if self._fp == NULL :
            raise NameError("cannot create fluid problem")
    def __dealloc__(self):
        if(self._fp != NULL) :
            fluid_problem_free(self._fp)
    def __getattr__(self, name) :
        if name == "mu" :
            return self._fp.mu
        elif name == "rho" :
            return self._fp.rho
        elif name == "epsilon" :
            return self._fp.epsilon
    def __setattr__(self, name, val) :
        if name == "mu" :
            self._fp.mu = val
        if name == "rho" :
            self._fp.rho = val
        elif name == "epsilon" :
            self._fp.epsilon = val
    def export(self, output_dir, double t, int it) :
        fluid_problem_export(self._fp, output_dir.encode(), t, it)
    def implicit_euler(self, double dt) :
        fluid_problem_implicit_euler(self._fp, dt)
74
75
76
77
78
        forces = numpy.ndarray([self._fp.n_particles,2],"d")
        for i in range(self._fp.n_particles) :
            forces[i,0] = self._fp.particle_force[i*2+0]
            forces[i,1] = self._fp.particle_force[i*2+1]
        return forces
79

80
81
82
83
84
85
    def set_particles(self, mass, volume, position, velocity):
        cdef double[:,::1] _mass = numpy.ascontiguousarray(mass)
        cdef double[:,::1] _volume = numpy.ascontiguousarray(volume)
        cdef double[:,::1] _position = numpy.ascontiguousarray(position)
        cdef double[:,::1] _velocity = numpy.ascontiguousarray(velocity)
        fluid_problem_set_particles(self._fp,mass.shape[0],&_mass[0][0],&_volume[0][0],&_position[0][0],&_velocity[0][0])