Changeset 3012

Timestamp:

May 30, 2006, 11:00:57 AM (19 years ago)

Author:

duncan

Message:

started fixing the black screen of death at the fit.py end.

Location:

inundation

Files:

• examples/bedslope_bsod_test.py (modified) (3 diffs)
• fit_interpolate/fit.py (modified) (4 diffs)
• fit_interpolate/test_fit.py (modified) (1 diff)
• pyvolution/mesh.py (modified) (2 diffs)
• pyvolution/test_mesh.py (modified) (1 diff)

inundation/examples/bedslope_bsod_test.py

@@ -30 +30 @@
 c3 = [331000, 6267000]
 bound_poly = [c0, c1, c2, c3]
-# duplicate point, this should fail, i.e. conjugate gradient doesn't converge
-#bound_poly = [c0, c1, c2, c3, c2]
 
 # set up internal region
@@ -40 +38 @@
 interior_polygon = [p0, p1, p2, p3]
 # duplicate point, this should fail, i.e. conjugate gradient doesn't converge
-#interior_polygon = [p0, p1, p2, p3, p0] 
+interior_polygon = [p0, p1, p2, p3, p0] 
 
-interior_regions = [[interior_polygon, 25000]]
+interior_region_area = 2500000000 #25000
+interior_regions = [[interior_polygon, interior_region_area]]
 
-meshname = 'bedslope_bsod.msh'
+meshname = 'bedslope_bsod.tsh'
 
+exterior_region_area = 10000000000000000 #100000
 # create a mesh
 _ = cache(create_mesh_from_regions,
@@ -51 +51 @@
           {'boundary_tags': {'bottom': [0], 'right': [1],
                              'top': [2], 'left': [3]},
-           'maximum_triangle_area': 100000,
+           'maximum_triangle_area': exterior_region_area,
            'filename': meshname,           
            'interior_regions': interior_regions},

inundation/fit_interpolate/fit.py

@@ -21 +21 @@
 """
 
-from Numeric import zeros, Float, ArrayType
+from Numeric import zeros, Float, ArrayType,take 
 
 from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
@@ -73 +73 @@
               a mesh origin, since geospatial has its own mesh origin.
         """
-
         # Initialise variabels
-        #self._A_can_be_reused = False
-        #self._point_coordinates = None
 
         if alpha is None:
@@ -307 +304 @@
                                                 geo_reference=point_origin)
             self.build_fit_subset(point_coordinates, z, verbose)
-       
+
+        #FIXME: Move this down
         #Check sanity
         m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
@@ -320 +318 @@
 
         self._build_coefficient_matrix_B(verbose)
+        indexes = range(0,m)
+        loners = self.mesh.get_lone_vertices()
+        #print "loners",loners
+        loners.sort()
+        loners.reverse()
+        #print "loners",loners
+        for lone_index in loners:
+            indexes.pop(lone_index)
+        #print "fit self.Atz",take(self.Atz,indexes)
+        
         return conjugate_gradient(self.B, self.Atz, self.Atz,
                                   imax=2*len(self.Atz) )

inundation/fit_interpolate/test_fit.py

@@ -804 +804 @@
         os.remove(point_file)
 
+
+    def test_smooth_att_to_mesh_with_excess_verts(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        d = [1.0, 1.0]
+        e = [18.0, 1000.0]
+        
+        points = [a, b, c, d, e]
+        triangles = [ [1,0,2] ]   #bac
+
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 2.0]
+        d3 = [3.0,1.0]
+        d4 = [2.0,3.0]
+        d5 = [2.0,2.0]
+        d6 = [1.0,3.0]
+        data_coords = [d1, d2, d3, d4, d5, d6]
+        z = linear_function(data_coords)
+        #print "z",z 
+
+        interp = Fit(points, triangles, alpha=0.0)
+        f = interp.fit(data_coords, z)
+        answer = linear_function(points)
+
+        #  Removing the bad verts that we don't care about
+        f = f[0:3]
+        answer = answer[0:3]
+
+        #print "f\n",f
+        #print "answer\n",answer
+
+        assert allclose(f, answer)
+        # Note: This test is iffy. Fit can fail if lone verts are inputted.
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Fit,'test')
-    #suite = unittest.makeSuite(Test_Fit,'test_smoothing_and_interpolation')
+    #suite = unittest.makeSuite(Test_Fit,'test_smooth_att_to_mesh_with_excess_verts')
     #suite = unittest.makeSuite(Test_Fit,'test_smooth_attributes_to_mesh_one_point')
     runner = unittest.TextTestRunner(verbosity=1)

inundation/pyvolution/mesh.py

@@ -616 +616 @@
                 assert x < epsilon, msg
 
-
+        self.lone_vertices = []
         #Check that all vertices have been registered
         for v_id, v in enumerate(self.vertexlist):
 
-            msg = 'Vertex %s does not belong to an element.'
+            #msg = 'Vertex %s does not belong to an element.'
             #assert v is not None, msg
             if v is None:
-                print msg%v_id
+                #print msg%v_id
+                self.lone_vertices.append(v_id)
 
         #Check integrity of neighbour structure
@@ -673 +674 @@
         #See domain.set_boundary
 
+    def get_lone_vertices(self):
+        """Return a list of vertices that are not connected to any triangles.
+
+        Precondition
+        FIXME(DSG - DSG) Pull the code out of check integrity that builds this
+                         structure.
+        check_integrity has to have been called. 
+        """
+        return self.lone_vertices
 
     def get_centroid_coordinates(self):

inundation/pyvolution/test_mesh.py

@@ -912 +912 @@
 
 
+
+    def test_lone_vertices(self):
+        a = [2.0, 1.0]
+        b = [6.0, 2.0]
+        c = [1.0, 3.0]
+        d = [2.0, 4.0]
+
+        points = [a, b, c, d]
+        vertices = [[0,1,2]]
+
+        mesh = Mesh(points, vertices)
+        mesh.check_integrity()
+        loners = mesh.get_lone_vertices()
+        self.failUnless(loners==[3],
+                        'FAILED!')
+
+        
+        a = [2.0, 1.0]
+        b = [6.0, 2.0]
+        c = [1.0, 3.0]
+        d = [2.0, 4.0]
+
+        points = [d, a, b, c]
+        vertices = [[3,1,2]]
+
+        mesh = Mesh(points, vertices)
+        mesh.check_integrity()
+        loners = mesh.get_lone_vertices()
+        self.failUnless(loners==[0],
+                        'FAILED!') 
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon')
-    #suite = unittest.makeSuite(Test_Mesh,'test')
+    #suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon')
+    suite = unittest.makeSuite(Test_Mesh,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)