Changeset 6226
- Timestamp:
- Jan 21, 2009, 5:28:57 PM (16 years ago)
- Location:
- anuga_core/source/anuga/abstract_2d_finite_volumes
- Files:
-
- 2 edited
-
domain.py (modified) (33 diffs)
-
quantity.py (modified) (70 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r6191 r6226 8 8 """ 9 9 10 import types 11 from time import time as walltime 12 10 13 from anuga.config import epsilon 11 14 from anuga.config import beta_euler, beta_rk2 12 13 15 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh 14 16 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ … … 22 24 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 23 25 import Transmissive_boundary 24 25 26 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain 26 27 from anuga.abstract_2d_finite_volumes.region\ 27 28 import Set_region as region_set_region 28 29 29 from anuga.utilities.polygon import inside_polygon 30 30 from anuga.abstract_2d_finite_volumes.util import get_textual_float 31 32 31 from quantity import Quantity 33 32 34 import types35 from time import time as walltime36 37 33 import Numeric as num 38 34 39 35 40 36 ## 41 # @brief Basic Domain class37 # @brief Generic Domain class 42 38 class Domain: 43 39 … … 49 45 # @param conserved_quantities List of names of quantities to be conserved. 50 46 # @param other_quantities List of names of other quantities. 51 # @param tagged_elements 52 # @param geo_reference 53 # @param use_inscribed_circle 54 # @param mesh_filename 55 # @param use_cache 56 # @param verbose 57 # @param full_send_dict 58 # @param ghost_recv_dict 59 # @param processor 60 # @param numproc 61 # @param number_of_full_nodes 62 # @param number_of_full_triangles 47 # @param tagged_elements ?? 48 # @param geo_reference ?? 49 # @param use_inscribed_circle ?? 50 # @param mesh_filename ?? 51 # @param use_cache ?? 52 # @param verbose True if this method is to be verbose. 53 # @param full_send_dict ?? 54 # @param ghost_recv_dict ?? 55 # @param processor ?? 56 # @param numproc ?? 57 # @param number_of_full_nodes ?? 58 # @param number_of_full_triangles ?? 63 59 def __init__(self, source=None, 64 60 triangles=None, … … 119 115 number_of_full_triangles=number_of_full_triangles, 120 116 verbose=verbose) 121 117 122 118 # Expose Mesh attributes (FIXME: Maybe turn into methods) 123 119 self.centroid_coordinates = self.mesh.centroid_coordinates 124 self.vertex_coordinates = self.mesh.vertex_coordinates 120 self.vertex_coordinates = self.mesh.vertex_coordinates 125 121 self.boundary = self.mesh.boundary 126 122 self.neighbours = self.mesh.neighbours 127 self.surrogate_neighbours = self.mesh.surrogate_neighbours 123 self.surrogate_neighbours = self.mesh.surrogate_neighbours 128 124 self.neighbour_edges = self.mesh.neighbour_edges 129 125 self.normals = self.mesh.normals 130 self.edgelengths = self.mesh.edgelengths 131 self.radii = self.mesh.radii 132 self.areas = self.mesh.areas 133 134 self.number_of_boundaries = self.mesh.number_of_boundaries 126 self.edgelengths = self.mesh.edgelengths 127 self.radii = self.mesh.radii 128 self.areas = self.mesh.areas 129 130 self.number_of_boundaries = self.mesh.number_of_boundaries 135 131 self.number_of_full_nodes = self.mesh.number_of_full_nodes 136 self.number_of_full_triangles = self.mesh.number_of_full_triangles 132 self.number_of_full_triangles = self.mesh.number_of_full_triangles 137 133 self.number_of_triangles_per_node = self.mesh.number_of_triangles_per_node 138 134 139 135 self.vertex_value_indices = self.mesh.vertex_value_indices 140 self.number_of_triangles = self.mesh.number_of_triangles 136 self.number_of_triangles = self.mesh.number_of_triangles 141 137 142 138 self.geo_reference = self.mesh.geo_reference 143 139 144 140 if verbose: print 'Initialising Domain' 145 141 … … 276 272 277 273 if verbose: print 'Domain: Done' 278 279 280 274 275 ###### 281 276 # Expose underlying Mesh functionality 277 ###### 278 282 279 def __len__(self): 283 280 return len(self.mesh) … … 285 282 def get_centroid_coordinates(self, *args, **kwargs): 286 283 return self.mesh.get_centroid_coordinates(*args, **kwargs) 287 284 288 285 def get_radii(self, *args, **kwargs): 289 return self.mesh.get_radii(*args, **kwargs) 290 286 return self.mesh.get_radii(*args, **kwargs) 287 291 288 def get_areas(self, *args, **kwargs): 292 return self.mesh.get_areas(*args, **kwargs) 289 return self.mesh.get_areas(*args, **kwargs) 293 290 294 291 def get_area(self, *args, **kwargs): … … 296 293 297 294 def get_vertex_coordinates(self, *args, **kwargs): 298 return self.mesh.get_vertex_coordinates(*args, **kwargs) 295 return self.mesh.get_vertex_coordinates(*args, **kwargs) 299 296 300 297 def get_triangles(self, *args, **kwargs): 301 return self.mesh.get_triangles(*args, **kwargs) 302 298 return self.mesh.get_triangles(*args, **kwargs) 299 303 300 def get_nodes(self, *args, **kwargs): 304 301 return self.mesh.get_nodes(*args, **kwargs) 305 302 306 303 def get_number_of_nodes(self, *args, **kwargs): 307 304 return self.mesh.get_number_of_nodes(*args, **kwargs) 308 305 309 306 def get_normal(self, *args, **kwargs): 310 return self.mesh.get_normal(*args, **kwargs) 311 307 return self.mesh.get_normal(*args, **kwargs) 308 312 309 def get_intersecting_segments(self, *args, **kwargs): 313 310 return self.mesh.get_intersecting_segments(*args, **kwargs) 314 311 315 312 def get_disconnected_triangles(self, *args, **kwargs): 316 313 return self.mesh.get_disconnected_triangles(*args, **kwargs) 317 314 318 315 def get_boundary_tags(self, *args, **kwargs): 319 316 return self.mesh.get_boundary_tags(*args, **kwargs) … … 321 318 def get_boundary_polygon(self, *args, **kwargs): 322 319 return self.mesh.get_boundary_polygon(*args, **kwargs) 323 320 324 321 def get_number_of_triangles_per_node(self, *args, **kwargs): 325 322 return self.mesh.get_number_of_triangles_per_node(*args, **kwargs) 326 323 327 324 def get_triangles_and_vertices_per_node(self, *args, **kwargs): 328 325 return self.mesh.get_triangles_and_vertices_per_node(*args, **kwargs) 329 326 330 327 def get_interpolation_object(self, *args, **kwargs): 331 return self.mesh.get_interpolation_object(*args, **kwargs) 332 328 return self.mesh.get_interpolation_object(*args, **kwargs) 329 333 330 def get_tagged_elements(self, *args, **kwargs): 334 return self.mesh.get_tagged_elements(*args, **kwargs) 335 331 return self.mesh.get_tagged_elements(*args, **kwargs) 332 336 333 def get_lone_vertices(self, *args, **kwargs): 337 return self.mesh.get_lone_vertices(*args, **kwargs) 338 334 return self.mesh.get_lone_vertices(*args, **kwargs) 335 339 336 def get_unique_vertices(self, *args, **kwargs): 340 return self.mesh.get_unique_vertices(*args, **kwargs) 337 return self.mesh.get_unique_vertices(*args, **kwargs) 341 338 342 339 def get_georeference(self, *args, **kwargs): 343 340 return self.mesh.get_georeference(*args, **kwargs) 344 341 345 342 def set_georeference(self, *args, **kwargs): 346 self.mesh.set_georeference(*args, **kwargs) 347 343 self.mesh.set_georeference(*args, **kwargs) 344 348 345 def build_tagged_elements_dictionary(self, *args, **kwargs): 349 346 self.mesh.build_tagged_elements_dictionary(*args, **kwargs) 350 347 351 348 def statistics(self, *args, **kwargs): 352 return self.mesh.statistics(*args, **kwargs) 353 354 355 349 return self.mesh.statistics(*args, **kwargs) 350 356 351 ## 357 352 # @brief Get conserved quantities for a volume. … … 365 360 vertex=None, 366 361 edge=None): 367 """Get conserved quantities at volume vol_id 362 """Get conserved quantities at volume vol_id. 368 363 369 364 If vertex is specified use it as index for vertex values … … 397 392 # @param time The new model time (seconds). 398 393 def set_time(self, time=0.0): 399 """Set the model time (seconds)""" 394 """Set the model time (seconds).""" 395 400 396 # FIXME: this is setting the relative time 401 397 # Note that get_time and set_time are now not symmetric … … 407 403 # @return The absolute model time (seconds). 408 404 def get_time(self): 409 """Get the absolute model time (seconds) """405 """Get the absolute model time (seconds).""" 410 406 411 407 return self.time + self.starttime … … 415 411 # @param beta The new beta value. 416 412 def set_beta(self, beta): 417 """Set default beta for limiting """413 """Set default beta for limiting.""" 418 414 419 415 self.beta = beta 420 416 for name in self.quantities: 421 #print 'setting beta for quantity ',name422 417 Q = self.quantities[name] 423 418 Q.set_beta(beta) … … 427 422 # @return The beta value used for limiting. 428 423 def get_beta(self): 429 """Get default beta for limiting """424 """Get default beta for limiting.""" 430 425 431 426 return self.beta … … 436 431 # @note If 'n' is not 1 or 2, raise exception. 437 432 def set_default_order(self, n): 438 """Set default (spatial) order to either 1 or 2 """433 """Set default (spatial) order to either 1 or 2.""" 439 434 440 435 msg = 'Default order must be either 1 or 2. I got %s' % n … … 443 438 self.default_order = n 444 439 self._order_ = self.default_order 445 446 if self.default_order == 1:447 pass448 #self.set_timestepping_method('euler')449 #self.set_all_limiters(beta_euler)450 451 if self.default_order == 2:452 pass453 #self.set_timestepping_method('rk2')454 #self.set_all_limiters(beta_rk2)455 440 456 441 ## … … 681 666 ## 682 667 # @brief Set quantities based on a regional tag. 683 # @param args 684 # @param kwargs 668 # @param args 669 # @param kwargs 685 670 def set_region(self, *args, **kwargs): 686 671 """Set quantities based on a regional tag. … … 1228 1213 self.starttime = float(time) 1229 1214 1230 # --------------------------1215 ################################################################################ 1231 1216 # Main components of evolve 1232 # --------------------------1217 ################################################################################ 1233 1218 1234 1219 ## … … 1275 1260 'evolving system, ' 1276 1261 'e.g. using the method set_boundary.\n' 1277 'This system has the boundary tags %s ' ) \1278 % self.get_boundary_tags() 1262 'This system has the boundary tags %s ' 1263 % self.get_boundary_tags()) 1279 1264 assert hasattr(self, 'boundary_objects'), msg 1280 1265 … … 1287 1272 1288 1273 if finaltime is not None and duration is not None: 1289 # print 'F', finaltime, duration1290 1274 msg = 'Only one of finaltime and duration may be specified' 1291 1275 raise Exception, msg … … 1296 1280 self.finaltime = self.starttime + float(duration) 1297 1281 1298 N = len(self) # Number of triangles1299 self.yieldtime = 0.0 # Track time between 'yields'1282 N = len(self) # Number of triangles 1283 self.yieldtime = 0.0 # Track time between 'yields' 1300 1284 1301 1285 # Initialise interval of timestep sizes (for reporting only) … … 1322 1306 1323 1307 if skip_initial_step is False: 1324 yield(self.time) # Yield initial values1308 yield(self.time) # Yield initial values 1325 1309 1326 1310 while True: … … 1347 1331 if self.time > finaltime: 1348 1332 # FIXME (Ole, 30 April 2006): Do we need this check? 1349 # Probably not (Ole, 18 September 2008). Now changed to1350 # Exception1333 # Probably not (Ole, 18 September 2008). 1334 # Now changed to Exception. 1351 1335 msg = ('WARNING (domain.py): time overshot finaltime. ' 1352 1336 'Contact Ole.Nielsen@ga.gov.au') … … 1417 1401 self.backup_conserved_quantities() 1418 1402 1419 # --------------------------------------1403 ###### 1420 1404 # First euler step 1421 # --------------------------------------1405 ###### 1422 1406 1423 1407 # Compute fluxes across each element edge … … 1442 1426 self.update_boundary() 1443 1427 1444 # ------------------------------------1428 ###### 1445 1429 # Second Euler step 1446 # ------------------------------------1430 ###### 1447 1431 1448 1432 # Compute fluxes across each element edge … … 1452 1436 self.update_conserved_quantities() 1453 1437 1454 # ------------------------------------1438 ###### 1455 1439 # Combine initial and final values 1456 1440 # of conserved quantities and cleanup 1457 # ------------------------------------1441 ###### 1458 1442 1459 1443 # Combine steps … … 1484 1468 initial_time = self.time 1485 1469 1486 # --------------------------------------1470 ###### 1487 1471 # First euler step 1488 # --------------------------------------1472 ###### 1489 1473 1490 1474 # Compute fluxes across each element edge … … 1509 1493 self.update_boundary() 1510 1494 1511 # ------------------------------------1495 ###### 1512 1496 # Second Euler step 1513 # ------------------------------------1497 ###### 1514 1498 1515 1499 # Compute fluxes across each element edge … … 1519 1503 self.update_conserved_quantities() 1520 1504 1521 # ------------------------------------1522 # Combine steps to obtain intermediate1523 # solution at time t^n + 0.5 h1524 # ------------------------------------1505 ###### 1506 # Combine steps to obtain intermediate 1507 # solution at time t^n + 0.5 h 1508 ###### 1525 1509 1526 1510 # Combine steps … … 1539 1523 self.update_boundary() 1540 1524 1541 # ------------------------------------1525 ###### 1542 1526 # Third Euler step 1543 # ------------------------------------1527 ###### 1544 1528 1545 1529 # Compute fluxes across each element edge … … 1549 1533 self.update_conserved_quantities() 1550 1534 1551 # ------------------------------------1535 ###### 1552 1536 # Combine final and initial values 1553 1537 # and cleanup 1554 # ------------------------------------1538 ###### 1555 1539 1556 1540 # Combine steps … … 1631 1615 1632 1616 ## 1633 # @brief 1634 # @param yieldstep 1635 # @param finaltime 1617 # @brief 1618 # @param yieldstep 1619 # @param finaltime 1636 1620 def update_timestep(self, yieldstep, finaltime): 1637 1621 from anuga.config import min_timestep, max_timestep … … 1779 1763 elif self._order_ == 2: 1780 1764 Q.extrapolate_second_order() 1781 #Q.limit()1782 1765 else: 1783 1766 raise Exception, 'Unknown order' 1784 #Q.interpolate_from_vertices_to_edges()1785 1767 1786 1768 ## 1787 1769 # @brief Calculate the norm of the centroid values of a specific quantity, 1788 1770 # using normfunc. 1789 # @param quantity 1790 # @param normfunc 1771 # @param quantity 1772 # @param normfunc 1791 1773 def centroid_norm(self, quantity, normfunc): 1792 """Calculate the norm of the centroid values 1793 of a specific quantity,using normfunc.1774 """Calculate the norm of the centroid values of a specific quantity, 1775 using normfunc. 1794 1776 1795 1777 normfunc should take a list to a float. … … 1801 1783 1802 1784 1803 # ------------------1785 ###### 1804 1786 # Initialise module 1805 # ------------------1787 ###### 1806 1788 1807 1789 # Optimisation with psyco 1808 1790 from anuga.config import use_psyco 1791 1809 1792 if use_psyco: 1810 1793 try: -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r6221 r6226 17 17 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar 18 18 from anuga.utilities.polygon import inside_polygon 19 20 19 from anuga.geospatial_data.geospatial_data import Geospatial_data 21 20 from anuga.fit_interpolate.fit import fit_to_mesh 22 21 from anuga.config import points_file_block_line_size as default_block_line_size 23 22 from anuga.config import epsilon 24 from anuga.caching import cache25 26 27 23 28 24 import Numeric as num 29 25 30 26 27 ## 28 # @brief Implement values at each triangular element. 31 29 class Quantity: 32 30 33 def __init__(self, domain, vertex_values=None): 34 31 ## 32 # @brief Construct values art each triangular element. 33 # @param domain ?? 34 # @param vertex_values ?? 35 def __init__(self, domain, 36 vertex_values=None): 35 37 from anuga.abstract_2d_finite_volumes.domain import Domain 38 36 39 msg = 'First argument in Quantity.__init__ ' 37 40 msg += 'must be of class Domain (or a subclass thereof). ' 38 msg += 'I got %s.' % str(domain.__class__)41 msg += 'I got %s.' % str(domain.__class__) 39 42 assert isinstance(domain, Domain), msg 40 43 41 44 if vertex_values is None: 42 N = len(domain) # number_of_elements45 N = len(domain) # number_of_elements 43 46 self.vertex_values = num.zeros((N, 3), num.Float) 44 47 else: … … 46 49 47 50 N, V = self.vertex_values.shape 48 assert V == 3,\ 49 'Three vertex values per element must be specified' 50 51 52 msg = 'Number of vertex values (%d) must be consistent with'\ 53 %N 54 msg += 'number of elements in specified domain (%d).'\ 55 %len(domain) 56 51 assert V == 3, 'Three vertex values per element must be specified' 52 53 msg = 'Number of vertex values (%d) must be consistent with' % N 54 msg += 'number of elements in specified domain (%d).' % len(domain) 57 55 assert N == len(domain), msg 58 56 … … 68 66 69 67 # Allocate space for Limiter Phi 70 self.phi = num.zeros(N, num.Float) 68 self.phi = num.zeros(N, num.Float) 71 69 72 70 # Intialise centroid and edge_values … … 87 85 self.set_beta(1.0) 88 86 89 90 87 # 91 88 # Methods for operator overloading 89 # 90 91 ## 92 # @brief Get length of object. 93 # @return Return axis length of array. 92 94 def __len__(self): 93 95 return self.centroid_values.shape[0] 94 96 95 97 ## 98 # @brief Negate all values in this quantity. 96 99 def __neg__(self): 97 100 """Negate all values in this quantity giving meaning to the … … 103 106 return Q 104 107 105 108 ## 109 # @brief Add two quantities. 110 # @param other The other quantity (to the right). 106 111 def __add__(self, other): 107 112 """Add to self anything that could populate a quantity … … 119 124 return result 120 125 126 ## 127 # @brief Add two quantities. 128 # @param other The other quantity (to the left). 121 129 def __radd__(self, other): 122 130 """Handle cases like 7+Q, where Q is an instance of class Quantity 123 131 """ 132 124 133 return self + other 125 134 126 135 ## 136 # @brief Subtract two quantities. 137 # @param other The other quantity (to the right). 127 138 def __sub__(self, other): 128 return self + -other #Invoke __neg__ 129 139 return self + -other # Invoke self.__neg__() 140 141 ## 142 # @brief Multiply two quantities. 143 # @param other The other quantity (to the right). 130 144 def __mul__(self, other): 131 145 """Multiply self with anything that could populate a quantity … … 138 152 if isinstance(other, Quantity): 139 153 Q = other 140 else: 154 else: 141 155 Q = Quantity(self.domain) 142 156 Q.set_values(other) … … 151 165 result.edge_values = self.edge_values * Q.edge_values 152 166 result.centroid_values = self.centroid_values * Q.centroid_values 153 167 154 168 return result 155 169 170 ## 171 # @brief Multiply two quantities. 172 # @param other The other quantity (to the left). 156 173 def __rmul__(self, other): 157 174 """Handle cases like 3*Q, where Q is an instance of class Quantity 158 175 """ 176 159 177 return self * other 160 178 179 ## 180 # @brief Divide two quantities. 181 # @param other The other quantity (to the right). 161 182 def __div__(self, other): 162 183 """Divide self with anything that could populate a quantity … … 172 193 if isinstance(other, Quantity): 173 194 Q = other 174 else: 195 else: 175 196 Q = Quantity(self.domain) 176 197 Q.set_values(other) … … 188 209 return result 189 210 211 ## 212 # @brief Divide two quantities. 213 # @param other The other quantity (to the left). 190 214 def __rdiv__(self, other): 191 215 """Handle cases like 3/Q, where Q is an instance of class Quantity 192 216 """ 217 193 218 return self / other 194 219 220 ## 221 # @brief Raise a quaintity to some power. 222 # @param other The power to raise quantity. 195 223 def __pow__(self, other): 196 224 """Raise quantity to (numerical) power … … 201 229 Example using __pow__: 202 230 Q = (Q1**2 + Q2**2)**0.5 203 204 231 """ 205 232 206 233 if isinstance(other, Quantity): 207 234 Q = other 208 else: 235 else: 209 236 Q = Quantity(self.domain) 210 237 Q.set_values(other) … … 222 249 return result 223 250 224 #def __sqrt__(self, other): 225 # """Define in terms of x**0.5 226 # """ 227 # pass 228 229 def set_beta(self,beta): 230 """Set default beta value for limiting 231 """ 251 ## 252 # @brief Set default beta value for limiting. 253 # @param beta ?? 254 def set_beta(self, beta): 255 """Set default beta value for limiting """ 232 256 233 257 if beta < 0.0: … … 235 259 if beta > 2.0: 236 260 print 'WARNING: setting beta > 2.0' 237 261 238 262 self.beta = beta 239 263 264 ## 265 # @brief Get the current beta value. 266 # @return The current beta value. 240 267 def get_beta(self): 241 """Get default beta value for limiting 242 """ 268 """Get default beta value for limiting""" 243 269 244 270 return self.beta 245 271 272 ## 273 # @brief Compute interpolated values at edges and centroid. 274 # @note vertex_values must be set before calling this. 246 275 def interpolate(self): 247 276 """Compute interpolated values at edges and centroid 248 277 Pre-condition: vertex_values have been set 249 278 """ 250 279 251 280 # FIXME (Ole): Maybe this function 252 281 # should move to the C-interface? 253 282 # However, it isn't called by validate_all.py, so it 254 283 # may not be that important to optimise it? 255 284 256 285 N = self.vertex_values.shape[0] 257 286 for i in range(N): … … 264 293 self.interpolate_from_vertices_to_edges() 265 294 266 295 ## 296 # @brief ?? 267 297 def interpolate_from_vertices_to_edges(self): 268 # Call correct module function 269 # (either from this module or C-extension) 298 # Call correct module function (either from this module or C-extension) 270 299 interpolate_from_vertices_to_edges(self) 271 300 301 ## 302 # @brief ?? 272 303 def interpolate_from_edges_to_vertices(self): 273 # Call correct module function 274 # (either from this module or C-extension) 304 # Call correct module function (either from this module or C-extension) 275 305 interpolate_from_edges_to_vertices(self) 276 277 278 279 306 280 307 #--------------------------------------------- 281 308 # Public interface for setting quantity values 282 309 #--------------------------------------------- 283 def set_values(self, 284 numeric=None, # List, numeric array or constant 285 quantity=None, # Another quantity 286 function=None, # Callable object: f(x,y) 287 geospatial_data=None, # Arbitrary dataset 288 filename=None, attribute_name=None, # Input from file 289 alpha=None, 290 location='vertices', 291 polygon=None, 292 indices=None, 293 smooth=False, 294 verbose=False, 295 use_cache=False): 296 310 311 ## 312 # @brief Set values for quantity based on different sources. 313 # @param numeric A num array, list or constant value. 314 # @param quantity Another Quantity. 315 # @param function Any callable object that takes two 1d arrays. 316 # @param geospatial_data Arbitrary instance of class Geospatial_data 317 # @param filename Path to a points file. 318 # @param attribute_name If specified any array using that name will be used. 319 # @param alpha Smoothing parameter to be used with fit_interpolate.fit. 320 # @param location Where to store values (vertices, edges, centroids). 321 # @param polygon Restrict update to locations that fall inside polygon. 322 # @param indices Restrict update to locations specified by this. 323 # @param smooth If True, smooth vertex values. 324 # @param verbose True if this method is to be verbose. 325 # @param use_cache If True cache results for fit_interpolate.fit. 326 # @note Exactly one of 'numeric', 'quantity', 'function', 'filename' 327 # must be present. 328 def set_values(self, numeric=None, # List, numeric array or constant 329 quantity=None, # Another quantity 330 function=None, # Callable object: f(x,y) 331 geospatial_data=None, # Arbitrary dataset 332 filename=None, 333 attribute_name=None, # Input from file 334 alpha=None, 335 location='vertices', 336 polygon=None, 337 indices=None, 338 smooth=False, 339 verbose=False, 340 use_cache=False): 297 341 """Set values for quantity based on different sources. 298 342 … … 366 410 constant values so far. 367 411 368 indices: Restrict update of quantity to locations that are 412 indices: Restrict update of quantity to locations that are 369 413 identified by indices (e.g. node ids if location 370 414 is 'unique vertices' or triangle ids otherwise). 371 415 372 416 verbose: True means that output to stdout is generated 373 417 … … 391 435 # perhaps the notion of location and indices simplified 392 436 393 # FIXME (Ole): Need to compute indices based on polygon 437 # FIXME (Ole): Need to compute indices based on polygon 394 438 # (and location) and use existing code after that. 395 439 396 440 # See ticket:275, ticket:250, ticeket:254 for refactoring plan 397 441 398 442 if polygon is not None: 399 443 if indices is not None: … … 405 449 msg += 'a constant.' 406 450 if numeric is None: 407 raise Exception, msg 451 raise Exception, msg 408 452 else: 409 453 # Check that numeric is as constant 410 454 assert type(numeric) in [FloatType, IntType, LongType], msg 411 455 412 413 456 location = 'centroids' 414 415 457 416 458 points = self.domain.get_centroid_coordinates(absolute=True) 417 459 indices = inside_polygon(points, polygon) 418 419 self.set_values_from_constant(numeric, 420 location, indices, verbose) 421 460 461 self.set_values_from_constant(numeric, location, indices, verbose) 422 462 423 463 self.extrapolate_first_order() 424 464 425 465 if smooth: 426 self.smooth_vertex_values(use_cache=use_cache, 427 verbose=verbose) 428 429 466 self.smooth_vertex_values() 467 430 468 return 431 432 433 434 435 436 469 437 470 # General input checks 438 471 L = [numeric, quantity, function, geospatial_data, filename] 439 msg = 'Exactly one of the arguments '+\ 440 'numeric, quantity, function, geospatial_data, '+\ 441 'or filename must be present.' 472 msg = ('Exactly one of the arguments numeric, quantity, function, ' 473 'geospatial_data, or filename must be present.') 442 474 assert L.count(None) == len(L)-1, msg 443 444 475 445 476 if location == 'edges': 446 477 msg = 'edges has been deprecated as valid location' 447 478 raise Exception, msg 448 479 449 480 if location not in ['vertices', 'centroids', 'unique vertices']: 450 msg = 'Invalid location: %s' % location481 msg = 'Invalid location: %s' % location 451 482 raise Exception, msg 452 453 483 454 484 msg = 'Indices must be a list or None' 455 485 assert type(indices) in [ListType, NoneType, num.ArrayType], msg 456 486 457 458 459 487 # Determine which 'set_values_from_...' to use 460 461 488 if numeric is not None: 462 489 if type(numeric) in [FloatType, IntType, LongType]: 463 self.set_values_from_constant(numeric, 464 location,indices, verbose)490 self.set_values_from_constant(numeric, location, 491 indices, verbose) 465 492 elif type(numeric) in [num.ArrayType, ListType]: 466 self.set_values_from_array(numeric, 467 location, indices, 468 use_cache=use_cache, 469 verbose=verbose) 493 self.set_values_from_array(numeric, location, indices, verbose) 470 494 elif callable(numeric): 471 self.set_values_from_function(numeric, 472 location, indices, 473 use_cache=use_cache, 474 verbose=verbose) 495 self.set_values_from_function(numeric, location, 496 indices, verbose) 475 497 elif isinstance(numeric, Quantity): 476 self.set_values_from_quantity(numeric, 477 location, indices, 478 verbose=verbose) 498 self.set_values_from_quantity(numeric, location, 499 indices, verbose) 479 500 elif isinstance(numeric, Geospatial_data): 480 self.set_values_from_geospatial_data(numeric, 481 alpha, 501 self.set_values_from_geospatial_data(numeric, alpha, location, 502 indices, verbose=verbose, 503 use_cache=use_cache) 504 else: 505 msg = 'Illegal type for argument numeric: %s' % str(numeric) 506 raise msg 507 elif quantity is not None: 508 self.set_values_from_quantity(quantity, location, indices, verbose) 509 elif function is not None: 510 msg = 'Argument function must be callable' 511 assert callable(function), msg 512 self.set_values_from_function(function, location, indices, verbose) 513 elif geospatial_data is not None: 514 self.set_values_from_geospatial_data(geospatial_data, alpha, 482 515 location, indices, 483 516 verbose=verbose, 484 517 use_cache=use_cache) 485 else:486 msg = 'Illegal type for argument numeric: %s' %str(numeric)487 raise msg488 489 elif quantity is not None:490 self.set_values_from_quantity(quantity,491 location, indices, verbose)492 elif function is not None:493 msg = 'Argument function must be callable'494 assert callable(function), msg495 self.set_values_from_function(function,496 location,497 indices,498 use_cache=use_cache,499 verbose=verbose)500 elif geospatial_data is not None:501 self.set_values_from_geospatial_data(geospatial_data,502 alpha,503 location, indices,504 verbose=verbose,505 use_cache=use_cache)506 507 518 elif filename is not None: 508 519 if hasattr(self.domain, 'points_file_block_line_size'): … … 510 521 else: 511 522 max_read_lines = default_block_line_size 512 self.set_values_from_file(filename, attribute_name, alpha, 513 location, indices, 514 verbose=verbose, 523 self.set_values_from_file(filename, attribute_name, alpha, location, 524 indices, verbose=verbose, 515 525 max_read_lines=max_read_lines, 516 526 use_cache=use_cache) 517 527 else: 518 raise Exception, 'This can\'t happen :-)' 519 520 528 raise Exception, "This can't happen :-)" 521 529 522 530 # Update all locations in triangles … … 531 539 532 540 533 #------------------------------------------------------------- 541 #--------------------------------------------------------------------------- 534 542 # Specific internal functions for setting values based on type 535 #------------------------------------------------------------- 536 537 def set_values_from_constant(self, X, 538 location, indices, verbose): 539 """Set quantity values from specified constant X 540 """ 543 #--------------------------------------------------------------------------- 544 545 ## 546 # @brief Set quantity values from specified constant. 547 # @param X The constant to set quantity values to. 548 # @param location 549 # @param indices 550 # @param verbose 551 def set_values_from_constant(self, X, location, indices, verbose): 552 """Set quantity values from specified constant X""" 541 553 542 554 # FIXME (Ole): Somehow indices refer to centroids 543 555 # rather than vertices as default. See unit test 544 556 # test_set_vertex_values_using_general_interface_with_subset(self): 545 546 557 547 558 if location == 'centroids': … … 552 563 for i in indices: 553 564 self.centroid_values[i] = X 554 555 #elif location == 'edges':556 # if indices is None:557 # self.edge_values[:] = X558 # else:559 # # Brute force560 # for i in indices:561 # self.edge_values[i] = X562 563 565 elif location == 'unique vertices': 564 566 if indices is None: 565 567 self.edge_values[:] = X #FIXME (Ole): Shouldn't this be vertex_values? 566 568 else: 567 568 569 # Go through list of unique vertices 569 570 for unique_vert_id in indices: 570 571 triangles =self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)572 571 triangles = \ 572 self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id) 573 573 574 # In case there are unused points 574 575 if len(triangles) == 0: 575 576 continue 576 577 577 578 # Go through all triangle, vertex pairs 578 579 # and set corresponding vertex value … … 590 591 self.vertex_values[i_vertex] = X 591 592 592 593 594 593 ## 594 # @brief Set values for a quantity. 595 # @param values Array of values. 596 # @param location Where values are to be stored. 597 # @param indices Limit update to these indices. 598 # @param verbose True if this method is to be verbose. 595 599 def set_values_from_array(self, values, 596 location='vertices', 597 indices=None, 598 use_cache=False, 599 verbose=False): 600 location='vertices', 601 indices=None, 602 verbose=False): 600 603 """Set values for quantity 601 604 … … 628 631 if indices is not None: 629 632 indices = num.array(indices, num.Int) 630 msg = 'Number of values must match number of indices:'631 msg += ' You specified %d values and %d indices'\632 % (values.shape[0], indices.shape[0])633 msg = ('Number of values must match number of indices: You ' 634 'specified %d values and %d indices' 635 % (values.shape[0], indices.shape[0])) 633 636 assert values.shape[0] == indices.shape[0], msg 634 637 … … 650 653 for i in range(len(indices)): 651 654 self.centroid_values[indices[i]] = values[i] 652 653 655 elif location == 'unique vertices': 654 assert len(values.shape) == 1 or num.allclose(values.shape[1:], 1),\ 655 'Values array must be 1d' 656 657 self.set_vertex_values(values.flat, 658 indices=indices, 659 use_cache=use_cache, 660 verbose=verbose) 661 656 assert (len(values.shape) == 1 or num.allclose(values.shape[1:], 1), 657 'Values array must be 1d') 658 659 self.set_vertex_values(values.flat, indices=indices) 662 660 else: 663 661 # Location vertices 664 662 if len(values.shape) == 1: 665 # This is the common case arising from fitted 666 # values (e.g. from pts file). 667 self.set_vertex_values(values, 668 indices=indices, 669 use_cache=use_cache, 670 verbose=verbose) 671 663 self.set_vertex_values(values, indices=indices) 672 664 elif len(values.shape) == 2: 673 665 # Vertex values are given as a triplet for each triangle 674 675 666 msg = 'Array must be N x 3' 676 667 assert values.shape[1] == 3, msg … … 683 674 else: 684 675 msg = 'Values array must be 1d or 2d' 685 raise msg 686 687 688 def set_values_from_quantity(self, q, 689 location, indices, verbose): 676 raise Exception, msg 677 678 ## 679 # @brief Set quantity values from a specified quantity instance. 680 # @param q The quantity instance to take values from. 681 # @param location IGNORED, 'vertices' ALWAYS USED! 682 # @param indices ?? 683 # @param verbose True if this method is to be verbose. 684 def set_values_from_quantity(self, q, location, indices, verbose): 690 685 """Set quantity values from specified quantity instance q 691 686 … … 702 697 703 698 self.set_values(A, location='vertices', 704 indices=indices, 705 verbose=verbose) 706 707 699 indices=indices, verbose=verbose) 700 701 ## 702 # @brief Set quantity values from a specified quantity instance. 703 # @param f Callable that takes two 1d array -> 1d array. 704 # @param location Where values are to be stored. 705 # @param indices ?? 706 # @param verbose True if this method is to be verbose. 708 707 def set_values_from_function(self, f, 709 location='vertices', 710 indices=None, 711 use_cache=False, 712 verbose=False): 708 location='vertices', 709 indices=None, 710 verbose=False): 713 711 """Set values for quantity using specified function 714 712 715 713 Input 716 717 714 f: x, y -> z Function where x, y and z are arrays 718 715 location: Where values are to be stored. … … 720 717 unique vertices 721 718 Default is "vertices" 722 indices: 723 724 719 indices: 725 720 """ 726 721 727 722 # FIXME: Should check that function returns something sensible and 728 # raise a meaningful lexception if it returns None for example723 # raise a meaningful exception if it returns None for example 729 724 730 725 # FIXME: Should supply absolute coordinates 731 732 726 733 727 # Compute the function values and call set_values again … … 735 729 if indices is None: 736 730 indices = range(len(self)) 737 731 738 732 V = num.take(self.domain.get_centroid_coordinates(), indices) 739 740 x = V[:,0]; y = V[:,1] 741 if use_cache is True: 742 res = cache(f, (x, y), 743 verbose=verbose) 744 else: 745 res = f(x, y) 746 747 self.set_values(res, 748 location=location, 749 indices=indices) 750 733 self.set_values(f(V[:,0], V[:,1]), 734 location=location, indices=indices) 751 735 elif location == 'vertices': 752 # This is the default branch taken by set_quantity753 754 736 M = self.domain.number_of_triangles 755 737 V = self.domain.get_vertex_coordinates() 756 738 757 x = V[:,0]; y = V[:,1] 758 if use_cache is True: 759 #print 'Caching function' 760 values = cache(f, (x, y), 761 verbose=verbose) 762 else: 763 if verbose is True: 764 print 'Evaluating function in set_values' 765 values = f(x, y) 766 767 #print 'value', min(x), max(x), min(y), max(y), max(values), len(values) 768 739 x = V[:,0]; 740 y = V[:,1]; 741 values = f(x, y) 769 742 770 743 # FIXME (Ole): This code should replace all the … … 773 746 # If that could be resolved this one will be 774 747 # more robust and simple. 775 776 #values = num.reshape(values, (M,3))777 #self.set_values(values,778 # location='vertices',779 # indices=indices)780 781 748 782 749 # This should be removed 783 750 if is_scalar(values): 784 751 # Function returned a constant value 785 self.set_values_from_constant(values, 786 location,indices, verbose)752 self.set_values_from_constant(values, location, 753 indices, verbose) 787 754 return 788 755 789 # This should be removed 756 # This should be removed 790 757 if indices is None: 791 758 for j in range(3): 792 self.vertex_values[:, j] = values[j::3]793 else: 759 self.vertex_values[:, j] = values[j::3] 760 else: 794 761 # Brute force 795 762 for i in indices: 796 763 for j in range(3): 797 self.vertex_values[i,j] = values[3*i+j] 798 799 764 self.vertex_values[i, j] = values[3*i + j] 800 765 else: 801 raise 'Not implemented: %s' %location 802 803 804 805 def set_values_from_geospatial_data(self, geospatial_data, alpha, 806 location, indices, 807 verbose=False, 808 use_cache=False): 809 """ Set values based on geo referenced geospatial data object. 810 """ 766 raise Exception, 'Not implemented: %s' % location 767 768 ## 769 # @brief Set values based on geo referenced geospatial data object. 770 # @param geospatial_data ?? 771 # @param alpha ?? 772 # @param location ?? 773 # @param indices ?? 774 # @param verbose ?? 775 # @param use_cache ?? 776 def set_values_from_geospatial_data(self, geospatial_data, 777 alpha, 778 location, 779 indices, 780 verbose=False, 781 use_cache=False): 782 """Set values based on geo referenced geospatial data object.""" 783 784 from anuga.coordinate_transforms.geo_reference import Geo_reference 811 785 812 786 points = geospatial_data.get_data_points(absolute=False) … … 814 788 data_georef = geospatial_data.get_geo_reference() 815 789 816 817 from anuga.coordinate_transforms.geo_reference import Geo_reference818 819 820 790 points = ensure_numeric(points, num.Float) 821 791 values = ensure_numeric(values, num.Float) 822 792 823 793 if location != 'vertices': 824 msg = 'set_values_from_points is only defined for '+\ 825 'location=\'vertices\'' 826 raise ms 827 828 #coordinates = self.domain.get_nodes() 829 #triangles = self.domain.get_triangles() 830 794 msg = ("set_values_from_points is only defined for " 795 "location='vertices'") 796 raise Exception, msg 831 797 832 798 # Take care of georeferencing … … 834 800 data_georef = Geo_reference() 835 801 836 837 802 mesh_georef = self.domain.geo_reference 838 803 839 840 804 # Call fit_interpolate.fit function 841 # args = (coordinates, triangles, points, values)842 805 args = (points, ) 843 806 kwargs = {'vertex_coordinates': None, … … 850 813 'verbose': verbose} 851 814 852 vertex_attributes = apply(fit_to_mesh, 853 args, kwargs) 815 vertex_attributes = apply(fit_to_mesh, args, kwargs) 854 816 855 817 # Call underlying method using array values 856 self.set_values_from_array(vertex_attributes, 857 location, indices, 858 use_cache=use_cache, 859 verbose=verbose) 860 861 862 863 def set_values_from_points(self, points, values, alpha, 864 location, indices, 865 data_georef=None, 866 verbose=False, 867 use_cache=False): 868 """ 869 Set quantity values from arbitray data points using 870 fit_interpolate.fit 871 """ 818 self.set_values_from_array(vertex_attributes, location, 819 indices, verbose) 820 821 ## 822 # @brief Set quantity values from arbitray data points. 823 # @param points ?? 824 # @param values ?? 825 # @param alpha ?? 826 # @param location ?? 827 # @param indices ?? 828 # @param data_georef ?? 829 # @param verbose True if this method is to be verbose. 830 # @param use_cache ?? 831 def set_values_from_points(self, points, 832 values, 833 alpha, 834 location, 835 indices, 836 data_georef=None, 837 verbose=False, 838 use_cache=False): 839 """Set quantity values from arbitray data points using fit_interpolate.fit""" 872 840 873 841 raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead' 874 875 876 def set_values_from_file(self, filename, attribute_name, alpha, 877 location, indices, 878 verbose=False, 879 use_cache=False, 880 max_read_lines=None): 881 """Set quantity based on arbitrary points in a points file 882 using attribute_name selects name of attribute 883 present in file. 842 843 ## 844 # @brief Set quantity based on arbitrary points in a points file. 845 # @param filename Path to the points file. 846 # @param attribute_name 847 # @param alpha 848 # @param location 849 # @param indices 850 # @param verbose True if this method is to be verbose. 851 # @param use_cache 852 # @param max_read_lines 853 def set_values_from_file(self, filename, 854 attribute_name, 855 alpha, 856 location, 857 indices, 858 verbose=False, 859 use_cache=False, 860 max_read_lines=None): 861 """Set quantity based on arbitrary points in a points file using 862 attribute_name selects name of attribute present in file. 884 863 If attribute_name is not specified, use first available attribute 885 864 as defined in geospatial_data. 886 865 """ 887 866 888 867 from types import StringType 868 889 869 msg = 'Filename must be a text string' 890 870 assert type(filename) == StringType, msg 891 871 892 893 872 if location != 'vertices': 894 msg = 'set_values_from_file is only defined for '+\ 895 'location=\'vertices\'' 896 raise msg 897 898 if True: 873 msg = "set_values_from_file is only defined for location='vertices'" 874 raise Exception, msg 875 876 if True: 899 877 # Use mesh as defined by domain 900 # This used to cause problems for caching 901 # due to quantities changing, but 902 # it now works using the appropriate Mesh object. 878 # This used to cause problems for caching due to quantities 879 # changing, but it now works using the appropriate Mesh object. 903 880 # This should close ticket:242 904 881 vertex_attributes = fit_to_mesh(filename, 905 mesh=self.domain.mesh, 882 mesh=self.domain.mesh, 906 883 alpha=alpha, 907 884 attribute_name=attribute_name, … … 911 888 else: 912 889 # This variant will cause Mesh object to be recreated 913 # in fit_to_mesh thus doubling up on the neighbour structure 890 # in fit_to_mesh thus doubling up on the neighbour structure 914 891 # FIXME(Ole): This is now obsolete 19 Jan 2009. 915 892 nodes = self.domain.get_nodes(absolute=True) 916 triangles = self.domain.get_triangles() 893 triangles = self.domain.get_triangles() 917 894 vertex_attributes = fit_to_mesh(filename, 918 nodes, triangles, 895 nodes, triangles, 919 896 mesh=None, 920 897 alpha=alpha, … … 923 900 verbose=verbose, 924 901 max_read_lines=max_read_lines) 925 902 926 903 # Call underlying method using array values 927 if verbose: 928 print 'Applying fitted data to domain' 929 self.set_values_from_array(vertex_attributes, 930 location, indices, 931 use_cache=use_cache, 932 verbose=verbose) 933 934 935 936 #----------------------------------------------------- 904 self.set_values_from_array(vertex_attributes, location, 905 indices, verbose) 906 907 ## 908 # @brief Get index for maximum or minimum value of quantity. 909 # @param mode Either 'max' or 'min'. 910 # @param indices Set of IDs of elements to work on. 937 911 def get_extremum_index(self, mode=None, indices=None): 938 912 """Return index for maximum or minimum value of quantity (on centroids) … … 958 932 if mode is None or mode == 'max': 959 933 i = num.argmax(V) 960 elif mode == 'min': 934 elif mode == 'min': 961 935 i = num.argmin(V) 962 963 936 else: 937 raise ValueError, 'Bad mode value, got: %s' % str(mode) 938 964 939 if indices is None: 965 940 return i … … 967 942 return indices[i] 968 943 969 944 ## 945 # @brief Get index for maximum value of quantity. 946 # @param indices Set of IDs of elements to work on. 970 947 def get_maximum_index(self, indices=None): 971 """See get extreme index for details 972 """ 973 974 return self.get_extremum_index(mode='max', 975 indices=indices) 976 977 978 948 """See get extreme index for details""" 949 950 return self.get_extremum_index(mode='max', indices=indices) 951 952 ## 953 # @brief Return maximum value of quantity (on centroids). 954 # @param indices Set of IDs of elements to work on. 979 955 def get_maximum_value(self, indices=None): 980 956 """Return maximum value of quantity (on centroids) … … 987 963 988 964 Note, we do not seek the maximum at vertices as each vertex can 989 have multiple values - one for each triangle sharing it 990 """ 991 965 have multiple values - one for each triangle sharing it 966 """ 992 967 993 968 i = self.get_maximum_index(indices) 994 V = self.get_values(location='centroids') #, indices=indices)995 969 V = self.get_values(location='centroids') #, indices=indices) 970 996 971 return V[i] 997 998 972 973 ## 974 # @brief Get location of maximum value of quantity (on centroids). 975 # @param indices Set of IDs of elements to work on. 999 976 def get_maximum_location(self, indices=None): 1000 977 """Return location of maximum value of quantity (on centroids) … … 1005 982 Usage: 1006 983 x, y = get_maximum_location() 1007 1008 984 1009 985 Notes: … … 1012 988 1013 989 If there are multiple cells with same maximum value, the 1014 first cell encountered in the triangle array is returned. 990 first cell encountered in the triangle array is returned. 1015 991 """ 1016 992 … … 1020 996 return x, y 1021 997 1022 998 ## 999 # @brief Get index for minimum value of quantity. 1000 # @param indices Set of IDs of elements to work on. 1023 1001 def get_minimum_index(self, indices=None): 1024 """See get extreme index for details 1025 """ 1026 1027 return self.get_extremum_index(mode='min', 1028 indices=indices)1029 1030 1002 """See get extreme index for details""" 1003 1004 return self.get_extremum_index(mode='min', indices=indices) 1005 1006 ## 1007 # @brief Return minimum value of quantity (on centroids). 1008 # @param indices Set of IDs of elements to work on. 1031 1009 def get_minimum_value(self, indices=None): 1032 1010 """Return minimum value of quantity (on centroids) … … 1038 1016 v = get_minimum_value() 1039 1017 1040 See get_maximum_value for more details. 1041 """ 1042 1018 See get_maximum_value for more details. 1019 """ 1043 1020 1044 1021 i = self.get_minimum_index(indices) 1045 1022 V = self.get_values(location='centroids') 1046 1023 1047 1024 return V[i] 1048 1049 1025 1026 1027 ## 1028 # @brief Get location of minimum value of quantity (on centroids). 1029 # @param indices Set of IDs of elements to work on. 1050 1030 def get_minimum_location(self, indices=None): 1051 1031 """Return location of minimum value of quantity (on centroids) … … 1056 1036 Usage: 1057 1037 x, y = get_minimum_location() 1058 1059 1038 1060 1039 Notes: … … 1063 1042 1064 1043 If there are multiple cells with same maximum value, the 1065 first cell encountered in the triangle array is returned. 1044 first cell encountered in the triangle array is returned. 1066 1045 """ 1067 1046 … … 1071 1050 return x, y 1072 1051 1073 1074 1052 ## 1053 # @brief Get values at interpolation points. 1054 # @param interpolation_points List of UTM coords or geospatial data object. 1055 # @param use_cache ?? 1056 # @param verbose True if this method is to be verbose. 1075 1057 def get_interpolated_values(self, interpolation_points, 1076 use_cache=False,1077 verbose=False):1078 """ 1079 1080 The argument interpolation points must be given as either a 1058 use_cache=False, 1059 verbose=False): 1060 """Get values at interpolation points 1061 1062 The argument interpolation points must be given as either a 1081 1063 list of absolute UTM coordinates or a geospatial data object. 1082 1064 """ 1083 1084 1065 1085 1066 # FIXME (Ole): Points might be converted to coordinates relative to mesh origin 1086 # This could all be refactored using the 1087 # 'change_points_geo_ref' method of Class geo_reference. 1067 # This could all be refactored using the 1068 # 'change_points_geo_ref' method of Class geo_reference. 1088 1069 # The purpose is to make interpolation points relative 1089 1070 # to the mesh origin. 1090 1071 # 1091 1072 # Speed is also a consideration here. 1092 1093 1094 # Ensure that interpolation points is either a list of 1073 1074 # Ensure that interpolation points is either a list of 1095 1075 # points, Nx2 array, or geospatial and convert to Numeric array 1096 if isinstance(interpolation_points, Geospatial_data): 1076 if isinstance(interpolation_points, Geospatial_data): 1097 1077 # Ensure interpolation points are in absolute UTM coordinates 1098 interpolation_points = interpolation_points.get_data_points(absolute=True) 1099 1078 interpolation_points = \ 1079 interpolation_points.get_data_points(absolute=True) 1080 1100 1081 # Reconcile interpolation points with georeference of domain 1101 interpolation_points = self.domain.geo_reference.get_relative(interpolation_points) 1082 interpolation_points = \ 1083 self.domain.geo_reference.get_relative(interpolation_points) 1102 1084 interpolation_points = ensure_numeric(interpolation_points) 1103 1085 1104 1086 1105 1087 # Get internal representation (disconnected) of vertex values 1106 1088 vertex_values, triangles = self.get_vertex_values(xy=False, 1107 smooth=False) 1108 1089 smooth=False) 1090 1109 1091 # Get possibly precomputed interpolation object 1110 1092 I = self.domain.get_interpolation_object() 1111 1093 1112 # Call interpolate method with interpolation points 1094 # Call interpolate method with interpolation points 1113 1095 result = I.interpolate_block(vertex_values, interpolation_points, 1114 use_cache=use_cache, 1115 verbose=verbose) 1116 1096 use_cache=use_cache, verbose=verbose) 1097 1117 1098 return result 1118 1119 1120 1121 1122 def get_values(self, 1123 interpolation_points=None, 1124 location='vertices', 1125 indices=None, 1126 use_cache=False, 1127 verbose=False): 1128 """get values for quantity 1129 1130 Extract values for quantity as a Numeric array. 1131 1099 1100 ## 1101 # @brief Get values as an array. 1102 # @param interpolation_points List of coords to get values at. 1103 # @param location Where to store results. 1104 # @param indices Set of IDs of elements to work on. 1105 # @param use_cache 1106 # @param verbose True if this method is to be verbose. 1107 def get_values(self, interpolation_points=None, 1108 location='vertices', 1109 indices=None, 1110 use_cache=False, 1111 verbose=False): 1112 """Get values for quantity 1113 1114 Extract values for quantity as a Numeric array. 1115 1132 1116 Inputs: 1133 1117 interpolation_points: List of x, y coordinates where value is 1134 sought (using interpolation). If points 1135 are given, values of location and indices 1118 sought (using interpolation). If points 1119 are given, values of location and indices 1136 1120 are ignored. Assume either absolute UTM 1137 1121 coordinates or geospatial data object. 1138 1122 1139 1123 location: Where values are to be stored. 1140 1124 Permissible options are: vertices, edges, centroids … … 1148 1132 values will be a list or a Numerical array of length N, N being 1149 1133 the number of elements. 1150 1134 1151 1135 In case of location == 'vertices' or 'edges' the dimension of 1152 1136 returned values will be of dimension Nx3 … … 1154 1138 In case of location == 'unique vertices' the average value at 1155 1139 each vertex will be returned and the dimension of returned values 1156 will be a 1d array of length "number of vertices" 1157 1140 will be a 1d array of length "number of vertices" 1141 1158 1142 Indices is the set of element ids that the operation applies to. 1159 1143 … … 1161 1145 internal ordering. 1162 1146 """ 1163 1164 1147 1165 1148 # FIXME (Ole): I reckon we should have the option of passing a 1166 1149 # polygon into get_values. The question becomes how 1167 1150 # resulting values should be ordered. 1168 1151 1169 1152 if verbose is True: 1170 print 'Getting values from %s' % location1153 print 'Getting values from %s' % location 1171 1154 1172 1155 if interpolation_points is not None: … … 1174 1157 use_cache=use_cache, 1175 1158 verbose=verbose) 1176 1177 1159 1178 1160 # FIXME (Ole): Consider deprecating 'edges' - but not if it is used 1179 # elsewhere in ANUGA. 1161 # elsewhere in ANUGA. 1180 1162 # Edges have already been deprecated in set_values, see changeset:5521, 1181 1163 # but *might* be useful in get_values. Any thoughts anyone? 1182 1183 if location not in ['vertices', 'centroids', 'edges',1184 ' unique vertices']:1185 msg = 'Invalid location: %s' % location1164 1165 if location not in ['vertices', 'centroids', 1166 'edges', 'unique vertices']: 1167 msg = 'Invalid location: %s' % location 1186 1168 raise msg 1187 1169 1188 1170 import types 1171 1189 1172 assert type(indices) in [types.ListType, types.NoneType, 1190 num.ArrayType], \1191 1173 num.ArrayType], \ 1174 'Indices must be a list or None' 1192 1175 1193 1176 if location == 'centroids': 1194 1177 if (indices == None): 1195 1178 indices = range(len(self)) 1196 return num.take(self.centroid_values, indices)1179 return num.take(self.centroid_values, indices) 1197 1180 elif location == 'edges': 1198 1181 if (indices == None): 1199 1182 indices = range(len(self)) 1200 return num.take(self.edge_values, indices)1183 return num.take(self.edge_values, indices) 1201 1184 elif location == 'unique vertices': 1202 1185 if (indices == None): … … 1207 1190 for unique_vert_id in indices: 1208 1191 triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id) 1209 1192 1210 1193 # In case there are unused points 1211 1194 if len(triangles) == 0: 1212 1195 msg = 'Unique vertex not associated with triangles' 1213 raise msg1196 raise Exception, msg 1214 1197 1215 1198 # Go through all triangle, vertex pairs 1216 1199 # Average the values 1217 1218 1200 # FIXME (Ole): Should we merge this with get_vertex_values 1219 1201 sum = 0 1220 1202 for triangle_id, vertex_id in triangles: 1221 1203 sum += self.vertex_values[triangle_id, vertex_id] 1222 vert_values.append(sum /len(triangles))1204 vert_values.append(sum / len(triangles)) 1223 1205 return num.array(vert_values, num.Float) 1224 1206 else: … … 1227 1209 return num.take(self.vertex_values, indices) 1228 1210 1229 1230 1231 def set_vertex_values(self, A, 1232 indices=None, 1233 use_cache=False, 1234 verbose=False): 1211 ## 1212 # @brief Set vertex values for all unique vertices based on array. 1213 # @param A Array to set values with. 1214 # @param indices Set of IDs of elements to work on. 1215 def set_vertex_values(self, A, indices=None): 1235 1216 """Set vertex values for all unique vertices based on input array A 1236 which has one entry per unique vertex, i.e. 1237 one value for each row inarray self.domain.nodes.1217 which has one entry per unique vertex, i.e. one value for each row in 1218 array self.domain.nodes. 1238 1219 1239 1220 indices is the list of vertex_id's that will be set. … … 1242 1223 """ 1243 1224 1244 1245 # Assert that A can be converted to a Numeric array of appropriate dim 1225 # Check that A can be converted to array and is of appropriate dim 1246 1226 A = ensure_numeric(A, num.Float) 1247 1248 # print 'SHAPE A', A.shape1249 1227 assert len(A.shape) == 1 1250 1228 … … 1256 1234 vertex_list = indices 1257 1235 1258 1259 #FIXME(Ole): This function ought to be faster. 1260 # We need to get the triangles_and_vertices list 1261 # from domain in one hit, then cache the computation of the 1262 # Nx3 array of vertex values that can then be assigned using 1263 # set_values_from_array. 1264 # 1265 # Alternatively, some C code would be handy 1266 # 1267 self._set_vertex_values(vertex_list, A) 1268 1269 1270 def _set_vertex_values(self, vertex_list, A): 1271 """Go through list of unique vertices 1272 This is the common case e.g. when values 1273 are obtained from a pts file through fitting 1274 """ 1275 1236 # Go through list of unique vertices 1276 1237 for i_index, unique_vert_id in enumerate(vertex_list): 1277 1278 1238 triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id) 1279 1239 1280 1240 # In case there are unused points 1281 if len(triangles) == 0: continue 1241 if len(triangles) == 0: 1242 continue 1282 1243 1283 1244 # Go through all triangle, vertex pairs … … 1289 1250 self.interpolate() 1290 1251 1291 1292 def smooth_vertex_values(self, 1293 use_cache=False, 1294 verbose=False): 1295 """ Smooths vertex values. 1296 """ 1297 1298 A,V = self.get_vertex_values(xy=False, smooth=True) 1299 self.set_vertex_values(A, 1300 use_cache=use_cache, 1301 verbose=verbose) 1302 1303 1304 1252 ## 1253 # @brief Smooth vertex values. 1254 def smooth_vertex_values(self): 1255 """Smooths vertex values.""" 1256 1257 A, V = self.get_vertex_values(xy=False, smooth=True) 1258 self.set_vertex_values(A) 1259 1260 ############################################################################ 1305 1261 # Methods for outputting model results 1306 def get_vertex_values(self, 1307 xy=True, 1308 smooth=None, 1309 precision=None): 1262 ############################################################################ 1263 1264 ## 1265 # @brief Get vertex values like an OBJ format i.e. one value per node. 1266 # @param xy True if we return X and Y as well as A and V. 1267 # @param smooth True if vertex values are to be smoothed. 1268 # @param precision The type of the result values (default float). 1269 def get_vertex_values(self, xy=True, 1270 smooth=None, 1271 precision=None): 1310 1272 """Return vertex values like an OBJ format i.e. one value per node. 1311 1273 … … 1316 1278 Mx3, where M is the number of triangles. Each row has three indices 1317 1279 defining the triangle and they correspond to elements in the arrays 1318 X, Y and A. 1319 1320 if smooth is True, vertex values corresponding to one common 1321 coordinate set will be smoothed by taking the average of vertex values for each node. 1322 In this case vertex coordinates will be 1323 de-duplicated corresponding to the original nodes as obtained from 1324 the method general_mesh.get_nodes() 1325 1326 If no smoothings is required, vertex coordinates and values will 1327 be aggregated as a concatenation of values at 1328 vertices 0, vertices 1 and vertices 2. This corresponds to 1329 the node coordinates obtained from the method 1330 general_mesh.get_vertex_coordinates() 1331 1280 X, Y and A. 1281 1282 If smooth is True, vertex values corresponding to one common coordinate 1283 set will be smoothed by taking the average of vertex values for each 1284 node. In this case vertex coordinates will be de-duplicated 1285 corresponding to the original nodes as obtained from the method 1286 general_mesh.get_nodes() 1287 1288 If no smoothings is required, vertex coordinates and values will be 1289 aggregated as a concatenation of values at vertices 0, vertices 1 and 1290 vertices 2. This corresponds to the node coordinates obtained from the 1291 method general_mesh.get_vertex_coordinates() 1332 1292 1333 1293 Calling convention 1334 1294 if xy is True: 1335 X, Y,A,V = get_vertex_values1295 X, Y, A, V = get_vertex_values 1336 1296 else: 1337 A,V = get_vertex_values 1338 1339 """ 1340 1297 A, V = get_vertex_values 1298 """ 1341 1299 1342 1300 if smooth is None: … … 1349 1307 if precision is None: 1350 1308 precision = num.Float 1351 1352 1309 1353 1310 if smooth is True: 1354 # Ensure continuous vertex values by averaging 1355 # values at each node 1356 1311 # Ensure continuous vertex values by averaging values at each node 1357 1312 V = self.domain.get_triangles() 1358 1313 N = self.domain.number_of_full_nodes # Ignore ghost nodes if any 1359 1314 A = num.zeros(N, num.Float) 1360 points = self.domain.get_nodes() 1361 1315 points = self.domain.get_nodes() 1316 1362 1317 if 1: 1363 1318 # Fast C version … … 1367 1322 A) 1368 1323 A = A.astype(precision) 1369 else: 1370 1324 else: 1371 1325 # Slow Python version 1372 1373 1326 current_node = 0 1374 1327 k = 0 # Track triangles touching on node … … 1376 1329 for index in self.domain.vertex_value_indices: 1377 1330 if current_node == N: 1378 msg = 'Current node exceeding number of nodes (%d) ' % (N)1331 msg = 'Current node exceeding number of nodes (%d) ' % N 1379 1332 raise msg 1380 1381 1382 1333 1383 1334 k += 1 1384 1335 1385 1336 volume_id = index / 3 1386 1337 vertex_id = index % 3 1387 1388 #assert V[volume_id, vertex_id] == current_node 1389 1338 1390 1339 v = self.vertex_values[volume_id, vertex_id] 1391 1340 total += v 1392 1341 1393 #print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)1394 1342 if self.domain.number_of_triangles_per_node[current_node] == k: 1395 1343 A[current_node] = total/k 1396 1397 1344 1398 1345 # Move on to next node 1399 1346 total = 0.0 1400 1347 k = 0 1401 1348 current_node += 1 1402 1403 1404 1405 1349 else: 1406 # Return disconnected internal vertex values 1350 # Return disconnected internal vertex values 1407 1351 V = self.domain.get_disconnected_triangles() 1408 1352 points = self.domain.get_vertex_coordinates() 1409 1353 A = self.vertex_values.flat.astype(precision) 1410 1354 1411 1412 # Return 1355 # Return 1413 1356 if xy is True: 1414 1357 X = points[:,0].astype(precision) 1415 1358 Y = points[:,1].astype(precision) 1416 1359 1417 1360 return X, Y, A, V 1418 1361 else: 1419 return A, V 1420 1421 1422 1362 return A, V 1363 1364 ## 1365 # @brief Extrapolate conserved quantities from centroid. 1423 1366 def extrapolate_first_order(self): 1424 """Extrapolate conserved quantities from centroid to 1425 vertices and edges for each volume using 1426 first order scheme. 1367 """Extrapolate conserved quantities from centroid to vertices and edges 1368 for each volume using first order scheme. 1427 1369 """ 1428 1370 … … 1438 1380 self.y_gradient *= 0.0 1439 1381 1440 1382 ## 1383 # @brief Compute the integral of quantity across entire domain. 1384 # @return The integral. 1441 1385 def get_integral(self): 1442 """Compute the integral of quantity across entire domain 1443 """ 1386 """Compute the integral of quantity across entire domain.""" 1387 1444 1388 areas = self.domain.get_areas() 1445 1389 integral = 0 … … 1451 1395 return integral 1452 1396 1397 ## 1398 # @brief get the gradients. 1453 1399 def get_gradients(self): 1454 """Provide gradients. Use compute_gradients first 1455 """ 1400 """Provide gradients. Use compute_gradients first.""" 1456 1401 1457 1402 return self.x_gradient, self.y_gradient 1458 1403 1459 1404 ## 1405 # @brief ?? 1406 # @param timestep ?? 1460 1407 def update(self, timestep): 1461 1408 # Call correct module function … … 1463 1410 return update(self, timestep) 1464 1411 1412 ## 1413 # @brief ?? 1465 1414 def compute_gradients(self): 1466 1415 # Call correct module function … … 1468 1417 return compute_gradients(self) 1469 1418 1419 ## 1420 # @brief ?? 1470 1421 def limit(self): 1471 1422 # Call correct module depending on whether … … 1473 1424 limit_old(self) 1474 1425 1426 ## 1427 # @brief ?? 1475 1428 def limit_vertices_by_all_neighbours(self): 1476 1429 # Call correct module function … … 1478 1431 limit_vertices_by_all_neighbours(self) 1479 1432 1433 ## 1434 # @brief ?? 1480 1435 def limit_edges_by_all_neighbours(self): 1481 1436 # Call correct module function … … 1483 1438 limit_edges_by_all_neighbours(self) 1484 1439 1440 ## 1441 # @brief ?? 1485 1442 def limit_edges_by_neighbour(self): 1486 1443 # Call correct module function 1487 1444 # (either from this module or C-extension) 1488 limit_edges_by_neighbour(self) 1489 1445 limit_edges_by_neighbour(self) 1446 1447 ## 1448 # @brief ?? 1490 1449 def extrapolate_second_order(self): 1491 1450 # Call correct module function … … 1493 1452 compute_gradients(self) 1494 1453 extrapolate_from_gradient(self) 1495 1454 1455 ## 1456 # @brief ?? 1496 1457 def extrapolate_second_order_and_limit_by_edge(self): 1497 1458 # Call correct module function … … 1499 1460 extrapolate_second_order_and_limit_by_edge(self) 1500 1461 1462 ## 1463 # @brief ?? 1501 1464 def extrapolate_second_order_and_limit_by_vertex(self): 1502 1465 # Call correct module function … … 1504 1467 extrapolate_second_order_and_limit_by_vertex(self) 1505 1468 1469 ## 1470 # @brief ?? 1471 # @param bound ?? 1506 1472 def bound_vertices_below_by_constant(self, bound): 1507 1473 # Call correct module function … … 1509 1475 bound_vertices_below_by_constant(self, bound) 1510 1476 1477 ## 1478 # @brief ?? 1479 # @param quantity ?? 1511 1480 def bound_vertices_below_by_quantity(self, quantity): 1512 1481 # Call correct module function … … 1515 1484 # check consistency 1516 1485 assert self.domain == quantity.domain 1517 bound_vertices_below_by_quantity(self, quantity) 1518 1486 bound_vertices_below_by_quantity(self, quantity) 1487 1488 ## 1489 # @brief ?? 1519 1490 def backup_centroid_values(self): 1520 1491 # Call correct module function … … 1522 1493 backup_centroid_values(self) 1523 1494 1524 def saxpy_centroid_values(self,a,b): 1495 ## 1496 # @brief ?? 1497 # @param a ?? 1498 # @param b ?? 1499 def saxpy_centroid_values(self, a, b): 1525 1500 # Call correct module function 1526 1501 # (either from this module or C-extension) 1527 saxpy_centroid_values(self,a,b) 1528 1529 #Conserved_quantity = Quantity 1530 1502 saxpy_centroid_values(self, a, b) 1503 1504 1505 ## 1506 # @brief OBSOLETE! 1531 1507 class Conserved_quantity(Quantity): 1532 """Class conserved quantity being removed, use Quantity 1533 1534 """ 1508 """Class conserved quantity being removed, use Quantity.""" 1535 1509 1536 1510 def __init__(self, domain, vertex_values=None): 1537 #Quantity.__init__(self, domain, vertex_values)1538 1539 1511 msg = 'ERROR: Use Quantity instead of Conserved_quantity' 1540 1541 1512 raise Exception, msg 1542 1513 1543 1514 1515 ###### 1516 # Prepare the C extensions. 1517 ###### 1544 1518 1545 1519 from anuga.utilities import compile 1546 if compile.can_use_C_extension('quantity_ext.c'): 1547 # Underlying C implementations can be accessed 1520 1521 if compile.can_use_C_extension('quantity_ext.c'): 1522 # Underlying C implementations can be accessed 1548 1523 1549 1524 from quantity_ext import \ … … 1564 1539 interpolate_from_vertices_to_edges,\ 1565 1540 interpolate_from_edges_to_vertices,\ 1566 update 1541 update 1567 1542 else: 1568 msg = 'C implementations could not be accessed by %s.\n ' % __file__1543 msg = 'C implementations could not be accessed by %s.\n ' % __file__ 1569 1544 msg += 'Make sure compile_all.py has been run as described in ' 1570 1545 msg += 'the ANUGA installation guide.' 1571 1546 raise Exception, msg 1572 1573 1574
Note: See TracChangeset
for help on using the changeset viewer.
