Changeset 328 for inundation/ga/storm_surge/pyvolution/least_squares.py

Timestamp:

Sep 21, 2004, 1:26:23 PM (20 years ago)

Author:

duncan

Message:

adding reading .tsh files to smoothing

File:

• inundation/ga/storm_surge/pyvolution/least_squares.py (modified) (5 diffs)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -15 +15 @@
 EPSILON = 1.0e-13
 
-
 from general_mesh import General_Mesh
 from Numeric import zeros, array, Float, Int, dot,transpose
 from LinearAlgebra import solve_linear_equations
+
+
+def smooth_attributes_to_mesh(vertex_coordinates,
+                              triangles,
+                              point_coordinates,
+                              point_attributes):
+    """
+    """
+    interp = Interpolation(vertex_coordinates,
+                              triangles,
+                              point_coordinates)
+    vertex_attributes = interp.smooth_attributes_to_mesh(point_attributes)
+    return vertex_attributes
 
 
@@ -52 +64 @@
         self.mesh = General_Mesh(vertex_coordinates, triangles)
 
-        self.A_m = zeros((1,1), Float) #Not initialised
+        self.A = zeros((1,1), Float) #Not initialised
 
         if point_coordinates:
-            self.build_A_m(point_coordinates)
+            self.build_A(point_coordinates)
         
-    def build_A_m(self, point_coordinates):
+    def build_A(self, point_coordinates):
+        """
+        "A" is the matrix of hat functions.  It has 1 row per data point and 1 column per vertex.
+
+        Post Condition
+        "A" has been initialised
+        """
 
         #Convert input to Numeric arrays
@@ -66 +84 @@
         n = point_coordinates.shape[0]               #Number of data points         
 
-        self.A_m = zeros((n,m), Float)
+        self.A = zeros((n,m), Float)
 
         #Compute matrix elements
@@ -99 +117 @@
                 if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
 
-                    #Assign values to A_m
+                    #Assign values to A
                     j = self.mesh.triangles[k,0] #Global vertex id
-                    self.A_m[i, j] = sigma0
+                    self.A[i, j] = sigma0
 
                     j = self.mesh.triangles[k,1] #Global vertex id
-                    self.A_m[i, j] = sigma1
+                    self.A[i, j] = sigma1
 
                     j = self.mesh.triangles[k,2] #Global vertex id
-                    self.A_m[i, j] = sigma2
+                    self.A[i, j] = sigma2
     
 
-    def smooth_att_to_mesh(self, z):
+    def smooth_attributes_to_mesh(self, z):
         """ Linearly interpolate many points with an attribute to the vertices.
         Pre Condition:
-          A_m has been initialised
+          "A" has been initialised
           
         Inputs:
@@ -121 +139 @@
         #Convert input to Numeric arrays
         z = array(z).astype(Float)
-        At_m = transpose(self.A_m)
+        At = transpose(self.A)
         #print "z", z
-        #print "At_m",At_m
-        self.AtA_m = dot(At_m, self.A_m) 
-        Atz_m = dot(At_m, z)
-        f = solve_linear_equations(self.AtA_m,Atz_m)
+        #print "At",At
+        self.AtA = dot(At, self.A) 
+        Atz = dot(At, z)
+        f = solve_linear_equations(self.AtA,Atz)
         return f
 
-    def smooth_atts_to_mesh(self, z_m):
-        """ Linearly interpolate many points with attributes to the vertices.
+        
+    def interpolate_attributes_to_points(self, f):
+        """ Linearly interpolate vertices with attributes to the points.
         Pre Condition:
-          A_m has been initialised
-          
-        Inputs:
-          z_m: A Matrix of data at the point_coordinates.
-               One attribute per colum.
-          
-        """
-        #Convert input to Numeric arrays
-        z_m = array(z_m).astype(Float)
-
-        
-        #Build n x m interpolation matrix        
-        m = self.mesh.coordinates.shape[0] #Number of vertices
-        n = z_m.shape[1]               #Number of data points         
-
-        f_m = zeros((n,m), Float)
-        
-        for i in range(z_m.shape[1]):
-            f_m[:,i] = smooth_att_to_mesh(z_m[:,i])
-        
-        
-    def smooth_att_to_points(self, f):
-        """ Linearly interpolate vertices with an attribute to the points.
-        Pre Condition:
-          A_m has been initialised
+          A has been initialised
         
         Inputs:
-          z: Vector of data at the point_coordinates.  One value per point.
-          
-        """
-        return dot(self.A_m,f)
-
-    def smooth_atts_to_points(self, f_m):
-        """ Linearly interpolate vertices with attributes to the points.
-        Pre Condition:
-          A_m has been initialised
-        
-        Inputs:
-          f_m: Matrix of data at the  vertices. 
+          f: Matrix of data at the  vertices. 
                One attribute per colum. 
           
         """
-        #Convert input to Numeric arrays
-        z_m = array(z_m).astype(Float)
+        return dot(self.A,f)
 
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    
+    import os, sys
+    from load_mesh.loadASCII import mesh_file_to_mesh_dictionary
+    usage = "usage: %s pmesh_file_name" %         os.path.basename(sys.argv[0])
+
+    if len(sys.argv) < 2:
+        print usage
+    else:
+        mesh_file_name = sys.argv[1]
+
+        meshdic = mesh_file_to_mesh_dictionary(mesh_file_name)
+        vertex_coordinates = meshdic['generatedpointlist']
+        triangles = meshdic['generatedtrianglelist']
         
-        #Build n x m interpolation matrix        
-        m = self.mesh.coordinates.shape[0] #Number of vertices
-        n = z_m.shape[1]               #Number of data points         
-
-        f_m = zeros((n,m), Float)
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 3.0]
+        d3 = [3.0,1.0]
+        z1 = 2
+        z2 = 4
+        z3 = 4
+        data_coords = [ d1, d2, d3]
+        z = [z1, z2, z3]
         
-        for i in range(z_m.shape[1]):
-            f_m[:,i] = smooth_att_to_mesh(z_m[:,i])
+        f = smooth_attributes_to_mesh(vertex_coordinates,
+                                      triangles,
+                                      data_coords,
+                                      z)
+        print f
         
+