Changeset 1632 for inundation/ga/storm_surge/pyvolution

Timestamp:

Jul 22, 2005, 5:59:19 PM (20 years ago)

Author:

ole

Message:

Interpolation work and some refactoring

Location:

inundation/ga/storm_surge/pyvolution

Files:

• general_mesh.py (modified) (4 diffs)
• interpolate_sww.py (modified) (2 diffs)
• least_squares.py (modified) (1 diff)
• quantity.py (modified) (1 diff)
• sparse.py (modified) (1 diff)
• test_least_squares.py (modified) (4 diffs)
• test_mesh.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -173 +173 @@
 
 
-    def get_vertex_coordinates(self):
+    def get_vertex_coordinates(self, obj = False):
         """Return all vertex coordinates.
         Return all vertex coordinates for all triangles as an Nx6 array
         (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
-        """
-        return self.vertex_coordinates
+
+        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
+        FIXME, we might make that the default.
+
+        
+        """
+
+        if obj is True:
+            from Numeric import concatenate, reshape
+            V = self.vertex_coordinates
+            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
+
+            N = V.shape[0]
+            return reshape(V, (3*N, 2))
+        else:    
+            return self.vertex_coordinates
 
 
@@ -193 +207 @@
         """
 
-        #FIXME: Perhaps they should be ordered as in obj files??
+        #FIXME (Ole): Perhaps they should be ordered as in obj files??
         #See quantity.get_vertex_values
+        #FIXME (Ole) - oh yes they should
 
         from Numeric import zeros, Float
@@ -210 +225 @@
         return vertex_coordinates
 
-    def get_vertices(self,  indexes=None):
+    def get_vertices(self, indexes=None):
         """Get connectivity
         indexes is the set of element ids of interest
@@ -221 +236 @@
 
         return  take(self.triangles, indexes)
+
+    #FIXME - merge these two
+    def get_triangles(self, obj = False):
+        """Get connetivity
+        Return triangles (triplets of indices into point coordinates)
+        
+        """
+
+        if obj is True:
+            from Numeric import array, reshape, Int        
+            m = len(self)  #Number of triangles
+            M = 3*m        #Total number of unique vertices
+            T = reshape(array(range(M)).astype(Int), (m,3))
+        else:
+            T = self.triangles
+
+        return T     
+
+    
 
     def get_unique_vertices(self,  indexes=None):

inundation/ga/storm_surge/pyvolution/interpolate_sww.py

@@ -84 +84 @@
         
             
-    def read_sww(self,file_name, quantity_name):
+    def read_sww(self, file_name, quantity_name):
         """
         Read in an sww file.
@@ -136 +136 @@
         The time list is defined
         """
-        from load_mesh.loadASCII import  export_points_file
+        from load_mesh.loadASCII import export_points_file
 
-        
         xya_dict = {}
         xya_dict['pointlist'] = self.point_coordinates
         xya_dict['attributelist'] = self.interpolated_quantity
         export_points_file(point_file, xya_dict)
+        
         
 #-------------------------------------------------------------

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -279 +279 @@
           vertex_coordinates: List of coordinate pairs [xi, eta] of
           points constituting mesh (or a an m x 2 Numeric array)
+          Points may appear multiple times
+          (e.g. if vertices have discontinuities)
 
           triangles: List of 3-tuples (or a Numeric array) of

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -538 +538 @@
         else:
             #Don't smooth
+            #obj machinery moved to general_mesh
 
             # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
             # These vert_id's will relate to the verts created below
-            m = len(self.domain)  #Number of volumes
-            M = 3*m        #Total number of unique vertices
-            V = reshape(array(range(M)).astype(Int), (m,3))
+            #m = len(self.domain)  #Number of volumes
+            #M = 3*m        #Total number of unique vertices
+            #V = reshape(array(range(M)).astype(Int), (m,3))
+            
+            V = self.domain.get_triangles(obj=True)
+            #FIXME use get_vertices, when ready
 
             A = self.vertex_values.flat

inundation/ga/storm_surge/pyvolution/sparse.py

@@ -115 +115 @@
             B = array(other)
         except:
-            print 'FIXME: Only Numeric types implemented so far'
+            msg = 'FIXME: Only Numeric types implemented so far'
+            raise msg
+            
 
         #Assume numeric types from now on

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -925 +925 @@
 
 
+    def test_interpolation_from_discontinuous_vertex_values(self):
+        """test_interpolation_from_discontinuous_vertex_values.
+        This will test the format used internally in pyvolution and also
+        interpolation from sww files
+        """
+        
+        from mesh import Mesh
+
+
+        #Setup mesh used to represent discontinuous function
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [b, a, c,
+                  b, c, e,
+                  e, c, f,
+                  d, b, e]
+        
+        #bac, bce, ecf, dbe
+        triangles = [[0,1,2], [3,4,5], [6,7,8], [9,10,11]]
+
+        
+        vertex_values = [0.,0.,0.,1.,1.,1.,2.,2.,2.,7.,3.,3.]
+                          
+        
+
+        #New datapoints where interpolated values are sought
+        data_points = [[0.0, 0.0],  #T0
+                       [0.5, 0.5],  #T0
+                       [1.5, 1.5],  #T1
+                       [2.5, 0.5],  #T2
+                       [0.0, 3.0],  #T3
+                       [1.0, 2.0],  #In between T1 and T3 (T1 is used) FIXME?
+                       [2.0, 1.0],  #In between T1 and T2 (T1 is used) FIXME?
+                       [1.0, 1.0]]  #In between T1 and T0 (T0 is used) FIXME?
+        
+
+
+
+        #Build interpolation matrix
+        interp = Interpolation(points, triangles, data_points)
+                               #, alpha=0.0, precrop = True)
+
+        #print interp.A.todense()
+        #print vertex_values
+
+        #Interpolate using fitted surface
+        z = interp.interpolate(vertex_values)
+
+        #print z
+
+        assert allclose(z, [0,0,1,2,5,1,1,0])
+
+
+
+
+
+
+
 
     def test_fit_and_interpolation_with_different_origins(self):
@@ -1167 +1230 @@
         mesh_dic = {}
         mesh_dic['vertices'] = [[0.0, 0.0],
-                                          [0.0, 5.0],
-                                          [5.0, 0.0]]
+                                [0.0, 5.0],
+                                [5.0, 0.0]]
         mesh_dic['triangles'] =  [[0, 2, 1]]
         mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
@@ -1210 +1273 @@
 
     def test_fit_to_msh_netcdf_fileII(self):
-        from load_mesh.loadASCII import import_mesh_file,export_mesh_file
+        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
         import tempfile
         import os
@@ -1217 +1280 @@
         mesh_dic = {}
         mesh_dic['vertices'] = [[0.0, 0.0],
-                                          [0.0, 5.0],
-                                          [5.0, 0.0]]
+                                [0.0, 5.0],
+                                [5.0, 0.0]]
         mesh_dic['triangles'] =  [[0, 2, 1]]
         mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]

inundation/ga/storm_surge/pyvolution/test_mesh.py

@@ -70 +70 @@
         V = mesh.get_vertex_coordinates()
         assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
+
+        V = mesh.get_vertex_coordinates(obj=True)
+        assert allclose(V, [ [0.0, 0.0],
+                             [4.0, 0.0],
+                             [0.0, 3.0] ])
 
         V0 = mesh.get_vertex_coordinate(0, 0)