Changeset 4780

Timestamp:

Nov 5, 2007, 1:11:52 PM (17 years ago)

Author:

duncan

Message:

Bug fix

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (3 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (1 diff)
• shallow_water/data_manager.py (modified) (2 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (4 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -179 +179 @@
         self.build_tagged_elements_dictionary(tagged_elements)
 
+        # Build a list of vertices that are not connected to any triangles
+        self.lone_vertices = []
+        #Check that all vertices have been registered
+        for node, count in enumerate(self.number_of_triangles_per_node):    
+            #msg = 'Node %d does not belong to an element.' %node
+            #assert count > 0, msg
+            if count == 0:
+                self.lone_vertices.append(node)
+                
         #Update boundary indices FIXME: OBSOLETE
         #self.build_boundary_structure()
@@ -667 +676 @@
                 assert x < epsilon, msg
 
-        self.lone_vertices = []
-        #Check that all vertices have been registered
-        for node, count in enumerate(self.number_of_triangles_per_node):
-        
-            #msg = 'Node %d does not belong to an element.' %node
-            #assert count > 0, msg
-            if count == 0:
-                self.lone_vertices.append(node)
 
 
@@ -762 +763 @@
         """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

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -693 +693 @@
                                         verbose = False,
                                         use_cache = False):
-
         #FIXME: Use this function for the time being. Later move code in here
 

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4774 +4774 @@
     # this could be a problem. I would prefer taking them from
     # the instantiation of Domain.
+    #
+    # (DSG) There is not always a Domain instance when Write_sww is used.
+    # Check to see if this is the same level of hardwiring as is in
+    # shallow water doamain.
+    
     sww_quantities = Domain.conserved_quantities
 
@@ -4876 +4881 @@
         outfile.createVariable('volumes', Int, ('number_of_volumes',
                                                 'number_of_vertices'))
-        
+
         outfile.createVariable('time', precision,
                                ('number_of_timesteps',))

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -3 +3 @@
 import unittest, os
 from math import sqrt, pi
+import tempfile
 
 from anuga.config import g, epsilon
@@ -8 +9 @@
 from anuga.utilities.numerical_tools import mean
 from anuga.utilities.polygon import is_inside_polygon
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+from anuga.geospatial_data.geospatial_data import Geospatial_data
 
 from shallow_water_domain import *
@@ -13 +17 @@
 # Get gateway to C implementation of flux function for direct testing
 from shallow_water_ext import flux_function_central as flux_function
+
+# For test_fitting_using_shallow_water_domain example
+def linear_function(point):
+    point = array(point)
+    return point[:,0]+point[:,1]
 
 class Weir:
@@ -5185 +5194 @@
 
         #-------------------------------------------------------------
+
+    def test_get_lone_vertices(self):
+        
+        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 = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        boundary = { (0, 0): 'Third',
+                     (0, 2): 'First',
+                     (2, 0): 'Second',
+                     (2, 1): 'Second',
+                     (3, 1): 'Second',
+                     (3, 2): 'Third'}
+
+
+        domain = Domain(points, vertices, boundary)
+        #domain.check_integrity()
+        domain.get_lone_vertices()
+
+        
+    def test_fitting_using_shallow_water_domain(self):
+        
+        #Mesh in zone 56 (absolute coords)
+
+        x0 = 314036.58727982
+        y0 = 6224951.2960092
+
+        a = [x0+0.0, y0+0.0]
+        b = [x0+0.0, y0+2.0]
+        c = [x0+2.0, y0+0.0]
+        d = [x0+0.0, y0+4.0]
+        e = [x0+2.0, y0+2.0]
+        f = [x0+4.0, y0+0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        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))
+        mesh4.check_integrity()
+        quantity = Quantity(mesh4)
+
+        #Get (enough) datapoints (relative to georef)
+        data_points_rel = [[ 0.66666667, 0.66666667],
+                       [ 1.33333333, 1.33333333],
+                       [ 2.66666667, 0.66666667],
+                       [ 0.66666667, 2.66666667],
+                       [ 0.0, 1.0],
+                       [ 0.0, 3.0],
+                       [ 1.0, 0.0],
+                       [ 1.0, 1.0],
+                       [ 1.0, 2.0],
+                       [ 1.0, 3.0],
+                       [ 2.0, 1.0],
+                       [ 3.0, 0.0],
+                       [ 3.0, 1.0]]
+
+        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'
+        
+        #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())
+
+        assert allclose(quantity.vertex_values.flat, answer)
+
+
+        #Check that values can be set from file using default attribute
+        quantity.set_values(filename = ptsfile, alpha = 0)
+        assert allclose(quantity.vertex_values.flat, answer)
+
+        #Cleanup
+        import os
+        os.remove(ptsfile)
+
+
         
 if __name__ == "__main__":
 
     suite = unittest.makeSuite(Test_Shallow_Water,'test')
-    #suite = unittest.makeSuite(Test_Shallow_Water,'test_extrema')    
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')