GRASS Programmer's Manual  6.4.2(2012)
test_gradient.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:      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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines