Changeset 260 for inundation/ga

Timestamp:

Sep 1, 2004, 3:31:13 AM (20 years ago)

Author:

ole

Message:

Cleaned up compute_gradients. Prepared for c extension

Location:

inundation/ga/storm_surge/pyvolution

Files:

• quantity.py (modified) (5 diffs)
• quantity_ext.c (modified) (1 diff)
• shallow_water.py (modified) (1 diff)
• shallow_water_ext.c (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -255 +255 @@
 
     def compute_gradients(self):
-        """Compute gradients of triangle surfaces defined by centroids of
-        neighbouring volumes.
-        If one face is on the boundary, use own centroid as neighbour centroid.
-        If two or more are on the boundary, fall back to first order scheme.
-        
-        Also return minimum and maximum of conserved quantities 
-        """
-
-        
-        from Numeric import array, zeros, Float
-        from util import gradient
-
-        N = self.centroid_values.shape[0]
-
-        a = zeros(N, Float)
-        b = zeros(N, Float)
-        
-        for k in range(N):
-        
-            number_of_boundaries = self.domain.number_of_boundaries[k]
-
-            if number_of_boundaries == 3:                 
-                #We have zero neighbouring volumes -
-                #Fall back to first order scheme
-                pass
-            
-            elif number_of_boundaries == 2:
-                #Special case where we have only one neighbouring volume.
-
-                #Find index of the one neighbour
-                for k0 in self.domain.neighbours[k,:]:
-                    if k0 >= 0:
-                        break
-
-                assert k0 != k
-                assert k0 >= 0
-
-                k1 = k  #Self
-
-                #Get data
-                q0 = self.centroid_values[k0]
-                q1 = self.centroid_values[k1]
-
-                x0, y0 = self.domain.centroids[k0] #V0 centroid
-                x1, y1 = self.domain.centroids[k1] #V1 centroid        
-
-                #Gradient
-                det = x0*y1 - x1*y0
-                if det != 0.0:
-                    a[k] = (y1*q0 - y0*q1)/det
-                    b[k] = (x0*q1 - x1*q0)/det
-
-            else:
-                #One or zero missing neighbours
-                #In case of one boundary - own centroid
-                #has been inserted as a surrogate for the one
-                #missing neighbour in neighbour_surrogates
-
-                #Get data        
-                k0 = self.domain.surrogate_neighbours[k,0]
-                k1 = self.domain.surrogate_neighbours[k,1]
-                k2 = self.domain.surrogate_neighbours[k,2]
-
-                q0 = self.centroid_values[k0]
-                q1 = self.centroid_values[k1]
-                q2 = self.centroid_values[k2]                    
-
-                x0, y0 = self.domain.centroids[k0] #V0 centroid
-                x1, y1 = self.domain.centroids[k1] #V1 centroid
-                x2, y2 = self.domain.centroids[k2] #V2 centroid
-
-                #Gradient
-                a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-        
-        return a, b
-        
+        #Call correct module function
+        #(either from this module or C-extension)
+        return compute_gradients(self)
+            
 
     def limit(self):
@@ -365 +293 @@
     a, b = self.compute_gradients()
 
-    V = self.domain.get_vertex_coordinates()
+    X = self.domain.get_vertex_coordinates()
     qc = self.centroid_values
     qv = self.vertex_values
@@ -375 +303 @@
         
         #vertex coordinates
-        x0, y0, x1, y1, x2, y2 = V[k,:]
+        x0, y0, x1, y1, x2, y2 = X[k,:]
         
         #Extrapolate
@@ -382 +310 @@
         qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
 
-            
+
+def compute_gradients(quantity):                    
+    """Compute gradients of triangle surfaces defined by centroids of
+    neighbouring volumes.
+    If one edge is on the boundary, use own centroid as neighbour centroid.
+    If two or more are on the boundary, fall back to first order scheme.
+    """
+
+    from Numeric import zeros, Float
+    from util import gradient
+    
+    centroids = quantity.domain.centroids
+    surrogate_neighbours = quantity.domain.surrogate_neighbours    
+    centroid_values = quantity.centroid_values    
+    number_of_boundaries = quantity.domain.number_of_boundaries     
+    
+    N = centroid_values.shape[0]
+
+    a = zeros(N, Float)
+    b = zeros(N, Float)
+    
+    for k in range(N):
+        if number_of_boundaries[k] < 2:
+            #Two or three true neighbours
+
+            #Get indices of neighbours (or self when used as surrogate)     
+            k0, k1, k2 = surrogate_neighbours[k,:]
+
+            q0 = centroid_values[k0]
+            q1 = centroid_values[k1]
+            q2 = centroid_values[k2]                    
+
+            x0, y0 = centroids[k0] #V0 centroid
+            x1, y1 = centroids[k1] #V1 centroid
+            x2, y2 = centroids[k2] #V2 centroid
+
+            #Gradient
+            a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+        
+        elif number_of_boundaries[k] == 2:
+            #One true neighbour
+
+            #Get index of the one neighbour
+            for k0 in surrogate_neighbours[k,:]:
+                if k0 != k: break
+            assert k0 != k
+
+            k1 = k  #self
+
+            #Get data
+            q0 = centroid_values[k0]
+            q1 = centroid_values[k1]
+
+            x0, y0 = centroids[k0] #V0 centroid
+            x1, y1 = centroids[k1] #V1 centroid        
+
+            #Gradient
+            det = x0*y1 - x1*y0
+            if det != 0.0:
+                a[k] = (y1*q0 - y0*q1)/det
+                b[k] = (x0*q1 - x1*q0)/det
+
+        else:
+            #No true neighbours -       
+            #Fall back to first order scheme
+            pass
+        
+    
+    return a, b
+        
+                    
 
 def limit(quantity):        
@@ -452 +450 @@
     #Replace python version with c implementations
         
-    from quantity_ext import limit #, extrapolate_second_order        
+    from quantity_ext import limit #, compute_gradients #, extrapolate_second_order        

inundation/ga/storm_surge/pyvolution/quantity_ext.c

@@ -20 +20 @@
 /////////////////////////////////////////////////
 // Gateways to Python
+
+
+
+
+
+/*
+
+def compute_gradients(quantity):                    
+    """Compute gradients of triangle surfaces defined by centroids of
+    neighbouring volumes.
+    If one edge is on the boundary, use own centroid as neighbour centroid.
+    If two or more are on the boundary, fall back to first order scheme.
+    """
+
+    from Numeric import zeros, Float
+    from util import gradient
+    
+    centroids = quantity.domain.centroids
+    surrogate_neighbours = quantity.domain.surrogate_neighbours    
+    centroid_values = quantity.centroid_values    
+    number_of_boundaries = quantity.domain.number_of_boundaries     
+    
+    N = centroid_values.shape[0]
+
+    a = zeros(N, Float)
+    b = zeros(N, Float)
+    
+    for k in range(N):
+        if number_of_boundaries[k] < 2:
+            #Two or three true neighbours
+
+            #Get indices of neighbours (or self when used as surrogate)     
+            k0, k1, k2 = surrogate_neighbours[k,:]
+
+            q0 = centroid_values[k0]
+            q1 = centroid_values[k1]
+            q2 = centroid_values[k2]                    
+
+            x0, y0 = centroids[k0] #V0 centroid
+            x1, y1 = centroids[k1] #V1 centroid
+            x2, y2 = centroids[k2] #V2 centroid
+
+            #Gradient
+            a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+        
+        elif number_of_boundaries[k] == 2:
+            #One true neighbour
+
+            #Get index of the one neighbour
+            for k0 in surrogate_neighbours[k,:]:
+                if k0 != k: break
+            assert k0 != k
+
+            k1 = k  #self
+
+            #Get data
+            q0 = centroid_values[k0]
+            q1 = centroid_values[k1]
+
+            x0, y0 = centroids[k0] #V0 centroid
+            x1, y1 = centroids[k1] #V1 centroid        
+
+            #Gradient
+            det = x0*y1 - x1*y0
+            if det != 0.0:
+                a[k] = (y1*q0 - y0*q1)/det
+                b[k] = (x0*q1 - x1*q0)/det
+
+        else:
+            #No true neighbours -       
+            #Fall back to first order scheme
+            pass
+        
+    
+    return a, b
+        
+*/
+
+
+
+
+
+
+
+
+
+
 
 /*

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -936 +936 @@
     #Replace python version with c implementations
     
-    pass
     from shallow_water_ext import rotate
     compute_fluxes = compute_fluxes_c

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -302 +302 @@
   _rotate((double *) r -> data, n1, n2);
 
-
+  //Release numeric arrays
   Py_DECREF(q);    
   Py_DECREF(normal);