Changeset 7685 for anuga_core/source/anuga/fit_interpolate/interpolate.py
- Timestamp:
- Apr 14, 2010, 6:53:12 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/interpolate.py
r7682 r7685 155 155 use_cache=use_cache, 156 156 verbose=verbose, 157 157 output_centroids=output_centroids) 158 158 159 159 return result … … 230 230 start_blocking_len=500000, 231 231 verbose=False, 232 232 output_centroids=False): 233 233 """Interpolate mesh data f to determine values, z, at points. 234 234 … … 324 324 # @param use_cache True if caching should be used to accelerate the calculations 325 325 # @param verbose True if this function is verbose. 326 # @return ??326 # @return interpolated f 327 327 def interpolate_block(self, f, point_coordinates, 328 328 use_cache=False, verbose=False, output_centroids=False): … … 407 407 ## 408 408 # @brief Get interpolated data at given points. 409 410 411 409 # Applies a transform to all points to calculate the 410 # interpolated values. Points outside the mesh are returned as NaN. 411 # @note self._A matrix must be valid 412 412 # @param f Array of arbitrary data 413 413 # @param verbose True if this function is to be verbose. … … 431 431 # @brief Build NxM interpolation matrix. 432 432 # @param point_coordinates Points to sample at 433 434 435 433 # @param output_centroids set to True to always sample from the centre 434 # of the intersected triangle, instead of the intersection 435 # point. 436 436 # @param verbose True if this function is to be verbose. 437 437 # @return Interpolation matrix A, plus lists of the points inside and outside the mesh 438 438 # and the list of centroids, if requested. 439 439 def _build_interpolation_matrix_A(self, 440 440 point_coordinates, … … 464 464 if verbose: log.critical('Getting indices inside mesh boundary') 465 465 466 inside_poly_indices, outside_poly_indices = \ 466 # Quick test against boundary, but will not deal with holes in the mesh 467 inside_boundary_indices, outside_poly_indices = \ 467 468 in_and_outside_polygon(point_coordinates, 468 469 self.mesh.get_boundary_polygon(), … … 474 475 475 476 # Since you can block, throw a warning, not an error. 476 if verbose and 0 == len(inside_ poly_indices):477 if verbose and 0 == len(inside_boundary_indices): 477 478 log.critical('WARNING: No points within the mesh!') 478 479 … … 485 486 A = Sparse(n,m) 486 487 487 n = len(inside_ poly_indices)488 n = len(inside_boundary_indices) 488 489 489 490 centroids = [] 490 491 492 inside_poly_indices = [] 493 491 494 # Compute matrix elements for points inside the mesh 492 495 if verbose: log.critical('Building interpolation matrix from %d points' 493 496 % n) 494 497 495 for d, i in enumerate(inside_ poly_indices):498 for d, i in enumerate(inside_boundary_indices): 496 499 # For each data_coordinate point 497 500 if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d' … … 501 504 element_found, sigma0, sigma1, sigma2, k = \ 502 505 search_tree_of_vertices(self.root, x) 506 507 # Update interpolation matrix A if necessary 508 if element_found is True: 509 510 if verbose: 511 print 'Point is within mesh:', d, i 503 512 504 # Update interpolation matrix A if necessary 505 if element_found is True:513 inside_poly_indices.append(i) 514 506 515 # Assign values to matrix A 507 516 j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0 … … 509 518 j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2 510 519 js = [j0, j1, j2] 511 520 512 521 if output_centroids is False: 513 522 # Weight each vertex according to its distance from x … … 516 525 A[i, j] = sigmas[j] 517 526 else: 518 527 # If centroids are needed, weight all 3 vertices equally 519 528 for j in js: 520 529 A[i, j] = 1.0/3.0 521 centroids.append(self.mesh.centroid_coordinates[k]) 530 centroids.append(self.mesh.centroid_coordinates[k]) 522 531 else: 523 # Mesh has a hole - move this point to the outside list524 # num.delete(inside_poly_indices,[i], axis=0)525 num.append(outside_poly_indices, [i], axis=0) 526 # msg = 'Could not find triangle for point', x 527 # raise Exception(msg)532 if verbose: 533 log.critical('Mesh has a hole - moving this point to outside list') 534 535 # This is a numpy arrays, so we need to do a slow transfer 536 outside_poly_indices = num.append(outside_poly_indices, [i], axis=0) 528 537 529 538 return A, inside_poly_indices, outside_poly_indices, centroids … … 775 784 verbose=False, 776 785 gauge_neighbour_id=None, 777 786 output_centroids=False): 778 787 """Initialise object and build spatial interpolation if required 779 788 … … 1003 1012 verbose=False, 1004 1013 output_centroids=output_centroids) 1005 self.centroids = interpol.centroids 1014 self.centroids = interpol.centroids 1006 1015 elif triangles is None and vertex_coordinates is not None: 1007 1016 result = interpolate_polyline(Q, … … 1012 1021 1013 1022 #assert len(result), len(interpolation_points) 1014 self.precomputed_values[name][i, :] = result 1015 1023 self.precomputed_values[name][i, :] = result 1024 1016 1025 # Report 1017 1026 if verbose: 1018 log.critical(self.statistics()) 1027 log.critical(self.statistics()) 1019 1028 else: 1020 1029 # Store quantitites as is
Note: See TracChangeset
for help on using the changeset viewer.
