Changeset 346

Timestamp:

Sep 28, 2004, 4:18:34 PM (21 years ago)

Author:

duncan

Message:

added test for fit_to_mesh_file, removed bugs from loadASCII.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• least_squares.py (modified) (6 diffs)
• load_mesh/loadASCII.py (modified) (8 diffs)
• test_least_squares.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -25 +25 @@
 from LinearAlgebra import solve_linear_equations
 
+ALPHA_INITIAL = 0.001
+
+def fit_to_mesh_file(mesh_file, point_file, mesh_output_file, alpha= ALPHA_INITIAL):
+    """
+    Given a mesh file (tsh) and a point attribute file (xya), fit
+    point attributes to the mesh and write a mesh file with the
+    results.
+    """
+    from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
+                 load_xya_file, export_trianglulation_file
+    # load in the .tsh file
+    mesh_dic = mesh_file_to_mesh_dictionary(mesh_file)
+    vertex_coordinates = mesh_dic['generatedpointlist']
+    triangles = mesh_dic['generatedtrianglelist']
+        
+    # load in the .xya file
+    point_dict = load_xya_file(point_file)
+    point_coordinates = point_dict['pointlist']
+    point_attributes = point_dict['pointattributelist']
+        
+        
+    f = fit_to_mesh(vertex_coordinates,
+                    triangles,
+                    point_coordinates,
+                    point_attributes,
+                    alpha = alpha)
+    # convert array to list of lists
+    mesh_dic['generatedpointattributelist'] = f.tolist()
+    export_trianglulation_file(mesh_output_file, mesh_dic)
+        
 
 def fit_to_mesh(vertex_coordinates,
@@ -30 +60 @@
                               point_coordinates,
                               point_attributes,
-                              alpha = 0.0):
+                              alpha = ALPHA_INITIAL):
     """
     Fit a smooth surface to a trianglulation,
@@ -66 +96 @@
                  triangles,
                  point_coordinates = None,
-                 alpha = 0.001):
+                 alpha = ALPHA_INITIAL):
 
         
@@ -155 +185 @@
 
                 #FIXME: Maybe move out to test or something
-                epsilon = 1.0e-13
+                epsilon = 1.0e-6
                 assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
 
@@ -329 +359 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    
+    """
+    Load in a mesh and data points with attributes.
+    Fit the attributes to the mesh.
+    Save a new mesh file.
+    """
     import os, sys
-    from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
-                 load_xya_file, export_trianglulation_file
-    usage = "usage: %s mesh_input.tsh point.xya  mesh_output.tsh" %         os.path.basename(sys.argv[0])
+    usage = "usage: %s mesh_input.tsh point.xya  mesh_output.tsh alpha" %         os.path.basename(sys.argv[0])
 
     if len(sys.argv) < 4:
@@ -341 +373 @@
         point_file = sys.argv[2]
         mesh_output_file = sys.argv[3]
-
-        # load in the .tsh file
-        mesh_dic = mesh_file_to_mesh_dictionary(mesh_file)
-        vertex_coordinates = mesh_dic['generatedpointlist']
-        triangles = mesh_dic['generatedtrianglelist']
-
-        # load in the .xya file
-        point_dict = load_xya_file(point_file)
-        point_coordinates = point_dict['pointlist']
-        point_attributes = point_dict['pointattributelist']
-        
-        
-        f = fit_to_mesh(vertex_coordinates,
-                                      triangles,
-                                      point_coordinates,
-                                      point_attributes)
-        # convert array to list of lists
-        mesh_dic['generatedpointattributelist'] = f.tolist()
-        export_trianglulation_file(mesh_output_file, mesh_dic)
-        #FIXME  do unit test
-        
-        
-
-
-
-
-
-
-
-
-
-
-
+        if len(sys.argv) > 4:
+            alpha = sys.argv[4]
+        else:
+            alpha = ALPHA_INITIAL
+        fit_to_mesh_file(mesh_file, point_file, mesh_output_file, alpha)
+        
+
+
+
+
+
+
+
+
+
+
+

inundation/ga/storm_surge/pyvolution/load_mesh/loadASCII.py

@@ -42 +42 @@
     #for fragment in fragments:
     #    print fragment
-    NumOfVertices = fragments[0]
-    NumOfVertAttributes = fragments[1]
+    if fragments ==[]:
+        NumOfVertices = 0
+        NumOfVertAttributes = 0
+    else:
+        NumOfVertices = fragments[0]
+        NumOfVertAttributes = fragments[1]
     points = []
     pointattributes = []
@@ -140 +144 @@
 
 def read_mesh(fd):
+    """
+    Note, if a file has no mesh info, it can still be read - the meshdic
+    returned will be 'empty'.
+    """
     delimiter = " " # warning: split() calls are using default whitespace
     
@@ -147 +155 @@
     fragments = line.split()
     #for fragment in fragments:
-    #    print fragment
-    NumOfVertices = fragments[0]
-    NumOfVertAttributes = fragments[1]
+    if fragments ==[]:
+        NumOfVertices = 0
+        NumOfVertAttributes = 0
+    else:
+        NumOfVertices = fragments[0]
+        NumOfVertAttributes = fragments[1]
     points = []
     pointattributes = []
@@ -173 +184 @@
     #for fragment in fragments:
     #    print fragment
-    NumOfSegments = fragments[0]
+    if fragments ==[]:
+        NumOfSegments = 0
+    else:
+        NumOfSegments = fragments[0]
     segments = []
     segmentmarkers = []
@@ -199 +213 @@
     #for fragment in fragments:
     #    print fragment
-    numOfHoles = fragments[0]
+    if fragments ==[]:
+        numOfHoles = 0
+    else:
+        numOfHoles = fragments[0]
     holes = []
     for index in range(int(numOfHoles)):        
@@ -216 +233 @@
     #for fragment in fragments:
     #    print fragment
-    numOfRegions = fragments[0]
+    if fragments ==[]:
+        numOfRegions = 0
+    else:
+        numOfRegions = fragments[0]
     regions = []
     regionattributes = []
@@ -321 +341 @@
                 else:
                     neighbors +=  "-1 "
-                
+        if triangles_attributes[index] == ['']:
+            att = ""
+        else:
+            att = str(triangles_attributes[index])
         fd.write(str(index) + " "
                  + str(tri[0]) + " " 
@@ -327 +350 @@
                  + str(tri[2]) + " " 
                  + neighbors + " "
-                 + str(triangles_attributes[index]) + "\n")
+                 + att + "\n")
             
     #One line:  <# of segments> 

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -489 +489 @@
         z1 = interp.interpolate(f)
         assert allclose(z, z1)
+
+    def test_fit_to_mesh_file(self):
+        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
+             export_trianglulation_file    
+        import tempfile
+        import os
+        
+        # create a .tsh file, no user outline
+        mesh_dic = {}
+        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
+                                          [0.0, 5.0],
+                                          [5.0, 0.0]]
+        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
+        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
+        mesh_dic['generatedtriangleattributelist'] = [['']]
+        mesh_dic['generatedpointattributelist'] = [[], [], []]
+        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
+        mesh_dic['generatedsegmentmarkerlist'] = ['external',
+                                                  'external',
+                                                  'external']       
+        mesh_file = tempfile.mktemp(".tsh")
+        export_trianglulation_file(mesh_file,mesh_dic)
+        
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("# demo \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
+        fd.close()
+        
+        mesh_output_file = "new_trianlge.tsh"
+        fit_to_mesh_file(mesh_file,
+                         point_file,
+                         mesh_output_file,
+                         alpha = 0.0)  
+        # load in the .tsh file we jusdt wrote
+        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
+
+        assert allclose(mesh_dic['generatedpointattributelist'],
+                        [[0.0, 0.0], [5.0, 10.0], [5.0,10.0]])
+        
+        #clean up
+        os.remove(mesh_file)
+        os.remove(point_file)
         
 #-------------------------------------------------------------