Changeset 4728
- Timestamp:
- Sep 11, 2007, 5:27:16 PM (18 years ago)
- Location:
- anuga_core/source/anuga/abstract_2d_finite_volumes
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r4723 r4728 150 150 Pre-condition: vertex_values have been set 151 151 """ 152 152 # FIXME (Ole): This function is not called 153 # in real model runs. Consider removing it. 154 153 155 N = self.vertex_values.shape[0] 154 156 for i in range(N): -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
r4712 r4728 36 36 37 37 if (number_of_boundaries[k] < 2) { 38 // Two or three true neighbours39 40 // Get indices of neighbours (or self when used as surrogate)41 // k0, k1, k2 = surrogate_neighbours[k,:]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 42 43 43 k0 = surrogate_neighbours[index3 + 0]; … … 48 48 if (k0 == k1 || k1 == k2) return -1; 49 49 50 // Get data50 // Get data 51 51 q0 = centroid_values[k0]; 52 52 q1 = centroid_values[k1]; … … 57 57 x2 = centroids[k2*2]; y2 = centroids[k2*2+1]; 58 58 59 // Gradient59 // Gradient 60 60 _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]); 61 61 62 62 } else if (number_of_boundaries[k] == 2) { 63 // One true neighbour64 65 // Get index of the one neighbour63 // One true neighbour 64 65 // Get index of the one neighbour 66 66 i=0; k0 = k; 67 67 while (i<3 && k0==k) { … … 73 73 k1 = k; //self 74 74 75 // Get data75 // Get data 76 76 q0 = centroid_values[k0]; 77 77 q1 = centroid_values[k1]; … … 80 80 x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; 81 81 82 // Two point gradient82 // Two point gradient 83 83 _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]); 84 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 85 } 93 86 // else: 94 87 // #No true neighbours - 95 88 // #Fall back to first order scheme 96 // pass97 89 } 98 90 return 0; … … 117 109 k2 = 2*k; 118 110 119 // Centroid coordinates111 // Centroid coordinates 120 112 x = centroids[k2]; y = centroids[k2+1]; 121 113 122 // vertex coordinates123 // x0, y0, x1, y1, x2, y2 = X[k,:]114 // vertex coordinates 115 // x0, y0, x1, y1, x2, y2 = X[k,:] 124 116 x0 = vertex_coordinates[k6 + 0]; 125 117 y0 = vertex_coordinates[k6 + 1]; … … 129 121 y2 = vertex_coordinates[k6 + 5]; 130 122 131 // Extrapolate123 // Extrapolate 132 124 vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y); 133 125 vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y); … … 156 148 q2 = vertex_values[k3 + 2]; 157 149 158 //printf("%f, %f, %f\n", q0, q1, q2);159 150 edge_values[k3 + 0] = 0.5*(q1+q2); 160 151 edge_values[k3 + 1] = 0.5*(q0+q2); … … 167 158 double* centroid_values, 168 159 double* centroid_backup_values) { 169 // Backup centroid values160 // Backup centroid values 170 161 171 162 … … 186 177 double* centroid_values, 187 178 double* centroid_backup_values) { 188 // saxby centroid values179 // Saxby centroid values 189 180 190 181 … … 206 197 double* explicit_update, 207 198 double* semi_implicit_update) { 208 // Update centroid values based on values stored in209 // explicit_update and semi_implicit_update as well as given timestep199 // Update centroid values based on values stored in 200 // explicit_update and semi_implicit_update as well as given timestep 210 201 211 202 … … 214 205 215 206 216 // Divide semi_implicit update by conserved quantity207 // Divide semi_implicit update by conserved quantity 217 208 for (k=0; k<N; k++) { 218 209 x = centroid_values[k]; … … 225 216 226 217 227 // Semi implicit updates218 // Semi implicit updates 228 219 for (k=0; k<N; k++) { 229 220 denominator = 1.0 - timestep*semi_implicit_update[k]; … … 236 227 } 237 228 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 229 230 // Explicit updates 244 231 for (k=0; k<N; k++) { 245 232 centroid_values[k] += timestep*explicit_update[k]; … … 247 234 248 235 249 // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py236 // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py 250 237 for (k=0;k<N;k++){ 251 238 semi_implicit_update[k]=0.0; … … 261 248 double* vertex_values, 262 249 double* A) { 263 // Average vertex values to obtain one value per node250 // Average vertex values to obtain one value per node 264 251 265 252 int i, index; … … 270 257 for (i=0; i<N; i++) { 271 258 272 // if (current_node == N) {273 // printf("Current node exceeding number of nodes (%d)", N);274 // return 1;259 // if (current_node == N) { 260 // printf("Current node exceeding number of nodes (%d)", N); 261 // return 1; 275 262 // } 276 263 … … 278 265 k += 1; 279 266 280 // volume_id = index / 3281 // vertex_id = index % 3282 // total += self.vertex_values[volume_id, vertex_id]267 // volume_id = index / 3 268 // vertex_id = index % 3 269 // total += self.vertex_values[volume_id, vertex_id] 283 270 total += vertex_values[index]; 284 271 285 // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);272 // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total); 286 273 if (number_of_triangles_per_node[current_node] == k) { 287 274 A[current_node] = total/k; … … 334 321 } 335 322 336 // Release and return323 // Release and return 337 324 Py_DECREF(centroid_values); 338 325 Py_DECREF(explicit_update); … … 368 355 369 356 370 // Release and return357 // Release and return 371 358 Py_DECREF(centroid_values); 372 359 Py_DECREF(centroid_backup_values); … … 401 388 402 389 403 // Release and return390 // Release and return 404 391 Py_DECREF(centroid_values); 405 392 Py_DECREF(centroid_backup_values); … … 438 425 } 439 426 440 // Release and return427 // Release and return 441 428 Py_DECREF(vertex_values); 442 429 Py_DECREF(edge_values); … … 514 501 } 515 502 516 // Get pertinent variables503 // Get pertinent variables 517 504 518 505 centroids = get_consecutive_array(domain, "centroid_coordinates"); … … 523 510 N = centroid_values -> dimensions[0]; 524 511 525 // Release512 // Release 526 513 Py_DECREF(domain); 527 514 528 // Allocate space for return vectors a and b (don't DECREF)515 // Allocate space for return vectors a and b (don't DECREF) 529 516 dimensions[0] = N; 530 517 a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); … … 546 533 } 547 534 548 // Release535 // Release 549 536 Py_DECREF(centroids); 550 537 Py_DECREF(centroid_values); … … 552 539 Py_DECREF(surrogate_neighbours); 553 540 554 // Build result, release and return541 // Build result, release and return 555 542 R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); 556 543 Py_DECREF(a); … … 591 578 } 592 579 593 // Get pertinent variables580 // Get pertinent variables 594 581 centroids = get_consecutive_array(domain, "centroid_coordinates"); 595 582 centroid_values = get_consecutive_array(quantity, "centroid_values"); … … 601 588 N = centroid_values -> dimensions[0]; 602 589 603 // Release590 // Release 604 591 Py_DECREF(domain); 605 592 … … 647 634 648 635 649 // Release636 // Release 650 637 Py_DECREF(centroids); 651 638 Py_DECREF(centroid_values); … … 692 679 neighbours = get_consecutive_array(domain, "neighbours"); 693 680 694 // Get safety factor beta_w681 // Get safety factor beta_w 695 682 Tmp = PyObject_GetAttrString(domain, "beta_w"); 696 683 if (!Tmp) { … … 713 700 N = qc -> dimensions[0]; 714 701 715 // Find min and max of this and neighbour's centroid values702 // Find min and max of this and neighbour's centroid values 716 703 qmin = malloc(N * sizeof(double)); 717 704 qmax = malloc(N * sizeof(double)); … … 765 752 Py_InitModule("quantity_ext", MethodTable); 766 753 767 import_array(); //Necessary for handling of NumPY structures768 } 754 import_array(); // Necessary for handling of NumPY structures 755 }
Note: See TracChangeset
for help on using the changeset viewer.