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* edge_values, |
---|
101 | double* a, |
---|
102 | double* b) { |
---|
103 | |
---|
104 | int k, k2, k3, k6; |
---|
105 | double x, y, x0, y0, x1, y1, x2, y2; |
---|
106 | |
---|
107 | for (k=0; k<N; k++){ |
---|
108 | k6 = 6*k; |
---|
109 | k3 = 3*k; |
---|
110 | k2 = 2*k; |
---|
111 | |
---|
112 | // Centroid coordinates |
---|
113 | x = centroids[k2]; y = centroids[k2+1]; |
---|
114 | |
---|
115 | // vertex coordinates |
---|
116 | // x0, y0, x1, y1, x2, y2 = X[k,:] |
---|
117 | x0 = vertex_coordinates[k6 + 0]; |
---|
118 | y0 = vertex_coordinates[k6 + 1]; |
---|
119 | x1 = vertex_coordinates[k6 + 2]; |
---|
120 | y1 = vertex_coordinates[k6 + 3]; |
---|
121 | x2 = vertex_coordinates[k6 + 4]; |
---|
122 | y2 = vertex_coordinates[k6 + 5]; |
---|
123 | |
---|
124 | // Extrapolate |
---|
125 | vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y); |
---|
126 | vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y); |
---|
127 | vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y); |
---|
128 | |
---|
129 | |
---|
130 | edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]); |
---|
131 | edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]); |
---|
132 | edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]); |
---|
133 | |
---|
134 | } |
---|
135 | return 0; |
---|
136 | } |
---|
137 | |
---|
138 | |
---|
139 | |
---|
140 | |
---|
141 | int _interpolate(int N, |
---|
142 | double* vertex_values, |
---|
143 | double* edge_values) { |
---|
144 | |
---|
145 | int k, k3; |
---|
146 | double q0, q1, q2; |
---|
147 | |
---|
148 | |
---|
149 | for (k=0; k<N; k++) { |
---|
150 | k3 = 3*k; |
---|
151 | |
---|
152 | q0 = vertex_values[k3 + 0]; |
---|
153 | q1 = vertex_values[k3 + 1]; |
---|
154 | q2 = vertex_values[k3 + 2]; |
---|
155 | |
---|
156 | edge_values[k3 + 0] = 0.5*(q1+q2); |
---|
157 | edge_values[k3 + 1] = 0.5*(q0+q2); |
---|
158 | edge_values[k3 + 2] = 0.5*(q0+q1); |
---|
159 | } |
---|
160 | return 0; |
---|
161 | } |
---|
162 | |
---|
163 | int _backup_centroid_values(int N, |
---|
164 | double* centroid_values, |
---|
165 | double* centroid_backup_values) { |
---|
166 | // Backup centroid values |
---|
167 | |
---|
168 | |
---|
169 | int k; |
---|
170 | |
---|
171 | for (k=0; k<N; k++) { |
---|
172 | centroid_backup_values[k] = centroid_values[k]; |
---|
173 | } |
---|
174 | |
---|
175 | |
---|
176 | return 0; |
---|
177 | } |
---|
178 | |
---|
179 | |
---|
180 | int _saxpy_centroid_values(int N, |
---|
181 | double a, |
---|
182 | double b, |
---|
183 | double* centroid_values, |
---|
184 | double* centroid_backup_values) { |
---|
185 | // Saxby centroid values |
---|
186 | |
---|
187 | |
---|
188 | int k; |
---|
189 | |
---|
190 | |
---|
191 | for (k=0; k<N; k++) { |
---|
192 | centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k]; |
---|
193 | } |
---|
194 | |
---|
195 | |
---|
196 | return 0; |
---|
197 | } |
---|
198 | |
---|
199 | |
---|
200 | int _update(int N, |
---|
201 | double timestep, |
---|
202 | double* centroid_values, |
---|
203 | double* explicit_update, |
---|
204 | double* semi_implicit_update) { |
---|
205 | // Update centroid values based on values stored in |
---|
206 | // explicit_update and semi_implicit_update as well as given timestep |
---|
207 | |
---|
208 | |
---|
209 | int k; |
---|
210 | double denominator, x; |
---|
211 | |
---|
212 | |
---|
213 | // Divide semi_implicit update by conserved quantity |
---|
214 | for (k=0; k<N; k++) { |
---|
215 | x = centroid_values[k]; |
---|
216 | if (x == 0.0) { |
---|
217 | semi_implicit_update[k] = 0.0; |
---|
218 | } else { |
---|
219 | semi_implicit_update[k] /= x; |
---|
220 | } |
---|
221 | } |
---|
222 | |
---|
223 | |
---|
224 | // Semi implicit updates |
---|
225 | for (k=0; k<N; k++) { |
---|
226 | denominator = 1.0 - timestep*semi_implicit_update[k]; |
---|
227 | if (denominator == 0.0) { |
---|
228 | return -1; |
---|
229 | } else { |
---|
230 | //Update conserved_quantities from semi implicit updates |
---|
231 | centroid_values[k] /= denominator; |
---|
232 | } |
---|
233 | } |
---|
234 | |
---|
235 | |
---|
236 | // Explicit updates |
---|
237 | for (k=0; k<N; k++) { |
---|
238 | centroid_values[k] += timestep*explicit_update[k]; |
---|
239 | } |
---|
240 | |
---|
241 | |
---|
242 | // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py |
---|
243 | for (k=0;k<N;k++){ |
---|
244 | semi_implicit_update[k]=0.0; |
---|
245 | } |
---|
246 | |
---|
247 | return 0; |
---|
248 | } |
---|
249 | |
---|
250 | |
---|
251 | int _average_vertex_values(int N, |
---|
252 | long* vertex_value_indices, |
---|
253 | long* number_of_triangles_per_node, |
---|
254 | double* vertex_values, |
---|
255 | double* A) { |
---|
256 | // Average vertex values to obtain one value per node |
---|
257 | |
---|
258 | int i, index; |
---|
259 | int k = 0; //Track triangles touching each node |
---|
260 | int current_node = 0; |
---|
261 | double total = 0.0; |
---|
262 | |
---|
263 | for (i=0; i<N; i++) { |
---|
264 | |
---|
265 | // if (current_node == N) { |
---|
266 | // printf("Current node exceeding number of nodes (%d)", N); |
---|
267 | // return 1; |
---|
268 | // } |
---|
269 | |
---|
270 | index = vertex_value_indices[i]; |
---|
271 | k += 1; |
---|
272 | |
---|
273 | // volume_id = index / 3 |
---|
274 | // vertex_id = index % 3 |
---|
275 | // total += self.vertex_values[volume_id, vertex_id] |
---|
276 | total += vertex_values[index]; |
---|
277 | |
---|
278 | // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total); |
---|
279 | if (number_of_triangles_per_node[current_node] == k) { |
---|
280 | A[current_node] = total/k; |
---|
281 | |
---|
282 | // Move on to next node |
---|
283 | total = 0.0; |
---|
284 | k = 0; |
---|
285 | current_node += 1; |
---|
286 | } |
---|
287 | } |
---|
288 | |
---|
289 | return 0; |
---|
290 | } |
---|
291 | |
---|
292 | |
---|
293 | ///////////////////////////////////////////////// |
---|
294 | // Gateways to Python |
---|
295 | PyObject *update(PyObject *self, PyObject *args) { |
---|
296 | // FIXME (Ole): It would be great to turn this text into a Python DOC string |
---|
297 | |
---|
298 | /*"""Update centroid values based on values stored in |
---|
299 | explicit_update and semi_implicit_update as well as given timestep |
---|
300 | |
---|
301 | Function implementing forcing terms must take on argument |
---|
302 | which is the domain and they must update either explicit |
---|
303 | or implicit updates, e,g,: |
---|
304 | |
---|
305 | def gravity(domain): |
---|
306 | .... |
---|
307 | domain.quantities['xmomentum'].explicit_update = ... |
---|
308 | domain.quantities['ymomentum'].explicit_update = ... |
---|
309 | |
---|
310 | |
---|
311 | |
---|
312 | Explicit terms must have the form |
---|
313 | |
---|
314 | G(q, t) |
---|
315 | |
---|
316 | and explicit scheme is |
---|
317 | |
---|
318 | q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t) |
---|
319 | |
---|
320 | |
---|
321 | Semi implicit forcing terms are assumed to have the form |
---|
322 | |
---|
323 | G(q, t) = H(q, t) q |
---|
324 | |
---|
325 | and the semi implicit scheme will then be |
---|
326 | |
---|
327 | q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1}) |
---|
328 | |
---|
329 | */ |
---|
330 | |
---|
331 | PyObject *quantity; |
---|
332 | PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update; |
---|
333 | |
---|
334 | double timestep; |
---|
335 | int N, err; |
---|
336 | |
---|
337 | |
---|
338 | // Convert Python arguments to C |
---|
339 | if (!PyArg_ParseTuple(args, "Od", &quantity, ×tep)) { |
---|
340 | PyErr_SetString(PyExc_RuntimeError, |
---|
341 | "quantity_ext.c: update could not parse input"); |
---|
342 | return NULL; |
---|
343 | } |
---|
344 | |
---|
345 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
346 | explicit_update = get_consecutive_array(quantity, "explicit_update"); |
---|
347 | semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update"); |
---|
348 | |
---|
349 | N = centroid_values -> dimensions[0]; |
---|
350 | |
---|
351 | err = _update(N, timestep, |
---|
352 | (double*) centroid_values -> data, |
---|
353 | (double*) explicit_update -> data, |
---|
354 | (double*) semi_implicit_update -> data); |
---|
355 | |
---|
356 | |
---|
357 | if (err != 0) { |
---|
358 | PyErr_SetString(PyExc_RuntimeError, |
---|
359 | "Zero division in semi implicit update - call Stephen :)"); |
---|
360 | return NULL; |
---|
361 | } |
---|
362 | |
---|
363 | // Release and return |
---|
364 | Py_DECREF(centroid_values); |
---|
365 | Py_DECREF(explicit_update); |
---|
366 | Py_DECREF(semi_implicit_update); |
---|
367 | |
---|
368 | return Py_BuildValue(""); |
---|
369 | } |
---|
370 | |
---|
371 | |
---|
372 | PyObject *backup_centroid_values(PyObject *self, PyObject *args) { |
---|
373 | |
---|
374 | PyObject *quantity; |
---|
375 | PyArrayObject *centroid_values, *centroid_backup_values; |
---|
376 | |
---|
377 | int N, err; |
---|
378 | |
---|
379 | |
---|
380 | // Convert Python arguments to C |
---|
381 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
382 | PyErr_SetString(PyExc_RuntimeError, |
---|
383 | "quantity_ext.c: backup_centroid_values could not parse input"); |
---|
384 | return NULL; |
---|
385 | } |
---|
386 | |
---|
387 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
388 | centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values"); |
---|
389 | |
---|
390 | N = centroid_values -> dimensions[0]; |
---|
391 | |
---|
392 | err = _backup_centroid_values(N, |
---|
393 | (double*) centroid_values -> data, |
---|
394 | (double*) centroid_backup_values -> data); |
---|
395 | |
---|
396 | |
---|
397 | // Release and return |
---|
398 | Py_DECREF(centroid_values); |
---|
399 | Py_DECREF(centroid_backup_values); |
---|
400 | |
---|
401 | return Py_BuildValue(""); |
---|
402 | } |
---|
403 | |
---|
404 | PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) { |
---|
405 | |
---|
406 | PyObject *quantity; |
---|
407 | PyArrayObject *centroid_values, *centroid_backup_values; |
---|
408 | |
---|
409 | double a,b; |
---|
410 | int N, err; |
---|
411 | |
---|
412 | |
---|
413 | // Convert Python arguments to C |
---|
414 | if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) { |
---|
415 | PyErr_SetString(PyExc_RuntimeError, |
---|
416 | "quantity_ext.c: saxpy_centroid_values could not parse input"); |
---|
417 | return NULL; |
---|
418 | } |
---|
419 | |
---|
420 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
421 | centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values"); |
---|
422 | |
---|
423 | N = centroid_values -> dimensions[0]; |
---|
424 | |
---|
425 | err = _saxpy_centroid_values(N,a,b, |
---|
426 | (double*) centroid_values -> data, |
---|
427 | (double*) centroid_backup_values -> data); |
---|
428 | |
---|
429 | |
---|
430 | // Release and return |
---|
431 | Py_DECREF(centroid_values); |
---|
432 | Py_DECREF(centroid_backup_values); |
---|
433 | |
---|
434 | return Py_BuildValue(""); |
---|
435 | } |
---|
436 | |
---|
437 | |
---|
438 | PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) { |
---|
439 | // |
---|
440 | //Compute edge values from vertex values using linear interpolation |
---|
441 | // |
---|
442 | |
---|
443 | PyObject *quantity; |
---|
444 | PyArrayObject *vertex_values, *edge_values; |
---|
445 | |
---|
446 | int N, err; |
---|
447 | |
---|
448 | // Convert Python arguments to C |
---|
449 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
450 | PyErr_SetString(PyExc_RuntimeError, |
---|
451 | "quantity_ext.c: interpolate_from_vertices_to_edges could not parse input"); |
---|
452 | return NULL; |
---|
453 | } |
---|
454 | |
---|
455 | vertex_values = get_consecutive_array(quantity, "vertex_values"); |
---|
456 | edge_values = get_consecutive_array(quantity, "edge_values"); |
---|
457 | |
---|
458 | N = vertex_values -> dimensions[0]; |
---|
459 | |
---|
460 | err = _interpolate(N, |
---|
461 | (double*) vertex_values -> data, |
---|
462 | (double*) edge_values -> data); |
---|
463 | |
---|
464 | if (err != 0) { |
---|
465 | PyErr_SetString(PyExc_RuntimeError, |
---|
466 | "Interpolate could not be computed"); |
---|
467 | return NULL; |
---|
468 | } |
---|
469 | |
---|
470 | // Release and return |
---|
471 | Py_DECREF(vertex_values); |
---|
472 | Py_DECREF(edge_values); |
---|
473 | |
---|
474 | return Py_BuildValue(""); |
---|
475 | } |
---|
476 | |
---|
477 | |
---|
478 | PyObject *average_vertex_values(PyObject *self, PyObject *args) { |
---|
479 | |
---|
480 | PyArrayObject |
---|
481 | *vertex_value_indices, |
---|
482 | *number_of_triangles_per_node, |
---|
483 | *vertex_values, |
---|
484 | *A; |
---|
485 | |
---|
486 | |
---|
487 | int N, err; |
---|
488 | |
---|
489 | // Convert Python arguments to C |
---|
490 | if (!PyArg_ParseTuple(args, "OOOO", |
---|
491 | &vertex_value_indices, |
---|
492 | &number_of_triangles_per_node, |
---|
493 | &vertex_values, |
---|
494 | &A)) { |
---|
495 | PyErr_SetString(PyExc_RuntimeError, |
---|
496 | "quantity_ext.c: average_vertex_values could not parse input"); |
---|
497 | return NULL; |
---|
498 | } |
---|
499 | |
---|
500 | N = vertex_value_indices -> dimensions[0]; |
---|
501 | // printf("Got parameters, N=%d\n", N); |
---|
502 | err = _average_vertex_values(N, |
---|
503 | (long*) vertex_value_indices -> data, |
---|
504 | (long*) number_of_triangles_per_node -> data, |
---|
505 | (double*) vertex_values -> data, |
---|
506 | (double*) A -> data); |
---|
507 | |
---|
508 | //printf("Error %d", err); |
---|
509 | if (err != 0) { |
---|
510 | PyErr_SetString(PyExc_RuntimeError, |
---|
511 | "average_vertex_values could not be computed"); |
---|
512 | return NULL; |
---|
513 | } |
---|
514 | |
---|
515 | return Py_BuildValue(""); |
---|
516 | } |
---|
517 | |
---|
518 | |
---|
519 | |
---|
520 | PyObject *compute_gradients(PyObject *self, PyObject *args) { |
---|
521 | //"""Compute gradients of triangle surfaces defined by centroids of |
---|
522 | //neighbouring volumes. |
---|
523 | //If one edge is on the boundary, use own centroid as neighbour centroid. |
---|
524 | //If two or more are on the boundary, fall back to first order scheme. |
---|
525 | //""" |
---|
526 | |
---|
527 | |
---|
528 | PyObject *quantity, *domain, *R; |
---|
529 | PyArrayObject |
---|
530 | *centroids, //Coordinates at centroids |
---|
531 | *centroid_values, //Values at centroids |
---|
532 | *number_of_boundaries, //Number of boundaries for each triangle |
---|
533 | *surrogate_neighbours, //True neighbours or - if one missing - self |
---|
534 | *a, *b; //Return values |
---|
535 | |
---|
536 | int dimensions[1], N, err; |
---|
537 | |
---|
538 | // Convert Python arguments to C |
---|
539 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
540 | PyErr_SetString(PyExc_RuntimeError, |
---|
541 | "quantity_ext.c: compute_gradients could not parse input"); |
---|
542 | return NULL; |
---|
543 | } |
---|
544 | |
---|
545 | domain = PyObject_GetAttrString(quantity, "domain"); |
---|
546 | if (!domain) { |
---|
547 | PyErr_SetString(PyExc_RuntimeError, |
---|
548 | "compute_gradients could not obtain domain object from quantity"); |
---|
549 | return NULL; |
---|
550 | } |
---|
551 | |
---|
552 | // Get pertinent variables |
---|
553 | |
---|
554 | centroids = get_consecutive_array(domain, "centroid_coordinates"); |
---|
555 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
556 | surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); |
---|
557 | number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); |
---|
558 | |
---|
559 | N = centroid_values -> dimensions[0]; |
---|
560 | |
---|
561 | // Release |
---|
562 | Py_DECREF(domain); |
---|
563 | |
---|
564 | // Allocate space for return vectors a and b (don't DECREF) |
---|
565 | dimensions[0] = N; |
---|
566 | a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
567 | b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
568 | |
---|
569 | |
---|
570 | |
---|
571 | err = _compute_gradients(N, |
---|
572 | (double*) centroids -> data, |
---|
573 | (double*) centroid_values -> data, |
---|
574 | (long*) number_of_boundaries -> data, |
---|
575 | (long*) surrogate_neighbours -> data, |
---|
576 | (double*) a -> data, |
---|
577 | (double*) b -> data); |
---|
578 | |
---|
579 | if (err != 0) { |
---|
580 | PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); |
---|
581 | return NULL; |
---|
582 | } |
---|
583 | |
---|
584 | // Release |
---|
585 | Py_DECREF(centroids); |
---|
586 | Py_DECREF(centroid_values); |
---|
587 | Py_DECREF(number_of_boundaries); |
---|
588 | Py_DECREF(surrogate_neighbours); |
---|
589 | |
---|
590 | // Build result, release and return |
---|
591 | R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); |
---|
592 | Py_DECREF(a); |
---|
593 | Py_DECREF(b); |
---|
594 | return R; |
---|
595 | } |
---|
596 | |
---|
597 | |
---|
598 | |
---|
599 | PyObject *extrapolate_second_order(PyObject *self, PyObject *args) { |
---|
600 | |
---|
601 | PyObject *quantity, *domain; |
---|
602 | PyArrayObject |
---|
603 | *centroids, //Coordinates at centroids |
---|
604 | *centroid_values, //Values at centroids |
---|
605 | *vertex_coordinates, //Coordinates at vertices |
---|
606 | *vertex_values, //Values at vertices |
---|
607 | *edge_values, //Values at edges |
---|
608 | *number_of_boundaries, //Number of boundaries for each triangle |
---|
609 | *surrogate_neighbours, //True neighbours or - if one missing - self |
---|
610 | *a, *b; //Gradients |
---|
611 | |
---|
612 | //int N, err; |
---|
613 | int dimensions[1], N, err; |
---|
614 | //double *a, *b; //Gradients |
---|
615 | |
---|
616 | // Convert Python arguments to C |
---|
617 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
618 | PyErr_SetString(PyExc_RuntimeError, |
---|
619 | "extrapolate_second_order could not parse input"); |
---|
620 | return NULL; |
---|
621 | } |
---|
622 | |
---|
623 | domain = PyObject_GetAttrString(quantity, "domain"); |
---|
624 | if (!domain) { |
---|
625 | PyErr_SetString(PyExc_RuntimeError, |
---|
626 | "extrapolate_second_order could not obtain domain object from quantity"); |
---|
627 | return NULL; |
---|
628 | } |
---|
629 | |
---|
630 | // Get pertinent variables |
---|
631 | centroids = get_consecutive_array(domain, "centroid_coordinates"); |
---|
632 | centroid_values = get_consecutive_array(quantity, "centroid_values"); |
---|
633 | surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); |
---|
634 | number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); |
---|
635 | vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates"); |
---|
636 | vertex_values = get_consecutive_array(quantity, "vertex_values"); |
---|
637 | edge_values = get_consecutive_array(quantity, "edge_values"); |
---|
638 | |
---|
639 | N = centroid_values -> dimensions[0]; |
---|
640 | |
---|
641 | // Release |
---|
642 | Py_DECREF(domain); |
---|
643 | |
---|
644 | //Allocate space for return vectors a and b (don't DECREF) |
---|
645 | dimensions[0] = N; |
---|
646 | a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
647 | b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); |
---|
648 | |
---|
649 | //FIXME: Odd that I couldn't use normal arrays |
---|
650 | //Allocate space for return vectors a and b (don't DECREF) |
---|
651 | //a = (double*) malloc(N * sizeof(double)); |
---|
652 | //if (!a) return NULL; |
---|
653 | //b = (double*) malloc(N * sizeof(double)); |
---|
654 | //if (!b) return NULL; |
---|
655 | |
---|
656 | |
---|
657 | err = _compute_gradients(N, |
---|
658 | (double*) centroids -> data, |
---|
659 | (double*) centroid_values -> data, |
---|
660 | (long*) number_of_boundaries -> data, |
---|
661 | (long*) surrogate_neighbours -> data, |
---|
662 | (double*) a -> data, |
---|
663 | (double*) b -> data); |
---|
664 | |
---|
665 | if (err != 0) { |
---|
666 | PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); |
---|
667 | return NULL; |
---|
668 | } |
---|
669 | |
---|
670 | err = _extrapolate(N, |
---|
671 | (double*) centroids -> data, |
---|
672 | (double*) centroid_values -> data, |
---|
673 | (double*) vertex_coordinates -> data, |
---|
674 | (double*) vertex_values -> data, |
---|
675 | (double*) edge_values -> data, |
---|
676 | (double*) a -> data, |
---|
677 | (double*) b -> data); |
---|
678 | |
---|
679 | |
---|
680 | if (err != 0) { |
---|
681 | PyErr_SetString(PyExc_RuntimeError, |
---|
682 | "Internal function _extrapolate failed"); |
---|
683 | return NULL; |
---|
684 | } |
---|
685 | |
---|
686 | |
---|
687 | |
---|
688 | // Release |
---|
689 | Py_DECREF(centroids); |
---|
690 | Py_DECREF(centroid_values); |
---|
691 | Py_DECREF(number_of_boundaries); |
---|
692 | Py_DECREF(surrogate_neighbours); |
---|
693 | Py_DECREF(vertex_coordinates); |
---|
694 | Py_DECREF(vertex_values); |
---|
695 | Py_DECREF(edge_values); |
---|
696 | Py_DECREF(a); |
---|
697 | Py_DECREF(b); |
---|
698 | |
---|
699 | return Py_BuildValue(""); |
---|
700 | } |
---|
701 | |
---|
702 | |
---|
703 | |
---|
704 | PyObject *limit_old(PyObject *self, PyObject *args) { |
---|
705 | //Limit slopes for each volume to eliminate artificial variance |
---|
706 | //introduced by e.g. second order extrapolator |
---|
707 | |
---|
708 | //This is an unsophisticated limiter as it does not take into |
---|
709 | //account dependencies among quantities. |
---|
710 | |
---|
711 | //precondition: |
---|
712 | // vertex values are estimated from gradient |
---|
713 | //postcondition: |
---|
714 | // vertex values are updated |
---|
715 | // |
---|
716 | |
---|
717 | PyObject *quantity, *domain, *Tmp; |
---|
718 | PyArrayObject |
---|
719 | *qv, //Conserved quantities at vertices |
---|
720 | *qc, //Conserved quantities at centroids |
---|
721 | *neighbours; |
---|
722 | |
---|
723 | int k, i, n, N, k3; |
---|
724 | double beta_w; //Safety factor |
---|
725 | double *qmin, *qmax, qn; |
---|
726 | |
---|
727 | // Convert Python arguments to C |
---|
728 | if (!PyArg_ParseTuple(args, "O", &quantity)) { |
---|
729 | PyErr_SetString(PyExc_RuntimeError, |
---|
730 | "quantity_ext.c: limit_old could not parse input"); |
---|
731 | return NULL; |
---|
732 | } |
---|
733 | |
---|
734 | domain = PyObject_GetAttrString(quantity, "domain"); |
---|
735 | if (!domain) { |
---|
736 | PyErr_SetString(PyExc_RuntimeError, |
---|
737 | "quantity_ext.c: limit_old could not obtain domain object from quantity"); |
---|
738 | |
---|
739 | return NULL; |
---|
740 | } |
---|
741 | |
---|
742 | //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours"); |
---|
743 | neighbours = get_consecutive_array(domain, "neighbours"); |
---|
744 | |
---|
745 | // Get safety factor beta_w |
---|
746 | Tmp = PyObject_GetAttrString(domain, "beta_w"); |
---|
747 | if (!Tmp) { |
---|
748 | PyErr_SetString(PyExc_RuntimeError, |
---|
749 | "quantity_ext.c: limit_old could not obtain beta_w object from domain"); |
---|
750 | |
---|
751 | return NULL; |
---|
752 | } |
---|
753 | |
---|
754 | beta_w = PyFloat_AsDouble(Tmp); |
---|
755 | |
---|
756 | Py_DECREF(Tmp); |
---|
757 | Py_DECREF(domain); |
---|
758 | |
---|
759 | |
---|
760 | qc = get_consecutive_array(quantity, "centroid_values"); |
---|
761 | qv = get_consecutive_array(quantity, "vertex_values"); |
---|
762 | |
---|
763 | |
---|
764 | N = qc -> dimensions[0]; |
---|
765 | |
---|
766 | // Find min and max of this and neighbour's centroid values |
---|
767 | qmin = malloc(N * sizeof(double)); |
---|
768 | qmax = malloc(N * sizeof(double)); |
---|
769 | for (k=0; k<N; k++) { |
---|
770 | k3=k*3; |
---|
771 | |
---|
772 | qmin[k] = ((double*) qc -> data)[k]; |
---|
773 | qmax[k] = qmin[k]; |
---|
774 | |
---|
775 | for (i=0; i<3; i++) { |
---|
776 | n = ((long*) neighbours -> data)[k3+i]; |
---|
777 | if (n >= 0) { |
---|
778 | qn = ((double*) qc -> data)[n]; //Neighbour's centroid value |
---|
779 | |
---|
780 | qmin[k] = min(qmin[k], qn); |
---|
781 | qmax[k] = max(qmax[k], qn); |
---|
782 | } |
---|
783 | //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]); |
---|
784 | //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]); |
---|
785 | } |
---|
786 | } |
---|
787 | |
---|
788 | // Call underlying routine |
---|
789 | _limit_old(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax); |
---|
790 | |
---|
791 | free(qmin); |
---|
792 | free(qmax); |
---|
793 | return Py_BuildValue(""); |
---|
794 | } |
---|
795 | |
---|
796 | |
---|
797 | |
---|
798 | // Method table for python module |
---|
799 | static struct PyMethodDef MethodTable[] = { |
---|
800 | {"limit_old", limit_old, METH_VARARGS, "Print out"}, |
---|
801 | {"update", update, METH_VARARGS, "Print out"}, |
---|
802 | {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"}, |
---|
803 | {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"}, |
---|
804 | {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"}, |
---|
805 | {"extrapolate_second_order", extrapolate_second_order, |
---|
806 | METH_VARARGS, "Print out"}, |
---|
807 | {"interpolate_from_vertices_to_edges", |
---|
808 | interpolate_from_vertices_to_edges, |
---|
809 | METH_VARARGS, "Print out"}, |
---|
810 | {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"}, |
---|
811 | {NULL, NULL, 0, NULL} // sentinel |
---|
812 | }; |
---|
813 | |
---|
814 | // Module initialisation |
---|
815 | void initquantity_ext(void){ |
---|
816 | Py_InitModule("quantity_ext", MethodTable); |
---|
817 | |
---|
818 | import_array(); // Necessary for handling of NumPY structures |
---|
819 | } |
---|