GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
N_gradient_calc.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: gradient management functions
9 * part of the gpde library
10 *
11 * COPYRIGHT: (C) 2000 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18 
19 #include <grass/N_pde.h>
20 
30 {
31  double minx, miny;
32  double maxx, maxy;
33  double sumx, sumy;
34  int nonullx, nonully;
35 
36  G_debug(3,
37  "N_calc_gradient_field_2d_stats: compute gradient field stats");
38 
39  N_calc_array_2d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
40  N_calc_array_2d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
41 
42  if (minx < miny)
43  field->min = minx;
44  else
45  field->min = miny;
46 
47  if (maxx > maxy)
48  field->max = maxx;
49  else
50  field->max = maxy;
51 
52  field->sum = sumx + sumy;
53  field->nonull = nonullx + nonully;
54  field->mean = field->sum / (double)field->nonull;
55 
56  return;
57 }
58 
106  N_array_2d * weight_x,
107  N_array_2d * weight_y,
108  N_geom_data * geom,
110  gradfield)
111 {
112  int i, j;
113  int rows, cols;
114  double dx, dy, p1, p2, r1, r2, mean, grad, res;
115  N_gradient_field_2d *field = gradfield;
116 
117 
118  if (pot->cols != weight_x->cols || pot->cols != weight_y->cols)
120  ("N_compute_gradient_field_2d: the arrays are not of equal size");
121 
122  if (pot->rows != weight_x->rows || pot->rows != weight_y->rows)
124  ("N_compute_gradient_field_2d: the arrays are not of equal size");
125 
126  if (pot->cols != geom->cols || pot->rows != geom->rows)
128  ("N_compute_gradient_field_2d: array sizes and geometry data are different");
129 
130 
131  G_debug(3, "N_compute_gradient_field_2d: compute gradient field");
132 
133  rows = pot->rows;
134  cols = pot->cols;
135  dx = geom->dx;
136  dy = geom->dy;
137 
138  if (field == NULL) {
139  field = N_alloc_gradient_field_2d(cols, rows);
140  }
141  else {
142  if (field->cols != geom->cols || field->rows != geom->rows)
144  ("N_compute_gradient_field_2d: gradient field sizes and geometry data are different");
145  }
146 
147 
148  for (j = 0; j < rows; j++)
149  for (i = 0; i < cols - 1; i++) {
150  grad = 0;
151  mean = 0;
152 
153  /* Only compute if the arrays are not null */
154  if (!N_is_array_2d_value_null(pot, i, j) &&
155  !N_is_array_2d_value_null(pot, i + 1, j)) {
156  p1 = N_get_array_2d_d_value(pot, i, j);
157  p2 = N_get_array_2d_d_value(pot, i + 1, j);
158  grad = (p1 - p2) / dx; /* gradient */
159  }
160  if (!N_is_array_2d_value_null(weight_x, i, j) &&
161  !N_is_array_2d_value_null(weight_x, i + 1, j)) {
162  r1 = N_get_array_2d_d_value(weight_x, i, j);
163  r2 = N_get_array_2d_d_value(weight_x, i + 1, j);
164  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
165  }
166 
167  res = mean * grad;
168 
169  N_put_array_2d_d_value(field->x_array, i + 1, j, res);
170 
171  }
172 
173  for (j = 0; j < rows - 1; j++)
174  for (i = 0; i < cols; i++) {
175  grad = 0;
176  mean = 0;
177 
178  /* Only compute if the arrays are not null */
179  if (!N_is_array_2d_value_null(pot, i, j) &&
180  !N_is_array_2d_value_null(pot, i, j + 1)) {
181  p1 = N_get_array_2d_d_value(pot, i, j);
182  p2 = N_get_array_2d_d_value(pot, i, j + 1);
183  grad = (p1 - p2) / dy; /* gradient */
184  }
185  if (!N_is_array_2d_value_null(weight_y, i, j) &&
186  !N_is_array_2d_value_null(weight_y, i, j + 1)) {
187  r1 = N_get_array_2d_d_value(weight_y, i, j);
188  r2 = N_get_array_2d_d_value(weight_y, i, j + 1);
189  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
190  }
191 
192  res = -1 * mean * grad;
193 
194  N_put_array_2d_d_value(field->y_array, i, j + 1, res);
195 
196  }
197 
198  /*Compute gradient field statistics */
200 
201  return field;
202 }
203 
242 void
244  N_array_2d * x_comp,
245  N_array_2d * y_comp)
246 {
247  int i, j;
248  int rows, cols;
249  double vx, vy;
250  N_array_2d *x = x_comp;
251  N_array_2d *y = y_comp;
252  N_gradient_2d grad;
253 
254 
255  if (!x)
256  G_fatal_error("N_compute_gradient_components_2d: x array is empty");
257  if (!y)
258  G_fatal_error("N_compute_gradient_components_2d: y array is empty");
259 
260  cols = field->x_array->cols;
261  rows = field->x_array->rows;
262 
263  /*Check the array sizes */
264  if (x->cols != cols || x->rows != rows)
266  ("N_compute_gradient_components_2d: the size of the x array doesn't fit the gradient field size");
267  if (y->cols != cols || y->rows != rows)
269  ("N_compute_gradient_components_2d: the size of the y array doesn't fit the gradient field size");
270 
271  for (j = 0; j < rows; j++)
272  for (i = 0; i < cols; i++) {
273  N_get_gradient_2d(field, &grad, i, j);
274 
275  /* in case a gradient is zero, we expect a no flow boundary */
276  if (grad.WC == 0.0 || grad.EC == 0.0)
277  vx = (grad.WC + grad.EC);
278  else
279  vx = (grad.WC + grad.EC) / 2;
280  if (grad.NC == 0.0 || grad.SC == 0.0)
281  vy = (grad.NC + grad.SC);
282  else
283  vy = (grad.NC + grad.SC) / 2;
284 
285  N_put_array_2d_d_value(x, i, j, vx);
286  N_put_array_2d_d_value(y, i, j, vy);
287  }
288 
289  return;
290 }
291 
301 {
302  double minx, miny, minz;
303  double maxx, maxy, maxz;
304  double sumx, sumy, sumz;
305  int nonullx, nonully, nonullz;
306 
307  G_debug(3,
308  "N_calc_gradient_field_3d_stats: compute gradient field stats");
309 
310  N_calc_array_3d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
311  N_calc_array_3d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
312  N_calc_array_3d_stats(field->z_array, &minz, &maxz, &sumz, &nonullz, 0);
313 
314  if (minx <= minz && minx <= miny)
315  field->min = minx;
316  if (miny <= minz && miny <= minx)
317  field->min = miny;
318  if (minz <= minx && minz <= miny)
319  field->min = minz;
320 
321  if (maxx >= maxz && maxx >= maxy)
322  field->max = maxx;
323  if (maxy >= maxz && maxy >= maxx)
324  field->max = maxy;
325  if (maxz >= maxx && maxz >= maxy)
326  field->max = maxz;
327 
328  field->sum = sumx + sumy + sumz;
329  field->nonull = nonullx + nonully + nonullz;
330  field->mean = field->sum / (double)field->nonull;
331 
332  return;
333 }
334 
335 
387  N_array_3d * weight_x,
388  N_array_3d * weight_y,
389  N_array_3d * weight_z,
390  N_geom_data * geom,
392  gradfield)
393 {
394  int i, j, k;
395  int cols, rows, depths;
396  double dx, dy, dz, p1, p2, r1, r2, mean, grad, res;
397  N_gradient_field_3d *field = gradfield;
398 
399 
400  if (pot->cols != weight_x->cols || pot->cols != weight_y->cols ||
401  pot->cols != weight_z->cols)
403  ("N_compute_gradient_field_3d: the arrays are not of equal size");
404 
405  if (pot->rows != weight_x->rows || pot->rows != weight_y->rows ||
406  pot->rows != weight_z->rows)
408  ("N_compute_gradient_field_3d: the arrays are not of equal size");
409 
410  if (pot->depths != weight_x->depths || pot->depths != weight_y->depths ||
411  pot->depths != weight_z->depths)
413  ("N_compute_gradient_field_3d: the arrays are not of equal size");
414 
415  if (pot->cols != geom->cols || pot->rows != geom->rows ||
416  pot->depths != geom->depths)
418  ("N_compute_gradient_field_3d: array sizes and geometry data are different");
419 
420  G_debug(3, "N_compute_gradient_field_3d: compute gradient field");
421 
422  cols = geom->cols;
423  rows = geom->rows;
424  depths = geom->depths;
425  dx = geom->dx;
426  dy = geom->dy;
427  dz = geom->dz;
428 
429  if (gradfield == NULL) {
430  field = N_alloc_gradient_field_3d(cols, rows, depths);
431  }
432  else {
433  if (field->cols != geom->cols || field->rows != geom->rows ||
434  field->depths != geom->depths)
436  ("N_compute_gradient_field_3d: gradient field sizes and geometry data are different");
437  }
438 
439  for (k = 0; k < depths; k++)
440  for (j = 0; j < rows; j++)
441  for (i = 0; i < cols - 1; i++) {
442  grad = 0;
443  mean = 0;
444 
445  /*Only compute if the arrays are not null */
446  if (!N_is_array_3d_value_null(pot, i, j, k) &&
447  !N_is_array_3d_value_null(pot, i + 1, j, k)) {
448  p1 = N_get_array_3d_d_value(pot, i, j, k);
449  p2 = N_get_array_3d_d_value(pot, i + 1, j, k);
450  grad = (p1 - p2) / dx; /* gradient */
451  }
452  if (!N_is_array_3d_value_null(weight_x, i, j, k) &&
453  !N_is_array_3d_value_null(weight_x, i + 1, j, k)) {
454  r1 = N_get_array_3d_d_value(weight_x, i, j, k);
455  r2 = N_get_array_3d_d_value(weight_x, i + 1, j, k);
456  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
457  }
458 
459  res = mean * grad;
460 
461  G_debug(6,
462  "N_compute_gradient_field_3d: X-direction insert value %6.5g at %i %i %i ",
463  res, k, j, i + 1);
464 
465  N_put_array_3d_d_value(field->x_array, i + 1, j, k, res);
466 
467  }
468 
469  for (k = 0; k < depths; k++)
470  for (j = 0; j < rows - 1; j++)
471  for (i = 0; i < cols; i++) {
472  grad = 0;
473  mean = 0;
474 
475  /* Only compute if the arrays are not null */
476  if (!N_is_array_3d_value_null(pot, i, j, k) &&
477  !N_is_array_3d_value_null(pot, i, j + 1, k)) {
478  p1 = N_get_array_3d_d_value(pot, i, j, k);
479  p2 = N_get_array_3d_d_value(pot, i, j + 1, k);
480  grad = (p1 - p2) / dy; /* gradient */
481  }
482  if (!N_is_array_3d_value_null(weight_y, i, j, k) &&
483  !N_is_array_3d_value_null(weight_y, i, j + 1, k)) {
484  r1 = N_get_array_3d_d_value(weight_y, i, j, k);
485  r2 = N_get_array_3d_d_value(weight_y, i, j + 1, k);
486  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
487  }
488 
489  res = -1 * mean * grad; /*invert the direction, because we count from north to south,
490  * but the gradient is defined in y direction */
491 
492  G_debug(6,
493  "N_compute_gradient_field_3d: Y-direction insert value %6.5g at %i %i %i ",
494  res, k, j + 1, i);
495 
496  N_put_array_3d_d_value(field->y_array, i, j + 1, k, res);
497 
498  }
499 
500  for (k = 0; k < depths - 1; k++)
501  for (j = 0; j < rows; j++)
502  for (i = 0; i < cols; i++) {
503  grad = 0;
504  mean = 0;
505 
506  /* Only compute if the arrays are not null */
507  if (!N_is_array_3d_value_null(pot, i, j, k) &&
508  !N_is_array_3d_value_null(pot, i, j, k + 1)) {
509  p1 = N_get_array_3d_d_value(pot, i, j, k);
510  p2 = N_get_array_3d_d_value(pot, i, j, k + 1);
511  grad = (p1 - p2) / dz; /* gradient */
512  }
513  if (!N_is_array_3d_value_null(weight_z, i, j, k) &&
514  !N_is_array_3d_value_null(weight_z, i, j, k + 1)) {
515  r1 = N_get_array_3d_d_value(weight_z, i, j, k);
516  r2 = N_get_array_3d_d_value(weight_z, i, j, k + 1);
517  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
518  }
519 
520  res = mean * grad;
521 
522  G_debug(6,
523  "N_compute_gradient_field_3d: Z-direction insert value %6.5g at %i %i %i ",
524  res, k + 1, j, i);
525 
526  N_put_array_3d_d_value(field->z_array, i, j, k + 1, res);
527 
528  }
529 
530  /*Compute gradient field statistics */
532 
533  return field;
534 }
535 
578 void
580  N_array_3d * x_comp,
581  N_array_3d * y_comp,
582  N_array_3d * z_comp)
583 {
584  int i, j, k;
585  int rows, cols, depths;
586  double vx, vy, vz;
587  N_array_3d *x = x_comp;
588  N_array_3d *y = y_comp;
589  N_array_3d *z = z_comp;
590  N_gradient_3d grad;
591 
592 
593  if (!x)
594  G_fatal_error("N_compute_gradient_components_3d: x array is empty");
595  if (!y)
596  G_fatal_error("N_compute_gradient_components_3d: y array is empty");
597  if (!z)
598  G_fatal_error("N_compute_gradient_components_3d: z array is empty");
599 
600  cols = field->x_array->cols;
601  rows = field->x_array->rows;
602  depths = field->x_array->depths;
603 
604  /*Check the array sizes */
605  if (x->cols != cols || x->rows != rows || x->depths != depths)
607  ("N_compute_gradient_components_3d: the size of the x array doesn't fit the gradient field size");
608  if (y->cols != cols || y->rows != rows || y->depths != depths)
610  ("N_compute_gradient_components_3d: the size of the y array doesn't fit the gradient field size");
611  if (z->cols != cols || z->rows != rows || z->depths != depths)
613  ("N_compute_gradient_components_3d: the size of the z array doesn't fit the gradient field size");
614 
615  for (k = 0; k < depths; k++)
616  for (j = 0; j < rows; j++)
617  for (i = 0; i < cols; i++) {
618  N_get_gradient_3d(field, &grad, i, j, k);
619  /* in case a gradient is zero, we expect a no flow boundary */
620  if (grad.WC == 0.0 || grad.EC == 0.0)
621  vx = (grad.WC + grad.EC);
622  else
623  vx = (grad.WC + grad.EC) / 2;
624  if (grad.NC == 0.0 || grad.SC == 0.0)
625  vy = (grad.NC + grad.SC);
626  else
627  vy = (grad.NC + grad.SC) / 2;
628  if (grad.TC == 0.0 || grad.BC == 0.0)
629  vz = (grad.TC + grad.BC);
630  else
631  vz = (grad.TC + grad.BC) / 2;
632 
633  N_put_array_3d_d_value(x, i, j, k, vx);
634  N_put_array_3d_d_value(y, i, j, k, vy);
635  N_put_array_3d_d_value(z, i, j, k, vz);
636  }
637 
638 
639  return;
640 }
double EC
Definition: N_pde.h:520
int N_is_array_3d_value_null(N_array_3d *data, int col, int row, int depth)
This function returns 1 if value of N_array_3d data at position col, row, depth is of type null...
Definition: N_arrays.c:878
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)
This function computes the gradient based on the input N_array_3d pot (that means potential)...
int rows
Definition: N_pde.h:181
void N_compute_gradient_field_components_2d(N_gradient_field_2d *field, N_array_2d *x_comp, N_array_2d *y_comp)
Calculate the x and y vector components from a gradient field for each cell and stores them in the pr...
Gradient between the cells in X and Y direction.
Definition: N_pde.h:509
int depths
Definition: N_pde.h:223
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
Definition: N_arrays.c:581
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)
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.
void N_calc_gradient_field_3d_stats(N_gradient_field_3d *field)
Calculate basic statistics of a gradient field.
int cols
Definition: N_pde.h:181
void N_calc_array_2d_stats(N_array_2d *a, double *min, double *max, double *sum, int *nonull, int withoffset)
Calculate basic statistics of the N_array_2d struct.
double WC
Definition: N_pde.h:520
int y
Definition: plot.c:34
int depths
Definition: N_pde.h:142
N_array_3d * y_array
Definition: N_pde.h:637
double SC
Definition: N_pde.h:520
double dy
Definition: N_pde.h:137
N_gradient_3d * N_get_gradient_3d(N_gradient_field_3d *field, N_gradient_3d *gradient, int col, int row, int depth)
Return a N_gradient_3d structure calculated from the input gradient field at position [depth][row][co...
Definition: N_gradient.c:248
double NC
Definition: N_pde.h:520
Geometric information about the structured grid.
Definition: N_pde.h:130
N_array_2d * y_array
Definition: N_pde.h:625
double dz
Definition: N_pde.h:138
N_gradient_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
Definition: N_gradient.c:920
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)
Calculate the x, y and z vector components from a gradient field for each cell and store them in the ...
int rows
Definition: N_pde.h:143
double TC
Definition: N_pde.h:520
N_array_3d * z_array
Definition: N_pde.h:638
int cols
Definition: N_pde.h:223
double BC
Definition: N_pde.h:520
double EC
Definition: N_pde.h:512
int cols
Definition: N_pde.h:144
Gradient between the cells in X, Y and Z direction.
Definition: N_pde.h:517
void N_calc_array_3d_stats(N_array_3d *a, double *min, double *max, double *sum, int *nonull, int withoffset)
Calculate basic statistics of the N_array_3d struct.
int N_is_array_2d_value_null(N_array_2d *data, int col, int row)
Returns 1 if the value of N_array_2d struct at postion col, row is of type null, otherwise 0...
Definition: N_arrays.c:228
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: N_arrays.c:987
N_gradient_2d * N_get_gradient_2d(N_gradient_field_2d *field, N_gradient_2d *gradient, int col, int row)
Return a N_gradient_2d structure calculated from the input gradient field at position [row][col]...
Definition: N_gradient.c:115
return NULL
Definition: dbfopen.c:1394
tuple cols
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
N_array_3d * x_array
Definition: N_pde.h:636
N_gradient_field_3d * N_alloc_gradient_field_3d(int cols, int rows, int depths)
Allocate a N_gradient_field_3d.
Definition: N_gradient.c:1018
double WC
Definition: N_pde.h:512
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition: N_arrays.c:375
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
Definition: N_arrays.c:1172
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition: N_tools.c:119
int G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
void N_calc_gradient_field_2d_stats(N_gradient_field_2d *field)
Calculate basic statistics of a gradient field.
double dx
Definition: N_pde.h:136
double SC
Definition: N_pde.h:512
double NC
Definition: N_pde.h:512
N_array_2d * x_array
Definition: N_pde.h:624
int rows
Definition: N_pde.h:223