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