1 /*******************************************************************************
2 * ALBERTA: an Adaptive multi Level finite element toolbox using
3 * Bisectioning refinement and Error control by Residual
4 * Techniques for scientific Applications
5 *
6 * file: alberta.h
7 *
8 * description: public header file of the ALBERTA package
9 *
10 *******************************************************************************
11 *
12 * authors: Alfred Schmidt
13 * Zentrum fuer Technomathematik
14 * Fachbereich 3 Mathematik/Informatik
15 * Universitaet Bremen
16 * Bibliothekstr. 2
17 * D-28359 Bremen, Germany
18 *
19 * Kunibert G. Siebert
20 * Institut fuer Mathematik
21 * Universitaet Augsburg
22 * Universitaetsstr. 14
23 * D-86159 Augsburg, Germany
24 *
25 * Daniel Koester
26 * Institut fuer Mathematik
27 * Universitaet Augsburg
28 * Universitaetsstr. 14
29 * D-86159 Augsburg, Germany
30 *
31 * Claus-Justus Heine
32 * Abteilung fuer Angewandte Mathematik
33 * Universitaet Freiburg
34 * Hermann-Herder-Strasse 10
35 * 79104 Freiburg, Germany
36 *
37 * http://www.alberta-fem.de
38 *
39 * (c) by A. Schmidt and K.G. Siebert (1996-2005),
40 * D. Koester (2002-2005),
41 * C.-J. Heine (2002-2009).
42 *
43 ******************************************************************************/
44
45 /*******************************************************************************
46 * Header-File for ALBERTA utilities
47 ******************************************************************************/
48 #ifndef _ALBERTA_H_
49 #define _ALBERTA_H_
50
51 #include <alberta/alberta_util.h>
52
53 #ifdef __cplusplus
54 extern "C" {
55 #elif 0
56 } /* some editors try to indent because of the brace above */
57 #endif
58
59 /* The version string used for disk-IO. Note the space after the
60 * version number, the version number is expected to be of the forn
61 * "X.Y " or "X.YZ". The IO version must not necessarily relate to the
62 * version of the library.
63 */
64
65 #define ALBERTA_MAGIC "ALBERTA: Version "
66 #define ALBERTA_VERSION ALBERTA_MAGIC"2.3 "
67
68 /*******************************************************************************
69 * Definition of the space dimension and of parameters depending on the
70 * space dimension:
71 *
72 * DIM_OF_WORLD: space dimension
73 *
74 * The ?D-suffix signals different simplex dimensions (formerly ==DIM).
75 ******************************************************************************/
76
77 #ifndef DIM_OF_WORLD
78 # error DIM_OF_WORLD UNDEFINED
79 #endif
80
81 #ifndef ALBERTA_DEBUG
82 /* #warning is a GNU extension, and there is no way to get at compiler
83 * switches like -Werror form here, so do NOT use it.
84 */
85 /* # warning ALBERTA_DEBUG WAS NOT DEFINED! DEFAULTING TO 0. */
86 # define ALBERTA_DEBUG 0
87 #endif
88
89 /* meshes of at most this dimension are supported. This macro can
90 * be use to lay-out static arrays.
91 */
92 #define DIM_LIMIT 3
93
94 /* The master dimension limit, meshes of higher dimension are not
95 * supported.
96 */
97 #define DIM_MAX MIN(DIM_OF_WORLD, DIM_LIMIT)
98
99 /* Various constants for dimension dependent geometrical quantities
100 *
101 * A note about the terminology: caused by historic reasons the
102 * terminology for simplexes is of naive nature: 1d and 2d simplexes
103 * do not have faces (unluckily). Instead, 1d simplexes have vertices
104 * and 2d simplexes have edges. Of course, 3d simplexes also have
105 * vertices and edges, but: those are _very_ different from 1d and 2d
106 * vertices and edges.
107 *
108 * The best solution would be to reserve the name "face" for a
109 * co-dimension one face-simplex. Unluckily, there does not seem to be
110 * consense about this. Therefore, the naive ALBERTA notion of "face"
111 * is maintained, but supplemented by the notion of a "wall", which
112 * denotes a co-dimension one (1) face-simplex. Whenever the notation
113 * "wall" is used somewhere in ALBERTA then it denotes stuff
114 * concerning co-dimension 1 simplexes (i.e. "faces"). Like WALL_QUAD,
115 * EL_INFO->wall_bound, wall-transformation (for periodic boundaries)
116 * etc.
117 *
118 * In my opionion this is a FIXME. "face" should be used everywhere to
119 * denote a co-dimension 1 face-simplex. In German "Face" means "Seite".
120 *
121 * cH.
122 */
123 #define N_VERTICES(DIM) ((DIM)+1)
124 #define N_EDGES(DIM) ((DIM)*((DIM)+1)/2)
125 #define N_WALLS(DIM) ((DIM)+1) /* number of codim 1 subsimplexes */
126 #define N_FACES(DIM) (((DIM) == 3) * N_WALLS(DIM))
127 #define N_NEIGH(DIM) (((DIM) != 0) * N_WALLS(DIM))
128 #define N_LAMBDA(DIM) N_VERTICES(DIM)
129 #define DIM_FAC(DIM) ((DIM) < 2 ? 1 : (DIM) == 2 ? 2 : 6)
130
131 #define VERTEX_OF_EDGE(DIM, EDGE) \
132 ((DIM) == 1 \
133 ? vertex_of_edge_1d[(EDGE)] \
134 : ((DIM == 2) \
135 ? vertex_of_edge_2d[(EDGE)] \
136 : vertex_of_edge_3d[(EDGE)]))
137
138 #define VERTEX_OF_WALL(DIM, WALL) \
139 ((DIM) == 1 \
140 ? vertex_of_wall_1d[(WALL)] \
141 : ((DIM == 2) \
142 ? vertex_of_wall_2d[(WALL)] \
143 : vertex_of_wall_3d[(WALL)]))
144
145 #define N_VERTICES_0D N_VERTICES(0)
146 #define N_EDGES_0D N_EDGES(0)
147 #define N_FACES_0D N_FACES(0)
148 #define N_NEIGH_0D N_NEIGH(0)
149 #define N_WALLS_0D N_WALLS(0)
150 #define N_LAMBDA_0D N_LAMBDA(0)
151 #define DIM_FAC_0D DIM_FAC(0)
152 #define VERTEX_OF_EDGE_0D NULL
153 #define VERTEX_OF_WALL_0D NULL
154
155 #define N_VERTICES_1D N_VERTICES(1)
156 #define N_EDGES_1D N_EDGES(1)
157 #define N_FACES_1D N_FACES(1)
158 #define N_NEIGH_1D N_NEIGH(1)
159 #define N_WALLS_1D N_WALLS(1)
160 #define N_LAMBDA_1D N_LAMBDA(1)
161 #define DIM_FAC_1D DIM_FAC(1)
162 #define VERTEX_OF_EDGE_1D(E) VERTEX_OF_EDGE(1, E)
163 #define VERTEX_OF_WALL_1D(W) VERTEX_OF_WALL(1, W)
164
165 #define N_VERTICES_2D N_VERTICES(2)
166 #define N_EDGES_2D N_EDGES(2)
167 #define N_FACES_2D N_FACES(2)
168 #define N_NEIGH_2D N_NEIGH(2)
169 #define N_WALLS_2D N_WALLS(2)
170 #define N_LAMBDA_2D N_LAMBDA(2)
171 #define DIM_FAC_2D DIM_FAC(2)
172 #define VERTEX_OF_EDGE_2D(E) VERTEX_OF_EDGE(2, E)
173 #define VERTEX_OF_WALL_2D(W) VERTEX_OF_WALL(2, W)
174
175 #define N_VERTICES_3D N_VERTICES(3)
176 #define N_EDGES_3D N_EDGES(3)
177 #define N_FACES_3D N_FACES(3)
178 #define N_NEIGH_3D N_NEIGH(3)
179 #define N_WALLS_3D N_WALLS(3)
180 #define N_LAMBDA_3D N_LAMBDA(3)
181 #define DIM_FAC_3D DIM_FAC(3)
182 #define VERTEX_OF_EDGE_3D(E) VETEX_OF_EDGE(3, E)
183 #define VERTEX_OF_WALL_3D(W) VERTEX_OF_WALL(3, W)
184
185 /* The maximal number for a given DIM_OF_WORLD.
186 */
187 #define N_LAMBDA_MAX N_VERTICES(DIM_MAX)
188 #define N_VERTICES_MAX N_VERTICES(DIM_MAX)
189 #define N_EDGES_MAX N_EDGES(DIM_MAX)
190 #define N_FACES_MAX N_FACES(DIM_MAX)
191 #define N_NEIGH_MAX N_NEIGH(DIM_MAX)
192 #define N_WALLS_MAX N_WALLS(DIM_MAX)
193 #define DIM_FAC_MAX DIM_FAC(DIM_MAX)
194
195 /* The maximal number which can possibly occur within ALBERTA. These
196 * macros can be used to layout static arrays. Of course, DIM_LIMIT is
197 * hard-wired to 3, and that will not change.
198 */
199 #define N_LAMBDA_LIMIT N_VERTICES(DIM_LIMIT)
200 #define N_VERTICES_LIMIT N_VERTICES(DIM_LIMIT)
201 #define N_EDGES_LIMIT N_EDGES(DIM_LIMIT)
202 #define N_FACES_LIMIT N_FACES(DIM_LIMIT)
203 #define N_NEIGH_LIMIT N_NEIGH(DIM_LIMIT)
204 #define N_WALLS_LIMIT N_WALLS(DIM_LIMIT)
205 #define DIM_FAC_LIMIT DIM_FAC(DIM_LIMIT)
206
207 /* As N_LAMBDA depends on DIM_MAX, we provide some convenience
208 * macros for initializing arrays etc.
209 */
210 #if DIM_MAX == 0
211 # define INIT_BARY_0D(a) { 1.0 }
212 # define INIT_BARY_1D(a, b) { 1.0 }
213 # define INIT_BARY_2D(a, b, c) { 1.0 }
214 # define INIT_BARY_3D(a, b, c, d) { 1.0 }
215 # define INIT_BARY_MAX(a, b, c, d) INIT_BARY_0D(a)
216 #elif DIM_MAX == 1
217 # define INIT_BARY_0D(a) { (a), 0.0 }
218 # define INIT_BARY_1D(a, b) { (a), (b) }
219 # define INIT_BARY_2D(a, b, c) { (a), (b) }
220 # define INIT_BARY_3D(a, b, c, d) { (a), (b) }
221 # define INIT_BARY_MAX(a, b, c, d) INIT_BARY_1D(a, b)
222 #elif DIM_MAX == 2
223 # define INIT_BARY_0D(a) { (a), 0.0, 0.0 }
224 # define INIT_BARY_1D(a, b) { (a), (b), 0.0 }
225 # define INIT_BARY_2D(a, b, c) { (a), (b), (c) }
226 # define INIT_BARY_3D(a, b, c, d) { (a), (b), (c) }
227 # define INIT_BARY_MAX(a, b, c, d) INIT_BARY_2D(a, b, c)
228 #elif DIM_MAX == 3
229 # define INIT_BARY_0D(a) { (a), 0.0, 0.0, 0.0 }
230 # define INIT_BARY_1D(a, b) { (a), (b), 0.0, 0.0 }
231 # define INIT_BARY_2D(a, b, c) { (a), (b), (c), 0.0 }
232 # define INIT_BARY_3D(a, b, c, d) { (a), (b), (c), (d) }
233 # define INIT_BARY_MAX(a, b, c, d) INIT_BARY_3D(a, b, c, d)
234 #else
235 # error Unsupported DIM_MAX
236 #endif
237
238 /* Various matrix and vector types. Please use them as apropriate.
239 *
240 * A famous fortune reading (grin):
241 *
242 * The primary purpose of the DATA statement is to give names to
243 * constants; instead of referring to pi as 3.141592653589793 at every
244 * appearance, the variable PI can be given that value with a DATA
245 * statement and used instead of the longer form of the constant.
246 * This also simplifies modifying the program, should the value of pi
247 * change.
248 *
249 * -- FORTRAN manual for Xerox Computers
250 *
251 */
252 typedef REAL REAL_B[N_LAMBDA_MAX];
253 typedef REAL_B REAL_BB[N_LAMBDA_MAX];
254 typedef REAL REAL_D[DIM_OF_WORLD];
255 typedef REAL_D REAL_DD[DIM_OF_WORLD];
256 typedef REAL_D REAL_BD[N_LAMBDA_MAX];
257 typedef REAL_BD REAL_BBD[N_LAMBDA_MAX];
258 typedef REAL_DD REAL_DDD[DIM_OF_WORLD];
259 typedef REAL_DD REAL_BDD[N_LAMBDA_MAX];
260 typedef REAL_BDD REAL_BBDD[N_LAMBDA_MAX];
261 typedef REAL_B REAL_DB[DIM_OF_WORLD];
262 typedef REAL_BB REAL_DBB[DIM_OF_WORLD];
263 typedef REAL_BB REAL_BBB[N_LAMBDA_MAX];
264 typedef REAL_BBB REAL_BBBB[N_LAMBDA_MAX];
265 typedef REAL_BBB REAL_DBBB[DIM_OF_WORLD];
266 typedef REAL_BBBB REAL_DBBBB[DIM_OF_WORLD];
267 typedef REAL_DB REAL_BDB[N_LAMBDA_MAX];
268 typedef REAL_DBB REAL_BDBB[N_LAMBDA_MAX];
269
270 /******************************************************************************/
271
272 /* The maximum number of quadrature points and local basis functions
273 * for each dimension. Useful to define C99 variable size arrays.
274 * Note that those fields are intentionlly are NOT labeled "const", it
275 * is possible to generate quadrature formulas of arbitrary degree
276 * using "get_product_quad(). The INIT_ELEMENT() frame-work provides
277 * means to introduce quadratures with per-element initialization,
278 * etc. An application may install new quadrature rules using
279 * "register_quadrature()" or "add_quadrature()". Similar thing hold
280 * for basis functions.
281 */
282 extern int n_quad_points_max[];
283 extern int n_bas_fcts_max[];
284
285 /*******************************************************************************
286 * access to element index via element or element_info structure
287 ******************************************************************************/
288
289 #if ALBERTA_DEBUG
290 #define INDEX(el) ((el) ? (el)->index : -1)
291 #else
292 #define INDEX(el) -1
293 #endif
294
295 /*******************************************************************************
296 * access to leaf data (only for leaf elements)
297 ******************************************************************************/
298
299 #define IS_LEAF_EL(el) (!(el)->child[0])
300 #define LEAF_DATA(el) ((void *)(el)->child[1])
301
302 /*******************************************************************************
303 * boundary types
304 ******************************************************************************/
305
306 #define INTERIOR 0
307 #define DIRICHLET 1
308 #define NEUMANN -1
309
310 #define IS_NEUMANN(bound) ((bound) <= NEUMANN)
311 #define IS_DIRICHLET(bound) ((bound) >= DIRICHLET)
312 #define IS_INTERIOR(bound) ((bound) == 0)
313
314 /*******************************************************************************
315 * node types (indices in n_dof[] vectors, e.g.)
316 ******************************************************************************/
317
318 /* Be careful: in 1D we have only VERTEX and CENTER nodes (although
319 * that violates the usual geometric meaning of VERTEX/EDGE/FACE:
320 * looking at the 2d/3d code one really would expect 1d CENTER DOFs to
321 * be 1d EDGE DOFs, but this is not the case).
322 *
323 * So:
324 *
325 * 1d: VERTEX and CENTER, EL_INFO->wall_bound refers to VERTEX boundary type
326 * 2d: VERTEX and CENTER and EDGE, EL_INFO->wall_bound refers to EDGEs
327 * 3d: FACE comes into play, EL_INFO->wall_bound refers to FACEs
328 */
329 enum node_types {
330 VERTEX = 0,
331 CENTER,
332 EDGE,
333 FACE,
334 N_NODE_TYPES
335 };
336
337 /*******************************************************************************
338 * basic types of the grid
339 ******************************************************************************/
340
341 typedef signed int DOF;
342 typedef enum node_types NODE_TYPES;
343 #define N_BNDRY_TYPES 256
344 typedef U_CHAR BNDRY_TYPE;
345 typedef BITS_256 BNDRY_FLAGS;
346 typedef struct el EL;
347 typedef struct macro_el MACRO_EL;
348 typedef struct el_info EL_INFO;
349 typedef struct el_geom_cache EL_GEOM_CACHE;
350 typedef struct rc_list_el RC_LIST_EL;
351 typedef struct mesh MESH;
352
353 typedef struct parametric PARAMETRIC;
354 typedef struct traverse_stack TRAVERSE_STACK;
355
356 typedef struct adapt_stat ADAPT_STAT;
357 typedef struct adapt_instat ADAPT_INSTAT;
358
359 #ifndef DOF_ADMIN_DEF
360 typedef struct dof_admin DOF_ADMIN;
361 typedef struct dof_int_vec DOF_INT_VEC;
362 typedef struct dof_dof_vec DOF_DOF_VEC;
363 typedef struct dof_uchar_vec DOF_UCHAR_VEC;
364 typedef struct dof_schar_vec DOF_SCHAR_VEC;
365 typedef struct dof_real_vec DOF_REAL_VEC;
366 typedef struct dof_real_d_vec DOF_REAL_D_VEC;
367 typedef struct dof_real_dd_vec DOF_REAL_DD_VEC;
368 typedef struct dof_ptr_vec DOF_PTR_VEC;
369 typedef struct dof_real_vec_d DOF_REAL_VEC_D;
370 typedef struct dof_matrix DOF_MATRIX;
371 typedef struct matrix_row MATRIX_ROW;
372 typedef struct matrix_row_real MATRIX_ROW_REAL;
373 typedef struct matrix_row_real_d MATRIX_ROW_REAL_D;
374 typedef struct matrix_row_real_dd MATRIX_ROW_REAL_DD;
375 #endif
376
377 typedef struct el_matrix EL_MATRIX;
378 typedef struct el_int_vec EL_INT_VEC;
379 typedef struct el_dof_vec EL_DOF_VEC;
380 typedef struct el_uchar_vec EL_UCHAR_VEC;
381 typedef struct el_schar_vec EL_SCHAR_VEC;
382 typedef struct el_bndry_vec EL_BNDRY_VEC;
383 typedef struct el_ptr_vec EL_PTR_VEC;
384 typedef struct el_real_vec EL_REAL_VEC;
385 typedef struct el_real_dd_vec EL_REAL_DD_VEC;
386 typedef struct el_real_d_vec EL_REAL_D_VEC;
387 typedef struct el_real_vec_d EL_REAL_VEC_D;
388
389 typedef struct bas_fcts BAS_FCTS;
390 typedef struct fe_space FE_SPACE;
391
392 typedef struct quadrature QUAD;
393 typedef struct quadrature QUADRATURE;
394 typedef struct quad_fast QUAD_FAST;
395
396 typedef struct quad_el_cache QUAD_EL_CACHE;
397
398 typedef struct wall_quadrature WALL_QUAD;
399 typedef struct wall_quad_fast WALL_QUAD_FAST;
400
401 typedef struct macro_data MACRO_DATA;
402 typedef struct node_projection NODE_PROJ;
403 typedef struct node_projection NODE_PROJECTION;
404 typedef struct aff_trafo AFF_TRAFO;
405
406 typedef struct dof_comp_hook DOF_COMP_HOOK;
407
408 typedef REAL (*LOC_FCT_AT_QP)(const EL_INFO *el_info,
409 const QUAD *quad, int iq,
410 void *ud);
411 typedef const REAL *(*LOC_FCT_D_AT_QP)(REAL_D result,
412 const EL_INFO *el_info,
413 const QUAD *quad, int iq,
414 void *ud);
415 typedef const REAL *(*GRD_LOC_FCT_AT_QP)(REAL_D res,
416 const EL_INFO *el_info,
417 const REAL_BD Lambda,
418 const QUAD *quad, int iq,
419 void *ud);
420 typedef const REAL_D *(*GRD_LOC_FCT_D_AT_QP)(REAL_DD res,
421 const EL_INFO *el_info,
422 const REAL_BD Lambda,
423 const QUAD *quad, int iq,
424 void *ud);
425
426 typedef REAL (*FCT_AT_X)(const REAL_D x);
427 typedef const REAL *(*GRD_FCT_AT_X)(const REAL_D x, REAL_D result);
428 typedef const REAL_D *(*D2_FCT_AT_X)(const REAL_D x, REAL_DD result);
429
430 typedef const REAL *(*FCT_D_AT_X)(const REAL_D x, REAL_D result);
431 typedef const REAL_D *(*GRD_FCT_D_AT_X)(const REAL_D x, REAL_DD result);
432 typedef const REAL_DD *(*D2_FCT_D_AT_X)(const REAL_D x, REAL_DDD result);
433
434
435 /*******************************************************************************
436 * Macro for calling STRUCT->init_element() routines, per element
437 * initializers for (at least) QUAD, BAS_FCTS, QUAD_FAST, WALL_QUAD,
438 * WALL_QUAD_FAST
439 ******************************************************************************/
440
441 /* INIT_ELEMENT(el_info, object)
442 *
443 * The following convention holds:
444 *
445 * a) This macro evaluates to INIT_EL_TAG_DFLT when no
446 * element-initializer is present.
447 *
448 * b) An init_element() method MUST allow a NULL pointer for the
449 * el_info argument. If called with el_info == NULL the
450 * init_element() method must restore its default behaviour. The
451 * "default case" is what the implementation defines as default;
452 * for performance reasons the default case should be the one which
453 * applies to the majority of mesh elements.
454 *
455 * c) The return value of the init_element() method must be
456 * INIT_EL_TAG_DFLT for the default case.
457 *
458 * d) The return value of the init_element() method must be
459 * INIT_EL_TAG_NULL for the NULL case, meaning, e.g., the number of
460 * basis functions is 0, or the number of quadrature points is
461 * zero. The application can assume that in the NULL case the
462 * structure does not contain any real data.
463 *
464 * e) In all other cases the return value is a tag which is used to
465 * efficiently cache values of intermediate computations, e.g. the
466 * values of basis functions at quadrature points. This tag should
467 * be locally unique, meaning that consecutive invocations of
468 * init_element() should return different tags for different
469 * simplexes. This can be used for optimizations: if the tag
470 * returned by an init_element() routine does not change, then the
471 * calling function may assume that the underlying object has not
472 * changed.
473 */
474
475 enum {
476 INIT_EL_TAG_NONE = 0, /* invalid tag */
477 INIT_EL_TAG_DFLT = 1, /* default case */
478 INIT_EL_TAG_NULL = 2 /* something is 0, e.g. no quad-points, basis
479 * functions are identically zero and so on.
480 */
481 };
482
483 typedef unsigned int INIT_EL_TAG;
484 typedef INIT_EL_TAG (*INIT_ELEMENT_FCT)(const EL_INFO *el_info, void *thisptr);
485
486 /* Tag context. */
487 typedef struct init_el_tag_ctx {
488 INIT_EL_TAG tag;
489 unsigned int cnt;
490 } INIT_EL_TAG_CTX;
491
492 #define INIT_EL_TAG_CTX_INIT(ctx) \
493 { \
494 (ctx)->tag = INIT_EL_TAG_DFLT; \
495 (ctx)->cnt = 0; \
496 }
497
498 /* Generate a new unique tag != NULL & DFLT */
499 #define INIT_EL_TAG_CTX_UNIQ(ctx) \
500 { \
501 (ctx)->tag = INIT_EL_TAG_NULL + (++((ctx)->cnt)); \
502 if ((ctx)->tag == INIT_EL_TAG_NONE) { \
503 (ctx)->cnt = 1; \
504 (ctx)->tag = INIT_EL_TAG_NULL + 1; \
505 } \
506 }
507 #define INIT_EL_TAG_CTX_NULL(ctx) (ctx)->tag = INIT_EL_TAG_NULL
508 #define INIT_EL_TAG_CTX_DFLT(ctx) (ctx)->tag = INIT_EL_TAG_DFLT
509 #define INIT_EL_TAG_CTX_TAG(ctx) (ctx)->tag
510
511 #define INIT_ELEMENT_METHOD(obj) (obj)->init_element
512 #define INIT_ELEMENT_FLAGS(obj) (obj)->fill_flags
513 #define INIT_ELEMENT_DEFUN(obj, init_el, flags) \
514 { \
515 INIT_ELEMENT_METHOD(obj) = init_el; \
516 INIT_ELEMENT_FLAGS(obj) = (flags); \
517 INIT_EL_TAG_CTX_INIT(&(obj)->tag_ctx); \
518 }
519
520 #define INIT_OBJECT(object) (void)INIT_ELEMENT(NULL, object)
521 #define INIT_ELEMENT_DECL \
522 INIT_ELEMENT_FCT init_element; \
523 FLAGS fill_flags; \
524 INIT_EL_TAG_CTX tag_ctx
525
526 #define INIT_ELEMENT_INITIALIZER(init_el, flags) \
527 (init_el), (flags), { INIT_EL_TAG_DFLT, 0 }
528
529 #define INIT_ELEMENT(el_info, object) \
530 (INIT_ELEMENT_NEEDED(object) \
531 ? INIT_ELEMENT_METHOD(object)(el_info, (void *)object) : INIT_EL_TAG_DFLT)
532
533 #define INIT_ELEMENT_NEEDED(object) (INIT_ELEMENT_METHOD(object) != NULL)
534 #define INIT_ELEMENT_SETUP(el_info, object, tagvar, null_action, chg_action) \
535 { \
536 INIT_EL_TAG AI_init_el_tag; \
537 \
538 AI_init_el_tag = INIT_ELEMENT(el_info, object); \
539 if (AI_init_el_tag == (INIT_EL_TAG)INIT_EL_TAG_NULL) { \
540 (tagvar) = INIT_EL_TAG_NULL; \
541 null_action; \
542 } else if (AI_init_el_tag != (tagvar)) { \
543 (tagvar) = AI_init_el_tag; \
544 chg_action; \
545 } \
546 }
547
548 /* Initialization of single items of DOF-chain */
549 #define INIT_OBJECT_SINGLE(object) (void)INIT_ELEMENT_SINGLE(NULL, object)
550 #define INIT_ELEMENT_SINGLE(el_info, object) \
551 (INIT_ELEMENT_NEEDED((object)->unchained) \
552 ? INIT_ELEMENT_METHOD((object)->unchained)(el_info, (void *)object) \
553 : INIT_EL_TAG_DFLT)
554
555 /*******************************************************************************
556 * node projection descriptor:
557 *
558 * a function pointer which calculates the projected location of a
559 * new vertex resulting from refinement.
560 *
561 ******************************************************************************/
562
563 struct node_projection
564 {
565 void (*func)(REAL_D old_coord, const EL_INFO *eli, const REAL_B lambda);
566 };
567
568 /*******************************************************************************
569 * The geometric connectivity of periodic meshes is described by wall-
570 * transformations, affine isometries which map the current mesh to its
571 * periodic neighbour accross the wall of an element.
572 * The transformation operates as usual: Mx + t
573 ******************************************************************************/
574
575 struct aff_trafo
576 {
577 REAL_DD M;
578 REAL_D t;
579 };
580
581 /*******************************************************************************
582 * one single element (triangle) of the grid:
583 *******************************************************************************
584 *
585 * position of the nodes in 1d:
586 *
587 * 0 _____ 1 or 0 _____ 1
588 * 2
589 *
590 * child[0] child[1]
591 * refinement: 0 _____ 1 0 ___ 1 0 ___ 1
592 * 2 2
593 *
594 *******************************************************************************
595 *
596 * position of the nodes in 2d
597 * 2 2 2 2
598 * /\ or /\ or /\ or /\
599 * / \ 4/ \ 3 / \ 4/ \ 3
600 * / \ / \ / 3 \ / 6 \
601 * 0/______\1 0/______\1 0/______\1 0/______\1
602 * 5 5
603 *
604 * refinement: 2 child[0] 0 1 child[1]
605 * /\ /| |\
606 * 4/ \ 3 --> 5/ |4 3| \ 5
607 * / 6 \ /6 | |6 \
608 * 0/______\1 1/___| |___\0
609 * 5 3 2 2 4
610 *
611 *******************************************************************************
612 *
613 * 3d refinement: vertex numbering after (Baensch +) Kossaczky
614 *
615 * edges:
616 * E0: between V0, V1
617 * E1: between V0, V2
618 * E2: between V0, V3
619 * E3: between V1, V2
620 * E4: between V1, V3
621 * E5: between V2, V3
622 *
623 * Always edge 0 (between vertices 0 and 1) is bisected.
624 *
625 * V1
626 * -+
627 * ****- ||
628 * E0 ****-- | |
629 * ****-- | | E3
630 * ****-- | |
631 * ****-- | |
632 * V0 +. . . . . . . . . . . . . . | . . . |
633 * --- (E1) | + V2
634 * --- | /
635 * --- |E4 /
636 * --- | /
637 * E2 --- | / E5
638 * --- | /
639 * --- | /
640 * ---|/
641 * +
642 * V3
643 *
644 *******************************************************************************
645 * child: pointers to the two children of the element
646 * if (child[0]==NULL) element is a leaf of the
647 * tree
648 * dof: vector of pointers to dof vectors :-)
649 * new_coord: in case of curved boundary, coords of ref.edge midpoint
650 * index: global element index (only for test purposes)
651 * mark: element is a leaf:
652 * mark == 0 do not refine/coarsen
653 * mark > 0 refine (mark times)
654 * mark < 0 may be coarsened (mark times)
655 ******************************************************************************/
656
657 struct el
658 {
659 EL *child[2];
660 DOF **dof;
661 S_CHAR mark;
662 REAL *new_coord;
663
664 #if ALBERTA_DEBUG
665 int index;
666 #endif
667 };
668
669 /*******************************************************************************
670 * child_vertex_3d[el_type][child][i] =
671 * parent's local vertex index of new vertex i
672 * 4 stands for the newly generated vertex
673 *******************************************************************************
674 * child_edge_3d[el_type][child][i] =
675 * parent's local edge index of new edge i
676 * new edge 2 is half of old edge 0,
677 * new edges 4,5 are really new edges, and value is different:
678 * child_edge_3d[][][4,5] = index of same edge in other child
679 *******************************************************************************
680 * vertex_of_edge_?d[edge][i], i = 1,2 are the two vertices of edge
681 ******************************************************************************/
682 /* The numbering of the vertices of all sub-simplices is given by the
683 * following fields.
684 *
685 * Note that we favour a lexicographical ordering of wall (face)
686 * indices in favour of a cyclic ordering in 3d, while we use a cyclic
687 * ordering in 2d.
688 *
689 * Note that all arrays contain "excess" element; this way it is
690 * possible to avoid modulo calculus just to keep the array index in
691 * range, i.e. it is legal to index like follows:
692 *
693 * w_3d = 3;
694 * v_2d = 2;
695 * for (w = 0; w < N_WALLS_3D; w++) {
696 * for (v = 0; v < N_VERTICES_2D; v++) {
697 * DO SOMETHING WITH vertex_of_wall_3d[w_3d + w][v_2d+v]
698 * }
699 * }
700 */
701 static const int vertex_of_edge_1d[N_EDGES_1D][2*N_VERTICES_1D-1] = {
702 {0, 1, 0}
703 };
704 static const int vertex_of_wall_1d[2*N_WALLS_1D-1][N_VERTICES_0D] = {
705 { 1 }, { 0 }, { 1 }
706 };
707 static const int vertex_of_edge_2d[2*N_EDGES_2D-1][2*N_VERTICES_1D-1] = {
708 {1, 2, 1}, {2, 0, 2}, {0, 1, 0}, {1, 2, 1}, {2, 0, 2}
709 };
710 #define vertex_of_wall_2d vertex_of_edge_2d
711 static const int vertex_of_edge_3d[N_EDGES_3D][2*N_VERTICES_1D-1] = {
712 {0, 1, 0}, {0, 2, 0}, {0, 3, 0}, {1, 2, 1}, {1, 3, 1}, {2, 3, 2}
713 };
714 static const int vertex_of_wall_3d[2*N_WALLS_3D-1][2*N_VERTICES_2D-1] = {
715 {1, 2, 3, 1, 2}, {0, 2, 3, 0, 2}, {0, 1, 3, 0, 1}, {0, 1, 2, 0, 1},
716 {1, 2, 3, 1, 2}, {0, 2, 3, 0, 2}, {0, 1, 3, 0, 1},
717 };
718
719 /*******************************************************************************
720 * edge_of_vertices_3d[i][j]: gives the local index of edge with vertices i, j
721 ******************************************************************************/
722 static const int edge_of_vertices_3d[N_VERTICES_3D][N_VERTICES_3D] = {
723 {-1, 0, 1, 2 },
724 { 0, -1, 3, 4 },
725 { 1, 3, -1, 5 },
726 { 2, 4, 5, -1 }
727 };
728
729 /*******************************************************************************
730 * face_of_edge_3d[e][0/1]: gives the local number of the two adjacent faces
731 ******************************************************************************/
732 static const int face_of_edge_3d[N_EDGES_3D][2] = {
733 { 2, 3 }, { 1, 3 }, { 1, 2 }, { 0, 3 }, { 0, 2 }, { 0, 1 }
734 };
735
736 /* defined in element_Xd.c: permuted ordering of wall vertices to be
737 * able to match quadrature points of neighbouring elements on the
738 * separating wall.
739 *
740 * wall_orientation[_rel]_Xd() returns an index into those arrays.
741 */
742 extern
743 const int sorted_wall_vertices_1d[N_WALLS_1D][DIM_FAC_1D][2*N_VERTICES_0D-1];
744 extern
745 const int sorted_wall_vertices_2d[N_WALLS_2D][DIM_FAC_2D][2*N_VERTICES_1D-1];
746 extern
747 const int sorted_wall_vertices_3d[N_WALLS_3D][DIM_FAC_3D][2*N_VERTICES_2D-1];
748
749 /*******************************************************************************
750 * PARAMETRIC structure, entry in MESH structure
751 * description of parametric meshes and elements
752 ******************************************************************************/
753
754 /*
755 * The values to pass as "flags" argument to
756 * use_lagrange_parametric(..., flags)
757 */
758 typedef enum param_strategy {
759 PARAM_ALL = 0, /* all elements are to be parametric */
760 PARAM_CURVED_CHILDS = 1, /* bisection along {\lambda_0 = \lambda_1} */
761 PARAM_STRAIGHT_CHILDS = 2 /* bisection along straight lines/planes */
762 } PARAM_STRATEGY;
763
764 #define PARAM_STRATEGY_MASK \
765 (PARAM_ALL|PARAM_CURVED_CHILDS|PARAM_STRAIGHT_CHILDS)
766 #define PARAM_PERIODIC_COORDS 0x04 /* The coordinates themselves are
767 * periodic, normally parametric
768 * coordinates of a periodic mesh
769 * are _NOT_ periodic.
770 */
771
772 /* General hook-structure for parametric
773 * meshes. use_lagrange_parametric() populates this struture for the
774 * standard conforming iso-parametric geometry approximations.
775 */
776 struct parametric
777 {
778 const char *name; /* textual description analogous to BAS_FCTS. */
779
780 /* true: some elements may be non-parametric */
781 bool not_all;
782
783 /* true: standard routines coord_to_world, etc. may be used to get
784 * data about the reference triangulation. Set to "false" by
785 * default.
786 */
787 bool use_reference_mesh;
788
789 /* init_element(el_info, param) == false : non-parametric element,
790 * init_element(el_info, param) == true : parametric element.
791 *
792 * NOTE: If PARAMETRIC::init_element(el_info, ...) returns false,
793 * then it is supposed to fill el_info->coord with the current
794 * element's coordinate information despite the fact the "el_info"
795 * is _CONST_. This way the normal per-element functions can be used
796 * (e.g. el_det(), el_grd_lambda() etc.) instead of the parametric
797 * ones. This simplifies the program flow (and source code) for
798 * partially parametric meshes a _LOT_.
799 */
800 bool (*init_element)(const EL_INFO *el_info, const PARAMETRIC *parametric);
801
802 void (*vertex_coords)(EL_INFO *info);
803
804 void (*coord_to_world)(const EL_INFO *info, const QUAD *quad,
805 int n, const REAL_B lambda[], REAL_D *world);
806
807 /* Be careful with this function, for some world coordinates there
808 * are no barycentric coordinates. Works best if world is not too
809 * far away from our simplex. */
810 void (*world_to_coord)(const EL_INFO *info, int N,
811 const REAL_D world[],
812 REAL_B lambda[], int k[]);
813 void (*det)(const EL_INFO *info, const QUAD *quad,
814 int n, const REAL_B lambda[], REAL dets[]);
815 void (*grd_lambda)(const EL_INFO *info, const QUAD *quad,
816 int n, const REAL_B lambda[],
817 REAL_BD Lambda[], REAL_BDD DLambda[], REAL dets[]);
818 void (*grd_world)(const EL_INFO *info, const QUAD *quad,
819 int n, const REAL_B lambda[],
820 REAL_BD grd_Xtr[], REAL_BDB D2_Xtr[], REAL_BDBB D3_Xtr[]);
821 void (*wall_normal)(const EL_INFO *el_info, int wall,
822 const QUAD *wall_quad,
823 int n, const REAL_B lambda[],
824 REAL_D nu[], REAL_DB grd_nu[], REAL_DBB D2_nu[],
825 REAL dets[]);
826
827 /* inherit_parametric is used by get_submesh(), unchain_parametric()
828 * is used by unchain_submesh(). Can be left out if the sub-mesh
829 * feature is not used.
830 */
831 void (*inherit_parametric)(MESH *slave);
832 void (*unchain_parametric)(MESH *slave);
833
834 void *data; /* private data for specific implementations */
835 };
836
837 /*******************************************************************************
838 * EL_GEOM_CACHE structure; geometric information for non-parametric
839 * meshes which is not generated during mesh-traversal, but needed in
840 * several places. Data like det, Lambda, wall-normals, wall-det.
841 * Access to this cache _MUST_ go through fill_el_geom_cache(). See
842 * also fill_quad_el_cache(el_info), especially for parametric meshes.
843 ******************************************************************************/
844
845 struct el_geom_cache
846 {
847 EL *current_el;
848 FLAGS fill_flag;
849 REAL det;
850 REAL_BD Lambda;
851 int orientation[N_WALLS_MAX][2];
852 int rel_orientation[N_WALLS_MAX];
853 REAL wall_det[N_WALLS_MAX];
854 REAL_D wall_normal[N_WALLS_MAX];
855 };
856
857 #define FILL_EL_DET (1 << 0)
858 #define FILL_EL_LAMBDA (1 << 1)
859
860 #define FILL_EL_WALL_SHIFT(wall) (2 + 4*(wall))
861 #define FILL_EL_WALL_MASK(wall) (0x7 << FILL_EL_WALL_SHIFT(wall))
862
863 #define FILL_EL_WALL_DET(wall) (1 << (FILL_EL_WALL_SHIFT(wall)+0))
864 #define FILL_EL_WALL_NORMAL(wall) (1 << (FILL_EL_WALL_SHIFT(wall)+1))
865 #define FILL_EL_WALL_ORIENTATION(wall) (1 << (FILL_EL_WALL_SHIFT(wall)+2))
866 #define FILL_EL_WALL_REL_ORIENTATION(wall) (1 << (FILL_EL_WALL_SHIFT(wall)+3))
867
868 #define FILL_EL_WALL_DETS \
869 (FILL_EL_WALL_DET(0)|FILL_EL_WALL_DET(1)| \
870 FILL_EL_WALL_DET(2)|FILL_EL_WALL_DET(3))
871
872 #define FILL_EL_WALL_NORMALS \
873 (FILL_EL_WALL_NORMAL(0)|FILL_EL_WALL_NORMAL(1)| \
874 FILL_EL_WALL_NORMAL(2)|FILL_EL_WALL_NORMAL(3))
875
876 #define FILL_EL_WALL_ORIENTATIONS \
877 (FILL_EL_WALL_ORIENTATION(0)|FILL_EL_WALL_ORIENTATION(1)| \
878 FILL_EL_WALL_ORIENTATION(2)|FILL_EL_WALL_ORIENTATION(3))
879
880 #define FILL_EL_WALL_REL_ORIENTATIONS \
881 (FILL_EL_WALL_REL_ORIENTATION(0)|FILL_EL_WALL_REL_ORIENTATION(1)| \
882 FILL_EL_WALL_REL_ORIENTATION(2)|FILL_EL_WALL_REL_ORIENTATION(3))
883
884 static inline const EL_GEOM_CACHE *
885 fill_el_geom_cache(const EL_INFO *el_info, FLAGS fill_flag);
886
887 /*******************************************************************************
888 * additional information to elements during hierarchy traversal
889 *******************************************************************************
890 *
891 * mesh: pointer to the mesh structure
892 * coord: world coordinates of the vertices. For curved
893 * parametric meshes the corresponding information is
894 * filled by the function hooks in the PARAMETRIC
895 * structur.
896 * macro_el: pointer to the macro-element we belong to.
897 * el: node in the mesh-tree.
898 * parent: pointer to an EL_INFO structure describing the parent of this
899 * element in the mesh tree.
900 * fill_flag: copy of the fill-flags used to generate the EL_INFO
901 * structure.
902 * level: the depth in the mesh-tree, level 0 means root-level (i.e.
903 * an element of the macro triangulation).
904 *
905 * macro_wall: Mapping of the boundary facets of the element to
906 * the boundary facets of the containing macro
907 * element. macro_wall[w] == -1 means that the
908 * boundary facet number "w" is located in the
909 * interior of the containing macro-element, i.e.
910 * EL_INFO::macro_el.
911 *
912 * wall_bound: Boundary type of the co-dim 1 facets (all
913 * dimensions). Boundary types range from 0 (interior
914 * faces) to 255. Boundary types are just markers
915 * without interpretation, used to group boundary
916 * facets which share common properties. Needs
917 * FILL_BOUND. Boundary types can also be access via
918 * the EL_INFO::macro_wall[] component which is filled
919 * all the time. There is also a support function
920 * wall_bound() for this purpose.
921 * vertex_bound: Boundary type of the vertices. This is a bit-field:
922 * bit N is set if any of the co-dim 1 facets the
923 * vertex belongs to has boundary-type N. Needs FILL_BOUND.
924 * edge_bound: Boundary type of the edges (only 3d). Also a bit field,
925 * obtained in the same manner as vertex_bound. Needs
926 * FILL_BOUND.
927 *
928 * active_projection: node projection function for the new vertex
929 * which would result from a refinement of the current
930 * element. Needs FILL_PROJECTION.
931 *
932 * neigh: pointer to the adjacent elements NULL-pointer for a
933 * part of the boundary. Needs FILL_NEIGH.
934 * opp_coord: world coordinates of opposite vertices. Needs
935 * FILL_NEIGH|FILL_COORD.
936 * opp_vertex: local indices of opposite vertices. Needs FILL_NEIGH.
937 *
938 * el_type: type of the element, 0, 1, or 2. Only meaningful in 3d.
939 * orientation: orientation of the tetrahedron relative to the macro
940 * element. This is only set for 3d, otherwise it is fixed
941 * at 0. For DIM == 3 this gives the orientation
942 * w.r.t. to the standard co-ordinate frame, for higher
943 * dimenstion "absolute" orientation makes no sense; so
944 * "orientation" will be 1 for all macro elements, regard-
945 * less of their actual relative orientations.
946 *
947 * el_geom_cache: A cache to store derived quantities which are not
948 * computed during mesh-traversal, but are derived
949 * from the co-ordinate information. Access to the
950 * cache _must_ go through fill_el_geom_cache(). This
951 * stuff needs -- of course -- the FILL_COORDS
952 * fill-flag.
953 *
954 ******************************************************************************/
955
956 struct el_info
957 {
958 MESH *mesh;
959 REAL_D coord[N_VERTICES_MAX];
960 const MACRO_EL *macro_el;
961 EL *el;
962 const EL_INFO *parent;
963 FLAGS fill_flag;
964 int level;
965
966 S_CHAR macro_wall[N_WALLS_MAX];
967
968 BNDRY_TYPE wall_bound[N_WALLS_MAX];
969 BNDRY_FLAGS vertex_bound[N_VERTICES_MAX];
970 BNDRY_FLAGS edge_bound[N_EDGES_MAX];
971 #if DIM_MAX > 1
972 BNDRY_TYPE face_bound[MAX(1, N_FACES_MAX)];
973 #endif
974
975 const NODE_PROJ *active_projection;
976
977 EL *neigh[N_NEIGH_MAX];
978 S_CHAR opp_vertex[N_NEIGH_MAX];
979 REAL_D opp_coord[N_NEIGH_MAX];
980
981 U_CHAR el_type;
982 S_CHAR orientation;
983
984 struct master_info {
985 EL *el;
986 int opp_vertex;
987 REAL_D opp_coord;
988 U_CHAR el_type;
989 S_CHAR orientation;
990 } master, mst_neigh;
991
992 EL_GEOM_CACHE el_geom_cache;
993 };
994
995 /* Some "standard" bit-field operations, meant to hide the
996 * N_BNDRY_TYPES argument.
997 */
998 #define BNDRY_FLAGS_INIT(flags) bitfield_zap((flags), N_BNDRY_TYPES)
999 #define BNDRY_FLAGS_ALL(flags) bitfield_fill((flags), N_BNDRY_TYPES)
1000 #define BNDRY_FLAGS_CPY(to, from) bitfield_cpy((to), (from), N_BNDRY_TYPES)
1001 #define BNDRY_FLAGS_AND(to, from) bitfield_and((to), (from), N_BNDRY_TYPES)
1002 #define BNDRY_FLAGS_OR(to, from) bitfield_or((to), (from), N_BNDRY_TYPES)
1003 #define BNDRY_FLAGS_XOR(to, from) bitfield_xor((to), (from), N_BNDRY_TYPES)
1004 #define BNDRY_FLAGS_CMP(a, b) bitfield_cmp((a), (b), N_BNDRY_TYPES)
1005
1006 /* bit 0 flags boundary segments, if not set we are in the interior */
1007 #define BNDRY_FLAGS_IS_INTERIOR(mask) (!bitfield_tst((mask), 0))
1008
1009 /* Set bit 0 to mark this as a boundary bit-mask. */
1010 #define BNDRY_FLAGS_MARK_BNDRY(flags) bitfield_set((flags), INTERIOR)
1011
1012 /* Return TRUE if SEGMENT has BIT set _and_ BIT != 0. */
1013 #define BNDRY_FLAGS_IS_AT_BNDRY(segment, bit) \
1014 ((bit) && bitfield_tst((segment), (bit)))
1015
1016 /* Set a bit in the boundary-type mask. The precise meaning of BIT:
1017 *
1018 * BIT == 0: clear the boundary mask (meaning: interior node)
1019 * BIT > 0: set bit BIT and also bit 0 (meaning: boundary node)
1020 */
1021 #define BNDRY_FLAGS_SET(flags, bit) \
1022 if ((bit) != INTERIOR) { \
1023 bitfield_set((flags), INTERIOR); \
1024 bitfield_set((flags), (bit)); \
1025 } else { \
1026 BNDRY_FLAGS_INIT(flags); \
1027 }
1028 /* return TRUE if SEGMENT and MASK have non-zero overlap */
1029 #define BNDRY_FLAGS_IS_PARTOF(segment, mask) \
1030 bitfield_andp((segment), (mask), 1 /* offset */, N_BNDRY_TYPES)
1031 /* FindFirstBoundaryBit, return INTERIOR for interior nodes, otherwise the
1032 * number of the first bit set in MASK.
1033 */
1034 #define BNDRY_FLAGS_FFBB(mask) bitfield_ffs(mask, 1 /* offset */, N_BNDRY_TYPES)
1035
1036 /*******************************************************************************
1037 * RC_LIST_EL structure to describe a refinement/coarsening patch.
1038 * el_info: contains information about the patch element. This is not
1039 * a pointer since EL_INFO structures are often overwritten
1040 * during mesh traversal.
1041 * no: index of the patch element in the patch.
1042 * flags: see the RCLE_... defines below for a description.
1043 * neigh: neighbours to the right/left in the orientation of the
1044 * edge, or NULL pointer for a boundary face. (dim == 3 only)
1045 * opp_vertex: the opposite vertex of neigh[0/1]. (dim == 3 only)
1046 ******************************************************************************/
1047
1048 struct rc_list_el
1049 {
1050 EL_INFO el_info;
1051 int no;
1052 FLAGS flags;
1053 RC_LIST_EL *neigh[2];
1054 int opp_vertex[2];
1055 };
1056
1057 /* Valid settings for RC_LIST_EL->flags. The "PERIODIC" flags can be
1058 * exploited in refine_interpol/coarse_restrict routines.
1059 */
1060 #define RCLE_NONE 0x0 /* just nothing special */
1061 #define RCLE_COARSE_EDGE_COMPAT (1 << 0) /* set if the coarsening edge
1062 * of the patch element is
1063 * the coarsening edge of the
1064 * patch. Only for internal
1065 * use.
1066 */
1067
1068 /*******************************************************************************
1069 * flags, which information should be present in the EL_INFO structure
1070 ******************************************************************************/
1071
1072 #define FILL_NOTHING 0x0000L
1073 #define FILL_COORDS 0x0001L
1074 #define FILL_BOUND 0x0002L
1075 #define FILL_NEIGH 0x0004L
1076 #define FILL_OPP_COORDS 0x0008L
1077 #define FILL_ORIENTATION 0x0010L
1078 #define FILL_PROJECTION 0x0020L
1079 #define FILL_MACRO_WALLS 0x0040L
1080 #define FILL_WALL_MAP FILL_MACRO_WALLS
1081 #define FILL_NON_PERIODIC 0x0080L
1082 #define FILL_MASTER_INFO 0x0100L
1083 #define FILL_MASTER_NEIGH 0x0200L
1084
1085 #define FILL_ANY \
1086 (FILL_COORDS|FILL_BOUND|FILL_NEIGH|FILL_OPP_COORDS| \
1087 FILL_ORIENTATION|FILL_PROJECTION|FILL_MACRO_WALLS| \
1088 FILL_NON_PERIODIC|FILL_MASTER_INFO|FILL_MASTER_NEIGH)
1089
1090 /*******************************************************************************
1091 * flags for mesh traversal
1092 ******************************************************************************/
1093
1094 #define CALL_EVERY_EL_PREORDER 0x010000L
1095 #define CALL_EVERY_EL_INORDER 0x020000L
1096 #define CALL_EVERY_EL_POSTORDER 0x040000L
1097 #define CALL_LEAF_EL 0x080000L
1098 #define CALL_LEAF_EL_LEVEL 0x100000L
1099 #define CALL_EL_LEVEL 0x200000L
1100 #define CALL_MG_LEVEL 0x400000L /* used in multigrid methods */
1101
1102 #define TEST_FLAG(flags, el_info) \
1103 TEST_EXIT(!((((el_info)->fill_flag)^(flags)) & (flags)), \
1104 "flag "#flags" not set\n")
1105
1106 #if ALBERTA_DEBUG==1
1107 # define DEBUG_TEST_FLAG(flags, el_info) \
1108 if((((el_info)->fill_flag)^(flags)) & (flags)) \
1109 print_error_funcname(funcName, __FILE__, __LINE__), \
1110 print_error_msg_exit("flag "#flags" not set\n")
1111 #else
1112 # define DEBUG_TEST_FLAG(flags, el_info) do { funcName = funcName; } while (0)
1113 #endif
1114 /*******************************************************************************
1115 * one single element of the macro triangulation:
1116 *******************************************************************************
1117 * el: pointer to the element data of the macro element
1118 * coord: world coordinates of the nodes on the macro element
1119 * wall_bound: Boundary type of the co-dim 1 facets (all
1120 * dimensions). Boundary types range from 0 (interior
1121 * faces) to 127. Boundary types are just a markers
1122 * without interpretation, used to group boundary
1123 * facets which share common properties.
1124 * vertex_bound: Boundary type of the vertices. This is a bit-field:
1125 * bit N is set if any of the co-dim 1 facets the
1126 * vertex belongs to has boundary-type N.
1127 * edge_bound: Boundary type of the edges (only 3d). Also a bit field,
1128 * obtained in the same manner as vertex_bound.
1129 * projection: possible node projection functions for all nodes [0]
1130 * or for specific edges or faces (dim > 1), which will
1131 * override entry [0].
1132 * index: unique global index of macro element
1133 * neigh: pointer to the adjacent macro elements
1134 * NULL-pointer for a part of the boundary
1135 * opp_vertex: local index of opposite vertex w.r.t. neighbour numbering
1136 * neigh_vertices: local indices of common vertices of the periodic
1137 * neighbour, this component is set only for the virtual
1138 * neighbours on periodic meshes.
1139 * neigh_vertices[wall][loc_idx] is the local vertex number
1140 * on the neighbour the vertex with local number
1141 * (wall + 1 + loc_idx) % N_VERTICES(MESH_DIM) on this
1142 * element is mapped to.
1143 * wall_trafo: only for periodic meshes: the affine transformation which
1144 * maps the mesh across the corresponding wall to the
1145 * neighbour facet. The wall transformation must be
1146 * affine isometries.
1147 * np_vertex_bound: boundary type of the vertices when treating a periodic
1148 * mesh as non-periodic
1149 * np_edge_bound: like np_vertex_bound
1150 * el_type: type of corresponding element.
1151 * orientation: orientation of corresponding element, used in 3d.
1152 *
1153 ******************************************************************************/
1154
1155 struct macro_el
1156 {
1157 EL *el;
1158 REAL_D *coord[N_VERTICES_MAX];
1159
1160 BNDRY_TYPE wall_bound[N_WALLS_MAX];
1161 BNDRY_FLAGS vertex_bound[N_VERTICES_MAX];
1162 #if DIM_MAX > 1
1163 BNDRY_FLAGS edge_bound[N_EDGES_MAX];
1164 #endif
1165 #if DIM_MAX > 2
1166 BNDRY_TYPE face_bound[N_FACES_MAX];
1167 #endif
1168
1169 NODE_PROJ *projection[N_NEIGH_MAX + 1];
1170
1171 int index;
1172
1173 MACRO_EL *neigh[N_NEIGH_MAX];
1174 S_CHAR opp_vertex[N_NEIGH_MAX];
1175 S_CHAR neigh_vertices[N_NEIGH_MAX][N_VERTICES(DIM_MAX-1)];
1176 AFF_TRAFO *wall_trafo[N_NEIGH_MAX];
1177 BNDRY_FLAGS np_vertex_bound[N_VERTICES_MAX];
1178 #if DIM_MAX > 1
1179 BNDRY_FLAGS np_edge_bound[N_EDGES_MAX];
1180 #endif
1181
1182 S_CHAR orientation;
1183
1184 U_CHAR el_type;
1185
1186 /* The chain to the master macro element if we belong to a trace-mesh */
1187 struct {
1188 MACRO_EL *macro_el;
1189 S_CHAR opp_vertex;
1190 BNDRY_FLAGS vertex_bound[MAX(1, N_VERTICES(DIM_MAX-1))];
1191 BNDRY_FLAGS np_vertex_bound[MAX(1, N_VERTICES(DIM_MAX-1))];
1192 #if DIM_MAX > 1
1193 BNDRY_FLAGS edge_bound[N_EDGES(MAX(1, DIM_MAX-1))];
1194 BNDRY_FLAGS np_edge_bound[N_EDGES(MAX(1, DIM_MAX-1))];
1195 #endif
1196 } master;
1197 };
1198
1199 /* Some support functions to access boundary-facet related data only
1200 * stored on the macro-element level.
1201 */
1202
wall_bound(const EL_INFO * el_info,int wall)1203 static inline BNDRY_TYPE wall_bound(const EL_INFO *el_info, int wall)
1204 {
1205 int mwall = el_info->macro_wall[wall];
1206
1207 if (mwall < 0) {
1208 return INTERIOR;
1209 }
1210
1211 if ((el_info->fill_flag & FILL_NON_PERIODIC)) {
1212 return el_info->macro_el->wall_bound[mwall];
1213 }
1214
1215 if (el_info->macro_el->neigh_vertices[mwall][0] < 0) {
1216 return el_info->macro_el->wall_bound[mwall];
1217 } else {
1218 return INTERIOR;
1219 }
1220 }
1221
wall_trafo(const EL_INFO * el_info,int wall)1222 static inline const AFF_TRAFO *wall_trafo(const EL_INFO *el_info, int wall)
1223 {
1224 int mwall;
1225
1226 if ((el_info->fill_flag & FILL_NON_PERIODIC)) {
1227 return NULL;
1228 }
1229
1230 mwall = el_info->macro_wall[wall];
1231
1232 return mwall < 0 ? NULL : el_info->macro_el->wall_trafo[mwall];
1233 }
1234
wall_proj(const EL_INFO * el_info,int wall)1235 static inline const NODE_PROJ *wall_proj(const EL_INFO *el_info, int wall)
1236 {
1237 if (wall < 0) {
1238 return el_info->macro_el->projection[0];
1239 } else {
1240 int mwall = el_info->macro_wall[wall];
1241
1242 return el_info->macro_el->projection[mwall+1];
1243 }
1244 }
1245
1246 /*******************************************************************************
1247 * index based storage of macro triangulations
1248 ******************************************************************************/
1249
1250 struct macro_data
1251 {
1252 int dim; /* dimension of the elements */
1253
1254 int n_total_vertices;
1255 int n_macro_elements;
1256
1257 REAL_D *coords; /* Length will be n_total_vertices */
1258
1259 int *mel_vertices; /* mel_vertices[i*N_VERTICES(dim)+j]:
1260 * global index of jth vertex of element i
1261 */
1262 int *neigh; /* neigh[i*N_NEIGH(dim)+j]:
1263 * neighbour j of element i or -1 at boundaries
1264 */
1265 int *opp_vertex; /* opp_vertex[i*N_NEIGH(dim)+j]: if set (need not
1266 * be) the local vertex number w.r.t. the neighbour
1267 * of the vertex opposit the separating wall.
1268 */
1269 BNDRY_TYPE *boundary; /* boundary[i*N_NEIGH(dim)+j]:
1270 * boundary type of jth co-dim 1 facet of element i
1271 *
1272 * WARNING: In 1D the local index corresponds
1273 * to vertex 1 & vice versa! (Consistent with
1274 * macro_data.neigh)
1275 */
1276 U_CHAR *el_type; /* el_type[i]: type of element i only used in 3d! */
1277 int (*wall_vtx_trafos)[N_VERTICES(DIM_MAX-1)][2]; /* the wall trafos */
1278 /* Wall transformations are in terms of mappings between
1279 * vertices. i-th wall trafo: global vertex number
1280 * wall_vtx_trafos[i][v][0] maps to wall_vtx_trafos[i][v][1], v loops
1281 * through the local vertex number of the respective wall.
1282 */
1283 int n_wall_vtx_trafos;/* for periodic meshes: number of
1284 * combinatorical wall trafos.
1285 */
1286 int *el_wall_vtx_trafos;
1287 /* el_wall_vtx_trafos[i*N_WALLS(dim)+j] number of the wall
1288 * transformation of the j-th wall for the i-th element. > 0:
1289 * #wall_trafo+1. < 0: inverse of -(#wall_trafo+1)
1290 */
1291 AFF_TRAFO *wall_trafos; /* The group generators of the space group
1292 * defining the periodic structure of the
1293 * mesh.
1294 */
1295 int n_wall_trafos;
1296 int *el_wall_trafos; /* N = el_wall_trafos[i*N_NEIGH(dim)+j]:
1297 *
1298 * number of the wall transformation mapping to
1299 * the neighbouring fundamental domain across
1300 * the given wall.
1301 *
1302 * If negative: inverse of generator -N-1
1303 * If positive: generator +N-1
1304 */
1305 #if ALBERTA_DEBUG
1306 char **mel_comment; /* for debugging */
1307 #endif
1308 };
1309
1310 #ifndef DOF_ADMIN_DEF
1311 # define DOF_ADMIN_DEF
1312
1313 /*******************************************************************************
1314 * dof handling
1315 ******************************************************************************/
1316 /* presumably the largest native integer type */
1317 # define DOF_FREE_UNIT_TYPE long
1318 typedef unsigned DOF_FREE_UNIT_TYPE DOF_FREE_UNIT;
1319 # define DOF_FREE_SIZE ((int)(8*sizeof(DOF_FREE_UNIT)))
1320 # define DOF_UNIT_ALL_FREE (~0UL)
1321 extern const DOF_FREE_UNIT dof_free_bit[DOF_FREE_SIZE]; /* in dof_admin.c */
1322 # define DOF_UNUSED (-1) /* el->dof[][] == DOF_UNUSED, mark unused DOFs */
1323
1324 # define FOR_ALL_DOFS(admin, todo) \
1325 if ((admin)->hole_count == 0) { \
1326 int dof; \
1327 \
1328 for (dof = 0; dof < (admin)->used_count; dof++) { \
1329 todo; \
1330 } \
1331 } else { \
1332 DOF_FREE_UNIT _dfu, *_dof_free = (admin)->dof_free; \
1333 int _i, _ibit, dof=0; \
1334 int _n= ((admin)->size_used + DOF_FREE_SIZE-1) / DOF_FREE_SIZE; \
1335 \
1336 for (_i = 0; _i < _n; _i++) { \
1337 if ((_dfu = _dof_free[_i])) { \
1338 if (_dfu == DOF_UNIT_ALL_FREE) { \
1339 dof += DOF_FREE_SIZE; \
1340 } else { \
1341 for (_ibit = 0; \
1342 _ibit < DOF_FREE_SIZE; \
1343 _ibit++, dof++, _dfu >>= 1) { \
1344 if ((_dfu & 1) == 0) { \
1345 todo; \
1346 } \
1347 } \
1348 } \
1349 } else { \
1350 for (_ibit = 0; _ibit < DOF_FREE_SIZE; _ibit++, dof++) { \
1351 todo; \
1352 } \
1353 } \
1354 } \
1355 }
1356
1357 # define FOR_ALL_FREE_DOFS(admin, todo) \
1358 if ((admin)->hole_count == 0) { \
1359 int dof; \
1360 for (dof = (admin)->used_count; dof < (admin)->size; dof++) { \
1361 todo; \
1362 } \
1363 } else { \
1364 DOF_FREE_UNIT _dfu, *_dof_free = (admin)->dof_free; \
1365 int _i, _ibit, dof=0; \
1366 int _n= ((admin)->size + DOF_FREE_SIZE-1) / DOF_FREE_SIZE; \
1367 \
1368 for (_i = 0; _i < _n; _i++) { \
1369 if ((_dfu = _dof_free[_i])) { \
1370 if (_dfu == DOF_UNIT_ALL_FREE) { \
1371 for (_ibit = 0 ; _ibit < DOF_FREE_SIZE; _ibit++, dof++) { \
1372 todo; \
1373 } \
1374 } else { \
1375 for (_ibit = 0; \
1376 _ibit < DOF_FREE_SIZE; \
1377 _ibit++, dof++, _dfu >>= 1) { \
1378 if ((_dfu & 1) != 0) { \
1379 todo; \
1380 } \
1381 } \
1382 } \
1383 } else { \
1384 dof += DOF_FREE_SIZE; \
1385 } \
1386 } \
1387 }
1388
1389 /* Stop if dof >= size_used */
1390 # define FOR_ALL_USED_FREE_DOFS(admin, todo) \
1391 FOR_ALL_FREE_DOFS(admin, \
1392 if (dof >= admin->size_used) { \
1393 break; \
1394 } \
1395 todo)
1396
1397 /* Possible values for DOF_ADMIN->flags */
1398 # define ADM_FLAGS_DFLT 0 /* nothing special */
1399 # define ADM_PRESERVE_COARSE_DOFS (1 << 0) /* preserve non-leaf DOFs */
1400 # define ADM_PERIODIC (1 << 1) /* periodic ADMIN on a
1401 * periodic mesh
1402 */
1403 #define ADM_FLAGS_MASK (ADM_PRESERVE_COARSE_DOFS | ADM_PERIODIC)
1404
1405 struct dof_admin
1406 {
1407 MESH *mesh;
1408 const char *name;
1409
1410 DOF_FREE_UNIT *dof_free; /* flag bit vector */
1411 unsigned int dof_free_size;/* flag bit vector size */
1412 unsigned int first_hole; /* index of first non-zero dof_free entry */
1413
1414 FLAGS flags;
1415
1416 DOF size; /* allocated size of dof_list vector */
1417 DOF used_count; /* number of used dof indices */
1418 DOF hole_count; /* number of FREED dof indices (NOT size-used)*/
1419 DOF size_used; /* > max. index of a used entry */
1420
1421 int n_dof[N_NODE_TYPES]; /* dofs from THIS dof_admin */
1422 int n0_dof[N_NODE_TYPES]; /* start of THIS admin's DOFs in the mesh. */
1423 /****************************************************************************/
1424 DOF_INT_VEC *dof_int_vec; /* linked list of int vectors */
1425 DOF_DOF_VEC *dof_dof_vec; /* linked list of dof vectors */
1426 DOF_DOF_VEC *int_dof_vec; /* linked list of dof vectors */
1427 DOF_UCHAR_VEC *dof_uchar_vec; /* linked list of u_char vectors */
1428 DOF_SCHAR_VEC *dof_schar_vec; /* linked list of s_char vectors */
1429 DOF_REAL_VEC *dof_real_vec; /* linked list of real vectors */
1430 DOF_REAL_D_VEC *dof_real_d_vec; /* linked list of real_d vectors */
1431 DOF_REAL_DD_VEC *dof_real_dd_vec; /* linked list of real_d vectors */
1432 DOF_PTR_VEC *dof_ptr_vec; /* linked list of void * vectors */
1433 DOF_MATRIX *dof_matrix; /* linked list of matrices */
1434
1435 DBL_LIST_NODE compress_hooks; /* linked list of custom compress
1436 * handlers.
1437 */
1438
1439 /*******************************************************************************
1440 * pointer for administration; don't touch!
1441 ******************************************************************************/
1442
1443 void *mem_info;
1444 };
1445
1446 /* DOF_COMP_HOOK is a linked list rooted in
1447 * DOF_ADMIN->compress_hooks. The user may install arbitrary many
1448 * custom compress-handlers via add_dof_compress_hook(),
1449 * del_dof_compress_hook().
1450 */
1451 struct dof_comp_hook
1452 {
1453 DBL_LIST_NODE node; /* our link to the compress_hooks list */
1454 void (*handler)(DOF first, DOF last, const DOF *new_dof, void *app_data);
1455 void *application_data;
1456 };
1457
1458 /*******************************************************************************
1459 * dof vector structures
1460 *******************************************************************************
1461 * next: pointer to next structure containing vector of same type
1462 * fe_space: pointer to fe_space structure
1463 * refine_interpol: dof interpolation during refinement
1464 * coarse_restrict: restriction of linear functionals evaluated on a finer
1465 * grid and stored in dof vector to the coarser grid
1466 * during coarsening
1467 * or dof interpolation during coarsening
1468 * size: allocated size of vector
1469 * vec[]: vector entries (entry is used if dof index is used)
1470 ******************************************************************************/
1471
1472 #define UCHAR_name uchar
1473 #define uchar_VECNAME SCHAR
1474 #define SCHAR_name schar
1475 #define schar_VECNAME SCHAR
1476 #define INT_name int
1477 #define int_VECNAME INT
1478 #define DOF_name dof
1479 #define dof_VECNAME INT
1480 #define PTR_name ptr
1481 #define ptr_VECNAME PTR
1482 #define REAL_name real
1483 #define real_VECNAME REAL
1484 #define REAL_D_name real_d
1485 #define real_d_VECNAME REAL_D
1486 #define REAL_DD_name real_dd
1487 #define real_dd_VECNAME REAL_DD
1488 #define BNDRY_name bndry
1489 #define bndry_VECNAME BNDRY
1490
1491 # define DECL_DOF_VEC(VECNAME, vectype) \
1492 struct CPP_CONCAT3(dof_, VECNAME##_name, _vec) \
1493 { \
1494 DOF_##VECNAME##_VEC *next; \
1495 const FE_SPACE *fe_space; \
1496 \
1497 const char *name; \
1498 \
1499 DOF size; \
1500 int reserved; \
1501 \
1502 vectype *vec; \
1503 \
1504 void (*refine_interpol)(DOF_##VECNAME##_VEC *, RC_LIST_EL *, int n); \
1505 void (*coarse_restrict)(DOF_##VECNAME##_VEC *, RC_LIST_EL *, int n); \
1506 void *user_data; \
1507 \
1508 DBL_LIST_NODE chain; \
1509 const DOF_##VECNAME##_VEC *unchained; \
1510 \
1511 EL_##VECNAME##_VEC *vec_loc; \
1512 \
1513 void *mem_info; /*pointer for administration; don't touch! */ \
1514 }; \
1515 typedef vectype VECNAME##_VEC_TYPE
1516
1517 DECL_DOF_VEC(INT, int);
1518 DECL_DOF_VEC(DOF, DOF);
1519 DECL_DOF_VEC(UCHAR, U_CHAR);
1520 DECL_DOF_VEC(SCHAR, S_CHAR);
1521 DECL_DOF_VEC(PTR, void *);
1522 DECL_DOF_VEC(REAL, REAL);
1523 DECL_DOF_VEC(REAL_D, REAL_D);
1524 DECL_DOF_VEC(REAL_DD, REAL_DD);
1525
1526 /* A finite element function/coefficient vector for DIM_OF_WORLD
1527 * problems. If the basis-functions are vector-valued themselves, then
1528 * the vector is actually REAL-valued (stride == 1), otherwise REAL_D
1529 * valued (stride == DIM_OF_WORLD).
1530 *
1531 * So: if stride == 1, then this is actually a DOF_REAL_VEC, if stride ==
1532 * DIM_OF_WORLD, then this is actually a DOF_REAL_D_VEC.
1533 */
1534 struct dof_real_vec_d
1535 {
1536 DOF_REAL_VEC_D *next;
1537 const FE_SPACE *fe_space;
1538
1539 const char *name;
1540
1541 DOF size;
1542 int stride; /* either 1 or DIM_OF_WORLD */
1543
1544 REAL *vec;
1545
1546 void (*refine_interpol)(DOF_REAL_VEC_D *, RC_LIST_EL *, int n);
1547 void (*coarse_restrict)(DOF_REAL_VEC_D *, RC_LIST_EL *, int n);
1548 void *user_data;
1549
1550 DBL_LIST_NODE chain;
1551 const DOF_REAL_VEC_D *unchained;
1552
1553 EL_REAL_VEC_D *vec_loc;
1554
1555 void *mem_info; /*pointer for administration; don't touch! */
1556 };
1557
1558 typedef REAL REAL_VEC_D_TYPE; /* needed?? */
1559
1560 /*******************************************************************************
1561 * sparse matrix with one row for each dof,
1562 * entries are either REAL or REAL_DD
1563 *******************************************************************************
1564 * next: pointer to next matrix (linked list in MESH)
1565 * matrix_row[]: pointers to row structures (or NULL if row index is unused)
1566 * size: currently allocated size of matrix_row[]
1567 ******************************************************************************/
1568
1569 /* "flag" values for "type" component */
1570 typedef enum matent_type {
1571 MATENT_NONE = -1,
1572 MATENT_REAL = 0,
1573 MATENT_REAL_D = 1,
1574 MATENT_REAL_DD = 2
1575 } MATENT_TYPE;
1576
1577 static const size_t matent_size[4] = {
1578 0, sizeof(REAL), sizeof(REAL_D), sizeof(REAL_DD)
1579 };
1580 # define MATENT_SIZE(type) matent_size[(type)+1]
1581
1582 struct dof_matrix
1583 {
1584 DOF_MATRIX *next;
1585 const FE_SPACE *row_fe_space;
1586 const FE_SPACE *col_fe_space;
1587 const char *name;
1588
1589 MATRIX_ROW **matrix_row; /* lists of matrix entries */
1590 DOF size; /* size of vector matrix_row */
1591 MATENT_TYPE type; /* type of matrix entries. */
1592 size_t n_entries; /* total number of entries in the
1593 * matrix, updated by
1594 * add_element_matrix(), set to 0 by
1595 * clear_dof_matrix().
1596 */
1597 bool is_diagonal;
1598 union {
1599 DOF_REAL_VEC *real;
1600 DOF_REAL_D_VEC *real_d;
1601 DOF_REAL_DD_VEC *real_dd;
1602 } diagonal; /* The diagonal entries, if is_diagonal == true */
1603 DOF_INT_VEC *diag_cols; /* The column indices of the diagonal entries */
1604 union {
1605 DOF_REAL_VEC *real;
1606 DOF_REAL_D_VEC *real_d;
1607 DOF_REAL_DD_VEC *real_dd;
1608 } inv_diag; /* a cache for the diagonal entries, may be NULL. */
1609
1610 BNDRY_FLAGS dirichlet_bndry; /* bndry-type bit-mask for
1611 * Dirichlet-boundary conditions built
1612 * into the matrix
1613 */
1614
1615 void (*refine_interpol)(DOF_MATRIX *, RC_LIST_EL *, int n);
1616 void (*coarse_restrict)(DOF_MATRIX *, RC_LIST_EL *, int n);
1617
1618 /* The list pointers for the block-matrix structure to support
1619 * direct sums of fe-spaces.
1620 */
1621 DBL_LIST_NODE row_chain;
1622 DBL_LIST_NODE col_chain;
1623 const DOF_MATRIX *unchained;
1624
1625 void *mem_info;
1626 };
1627
1628 /*******************************************************************************
1629 * row structure for sparse matrix, with either REAL or REAL_DD entries.
1630 *******************************************************************************
1631 * next: pointer to next structure containing entries of same row
1632 * col[]: column indices of entries (if >= 0; else unused)
1633 * entry[]: matrix entries
1634 ******************************************************************************/
1635
1636 # define ROW_LENGTH 9
1637
1638 /* The actual size of this structure is determined by the type of the
1639 * matrix entries. The correct length is allocated in
1640 * get_matrix_row().
1641 */
1642 # define SIZEOF_MATRIX_ROW(type) \
1643 (sizeof(MATRIX_ROW) - sizeof(REAL_DD) + ROW_LENGTH*sizeof(type))
1644
1645 struct matrix_row
1646 {
1647 MATRIX_ROW *next;
1648 MATENT_TYPE type;
1649 DOF col[ROW_LENGTH]; /* column indices */
1650 union {
1651 REAL real[1];
1652 REAL_D real_d[1];
1653 REAL_DD real_dd[1];
1654 } entry;
1655 };
1656
1657 struct matrix_row_real
1658 {
1659 MATRIX_ROW_REAL *next;
1660 MATENT_TYPE type;
1661 DOF col[ROW_LENGTH]; /* column indices */
1662 REAL entry[ROW_LENGTH]; /* matrix entries */
1663 };
1664
1665 struct matrix_row_real_d
1666 {
1667 MATRIX_ROW_REAL_D *next;
1668 MATENT_TYPE type;
1669 DOF col[ROW_LENGTH]; /* column indices */
1670 REAL_D entry[ROW_LENGTH]; /* matrix entries */
1671 };
1672
1673 struct matrix_row_real_dd
1674 {
1675 MATRIX_ROW_REAL_DD *next;
1676 MATENT_TYPE type;
1677 DOF col[ROW_LENGTH]; /* column indices */
1678 REAL_DD entry[ROW_LENGTH]; /* matrix entries */
1679 };
1680
1681 # define ENTRY_USED(col) ((col) >= 0)
1682 # define ENTRY_NOT_USED(col) ((col) < 0)
1683 # define UNUSED_ENTRY -1
1684 # define NO_MORE_ENTRIES -2
1685
1686 # ifndef __CBLAS_H__
1687 typedef enum { NoTranspose,
1688 Transpose,
1689 ConjugateTranspose } MatrixTranspose;
1690 # endif
1691
1692 /* In C++ we would call this construct an iterator ... */
1693 # define FOR_ALL_MAT_COLS(type, matrow, what) \
1694 { \
1695 MATRIX_ROW_##type *row; \
1696 int col_idx; \
1697 DOF col_dof; \
1698 \
1699 for (row = (MATRIX_ROW_##type *)(matrow); row; row = row->next) { \
1700 for (col_idx = 0; col_idx < ROW_LENGTH; col_idx++) { \
1701 col_dof = row->col[col_idx]; \
1702 if (ENTRY_USED(col_dof)) { \
1703 what; \
1704 } else if (col_dof == NO_MORE_ENTRIES) { \
1705 break; \
1706 } \
1707 } \
1708 if (col_dof == NO_MORE_ENTRIES) { \
1709 break; \
1710 } \
1711 } \
1712 }
1713
1714 #endif /* DOF_ADMIN_DEF */
1715
1716 typedef struct el_vec_head
1717 {
1718 int n_components;
1719 int n_components_max;
1720 DBL_LIST_NODE chain;
1721 int reserved;
1722 } EL_VEC_HEAD;
1723
1724 #define DECL_DOF_EL_VEC(VECNAME, vectype) \
1725 struct CPP_CONCAT3(el_,VECNAME##_name, _vec) \
1726 { \
1727 int n_components; \
1728 int n_components_max; \
1729 DBL_LIST_NODE chain; \
1730 int reserved; \
1731 vectype vec[1]; \
1732 }; \
1733 typedef vectype EL_##VECNAME##_VEC_TYPE
1734
1735 DECL_DOF_EL_VEC(INT, int);
1736 DECL_DOF_EL_VEC(DOF, DOF);
1737 DECL_DOF_EL_VEC(UCHAR, U_CHAR);
1738 DECL_DOF_EL_VEC(SCHAR, S_CHAR);
1739 DECL_DOF_EL_VEC(BNDRY, BNDRY_FLAGS);
1740 DECL_DOF_EL_VEC(PTR, void *);
1741 DECL_DOF_EL_VEC(REAL, REAL);
1742 DECL_DOF_EL_VEC(REAL_D, REAL_D);
1743 DECL_DOF_EL_VEC(REAL_DD, REAL_DD);
1744
1745 struct el_real_vec_d
1746 {
1747 int n_components;
1748 int n_components_max;
1749 DBL_LIST_NODE chain;
1750 int stride; /* either 1 or DIM_OF_WORLD */
1751 REAL vec[1];
1752 };
1753 typedef REAL EL_REAL_VEC_D_TYPE;
1754
1755 #undef DECL_DOF_EL_VEC
1756
1757 /*******************************************************************************
1758 * Here comes the MESH (giving access to the whole triangulation)
1759 ******************************************************************************/
1760
1761 struct mesh
1762 {
1763 const char *name;
1764
1765 int dim;
1766
1767 int n_vertices;
1768 int n_elements;
1769 int n_hier_elements;
1770
1771 int n_edges; /* Only used for dim > 1 */
1772 int n_faces; /* Only used for dim == 3 */
1773 int max_edge_neigh; /* Only used for dim == 3 */
1774
1775 bool is_periodic; /* true if it is possible to define periodic*/
1776 int per_n_vertices; /* DOF_ADMINS on this mesh. The per_n_... */
1777 int per_n_edges; /* entries count the number of quantities on*/
1778 int per_n_faces; /* the periodic mesh (i.e. n_faces counts */
1779 /* periodic faces twice, n_per_faces not). */
1780 AFF_TRAFO *const*wall_trafos;
1781 int n_wall_trafos;
1782
1783 int n_macro_el;
1784 MACRO_EL *macro_els;
1785
1786 REAL_D bbox[2]; /* bounding box for the mesh */
1787 REAL_D diam; /* bbox[1] - bbox[0] */
1788
1789 PARAMETRIC *parametric;
1790
1791 DOF_ADMIN **dof_admin;
1792 int n_dof_admin;
1793
1794 int n_dof_el; /* sum of all dofs from all admins */
1795 int n_dof[N_NODE_TYPES]; /* sum of vertex/edge/... dofs from
1796 * all admins
1797 */
1798 int n_node_el; /* number of used nodes on each element */
1799 int node[N_NODE_TYPES]; /* index of first vertex/edge/... node*/
1800
1801 unsigned int cookie; /* changed on each refine/coarsen. Use
1802 * this to check consistency of meshes
1803 * and DOF vectors when reading from
1804 * files.
1805 */
1806 int trace_id; /* if this is a trace-mesh (aka sub-mesh)
1807 * then this is a unique id identifying
1808 * it among all the other trace-meshes
1809 * chained to the parent mesh.
1810 */
1811
1812 void *mem_info; /* pointer for administration; don't touch! */
1813 };
1814
1815 /*******************************************************************************
1816 * data structure for basis function representation
1817 ******************************************************************************/
1818
1819 typedef REAL
1820 (*BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
1821 typedef const REAL *
1822 (*GRD_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
1823 typedef const REAL_B *
1824 (*D2_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
1825 typedef const REAL_BB *
1826 (*D3_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
1827 typedef const REAL_BBB *
1828 (*D4_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
1829
1830 typedef const REAL *
1831 (*BAS_FCT_D)(const REAL_B lambda, const BAS_FCTS *thisptr);
1832 typedef const REAL_B *
1833 (*GRD_BAS_FCT_D)(const REAL_B lambda, const BAS_FCTS *thisptr);
1834 typedef const REAL_BB *
1835 (*D2_BAS_FCT_D)(const REAL_B lambda, const BAS_FCTS *thisptr);
1836
1837 typedef void (*REF_INTER_FCT)(DOF_REAL_VEC *, RC_LIST_EL *, int);
1838 typedef void (*REF_INTER_D_FCT)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);
1839 typedef void (*REF_INTER_FCT_D)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);
1840
1841 #define PHI(bfcts, i, lambda) (bfcts)->phi[i](lambda, bfcts)
1842 #define GRD_PHI(bfcts, i, lambda) (bfcts)->grd_phi[i](lambda, bfcts)
1843 #define D2_PHI(bfcts, i, lambda) (bfcts)->D2_phi[i](lambda, bfcts)
1844 #define D3_PHI(bfcts, i, lambda) (bfcts)->D3_phi[i](lambda, bfcts)
1845 #define D4_PHI(bfcts, i, lambda) (bfcts)->D4_phi[i](lambda, bfcts)
1846
1847 #define PHI_D(bfcts, i, lambda) (bfcts)->phi_d[i](lambda, bfcts)
1848 #define GRD_PHI_D(bfcts, i, lambda) (bfcts)->grd_phi_d[i](lambda, bfcts)
1849 #define D2_PHI_D(bfcts, i, lambda) (bfcts)->D2_phi_d[i](lambda, bfcts)
1850
1851 #define GET_DOF_INDICES(bfcts, el, admin, dof) \
1852 (bfcts)->get_dof_indices(dof, el, admin, bfcts)
1853 #define INTERPOL(bfcts, coeff, el_info, wall, n, indices, f, ud) \
1854 (bfcts)->interpol(coeff, el_info, wall, n, indices, f, ud, bfcts)
1855 #define INTERPOL_D(bfcts, coeff, el_info, wall, n, indices, f, ud) \
1856 (bfcts)->interpol_d(coeff, el_info, wall, n, indices, f, ud, bfcts)
1857 #define INTERPOL_DOW(bfcts, coeff, el_info, wall, n, indices, f, ud) \
1858 (bfcts)->interpol_dow(coeff, el_info, wall, n, indices, f, ud, bfcts)
1859 #define GET_BOUND(bfcts, el_info, bound) \
1860 (bfcts)->get_bound(bound, el_info, bfcts)
1861
1862 struct bas_fcts
1863 {
1864 const char *name; /* textual description */
1865 int dim; /* dimension of the corresponding mesh. */
1866 int rdim; /* dimension of the range, 1 or DIM_OF_WORLD */
1867 int n_bas_fcts; /* nu_mber of basisfunctions on one el */
1868 int n_bas_fcts_max;/* max. number in presence of init_element() */
1869 int degree; /* maximal degree of the basis functions,
1870 * may vary on a per-element basis if
1871 * init_element() is != NULL.
1872 */
1873 int n_dof[N_NODE_TYPES]; /* dofs from these bas_fcts */
1874 int trace_admin; /* If >= 0, then the basis function set
1875 * needs a DOF_ADMIN living on a trace
1876 * mesh with id TRACE_ADMIN.
1877 */
1878
1879 /************** link to next set of bfcts in a direct sum *******************/
1880 DBL_LIST_NODE chain;
1881
1882 /* A pointer to the unchained version. It simply points back to the
1883 * same structure if this is an unchained basis-function
1884 * structure.
1885 */
1886 const BAS_FCTS *unchained;
1887
1888 /*************** per-element initializer (maybe NULL) ***********************/
1889
1890 INIT_ELEMENT_DECL;
1891
1892 /*************** the basis functions themselves ***********************/
1893 const BAS_FCT *phi;
1894 const GRD_BAS_FCT *grd_phi;
1895 const D2_BAS_FCT *D2_phi;
1896 const D3_BAS_FCT *D3_phi; /* Optional, implemented for Lagrange bfcts. */
1897 const D4_BAS_FCT *D4_phi; /* Optional, implemented for Lagrange bfcts. */
1898
1899 /* Vector valued basis functions are always factored as phi[i]() *
1900 * phi_d[i](). If phi_d[i]() is piece-wise constant, then
1901 * dir_pw_const should be true. The directions are never cached in
1902 * QUAD_FAST, only the scalar factor.
1903 */
1904 const BAS_FCT_D *phi_d;
1905 const GRD_BAS_FCT_D *grd_phi_d;
1906 const D2_BAS_FCT_D *D2_phi_d;
1907
1908 bool dir_pw_const; /* Direction is p.w. constant on the reference element. */
1909
1910 /*************** the trace space on the wall (e.g. Mini-element) ************/
1911 const BAS_FCTS *trace_bas_fcts; /* The trace space */
1912
1913 /* The local DOF mapping for the trace spaces,
1914 * < 3d:
1915 * [0][0][wall][slave local dof] == master local dof,
1916 * 3d:
1917 * [type > 0][orient < 0][wall][slave local dof] == master local dof.
1918 */
1919 const int *trace_dof_map[2][2][N_WALLS_MAX];
1920
1921 /* This obscure component can vary from wall to wall in the presence
1922 * of an INIT_ELEMENT() method. It is _always_ equal to
1923 * trace_bas_fcts->n_bas_fcts ... BUT ONLY after the respective
1924 * element initializer has been called for trace_bas_fcts on the
1925 * trace mesh. If an INIT_ELEMENT() method is present then it _MUST_
1926 * initialize trace_dof_map _AND_ n_trace_bas_fcts. Of course, in 3D
1927 * only the components corresponding to type and orientation of the
1928 * current EL_INFO object have to be taken care of by the
1929 * INIT_ELEMENT() method.
1930 */
1931 int n_trace_bas_fcts[N_WALLS_MAX];
1932
1933 /*************** interconnection to DOF_ADMIN and mesh ********************/
1934 const EL_DOF_VEC *(*get_dof_indices)(DOF *result,
1935 const EL *, const DOF_ADMIN *,
1936 const BAS_FCTS *thisptr);
1937 const EL_BNDRY_VEC *(*get_bound)(BNDRY_FLAGS *bndry_bits,
1938 const EL_INFO *eli,
1939 const BAS_FCTS *thisptr);
1940
1941 /*************** entries must be set for interpolation ********************/
1942
1943 void (*interpol)(EL_REAL_VEC *coeff,
1944 const EL_INFO *el_info, int wall,
1945 int n, const int *indices,
1946 LOC_FCT_AT_QP f, void *ud,
1947 const BAS_FCTS *thisptr);
1948 void (*interpol_d)(EL_REAL_D_VEC *coeff,
1949 const EL_INFO *el_info, int wall,
1950 int n, const int *indices,
1951 LOC_FCT_D_AT_QP f, void *ud,
1952 const BAS_FCTS *thisptr);
1953 void (*interpol_dow)(EL_REAL_VEC_D *coeff,
1954 const EL_INFO *el_info, int wall,
1955 int n, const int *indices,
1956 LOC_FCT_D_AT_QP f, void *ud,
1957 const BAS_FCTS *thisptr);
1958
1959 /******************** optional entries ***********************************/
1960
1961 const EL_INT_VEC *(*get_int_vec)(int result[],
1962 const EL *, const DOF_INT_VEC *);
1963 const EL_REAL_VEC *(*get_real_vec)(REAL result[],
1964 const EL *, const DOF_REAL_VEC *);
1965 const EL_REAL_D_VEC *(*get_real_d_vec)(REAL_D result[],
1966 const EL *, const DOF_REAL_D_VEC *);
1967 const EL_REAL_VEC_D *(*get_real_vec_d)(REAL result[],
1968 const EL *, const DOF_REAL_VEC_D *);
1969 const EL_UCHAR_VEC *(*get_uchar_vec)(U_CHAR result[],
1970 const EL *, const DOF_UCHAR_VEC *);
1971 const EL_SCHAR_VEC *(*get_schar_vec)(S_CHAR result[],
1972 const EL *, const DOF_SCHAR_VEC *);
1973 const EL_PTR_VEC *(*get_ptr_vec)(void *result[],
1974 const EL *, const DOF_PTR_VEC *);
1975 const EL_REAL_DD_VEC *(*get_real_dd_vec)(REAL_DD result[],
1976 const EL *, const DOF_REAL_DD_VEC *);
1977
1978 void (*real_refine_inter)(DOF_REAL_VEC *, RC_LIST_EL *, int);
1979 void (*real_coarse_inter)(DOF_REAL_VEC *, RC_LIST_EL *, int);
1980 void (*real_coarse_restr)(DOF_REAL_VEC *, RC_LIST_EL *, int);
1981
1982 void (*real_d_refine_inter)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);
1983 void (*real_d_coarse_inter)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);
1984 void (*real_d_coarse_restr)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);
1985
1986 void (*real_refine_inter_d)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);
1987 void (*real_coarse_inter_d)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);
1988 void (*real_coarse_restr_d)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);
1989
1990 void *ext_data; /* Implementation dependent extra data */
1991 };
1992
1993 /* Barycentric coordinates of Lagrange nodes. */
1994 #define LAGRANGE_NODES(bfcts) \
1995 ((const REAL_B *)(*(void **)(bfcts)->ext_data))
1996
1997 /******************************************************************************
1998 * FE spaces are a triple of DOFs and BAS_FCTs on a MESH
1999 *
2000 * Fe-Spaces may be vector-valued, if either the basis functions are
2001 * vector-valued, or on request. rdim codes the dimension of the
2002 * range-space, it may be either 1 or DIM_OF_WORLD, so it is rather a
2003 * flag-value. Scalar fe-space should in principle not be attached to
2004 * DOF_REAL_D_VECs (but can be).
2005 *
2006 * Fe-spaces may be chained to constitute the Cartesian product space
2007 * of a couple of distinct fe-spaces, e.g. to form the Cartesian
2008 * product of some "trivially" vector valued fe-space with another
2009 * space spanned by "really" vector-valued basis functions, e.g. to
2010 * add edge respectively face bubbles.
2011 *
2012 *****************************************************************************/
2013
2014 struct fe_space
2015 {
2016 const char *name;
2017 const DOF_ADMIN *admin;
2018 const BAS_FCTS *bas_fcts;
2019 MESH *mesh;
2020 int rdim;
2021 int ref_cnt;
2022 DBL_LIST_NODE chain;
2023 const FE_SPACE *unchained;
2024 };
2025
2026 /* How to check whether two FE_SPACE objects are essentially the same */
2027 #define FE_SPACE_EQ_P(fe1, fe2) \
2028 ((fe1) == (fe2) || \
2029 ((fe1)->admin == (fe2)->admin && \
2030 (fe1)->bas_fcts == (fe2)->bas_fcts && \
2031 (fe1)->mesh == (fe2)->mesh && \
2032 (fe1)->rdim == (fe2)->rdim))
2033
fe_space_is_eq(const FE_SPACE * fe1,const FE_SPACE * fe2)2034 static inline bool fe_space_is_eq(const FE_SPACE *fe1, const FE_SPACE *fe2)
2035 {
2036 return FE_SPACE_EQ_P(fe1, fe2);
2037 }
2038
2039 /*******************************************************************************
2040 * data structures for numerical integration
2041 ******************************************************************************/
2042
2043 struct quadrature
2044 {
2045 const char *name;
2046 int degree;
2047
2048 int dim; /* barycentric coords have (dim+1) components */
2049 int codim; /* the co-dimension */
2050 int subsplx; /* co-dim 1: face number, co-dim 2: edge number */
2051
2052 int n_points;
2053 int n_points_max; /* max. number in presence of INIT_ELEMENT() */
2054 const REAL_B *lambda;
2055 const REAL *w;
2056
2057 void *metadata; /* for internally kept per element caches etc. */
2058
2059 INIT_ELEMENT_DECL;
2060 };
2061
2062 /*******************************************************************************
2063 * per-element quadrature cache for co-ordinates, det, Lambda, DLambda
2064 * uch a cache structure is associated with each quadrature and can be
2065 * filled by fill_quad_el_cache(el_info, quad, fill_flag). The cache
2066 * is incremental, repeated calls will fill in data as needed.
2067 ******************************************************************************/
2068
2069 struct quad_el_cache
2070 {
2071 EL *current_el;
2072 FLAGS fill_flag;
2073 REAL_D *world;
2074 struct {
2075 REAL *det;
2076 REAL_BD *Lambda;
2077 REAL_BDD *DLambda;
2078 REAL_BD *grd_world;
2079 REAL_BDB *D2_world;
2080 REAL_BDBB *D3_world;
2081 REAL *wall_det; /* for co-dim 1 */
2082 REAL_D *wall_normal; /* for co-dim 1 */
2083 REAL_DB *grd_normal; /* for co-dim 1 */
2084 REAL_DBB *D2_normal; /* for co-dim 1 */
2085 } param;
2086 };
2087
2088 #define FILL_EL_QUAD_WORLD 0x0001
2089 #define FILL_EL_QUAD_DET 0x0002
2090 #define FILL_EL_QUAD_LAMBDA 0x0004
2091 #define FILL_EL_QUAD_DLAMBDA 0x0008
2092 #define FILL_EL_QUAD_GRD_WORLD 0x0010
2093 #define FILL_EL_QUAD_D2_WORLD 0x0020
2094 #define FILL_EL_QUAD_D3_WORLD 0x0040
2095 #define FILL_EL_QUAD_WALL_DET 0x0100
2096 #define FILL_EL_QUAD_WALL_NORMAL 0x0200
2097 #define FILL_EL_QUAD_GRD_NORMAL 0x0400
2098 #define FILL_EL_QUAD_D2_NORMAL 0x0800
2099
2100 static inline const QUAD_EL_CACHE *fill_quad_el_cache(const EL_INFO *el_info,
2101 const QUAD *quad,
2102 FLAGS fill);
2103
2104 /*******************************************************************************
2105 * data structure with precomputed values of basis functions at
2106 * quadrature nodes on the standard element
2107 ******************************************************************************/
2108
2109 #define INIT_PHI 0x01
2110 #define INIT_GRD_PHI 0x02
2111 #define INIT_D2_PHI 0x04
2112 #define INIT_D3_PHI 0x08
2113 #define INIT_D4_PHI 0x10
2114 #define INIT_TANGENTIAL 0x80 /* derivatives are tangential, for co-dim 1
2115 * quadratures.
2116 */
2117
2118 struct quad_fast
2119 {
2120 const QUAD *quad;
2121 const BAS_FCTS *bas_fcts;
2122
2123 FLAGS init_flag;
2124
2125 int dim;
2126 int n_points;
2127 int n_bas_fcts;
2128 int n_points_max;
2129 int n_bas_fcts_max;
2130 const REAL *w; /* shallow copy of quad->w */
2131 const REAL (*const*phi); /* [qp][bf] */
2132 const REAL_B (*const*grd_phi);
2133 const REAL_BB (*const*D2_phi);
2134 const REAL_BBB (*const*D3_phi);
2135 const REAL_BBBB (*const*D4_phi);
2136
2137 /* For vector valued basis functions with a p.w. constant
2138 * directional derivative we cache that direction and make it
2139 * available for applications. The component is initialized by the
2140 * INIT_ELEMENT() method.
2141 *
2142 * So: phi_d[i] gives the value of the directional factor for the
2143 * i-th basis function. If (!bas_fcts->dir_pw_const), then phi_d is
2144 * NULL.
2145 */
2146 const REAL_D *phi_d;
2147
2148 /* chain to next structure, if bas_fcts->chain is non-empty */
2149 DBL_LIST_NODE chain;
2150
2151 /* a clone of this structure, but as single item. */
2152 const QUAD_FAST *unchained;
2153
2154 INIT_ELEMENT_DECL;
2155
2156 void *internal;
2157 };
2158
2159 /*******************************************************************************
2160 * data structure for adaptive methods
2161 ******************************************************************************/
2162
2163 typedef enum adaptation_strategy {
2164 NoStrategy = 0,
2165 GlobalRefinement = 1,
2166 GR = GlobalRefinement,
2167 MinimumStrategy = 2,
2168 MS = MinimumStrategy,
2169 EqualDistributionStrategy = 3,
2170 ES = EqualDistributionStrategy,
2171 GuaranteedErrorReductionStrategy = 4,
2172 GERS = GuaranteedErrorReductionStrategy
2173 } ADAPTATION_STRATEGY;
2174
2175 struct adapt_stat
2176 {
2177 const char *name;
2178 REAL tolerance;
2179 REAL p; /* power in estimator norm */
2180 int max_iteration;
2181 int info;
2182
2183 REAL (*estimate)(MESH *mesh, ADAPT_STAT *adapt);
2184 REAL (*get_el_est)(EL *el); /* local error estimate */
2185 REAL (*get_el_estc)(EL *el); /* local coarsening error estimate*/
2186 U_CHAR (*marking)(MESH *mesh, ADAPT_STAT *adapt);
2187
2188 void *est_info; /* estimator parameters */
2189 REAL err_sum, err_max; /* sum and max of el_est */
2190
2191 void (*build_before_refine)(MESH *mesh, U_CHAR flag);
2192 void (*build_before_coarsen)(MESH *mesh, U_CHAR flag);
2193 void (*build_after_coarsen)(MESH *mesh, U_CHAR flag);
2194 void (*solve)(MESH *mesh);
2195
2196 int refine_bisections;
2197 bool coarsen_allowed; /* 0 : 1 */
2198 int coarse_bisections;
2199 FLAGS adaptation_fill_flags; /* Fill-flags used during adaptation */
2200
2201 ADAPTATION_STRATEGY strategy; /* 1=GR, 2=MS, 3=ES, 4=GERS */
2202 REAL MS_gamma, MS_gamma_c; /* maximum strategy */
2203 REAL ES_theta, ES_theta_c; /* equidistribution strategy */
2204 REAL GERS_theta_star, GERS_nu, GERS_theta_c; /* willy's strategy */
2205 };
2206
2207
2208 struct adapt_instat
2209 {
2210 const char *name;
2211
2212 ADAPT_STAT adapt_initial[1];
2213 ADAPT_STAT adapt_space[1];
2214
2215 REAL time;
2216 REAL start_time, end_time;
2217 REAL timestep;
2218
2219 void (*init_timestep)(MESH *mesh, ADAPT_INSTAT *adapt);
2220 void (*set_time)(MESH *mesh, ADAPT_INSTAT *adapt);
2221 void (*one_timestep)(MESH *mesh, ADAPT_INSTAT *adapt);
2222 REAL (*get_time_est)(MESH *mesh, ADAPT_INSTAT *adapt);
2223 void (*close_timestep)(MESH *mesh, ADAPT_INSTAT *adapt);
2224
2225 int strategy;
2226 int max_iteration;
2227
2228 REAL tolerance;
2229 REAL rel_initial_error;
2230 REAL rel_space_error;
2231 REAL rel_time_error;
2232 REAL time_theta_1;
2233 REAL time_theta_2;
2234 REAL time_delta_1;
2235 REAL time_delta_2;
2236 int info;
2237 };
2238
2239 #define MESH_REFINED 1
2240 #define MESH_COARSENED 2
2241
2242 typedef enum norm
2243 {
2244 NO_NORM = 0, /* uninitialized */
2245 H1_NORM = 1, /* H1-half-norm */
2246 L2_NORM = 2, /* L2-norm */
2247 L2H1_NORM = H1_NORM|L2_NORM /* full H1-norm */
2248 } NORM;
2249
2250 /*******************************************************************************
2251 * data structures for matrix and vector update
2252 ******************************************************************************/
2253
2254 struct el_matrix
2255 {
2256 MATENT_TYPE type;
2257 int n_row, n_col;
2258 int n_row_max, n_col_max;
2259 union {
2260 REAL *const*real;
2261 REAL_D *const*real_d;
2262 REAL_DD *const*real_dd;
2263 } data;
2264 DBL_LIST_NODE row_chain;
2265 DBL_LIST_NODE col_chain;
2266 };
2267
2268 typedef const EL_MATRIX *
2269 (*EL_MATRIX_FCT)(const EL_INFO *el_info, void *fill_info);
2270
2271 typedef struct el_matrix_info EL_MATRIX_INFO;
2272 struct el_matrix_info
2273 {
2274 const FE_SPACE *row_fe_space;
2275 const FE_SPACE *col_fe_space;
2276
2277 MATENT_TYPE krn_blk_type;
2278
2279 BNDRY_FLAGS dirichlet_bndry;
2280 REAL factor;
2281
2282 EL_MATRIX_FCT el_matrix_fct;
2283 void *fill_info;
2284
2285 const EL_MATRIX_FCT *neigh_el_mat_fcts;
2286 void *neigh_fill_info;
2287
2288 FLAGS fill_flag;
2289 };
2290
2291 typedef const EL_REAL_VEC *
2292 (*EL_VEC_FCT)(const EL_INFO *el_info, void *fill_info);
2293
2294 typedef struct el_vec_info EL_VEC_INFO;
2295 struct el_vec_info
2296 {
2297 const FE_SPACE *fe_space;
2298
2299 BNDRY_FLAGS dirichlet_bndry;
2300 REAL factor;
2301
2302 EL_VEC_FCT el_vec_fct;
2303 void *fill_info;
2304
2305 FLAGS fill_flag;
2306 };
2307
2308 typedef const EL_REAL_D_VEC *
2309 (*EL_VEC_D_FCT)(const EL_INFO *el_info, void *fill_info);
2310
2311 typedef struct el_vec_d_info EL_VEC_D_INFO;
2312 struct el_vec_d_info
2313 {
2314 const FE_SPACE *fe_space;
2315
2316 BNDRY_FLAGS dirichlet_bndry;
2317 REAL factor;
2318
2319 EL_VEC_D_FCT el_vec_fct;
2320 void *fill_info;
2321
2322 FLAGS fill_flag;
2323 };
2324
2325 typedef const EL_REAL_VEC_D *
2326 (*EL_VEC_FCT_D)(const EL_INFO *el_info, void *fill_info);
2327
2328 typedef struct el_vec_info_d EL_VEC_INFO_D;
2329 struct el_vec_info_d
2330 {
2331 const FE_SPACE *fe_space;
2332
2333 BNDRY_FLAGS dirichlet_bndry;
2334 REAL factor;
2335
2336 EL_VEC_FCT_D el_vec_fct;
2337 void *fill_info;
2338
2339 FLAGS fill_flag;
2340 };
2341
2342 /*******************************************************************************
2343 * a matrix of quadratures to use for each block of a block-operator
2344 * acting on a direct sum of FE-space.
2345 ******************************************************************************/
2346 typedef struct quad_tensor
2347 {
2348 const QUAD *quad; /* the quadrature rules to use for this block */
2349 DBL_LIST_NODE row_chain; /* cyclic list link for one row */
2350 DBL_LIST_NODE col_chain; /* cyclic list link for one column */
2351 DBL_LIST_NODE dep_chain; /* cyclic list link for the third index */
2352 } QUAD_TENSOR;
2353
2354 typedef struct wall_quad_tensor
2355 {
2356 const WALL_QUAD *quad; /* the quadrature rule to use for this block */
2357 DBL_LIST_NODE row_chain; /* cyclic list link for one row */
2358 DBL_LIST_NODE col_chain; /* cyclic list link for on e column */
2359 DBL_LIST_NODE dep_chain; /* cyclic list link for the third index */
2360 } WALL_QUAD_TENSOR;
2361
2362 /*******************************************************************************
2363 * data structure about the differential operator for matrix assemblage
2364 ******************************************************************************/
2365
2366 typedef struct operator_info OPERATOR_INFO;
2367 struct operator_info
2368 {
2369 const FE_SPACE *row_fe_space; /* range fe-space */
2370 const FE_SPACE *col_fe_space; /* domain fe-space */
2371
2372 const QUAD *quad[3];
2373 const QUAD_TENSOR *quad_tensor[3];
2374
2375 bool (*init_element)(const EL_INFO *el_info, const QUAD *quad[3], void *apd);
2376
2377 union {
2378 const REAL_B *(*real)(const EL_INFO *el_info,
2379 const QUAD *quad, int iq, void *apd);
2380 const REAL_BD *(*real_d)(const EL_INFO *el_info,
2381 const QUAD *quad, int iq, void *apd);
2382 const REAL_BDD *(*real_dd)(const EL_INFO *el_info,
2383 const QUAD *quad, int iq, void *apd);
2384 } LALt;
2385 MATENT_TYPE LALt_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
2386 bool LALt_pw_const;
2387 bool LALt_symmetric;
2388 int LALt_degree; /* quadrature degree of the LALt() kernel */
2389
2390 union {
2391 const REAL *(*real)(const EL_INFO *el_info,
2392 const QUAD *quad, int iq, void *apd);
2393 const REAL_D *(*real_d)(const EL_INFO *el_info,
2394 const QUAD *quad, int iq, void *apd);
2395 const REAL_DD *(*real_dd)(const EL_INFO *el_info,
2396 const QUAD *quad, int iq, void *apd);
2397 const REAL_DDD *(*real_ddd)(const EL_INFO *el_info,
2398 const QUAD *quad, int iq, void *apd);
2399 } Lb0;
2400 bool Lb0_pw_const;
2401 union {
2402 const REAL *(*real)(const EL_INFO *el_info,
2403 const QUAD *quad, int iq, void *apd);
2404 const REAL_D *(*real_d)(const EL_INFO *el_info,
2405 const QUAD *quad, int iq, void *apd);
2406 const REAL_DD *(*real_dd)(const EL_INFO *el_info,
2407 const QUAD *quad, int iq, void *apd);
2408 const REAL_DDD *(*real_ddd)(const EL_INFO *el_info,
2409 const QUAD *quad, int iq, void *apd);
2410 } Lb1;
2411 bool Lb1_pw_const;
2412 MATENT_TYPE Lb_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
2413 bool Lb0_Lb1_anti_symmetric;
2414 int Lb_degree; /* quadrature degree for the Lb0() & Lb1() kernel */
2415 const EL_REAL_VEC_D *(*advection_field)(const EL_INFO *el_info, void *apd);
2416 const FE_SPACE *adv_fe_space;
2417
2418 union {
2419 REAL (*real)(const EL_INFO *el_info,
2420 const QUAD *quad, int iq, void *apd);
2421 const REAL *(*real_d)(const EL_INFO *el_info,
2422 const QUAD *quad, int iq, void *apd);
2423 const REAL_D *(*real_dd)(const EL_INFO *el_info,
2424 const QUAD *quad, int iq, void *apd);
2425 } c;
2426 bool c_pw_const;
2427 MATENT_TYPE c_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
2428 int c_degree; /* quadrature degree for the c()-kernel */
2429
2430 BNDRY_FLAGS dirichlet_bndry; /* bndry-type bit-mask for
2431 * Dirichlet-boundary conditions
2432 * built into the matrix
2433 */
2434 FLAGS fill_flag;
2435 void *user_data; /* application data, passed to init_element */
2436 };
2437
2438 typedef struct bndry_operator_info BNDRY_OPERATOR_INFO;
2439
2440 struct bndry_operator_info
2441 {
2442 const FE_SPACE *row_fe_space;
2443 const FE_SPACE *col_fe_space;
2444
2445 const WALL_QUAD *quad[3];
2446 const WALL_QUAD_TENSOR *quad_tensor[3];
2447
2448 bool (*init_element)(const EL_INFO *el_info, int wall,
2449 const WALL_QUAD *quad[3], void *ud);
2450
2451 union {
2452 const REAL_B *(*real)(const EL_INFO *el_info,
2453 const QUAD *quad, int iq, void *apd);
2454 const REAL_BD *(*real_d)(const EL_INFO *el_info,
2455 const QUAD *quad, int iq, void *apd);
2456 const REAL_BDD *(*real_dd)(const EL_INFO *el_info,
2457 const QUAD *quad, int iq, void *apd);
2458 } LALt;
2459 MATENT_TYPE LALt_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
2460 bool LALt_pw_const;
2461 bool LALt_symmetric;
2462 int LALt_degree; /* quad-deg for the LALt() kernel */
2463
2464 union {
2465 const REAL *(*real)(const EL_INFO *el_info,
2466 const QUAD *quad, int iq, void *apd);
2467 const REAL_D *(*real_d)(const EL_INFO *el_info,
2468 const QUAD *quad, int iq, void *apd);
2469 const REAL_DD *(*real_dd)(const EL_INFO *el_info,
2470 const QUAD *quad, int iq, void *apd);
2471 const REAL_DDD *(*real_ddd)(const EL_INFO *el_info,
2472 const QUAD *quad, int iq, void *apd);
2473 } Lb0;
2474 bool Lb0_pw_const;
2475 union {
2476 const REAL *(*real)(const EL_INFO *el_info,
2477 const QUAD *quad, int iq, void *apd);
2478 const REAL_D *(*real_d)(const EL_INFO *el_info,
2479 const QUAD *quad, int iq, void *apd);
2480 const REAL_DD *(*real_dd)(const EL_INFO *el_info,
2481 const QUAD *quad, int iq, void *apd);
2482 const REAL_DDD *(*real_ddd)(const EL_INFO *el_info,
2483 const QUAD *quad, int iq, void *apd);
2484 } Lb1;
2485 bool Lb1_pw_const;
2486 MATENT_TYPE Lb_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
2487 bool Lb0_Lb1_anti_symmetric;
2488 int Lb_degree; /* quad-deg for the Lb0() and Lb1() kernel */
2489
2490 /* The following two entries are not used yet. */
2491 const EL_REAL_VEC_D *(*advection_field)(const EL_INFO *el_info, void *apd);
2492 const FE_SPACE *adv_fe_space;
2493
2494 union {
2495 REAL (*real)(const EL_INFO *el_info,
2496 const QUAD *quad, int iq, void *apd);
2497 const REAL *(*real_d)(const EL_INFO *el_info,
2498 const QUAD *quad, int iq, void *apd);
2499 const REAL_D *(*real_dd)(const EL_INFO *el_info,
2500 const QUAD *quad, int iq, void *apd);
2501 } c;
2502 bool c_pw_const;
2503 MATENT_TYPE c_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
2504 int c_degree; /* quad-deg of the c() kernel */
2505
2506 /* boundary segment(s) we belong to; if
2507 * BNDRY_FLAGS_IS_INTERIOR(bndry_type), then the operator is invoked
2508 * on all interior faces, e.g. to implement a DG-method.
2509 */
2510 BNDRY_FLAGS bndry_type;
2511
2512 bool discontinuous; /* assemble jumps w.r.t. the neighbour */
2513 bool tangential; /* use tangential gradients */
2514
2515 FLAGS fill_flag;
2516 void *user_data;
2517 };
2518
2519 /*******************************************************************************
2520 *
2521 * A structure describing the "closure" of an differential operator by
2522 * Dirichlet, Neumann or Robin boundary conditions.
2523 *
2524 * Note that the quadrature passed to neumann() and robin() is a
2525 * co-dim 1 quadrature, so quad->subsplx is the number of the boundary
2526 * face, and el_info->wall_bound[quad->subsplx] contains the boundary
2527 * classification.
2528 *
2529 ******************************************************************************/
2530
2531 typedef struct bndry_cond_info BNDRY_COND_INFO;
2532 struct bndry_cond_info
2533 {
2534 const FE_SPACE *fe_space;
2535
2536 REAL (*dirichlet)(const EL_INFO *el_info,
2537 const REAL_B lambda,
2538 BNDRY_FLAGS bndry,
2539 void *apd);
2540 BNDRY_FLAGS dirichlet_bndry;
2541
2542 REAL (*neumann)(const EL_INFO *el_info, const QUAD *quad, int iq, void *apd);
2543 const WALL_QUAD *neumann_quad;
2544 BNDRY_FLAGS neumann_bndry;
2545
2546 REAL (*robin)(const EL_INFO *el_info, const QUAD *quad, int iq, void *apd);
2547 const WALL_QUAD *robin_quad;
2548 REAL robin_const; /* if robin == NULL, use this constant value */
2549 BNDRY_FLAGS robin_bndry;
2550
2551 S_CHAR (*bndry_type)(const BNDRY_FLAGS bndry_bits);
2552
2553 void *user_data;
2554 };
2555
2556 typedef struct bndry_cond_info_d BNDRY_COND_INFO_D;
2557 struct bndry_cond_info_d
2558 {
2559 const FE_SPACE *fe_space;
2560
2561 const REAL *(*dirichlet)(REAL_D result,
2562 const EL_INFO *el_info,
2563 const REAL_B lambda,
2564 void *apd);
2565 BNDRY_FLAGS dirichlet_bndry;
2566
2567 const REAL *(*neumann)(REAL_D result,
2568 const EL_INFO *el_info,
2569 const QUAD *quad, int iq, void *apd);
2570 const WALL_QUAD *neumann_quad;
2571 BNDRY_FLAGS neumann_bound;
2572
2573 const REAL *(*robin)(REAL_D result,
2574 const EL_INFO *el_info,
2575 const QUAD *quad, int iq, void *apd);
2576 const WALL_QUAD *robin_quad;
2577 REAL_D robin_const; /* if robin == NULL, use this constant value */
2578 BNDRY_FLAGS robin_bound;
2579
2580 S_CHAR (*bndry_type)(const BNDRY_FLAGS bndry_bits);
2581
2582 void *user_data;
2583 };
2584
2585 /*******************************************************************************
2586 * calculate element stiffness matrices by preevaluated integrals
2587 * over the the reference element. In this notation, "PSI" means the
2588 * row-space, i.e. the space of test-functions, "PHI" means the
2589 * column-space, i.e. the space of ansatz functions.
2590 *
2591 * The special tri-linear caches Q001, Q010 and Q100 are meant to
2592 * support assembling tri-linear forms, e.g. for material
2593 * derivatives. The caches are present, but currently there is not
2594 * further support for tri-linear forms inside ALBERTA. The "1"
2595 * denotes which of the three factors is differentiated.
2596 ******************************************************************************/
2597
2598 typedef struct q11_psi_phi Q11_PSI_PHI;
2599 typedef struct q01_psi_phi Q01_PSI_PHI;
2600 typedef struct q10_psi_phi Q10_PSI_PHI;
2601 typedef struct q00_psi_phi Q00_PSI_PHI;
2602 typedef struct q001_eta_psi_phi Q001_ETA_PSI_PHI;
2603 typedef struct q010_eta_psi_phi Q010_ETA_PSI_PHI;
2604 typedef struct q100_eta_psi_phi Q100_ETA_PSI_PHI;
2605
2606 /* for i = 0 ... n_psi-1
2607 * for j = 0 ... n_phi-1
2608 * for m = 0 ... n_entries[i][j]
2609 *
2610 * Then we have: values[m][i][j] is the product of
2611 * (\nabla\psi_i)_{k[i][j][m]} with (\nabla\phi_j)_{l[i][j][m]}
2612 *
2613 * end for
2614 * end for
2615 * end for
2616 *
2617 * Products yielding zero are left out.
2618 */
2619 typedef struct q11_psi_phi_cache
2620 {
2621 int n_psi;
2622 int n_phi;
2623
2624 const int *const*n_entries;
2625 const REAL *const*const*values;
2626 const int *const*const*k;
2627 const int *const*const*l;
2628 } Q11_PSI_PHI_CACHE;
2629
2630 struct q11_psi_phi
2631 {
2632 const BAS_FCTS *psi;
2633 const BAS_FCTS *phi;
2634 const QUAD *quad;
2635
2636 const Q11_PSI_PHI_CACHE *cache;
2637
2638 INIT_ELEMENT_DECL;
2639 };
2640
2641 /* for i = 0 ... n_psi-1
2642 * for j = 0 ... n_phi-1
2643 * for m = 0 ... n_entries[i][j]
2644 *
2645 * Then we have: values[m][i][j] is the product of
2646 * \psi_i with (\nabla\phi_j)_{l[i][j][m]}
2647 *
2648 * end for
2649 * end for
2650 * end for
2651 *
2652 * Products yielding zero are left out.
2653 */
2654 typedef struct q01_psi_phi_cache
2655 {
2656 int n_psi;
2657 int n_phi;
2658
2659 const int *const*n_entries;
2660 const REAL *const*const*values;
2661 const int *const*const*l;
2662 } Q01_PSI_PHI_CACHE;
2663
2664 struct q01_psi_phi
2665 {
2666 const BAS_FCTS *psi;
2667 const BAS_FCTS *phi;
2668 const QUAD *quad;
2669
2670 const Q01_PSI_PHI_CACHE *cache;
2671
2672 INIT_ELEMENT_DECL;
2673 };
2674
2675 /* for i = 0 ... n_psi-1
2676 * for j = 0 ... n_phi-1
2677 * for m = 0 ... n_entries[i][j]
2678 *
2679 * Then we have: values[m][i][j] is the product of
2680 * (\nabla\psi_i)_{k[i][j][m]} with \phi_j
2681 *
2682 * end for
2683 * end for
2684 * end for
2685 *
2686 * Products yielding zero are left out.
2687 */
2688 typedef struct q10_psi_phi_cache
2689 {
2690 int n_psi;
2691 int n_phi;
2692
2693 const int *const*n_entries;
2694 const REAL *const*const*values;
2695 const int *const*const*k;
2696 } Q10_PSI_PHI_CACHE;
2697
2698 struct q10_psi_phi
2699 {
2700 const BAS_FCTS *psi;
2701 const BAS_FCTS *phi;
2702 const QUAD *quad;
2703
2704 const Q10_PSI_PHI_CACHE *cache;
2705
2706 INIT_ELEMENT_DECL;
2707 };
2708
2709 typedef struct q00_psi_phi_cache
2710 {
2711 int n_psi;
2712 int n_phi;
2713
2714 const REAL *const*values;
2715 } Q00_PSI_PHI_CACHE;
2716
2717 struct q00_psi_phi
2718 {
2719 const BAS_FCTS *psi;
2720 const BAS_FCTS *phi;
2721 const QUAD *quad;
2722
2723 const Q00_PSI_PHI_CACHE *cache;
2724
2725 INIT_ELEMENT_DECL;
2726 };
2727
2728 /* for i = 0 ... n_psi-1
2729 * for j = 0 ... n_phi-1
2730 * for k = 0 ... n_eta-1
2731 * for m = 0 ... n_entries[i][j][k]
2732 *
2733 * Then we have: values[m][i][j][k] is the product of
2734 * \psi_i, \phi_j with (\nabla\eta_k)_{l[i][j][k][m]}
2735 * end for
2736 * end for
2737 * end for
2738 * end for
2739 *
2740 * Products yielding zero are left out. The analogue holds for the 010
2741 * and 100 structures, where the position of the "1" marks which
2742 * factor is differentitated.
2743 */
2744 typedef struct q001_eta_psi_phi_cache
2745 {
2746 int n_eta;
2747 int n_psi;
2748 int n_phi;
2749
2750 const int *const*const*n_entries;
2751 const REAL *const*const*const*values;
2752 const int *const*const*const*l;
2753 } Q001_ETA_PSI_PHI_CACHE;
2754
2755 struct q001_eta_psi_phi
2756 {
2757 const BAS_FCTS *eta;
2758 const BAS_FCTS *psi;
2759 const BAS_FCTS *phi;
2760 const QUAD *quad;
2761
2762 const Q001_ETA_PSI_PHI_CACHE *cache;
2763
2764 INIT_ELEMENT_DECL;
2765 };
2766
2767 typedef struct q010_eta_psi_phi_cache
2768 {
2769 int n_eta;
2770 int n_psi;
2771 int n_phi;
2772
2773 const int *const*const*n_entries;
2774 const REAL *const*const*const*values;
2775 const int *const*const*const*l;
2776 } Q010_ETA_PSI_PHI_CACHE;
2777
2778 struct q010_eta_psi_phi
2779 {
2780 const BAS_FCTS *eta;
2781 const BAS_FCTS *psi;
2782 const BAS_FCTS *phi;
2783 const QUAD *quad;
2784
2785 const Q010_ETA_PSI_PHI_CACHE *cache;
2786
2787 INIT_ELEMENT_DECL;
2788 };
2789
2790 typedef struct q100_eta_psi_phi_cache
2791 {
2792 int n_eta;
2793 int n_psi;
2794 int n_phi;
2795
2796 const int *const*const*n_entries;
2797 const REAL *const*const*const*values;
2798 const int *const*const*const*l;
2799 } Q100_ETA_PSI_PHI_CACHE;
2800
2801 struct q100_eta_psi_phi
2802 {
2803 const BAS_FCTS *eta;
2804 const BAS_FCTS *psi;
2805 const BAS_FCTS *phi;
2806 const QUAD *quad;
2807
2808 const Q100_ETA_PSI_PHI_CACHE *cache;
2809
2810 INIT_ELEMENT_DECL;
2811 };
2812
2813 /******************************************************************************/
2814
2815 /* Data for assembling a theta splitting scheme. */
2816 typedef struct el_sys_info_instat EL_SYS_INFO_INSTAT;
2817 struct el_sys_info_instat
2818 {
2819 const FE_SPACE *row_fe_space;
2820 const FE_SPACE *col_fe_space;
2821
2822 INIT_EL_TAG (*el_update_fct)(const EL_INFO *el_info,
2823 REAL tau, REAL theta,
2824 EL_SYS_INFO_INSTAT *thisptr);
2825 EL_MATRIX *el_matrix;
2826 EL_REAL_VEC *el_load;
2827
2828 const EL_MATRIX *el_stiff;
2829 const EL_MATRIX *el_mass;
2830 const EL_REAL_VEC *u_h_loc;
2831
2832 FLAGS fill_flag;
2833 BNDRY_FLAGS dirichlet_bndry; /* bndry-type bit-mask for
2834 * Dirichlet-boundary conditions
2835 * built into the matrix
2836 */
2837 MATENT_TYPE krn_blk_type; /* MATENT_REAL for scalar problems */
2838 };
2839
2840 typedef struct el_sys_info_dow_instat EL_SYS_INFO_DOW_INSTAT;
2841 struct el_sys_info_dow_instat
2842 {
2843 const FE_SPACE *row_fe_space;
2844 const FE_SPACE *col_fe_space;
2845
2846 INIT_EL_TAG (*el_update_fct)(const EL_INFO *el_info,
2847 REAL tau, REAL theta,
2848 EL_SYS_INFO_DOW_INSTAT *thisptr);
2849 EL_MATRIX *el_matrix;
2850 EL_REAL_VEC_D *el_load;
2851
2852 const EL_MATRIX *el_stiff;
2853 const EL_MATRIX *el_mass;
2854 const EL_REAL_VEC_D *u_h_loc;
2855
2856 FLAGS fill_flag;
2857 BNDRY_FLAGS dirichlet_bndry; /* bndry-type bit-mask for
2858 * Dirichlet-boundary conditions
2859 * built into the matrix
2860 */
2861 MATENT_TYPE krn_blk_type;
2862 };
2863
2864 typedef struct el_sys_info_d_instat EL_SYS_INFO_D_INSTAT;
2865 struct el_sys_info_d_instat
2866 {
2867 const FE_SPACE *row_fe_space;
2868 const FE_SPACE *col_fe_space;
2869
2870 INIT_EL_TAG (*el_update_fct)(const EL_INFO *el_info,
2871 REAL tau, REAL theta,
2872 EL_SYS_INFO_D_INSTAT *thisptr);
2873 EL_MATRIX *el_matrix;
2874 EL_REAL_D_VEC *el_load;
2875
2876 const EL_MATRIX *el_stiff;
2877 const EL_MATRIX *el_mass;
2878 const EL_REAL_D_VEC *u_h_loc;
2879
2880 FLAGS fill_flag;
2881 BNDRY_FLAGS dirichlet_bndry; /* bndry-type bit-mask for
2882 * Dirichlet-boundary conditions
2883 * built into the matrix
2884 */
2885 MATENT_TYPE krn_blk_type;
2886 };
2887
2888
2889
2890 /*******************************************************************************
2891 * preconditioner types.
2892 ******************************************************************************/
2893 typedef enum {
2894 PreconEnd = -1, /* Terminator for variable argument precon functions */
2895 PreconRepeat = PreconEnd,
2896 NoPrecon = 0,
2897 DiagPrecon = 1,
2898 HBPrecon = 2,
2899 BPXPrecon = 3,
2900 SSORPrecon = 4, /* omega == 1, n_iter = 1 */
2901 __SSORPrecon = 5, /* SSOR, but with variable omega and n_iter */
2902 ILUkPrecon = 6, /* combinatorical ILU(k) */
2903 BlkDiagPrecon = 512,
2904 BlkSSORPrecon = 513,
2905 } OEM_PRECON;
2906
2907 #define N_BLOCK_PRECON_MAX 10
2908
2909 /*******************************************************************************
2910 * A data structure which can be use to define more complex
2911 * preconditioners. The purpose of this structure is to avoid defining
2912 * functions with an endless number of arguments. This
2913 * "parameter-transport-structure" can be passed to
2914 * init_precon_from_type(), instead of calling init_oem_precon().
2915 *
2916 * The general idea is:
2917 *
2918 * type -- one of the preconditoner types defined above.
2919 *
2920 * param -- if the preconditioner defined by "type" needs additional
2921 * parameters, then the corresponding section in the "param"
2922 * component has to be filled.
2923 *
2924 * Examples (using a C99-compliant C-compiler):
2925 *
2926 * type == __SSORPrecon:
2927 *
2928 * PRECON_TYPE prec = {
2929 * __SSORPrecon,
2930 * { .__SSOR = { 1.5, 2 } }
2931 * };
2932 *
2933 * type == BlkDiagPrecon, the FE-space is a direct sum of 3
2934 * components, e.g. the Crouzeix-Raviart-Mansfield Stokes
2935 * discretisation in 3d (Lagrange-2 + face-bubble + wall-bubbles):
2936 *
2937 * PRECON_TYPE = {
2938 * BlkDiagPrecon,
2939 * { .BlkDiag = {
2940 * { __SSORPrecon, { 1.0, 1 } },
2941 * { DiagPrecon },
2942 * { DIagPrecon }, }
2943 * }
2944 * }
2945 *
2946 *
2947 ******************************************************************************/
2948
2949 /* Helper struct for BlkPrecon parameters */
2950 struct __precon_type {
2951 OEM_PRECON type;
2952 union {
2953 struct {
2954 REAL omega;
2955 int n_iter;
2956 } __SSOR;
2957 struct {
2958 int level;
2959 } ILUk;
2960 } param;
2961 };
2962
2963 typedef struct precon_type
2964 {
2965 OEM_PRECON type;
2966 union {
2967 struct {
2968 REAL omega;
2969 int n_iter;
2970 } __SSOR;
2971 struct {
2972 int level;
2973 } ILUk;
2974 struct {
2975 struct __precon_type precon[N_BLOCK_PRECON_MAX];
2976 } BlkDiag;
2977 struct {
2978 struct __precon_type precon[N_BLOCK_PRECON_MAX];
2979 REAL omega;
2980 int n_iter;
2981 } BlkSSOR;
2982 } param;
2983 } PRECON_TYPE;
2984
2985 /*******************************************************************************
2986 * The precon "class".
2987 *
2988 * precon_data: Opaque data pointer.
2989 *
2990 * init_precon(): Has to be called after the operator has changed,
2991 * e.g. to update the inverse of the diagonal after the
2992 * matrix has changed.
2993 *
2994 * exit_precon(): Destroys the preconditioner, include the precon
2995 * structure itself.
2996 *
2997 * precon(): The preconditioner itself.
2998 *
2999 ******************************************************************************/
3000
3001 typedef struct precon PRECON;
3002 struct precon
3003 {
3004 void *precon_data;
3005
3006 bool (*init_precon)(void *precon_data);
3007 void (*precon)(void *precon_data, int n, REAL *vec);
3008 void (*exit_precon)(void *precon_data);
3009 };
3010
3011 extern const PRECON *get_diag_precon(const DOF_MATRIX *A,
3012 const DOF_SCHAR_VEC *bound);
3013 extern const PRECON *get_HB_precon(const DOF_MATRIX *matrix,
3014 const DOF_SCHAR_VEC *bound,
3015 int info);
3016 extern const PRECON *get_BPX_precon(const DOF_MATRIX *matrix,
3017 const DOF_SCHAR_VEC *bound,
3018 int info);
3019 extern const PRECON *get_SSOR_precon(const DOF_MATRIX *A,
3020 const DOF_SCHAR_VEC *bound,
3021 REAL omega,
3022 int n_iter);
3023 extern const PRECON *get_ILUk_precon(const DOF_MATRIX *A,
3024 const DOF_SCHAR_VEC *mask,
3025 int ilu_level, int info);
3026
3027 /*******************************************************************************
3028 * abstract multigrid
3029 ******************************************************************************/
3030
3031 typedef struct multi_grid_info MULTI_GRID_INFO;
3032
3033 struct multi_grid_info
3034 {
3035 REAL tolerance; /* tol. for resid */
3036 REAL exact_tolerance; /* tol. for exact_solver */
3037
3038 int cycle; /* 1=V-cycle, 2=W-cycle */
3039 int n_pre_smooth, n_in_smooth; /* no of smoothing loops */
3040 int n_post_smooth; /* no of smoothing loops */
3041 int mg_levels; /* current no. of levels */
3042 int exact_level; /* level for exact_solver */
3043 int max_iter; /* max. no of MG iter's */
3044 int info;
3045
3046 int (*init_multi_grid)(MULTI_GRID_INFO *mg_info);
3047 void (*pre_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
3048 void (*in_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
3049 void (*post_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
3050 void (*mg_restrict)(MULTI_GRID_INFO *mg_info, int level);
3051 void (*mg_prolongate)(MULTI_GRID_INFO *mg_info, int level);
3052 void (*exact_solver)(MULTI_GRID_INFO *mg_info, int level);
3053 REAL (*mg_resid)(MULTI_GRID_INFO *mg_info, int level);
3054 void (*exit_multi_grid)(MULTI_GRID_INFO *mg_info);
3055
3056 void *data; /* application dep. data */
3057
3058 };
3059 int MG(MULTI_GRID_INFO *mg_info);
3060
3061 /*******************************************************************************
3062 * concrete multigrid
3063 ******************************************************************************/
3064
3065 typedef struct mg_s_info MG_S_INFO;
3066 struct mg_s_info
3067 {
3068 MULTI_GRID_INFO *mg_info; /* abstract MG info */
3069
3070 const FE_SPACE *fe_space;
3071 const DOF_ADMIN *vertex_admin;
3072 DOF_MATRIX *mat;
3073 const DOF_REAL_VEC *f;
3074 DOF_REAL_VEC *u;
3075 const DOF_SCHAR_VEC *bound;
3076
3077 int smoother, exact_solver;
3078 REAL smooth_omega, exact_omega;
3079
3080 int size; /* current size of vectors*/
3081 DOF_MATRIX **matrix; /* one for each level */
3082 REAL **f_h; /* one for each level */
3083 REAL **u_h; /* one for each level */
3084 REAL **r_h; /* one for each level */
3085 int *dofs_per_level; /* count dofs per level */
3086
3087 int sort_size; /* size of sort vectors */
3088 DOF *sort_dof; /* dofs in order of levels*/
3089 DOF *(dof_parent[2]); /* (for linear elements) */
3090 U_CHAR *dof_level;
3091 S_CHAR *sort_bound; /* sorted bound */
3092
3093 int sort_invers_size; /* size of inv. sort list */
3094 int *sort_dof_invers; /* inverse sort list */
3095 };
3096 /*******************************************************************************
3097 * sort_dof[ sorted dof ] = unsorted dof
3098 * sort_dof_invers[ unsorted dof ] = sorted dof
3099 ******************************************************************************/
3100
3101 /* file MG_s1.c DOF_sort routines *********************************************/
3102 void MG_s_setup_levels(MG_S_INFO *mg_s_info);
3103 void MG_s_setup_mat_b(MG_S_INFO *mg_s_info,
3104 DOF_MATRIX *mat, const DOF_SCHAR_VEC *bound);
3105 void MG_s_dof_copy_to_sparse(MG_S_INFO *mg_s_info,
3106 const DOF_REAL_VEC *x, REAL *y);
3107 void MG_s_dof_copy_from_sparse(MG_S_INFO *mg_s_info,
3108 const REAL *x, DOF_REAL_VEC *y);
3109 void MG_s_reset_mat(MG_S_INFO *mg_s_info);
3110 void MG_s_sort_mat(MG_S_INFO *mg_s_info);
3111 void MG_s_free_mem(MG_S_INFO *mg_s_info);
3112
3113 /* file MG_s2.c: DOF_sort independent routines ********************************/
3114 void MG_s_restrict_mg_matrices(MG_S_INFO *mg_s_info);
3115 void MG_s_restrict(MULTI_GRID_INFO *mg_info, int mg_level);
3116 void MG_s_prolongate(MULTI_GRID_INFO *mg_info, int mg_level);
3117 REAL MG_s_resid(MULTI_GRID_INFO *mg_info, int mg_level);
3118 void MG_s_smoother(MULTI_GRID_INFO *mg_info, int mg_level, int n);
3119 void MG_s_exact_solver(MULTI_GRID_INFO *mg_info, int mg_level);
3120 void MG_s_gemv(MG_S_INFO *mg_s_info, int mg_level, MatrixTranspose transpose,
3121 REAL alpha, DOF_MATRIX *a, REAL *x, REAL beta, REAL *y);
3122
3123
3124 /* file MG_s.c: ***************************************************************/
3125 int mg_s(DOF_MATRIX *matrix, DOF_REAL_VEC *u, const DOF_REAL_VEC *f,
3126 const DOF_SCHAR_VEC *bound,
3127 REAL tol, int max_iter, int info, char *prefix);
3128 MG_S_INFO *mg_s_init(DOF_MATRIX *matrix, const DOF_SCHAR_VEC *bound,
3129 int info, char *prefix);
3130 int mg_s_solve(MG_S_INFO *mg_s_info,
3131 DOF_REAL_VEC *u, const DOF_REAL_VEC *f, REAL tol, int max_iter);
3132 void mg_s_exit(MG_S_INFO *mg_s_info);
3133
3134 /*******************************************************************************
3135 * Graphic output Definitions
3136 ******************************************************************************/
3137
3138 typedef void * GRAPH_WINDOW;
3139 typedef float GRAPH_RGBCOLOR[3];
3140
3141 /** flags used by graph_mesh(): ****/
3142 #define GRAPH_MESH_BOUNDARY 1
3143 #define GRAPH_MESH_ELEMENT_MARK 2
3144 #define GRAPH_MESH_VERTEX_DOF 4
3145 #define GRAPH_MESH_ELEMENT_INDEX 8
3146
3147 /*******************************************************************************
3148 * very useful macro definitons
3149 ******************************************************************************/
3150
3151 #define GET_MESH(dim, name, macro_data, init_node_proj, init_wall_trafo) \
3152 check_and_get_mesh((dim), DIM_OF_WORLD, ALBERTA_DEBUG, \
3153 ALBERTA_VERSION, (name), (macro_data), \
3154 (init_node_proj), (init_wall_trafo))
3155
3156 #define GET_DOF_VEC(ptr, dof_vec) \
3157 { \
3158 DEBUG_TEST_EXIT((dof_vec) && (dof_vec)->vec, \
3159 "%s == NULL\n", (dof_vec) ? NAME(dof_vec) : #dof_vec); \
3160 (ptr) = (dof_vec)->vec; \
3161 }
3162
3163 /*******************************************************************************
3164 * defined in graphXO.c
3165 ******************************************************************************/
3166 extern const GRAPH_RGBCOLOR rgb_black;
3167 extern const GRAPH_RGBCOLOR rgb_white;
3168 extern const GRAPH_RGBCOLOR rgb_red;
3169 extern const GRAPH_RGBCOLOR rgb_green;
3170 extern const GRAPH_RGBCOLOR rgb_blue;
3171 extern const GRAPH_RGBCOLOR rgb_yellow;
3172 extern const GRAPH_RGBCOLOR rgb_magenta;
3173 extern const GRAPH_RGBCOLOR rgb_cyan;
3174 extern const GRAPH_RGBCOLOR rgb_grey50;
3175
3176 extern const GRAPH_RGBCOLOR rgb_albert;
3177 extern const GRAPH_RGBCOLOR rgb_alberta;
3178
3179
3180 /*******************************************************************************
3181 * used in wall_quad
3182 ******************************************************************************/
3183
3184 /* A collection of quadrature rules for the integration over walls
3185 * (3d: faces, 2d: edges) of a simplex.
3186 *
3187 * Each of the quadrature rules WALL_QUAD::quad[wall] may have its
3188 * own INIT_ELEMENT method.
3189 *
3190 * INIT_ELEMENT(el_info, WALL_QUAD) may or may not be called: it is
3191 * legal to only call INIT_ELEMENT(el_info, WALL_QUAD::quad[wall)
3192 * individually.
3193 *
3194 * If INIT_ELEMENT(el_info, WALL_QUAD) is called, then it has to
3195 * initialize all quadrature rules for all walls, so the sub-ordinate
3196 * initializers need not be called in this case.
3197 */
3198 struct wall_quadrature
3199 {
3200 const char *name;
3201 int degree;
3202 int dim;
3203 int n_points_max;
3204 QUAD quad[N_WALLS_MAX];
3205
3206 INIT_ELEMENT_DECL;
3207
3208 void *metadata;
3209 };
3210
3211 /* Convenience structure for WALL_QUAD: its is legal to call
3212 *
3213 * get_quad_fast(bas_fcts, WALL_QUAD::quad[wall], ...)
3214 *
3215 * individually, however
3216 *
3217 * get_wall_quad_fast()
3218 *
3219 * does this in a single run.
3220 *
3221 * If INIT_ELEMENT(el_info, WALL_QUAD_FAST) is called, then the
3222 * sub-ordinate initializers
3223 * INIT_ELEMENT(el_info,WALL_QUAD_FAST::quad_fast[wall]) need not be
3224 * called.
3225 */
3226 struct wall_quad_fast
3227 {
3228 const WALL_QUAD *wall_quad;
3229 const BAS_FCTS *bas_fcts;
3230
3231 FLAGS init_flag;
3232 const QUAD_FAST *quad_fast[N_WALLS_MAX];
3233
3234 INIT_ELEMENT_DECL;
3235 };
3236
3237 /* initialize the meta-data for the given WALL_QUAD, no need to call
3238 * this if the WALL_QUAD has been aquired by get_wall_quad(), only
3239 * needed for externally defined extension quadrature rules.
3240 */
3241 extern void register_wall_quadrature(WALL_QUAD *wall_quad);
3242
3243 /* Return a suitable quadrature for integrating over the given wall
3244 * (neigh number), but the barycentric co-ordinates of QUAD->lambda
3245 * are relative to the neighbour element.
3246 */
3247 extern const QUAD *get_neigh_quad(const EL_INFO *el_info,
3248 const WALL_QUAD *wall_quad,
3249 int neigh);
3250
3251 /* Return a suitable QUAD_FAST structure for integrating over the
3252 * given wall, but relative to the neighbour element. If the returned
3253 * QUAD_FAST object has a per-element initializer, then it must be
3254 * called with an EL_INFO structure for the neighbour element.
3255 *
3256 * It is also legal to just call
3257 *
3258 * get_quad_fast(bas_fcts, get_neigh_quad(el_info, wall_quad, neigh), ...)
3259 *
3260 * but get_neigh_quad_fast() is slightly more efficient.
3261 */
3262 const QUAD_FAST *get_neigh_quad_fast(const EL_INFO *el_info,
3263 const WALL_QUAD_FAST *wqfast,
3264 int neigh);
3265
3266 /*******************************************************************************
3267 * functions supplied by ALBERTA
3268 ******************************************************************************/
3269
3270 /*** file coarsen.c *******************************************************/
3271 extern U_CHAR coarsen(MESH *mesh, FLAGS fill_flags);
3272 extern U_CHAR global_coarsen(MESH *mesh, int no, FLAGS fill_flags);
3273 extern int get_max_level(MESH *mesh);
3274
3275 /*** file dof_admin.c *****************************************************/
3276 extern void free_dof_index(DOF_ADMIN *admin, int dof);
3277 extern DOF get_dof_index(DOF_ADMIN *admin);
3278 extern void enlarge_dof_lists(DOF_ADMIN *admin, int minsize);
3279
3280 extern const DOF_ADMIN *get_vertex_admin(MESH *mesh, FLAGS flags);
3281
3282 extern void test_dof_matrix(DOF_MATRIX *matrix);
3283 extern void dof_matrix_set_diagonal(DOF_MATRIX *matrix, bool diag);
3284 extern void dof_matrix_try_diagonal(DOF_MATRIX *matrix);
3285 extern void add_element_matrix(DOF_MATRIX *matrix,
3286 REAL factor,
3287 const EL_MATRIX *el_matrix,
3288 MatrixTranspose transpose,
3289 const EL_DOF_VEC *row_dof,
3290 const EL_DOF_VEC *col_dof,
3291 const EL_SCHAR_VEC *bound);
3292 extern void add_element_vec(DOF_REAL_VEC *drv, REAL factor,
3293 const EL_REAL_VEC *el_vec,
3294 const EL_DOF_VEC *dof,
3295 const EL_SCHAR_VEC *bound);
3296 extern void add_element_d_vec(DOF_REAL_D_VEC *drdv, REAL factor,
3297 const EL_REAL_D_VEC *el_vec,
3298 const EL_DOF_VEC *dof,
3299 const EL_SCHAR_VEC *bound);
3300 extern void add_element_vec_dow(DOF_REAL_VEC_D *drdv, REAL factor,
3301 const EL_REAL_VEC_D *el_vec,
3302 const EL_DOF_VEC *dof,
3303 const EL_SCHAR_VEC *bound);
3304
3305 extern void update_matrix(DOF_MATRIX *dof_matrix,
3306 const EL_MATRIX_INFO *minfo,
3307 MatrixTranspose transpose);
3308 void update_real_vec(DOF_REAL_VEC *drv, const EL_VEC_INFO *vec_info);
3309 void update_real_d_vec(DOF_REAL_D_VEC *drdv, const EL_VEC_D_INFO *vecd_info);
3310 void update_real_vec_dow(DOF_REAL_VEC_D *drdv, const EL_VEC_INFO_D *vecd_info);
3311
3312 extern void dof_compress(MESH *mesh);
3313 extern void add_dof_compress_hook(const DOF_ADMIN *admin, DOF_COMP_HOOK *hook);
3314 extern void del_dof_compress_hook(DOF_COMP_HOOK *hook);
3315 extern void clear_dof_matrix(DOF_MATRIX *matrix);
3316
3317 extern void print_dof_matrix(const DOF_MATRIX *matrix);
3318 extern void print_dof_real_vec(const DOF_REAL_VEC *drv);
3319 extern void print_dof_real_d_vec(const DOF_REAL_D_VEC *drdv);
3320 extern void print_dof_real_dd_vec(const DOF_REAL_DD_VEC *drdv);
3321 extern void print_dof_real_vec_dow(const DOF_REAL_VEC_D *drvd);
3322 extern void print_real_vec_maple(REAL *vector, int size, const char *vec_name);
3323 extern void print_dof_real_vec_maple(const DOF_REAL_VEC *drv,
3324 const char *vec_name);
3325 extern void print_dof_real_d_vec_maple(const DOF_REAL_D_VEC *drdv,
3326 const char *vec_name);
3327 extern void print_dof_real_vec_dow_maple(const DOF_REAL_VEC_D *drvd,
3328 const char *vec_name);
3329 extern void print_dof_matrix_maple(const DOF_MATRIX *matrix,
3330 const char *matrix_name);
3331 extern void fprint_real_vec_maple(FILE *fp,
3332 REAL *vector, int size, const char *vec_name);
3333 extern void fprint_dof_real_vec_maple(FILE *fp,
3334 const DOF_REAL_VEC *drv,
3335 const char *vec_name);
3336 extern void fprint_dof_real_d_vec_maple(FILE *fp,
3337 const DOF_REAL_D_VEC *drdv,
3338 const char *vec_name);
3339 extern void fprint_dof_real_vec_dow_maple(FILE *fp,
3340 const DOF_REAL_VEC_D *drvd,
3341 const char *vec_name);
3342 extern void fprint_dof_matrix_maple(FILE *fp,
3343 const DOF_MATRIX *matrix,
3344 const char *matrix_name);
3345 extern void file_print_real_vec_maple(const char *file_name, const char *mode,
3346 REAL *vector, int size,
3347 const char *vec_name);
3348 extern void file_print_dof_real_vec_maple(const char *file_name,
3349 const char *mode,
3350 const DOF_REAL_VEC *drv,
3351 const char *vec_name);
3352 extern void file_print_dof_real_d_vec_maple(const char *file_name,
3353 const char *mode,
3354 const DOF_REAL_D_VEC *drdv,
3355 const char *vec_name);
3356 extern void file_print_dof_real_vec_dow_maple(const char *file_name,
3357 const char *mode,
3358 const DOF_REAL_VEC_D *drvd,
3359 const char *vec_name);
3360 extern void file_print_dof_matrix_maple(const char *file_name,
3361 const char *mode,
3362 const DOF_MATRIX *matrix,
3363 const char *matrix_name);
3364 extern void print_dof_ptr_vec(const DOF_PTR_VEC *dpv);
3365 extern void print_dof_int_vec(const DOF_INT_VEC *div);
3366 extern void print_dof_uchar_vec(const DOF_UCHAR_VEC *div);
3367 extern void print_dof_schar_vec(const DOF_SCHAR_VEC *div);
3368 /* BLAS 1 */
3369 extern REAL dof_nrm2(const DOF_REAL_VEC *x);
3370 extern REAL dof_asum(const DOF_REAL_VEC *x);
3371 extern void dof_set(REAL alpha, DOF_REAL_VEC *x);
3372 extern void dof_scal(REAL alpha, DOF_REAL_VEC *x);
3373 extern REAL dof_dot(const DOF_REAL_VEC *x, const DOF_REAL_VEC *y);
3374 extern void dof_copy(const DOF_REAL_VEC *x, DOF_REAL_VEC *y);
3375 extern void dof_axpy(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y);
3376 /* some non BLAS */
3377 extern void dof_xpay(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y);
3378 extern REAL dof_min(const DOF_REAL_VEC *x);
3379 extern REAL dof_max(const DOF_REAL_VEC *x);
3380 /* BLAS 2 */
3381 extern void dof_gemv(MatrixTranspose transpose, REAL alpha,
3382 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3383 const DOF_REAL_VEC *x,
3384 REAL beta, DOF_REAL_VEC *y);
3385 extern void dof_mv(MatrixTranspose transpose,
3386 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3387 const DOF_REAL_VEC *x, DOF_REAL_VEC *y);
3388
3389 /* now the same for REAL_D */
3390 extern void dof_axpy_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
3391 extern void dof_copy_d(const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
3392 extern REAL dof_dot_d(const DOF_REAL_D_VEC *x, const DOF_REAL_D_VEC *y);
3393 extern REAL dof_nrm2_d(const DOF_REAL_D_VEC *x);
3394 extern REAL dof_asum_d(const DOF_REAL_D_VEC *x);
3395 extern void dof_scal_d(REAL alpha, DOF_REAL_D_VEC *x);
3396 extern void dof_set_d(REAL alpha, DOF_REAL_D_VEC *x);
3397 extern void dof_xpay_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
3398 extern REAL dof_min_d(const DOF_REAL_D_VEC *x);
3399 extern REAL dof_max_d(const DOF_REAL_D_VEC *x);
3400
3401 extern void
3402 dof_axpy_dow(REAL alpha, const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
3403 extern void dof_copy_dow(const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
3404 extern REAL dof_dot_dow(const DOF_REAL_VEC_D *x, const DOF_REAL_VEC_D *y);
3405 extern REAL dof_nrm2_dow(const DOF_REAL_VEC_D *x);
3406 extern REAL dof_asum_dow(const DOF_REAL_VEC_D *x);
3407 extern void dof_scal_dow(REAL alpha, DOF_REAL_VEC_D *x);
3408 extern void dof_set_dow(REAL alpha, DOF_REAL_VEC_D *x);
3409 extern void
3410 dof_xpay_dow(REAL alpha, const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
3411 extern REAL dof_min_dow(const DOF_REAL_VEC_D *x);
3412 extern REAL dof_max_dow(const DOF_REAL_VEC_D *x);
3413 /* BLAS 2 for REAL_D */
3414 extern void dof_gemv_d(MatrixTranspose transpose, REAL alpha,
3415 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3416 const DOF_REAL_D_VEC *x,
3417 REAL beta, DOF_REAL_D_VEC *y);
3418 extern void dof_gemv_rdr(MatrixTranspose transpose, REAL alpha,
3419 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3420 const DOF_REAL_VEC *x,
3421 REAL beta, DOF_REAL_D_VEC *y);
3422 extern void dof_gemv_rrd(MatrixTranspose transpose, REAL alpha,
3423 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3424 const DOF_REAL_D_VEC *x,
3425 REAL beta, DOF_REAL_VEC *y);
3426 extern void dof_mv_d(MatrixTranspose transpose,
3427 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3428 const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
3429 extern void dof_mv_rdr(MatrixTranspose transpose,
3430 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3431 const DOF_REAL_VEC *x,
3432 DOF_REAL_D_VEC *y);
3433 extern void dof_mv_rrd(MatrixTranspose transpose,
3434 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3435 const DOF_REAL_D_VEC *x,
3436 DOF_REAL_VEC *y);
3437
3438 extern void dof_gemv_dow(MatrixTranspose transpose, REAL alpha,
3439 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3440 const DOF_REAL_VEC_D *x,
3441 REAL beta, DOF_REAL_VEC_D *y);
3442 extern void dof_gemv_dow_scl(MatrixTranspose transpose, REAL alpha,
3443 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3444 const DOF_REAL_VEC *x,
3445 REAL beta, DOF_REAL_VEC_D *y);
3446 extern void dof_gemv_scl_dow(MatrixTranspose transpose, REAL alpha,
3447 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3448 const DOF_REAL_VEC_D *x,
3449 REAL beta, DOF_REAL_VEC *y);
3450 extern void dof_mv_dow(MatrixTranspose transpose,
3451 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3452 const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
3453 extern void dof_mv_dow_scl(MatrixTranspose transpose,
3454 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3455 const DOF_REAL_VEC *x,
3456 DOF_REAL_VEC_D *y);
3457 extern void dof_mv_scl_dow(MatrixTranspose transpose,
3458 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3459 const DOF_REAL_VEC_D *x,
3460 DOF_REAL_VEC *y);
3461
3462 /* copy operation for DOF_MATRIXes */
3463 extern void dof_matrix_copy(DOF_MATRIX *dst, const DOF_MATRIX *src);
3464
3465 /* low-level administration routines, use with caution */
3466 extern void add_dof_matrix_to_admin(DOF_MATRIX *, DOF_ADMIN *);
3467 extern void add_dof_int_vec_to_admin(DOF_INT_VEC *, DOF_ADMIN *);
3468 extern void add_int_dof_vec_to_admin(DOF_DOF_VEC *, DOF_ADMIN *);
3469 extern void add_dof_dof_vec_to_admin(DOF_DOF_VEC *, DOF_ADMIN *);
3470 extern void add_dof_uchar_vec_to_admin(DOF_UCHAR_VEC *, DOF_ADMIN *);
3471 extern void add_dof_schar_vec_to_admin(DOF_SCHAR_VEC *, DOF_ADMIN *);
3472 extern void add_dof_real_vec_to_admin(DOF_REAL_VEC *, DOF_ADMIN *);
3473 extern void add_dof_real_d_vec_to_admin(DOF_REAL_D_VEC *, DOF_ADMIN *);
3474 extern void add_dof_real_dd_vec_to_admin(DOF_REAL_DD_VEC *, DOF_ADMIN *);
3475 extern void add_dof_ptr_vec_to_admin(DOF_PTR_VEC *, DOF_ADMIN *);
3476 extern void remove_dof_matrix_from_admin(DOF_MATRIX *);
3477 extern void remove_dof_int_vec_from_admin(DOF_INT_VEC *);
3478 extern void remove_dof_dof_vec_from_admin(DOF_DOF_VEC *);
3479 extern void remove_int_dof_vec_from_admin(DOF_DOF_VEC *);
3480 extern void remove_dof_uchar_vec_from_admin(DOF_UCHAR_VEC *);
3481 extern void remove_dof_schar_vec_from_admin(DOF_SCHAR_VEC *);
3482 extern void remove_dof_real_vec_from_admin(DOF_REAL_VEC *);
3483 extern void remove_dof_real_vec_from_admin(DOF_REAL_VEC *);
3484 extern void remove_dof_real_d_vec_from_admin(DOF_REAL_D_VEC *);
3485 extern void remove_dof_real_dd_vec_from_admin(DOF_REAL_DD_VEC *);
3486 extern void remove_dof_ptr_vec_from_admin(DOF_PTR_VEC *);
3487
3488 /*** file wall_quad.c *****************************************************/
3489
3490 extern const WALL_QUAD *wall_quad_from_quad(const QUAD *quad);
3491 extern const WALL_QUAD *get_wall_quad(int dim, int degree);
3492 extern const WALL_QUAD_FAST *
3493 get_wall_quad_fast(const BAS_FCTS *, const WALL_QUAD *, FLAGS init_flag);
3494 extern const QUAD *get_neigh_quad(const EL_INFO *el_info,
3495 const WALL_QUAD *wall_quad,
3496 int wall);
3497 extern const QUAD_FAST *get_neigh_quad_fast(const EL_INFO *el_info,
3498 const WALL_QUAD_FAST *bndry_qfast,
3499 int wall);
3500
3501 /*** file macro.c *********************************************************/
3502 extern void macro_test(MACRO_DATA *data, const char *new_filename);
3503
3504 extern MACRO_DATA *read_macro(const char *name);
3505 extern MACRO_DATA *read_macro_bin(const char *name);
3506 extern MACRO_DATA *read_macro_xdr(const char *name);
3507
3508 extern bool write_macro(MESH *mesh, const char *name);
3509 extern bool write_macro_bin(MESH *mesh, const char *name);
3510 extern bool write_macro_xdr(MESH *mesh, const char *name);
3511
3512 extern bool write_macro_data(MACRO_DATA *data, const char *name);
3513 extern bool write_macro_data_bin(MACRO_DATA *data, const char *name);
3514 extern bool write_macro_data_xdr(MACRO_DATA *data, const char *name);
3515
3516 extern MACRO_DATA *alloc_macro_data(int dim, int nv, int ne);
3517 extern void free_macro_data(MACRO_DATA *data);
3518 extern void compute_neigh_fast(MACRO_DATA *data);
3519 extern void default_boundary(MACRO_DATA *data, U_CHAR type, bool overwrite);
3520
3521 extern MACRO_DATA *mesh2macro_data(MESH *mesh);
3522 extern void
3523 macro_data2mesh(MESH *mesh, const MACRO_DATA *data,
3524 NODE_PROJECTION *(*n_proj)(MESH *,MACRO_EL *,int),
3525 AFF_TRAFO *(*init_wall_trafos)(MESH *, MACRO_EL *, int wall));
3526
3527 /*** file memory.c ********************************************************/
3528 extern MESH *
3529 check_and_get_mesh(int dim, int dow, int neigh,
3530 const char *version, const char *name,
3531 const MACRO_DATA *macro_data,
3532 NODE_PROJ *(*init_node_proj)(MESH *, MACRO_EL *, int),
3533 AFF_TRAFO *(*init_wall_trafo)(MESH *, MACRO_EL *, int wall));
3534 extern void free_dof_admin(DOF_ADMIN *admin, MESH *mesh);
3535 extern void free_int_dof_vec(DOF_DOF_VEC *vec);
3536 extern void free_dof_int_vec(DOF_INT_VEC *vec);
3537 extern void free_dof_dof_vec(DOF_DOF_VEC *vec);
3538 extern void free_dof_matrix(DOF_MATRIX *mat);
3539 extern void free_dof_real_vec(DOF_REAL_VEC *vec);
3540 extern void free_dof_real_d_vec(DOF_REAL_D_VEC *vec);
3541 extern void free_dof_real_dd_vec(DOF_REAL_DD_VEC *vec);
3542 extern void free_dof_real_vec_d(DOF_REAL_VEC_D *vec);
3543 extern void free_dof_schar_vec(DOF_SCHAR_VEC *vec);
3544 extern void free_dof_uchar_vec(DOF_UCHAR_VEC *vec);
3545 extern void free_dof_ptr_vec(DOF_PTR_VEC *vec);
3546 extern void free_fe_space(const FE_SPACE *fe_space);
3547 extern void free_real_d(MESH *mesh, REAL *ptr);
3548 extern void free_matrix_row(const FE_SPACE *, MATRIX_ROW *);
3549 extern void free_real_matrix_row(const FE_SPACE *, MATRIX_ROW_REAL *);
3550 extern void free_real_d_matrix_row(const FE_SPACE *, MATRIX_ROW_REAL_D *);
3551 extern void free_real_dd_matrix_row(const FE_SPACE *, MATRIX_ROW_REAL_DD *);
3552 extern void free_element(EL *el, MESH *mesh);
3553 extern void free_rc_list(MESH *mesh, RC_LIST_EL *list); /* only for 3D */
3554 extern void free_mesh(MESH *);
3555 extern void free_dof(DOF *dof, MESH *mesh, int position, FLAGS flags);
3556
3557 extern DOF *get_dof(MESH *mesh, int position);
3558 extern DOF *get_periodic_dof(MESH *mesh, int position,
3559 const DOF *twin);
3560 extern const FE_SPACE *copy_fe_space(const FE_SPACE *fe_space);
3561 extern const FE_SPACE *clone_fe_space(const FE_SPACE *fe_space, int rdim);
3562 extern const FE_SPACE *get_fe_space(MESH *mesh,
3563 const char *name,
3564 const BAS_FCTS *bas_fcts,
3565 int rdim,
3566 FLAGS adm_flags);
3567 extern const FE_SPACE *get_dof_space(MESH *mesh, const char *name,
3568 const int ndof[N_NODE_TYPES],
3569 FLAGS adm_flags);
3570 extern DOF_INT_VEC *get_dof_int_vec(const char *name, const FE_SPACE *);
3571 extern DOF_DOF_VEC *get_int_dof_vec(const char *name, const FE_SPACE *);
3572 extern DOF_DOF_VEC *get_dof_dof_vec(const char *name, const FE_SPACE *);
3573 extern DOF_MATRIX *get_dof_matrix(const char *name,
3574 const FE_SPACE *row_fe_space,
3575 const FE_SPACE *col_fe_space);
3576 extern DOF_REAL_VEC *get_dof_real_vec(const char *name, const FE_SPACE *);
3577 extern DOF_REAL_D_VEC *get_dof_real_d_vec(const char *name, const FE_SPACE *);
3578 extern DOF_REAL_DD_VEC *get_dof_real_dd_vec(const char *name, const FE_SPACE *);
3579 extern DOF_REAL_VEC_D *get_dof_real_vec_d(const char *name, const FE_SPACE *);
3580 extern DOF_SCHAR_VEC *get_dof_schar_vec(const char *name, const FE_SPACE *);
3581 extern DOF_UCHAR_VEC *get_dof_uchar_vec(const char *name, const FE_SPACE *);
3582 extern DOF_PTR_VEC *get_dof_ptr_vec(const char *name, const FE_SPACE *);
3583 extern REAL *get_real_d(MESH *mesh);
3584 extern MATRIX_ROW *get_matrix_row(const FE_SPACE *, MATENT_TYPE type);
3585 extern EL *get_element(MESH *mesh);
3586 extern RC_LIST_EL *get_rc_list(MESH *mesh); /* only for 3D */
3587 extern size_t init_leaf_data(MESH *mesh, size_t size,
3588 void (*refine_leaf_data)(EL *parent,
3589 EL *child[2]),
3590 void (*coarsen_leaf_data)(EL *parent,
3591 EL *child[2]));
3592 extern EL_MATRIX *get_el_matrix(const FE_SPACE *row_fe_space,
3593 const FE_SPACE *col_fe_space,
3594 MATENT_TYPE op_type);
3595 extern void free_el_matrix(EL_MATRIX *el_mat);
3596 extern void print_el_matrix(const EL_MATRIX *el_mat);
3597
3598 extern EL_INT_VEC *get_el_int_vec(const BAS_FCTS *bas_fcts);
3599 extern void free_el_int_vec(EL_INT_VEC *el_vec);
3600 extern EL_DOF_VEC *get_el_dof_vec(const BAS_FCTS *bas_fcts);
3601 extern void free_el_dof_vec(EL_DOF_VEC *el_vec);
3602 extern EL_UCHAR_VEC *get_el_uchar_vec(const BAS_FCTS *bas_fcts);
3603 extern void free_el_uchar_vec(EL_UCHAR_VEC *el_vec);
3604 extern EL_SCHAR_VEC *get_el_schar_vec(const BAS_FCTS *bas_fcts);
3605 extern void free_el_schar_vec(EL_SCHAR_VEC *el_vec);
3606 extern EL_BNDRY_VEC *get_el_bndry_vec(const BAS_FCTS *bas_fcts);
3607 extern void free_el_bndry_vec(EL_BNDRY_VEC *el_vec);
3608 extern EL_PTR_VEC *get_el_ptr_vec(const BAS_FCTS *bas_fcts);
3609 extern void free_el_ptr_vec(EL_PTR_VEC *el_vec);
3610 extern EL_REAL_VEC *get_el_real_vec(const BAS_FCTS *bas_fcts);
3611 extern void free_el_real_vec(EL_REAL_VEC *el_vec);
3612 extern EL_REAL_D_VEC *get_el_real_d_vec(const BAS_FCTS *bas_fcts);
3613 extern void free_el_real_d_vec(EL_REAL_D_VEC *el_vec);
3614 extern EL_REAL_DD_VEC *get_el_real_dd_vec(const BAS_FCTS *bas_fcts);
3615 extern void free_el_real_dd_vec(EL_REAL_DD_VEC *el_vec);
3616 extern EL_REAL_VEC_D *get_el_real_vec_d(const BAS_FCTS *bas_fcts);
3617 extern void free_el_real_vec_d(EL_REAL_VEC_D *el_vec);
3618
3619 /*** file submesh.c ********************************************************/
3620 extern MESH *get_submesh(MESH *master, const char *name,
3621 bool (*binding_method)(MESH *master, MACRO_EL *el,
3622 int wall, void *data),
3623 void *data);
3624 extern MESH *get_bndry_submesh(MESH *master, const char *name);
3625 extern MESH *read_bndry_submesh_xdr(MESH *master, const char *slave_filename);
3626 extern MESH *get_bndry_submesh_by_type(MESH *master, const char *name,
3627 BNDRY_TYPE type);
3628 extern MESH *read_bndry_submesh_by_type_xdr(MESH *master,
3629 const char *slave_filename,
3630 BNDRY_TYPE type);
3631
3632 extern MESH *get_bndry_submesh_by_segment(MESH *master,
3633 const char *name,
3634 const BNDRY_FLAGS segment);
3635 extern MESH *read_bndry_submesh_by_segment(MESH *master,
3636 const char *slave_filename,
3637 const BNDRY_FLAGS segment);
3638 extern MESH *read_bndry_submesh_by_segment_xdr(MESH *master,
3639 const char *slave_filename,
3640 const BNDRY_FLAGS segment);
3641
3642 extern
3643 MESH *lookup_submesh_by_binding(MESH *master,
3644 bool (*binding_method)(MESH *master,
3645 MACRO_EL *el, int wall,
3646 void *data),
3647 void *data);
3648 extern
3649 MESH *lookup_bndry_submesh_by_type(MESH *master, BNDRY_TYPE type);
3650 extern
3651 MESH *lookup_bndry_submesh_by_segment(MESH *master, const BNDRY_FLAGS segment);
3652 extern
3653 MESH *lookup_bndry_submesh(MESH *master);
3654
3655 extern void unchain_submesh(MESH *slave);
3656
3657 extern void bind_submesh(MESH *master,
3658 MESH *slave,
3659 bool (*binding_method)(MESH *master, MACRO_EL *el,
3660 int wall, void *data),
3661 void *data);
3662
3663 extern MESH *read_submesh(MESH *master,
3664 const char *slave_filename,
3665 bool (*binding_method)(MESH *master, MACRO_EL *el,
3666 int wall, void *data),
3667 void *data);
3668
3669 extern MESH *read_submesh_xdr(MESH *master,
3670 const char *slave_filename,
3671 bool (*binding_method)(MESH *master, MACRO_EL *el,
3672 int wall, void *data),
3673 void *data);
3674 extern MESH *lookup_submesh_by_id(MESH *mesh, int id);
3675 extern MESH *lookup_submesh_by_name(MESH *mesh, const char *name);
3676
3677 extern void get_slave_dof_mapping(const FE_SPACE *m_fe_space,
3678 DOF_INT_VEC *s_map);
3679 extern MESH *get_master(MESH *slave);
3680 extern const EL_DOF_VEC *get_master_dof_indices(EL_DOF_VEC *result,
3681 const EL_INFO *s_el_info,
3682 const FE_SPACE *m_fe_space);
3683 extern const EL_BNDRY_VEC *get_master_bound(EL_BNDRY_VEC *result,
3684 const EL_INFO *s_el_info,
3685 const BAS_FCTS *m_bas_fcts);
3686 extern void
3687 fill_master_el_info(EL_INFO *mst_el_info,
3688 const EL_INFO *slv_el_info, FLAGS fill_flags);
3689 extern const EL *get_slave_el(const EL *el, int wall, MESH *trace_mesh);
3690 extern void fill_slave_el_info(EL_INFO *slv_el_info,
3691 const EL_INFO *el_info, int wall,
3692 MESH *trace_mesh);
3693
3694 void trace_to_bulk_coords_2d(REAL_B result,
3695 const REAL_B lambda,
3696 const EL_INFO *el_info);
3697 void trace_to_bulk_coords_1d(REAL_B result,
3698 const REAL_B lambda,
3699 const EL_INFO *el_info);
3700 void trace_to_bulk_coords_0d(REAL_B result,
3701 const REAL_B lambda,
3702 const EL_INFO *el_info);
3703
3704 void bulk_to_trace_coords_2d(REAL_B result,
3705 const REAL_B lambda,
3706 const EL_INFO *el_info);
3707 void bulk_to_trace_coords_1d(REAL_B result,
3708 const REAL_B lambda,
3709 const EL_INFO *el_info);
3710 void bulk_to_trace_coords_0d(REAL_B result,
3711 const REAL_B lambda,
3712 const EL_INFO *el_info);
3713
3714 static inline
trace_to_bulk_coords(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)3715 void trace_to_bulk_coords(REAL_B result,
3716 const REAL_B lambda,
3717 const EL_INFO *el_info)
3718 {
3719 FUNCNAME("trace_to_bulk_coords");
3720 switch (el_info->mesh->dim) {
3721 case 2: trace_to_bulk_coords_2d(result, lambda, el_info); break;
3722 case 1: trace_to_bulk_coords_1d(result, lambda, el_info); break;
3723 case 0: trace_to_bulk_coords_0d(result, lambda, el_info); break;
3724 default:
3725 ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
3726 break;
3727 }
3728 }
3729
3730 static inline
trace_to_bulk_coords_dim(int dim,REAL_B result,const REAL_B lambda,const EL_INFO * el_info)3731 void trace_to_bulk_coords_dim(int dim,
3732 REAL_B result,
3733 const REAL_B lambda,
3734 const EL_INFO *el_info)
3735 {
3736 FUNCNAME("trace_to_bulk_coords_dim");
3737 switch (dim) {
3738 case 2: trace_to_bulk_coords_2d(result, lambda, el_info); break;
3739 case 1: trace_to_bulk_coords_1d(result, lambda, el_info); break;
3740 case 0: trace_to_bulk_coords_0d(result, lambda, el_info); break;
3741 default:
3742 ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
3743 break;
3744 }
3745 }
3746
3747 static inline
bulk_to_trace_coords(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)3748 void bulk_to_trace_coords(REAL_B result,
3749 const REAL_B lambda,
3750 const EL_INFO *el_info)
3751 {
3752 FUNCNAME("bulk_to_trace_coords");
3753 switch (el_info->mesh->dim) {
3754 case 2: bulk_to_trace_coords_2d(result, lambda, el_info); break;
3755 case 1: bulk_to_trace_coords_1d(result, lambda, el_info); break;
3756 case 0: bulk_to_trace_coords_0d(result, lambda, el_info); break;
3757 default:
3758 ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
3759 break;
3760 }
3761 }
3762
3763 static inline
bulk_to_trace_coords_dim(int dim,REAL_B result,const REAL_B lambda,const EL_INFO * el_info)3764 void bulk_to_trace_coords_dim(int dim,
3765 REAL_B result,
3766 const REAL_B lambda,
3767 const EL_INFO *el_info)
3768 {
3769 FUNCNAME("bulk_to_trace_coords_dim");
3770 switch (dim) {
3771 case 2: bulk_to_trace_coords_2d(result, lambda, el_info); break;
3772 case 1: bulk_to_trace_coords_1d(result, lambda, el_info); break;
3773 case 0: bulk_to_trace_coords_0d(result, lambda, el_info); break;
3774 default:
3775 ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
3776 break;
3777 }
3778 }
3779
3780 #define TRACE_DOF_VEC_PROTO(TYPE, typename) \
3781 extern void trace_##typename##_vec(DOF_##TYPE##_VEC *svec, \
3782 const DOF_##TYPE##_VEC *mvec)
3783 TRACE_DOF_VEC_PROTO(REAL, dof_real);
3784 TRACE_DOF_VEC_PROTO(REAL_D, dof_real_d);
3785 TRACE_DOF_VEC_PROTO(INT, dof_int);
3786 TRACE_DOF_VEC_PROTO(DOF, dof_dof);
3787 TRACE_DOF_VEC_PROTO(DOF, int_dof);
3788 TRACE_DOF_VEC_PROTO(UCHAR, dof_uchar);
3789 TRACE_DOF_VEC_PROTO(SCHAR, dof_schar);
3790 TRACE_DOF_VEC_PROTO(PTR, dof_ptr);
3791 extern
3792 void trace_dof_real_vec_d(DOF_REAL_VEC_D *svec, const DOF_REAL_VEC_D *mvec);
3793
3794 extern void update_master_matrix(DOF_MATRIX *m_dof_matrix,
3795 const EL_MATRIX_INFO *s_minfo,
3796 MatrixTranspose transpose);
3797 extern void update_master_real_vec(DOF_REAL_VEC *m_drv,
3798 const EL_VEC_INFO *s_vec_info);
3799 extern void update_master_real_d_vec(DOF_REAL_D_VEC *m_drdv,
3800 const EL_VEC_D_INFO *s_vecd_info);
3801
3802 /*** file level.c ********************************************************/
3803 extern REAL level_element_det_2d(const REAL_D coord[]);
3804 extern void level_coord_to_world_2d(const REAL_D coord[],
3805 const REAL_B lambda,
3806 REAL_D world);
3807 extern void level_coord_to_el_coord_2d(const REAL_B v_lambda[],
3808 const REAL_B lambda,
3809 REAL_B el_lambda);
3810 extern REAL level_element_det_3d(const REAL_D coord[]);
3811 extern void level_coord_to_world_3d(const REAL_D coord[],
3812 const REAL_B lambda,
3813 REAL_D world);
3814 extern void level_coord_to_el_coord_3d(const REAL_B v_lambda[],
3815 const REAL_B lambda,
3816 REAL_B el_lambda);
3817
3818 extern int find_level(MESH *mesh, FLAGS fill_flag, const DOF_REAL_VEC *Level,
3819 REAL value,
3820 int (*init)(const EL_INFO *el_info,
3821 REAL v[],
3822 int N, int wall, const REAL_B lambda[]),
3823 void (*cal)(const EL_INFO *el_info,
3824 REAL v[],
3825 int i,
3826 int wall, const REAL_B lambda[],
3827 const REAL_D coord[]));
3828 extern void set_element_mark(MESH *mesh, FLAGS fill_flag, S_CHAR mark);
3829
3830 /*** file numint.c ********************************************************/
3831 const QUAD *get_quadrature(int dim, int degree);
3832 void register_quadrature(QUAD *quad);
3833 bool new_quadrature(const QUAD *quad);
3834 const QUAD *get_product_quad(const QUAD *quad);
3835 const QUAD *get_lumping_quadrature(int dim);
3836 static inline const QUAD_EL_CACHE *fill_quad_el_cache(const EL_INFO *el_info,
3837 const QUAD *quad,
3838 FLAGS need);
3839 void print_quadrature(const QUAD *quad);
3840 void check_quadrature(const QUAD *quad);
3841
3842 REAL integrate_std_simp(const QUAD *quad, REAL (*f)(const REAL_B lambda));
3843
3844 #ifndef __cplusplus
3845 /* These are some functions defined in evaluate.h. Due to the fact
3846 * that C++ still lacks C99 variable length array support (at least in
3847 * general) evaluate.h is not included, and anything depending on it
3848 * is disabled.
3849 */
3850 static inline
3851 const REAL *f_at_qp(REAL quad_vec[],
3852 const QUAD *quad, REAL (*f)(const REAL_B lambda));
3853 static inline
3854 const REAL_D *f_d_at_qp(REAL_D quad_vec[],
3855 const QUAD *quad,
3856 const REAL *(*f)(const REAL_B lambda));
3857 static inline
3858 const REAL_D *grd_f_at_qp(REAL_D quad_vec[],
3859 const QUAD *quad,
3860 const REAL *(*grd_f)(const REAL_B lambda));
3861 static inline
3862 const REAL_DD *grd_f_d_at_qp(REAL_DD quad_vec[],
3863 const QUAD *quad,
3864 const REAL_D *(*grd_f)(const REAL_B lambda));
3865
3866 static inline
3867 const REAL *fx_at_qp(REAL quad_vec[],
3868 const EL_INFO *el_info,
3869 const QUAD *quad, FCT_AT_X f);
3870 static inline
3871 const REAL_D *fx_d_at_qp(REAL_D quad_vec[],
3872 const EL_INFO *el_info,
3873 const QUAD *quad,
3874 FCT_D_AT_X f);
3875 static inline
3876 const REAL_D *grd_fx_at_qp(REAL_D quad_vec[],
3877 const EL_INFO *el_info,
3878 const QUAD *quad,
3879 GRD_FCT_AT_X grd_f);
3880 static inline
3881 const REAL_DD *grd_fx_d_at_qp(REAL_DD quad_vec[],
3882 const EL_INFO *el_info,
3883 const QUAD *quad,
3884 GRD_FCT_D_AT_X grd_f);
3885
3886 static inline
3887 const REAL *f_loc_at_qp(REAL quad_vec[],
3888 const EL_INFO *el_info, const QUAD *quad,
3889 LOC_FCT_AT_QP f_at_qp, void *ud);
3890 static inline
3891 const REAL_D *f_loc_d_at_qp(REAL_D quad_vec[],
3892 const EL_INFO *el_info, const QUAD *quad,
3893 LOC_FCT_D_AT_QP f_at_qp, void *ud);
3894 static inline
3895 const REAL_D *grd_f_loc_at_qp(REAL_D quad_vec[],
3896 const EL_INFO *el_info, const QUAD *quad,
3897 const REAL_BD Lambda,
3898 GRD_LOC_FCT_AT_QP grd_f_at_qp, void *ud);
3899 static inline
3900 const REAL_DD *grd_f_loc_d_at_qp(REAL_DD quad_vec[],
3901 const EL_INFO *el_info, const QUAD *quad,
3902 const REAL_BD Lambda,
3903 GRD_LOC_FCT_D_AT_QP grd_f_at_qp, void *ud);
3904 static inline
3905 const REAL_D *
3906 param_grd_f_loc_at_qp(REAL_D quad_vec[],
3907 const EL_INFO *el_info,
3908 const QUAD *quad, const REAL_BD Lambda[],
3909 GRD_LOC_FCT_AT_QP grd_f_at_qp, void *ud);
3910 static inline
3911 const REAL_DD *
3912 param_grd_f_loc_d_at_qp(REAL_DD quad_vec[],
3913 const EL_INFO *el_info,
3914 const QUAD *quad, const REAL_BD Lambda[],
3915 GRD_LOC_FCT_D_AT_QP grd_f_at_qp, void *ud);
3916 #endif /* __cplusplus */
3917
3918 const QUAD_FAST *get_quad_fast(const BAS_FCTS *, const QUAD *, FLAGS init_flag);
3919
3920 const REAL_D *const*get_quad_fast_phi_dow(const QUAD_FAST *cache);
3921 const REAL_DB *const*get_quad_fast_grd_phi_dow(const QUAD_FAST *cache);
3922 const REAL_DBB *const*get_quad_fast_D2_phi_dow(const QUAD_FAST *cache);
3923
3924 /*** file refine.c ********************************************************/
3925 extern U_CHAR refine(MESH *mesh, FLAGS fill_flags);
3926 extern U_CHAR global_refine(MESH *mesh, int mark, FLAGS fill_flags);
3927
3928 /******************************************************************************/
3929
3930 /*** file adapt.c *********************************************************/
3931 extern void adapt_method_stat(MESH *mesh, ADAPT_STAT *adapt);
3932 extern void adapt_method_instat(MESH *mesh, ADAPT_INSTAT *adapt);
3933 extern int marking(MESH *mesh, ADAPT_STAT *adapt);
3934 extern ADAPT_INSTAT *get_adapt_instat(int dim, const char *name,
3935 const char *prefix,
3936 int info, ADAPT_INSTAT *adapt_instat);
3937 extern ADAPT_STAT *get_adapt_stat(int dim, const char *name,
3938 const char *prefix,
3939 int info, ADAPT_STAT *adapt_stat);
3940
3941 /*** file quad_cache.c ******************************************************/
3942 extern const Q00_PSI_PHI *get_q00_psi_phi(const BAS_FCTS *psi,
3943 const BAS_FCTS *phi,
3944 const QUAD *quad);
3945 extern const Q01_PSI_PHI *get_q01_psi_phi(const BAS_FCTS *psi,
3946 const BAS_FCTS *phi,
3947 const QUAD *quad);
3948 extern const Q10_PSI_PHI *get_q10_psi_phi(const BAS_FCTS *psi,
3949 const BAS_FCTS *phi,
3950 const QUAD *quad);
3951 extern const Q11_PSI_PHI *get_q11_psi_phi(const BAS_FCTS *psi,
3952 const BAS_FCTS *phi,
3953 const QUAD *quad);
3954
3955 extern const Q001_ETA_PSI_PHI *get_q001_eta_psi_phi(const BAS_FCTS *eta,
3956 const BAS_FCTS *psi,
3957 const BAS_FCTS *phi,
3958 const QUAD *quad);
3959 extern const Q010_ETA_PSI_PHI *get_q010_eta_psi_phi(const BAS_FCTS *eta,
3960 const BAS_FCTS *psi,
3961 const BAS_FCTS *phi,
3962 const QUAD *quad);
3963 extern const Q100_ETA_PSI_PHI *get_q100_eta_psi_phi(const BAS_FCTS *eta,
3964 const BAS_FCTS *psi,
3965 const BAS_FCTS *phi,
3966 const QUAD *quad);
3967
3968 /*** file assemble.c ******************************************************/
3969 extern const EL_MATRIX_INFO *fill_matrix_info(const OPERATOR_INFO *,
3970 EL_MATRIX_INFO *res);
3971 extern const EL_MATRIX_INFO *fill_matrix_info_ext(EL_MATRIX_INFO *res,
3972 const OPERATOR_INFO *,
3973 const BNDRY_OPERATOR_INFO *,
3974 ...);
3975 extern const QUAD_TENSOR *get_quad_matrix(const FE_SPACE *row_fe_space,
3976 const FE_SPACE *col_fe_space,
3977 int krn_degree,
3978 int n_derivatives);
3979 extern const QUAD_TENSOR *get_quad_tensor(const FE_SPACE *row_fe_space,
3980 const FE_SPACE *col_fe_space,
3981 const FE_SPACE *depth_fe_space,
3982 int krn_degree,
3983 int n_derivatives);
3984 extern void free_quad_tensor(const QUAD_TENSOR *qtensor);
3985
3986 /*** file assemble-instat.c ***********************************************/
3987
3988 extern EL_SYS_INFO_INSTAT *
3989 fill_sys_info_instat(const OPERATOR_INFO *stiff_info,
3990 const OPERATOR_INFO *mass_info,
3991 const DOF_REAL_VEC *u_h);
3992 extern EL_SYS_INFO_DOW_INSTAT *
3993 fill_sys_info_instat_dow(const OPERATOR_INFO *stiff_info,
3994 const OPERATOR_INFO *mass_info,
3995 const DOF_REAL_VEC_D *u_h);
3996 static inline EL_SYS_INFO_D_INSTAT *
fill_sys_info_instat_d(const OPERATOR_INFO * stiff_info,const OPERATOR_INFO * mass_info,const DOF_REAL_D_VEC * u_h)3997 fill_sys_info_instat_d(const OPERATOR_INFO *stiff_info,
3998 const OPERATOR_INFO *mass_info,
3999 const DOF_REAL_D_VEC *u_h)
4000 {
4001 return (EL_SYS_INFO_D_INSTAT *)
4002 fill_sys_info_instat_dow(
4003 stiff_info, mass_info, (const DOF_REAL_VEC_D *)u_h);
4004 }
4005 extern void free_sys_info_instat(EL_SYS_INFO_INSTAT *elsii);
4006 extern void free_sys_info_d_instat(EL_SYS_INFO_D_INSTAT *elsii);
4007 extern void update_system_instat(DOF_MATRIX *dof_matrix,
4008 DOF_REAL_VEC *f_h,
4009 REAL tau,
4010 REAL theta,
4011 EL_SYS_INFO_INSTAT *elsii);
4012 extern void update_system_instat_dow(DOF_MATRIX *dof_matrix,
4013 DOF_REAL_VEC_D *f_h,
4014 REAL tau,
4015 REAL theta,
4016 EL_SYS_INFO_DOW_INSTAT *elsii);
4017 static inline
update_system_instat_d(DOF_MATRIX * dof_matrix,DOF_REAL_D_VEC * f_h,REAL tau,REAL theta,EL_SYS_INFO_D_INSTAT * elsii)4018 void update_system_instat_d(DOF_MATRIX *dof_matrix,
4019 DOF_REAL_D_VEC *f_h,
4020 REAL tau,
4021 REAL theta,
4022 EL_SYS_INFO_D_INSTAT *elsii)
4023 {
4024 update_system_instat_dow(dof_matrix, (DOF_REAL_VEC_D *)f_h, tau, theta,
4025 (EL_SYS_INFO_DOW_INSTAT *)elsii);
4026 }
4027
4028 /*** file bas_fct.c *******************************************************/
4029 extern const BAS_FCTS *get_bas_fcts(int dim, const char *name);
4030 extern const BAS_FCTS *get_discontinuous_lagrange(int dim, int degree);
4031 extern const BAS_FCTS *get_lagrange(int dim, int degree);
4032 extern const BAS_FCTS *get_disc_ortho_poly(int dim, int degree);
4033 extern const BAS_FCTS *new_bas_fcts(const BAS_FCTS * bas_fcts);
4034 extern BAS_FCTS *chain_bas_fcts(const BAS_FCTS *head, BAS_FCTS *tail);
4035
4036 typedef const BAS_FCTS *
4037 (*BAS_FCTS_INIT_FCT)(int dim, int dow, const char *name);
4038
4039 extern void add_bas_fcts_plugin(BAS_FCTS_INIT_FCT init_fct);
4040
4041 extern const EL_INT_VEC *
4042 default_get_int_vec(int rvec[], const EL *el, const DOF_INT_VEC *vec);
4043 extern const EL_REAL_VEC *
4044 default_get_real_vec(REAL rvec[], const EL *el, const DOF_REAL_VEC *vec);
4045 extern const EL_REAL_D_VEC *
4046 default_get_real_d_vec(REAL_D rvec[],
4047 const EL *el, const DOF_REAL_D_VEC *vec);
4048 extern const EL_REAL_DD_VEC *
4049 default_get_real_dd_vec(REAL_DD rvec[],
4050 const EL *el, const DOF_REAL_DD_VEC *vec);
4051 extern const EL_REAL_VEC_D *
4052 default_get_real_vec_d(REAL rvec[],
4053 const EL *el, const DOF_REAL_VEC_D *vec);
4054 extern const EL_UCHAR_VEC *
4055 default_get_uchar_vec(U_CHAR rvec[],
4056 const EL *el, const DOF_UCHAR_VEC *vec);
4057 extern const EL_SCHAR_VEC *
4058 default_get_schar_vec(S_CHAR rvec[],
4059 const EL *el, const DOF_SCHAR_VEC *vec);
4060 extern const EL_PTR_VEC *
4061 default_get_ptr_vec(void *rvec[],
4062 const EL *el, const DOF_PTR_VEC *vec);
4063
4064 /*** file check.c *********************************************************/
4065 extern void check_mesh(MESH *mesh);
4066
4067 /*** file element.c *******************************************************/
4068
4069 /* These routines are partially available as _?d-versions to avoid looking
4070 * up the dimension. This should be a small efficiency bonus.
4071 */
4072 extern void
4073 fill_neigh_el_info(EL_INFO *neigh_info,
4074 const EL_INFO *el_info, int wall, int rel_perm);
4075 static inline const EL_GEOM_CACHE *
4076 fill_el_geom_cache(const EL_INFO *el_info, FLAGS fill_flag);
4077
4078 #if 0
4079 /* implemented as inline-functions further below */
4080 extern int wall_orientation(int dim, const EL *el, int wall);
4081 extern int world_to_coord(const EL_INFO *el_info, const REAL *,
4082 REAL_B);
4083 extern const REAL *coord_to_world(const EL_INFO *, const REAL_B, REAL_D);
4084 extern REAL el_det(const EL_INFO *el_info);
4085 extern REAL el_volume(const EL_INFO *el_info);
4086 extern REAL el_grd_lambda(const EL_INFO *el_info,
4087 REAL_BD grd_lam);
4088 #endif /* inlines further below */
4089
4090 /* Dimension dependent routines, 0d, just dummies in most cases. */
4091 extern int wall_orientation_0d(const EL *el, int wall);
4092 extern int wall_rel_orientation_0d(const EL *el, const EL *neigh,
4093 int wall, int ov);
4094 extern int world_to_coord_0d(const EL_INFO *el_info, const REAL *, REAL_B);
4095 extern const REAL *coord_to_world_0d(const EL_INFO *, const REAL_B, REAL_D);
4096 extern REAL el_det_0d(const EL_INFO *);
4097 extern REAL el_volume_0d(const EL_INFO *el_info);
4098 extern REAL el_grd_lambda_0d(const EL_INFO *el_info,
4099 REAL_BD grd_lam);
4100 extern REAL get_wall_normal_0d(const EL_INFO *el_info, int wall, REAL *normal);
4101
4102 /* Dimension dependent routines, 1d */
4103 extern int wall_orientation_1d(const EL *el, int wall);
4104 extern int wall_rel_orientation_1d(const EL *el, const EL *neigh,
4105 int wall, int ov);
4106 extern int world_to_coord_1d(const EL_INFO *el_info, const REAL *,
4107 REAL_B);
4108 extern const REAL *coord_to_world_1d(const EL_INFO *, const REAL_B, REAL_D);
4109 extern REAL el_det_1d(const EL_INFO *);
4110 extern REAL el_volume_1d(const EL_INFO *el_info);
4111 extern REAL el_grd_lambda_1d(const EL_INFO *,
4112 REAL_BD grd_lam);
4113 extern REAL get_wall_normal_1d(const EL_INFO *el_info, int wall, REAL *normal);
4114
4115 #if DIM_MAX > 1
4116 /* Dimension dependent routines, 2d */
4117 extern int wall_orientation_2d(const EL *el, int wall);
4118 extern int wall_rel_orientation_2d(const EL *el, const EL *neigh,
4119 int wall, int ov);
4120 extern int world_to_coord_2d(const EL_INFO *el_info, const REAL *,
4121 REAL_B);
4122 extern const REAL *coord_to_world_2d(const EL_INFO *, const REAL_B, REAL_D);
4123 extern REAL el_det_2d(const EL_INFO *);
4124 extern REAL el_volume_2d(const EL_INFO *el_info);
4125 extern REAL el_grd_lambda_2d(const EL_INFO *,
4126 REAL_BD grd_lam);
4127 extern REAL get_wall_normal_2d(const EL_INFO *el_info, int wall, REAL *normal);
4128
4129 #if DIM_MAX > 2
4130 /* Dimension dependent routines, 3d */
4131 extern int wall_orientation_3d(const EL *el, int wall);
4132 extern int wall_rel_orientation_3d(const EL *el, const EL *neigh,
4133 int wall, int oppv);
4134 extern int world_to_coord_3d(const EL_INFO *el_info, const REAL *, REAL_B);
4135 extern const REAL *coord_to_world_3d(const EL_INFO *, const REAL_B, REAL_D);
4136 extern REAL el_det_3d(const EL_INFO *);
4137 extern REAL el_volume_3d(const EL_INFO *el_info);
4138 extern REAL el_grd_lambda_3d(const EL_INFO *,
4139 REAL_BD grd_lam);
4140 extern REAL get_wall_normal_3d(const EL_INFO *el_info, int wall, REAL *normal);
4141 #endif
4142 #endif
4143
4144 /* Below we provide wrapper functions which distinguish the dimension
4145 * dependent routines by the co-dimension rather than by the dimension
4146 * of the underlying mesh. We start by defining a preprocessor macro
4147 * which spares us some typing and especially typos.
4148 *
4149 * In addition, we provide wrapper functions which decide by looking
4150 * at el_info->mesh->dim what to do.
4151 *
4152 */
4153 #if DIM_OF_WORLD == 1
4154
4155 # define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames) \
4156 static inline ret name##suf argtypes \
4157 { \
4158 FUNCNAME(#name); \
4159 \
4160 switch (dim) { \
4161 case 0: return name##_0d argnames; \
4162 case 1: return name##_1d argnames; \
4163 default: \
4164 ERROR_EXIT("Illegal dim!\n"); \
4165 return (ret)0; /* shut-off a compiler warning */ \
4166 } \
4167 } \
4168 struct _AI_semicolon_dummy
4169
4170 # define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames) \
4171 static inline ret name##_0cd argtypes { return name##_1d argnames; } \
4172 static inline ret name##_1cd argtypes { return name##_0d argnames; } \
4173 struct _AI_semicolon_dummy
4174
4175 /* Variants which start at DOW == 2 and thus are empty here */
4176 # define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames) \
4177 struct _AI_semicolon_dummy
4178 # define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames) \
4179 struct _AI_semicolon_dummy
4180
4181 #elif DIM_OF_WORLD == 2
4182
4183 # define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames) \
4184 static inline ret name##suf argtypes \
4185 { \
4186 FUNCNAME(#name); \
4187 \
4188 switch (dim) { \
4189 case 0: return name##_0d argnames; \
4190 case 1: return name##_1d argnames; \
4191 case 2: return name##_2d argnames; \
4192 default: \
4193 ERROR_EXIT("Illegal dim!\n"); \
4194 return (ret)0L; /* shut-off a compiler warning */ \
4195 } \
4196 } \
4197 struct _AI_semicolon_dummy
4198
4199 # define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames) \
4200 static inline ret name##_0cd argtypes { return name##_2d argnames; } \
4201 static inline ret name##_1cd argtypes { return name##_1d argnames; } \
4202 static inline ret name##_2cd argtypes { return name##_0d argnames; } \
4203 struct _AI_semicolon_dummy
4204
4205 /* Variants which start at DOW == 2 */
4206 # define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames) \
4207 static inline ret name##_0cd argtypes { return name##_2d argnames; } \
4208 struct _AI_semicolon_dummy
4209 # define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames) \
4210 static inline void name##_0cd argtypes { name##_2d argnames; } \
4211 struct _AI_semicolon_dummy
4212
4213 #elif DIM_OF_WORLD == 3
4214
4215 # define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames) \
4216 static inline ret name##suf argtypes \
4217 { \
4218 FUNCNAME(#name); \
4219 \
4220 switch (dim) { \
4221 case 0: return name##_0d argnames; \
4222 case 1: return name##_1d argnames; \
4223 case 2: return name##_2d argnames; \
4224 case 3: return name##_3d argnames; \
4225 default: \
4226 ERROR_EXIT("Illegal dim!\n"); \
4227 return (ret)0L; /* shut-off a compiler warning */ \
4228 } \
4229 } \
4230 struct _AI_semicolon_dummy
4231
4232 # define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames) \
4233 static inline ret name##_0cd argtypes { return name##_3d argnames; } \
4234 static inline ret name##_1cd argtypes { return name##_2d argnames; } \
4235 static inline ret name##_2cd argtypes { return name##_1d argnames; } \
4236 static inline ret name##_3cd argtypes { return name##_0d argnames; } \
4237 struct _AI_semicolon_dummy
4238
4239 /* Variants which start at DOW == 2 */
4240 # define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames) \
4241 static inline ret name##_0cd argtypes { return name##_3d argnames; } \
4242 static inline ret name##_1cd argtypes { return name##_2d argnames; } \
4243 struct _AI_semicolon_dummy
4244 # define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames) \
4245 static inline void name##_0cd argtypes { name##_3d argnames; } \
4246 static inline void name##_1cd argtypes { name##_2d argnames; } \
4247 struct _AI_semicolon_dummy
4248
4249 #elif DIM_OF_WORLD > 3
4250
4251 # define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames) \
4252 static inline ret name##suf argtypes \
4253 { \
4254 FUNCNAME(#name); \
4255 \
4256 switch (dim) { \
4257 case 0: return name##_0d argnames; \
4258 case 1: return name##_1d argnames; \
4259 case 2: return name##_2d argnames; \
4260 case 3: return name##_3d argnames; \
4261 default: \
4262 ERROR_EXIT("Illegal dim!\n"); \
4263 return (ret)0L; /* shut-off a compiler warning */ \
4264 } \
4265 } \
4266 struct _AI_semicolon_dummy
4267
4268 # define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames) \
4269 struct _AI_semicolon_dummy
4270 # define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames) \
4271 struct _AI_semicolon_dummy
4272 # define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames) \
4273 struct _AI_semicolon_dummy
4274
4275 #endif
4276
4277 /* ..._Xcd() alias definitions */
4278 ALBERTA_CODIM_ALIAS(int, world_to_coord,
4279 (const EL_INFO *el_info,
4280 const REAL *xy,
4281 REAL_B lambda),
4282 (el_info, xy, lambda));
4283 ALBERTA_CODIM_ALIAS(const REAL *, coord_to_world,
4284 (const EL_INFO *el_info, const REAL_B l, REAL_D w),
4285 (el_info, l, w));
4286 ALBERTA_CODIM_ALIAS(REAL, el_volume, (const EL_INFO *el_info), (el_info));
4287 ALBERTA_CODIM_ALIAS(REAL, el_det, (const EL_INFO *el_info), (el_info));
4288 ALBERTA_CODIM_ALIAS(REAL, el_grd_lambda,
4289 (const EL_INFO *el_info,
4290 REAL_BD grd_lam),
4291 (el_info, grd_lam));
4292 ALBERTA_CODIM_ALIAS(REAL, get_wall_normal,
4293 (const EL_INFO *el_info, int i0, REAL *normal),
4294 (el_info, i0, normal));
4295 ALBERTA_CODIM_ALIAS(int, wall_orientation,
4296 (const EL *el, int wall),
4297 (el, wall));
4298 ALBERTA_CODIM_ALIAS(int, wall_rel_orientation,
4299 (const EL *el, const EL *neigh, int wall, int oppv),
4300 (el, neigh, wall, oppv));
4301 static const int sorted_wall_vertices_0d[1][1][1] = {{{ 0 }}}; /* dummy */
4302 ALBERTA_CODIM_ALIAS(const int *, sorted_wall_vertices,
4303 (int wall, int permno),
4304 [wall][permno]);
4305 static const int vertex_of_wall_0d[1][1] = {{ 0 }}; /* dummy */
4306 ALBERTA_CODIM_ALIAS(const int *, vertex_of_wall,
4307 (int wall),
4308 [wall]);
4309 static const int vertex_of_edge_0d[1][1] = {{ 0 }}; /* dummy */
4310 ALBERTA_CODIM_ALIAS(const int *, vertex_of_edge,
4311 (int edge),
4312 [edge]);
4313
4314 /* Wrappers which look at el_info->mesh->dim */
4315 ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
4316 int, world_to_coord, /**/,
4317 (const EL_INFO *el_info, const REAL *x, REAL_B lambda),
4318 (el_info, x, lambda));
4319 ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
4320 const REAL *, coord_to_world, /**/,
4321 (const EL_INFO *el_info, const REAL_B lambda, REAL_D x),
4322 (el_info, lambda, x));
4323 ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
4324 REAL, el_volume, /**/,
4325 (const EL_INFO *el_info), (el_info));
4326 ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
4327 REAL, el_det, /**/,
4328 (const EL_INFO *el_info), (el_info));
4329 ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
4330 REAL, el_grd_lambda, /**/,
4331 (const EL_INFO *el_info,
4332 REAL_BD grd_lam),
4333 (el_info, grd_lam));
4334 ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
4335 REAL, get_wall_normal, /**/,
4336 (const EL_INFO *el_info, int wall, REAL *normal),
4337 (el_info, wall, normal));
4338
4339 /* Wrappers with addtional "dim" as argument */
4340 ALBERTA_CODIM_WRAPPER(dim, int, wall_orientation, /**/,
4341 (int dim, const EL *el, int wall),
4342 (el, wall));
4343 ALBERTA_CODIM_WRAPPER(dim, int, wall_rel_orientation, /**/,
4344 (int dim,
4345 const EL *el, const EL *neigh, int wall, int oppv),
4346 (el, neigh, wall, oppv));
4347 ALBERTA_CODIM_WRAPPER(dim, const int *, sorted_wall_vertices, /**/,
4348 (int dim, int wall, int permno), [wall][permno]);
4349 ALBERTA_CODIM_WRAPPER(dim, const int *, vertex_of_wall, /**/,
4350 (int dim, int wall), [wall]);
4351 ALBERTA_CODIM_WRAPPER(dim, const int *, vertex_of_edge, /**/,
4352 (int dim, int edge), [edge]);
4353 ALBERTA_CODIM_WRAPPER(dim, int, world_to_coord, _dim,
4354 (int dim,
4355 const EL_INFO *el_info, const REAL *x, REAL_B lambda),
4356 (el_info, x, lambda));
4357 ALBERTA_CODIM_WRAPPER(dim, const REAL *, coord_to_world, _dim,
4358 (int dim,
4359 const EL_INFO *el_info, const REAL_B lambda, REAL_D x),
4360 (el_info, lambda, x));
4361 ALBERTA_CODIM_WRAPPER(dim, REAL, el_volume, _dim,
4362 (int dim, const EL_INFO *el_info), (el_info));
4363 ALBERTA_CODIM_WRAPPER(dim, REAL, el_det, _dim,
4364 (int dim, const EL_INFO *el_info), (el_info));
4365 ALBERTA_CODIM_WRAPPER(dim, REAL, el_grd_lambda, _dim,
4366 (int dim, const EL_INFO *el_info, REAL_BD grd_lam),
4367 (el_info, grd_lam));
4368 ALBERTA_CODIM_WRAPPER(dim, REAL, get_wall_normal, _dim,
4369 (int dim, const EL_INFO *el_info, int wall, REAL *normal),
4370 (el_info, wall, normal));
4371
4372 /* Some special wrapper functions, used for some stuff defined in
4373 * level.c
4374 */
4375 ALBERTA_CODIM_ALIAS_2(REAL, level_element_det, (const REAL_D coord[]), (coord));
4376 ALBERTA_VOID_CODIM_ALIAS_2(level_coord_to_world,
4377 (const REAL_D coord[],
4378 const REAL_B lambda,
4379 REAL_D world),
4380 (coord, lambda, world));
4381 ALBERTA_VOID_CODIM_ALIAS_2(level_coord_to_el_coord,
4382 (const REAL_B v_lambda[],
4383 const REAL_B lambda,
4384 REAL_B el_lambda),
4385 (v_lambda, lambda, el_lambda));
4386
4387 /*** file estimator{_dowb}.c **********************************************/
4388
4389 /* The values accepted by the f_flags arguments of ellipt_est() &
4390 * friends.
4391 */
4392 #define INIT_UH 1
4393 #define INIT_GRD_UH 2
4394
4395 /*
4396 {
4397 const void *est_handle;
4398 REAL est_el;
4399 const EL_GEOM_CACHE *elgc;
4400
4401 est_handle = estimator_init(...);
4402 TRAVERSE_FIRST(mesh, -1, fill_flag) {
4403 est_el = element_est(el_info, est_handle);
4404 #if needed
4405 uh = element_est_uh(est_handle);
4406 grd_uh = element_est_grd_uh(est_handle);
4407 el_gc = fill_el_geom_cache(el_info, FILL_EL_DET);
4408 est_el += el_gc->det * additional_stuff(el_cache,...);
4409 #endif
4410 element_est_finish(est_el, est_handle);
4411 } TRAVERSE_NEXT();
4412
4413 estimate = estimator_finish(..., est_handle)
4414 }
4415 */
4416
4417 extern const void *ellipt_est_init(const DOF_REAL_VEC *uh,
4418 ADAPT_STAT *adapt,
4419 REAL *(*rw_est)(EL *),
4420 REAL *(*rw_estc)(EL *),
4421 const QUAD *quad,
4422 const WALL_QUAD *wall_quad,
4423 NORM norm,
4424 REAL C[3],
4425 const REAL_DD A,
4426 const BNDRY_FLAGS dirichlet_bndry,
4427 REAL (*f)(const EL_INFO *el_info,
4428 const QUAD *quad,
4429 int qp,
4430 REAL uh_qp,
4431 const REAL_D grd_uh_gp),
4432 FLAGS f_flags,
4433 REAL (*gn)(const EL_INFO *el_info,
4434 const QUAD *quad,
4435 int qp,
4436 REAL uh_qp,
4437 const REAL_D normal),
4438 FLAGS gn_flags);
4439 extern const void *heat_est_init(const DOF_REAL_VEC *uh,
4440 const DOF_REAL_VEC *uh_old,
4441 ADAPT_INSTAT *adapt,
4442 REAL *(*rw_est)(EL *),
4443 REAL *(*rw_estc)(EL *),
4444 const QUAD *quad,
4445 const WALL_QUAD *wall_quad,
4446 REAL C[4],
4447 const REAL_DD A,
4448 const BNDRY_FLAGS dirichlet_bndry,
4449 REAL (*f)(const EL_INFO *el_info,
4450 const QUAD *quad,
4451 int qp,
4452 REAL uh_qp,
4453 const REAL_D grd_uh_gp,
4454 REAL time),
4455 FLAGS f_flags,
4456 REAL (*gn)(const EL_INFO *el_info,
4457 const QUAD *quad,
4458 int qp,
4459 REAL uh_qp,
4460 const REAL_D normal,
4461 REAL time),
4462 FLAGS gn_flags);
4463
4464 extern REAL element_est(const EL_INFO *el_info, const void *est_handle);
4465 extern void element_est_finish(const EL_INFO *el_info,
4466 REAL est_el, const void *est_handle);
4467 extern REAL ellipt_est_finish(ADAPT_STAT *adapt, const void *est_handle);
4468 extern REAL heat_est_finish(ADAPT_INSTAT *adapt, const void *est_handle);
4469 extern const REAL *element_est_uh(const void *est_handle);
4470 extern const REAL_D *element_est_grd_uh(const void *est_handle);
4471 extern const EL_REAL_VEC *element_est_uh_loc(const void *est_handle);
4472 extern const EL_REAL_VEC *element_est_uh_old_loc(const void *est_handle);
4473
4474 extern REAL ellipt_est(const DOF_REAL_VEC *uh,
4475 ADAPT_STAT *adapt,
4476 REAL *(*rw_est)(EL *),
4477 REAL *(*rw_estc)(EL *),
4478 int quad_degree,
4479 NORM norm,
4480 REAL C[3],
4481 const REAL_DD A,
4482 const BNDRY_FLAGS dirichlet_bndry,
4483 REAL (*f)(const EL_INFO *el_info,
4484 const QUAD *quad,
4485 int qp,
4486 REAL uh_qp,
4487 const REAL_D grd_uh_gp),
4488 FLAGS f_flags,
4489 REAL (*gn)(const EL_INFO *el_info,
4490 const QUAD *quad,
4491 int qp,
4492 REAL uh_qp,
4493 const REAL_D normal),
4494 FLAGS gn_flags);
4495 extern REAL heat_est(const DOF_REAL_VEC *uh,
4496 const DOF_REAL_VEC *uh_old,
4497 ADAPT_INSTAT *adapt,
4498 REAL *(*rw_est)(EL *),
4499 REAL *(*rw_estc)(EL *),
4500 int quad_degree,
4501 REAL C[4],
4502 const REAL_DD A,
4503 const BNDRY_FLAGS dirichlet_bndry,
4504 REAL (*f)(const EL_INFO *el_info,
4505 const QUAD *quad,
4506 int qp,
4507 REAL uh_qp,
4508 const REAL_D grd_uh_gp,
4509 REAL time),
4510 FLAGS f_flags,
4511 REAL (*gn)(const EL_INFO *el_info,
4512 const QUAD *quad,
4513 int qp,
4514 REAL uh_qp,
4515 const REAL_D normal,
4516 REAL time),
4517 FLAGS gn_flags);
4518
4519 /*** file estimator_dowb.c **************************************************/
4520 extern const void *ellipt_est_dow_init(const DOF_REAL_VEC_D *uh,
4521 ADAPT_STAT *adapt,
4522 REAL *(*rw_est)(EL *),
4523 REAL *(*rw_estc)(EL *),
4524 const QUAD *quad,
4525 const WALL_QUAD *wall_quad,
4526 NORM norm,
4527 REAL C[3],
4528 const void *A,
4529 MATENT_TYPE A_type,
4530 MATENT_TYPE A_blocktype,
4531 bool sym_grad,
4532 const BNDRY_FLAGS dirichlet_bndry,
4533 const REAL *(*f)(REAL_D result,
4534 const EL_INFO *el_info,
4535 const QUAD *quad,
4536 int qp,
4537 const REAL_D uh_qp,
4538 const REAL_DD grd_uh_gp),
4539 FLAGS f_flags,
4540 const REAL *(*gn)(REAL_D result,
4541 const EL_INFO *el_info,
4542 const QUAD *quad,
4543 int qp,
4544 const REAL_D uh_qp,
4545 const REAL_D normal),
4546 FLAGS gn_flags);
4547 extern const void *heat_est_dow_init(const DOF_REAL_VEC_D *uh,
4548 const DOF_REAL_VEC_D *uh_old,
4549 ADAPT_INSTAT *adapt,
4550 REAL *(*rw_est)(EL *),
4551 REAL *(*rw_estc)(EL *),
4552 const QUAD *quad,
4553 const WALL_QUAD *wall_quad,
4554 REAL C[4],
4555 const void *A,
4556 MATENT_TYPE A_type,
4557 MATENT_TYPE A_blocktype,
4558 bool sym_grad,
4559 const BNDRY_FLAGS dirichlet_bndry,
4560 const REAL *(*f)(REAL_D result,
4561 const EL_INFO *el_info,
4562 const QUAD *quad,
4563 int qp,
4564 const REAL_D uh_qp,
4565 const REAL_DD grd_uh_gp,
4566 REAL time),
4567 FLAGS f_flags,
4568 const REAL *(*gn)(REAL_D result,
4569 const EL_INFO *el_info,
4570 const QUAD *quad,
4571 int qp,
4572 const REAL_D uh_qp,
4573 const REAL_D normal,
4574 REAL time),
4575 FLAGS gn_flags);
4576
4577 extern REAL element_est_dow(const EL_INFO *el_info, const void *est_handle);
4578 extern void element_est_dow_finish(const EL_INFO *el_info,
4579 REAL est_el, const void *est_handle);
4580 extern REAL ellipt_est_dow_finish(ADAPT_STAT *adapt, const void *est_handle);
4581 extern REAL heat_est_dow_finish(ADAPT_INSTAT *adapt, const void *est_handle);
4582 extern const REAL_D *element_est_uh_dow(const void *est_handle);
4583 extern const REAL_DD *element_est_grd_uh_dow(const void *est_handle);
4584 extern const EL_REAL_VEC_D *element_est_uh_loc_dow(const void *est_handle);
4585 extern const EL_REAL_VEC_D *element_est_uh_old_loc_dow(const void *est_handle);
4586
4587 extern REAL ellipt_est_dow(const DOF_REAL_VEC_D *uh,
4588 ADAPT_STAT *adapt,
4589 REAL *(*rw_est)(EL *),
4590 REAL *(*rw_estc)(EL *),
4591 int quad_degree,
4592 NORM norm,
4593 REAL C[3],
4594 const void *A,
4595 MATENT_TYPE A_type,
4596 MATENT_TYPE A_blocktype,
4597 bool sym_grad,
4598 const BNDRY_FLAGS dirichlet_bndry,
4599 const REAL *(*f)(REAL_D result,
4600 const EL_INFO *el_info,
4601 const QUAD *quad,
4602 int qp,
4603 const REAL_D uh_qp,
4604 const REAL_DD grd_uh_gp),
4605 FLAGS f_flags,
4606 const REAL *(*gn)(REAL_D result,
4607 const EL_INFO *el_info,
4608 const QUAD *quad,
4609 int qp,
4610 const REAL_D uh_qp,
4611 const REAL_D normal),
4612 FLAGS gn_flags);
4613 extern REAL heat_est_dow(const DOF_REAL_VEC_D *uh,
4614 const DOF_REAL_VEC_D *uh_old,
4615 ADAPT_INSTAT *adapt,
4616 REAL *(*rw_est)(EL *),
4617 REAL *(*rw_estc)(EL *),
4618 int quad_degree,
4619 REAL C[4],
4620 const void *A,
4621 MATENT_TYPE A_type,
4622 MATENT_TYPE A_blocktype,
4623 bool sym_grad,
4624 const BNDRY_FLAGS dirichlet_bndry,
4625 const REAL *(*f)(REAL_D result,
4626 const EL_INFO *el_info,
4627 const QUAD *quad,
4628 int qp,
4629 const REAL_D uh_qp,
4630 const REAL_DD grd_uh_gp,
4631 REAL time),
4632 FLAGS f_flags,
4633 const REAL *(*gn)(REAL_D result,
4634 const EL_INFO *el_info,
4635 const QUAD *quad,
4636 int qp,
4637 const REAL_D uh_qp,
4638 const REAL_D normal,
4639 REAL time),
4640 FLAGS gn_flags);
4641
4642 /* DOF_REAL_D_VEC versions */
4643 static inline
ellipt_est_d_init(const DOF_REAL_D_VEC * uh,ADAPT_STAT * adapt,REAL * (* rw_est)(EL *),REAL * (* rw_estc)(EL *),const QUAD * quad,const WALL_QUAD * wall_quad,NORM norm,REAL C[3],const void * A,MATENT_TYPE A_type,MATENT_TYPE A_blocktype,bool sym_grad,const BNDRY_FLAGS dirichlet_bndry,const REAL * (* f)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_DD grd_uh_gp),FLAGS f_flags,const REAL * (* gn)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_D normal),FLAGS gn_flags)4644 const void *ellipt_est_d_init(const DOF_REAL_D_VEC *uh,
4645 ADAPT_STAT *adapt,
4646 REAL *(*rw_est)(EL *),
4647 REAL *(*rw_estc)(EL *),
4648 const QUAD *quad,
4649 const WALL_QUAD *wall_quad,
4650 NORM norm,
4651 REAL C[3],
4652 const void *A,
4653 MATENT_TYPE A_type,
4654 MATENT_TYPE A_blocktype,
4655 bool sym_grad,
4656 const BNDRY_FLAGS dirichlet_bndry,
4657 const REAL *(*f)(REAL_D result,
4658 const EL_INFO *el_info,
4659 const QUAD *quad,
4660 int qp,
4661 const REAL_D uh_qp,
4662 const REAL_DD grd_uh_gp),
4663 FLAGS f_flags,
4664 const REAL *(*gn)(REAL_D result,
4665 const EL_INFO *el_info,
4666 const QUAD *quad,
4667 int qp,
4668 const REAL_D uh_qp,
4669 const REAL_D normal),
4670 FLAGS gn_flags)
4671 {
4672 return ellipt_est_dow_init(
4673 (const DOF_REAL_VEC_D *)uh,
4674 adapt, rw_est, rw_estc, quad, wall_quad, norm, C, A, A_type, A_blocktype,
4675 sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
4676 }
4677
4678 static inline
ellipt_est_d(const DOF_REAL_D_VEC * uh,ADAPT_STAT * adapt,REAL * (* rw_est)(EL *),REAL * (* rw_estc)(EL *),int quad_degree,NORM norm,REAL C[3],const void * A,MATENT_TYPE A_type,MATENT_TYPE A_blocktype,bool sym_grad,const BNDRY_FLAGS dirichlet_bndry,const REAL * (* f)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_DD grd_uh_gp),FLAGS f_flags,const REAL * (* gn)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_D normal),FLAGS gn_flags)4679 REAL ellipt_est_d(const DOF_REAL_D_VEC *uh,
4680 ADAPT_STAT *adapt,
4681 REAL *(*rw_est)(EL *),
4682 REAL *(*rw_estc)(EL *),
4683 int quad_degree,
4684 NORM norm,
4685 REAL C[3],
4686 const void *A,
4687 MATENT_TYPE A_type,
4688 MATENT_TYPE A_blocktype,
4689 bool sym_grad,
4690 const BNDRY_FLAGS dirichlet_bndry,
4691 const REAL *(*f)(REAL_D result,
4692 const EL_INFO *el_info,
4693 const QUAD *quad,
4694 int qp,
4695 const REAL_D uh_qp,
4696 const REAL_DD grd_uh_gp),
4697 FLAGS f_flags,
4698 const REAL *(*gn)(REAL_D result,
4699 const EL_INFO *el_info,
4700 const QUAD *quad,
4701 int qp,
4702 const REAL_D uh_qp,
4703 const REAL_D normal),
4704 FLAGS gn_flags)
4705 {
4706 return ellipt_est_dow(
4707 (const DOF_REAL_VEC_D *)uh,
4708 adapt, rw_est, rw_estc, quad_degree, norm, C, A, A_type, A_blocktype,
4709 sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
4710 }
4711
4712 static inline
heat_est_d_init(const DOF_REAL_D_VEC * uh,const DOF_REAL_D_VEC * uh_old,ADAPT_INSTAT * adapt,REAL * (* rw_est)(EL *),REAL * (* rw_estc)(EL *),const QUAD * quad,const WALL_QUAD * wall_quad,REAL C[4],const void * A,MATENT_TYPE A_type,MATENT_TYPE A_blocktype,bool sym_grad,const BNDRY_FLAGS dirichlet_bndry,const REAL * (* f)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_DD grd_uh_gp,REAL time),FLAGS f_flags,const REAL * (* gn)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_D normal,REAL time),FLAGS gn_flags)4713 const void *heat_est_d_init(const DOF_REAL_D_VEC *uh,
4714 const DOF_REAL_D_VEC *uh_old,
4715 ADAPT_INSTAT *adapt,
4716 REAL *(*rw_est)(EL *),
4717 REAL *(*rw_estc)(EL *),
4718 const QUAD *quad,
4719 const WALL_QUAD *wall_quad,
4720 REAL C[4],
4721 const void *A,
4722 MATENT_TYPE A_type,
4723 MATENT_TYPE A_blocktype,
4724 bool sym_grad,
4725 const BNDRY_FLAGS dirichlet_bndry,
4726 const REAL *(*f)(REAL_D result,
4727 const EL_INFO *el_info,
4728 const QUAD *quad,
4729 int qp,
4730 const REAL_D uh_qp,
4731 const REAL_DD grd_uh_gp,
4732 REAL time),
4733 FLAGS f_flags,
4734 const REAL *(*gn)(REAL_D result,
4735 const EL_INFO *el_info,
4736 const QUAD *quad,
4737 int qp,
4738 const REAL_D uh_qp,
4739 const REAL_D normal,
4740 REAL time),
4741 FLAGS gn_flags)
4742 {
4743 return heat_est_dow_init(
4744 (const DOF_REAL_VEC_D *)uh, (const DOF_REAL_VEC_D *)uh_old,
4745 adapt, rw_est, rw_estc, quad, wall_quad, C, A, A_type, A_blocktype,
4746 sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
4747 }
4748
4749 static inline
heat_est_d(const DOF_REAL_D_VEC * uh,const DOF_REAL_D_VEC * uh_old,ADAPT_INSTAT * adapt,REAL * (* rw_est)(EL *),REAL * (* rw_estc)(EL *),int quad_degree,REAL C[4],const void * A,MATENT_TYPE A_type,MATENT_TYPE A_blocktype,bool sym_grad,const BNDRY_FLAGS dirichlet_bndry,const REAL * (* f)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_DD grd_uh_gp,REAL time),FLAGS f_flags,const REAL * (* gn)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_D normal,REAL time),FLAGS gn_flags)4750 REAL heat_est_d(const DOF_REAL_D_VEC *uh,
4751 const DOF_REAL_D_VEC *uh_old,
4752 ADAPT_INSTAT *adapt,
4753 REAL *(*rw_est)(EL *),
4754 REAL *(*rw_estc)(EL *),
4755 int quad_degree,
4756 REAL C[4],
4757 const void *A,
4758 MATENT_TYPE A_type,
4759 MATENT_TYPE A_blocktype,
4760 bool sym_grad,
4761 const BNDRY_FLAGS dirichlet_bndry,
4762 const REAL *(*f)(REAL_D result,
4763 const EL_INFO *el_info,
4764 const QUAD *quad,
4765 int qp,
4766 const REAL_D uh_qp,
4767 const REAL_DD grd_uh_gp,
4768 REAL time),
4769 FLAGS f_flags,
4770 const REAL *(*gn)(REAL_D result,
4771 const EL_INFO *el_info,
4772 const QUAD *quad,
4773 int qp,
4774 const REAL_D uh_qp,
4775 const REAL_D normal,
4776 REAL time),
4777 FLAGS gn_flags)
4778 {
4779 return heat_est_dow(
4780 (const DOF_REAL_VEC_D *)uh, (const DOF_REAL_VEC_D *)uh_old,
4781 adapt, rw_est, rw_estc, quad_degree, C, A, A_type, A_blocktype,
4782 sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
4783 }
4784
4785 static inline
element_est_d(const EL_INFO * el_info,const void * est_handle)4786 REAL element_est_d(const EL_INFO *el_info, const void *est_handle)
4787 {
4788 return element_est_dow(el_info, est_handle);
4789 }
4790
4791 static inline
element_est_d_finish(const EL_INFO * el_info,REAL est_el,const void * est_handle)4792 void element_est_d_finish(const EL_INFO *el_info,
4793 REAL est_el, const void *est_handle)
4794 {
4795 element_est_dow_finish(el_info, est_el, est_handle);
4796 }
4797
4798 static inline
ellipt_est_d_finish(ADAPT_STAT * adapt,const void * est_handle)4799 REAL ellipt_est_d_finish(ADAPT_STAT *adapt, const void *est_handle)
4800 {
4801 return ellipt_est_dow_finish(adapt, est_handle);
4802 }
4803
4804 static inline
heat_est_d_finish(ADAPT_INSTAT * adapt,const void * est_handle)4805 REAL heat_est_d_finish(ADAPT_INSTAT *adapt, const void *est_handle)
4806 {
4807 return heat_est_dow_finish(adapt, est_handle);
4808 }
4809
4810 static inline
element_est_uh_d(const void * est_handle)4811 const REAL_D *element_est_uh_d(const void *est_handle)
4812 {
4813 return element_est_uh_dow(est_handle);
4814 }
4815
4816 static inline
element_est_grd_uh_d(const void * est_handle)4817 const REAL_DD *element_est_grd_uh_d(const void *est_handle) {
4818 return element_est_grd_uh_dow(est_handle);
4819 }
4820
4821 static inline
element_est_uh_loc_d(const void * est_handle)4822 const EL_REAL_D_VEC *element_est_uh_loc_d(const void *est_handle)
4823 {
4824 return (const EL_REAL_D_VEC *)element_est_uh_loc_dow(est_handle);
4825 }
4826
4827 static inline
element_est_uh_old_loc_d(const void * est_handle)4828 const EL_REAL_D_VEC *element_est_uh_old_loc_d(const void *est_handle)
4829 {
4830 return (const EL_REAL_D_VEC *)element_est_uh_old_loc_dow(est_handle);
4831 }
4832
4833 /*** file error.c *********************************************************/
4834 REAL max_err_at_qp(FCT_AT_X u, const DOF_REAL_VEC *, const QUAD *);
4835 REAL max_err_dow_at_qp(FCT_D_AT_X u, const DOF_REAL_VEC_D *, const QUAD *);
4836 REAL max_err_at_vert(FCT_AT_X f, const DOF_REAL_VEC *);
4837 REAL max_err_dow_at_vert(FCT_D_AT_X u, const DOF_REAL_VEC_D *uh);
4838
4839 REAL L2_err(FCT_AT_X u, const DOF_REAL_VEC *uh,
4840 const QUAD *quad,
4841 bool rel_err, bool mean_value_adjust,
4842 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
4843 REAL L2_err_dow(FCT_D_AT_X u, const DOF_REAL_VEC_D *uh,
4844 const QUAD *quad,
4845 bool rel_err, bool mean_value_adjust,
4846 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
4847 REAL L2_err_weighted(FCT_AT_X weight, FCT_AT_X u, const DOF_REAL_VEC *uh,
4848 const QUAD *quad,
4849 bool rel_err, bool mean_value_adjust,
4850 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
4851 REAL L2_err_dow_weighted(FCT_AT_X wieght, FCT_D_AT_X u,
4852 const DOF_REAL_VEC_D *uh,
4853 const QUAD *quad,
4854 bool rel_err, bool mean_value_adjust,
4855 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
4856
4857 REAL H1_err(GRD_FCT_AT_X grd_u, const DOF_REAL_VEC *u_h,
4858 const QUAD *quad, bool rel_err,
4859 REAL *(*rw_err_el)(EL *el), REAL *max_h1_err2);
4860 REAL H1_err_dow(GRD_FCT_D_AT_X grd_u, const DOF_REAL_VEC_D *uh,
4861 const QUAD *quad, bool rel_err,
4862 REAL *(*rw_err_esl)(EL *el), REAL *max_h1_err2);
4863 REAL deform_err(GRD_FCT_D_AT_X grd_u, const DOF_REAL_VEC_D *uh,
4864 const QUAD *quad,
4865 bool rel_err, REAL *(*rw_err_el)(EL *), REAL *max_el_err2);
4866 REAL H1_err_weighted(FCT_AT_X weight, GRD_FCT_AT_X grd_u,
4867 const DOF_REAL_VEC *u_h,
4868 const QUAD *quad, bool rel_err,
4869 REAL *(*rw_err_el)(EL *el), REAL *max_h1_err2);
4870 REAL H1_err_dow_weighted(FCT_AT_X weight, GRD_FCT_D_AT_X grd_u,
4871 const DOF_REAL_VEC_D *uh,
4872 const QUAD *quad, bool rel_err,
4873 REAL *(*rw_err_esl)(EL *el), REAL *max_h1_err2);
4874 REAL deform_err_weighted(FCT_AT_X weight, GRD_FCT_D_AT_X grd_u,
4875 const DOF_REAL_VEC_D *uh,
4876 const QUAD *quad,
4877 bool rel_err, REAL *(*rw_err_el)(EL *),
4878 REAL *max_el_err2);
4879
4880 REAL max_err_at_qp_loc(LOC_FCT_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4881 const DOF_REAL_VEC *uh, const QUAD *quad);
4882 REAL max_err_dow_at_qp_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4883 const DOF_REAL_VEC_D *uh, const QUAD *quad);
4884 REAL max_err_at_vert_loc(LOC_FCT_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4885 const DOF_REAL_VEC *uh);
4886 REAL max_err_dow_at_vert_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4887 const DOF_REAL_VEC_D *uh);
4888
4889 REAL L2_err_loc(LOC_FCT_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4890 const DOF_REAL_VEC *uh,
4891 const QUAD *quad,
4892 bool rel_err, bool mean_value_adjust,
4893 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
4894 REAL H1_err_loc(GRD_LOC_FCT_AT_QP grd_u_at_qp, void *ud, FLAGS fill_flag,
4895 const DOF_REAL_VEC *uh,
4896 const QUAD *quad, bool rel_err, REAL *(*rw_err_el)(EL *),
4897 REAL *max_h1_err2);
4898 REAL L2_err_loc_dow(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4899 const DOF_REAL_VEC_D *uh,
4900 const QUAD *quad,
4901 bool rel_err, bool mean_value_adjust,
4902 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
4903 REAL H1_err_loc_dow(GRD_LOC_FCT_D_AT_QP grd_u_at_qp, void *ud, FLAGS fill_flag,
4904 const DOF_REAL_VEC_D *uh, const QUAD *quad,
4905 bool rel_err,
4906 REAL *(*rw_err_el)(EL *), REAL *max_h1_err2);
4907 REAL deform_err_loc(GRD_LOC_FCT_D_AT_QP grd_u_loc, void *ud, FLAGS fill_flag,
4908 const DOF_REAL_VEC_D *uh, const QUAD *quad,
4909 bool rel_err,
4910 REAL *(*rw_err_el)(EL *), REAL *max_el_err2);
4911
4912 /* DOF_REAL_D_VEC versions */
4913
4914 static inline
max_err_d_at_qp(FCT_D_AT_X u,const DOF_REAL_D_VEC * uh,const QUAD * quad)4915 REAL max_err_d_at_qp(FCT_D_AT_X u, const DOF_REAL_D_VEC *uh, const QUAD *quad)
4916 {
4917 return max_err_dow_at_qp(u, (const DOF_REAL_VEC_D *)uh, quad);
4918 }
4919 static inline
max_err_d_at_vert(FCT_D_AT_X u,const DOF_REAL_D_VEC * uh)4920 REAL max_err_d_at_vert(FCT_D_AT_X u, const DOF_REAL_D_VEC *uh)
4921 {
4922 return max_err_dow_at_vert(u, (const DOF_REAL_VEC_D *)uh);
4923 }
4924 static inline
L2_err_d(FCT_D_AT_X u,const DOF_REAL_D_VEC * uh,const QUAD * quad,bool rel_err,bool mean_value_adjust,REAL * (* rw_err_el)(EL * el),REAL * max_l2_err2)4925 REAL L2_err_d(FCT_D_AT_X u, const DOF_REAL_D_VEC *uh,
4926 const QUAD *quad,
4927 bool rel_err, bool mean_value_adjust,
4928 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2)
4929 {
4930 return L2_err_dow(u, (const DOF_REAL_VEC_D *)uh,
4931 quad, rel_err, mean_value_adjust, rw_err_el, max_l2_err2);
4932 }
4933 static inline
H1_err_d(GRD_FCT_D_AT_X grd_u,const DOF_REAL_D_VEC * uh,const QUAD * quad,bool rel_err,REAL * (* rw_err_el)(EL * el),REAL * max_h1_err2)4934 REAL H1_err_d(GRD_FCT_D_AT_X grd_u, const DOF_REAL_D_VEC *uh,
4935 const QUAD *quad, bool rel_err,
4936 REAL *(*rw_err_el)(EL *el), REAL *max_h1_err2)
4937 {
4938 return H1_err_dow(grd_u, (const DOF_REAL_VEC_D *)uh,
4939 quad, rel_err, rw_err_el, max_h1_err2);
4940 }
4941 static inline
max_err_d_at_qp_loc(LOC_FCT_D_AT_QP u_at_qp,void * ud,FLAGS fill_flag,const DOF_REAL_D_VEC * uh,const QUAD * quad)4942 REAL max_err_d_at_qp_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4943 const DOF_REAL_D_VEC *uh, const QUAD *quad)
4944 {
4945 return max_err_dow_at_qp_loc(u_at_qp, ud, fill_flag,
4946 (const DOF_REAL_VEC_D *)uh, quad);
4947 }
4948 static inline
max_err_d_at_vert_loc(LOC_FCT_D_AT_QP u_at_qp,void * ud,FLAGS fill_flag,const DOF_REAL_D_VEC * uh)4949 REAL max_err_d_at_vert_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4950 const DOF_REAL_D_VEC *uh)
4951 {
4952 return max_err_dow_at_vert_loc(u_at_qp, ud, fill_flag,
4953 (const DOF_REAL_VEC_D *)uh);
4954 }
4955 static inline
L2_err_loc_d(LOC_FCT_D_AT_QP u_at_qp,void * ud,FLAGS fill_flag,const DOF_REAL_D_VEC * uh,const QUAD * quad,bool rel_err,bool mean_value_adjust,REAL * (* rw_err_el)(EL * el),REAL * max_l2_err2)4956 REAL L2_err_loc_d(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
4957 const DOF_REAL_D_VEC *uh,
4958 const QUAD *quad,
4959 bool rel_err, bool mean_value_adjust,
4960 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2)
4961 {
4962 return L2_err_loc_dow(u_at_qp, ud, fill_flag,
4963 (const DOF_REAL_VEC_D *)uh,
4964 quad,
4965 rel_err, mean_value_adjust, rw_err_el, max_l2_err2);
4966 }
4967 static inline
H1_err_loc_d(GRD_LOC_FCT_D_AT_QP grd_u_at_qp,void * ud,FLAGS fill_flag,const DOF_REAL_D_VEC * uh,const QUAD * quad,bool rel_err,REAL * (* rw_err_el)(EL *),REAL * max_h1_err2)4968 REAL H1_err_loc_d(GRD_LOC_FCT_D_AT_QP grd_u_at_qp, void *ud, FLAGS fill_flag,
4969 const DOF_REAL_D_VEC *uh, const QUAD *quad,
4970 bool rel_err,
4971 REAL *(*rw_err_el)(EL *), REAL *max_h1_err2)
4972 {
4973 return H1_err_loc_dow(grd_u_at_qp, ud, fill_flag,
4974 (const DOF_REAL_VEC_D *)uh,
4975 quad, rel_err, rw_err_el, max_h1_err2);
4976 }
4977
4978 /*** file eval.c **********************************************************/
4979
4980 /* evaluation routines are defined as inline functions in the file evaluate.h */
4981
4982 REAL H1_norm_uh(const QUAD *quad, const DOF_REAL_VEC *u_h);
4983 REAL L2_norm_uh(const QUAD *quad, const DOF_REAL_VEC *u_h);
4984 REAL L8_uh_at_qp(REAL *minp, REAL *maxp,
4985 const QUAD *quad, const DOF_REAL_VEC *u_h);
4986 REAL H1_norm_uh_dow(const QUAD *quad, const DOF_REAL_VEC_D *u_h);
4987 REAL L2_norm_uh_dow(const QUAD *quad, const DOF_REAL_VEC_D *u_h);
4988 REAL L8_uh_at_qp_dow(REAL *minp, REAL *maxp,
4989 const QUAD *quad, const DOF_REAL_VEC_D *u_h);
L2_norm_uh_d(const QUAD * quad,const DOF_REAL_D_VEC * u_h)4990 static inline REAL L2_norm_uh_d(const QUAD *quad, const DOF_REAL_D_VEC *u_h)
4991 {
4992 return L2_norm_uh_dow(quad, (const DOF_REAL_VEC_D *)u_h);
4993 }
L8_uh_at_qp_d(REAL * minp,REAL * maxp,const QUAD * quad,const DOF_REAL_D_VEC * u_h)4994 static inline REAL L8_uh_at_qp_d(REAL *minp, REAL *maxp,
4995 const QUAD *quad, const DOF_REAL_D_VEC *u_h)
4996 {
4997 return L8_uh_at_qp_dow(minp, maxp, quad, (const DOF_REAL_VEC_D *)u_h);
4998 }
H1_norm_uh_d(const QUAD * quad,const DOF_REAL_D_VEC * u_h)4999 static inline REAL H1_norm_uh_d(const QUAD *quad, const DOF_REAL_D_VEC *u_h)
5000 {
5001 return H1_norm_uh_dow(quad, (const DOF_REAL_VEC_D *)u_h);
5002 }
5003
5004 void interpol(FCT_AT_X fct, DOF_REAL_VEC *u_h);
5005 void interpol_dow(FCT_D_AT_X fct, DOF_REAL_VEC_D *uh);
5006 void interpol_loc(DOF_REAL_VEC *vec,
5007 LOC_FCT_AT_QP fct_at_qp, void *app_data,
5008 FLAGS fill_flags);
5009 void interpol_loc_dow(DOF_REAL_VEC_D *vec,
5010 LOC_FCT_D_AT_QP fct_at_qp, void *app_data,
5011 FLAGS fill_flags);
5012
5013 static inline
interpol_d(FCT_D_AT_X f,DOF_REAL_D_VEC * u_h)5014 void interpol_d(FCT_D_AT_X f, DOF_REAL_D_VEC *u_h)
5015 {
5016 interpol_dow(f, (DOF_REAL_VEC_D *)u_h);
5017 }
5018
5019 static inline
interpol_loc_d(DOF_REAL_D_VEC * vec,LOC_FCT_D_AT_QP fct_at_qp,void * app_data,FLAGS fill_flags)5020 void interpol_loc_d(DOF_REAL_D_VEC *vec,
5021 LOC_FCT_D_AT_QP fct_at_qp, void *app_data,
5022 FLAGS fill_flags)
5023 {
5024 interpol_loc_dow((DOF_REAL_VEC_D *)vec, fct_at_qp, app_data, fill_flags);
5025 }
5026
5027 /*** file graphXO.c *******************************************************/
5028 GRAPH_WINDOW graph_open_window(const char *title, const char *geometry,
5029 REAL *world, MESH *mesh);
5030 void graph_close_window(GRAPH_WINDOW win);
5031 void graph_clear_window(GRAPH_WINDOW win, const GRAPH_RGBCOLOR c);
5032
5033 void graph_mesh(GRAPH_WINDOW win, MESH *mesh, const GRAPH_RGBCOLOR c,
5034 FLAGS flag);
5035 void graph_drv(GRAPH_WINDOW win, const DOF_REAL_VEC *uh,
5036 REAL min, REAL max, int refine);
5037 void graph_drv_d(GRAPH_WINDOW win, const DOF_REAL_D_VEC *uh,
5038 REAL min, REAL max, int refine);
5039 void graph_el_est(GRAPH_WINDOW win, MESH *mesh, REAL (*get_el_est)(EL *el),
5040 REAL min, REAL max);
5041
5042 void graph_point(GRAPH_WINDOW, const REAL [2], const GRAPH_RGBCOLOR, float);
5043 void graph_points(GRAPH_WINDOW win, int np, REAL (*p)[2],
5044 const GRAPH_RGBCOLOR c, float ps);
5045 void graph_line(GRAPH_WINDOW, const REAL [2], const REAL [2],
5046 const GRAPH_RGBCOLOR, float);
5047
5048 void graph_fvalues_2d(GRAPH_WINDOW win, MESH *mesh,
5049 REAL(*fct)(const EL_INFO *el_info, const REAL *lambda),
5050 FLAGS flags, REAL min, REAL max, int refine);
5051 void graph_level_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v, REAL level,
5052 const GRAPH_RGBCOLOR c, int refine);
5053 void graph_levels_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v,
5054 int n, REAL const *levels, const GRAPH_RGBCOLOR *color,
5055 int refine);
5056 void graph_level_d_2d(GRAPH_WINDOW, const DOF_REAL_D_VEC *,
5057 REAL, const GRAPH_RGBCOLOR, int);
5058 void graph_levels_d_2d(GRAPH_WINDOW, const DOF_REAL_D_VEC *,
5059 int, const REAL *, const GRAPH_RGBCOLOR *, int);
5060
5061 /* multigrid level display routines: */
5062 void graph_mesh_mg_2d(GRAPH_WINDOW win, MESH *mesh, const GRAPH_RGBCOLOR c,
5063 FLAGS flags, int mg_level);
5064 void graph_values_mg_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v,
5065 REAL min, REAL max, int refine,
5066 int mg_level, const FE_SPACE *fe_space,
5067 const int *sort_dof_invers);
5068
5069 /*** file l2scp.c *********************************************************/
5070 void L2scp_fct_bas(FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh);
5071 void L2scp_fct_bas_dow(FCT_D_AT_X, const QUAD *quad, DOF_REAL_VEC_D *fhd);
5072 void L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
5073 LOC_FCT_AT_QP f_at_qp, void *fct_data, FLAGS fill_flag,
5074 const QUAD *quad);
5075 void L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
5076 LOC_FCT_D_AT_QP f_at_qp, void *ud, FLAGS fill_flag,
5077 const QUAD *quad);
5078
5079 void H1scp_fct_bas(GRD_FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh);
5080 void H1scp_fct_bas_dow(GRD_FCT_D_AT_X f, const QUAD *quad, DOF_REAL_VEC_D *fh);
5081 void H1scp_fct_bas_loc(DOF_REAL_VEC *fh,
5082 GRD_LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
5083 const QUAD *quad);
5084 void H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
5085 GRD_LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
5086 const QUAD *quad);
5087
5088 bool bndry_L2scp_fct_bas(DOF_REAL_VEC *fh,
5089 REAL (*f)(const REAL_D x, const REAL_D normal),
5090 const BNDRY_FLAGS bndry_seg,
5091 const WALL_QUAD *quad);
5092 bool bndry_L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
5093 LOC_FCT_AT_QP f_at_qp, void *ud, FLAGS fill_flag,
5094 const BNDRY_FLAGS bndry_seg,
5095 const WALL_QUAD *quad);
5096 bool bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
5097 const REAL *(*f)(const REAL_D x,
5098 const REAL_D normal,
5099 REAL_D result),
5100 const BNDRY_FLAGS bndry_seg,
5101 const WALL_QUAD *quad);
5102 bool bndry_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
5103 LOC_FCT_D_AT_QP f_at_qp, void *ud,
5104 FLAGS fill_flag,
5105 const BNDRY_FLAGS bndry_seg,
5106 const WALL_QUAD *quad);
5107
5108 void trace_L2scp_fct_bas(DOF_REAL_VEC *fh, FCT_AT_X f,
5109 MESH *trace_mesh, const QUAD *quad);
5110 void trace_L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
5111 LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
5112 MESH *trace_mesh,
5113 const QUAD *quad);
5114 void trace_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
5115 FCT_D_AT_X f,
5116 MESH *trace_mesh,
5117 const QUAD *quad);
5118 void trace_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
5119 LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
5120 MESH *trace_mesh,
5121 const QUAD *quad);
5122
5123 bool bndry_H1scp_fct_bas(DOF_REAL_VEC *fh,
5124 const REAL *(*f)(REAL_D result,
5125 const REAL_D x, const REAL_D normal),
5126 const BNDRY_FLAGS scp_segment,
5127 const WALL_QUAD *quad);
5128 bool bndry_H1scp_fct_bas_loc(DOF_REAL_VEC *fh,
5129 GRD_LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
5130 const BNDRY_FLAGS scp_segment,
5131 const WALL_QUAD *quad);
5132 bool bndry_H1scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
5133 const REAL_D *(*f)(REAL_DD result,
5134 const REAL_D x,
5135 const REAL_D normal),
5136 const BNDRY_FLAGS bndry_seg,
5137 const WALL_QUAD *quad);
5138 bool bndry_H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
5139 GRD_LOC_FCT_D_AT_QP f_loc,
5140 void *fd, FLAGS fill_flag,
5141 const BNDRY_FLAGS bndry_seg,
5142 const WALL_QUAD *quad);
5143
5144 extern bool
5145 dirichlet_bound(DOF_REAL_VEC *fh, DOF_REAL_VEC *uh,
5146 DOF_SCHAR_VEC *bound,
5147 const BNDRY_FLAGS dirichlet_segment,
5148 REAL (*g)(const REAL_D));
5149 extern bool
5150 dirichlet_bound_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
5151 DOF_SCHAR_VEC *bound,
5152 const BNDRY_FLAGS dirichlet_segment,
5153 const REAL *(*g)(const REAL_D, REAL_D));
5154 extern bool
5155 dirichlet_bound_loc(DOF_REAL_VEC *fh,
5156 DOF_REAL_VEC *uh,
5157 DOF_SCHAR_VEC *bound,
5158 const BNDRY_FLAGS dirichlet_segment,
5159 LOC_FCT_AT_QP g_at_qp, void *ud, FLAGS fill_flags);
5160 extern bool
5161 dirichlet_bound_loc_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
5162 DOF_SCHAR_VEC *bound,
5163 const BNDRY_FLAGS dirichlet_segment,
5164 LOC_FCT_D_AT_QP g_at_qp, void *ud, FLAGS fill_flags);
5165
5166 REAL mean_value(MESH *mesh,
5167 REAL (*f)(const REAL_D),
5168 const DOF_REAL_VEC *fh, const QUAD *quad);
5169 const REAL *mean_value_dow(MESH *mesh,
5170 FCT_D_AT_X f,
5171 const DOF_REAL_VEC_D *fh,
5172 const QUAD *quad,
5173 REAL_D mean);
5174 REAL mean_value_loc(MESH *mesh,
5175 LOC_FCT_AT_QP f_at_qp, void *ud, FLAGS fill_flags,
5176 const DOF_REAL_VEC *fh, const QUAD *quad);
5177 const REAL *mean_value_loc_dow(REAL_D mean,
5178 MESH *mesh,
5179 LOC_FCT_D_AT_QP f_at_qp,
5180 void *ud, FLAGS fill_flag,
5181 const DOF_REAL_VEC_D *fh,
5182 const QUAD *quad);
5183
5184 void robin_bound_matrix_info(EL_MATRIX_INFO *robin_info,
5185 const FE_SPACE *row_fe_space,
5186 const FE_SPACE *col_fe_space,
5187 const BNDRY_FLAGS robin_segment,
5188 REAL alpha_r,
5189 const WALL_QUAD *wall_quad,
5190 REAL exponent);
5191 void robin_bound(DOF_MATRIX *matrix,
5192 const BNDRY_FLAGS robin_seg,
5193 REAL alpha_r,
5194 const WALL_QUAD *wall_quad,
5195 REAL exponent);
5196
5197 bool boundary_conditions(DOF_MATRIX *matrix,
5198 DOF_REAL_VEC *fh,
5199 DOF_REAL_VEC *uh,
5200 DOF_SCHAR_VEC *bound,
5201 const BNDRY_FLAGS dirichlet_segment,
5202 REAL (*g)(const REAL_D x),
5203 REAL (*gn)(const REAL_D x, const REAL_D normal),
5204 REAL alpha_r,
5205 const WALL_QUAD *wall_quad);
5206 bool boundary_conditions_loc(DOF_MATRIX *matrix,
5207 DOF_REAL_VEC *fh,
5208 DOF_REAL_VEC *uh,
5209 DOF_SCHAR_VEC *bound,
5210 const BNDRY_FLAGS dirichlet_segment,
5211 LOC_FCT_AT_QP g_at_qp,
5212 LOC_FCT_AT_QP gn_at_qp,
5213 void *app_data, FLAGS fill_flags,
5214 REAL alpha_r,
5215 const WALL_QUAD *wall_quad);
5216 bool boundary_conditions_dow(DOF_MATRIX *matrix,
5217 DOF_REAL_VEC_D *fh,
5218 DOF_REAL_VEC_D *uh,
5219 DOF_SCHAR_VEC *bound,
5220 const BNDRY_FLAGS dirichlet_segment,
5221 const REAL *(*g)(const REAL_D x, REAL_D res),
5222 const REAL *(*gn)(const REAL_D x,
5223 const REAL_D normal,
5224 REAL_D res),
5225 REAL alpha_r,
5226 const WALL_QUAD *wall_quad);
5227 bool boundary_conditions_loc_dow(DOF_MATRIX *matrix,
5228 DOF_REAL_VEC_D *fh,
5229 DOF_REAL_VEC_D *uh,
5230 DOF_SCHAR_VEC *bound,
5231 const BNDRY_FLAGS dirichlet_segment,
5232 LOC_FCT_D_AT_QP g_at_qp,
5233 LOC_FCT_D_AT_QP gn_at_qp,
5234 void *app_data, FLAGS fill_flags,
5235 REAL alpha_r,
5236 const WALL_QUAD *wall_quad);
5237
5238 /* ... _d versions */
5239 static inline
L2scp_fct_bas_d(FCT_D_AT_X fct,const QUAD * quad,DOF_REAL_D_VEC * fhd)5240 void L2scp_fct_bas_d(FCT_D_AT_X fct, const QUAD *quad, DOF_REAL_D_VEC *fhd)
5241 {
5242 L2scp_fct_bas_dow(fct, quad, (DOF_REAL_VEC_D *)fhd);
5243 }
5244 static inline
L2scp_fct_bas_loc_d(DOF_REAL_D_VEC * fh,LOC_FCT_D_AT_QP f_at_qp,void * fd,FLAGS fill_flag,const QUAD * quad)5245 void L2scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh,
5246 LOC_FCT_D_AT_QP f_at_qp, void *fd, FLAGS fill_flag,
5247 const QUAD *quad)
5248 {
5249 L2scp_fct_bas_loc_dow((DOF_REAL_VEC_D *)fh, f_at_qp, fd, fill_flag, quad);
5250 }
5251 static inline
H1scp_fct_bas_d(GRD_FCT_D_AT_X f,const QUAD * quad,DOF_REAL_D_VEC * fh)5252 void H1scp_fct_bas_d(GRD_FCT_D_AT_X f, const QUAD *quad, DOF_REAL_D_VEC *fh)
5253 {
5254 H1scp_fct_bas_dow(f, quad, (DOF_REAL_VEC_D *)fh);
5255 }
5256 static inline
H1scp_fct_bas_loc_d(DOF_REAL_D_VEC * fh,GRD_LOC_FCT_D_AT_QP f,void * fd,FLAGS fill_flag,const QUAD * quad)5257 void H1scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh,
5258 GRD_LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
5259 const QUAD *quad)
5260 {
5261 H1scp_fct_bas_loc_dow((DOF_REAL_VEC_D *)fh, f, fd, fill_flag, quad);
5262 }
5263 static inline
bndry_L2scp_fct_bas_d(DOF_REAL_D_VEC * fh,const REAL * (* f)(const REAL_D x,const REAL_D normal,REAL_D result),const BNDRY_FLAGS bndry_seg,const WALL_QUAD * quad)5264 bool bndry_L2scp_fct_bas_d(DOF_REAL_D_VEC *fh,
5265 const REAL *(*f)(const REAL_D x,
5266 const REAL_D normal,
5267 REAL_D result),
5268 const BNDRY_FLAGS bndry_seg,
5269 const WALL_QUAD *quad)
5270 {
5271 return bndry_L2scp_fct_bas_dow((DOF_REAL_VEC_D *)fh, f, bndry_seg, quad);
5272 }
5273 static inline
bndry_L2scp_fct_bas_loc_d(DOF_REAL_D_VEC * fh,LOC_FCT_D_AT_QP f_at_qp,void * fd,FLAGS fill_flag,const BNDRY_FLAGS bndry_seg,const WALL_QUAD * quad)5274 bool bndry_L2scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh,
5275 LOC_FCT_D_AT_QP f_at_qp, void *fd,
5276 FLAGS fill_flag,
5277 const BNDRY_FLAGS bndry_seg,
5278 const WALL_QUAD *quad)
5279 {
5280 return bndry_L2scp_fct_bas_loc_dow((DOF_REAL_VEC_D *)fh,
5281 f_at_qp, fd, fill_flag, bndry_seg, quad);
5282 }
5283 static inline
dirichlet_bound_d(DOF_REAL_D_VEC * fh,DOF_REAL_D_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,const REAL * (* g)(const REAL_D,REAL_D))5284 bool dirichlet_bound_d(DOF_REAL_D_VEC *fh, DOF_REAL_D_VEC *uh,
5285 DOF_SCHAR_VEC *bound,
5286 const BNDRY_FLAGS dirichlet_segment,
5287 const REAL *(*g)(const REAL_D, REAL_D))
5288 {
5289 return dirichlet_bound_dow((DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
5290 bound, dirichlet_segment, g);
5291 }
5292 static inline
dirichlet_bound_loc_d(DOF_REAL_D_VEC * fh,DOF_REAL_D_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,LOC_FCT_D_AT_QP g_at_qp,void * gd,FLAGS fill_flags)5293 bool dirichlet_bound_loc_d(DOF_REAL_D_VEC *fh, DOF_REAL_D_VEC *uh,
5294 DOF_SCHAR_VEC *bound,
5295 const BNDRY_FLAGS dirichlet_segment,
5296 LOC_FCT_D_AT_QP g_at_qp, void *gd, FLAGS fill_flags)
5297 {
5298 return dirichlet_bound_loc_dow((DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
5299 bound, dirichlet_segment,
5300 g_at_qp, gd, fill_flags);
5301 }
5302 static inline
mean_value_d(MESH * mesh,FCT_D_AT_X f,const DOF_REAL_D_VEC * fh,const QUAD * quad,REAL_D mean)5303 const REAL *mean_value_d(MESH *mesh,
5304 FCT_D_AT_X f,
5305 const DOF_REAL_D_VEC *fh,
5306 const QUAD *quad,
5307 REAL_D mean)
5308 {
5309 return mean_value_dow(mesh, f, (const DOF_REAL_VEC_D *)fh, quad, mean);
5310 }
5311 static inline
mean_value_loc_d(REAL_D mean,MESH * mesh,LOC_FCT_D_AT_QP f_at_qp,void * fd,FLAGS fill_flag,const DOF_REAL_D_VEC * fh,const QUAD * quad)5312 const REAL *mean_value_loc_d(REAL_D mean,
5313 MESH *mesh,
5314 LOC_FCT_D_AT_QP f_at_qp, void *fd, FLAGS fill_flag,
5315 const DOF_REAL_D_VEC *fh,
5316 const QUAD *quad)
5317 {
5318 return mean_value_loc_dow(
5319 mean, mesh, f_at_qp, fd, fill_flag, (const DOF_REAL_VEC_D *)fh, quad);
5320 }
5321 static inline
boundary_conditions_d(DOF_MATRIX * matrix,DOF_REAL_D_VEC * fh,DOF_REAL_D_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,const REAL * (* g)(const REAL_D x,REAL_D res),const REAL * (* gn)(const REAL_D x,const REAL_D normal,REAL_D res),REAL alpha_r,const WALL_QUAD * wall_quad)5322 bool boundary_conditions_d(DOF_MATRIX *matrix,
5323 DOF_REAL_D_VEC *fh,
5324 DOF_REAL_D_VEC *uh,
5325 DOF_SCHAR_VEC *bound,
5326 const BNDRY_FLAGS dirichlet_segment,
5327 const REAL *(*g)(const REAL_D x, REAL_D res),
5328 const REAL *(*gn)(const REAL_D x,
5329 const REAL_D normal,
5330 REAL_D res),
5331 REAL alpha_r,
5332 const WALL_QUAD *wall_quad)
5333 {
5334 return boundary_conditions_dow(matrix,
5335 (DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
5336 bound, dirichlet_segment, g, gn, alpha_r,
5337 wall_quad);
5338 }
5339 static inline
boundary_conditions_loc_d(DOF_MATRIX * matrix,DOF_REAL_D_VEC * fh,DOF_REAL_D_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,LOC_FCT_D_AT_QP g_at_qp,LOC_FCT_D_AT_QP gn_at_qp,void * gdata,FLAGS fill_flags,REAL alpha_r,const WALL_QUAD * wall_quad)5340 bool boundary_conditions_loc_d(DOF_MATRIX *matrix,
5341 DOF_REAL_D_VEC *fh,
5342 DOF_REAL_D_VEC *uh,
5343 DOF_SCHAR_VEC *bound,
5344 const BNDRY_FLAGS dirichlet_segment,
5345 LOC_FCT_D_AT_QP g_at_qp,
5346 LOC_FCT_D_AT_QP gn_at_qp,
5347 void *gdata, FLAGS fill_flags,
5348 REAL alpha_r,
5349 const WALL_QUAD *wall_quad)
5350 {
5351 return boundary_conditions_loc_dow(matrix,
5352 (DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
5353 bound, dirichlet_segment,
5354 g_at_qp, gn_at_qp, gdata, fill_flags,
5355 alpha_r, wall_quad);
5356 }
5357
5358 /* file oem_solve.c *******************************************************/
5359 extern OEM_MV_FCT get_oem_solver(OEM_SOLVER solver);
5360 extern OEM_DATA *init_oem_solve(const DOF_MATRIX *A,
5361 const DOF_SCHAR_VEC *bound,
5362 REAL tol, const PRECON *precon,
5363 int restart, int max_iter, int info);
5364 extern const PRECON *init_oem_precon(const DOF_MATRIX *A,
5365 const DOF_SCHAR_VEC *mask,
5366 int info, OEM_PRECON precon,
5367 ... /* ssor_omega, ssor_n_iter etc. */);
5368 extern const PRECON *vinit_oem_precon(const DOF_MATRIX *A,
5369 const DOF_SCHAR_VEC *mask,
5370 int info, OEM_PRECON,
5371 va_list ap);
5372 const PRECON *init_precon_from_type(const DOF_MATRIX *A,
5373 const DOF_SCHAR_VEC *mask,
5374 int info,
5375 const PRECON_TYPE *prec_type);
5376
5377 extern void release_oem_solve(const OEM_DATA *oem);
5378 extern int oem_mat_vec(void *ud, int dim, const REAL *x, REAL *y);
5379 extern OEM_MV_FCT
5380 init_oem_mat_vec(void **datap,
5381 MatrixTranspose transpose, const DOF_MATRIX *A,
5382 const DOF_SCHAR_VEC *bound);
5383 extern void exit_oem_mat_vec(void *);
5384
5385 extern int call_oem_solve_s(const OEM_DATA *oem, OEM_SOLVER solver,
5386 const DOF_REAL_VEC *f, DOF_REAL_VEC *u);
5387 extern int
5388 oem_solve_s(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
5389 const DOF_REAL_VEC *f, DOF_REAL_VEC *u, OEM_SOLVER solver,
5390 REAL tol, const PRECON *precon,
5391 int restart, int max_iter, int info);
5392
5393 extern int call_oem_solve_dow(const OEM_DATA *oem, OEM_SOLVER solver,
5394 const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u);
5395 extern int
5396 oem_solve_dow(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
5397 const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u, OEM_SOLVER solver,
5398 REAL tol, const PRECON *precon,
5399 int restart, int max_iter, int info);
5400
5401 static inline
call_oem_solve_d(const OEM_DATA * oem,OEM_SOLVER solver,const DOF_REAL_D_VEC * f,DOF_REAL_D_VEC * u)5402 int call_oem_solve_d(const OEM_DATA *oem, OEM_SOLVER solver,
5403 const DOF_REAL_D_VEC *f, DOF_REAL_D_VEC *u)
5404 {
5405 return call_oem_solve_dow(oem, solver,
5406 (const DOF_REAL_VEC_D *)f, (DOF_REAL_VEC_D *)u);
5407 }
5408 static inline
oem_solve_d(const DOF_MATRIX * A,const DOF_SCHAR_VEC * bound,const DOF_REAL_D_VEC * f,DOF_REAL_D_VEC * u,OEM_SOLVER solver,REAL tol,const PRECON * precon,int restart,int max_iter,int info)5409 int oem_solve_d(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
5410 const DOF_REAL_D_VEC *f, DOF_REAL_D_VEC *u, OEM_SOLVER solver,
5411 REAL tol, const PRECON *precon,
5412 int restart, int max_iter, int info)
5413 {
5414 return oem_solve_dow(A, bound,
5415 (const DOF_REAL_VEC_D *)f, (DOF_REAL_VEC_D *)u, solver,
5416 tol, precon, restart, max_iter, info);
5417 }
5418
5419 /* file oem_sp_solve.c ******************************************************/
5420 typedef struct sp_constraint
5421 {
5422 const DOF_MATRIX *B, *Bt;
5423 const DOF_SCHAR_VEC *bound;
5424 OEM_MV_FCT project;
5425 void *project_data; /* possibly "OEM_DATA *" */
5426 OEM_MV_FCT precon;
5427 void *precon_data; /* possibly "OEM_DATA *" */
5428 REAL proj_factor, prec_factor;
5429 } SP_CONSTRAINT;
5430
5431 SP_CONSTRAINT *init_sp_constraint(const DOF_MATRIX *B,
5432 const DOF_MATRIX *Bt,
5433 const DOF_SCHAR_VEC *bound,
5434 REAL tol, int info,
5435 const DOF_MATRIX *Yproj,
5436 OEM_SOLVER Yproj_solver,
5437 int Yproj_max_iter, const PRECON *Yproj_prec,
5438 const DOF_MATRIX *Yprec,
5439 OEM_SOLVER Yprec_solver,
5440 int Yprec_max_iter, const PRECON *Yprec_prec,
5441 REAL Yproj_frac, REAL Yprec_frac);
5442 void release_sp_constraint(SP_CONSTRAINT *constr);
5443 int oem_sp_schur_solve(OEM_SOLVER sp_solver,
5444 REAL sp_tol, int sp_max_iter, int sp_info,
5445 OEM_MV_FCT principal_inverse,
5446 OEM_DATA *principal_data,
5447 const DOF_REAL_VEC_D *f,
5448 DOF_REAL_VEC_D *u,
5449 SP_CONSTRAINT *constraint,
5450 const DOF_REAL_VEC *g,
5451 DOF_REAL_VEC *p,
5452 ...);
5453 int oem_sp_solve_dow_scl(OEM_SOLVER sp_solver,
5454 REAL sp_tol, REAL tol_incr,
5455 int sp_max_iter, int sp_info,
5456 const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
5457 OEM_SOLVER A_solver,
5458 int A_max_iter, const PRECON *A_precon,
5459 DOF_MATRIX *B,
5460 DOF_MATRIX *Bt,
5461 DOF_MATRIX *Yproj,
5462 OEM_SOLVER Yproj_solver,
5463 int Yproj_max_iter, const PRECON *Yproj_precon,
5464 DOF_MATRIX *Yprec,
5465 OEM_SOLVER Yprec_solver,
5466 int Yprec_max_iter, const PRECON *Yprec_precon,
5467 REAL Yproj_frac, REAL Ypre_frac,
5468 const DOF_REAL_VEC_D *f,
5469 const DOF_REAL_VEC *g,
5470 DOF_REAL_VEC_D *x,
5471 DOF_REAL_VEC *y);
5472 REAL sp_dirichlet_bound_dow_scl(MatrixTranspose transpose,
5473 const DOF_MATRIX *Bt,
5474 const DOF_SCHAR_VEC *bound,
5475 const DOF_REAL_VEC_D *u_h,
5476 DOF_REAL_VEC *g_h);
5477 REAL sp_flux_adjust_dow_scl(MatrixTranspose transpose,
5478 const DOF_MATRIX *Bt,
5479 const DOF_SCHAR_VEC *bound,
5480 const DOF_REAL_VEC_D *u_h,
5481 DOF_REAL_VEC *g_h,
5482 REAL initial_flux,
5483 bool do_adjust);
5484
5485 static inline
oem_sp_solve_ds(OEM_SOLVER sp_solver,REAL sp_tol,REAL tol_incr,int sp_max_iter,int sp_info,const DOF_MATRIX * A,const DOF_SCHAR_VEC * bound,OEM_SOLVER A_solver,int A_max_iter,const PRECON * A_prec,DOF_MATRIX * B,DOF_MATRIX * Bt,DOF_MATRIX * Yproj,OEM_SOLVER Yproj_solver,int Yproj_max_iter,const PRECON * Yproj_prec,DOF_MATRIX * Yprec,OEM_SOLVER Yprec_solver,int Yprec_max_iter,const PRECON * Yprec_prec,REAL Yproj_frac,REAL Yprec_frac,const DOF_REAL_D_VEC * f,const DOF_REAL_VEC * g,DOF_REAL_D_VEC * x,DOF_REAL_VEC * y)5486 int oem_sp_solve_ds(
5487 OEM_SOLVER sp_solver,
5488 REAL sp_tol, REAL tol_incr,
5489 int sp_max_iter, int sp_info,
5490 const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
5491 OEM_SOLVER A_solver, int A_max_iter, const PRECON *A_prec,
5492 DOF_MATRIX *B,
5493 DOF_MATRIX *Bt,
5494 DOF_MATRIX *Yproj,
5495 OEM_SOLVER Yproj_solver, int Yproj_max_iter, const PRECON *Yproj_prec,
5496 DOF_MATRIX *Yprec,
5497 OEM_SOLVER Yprec_solver, int Yprec_max_iter, const PRECON *Yprec_prec,
5498 REAL Yproj_frac, REAL Yprec_frac,
5499 const DOF_REAL_D_VEC *f,
5500 const DOF_REAL_VEC *g,
5501 DOF_REAL_D_VEC *x,
5502 DOF_REAL_VEC *y)
5503 {
5504 return oem_sp_solve_dow_scl(sp_solver, sp_tol, tol_incr, sp_max_iter, sp_info,
5505 A, bound, A_solver, A_max_iter, A_prec,
5506 B, Bt,
5507 Yproj, Yproj_solver, Yproj_max_iter, Yproj_prec,
5508 Yprec, Yprec_solver, Yprec_max_iter, Yprec_prec,
5509 Yproj_frac, Yprec_frac,
5510 (const DOF_REAL_VEC_D *)f, g,
5511 (DOF_REAL_VEC_D *)x, y);
5512 }
5513 static inline
sp_dirichlet_bound_ds(MatrixTranspose transpose,const DOF_MATRIX * Bt,const DOF_SCHAR_VEC * bound,const DOF_REAL_D_VEC * u_h,DOF_REAL_VEC * g_h)5514 REAL sp_dirichlet_bound_ds(MatrixTranspose transpose,
5515 const DOF_MATRIX *Bt,
5516 const DOF_SCHAR_VEC *bound,
5517 const DOF_REAL_D_VEC *u_h,
5518 DOF_REAL_VEC *g_h)
5519 {
5520 return sp_dirichlet_bound_dow_scl(transpose, Bt, bound,
5521 (const DOF_REAL_VEC_D *)u_h, g_h);
5522 }
5523
5524
5525 /* file parametric.c ********************************************************/
5526 void use_lagrange_parametric(MESH *mesh, int degree,
5527 NODE_PROJ *n_proj,
5528 FLAGS flags);
5529 DOF_REAL_D_VEC *get_lagrange_coords(MESH *mesh);
5530 DOF_PTR_VEC *get_lagrange_edge_projections(MESH *mesh);
5531
5532 typedef enum param_copy_direction {
5533 COPY_FROM_MESH = false,
5534 COPY_TO_MESH = true
5535 } PARAM_COPY_DIRECTION;
5536
5537 void copy_lagrange_coords(MESH *mesh, DOF_REAL_D_VEC *coords, bool tomesh);
5538
5539 /*-- file sor.c *************************************************************/
5540 int sor_d(DOF_MATRIX *a, const DOF_REAL_D_VEC *f, const DOF_SCHAR_VEC *b,
5541 DOF_REAL_D_VEC *u, REAL omega, REAL tol, int max_iter, int info);
5542 int sor_s(DOF_MATRIX *a, const DOF_REAL_VEC *f, const DOF_SCHAR_VEC *b,
5543 DOF_REAL_VEC *u, REAL omega, REAL tol, int max_iter, int info);
5544
5545 /*** file ssor.c ************************************************************/
5546 int ssor_d(DOF_MATRIX *a, const DOF_REAL_D_VEC *f, const DOF_SCHAR_VEC *b,
5547 DOF_REAL_D_VEC *u, REAL omega, REAL tol, int max_iter, int info);
5548 int ssor_s(DOF_MATRIX *a, const DOF_REAL_VEC *f, const DOF_SCHAR_VEC *b,
5549 DOF_REAL_VEC *u, REAL omega, REAL tol, int max_iter, int info);
5550
5551 /*** file traverse_r.c *****************************************************/
5552 extern void mesh_traverse(MESH *mesh, int level, FLAGS fill_flag,
5553 void (*el_fct)(const EL_INFO *, void *data),
5554 void *data);
5555 extern void fill_macro_info(MESH *mesh, const MACRO_EL *mel, EL_INFO *elinfo);
5556 extern void fill_elinfo(int ichild, FLAGS mask,
5557 const EL_INFO *parent_info, EL_INFO *elinfo);
5558
5559 /*** file traverse_nr.c *****************************************************/
5560 extern TRAVERSE_STACK *get_traverse_stack(void);
5561 extern void free_traverse_stack(TRAVERSE_STACK *stack);
5562 extern const EL_INFO *traverse_first(TRAVERSE_STACK *stack,
5563 MESH *mesh, int level, FLAGS fill_flag);
5564 extern const EL_INFO *traverse_next(TRAVERSE_STACK *stack, const EL_INFO *);
5565 extern const EL_INFO *traverse_neighbour(TRAVERSE_STACK *stack, const EL_INFO *,
5566 int neighbour);
5567 extern const EL_INFO *traverse_parent(const TRAVERSE_STACK *stack,
5568 const EL_INFO *child);
5569 extern const EL_INFO *subtree_traverse_first(TRAVERSE_STACK *stack,
5570 const EL_INFO *local_root,
5571 int level, FLAGS fill_flag);
5572 void clear_traverse_mark(TRAVERSE_STACK *stack);
5573
5574 #define TRAVERSE_FIRST(mesh, level, fill_flag) \
5575 { \
5576 TRAVERSE_STACK *stack = get_traverse_stack(); \
5577 const EL_INFO *el_info; \
5578 if ((el_info = traverse_first(stack, (mesh), (level), (fill_flag)))) do
5579 #define TRAVERSE_NEXT(/**/) \
5580 while ((el_info = traverse_next(stack, el_info))); \
5581 free_traverse_stack(stack); \
5582 }
5583
5584 /* Compatibility */
5585 #define TRAVERSE_START(mesh, level, fill_flag) \
5586 TRAVERSE_FIRST(mesh, level, fill_flag)
5587 #define TRAVERSE_STOP(/**/) TRAVERSE_NEXT()
5588
5589 #define TRAVERSE_NEIGHBOUR(el_info, neighbour) \
5590 traverse_neighbour(stack, el_info, neighbour)
5591
5592 /* file trav_xy.c ***********************************************************/
5593 extern int find_el_at_pt(MESH *mesh, const REAL_D xy,
5594 EL_INFO **el_info_p, FLAGS flag, REAL_B bary,
5595 const MACRO_EL *start_mel,
5596 const REAL_D xy0, REAL *sp);
5597
5598 /*** file read_mesh.c ******************************************************/
5599 extern MESH *
5600 read_mesh(const char *fn, REAL *timeptr,
5601 NODE_PROJ *(*n_proj)(MESH *, MACRO_EL *, int),
5602 MESH *master);
5603 DOF_REAL_VEC *read_dof_real_vec(const char *, MESH *, FE_SPACE *);
5604 DOF_REAL_D_VEC *read_dof_real_d_vec(const char *, MESH *, FE_SPACE *);
5605 DOF_REAL_VEC_D *read_dof_real_vec_d(const char *, MESH *, FE_SPACE *);
5606 DOF_INT_VEC *read_dof_int_vec(const char *, MESH *, FE_SPACE *);
5607 DOF_SCHAR_VEC *read_dof_schar_vec(const char *, MESH *, FE_SPACE *);
5608 DOF_UCHAR_VEC *read_dof_uchar_vec(const char *, MESH *, FE_SPACE *);
5609
5610 extern MESH *
5611 fread_mesh(FILE *fp, REAL *timeptr,
5612 NODE_PROJ *(*n_proj)(MESH *, MACRO_EL *, int),
5613 MESH *master);
5614 DOF_REAL_VEC *fread_dof_real_vec(FILE *fp, MESH *, FE_SPACE *);
5615 DOF_REAL_D_VEC *fread_dof_real_d_vec(FILE *fp, MESH *, FE_SPACE *);
5616 DOF_REAL_VEC_D *fread_dof_real_vec_d(FILE *fp, MESH *, FE_SPACE *);
5617 DOF_INT_VEC *fread_dof_int_vec(FILE *fp, MESH *, FE_SPACE *);
5618 DOF_SCHAR_VEC *fread_dof_schar_vec(FILE *fp, MESH *, FE_SPACE *);
5619 DOF_UCHAR_VEC *fread_dof_uchar_vec(FILE *fp, MESH *, FE_SPACE *);
5620
5621 extern MESH *
5622 read_mesh_xdr(const char *file, REAL *timeptr,
5623 NODE_PROJ *(*)(MESH *, MACRO_EL *, int),
5624 MESH *master);
5625 DOF_REAL_VEC *read_dof_real_vec_xdr(const char *, MESH *, FE_SPACE *);
5626 DOF_REAL_D_VEC *read_dof_real_d_vec_xdr(const char *, MESH *, FE_SPACE *);
5627 DOF_REAL_VEC_D *read_dof_real_vec_d_xdr(const char *, MESH *, FE_SPACE *);
5628 DOF_INT_VEC *read_dof_int_vec_xdr(const char *, MESH *, FE_SPACE *);
5629 DOF_SCHAR_VEC *read_dof_schar_vec_xdr(const char *, MESH *, FE_SPACE *);
5630 DOF_UCHAR_VEC *read_dof_uchar_vec_xdr(const char *, MESH *, FE_SPACE *);
5631
5632 extern MESH *
5633 fread_mesh_xdr(FILE *fp, REAL *timeptr,
5634 NODE_PROJ *(*)(MESH *, MACRO_EL *, int),
5635 MESH *master);
5636 DOF_REAL_VEC *fread_dof_real_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
5637 DOF_REAL_D_VEC *fread_dof_real_d_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
5638 DOF_REAL_VEC_D *fread_dof_real_vec_d_xdr(FILE *fp, MESH *, FE_SPACE *);
5639 DOF_INT_VEC *fread_dof_int_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
5640 DOF_SCHAR_VEC *fread_dof_schar_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
5641 DOF_UCHAR_VEC *fread_dof_uchar_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
5642
5643 #if 0
5644 /* IFF format (multiple objects in one file) */
5645
5646 /* FORM
5647 length
5648 AFEM
5649 MESH
5650 length
5651 data
5652 AVEC
5653 length
5654 data
5655 AVEC
5656 length
5657 data
5658
5659 where AVEC is one of DRV DRDV DUCV DSCV DINV
5660 */
5661
5662 #define IFF_TAG_ALBERTA "AFEM"
5663 #define IFF_TAG_MESH "MESH"
5664 #define IFF_TAG_REAL_VEC "DRV "
5665 #define IFF_TAG_REAL_D_VEC "DRDV"
5666 #define IFF_TAG_INT_VEC "DINV"
5667 #define IFF_TAG_UCHAR_VEC "DUCV"
5668 #define IFF_TAG_SCHAR_VEC "DSCV"
5669
5670 static inline bool is_iff_tag(const char tag1[4], const char tag2[4])
5671 {
5672 return memcmp(tag1, tag2, 4) == 0;
5673 }
5674
5675 /* fread_iff() just read in the tag and the size, afterwards the
5676 * normal read_xdr() routines can be used to suck in the following
5677 * data, depending on what TAG contains.
5678 */
5679 bool fread_iff(FILE *fp, char tag[4], unsigned int *size);
5680 FILE *read_iff(const char *filename, char tag[4], unsigned int *size);
5681 #endif
5682
5683 /*** file write_mesh.c *****************************************************/
5684 bool write_mesh(MESH *, const char *, REAL);
5685 bool write_dof_real_vec(const DOF_REAL_VEC *, const char *);
5686 bool write_dof_real_vec_d(const DOF_REAL_VEC_D *, const char *);
5687 bool write_dof_real_d_vec(const DOF_REAL_D_VEC *, const char *);
5688 bool write_dof_int_vec(const DOF_INT_VEC *, const char *);
5689 bool write_dof_schar_vec(const DOF_SCHAR_VEC *, const char *);
5690 bool write_dof_uchar_vec(const DOF_UCHAR_VEC *, const char *);
5691
5692 bool fwrite_mesh(MESH *, FILE *fp, REAL);
5693 bool fwrite_dof_real_vec(const DOF_REAL_VEC *, FILE *fp);
5694 bool fwrite_dof_real_d_vec(const DOF_REAL_D_VEC *, FILE *fp);
5695 bool fwrite_dof_int_vec(const DOF_INT_VEC *, FILE *fp);
5696 bool fwrite_dof_schar_vec(const DOF_SCHAR_VEC *, FILE *fp);
5697 bool fwrite_dof_uchar_vec(const DOF_UCHAR_VEC *, FILE *fp);
5698
5699 bool write_mesh_xdr(MESH *, const char *, REAL);
5700 bool write_dof_real_vec_xdr(const DOF_REAL_VEC *, const char *);
5701 bool write_dof_real_vec_d_xdr(const DOF_REAL_VEC_D *, const char *);
5702 bool write_dof_real_d_vec_xdr(const DOF_REAL_D_VEC *, const char *);
5703 bool write_dof_int_vec_xdr(const DOF_INT_VEC *, const char *);
5704 bool write_dof_schar_vec_xdr(const DOF_SCHAR_VEC *, const char *);
5705 bool write_dof_uchar_vec_xdr(const DOF_UCHAR_VEC *, const char *);
5706
5707 bool fwrite_mesh_xdr(MESH *, FILE *fp, REAL time);
5708 bool fwrite_dof_real_vec_xdr(const DOF_REAL_VEC *, FILE *fp);
5709 bool fwrite_dof_real_vec_d_xdr(const DOF_REAL_VEC_D *, FILE *fp);
5710 bool fwrite_dof_real_d_vec_xdr(const DOF_REAL_D_VEC *, FILE *fp);
5711 bool fwrite_dof_int_vec_xdr(const DOF_INT_VEC *, FILE *fp);
5712 bool fwrite_dof_schar_vec_xdr(const DOF_SCHAR_VEC *, FILE *fp);
5713 bool fwrite_dof_uchar_vec_xdr(const DOF_UCHAR_VEC *, FILE *fp);
5714
5715 bool write_dof_matrix_pbm(const DOF_MATRIX *matrix, const char *filename);
5716 bool fwrite_dof_matrix_pbm(const DOF_MATRIX *matrix, FILE *file);
5717
5718 /* Writing of IFF-files, writing requires file-positioning support. */
5719 #if 0
5720 bool fwrite_iff(FILE *fp);
5721 bool fterm_iff(FILE *fp);
5722 FILE *fopen_iff(const char *filename, bool append);
5723 bool fclose_iff(FILE *fp);
5724 bool fwrite_mesh_iff(MESH *mesh, REAL time, FILE *fp);
5725 bool fwrite_dof_real_vec_iff(const DOF_REAL_VEC *v, FILE *fp);
5726 bool fwrite_dof_real_d_vec_iff(const DOF_REAL_D_VEC *v, FILE *fp);
5727 bool fwrite_dof_int_vec_iff(const DOF_INT_VEC *v, FILE *fp);
5728 bool fwrite_dof_schar_vec_iff(const DOF_SCHAR_VEC *v, FILE *fp);
5729 bool fwrite_dof_uchar_vec_iff(const DOF_UCHAR_VEC *v, FILE *fp);
5730 #endif
5731
5732 /*** file write_mesh_gmv.c *************************************************/
5733
5734 bool write_mesh_gmv(MESH *mesh, const char *file_name, bool write_ascii,
5735 bool use_refined_grid,
5736 const int n_drv,
5737 DOF_REAL_VEC **drv_ptr,
5738 const int n_drdv,
5739 DOF_REAL_VEC_D **drv_dow_ptr,
5740 DOF_REAL_VEC_D *velocity,
5741 REAL sim_time);
5742
5743 bool write_dof_vec_gmv(MESH *mesh,
5744 const char *mesh_file,
5745 const char *file_name, bool write_ascii,
5746 bool use_refined_grid,
5747 const int n_drv,
5748 DOF_REAL_VEC **drv_ptr,
5749 const int n_drdv,
5750 DOF_REAL_VEC_D **drv_dow_ptr,
5751 DOF_REAL_VEC_D *velocity,
5752 REAL sim_time);
5753
5754 /*-- file write_mesh_ps.c ***************************************************/
5755 void write_mesh_ps(MESH *mesh, const char *filename, const char *title,
5756 const REAL x[2], const REAL y[2], bool keepaspect,
5757 bool draw_bound);
5758
5759 /*******************************************************************************
5760 * interface for Lagrange elements for the gltools
5761 * file gltools.c
5762 ******************************************************************************/
5763 typedef void* GLTOOLS_WINDOW;
5764
5765 GLTOOLS_WINDOW open_gltools_window(const char *, const char *, const REAL *,
5766 MESH *, int);
5767 void close_gltools_window(GLTOOLS_WINDOW);
5768
5769 extern int gltools_get_next_dialog(void);
5770 extern void gltools_set_next_dialog(int dialog);
5771 extern void gltools_est(GLTOOLS_WINDOW, MESH *,
5772 REAL (*rw_est)(EL *el), REAL min, REAL max, REAL time);
5773 extern void gltools_disp_mesh(GLTOOLS_WINDOW, MESH *,
5774 bool mark, const DOF_REAL_VEC_D *, REAL time);
5775 extern void gltools_mesh(GLTOOLS_WINDOW win, MESH *, bool mark, REAL time);
5776 extern void gltools_disp_drv(GLTOOLS_WINDOW, const DOF_REAL_VEC *,
5777 REAL min, REAL max, const DOF_REAL_VEC_D *,
5778 REAL time);
5779 extern void gltools_drv(GLTOOLS_WINDOW, const DOF_REAL_VEC *,
5780 REAL min, REAL max, REAL time);
5781 extern void gltools_disp_drv_d(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
5782 REAL min, REAL max, const DOF_REAL_VEC_D *,
5783 REAL time);
5784 extern void gltools_drv_d(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
5785 REAL min, REAL max, REAL time);
5786 extern void gltools_disp_vec(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
5787 REAL min, REAL max, const DOF_REAL_VEC_D *,
5788 REAL time);
5789 extern void gltools_vec(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
5790 REAL min, REAL max, REAL time);
5791
5792 /*******************************************************************************
5793 * interface for Lagrange elements for the dxtools
5794 * file dxtools.c
5795 ******************************************************************************/
5796
5797 typedef struct dxtools_window DXTOOLS_WINDOW;
5798
5799 extern DXTOOLS_WINDOW *open_dxtools_window(const char *title,
5800 const char *geometry);
5801
5802 extern void close_dxtools_window(DXTOOLS_WINDOW *win);
5803 extern void dxtools_mesh(DXTOOLS_WINDOW *win, MESH *mesh);
5804 extern void dxtools_drv(DXTOOLS_WINDOW *win, const DOF_REAL_VEC *u);
5805 extern void dxtools_drdv(DXTOOLS_WINDOW *win, const DOF_REAL_D_VEC *u);
5806
5807 #ifdef __cplusplus
5808 } /* extern "C" */
5809 #endif
5810
5811 /*******************************************************************************
5812 *
5813 * A couple of header files containing inline functions for various purposes
5814 *
5815 ******************************************************************************/
5816
5817 #include "alberta_inlines.h" /* DIM_OF_WORLD blas */
5818
5819 #include "dof_chains.h" /* support for chains of objects */
5820
5821 #include "el_vec.h" /* element vectors and matrices */
5822
5823 #ifndef __cplusplus
5824 # include "evaluate.h" /* evaluation of finite element functions */
5825 #endif
5826
5827 #endif /* !_ALBERTA_H_ */
5828