Changeset 265

Timestamp:

Sep 1, 2004, 2:54:04 PM (20 years ago)

Author:

ole

Message:

Interpoalte_from_vertices_to_edges - implemented in C

Location:

inundation/ga/storm_surge/pyvolution

Files:

• quantity.py (modified) (5 diffs)
• quantity_ext.c (modified) (4 diffs)
• show_balanced_limiters.py (modified) (1 diff)
• test_quantity.py (modified) (2 diffs)
• util_ext.h (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -31 +31 @@
             self.vertex_values = zeros((N, 3), Float)
         else:    
-            self.vertex_values = array(vertex_values)
+            self.vertex_values = array(vertex_values).astype(Float)
 
             N, V = self.vertex_values.shape
@@ -72 +72 @@
 
     def interpolate_from_vertices_to_edges(self):
-        for k in range(self.vertex_values.shape[0]):
-            q0 = self.vertex_values[k, 0]
-            q1 = self.vertex_values[k, 1]
-            q2 = self.vertex_values[k, 2]
-            
-            self.edge_values[k, 0] = 0.5*(q1+q2)
-            self.edge_values[k, 1] = 0.5*(q0+q2)  
-            self.edge_values[k, 2] = 0.5*(q0+q1)
-
+        #Call correct module function
+        #(either from this module or C-extension)
+        interpolate_from_vertices_to_edges(self)
             
             
@@ -285 +279 @@
 
 
-def extrapolate_second_order(self):        
+
+def interpolate_from_vertices_to_edges(quantity):
+    """Compute edge values from vertex values using linear interpolation
+    """
+
+    for k in range(quantity.vertex_values.shape[0]):
+        q0 = quantity.vertex_values[k, 0]
+        q1 = quantity.vertex_values[k, 1]
+        q2 = quantity.vertex_values[k, 2]
+            
+        quantity.edge_values[k, 0] = 0.5*(q1+q2)
+        quantity.edge_values[k, 1] = 0.5*(q0+q2)  
+        quantity.edge_values[k, 2] = 0.5*(q0+q1)
+
+
+
+def extrapolate_second_order(quantity):        
     """Extrapolate conserved quantities from centroid to
     vertices for each volume using
@@ -291 +301 @@
     """
         
-    a, b = self.compute_gradients()
-
-    X = self.domain.get_vertex_coordinates()
-    qc = self.centroid_values
-    qv = self.vertex_values
+    a, b = quantity.compute_gradients()
+
+    X = quantity.domain.get_vertex_coordinates()
+    qc = quantity.centroid_values
+    qv = quantity.vertex_values
     
     #Check each triangle
-    for k in range(self.domain.number_of_elements):
+    for k in range(quantity.domain.number_of_elements):
         #Centroid coordinates            
-        x, y = self.domain.centroids[k]
+        x, y = quantity.domain.centroids[k]
         
         #vertex coordinates
@@ -452 +462 @@
         
     from quantity_ext import limit, compute_gradients,\
-    extrapolate_second_order        
+    extrapolate_second_order, interpolate_from_vertices_to_edges 
+

inundation/ga/storm_surge/pyvolution/quantity_ext.c

@@ -127 +127 @@
     vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
     
-    /*
-    printf("V %f, %f, %f\n", vertex_values[k3+0], 
-           vertex_values[k3+1], 
-           vertex_values[k3+2]); 
-    */
-           
-    //printf("k=%d, a=%f, b=%f\n", k, a[k], b[k]); 
-    
-           
-    //printf("C %f, %f, %f\n", centroid_values[k], 
-    //     centroid_values[k], 
-    //     centroid_values[k]);            
-           
-    //printf("X %f, %f, %f, %f, %f, %f\n", x0, y0, x1, y1, x2, y2);
   }
   return 0;
@@ -146 +132 @@
 
 
+
+
+int _interpolate(int N,
+                 double* vertex_values,
+                 double* edge_values) {
+  
+  int k, k3;
+  double q0, q1, q2;
+
+                     
+  for (k=0; k<N; k++) {
+    k3 = 3*k;    
+    
+    q0 = vertex_values[k3 + 0];
+    q1 = vertex_values[k3 + 1];
+    q2 = vertex_values[k3 + 2];
+            
+    //printf("%f, %f, %f\n", q0, q1, q2);
+    edge_values[k3 + 0] = 0.5*(q1+q2);
+    edge_values[k3 + 1] = 0.5*(q0+q2);  
+    edge_values[k3 + 2] = 0.5*(q0+q1);
+  }  
+  return 0;
+}
+                     
   
 /////////////////////////////////////////////////
 // Gateways to Python
+
+
+PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
+  
+  PyObject *quantity;
+  PyArrayObject *vertex_values, *edge_values;
+    
+  int N, err;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "O", &quantity)) 
+    return NULL;
+
+  vertex_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "vertex_values");           
+  if (!vertex_values) return NULL;                
+
+  edge_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "edge_values");   
+  if (!edge_values) return NULL;          
+
+  N = vertex_values -> dimensions[0];  
+  
+  err = _interpolate(N,
+                     (double*) vertex_values -> data,
+                     (double*) edge_values -> data);
+                     
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
+    return NULL;
+  }                  
+  
+  //Release and return               
+  Py_DECREF(vertex_values);    
+  Py_DECREF(edge_values); 
+
+  return Py_BuildValue("");
+}
+
+
 PyObject *compute_gradients(PyObject *self, PyObject *args) {
   
@@ -272 +323 @@
     PyObject_GetAttrString(quantity, "vertex_values");           
   if (!vertex_values) return NULL;                
+
+  
+  /*
+  printf("In extrapolate C routine\n");  
+  printf("d0=%d, d1=%d\n", 
+         vertex_values -> dimensions[0],
+         vertex_values -> dimensions[1]);
+  */
+  
+  vertex_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "vertex_values");           
+  if (!vertex_values) return NULL;                
   
   N = centroid_values -> dimensions[0];
@@ -411 +474 @@
   {"extrapolate_second_order", extrapolate_second_order, 
    METH_VARARGS, "Print out"},        
+  {"interpolate_from_vertices_to_edges", 
+   interpolate_from_vertices_to_edges,
+   METH_VARARGS, "Print out"},           
   {NULL, NULL, 0, NULL}   /* sentinel */
 };

inundation/ga/storm_surge/pyvolution/show_balanced_limiters.py

@@ -80 +80 @@
 
 #Evolve
-for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
+for t in domain.evolve(yieldstep = 0.1, finaltime = 50):
     domain.write_time()
     #print domain.quantities['level'].centroid_values[:6]

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -77 +77 @@
 
     def test_interpolation2(self):
-        quantity = Quantity(self.mesh4,
+        quantity = Conserved_quantity(self.mesh4,
                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
 
+
+        quantity.extrapolate_second_order()
+
+        assert allclose(quantity.vertex_values, [[2., 2., 2.],
+                                               [3.+2./3, 6.+2./3, 4.+2./3],
+                                               [7.5, 0.5, 1.],
+                                               [-5, -2.5, 7.5]])
+        
+        #print quantity.edge_values        
         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
                                                [5., 5., 5.],
@@ -236 +245 @@
             assert allclose(b[i], 1.0)
 
-
         quantity.extrapolate_second_order()
-
+        
+        #print quantity.vertex_values
         assert allclose(quantity.vertex_values[1,0], 2.0)
         assert allclose(quantity.vertex_values[1,1], 6.0)        

inundation/ga/storm_surge/pyvolution/util_ext.h

@@ -10 +10 @@
 // Ole Nielsen, GA 2004
         
-
+#include "Python.h"     
+#include "Numeric/arrayobject.h"
 #include "math.h"
 
@@ -49 +50 @@
 }
 
+
+void print_numeric_array(PyArrayObject *x) {  
+  int i, j;
+  for (i=0; i<x->dimensions[0]; i++) {  
+    for (j=0; j<x->dimensions[1]; j++) {
+      printf("%f ", *(double*) (x->data + i*x->strides[0] + j*x->strides[1]));
+    }
+    printf("\n");  
+  }
+  printf("\n");    
+}
+
+void print_numeric_vector(PyArrayObject *x) {  
+  int i;
+  for (i=0; i<x->dimensions[0]; i++) {
+    printf("%f ", *(double*) (x->data + i*x->strides[0]));  
+  }
+  printf("\n");  
+}