Basilisk source code (http://basilisk.fr/src/)

root / src / common.i

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
%{
  #define SWIG_FILE_WITH_INIT
  #define face
  #define vertex

  typedef int scalar;
  typedef struct {
    scalar x, y;
  } vector;
  typedef struct {
    vector x, y;
  } tensor;

  extern void _init_solver();
%}

%rename(_scalar) scalar;
typedef int scalar;

%rename(_vector) vector;
typedef struct {
  scalar x, y;
} vector;

%rename(_tensor) tensor;
typedef struct {
  vector x, y;
} tensor;

// scalar lists
%typemap(in) scalar * {
  /* Check if is a list */
  if (PyList_Check($input)) {
    int size = PyList_Size($input);
    int i = 0;
    $1 = (scalar *) malloc((size+1)*sizeof(scalar));
    for (i = 0; i < size; i++) {
      PyObject *o = PyList_GetItem($input,i);
      if (PyInt_Check(o))
	$1[i] = PyInt_AsLong(PyList_GetItem($input,i));
      else {
	PyErr_SetString(PyExc_TypeError,"list must contain scalars");
	free($1);
	return NULL;
      }
    }
    $1[i] = -1;
  } else {
    PyErr_SetString(PyExc_TypeError,"not a list");
    return NULL;
  }
}

// This cleans up the scalar * array we malloc'd before the function call
%typemap(freearg) scalar * {
  free((char *) $1);
}

%exception py_scalar_init {
  $action
  if (result) {
     PyErr_SetString(PyExc_RuntimeError, "initialization failed");
     return NULL;
  }
}
%inline %{
  extern void init_grid (int n);
  extern void free_grid (void);
  struct _origin { double x, y, z; };
  extern void origin (struct _origin p);
  extern void size (double L);

  extern int py_scalar_init (scalar s, PyObject * f);
  extern int py_register_event (PyObject * action, PyObject * i, 
				PyObject * t);
%}

%include "numpy.i"

%init %{
  import_array();
  _init_solver();
%}

%pythoncode %{
from numpy import empty_like
class scalar(int):
    def __new__(cls,i=None):
        if i == None:
            i = new_scalar("python")
        return int.__new__(cls,i)
    def __del__(self):
        if callable(_delete):
            _delete([self])
    def __setattr__(self, name, value):
        if name == "f":
            py_scalar_init(self,value)
        else:
            self.__dict__[name] = value
    def f(self,x,y=0):
        try:
            ndim = x.ndim
        except AttributeError:
            return interpolate(self,x,y)
        if ndim == 1:
            return _interpolate1D(self,x,x.size)
        elif ndim == 2:
            z = empty_like(x)
            _interpolate2D(self,x,y,z)
            return z
    def norm(self):
        return normf(self)
    def stats(self):
        return statsf(self)

class vector:
    def __init__(self,v):
        self.x = scalar(v.x)
        self.y = scalar(v.y)

class tensor:
    def __init__(self,t):
        self.x = vector(t.x)
        self.y = vector(t.y)

def event (action, i = None, t = None):
    py_register_event(action, i, t)

from random import uniform
def noise():
    return uniform(-1.,1.)
%}