GRASS Programmer's Manual  6.4.2(2012)
N_gradient_calc.c
Go to the documentation of this file.
00001 
00002 /*****************************************************************************
00003 *
00004 * MODULE:       Grass PDE Numerical Library
00005 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
00006 *               soerengebbert <at> gmx <dot> de
00007 *               
00008 * PURPOSE:      gradient management functions 
00009 *               part of the gpde library
00010 *
00011 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
00012 *
00013 *               This program is free software under the GNU General Public
00014 *               License (>=v2). Read the file COPYING that comes with GRASS
00015 *               for details.
00016 *
00017 *****************************************************************************/
00018 
00019 #include <grass/N_pde.h>
00020 
00029 void N_calc_gradient_field_2d_stats(N_gradient_field_2d * field)
00030 {
00031     double minx, miny;
00032     double maxx, maxy;
00033     double sumx, sumy;
00034     int nonullx, nonully;
00035 
00036     G_debug(3,
00037             "N_calc_gradient_field_2d_stats: compute gradient field stats");
00038 
00039     N_calc_array_2d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
00040     N_calc_array_2d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
00041 
00042     if (minx < miny)
00043         field->min = minx;
00044     else
00045         field->min = miny;
00046 
00047     if (maxx > maxy)
00048         field->max = maxx;
00049     else
00050         field->max = maxy;
00051 
00052     field->sum = sumx + sumy;
00053     field->nonull = nonullx + nonully;
00054     field->mean = field->sum / (double)field->nonull;
00055 
00056     return;
00057 }
00058 
00105 N_gradient_field_2d *N_compute_gradient_field_2d(N_array_2d * pot,
00106                                                  N_array_2d * weight_x,
00107                                                  N_array_2d * weight_y,
00108                                                  N_geom_data * geom,
00109                                                  N_gradient_field_2d *
00110                                                  gradfield)
00111 {
00112     int i, j;
00113     int rows, cols;
00114     double dx, dy, p1, p2, r1, r2, mean, grad, res;
00115     N_gradient_field_2d *field = gradfield;
00116 
00117 
00118     if (pot->cols != weight_x->cols || pot->cols != weight_y->cols)
00119         G_fatal_error
00120             ("N_compute_gradient_field_2d: the arrays are not of equal size");
00121 
00122     if (pot->rows != weight_x->rows || pot->rows != weight_y->rows)
00123         G_fatal_error
00124             ("N_compute_gradient_field_2d: the arrays are not of equal size");
00125 
00126     if (pot->cols != geom->cols || pot->rows != geom->rows)
00127         G_fatal_error
00128             ("N_compute_gradient_field_2d: array sizes and geometry data are different");
00129 
00130 
00131     G_debug(3, "N_compute_gradient_field_2d: compute gradient field");
00132 
00133     rows = pot->rows;
00134     cols = pot->cols;
00135     dx = geom->dx;
00136     dy = geom->dy;
00137 
00138     if (field == NULL) {
00139         field = N_alloc_gradient_field_2d(cols, rows);
00140     }
00141     else {
00142         if (field->cols != geom->cols || field->rows != geom->rows)
00143             G_fatal_error
00144                 ("N_compute_gradient_field_2d: gradient field sizes and geometry data are different");
00145     }
00146 
00147 
00148     for (j = 0; j < rows; j++)
00149         for (i = 0; i < cols - 1; i++) {
00150             grad = 0;
00151             mean = 0;
00152 
00153             /* Only compute if the arrays are not null */
00154             if (!N_is_array_2d_value_null(pot, i, j) &&
00155                 !N_is_array_2d_value_null(pot, i + 1, j)) {
00156                 p1 = N_get_array_2d_d_value(pot, i, j);
00157                 p2 = N_get_array_2d_d_value(pot, i + 1, j);
00158                 grad = (p1 - p2) / dx;  /* gradient */
00159             }
00160             if (!N_is_array_2d_value_null(weight_x, i, j) &&
00161                 !N_is_array_2d_value_null(weight_x, i + 1, j)) {
00162                 r1 = N_get_array_2d_d_value(weight_x, i, j);
00163                 r2 = N_get_array_2d_d_value(weight_x, i + 1, j);
00164                 mean = N_calc_harmonic_mean(r1, r2);    /*harmonical mean */
00165             }
00166 
00167             res = mean * grad;
00168 
00169             N_put_array_2d_d_value(field->x_array, i + 1, j, res);
00170 
00171         }
00172 
00173     for (j = 0; j < rows - 1; j++)
00174         for (i = 0; i < cols; i++) {
00175             grad = 0;
00176             mean = 0;
00177 
00178             /* Only compute if the arrays are not null */
00179             if (!N_is_array_2d_value_null(pot, i, j) &&
00180                 !N_is_array_2d_value_null(pot, i, j + 1)) {
00181                 p1 = N_get_array_2d_d_value(pot, i, j);
00182                 p2 = N_get_array_2d_d_value(pot, i, j + 1);
00183                 grad = (p1 - p2) / dy;  /* gradient */
00184             }
00185             if (!N_is_array_2d_value_null(weight_y, i, j) &&
00186                 !N_is_array_2d_value_null(weight_y, i, j + 1)) {
00187                 r1 = N_get_array_2d_d_value(weight_y, i, j);
00188                 r2 = N_get_array_2d_d_value(weight_y, i, j + 1);
00189                 mean = N_calc_harmonic_mean(r1, r2);    /*harmonical mean */
00190             }
00191 
00192             res = -1 * mean * grad;
00193 
00194             N_put_array_2d_d_value(field->y_array, i, j + 1, res);
00195 
00196         }
00197 
00198     /*Compute gradient field statistics */
00199     N_calc_gradient_field_2d_stats(field);
00200 
00201     return field;
00202 }
00203 
00242 void
00243 N_compute_gradient_field_components_2d(N_gradient_field_2d * field,
00244                                        N_array_2d * x_comp,
00245                                        N_array_2d * y_comp)
00246 {
00247     int i, j;
00248     int rows, cols;
00249     double vx, vy;
00250     N_array_2d *x = x_comp;
00251     N_array_2d *y = y_comp;
00252     N_gradient_2d grad;
00253 
00254 
00255     if (!x)
00256         G_fatal_error("N_compute_gradient_components_2d: x array is empty");
00257     if (!y)
00258         G_fatal_error("N_compute_gradient_components_2d: y array is empty");
00259 
00260     cols = field->x_array->cols;
00261     rows = field->x_array->rows;
00262 
00263     /*Check the array sizes */
00264     if (x->cols != cols || x->rows != rows)
00265         G_fatal_error
00266             ("N_compute_gradient_components_2d: the size of the x array doesn't fit the gradient field size");
00267     if (y->cols != cols || y->rows != rows)
00268         G_fatal_error
00269             ("N_compute_gradient_components_2d: the size of the y array doesn't fit the gradient field size");
00270 
00271     for (j = 0; j < rows; j++)
00272         for (i = 0; i < cols; i++) {
00273             N_get_gradient_2d(field, &grad, i, j);
00274 
00275             /* in case a gradient is zero, we expect a no flow boundary */
00276             if (grad.WC == 0.0 || grad.EC == 0.0)
00277                 vx = (grad.WC + grad.EC);
00278             else
00279                 vx = (grad.WC + grad.EC) / 2;
00280             if (grad.NC == 0.0 || grad.SC == 0.0)
00281                 vy = (grad.NC + grad.SC);
00282             else
00283                 vy = (grad.NC + grad.SC) / 2;
00284 
00285             N_put_array_2d_d_value(x, i, j, vx);
00286             N_put_array_2d_d_value(y, i, j, vy);
00287         }
00288 
00289     return;
00290 }
00291 
00300 void N_calc_gradient_field_3d_stats(N_gradient_field_3d * field)
00301 {
00302     double minx, miny, minz;
00303     double maxx, maxy, maxz;
00304     double sumx, sumy, sumz;
00305     int nonullx, nonully, nonullz;
00306 
00307     G_debug(3,
00308             "N_calc_gradient_field_3d_stats: compute gradient field stats");
00309 
00310     N_calc_array_3d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
00311     N_calc_array_3d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
00312     N_calc_array_3d_stats(field->z_array, &minz, &maxz, &sumz, &nonullz, 0);
00313 
00314     if (minx <= minz && minx <= miny)
00315         field->min = minx;
00316     if (miny <= minz && miny <= minx)
00317         field->min = miny;
00318     if (minz <= minx && minz <= miny)
00319         field->min = minz;
00320 
00321     if (maxx >= maxz && maxx >= maxy)
00322         field->max = maxx;
00323     if (maxy >= maxz && maxy >= maxx)
00324         field->max = maxy;
00325     if (maxz >= maxx && maxz >= maxy)
00326         field->max = maxz;
00327 
00328     field->sum = sumx + sumy + sumz;
00329     field->nonull = nonullx + nonully + nonullz;
00330     field->mean = field->sum / (double)field->nonull;
00331 
00332     return;
00333 }
00334 
00335 
00386 N_gradient_field_3d *N_compute_gradient_field_3d(N_array_3d * pot,
00387                                                  N_array_3d * weight_x,
00388                                                  N_array_3d * weight_y,
00389                                                  N_array_3d * weight_z,
00390                                                  N_geom_data * geom,
00391                                                  N_gradient_field_3d *
00392                                                  gradfield)
00393 {
00394     int i, j, k;
00395     int cols, rows, depths;
00396     double dx, dy, dz, p1, p2, r1, r2, mean, grad, res;
00397     N_gradient_field_3d *field = gradfield;
00398 
00399 
00400     if (pot->cols != weight_x->cols || pot->cols != weight_y->cols ||
00401         pot->cols != weight_z->cols)
00402         G_fatal_error
00403             ("N_compute_gradient_field_3d: the arrays are not of equal size");
00404 
00405     if (pot->rows != weight_x->rows || pot->rows != weight_y->rows ||
00406         pot->rows != weight_z->rows)
00407         G_fatal_error
00408             ("N_compute_gradient_field_3d: the arrays are not of equal size");
00409 
00410     if (pot->depths != weight_x->depths || pot->depths != weight_y->depths ||
00411         pot->depths != weight_z->depths)
00412         G_fatal_error
00413             ("N_compute_gradient_field_3d: the arrays are not of equal size");
00414 
00415     if (pot->cols != geom->cols || pot->rows != geom->rows ||
00416         pot->depths != geom->depths)
00417         G_fatal_error
00418             ("N_compute_gradient_field_3d: array sizes and geometry data are different");
00419 
00420     G_debug(3, "N_compute_gradient_field_3d: compute gradient field");
00421 
00422     cols = geom->cols;
00423     rows = geom->rows;
00424     depths = geom->depths;
00425     dx = geom->dx;
00426     dy = geom->dy;
00427     dz = geom->dz;
00428 
00429     if (gradfield == NULL) {
00430         field = N_alloc_gradient_field_3d(cols, rows, depths);
00431     }
00432     else {
00433         if (field->cols != geom->cols || field->rows != geom->rows ||
00434             field->depths != geom->depths)
00435             G_fatal_error
00436                 ("N_compute_gradient_field_3d: gradient field sizes and geometry data are different");
00437     }
00438 
00439     for (k = 0; k < depths; k++)
00440         for (j = 0; j < rows; j++)
00441             for (i = 0; i < cols - 1; i++) {
00442                 grad = 0;
00443                 mean = 0;
00444 
00445                 /*Only compute if the arrays are not null */
00446                 if (!N_is_array_3d_value_null(pot, i, j, k) &&
00447                     !N_is_array_3d_value_null(pot, i + 1, j, k)) {
00448                     p1 = N_get_array_3d_d_value(pot, i, j, k);
00449                     p2 = N_get_array_3d_d_value(pot, i + 1, j, k);
00450                     grad = (p1 - p2) / dx;      /* gradient */
00451                 }
00452                 if (!N_is_array_3d_value_null(weight_x, i, j, k) &&
00453                     !N_is_array_3d_value_null(weight_x, i + 1, j, k)) {
00454                     r1 = N_get_array_3d_d_value(weight_x, i, j, k);
00455                     r2 = N_get_array_3d_d_value(weight_x, i + 1, j, k);
00456                     mean = N_calc_harmonic_mean(r1, r2);        /*harmonical mean */
00457                 }
00458 
00459                 res = mean * grad;
00460 
00461                 G_debug(6,
00462                         "N_compute_gradient_field_3d: X-direction insert value %6.5g at %i %i %i ",
00463                         res, k, j, i + 1);
00464 
00465                 N_put_array_3d_d_value(field->x_array, i + 1, j, k, res);
00466 
00467             }
00468 
00469     for (k = 0; k < depths; k++)
00470         for (j = 0; j < rows - 1; j++)
00471             for (i = 0; i < cols; i++) {
00472                 grad = 0;
00473                 mean = 0;
00474 
00475                 /* Only compute if the arrays are not null */
00476                 if (!N_is_array_3d_value_null(pot, i, j, k) &&
00477                     !N_is_array_3d_value_null(pot, i, j + 1, k)) {
00478                     p1 = N_get_array_3d_d_value(pot, i, j, k);
00479                     p2 = N_get_array_3d_d_value(pot, i, j + 1, k);
00480                     grad = (p1 - p2) / dy;      /* gradient */
00481                 }
00482                 if (!N_is_array_3d_value_null(weight_y, i, j, k) &&
00483                     !N_is_array_3d_value_null(weight_y, i, j + 1, k)) {
00484                     r1 = N_get_array_3d_d_value(weight_y, i, j, k);
00485                     r2 = N_get_array_3d_d_value(weight_y, i, j + 1, k);
00486                     mean = N_calc_harmonic_mean(r1, r2);        /*harmonical mean */
00487                 }
00488 
00489                 res = -1 * mean * grad; /*invert the direction, because we count from north to south,
00490                                          * but the gradient is defined in y direction */
00491 
00492                 G_debug(6,
00493                         "N_compute_gradient_field_3d: Y-direction insert value %6.5g at %i %i %i ",
00494                         res, k, j + 1, i);
00495 
00496                 N_put_array_3d_d_value(field->y_array, i, j + 1, k, res);
00497 
00498             }
00499 
00500     for (k = 0; k < depths - 1; k++)
00501         for (j = 0; j < rows; j++)
00502             for (i = 0; i < cols; i++) {
00503                 grad = 0;
00504                 mean = 0;
00505 
00506                 /* Only compute if the arrays are not null */
00507                 if (!N_is_array_3d_value_null(pot, i, j, k) &&
00508                     !N_is_array_3d_value_null(pot, i, j, k + 1)) {
00509                     p1 = N_get_array_3d_d_value(pot, i, j, k);
00510                     p2 = N_get_array_3d_d_value(pot, i, j, k + 1);
00511                     grad = (p1 - p2) / dz;      /* gradient */
00512                 }
00513                 if (!N_is_array_3d_value_null(weight_z, i, j, k) &&
00514                     !N_is_array_3d_value_null(weight_z, i, j, k + 1)) {
00515                     r1 = N_get_array_3d_d_value(weight_z, i, j, k);
00516                     r2 = N_get_array_3d_d_value(weight_z, i, j, k + 1);
00517                     mean = N_calc_harmonic_mean(r1, r2);        /*harmonical mean */
00518                 }
00519 
00520                 res = mean * grad;
00521 
00522                 G_debug(6,
00523                         "N_compute_gradient_field_3d: Z-direction insert value %6.5g at %i %i %i ",
00524                         res, k + 1, j, i);
00525 
00526                 N_put_array_3d_d_value(field->z_array, i, j, k + 1, res);
00527 
00528             }
00529 
00530     /*Compute gradient field statistics */
00531     N_calc_gradient_field_3d_stats(field);
00532 
00533     return field;
00534 }
00535 
00578 void
00579 N_compute_gradient_field_components_3d(N_gradient_field_3d * field,
00580                                        N_array_3d * x_comp,
00581                                        N_array_3d * y_comp,
00582                                        N_array_3d * z_comp)
00583 {
00584     int i, j, k;
00585     int rows, cols, depths;
00586     double vx, vy, vz;
00587     N_array_3d *x = x_comp;
00588     N_array_3d *y = y_comp;
00589     N_array_3d *z = z_comp;
00590     N_gradient_3d grad;
00591 
00592 
00593     if (!x)
00594         G_fatal_error("N_compute_gradient_components_3d: x array is empty");
00595     if (!y)
00596         G_fatal_error("N_compute_gradient_components_3d: y array is empty");
00597     if (!z)
00598         G_fatal_error("N_compute_gradient_components_3d: z array is empty");
00599 
00600     cols = field->x_array->cols;
00601     rows = field->x_array->rows;
00602     depths = field->x_array->depths;
00603 
00604     /*Check the array sizes */
00605     if (x->cols != cols || x->rows != rows || x->depths != depths)
00606         G_fatal_error
00607             ("N_compute_gradient_components_3d: the size of the x array doesn't fit the gradient field size");
00608     if (y->cols != cols || y->rows != rows || y->depths != depths)
00609         G_fatal_error
00610             ("N_compute_gradient_components_3d: the size of the y array doesn't fit the gradient field size");
00611     if (z->cols != cols || z->rows != rows || z->depths != depths)
00612         G_fatal_error
00613             ("N_compute_gradient_components_3d: the size of the z array doesn't fit the gradient field size");
00614 
00615     for (k = 0; k < depths; k++)
00616         for (j = 0; j < rows; j++)
00617             for (i = 0; i < cols; i++) {
00618                 N_get_gradient_3d(field, &grad, i, j, k);
00619                 /* in case a gradient is zero, we expect a no flow boundary */
00620                 if (grad.WC == 0.0 || grad.EC == 0.0)
00621                     vx = (grad.WC + grad.EC);
00622                 else
00623                     vx = (grad.WC + grad.EC) / 2;
00624                 if (grad.NC == 0.0 || grad.SC == 0.0)
00625                     vy = (grad.NC + grad.SC);
00626                 else
00627                     vy = (grad.NC + grad.SC) / 2;
00628                 if (grad.TC == 0.0 || grad.BC == 0.0)
00629                     vz = (grad.TC + grad.BC);
00630                 else
00631                     vz = (grad.TC + grad.BC) / 2;
00632 
00633                 N_put_array_3d_d_value(x, i, j, k, vx);
00634                 N_put_array_3d_d_value(y, i, j, k, vy);
00635                 N_put_array_3d_d_value(z, i, j, k, vz);
00636             }
00637 
00638 
00639     return;
00640 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines