1 | // Python - C extension for quantity module. |
---|
2 | // |
---|
3 | // To compile (Python2.3): |
---|
4 | // gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O |
---|
5 | // gcc -shared util_ext.o -o util_ext.so |
---|
6 | // |
---|
7 | // See the module quantity.py |
---|
8 | // |
---|
9 | // |
---|
10 | // Ole Nielsen, GA 2004 |
---|
11 | |
---|
12 | #include "Python.h" |
---|
13 | #include "Numeric/arrayobject.h" |
---|
14 | #include "math.h" |
---|
15 | |
---|
16 | //Shared code snippets |
---|
17 | #include "util_ext.h" |
---|
18 | //#include "quantity_ext.h" //Obsolete |
---|
19 | |
---|
20 | |
---|
21 | |
---|
22 | int _compute_gradients(int N, |
---|
23 | double* centroids, |
---|
24 | double* centroid_values, |
---|
25 | long* number_of_boundaries, |
---|
26 | long* surrogate_neighbours, |
---|
27 | double* a, |
---|
28 | double* b){ |
---|
29 | |
---|
30 | int i, k, k0, k1, k2, index3; |
---|
31 | double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det; |
---|
32 | |
---|
33 | |
---|
34 | for (k=0; k<N; k++) { |
---|
35 | index3 = 3*k; |
---|
36 | |
---|
37 | if (number_of_boundaries[k] < 2) { |
---|
38 | // Two or three true neighbours |
---|
39 | |
---|
40 | // Get indices of neighbours (or self when used as surrogate) |
---|
41 | // k0, k1, k2 = surrogate_neighbours[k,:] |
---|
42 | |
---|
43 | k0 = surrogate_neighbours[index3 + 0]; |
---|
44 | k1 = surrogate_neighbours[index3 + 1]; |
---|
45 | k2 = surrogate_neighbours[index3 + 2]; |
---|
46 | |
---|
47 | |
---|
48 | if (k0 == k1 || k1 == k2) return -1; |
---|
49 | |
---|
50 | // Get data |
---|
51 | q0 = centroid_values[k0]; |
---|
52 | q1 = centroid_values[k1]; |
---|
53 | q2 = centroid_values[k2]; |
---|
54 | |
---|
55 | x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; |
---|
56 | x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; |
---|
57 | x2 = centroids[k2*2]; y2 = centroids[k2*2+1]; |
---|
58 | |
---|
59 | // Gradient |
---|
60 | _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]); |
---|
61 | |
---|
62 | } else if (number_of_boundaries[k] == 2) { |
---|
63 | // One true neighbour |
---|
64 | |
---|
65 | // Get index of the one neighbour |
---|
66 | i=0; k0 = k; |
---|
67 | while (i<3 && k0==k) { |
---|
68 | k0 = surrogate_neighbours[index3 + i]; |
---|
69 | i++; |
---|
70 | } |
---|
71 | if (k0 == k) return -1; |
---|
72 | |
---|
73 | k1 = k; //self |
---|
74 | |
---|
75 | // Get data |
---|
76 | q0 = centroid_values[k0]; |
---|
77 | q1 = centroid_values[k1]; |
---|
78 | |
---|
79 | x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; |
---|
80 | x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; |
---|
81 | |
---|
82 | // Two point gradient |
---|
83 | _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]); |
---|
84 | |
---|
85 | } |
---|
86 | // else: |
---|
87 | // #No true neighbours - |
---|
88 | // #Fall back to first order scheme |
---|
89 | } |
---|
90 | return 0; |
---|
91 | } |
---|
92 | |
---|
93 | |
---|
94 | |
---|
95 | int _extrapolate(int N, |
---|
96 | double* centroids, |
---|
97 | double* centroid_values, |
---|
98 | double* vertex_coordinates, |
---|
99 | double* vertex_values, |
---|
100 | double* a, |
---|
101 | double* b) { |
---|
102 | |
---|
103 | int k, k2, k3, k6; |
---|
104 | double x, y, x0, y0, x1, y1, x2, y2; |
---|
105 | |
---|
106 | for (k=0; k<N; k++){ |
---|
107 | k6 = 6*k; |
---|
108 | k3 = 3*k; |
---|
109 | k2 = 2*k; |
---|
110 | |
---|
111 | // Centroid coordinates |
---|
112 | x = centroids[k2]; y = centroids[k2+1]; |
---|
113 | |
---|
114 | // vertex coordinates |
---|
115 | // x0, y0, x1, y1, x2, y2 = X[k,:] |
---|
116 | x0 = vertex_coordinates[k6 + 0]; |
---|
117 | y0 = vertex_coordinates[k6 + 1]; |
---|
118 | x1 = vertex_coordinates[k6 + 2]; |
---|
119 | y1 = vertex_coordinates[k6 + 3]; |
---|
120 | x2 = vertex_coordinates[k6 + 4]; |
---|
121 | y2 = vertex_coordinates[k6 + 5]; |
---|
122 | |
---|
123 | // Extrapolate |
---|
124 | vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y); |
---|
125 | vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y); |
---|
126 | vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y); |
---|
127 | |
---|
128 | } |
---|
129 | return 0; |
---|
130 | } |
---|
131 | |
---|
132 | |
---|
133 | |
---|
134 | |
---|
135 | int _interpolate(int N, |
---|
136 | double* vertex_values, |
---|
137 | double* edge_values) { |
---|
138 | |
---|
139 | int k, k3; |
---|
140 | double q0, q1, q2; |
---|
141 | |
---|
142 | |
---|
143 | for (k=0; k<N; k++) { |
---|
144 | k3 = 3*k; |
---|
145 | |
---|
146 | q0 = vertex_values[k3 + 0]; |
---|
147 | q1 = vertex_values[k3 + 1]; |
---|
148 | q2 = vertex_values[k3 + 2]; |
---|
149 | |
---|
150 | edge_values[k3 + 0] = 0.5*(q1+q2); |
---|
151 | edge_values[k3 + 1] = 0.5*(q0+q2); |
---|
152 | edge_values[k3 + 2] = 0.5*(q0+q1); |
---|
153 | } |
---|
154 | return 0; |
---|
155 | } |
---|
156 | |
---|
157 | int _backup_centroid_values(int N, |
---|
158 | double* centroid_values, |
---|
159 | double* centroid_backup_values) { |
---|
160 | // Backup centroid values |
---|
161 | |
---|
162 | |
---|
163 | int k; |
---|
164 | |
---|
165 | for (k=0; k<N; k++) { |
---|
166 | centroid_backup_values[k] = centroid_values[k]; |
---|
167 | } |
---|
168 | |
---|
169 | |
---|
170 | return 0; |
---|
171 | } |
---|
172 | |
---|
173 | |
---|
174 | int _saxpy_centroid_values(int N, |
---|
175 | double a, |
---|
176 | double b, |
---|
177 | double* centroid_values, |
---|
178 | double* centroid_backup_values) { |
---|
179 | // Saxby centroid values |
---|
180 | |
---|
181 | |
---|
182 | int k; |
---|
183 | |
---|
184 | |
---|
185 | for (k=0; k<N; k++) { |
---|
186 | centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k]; |
---|
187 | } |
---|
188 | |
---|
189 | |
---|
190 | return 0; |
---|
191 | } |
---|
192 | |
---|
193 | |
---|
194 | int _update(int N, |
---|
195 | double timestep, |
---|
196 | double* centroid_values, |
---|
197 | double* explicit_update, |
---|
198 | double* semi_implicit_update) { |
---|
199 | // Update centroid values based on values stored in |
---|
200 | // explicit_update and semi_implicit_update as well as given timestep |
---|
201 | |
---|
202 | |
---|
203 | int k; |
---|
204 | double denominator, x; |
---|
205 | |
---|
206 | |
---|
207 | // Divide semi_implicit update by conserved quantity |
---|
208 | for (k=0; k<N; k++) { |
---|
209 | x = centroid_values[k]; |
---|
210 | if (x == 0.0) { |
---|
211 | semi_implicit_update[k] = 0.0; |
---|
212 | } else { |
---|
213 | semi_implicit_update[k] /= x; |
---|
214 | } |
---|
215 | } |
---|
216 | |
---|
217 | |
---|
218 | // Semi implicit updates |
---|
219 | for (k=0; k<N; k++) { |
---|
220 | denominator = 1.0 - timestep*semi_implicit_update[k]; |
---|
221 | if (denominator == 0.0) { |
---|
222 | return -1; |
---|
223 | } else { |
---|
224 | //Update conserved_quantities from semi implicit updates |
---|
225 | centroid_values[k] /= denominator; |
---|
226 | } |
---|
227 | } |
---|
228 | |
---|
229 | |
---|
230 | // Explicit updates |
---|
231 | for (k=0; k<N; k++) { |
---|
232 | centroid_values[k] += timestep*explicit_update[k]; |
---|
233 | } |
---|
234 | |
---|
235 | |
---|
236 | // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py |
---|
237 | for (k=0;k<N;k++){ |
---|
238 | semi_implicit_update[k]=0.0; |
---|
239 | } |
---|
240 | |
---|
241 | return 0; |
---|
242 | } |
---|
243 | |
---|
244 | |
---|
245 | int _average_vertex_values(int N, |
---|
246 | long* vertex_value_indices, |
---|
247 | long* number_of_triangles_per_node, |
---|
248 | double* vertex_values, |
---|
249 | double* A) { |
---|
250 | // Average vertex values to obtain one value per node |
---|
251 | |
---|
252 | int i, index; |
---|
253 | int k = 0; //Track triangles touching each node |
---|
254 | int current_node = 0; |
---|
255 | double total = 0.0; |
---|
256 | |
---|
257 | for (i=0; i<N; i++) { |
---|
258 | |
---|
259 | // if (current_node == N) { |
---|
260 | // printf("Current node exceeding number of nodes (%d)", N); |
---|
261 | // return 1; |
---|
262 | // } |
---|
263 | |
---|
264 | index = vertex_value_indices[i]; |
---|
265 | k += 1; |
---|
266 | |
---|
267 | // volume_id = index / 3 |
---|
268 | // vertex_id = index % 3 |
---|
269 | // total += self.vertex_values[volume_id, vertex_id] |
---|
270 | total += vertex_values[index]; |
---|
271 | |
---|
272 | // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total); |
---|
273 | if (number_of_triangles_per_node[current_node] == k) { |
---|
274 | A[current_node] = total/k; |
---|
275 | |
---|
276 | // Move on to next node |
---|
277 | total = 0.0; |
---|
278 | k = 0; |
---|
279 | current_node += 1; |
---|
280 | } |
---|
281 | } |
---|
282 | |
---|
283 | return 0; |
---|
284 | } |
---|
285 | |
---|
286 | |
---|
287 | ///////////////////////////////////////////////// |
---|
288 | // Gateways to Python |
---|
289 | PyObject *update(PyObject *self, PyObject *args) { |
---|
290 | |
---|
291 | PyObject *quantity; |
---|
292 | PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update; |
---|
293 | |
---|
294 | double timestep; |
---|
295 | int N, err; |
---|
296 | |
---|
297 | |
---|
298 | // Convert Python arguments to C |
---|
299 | if (!PyArg_ParseTuple(args, "Od", &quantity, ×tep)) { |
---|
300 | PyErr_SetString(PyExc_RuntimeError, |
---|
301 | "quantity_ext.c: update could not parse input"); |
---|
302 | return NULL; |
---|
303 | } |
---|
304 | |
---|
305 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
306 | explicit_update = get_consecutive_array(quantity, "explicit_update"); |
---|
307 | semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update"); |
---|
308 | |
---|
309 | N = centroid_values -> dimensions[0]; |
---|
310 | |
---|
311 | err = _update(N, timestep, |
---|
312 | (double*) centroid_values -> data, |
---|
313 | (double*) explicit_update -> data, |
---|
314 | (double*) semi_implicit_update -> data); |
---|
315 | |
---|
316 | |
---|
317 | if (err != 0) { |
---|
318 | PyErr_SetString(PyExc_RuntimeError, |
---|
319 | "Zero division in semi implicit update - call Stephen :)"); |
---|
320 | return NULL; |
---|
321 | } |
---|
322 | |
---|
323 | // Release and return |
---|
324 | Py_DECREF(centroid_values); |
---|
325 | Py_DECREF(explicit_update); |
---|
326 | Py_DECREF(semi_implicit_update); |
---|
327 | |
---|
328 | return Py_BuildValue(""); |
---|
329 | } |
---|
330 | |
---|
331 | |
---|
332 | PyObject *backup_centroid_values(PyObject *self, PyObject *args) { |
---|
333 | |
---|
334 | PyObject *quantity; |
---|
335 | PyArrayObject *centroid_values, *centroid_backup_values; |
---|
336 | |
---|
337 | int N, err; |
---|
338 | |
---|
339 | |
---|
340 | // Convert Python arguments to C |
---|
341 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
342 | PyErr_SetString(PyExc_RuntimeError, |
---|
343 | "quantity_ext.c: backup_centroid_values could not parse input"); |
---|
344 | return NULL; |
---|
345 | } |
---|
346 | |
---|
347 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
348 | centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values"); |
---|
349 | |
---|
350 | N = centroid_values -> dimensions[0]; |
---|
351 | |
---|
352 | err = _backup_centroid_values(N, |
---|
353 | (double*) centroid_values -> data, |
---|
354 | (double*) centroid_backup_values -> data); |
---|
355 | |
---|
356 | |
---|
357 | // Release and return |
---|
358 | Py_DECREF(centroid_values); |
---|
359 | Py_DECREF(centroid_backup_values); |
---|
360 | |
---|
361 | return Py_BuildValue(""); |
---|
362 | } |
---|
363 | |
---|
364 | PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) { |
---|
365 | |
---|
366 | PyObject *quantity; |
---|
367 | PyArrayObject *centroid_values, *centroid_backup_values; |
---|
368 | |
---|
369 | double a,b; |
---|
370 | int N, err; |
---|
371 | |
---|
372 | |
---|
373 | // Convert Python arguments to C |
---|
374 | if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) { |
---|
375 | PyErr_SetString(PyExc_RuntimeError, |
---|
376 | "quantity_ext.c: saxpy_centroid_values could not parse input"); |
---|
377 | return NULL; |
---|
378 | } |
---|
379 | |
---|
380 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
381 | centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values"); |
---|
382 | |
---|
383 | N = centroid_values -> dimensions[0]; |
---|
384 | |
---|
385 | err = _saxpy_centroid_values(N,a,b, |
---|
386 | (double*) centroid_values -> data, |
---|
387 | (double*) centroid_backup_values -> data); |
---|
388 | |
---|
389 | |
---|
390 | // Release and return |
---|
391 | Py_DECREF(centroid_values); |
---|
392 | Py_DECREF(centroid_backup_values); |
---|
393 | |
---|
394 | return Py_BuildValue(""); |
---|
395 | } |
---|
396 | |
---|
397 | |
---|
398 | PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) { |
---|
399 | |
---|
400 | PyObject *quantity; |
---|
401 | PyArrayObject *vertex_values, *edge_values; |
---|
402 | |
---|
403 | int N, err; |
---|
404 | |
---|
405 | // Convert Python arguments to C |
---|
406 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
407 | PyErr_SetString(PyExc_RuntimeError, |
---|
408 | "quantity_ext.c: interpolate_from_vertices_to_edges could not parse input"); |
---|
409 | return NULL; |
---|
410 | } |
---|
411 | |
---|
412 | vertex_values = get_consecutive_array(quantity, "vertex_values"); |
---|
413 | edge_values = get_consecutive_array(quantity, "edge_values"); |
---|
414 | |
---|
415 | N = vertex_values -> dimensions[0]; |
---|
416 | |
---|
417 | err = _interpolate(N, |
---|
418 | (double*) vertex_values -> data, |
---|
419 | (double*) edge_values -> data); |
---|
420 | |
---|
421 | if (err != 0) { |
---|
422 | PyErr_SetString(PyExc_RuntimeError, |
---|
423 | "Interpolate could not be computed"); |
---|
424 | return NULL; |
---|
425 | } |
---|
426 | |
---|
427 | // Release and return |
---|
428 | Py_DECREF(vertex_values); |
---|
429 | Py_DECREF(edge_values); |
---|
430 | |
---|
431 | return Py_BuildValue(""); |
---|
432 | } |
---|
433 | |
---|
434 | |
---|
435 | PyObject *average_vertex_values(PyObject *self, PyObject *args) { |
---|
436 | |
---|
437 | PyArrayObject |
---|
438 | *vertex_value_indices, |
---|
439 | *number_of_triangles_per_node, |
---|
440 | *vertex_values, |
---|
441 | *A; |
---|
442 | |
---|
443 | |
---|
444 | int N, err; |
---|
445 | |
---|
446 | // Convert Python arguments to C |
---|
447 | if (!PyArg_ParseTuple(args, "OOOO", |
---|
448 | &vertex_value_indices, |
---|
449 | &number_of_triangles_per_node, |
---|
450 | &vertex_values, |
---|
451 | &A)) { |
---|
452 | PyErr_SetString(PyExc_RuntimeError, |
---|
453 | "quantity_ext.c: average_vertex_values could not parse input"); |
---|
454 | return NULL; |
---|
455 | } |
---|
456 | |
---|
457 | N = vertex_value_indices -> dimensions[0]; |
---|
458 | // printf("Got parameters, N=%d\n", N); |
---|
459 | err = _average_vertex_values(N, |
---|
460 | (long*) vertex_value_indices -> data, |
---|
461 | (long*) number_of_triangles_per_node -> data, |
---|
462 | (double*) vertex_values -> data, |
---|
463 | (double*) A -> data); |
---|
464 | |
---|
465 | //printf("Error %d", err); |
---|
466 | if (err != 0) { |
---|
467 | PyErr_SetString(PyExc_RuntimeError, |
---|
468 | "average_vertex_values could not be computed"); |
---|
469 | return NULL; |
---|
470 | } |
---|
471 | |
---|
472 | return Py_BuildValue(""); |
---|
473 | } |
---|
474 | |
---|
475 | |
---|
476 | |
---|
477 | PyObject *compute_gradients(PyObject *self, PyObject *args) { |
---|
478 | |
---|
479 | PyObject *quantity, *domain, *R; |
---|
480 | PyArrayObject |
---|
481 | *centroids, //Coordinates at centroids |
---|
482 | *centroid_values, //Values at centroids |
---|
483 | *number_of_boundaries, //Number of boundaries for each triangle |
---|
484 | *surrogate_neighbours, //True neighbours or - if one missing - self |
---|
485 | *a, *b; //Return values |
---|
486 | |
---|
487 | int dimensions[1], N, err; |
---|
488 | |
---|
489 | // Convert Python arguments to C |
---|
490 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
491 | PyErr_SetString(PyExc_RuntimeError, |
---|
492 | "quantity_ext.c: compute_gradients could not parse input"); |
---|
493 | return NULL; |
---|
494 | } |
---|
495 | |
---|
496 | domain = PyObject_GetAttrString(quantity, "domain"); |
---|
497 | if (!domain) { |
---|
498 | PyErr_SetString(PyExc_RuntimeError, |
---|
499 | "compute_gradients could not obtain domain object from quantity"); |
---|
500 | return NULL; |
---|
501 | } |
---|
502 | |
---|
503 | // Get pertinent variables |
---|
504 | |
---|
505 | centroids = get_consecutive_array(domain, "centroid_coordinates"); |
---|
506 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
507 | surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); |
---|
508 | number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); |
---|
509 | |
---|
510 | N = centroid_values -> dimensions[0]; |
---|
511 | |
---|
512 | // Release |
---|
513 | Py_DECREF(domain); |
---|
514 | |
---|
515 | // Allocate space for return vectors a and b (don't DECREF) |
---|
516 | dimensions[0] = N; |
---|
517 | a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
518 | b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
519 | |
---|
520 | |
---|
521 | |
---|
522 | err = _compute_gradients(N, |
---|
523 | (double*) centroids -> data, |
---|
524 | (double*) centroid_values -> data, |
---|
525 | (long*) number_of_boundaries -> data, |
---|
526 | (long*) surrogate_neighbours -> data, |
---|
527 | (double*) a -> data, |
---|
528 | (double*) b -> data); |
---|
529 | |
---|
530 | if (err != 0) { |
---|
531 | PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); |
---|
532 | return NULL; |
---|
533 | } |
---|
534 | |
---|
535 | // Release |
---|
536 | Py_DECREF(centroids); |
---|
537 | Py_DECREF(centroid_values); |
---|
538 | Py_DECREF(number_of_boundaries); |
---|
539 | Py_DECREF(surrogate_neighbours); |
---|
540 | |
---|
541 | // Build result, release and return |
---|
542 | R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); |
---|
543 | Py_DECREF(a); |
---|
544 | Py_DECREF(b); |
---|
545 | return R; |
---|
546 | } |
---|
547 | |
---|
548 | |
---|
549 | |
---|
550 | PyObject *extrapolate_second_order(PyObject *self, PyObject *args) { |
---|
551 | |
---|
552 | PyObject *quantity, *domain; |
---|
553 | PyArrayObject |
---|
554 | *centroids, //Coordinates at centroids |
---|
555 | *centroid_values, //Values at centroids |
---|
556 | *vertex_coordinates, //Coordinates at vertices |
---|
557 | *vertex_values, //Values at vertices |
---|
558 | *number_of_boundaries, //Number of boundaries for each triangle |
---|
559 | *surrogate_neighbours, //True neighbours or - if one missing - self |
---|
560 | *a, *b; //Gradients |
---|
561 | |
---|
562 | //int N, err; |
---|
563 | int dimensions[1], N, err; |
---|
564 | //double *a, *b; //Gradients |
---|
565 | |
---|
566 | // Convert Python arguments to C |
---|
567 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
568 | PyErr_SetString(PyExc_RuntimeError, |
---|
569 | "extrapolate_second_order could not parse input"); |
---|
570 | return NULL; |
---|
571 | } |
---|
572 | |
---|
573 | domain = PyObject_GetAttrString(quantity, "domain"); |
---|
574 | if (!domain) { |
---|
575 | PyErr_SetString(PyExc_RuntimeError, |
---|
576 | "extrapolate_second_order could not obtain domain object from quantity"); |
---|
577 | return NULL; |
---|
578 | } |
---|
579 | |
---|
580 | // Get pertinent variables |
---|
581 | centroids = get_consecutive_array(domain, "centroid_coordinates"); |
---|
582 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
583 | surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); |
---|
584 | number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); |
---|
585 | vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates"); |
---|
586 | vertex_values = get_consecutive_array(quantity, "vertex_values"); |
---|
587 | |
---|
588 | N = centroid_values -> dimensions[0]; |
---|
589 | |
---|
590 | // Release |
---|
591 | Py_DECREF(domain); |
---|
592 | |
---|
593 | //Allocate space for return vectors a and b (don't DECREF) |
---|
594 | dimensions[0] = N; |
---|
595 | a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
596 | b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
597 | |
---|
598 | //FIXME: Odd that I couldn't use normal arrays |
---|
599 | //Allocate space for return vectors a and b (don't DECREF) |
---|
600 | //a = (double*) malloc(N * sizeof(double)); |
---|
601 | //if (!a) return NULL; |
---|
602 | //b = (double*) malloc(N * sizeof(double)); |
---|
603 | //if (!b) return NULL; |
---|
604 | |
---|
605 | |
---|
606 | err = _compute_gradients(N, |
---|
607 | (double*) centroids -> data, |
---|
608 | (double*) centroid_values -> data, |
---|
609 | (long*) number_of_boundaries -> data, |
---|
610 | (long*) surrogate_neighbours -> data, |
---|
611 | (double*) a -> data, |
---|
612 | (double*) b -> data); |
---|
613 | |
---|
614 | if (err != 0) { |
---|
615 | PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); |
---|
616 | return NULL; |
---|
617 | } |
---|
618 | |
---|
619 | err = _extrapolate(N, |
---|
620 | (double*) centroids -> data, |
---|
621 | (double*) centroid_values -> data, |
---|
622 | (double*) vertex_coordinates -> data, |
---|
623 | (double*) vertex_values -> data, |
---|
624 | (double*) a -> data, |
---|
625 | (double*) b -> data); |
---|
626 | |
---|
627 | |
---|
628 | if (err != 0) { |
---|
629 | PyErr_SetString(PyExc_RuntimeError, |
---|
630 | "Internal function _extrapolate failed"); |
---|
631 | return NULL; |
---|
632 | } |
---|
633 | |
---|
634 | |
---|
635 | |
---|
636 | // Release |
---|
637 | Py_DECREF(centroids); |
---|
638 | Py_DECREF(centroid_values); |
---|
639 | Py_DECREF(number_of_boundaries); |
---|
640 | Py_DECREF(surrogate_neighbours); |
---|
641 | Py_DECREF(vertex_coordinates); |
---|
642 | Py_DECREF(vertex_values); |
---|
643 | Py_DECREF(a); |
---|
644 | Py_DECREF(b); |
---|
645 | |
---|
646 | return Py_BuildValue(""); |
---|
647 | } |
---|
648 | |
---|
649 | |
---|
650 | |
---|
651 | PyObject *limit(PyObject *self, PyObject *args) { |
---|
652 | |
---|
653 | PyObject *quantity, *domain, *Tmp; |
---|
654 | PyArrayObject |
---|
655 | *qv, //Conserved quantities at vertices |
---|
656 | *qc, //Conserved quantities at centroids |
---|
657 | *neighbours; |
---|
658 | |
---|
659 | int k, i, n, N, k3; |
---|
660 | double beta_w; //Safety factor |
---|
661 | double *qmin, *qmax, qn; |
---|
662 | |
---|
663 | // Convert Python arguments to C |
---|
664 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
665 | PyErr_SetString(PyExc_RuntimeError, |
---|
666 | "quantity_ext.c: limit could not parse input"); |
---|
667 | return NULL; |
---|
668 | } |
---|
669 | |
---|
670 | domain = PyObject_GetAttrString(quantity, "domain"); |
---|
671 | if (!domain) { |
---|
672 | PyErr_SetString(PyExc_RuntimeError, |
---|
673 | "quantity_ext.c: limit could not obtain domain object from quantity"); |
---|
674 | |
---|
675 | return NULL; |
---|
676 | } |
---|
677 | |
---|
678 | //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours"); |
---|
679 | neighbours = get_consecutive_array(domain, "neighbours"); |
---|
680 | |
---|
681 | // Get safety factor beta_w |
---|
682 | Tmp = PyObject_GetAttrString(domain, "beta_w"); |
---|
683 | if (!Tmp) { |
---|
684 | PyErr_SetString(PyExc_RuntimeError, |
---|
685 | "quantity_ext.c: limit could not obtain beta_w object from domain"); |
---|
686 | |
---|
687 | return NULL; |
---|
688 | } |
---|
689 | |
---|
690 | beta_w = PyFloat_AsDouble(Tmp); |
---|
691 | |
---|
692 | Py_DECREF(Tmp); |
---|
693 | Py_DECREF(domain); |
---|
694 | |
---|
695 | |
---|
696 | qc = get_consecutive_array(quantity, "centroid_values"); |
---|
697 | qv = get_consecutive_array(quantity, "vertex_values"); |
---|
698 | |
---|
699 | |
---|
700 | N = qc -> dimensions[0]; |
---|
701 | |
---|
702 | // Find min and max of this and neighbour's centroid values |
---|
703 | qmin = malloc(N * sizeof(double)); |
---|
704 | qmax = malloc(N * sizeof(double)); |
---|
705 | for (k=0; k<N; k++) { |
---|
706 | k3=k*3; |
---|
707 | |
---|
708 | qmin[k] = ((double*) qc -> data)[k]; |
---|
709 | qmax[k] = qmin[k]; |
---|
710 | |
---|
711 | for (i=0; i<3; i++) { |
---|
712 | n = ((long*) neighbours -> data)[k3+i]; |
---|
713 | if (n >= 0) { |
---|
714 | qn = ((double*) qc -> data)[n]; //Neighbour's centroid value |
---|
715 | |
---|
716 | qmin[k] = min(qmin[k], qn); |
---|
717 | qmax[k] = max(qmax[k], qn); |
---|
718 | } |
---|
719 | //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]); |
---|
720 | //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]); |
---|
721 | } |
---|
722 | } |
---|
723 | |
---|
724 | // Call underlying routine |
---|
725 | _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax); |
---|
726 | |
---|
727 | free(qmin); |
---|
728 | free(qmax); |
---|
729 | return Py_BuildValue(""); |
---|
730 | } |
---|
731 | |
---|
732 | |
---|
733 | |
---|
734 | // Method table for python module |
---|
735 | static struct PyMethodDef MethodTable[] = { |
---|
736 | {"limit", limit, METH_VARARGS, "Print out"}, |
---|
737 | {"update", update, METH_VARARGS, "Print out"}, |
---|
738 | {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"}, |
---|
739 | {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"}, |
---|
740 | {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"}, |
---|
741 | {"extrapolate_second_order", extrapolate_second_order, |
---|
742 | METH_VARARGS, "Print out"}, |
---|
743 | {"interpolate_from_vertices_to_edges", |
---|
744 | interpolate_from_vertices_to_edges, |
---|
745 | METH_VARARGS, "Print out"}, |
---|
746 | {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"}, |
---|
747 | {NULL, NULL, 0, NULL} // sentinel |
---|
748 | }; |
---|
749 | |
---|
750 | // Module initialisation |
---|
751 | void initquantity_ext(void){ |
---|
752 | Py_InitModule("quantity_ext", MethodTable); |
---|
753 | |
---|
754 | import_array(); // Necessary for handling of NumPY structures |
---|
755 | } |
---|