GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
n_les_assemble.c
Go to the documentation of this file.
1/*****************************************************************************
2 *
3 * MODULE: Grass PDE Numerical Library
4 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
5 * soerengebbert <at> gmx <dot> de
6 *
7 * PURPOSE: functions to assemble a linear equation system
8 * part of the gpde library
9 *
10 * COPYRIGHT: (C) 2000 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General Public
13 * License (>=v2). Read the file COPYING that comes with GRASS
14 * for details.
15 *
16 *****************************************************************************/
17
18#include <math.h>
19#include <grass/N_pde.h>
20
21/* local protos */
22static int make_les_entry_2d(int i, int j, int offset_i, int offset_j,
23 int count, int pos, N_les *les,
24 G_math_spvector *spvect, N_array_2d *cell_count,
25 N_array_2d *status, N_array_2d *start_val,
26 double entry, int cell_type);
27
28static int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
29 int offset_k, int count, int pos, N_les *les,
30 G_math_spvector *spvect, N_array_3d *cell_count,
31 N_array_3d *status, N_array_3d *start_val,
32 double entry, int cell_type);
33
34/* *************************************************************** *
35 * ********************** N_alloc_5star ************************** *
36 * *************************************************************** */
37/*!
38 * \brief allocate a 5 point star data structure
39 *
40 * \return N_data_star *
41 * */
43{
44 N_data_star *star = (N_data_star *)G_calloc(1, sizeof(N_data_star));
45
46 star->type = N_5_POINT_STAR;
47 star->count = 5;
48 return star;
49}
50
51/* *************************************************************** *
52 * ********************* N_alloc_7star *************************** *
53 * *************************************************************** */
54/*!
55 * \brief allocate a 7 point star data structure
56 *
57 * \return N_data_star *
58 * */
60{
61 N_data_star *star = (N_data_star *)G_calloc(1, sizeof(N_data_star));
62
63 star->type = N_7_POINT_STAR;
64 star->count = 7;
65 return star;
66}
67
68/* *************************************************************** *
69 * ********************* N_alloc_9star *************************** *
70 * *************************************************************** */
71/*!
72 * \brief allocate a 9 point star data structure
73 *
74 * \return N_data_star *
75 *
76 * \attention The 9 point start is not yet implemented in the matrix assembling
77 * function
78 *
79 * */
81{
82 N_data_star *star = (N_data_star *)G_calloc(1, sizeof(N_data_star));
83
84 star->type = N_9_POINT_STAR;
85 star->count = 9;
86 return star;
87}
88
89/* *************************************************************** *
90 * ********************* N_alloc_27star ************************** *
91 * *************************************************************** */
92/*!
93 * \brief allocate a 27 point star data structure
94 *
95 * \return N_data_star *
96 *
97 * \attention The 27 point start is not yet implemented in the matrix assembling
98 * function
99 *
100 * */
102{
103 N_data_star *star = (N_data_star *)G_calloc(1, sizeof(N_data_star));
104
105 star->type = N_27_POINT_STAR;
106 star->count = 27;
107 return star;
108}
109
110/* *************************************************************** *
111 * ********************** N_create_5star ************************* *
112 * *************************************************************** */
113/*!
114 * \brief allocate and initialize a 5 point star data structure
115 *
116 * \param C double
117 * \param W double
118 * \param E double
119 * \param N double
120 * \param S double
121 * \param V double
122 * \return N_data_star *
123 * */
124N_data_star *N_create_5star(double C, double W, double E, double N, double S,
125 double V)
126{
127 N_data_star *star = N_alloc_5star();
128
129 star->C = C;
130 star->W = W;
131 star->E = E;
132 star->N = N;
133 star->S = S;
134
135 star->V = V;
136
137 G_debug(5, "N_create_5star: w %g e %g n %g s %g c %g v %g\n", star->W,
138 star->E, star->N, star->S, star->C, star->V);
139
140 return star;
141}
142
143/* *************************************************************** *
144 * ************************* N_create_7star ********************** *
145 * *************************************************************** */
146/*!
147 * \brief allocate and initialize a 7 point star data structure
148 *
149 * \param C double
150 * \param W double
151 * \param E double
152 * \param N double
153 * \param S double
154 * \param T double
155 * \param B double
156 * \param V double
157 * \return N_data_star *
158 * */
159N_data_star *N_create_7star(double C, double W, double E, double N, double S,
160 double T, double B, double V)
161{
162 N_data_star *star = N_alloc_7star();
163
164 star->C = C;
165 star->W = W;
166 star->E = E;
167 star->N = N;
168 star->S = S;
169
170 star->T = T;
171 star->B = B;
172
173 star->V = V;
174
175 G_debug(5, "N_create_7star: w %g e %g n %g s %g t %g b %g c %g v %g\n",
176 star->W, star->E, star->N, star->S, star->T, star->B, star->C,
177 star->V);
178
179 return star;
180}
181
182/* *************************************************************** *
183 * ************************ N_create_9star *********************** *
184 * *************************************************************** */
185/*!
186 * \brief allocate and initialize a 9 point star data structure
187 *
188 * \param C double
189 * \param W double
190 * \param E double
191 * \param N double
192 * \param S double
193 * \param NW double
194 * \param SW double
195 * \param NE double
196 * \param SE double
197 * \param V double
198 * \return N_data_star *
199 * */
200N_data_star *N_create_9star(double C, double W, double E, double N, double S,
201 double NW, double SW, double NE, double SE,
202 double V)
203{
204 N_data_star *star = N_alloc_9star();
205
206 star->C = C;
207 star->W = W;
208 star->E = E;
209 star->N = N;
210 star->S = S;
211
212 star->NW = NW;
213 star->SW = SW;
214 star->NE = NE;
215 star->SE = SE;
216
217 star->V = V;
218
219 G_debug(5,
220 "N_create_9star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g "
221 "v %g\n",
222 star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE,
223 star->SE, star->C, star->V);
224
225 return star;
226}
227
228/* *************************************************************** *
229 * ************************ N_create_27star *********************** *
230 * *************************************************************** */
231/*!
232 * \brief allocate and initialize a 27 point star data structure
233 *
234 * \param C double
235 * \param W double
236 * \param E double
237 * \param N double
238 * \param S double
239 * \param NW double
240 * \param SW double
241 * \param NE double
242 * \param SE double
243 * \param T double
244 * \param W_T double
245 * \param E_T double
246 * \param N_T double
247 * \param S_T double
248 * \param NW_T double
249 * \param SW_T double
250 * \param NE_T double
251 * \param SE_T double
252 * \param B double
253 * \param W_B double
254 * \param E_B double
255 * \param N_B double
256 * \param S_B double
257 * \param NW_B double
258 * \param SW_B double
259 * \param NE_B double
260 * \param SE_B double
261 * \param V double
262 * \return N_data_star *
263 * */
264N_data_star *N_create_27star(double C, double W, double E, double N, double S,
265 double NW, double SW, double NE, double SE,
266 double T, double W_T, double E_T, double N_T,
267 double S_T, double NW_T, double SW_T, double NE_T,
268 double SE_T, double B, double W_B, double E_B,
269 double N_B, double S_B, double NW_B, double SW_B,
270 double NE_B, double SE_B, double V)
271{
272 N_data_star *star = N_alloc_27star();
273
274 star->C = C;
275 star->W = W;
276 star->E = E;
277 star->N = N;
278 star->S = S;
279
280 star->NW = NW;
281 star->SW = SW;
282 star->NE = NE;
283 star->SE = SE;
284
285 star->T = T;
286 star->W_T = W_T;
287 star->E_T = E_T;
288 star->N_T = N_T;
289 star->S_T = S_T;
290
291 star->NW_T = NW_T;
292 star->SW_T = SW_T;
293 star->NE_T = NE_T;
294 star->SE_T = SE_T;
295
296 star->B = B;
297 star->W_B = W_B;
298 star->E_B = E_B;
299 star->N_B = N_B;
300 star->S_B = S_B;
301
302 star->NW_B = NW_B;
303 star->SW_B = SW_B;
304 star->NE_B = NE_B;
305 star->SE_B = SE_B;
306
307 star->V = V;
308
309 G_debug(5,
310 "N_create_27star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c "
311 "%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);
314
315 G_debug(5,
316 "N_create_27star: w_t %g e_t %g n_t %g s_t %g nw_t %g sw_t %g "
317 "ne_t %g se_t %g t %g \n",
318 star->W_T, star->E_T, star->N_T, star->S_T, star->NW_T, star->SW_T,
319 star->NE_T, star->SE_T, star->T);
320
321 G_debug(5,
322 "N_create_27star: w_b %g e_b %g n_b %g s_b %g nw_b %g sw_b %g "
323 "ne_b %g se_B %g b %g\n",
324 star->W_B, star->E_B, star->N_B, star->S_B, star->NW_B, star->SW_B,
325 star->NE_B, star->SE_B, star->B);
326
327 return star;
328}
329
330/* *************************************************************** *
331 * ****************** N_set_les_callback_3d_func ***************** *
332 * *************************************************************** */
333/*!
334 * \brief Set the callback function which is called while assembling the les in
335 * 3d
336 *
337 * \param data N_les_callback_3d *
338 * \param callback_func_3d N_data_star *
339 * \return void
340 * */
342 N_data_star *(*callback_func_3d)(void *,
343 N_geom_data *,
344 int, int, int))
345{
346 data->callback = callback_func_3d;
347}
348
349/* *************************************************************** *
350 * *************** N_set_les_callback_2d_func ******************** *
351 * *************************************************************** */
352/*!
353 * \brief Set the callback function which is called while assembling the les in
354 * 2d
355 *
356 * \param data N_les_callback_2d *
357 * \param callback_func_2d N_data_star *
358 * \return void
359 * */
361 N_data_star *(*callback_func_2d)(void *,
362 N_geom_data *,
363 int, int))
364{
365 data->callback = callback_func_2d;
366}
367
368/* *************************************************************** *
369 * ************** N_alloc_les_callback_3d ************************ *
370 * *************************************************************** */
371/*!
372 * \brief Allocate the structure holding the callback function
373 *
374 * A template callback is set. Use N_set_les_callback_3d_func
375 * to set up a specific function.
376 *
377 * \return N_les_callback_3d *
378 * */
380{
381 N_les_callback_3d *call;
382
383 call = (N_les_callback_3d *)G_calloc(1, sizeof(N_les_callback_3d *));
385
386 return call;
387}
388
389/* *************************************************************** *
390 * *************** N_alloc_les_callback_2d *********************** *
391 * *************************************************************** */
392/*!
393 * \brief Allocate the structure holding the callback function
394 *
395 * A template callback is set. Use N_set_les_callback_2d_func
396 * to set up a specific function.
397 *
398 * \return N_les_callback_2d *
399 * */
401{
402 N_les_callback_2d *call;
403
404 call = (N_les_callback_2d *)G_calloc(1, sizeof(N_les_callback_2d *));
406
407 return call;
408}
409
410/* *************************************************************** *
411 * ******************** N_callback_template_3d ******************* *
412 * *************************************************************** */
413/*!
414 * \brief A callback template creates a 7 point star structure
415 *
416 * This is a template callback for mass balance calculation with 7 point stars
417 * based on 3d data (g3d).
418 *
419 * \param data void * (unused)
420 * \param geom N_geom_data *
421 * \param depth int (unused)
422 * \param row int (unused)
423 * \param col int (unused)
424 * \return N_data_star *
425 *
426 * */
428 int col UNUSED, int row UNUSED,
429 int depth UNUSED)
430{
431 N_data_star *star = N_alloc_7star();
432
433 star->E = 1 / geom->dx;
434 star->W = 1 / geom->dx;
435 star->N = 1 / geom->dy;
436 star->S = 1 / geom->dy;
437 star->T = 1 / geom->dz;
438 star->B = 1 / geom->dz;
439 star->C = -1 * (2 / geom->dx + 2 / geom->dy + 2 / geom->dz);
440 star->V = -1;
441
442 G_debug(
443 5, "N_callback_template_3d: w %g e %g n %g s %g t %g b %g c %g v %g\n",
444 star->W, star->E, star->N, star->S, star->T, star->B, star->C, star->V);
445
446 return star;
447}
448
449/* *************************************************************** *
450 * ********************* N_callback_template_2d ****************** *
451 * *************************************************************** */
452/*!
453 * \brief A callback template creates a 9 point star structure
454 *
455 * This is a template callback for mass balance calculation with 9 point stars
456 * based on 2d data (raster).
457 *
458 * \param data void * (unused)
459 * \param geom N_geom_data *
460 * \param row int (unused)
461 * \param col int (unused)
462 * \return N_data_star *
463 *
464 * */
466 int col UNUSED, int row UNUSED)
467{
468 N_data_star *star = N_alloc_9star();
469
470 star->E = 1 / geom->dx;
471 star->NE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
472 star->SE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
473 star->W = 1 / geom->dx;
474 star->NW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
475 star->SW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
476 star->N = 1 / geom->dy;
477 star->S = 1 / geom->dy;
478 star->C = -1 * (star->E + star->NE + star->SE + star->W + star->NW +
479 star->SW + star->N + star->S);
480 star->V = 0;
481
482 return star;
483}
484
485/* *************************************************************** *
486 * ******************** N_assemble_les_2d ************************ *
487 * *************************************************************** */
488/*!
489 * \brief Assemble a linear equation system (les) based on 2d location data
490 * (raster) and active cells
491 *
492 * This function calls #N_assemble_les_2d_param
493 *
494 */
495N_les *N_assemble_les_2d(int les_type, N_geom_data *geom, N_array_2d *status,
496 N_array_2d *start_val, void *data,
497 N_les_callback_2d *call)
498{
499 return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
500 call, N_CELL_ACTIVE);
501}
502
503/*!
504 * \brief Assemble a linear equation system (les) based on 2d location data
505 * (raster) and active cells
506 *
507 * This function calls #N_assemble_les_2d_param
508 *
509 */
511 N_array_2d *status, N_array_2d *start_val,
512 void *data, N_les_callback_2d *call)
513{
514 return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
515 call, N_CELL_ACTIVE);
516}
517
518/*!
519 * \brief Assemble a linear equation system (les) based on 2d location data
520 * (raster) and active and dirichlet cells
521 *
522 * This function calls #N_assemble_les_2d_param
523 *
524 */
526 N_array_2d *status, N_array_2d *start_val,
527 void *data, N_les_callback_2d *call)
528{
529 return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
530 call, N_CELL_DIRICHLET);
531}
532
533/*!
534 * \brief Assemble a linear equation system (les) based on 2d location data
535 * (raster)
536 *
537 *
538 * The linear equation system type can be set to N_NORMAL_LES to create a
539 * regular matrix, or to N_SPARSE_LES to create a sparse matrix. This function
540 * returns a new created linear equation system which can be solved with linear
541 * equation solvers. An 2d array with start values and an 2d status array must
542 * be provided as well as the location geometry and a void pointer to data
543 * passed to the callback which creates the les row entries. This callback
544 * must be defined in the N_les_callback_2d structure.
545 *
546 * The creation of the les is parallelized with OpenMP.
547 * If you implement new callbacks, please make sure that the
548 * function calls are thread safe.
549 *
550 *
551 * the les can be created in two ways, with dirichlet and similar cells and
552 * without them, to spare some memory. If the les is created with dirichlet
553 * cell, the dirichlet boundary condition must be added.
554 *
555 * \param les_type int
556 * \param geom N_geom_data*
557 * \param status N_array_2d *
558 * \param start_val N_array_2d *
559 * \param data void *
560 * \param cell_type int -- les assemble based on N_CELL_ACTIVE or
561 * N_CELL_DIRICHLET \param call N_les_callback_2d * \return N_les *
562 * */
564 N_array_2d *status, N_array_2d *start_val,
565 void *data, N_les_callback_2d *call,
566 int cell_type)
567{
568 int i, j, count = 0, pos = 0;
569 int cell_type_count = 0;
570 int **index_ij;
571 N_array_2d *cell_count;
572 N_les *les = NULL;
573
574 G_debug(
575 2,
576 "N_assemble_les_2d: starting to assemble the linear equation system");
577
578 /* At first count the number of valid cells and save
579 * each number in a new 2d array. Those numbers are used
580 * to create the linear equation system.
581 * */
582
583 cell_count = N_alloc_array_2d(geom->cols, geom->rows, 1, CELL_TYPE);
584
585 /* include dirichlet cells in the les */
586 if (cell_type == N_CELL_DIRICHLET) {
587 for (j = 0; j < geom->rows; j++) {
588 for (i = 0; i < geom->cols; i++) {
589 /*use all non-inactive cells for les creation */
590 if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) &&
592 cell_type_count++;
593 }
594 }
595 }
596 /*use only active cell in the les */
597 if (cell_type == N_CELL_ACTIVE) {
598 for (j = 0; j < geom->rows; j++) {
599 for (i = 0; i < geom->cols; i++) {
600 /*count only active cells */
601 if (N_CELL_ACTIVE == N_get_array_2d_d_value(status, i, j))
602 cell_type_count++;
603 }
604 }
605 }
606
607 G_debug(2, "N_assemble_les_2d: number of used cells %i\n", cell_type_count);
608
609 if (cell_type_count == 0)
610 G_fatal_error("Not enough cells [%i] to create the linear equation "
611 "system. Check the cell status. Only active cells (value "
612 "= 1) are used to create the equation system.",
613 cell_type_count);
614
615 /* Then allocate the memory for the linear equation system (les).
616 * Only valid cells are used to create the les. */
617 index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
618 for (i = 0; i < cell_type_count; i++)
619 index_ij[i] = (int *)G_calloc(2, sizeof(int));
620
621 les = N_alloc_les_Ax_b(cell_type_count, les_type);
622
623 count = 0;
624
625 /*count the number of cells which should be used to create the linear
626 * equation system */
627 /*save the i and j indices and create a ordered numbering */
628 for (j = 0; j < geom->rows; j++) {
629 for (i = 0; i < geom->cols; i++) {
630 /*count every non-inactive cell */
631 if (cell_type == N_CELL_DIRICHLET) {
632 if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) &&
634 N_put_array_2d_c_value(cell_count, i, j, count);
635 index_ij[count][0] = i;
636 index_ij[count][1] = j;
637 count++;
638 G_debug(5,
639 "N_assemble_les_2d: non-inactive cells count %i at "
640 "pos x[%i] y[%i]\n",
641 count, i, j);
642 }
643 /*count every active cell */
644 }
645 else if (N_CELL_ACTIVE == N_get_array_2d_c_value(status, i, j)) {
646 N_put_array_2d_c_value(cell_count, i, j, count);
647 index_ij[count][0] = i;
648 index_ij[count][1] = j;
649 count++;
650 G_debug(5,
651 "N_assemble_les_2d: active cells count %i at pos x[%i] "
652 "y[%i]\n",
653 count, i, j);
654 }
655 }
656 }
657
658 G_debug(2, "N_assemble_les_2d: starting the parallel assemble loop");
659
660 /* Assemble the matrix in parallel */
661#pragma omp parallel for private(i, j, pos, count) schedule(static)
662 for (count = 0; count < cell_type_count; count++) {
663 i = index_ij[count][0];
664 j = index_ij[count][1];
665
666 /*create the entries for the */
667 N_data_star *items = call->callback(data, geom, i, j);
668
669 /* we need a sparse vector pointer anytime */
670 G_math_spvector *spvect = NULL;
671
672 /*allocate a sprase vector */
673 if (les_type == N_SPARSE_LES) {
674 spvect = G_math_alloc_spvector(items->count);
675 }
676 /* initial conditions */
677 les->x[count] = N_get_array_2d_d_value(start_val, i, j);
678
679 /* the entry in the vector b */
680 les->b[count] = items->V;
681
682 /* pos describes the position in the sparse vector.
683 * the first entry is always the diagonal entry of the matrix*/
684 pos = 0;
685
686 if (les_type == N_SPARSE_LES) {
687 spvect->index[pos] = count;
688 spvect->values[pos] = items->C;
689 }
690 else {
691 les->A[count][count] = items->C;
692 }
693 /* western neighbour, entry is col - 1 */
694 if (i > 0) {
695 pos = make_les_entry_2d(i, j, -1, 0, count, pos, les, spvect,
696 cell_count, status, start_val, items->W,
697 cell_type);
698 }
699 /* eastern neighbour, entry col + 1 */
700 if (i < geom->cols - 1) {
701 pos = make_les_entry_2d(i, j, 1, 0, count, pos, les, spvect,
702 cell_count, status, start_val, items->E,
703 cell_type);
704 }
705 /* northern neighbour, entry row - 1 */
706 if (j > 0) {
707 pos = make_les_entry_2d(i, j, 0, -1, count, pos, les, spvect,
708 cell_count, status, start_val, items->N,
709 cell_type);
710 }
711 /* southern neighbour, entry row + 1 */
712 if (j < geom->rows - 1) {
713 pos = make_les_entry_2d(i, j, 0, 1, count, pos, les, spvect,
714 cell_count, status, start_val, items->S,
715 cell_type);
716 }
717 /*in case of a nine point star, we have additional entries */
718 if (items->type == N_9_POINT_STAR) {
719 /* north-western neighbour, entry is col - 1 row - 1 */
720 if (i > 0 && j > 0) {
721 pos = make_les_entry_2d(i, j, -1, -1, count, pos, les, spvect,
722 cell_count, status, start_val,
723 items->NW, cell_type);
724 }
725 /* north-eastern neighbour, entry col + 1 row - 1 */
726 if (i < geom->cols - 1 && j > 0) {
727 pos = make_les_entry_2d(i, j, 1, -1, count, pos, les, spvect,
728 cell_count, status, start_val,
729 items->NE, cell_type);
730 }
731 /* south-western neighbour, entry is col - 1 row + 1 */
732 if (i > 0 && 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->SW, cell_type);
736 }
737 /* south-eastern neighbour, entry col + 1 row + 1 */
738 if (i < geom->cols - 1 && j < geom->rows - 1) {
739 pos = make_les_entry_2d(i, j, 1, 1, count, pos, les, spvect,
740 cell_count, status, start_val,
741 items->SE, cell_type);
742 }
743 }
744
745 /*How many entries in the les */
746 if (les->type == N_SPARSE_LES) {
747 spvect->cols = pos + 1;
748 G_math_add_spvector(les->Asp, spvect, count);
749 }
750
751 if (items)
752 G_free(items);
753 }
754
755 /*release memory */
756 N_free_array_2d(cell_count);
757
758 for (i = 0; i < cell_type_count; i++)
759 G_free(index_ij[i]);
760
761 G_free(index_ij);
762
763 return les;
764}
765
766/*!
767 * \brief Integrate Dirichlet or Transmission boundary conditions into the les
768 * (2s)
769 *
770 * Dirichlet and Transmission boundary conditions will be integrated into
771 * the provided linear equation system. This is meaningful if
772 * the les was created with #N_assemble_les_2d_dirichlet, because in
773 * this case Dirichlet boundary conditions are not automatically included.
774 *
775 * The provided les will be modified:
776 *
777 * Ax = b will be split into Ax_u + Ax_d = b
778 *
779 * x_u - the unknowns
780 * x_d - the Dirichlet cells
781 *
782 * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
783 *
784 * | A_u 0 | x_u
785 * | 0 I | x_d
786 *
787 * \param les N_les* -- the linear equation system
788 * \param geom N_geom_data* -- geometrical data information
789 * \param status N_array_2d* -- the status array containing the cell types
790 * \param start_val N_array_2d* -- an array with start values
791 * \return int -- 1 = success, 0 = failure
792 * */
794 N_array_2d *status, N_array_2d *start_val)
795{
796 int rows, cols;
797 int count = 0;
798 int i, j, x, y, stat;
799 double *dvect1;
800 double *dvect2;
801
802 G_debug(2, "N_les_integrate_dirichlet_2d: integrating the dirichlet "
803 "boundary condition");
804
805 rows = geom->rows;
806 cols = geom->cols;
807
808 /*we need to additional vectors */
809 dvect1 = (double *)G_calloc(les->cols, sizeof(double));
810 dvect2 = (double *)G_calloc(les->cols, sizeof(double));
811
812 /*fill the first one with the x vector data of Dirichlet cells */
813 count = 0;
814 for (y = 0; y < rows; y++) {
815 for (x = 0; x < cols; x++) {
816 stat = N_get_array_2d_c_value(status, x, y);
817 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
818 dvect1[count] = N_get_array_2d_d_value(start_val, x, y);
819 count++;
820 }
821 else if (stat == N_CELL_ACTIVE) {
822 dvect1[count] = 0.0;
823 count++;
824 }
825 }
826 }
827
828#pragma omp parallel default(shared)
829 {
830 /*perform the matrix vector product and */
831 if (les->type == N_SPARSE_LES)
832 G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows);
833 else
834 G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols);
835#pragma omp for schedule(static) private(i)
836 for (i = 0; i < les->cols; i++)
837 les->b[i] = les->b[i] - dvect2[i];
838 }
839
840 /*now set the Dirichlet cell rows and cols to zero and the
841 * diagonal entry to 1*/
842 count = 0;
843 for (y = 0; y < rows; y++) {
844 for (x = 0; x < cols; x++) {
845 stat = N_get_array_2d_c_value(status, x, y);
846 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
847 if (les->type == N_SPARSE_LES) {
848 /*set the rows to zero */
849 for (i = 0; (unsigned int)i < les->Asp[count]->cols; i++)
850 les->Asp[count]->values[i] = 0.0;
851 /*set the cols to zero */
852 for (i = 0; i < les->rows; i++) {
853 for (j = 0; (unsigned int)j < les->Asp[i]->cols; j++) {
854 if (les->Asp[i]->index[j] == (unsigned int)count)
855 les->Asp[i]->values[j] = 0.0;
856 }
857 }
858
859 /*entry on the diagonal */
860 les->Asp[count]->values[0] = 1.0;
861 }
862 else {
863 /*set the rows to zero */
864 for (i = 0; i < les->cols; i++)
865 les->A[count][i] = 0.0;
866 /*set the cols to zero */
867 for (i = 0; i < les->rows; i++)
868 les->A[i][count] = 0.0;
869
870 /*entry on the diagonal */
871 les->A[count][count] = 1.0;
872 }
873 }
874 if (stat >= N_CELL_ACTIVE)
875 count++;
876 }
877 }
878
879 return 0;
880}
881
882/* **************************************************************** */
883/* **** make an entry in the les (2d) ***************************** */
884/* **************************************************************** */
885int make_les_entry_2d(int i, int j, int offset_i, int offset_j, int count,
886 int pos, N_les *les, G_math_spvector *spvect,
887 N_array_2d *cell_count, N_array_2d *status,
888 N_array_2d *start_val, double entry, int cell_type)
889{
890 int K;
891 int di = offset_i;
892 int dj = offset_j;
893
894 K = N_get_array_2d_c_value(cell_count, i + di, j + dj) -
895 N_get_array_2d_c_value(cell_count, i, j);
896
897 /* active cells build the linear equation system */
898 if (cell_type == N_CELL_ACTIVE) {
899 /* dirichlet or transmission cells must be handled like this */
900 if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_ACTIVE &&
901 N_get_array_2d_c_value(status, i + di, j + dj) < N_MAX_CELL_STATE)
902 les->b[count] -=
903 N_get_array_2d_d_value(start_val, i + di, j + dj) * entry;
904 else if (N_get_array_2d_c_value(status, i + di, j + dj) ==
906 if ((count + K) >= 0 && (count + K) < les->cols) {
907 G_debug(5,
908 " make_les_entry_2d: (N_CELL_ACTIVE) create matrix "
909 "entry at row[%i] col[%i] value %g\n",
910 count, count + K, entry);
911 pos++;
912 if (les->type == N_SPARSE_LES) {
913 spvect->index[pos] = count + K;
914 spvect->values[pos] = entry;
915 }
916 else {
917 les->A[count][count + K] = entry;
918 }
919 }
920 }
921 } /* if dirichlet cells should be used then check for all valid cell
922 neighbours */
923 else if (cell_type == N_CELL_DIRICHLET) {
924 /* all valid cells */
925 if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_INACTIVE &&
926 N_get_array_2d_c_value(status, i + di, j + dj) < N_MAX_CELL_STATE) {
927 if ((count + K) >= 0 && (count + K) < les->cols) {
928 G_debug(5,
929 " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix "
930 "entry at row[%i] col[%i] value %g\n",
931 count, count + K, entry);
932 pos++;
933 if (les->type == N_SPARSE_LES) {
934 spvect->index[pos] = count + K;
935 spvect->values[pos] = entry;
936 }
937 else {
938 les->A[count][count + K] = entry;
939 }
940 }
941 }
942 }
943
944 return pos;
945}
946
947/* *************************************************************** *
948 * ******************** N_assemble_les_3d ************************ *
949 * *************************************************************** */
950/*!
951 * \brief Assemble a linear equation system (les) based on 3d location data
952 * (g3d) active cells
953 *
954 * This function calls #N_assemble_les_3d_param
955 * */
956N_les *N_assemble_les_3d(int les_type, N_geom_data *geom, N_array_3d *status,
957 N_array_3d *start_val, void *data,
958 N_les_callback_3d *call)
959{
960 return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
961 call, N_CELL_ACTIVE);
962}
963
964/*!
965 * \brief Assemble a linear equation system (les) based on 3d location data
966 * (g3d) active cells
967 *
968 * This function calls #N_assemble_les_3d_param
969 * */
971 N_array_3d *status, N_array_3d *start_val,
972 void *data, N_les_callback_3d *call)
973{
974 return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
975 call, N_CELL_ACTIVE);
976}
977
978/*!
979 * \brief Assemble a linear equation system (les) based on 3d location data
980 * (g3d) active and dirichlet cells
981 *
982 * This function calls #N_assemble_les_3d_param
983 * */
985 N_array_3d *status, N_array_3d *start_val,
986 void *data, N_les_callback_3d *call)
987{
988 return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
989 call, N_CELL_DIRICHLET);
990}
991
992/*!
993 * \brief Assemble a linear equation system (les) based on 3d location data
994 * (g3d)
995 *
996 * The linear equation system type can be set to N_NORMAL_LES to create a
997 * regular matrix, or to N_SPARSE_LES to create a sparse matrix. This function
998 * returns a new created linear equation system which can be solved with linear
999 * equation solvers. An 3d array with start values and an 3d status array must
1000 * be provided as well as the location geometry and a void pointer to data
1001 * passed to the callback which creates the les row entries. This callback
1002 * must be defined in the N_les_callback_3d structure.
1003 *
1004 * The creation of the les is parallelized with OpenMP.
1005 * If you implement new callbacks, please make sure that the
1006 * function calls are thread safe.
1007 *
1008 * the les can be created in two ways, with dirichlet and similar cells and
1009 * without them, to spare some memory. If the les is created with dirichlet
1010 * cell, the dirichlet boundary condition must be added.
1011 *
1012 * \param les_type int
1013 * \param geom N_geom_data*
1014 * \param status N_array_3d *
1015 * \param start_val N_array_3d *
1016 * \param data void *
1017 * \param call N_les_callback_3d *
1018 * \param cell_type int -- les assemble based on N_CELL_ACTIVE or
1019 * N_CELL_DIRICHLET \return N_les *
1020 * */
1022 N_array_3d *status, N_array_3d *start_val,
1023 void *data, N_les_callback_3d *call,
1024 int cell_type)
1025{
1026 int i, j, k, count = 0, pos = 0;
1027 int cell_type_count = 0;
1028 N_array_3d *cell_count;
1029 N_les *les = NULL;
1030 int **index_ij;
1031
1032 G_debug(
1033 2,
1034 "N_assemble_les_3d: starting to assemble the linear equation system");
1035
1036 cell_count =
1037 N_alloc_array_3d(geom->cols, geom->rows, geom->depths, 1, DCELL_TYPE);
1038
1039 /* First count the number of valid cells and save
1040 * each number in a new 3d array. Those numbers are used
1041 * to create the linear equation system.*/
1042
1043 if (cell_type == N_CELL_DIRICHLET) {
1044 /* include dirichlet cells in the les */
1045 for (k = 0; k < geom->depths; k++) {
1046 for (j = 0; j < geom->rows; j++) {
1047 for (i = 0; i < geom->cols; i++) {
1048 /*use all non-inactive cells for les creation */
1049 if (N_CELL_INACTIVE <
1050 (int)N_get_array_3d_d_value(status, i, j, k) &&
1051 (int)N_get_array_3d_d_value(status, i, j, k) <
1053 cell_type_count++;
1054 }
1055 }
1056 }
1057 }
1058 else {
1059 /*use only active cell in the les */
1060 for (k = 0; k < geom->depths; k++) {
1061 for (j = 0; j < geom->rows; j++) {
1062 for (i = 0; i < geom->cols; i++) {
1063 /*count only active cells */
1064 if (N_CELL_ACTIVE ==
1065 (int)N_get_array_3d_d_value(status, i, j, k))
1066 cell_type_count++;
1067 }
1068 }
1069 }
1070 }
1071
1072 G_debug(2, "N_assemble_les_3d: number of used cells %i\n",
1073 cell_type_count);
1074
1075 if (cell_type_count == 0.0)
1077 "Not enough active cells [%i] to create the linear equation "
1078 "system. Check the cell status. Only active cells (value = 1) are "
1079 "used to create the equation system.",
1080 cell_type_count);
1081
1082 /* allocate the memory for the linear equation system (les).
1083 * Only valid cells are used to create the les. */
1084 les = N_alloc_les_Ax_b(cell_type_count, les_type);
1085
1086 index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
1087 for (i = 0; i < cell_type_count; i++)
1088 index_ij[i] = (int *)G_calloc(3, sizeof(int));
1089
1090 count = 0;
1091 /*count the number of cells which should be used to create the linear
1092 * equation system */
1093 /*save the k, i and j indices and create a ordered numbering */
1094 for (k = 0; k < geom->depths; k++) {
1095 for (j = 0; j < geom->rows; j++) {
1096 for (i = 0; i < geom->cols; i++) {
1097 if (cell_type == N_CELL_DIRICHLET) {
1098 if (N_CELL_INACTIVE <
1099 (int)N_get_array_3d_d_value(status, i, j, k) &&
1100 (int)N_get_array_3d_d_value(status, i, j, k) <
1102 N_put_array_3d_d_value(cell_count, i, j, k, count);
1103 index_ij[count][0] = i;
1104 index_ij[count][1] = j;
1105 index_ij[count][2] = k;
1106 count++;
1107 G_debug(5,
1108 "N_assemble_les_3d: non-inactive cells count "
1109 "%i at pos x[%i] y[%i] z[%i]\n",
1110 count, i, j, k);
1111 }
1112 }
1113 else if (N_CELL_ACTIVE ==
1114 (int)N_get_array_3d_d_value(status, i, j, k)) {
1115 N_put_array_3d_d_value(cell_count, i, j, k, count);
1116 index_ij[count][0] = i;
1117 index_ij[count][1] = j;
1118 index_ij[count][2] = k;
1119 count++;
1120 G_debug(5,
1121 "N_assemble_les_3d: active cells count %i at pos "
1122 "x[%i] y[%i] z[%i]\n",
1123 count, i, j, k);
1124 }
1125 }
1126 }
1127 }
1128
1129 G_debug(2, "N_assemble_les_3d: starting the parallel assemble loop");
1130
1131#pragma omp parallel for private(i, j, k, pos, count) schedule(static)
1132 for (count = 0; count < cell_type_count; count++) {
1133 i = index_ij[count][0];
1134 j = index_ij[count][1];
1135 k = index_ij[count][2];
1136
1137 /*create the entries for the */
1138 N_data_star *items = call->callback(data, geom, i, j, k);
1139
1140 G_math_spvector *spvect = NULL;
1141
1142 /*allocate a sprase vector */
1143 if (les_type == N_SPARSE_LES)
1144 spvect = G_math_alloc_spvector(items->count);
1145 /* initial conditions */
1146
1147 les->x[count] = N_get_array_3d_d_value(start_val, i, j, k);
1148
1149 /* the entry in the vector b */
1150 les->b[count] = items->V;
1151
1152 /* pos describes the position in the sparse vector.
1153 * the first entry is always the diagonal entry of the matrix*/
1154 pos = 0;
1155
1156 if (les_type == N_SPARSE_LES) {
1157 spvect->index[pos] = count;
1158 spvect->values[pos] = items->C;
1159 }
1160 else {
1161 les->A[count][count] = items->C;
1162 }
1163 /* western neighbour, entry is col - 1 */
1164 if (i > 0) {
1165 pos = make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect,
1166 cell_count, status, start_val, items->W,
1167 cell_type);
1168 }
1169 /* eastern neighbour, entry col + 1 */
1170 if (i < geom->cols - 1) {
1171 pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect,
1172 cell_count, status, start_val, items->E,
1173 cell_type);
1174 }
1175 /* northern neighbour, entry row -1 */
1176 if (j > 0) {
1177 pos = make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect,
1178 cell_count, status, start_val, items->N,
1179 cell_type);
1180 }
1181 /* southern neighbour, entry row +1 */
1182 if (j < geom->rows - 1) {
1183 pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect,
1184 cell_count, status, start_val, items->S,
1185 cell_type);
1186 }
1187 /*only for a 7 star entry needed */
1188 if (items->type == N_7_POINT_STAR || items->type == N_27_POINT_STAR) {
1189 /* the upper cell (top), entry depth + 1 */
1190 if (k < geom->depths - 1) {
1191 pos = make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les,
1192 spvect, cell_count, status, start_val,
1193 items->T, cell_type);
1194 }
1195 /* the lower cell (bottom), entry depth - 1 */
1196 if (k > 0) {
1197 pos = make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les,
1198 spvect, cell_count, status, start_val,
1199 items->B, cell_type);
1200 }
1201 }
1202
1203 /*How many entries in the les */
1204 if (les->type == N_SPARSE_LES) {
1205 spvect->cols = pos + 1;
1206 G_math_add_spvector(les->Asp, spvect, count);
1207 }
1208
1209 if (items)
1210 G_free(items);
1211 }
1212
1213 N_free_array_3d(cell_count);
1214
1215 for (i = 0; i < cell_type_count; i++)
1216 G_free(index_ij[i]);
1217
1218 G_free(index_ij);
1219
1220 return les;
1221}
1222
1223/*!
1224 * \brief Integrate Dirichlet or Transmission boundary conditions into the les
1225 * (3d)
1226 *
1227 * Dirichlet and Transmission boundary conditions will be integrated into
1228 * the provided linear equation system. This is meaningful if
1229 * the les was created with #N_assemble_les_2d_dirichlet, because in
1230 * this case Dirichlet boundary conditions are not automatically included.
1231 *
1232 * The provided les will be modified:
1233 *
1234 * Ax = b will be split into Ax_u + Ax_d = b
1235 *
1236 * x_u - the unknowns
1237 * x_d - the Dirichlet cells
1238 *
1239 * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
1240 *
1241 * | A_u 0 | x_u
1242 * | 0 I | x_d
1243 *
1244 * \param les N_les* -- the linear equation system
1245 * \param geom N_geom_data* -- geometrical data information
1246 * \param status N_array_2d* -- the status array containing the cell types
1247 * \param start_val N_array_2d* -- an array with start values
1248 * \return int -- 1 = success, 0 = failure
1249 * */
1251 N_array_3d *status, N_array_3d *start_val)
1252{
1253 int rows, cols, depths;
1254 int count = 0;
1255 int i, j, x, y, z, stat;
1256 double *dvect1;
1257 double *dvect2;
1258
1259 G_debug(2, "N_les_integrate_dirichlet_3d: integrating the dirichlet "
1260 "boundary condition");
1261
1262 rows = geom->rows;
1263 cols = geom->cols;
1264 depths = geom->depths;
1265
1266 /*we need to additional vectors */
1267 dvect1 = (double *)G_calloc(les->cols, sizeof(double));
1268 dvect2 = (double *)G_calloc(les->cols, sizeof(double));
1269
1270 /*fill the first one with the x vector data of Dirichlet cells */
1271 count = 0;
1272 for (z = 0; z < depths; z++) {
1273 for (y = 0; y < rows; y++) {
1274 for (x = 0; x < cols; x++) {
1275 stat = (int)N_get_array_3d_d_value(status, x, y, z);
1276 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
1277 dvect1[count] = N_get_array_3d_d_value(start_val, x, y, z);
1278 count++;
1279 }
1280 else if (stat == N_CELL_ACTIVE) {
1281 dvect1[count] = 0.0;
1282 count++;
1283 }
1284 }
1285 }
1286 }
1287
1288#pragma omp parallel default(shared)
1289 {
1290 /*perform the matrix vector product and */
1291 if (les->type == N_SPARSE_LES)
1292 G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows);
1293 else
1294 G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols);
1295#pragma omp for schedule(static) private(i)
1296 for (i = 0; i < les->cols; i++)
1297 les->b[i] = les->b[i] - dvect2[i];
1298 }
1299
1300 /*now set the Dirichlet cell rows and cols to zero and the
1301 * diagonal entry to 1*/
1302 count = 0;
1303 for (z = 0; z < depths; z++) {
1304 for (y = 0; y < rows; y++) {
1305 for (x = 0; x < cols; x++) {
1306 stat = (int)N_get_array_3d_d_value(status, x, y, z);
1307 if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
1308 if (les->type == N_SPARSE_LES) {
1309 /*set the rows to zero */
1310 for (i = 0; (unsigned int)i < les->Asp[count]->cols;
1311 i++)
1312 les->Asp[count]->values[i] = 0.0;
1313 /*set the cols to zero */
1314 for (i = 0; i < les->rows; i++) {
1315 for (j = 0; (unsigned int)j < les->Asp[i]->cols;
1316 j++) {
1317 if (les->Asp[i]->index[j] ==
1318 (unsigned int)count)
1319 les->Asp[i]->values[j] = 0.0;
1320 }
1321 }
1322
1323 /*entry on the diagonal */
1324 les->Asp[count]->values[0] = 1.0;
1325 }
1326 else {
1327 /*set the rows to zero */
1328 for (i = 0; i < les->cols; i++)
1329 les->A[count][i] = 0.0;
1330 /*set the cols to zero */
1331 for (i = 0; i < les->rows; i++)
1332 les->A[i][count] = 0.0;
1333
1334 /*entry on the diagonal */
1335 les->A[count][count] = 1.0;
1336 }
1337 }
1338 count++;
1339 }
1340 }
1341 }
1342
1343 return 0;
1344}
1345
1346/* **************************************************************** */
1347/* **** make an entry in the les (3d) ***************************** */
1348/* **************************************************************** */
1349int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
1350 int offset_k, int count, int pos, N_les *les,
1351 G_math_spvector *spvect, N_array_3d *cell_count,
1352 N_array_3d *status, N_array_3d *start_val, double entry,
1353 int cell_type)
1354{
1355 int K;
1356 int di = offset_i;
1357 int dj = offset_j;
1358 int dk = offset_k;
1359
1360 K = (int)N_get_array_3d_d_value(cell_count, i + di, j + dj, k + dk) -
1361 (int)N_get_array_3d_d_value(cell_count, i, j, k);
1362
1363 if (cell_type == N_CELL_ACTIVE) {
1364 if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) >
1365 N_CELL_ACTIVE &&
1366 (int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) <
1368 les->b[count] -=
1369 N_get_array_3d_d_value(start_val, i + di, j + dj, k + dk) *
1370 entry;
1371 else if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) ==
1372 N_CELL_ACTIVE) {
1373 if ((count + K) >= 0 && (count + K) < les->cols) {
1374 G_debug(5,
1375 " make_les_entry_3d: (N_CELL_ACTIVE) create matrix "
1376 "entry at row[%i] col[%i] value %g\n",
1377 count, count + K, entry);
1378 pos++;
1379 if (les->type == N_SPARSE_LES) {
1380 spvect->index[pos] = count + K;
1381 spvect->values[pos] = entry;
1382 }
1383 else {
1384 les->A[count][count + K] = entry;
1385 }
1386 }
1387 }
1388 }
1389 else if (cell_type == N_CELL_DIRICHLET) {
1390 if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) !=
1392 if ((count + K) >= 0 && (count + K) < les->cols) {
1393 G_debug(5,
1394 " make_les_entry_3d: (N_CELL_DIRICHLET) create matrix "
1395 "entry at row[%i] col[%i] value %g\n",
1396 count, count + K, entry);
1397 pos++;
1398 if (les->type == N_SPARSE_LES) {
1399 spvect->index[pos] = count + K;
1400 spvect->values[pos] = entry;
1401 }
1402 else {
1403 les->A[count][count + K] = entry;
1404 }
1405 }
1406 }
1407 }
1408
1409 return pos;
1410}
#define N_27_POINT_STAR
Definition N_pde.h:43
#define N_5_POINT_STAR
Definition N_pde.h:40
#define N_9_POINT_STAR
Definition N_pde.h:42
#define N_CELL_DIRICHLET
Definition N_pde.h:32
#define N_CELL_INACTIVE
Definition N_pde.h:30
#define N_CELL_ACTIVE
Definition N_pde.h:31
#define N_MAX_CELL_STATE
the maximum number of available cell states (eg: boundary condition, inactiven active)
Definition N_pde.h:38
#define N_7_POINT_STAR
Definition N_pde.h:41
#define N_SPARSE_LES
Definition N_pde.h:26
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
Compute the matrix - vector product of matrix A and vector x.
#define NULL
Definition ccmath.h:32
#define SE
Definition dataquad.h:31
#define SW
Definition dataquad.h:30
#define NE
Definition dataquad.h:29
#define NW
Definition dataquad.h:28
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
int count
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.
Definition n_arrays.c:719
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.
Definition n_arrays.c:314
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition n_arrays.c:774
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.
Definition n_arrays.c:380
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition n_arrays.c:132
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.
Definition n_arrays.c:1148
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition n_arrays.c:75
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.
Definition n_arrays.c:979
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.
Definition n_arrays.c:516
N_les * N_alloc_les_Ax_b(int rows, int type)
Allocate memory for a quadratic linear equation system which includes the Matrix A,...
Definition n_les.c:149
N_data_star * N_alloc_27star(void)
allocate a 27 point star data structure
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...
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_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
N_data_star * N_alloc_7star(void)
allocate a 7 point star data structure
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.
N_data_star * N_alloc_9star(void)
allocate a 9 point star data structure
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_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.
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)
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
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.
N_data_star * N_callback_template_3d(void *data UNUSED, N_geom_data *geom, int col UNUSED, int row UNUSED, int depth UNUSED)
A callback template creates a 7 point star structure.
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_5star(void)
allocate a 5 point star data structure
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_data_star * N_callback_template_2d(void *data UNUSED, N_geom_data *geom, int col UNUSED, int row UNUSED)
A callback template creates a 9 point star structure.
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_callback_2d * N_alloc_les_callback_2d(void)
Allocate the structure holding the callback function.
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.
N_les_callback_3d * N_alloc_les_callback_3d(void)
Allocate the structure holding the callback function.
void N_set_les_callback_3d_func(N_les_callback_3d *data, N_data_star *(*callback_func_3d)(void *, N_geom_data *, int, int, int))
Set the callback function which is called while assembling the les in 3d.
void N_set_les_callback_2d_func(N_les_callback_2d *data, N_data_star *(*callback_func_2d)(void *, N_geom_data *, int, int))
Set the callback function which is called while assembling the les in 2d.
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_Ax_sparse(G_math_spvector **Asp, double *x, double *y, int rows)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
G_math_spvector * G_math_alloc_spvector(int cols)
Allocate memory for a sparse vector.
int G_math_add_spvector(G_math_spvector **Asp, G_math_spvector *spvector, int row)
Adds a sparse vector to a sparse matrix at position row.
Matrix entries for a mass balance 5/7/9 star system.
Definition N_pde.h:295
double W_B
Definition N_pde.h:302
double E
Definition N_pde.h:298
double S
Definition N_pde.h:298
double NE
Definition N_pde.h:298
double E_B
Definition N_pde.h:302
double N_T
Definition N_pde.h:300
double SE
Definition N_pde.h:298
double N_B
Definition N_pde.h:302
double NW
Definition N_pde.h:298
double W
Definition N_pde.h:298
double S_T
Definition N_pde.h:300
double NW_T
Definition N_pde.h:300
double E_T
Definition N_pde.h:300
double NW_B
Definition N_pde.h:302
double S_B
Definition N_pde.h:302
double W_T
Definition N_pde.h:300
double SW
Definition N_pde.h:298
double SW_T
Definition N_pde.h:300
double C
Definition N_pde.h:298
int type
Definition N_pde.h:296
double T
Definition N_pde.h:300
double SW_B
Definition N_pde.h:302
double NE_T
Definition N_pde.h:300
double N
Definition N_pde.h:298
double B
Definition N_pde.h:302
double SE_T
Definition N_pde.h:300
int count
Definition N_pde.h:297
double SE_B
Definition N_pde.h:302
double NE_B
Definition N_pde.h:302
double V
Definition N_pde.h:298
Geometric information about the structured grid.
Definition N_pde.h:101
double dx
Definition N_pde.h:108
int depths
Definition N_pde.h:114
double dy
Definition N_pde.h:109
double dz
Definition N_pde.h:110
int rows
Definition N_pde.h:115
int cols
Definition N_pde.h:116
callback structure for 2d matrix assembling
Definition N_pde.h:315
N_data_star *(* callback)(void *, N_geom_data *, int, int)
Definition N_pde.h:316
callback structure for 3d matrix assembling
Definition N_pde.h:308
N_data_star *(* callback)(void *, N_geom_data *, int, int, int)
Definition N_pde.h:309
The linear equation system (les) structure.
Definition N_pde.h:71
int cols
Definition N_pde.h:77
int rows
Definition N_pde.h:76
double * x
Definition N_pde.h:72
double ** A
Definition N_pde.h:74
double * b
Definition N_pde.h:73
int type
Definition N_pde.h:79
G_math_spvector ** Asp
Definition N_pde.h:75
#define x