21 #include <grass/N_pde.h>    24 static int make_les_entry_2d(
int i, 
int j, 
int offset_i, 
int offset_j,
    31 static int make_les_entry_3d(
int i, 
int j, 
int k, 
int offset_i, 
int offset_j,
   139     G_debug(5, 
"N_create_5star:  w %g e %g n %g s %g c %g v %g\n", star->
W,
   140             star->
E, star->
N, star->
S, star->
C, star->
V);
   162                             double S, 
double T, 
double B, 
double V)
   177     G_debug(5, 
"N_create_7star:  w %g e %g n %g s %g t %g b %g c %g v %g\n",
   178             star->
W, star->
E, star->
N, star->
S, star->
T, star->
B, star->
C,
   203                             double S, 
double NW, 
double SW, 
double NE,
   222             "N_create_9star:  w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
   223             star->
W, star->
E, star->
N, star->
S, star->
NW, star->
SW, star->
NE,
   224             star->
SE, star->
C, star->
V);
   266                              double NW, 
double SW, 
double NE, 
double SE,
   267                              double T, 
double W_T, 
double E_T, 
double N_T,
   268                              double S_T, 
double NW_T, 
double SW_T,
   269                              double NE_T, 
double SE_T, 
double B, 
double W_B,
   270                              double E_B, 
double N_B, 
double S_B, 
double NW_B,
   271                              double SW_B, 
double NE_B, 
double SE_B, 
double V)
   311             "N_create_27star:  w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
   312             star->
W, star->
E, star->
N, star->
S, star->
NW, star->
SW, star->
NE,
   313             star->
SE, star->
C, star->
V);
   316             "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",
   321             "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",
   429     star->
E = 1 / geom->
dx;
   430     star->
W = 1 / geom->
dx;
   431     star->
N = 1 / geom->
dy;
   432     star->
S = 1 / geom->
dy;
   433     star->
T = 1 / geom->
dz;
   434     star->
B = 1 / geom->
dz;
   435     star->
C = -1 * (2 / geom->
dx + 2 / geom->
dy + 2 / geom->
dz);
   439             "N_callback_template_3d:  w %g e %g n %g s %g t %g b %g c %g v %g\n",
   440             star->
W, star->
E, star->
N, star->
S, star->
T, star->
B, star->
C,
   468     star->
E = 1 / geom->
dx;
   469     star->
NE = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
   470     star->
SE = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
   471     star->
W = 1 / geom->
dx;
   472     star->
NW = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
   473     star->
SW = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
   474     star->
N = 1 / geom->
dy;
   475     star->
S = 1 / geom->
dy;
   477         -1 * (star->
E + star->
NE + star->
SE + star->
W + star->
NW + star->
SW +
   565     int i, j, 
count = 0, pos = 0;
   566     int cell_type_count = 0;
   572             "N_assemble_les_2d: starting to assemble the linear equation system");
   583         for (j = 0; j < geom->
rows; j++) {
   584             for (i = 0; i < geom->
cols; i++) {
   594         for (j = 0; j < geom->
rows; j++) {
   595             for (i = 0; i < geom->
cols; i++) {
   603     G_debug(2, 
"N_assemble_les_2d: number of used cells %i\n",
   606     if (cell_type_count == 0)
   608             (
"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.",
   613     index_ij = (
int **)
G_calloc(cell_type_count, 
sizeof(
int *));
   614     for (i = 0; i < cell_type_count; i++)
   615         index_ij[i] = (
int *)
G_calloc(2, 
sizeof(
int));
   623     for (j = 0; j < geom->
rows; j++) {
   624         for (i = 0; i < geom->
cols; i++) {
   630                     index_ij[
count][0] = i;
   631                     index_ij[
count][1] = j;
   634                             "N_assemble_les_2d: non-inactive cells count %i at pos x[%i] y[%i]\n",
   641                 index_ij[
count][0] = i;
   642                 index_ij[
count][1] = j;
   645                         "N_assemble_les_2d: active cells count %i at pos x[%i] y[%i]\n",
   651     G_debug(2, 
"N_assemble_les_2d: starting the parallel assemble loop");
   654 #pragma omp parallel for private(i, j, pos, count) schedule(static)   655     for (count = 0; count < cell_type_count; count++) {
   656         i = index_ij[
count][0];
   657         j = index_ij[
count][1];
   681             spvect->
values[pos] = items->
C;
   688             pos = make_les_entry_2d(i, j, -1, 0, count, pos, les, spvect,
   689                                     cell_count, status, start_val, items->
W,
   693         if (i < geom->
cols - 1) {
   694             pos = make_les_entry_2d(i, j, 1, 0, count, pos, les, spvect,
   695                                     cell_count, status, start_val, items->
E,
   701                 make_les_entry_2d(i, j, 0, -1, count, pos, les, spvect,
   702                                   cell_count, status, start_val, items->
N,
   706         if (j < geom->
rows - 1) {
   707             pos = make_les_entry_2d(i, j, 0, 1, count, pos, les, spvect,
   708                                     cell_count, status, start_val, items->
S,
   714             if (i > 0 && j > 0) {
   715                 pos = make_les_entry_2d(i, j, -1, -1, count, pos, les, spvect,
   716                                         cell_count, status, start_val,
   717                                         items->
NW, cell_type);
   720             if (i < geom->
cols - 1 && j > 0) {
   721                 pos = make_les_entry_2d(i, j, 1, -1, count, pos, les, spvect,
   722                                         cell_count, status, start_val,
   723                                         items->
NE, cell_type);
   726             if (i > 0 && j < geom->
rows - 1) {
   727                 pos = make_les_entry_2d(i, j, -1, 1, count, pos, les, spvect,
   728                                         cell_count, status, start_val,
   729                                         items->
SW, cell_type);
   732             if (i < geom->
cols - 1 && j < geom->
rows - 1) {
   733                 pos = make_les_entry_2d(i, j, 1, 1, count, pos, les, spvect,
   734                                         cell_count, status, start_val,
   735                                         items->
SE, cell_type);
   741             spvect->
cols = pos + 1;
   752     for (i = 0; i < cell_type_count; i++)
   791     int i, j, 
x, y, stat;
   796             "N_les_integrate_dirichlet_2d: integrating the dirichlet boundary condition");
   807     for (y = 0; y < 
rows; y++) {
   808         for (x = 0; x < 
cols; x++) {
   821 #pragma omp parallel default(shared)   828 #pragma omp for schedule (static) private(i)   829         for (i = 0; i < les->
cols; i++)
   830             les->
b[i] = les->
b[i] - dvect2[i];
   836     for (y = 0; y < 
rows; y++) {
   837         for (x = 0; x < 
cols; x++) {
   845                     for (i = 0; i < les->
rows; i++) {
   846                         for (j = 0; j < les->
Asp[i]->
cols; j++) {
   858                     for (i = 0; i < les->
cols; i++)
   859                         les->
A[count][i] = 0.0;
   861                     for (i = 0; i < les->
rows; i++)
   862                         les->
A[i][count] = 0.0;
   865                     les->
A[count][count] = 1.0;
   880 int make_les_entry_2d(
int i, 
int j, 
int offset_i, 
int offset_j, 
int count,
   883                       N_array_2d * start_val, 
double entry, 
int cell_type)
   901             if ((count + K) >= 0 && (count + K) < les->
cols) {
   903                         " make_les_entry_2d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
   904                         count, count + K, entry);
   907                     spvect->
index[pos] = count + K;
   908                     spvect->
values[pos] = entry;
   911                     les->
A[
count][count + K] = entry;
   921             if ((count + K) >= 0 && (count + K) < les->
cols) {
   923                         " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
   924                         count, count + K, entry);
   927                     spvect->
index[pos] = count + K;
   928                     spvect->
values[pos] = entry;
   931                     les->
A[
count][count + K] = entry;
  1017     int i, j, k, count = 0, pos = 0;
  1018     int cell_type_count = 0;
  1024             "N_assemble_les_3d: starting to assemble the linear equation system");
  1035         for (k = 0; k < geom->
depths; k++) {
  1036             for (j = 0; j < geom->
rows; j++) {
  1037                 for (i = 0; i < geom->
cols; i++) {
  1050         for (k = 0; k < geom->
depths; k++) {
  1051             for (j = 0; j < geom->
rows; j++) {
  1052                 for (i = 0; i < geom->
cols; i++) {
  1064             "N_assemble_les_3d: number of  used cells %i\n", cell_type_count);
  1066     if (cell_type_count == 0.0)
  1068             (
"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.",
  1075     index_ij = (
int **)
G_calloc(cell_type_count, 
sizeof(
int *));
  1076     for (i = 0; i < cell_type_count; i++)
  1077         index_ij[i] = (
int *)
G_calloc(3, 
sizeof(
int));
  1082     for (k = 0; k < geom->
depths; k++) {
  1083         for (j = 0; j < geom->
rows; j++) {
  1084             for (i = 0; i < geom->
cols; i++) {
  1091                         index_ij[
count][0] = i;
  1092                         index_ij[
count][1] = j;
  1093                         index_ij[
count][2] = k;
  1096                                 "N_assemble_les_3d: non-inactive cells count %i at pos x[%i] y[%i] z[%i]\n",
  1103                     index_ij[
count][0] = i;
  1104                     index_ij[
count][1] = j;
  1105                     index_ij[
count][2] = k;
  1108                             "N_assemble_les_3d: active cells count %i at pos x[%i] y[%i] z[%i]\n",
  1115     G_debug(2, 
"N_assemble_les_3d: starting the parallel assemble loop");
  1117 #pragma omp parallel for private(i, j, k, pos, count) schedule(static)  1118     for (count = 0; count < cell_type_count; count++) {
  1119         i = index_ij[
count][0];
  1120         j = index_ij[
count][1];
  1121         k = index_ij[
count][2];
  1144             spvect->
values[pos] = items->
C;
  1152                 make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect,
  1153                                   cell_count, status, start_val, items->
W,
  1157         if (i < geom->
cols - 1) {
  1158             pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect,
  1159                                     cell_count, status, start_val, items->
E,
  1165                 make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect,
  1166                                   cell_count, status, start_val, items->
N,
  1170         if (j < geom->
rows - 1) {
  1171             pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect,
  1172                                     cell_count, status, start_val, items->
S,
  1178             if (k < geom->
depths - 1) {
  1180                     make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les,
  1181                                       spvect, cell_count, status, start_val,
  1182                                       items->
T, cell_type);
  1187                     make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les,
  1188                                       spvect, cell_count, status, start_val,
  1189                                       items->
B, cell_type);
  1195             spvect->
cols = pos + 1;
  1205     for (i = 0; i < cell_type_count; i++)
  1244     int i, j, 
x, y, z, stat;
  1249             "N_les_integrate_dirichlet_3d: integrating the dirichlet boundary condition");
  1256     dvect1 = (
double *)
G_calloc(les->
cols, 
sizeof(
double));
  1257     dvect2 = (
double *)
G_calloc(les->
cols, 
sizeof(
double));
  1261     for (z = 0; z < 
depths; z++) {
  1262         for (y = 0; y < 
rows; y++) {
  1263             for (x = 0; x < 
cols; x++) {
  1271                     dvect1[
count] = 0.0;
  1278 #pragma omp parallel default(shared)  1285 #pragma omp for schedule (static) private(i)  1286         for (i = 0; i < les->
cols; i++)
  1287             les->
b[i] = les->
b[i] - dvect2[i];
  1293     for (z = 0; z < 
depths; z++) {
  1294         for (y = 0; y < 
rows; y++) {
  1295             for (x = 0; x < 
cols; x++) {
  1303                         for (i = 0; i < les->
rows; i++) {
  1304                             for (j = 0; j < les->
Asp[i]->
cols; j++) {
  1305                                 if (les->
Asp[i]->
index[j] == count)
  1316                         for (i = 0; i < les->
cols; i++)
  1317                             les->
A[count][i] = 0.0;
  1319                         for (i = 0; i < les->
rows; i++)
  1320                             les->
A[i][count] = 0.0;
  1323                         les->
A[count][count] = 1.0;
  1338 int make_les_entry_3d(
int i, 
int j, 
int k, 
int offset_i, 
int offset_j,
  1339                       int offset_k, 
int count, 
int pos, 
N_les * les,
  1342                       double entry, 
int cell_type)
  1362             if ((count + K) >= 0 && (count + K) < les->
cols) {
  1364                         " make_les_entry_3d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
  1365                         count, count + K, entry);
  1368                     spvect->
index[pos] = count + K;
  1369                     spvect->
values[pos] = entry;
  1372                     les->
A[
count][count + K] = entry;
  1380             if ((count + K) >= 0 && (count + K) < les->
cols) {
  1382                         " make_les_entry_3d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
  1383                         count, count + K, entry);
  1386                     spvect->
index[pos] = count + K;
  1387                     spvect->
values[pos] = entry;
  1390                     les->
A[
count][count + K] = entry;
 
N_data_star *(* callback)()
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
N_data_star * N_create_9star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double V)
allocate and initialize a 9 point star data structure 
Matrix entries for a mass balance 5/7/9 star system. 
N_les * N_alloc_les_Ax_b(int rows, int type)
Allocate memory for a quadratic linear equation system which includes the Matrix A, vector x and vector b. 
N_les * N_assemble_les_2d_dirichlet(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active and dirichlet c...
G_math_spvector * G_math_alloc_spvector(int)
Allocate memory for a sparse vector. 
The row vector of the sparse matrix. 
int N_les_integrate_dirichlet_2d(N_les *les, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val)
Integrate Dirichlet or Transmission boundary conditions into the les (2s) 
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth. 
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure. 
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure. 
N_les_callback_3d * N_alloc_les_callback_3d(void)
Allocate the structure holding the callback function. 
N_data_star * N_alloc_5star(void)
allocate a 5 point star data structure 
void G_free(void *)
Free allocated memory. 
#define N_MAX_CELL_STATE
the maximum number of available cell states (eg: boundary condition, inactiven active) ...
int N_les_integrate_dirichlet_3d(N_les *les, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val)
Integrate Dirichlet or Transmission boundary conditions into the les (3d) 
N_les * N_assemble_les_2d(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active cells...
void N_put_array_2d_c_value(N_array_2d *data, int col, int row, CELL value)
Writes a CELL value to the N_array_2d struct at position col, row. 
callback structure for 2d matrix assembling 
N_data_star * N_callback_template_2d(void *data, N_geom_data *geom, int col, int row)
A callback template creates a 9 point star structure. 
void G_math_Ax_sparse(G_math_spvector **, double *, double *, int)
Compute the matrix - vector product of sparse matrix **Asp and vector x. 
N_les * N_assemble_les_2d_param(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call, int cell_type)
Assemble a linear equation system (les) based on 2d location data (raster) 
N_les * N_assemble_les_3d_dirichlet(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active and dirichlet cells...
Geometric information about the structured grid. 
N_les * N_assemble_les_3d(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active cells. 
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure. 
N_data_star * N_alloc_9star(void)
allocate a 9 point star data structure 
N_les_callback_2d * N_alloc_les_callback_2d(void)
Allocate the structure holding the callback function. 
N_data_star * N_create_27star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double T, double W_T, double E_T, double N_T, double S_T, double NW_T, double SW_T, double NE_T, double SE_T, double B, double W_B, double E_B, double N_B, double S_B, double NW_B, double SW_B, double NE_B, double SE_B, double V)
allocate and initialize a 27 point star data structure 
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth. 
N_data_star * N_alloc_7star(void)
allocate a 7 point star data structure 
void N_set_les_callback_2d_func(N_les_callback_2d *data, N_data_star *(*callback_func_2d)())
Set the callback function which is called while assembling the les in 2d. 
N_data_star * N_create_7star(double C, double W, double E, double N, double S, double T, double B, double V)
allocate and initialize a 7 point star data structure 
N_data_star * N_alloc_27star(void)
allocate a 27 point star data structure 
int depths
number of depths for 3D data 
N_les * N_assemble_les_3d_active(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active cells. 
int cols
Number of columns for 2D data. 
CELL N_get_array_2d_c_value(N_array_2d *data, int col, int row)
Returns the value of type CELL at position col, row. 
void N_set_les_callback_3d_func(N_les_callback_3d *data, N_data_star *(*callback_func_3d)())
Set the callback function which is called while assembling the les in 3d. 
N_data_star *(* callback)()
N_les * N_assemble_les_2d_active(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active cells...
N_les * N_assemble_les_3d_param(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call, int cell_type)
Assemble a linear equation system (les) based on 3d location data (g3d) 
void G_math_d_Ax(double **, double *, double *, int, int)
Compute the matrix - vector product of matrix A and vector x. 
callback structure for 3d matrix assembling 
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row. 
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d. 
int rows
Number of rows for 2D data. 
int G_math_add_spvector(G_math_spvector **, G_math_spvector *, int)
Adds a sparse vector to a sparse matrix at position row. 
int G_debug(int, const char *,...) __attribute__((format(printf
N_data_star * N_create_5star(double C, double W, double E, double N, double S, double V)
allocate and initialize a 5 point star data structure 
N_data_star * N_callback_template_3d(void *data, N_geom_data *geom, int col, int row, int depth)
A callback template creates a 7 point star structure. 
The linear equation system (les) structure.