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 matrix assembling 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_matrix_assemble_2d(void); 00029 static int test_matrix_assemble_3d(void); 00030 static N_array_2d *create_status_array_2d(void); 00031 static N_array_3d *create_status_array_3d(void); 00032 static N_array_2d *create_value_array_2d(void); 00033 static N_array_3d *create_value_array_3d(void); 00034 00035 /* *************************************************************** */ 00036 /* Performe the les assmbling tests ****************************** */ 00037 /* *************************************************************** */ 00038 int unit_test_assemble(void) 00039 { 00040 int sum = 0; 00041 00042 G_message(_("\n++ Running assembling unit tests ++")); 00043 00044 G_message(_("\t 1. testing 2d assembling")); 00045 sum += test_matrix_assemble_2d(); 00046 00047 G_message(_("\t 2. testing 3d assembling")); 00048 sum += test_matrix_assemble_3d(); 00049 00050 if (sum > 0) 00051 G_warning(_("\n-- Assembling unit tests failure --")); 00052 else 00053 G_message(_("\n-- Assembling 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_status_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 00072 if (j == 1) { 00073 N_put_array_2d_c_value(data, i, j, 2); 00074 } 00075 else { 00076 N_put_array_2d_c_value(data, i, j, 1); 00077 } 00078 } 00079 } 00080 return data; 00081 } 00082 00083 /* *************************************************************** */ 00084 /* Create a value array ****************************************** */ 00085 /* *************************************************************** */ 00086 N_array_2d *create_value_array_2d(void) 00087 { 00088 N_array_2d *data; 00089 int i, j; 00090 00091 data = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, DCELL_TYPE); 00092 00093 #pragma omp parallel for private (i, j) shared (data) 00094 for (j = 0; j < TEST_N_NUM_ROWS; j++) { 00095 for (i = 0; i < TEST_N_NUM_COLS; i++) { 00096 00097 if (j == 1) { 00098 N_put_array_2d_d_value(data, i, j, 50); 00099 } 00100 else { 00101 N_put_array_2d_d_value(data, i, j, 1); 00102 } 00103 } 00104 } 00105 return data; 00106 } 00107 00108 /* *************************************************************** */ 00109 /* Create the status array with values of 1 and 2 **************** */ 00110 /* *************************************************************** */ 00111 N_array_3d *create_status_array_3d(void) 00112 { 00113 N_array_3d *data; 00114 int i, j, k; 00115 00116 data = 00117 N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 00118 1, FCELL_TYPE); 00119 00120 #pragma omp parallel for private (i, j, k) shared (data) 00121 for (k = 0; k < TEST_N_NUM_DEPTHS; k++) 00122 for (j = 0; j < TEST_N_NUM_ROWS; j++) { 00123 for (i = 0; i < TEST_N_NUM_COLS; i++) { 00124 00125 if (i == 0 && j == 1) { 00126 N_put_array_3d_f_value(data, i, j, k, 2.0); 00127 } 00128 else { 00129 00130 N_put_array_3d_f_value(data, i, j, k, 1.0); 00131 } 00132 } 00133 } 00134 00135 return data; 00136 00137 } 00138 00139 /* *************************************************************** */ 00140 /* Create a value array ****************************************** */ 00141 /* *************************************************************** */ 00142 N_array_3d *create_value_array_3d(void) 00143 { 00144 N_array_3d *data; 00145 int i, j, k; 00146 00147 data = 00148 N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 00149 1, DCELL_TYPE); 00150 00151 #pragma omp parallel for private (i, j, k) shared (data) 00152 for (k = 0; k < TEST_N_NUM_DEPTHS; k++) 00153 for (j = 0; j < TEST_N_NUM_ROWS; j++) { 00154 for (i = 0; i < TEST_N_NUM_COLS; i++) { 00155 00156 00157 if (i == 0 && j == 1) { 00158 N_put_array_3d_f_value(data, i, j, k, 50); 00159 } 00160 else { 00161 00162 N_put_array_3d_f_value(data, i, j, k, 1); 00163 } 00164 } 00165 } 00166 00167 return data; 00168 00169 } 00170 00171 /* *************************************************************** */ 00172 /* Test the matrix assembling with 3d array data ***************** */ 00173 /* *************************************************************** */ 00174 int test_matrix_assemble_3d(void) 00175 { 00176 N_geom_data *geom; 00177 N_les *les; 00178 N_les_callback_3d *call; 00179 N_array_3d *status; 00180 N_array_3d *start_val; 00181 00182 00183 /*set the callback to default */ 00184 call = N_alloc_les_callback_3d(); 00185 00186 status = create_status_array_3d(); 00187 start_val = create_value_array_3d(); 00188 00189 geom = N_alloc_geom_data(); 00190 00191 geom->dx = 1; 00192 geom->dy = 1; 00193 geom->dz = 1; 00194 00195 geom->Az = 1; 00196 00197 geom->depths = TEST_N_NUM_DEPTHS; 00198 geom->rows = TEST_N_NUM_ROWS; 00199 geom->cols = TEST_N_NUM_COLS; 00200 00201 /*Assemble the matrix */ 00202 les = 00203 N_assemble_les_3d(N_SPARSE_LES, geom, status, start_val, NULL, call); 00204 N_free_les(les); 00205 les = 00206 N_assemble_les_3d_active(N_SPARSE_LES, geom, status, start_val, NULL, 00207 call); 00208 N_free_les(les); 00209 les = 00210 N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, status, start_val, 00211 NULL, call); 00212 N_les_integrate_dirichlet_3d(les, geom, status, start_val); 00213 N_free_les(les); 00214 00215 les = 00216 N_assemble_les_3d(N_NORMAL_LES, geom, status, start_val, NULL, call); 00217 N_free_les(les); 00218 les = 00219 N_assemble_les_3d_active(N_NORMAL_LES, geom, status, start_val, NULL, 00220 call); 00221 N_free_les(les); 00222 les = 00223 N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, status, start_val, 00224 NULL, call); 00225 N_les_integrate_dirichlet_3d(les, geom, status, start_val); 00226 N_free_les(les); 00227 00228 G_free(geom); 00229 G_free(call); 00230 00231 return 0; 00232 } 00233 00234 00235 /* *************************************************************** */ 00236 /* Test the matrix assembling with 2d array data ***************** */ 00237 /* *************************************************************** */ 00238 int test_matrix_assemble_2d(void) 00239 { 00240 00241 N_geom_data *geom; 00242 N_les *les; 00243 N_les_callback_2d *call; 00244 N_array_2d *status; 00245 N_array_2d *start_val; 00246 00247 /*set the callback to default */ 00248 call = N_alloc_les_callback_2d(); 00249 00250 status = create_status_array_2d(); 00251 start_val = create_value_array_2d(); 00252 00253 geom = N_alloc_geom_data(); 00254 00255 geom->dx = 1; 00256 geom->dy = 1; 00257 00258 geom->Az = 1; 00259 00260 geom->rows = TEST_N_NUM_ROWS; 00261 geom->cols = TEST_N_NUM_COLS; 00262 00263 /*Assemble the matrix */ 00264 les = 00265 N_assemble_les_2d(N_SPARSE_LES, geom, status, start_val, NULL, call); 00266 N_free_les(les); 00267 les = 00268 N_assemble_les_2d_active(N_SPARSE_LES, geom, status, start_val, NULL, 00269 call); 00270 N_free_les(les); 00271 les = 00272 N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, status, start_val, 00273 NULL, call); 00274 N_les_integrate_dirichlet_2d(les, geom, status, start_val); 00275 N_free_les(les); 00276 00277 les = 00278 N_assemble_les_2d(N_NORMAL_LES, geom, status, start_val, NULL, call); 00279 N_free_les(les); 00280 les = 00281 N_assemble_les_2d_active(N_NORMAL_LES, geom, status, start_val, NULL, 00282 call); 00283 N_free_les(les); 00284 les = 00285 N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, status, start_val, 00286 NULL, call); 00287 N_les_integrate_dirichlet_2d(les, geom, status, start_val); 00288 N_free_les(les); 00289 00290 G_free(geom); 00291 G_free(call); 00292 00293 return 0; 00294 }