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