Changeset 5866


Ignore:
Timestamp:
Oct 27, 2008, 3:20:26 PM (16 years ago)
Author:
ole
Message:

Optimised 'get_interpolated_values' in a more general and robust way than
what was done in changeset:5864

Now, each mesh precomputes and maintains a local copy of
its interpolation object.

Location:
anuga_core/source/anuga
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r5862 r5866  
    968968       
    969969
     970    def get_interpolation_object(self):
     971        """Get object I that will allow linear interpolation using this mesh
     972       
     973        This is a time consuming process but it needs only to be
     974        once for the mesh.
     975       
     976        Interpolation can then be done using
     977       
     978        result = I.interpolate_block(vertex_values, interpolation_points)       
     979       
     980        where vertex values have been obtained from a quantity using
     981        vertex_values, triangles = self.get_vertex_values()
     982        """
     983
     984        if hasattr(self, 'interpolation_object'):
     985            I = self.interpolation_object
     986        else:
     987            from anuga.fit_interpolate.interpolate import Interpolate
     988                       
     989            # Get discontinuous mesh - this will match internal
     990            # representation of vertex values
     991            triangles = self.get_disconnected_triangles()
     992            vertex_coordinates = self.get_vertex_coordinates()
     993
     994            I = Interpolate(vertex_coordinates, triangles)
     995            self.interpolation_object = I
     996       
     997        return I               
    970998       
    971999

  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r5864 r5866  
    10291029
    10301030
    1031     def obsolete_get_interpolated_values(self, interpolation_points,
     1031    def get_interpolated_values(self, interpolation_points,
     1032                                use_cache=False,
    10321033                                verbose=False):
    10331034        """ Get values at interpolation points
    10341035       
    1035         If interpolation points have been given previously, the
    1036         associated matrices will be reused to save time.
    1037        
    10381036        The argument interpolation points must be given as either a
    10391037        list of absolute UTM coordinates or a geospatial data object.
    10401038        """
    1041 
    1042         # FIXME (Ole): Could do with an input check (should be generalised
    1043         # and used widely)
    1044         # That will check that interpolation points is either a list of
    1045         # points, Nx2 array, or geospatial
    1046        
    1047         # Ensure points are converted to coordinates relative to mesh origin
    1048         # FIXME (Ole): This could all be refactored using the
     1039       
     1040
     1041        # FIXME (Ole): Points might be converted to coordinates relative to mesh origin
     1042        # This could all be refactored using the
    10491043        # 'change_points_geo_ref' method of Class geo_reference.
    10501044        # The purpose is to make interpolation points relative
    … …  
    10541048       
    10551049       
    1056         if isinstance(interpolation_points, Geospatial_data):       
    1057             # Ensure interpolation points are in absolute UTM coordinates
    1058             interpolation_points = interpolation_points.get_data_points(absolute=True)
    1059            
    1060         # Reconcile interpolation points with georeference of domain
    1061         interpolation_points = self.domain.geo_reference.get_relative(interpolation_points)
    1062         interpolation_points = ensure_numeric(interpolation_points)
    1063        
    1064        
    1065        
    1066         # Get internal (discontinuous) triangles for use with the
    1067         # interpolation object.
    1068         x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
    1069                                                                 smooth=False)
    1070         # FIXME (Ole): This concat should roll into get_vertex_values
    1071         vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
    1072                                          axis=1)
    1073 
    1074         # Use caching to reuse interpolation information
    1075         from anuga.fit_interpolate.interpolate import interpolate
    1076        
    1077         return interpolate(vertex_coordinates,
    1078                            triangles,
    1079                            vertex_values,
    1080                            interpolation_points,
    1081                            use_cache=True,
    1082                            verbose=verbose)
    1083        
    1084        
    1085 
    1086     def get_interpolated_values(self, interpolation_points,
    1087                                 verbose=False):
    1088         """ Get values at interpolation points
    1089        
    1090         If interpolation points have been given previously, the
    1091         associated matrices will be reused to save time.
    1092        
    1093         The argument interpolation points must be given as either a
    1094         list of absolute UTM coordinates or a geospatial data object.
    1095         """
    1096 
    1097         # FIXME (Ole): Could do with an input check (should be generalised
    1098         # and used widely)
    1099         # That will check that interpolation points is either a list of
    1100         # points, Nx2 array, or geospatial
    1101        
    1102         # Ensure points are converted to coordinates relative to mesh origin
    1103         # FIXME (Ole): This could all be refactored using the
    1104         # 'change_points_geo_ref' method of Class geo_reference.
    1105         # The purpose is to make interpolation points relative
    1106         # to the mesh origin.
    1107         #
    1108         # Speed is also a consideration here.
    1109        
    1110         # FIXME (Ole): The reuse of I is a little rough, but it sure does speed up
    1111         # frequent interpolations such as those found in culvert_class.
    1112        
    1113         use_cache = True
    1114        
    1115         # Get internal (discontinuous) triangles for use with the
    1116         # interpolation object.
    1117         x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
    1118                                                                 smooth=False)
    1119         # FIXME (Ole): This concat should roll into get_vertex_values
    1120         vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
    1121                                          axis=1)                                                               
    1122            
    1123            
     1050        # Ensure that interpolation points is either a list of
     1051        # points, Nx2 array, or geospatial and convert to Numeric array
    11241052        if isinstance(interpolation_points, Geospatial_data):       
    11251053            # Ensure interpolation points are in absolute UTM coordinates
    … …  
    11291057        interpolation_points = self.domain.geo_reference.get_relative(interpolation_points)
    11301058        interpolation_points = ensure_numeric(interpolation_points)
    1131        
    1132         # Check if same interpolation points have been used before
    1133         if hasattr(self, 'stored_I') and\
    1134                 alltrue(interpolation_points == self.stored_interpolation_points):
    1135             I = self.stored_I
    1136         else:   
    1137        
    1138             # Use caching to reuse interpolation information
    1139             from anuga.fit_interpolate.interpolate import Interpolate
    1140             from anuga.caching import cache
     1059
     1060       
     1061        # Get internal representation (disconnected) of vertex values
     1062        vertex_values, triangles = self.get_vertex_values(xy=False,
     1063                                                          smooth=False)               
    11411064   
    1142 
    1143             # Create interpolation object with matrix
    1144             args = (ensure_numeric(vertex_coordinates, Float),
    1145                     ensure_numeric(triangles))
    1146                      
    1147             if use_cache is True:
    1148                 I = cache(Interpolate, args,
    1149                           verbose=verbose)
    1150             else:
    1151                 I = apply(Interpolate, args)
    1152                          
    1153             self.stored_I = I   
    1154             self.stored_interpolation_points = interpolation_points
    1155            
    1156            
    1157    
    1158         # Call interpolate method with interpolation points
     1065        # Get possibly precomputed interpolation object
     1066        I = self.domain.get_interpolation_object()
     1067
     1068        # Call interpolate method with interpolation points               
    11591069        result = I.interpolate_block(vertex_values, interpolation_points,
    11601070                                     use_cache=use_cache,
    … …  
    11641074       
    11651075       
    1166         return interpolate(vertex_coordinates,
    1167                            triangles,
    1168                            vertex_values,
    1169                            interpolation_points,
    1170                            use_cache=True,
    1171                            verbose=verbose)
    1172                            
    11731076
    11741077
    … …  
    11771080                   location='vertices',
    11781081                   indices=None,
     1082                   use_cache=False,
    11791083                   verbose=False):
    11801084        """get values for quantity
    … …  
    12251129        if interpolation_points is not None:
    12261130            return self.get_interpolated_values(interpolation_points,
     1131                                                use_cache=use_cache,
    12271132                                                verbose=verbose)
    12281133       
    … …  
    14351340
    14361341        else:
    1437             # Allow discontinuous vertex values
     1342            # Return disconnected internal vertex values
    14381343            V = self.domain.get_disconnected_triangles()
    14391344            points = self.domain.get_vertex_coordinates()

  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5862 r5866  
    417417        ymomentum = self.get_quantity('ymomentum')       
    418418       
    419         uh = xmomentum.get_values(interpolation_points=midpoints)
    420         vh = ymomentum.get_values(interpolation_points=midpoints)       
     419        uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
     420        vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
    421421       
    422422        # Compute and sum flows across each segment
    … …  
    488488        ymomentum = self.get_quantity('ymomentum')       
    489489
    490         w = stage.get_values(interpolation_points=midpoints)
    491         z = elevation.get_values(interpolation_points=midpoints)       
    492         uh = xmomentum.get_values(interpolation_points=midpoints)
    493         vh = ymomentum.get_values(interpolation_points=midpoints)       
     490        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
     491        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)       
     492        uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
     493        vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
    494494        h = w-z # Depth
    495495       
Note: See TracChangeset for help on using the changeset viewer.