|
GRASS Programmer's Manual
6.5.svn(2012)-r51648
|
|
GRASS Programmer's Manual
6.5.svn(2012)-r51648
|
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: functions to assemble a linear equation system 00009 * part of the gpde library 00010 * 00011 * COPYRIGHT: (C) 2000 by the GRASS Development Team 00012 * 00013 * This program is free software under the GNU General Public 00014 * License (>=v2). Read the file COPYING that comes with GRASS 00015 * for details. 00016 * 00017 *****************************************************************************/ 00018 00019 00020 #include <math.h> 00021 #include "grass/N_pde.h" 00022 00023 /* local protos */ 00024 static int make_les_entry_2d(int i, int j, int offset_i, int offset_j, 00025 int count, int pos, N_les * les, 00026 G_math_spvector * spvect, 00027 N_array_2d * cell_count, N_array_2d * status, 00028 N_array_2d * start_val, double entry, 00029 int cell_type); 00030 00031 static int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j, 00032 int offset_k, int count, int pos, N_les * les, 00033 G_math_spvector * spvect, 00034 N_array_3d * cell_count, N_array_3d * status, 00035 N_array_3d * start_val, double entry, 00036 int cell_type); 00037 00038 /* *************************************************************** * 00039 * ********************** N_alloc_5star ************************** * 00040 * *************************************************************** */ 00046 N_data_star *N_alloc_5star(void) 00047 { 00048 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star)); 00049 00050 star->type = N_5_POINT_STAR; 00051 star->count = 5; 00052 return star; 00053 } 00054 00055 /* *************************************************************** * 00056 * ********************* N_alloc_7star *************************** * 00057 * *************************************************************** */ 00063 N_data_star *N_alloc_7star(void) 00064 { 00065 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star)); 00066 00067 star->type = N_7_POINT_STAR; 00068 star->count = 7; 00069 return star; 00070 } 00071 00072 /* *************************************************************** * 00073 * ********************* N_alloc_9star *************************** * 00074 * *************************************************************** */ 00083 N_data_star *N_alloc_9star(void) 00084 { 00085 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star)); 00086 00087 star->type = N_9_POINT_STAR; 00088 star->count = 9; 00089 return star; 00090 } 00091 00092 /* *************************************************************** * 00093 * ********************* N_alloc_27star ************************** * 00094 * *************************************************************** */ 00103 N_data_star *N_alloc_27star(void) 00104 { 00105 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star)); 00106 00107 star->type = N_27_POINT_STAR; 00108 star->count = 27; 00109 return star; 00110 } 00111 00112 /* *************************************************************** * 00113 * ********************** N_create_5star ************************* * 00114 * *************************************************************** */ 00126 N_data_star *N_create_5star(double C, double W, double E, double N, 00127 double S, double V) 00128 { 00129 N_data_star *star = N_alloc_5star(); 00130 00131 star->C = C; 00132 star->W = W; 00133 star->E = E; 00134 star->N = N; 00135 star->S = S; 00136 00137 star->V = V; 00138 00139 G_debug(5, "N_create_5star: w %g e %g n %g s %g c %g v %g\n", star->W, 00140 star->E, star->N, star->S, star->C, star->V); 00141 00142 return star; 00143 } 00144 00145 /* *************************************************************** * 00146 * ************************* N_create_7star ********************** * 00147 * *************************************************************** */ 00161 N_data_star *N_create_7star(double C, double W, double E, double N, 00162 double S, double T, double B, double V) 00163 { 00164 N_data_star *star = N_alloc_7star(); 00165 00166 star->C = C; 00167 star->W = W; 00168 star->E = E; 00169 star->N = N; 00170 star->S = S; 00171 00172 star->T = T; 00173 star->B = B; 00174 00175 star->V = V; 00176 00177 G_debug(5, "N_create_7star: w %g e %g n %g s %g t %g b %g c %g v %g\n", 00178 star->W, star->E, star->N, star->S, star->T, star->B, star->C, 00179 star->V); 00180 00181 return star; 00182 } 00183 00184 /* *************************************************************** * 00185 * ************************ N_create_9star *********************** * 00186 * *************************************************************** */ 00202 N_data_star *N_create_9star(double C, double W, double E, double N, 00203 double S, double NW, double SW, double NE, 00204 double SE, double V) 00205 { 00206 N_data_star *star = N_alloc_9star(); 00207 00208 star->C = C; 00209 star->W = W; 00210 star->E = E; 00211 star->N = N; 00212 star->S = S; 00213 00214 star->NW = NW; 00215 star->SW = SW; 00216 star->NE = NE; 00217 star->SE = SE; 00218 00219 star->V = V; 00220 00221 G_debug(5, 00222 "N_create_9star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n", 00223 star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE, 00224 star->SE, star->C, star->V); 00225 00226 return star; 00227 } 00228 00229 /* *************************************************************** * 00230 * ************************ N_create_27star *********************** * 00231 * *************************************************************** */ 00265 N_data_star *N_create_27star(double C, double W, double E, double N, double S, 00266 double NW, double SW, double NE, double SE, 00267 double T, double W_T, double E_T, double N_T, 00268 double S_T, double NW_T, double SW_T, 00269 double NE_T, double SE_T, double B, double W_B, 00270 double E_B, double N_B, double S_B, double NW_B, 00271 double SW_B, double NE_B, double SE_B, double V) 00272 { 00273 N_data_star *star = N_alloc_27star(); 00274 00275 star->C = C; 00276 star->W = W; 00277 star->E = E; 00278 star->N = N; 00279 star->S = S; 00280 00281 star->NW = NW; 00282 star->SW = SW; 00283 star->NE = NE; 00284 star->SE = SE; 00285 00286 star->T = T; 00287 star->W_T = W_T; 00288 star->E_T = E_T; 00289 star->N_T = N_T; 00290 star->S_T = S_T; 00291 00292 star->NW_T = NW_T; 00293 star->SW_T = SW_T; 00294 star->NE_T = NE_T; 00295 star->SE_T = SE_T; 00296 00297 star->B = B; 00298 star->W_B = W_B; 00299 star->E_B = E_B; 00300 star->N_B = N_B; 00301 star->S_B = S_B; 00302 00303 star->NW_B = NW_B; 00304 star->SW_B = SW_B; 00305 star->NE_B = NE_B; 00306 star->SE_B = SE_B; 00307 00308 star->V = V; 00309 00310 G_debug(5, 00311 "N_create_27star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n", 00312 star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE, 00313 star->SE, star->C, star->V); 00314 00315 G_debug(5, 00316 "N_create_27star: w_t %g e_t %g n_t %g s_t %g nw_t %g sw_t %g ne_t %g se_t %g t %g \n", 00317 star->W_T, star->E_T, star->N_T, star->S_T, star->NW_T, 00318 star->SW_T, star->NE_T, star->SE_T, star->T); 00319 00320 G_debug(5, 00321 "N_create_27star: w_b %g e_b %g n_b %g s_b %g nw_b %g sw_b %g ne_b %g se_B %g b %g\n", 00322 star->W_B, star->E_B, star->N_B, star->S_B, star->NW_B, 00323 star->SW_B, star->NE_B, star->SE_B, star->B); 00324 00325 00326 00327 return star; 00328 } 00329 00330 00331 /* *************************************************************** * 00332 * ****************** N_set_les_callback_3d_func ***************** * 00333 * *************************************************************** */ 00341 void 00342 N_set_les_callback_3d_func(N_les_callback_3d * data, 00343 N_data_star * (*callback_func_3d) ()) 00344 { 00345 data->callback = callback_func_3d; 00346 } 00347 00348 /* *************************************************************** * 00349 * *************** N_set_les_callback_2d_func ******************** * 00350 * *************************************************************** */ 00358 void 00359 N_set_les_callback_2d_func(N_les_callback_2d * data, 00360 N_data_star * (*callback_func_2d) ()) 00361 { 00362 data->callback = callback_func_2d; 00363 } 00364 00365 /* *************************************************************** * 00366 * ************** N_alloc_les_callback_3d ************************ * 00367 * *************************************************************** */ 00376 N_les_callback_3d *N_alloc_les_callback_3d(void) 00377 { 00378 N_les_callback_3d *call; 00379 00380 call = (N_les_callback_3d *) G_calloc(1, sizeof(N_les_callback_3d *)); 00381 call->callback = N_callback_template_3d; 00382 00383 return call; 00384 } 00385 00386 /* *************************************************************** * 00387 * *************** N_alloc_les_callback_2d *********************** * 00388 * *************************************************************** */ 00397 N_les_callback_2d *N_alloc_les_callback_2d(void) 00398 { 00399 N_les_callback_2d *call; 00400 00401 call = (N_les_callback_2d *) G_calloc(1, sizeof(N_les_callback_2d *)); 00402 call->callback = N_callback_template_2d; 00403 00404 return call; 00405 } 00406 00407 /* *************************************************************** * 00408 * ******************** N_callback_template_3d ******************* * 00409 * *************************************************************** */ 00424 N_data_star *N_callback_template_3d(void *data, N_geom_data * geom, int col, 00425 int row, int depth) 00426 { 00427 N_data_star *star = N_alloc_7star(); 00428 00429 star->E = 1 / geom->dx; 00430 star->W = 1 / geom->dx; 00431 star->N = 1 / geom->dy; 00432 star->S = 1 / geom->dy; 00433 star->T = 1 / geom->dz; 00434 star->B = 1 / geom->dz; 00435 star->C = -1 * (2 / geom->dx + 2 / geom->dy + 2 / geom->dz); 00436 star->V = -1; 00437 00438 G_debug(5, 00439 "N_callback_template_3d: w %g e %g n %g s %g t %g b %g c %g v %g\n", 00440 star->W, star->E, star->N, star->S, star->T, star->B, star->C, 00441 star->V); 00442 00443 00444 return star; 00445 } 00446 00447 /* *************************************************************** * 00448 * ********************* N_callback_template_2d ****************** * 00449 * *************************************************************** */ 00463 N_data_star *N_callback_template_2d(void *data, N_geom_data * geom, int col, 00464 int row) 00465 { 00466 N_data_star *star = N_alloc_9star(); 00467 00468 star->E = 1 / geom->dx; 00469 star->NE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy); 00470 star->SE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy); 00471 star->W = 1 / geom->dx; 00472 star->NW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy); 00473 star->SW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy); 00474 star->N = 1 / geom->dy; 00475 star->S = 1 / geom->dy; 00476 star->C = 00477 -1 * (star->E + star->NE + star->SE + star->W + star->NW + star->SW + 00478 star->N + star->S); 00479 star->V = 0; 00480 00481 return star; 00482 } 00483 00484 /* *************************************************************** * 00485 * ******************** N_assemble_les_2d ************************ * 00486 * *************************************************************** */ 00493 N_les *N_assemble_les_2d(int les_type, N_geom_data * geom, 00494 N_array_2d * status, N_array_2d * start_val, 00495 void *data, N_les_callback_2d * call) 00496 { 00497 return N_assemble_les_2d_param(les_type, geom, status, start_val, data, 00498 call, N_CELL_ACTIVE); 00499 } 00500 00507 N_les *N_assemble_les_2d_active(int les_type, N_geom_data * geom, 00508 N_array_2d * status, N_array_2d * start_val, 00509 void *data, N_les_callback_2d * call) 00510 { 00511 return N_assemble_les_2d_param(les_type, geom, status, start_val, data, 00512 call, N_CELL_ACTIVE); 00513 } 00514 00521 N_les *N_assemble_les_2d_dirichlet(int les_type, N_geom_data * geom, 00522 N_array_2d * status, 00523 N_array_2d * start_val, void *data, 00524 N_les_callback_2d * call) 00525 { 00526 return N_assemble_les_2d_param(les_type, geom, status, start_val, data, 00527 call, N_CELL_DIRICHLET); 00528 } 00529 00560 N_les *N_assemble_les_2d_param(int les_type, N_geom_data * geom, 00561 N_array_2d * status, N_array_2d * start_val, 00562 void *data, N_les_callback_2d * call, 00563 int cell_type) 00564 { 00565 int i, j, count = 0, pos = 0; 00566 00567 int cell_type_count = 0; 00568 00569 int **index_ij; 00570 00571 N_array_2d *cell_count; 00572 00573 N_les *les = NULL; 00574 00575 G_debug(2, 00576 "N_assemble_les_2d: starting to assemble the linear equation system"); 00577 00578 /* At first count the number of valid cells and save the 00579 * each number in a new 2d array. Those numbers are used 00580 * to create the linear equation system. 00581 * */ 00582 00583 cell_count = N_alloc_array_2d(geom->cols, geom->rows, 1, CELL_TYPE); 00584 00585 /* include dirichlet cells in the les */ 00586 if (cell_type == N_CELL_DIRICHLET) { 00587 for (j = 0; j < geom->rows; j++) { 00588 for (i = 0; i < geom->cols; i++) { 00589 /*use all non-inactive cells for les creation */ 00590 if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) && 00591 N_get_array_2d_c_value(status, i, j) < N_MAX_CELL_STATE) 00592 cell_type_count++; 00593 } 00594 } 00595 } 00596 /*use only active cell in the les */ 00597 if (cell_type == N_CELL_ACTIVE) { 00598 for (j = 0; j < geom->rows; j++) { 00599 for (i = 0; i < geom->cols; i++) { 00600 /*count only active cells */ 00601 if (N_CELL_ACTIVE == N_get_array_2d_d_value(status, i, j)) 00602 cell_type_count++; 00603 } 00604 } 00605 } 00606 00607 G_debug(2, "N_assemble_les_2d: number of used cells %i\n", 00608 cell_type_count); 00609 00610 if (cell_type_count == 0) 00611 G_fatal_error 00612 ("Not enough cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.", 00613 cell_type_count); 00614 00615 /* Then allocate the memory for the linear equation system (les). 00616 * Only valid cells are used to create the les. */ 00617 index_ij = (int **)G_calloc(cell_type_count, sizeof(int *)); 00618 for (i = 0; i < cell_type_count; i++) 00619 index_ij[i] = (int *)G_calloc(2, sizeof(int)); 00620 00621 les = N_alloc_les_Ax_b(cell_type_count, les_type); 00622 00623 count = 0; 00624 00625 /*count the number of cells which should be used to create the linear equation system */ 00626 /*save the i and j indices and create a ordered numbering */ 00627 for (j = 0; j < geom->rows; j++) { 00628 for (i = 0; i < geom->cols; i++) { 00629 /*count every non-inactive cell */ 00630 if (cell_type == N_CELL_DIRICHLET) { 00631 if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) && 00632 N_get_array_2d_c_value(status, i, j) < N_MAX_CELL_STATE) { 00633 N_put_array_2d_c_value(cell_count, i, j, count); 00634 index_ij[count][0] = i; 00635 index_ij[count][1] = j; 00636 count++; 00637 G_debug(5, 00638 "N_assemble_les_2d: non-inactive cells count %i at pos x[%i] y[%i]\n", 00639 count, i, j); 00640 } 00641 /*count every active cell */ 00642 } 00643 else if (N_CELL_ACTIVE == N_get_array_2d_c_value(status, i, j)) { 00644 N_put_array_2d_c_value(cell_count, i, j, count); 00645 index_ij[count][0] = i; 00646 index_ij[count][1] = j; 00647 count++; 00648 G_debug(5, 00649 "N_assemble_les_2d: active cells count %i at pos x[%i] y[%i]\n", 00650 count, i, j); 00651 } 00652 } 00653 } 00654 00655 G_debug(2, "N_assemble_les_2d: starting the parallel assemble loop"); 00656 00657 /* Assemble the matrix in parallel */ 00658 #pragma omp parallel for private(i, j, pos, count) schedule(static) 00659 for (count = 0; count < cell_type_count; count++) { 00660 i = index_ij[count][0]; 00661 j = index_ij[count][1]; 00662 00663 /*create the entries for the */ 00664 N_data_star *items = call->callback(data, geom, i, j); 00665 00666 /* we need a sparse vector pointer anytime */ 00667 G_math_spvector *spvect = NULL; 00668 00669 /*allocate a sprase vector */ 00670 if (les_type == N_SPARSE_LES) { 00671 spvect = G_math_alloc_spvector(items->count); 00672 } 00673 /* initial conditions */ 00674 les->x[count] = N_get_array_2d_d_value(start_val, i, j); 00675 00676 /* the entry in the vector b */ 00677 les->b[count] = items->V; 00678 00679 /* pos describes the position in the sparse vector. 00680 * the first entry is always the diagonal entry of the matrix*/ 00681 pos = 0; 00682 00683 if (les_type == N_SPARSE_LES) { 00684 spvect->index[pos] = count; 00685 spvect->values[pos] = items->C; 00686 } 00687 else { 00688 les->A[count][count] = items->C; 00689 } 00690 /* western neighbour, entry is col - 1 */ 00691 if (i > 0) { 00692 pos = make_les_entry_2d(i, j, -1, 0, count, pos, les, spvect, 00693 cell_count, status, start_val, items->W, 00694 cell_type); 00695 } 00696 /* eastern neighbour, entry col + 1 */ 00697 if (i < geom->cols - 1) { 00698 pos = make_les_entry_2d(i, j, 1, 0, count, pos, les, spvect, 00699 cell_count, status, start_val, items->E, 00700 cell_type); 00701 } 00702 /* northern neighbour, entry row - 1 */ 00703 if (j > 0) { 00704 pos = 00705 make_les_entry_2d(i, j, 0, -1, count, pos, les, spvect, 00706 cell_count, status, start_val, items->N, 00707 cell_type); 00708 } 00709 /* southern neighbour, entry row + 1 */ 00710 if (j < geom->rows - 1) { 00711 pos = make_les_entry_2d(i, j, 0, 1, count, pos, les, spvect, 00712 cell_count, status, start_val, items->S, 00713 cell_type); 00714 } 00715 /*in case of a nine point star, we have additional entries */ 00716 if (items->type == N_9_POINT_STAR) { 00717 /* north-western neighbour, entry is col - 1 row - 1 */ 00718 if (i > 0 && j > 0) { 00719 pos = make_les_entry_2d(i, j, -1, -1, count, pos, les, spvect, 00720 cell_count, status, start_val, 00721 items->NW, cell_type); 00722 } 00723 /* north-eastern neighbour, entry col + 1 row - 1 */ 00724 if (i < geom->cols - 1 && j > 0) { 00725 pos = make_les_entry_2d(i, j, 1, -1, count, pos, les, spvect, 00726 cell_count, status, start_val, 00727 items->NE, cell_type); 00728 } 00729 /* south-western neighbour, entry is col - 1 row + 1 */ 00730 if (i > 0 && j < geom->rows - 1) { 00731 pos = make_les_entry_2d(i, j, -1, 1, count, pos, les, spvect, 00732 cell_count, status, start_val, 00733 items->SW, cell_type); 00734 } 00735 /* south-eastern neighbour, entry col + 1 row + 1 */ 00736 if (i < geom->cols - 1 && j < geom->rows - 1) { 00737 pos = make_les_entry_2d(i, j, 1, 1, count, pos, les, spvect, 00738 cell_count, status, start_val, 00739 items->SE, cell_type); 00740 } 00741 } 00742 00743 /*How many entries in the les */ 00744 if (les->type == N_SPARSE_LES) { 00745 spvect->cols = pos + 1; 00746 G_math_add_spvector(les->Asp, spvect, count); 00747 } 00748 00749 if (items) 00750 G_free(items); 00751 } 00752 00753 /*release memory */ 00754 N_free_array_2d(cell_count); 00755 00756 for (i = 0; i < cell_type_count; i++) 00757 G_free(index_ij[i]); 00758 00759 G_free(index_ij); 00760 00761 return les; 00762 } 00763 00790 int N_les_integrate_dirichlet_2d(N_les * les, N_geom_data * geom, 00791 N_array_2d * status, N_array_2d * start_val) 00792 { 00793 int rows, cols; 00794 00795 int count = 0; 00796 00797 int i, j, x, y, stat; 00798 00799 double *dvect1; 00800 00801 double *dvect2; 00802 00803 G_debug(2, 00804 "N_les_integrate_dirichlet_2d: integrating the dirichlet boundary condition"); 00805 00806 rows = geom->rows; 00807 cols = geom->cols; 00808 00809 /*we nned to additional vectors */ 00810 dvect1 = (double *)G_calloc(les->cols, sizeof(double)); 00811 dvect2 = (double *)G_calloc(les->cols, sizeof(double)); 00812 00813 /*fill the first one with the x vector data of Dirichlet cells */ 00814 count = 0; 00815 for (y = 0; y < rows; y++) { 00816 for (x = 0; x < cols; x++) { 00817 stat = N_get_array_2d_c_value(status, x, y); 00818 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) { 00819 dvect1[count] = N_get_array_2d_d_value(start_val, x, y); 00820 count++; 00821 } 00822 else if (stat == N_CELL_ACTIVE) { 00823 dvect1[count] = 0.0; 00824 count++; 00825 } 00826 } 00827 } 00828 00829 #pragma omp parallel default(shared) 00830 { 00831 /*performe the matrix vector product and */ 00832 if (les->type == N_SPARSE_LES) 00833 G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows); 00834 else 00835 G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols); 00836 #pragma omp for schedule (static) private(i) 00837 for (i = 0; i < les->cols; i++) 00838 les->b[i] = les->b[i] - dvect2[i]; 00839 } 00840 00841 /*now set the Dirichlet cell rows and cols to zero and the 00842 * diagonal entry to 1*/ 00843 count = 0; 00844 for (y = 0; y < rows; y++) { 00845 for (x = 0; x < cols; x++) { 00846 stat = N_get_array_2d_c_value(status, x, y); 00847 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) { 00848 if (les->type == N_SPARSE_LES) { 00849 /*set the rows to zero */ 00850 for (i = 0; i < les->Asp[count]->cols; i++) 00851 les->Asp[count]->values[i] = 0.0; 00852 /*set the cols to zero */ 00853 for (i = 0; i < les->rows; i++) { 00854 for (j = 0; j < les->Asp[i]->cols; j++) { 00855 if (les->Asp[i]->index[j] == count) 00856 les->Asp[i]->values[j] = 0.0; 00857 } 00858 } 00859 00860 /*entry on the diagonal */ 00861 les->Asp[count]->values[0] = 1.0; 00862 00863 } 00864 else { 00865 /*set the rows to zero */ 00866 for (i = 0; i < les->cols; i++) 00867 les->A[count][i] = 0.0; 00868 /*set the cols to zero */ 00869 for (i = 0; i < les->rows; i++) 00870 les->A[i][count] = 0.0; 00871 00872 /*entry on the diagonal */ 00873 les->A[count][count] = 1.0; 00874 } 00875 } 00876 if (stat >= N_CELL_ACTIVE) 00877 count++; 00878 } 00879 } 00880 00881 return 0; 00882 00883 } 00884 00885 /* **************************************************************** */ 00886 /* **** make an entry in the les (2d) ***************************** */ 00887 /* **************************************************************** */ 00888 int make_les_entry_2d(int i, int j, int offset_i, int offset_j, int count, 00889 int pos, N_les * les, G_math_spvector * spvect, 00890 N_array_2d * cell_count, N_array_2d * status, 00891 N_array_2d * start_val, double entry, int cell_type) 00892 { 00893 int K; 00894 00895 int di = offset_i; 00896 00897 int dj = offset_j; 00898 00899 K = N_get_array_2d_c_value(cell_count, i + di, j + dj) - 00900 N_get_array_2d_c_value(cell_count, i, j); 00901 00902 /* active cells build the linear equation system */ 00903 if (cell_type == N_CELL_ACTIVE) { 00904 /* dirichlet or transmission cells must be handled like this */ 00905 if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_ACTIVE && 00906 N_get_array_2d_c_value(status, i + di, j + dj) < N_MAX_CELL_STATE) 00907 les->b[count] -= 00908 N_get_array_2d_d_value(start_val, i + di, j + dj) * entry; 00909 else if (N_get_array_2d_c_value(status, i + di, j + dj) == 00910 N_CELL_ACTIVE) { 00911 if ((count + K) >= 0 && (count + K) < les->cols) { 00912 G_debug(5, 00913 " make_les_entry_2d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n", 00914 count, count + K, entry); 00915 pos++; 00916 if (les->type == N_SPARSE_LES) { 00917 spvect->index[pos] = count + K; 00918 spvect->values[pos] = entry; 00919 } 00920 else { 00921 les->A[count][count + K] = entry; 00922 } 00923 } 00924 } 00925 } /* if dirichlet cells should be used then check for all valid cell neighbours */ 00926 else if (cell_type == N_CELL_DIRICHLET) { 00927 /* all valid cells */ 00928 if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_INACTIVE 00929 && N_get_array_2d_c_value(status, i + di, 00930 j + dj) < N_MAX_CELL_STATE) { 00931 if ((count + K) >= 0 && (count + K) < les->cols) { 00932 G_debug(5, 00933 " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n", 00934 count, count + K, entry); 00935 pos++; 00936 if (les->type == N_SPARSE_LES) { 00937 spvect->index[pos] = count + K; 00938 spvect->values[pos] = entry; 00939 } 00940 else { 00941 les->A[count][count + K] = entry; 00942 } 00943 } 00944 } 00945 } 00946 00947 return pos; 00948 } 00949 00950 00951 /* *************************************************************** * 00952 * ******************** N_assemble_les_3d ************************ * 00953 * *************************************************************** */ 00959 N_les *N_assemble_les_3d(int les_type, N_geom_data * geom, 00960 N_array_3d * status, N_array_3d * start_val, 00961 void *data, N_les_callback_3d * call) 00962 { 00963 return N_assemble_les_3d_param(les_type, geom, status, start_val, data, 00964 call, N_CELL_ACTIVE); 00965 } 00966 00972 N_les *N_assemble_les_3d_active(int les_type, N_geom_data * geom, 00973 N_array_3d * status, N_array_3d * start_val, 00974 void *data, N_les_callback_3d * call) 00975 { 00976 return N_assemble_les_3d_param(les_type, geom, status, start_val, data, 00977 call, N_CELL_ACTIVE); 00978 } 00979 00985 N_les *N_assemble_les_3d_dirichlet(int les_type, N_geom_data * geom, 00986 N_array_3d * status, 00987 N_array_3d * start_val, void *data, 00988 N_les_callback_3d * call) 00989 { 00990 return N_assemble_les_3d_param(les_type, geom, status, start_val, data, 00991 call, N_CELL_DIRICHLET); 00992 } 00993 01022 N_les *N_assemble_les_3d_param(int les_type, N_geom_data * geom, 01023 N_array_3d * status, N_array_3d * start_val, 01024 void *data, N_les_callback_3d * call, 01025 int cell_type) 01026 { 01027 int i, j, k, count = 0, pos = 0; 01028 01029 int cell_type_count = 0; 01030 01031 N_array_3d *cell_count; 01032 01033 N_les *les = NULL; 01034 01035 int **index_ij; 01036 01037 G_debug(2, 01038 "N_assemble_les_3d: starting to assemble the linear equation system"); 01039 01040 cell_count = 01041 N_alloc_array_3d(geom->cols, geom->rows, geom->depths, 1, DCELL_TYPE); 01042 01043 /* First count the number of valid cells and save 01044 * each number in a new 3d array. Those numbers are used 01045 * to create the linear equation system.*/ 01046 01047 if (cell_type == N_CELL_DIRICHLET) { 01048 /* include dirichlet cells in the les */ 01049 for (k = 0; k < geom->depths; k++) { 01050 for (j = 0; j < geom->rows; j++) { 01051 for (i = 0; i < geom->cols; i++) { 01052 /*use all non-inactive cells for les creation */ 01053 if (N_CELL_INACTIVE < 01054 (int)N_get_array_3d_d_value(status, i, j, k) && 01055 (int)N_get_array_3d_d_value(status, i, j, 01056 k) < N_MAX_CELL_STATE) 01057 cell_type_count++; 01058 } 01059 } 01060 } 01061 } 01062 else { 01063 /*use only active cell in the les */ 01064 for (k = 0; k < geom->depths; k++) { 01065 for (j = 0; j < geom->rows; j++) { 01066 for (i = 0; i < geom->cols; i++) { 01067 /*count only active cells */ 01068 if (N_CELL_ACTIVE 01069 == (int)N_get_array_3d_d_value(status, i, j, k)) 01070 cell_type_count++; 01071 01072 } 01073 } 01074 } 01075 } 01076 01077 G_debug(2, 01078 "N_assemble_les_3d: number of used cells %i\n", cell_type_count); 01079 01080 if (cell_type_count == 0.0) 01081 G_fatal_error 01082 ("Not enough active cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.", 01083 cell_type_count); 01084 01085 /* allocate the memory for the linear equation system (les). 01086 * Only valid cells are used to create the les. */ 01087 les = N_alloc_les_Ax_b(cell_type_count, les_type); 01088 01089 index_ij = (int **)G_calloc(cell_type_count, sizeof(int *)); 01090 for (i = 0; i < cell_type_count; i++) 01091 index_ij[i] = (int *)G_calloc(3, sizeof(int)); 01092 01093 count = 0; 01094 /*count the number of cells which should be used to create the linear equation system */ 01095 /*save the k, i and j indices and create a ordered numbering */ 01096 for (k = 0; k < geom->depths; k++) { 01097 for (j = 0; j < geom->rows; j++) { 01098 for (i = 0; i < geom->cols; i++) { 01099 if (cell_type == N_CELL_DIRICHLET) { 01100 if (N_CELL_INACTIVE < 01101 (int)N_get_array_3d_d_value(status, i, j, k) && 01102 (int)N_get_array_3d_d_value(status, i, j, 01103 k) < N_MAX_CELL_STATE) { 01104 N_put_array_3d_d_value(cell_count, i, j, k, count); 01105 index_ij[count][0] = i; 01106 index_ij[count][1] = j; 01107 index_ij[count][2] = k; 01108 count++; 01109 G_debug(5, 01110 "N_assemble_les_3d: non-inactive cells count %i at pos x[%i] y[%i] z[%i]\n", 01111 count, i, j, k); 01112 } 01113 } 01114 else if (N_CELL_ACTIVE == 01115 (int)N_get_array_3d_d_value(status, i, j, k)) { 01116 N_put_array_3d_d_value(cell_count, i, j, k, count); 01117 index_ij[count][0] = i; 01118 index_ij[count][1] = j; 01119 index_ij[count][2] = k; 01120 count++; 01121 G_debug(5, 01122 "N_assemble_les_3d: active cells count %i at pos x[%i] y[%i] z[%i]\n", 01123 count, i, j, k); 01124 } 01125 } 01126 } 01127 } 01128 01129 G_debug(2, "N_assemble_les_3d: starting the parallel assemble loop"); 01130 01131 #pragma omp parallel for private(i, j, k, pos, count) schedule(static) 01132 for (count = 0; count < cell_type_count; count++) { 01133 i = index_ij[count][0]; 01134 j = index_ij[count][1]; 01135 k = index_ij[count][2]; 01136 01137 /*create the entries for the */ 01138 N_data_star *items = call->callback(data, geom, i, j, k); 01139 01140 G_math_spvector *spvect = NULL; 01141 01142 /*allocate a sprase vector */ 01143 if (les_type == N_SPARSE_LES) 01144 spvect = G_math_alloc_spvector(items->count); 01145 /* initial conditions */ 01146 01147 les->x[count] = N_get_array_3d_d_value(start_val, i, j, k); 01148 01149 /* the entry in the vector b */ 01150 les->b[count] = items->V; 01151 01152 /* pos describes the position in the sparse vector. 01153 * the first entry is always the diagonal entry of the matrix*/ 01154 pos = 0; 01155 01156 if (les_type == N_SPARSE_LES) { 01157 spvect->index[pos] = count; 01158 spvect->values[pos] = items->C; 01159 } 01160 else { 01161 les->A[count][count] = items->C; 01162 } 01163 /* western neighbour, entry is col - 1 */ 01164 if (i > 0) { 01165 pos = 01166 make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect, 01167 cell_count, status, start_val, items->W, 01168 cell_type); 01169 } 01170 /* eastern neighbour, entry col + 1 */ 01171 if (i < geom->cols - 1) { 01172 pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect, 01173 cell_count, status, start_val, items->E, 01174 cell_type); 01175 } 01176 /* northern neighbour, entry row -1 */ 01177 if (j > 0) { 01178 pos = 01179 make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect, 01180 cell_count, status, start_val, items->N, 01181 cell_type); 01182 } 01183 /* southern neighbour, entry row +1 */ 01184 if (j < geom->rows - 1) { 01185 pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect, 01186 cell_count, status, start_val, items->S, 01187 cell_type); 01188 } 01189 /*only for a 7 star entry needed */ 01190 if (items->type == N_7_POINT_STAR || items->type == N_27_POINT_STAR) { 01191 /* the upper cell (top), entry depth + 1 */ 01192 if (k < geom->depths - 1) { 01193 pos = 01194 make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les, 01195 spvect, cell_count, status, start_val, 01196 items->T, cell_type); 01197 } 01198 /* the lower cell (bottom), entry depth - 1 */ 01199 if (k > 0) { 01200 pos = 01201 make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les, 01202 spvect, cell_count, status, start_val, 01203 items->B, cell_type); 01204 } 01205 } 01206 01207 /*How many entries in the les */ 01208 if (les->type == N_SPARSE_LES) { 01209 spvect->cols = pos + 1; 01210 G_math_add_spvector(les->Asp, spvect, count); 01211 } 01212 01213 if (items) 01214 G_free(items); 01215 } 01216 01217 N_free_array_3d(cell_count); 01218 01219 for (i = 0; i < cell_type_count; i++) 01220 G_free(index_ij[i]); 01221 01222 G_free(index_ij); 01223 01224 return les; 01225 } 01226 01253 int N_les_integrate_dirichlet_3d(N_les * les, N_geom_data * geom, 01254 N_array_3d * status, N_array_3d * start_val) 01255 { 01256 int rows, cols, depths; 01257 01258 int count = 0; 01259 01260 int i, j, x, y, z, stat; 01261 01262 double *dvect1; 01263 01264 double *dvect2; 01265 01266 G_debug(2, 01267 "N_les_integrate_dirichlet_3d: integrating the dirichlet boundary condition"); 01268 01269 rows = geom->rows; 01270 cols = geom->cols; 01271 depths = geom->depths; 01272 01273 /*we nned to additional vectors */ 01274 dvect1 = (double *)G_calloc(les->cols, sizeof(double)); 01275 dvect2 = (double *)G_calloc(les->cols, sizeof(double)); 01276 01277 /*fill the first one with the x vector data of Dirichlet cells */ 01278 count = 0; 01279 for (z = 0; z < depths; z++) { 01280 for (y = 0; y < rows; y++) { 01281 for (x = 0; x < cols; x++) { 01282 stat = (int)N_get_array_3d_d_value(status, x, y, z); 01283 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) { 01284 dvect1[count] = 01285 N_get_array_3d_d_value(start_val, x, y, z); 01286 count++; 01287 } 01288 else if (stat == N_CELL_ACTIVE) { 01289 dvect1[count] = 0.0; 01290 count++; 01291 } 01292 } 01293 } 01294 } 01295 01296 #pragma omp parallel default(shared) 01297 { 01298 /*performe the matrix vector product and */ 01299 if (les->type == N_SPARSE_LES) 01300 G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows); 01301 else 01302 G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols); 01303 #pragma omp for schedule (static) private(i) 01304 for (i = 0; i < les->cols; i++) 01305 les->b[i] = les->b[i] - dvect2[i]; 01306 } 01307 01308 /*now set the Dirichlet cell rows and cols to zero and the 01309 * diagonal entry to 1*/ 01310 count = 0; 01311 for (z = 0; z < depths; z++) { 01312 for (y = 0; y < rows; y++) { 01313 for (x = 0; x < cols; x++) { 01314 stat = (int)N_get_array_3d_d_value(status, x, y, z); 01315 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) { 01316 if (les->type == N_SPARSE_LES) { 01317 /*set the rows to zero */ 01318 for (i = 0; i < les->Asp[count]->cols; i++) 01319 les->Asp[count]->values[i] = 0.0; 01320 /*set the cols to zero */ 01321 for (i = 0; i < les->rows; i++) { 01322 for (j = 0; j < les->Asp[i]->cols; j++) { 01323 if (les->Asp[i]->index[j] == count) 01324 les->Asp[i]->values[j] = 0.0; 01325 } 01326 } 01327 01328 /*entry on the diagonal */ 01329 les->Asp[count]->values[0] = 1.0; 01330 01331 } 01332 else { 01333 /*set the rows to zero */ 01334 for (i = 0; i < les->cols; i++) 01335 les->A[count][i] = 0.0; 01336 /*set the cols to zero */ 01337 for (i = 0; i < les->rows; i++) 01338 les->A[i][count] = 0.0; 01339 01340 /*entry on the diagonal */ 01341 les->A[count][count] = 1.0; 01342 } 01343 } 01344 count++; 01345 } 01346 } 01347 } 01348 01349 return 0; 01350 01351 } 01352 01353 /* **************************************************************** */ 01354 /* **** make an entry in the les (3d) ***************************** */ 01355 /* **************************************************************** */ 01356 int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j, 01357 int offset_k, int count, int pos, N_les * les, 01358 G_math_spvector * spvect, N_array_3d * cell_count, 01359 N_array_3d * status, N_array_3d * start_val, 01360 double entry, int cell_type) 01361 { 01362 int K; 01363 01364 int di = offset_i; 01365 01366 int dj = offset_j; 01367 01368 int dk = offset_k; 01369 01370 K = (int)N_get_array_3d_d_value(cell_count, i + di, j + dj, k + dk) - 01371 (int)N_get_array_3d_d_value(cell_count, i, j, k); 01372 01373 if (cell_type == N_CELL_ACTIVE) { 01374 if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) > 01375 N_CELL_ACTIVE && 01376 (int)N_get_array_3d_d_value(status, i + di, j + dj, 01377 k + dk) < N_MAX_CELL_STATE) 01378 les->b[count] -= 01379 N_get_array_3d_d_value(start_val, i + di, j + dj, 01380 k + dk) * entry; 01381 else if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) 01382 == N_CELL_ACTIVE) { 01383 if ((count + K) >= 0 && (count + K) < les->cols) { 01384 G_debug(5, 01385 " make_les_entry_3d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n", 01386 count, count + K, entry); 01387 pos++; 01388 if (les->type == N_SPARSE_LES) { 01389 spvect->index[pos] = count + K; 01390 spvect->values[pos] = entry; 01391 } 01392 else { 01393 les->A[count][count + K] = entry; 01394 } 01395 } 01396 } 01397 } 01398 else if (cell_type == N_CELL_DIRICHLET) { 01399 if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) 01400 != N_CELL_INACTIVE) { 01401 if ((count + K) >= 0 && (count + K) < les->cols) { 01402 G_debug(5, 01403 " make_les_entry_3d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n", 01404 count, count + K, entry); 01405 pos++; 01406 if (les->type == N_SPARSE_LES) { 01407 spvect->index[pos] = count + K; 01408 spvect->values[pos] = entry; 01409 } 01410 else { 01411 les->A[count][count + K] = entry; 01412 } 01413 } 01414 } 01415 } 01416 01417 return pos; 01418 }
1.8.0 
