Changeset 262 for inundation/ga/storm_surge

Timestamp:

Sep 1, 2004, 5:26:59 AM (20 years ago)

Author:

ole

Message:

extrapolate 2 order - c implementation

Location:

inundation/ga/storm_surge/pyvolution

Files:

• mesh.py (modified) (1 diff)
• quantity.py (modified) (1 diff)
• quantity_ext.c (modified) (8 diffs)
• test_quantity.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -81 +81 @@
         #Allocate space for geometric quantities
         self.centroids = zeros((N, 2), Float)
+        #FIXME: Should be renamed to centroid_coordinates
+        
         self.areas = zeros(N, Float)
         self.radii = zeros(N, Float)

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -451 +451 @@
     #Replace python version with c implementations
         
-    from quantity_ext import limit, compute_gradients #, extrapolate_second_order        
+    from quantity_ext import limit, compute_gradients,\
+    extrapolate_second_order        

inundation/ga/storm_surge/pyvolution/quantity_ext.c

@@ -94 +94 @@
 
 
-
-
+int _extrapolate(int N,
+                 double* centroids, 
+                 double* centroid_values,                
+                 double* vertex_coordinates,
+                 double* vertex_values,
+                 double* a,
+                 double* b) {
+  
+  int k, k2, k3, k6;
+  double x, y, x0, y0, x1, y1, x2, y2;
+                     
+  for (k=0; k<N; k++) {
+    k6 = 6*k;
+    k3 = 3*k;    
+    k2 = 2*k;
+        
+    //Centroid coordinates            
+    x = centroids[k2]; y = centroids[k2+1]; 
+
+    //vertex coordinates
+    //x0, y0, x1, y1, x2, y2 = X[k,:]
+    x0 = vertex_coordinates[k6 + 0];
+    y0 = vertex_coordinates[k6 + 1];    
+    x1 = vertex_coordinates[k6 + 2];
+    y1 = vertex_coordinates[k6 + 3];        
+    x2 = vertex_coordinates[k6 + 4];
+    y2 = vertex_coordinates[k6 + 5];           
+
+    //Extrapolate
+    vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
+    vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
+    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;
+}
+
+
+  
 /////////////////////////////////////////////////
 // Gateways to Python
-
-
-
 PyObject *compute_gradients(PyObject *self, PyObject *args) {
   
-  PyObject *quantity, *domain;
+  PyObject *quantity, *domain, *R;
   PyArrayObject 
     *centroids,            //Coordinates at centroids
@@ -162 +210 @@
   }
                      
-  
-//Return values
-  return Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); 
+  //Release                  
+  Py_DECREF(centroids);    
+  Py_DECREF(centroid_values);    
+  Py_DECREF(number_of_boundaries); 
+  Py_DECREF(surrogate_neighbours);           
+         
+  //Build result, release and return
+  R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); 
+  Py_DECREF(a);      
+  Py_DECREF(b);        
+  return R;
 }
 
 
-/*
+
 PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
-
-
-
-PyObject *compute_gradients(PyObject *self, PyObject *args) {
   
   PyObject *quantity, *domain;
@@ -179 +231 @@
     *centroids,            //Coordinates at centroids
     *centroid_values,      //Values at centroids    
+    *vertex_coordinates,   //Coordinates at vertices            
+    *vertex_values,        //Values at vertices        
     *number_of_boundaries, //Number of boundaries for each triangle
     *surrogate_neighbours, //True neighbours or - if one missing - self
-    *a, *b;                //Return values
+    *a, *b;                //Gradients
    
-  int dimensions[1], N, err;
+  //int N, err;
+  int dimensions[1], N, err;  
+  //double *a, *b;  //Gradients
   
   // Convert Python arguments to C  
@@ -208 +264 @@
     PyObject_GetAttrString(domain, "number_of_boundaries");           
   if (!number_of_boundaries) return NULL;            
+  
+  vertex_coordinates = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "vertex_coordinates");           
+  if (!vertex_coordinates) return NULL;              
+  
+  vertex_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "vertex_values");           
+  if (!vertex_values) return NULL;                
   
   N = centroid_values -> dimensions[0];
@@ -218 +282 @@
   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);  
-        
+
+  //FIXME: Odd that I couldn't use normal arrays          
+  //Allocate space for return vectors a and b (don't DECREF) 
+  //a = (double*) malloc(N * sizeof(double));
+  //if (!a) return NULL;  
+  //b = (double*) malloc(N * sizeof(double));  
+  //if (!b) return NULL;
 
   
@@ -227 +297 @@
                            (int*) surrogate_neighbours -> data,
                            (double*) a -> data,
-                           (double*) b -> data);     
+                           (double*) b -> data);
                            
   if (err != 0) {
@@ -233 +303 @@
     return NULL;
   }
-
-
-  PyObject *quantity, *domain, *Tmp;
-  PyArrayObject 
-    *qv, //Conserved quantities at vertices 
-    *qc, //Conserved quantities at centroids
-    *neighbours;
-    
-  int k, i, n, N, k3;
-  double beta; //Safety factor
-  double *qmin, *qmax, qn;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "O", &quantity)) 
-    return NULL;
-
-  domain = PyObject_GetAttrString(quantity, "domain");     
-  if (!domain) 
-    return NULL;  
-    
-  neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
-
-  qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); 
-  qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
-  N = qc -> dimensions[0];
-*/
-  /*
-def extrapolate_second_order(self):        
-    """Extrapolate conserved quantities from centroid to
-    vertices for each volume using
-    second order scheme.
-    """
-        
-    a, b = self.compute_gradients()
-
-    V = self.domain.get_vertex_coordinates()
-    qc = self.centroid_values
-    qv = self.vertex_values
-    
-    #Check each triangle
-    for k in range(self.domain.number_of_elements):
-        #Centroid coordinates            
-        x, y = self.domain.centroids[k]
-        
-        #vertex coordinates
-        x0, y0, x1, y1, x2, y2 = V[k,:]
-        
-        #Extrapolate
-        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
-        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
-        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
-  */
-//}    
+  
+  err = _extrapolate(N, 
+                     (double*) centroids -> data, 
+                     (double*) centroid_values -> data,              
+                     (double*) vertex_coordinates -> data,
+                     (double*) vertex_values -> data,
+                     (double*) a -> data,
+                     (double*) b -> data);                   
+  //a, b);
+                     
+                     
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, 
+                    "Internal function _extrapolate failed");
+    return NULL;
+  }                  
+                
+  
+
+  //Release  
+  Py_DECREF(centroids);    
+  Py_DECREF(centroid_values);    
+  Py_DECREF(number_of_boundaries); 
+  Py_DECREF(surrogate_neighbours);
+  Py_DECREF(vertex_coordinates);                   
+  Py_DECREF(vertex_values);                 
+  Py_DECREF(a);                   
+  Py_DECREF(b);                     
+  //free(a);      
+  //free(b);        
+  
+  return Py_BuildValue("");    
+}    
 
 
@@ -359 +409 @@
   {"limit", limit, METH_VARARGS, "Print out"},    
   {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},      
+  {"extrapolate_second_order", extrapolate_second_order, 
+   METH_VARARGS, "Print out"},        
   {NULL, NULL, 0, NULL}   /* sentinel */
 };

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -194 +194 @@
         quantity.set_values([2.0, 4.0, 8.0, 2.0],
                             location = 'centroids')
-        
+
+
+
         #Gradients
         a, b = quantity.compute_gradients()
-
 
         #gradient bewteen t0 and t1 is undefined as det==0
@@ -209 +210 @@
 
         quantity.extrapolate_second_order()
-        
+
+        #print quantity.vertex_values
         assert allclose(quantity.vertex_values, [[2., 2.,  2.],
                                                  [0., 6.,  6.],
@@ -215 +217 @@
                                                  [0., 0.,  6.]])
 
-
-
+                
     def test_second_order_extrapolation2(self):
         quantity = Conserved_quantity(self.mesh4)