Changeset 8123


Ignore:
Timestamp:
Feb 24, 2011, 5:53:54 PM (13 years ago)
Author:
steve
Message:

Added in a procedure to calculate local gradients on each triangle using
vectex values

Location:
trunk/anuga_core/source/anuga/abstract_2d_finite_volumes
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r8017 r8123  
    14721472        return compute_gradients(self)
    14731473
     1474
     1475    ##
     1476    # @brief ??
     1477    def compute_local_gradients(self):
     1478        # Call correct module function
     1479        # (either from this module or C-extension)
     1480        return compute_local_gradients(self)
     1481
     1482
     1483
     1484
    14741485    ##
    14751486    # @brief ??
     
    15821593         saxpy_centroid_values,\
    15831594         compute_gradients,\
     1595         compute_local_gradients,\
    15841596         limit_old,\
    15851597         limit_vertices_by_all_neighbours,\
  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r7616 r8123  
    9393}
    9494
     95
     96int _compute_local_gradients(int N,
     97                               double* vertex_coordinates,
     98                               double* vertex_values,
     99                               double* a,
     100                               double* b) {
     101
     102  int k, k2, k3, k6;
     103  double x0, y0, x1, y1, x2, y2, v0, v1, v2;
     104 
     105  for (k=0; k<N; k++) {
     106    k6 = 6*k;
     107    k3 = 3*k;
     108    k2 = 2*k;
     109     
     110    // vertex coordinates
     111    // x0, y0, x1, y1, x2, y2 = X[k,:]
     112    x0 = vertex_coordinates[k6 + 0];
     113    y0 = vertex_coordinates[k6 + 1];
     114    x1 = vertex_coordinates[k6 + 2];
     115    y1 = vertex_coordinates[k6 + 3];
     116    x2 = vertex_coordinates[k6 + 4];
     117    y2 = vertex_coordinates[k6 + 5];
     118     
     119    v0 = vertex_values[k3+0];
     120    v1 = vertex_values[k3+1];
     121    v2 = vertex_values[k3+2];
     122   
     123    // Gradient
     124    _gradient(x0, y0, x1, y1, x2, y2, v0, v1, v2, &a[k], &b[k]);
     125
     126   
     127    }
     128    return 0;
     129}
    95130
    96131int _extrapolate_from_gradient(int N,
     
    11621197
    11631198
     1199
     1200PyObject *compute_local_gradients(PyObject *self, PyObject *args) {
     1201
     1202        PyObject *quantity, *domain;
     1203        PyArrayObject
     1204            *vertex_coordinates,   //Coordinates at vertices
     1205            *vertex_values,        //Values at vertices
     1206            *x_gradient,           //x gradient
     1207            *y_gradient;           //y gradient
     1208
     1209        //int N, err;
     1210        //int dimensions[1];
     1211        int N, err;
     1212        //double *a, *b;  //Gradients
     1213
     1214        // Convert Python arguments to C
     1215        if (!PyArg_ParseTuple(args, "O", &quantity)) {
     1216          PyErr_SetString(PyExc_RuntimeError,
     1217                          "compute_local_gradient could not parse input");     
     1218          return NULL;
     1219        }
     1220
     1221        domain = PyObject_GetAttrString(quantity, "domain");
     1222        if (!domain) {
     1223          PyErr_SetString(PyExc_RuntimeError,
     1224                          "compute_local_gradient could not obtain domain object from quantity");       
     1225          return NULL;
     1226        }
     1227
     1228        // Get pertinent variables
     1229        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
     1230        vertex_values        = get_consecutive_array(quantity, "vertex_values");
     1231        x_gradient           = get_consecutive_array(quantity, "x_gradient");
     1232        y_gradient           = get_consecutive_array(quantity, "y_gradient");
     1233
     1234        N = vertex_values -> dimensions[0];
     1235
     1236        // Release
     1237        Py_DECREF(domain);
     1238
     1239        err = _compute_local_gradients(N,
     1240                        (double*) vertex_coordinates -> data,
     1241                        (double*) vertex_values -> data,
     1242                        (double*) x_gradient -> data,
     1243                        (double*) y_gradient -> data);
     1244
     1245
     1246        if (err != 0) {
     1247          PyErr_SetString(PyExc_RuntimeError,
     1248                          "Internal function _compute_local_gradient failed");
     1249          return NULL;
     1250        }
     1251
     1252
     1253
     1254        // Release
     1255        Py_DECREF(vertex_coordinates);
     1256        Py_DECREF(vertex_values);
     1257        Py_DECREF(x_gradient);
     1258        Py_DECREF(y_gradient);
     1259
     1260        return Py_BuildValue("");
     1261}
     1262
     1263
    11641264PyObject *extrapolate_second_order_and_limit_by_edge(PyObject *self, PyObject *args) {
    11651265    /* Compute edge values using second order approximation and limit values
     
    20932193        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
    20942194        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
     2195    {"compute_local_gradients", compute_gradients, METH_VARARGS, "Print out"},
    20952196        {"extrapolate_from_gradient", extrapolate_from_gradient,
    20962197                METH_VARARGS, "Print out"},
Note: See TracChangeset for help on using the changeset viewer.