n_gradient_calc.c

Go to the documentation of this file.

/*****************************************************************************
 *
 * MODULE:       Grass PDE Numerical Library
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
 *                 soerengebbert <at> gmx <dot> de
 *
 * PURPOSE:      gradient management functions
 *               part of the gpde library
 *
 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *
 *****************************************************************************/
 
#include <grass/N_pde.h>
 
/*! \brief Calculate basic statistics of a gradient field
 *
 * The statistic is stored in the gradient field struct
 *
 * \param field N_gradient_2d_field *
 * \return void
 *
 * */
void N_calc_gradient_field_2d_stats(N_gradient_field_2d *field)
{
    double minx, miny;
    double maxx, maxy;
    double sumx, sumy;
    int nonullx, nonully;
 
    G_debug(3, "N_calc_gradient_field_2d_stats: compute gradient field stats");
 
    N_calc_array_2d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
    N_calc_array_2d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
 
    if (minx < miny)
        field->min = minx;
    else
        field->min = miny;
 
    if (maxx > maxy)
        field->max = maxx;
    else
        field->max = maxy;
 
    field->sum = sumx + sumy;
    field->nonull = nonullx + nonully;
    field->mean = field->sum / (double)field->nonull;
 
    return;
}
 
/*!
 * \brief This function computes the gradient based on the input N_array_2d pot
 * (potential), a weighting factor N_array_2d named weight and the distance
 between two cells
 * saved in the N_geom_data struct.
 *
 * The gradient is calculated between cells for each cell and direction.
 * An existing gradient field can be filled with new data or, if a NULL pointer
 is
 * given, a new gradient field will be allocated with the appropriate size.
 *
 *
 \verbatim
 ______________
 |    |    |    |
 |    |    |    |
 |----|-NC-|----|
 |    |    |    |
 |   WC    EC   |
 |    |    |    |
 |----|-SC-|----|
 |    |    |    |
 |____|____|____|
 
 
 x - direction:
 
 r = 2 * weight[row][col]*weight[row][col + 1] /
 (weight[row][col]*weight[row][col + 1]) EC = r * (pot[row][col] - pot[row][col
 + 1])/dx
 
 y - direction:
 
 r = 2 * weight[row][col]*weight[row + 1][col] / (weight[row][col]*weight[row +
 1][col]) SC = r * (pot[row][col] - pot[row + 1][col])/dy
 
 the values SC and EC are the values of the next row/col
 
 
 \endverbatim
 * \param pot N_array_2d * - the potential N_array_2d
 * \param weight_x N_array_2d * - the weighting factor N_array_2d used to modify
 the gradient in x-direction
 * \param weight_y N_array_2d * - the weighting factor N_array_2d used to modify
 the gradient in y-direction
 * \param geom N_geom_data * - geometry data structure
 * \param gradfield N_gradient_field_2d * - a gradient field of the correct
 size, if a NULL pointer is provided this gradient field will be new allocated
 * \return N_gradient_field_2d * - the pointer to the computed gradient field
 
 *
 * */
N_gradient_field_2d *N_compute_gradient_field_2d(N_array_2d *pot,
                                                 N_array_2d *weight_x,
                                                 N_array_2d *weight_y,
                                                 N_geom_data *geom,
                                                 N_gradient_field_2d *gradfield)
{
    int i, j;
    int rows, cols;
    double dx, dy, p1, p2, r1, r2, mean, grad, res;
    N_gradient_field_2d *field = gradfield;
 
    if (pot->cols != weight_x->cols || pot->cols != weight_y->cols)
        G_fatal_error(
            "N_compute_gradient_field_2d: the arrays are not of equal size");
 
    if (pot->rows != weight_x->rows || pot->rows != weight_y->rows)
        G_fatal_error(
            "N_compute_gradient_field_2d: the arrays are not of equal size");
 
    if (pot->cols != geom->cols || pot->rows != geom->rows)
        G_fatal_error("N_compute_gradient_field_2d: array sizes and geometry "
                      "data are different");
 
    G_debug(3, "N_compute_gradient_field_2d: compute gradient field");
 
    rows = pot->rows;
    cols = pot->cols;
    dx = geom->dx;
    dy = geom->dy;
 
    if (field == NULL) {
        field = N_alloc_gradient_field_2d(cols, rows);
    }
    else {
        if (field->cols != geom->cols || field->rows != geom->rows)
            G_fatal_error("N_compute_gradient_field_2d: gradient field sizes "
                          "and geometry data are different");
    }
 
    for (j = 0; j < rows; j++)
        for (i = 0; i < cols - 1; i++) {
            grad = 0;
            mean = 0;
 
            /* Only compute if the arrays are not null */
            if (!N_is_array_2d_value_null(pot, i, j) &&
                !N_is_array_2d_value_null(pot, i + 1, j)) {
                p1 = N_get_array_2d_d_value(pot, i, j);
                p2 = N_get_array_2d_d_value(pot, i + 1, j);
                grad = (p1 - p2) / dx; /* gradient */
            }
            if (!N_is_array_2d_value_null(weight_x, i, j) &&
                !N_is_array_2d_value_null(weight_x, i + 1, j)) {
                r1 = N_get_array_2d_d_value(weight_x, i, j);
                r2 = N_get_array_2d_d_value(weight_x, i + 1, j);
                mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
            }
 
            res = mean * grad;
 
            N_put_array_2d_d_value(field->x_array, i + 1, j, res);
        }
 
    for (j = 0; j < rows - 1; j++)
        for (i = 0; i < cols; i++) {
            grad = 0;
            mean = 0;
 
            /* Only compute if the arrays are not null */
            if (!N_is_array_2d_value_null(pot, i, j) &&
                !N_is_array_2d_value_null(pot, i, j + 1)) {
                p1 = N_get_array_2d_d_value(pot, i, j);
                p2 = N_get_array_2d_d_value(pot, i, j + 1);
                grad = (p1 - p2) / dy; /* gradient */
            }
            if (!N_is_array_2d_value_null(weight_y, i, j) &&
                !N_is_array_2d_value_null(weight_y, i, j + 1)) {
                r1 = N_get_array_2d_d_value(weight_y, i, j);
                r2 = N_get_array_2d_d_value(weight_y, i, j + 1);
                mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
            }
 
            res = -1 * mean * grad;
 
            N_put_array_2d_d_value(field->y_array, i, j + 1, res);
        }
 
    /*Compute gradient field statistics */
    N_calc_gradient_field_2d_stats(field);
 
    return field;
}
 
/*!
 * \brief Calculate the x and y vector components from a gradient field for each
 * cell and stores them in the provided N_array_2d structures
 *
 * The arrays must have the same size as the gradient field.
 
 \verbatim
 
 Based on this storages scheme the gradient vector for each cell is
 calculated and stored in the provided  N_array_2d structures
 
 ______________
 |    |    |    |
 |    |    |    |
 |----|-NC-|----|
 |    |    |    |
 |   WC    EC   |
 |    |    |    |
 |----|-SC-|----|
 |    |    |    |
 |____|____|____|
 
 x vector component:
 
 x = (WC + EC) / 2
 
 y vector component:
 
 y = (NC + SC) / 2
 
 \endverbatim
 *
 * \param field N_gradient_field_2d *
 * \param x_comp N_array_2d * - the array in which the x component will be
 written
 * \param y_comp N_array_2d * - the array in which the y component will be
 written
 *
 * \return void
 * */
void N_compute_gradient_field_components_2d(N_gradient_field_2d *field,
                                            N_array_2d *x_comp,
                                            N_array_2d *y_comp)
{
    int i, j;
 
    int rows, cols;
 
    double vx, vy;
 
    N_array_2d *x = x_comp;
 
    N_array_2d *y = y_comp;
 
    N_gradient_2d grad;
 
    if (!x)
        G_fatal_error("N_compute_gradient_components_2d: x array is empty");
    if (!y)
        G_fatal_error("N_compute_gradient_components_2d: y array is empty");
 
    cols = field->x_array->cols;
    rows = field->x_array->rows;
 
    /*Check the array sizes */
    if (x->cols != cols || x->rows != rows)
        G_fatal_error("N_compute_gradient_components_2d: the size of the x "
                      "array doesn't fit the gradient field size");
    if (y->cols != cols || y->rows != rows)
        G_fatal_error("N_compute_gradient_components_2d: the size of the y "
                      "array doesn't fit the gradient field size");
 
    for (j = 0; j < rows; j++)
        for (i = 0; i < cols; i++) {
            N_get_gradient_2d(field, &grad, i, j);
 
            /* in case a gradient is zero, we expect a no flow boundary */
            if (grad.WC == 0.0 || grad.EC == 0.0)
                vx = (grad.WC + grad.EC);
            else
                vx = (grad.WC + grad.EC) / 2;
            if (grad.NC == 0.0 || grad.SC == 0.0)
                vy = (grad.NC + grad.SC);
            else
                vy = (grad.NC + grad.SC) / 2;
 
            N_put_array_2d_d_value(x, i, j, vx);
            N_put_array_2d_d_value(y, i, j, vy);
        }
 
    return;
}
 
/*! \brief Calculate basic statistics of a gradient field
 *
 * The statistic is stored in the gradient field struct
 *
 * \param field N_gradient_3d_field *
 * \return void
 *
 * */
void N_calc_gradient_field_3d_stats(N_gradient_field_3d *field)
{
    double minx, miny, minz;
 
    double maxx, maxy, maxz;
 
    double sumx, sumy, sumz;
 
    int nonullx, nonully, nonullz;
 
    G_debug(3, "N_calc_gradient_field_3d_stats: compute gradient field stats");
 
    N_calc_array_3d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
    N_calc_array_3d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
    N_calc_array_3d_stats(field->z_array, &minz, &maxz, &sumz, &nonullz, 0);
 
    if (minx <= minz && minx <= miny)
        field->min = minx;
    if (miny <= minz && miny <= minx)
        field->min = miny;
    if (minz <= minx && minz <= miny)
        field->min = minz;
 
    if (maxx >= maxz && maxx >= maxy)
        field->max = maxx;
    if (maxy >= maxz && maxy >= maxx)
        field->max = maxy;
    if (maxz >= maxx && maxz >= maxy)
        field->max = maxz;
 
    field->sum = sumx + sumy + sumz;
    field->nonull = nonullx + nonully + nonullz;
    field->mean = field->sum / (double)field->nonull;
 
    return;
}
 
/*!
 * \brief This function computes the gradient based on the input N_array_3d pot
 * (that means potential), a weighting factor N_array_3d named weight and the
 distance between two cells
 * saved in the N_geom_data struct.
 *
 * The gradient is calculated between cells for each cell and direction.
 * An existing gradient field can be filled with new data or, if a NULL pointer
 is
 * given, a new gradient field will be allocated with the appropriate size.
 *
 *
 *
 *
 \verbatim
 
 |  /
 TC NC
 |/
 --WC-----EC--
 /|
 SC BC
 /  |
 
 x - direction:
 
 r = 2 * weight_x[depth][row][col]*weight_x[depth][row][col + 1] /
 (weight_X[depth][row][col]*weight_x[depth][row][col + 1]) EC = r *
 (pot[depth][row][col] - pot[depth][row][col + 1])/dx
 
 y - direction:
 
 r = 2 * weight_y[depth][row][col]*weight_y[depth][row + 1][col] /
 (weight_y[depth][row][col]*weight_y[depth][row + 1][col]) SC = r *
 (pot[depth][row][col] - pot[depth][row + 1][col])/dy
 
 z - direction:
 
 r = 2 * weight_z[depth][row][col]*weight_z[depth + 1][row][col] /
 (weight_z[depth][row][col]*weight_z[depth + 1][row][col]) TC = r *
 (pot[depth][row][col] - pot[depth + 1][row][col])/dy
 
 the values BC, NC, WC are the values of the next depth/row/col
 
 
 \endverbatim
 * \param pot N_array_3d * - the potential N_array_2d
 * \param weight_x N_array_3d * - the weighting factor N_array_3d used to modify
 the gradient in x-direction
 * \param weight_y N_array_3d * - the weighting factor N_array_3d used to modify
 the gradient in y-direction
 * \param weight_z N_array_3d * - the weighting factor N_array_3d used to modify
 the gradient in z-direction
 * \param geom N_geom_data * - geometry data structure
 * \param gradfield N_gradient_field_3d * - a gradient field of the correct
 size, if a NULL pointer is provided this gradient field will be new allocated
 * \return N_gradient_field_3d * - the pointer to the computed gradient field
 *
 * */
N_gradient_field_3d *
N_compute_gradient_field_3d(N_array_3d *pot, N_array_3d *weight_x,
                            N_array_3d *weight_y, N_array_3d *weight_z,
                            N_geom_data *geom, N_gradient_field_3d *gradfield)
{
    int i, j, k;
 
    int cols, rows, depths;
 
    double dx, dy, dz, p1, p2, r1, r2, mean, grad, res;
 
    N_gradient_field_3d *field = gradfield;
 
    if (pot->cols != weight_x->cols || pot->cols != weight_y->cols ||
        pot->cols != weight_z->cols)
        G_fatal_error(
            "N_compute_gradient_field_3d: the arrays are not of equal size");
 
    if (pot->rows != weight_x->rows || pot->rows != weight_y->rows ||
        pot->rows != weight_z->rows)
        G_fatal_error(
            "N_compute_gradient_field_3d: the arrays are not of equal size");
 
    if (pot->depths != weight_x->depths || pot->depths != weight_y->depths ||
        pot->depths != weight_z->depths)
        G_fatal_error(
            "N_compute_gradient_field_3d: the arrays are not of equal size");
 
    if (pot->cols != geom->cols || pot->rows != geom->rows ||
        pot->depths != geom->depths)
        G_fatal_error("N_compute_gradient_field_3d: array sizes and geometry "
                      "data are different");
 
    G_debug(3, "N_compute_gradient_field_3d: compute gradient field");
 
    cols = geom->cols;
    rows = geom->rows;
    depths = geom->depths;
    dx = geom->dx;
    dy = geom->dy;
    dz = geom->dz;
 
    if (gradfield == NULL) {
        field = N_alloc_gradient_field_3d(cols, rows, depths);
    }
    else {
        if (field->cols != geom->cols || field->rows != geom->rows ||
            field->depths != geom->depths)
            G_fatal_error("N_compute_gradient_field_3d: gradient field sizes "
                          "and geometry data are different");
    }
 
    for (k = 0; k < depths; k++)
        for (j = 0; j < rows; j++)
            for (i = 0; i < cols - 1; i++) {
                grad = 0;
                mean = 0;
 
                /*Only compute if the arrays are not null */
                if (!N_is_array_3d_value_null(pot, i, j, k) &&
                    !N_is_array_3d_value_null(pot, i + 1, j, k)) {
                    p1 = N_get_array_3d_d_value(pot, i, j, k);
                    p2 = N_get_array_3d_d_value(pot, i + 1, j, k);
                    grad = (p1 - p2) / dx; /* gradient */
                }
                if (!N_is_array_3d_value_null(weight_x, i, j, k) &&
                    !N_is_array_3d_value_null(weight_x, i + 1, j, k)) {
                    r1 = N_get_array_3d_d_value(weight_x, i, j, k);
                    r2 = N_get_array_3d_d_value(weight_x, i + 1, j, k);
                    mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
                }
 
                res = mean * grad;
 
                G_debug(6,
                        "N_compute_gradient_field_3d: X-direction insert value "
                        "%6.5g at %i %i %i ",
                        res, k, j, i + 1);
 
                N_put_array_3d_d_value(field->x_array, i + 1, j, k, res);
            }
 
    for (k = 0; k < depths; k++)
        for (j = 0; j < rows - 1; j++)
            for (i = 0; i < cols; i++) {
                grad = 0;
                mean = 0;
 
                /* Only compute if the arrays are not null */
                if (!N_is_array_3d_value_null(pot, i, j, k) &&
                    !N_is_array_3d_value_null(pot, i, j + 1, k)) {
                    p1 = N_get_array_3d_d_value(pot, i, j, k);
                    p2 = N_get_array_3d_d_value(pot, i, j + 1, k);
                    grad = (p1 - p2) / dy; /* gradient */
                }
                if (!N_is_array_3d_value_null(weight_y, i, j, k) &&
                    !N_is_array_3d_value_null(weight_y, i, j + 1, k)) {
                    r1 = N_get_array_3d_d_value(weight_y, i, j, k);
                    r2 = N_get_array_3d_d_value(weight_y, i, j + 1, k);
                    mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
                }
 
                res = -1 * mean * grad; /*invert the direction, because we count
                                         * from north to south, but the gradient
                                         * is defined in y direction */
 
                G_debug(6,
                        "N_compute_gradient_field_3d: Y-direction insert value "
                        "%6.5g at %i %i %i ",
                        res, k, j + 1, i);
 
                N_put_array_3d_d_value(field->y_array, i, j + 1, k, res);
            }
 
    for (k = 0; k < depths - 1; k++)
        for (j = 0; j < rows; j++)
            for (i = 0; i < cols; i++) {
                grad = 0;
                mean = 0;
 
                /* Only compute if the arrays are not null */
                if (!N_is_array_3d_value_null(pot, i, j, k) &&
                    !N_is_array_3d_value_null(pot, i, j, k + 1)) {
                    p1 = N_get_array_3d_d_value(pot, i, j, k);
                    p2 = N_get_array_3d_d_value(pot, i, j, k + 1);
                    grad = (p1 - p2) / dz; /* gradient */
                }
                if (!N_is_array_3d_value_null(weight_z, i, j, k) &&
                    !N_is_array_3d_value_null(weight_z, i, j, k + 1)) {
                    r1 = N_get_array_3d_d_value(weight_z, i, j, k);
                    r2 = N_get_array_3d_d_value(weight_z, i, j, k + 1);
                    mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
                }
 
                res = mean * grad;
 
                G_debug(6,
                        "N_compute_gradient_field_3d: Z-direction insert value "
                        "%6.5g at %i %i %i ",
                        res, k + 1, j, i);
 
                N_put_array_3d_d_value(field->z_array, i, j, k + 1, res);
            }
 
    /*Compute gradient field statistics */
    N_calc_gradient_field_3d_stats(field);
 
    return field;
}
 
/*!
 * \brief Calculate the x, y and z vector components from a gradient field for
 each cell
 *  and store them in the provided N_array_3d structures
 *
 * The arrays must have the same size as the gradient field.
 *
 \verbatim
 
 Based on this storages scheme the gradient vector for each cell is
 calculated and stored in the provided  N_array_3d structures
 
 
 |  /
 TC NC
 |/
 --WC-----EC--
 /|
 SC BC
 /  |
 
 
 x vector component:
 
 x = (WC + EC) / 2
 
 y vector component:
 
 y = (NC + SC) / 2
 
 z vector component:
 
 z = (TC + BC) / 2
 
 \endverbatim
 
 * \param field N_gradient_field_3d *
 * \param x_comp N_array_3d * - the array in which the x component will be
 written
 * \param y_comp N_array_3d * - the array in which the y component will be
 written
 * \param z_comp N_array_3d * - the array in which the z component will be
 written
 *
 * \return void
 * */
void N_compute_gradient_field_components_3d(N_gradient_field_3d *field,
                                            N_array_3d *x_comp,
                                            N_array_3d *y_comp,
                                            N_array_3d *z_comp)
{
    int i, j, k;
 
    int rows, cols, depths;
 
    double vx, vy, vz;
 
    N_array_3d *x = x_comp;
 
    N_array_3d *y = y_comp;
 
    N_array_3d *z = z_comp;
 
    N_gradient_3d grad;
 
    if (!x)
        G_fatal_error("N_compute_gradient_components_3d: x array is empty");
    if (!y)
        G_fatal_error("N_compute_gradient_components_3d: y array is empty");
    if (!z)
        G_fatal_error("N_compute_gradient_components_3d: z array is empty");
 
    cols = field->x_array->cols;
    rows = field->x_array->rows;
    depths = field->x_array->depths;
 
    /*Check the array sizes */
    if (x->cols != cols || x->rows != rows || x->depths != depths)
        G_fatal_error("N_compute_gradient_components_3d: the size of the x "
                      "array doesn't fit the gradient field size");
    if (y->cols != cols || y->rows != rows || y->depths != depths)
        G_fatal_error("N_compute_gradient_components_3d: the size of the y "
                      "array doesn't fit the gradient field size");
    if (z->cols != cols || z->rows != rows || z->depths != depths)
        G_fatal_error("N_compute_gradient_components_3d: the size of the z "
                      "array doesn't fit the gradient field size");
 
    for (k = 0; k < depths; k++)
        for (j = 0; j < rows; j++)
            for (i = 0; i < cols; i++) {
                N_get_gradient_3d(field, &grad, i, j, k);
                /* in case a gradient is zero, we expect a no flow boundary */
                if (grad.WC == 0.0 || grad.EC == 0.0)
                    vx = (grad.WC + grad.EC);
                else
                    vx = (grad.WC + grad.EC) / 2;
                if (grad.NC == 0.0 || grad.SC == 0.0)
                    vy = (grad.NC + grad.SC);
                else
                    vy = (grad.NC + grad.SC) / 2;
                if (grad.TC == 0.0 || grad.BC == 0.0)
                    vz = (grad.TC + grad.BC);
                else
                    vz = (grad.TC + grad.BC) / 2;
 
                N_put_array_3d_d_value(x, i, j, k, vx);
                N_put_array_3d_d_value(y, i, j, k, vy);
                N_put_array_3d_d_value(z, i, j, k, vz);
            }
 
    return;
}