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