Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Newer
Older
100644 462 lines (430 sloc) 20.432 kb
cf151c41 »
2011-06-20 SRC: sync.
1 #include <Python.h>
2 #include <numpy/ndarrayobject.h>
3
4 void
5 weights_gauss_radau006004 (const double *restrict sigma, int n, int ssi, int ssr,
6 double *restrict omega, int wsi, int wsl, int wsr)
7 {
8 int i;
9 double acc, sigma0, sigma1, sigma2, sigma3, sigma4, sigma5, omega9, omega18, omega15p, omega16m,
10 omega14p, omega3, omega8, omega16p, omega19, omega21, omega15m, omega11, omega12m, omega4,
11 omega12p, omega20, omega0, omega17p, omega10, omega17m, omega5, omega6, omega23, omega1,
12 omega13p, omega14m, omega13m, omega7, omega22, omega2;
13 for (i = 6; i < n - 6; i++)
14 {
15 sigma0 = sigma[i * ssi + 0 * ssr];
16 sigma1 = sigma[i * ssi + 1 * ssr];
17 sigma2 = sigma[i * ssi + 2 * ssr];
18 sigma3 = sigma[i * ssi + 3 * ssr];
19 sigma4 = sigma[i * ssi + 4 * ssr];
20 sigma5 = sigma[i * ssi + 5 * ssr];
21 acc = 0.0;
22 omega0 = (+0.00216450216450216) / ((sigma0 + 1.0e-6) * (sigma0 + 1.0e-6));
23 acc = acc + omega0;
24 omega1 = (+0.0649350649350649) / ((sigma1 + 1.0e-6) * (sigma1 + 1.0e-6));
25 acc = acc + omega1;
26 omega2 = (+0.324675324675325) / ((sigma2 + 1.0e-6) * (sigma2 + 1.0e-6));
27 acc = acc + omega2;
28 omega3 = (+0.432900432900433) / ((sigma3 + 1.0e-6) * (sigma3 + 1.0e-6));
29 acc = acc + omega3;
30 omega4 = (+0.162337662337662) / ((sigma4 + 1.0e-6) * (sigma4 + 1.0e-6));
31 acc = acc + omega4;
32 omega5 = (+0.012987012987013) / ((sigma5 + 1.0e-6) * (sigma5 + 1.0e-6));
33 acc = acc + omega5;
34 omega0 = (omega0) / (acc);
35 omega1 = (omega1) / (acc);
36 omega2 = (omega2) / (acc);
37 omega3 = (omega3) / (acc);
38 omega4 = (omega4) / (acc);
39 omega5 = (omega5) / (acc);
40 acc = 0.0;
41 omega6 = (+0.00677201832660494) / ((sigma0 + 1.0e-6) * (sigma0 + 1.0e-6));
42 acc = acc + omega6;
43 omega7 = (+0.107265913289826) / ((sigma1 + 1.0e-6) * (sigma1 + 1.0e-6));
44 acc = acc + omega7;
45 omega8 = (+0.380324261915151) / ((sigma2 + 1.0e-6) * (sigma2 + 1.0e-6));
46 acc = acc + omega8;
47 omega9 = (+0.386977736276944) / ((sigma3 + 1.0e-6) * (sigma3 + 1.0e-6));
48 acc = acc + omega9;
49 omega10 = (+0.112027858604345) / ((sigma4 + 1.0e-6) * (sigma4 + 1.0e-6));
50 acc = acc + omega10;
51 omega11 = (+0.00663221158712994) / ((sigma5 + 1.0e-6) * (sigma5 + 1.0e-6));
52 acc = acc + omega11;
53 omega6 = (omega6) / (acc);
54 omega7 = (omega7) / (acc);
55 omega8 = (omega8) / (acc);
56 omega9 = (omega9) / (acc);
57 omega10 = (omega10) / (acc);
58 omega11 = (omega11) / (acc);
59 acc = 0.0;
60 omega12p =
61 ((+0.00402560307167669) / (+2.02092329375878)) / ((sigma0 + 1.0e-6) * (sigma0 + 1.0e-6));
62 acc = acc + omega12p;
63 omega13p =
64 ((+0.0948077070036824) / (+2.02092329375878)) / ((sigma1 + 1.0e-6) * (sigma1 + 1.0e-6));
65 acc = acc + omega13p;
66 omega14p =
67 ((+0.432572020266995) / (+2.02092329375878)) / ((sigma2 + 1.0e-6) * (sigma2 + 1.0e-6));
68 acc = acc + omega14p;
69 omega15p =
70 ((+0.0239353628232871) / (+2.02092329375878)) / ((sigma3 + 1.0e-6) * (sigma3 + 1.0e-6));
71 acc = acc + omega15p;
72 omega16p =
73 ((+1.45860816934021) / (+2.02092329375878)) / ((sigma4 + 1.0e-6) * (sigma4 + 1.0e-6));
74 acc = acc + omega16p;
75 omega17p =
76 ((+0.00697443125292721) / (+2.02092329375878)) / ((sigma5 + 1.0e-6) * (sigma5 + 1.0e-6));
77 acc = acc + omega17p;
78 omega12p = (omega12p) / (acc);
79 omega13p = (omega13p) / (acc);
80 omega14p = (omega14p) / (acc);
81 omega15p = (omega15p) / (acc);
82 omega16p = (omega16p) / (acc);
83 omega17p = (omega17p) / (acc);
84 acc = 0.0;
85 omega12m =
86 ((+0.00201280153583835) / (+1.02092329375878)) / ((sigma0 + 1.0e-6) * (sigma0 + 1.0e-6));
87 acc = acc + omega12m;
88 omega13m =
89 ((+0.0474038535018412) / (+1.02092329375878)) / ((sigma1 + 1.0e-6) * (sigma1 + 1.0e-6));
90 acc = acc + omega13m;
91 omega14m =
92 ((+0.216286010133497) / (+1.02092329375878)) / ((sigma2 + 1.0e-6) * (sigma2 + 1.0e-6));
93 acc = acc + omega14m;
94 omega15m =
95 ((+0.0119676814116435) / (+1.02092329375878)) / ((sigma3 + 1.0e-6) * (sigma3 + 1.0e-6));
96 acc = acc + omega15m;
97 omega16m =
98 ((+0.729304084670107) / (+1.02092329375878)) / ((sigma4 + 1.0e-6) * (sigma4 + 1.0e-6));
99 acc = acc + omega16m;
100 omega17m =
101 ((+0.0139488625058544) / (+1.02092329375878)) / ((sigma5 + 1.0e-6) * (sigma5 + 1.0e-6));
102 acc = acc + omega17m;
103 omega12m = (omega12m) / (acc);
104 omega13m = (omega13m) / (acc);
105 omega14m = (omega14m) / (acc);
106 omega15m = (omega15m) / (acc);
107 omega16m = (omega16m) / (acc);
108 omega17m = (omega17m) / (acc);
109 acc = 0.0;
110 omega18 = (+0.00999709318423371) / ((sigma0 + 1.0e-6) * (sigma0 + 1.0e-6));
111 acc = acc + omega18;
112 omega19 = (+0.141307791850495) / ((sigma1 + 1.0e-6) * (sigma1 + 1.0e-6));
113 acc = acc + omega19;
114 omega20 = (+0.417983745381249) / ((sigma2 + 1.0e-6) * (sigma2 + 1.0e-6));
115 acc = acc + omega20;
116 omega21 = (+0.348255604187079) / ((sigma3 + 1.0e-6) * (sigma3 + 1.0e-6));
117 acc = acc + omega21;
118 omega22 = (+0.0791777371169748) / ((sigma4 + 1.0e-6) * (sigma4 + 1.0e-6));
119 acc = acc + omega22;
120 omega23 = (+0.003278028279983) / ((sigma5 + 1.0e-6) * (sigma5 + 1.0e-6));
121 acc = acc + omega23;
122 omega18 = (omega18) / (acc);
123 omega19 = (omega19) / (acc);
124 omega20 = (omega20) / (acc);
125 omega21 = (omega21) / (acc);
126 omega22 = (omega22) / (acc);
127 omega23 = (omega23) / (acc);
128 omega[i * wsi + 0 * wsl + 0 * wsr + 0] = omega0;
129 omega[i * wsi + 0 * wsl + 1 * wsr + 0] = omega1;
130 omega[i * wsi + 0 * wsl + 2 * wsr + 0] = omega2;
131 omega[i * wsi + 0 * wsl + 3 * wsr + 0] = omega3;
132 omega[i * wsi + 0 * wsl + 4 * wsr + 0] = omega4;
133 omega[i * wsi + 0 * wsl + 5 * wsr + 0] = omega5;
134 omega[i * wsi + 1 * wsl + 0 * wsr + 0] = omega6;
135 omega[i * wsi + 1 * wsl + 1 * wsr + 0] = omega7;
136 omega[i * wsi + 1 * wsl + 2 * wsr + 0] = omega8;
137 omega[i * wsi + 1 * wsl + 3 * wsr + 0] = omega9;
138 omega[i * wsi + 1 * wsl + 4 * wsr + 0] = omega10;
139 omega[i * wsi + 1 * wsl + 5 * wsr + 0] = omega11;
140 omega[i * wsi + 2 * wsl + 0 * wsr + 0] = omega12p;
141 omega[i * wsi + 2 * wsl + 0 * wsr + 1] = omega12m;
142 omega[i * wsi + 2 * wsl + 1 * wsr + 0] = omega13p;
143 omega[i * wsi + 2 * wsl + 1 * wsr + 1] = omega13m;
144 omega[i * wsi + 2 * wsl + 2 * wsr + 0] = omega14p;
145 omega[i * wsi + 2 * wsl + 2 * wsr + 1] = omega14m;
146 omega[i * wsi + 2 * wsl + 3 * wsr + 0] = omega15p;
147 omega[i * wsi + 2 * wsl + 3 * wsr + 1] = omega15m;
148 omega[i * wsi + 2 * wsl + 4 * wsr + 0] = omega16p;
149 omega[i * wsi + 2 * wsl + 4 * wsr + 1] = omega16m;
150 omega[i * wsi + 2 * wsl + 5 * wsr + 0] = omega17p;
151 omega[i * wsi + 2 * wsl + 5 * wsr + 1] = omega17m;
152 omega[i * wsi + 3 * wsl + 0 * wsr + 0] = omega18;
153 omega[i * wsi + 3 * wsl + 1 * wsr + 0] = omega19;
154 omega[i * wsi + 3 * wsl + 2 * wsr + 0] = omega20;
155 omega[i * wsi + 3 * wsl + 3 * wsr + 0] = omega21;
156 omega[i * wsi + 3 * wsl + 4 * wsr + 0] = omega22;
157 omega[i * wsi + 3 * wsl + 5 * wsr + 0] = omega23;
158 }
159 }
160
161 PyObject *
162 py_weights_gauss_radau006004 (PyObject * self, PyObject * args)
163 {
164 double *sigma, *omega;
165 PyArrayObject *sigma_py, *omega_py;
166
167 long int n;
168 int ssi, ssr, wsi, wsl, wsr;
169
170 /* parse options */
171
172 if (!PyArg_ParseTuple (args, "OO", &sigma_py, &omega_py))
173 return NULL;
174
175 if (sigma_py->nd != 2 || sigma_py->descr->type_num != PyArray_DOUBLE)
176 {
177 PyErr_SetString (PyExc_ValueError, "sigma must be two-dimensional and of type float");
178 return NULL;
179 }
180
181 if (omega_py->descr->type_num != PyArray_DOUBLE)
182 {
183 PyErr_SetString (PyExc_ValueError, "omega must be of type float");
184 return NULL;
185 }
186
187 if (!(omega_py->nd >= 2 && omega_py->nd <= 4))
188 {
189 PyErr_SetString (PyExc_ValueError, "omega must be two, three, or four dimensional");
190 return NULL;
191 }
192
193 /* get data, n, strides */
194 sigma = (double *) PyArray_DATA (sigma_py);
195 omega = (double *) PyArray_DATA (omega_py);
196
197 n = PyArray_DIM (omega_py, 0);
198
199 ssi = sigma_py->strides[0] / sizeof (double);
200 ssr = sigma_py->strides[1] / sizeof (double);
201
202 wsi = omega_py->strides[0] / sizeof (double);
203 if (omega_py->nd == 3)
204 {
205 wsl = omega_py->strides[1] / sizeof (double);
206 wsr = omega_py->strides[2] / sizeof (double);
207 }
208 else
209 {
210 wsl = 0;
211 wsr = omega_py->strides[1] / sizeof (double);
212 }
213
214 weights_gauss_radau006004 (sigma, n, ssi, ssr, omega, wsi, wsl, wsr);
215
216 Py_INCREF (Py_None);
217 return Py_None;
218 }
219
220 void
221 reconstruct_gauss_radau006004 (const double *restrict f, int n, int fsi,
222 const double *restrict omega, int wsi, int wsl, int wsr,
223 double *restrict fr, int frsi, int frsl)
224 {
225 int i;
226 double omega9, omega18, omega15p, omega16m, omega14p, omega3, omega8, omega16p, omega19, omega21,
227 omega15m, omega11, omega12m, omega4, omega12p, omega20, omega0, omega17p, omega10, omega17m,
228 omega5, omega6, omega23, omega1, omega13p, omega14m, omega13m, omega7, omega22, omega2, fs0,
229 fs1, fs2, fs3, fr9, fr18, fr13, fr17, fr3, fr8, fr21, fr11, fr14, fr7, fr20, fr0, fr4, fr10,
230 fr5, fr6, fr23, fr1, fr19, fr2, fr12, fr15, fr22, fr16;
231 for (i = 6; i < n - 6; i++)
232 {
233 omega0 = omega[i * wsi + 0 * wsl + 0 * wsr + 0];
234 omega1 = omega[i * wsi + 0 * wsl + 1 * wsr + 0];
235 omega2 = omega[i * wsi + 0 * wsl + 2 * wsr + 0];
236 omega3 = omega[i * wsi + 0 * wsl + 3 * wsr + 0];
237 omega4 = omega[i * wsi + 0 * wsl + 4 * wsr + 0];
238 omega5 = omega[i * wsi + 0 * wsl + 5 * wsr + 0];
239 omega6 = omega[i * wsi + 1 * wsl + 0 * wsr + 0];
240 omega7 = omega[i * wsi + 1 * wsl + 1 * wsr + 0];
241 omega8 = omega[i * wsi + 1 * wsl + 2 * wsr + 0];
242 omega9 = omega[i * wsi + 1 * wsl + 3 * wsr + 0];
243 omega10 = omega[i * wsi + 1 * wsl + 4 * wsr + 0];
244 omega11 = omega[i * wsi + 1 * wsl + 5 * wsr + 0];
245 omega12p = omega[i * wsi + 2 * wsl + 0 * wsr + 0];
246 omega12m = omega[i * wsi + 2 * wsl + 0 * wsr + 1];
247 omega13p = omega[i * wsi + 2 * wsl + 1 * wsr + 0];
248 omega13m = omega[i * wsi + 2 * wsl + 1 * wsr + 1];
249 omega14p = omega[i * wsi + 2 * wsl + 2 * wsr + 0];
250 omega14m = omega[i * wsi + 2 * wsl + 2 * wsr + 1];
251 omega15p = omega[i * wsi + 2 * wsl + 3 * wsr + 0];
252 omega15m = omega[i * wsi + 2 * wsl + 3 * wsr + 1];
253 omega16p = omega[i * wsi + 2 * wsl + 4 * wsr + 0];
254 omega16m = omega[i * wsi + 2 * wsl + 4 * wsr + 1];
255 omega17p = omega[i * wsi + 2 * wsl + 5 * wsr + 0];
256 omega17m = omega[i * wsi + 2 * wsl + 5 * wsr + 1];
257 omega18 = omega[i * wsi + 3 * wsl + 0 * wsr + 0];
258 omega19 = omega[i * wsi + 3 * wsl + 1 * wsr + 0];
259 omega20 = omega[i * wsi + 3 * wsl + 2 * wsr + 0];
260 omega21 = omega[i * wsi + 3 * wsl + 3 * wsr + 0];
261 omega22 = omega[i * wsi + 3 * wsl + 4 * wsr + 0];
262 omega23 = omega[i * wsi + 3 * wsl + 5 * wsr + 0];
263 fr0 =
264 (+2.45) * (f[(i + 0) * fsi]) + (-3.55) * (f[(i + 1) * fsi]) + (+3.95) * (f[(i + 2) * fsi]) +
265 (-2.71666666666667) * (f[(i + 3) * fsi]) + (+1.03333333333333) * (f[(i + 4) * fsi]) +
266 (-0.166666666666667) * (f[(i + 5) * fsi]);
267 fr1 =
268 (+0.166666666666667) * (f[(i - 1) * fsi]) + (+1.45) * (f[(i + 0) * fsi]) +
269 (-1.05) * (f[(i + 1) * fsi]) + (+0.616666666666667) * (f[(i + 2) * fsi]) +
270 (-0.216666666666667) * (f[(i + 3) * fsi]) + (+0.0333333333333333) * (f[(i + 4) * fsi]);
271 fr2 =
272 (-0.0333333333333333) * (f[(i - 2) * fsi]) + (+0.366666666666667) * (f[(i - 1) * fsi]) +
273 (+0.95) * (f[(i + 0) * fsi]) + (-0.383333333333333) * (f[(i + 1) * fsi]) +
274 (+0.116666666666667) * (f[(i + 2) * fsi]) + (-0.0166666666666667) * (f[(i + 3) * fsi]);
275 fr3 =
276 (+0.0166666666666667) * (f[(i - 3) * fsi]) + (-0.133333333333333) * (f[(i - 2) * fsi]) +
277 (+0.616666666666667) * (f[(i - 1) * fsi]) + (+0.616666666666667) * (f[(i + 0) * fsi]) +
278 (-0.133333333333333) * (f[(i + 1) * fsi]) + (+0.0166666666666667) * (f[(i + 2) * fsi]);
279 fr4 =
280 (-0.0166666666666667) * (f[(i - 4) * fsi]) + (+0.116666666666667) * (f[(i - 3) * fsi]) +
281 (-0.383333333333333) * (f[(i - 2) * fsi]) + (+0.95) * (f[(i - 1) * fsi]) +
282 (+0.366666666666667) * (f[(i + 0) * fsi]) + (-0.0333333333333333) * (f[(i + 1) * fsi]);
283 fr5 =
284 (+0.0333333333333333) * (f[(i - 5) * fsi]) + (-0.216666666666667) * (f[(i - 4) * fsi]) +
285 (+0.616666666666667) * (f[(i - 3) * fsi]) + (-1.05) * (f[(i - 2) * fsi]) +
286 (+1.45) * (f[(i - 1) * fsi]) + (+0.166666666666667) * (f[(i + 0) * fsi]);
287 fr6 =
288 (+1.62117649102017) * (f[(i + 0) * fsi]) + (-1.29008774199659) * (f[(i + 1) * fsi]) +
289 (+1.19272916844303) * (f[(i + 2) * fsi]) + (-0.755251069323849) * (f[(i + 3) * fsi]) +
290 (+0.274442607964963) * (f[(i + 4) * fsi]) + (-0.0430094561077231) * (f[(i + 5) * fsi]);
291 fr7 =
292 (+0.0430094561077231) * (f[(i - 1) * fsi]) + (+1.36311975437383) * (f[(i + 0) * fsi]) +
293 (-0.644945900380743) * (f[(i + 1) * fsi]) + (+0.332540046288564) * (f[(i + 2) * fsi]) +
294 (-0.110109227708003) * (f[(i + 3) * fsi]) + (+0.0163858713186246) * (f[(i + 4) * fsi]);
295 fr8 =
296 (-0.0163858713186246) * (f[(i - 2) * fsi]) + (+0.141324684019471) * (f[(i - 1) * fsi]) +
297 (+1.11733168459446) * (f[(i + 0) * fsi]) + (-0.31722847400825) * (f[(i + 1) * fsi]) +
298 (+0.0867519765091948) * (f[(i + 2) * fsi]) + (-0.0117939997962549) * (f[(i + 3) * fsi]);
299 fr9 =
300 (+0.0117939997962549) * (f[(i - 3) * fsi]) + (-0.0871498700961542) * (f[(i - 2) * fsi]) +
301 (+0.318234680963295) * (f[(i - 1) * fsi]) + (+0.881451688669366) * (f[(i + 0) * fsi]) +
302 (-0.140318477064426) * (f[(i + 1) * fsi]) + (+0.0159879777316652) * (f[(i + 2) * fsi]);
303 fr10 =
304 (-0.0159879777316652) * (f[(i - 4) * fsi]) + (+0.107721866186246) * (f[(i - 3) * fsi]) +
305 (-0.326969536071132) * (f[(i - 2) * fsi]) + (+0.637994235596599) * (f[(i - 1) * fsi]) +
306 (+0.641632022694388) * (f[(i + 0) * fsi]) + (-0.0443906106744351) * (f[(i + 1) * fsi]);
307 fr11 =
308 (+0.0443906106744351) * (f[(i - 5) * fsi]) + (-0.282331641778276) * (f[(i - 4) * fsi]) +
309 (+0.773581026302773) * (f[(i - 3) * fsi]) + (-1.21478174955983) * (f[(i - 2) * fsi]) +
310 (+1.30385339571313) * (f[(i - 1) * fsi]) + (+0.375288358647777) * (f[(i + 0) * fsi]);
311 fr12 =
312 (+0.67094398787024) * (f[(i + 0) * fsi]) + (+0.855895783827608) * (f[(i + 1) * fsi]) +
313 (-0.988824869991825) * (f[(i + 2) * fsi]) + (+0.67669237594118) * (f[(i + 3) * fsi]) +
314 (-0.25574635212777) * (f[(i + 4) * fsi]) + (+0.0410390744805702) * (f[(i + 5) * fsi]);
315 fr13 =
316 (-0.0410390744805702) * (f[(i - 1) * fsi]) + (+0.91717843475366) * (f[(i + 0) * fsi]) +
317 (+0.240309666619054) * (f[(i + 1) * fsi]) + (-0.168043380380424) * (f[(i + 2) * fsi]) +
318 (+0.0611062587326286) * (f[(i + 3) * fsi]) + (-0.00951190524434911) * (f[(i + 4) * fsi]);
319 fr14 =
320 (+0.0095119052443491) * (f[(i - 2) * fsi]) + (-0.0981105059466648) * (f[(i - 1) * fsi]) +
321 (+1.0598570134189) * (f[(i + 0) * fsi]) + (+0.0500715617320723) * (f[(i + 1) * fsi]) +
322 (-0.0253648017151878) * (f[(i + 2) * fsi]) + (+0.00403482726653401) * (f[(i + 3) * fsi]);
323 fr15 =
324 (-0.00403482726653401) * (f[(i - 3) * fsi]) + (+0.0337208688435532) * (f[(i - 2) * fsi]) +
325 (-0.158632914944675) * (f[(i - 1) * fsi]) + (+1.14055355874958) * (f[(i + 0) * fsi]) +
326 (-0.0104508472659379) * (f[(i + 1) * fsi]) + (-0.00115583811598369) * (f[(i + 2) * fsi]);
327 fr16 =
328 (+0.00115583811598369) * (f[(i - 4) * fsi]) + (-0.0109698559624361) * (f[(i - 3) * fsi]) +
329 (+0.0510584405833085) * (f[(i - 2) * fsi]) + (-0.181749677264349) * (f[(i - 1) * fsi]) +
330 (+1.15789113048933) * (f[(i + 0) * fsi]) + (-0.01738587596184) * (f[(i + 1) * fsi]);
331 fr17 =
332 (+0.01738587596184) * (f[(i - 5) * fsi]) + (-0.103159417655056) * (f[(i - 4) * fsi]) +
333 (+0.249818283465164) * (f[(i - 3) * fsi]) + (-0.296659078653492) * (f[(i - 2) * fsi]) +
334 (+0.0790384621632518) * (f[(i - 1) * fsi]) + (+1.05357587471829) * (f[(i + 0) * fsi]);
335 fr18 =
336 (+0.241784045018277) * (f[(i + 0) * fsi]) + (+1.423644765531) * (f[(i + 1) * fsi]) +
337 (-1.15775785153424) * (f[(i + 2) * fsi]) + (+0.705173595903445) * (f[(i + 3) * fsi]) +
338 (-0.251984964906828) * (f[(i + 4) * fsi]) + (+0.0391404099883515) * (f[(i + 5) * fsi]);
339 fr19 =
340 (-0.0391404099883517) * (f[(i - 1) * fsi]) + (+0.476626504948387) * (f[(i + 0) * fsi]) +
341 (+0.836538615705721) * (f[(i + 1) * fsi]) + (-0.374949651767207) * (f[(i + 2) * fsi]) +
342 (+0.118067446078169) * (f[(i + 3) * fsi]) + (-0.017142504976717) * (f[(i + 4) * fsi]);
343 fr20 =
344 (+0.017142504976717) * (f[(i - 2) * fsi]) + (-0.141995439848654) * (f[(i - 1) * fsi]) +
345 (+0.733764079599142) * (f[(i + 0) * fsi]) + (+0.49368851617138) * (f[(i + 1) * fsi]) +
346 (-0.117812077116452) * (f[(i + 2) * fsi]) + (+0.0152124162178667) * (f[(i + 3) * fsi]);
347 fr21 =
348 (-0.0152124162178667) * (f[(i - 3) * fsi]) + (+0.108417002283917) * (f[(i - 2) * fsi]) +
349 (-0.370181683116655) * (f[(i - 1) * fsi]) + (+1.03801240395648) * (f[(i + 0) * fsi]) +
350 (+0.265502272903379) * (f[(i + 1) * fsi]) + (-0.0265375798092516) * (f[(i + 2) * fsi]);
351 fr22 =
352 (+0.0265375798092516) * (f[(i - 4) * fsi]) + (-0.174437895073377) * (f[(i - 3) * fsi]) +
353 (+0.506480699422692) * (f[(i - 2) * fsi]) + (-0.900933279301688) * (f[(i - 1) * fsi]) +
354 (+1.43607610109525) * (f[(i + 0) * fsi]) + (+0.10627679404787) * (f[(i + 1) * fsi]);
355 fr23 =
356 (-0.10627679404787) * (f[(i - 5) * fsi]) + (+0.664198344096469) * (f[(i - 4) * fsi]) +
357 (-1.76858980579142) * (f[(i - 3) * fsi]) + (+2.63201658038008) * (f[(i - 2) * fsi]) +
358 (-2.49508519001973) * (f[(i - 1) * fsi]) + (+2.07373686538247) * (f[(i + 0) * fsi]);
359 fs0 =
360 (omega0) * (fr0) + (omega1) * (fr1) + (omega2) * (fr2) + (omega3) * (fr3) +
361 (omega4) * (fr4) + (omega5) * (fr5);
362 fs1 =
363 (omega6) * (fr6) + (omega7) * (fr7) + (omega8) * (fr8) + (omega9) * (fr9) +
364 (omega10) * (fr10) + (omega11) * (fr11);
365 fs2 =
366 ((+2.02092329375878) *
367 ((omega12p) * (fr12) + (omega13p) * (fr13) + (omega14p) * (fr14) + (omega15p) * (fr15) +
368 (omega16p) * (fr16) + (omega17p) * (fr17))) -
369 ((+1.02092329375878) *
370 ((omega12m) * (fr12) + (omega13m) * (fr13) + (omega14m) * (fr14) + (omega15m) * (fr15) +
371 (omega16m) * (fr16) + (omega17m) * (fr17)));
372 fs3 =
373 (omega18) * (fr18) + (omega19) * (fr19) + (omega20) * (fr20) + (omega21) * (fr21) +
374 (omega22) * (fr22) + (omega23) * (fr23);
375 fr[i * frsi + 0 * frsl] = fs0;
376 fr[i * frsi + 1 * frsl] = fs1;
377 fr[i * frsi + 2 * frsl] = fs2;
378 fr[i * frsi + 3 * frsl] = fs3;
379 }
380 }
381
382 PyObject *
383 py_reconstruct_gauss_radau006004 (PyObject * self, PyObject * args)
384 {
385 double *f, *omega, *fr;
386 PyArrayObject *f_py, *omega_py, *fr_py;
387
388 long int n;
389 int fsi, frsi, frsl, wsi, wsl, wsr;
390
391 /* parse options */
392
393 if (!PyArg_ParseTuple (args, "OOO", &f_py, &omega_py, &fr_py))
394 return NULL;
395
396 if (f_py->nd != 1 || f_py->descr->type_num != PyArray_DOUBLE)
397 {
398 PyErr_SetString (PyExc_ValueError, "f must be one-dimensional and of type float");
399 return NULL;
400 }
401
402 if (fr_py->descr->type_num != PyArray_DOUBLE)
403 {
404 PyErr_SetString (PyExc_ValueError, "fr must be of type float");
405 return NULL;
406 }
407
408 if (!(fr_py->nd == 1 || fr_py->nd == 2))
409 {
410 PyErr_SetString (PyExc_ValueError, "fr must be one or two dimensional");
411 return NULL;
412 }
413
414 if (omega_py->descr->type_num != PyArray_DOUBLE)
415 {
416 PyErr_SetString (PyExc_ValueError, "omega must be of type float");
417 return NULL;
418 }
419
420 if (!(omega_py->nd >= 2 && omega_py->nd <= 4))
421 {
422 PyErr_SetString (PyExc_ValueError, "omega must be two, three, or four dimensional");
423 return NULL;
424 }
425
426 /* get data, n, strides */
427 f = (double *) PyArray_DATA (f_py);
428 fr = (double *) PyArray_DATA (fr_py);
429 omega = (double *) PyArray_DATA (omega_py);
430
431 n = PyArray_DIM (omega_py, 0);
432
433 fsi = f_py->strides[0] / sizeof (double);
434 frsi = fr_py->strides[0] / sizeof (double);
435
436 if (n == 1)
437 {
438 frsl = 0;
439 }
440 else
441 {
442 frsl = fr_py->strides[1] / sizeof (double);
443 }
444
445 wsi = omega_py->strides[0] / sizeof (double);
446 if (omega_py->nd == 3)
447 {
448 wsl = omega_py->strides[1] / sizeof (double);
449 wsr = omega_py->strides[2] / sizeof (double);
450 }
451 else
452 {
453 wsl = 0;
454 wsr = omega_py->strides[1] / sizeof (double);
455 }
456
457 reconstruct_gauss_radau006004 (f, n, fsi, omega, wsi, wsl, wsr, fr, frsi, frsl);
458
459 Py_INCREF (Py_None);
460 return Py_None;
461 }
Something went wrong with that request. Please try again.