Changeset 4130

Timestamp:

Jan 5, 2007, 2:14:16 PM (18 years ago)

Author:

duncan

Message:

Added first iteration of blocking to quantity.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (3 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (9 diffs)
• fit_interpolate/fit.py (modified) (2 diffs)
• fit_interpolate/test_fit.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -17 +17 @@
 from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
      argmax, allclose, take, reshape
+
 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
 from anuga.geospatial_data.geospatial_data import Geospatial_data
+from anuga.fit_interpolate.fit import fit_to_mesh
 
 class Quantity:
@@ -655 +657 @@
 
 
-        from anuga.fit_interpolate.fit import fit_to_mesh
         from anuga.coordinate_transforms.geo_reference import Geo_reference
 
@@ -727 +728 @@
         """
 
-        from load_mesh.loadASCII import import_points_file
-        from anuga.geospatial_data.geospatial_data import\
-             points_dictionary2geospatial_data
-
         from types import StringType
         msg = 'Filename must be a text string'
         assert type(filename) == StringType, msg
 
-        geospatial_data = Geospatial_data(filename)
+
+        if location != 'vertices':
+            msg = 'set_values_from_points is only defined for '+\
+                  'location=\'vertices\''
+            raise ms
+
+        coordinates = self.domain.get_nodes(absolute=True)
+        triangles = self.domain.triangles      #FIXME
+
         
-        self.set_values_from_geospatial_data(geospatial_data,
-                                             alpha,
-                                             location, indices,
-                                             verbose = verbose,
-                                             use_cache = use_cache)
+        # FIXME handle attribute name 
+        vertex_attributes = fit_to_mesh(coordinates, triangles,filename,
+                                alpha=alpha) #, max_read_lines=max_read_lines)
+        
+        #Call underlying method using array values
+        self.set_values_from_array(vertex_attributes,
+                                   location, indices, verbose)
+
+        #geospatial_data = Geospatial_data(filename)
+        
+        #self.set_values_from_geospatial_data(geospatial_data,
+        #                                     alpha,
+         #                                    location, indices,
+          #                                   verbose = verbose,
+           #                                  use_cache = use_cache)
 
    

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -3 +3 @@
 import unittest
 from math import sqrt, pi
-
+import tempfile
 
 from quantity import *
@@ -499 +499 @@
                        [ 3.0, 1.0]]
 
-        z = linear_function(data_points)
-
-
-        #Create pts file
-        from load_mesh.loadASCII import export_points_file
-        ptsfile = 'testptsfile.pts'
+        data_geo_spatial = Geospatial_data(data_points,
+                         geo_reference = Geo_reference(56, 0, 0))
+        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
+        attributes = linear_function(data_points_absolute)
         att = 'spam_and_eggs'
-        points_dict = {'pointlist': data_points,
-                       'attributelist': {att: z}}
-
-        export_points_file(ptsfile, points_dict)
+        
+        #Create .txt file
+        ptsfile = tempfile.mktemp(".txt")
+        file = open(ptsfile,"w")
+        file.write(" x,y," + att + " \n")
+        for data_point, attribute in map(None, data_points_absolute
+                                         ,attributes):
+            row = str(data_point[0]) + ',' + str(data_point[1]) \
+                  + ',' + str(attribute)
+            file.write(row + "\n")
+        file.close()
+
 
         #Check that values can be set from file
@@ -551 +557 @@
         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
 
+        #absolute going in ..
         mesh4 = Domain(points, elements,
                        geo_reference = Geo_reference(56, 0, 0))
@@ -557 +564 @@
 
         #Get (enough) datapoints (relative to georef)
-        data_points = [[ 0.66666667, 0.66666667],
+        data_points_rel = [[ 0.66666667, 0.66666667],
                        [ 1.33333333, 1.33333333],
                        [ 2.66666667, 0.66666667],
@@ -571 +578 @@
                        [ 3.0, 1.0]]
 
-        z = linear_function(data_points)
-
-
-        #Create pts file
-        from load_mesh.loadASCII import export_points_file
-
-        ptsfile = 'testptsfile.pts'
+        data_geo_spatial = Geospatial_data(data_points_rel,
+                         geo_reference = Geo_reference(56, x0, y0))
+        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
+        attributes = linear_function(data_points_absolute)
         att = 'spam_and_eggs'
-
-        points_dict = {'pointlist': data_points,
-                       'attributelist': {att: z},
-                       'geo_reference': Geo_reference(zone = 56,
-                                                      xllcorner = x0,
-                                                      yllcorner = y0)}
-
-        export_points_file(ptsfile, points_dict)
-
+        
+        #Create .txt file
+        ptsfile = tempfile.mktemp(".txt")
+        file = open(ptsfile,"w")
+        file.write(" x,y," + att + " \n")
+        for data_point, attribute in map(None, data_points_absolute
+                                         ,attributes):
+            row = str(data_point[0]) + ',' + str(data_point[1]) \
+                  + ',' + str(attribute)
+            file.write(row + "\n")
+        file.close()
+
+        #file = open(ptsfile, 'r')
+        #lines = file.readlines()
+        #file.close()
+     
 
         #Check that values can be set from file
         quantity.set_values(filename = ptsfile,
                             attribute_name = att, alpha = 0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates() - [x0, y0])
+        answer = linear_function(quantity.domain.get_vertex_coordinates())
 
         assert allclose(quantity.vertex_values.flat, answer)
@@ -630 +641 @@
         quantity = Quantity(mesh4)
 
-        #Get (enough) datapoints (relative to georef)
+        #Get (enough) datapoints
         data_points = [[ x0+0.66666667, y0+0.66666667],
                        [ x0+1.33333333, y0+1.33333333],
@@ -645 +656 @@
                        [ x0+3.0, y0+1.0]]
 
-        z = linear_function(data_points)
-
-
-        #Create pts file
-        from load_mesh.loadASCII import export_points_file
-
-        ptsfile = 'testptsfile.pts'
+
+        data_geo_spatial = Geospatial_data(data_points,
+                         geo_reference = Geo_reference(56, 0, 0))
+        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
+        attributes = linear_function(data_points_absolute)
         att = 'spam_and_eggs'
-
-        points_dict = {'pointlist': data_points,
-                       'attributelist': {att: z},
-                       'geo_reference': Geo_reference(zone = 56,
-                                                      xllcorner = 0,
-                                                      yllcorner = 0)}
-
-        export_points_file(ptsfile, points_dict)
+        
+        #Create .txt file
+        ptsfile = tempfile.mktemp(".txt")
+        file = open(ptsfile,"w")
+        file.write(" x,y," + att + " \n")
+        for data_point, attribute in map(None, data_points_absolute
+                                         ,attributes):
+            row = str(data_point[0]) + ',' + str(data_point[1]) \
+                  + ',' + str(attribute)
+            file.write(row + "\n")
+        file.close()
 
 
@@ -666 +678 @@
         quantity.set_values(filename = ptsfile,
                             attribute_name = att, alpha = 0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates() + [x0, y0])
+        answer = linear_function(quantity.domain. \
+                                 get_vertex_coordinates(absolute=True))
 
 
@@ -1690 +1703 @@
     suite = unittest.makeSuite(Test_Quantity, 'test')
     #print "restricted test"
-    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
+    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_from_file_with_georef2')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/fit_interpolate/fit.py

@@ -415 +415 @@
 def fit_to_mesh(vertex_coordinates,
                 triangles,
-                point_coordinates=None,
+                point_coordinates, # this can also be a .csv/.txt file name
                 point_attributes=None,
-                alpha = DEFAULT_ALPHA,
-                verbose = False,
-                acceptable_overshoot = 1.01,
-                mesh_origin = None,
-                data_origin = None,
+                alpha=DEFAULT_ALPHA,
+                verbose=False,
+                acceptable_overshoot=1.01,
+                mesh_origin=None,
+                data_origin=None,
+                max_read_lines=None,
                 use_cache = False):
     """
@@ -477 +478 @@
     vertex_attributes = interp.fit(point_coordinates,
                                    point_attributes,
-                                   point_origin = data_origin,
-                                   verbose = verbose)
+                                   point_origin=data_origin,
+                                   max_read_lines=max_read_lines,
+                                   verbose=verbose)
 
         

anuga_core/source/anuga/fit_interpolate/test_fit.py

@@ -235 +235 @@
         
         f = interp.fit(fileName, max_read_lines=2)
+        answer = linear_function(vertices)
+        #print "f\n",f
+        #print "answer\n",answer
+        assert allclose(f, answer)
+
+    def test_fit_2_mesh(self):
+
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0,1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0,-2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
+
+
+        fileName = tempfile.mktemp(".ddd")
+        file = open(fileName,"w")
+        file.write(" x,y, elevation \n\
+-2.0, 2.0, 0.\n\
+-1.0, 1.0, 0.\n\
+0.0, 2.0 , 2.\n\
+1.0, 1.0 , 2.\n\
+ 2.0,  1.0 ,3. \n\
+ 0.0,  0.0 , 0.\n\
+ 1.0,  0.0 , 1.\n\
+ 0.0,  -1.0, -1.\n\
+ -0.2, -0.5, -0.7\n\
+ -0.9, -1.5, -2.4\n\
+ 0.5,  -1.9, -1.4\n\
+ 3.0,  1.0 , 4.\n")
+        file.close()
+        
+        f = fit_to_mesh(vertices, triangles,fileName,
+                                alpha=0.0, max_read_lines=2)
         answer = linear_function(vertices)
         #print "f\n",f