Changeset 8726

Timestamp:

Mar 5, 2013, 2:40:28 PM (12 years ago)

Author:

steve

Message:

Moving the documentation to a higher level

Location:

trunk

Files:

• anuga_core/source/anuga/file/sww.py (modified) (4 diffs)
• anuga_core/source/anuga/operators/kinematic_viscosity_operator.py (modified) (2 diffs)
• anuga_core/source/anuga/operators/kinematic_viscosity_operator_ext.c (modified) (12 diffs)
• anuga_core/source/anuga/utilities/sww_merge.py (modified) (4 diffs)
• anuga_documentation (moved) (moved from trunk/anuga_core/documentation)

trunk/anuga_core/source/anuga/file/sww.py

@@ -300 +300 @@
                 #
                 # In this branch it is assumed that elevation
-                # is also available as a quantity            
+                # is also available as a quantity
+
+
+                # Smoothing for the get_vertex_values will be obtained
+                # from the smooth setting in domain
             
                 Q = domain.quantities['stage']
@@ -543 +547 @@
         # dimension definitions
         outfile.createDimension('number_of_volumes', number_of_volumes)
+        outfile.createDimension('number_of_triangle_vertices', number_of_points)
         outfile.createDimension('number_of_vertices', 3)
         outfile.createDimension('numbers_in_range', 2)
@@ -708 +713 @@
         # variable definitions
         outfile.createVariable('tri_l2g',  netcdf_int, ('number_of_volumes',))
-        outfile.createVariable('node_l2g', netcdf_int, ('number_of_points',))
+        outfile.createVariable('node_l2g', netcdf_int, ('number_of_triangle_vertices',))
         outfile.createVariable('tri_full_flag', netcdf_int, ('number_of_volumes',))
 
@@ -717 +722 @@
         outfile.variables['tri_l2g'][:] = tri_l2g.astype(num.int32)
 
-        #print node_l2g.shape
+        print node_l2g.shape
         #print node_l2g
-        #print outfile.variables['node_l2g'].shape
+        print outfile.variables['node_l2g'].shape
 
         outfile.variables['node_l2g'][:] = node_l2g.astype(num.int32)

trunk/anuga_core/source/anuga/operators/kinematic_viscosity_operator.py

@@ -54 +54 @@
             self.diffusivity = self.domain.get_quantity(diffusivity)
             
-    
+
+        self.smooth = 0.001
+        
         assert isinstance(self.diffusivity, Quantity)
         
@@ -114 +116 @@
         """ Parabolic update of x and y velocity
 
-        (I + dt (div h grad) ) U^{n+1} = U^{n}
+        (I + dt (div d grad) ) U^{n+1} = U^{n}
 
         """

trunk/anuga_core/source/anuga/operators/kinematic_viscosity_operator_ext.c

@@ -9 +9 @@
 // taken from http://cprogramminglanguage.net/quicksort-algorithm-c-source-code.aspx
 
-void swap(int *x,int *y)
-{
-   int temp;
-   temp = *x;
-   *x = *y;
-   *y = temp;
-}
-
-int choose_pivot(int i,int j )
-{
-   return((i+j) /2);
-}
-
-void quicksort(int list[],int m,int n)
-{
-   int key,i,j,k;
-   if( m < n)
-   {
-      k = choose_pivot(m,n);
-      swap(&list[m],&list[k]);
-      key = list[m];
-      i = m+1;
-      j = n;
-      while(i <= j)
-      {
-         while((i <= n) && (list[i] <= key))
+void swap(int *x, int *y) {
+    int temp;
+    temp = *x;
+    *x = *y;
+    *y = temp;
+}
+
+int choose_pivot(int i, int j) {
+    return ((i + j) / 2);
+}
+
+void quicksort(int list[], int m, int n) {
+    int key, i, j, k;
+    if (m < n) {
+        k = choose_pivot(m, n);
+        swap(&list[m], &list[k]);
+        key = list[m];
+        i = m + 1;
+        j = n;
+        while (i <= j) {
+            while ((i <= n) && (list[i] <= key))
                 i++;
-         while((j >= m) && (list[j] > key))
+            while ((j >= m) && (list[j] > key))
                 j--;
-         if( i < j)
-                swap(&list[i],&list[j]);
-      }
-          // swap two elements
-      swap(&list[m],&list[j]);
-          // recursively sort the lesser list
-      quicksort(list,m,j-1);
-      quicksort(list,j+1,n);
-   }
-}
-
-int build_geo_structure(int n, 
-                        int tot_len, 
-                        double *centroids, 
-                        long *neighbours, 
-                        double *edgelengths, 
-                        double *edge_midpoints, 
-                        long *geo_indices, 
-                        double *geo_values) {
+            if (i < j)
+                swap(&list[i], &list[j]);
+        }
+        // swap two elements
+        swap(&list[m], &list[j]);
+        // recursively sort the lesser list
+        quicksort(list, m, j - 1);
+        quicksort(list, j + 1, n);
+    }
+}
+
+int build_geo_structure(int n,
+        int tot_len,
+        double *centroids,
+        long *neighbours,
+        double *edgelengths,
+        double *edge_midpoints,
+        long *geo_indices,
+        double *geo_values) {
     int i, edge, edge_counted, j, m;
-        double dist, this_x, this_y, other_x, other_y, edge_length;
-        edge_counted = 0;
-        for (i=0; i<n; i++) {
-                //The centroid coordinates of triangle i
-                this_x = centroids[2*i];
-                this_y = centroids[2*i+1];
-                for (edge=0; edge<3; edge++) {
-                
-                        j = neighbours[3*i+edge];
-                        
-                        //Get the index and the coordinates of the interacting point
-
-                        // Edge
-                        if (j < 0 ) {
+    double dist, this_x, this_y, other_x, other_y, edge_length;
+    edge_counted = 0;
+    for (i = 0; i < n; i++) {
+        //The centroid coordinates of triangle i
+        this_x = centroids[2 * i];
+        this_y = centroids[2 * i + 1];
+        for (edge = 0; edge < 3; edge++) {
+
+            j = neighbours[3 * i + edge];
+
+            //Get the index and the coordinates of the interacting point
+
+            // Edge
+            if (j < 0) {
                 m = -j - 1;
-                                geo_indices[3*i+edge] = n + m;
-                                edge_counted++;
-                                
-                                other_x = edge_midpoints[2*(3*i+edge)];
-                                other_y = edge_midpoints[2*(3*i+edge)+1];
-                        } else {
-                                geo_indices[3*i+edge] = j;
-                                
-                                other_x = centroids[2*j];
-                                other_y = centroids[2*j+1];
-                        }
-                        
-                        //Compute the interaction
-                        edge_length = edgelengths[3*i+edge];
-                        dist = sqrt((this_x-other_x)*(this_x-other_x) + (this_y-other_y)*(this_y-other_y));
-                        geo_values[3*i+edge] = - edge_length / dist;
-                }
-        }
-        return 0;
-}
-
+                geo_indices[3 * i + edge] = n + m;
+                edge_counted++;
+
+                other_x = edge_midpoints[2 * (3 * i + edge)];
+                other_y = edge_midpoints[2 * (3 * i + edge) + 1];
+            } else {
+                geo_indices[3 * i + edge] = j;
+
+                other_x = centroids[2 * j];
+                other_y = centroids[2 * j + 1];
+            }
+
+            //Compute the interaction
+            edge_length = edgelengths[3 * i + edge];
+            dist = sqrt((this_x - other_x)*(this_x - other_x) + (this_y - other_y)*(this_y - other_y));
+            geo_values[3 * i + edge] = -edge_length / dist;
+        }
+    }
+    return 0;
+}
 
 int build_elliptic_matrix_not_symmetric(int n,
-                          int tot_len,
-                          long *geo_indices,
-                          double *geo_values,
-                          double *cell_data,
-                          double *bdry_data,
-                          double *data,
-                          long *colind) {
-        int i, k, edge, j[4], sorted_j[4], this_index;
-        double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
-        for (i=0; i<n; i++) {
-                v_i = 0.0;
-                j[3] = i;
-                //Get the values of each interaction, and the column index at which they occur
-        for (edge=0; edge<3; edge++) {
-            j[edge] = geo_indices[3*i+edge];
-            //if (j[edge]<n) { //interior
-            //    h_j = cell_data[j[edge]];
-            //} else { //boundary
-            //    h_j = bdry_data[j[edge]-n];
-            //}
-            v[edge] = -cell_data[i]*geo_values[3*i+edge]; //the negative of the individual interaction
-            v_i += cell_data[i]*geo_values[3*i+edge]; //sum the three interactions
-        }
-        if (cell_data[i]<=0.0) {
-            v_i  = 0.0;
-            v[0] = 0.0;
-            v[1] = 0.0;
-            v[2] = 0.0;
-        }
-                //Organise the set of 4 values/indices into the data and colind arrays
-                for (k=0; k<4; k++) sorted_j[k] = j[k];
-                quicksort(sorted_j,0,3);
-                for (k=0; k<4; k++) { //loop through the nonzero indices
-                        this_index = sorted_j[k];
-                        if (this_index == i) {
-                                data[4*i+k] = v_i;
-                                colind[4*i+k] = i;
-                        } else if (this_index == j[0]) {
-                                data[4*i+k] = v[0];
-                                colind[4*i+k] = j[0];
-                        } else if (this_index == j[1]) {
-                                data[4*i+k] = v[1];
-                                colind[4*i+k] = j[1];
-                        } else { //this_index == j[2]
-                                data[4*i+k] = v[2];
-                                colind[4*i+k] = j[2];
-                        }
-                }
-        }
-        return 0;
-}
-
-int build_elliptic_matrix(int n,
-                          int tot_len,
-                          long *geo_indices,
-                          double *geo_values,
-                          double *cell_data,
-                          double *bdry_data,
-                          double *data,
-                          long *colind) {
-        int i, k, edge, j[4], sorted_j[4], this_index;
-        double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
-        for (i=0; i<n; i++) {
-                v_i = 0.0;
-                j[3] = i;
-                //Get the values of each interaction, and the column index at which they occur
-        for (edge=0; edge<3; edge++) {
-            j[edge] = geo_indices[3*i+edge];
-            if (j[edge]<n) { //interior
-                h_j = cell_data[j[edge]];
-            } else { //boundary
-                h_j = bdry_data[j[edge]-n];
-            }
-            v[edge] = -0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
-            v_i += 0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
-        }
-        if (cell_data[i]<=0.0) {
-            v_i  = 0.0;
-            v[0] = 0.0;
-            v[1] = 0.0;
-            v[2] = 0.0;
-        }
-                //Organise the set of 4 values/indices into the data and colind arrays
-                for (k=0; k<4; k++) sorted_j[k] = j[k];
-                quicksort(sorted_j,0,3);
-                for (k=0; k<4; k++) { //loop through the nonzero indices
-                        this_index = sorted_j[k];
-                        if (this_index == i) {
-                                data[4*i+k] = v_i;
-                                colind[4*i+k] = i;
-                        } else if (this_index == j[0]) {
-                                data[4*i+k] = v[0];
-                                colind[4*i+k] = j[0];
-                        } else if (this_index == j[1]) {
-                                data[4*i+k] = v[1];
-                                colind[4*i+k] = j[1];
-                        } else { //this_index == j[2]
-                                data[4*i+k] = v[2];
-                                colind[4*i+k] = j[2];
-                        }
-                }
-        }
-        return 0;
-}
-
-
-int update_elliptic_matrix_not_symmetric(int n,
         int tot_len,
         long *geo_indices,
@@ -215 +102 @@
         v_i = 0.0;
         j[3] = i;
-
         //Get the values of each interaction, and the column index at which they occur
         for (edge = 0; edge < 3; edge++) {
             j[edge] = geo_indices[3 * i + edge];
-            if (j[edge] < n) { //interior
-                h_j = cell_data[j[edge]];
-            } else { //boundary
-                h_j = bdry_data[j[edge] - n];
-            }
-            v[edge] = - cell_data[i] * geo_values[3 * i + edge]; //the negative of the individual interaction
-            v_i +=  cell_data[i] * geo_values[3 * i + edge]; //sum the three interactions
+            //if (j[edge]<n) { //interior
+            //    h_j = cell_data[j[edge]];
+            //} else { //boundary
+            //    h_j = bdry_data[j[edge]-n];
+            //}
+            v[edge] = -cell_data[i] * geo_values[3 * i + edge]; //the negative of the individual interaction
+            v_i += cell_data[i] * geo_values[3 * i + edge]; //sum the three interactions
         }
         if (cell_data[i] <= 0.0) {
@@ -256 +142 @@
 }
 
-int update_elliptic_matrix(int n,
+int build_elliptic_matrix(int n,
         int tot_len,
         long *geo_indices,
@@ -269 +155 @@
         v_i = 0.0;
         j[3] = i;
-
         //Get the values of each interaction, and the column index at which they occur
         for (edge = 0; edge < 3; edge++) {
@@ -310 +195 @@
 }
 
-
+int update_elliptic_matrix_not_symmetric(int n,
+        int tot_len,
+        long *geo_indices,
+        double *geo_values,
+        double *cell_data,
+        double *bdry_data,
+        double *data,
+        long *colind) {
+    int i, k, edge, j[4], sorted_j[4], this_index;
+    double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
+    for (i = 0; i < n; i++) {
+        v_i = 0.0;
+        j[3] = i;
+
+        //Get the values of each interaction, and the column index at which they occur
+        for (edge = 0; edge < 3; edge++) {
+            j[edge] = geo_indices[3 * i + edge];
+            if (j[edge] < n) { //interior
+                h_j = cell_data[j[edge]];
+            } else { //boundary
+                h_j = bdry_data[j[edge] - n];
+            }
+            v[edge] = -cell_data[i] * geo_values[3 * i + edge]; //the negative of the individual interaction
+            v_i += cell_data[i] * geo_values[3 * i + edge]; //sum the three interactions
+        }
+        if (cell_data[i] <= 0.0) {
+            v_i = 0.0;
+            v[0] = 0.0;
+            v[1] = 0.0;
+            v[2] = 0.0;
+        }
+        //Organise the set of 4 values/indices into the data and colind arrays
+        for (k = 0; k < 4; k++) sorted_j[k] = j[k];
+        quicksort(sorted_j, 0, 3);
+        for (k = 0; k < 4; k++) { //loop through the nonzero indices
+            this_index = sorted_j[k];
+            if (this_index == i) {
+                data[4 * i + k] = v_i;
+                colind[4 * i + k] = i;
+            } else if (this_index == j[0]) {
+                data[4 * i + k] = v[0];
+                colind[4 * i + k] = j[0];
+            } else if (this_index == j[1]) {
+                data[4 * i + k] = v[1];
+                colind[4 * i + k] = j[1];
+            } else { //this_index == j[2]
+                data[4 * i + k] = v[2];
+                colind[4 * i + k] = j[2];
+            }
+        }
+    }
+    return 0;
+}
+
+int update_elliptic_matrix(int n,
+        int tot_len,
+        long *geo_indices,
+        double *geo_values,
+        double *cell_data,
+        double *bdry_data,
+        double *data,
+        long *colind) {
+    int i, k, edge, j[4], sorted_j[4], this_index;
+    double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
+    for (i = 0; i < n; i++) {
+        v_i = 0.0;
+        j[3] = i;
+
+        //Get the values of each interaction, and the column index at which they occur
+        for (edge = 0; edge < 3; edge++) {
+            j[edge] = geo_indices[3 * i + edge];
+            if (j[edge] < n) { //interior
+                h_j = cell_data[j[edge]];
+            } else { //boundary
+                h_j = bdry_data[j[edge] - n];
+            }
+            v[edge] = -0.5 * (cell_data[i] + h_j) * geo_values[3 * i + edge]; //the negative of the individual interaction
+            v_i += 0.5 * (cell_data[i] + h_j) * geo_values[3 * i + edge]; //sum the three interactions
+        }
+        if (cell_data[i] <= 0.0) {
+            v_i = 0.0;
+            v[0] = 0.0;
+            v[1] = 0.0;
+            v[2] = 0.0;
+        }
+        //Organise the set of 4 values/indices into the data and colind arrays
+        for (k = 0; k < 4; k++) sorted_j[k] = j[k];
+        quicksort(sorted_j, 0, 3);
+        for (k = 0; k < 4; k++) { //loop through the nonzero indices
+            this_index = sorted_j[k];
+            if (this_index == i) {
+                data[4 * i + k] = v_i;
+                colind[4 * i + k] = i;
+            } else if (this_index == j[0]) {
+                data[4 * i + k] = v[0];
+                colind[4 * i + k] = j[0];
+            } else if (this_index == j[1]) {
+                data[4 * i + k] = v[1];
+                colind[4 * i + k] = j[1];
+            } else { //this_index == j[2]
+                data[4 * i + k] = v[2];
+                colind[4 * i + k] = j[2];
+            }
+        }
+    }
+    return 0;
+}
 
 /**
-* Python wrapper methods
-*/
+ * Python wrapper methods
+ */
 
 static PyObject *py_build_geo_structure(PyObject *self, PyObject *args) {
@@ -324 +315 @@
             *edgelengths,
             *edge_midpoint_coordinates,
-                        *geo_indices,
-                        *geo_values;
-    
+            *geo_indices,
+            *geo_values;
+
     //Convert Python arguments to C
     if (!PyArg_ParseTuple(args, "O", &kv_operator)) {
@@ -339 +330 @@
 
     //Extract parameters
-    n       = get_python_integer(kv_operator,"n");
-    tot_len = get_python_integer(kv_operator,"tot_len");
-
-    centroid_coordinates = get_consecutive_array(mesh,"centroid_coordinates");
-    neighbours = get_consecutive_array(mesh,"neighbours");
-    edgelengths = get_consecutive_array(mesh,"edgelengths");
-    edge_midpoint_coordinates = get_consecutive_array(mesh,"edge_midpoint_coordinates");
-    geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
-    geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
-    
+    n = get_python_integer(kv_operator, "n");
+    tot_len = get_python_integer(kv_operator, "tot_len");
+
+    centroid_coordinates = get_consecutive_array(mesh, "centroid_coordinates");
+    neighbours = get_consecutive_array(mesh, "neighbours");
+    edgelengths = get_consecutive_array(mesh, "edgelengths");
+    edge_midpoint_coordinates = get_consecutive_array(mesh, "edge_midpoint_coordinates");
+    geo_indices = get_consecutive_array(kv_operator, "geo_structure_indices");
+    geo_values = get_consecutive_array(kv_operator, "geo_structure_values");
+
     //Release
     Py_DECREF(mesh);
-    
-    err = build_geo_structure(n,tot_len, 
-                             (double *)centroid_coordinates -> data, 
-                             (long *)neighbours -> data, 
-                             (double *)edgelengths->data, 
-                             (double *)edge_midpoint_coordinates -> data, 
-                             (long *)geo_indices -> data, 
-                             (double *)geo_values -> data);
+
+    err = build_geo_structure(n, tot_len,
+            (double *) centroid_coordinates -> data,
+            (long *) neighbours -> data,
+            (double *) edgelengths->data,
+            (double *) edge_midpoint_coordinates -> data,
+            (long *) geo_indices -> data,
+            (double *) geo_values -> data);
     if (err != 0) {
         PyErr_SetString(PyExc_RuntimeError, "Could not build geo structure");
         return NULL;
     }
-    
+
     //Release the arrays
     Py_DECREF(centroid_coordinates);
@@ -371 +362 @@
     Py_DECREF(geo_indices);
     Py_DECREF(geo_values);
-    
+
     return Py_BuildValue("");
 }
 
 static PyObject *py_build_elliptic_matrix(PyObject *self, PyObject *args) {
-        PyObject *kv_operator;
-        int n, tot_len, err;
-        PyArrayObject
-                *cell_data,
-                *bdry_data,
-                *geo_indices,
-                *geo_values,
-                *_data,
-                *colind;
-        
-        //Convert Python arguments to C
+    PyObject *kv_operator;
+    int n, tot_len, err;
+    PyArrayObject
+            *cell_data,
+            *bdry_data,
+            *geo_indices,
+            *geo_values,
+            *_data,
+            *colind;
+
+    //Convert Python arguments to C
     if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) {
         PyErr_SetString(PyExc_RuntimeError, "build_elliptic_matrix could not parse input");
@@ -393 +384 @@
 
 
-    n       = get_python_integer(kv_operator,"n");
-    tot_len = get_python_integer(kv_operator,"tot_len");
-
-
-        geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
-        geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
-        _data = get_consecutive_array(kv_operator,"operator_data");
-        colind = get_consecutive_array(kv_operator,"operator_colind");
-        
-        err = build_elliptic_matrix(n,tot_len,
-                                                                (long *)geo_indices -> data, 
-                                                                (double *)geo_values -> data, 
-                                                                (double *)cell_data -> data,
-                                                                (double *)bdry_data -> data,
-                                                                (double *)_data -> data, 
-                                                                (long *)colind -> data);
+    n = get_python_integer(kv_operator, "n");
+    tot_len = get_python_integer(kv_operator, "tot_len");
+
+
+    geo_indices = get_consecutive_array(kv_operator, "geo_structure_indices");
+    geo_values = get_consecutive_array(kv_operator, "geo_structure_values");
+    _data = get_consecutive_array(kv_operator, "operator_data");
+    colind = get_consecutive_array(kv_operator, "operator_colind");
+
+    err = build_elliptic_matrix(n, tot_len,
+            (long *) geo_indices -> data,
+            (double *) geo_values -> data,
+            (double *) cell_data -> data,
+            (double *) bdry_data -> data,
+            (double *) _data -> data,
+            (long *) colind -> data);
     if (err != 0) {
         PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
         return NULL;
     }
-        
-        Py_DECREF(geo_indices);
-        Py_DECREF(geo_values);
-        Py_DECREF(_data);
-        Py_DECREF(colind);
-    
+
+    Py_DECREF(geo_indices);
+    Py_DECREF(geo_values);
+    Py_DECREF(_data);
+    Py_DECREF(colind);
+
     return Py_BuildValue("");
 }
 
-
 static PyObject *py_update_elliptic_matrix(PyObject *self, PyObject *args) {
-        PyObject *kv_operator;
-        int n, tot_len, err;
-        PyArrayObject
-                *cell_data,
-                *bdry_data,
-                *geo_indices,
-                *geo_values,
-                *_data,
-                *colind;
-
-        //Convert Python arguments to C
+    PyObject *kv_operator;
+    int n, tot_len, err;
+    PyArrayObject
+            *cell_data,
+            *bdry_data,
+            *geo_indices,
+            *geo_values,
+            *_data,
+            *colind;
+
+    //Convert Python arguments to C
     if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) {
         PyErr_SetString(PyExc_RuntimeError, "update_elliptic_matrix could not parse input");
@@ -441 +431 @@
 
 
-    n       = get_python_integer(kv_operator,"n");
-    tot_len = get_python_integer(kv_operator,"tot_len");
-
-
-        geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
-        geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
-        _data = get_consecutive_array(kv_operator,"operator_data");
-        colind = get_consecutive_array(kv_operator,"operator_colind");
-
-        err = update_elliptic_matrix(n,tot_len,
-                                                                (long *)geo_indices -> data,
-                                                                (double *)geo_values -> data,
-                                                                (double *)cell_data -> data,
-                                                                (double *)bdry_data -> data,
-                                                                (double *)_data -> data,
-                                                                (long *)colind -> data);
+    n = get_python_integer(kv_operator, "n");
+    tot_len = get_python_integer(kv_operator, "tot_len");
+
+
+    geo_indices = get_consecutive_array(kv_operator, "geo_structure_indices");
+    geo_values = get_consecutive_array(kv_operator, "geo_structure_values");
+    _data = get_consecutive_array(kv_operator, "operator_data");
+    colind = get_consecutive_array(kv_operator, "operator_colind");
+
+    err = update_elliptic_matrix(n, tot_len,
+            (long *) geo_indices -> data,
+            (double *) geo_values -> data,
+            (double *) cell_data -> data,
+            (double *) bdry_data -> data,
+            (double *) _data -> data,
+            (long *) colind -> data);
     if (err != 0) {
         PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
@@ -462 +452 @@
     }
 
-        Py_DECREF(geo_indices);
-        Py_DECREF(geo_values);
-        Py_DECREF(_data);
-        Py_DECREF(colind);
+    Py_DECREF(geo_indices);
+    Py_DECREF(geo_values);
+    Py_DECREF(_data);
+    Py_DECREF(colind);
 
     return Py_BuildValue("");
@@ -472 +462 @@
 
 static struct PyMethodDef MethodTable[] = {
-        {"build_geo_structure",py_build_geo_structure,METH_VARARGS,"Print out"},
-                {"build_elliptic_matrix",py_build_elliptic_matrix,METH_VARARGS,"Print out"},
-        {"update_elliptic_matrix",py_update_elliptic_matrix,METH_VARARGS,"Print out"},
-        {NULL,NULL,0,NULL} // sentinel
+    {"build_geo_structure", py_build_geo_structure, METH_VARARGS, "Print out"},
+    {"build_elliptic_matrix", py_build_elliptic_matrix, METH_VARARGS, "Print out"},
+    {"update_elliptic_matrix", py_update_elliptic_matrix, METH_VARARGS, "Print out"},
+    {NULL, NULL, 0, NULL} // sentinel
 };
-void initkinematic_viscosity_operator_ext(){
+
+void initkinematic_viscosity_operator_ext() {
     (void) Py_InitModule("kinematic_viscosity_operator_ext", MethodTable);
     import_array();

trunk/anuga_core/source/anuga/utilities/sww_merge.py

@@ -24 +24 @@
     swwfiles = [ domain_global_name+"_P"+str(np)+"_"+str(v)+".sww" for v in range(np)]
 
-    _sww_merge_parallel(swwfiles, output, verbose, delete_old)
+    _sww_merge_parallel_non_smooth(swwfiles, output, verbose, delete_old)
 
 
@@ -176 +176 @@
 
 
-def _sww_merge_parallel(swwfiles, output,  verbose=False, delete_old=False):
+def _sww_merge_parallel_smooth(swwfiles, output,  verbose=False, delete_old=False):
     """
         Merge a list of sww files into a single file.
@@ -183 +183 @@
 
         The sww files to be merged must have exactly the same timesteps.
+
+        It is assumed that the separate sww files have been stored in non_smooth
+        format.
 
         Note that some advanced information and custom quantities may not be
@@ -429 +432 @@
             os.remove(filename)
 
+def _sww_merge_parallel_non_smooth(swwfiles, output,  verbose=False, delete_old=False):
+    """
+        Merge a list of sww files into a single file.
+
+        Used to merge files created by parallel runs.
+
+        The sww files to be merged must have exactly the same timesteps.
+
+        It is assumed that the separate sww files have been stored in non_smooth
+        format.
+
+        Note that some advanced information and custom quantities may not be
+        exported.
+
+        swwfiles is a list of .sww files to merge.
+        output is the output filename, including .sww extension.
+        verbose True to log output information
+    """
+
+    if verbose:
+        print "MERGING SWW Files"
+
+
+    first_file = True
+    tri_offset = 0
+    for filename in swwfiles:
+        if verbose:
+            print 'Reading file ', filename, ':'
+
+        fid = NetCDFFile(filename, netcdf_mode_r)
+
+        if first_file:
+
+            times    = fid.variables['time'][:]
+            n_steps = len(times)
+            number_of_timesteps = fid.dimensions['number_of_timesteps']
+            #print n_steps, number_of_timesteps
+            starttime = int(fid.starttime)
+
+            out_s_quantities = {}
+            out_d_quantities = {}
+
+
+            xllcorner = fid.xllcorner
+            yllcorner = fid.yllcorner
+
+            number_of_global_triangles = int(fid.number_of_global_triangles)
+            number_of_global_nodes     = int(fid.number_of_global_nodes)
+            number_of_global_triangle_vertices = 3*number_of_global_triangles
+
+            order      = fid.order
+            xllcorner  = fid.xllcorner;
+            yllcorner  = fid.yllcorner ;
+            zone       = fid.zone;
+            false_easting  = fid.false_easting;
+            false_northing = fid.false_northing;
+            datum      = fid.datum;
+            projection = fid.projection;
+
+            g_volumes = num.zeros((number_of_global_triangles,3),num.int)
+            g_x = num.zeros((number_of_global_nodes,),num.float32)
+            g_y = num.zeros((number_of_global_nodes,),num.float32)
+
+            g_points = num.zeros((number_of_global_nodes,2),num.float32)
+
+            quantities = ['elevation', 'stage', 'xmomentum', 'ymomentum']
+            static_quantities = []
+            dynamic_quantities = []
+
+            for quantity in quantities:
+                # Test if elevation is static
+                if n_steps == fid.variables[quantity].shape[0]:
+                    dynamic_quantities.append(quantity)
+                else:
+                    static_quantities.append(quantity)
+
+            for quantity in static_quantities:
+                out_s_quantities[quantity] = num.zeros((3*number_of_global_triangles,),num.float32)
+
+            # Quantities are stored as a 2D array of timesteps x data.
+            for quantity in dynamic_quantities:
+                out_d_quantities[quantity] = \
+                      num.zeros((n_steps,3*number_of_global_triangles),num.float32)
+
+            description = 'merged:' + getattr(fid, 'description')
+            first_file = False
+
+
+        # Read in from files and add to global arrays
+
+        tri_l2g  = fid.variables['tri_l2g'][:]
+        node_l2g = fid.variables['node_l2g'][:]
+        tri_full_flag = fid.variables['tri_full_flag'][:]
+
+        volumes = num.array(fid.variables['volumes'][:],dtype=num.int)
+
+
+
+        #l_volumes = num.zeros_like(volumes)
+        #l_old_volumes = num.zeros_like(volumes)
+
+
+        # Change the local node ids to global id in the
+        # volume array
+
+
+        print volumes
+        print volumes.shape
+
+
+        print fid.variables['x'][0:6]
+        print fid.variables['y'][0:6]
+
+        print fid.variables['x'].shape
+
+
+        g_n0 = node_l2g[volumes[:,0]].reshape(-1,1)
+        g_n1 = node_l2g[volumes[:,1]].reshape(-1,1)
+        g_n2 = node_l2g[volumes[:,2]].reshape(-1,1)
+
+        #print g_n0.shape
+        l_volumes = num.hstack((g_n0,g_n1,g_n2))
+
+        #assert num.allclose(l_volumes, l_old_volumes)
+
+        # Just pick out the full triangles
+        #ftri_l2g = num.compress(tri_full_flag, tri_l2g)
+
+        #print l_volumes
+        #print tri_full_flag
+        print l_volumes.shape
+        print g_volumes.shape
+        print tri_l2g.shape
+
+        #fg_volumes = num.compress(tri_full_flag,l_volumes,axis=0)
+        g_volumes[tri_l2g] = l_volumes
+
+
+
+
+        #g_x[node_l2g] = fid.variables['x']
+        #g_y[node_l2g] = fid.variables['y']
+
+        g_points[tri_l2g,0] = fid.variables['x']
+        g_points[tri_l2g,1] = fid.variables['y']
+
+
+        #print number_of_timesteps
+
+
+        # FIXME SR: It seems that some of the "ghost" node quantity values
+        # are being storded. We should only store those nodes which are associated with
+        # full triangles. So we need an index array of "full" nodes, ie those in
+        # full triangles
+
+        #use numpy.compress and numpy.unique to get "full nodes
+
+        #f_volumes = num.compress(tri_full_flag,volumes,axis=0)
+        #fl_nodes = num.unique(f_volumes)
+        #f_node_l2g = node_l2g[fl_nodes]
+
+        #print len(node_l2g)
+        #print len(fl_nodes)
+
+        # Read in static quantities
+        for quantity in static_quantities:
+            #out_s_quantities[quantity][node_l2g] = \
+            #             num.array(fid.variables[quantity],dtype=num.float32)
+            q = fid.variables[quantity]
+            #print quantity, q.shape
+            out_s_quantities[quantity][f_node_l2g] = \
+                         num.array(q,dtype=num.float32)[fl_nodes]
+
+
+        #Collate all dynamic quantities according to their timestep
+        for quantity in dynamic_quantities:
+            q = fid.variables[quantity]
+            #print q.shape
+            for i in range(n_steps):
+                #out_d_quantities[quantity][i][node_l2g] = \
+                #           num.array(q[i],dtype=num.float32)
+                out_d_quantities[quantity][i][f_node_l2g] = \
+                           num.array(q[i],dtype=num.float32)[fl_nodes]
+
+
+
+
+        fid.close()
+
+
+    #---------------------------
+    # Write out the SWW file
+    #---------------------------
+    #print g_points.shape
+
+    #print number_of_global_triangles
+    #print number_of_global_nodes
+
+
+    if verbose:
+            print 'Writing file ', output, ':'
+    fido = NetCDFFile(output, netcdf_mode_w)
+    sww = Write_sww(static_quantities, dynamic_quantities)
+    sww.store_header(fido, starttime,
+                             number_of_global_triangles,
+                             number_of_global_nodes,
+                             description=description,
+                             sww_precision=netcdf_float32)
+
+
+
+    from anuga.coordinate_transforms.geo_reference import Geo_reference
+    geo_reference = Geo_reference()
+
+    sww.store_triangulation(fido, g_points, g_volumes, points_georeference=geo_reference)
+
+    fido.order      = order
+    fido.xllcorner  = xllcorner;
+    fido.yllcorner  = yllcorner ;
+    fido.zone       = zone;
+    fido.false_easting  = false_easting;
+    fido.false_northing = false_northing;
+    fido.datum      = datum;
+    fido.projection = projection;
+
+    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
+
+    # Write out all the dynamic quantities for each timestep
+
+    for i in range(n_steps):
+        fido.variables['time'][i] = times[i]
+
+    for q in dynamic_quantities:
+        q_values = out_d_quantities[q]
+        for i in range(n_steps):
+            fido.variables[q][i] = q_values[i]
+
+        # This updates the _range values
+        q_range = fido.variables[q + Write_sww.RANGE][:]
+        q_values_min = num.min(q_values)
+        if q_values_min < q_range[0]:
+            fido.variables[q + Write_sww.RANGE][0] = q_values_min
+        q_values_max = num.max(q_values)
+        if q_values_max > q_range[1]:
+            fido.variables[q + Write_sww.RANGE][1] = q_values_max
+
+
+    #print out_s_quantities
+    #print out_d_quantities
+
+    #print g_x
+    #print g_y
+
+    #print g_volumes
+
+    fido.close()
+
+    if delete_old:
+        import os
+        for filename in swwfiles:
+
+            if verbose:
+                print 'Deleting file ', filename, ':'
+            os.remove(filename)
+
+
+
+
+
+
 if __name__ == "__main__":