GRASS Programmer's Manual
6.4.2(2012)
|
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: Unit tests for gradient calculation 00009 * 00010 * COPYRIGHT: (C) 2000 by the GRASS Development Team 00011 * 00012 * This program is free software under the GNU General Public 00013 * License (>=v2). Read the file COPYING that comes with GRASS 00014 * for details. 00015 * 00016 *****************************************************************************/ 00017 00018 00019 00020 #include <stdio.h> 00021 #include <stdlib.h> 00022 #include <string.h> 00023 #include <grass/glocale.h> 00024 #include <grass/N_pde.h> 00025 #include "test_gpde_lib.h" 00026 00027 /* prototypes */ 00028 static int test_gradient_2d(void); 00029 static int test_gradient_3d(void); 00030 static N_array_2d *create_relax_array_2d(void); 00031 static N_array_3d *create_relax_array_3d(void); 00032 static N_array_2d *create_potential_array_2d(void); 00033 static N_array_3d *create_potential_array_3d(void); 00034 00035 /* *************************************************************** */ 00036 /* Performe the gradient tests *********************************** */ 00037 /* *************************************************************** */ 00038 int unit_test_gradient(void) 00039 { 00040 int sum = 0; 00041 00042 G_message(_("\n++ Running gradient unit tests ++")); 00043 00044 G_message(_("\t 1. testing 2d gradient")); 00045 sum += test_gradient_2d(); 00046 00047 G_message(_("\t 2. testing 3d gradient")); 00048 sum += test_gradient_3d(); 00049 00050 if (sum > 0) 00051 G_warning(_("\n-- Gradient unit tests failure --")); 00052 else 00053 G_message(_("\n-- Gradient unit tests finished successfully --")); 00054 00055 return sum; 00056 } 00057 00058 /* *************************************************************** */ 00059 /* Create the status array with values of 1 and 2 **************** */ 00060 /* *************************************************************** */ 00061 N_array_2d *create_relax_array_2d(void) 00062 { 00063 N_array_2d *data; 00064 int i, j; 00065 00066 data = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, CELL_TYPE); 00067 00068 #pragma omp parallel for private (i, j) shared (data) 00069 for (j = 0; j < TEST_N_NUM_ROWS; j++) { 00070 for (i = 0; i < TEST_N_NUM_COLS; i++) { 00071 N_put_array_2d_c_value(data, i, j, 1); 00072 } 00073 } 00074 return data; 00075 } 00076 00077 /* *************************************************************** */ 00078 /* Create a value array ****************************************** */ 00079 /* *************************************************************** */ 00080 N_array_2d *create_potential_array_2d(void) 00081 { 00082 N_array_2d *data; 00083 int i, j; 00084 00085 data = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, DCELL_TYPE); 00086 00087 #pragma omp parallel for private (i, j) shared (data) 00088 for (j = 0; j < TEST_N_NUM_ROWS; j++) { 00089 for (i = 0; i < TEST_N_NUM_COLS; i++) { 00090 N_put_array_2d_d_value(data, i, j, (double)i * j); 00091 } 00092 } 00093 return data; 00094 } 00095 00096 /* *************************************************************** */ 00097 /* Create the status array with values of 1 and 2 **************** */ 00098 /* *************************************************************** */ 00099 N_array_3d *create_relax_array_3d(void) 00100 { 00101 N_array_3d *data; 00102 int i, j, k; 00103 00104 data = 00105 N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 00106 1, FCELL_TYPE); 00107 00108 #pragma omp parallel for private (i, j, k) shared (data) 00109 for (k = 0; k < TEST_N_NUM_DEPTHS; k++) { 00110 for (j = 0; j < TEST_N_NUM_ROWS; j++) { 00111 for (i = 0; i < TEST_N_NUM_COLS; i++) { 00112 N_put_array_3d_f_value(data, i, j, k, 1.0); 00113 } 00114 } 00115 } 00116 00117 return data; 00118 00119 } 00120 00121 /* *************************************************************** */ 00122 /* Create a value array ****************************************** */ 00123 /* *************************************************************** */ 00124 N_array_3d *create_potential_array_3d(void) 00125 { 00126 N_array_3d *data; 00127 int i, j, k; 00128 00129 data = 00130 N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 00131 1, DCELL_TYPE); 00132 00133 #pragma omp parallel for private (i, j, k) shared (data) 00134 for (k = 0; k < TEST_N_NUM_DEPTHS; k++) 00135 for (j = 0; j < TEST_N_NUM_ROWS; j++) { 00136 for (i = 0; i < TEST_N_NUM_COLS; i++) { 00137 N_put_array_3d_f_value(data, i, j, k, (double)i * j * k); 00138 } 00139 } 00140 00141 return data; 00142 00143 } 00144 00145 /* *************************************************************** */ 00146 /* Test the gradient calculation with 3d array data ************** */ 00147 /* *************************************************************** */ 00148 int test_gradient_3d(void) 00149 { 00150 N_array_3d *relax = NULL; 00151 N_array_3d *pot = NULL; 00152 N_array_3d *xcomp = NULL; 00153 N_array_3d *ycomp = NULL; 00154 N_array_3d *zcomp = NULL; 00155 N_gradient_field_3d *field = NULL; 00156 N_gradient_3d *grad = NULL; 00157 N_geom_data *geom = NULL; 00158 00159 geom = N_alloc_geom_data(); 00160 00161 geom->dx = 1; 00162 geom->dy = 1; 00163 geom->dz = 1; 00164 00165 geom->Az = 1; 00166 geom->planimetric = 1; 00167 00168 geom->depths = TEST_N_NUM_DEPTHS; 00169 geom->rows = TEST_N_NUM_ROWS; 00170 geom->cols = TEST_N_NUM_COLS; 00171 00172 00173 relax = create_relax_array_3d(); 00174 pot = create_potential_array_3d(); 00175 00176 field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, NULL); 00177 field = 00178 N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field); 00179 00180 /*compute stats */ 00181 N_calc_gradient_field_3d_stats(field); 00182 N_print_gradient_field_3d_info(field); 00183 00184 N_free_gradient_field_3d(field); 00185 00186 N_free_array_3d(relax); 00187 N_free_array_3d(pot); 00188 00189 relax = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE); 00190 pot = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE); 00191 xcomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE); 00192 ycomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE); 00193 zcomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE); 00194 00195 N_put_array_3d_d_value(relax, 0, 0, 0, 1); 00196 N_put_array_3d_d_value(relax, 0, 1, 0, 1); 00197 N_put_array_3d_d_value(relax, 0, 2, 0, 1); 00198 N_put_array_3d_d_value(relax, 1, 0, 0, 1); 00199 N_put_array_3d_d_value(relax, 1, 1, 0, 1); 00200 N_put_array_3d_d_value(relax, 1, 2, 0, 1); 00201 N_put_array_3d_d_value(relax, 2, 0, 0, 1); 00202 N_put_array_3d_d_value(relax, 2, 1, 0, 1); 00203 N_put_array_3d_d_value(relax, 2, 2, 0, 1); 00204 00205 N_put_array_3d_d_value(relax, 0, 0, 1, 1); 00206 N_put_array_3d_d_value(relax, 0, 1, 1, 1); 00207 N_put_array_3d_d_value(relax, 0, 2, 1, 1); 00208 N_put_array_3d_d_value(relax, 1, 0, 1, 1); 00209 N_put_array_3d_d_value(relax, 1, 1, 1, 1); 00210 N_put_array_3d_d_value(relax, 1, 2, 1, 1); 00211 N_put_array_3d_d_value(relax, 2, 0, 1, 1); 00212 N_put_array_3d_d_value(relax, 2, 1, 1, 1); 00213 N_put_array_3d_d_value(relax, 2, 2, 1, 1); 00214 00215 00216 N_put_array_3d_d_value(relax, 0, 0, 2, 1); 00217 N_put_array_3d_d_value(relax, 0, 1, 2, 1); 00218 N_put_array_3d_d_value(relax, 0, 2, 2, 1); 00219 N_put_array_3d_d_value(relax, 1, 0, 2, 1); 00220 N_put_array_3d_d_value(relax, 1, 1, 2, 1); 00221 N_put_array_3d_d_value(relax, 1, 2, 2, 1); 00222 N_put_array_3d_d_value(relax, 2, 0, 2, 1); 00223 N_put_array_3d_d_value(relax, 2, 1, 2, 1); 00224 N_put_array_3d_d_value(relax, 2, 2, 2, 1); 00225 00226 00233 N_put_array_3d_d_value(pot, 0, 0, 0, 1.0); 00234 N_put_array_3d_d_value(pot, 1, 0, 0, 2.0); 00235 N_put_array_3d_d_value(pot, 2, 0, 0, 6.0); 00236 N_put_array_3d_d_value(pot, 0, 1, 0, 3.0); 00237 N_put_array_3d_d_value(pot, 1, 1, 0, 7.0); 00238 N_put_array_3d_d_value(pot, 2, 1, 0, 10.0); 00239 N_put_array_3d_d_value(pot, 0, 2, 0, 8.0); 00240 N_put_array_3d_d_value(pot, 1, 2, 0, 15.0); 00241 N_put_array_3d_d_value(pot, 2, 2, 0, 25.0); 00242 00243 N_put_array_3d_d_value(pot, 0, 0, 1, 1.2); 00244 N_put_array_3d_d_value(pot, 1, 0, 1, 2.2); 00245 N_put_array_3d_d_value(pot, 2, 0, 1, 6.2); 00246 N_put_array_3d_d_value(pot, 0, 1, 1, 3.2); 00247 N_put_array_3d_d_value(pot, 1, 1, 1, 7.2); 00248 N_put_array_3d_d_value(pot, 2, 1, 1, 10.2); 00249 N_put_array_3d_d_value(pot, 0, 2, 1, 8.2); 00250 N_put_array_3d_d_value(pot, 1, 2, 1, 15.2); 00251 N_put_array_3d_d_value(pot, 2, 2, 1, 25.2); 00252 00253 00254 N_put_array_3d_d_value(pot, 0, 0, 2, 1.5); 00255 N_put_array_3d_d_value(pot, 1, 0, 2, 2.5); 00256 N_put_array_3d_d_value(pot, 2, 0, 2, 6.5); 00257 N_put_array_3d_d_value(pot, 0, 1, 2, 3.5); 00258 N_put_array_3d_d_value(pot, 1, 1, 2, 7.5); 00259 N_put_array_3d_d_value(pot, 2, 1, 2, 10.5); 00260 N_put_array_3d_d_value(pot, 0, 2, 2, 8.5); 00261 N_put_array_3d_d_value(pot, 1, 2, 2, 15.5); 00262 N_put_array_3d_d_value(pot, 2, 2, 2, 25.5); 00263 00264 00265 geom->depths = 3; 00266 geom->rows = 3; 00267 geom->cols = 3; 00268 00269 field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, NULL); 00270 field = 00271 N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field); 00272 N_print_gradient_field_3d_info(field); 00273 00274 grad = N_get_gradient_3d(field, NULL, 0, 0, 0); 00275 printf 00276 ("Gradient 3d: NC %g == 0 ; SC %g == 2 ; WC %g == 0 ; EC %g == -1 BC %g == 0 TC %g == -0.2\n", 00277 grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC); 00278 00279 grad = N_get_gradient_3d(field, grad, 1, 0, 0); 00280 printf 00281 ("Gradient 3d: NC %g == 0 ; SC %g == 5 ; WC %g == -1 ; EC %g == -4 BC %g == 0 TC %g == -0.2\n", 00282 grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC); 00283 N_free_gradient_3d(grad); 00284 00285 grad = N_get_gradient_3d(field, NULL, 1, 1, 1); 00286 printf 00287 ("Gradient 3d: NC %g == 5 ; SC %g == 8 ; WC %g == -4 ; EC %g == -3 BC %g == -0.2 TC %g == -0.3\n", 00288 grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC); 00289 00290 grad = N_get_gradient_3d(field, grad, 1, 2, 2); 00291 printf 00292 ("Gradient 3d: NC %g == 8 ; SC %g == 0 ; WC %g == -7 ; EC %g == -10 BC %g == -0.3 TC %g == 0\n", 00293 grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC); 00294 N_free_gradient_3d(grad); 00295 00296 grad = N_get_gradient_3d(field, NULL, 2, 2, 2); 00297 printf 00298 ("Gradient 3d: NC %g ==15 ; SC %g == 0 ; WC %g == -10 ; EC %g == 0 BC %g == -0.3 TC %g == 0\n", 00299 grad->NC, grad->SC, grad->WC, grad->EC, grad->BC, grad->TC); 00300 N_free_gradient_3d(grad); 00301 00302 N_compute_gradient_field_components_3d(field, xcomp, ycomp, zcomp); 00303 00304 /* 00305 N_print_array_3d(xcomp); 00306 N_print_array_3d(ycomp); 00307 N_print_array_3d(zcomp); 00308 */ 00309 00310 N_free_gradient_field_3d(field); 00311 G_free(geom); 00312 00313 N_free_array_3d(xcomp); 00314 N_free_array_3d(ycomp); 00315 N_free_array_3d(zcomp); 00316 N_free_array_3d(relax); 00317 N_free_array_3d(pot); 00318 00319 return 0; 00320 } 00321 00322 00323 /* *************************************************************** */ 00324 /* Test the gradient calculation with 2d array data ************** */ 00325 /* *************************************************************** */ 00326 int test_gradient_2d(void) 00327 { 00328 N_array_2d *relax = NULL; 00329 N_array_2d *pot = NULL; 00330 N_array_2d *xcomp = NULL; 00331 N_array_2d *ycomp = NULL; 00332 N_gradient_field_2d *field = NULL; 00333 N_geom_data *geom = NULL; 00334 N_gradient_2d *grad = NULL; 00335 N_gradient_neighbours_2d *grad_2d = NULL; 00336 00337 geom = N_alloc_geom_data(); 00338 00339 geom->dx = 1; 00340 geom->dy = 1; 00341 geom->dz = 1; 00342 00343 geom->Az = 1; 00344 geom->planimetric = 1; 00345 00346 geom->rows = TEST_N_NUM_ROWS; 00347 geom->cols = TEST_N_NUM_COLS; 00348 00349 00350 relax = create_relax_array_2d(); 00351 pot = create_potential_array_2d(); 00352 00353 00354 field = N_compute_gradient_field_2d(pot, relax, relax, geom, field); 00355 field = N_compute_gradient_field_2d(pot, relax, relax, geom, field); 00356 00357 /*compute stats */ 00358 N_calc_gradient_field_2d_stats(field); 00359 N_print_gradient_field_2d_info(field); 00360 00361 N_free_gradient_field_2d(field); 00362 00363 N_free_array_2d(relax); 00364 N_free_array_2d(pot); 00365 00366 relax = N_alloc_array_2d(3, 3, 0, DCELL_TYPE); 00367 pot = N_alloc_array_2d(3, 3, 0, DCELL_TYPE); 00368 xcomp = N_alloc_array_2d(3, 3, 0, DCELL_TYPE); 00369 ycomp = N_alloc_array_2d(3, 3, 0, DCELL_TYPE); 00370 00371 N_put_array_2d_d_value(relax, 0, 0, 1.0); 00372 N_put_array_2d_d_value(relax, 0, 1, 1.0); 00373 N_put_array_2d_d_value(relax, 0, 2, 1.0); 00374 N_put_array_2d_d_value(relax, 1, 0, 1.0); 00375 N_put_array_2d_d_value(relax, 1, 1, 1.0); 00376 N_put_array_2d_d_value(relax, 1, 2, 1.0); 00377 N_put_array_2d_d_value(relax, 2, 0, 1.0); 00378 N_put_array_2d_d_value(relax, 2, 1, 1.0); 00379 N_put_array_2d_d_value(relax, 2, 2, 1.0); 00380 00387 N_put_array_2d_d_value(pot, 0, 0, 1.0); 00388 N_put_array_2d_d_value(pot, 1, 0, 2.0); 00389 N_put_array_2d_d_value(pot, 2, 0, 6.0); 00390 N_put_array_2d_d_value(pot, 0, 1, 3.0); 00391 N_put_array_2d_d_value(pot, 1, 1, 7.0); 00392 N_put_array_2d_d_value(pot, 2, 1, 10.0); 00393 N_put_array_2d_d_value(pot, 0, 2, 8.0); 00394 N_put_array_2d_d_value(pot, 1, 2, 15.0); 00395 N_put_array_2d_d_value(pot, 2, 2, 25.0); 00396 00397 geom->rows = 3; 00398 geom->cols = 3; 00399 00400 field = N_compute_gradient_field_2d(pot, relax, relax, geom, NULL); 00401 field = N_compute_gradient_field_2d(pot, relax, relax, geom, field); 00402 N_print_gradient_field_2d_info(field); 00403 00404 /*common gradient calculation */ 00405 grad = N_get_gradient_2d(field, NULL, 0, 0); 00406 G_message 00407 ("Gradient 2d: pos 0,0 NC %g == 0 ; SC %g == 2 ; WC %g == 0 ; EC %g == -1\n", 00408 grad->NC, grad->SC, grad->WC, grad->EC); 00409 00410 grad = N_get_gradient_2d(field, grad, 1, 0); 00411 G_message 00412 ("Gradient 2d: pos 1,0 NC %g == 0 ; SC %g == 5 ; WC %g == -1 ; EC %g == -4\n", 00413 grad->NC, grad->SC, grad->WC, grad->EC); 00414 N_free_gradient_2d(grad); 00415 00416 grad = N_get_gradient_2d(field, NULL, 1, 1); 00417 G_message 00418 ("Gradient 2d: pos 1,1 NC %g == 5 ; SC %g == 8 ; WC %g == -4 ; EC %g == -3\n", 00419 grad->NC, grad->SC, grad->WC, grad->EC); 00420 00421 grad = N_get_gradient_2d(field, grad, 1, 2); 00422 G_message 00423 ("Gradient 2d: pos 1,2 NC %g == 8 ; SC %g == 0 ; WC %g == -7 ; EC %g == -10\n", 00424 grad->NC, grad->SC, grad->WC, grad->EC); 00425 N_free_gradient_2d(grad); 00426 00427 grad = N_get_gradient_2d(field, NULL, 2, 2); 00428 G_message 00429 ("Gradient 2d: pos 2,2 NC %g ==15 ; SC %g == 0 ; WC %g == -10 ; EC %g == 0\n", 00430 grad->NC, grad->SC, grad->WC, grad->EC); 00431 N_free_gradient_2d(grad); 00432 00433 N_compute_gradient_field_components_2d(field, xcomp, ycomp); 00434 00435 /*gradient neighbour calculation */ 00436 grad_2d = N_get_gradient_neighbours_2d(field, NULL, 1, 1); 00437 grad_2d = N_get_gradient_neighbours_2d(field, grad_2d, 1, 1); 00438 G_message 00439 ("N_gradient_neighbours_x; pos 1,1 NWN %g NEN %g WC %g EC %g SWS %g SES %g\n", 00440 grad_2d->x->NWN, grad_2d->x->NEN, grad_2d->x->WC, grad_2d->x->EC, 00441 grad_2d->x->SWS, grad_2d->x->SES); 00442 G_message 00443 ("N_gradient_neighbours_y: pos 1,1 NWW %g NEE %g NC %g SC %g SWW %g SEE %g\n", 00444 grad_2d->y->NWW, grad_2d->y->NEE, grad_2d->y->NC, grad_2d->y->SC, 00445 grad_2d->y->SWW, grad_2d->y->SEE); 00446 00447 N_free_gradient_neighbours_2d(grad_2d); 00448 00449 /* 00450 N_print_array_2d(xcomp); 00451 N_print_array_2d(ycomp); 00452 */ 00453 00454 N_free_gradient_field_2d(field); 00455 G_free(geom); 00456 00457 N_free_array_2d(xcomp); 00458 N_free_array_2d(ycomp); 00459 N_free_array_2d(relax); 00460 N_free_array_2d(pot); 00461 00462 00463 return 0; 00464 }