Changeset 9650 for trunk/anuga_core/anuga/fit_interpolate/fitsmooth.c

Timestamp:

Feb 9, 2015, 10:13:02 PM (10 years ago)

Author:

steve

Message:

Reverted back to a working verions of fitsmooth.c

File:

• trunk/anuga_core/anuga/fit_interpolate/fitsmooth.c (modified) (27 diffs)

trunk/anuga_core/anuga/fit_interpolate/fitsmooth.c

@@ -152 +152 @@
     int k,k6;
     double x0,y0,x1,y1,x2,y2;
-    triangle *T;
 
     // set up quad tree and allocate memory
@@ -168 +167 @@
         x2 = vertex_coordinates[k6 + 4];
         y2 = vertex_coordinates[k6 + 5];
-        T = new_triangle(k,x0,y0,x1,y1,x2,y2);
+        triangle * T = new_triangle(k,x0,y0,x1,y1,x2,y2);
         quad_tree_insert_triangle(tree,T);
     }
@@ -188 +187 @@
               {
 
+
+
     int k;
+
+
+
     int i,w;
-    edge_key_t key;
-    double x;
-    double y;
-    triangle * T;
-    double * sigma;
-    int js[3];
-
-
     for(w=0;w<zdims;w++){
         for(i=0;i<N;i++){
@@ -204 +200 @@
     }
 
+    edge_key_t key;
+
+
+
 
     #pragma omp parallel for private(k,i,key,w)
@@ -209 +209 @@
 
 
-        x = point_coordinates[2*k];
-        y = point_coordinates[2*k+1];
-        T = search(quadtree,x,y);
+        double x = point_coordinates[2*k];
+        double y = point_coordinates[2*k+1];
+        triangle * T = search(quadtree,x,y);
 
         if(T!=NULL){
-            sigma = calculate_sigma(T,x,y);
-
+            double * sigma = calculate_sigma(T,x,y);
+            int js[3];
             for(i=0;i<3;i++){
                 js[i]=triangles[3*(T->index)+i];
@@ -251 +251 @@
                              double* Atz2,
                              int n, int zdim){
+
+    add_sparse_dok(dok_AtA1,1,dok_AtA2,1);
+
     int i;
-
-    add_sparse_dok(dok_AtA1,1,dok_AtA2,1);
-
     for(i=0;i<n*zdim;i++){
         Atz1[i]+=Atz2[i];
@@ -369 +369 @@
     PyArrayObject *vertex_coordinates;
     PyArrayObject *extents;
-    int n;
 
 
@@ -387 +386 @@
     CHECK_C_CONTIG(extents);
 
-    n = triangles->dimensions[0];
+    int n = triangles->dimensions[0];
 
     #ifdef PYVERSION273
@@ -417 +416 @@
     PyArrayObject *areas;
     PyArrayObject *vertex_coordinates;
-    sparse_dok * smoothing_mat; // Should be an input argument?
     int err;
-    int n;
 
     // Convert Python arguments to C
@@ -436 +433 @@
     CHECK_C_CONTIG(vertex_coordinates);
 
-    n = triangles->dimensions[0];
-
+    int n = triangles->dimensions[0];
+
+
+    sparse_dok * smoothing_mat; // Should be an input argument?
     smoothing_mat = make_dok();
 
@@ -478 +477 @@
     PyArrayObject *point_coordinates;
     PyArrayObject *z;
-    sparse_dok * dok_AtA; // Should be an input argument?
     int N; // Number of triangles
     int err;
     int npts;
     int zdims;
-    int i;
     PyObject *tree;
-    quad_tree *quadtree;
-    double ** Atz;
 
     // Convert Python arguments to C
@@ -507 +502 @@
 
     #ifdef PYVERSION273
-    quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
-    #else
-    quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
-    #endif
-
-
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    #else
+    quad_tree * quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
+    #endif
+
+    sparse_dok * dok_AtA; // Should be an input argument?
     dok_AtA = make_dok();
-    Atz = malloc(sizeof(double*)*zdims);
+    double ** Atz = malloc(sizeof(double*)*zdims);
+    int i;
     for(i=0;i<zdims;i++){
         Atz[i] = malloc(sizeof(double)*N);
@@ -566 +562 @@
     PyObject *tree;
     PyArrayObject *point;
-    quad_tree *quadtree;
-    double xp,yp;
-    int *ipointd;
-    double *pointd;
-    triangle *T;
-    PyObject *sigmalist;
-    long found;
-    long index;
-    double *sigma;
-    double ssigma[3];
-    PyObject *retlist;
-
 
     // Convert Python arguments to C
@@ -590 +574 @@
 
     #ifdef PYVERSION273
-    quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
-    #else
-    quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
-    #endif
-
-
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    #else
+    quad_tree * quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
+    #endif
+
+    double xp,yp;
     if(PyArray_TYPE(point)==7){
-        ipointd = (int*)point->data;
-        xp = (double)(ipointd[0]);
-        yp = (double)(ipointd[1]);
+        int *pointd = (int*)point->data;
+        xp = (double)(pointd[0]);
+        yp = (double)(pointd[1]);
     }
     else{
-        pointd = (double*)point->data;
+        double *pointd = (double*)point->data;
         xp = pointd[0];
         yp = pointd[1];
     }
 
-    T = search(quadtree,xp,yp);
-    sigmalist = PyList_New(3);
-
+    triangle * T = search(quadtree,xp,yp);
+    PyObject * sigmalist = PyList_New(3);
+    long found;
+    long index;
     
     if(T!=NULL){
-            sigma = calculate_sigma(T,xp,yp);
+            double * sigma = calculate_sigma(T,xp,yp);
             sigmalist = c_double_array_to_list(sigma,3);
             free(sigma);
@@ -618 +603 @@
             index = (long)T->index;
     }else{
-
-            ssigma[0]=ssigma[1]=ssigma[2]=-1;
-            sigmalist = c_double_array_to_list(ssigma,3);
+            double sigma[3];
+            sigma[0]=sigma[1]=sigma[2]=-1;
+            sigmalist = c_double_array_to_list(sigma,3);
             index = -10;
             found = 0;
     }
-    retlist = PyList_New(3);
+    PyObject * retlist = PyList_New(3);
     PyList_SET_ITEM(retlist,0,PyInt_FromLong(found));
     PyList_SET_ITEM(retlist,1,sigmalist);
@@ -641 +626 @@
     // Setting up variables to parse input
     PyObject *tree; // capsule to hold quad_tree pointer
-    quad_tree *quadtree;
 
     // Convert Python arguments to C
@@ -653 +637 @@
     // Extract quad_tree pointer from capsule
     #ifdef PYVERSION273
-    quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
-    #else
-    quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    #else
+    quad_tree * quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
     #endif
     
@@ -669 +653 @@
     // Setting up variables to parse input
     PyObject *tree; // capsule to hold quad_tree pointer
-    quad_tree *quadtree;
 
     // Convert Python arguments to C
@@ -681 +664 @@
     // Extract quad_tree pointer from capsule
     #ifdef PYVERSION273
-    quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
-    #else
-    quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    #else
+    quad_tree * quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
     #endif
     
@@ -707 +690 @@
     PyArrayObject *Atz1, *Atz2;
     int n,zdim; // number of nodes
-    sparse_dok *dok_AtA1;
-    sparse_dok *dok_AtA2;
 
     // Convert Python arguments to C
@@ -721 +702 @@
     // Get pointers to sparse_dok objects from capsules
     #ifdef PYVERSION273
-    dok_AtA1 = (sparse_dok*) PyCapsule_GetPointer(AtA_cap1,"sparse dok");
-    dok_AtA2 = (sparse_dok*) PyCapsule_GetPointer(AtA_cap2,"sparse dok");
-    #else
-    dok_AtA1 = (sparse_dok*) PyCObject_AsVoidPtr(AtA_cap1);
-    dok_AtA2 = (sparse_dok*) PyCObject_AsVoidPtr(AtA_cap2);
+    sparse_dok * dok_AtA1 = (sparse_dok*) PyCapsule_GetPointer(AtA_cap1,"sparse dok");
+    sparse_dok * dok_AtA2 = (sparse_dok*) PyCapsule_GetPointer(AtA_cap2,"sparse dok");
+    #else
+    sparse_dok * dok_AtA1 = (sparse_dok*) PyCObject_AsVoidPtr(AtA_cap1);
+    sparse_dok * dok_AtA2 = (sparse_dok*) PyCObject_AsVoidPtr(AtA_cap2);
     #endif
 
@@ -750 +731 @@
     PyObject *D_cap; // capsule object holding sparse_dok pointer
     int n; // number of matrix columns/rows
-    sparse_dok * D_mat;
-    PyObject *ret_D;
-    int i,j;
-    edge_key_t key;
-    edge_t *s;
 
     // Convert Python arguments to C
@@ -766 +742 @@
     // Get pointer to spars_dok struct
     #ifdef PYVERSION273
-    D_mat = (sparse_dok*) PyCapsule_GetPointer(D_cap,"sparse dok");
-    #else
-    D_mat = (sparse_dok*) PyCObject_AsVoidPtr(D_cap);
+    sparse_dok * D_mat = (sparse_dok*) PyCapsule_GetPointer(D_cap,"sparse dok");
+    #else
+    sparse_dok * D_mat = (sparse_dok*) PyCObject_AsVoidPtr(D_cap);
     #endif
 
@@ -780 +756 @@
 
     // Build new python list containing the full matrix to return
-    ret_D = PyList_New(n);
-
+    PyObject *ret_D = PyList_New(n);
+    int i,j;
+    edge_key_t key;
+    edge_t *s;
     for(i=0;i<n;i++)
     {
@@ -816 +794 @@
     PyObject *smoothing_mat_cap; // Capsule storing D pointer (sparse_dok struct)
     PyObject *AtA_cap; // Capsule storing AtA pointer (sparse_dok struct)
-    sparse_dok *smoothing_mat;
-    sparse_dok *dok_AtA;
-    sparse_csr *B;
-    PyObject *data;
-    PyObject *colind;
-    PyObject *row_ptr;
-    PyObject *lst;
 
     // Convert Python arguments to C
@@ -834 +805 @@
     // Extract pointers to c structs from capsules
     #ifdef PYVERSION273
-    smoothing_mat = (sparse_dok*) PyCapsule_GetPointer(smoothing_mat_cap,"sparse dok");
-    dok_AtA = (sparse_dok*) PyCapsule_GetPointer(AtA_cap,"sparse dok");
-    #else
-    smoothing_mat = (sparse_dok*) PyCObject_AsVoidPtr(smoothing_mat_cap);
-    dok_AtA = (sparse_dok*) PyCObject_AsVoidPtr(AtA_cap);
+    sparse_dok * smoothing_mat = (sparse_dok*) PyCapsule_GetPointer(smoothing_mat_cap,"sparse dok");
+    sparse_dok * dok_AtA = (sparse_dok*) PyCapsule_GetPointer(AtA_cap,"sparse dok");
+    #else
+    sparse_dok * smoothing_mat = (sparse_dok*) PyCObject_AsVoidPtr(smoothing_mat_cap);
+    sparse_dok * dok_AtA = (sparse_dok*) PyCObject_AsVoidPtr(AtA_cap);
     #endif
 
@@ -845 +816 @@
     
     // Create sparse_csr matrix and convert result to this format
+    sparse_csr * B;
     B = make_csr();
     convert_to_csr_ptr(B,smoothing_mat);
     
     // Extract the sparse_csr data to be returned as python lists
-    data = c_double_array_to_list(B->data,
+    PyObject *data = c_double_array_to_list(B->data,
                                 B->num_entries);
-    colind = c_int_array_to_list(B->colind,
+    PyObject *colind = c_int_array_to_list(B->colind,
                                 B->num_entries);
-    row_ptr = c_int_array_to_list(B->row_ptr,
+    PyObject *row_ptr = c_int_array_to_list(B->row_ptr,
                                 B->num_rows);
 
     // Build python list for return
-    lst = PyList_New(3);
+    PyObject *lst = PyList_New(3);
     PyList_SET_ITEM(lst, 0, data);
     PyList_SET_ITEM(lst, 1, colind);