Changeset 275

Timestamp:

Sep 6, 2004, 5:49:21 PM (20 years ago)

Author:

ole

Message:

get_vertex_values for (smoothed) data output

Location:

inundation/ga/storm_surge/pyvolution

Files:

• mesh.py (modified) (3 diffs)
• quantity.py (modified) (4 diffs)
• test_quantity.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -60 +60 @@
         from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
 
-        #Only one of them
         self.vertices = array(vertices).astype(Int)
         self.coordinates = array(coordinates)
-
 
         #Input checks
@@ -270 +268 @@
             c = self.vertices[i, 2]            
 
-            #Register the vertices as keys mapping to 
-            #(triangle, edge) tuples associated with them
+            #Register the vertices v as lists of
+            #(triangle_id, vertex_id) tuples associated with them
             #This is used for smoothing
             for vertex_id, v in enumerate([a,b,c]):
@@ -526 +524 @@
         """
 
+        #FIXME: Perhaps they should be ordered as in obj files??
+        #See quantity.get_vertex_values
+        
         from Numeric import zeros, Float
 

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -54 +54 @@
         self.interpolate()
 
-
+    def __len__(self):
+        return self.centroid_values.shape[0]
+        
     def interpolate(self):
         """Compute interpolated values at edges and centroid
@@ -206 +208 @@
         TODO: be able to smooth individual fields
         NOTE:  This function does not have a test.
+        FIXME: NOT DONE
         """
 
@@ -276 +279 @@
 
         if smooth == True:
-            M = len(self.domain.vertexlist) #Number of unique vertices
-            
-            A = zeros((M,N), precision)
-            V = zeros((len(self),3), Int)            
-
-            #Rebuild triangle structure
-            for k, v in enumerate(self.vertexdict.keys()):            
-                for volume_id, vertex_id in self.vertexdict[v]:                
-                    V[volume_id, vertex_id] = k
-
-
+
+            N = len(self.domain.vertexlist)
+            A = zeros(N, precision)
+            V = self.domain.vertices
+            
             #Smoothing loop
-            for j, i in enumerate(indices):
-                #For each quantity
-                #
-                #j is an index into returned quantities (0 <= j < N)
-                #i the given index into internal quantities
-                #(e.g. 0: bed elevation or 1: friction)
-
-                for k, v in enumerate(self.vertexdict.keys()):
-                    #Go through all contributions to vertex v
-                    #k is an index into returned vertices (0 <= k < M)
-                    #v is an index into internal vertices
-                    
-                    L = len(self.vertexdict[v])
-
-                    ulist = []
-                    for volume_id, vertex_id in self.vertexdict[v]:
-                        #For all contributing triange, vertex pairs
-                        
-                        cmd = 'u = Volume.%s_vertex%d[%d,%d]'\
-                              %(value_array, vertex_id, volume_id, i)
-                        exec(cmd)
-                        ulist.append(u)
-
-                    A[k,j] = reduction(ulist)    
+            for k in range(N):
+                L = self.domain.vertexlist[k]
+                
+                #Go through all triangle, vertex pairs
+                #contributing to vertex k and register vertex value
+                contributions = []                
+                for volume_id, vertex_id in L:
+
+                    v = self.vertex_values[volume_id, vertex_id]
+                    contributions.append(v)
+
+                A[k] = reduction(contributions)    
 
 
             if xy is True:
-                X = zeros(M, precision)
-                Y = zeros(M, precision)
-
-                for k, v in enumerate(self.vertexdict.keys()):            
-                    X[k], Y[k] = Point.coordinates[v].astype(precision)
-
+                X = self.domain.coordinates[:,0]
+                Y = self.domain.coordinates[:,1]
                 return X, Y, A, V
             else:
@@ -338 +319 @@
             #Do vertex coordinates    
             if xy is True:                
-                V = self.domain.get_vertex_coordinates()
-                X = concatenate((V[:,0], V[:,2], V[:,4]), axis=0).\
-                    astype(precision)
-                Y = concatenate((V[:,1], V[:,3], V[:,5]), axis=0).\
-                    astype(precision)
+                C = self.domain.get_vertex_coordinates()
+
+                X = C[:,0:6:2].copy()
+                Y = C[:,1:6:2].copy()                
+
                 
-                return X, Y, A , V           
+                return X.flat, Y.flat, A, V           
             else:
                 return A, V

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -407 +407 @@
 
         
+
+
+    #Test smoothing    
+    def test_smoothing(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Transmissive_boundary
+        from Numeric import zeros, Float
+        from util import mean
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(2, 2)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order=2
+        domain.reduction = mean
+
+        
+        #Set some field values
+        domain.set_quantity('elevation', lambda x,y: x)        
+        domain.set_quantity('friction', 0.03)
+
+
+        ###################### 
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+        
+
+        ###################### 
+        #Initial condition - with jumps
+
+        bed = domain.quantities['elevation'].vertex_values
+        level = zeros(bed.shape, Float)
+
+        h = 0.03
+        for i in range(level.shape[0]):
+            if i % 2 == 0:            
+                level[i,:] = bed[i,:] + h
+            else:
+                level[i,:] = bed[i,:]
+                
+        domain.set_quantity('level', level)
+
+        level = domain.quantities['level']
+
+        #Get smoothed level
+        A, V = level.get_vertex_values(xy=False, smooth=True)
+        Q = level.vertex_values
+
+        #First four points 
+        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)  
+        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
+        assert allclose(A[2], Q[3,0])
+        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
+
+        #Center point
+        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
+                               Q[5,0] + Q[6,2] + Q[7,1])/6)
+                                 
+
+        #Check V
+        assert allclose(V[0,:], [3,4,0])
+        assert allclose(V[1,:], [1,0,4])
+        assert allclose(V[2,:], [4,5,1])                        
+        assert allclose(V[3,:], [2,1,5])
+        assert allclose(V[4,:], [6,7,3])        
+        assert allclose(V[5,:], [4,3,7])
+        assert allclose(V[6,:], [7,8,4])
+        assert allclose(V[7,:], [5,4,8])                        
+
+        #Get smoothed level with XY
+        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=True)
+
+        assert allclose(A, A1)
+        assert allclose(V, V1)        
+
+        #Check XY
+        assert allclose(X[4], 0.5)
+        assert allclose(Y[4], 0.5)
+
+        assert allclose(X[7], 1.0)
+        assert allclose(Y[7], 0.5)                 
+
+
+
+
+    def test_vertex_values_no_smoothing(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Transmissive_boundary
+        from Numeric import zeros, Float
+        from util import mean
+
+        
+        #Create basic mesh
+        points, vertices, boundary = rectangular(2, 2)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order=2
+        domain.reduction = mean
+
+        
+        #Set some field values
+        domain.set_quantity('elevation', lambda x,y: x)        
+        domain.set_quantity('friction', 0.03)
+
+        #Get level
+        level = domain.quantities['level']        
+        A, V = level.get_vertex_values(xy=False, smooth=False)
+        Q = level.vertex_values
+
+        for k in range(8):
+            assert allclose(A[k], Q[k])
+
+        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=False)
+
+        
+        assert allclose(A, A1)
+        assert allclose(V, V1)        
+
+        #Check XY
+        assert allclose(X[1], 0.5)
+        assert allclose(Y[1], 0.5)
+        
+        assert allclose(X[4], 0.0)
+        assert allclose(Y[4], 0.0)
+
+        assert allclose(X[12], 1.0)
+        assert allclose(Y[12], 0.0)                
+        
+        
+                                 
+
         
 #-------------------------------------------------------------