GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
test_gradient.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: Unit tests for gradient calculation
9 *
10 * COPYRIGHT: (C) 2000 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General Public
13 * License (>=v2). Read the file COPYING that comes with GRASS
14 * for details.
15 *
16 *****************************************************************************/
17 
18 
19 
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <grass/glocale.h>
24 #include <grass/N_pde.h>
25 #include "test_gpde_lib.h"
26 
27 /* prototypes */
28 static int test_gradient_2d(void);
29 static int test_gradient_3d(void);
30 static N_array_2d *create_relax_array_2d(void);
31 static N_array_3d *create_relax_array_3d(void);
32 static N_array_2d *create_potential_array_2d(void);
33 static N_array_3d *create_potential_array_3d(void);
34 
35 /* *************************************************************** */
36 /* Performe the gradient tests *********************************** */
37 /* *************************************************************** */
39 {
40  int sum = 0;
41 
42  G_message(_("\n++ Running gradient unit tests ++"));
43 
44  G_message(_("\t 1. testing 2d gradient"));
45  sum += test_gradient_2d();
46 
47  G_message(_("\t 2. testing 3d gradient"));
48  sum += test_gradient_3d();
49 
50  if (sum > 0)
51  G_warning(_("\n-- Gradient unit tests failure --"));
52  else
53  G_message(_("\n-- Gradient unit tests finished successfully --"));
54 
55  return sum;
56 }
57 
58 /* *************************************************************** */
59 /* Create the status array with values of 1 and 2 **************** */
60 /* *************************************************************** */
61 N_array_2d *create_relax_array_2d(void)
62 {
64  int i, j;
65 
66  data = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, CELL_TYPE);
67 
68 #pragma omp parallel for private (i, j) shared (data)
69  for (j = 0; j < TEST_N_NUM_ROWS; j++) {
70  for (i = 0; i < TEST_N_NUM_COLS; i++) {
71  N_put_array_2d_c_value(data, i, j, 1);
72  }
73  }
74  return data;
75 }
76 
77 /* *************************************************************** */
78 /* Create a value array ****************************************** */
79 /* *************************************************************** */
80 N_array_2d *create_potential_array_2d(void)
81 {
83  int i, j;
84 
85  data = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, DCELL_TYPE);
86 
87 #pragma omp parallel for private (i, j) shared (data)
88  for (j = 0; j < TEST_N_NUM_ROWS; j++) {
89  for (i = 0; i < TEST_N_NUM_COLS; i++) {
90  N_put_array_2d_d_value(data, i, j, (double)i * j);
91  }
92  }
93  return data;
94 }
95 
96 /* *************************************************************** */
97 /* Create the status array with values of 1 and 2 **************** */
98 /* *************************************************************** */
99 N_array_3d *create_relax_array_3d(void)
100 {
101  N_array_3d *data;
102  int i, j, k;
103 
104  data =
106  1, FCELL_TYPE);
107 
108 #pragma omp parallel for private (i, j, k) shared (data)
109  for (k = 0; k < TEST_N_NUM_DEPTHS; k++) {
110  for (j = 0; j < TEST_N_NUM_ROWS; j++) {
111  for (i = 0; i < TEST_N_NUM_COLS; i++) {
112  N_put_array_3d_f_value(data, i, j, k, 1.0);
113  }
114  }
115  }
116 
117  return data;
118 
119 }
120 
121 /* *************************************************************** */
122 /* Create a value array ****************************************** */
123 /* *************************************************************** */
124 N_array_3d *create_potential_array_3d(void)
125 {
126  N_array_3d *data;
127  int i, j, k;
128 
129  data =
130  N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
131  1, DCELL_TYPE);
132 
133 #pragma omp parallel for private (i, j, k) shared (data)
134  for (k = 0; k < TEST_N_NUM_DEPTHS; k++)
135  for (j = 0; j < TEST_N_NUM_ROWS; j++) {
136  for (i = 0; i < TEST_N_NUM_COLS; i++) {
137  N_put_array_3d_f_value(data, i, j, k, (double)i * j * k);
138  }
139  }
140 
141  return data;
142 
143 }
144 
145 /* *************************************************************** */
146 /* Test the gradient calculation with 3d array data ************** */
147 /* *************************************************************** */
148 int test_gradient_3d(void)
149 {
150  N_array_3d *relax = NULL;
151  N_array_3d *pot = NULL;
152  N_array_3d *xcomp = NULL;
153  N_array_3d *ycomp = NULL;
154  N_array_3d *zcomp = NULL;
155  N_gradient_field_3d *field = NULL;
156  N_gradient_3d *grad = NULL;
157  N_geom_data *geom = NULL;
158 
159  geom = N_alloc_geom_data();
160 
161  geom->dx = 1;
162  geom->dy = 1;
163  geom->dz = 1;
164 
165  geom->Az = 1;
166  geom->planimetric = 1;
167 
168  geom->depths = TEST_N_NUM_DEPTHS;
169  geom->rows = TEST_N_NUM_ROWS;
170  geom->cols = TEST_N_NUM_COLS;
171 
172 
173  relax = create_relax_array_3d();
174  pot = create_potential_array_3d();
175 
176  field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, NULL);
177  field =
178  N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field);
179 
180  /*compute stats */
183 
185 
186  N_free_array_3d(relax);
187  N_free_array_3d(pot);
188 
189  relax = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
190  pot = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
191  xcomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
192  ycomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
193  zcomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
194 
195  N_put_array_3d_d_value(relax, 0, 0, 0, 1);
196  N_put_array_3d_d_value(relax, 0, 1, 0, 1);
197  N_put_array_3d_d_value(relax, 0, 2, 0, 1);
198  N_put_array_3d_d_value(relax, 1, 0, 0, 1);
199  N_put_array_3d_d_value(relax, 1, 1, 0, 1);
200  N_put_array_3d_d_value(relax, 1, 2, 0, 1);
201  N_put_array_3d_d_value(relax, 2, 0, 0, 1);
202  N_put_array_3d_d_value(relax, 2, 1, 0, 1);
203  N_put_array_3d_d_value(relax, 2, 2, 0, 1);
204 
205  N_put_array_3d_d_value(relax, 0, 0, 1, 1);
206  N_put_array_3d_d_value(relax, 0, 1, 1, 1);
207  N_put_array_3d_d_value(relax, 0, 2, 1, 1);
208  N_put_array_3d_d_value(relax, 1, 0, 1, 1);
209  N_put_array_3d_d_value(relax, 1, 1, 1, 1);
210  N_put_array_3d_d_value(relax, 1, 2, 1, 1);
211  N_put_array_3d_d_value(relax, 2, 0, 1, 1);
212  N_put_array_3d_d_value(relax, 2, 1, 1, 1);
213  N_put_array_3d_d_value(relax, 2, 2, 1, 1);
214 
215 
216  N_put_array_3d_d_value(relax, 0, 0, 2, 1);
217  N_put_array_3d_d_value(relax, 0, 1, 2, 1);
218  N_put_array_3d_d_value(relax, 0, 2, 2, 1);
219  N_put_array_3d_d_value(relax, 1, 0, 2, 1);
220  N_put_array_3d_d_value(relax, 1, 1, 2, 1);
221  N_put_array_3d_d_value(relax, 1, 2, 2, 1);
222  N_put_array_3d_d_value(relax, 2, 0, 2, 1);
223  N_put_array_3d_d_value(relax, 2, 1, 2, 1);
224  N_put_array_3d_d_value(relax, 2, 2, 2, 1);
225 
226 
233  N_put_array_3d_d_value(pot, 0, 0, 0, 1.0);
234  N_put_array_3d_d_value(pot, 1, 0, 0, 2.0);
235  N_put_array_3d_d_value(pot, 2, 0, 0, 6.0);
236  N_put_array_3d_d_value(pot, 0, 1, 0, 3.0);
237  N_put_array_3d_d_value(pot, 1, 1, 0, 7.0);
238  N_put_array_3d_d_value(pot, 2, 1, 0, 10.0);
239  N_put_array_3d_d_value(pot, 0, 2, 0, 8.0);
240  N_put_array_3d_d_value(pot, 1, 2, 0, 15.0);
241  N_put_array_3d_d_value(pot, 2, 2, 0, 25.0);
242 
243  N_put_array_3d_d_value(pot, 0, 0, 1, 1.2);
244  N_put_array_3d_d_value(pot, 1, 0, 1, 2.2);
245  N_put_array_3d_d_value(pot, 2, 0, 1, 6.2);
246  N_put_array_3d_d_value(pot, 0, 1, 1, 3.2);
247  N_put_array_3d_d_value(pot, 1, 1, 1, 7.2);
248  N_put_array_3d_d_value(pot, 2, 1, 1, 10.2);
249  N_put_array_3d_d_value(pot, 0, 2, 1, 8.2);
250  N_put_array_3d_d_value(pot, 1, 2, 1, 15.2);
251  N_put_array_3d_d_value(pot, 2, 2, 1, 25.2);
252 
253 
254  N_put_array_3d_d_value(pot, 0, 0, 2, 1.5);
255  N_put_array_3d_d_value(pot, 1, 0, 2, 2.5);
256  N_put_array_3d_d_value(pot, 2, 0, 2, 6.5);
257  N_put_array_3d_d_value(pot, 0, 1, 2, 3.5);
258  N_put_array_3d_d_value(pot, 1, 1, 2, 7.5);
259  N_put_array_3d_d_value(pot, 2, 1, 2, 10.5);
260  N_put_array_3d_d_value(pot, 0, 2, 2, 8.5);
261  N_put_array_3d_d_value(pot, 1, 2, 2, 15.5);
262  N_put_array_3d_d_value(pot, 2, 2, 2, 25.5);
263 
264 
265  geom->depths = 3;
266  geom->rows = 3;
267  geom->cols = 3;
268 
269  field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, NULL);
270  field =
271  N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field);
273 
274  grad = N_get_gradient_3d(field, NULL, 0, 0, 0);
275  printf
276  ("Gradient 3d: NC %g == 0 ; SC %g == 2 ; WC %g == 0 ; EC %g == -1 BC %g == 0 TC %g == -0.2\n",
277  grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC);
278 
279  grad = N_get_gradient_3d(field, grad, 1, 0, 0);
280  printf
281  ("Gradient 3d: NC %g == 0 ; SC %g == 5 ; WC %g == -1 ; EC %g == -4 BC %g == 0 TC %g == -0.2\n",
282  grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC);
283  N_free_gradient_3d(grad);
284 
285  grad = N_get_gradient_3d(field, NULL, 1, 1, 1);
286  printf
287  ("Gradient 3d: NC %g == 5 ; SC %g == 8 ; WC %g == -4 ; EC %g == -3 BC %g == -0.2 TC %g == -0.3\n",
288  grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC);
289 
290  grad = N_get_gradient_3d(field, grad, 1, 2, 2);
291  printf
292  ("Gradient 3d: NC %g == 8 ; SC %g == 0 ; WC %g == -7 ; EC %g == -10 BC %g == -0.3 TC %g == 0\n",
293  grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC);
294  N_free_gradient_3d(grad);
295 
296  grad = N_get_gradient_3d(field, NULL, 2, 2, 2);
297  printf
298  ("Gradient 3d: NC %g ==15 ; SC %g == 0 ; WC %g == -10 ; EC %g == 0 BC %g == -0.3 TC %g == 0\n",
299  grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC);
300  N_free_gradient_3d(grad);
301 
302  N_compute_gradient_field_components_3d(field, xcomp, ycomp, zcomp);
303 
304  /*
305  N_print_array_3d(xcomp);
306  N_print_array_3d(ycomp);
307  N_print_array_3d(zcomp);
308  */
309 
311  G_free(geom);
312 
313  N_free_array_3d(xcomp);
314  N_free_array_3d(ycomp);
315  N_free_array_3d(zcomp);
316  N_free_array_3d(relax);
317  N_free_array_3d(pot);
318 
319  return 0;
320 }
321 
322 
323 /* *************************************************************** */
324 /* Test the gradient calculation with 2d array data ************** */
325 /* *************************************************************** */
326 int test_gradient_2d(void)
327 {
328  N_array_2d *relax = NULL;
329  N_array_2d *pot = NULL;
330  N_array_2d *xcomp = NULL;
331  N_array_2d *ycomp = NULL;
332  N_gradient_field_2d *field = NULL;
333  N_geom_data *geom = NULL;
334  N_gradient_2d *grad = NULL;
335  N_gradient_neighbours_2d *grad_2d = NULL;
336 
337  geom = N_alloc_geom_data();
338 
339  geom->dx = 1;
340  geom->dy = 1;
341  geom->dz = 1;
342 
343  geom->Az = 1;
344  geom->planimetric = 1;
345 
346  geom->rows = TEST_N_NUM_ROWS;
347  geom->cols = TEST_N_NUM_COLS;
348 
349 
350  relax = create_relax_array_2d();
351  pot = create_potential_array_2d();
352 
353 
354  field = N_compute_gradient_field_2d(pot, relax, relax, geom, field);
355  field = N_compute_gradient_field_2d(pot, relax, relax, geom, field);
356 
357  /*compute stats */
360 
362 
363  N_free_array_2d(relax);
364  N_free_array_2d(pot);
365 
366  relax = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
367  pot = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
368  xcomp = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
369  ycomp = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
370 
371  N_put_array_2d_d_value(relax, 0, 0, 1.0);
372  N_put_array_2d_d_value(relax, 0, 1, 1.0);
373  N_put_array_2d_d_value(relax, 0, 2, 1.0);
374  N_put_array_2d_d_value(relax, 1, 0, 1.0);
375  N_put_array_2d_d_value(relax, 1, 1, 1.0);
376  N_put_array_2d_d_value(relax, 1, 2, 1.0);
377  N_put_array_2d_d_value(relax, 2, 0, 1.0);
378  N_put_array_2d_d_value(relax, 2, 1, 1.0);
379  N_put_array_2d_d_value(relax, 2, 2, 1.0);
380 
387  N_put_array_2d_d_value(pot, 0, 0, 1.0);
388  N_put_array_2d_d_value(pot, 1, 0, 2.0);
389  N_put_array_2d_d_value(pot, 2, 0, 6.0);
390  N_put_array_2d_d_value(pot, 0, 1, 3.0);
391  N_put_array_2d_d_value(pot, 1, 1, 7.0);
392  N_put_array_2d_d_value(pot, 2, 1, 10.0);
393  N_put_array_2d_d_value(pot, 0, 2, 8.0);
394  N_put_array_2d_d_value(pot, 1, 2, 15.0);
395  N_put_array_2d_d_value(pot, 2, 2, 25.0);
396 
397  geom->rows = 3;
398  geom->cols = 3;
399 
400  field = N_compute_gradient_field_2d(pot, relax, relax, geom, NULL);
401  field = N_compute_gradient_field_2d(pot, relax, relax, geom, field);
403 
404  /*common gradient calculation */
405  grad = N_get_gradient_2d(field, NULL, 0, 0);
406  G_message
407  ("Gradient 2d: pos 0,0 NC %g == 0 ; SC %g == 2 ; WC %g == 0 ; EC %g == -1\n",
408  grad->NC, grad->SC, grad->WC, grad->EC);
409 
410  grad = N_get_gradient_2d(field, grad, 1, 0);
411  G_message
412  ("Gradient 2d: pos 1,0 NC %g == 0 ; SC %g == 5 ; WC %g == -1 ; EC %g == -4\n",
413  grad->NC, grad->SC, grad->WC, grad->EC);
414  N_free_gradient_2d(grad);
415 
416  grad = N_get_gradient_2d(field, NULL, 1, 1);
417  G_message
418  ("Gradient 2d: pos 1,1 NC %g == 5 ; SC %g == 8 ; WC %g == -4 ; EC %g == -3\n",
419  grad->NC, grad->SC, grad->WC, grad->EC);
420 
421  grad = N_get_gradient_2d(field, grad, 1, 2);
422  G_message
423  ("Gradient 2d: pos 1,2 NC %g == 8 ; SC %g == 0 ; WC %g == -7 ; EC %g == -10\n",
424  grad->NC, grad->SC, grad->WC, grad->EC);
425  N_free_gradient_2d(grad);
426 
427  grad = N_get_gradient_2d(field, NULL, 2, 2);
428  G_message
429  ("Gradient 2d: pos 2,2 NC %g ==15 ; SC %g == 0 ; WC %g == -10 ; EC %g == 0\n",
430  grad->NC, grad->SC, grad->WC, grad->EC);
431  N_free_gradient_2d(grad);
432 
433  N_compute_gradient_field_components_2d(field, xcomp, ycomp);
434 
435  /*gradient neighbour calculation */
436  grad_2d = N_get_gradient_neighbours_2d(field, NULL, 1, 1);
437  grad_2d = N_get_gradient_neighbours_2d(field, grad_2d, 1, 1);
438  G_message
439  ("N_gradient_neighbours_x; pos 1,1 NWN %g NEN %g WC %g EC %g SWS %g SES %g\n",
440  grad_2d->x->NWN, grad_2d->x->NEN, grad_2d->x->WC, grad_2d->x->EC,
441  grad_2d->x->SWS, grad_2d->x->SES);
442  G_message
443  ("N_gradient_neighbours_y: pos 1,1 NWW %g NEE %g NC %g SC %g SWW %g SEE %g\n",
444  grad_2d->y->NWW, grad_2d->y->NEE, grad_2d->y->NC, grad_2d->y->SC,
445  grad_2d->y->SWW, grad_2d->y->SEE);
446 
448 
449  /*
450  N_print_array_2d(xcomp);
451  N_print_array_2d(ycomp);
452  */
453 
455  G_free(geom);
456 
457  N_free_array_2d(xcomp);
458  N_free_array_2d(ycomp);
459  N_free_array_2d(relax);
460  N_free_array_2d(pot);
461 
462 
463  return 0;
464 }