1 #ifndef _ALBERTA_INTERN_H_
2 #define _ALBERTA_INTERN_H_
3 /*--------------------------------------------------------------------------*/
4 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
5 /* Bisectioning refinement and Error control by Residual */
6 /* Techniques for scientific Applications */
7 /* */
8 /* file: alberta_intern.h */
9 /* */
10 /* */
11 /* description: private header file of the ALBERTA package */
12 /* */
13 /*--------------------------------------------------------------------------*/
14 /* */
15 /* authors: Alfred Schmidt */
16 /* Zentrum fuer Technomathematik */
17 /* Fachbereich 3 Mathematik/Informatik */
18 /* Universitaet Bremen */
19 /* Bibliothekstr. 2 */
20 /* D-28359 Bremen, Germany */
21 /* */
22 /* Kunibert G. Siebert */
23 /* Institut fuer Mathematik */
24 /* Universitaet Augsburg */
25 /* Universitaetsstr. 14 */
26 /* D-86159 Augsburg, Germany */
27 /* */
28 /* Daniel Koester */
29 /* Institut fuer Mathematik */
30 /* Universitaet Augsburg */
31 /* Universitaetsstr. 14 */
32 /* D-86159 Augsburg, Germany */
33 /* */
34 /* Claus-Justus Heine */
35 /* Abteilung fuer Angewandte Mathematik */
36 /* Albert-Ludwigs-Universitaet Freiburg */
37 /* Hermann-Herder-Str. 10 */
38 /* D-79104 Freiburg im Breisgau, Germany */
39 /* */
40 /* http://www.alberta-fem.de */
41 /* */
42 /* (c) by A. Schmidt and K.G. Siebert (1996-2005) */
43 /* (c) by D. Koester (2002-2005), C.-J. Heine (2002-2009) */
44 /* */
45 /*--------------------------------------------------------------------------*/
46
47 #ifdef HAVE_CONFIG_H
48 # include "config.h" /* probably a good idea to pull this one in here ... */
49 #endif
50
51 #include <sys/types.h>
52 #include <ctype.h>
53 #include <stdio.h>
54 #include <rpc/types.h>
55 #include <rpc/xdr.h>
56 #include <stdlib.h>
57 #include <string.h>
58 #include <time.h>
59
60 #include "alberta.h"
61
62 #ifdef __cplusplus
63 extern "C" {
64 #elif 0
65 } /* some editors try to indent because of the brace above */
66 #endif
67
68 #if !ALBERTA_DEBUG && GCC_ATTRIBUTE_FLATTEN
69 /* "always_inline" causes an internal compiler error, but the
70 * "flatten" seems to do the correct thing. We really want to do
71 * inlining of large functions here because of some constant boolean
72 * arguments which select the respective behaviour. Kind of C++
73 * templates for the poor C programmer.
74 */
75 # define FLATTEN_ATTR __FLATTEN_ATTRIBUTE__
76 # define FORCE_INLINE_ATTR __FORCE_INLINE_ATTRIBUTE__
77 #else
78 # define FLATTEN_ATTR /* nothing */
79 # define FORCE_INLINE_ATTR /* nothing */
80 #endif
81
82 #define LAGRANGE_DEG_MAX 4
83
84 #define FACULTY(N) \
85 (1 * MAX(1, N) * MAX(1, (N)-1) * MAX(1, (N)-2) * MAX(1, (N)-3) * \
86 MAX(1, (N)-4) * MAX(1, (N)-5) * MAX(1, (N)-6) * MAX(1, (N)-7) * \
87 MAX(1, (N)-8) * MAX(1, (N)-9) * MAX(1, (N)-10))
88 #define BINOMIAL(N, K) (FACULTY(N)/(FACULTY(K)*FACULTY((N)-(K))))
89
90 #define N_BAS_LAGRANGE(DEG, DIM) BINOMIAL(((long)((DEG)+(DIM))), ((long)(DEG)))
91
92 /* #define N_BAS_LAG_0D(DEG) N_BAS_LAGRANGE(0, 0) */
93
94 #define N_BAS_LAG_1D(DEG) N_BAS_LAGRANGE(DEG, 1)
95 #define N_BAS_LAG_2D(DEG) N_BAS_LAGRANGE(DEG, 2)
96 #define N_BAS_LAG_3D(DEG) N_BAS_LAGRANGE(DEG, 3)
97
98 typedef struct quad_metadata
99 {
100 /* The per-element cache, this _must_ come first. DO NOT CHANGE. */
101 QUAD_EL_CACHE el_cache;
102 int n_points;
103
104 struct {
105 void *ordinary; /* ordinary quad-fast objects codim 0 & 1 */
106 void *tangential; /* tangential quad-fast objects codim 1 */
107 } quad_fast_head; /* All QUAD_FAST objects for this quadrature */
108
109 /* cached values for parametric meshes, generic part */
110 void *param_data[LAGRANGE_DEG_MAX+1];
111 void (*delete_param_data)(void *param_data);
112 } QUAD_METADATA;
113
114 #define ___AI_STRINGIZE(ARG) #ARG
115 #define __AI_STRINGIZE(ARG) ___AI_STRINGIZE(ARG)
116 #define _AI_STRINGIZE(ARG) __AI_STRINGIZE(ARG)
117
118 #define ___AI_CONCAT(ARG1, ARG2) ARG1##ARG2
119 #define __AI_CONCAT(ARG1, ARG2) ___AI_CONCAT(ARG1, ARG2)
120 #define _AI_CONCAT(ARG1, ARG2) __AI_CONCAT(ARG1, ARG2)
121
122 /*--------------------------------------------------------------------------*/
123 /* information about size of leaf data, which will be hidden at child[1] and*/
124 /* pointers to functions which will be used for transferring such */
125 /* information during refinement/coarsening */
126 /*--------------------------------------------------------------------------*/
127 typedef struct leaf_data_info LEAF_DATA_INFO;
128
129 struct leaf_data_info
130 {
131 size_t leaf_data_size;
132 void (*refine_leaf_data)(EL *parent, EL *child[2]);
133 void (*coarsen_leaf_data)(EL *parent, EL *child[2]);
134 };
135
136 /*--------------------------------------------------------------------------*/
137 /* interpolation of dof vectors during refinement */
138 /* using DOF_VEC_LIST structure */
139 /*--------------------------------------------------------------------------*/
140
141 typedef struct dof_vec_list DOF_VEC_LIST;
142
143 struct dof_vec_list
144 {
145 int size; /* current size of total list */
146 void **list; /* total list, distributed to ones below */
147
148 int n_dof_int_vec, n_dof_dof_vec;
149 int n_dof_uchar_vec, n_dof_schar_vec;
150 int n_dof_real_vec, n_dof_real_d_vec;
151 int n_dof_ptr_vec;
152 int n_dof_matrix;
153
154 DOF_INT_VEC **dof_int_vec;
155 DOF_DOF_VEC **dof_dof_vec;
156 DOF_UCHAR_VEC **dof_uchar_vec;
157 DOF_SCHAR_VEC **dof_schar_vec;
158 DOF_REAL_VEC **dof_real_vec;
159 DOF_REAL_D_VEC **dof_real_d_vec;
160 DOF_PTR_VEC **dof_ptr_vec;
161 DOF_MATRIX **dof_matrix;
162 };
163
164 /*--------------------------------------------------------------------------*/
165 /* the structure stored in the entry mem_info of mesh. */
166 /* It stores information for memory management and macro elements. */
167 /*--------------------------------------------------------------------------*/
168 typedef struct mesh_mem_info MESH_MEM_INFO;
169
170 struct mesh_mem_info
171 {
172 void *dof_ptrs;
173 void *dofs[N_NODE_TYPES];
174 void *element;
175
176 void *rc_list;
177 void *real_d;
178
179 DOF_VEC_LIST *dvlist;
180 DOF_VEC_LIST *dvlist_np; /* non-periodic dvlist on periodic meshes */
181
182 void *leaf_data;
183 LEAF_DATA_INFO leaf_data_info[1];
184
185 /*--------------------------------------------------------------------------*/
186 /* special entries for the use of master/slave grids. */
187 /*--------------------------------------------------------------------------*/
188
189 /* pointers on trace grids */
190 MESH *master; /* pointer to master mesh */
191 /* binding of master mesh */
192 DOF_PTR_VEC *master_binding; /* pointers to master mesh simplices */
193 DOF_PTR_VEC *slave_binding; /* pointers to slave mesh simplices */
194 int next_trace_id;
195
196 /* pointers on master grids */
197 int n_slaves; /* no. of slaves of this mesh */
198 MESH **slaves; /* vector of pointers to slaves */
199
200 /*---8<---------------------------------------------------------------------*/
201 /*--- These last two entries are not really for memory management. ---*/
202 /*--- They describe the array used to store macro element coordinate ---*/
203 /*--- information, which is filled in macro.c and read_mesh.c. ---*/
204 /*--------------------------------------------------------------------->8---*/
205 int count;
206 REAL_D *coords;
207 };
208
_AI_MCOPY_DOW(REAL_DD src,REAL_DD dst)209 static inline REAL_D *_AI_MCOPY_DOW(REAL_DD src, REAL_DD dst)
210 {
211 return MCOPY_DOW((const REAL_D *)src, dst);
212 }
213
214 /* Determine whether the neighbour in an RC_LIST sits on the other
215 * side of a periodic wall. elem and neigh must not be NULL.
216 */
_AI_rc_list_periodic_neigh_p(RC_LIST_EL * rc_el,RC_LIST_EL * rc_neigh)217 static inline int _AI_rc_list_periodic_neigh_p(RC_LIST_EL *rc_el,
218 RC_LIST_EL *rc_neigh)
219 {
220 return
221 rc_el->el_info.el->dof[0] != rc_neigh->el_info.el->dof[0] &&
222 rc_el->el_info.el->dof[1] != rc_neigh->el_info.el->dof[0];
223 }
224
225 /* Same as above, but with a 2-element array of edge DOF pointers.
226 *
227 * BIG FAT NOTE: SIDE-EFFECT:
228 * The edge DOF-pointers are updated; this way the test can be used in
229 * a "for" loop around the refinement edge without having to update
230 * "edge". The orientation of "edge" is maintained.
231 */
_AI_rc_list_periodic_wall_p(RC_LIST_EL * rc_el,DOF * edge[2])232 static inline int _AI_rc_list_periodic_wall_p(RC_LIST_EL *rc_el, DOF *edge[2])
233 {
234 if (rc_el->el_info.el->dof[0] != edge[0] &&
235 rc_el->el_info.el->dof[0] != edge[1]) {
236 if (rc_el->el_info.el->dof[0][0] == edge[0][0]) {
237 edge[0] = rc_el->el_info.el->dof[0];
238 edge[1] = rc_el->el_info.el->dof[1];
239 } else {
240 edge[0] = rc_el->el_info.el->dof[1];
241 edge[1] = rc_el->el_info.el->dof[0];
242 }
243 DEBUG_TEST_EXIT(edge[0][0] < edge[1][0],
244 "Wrong orientation of refinement edge.");
245 return true;
246 } else {
247 return false;
248 }
249 }
250
251 /* How to compare two OPERATOR_INFO structures. the DOWB_ and RDR_
252 * versions need additional tests.
253 */
254 #define OP_INFO_EQ_P(oi1, oi2) \
255 (FE_SPACE_EQ_P((oi1)->row_fe_space, (oi2)->row_fe_space) && \
256 FE_SPACE_EQ_P((oi1)->col_fe_space, (oi2)->col_fe_space) && \
257 (oi1)->quad[2] == (oi2)->quad[2] && \
258 (oi1)->quad[1] == (oi2)->quad[1] && \
259 (oi1)->quad[0] == (oi2)->quad[0] && \
260 (oi1)->quad_tensor[0] == (oi2)->quad_tensor[0] && \
261 (oi1)->quad_tensor[1] == (oi2)->quad_tensor[1] && \
262 (oi1)->quad_tensor[2] == (oi2)->quad_tensor[2] && \
263 (oi1)->init_element == (oi2)->init_element && \
264 \
265 (oi1)->LALt.real == (oi2)->LALt.real && \
266 (oi1)->LALt_type == (oi2)->LALt_type && \
267 (oi1)->LALt_symmetric == (oi2)->LALt_symmetric && \
268 (oi1)->LALt_pw_const == (oi2)->LALt_pw_const && \
269 (oi1)->LALt_degree == (oi2)->LALt_degree && \
270 \
271 (oi1)->Lb0.real == (oi2)->Lb0.real && \
272 (oi1)->Lb0_pw_const == (oi2)->Lb0_pw_const && \
273 \
274 (oi1)->Lb1.real == (oi2)->Lb1.real && \
275 (oi1)->Lb1_pw_const == (oi2)->Lb1_pw_const && \
276 \
277 (oi1)->Lb0_Lb1_anti_symmetric == (oi2)->Lb0_Lb1_anti_symmetric && \
278 (oi1)->Lb_type == (oi2)->Lb_type && \
279 (oi1)->Lb_degree == (oi2)->Lb_degree && \
280 (oi1)->advection_field == (oi2)->advection_field && \
281 (oi1)->adv_fe_space == (oi2)->adv_fe_space && \
282 \
283 (oi1)->c.real == (oi2)->c.real && \
284 (oi1)->c_type == (oi2)->c_type && \
285 (oi1)->c_pw_const == (oi2)->c_pw_const && \
286 (oi1)->c_degree == (oi2)->c_degree && \
287 \
288 (oi1)->fill_flag == (oi2)->fill_flag /*&& \
289 (oi1)->user_data == (oi2)->user_data*/)
290
291 /* block-matrix/block-operator types for DIM_OF_WORLD/scalar combinations */
292
293 typedef enum op_block_type {
294 OP_TYPE_SS = 0,
295 OP_TYPE_SV,
296 OP_TYPE_VS,
297 OP_TYPE_CV,
298 OP_TYPE_VC,
299 OP_TYPE_VV,
300 N_OP_BLOCK_TYPES
301 } OP_BLOCK_TYPE;
302
operator_type(const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space)303 static inline OP_BLOCK_TYPE operator_type(const FE_SPACE *row_fe_space,
304 const FE_SPACE *col_fe_space)
305 {
306 FUNCNAME("operator_type");
307 const BAS_FCTS *row_bfcts, *col_bfcts;
308 OP_BLOCK_TYPE op_type;
309
310 row_bfcts = row_fe_space->bas_fcts;
311 col_bfcts = col_fe_space->bas_fcts;
312
313 if (row_fe_space->rdim == 1) {
314 op_type = col_bfcts->rdim == 1 ? OP_TYPE_SS : OP_TYPE_SV;
315 } else if (row_bfcts->rdim == 1) {
316 op_type = col_bfcts->rdim == 1 ? OP_TYPE_SS : OP_TYPE_CV;
317 } else if (col_fe_space->rdim == 1) {
318 op_type = OP_TYPE_VS;
319 } else {
320 op_type = col_bfcts->rdim == 1 ? OP_TYPE_VC : OP_TYPE_VV;
321 }
322 TEST_EXIT(VECTOR_BASIS_FUNCTIONS || op_type == OP_TYPE_SS,
323 "DIM_OF_WORLD-valued basis functions not supported in this "
324 "version of ALBERTA. re-configure and re-compile the library "
325 "without \"--disable-vector-basis-functions\".\n");
326
327 return op_type;
328 }
329
matent_type(const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space,MATENT_TYPE op_mat_type)330 static inline MATENT_TYPE matent_type(const FE_SPACE *row_fe_space,
331 const FE_SPACE *col_fe_space,
332 MATENT_TYPE op_mat_type)
333 {
334 switch (operator_type(row_fe_space, col_fe_space)) {
335 case OP_TYPE_VV:
336 case OP_TYPE_VS:
337 case OP_TYPE_SV:
338 return MATENT_REAL;
339 case OP_TYPE_SS:
340 return op_mat_type;
341 case OP_TYPE_VC:
342 case OP_TYPE_CV:
343 return MATENT_REAL_D;
344 default:
345 return MATENT_NONE;
346 }
347 }
348
349 /* Below are some template like constructs to emit variants of small
350 * functions for constant dimensions. The assemble stuff makes
351 * extensive use of this.
352 */
353
354 /* Make sure n_lambda is a constant by emitting several copies of the
355 * same function (the compiler cannot know that n_lambda only takes
356 * the values 1, 2, 3 or 4, otherwise it could decide by itself to
357 * emit several copies of the same function for different values of
358 * n_lambda).
359 */
360 #define DIM_NAME(NAMEBASE, DIM) NAMEBASE##_##DIM##D
361
362 #define EMIT_DIM_VERSION(CLASS, NAMEBASE, DIM, ARG_DECL, ARG_CALL) \
363 CLASS void FLATTEN_ATTR DIM_NAME(NAMEBASE, DIM) ARG_DECL \
364 { \
365 NAMEBASE ARG_CALL(N_LAMBDA(DIM)); \
366 } \
367 struct _AI_semicolon_dummy
368
369 #if DIM_MAX == 1
370 # define EMIT_DIM_VERSIONS(CLASS, NAMEBASE, ARG_DECL, ARG_CALL) \
371 EMIT_DIM_VERSION(CLASS, NAMEBASE, 1, ARG_DECL, ARG_CALL)
372 # define DECLARE_DIM_VERSIONS(CLASS, NAMEBASE, ARG_DECL) \
373 CLASS void DIM_NAME(NAMEBASE, 1) ARG_DECL
374 # define DEF_DIM_DUMMIES(NAMEBASE, FTYPE) \
375 static const FTYPE DIM_NAME(NAMEBASE, 1) = NULL
376 #elif DIM_MAX == 2
377 # define EMIT_DIM_VERSIONS(CLASS, NAMEBASE, ARG_DECL, ARG_CALL) \
378 EMIT_DIM_VERSION(CLASS, NAMEBASE, 1, ARG_DECL, ARG_CALL); \
379 EMIT_DIM_VERSION(CLASS, NAMEBASE, 2, ARG_DECL, ARG_CALL)
380 # define DECLARE_DIM_VERSIONS(CLASS, NAMEBASE, ARG_DECL) \
381 CLASS void DIM_NAME(NAMEBASE, 1) ARG_DECL; \
382 CLASS void DIM_NAME(NAMEBASE, 2) ARG_DECL
383 # define DEF_DIM_DUMMIES(NAMEBASE, FTYPE) \
384 static const FTYPE DIM_NAME(NAMEBASE, 1) = NULL; \
385 static const FTYPE DIM_NAME(NAMEBASE, 2) = NULL
386 #elif DIM_MAX == 3
387 # define EMIT_DIM_VERSIONS(CLASS, NAMEBASE, ARG_DECL, ARG_CALL) \
388 EMIT_DIM_VERSION(CLASS, NAMEBASE, 1, ARG_DECL, ARG_CALL); \
389 EMIT_DIM_VERSION(CLASS, NAMEBASE, 2, ARG_DECL, ARG_CALL); \
390 EMIT_DIM_VERSION(CLASS, NAMEBASE, 3, ARG_DECL, ARG_CALL)
391 # define DECLARE_DIM_VERSIONS(CLASS, NAMEBASE, ARG_DECL) \
392 CLASS void DIM_NAME(NAMEBASE, 1) ARG_DECL; \
393 CLASS void DIM_NAME(NAMEBASE, 2) ARG_DECL; \
394 CLASS void DIM_NAME(NAMEBASE, 3) ARG_DECL
395 # define DEF_DIM_DUMMIES(NAMEBASE, FTYPE) \
396 static const FTYPE DIM_NAME(NAMEBASE, 1) = NULL; \
397 static const FTYPE DIM_NAME(NAMEBASE, 2) = NULL; \
398 static const FTYPE DIM_NAME(NAMEBASE, 3) = NULL
399 #else
400 # error unsupported DIM_MAX
401 #endif
402
403 /* Generate the proper function name for a given dimension. */
404 #if DIM_MAX == 1
405 # define DIM_VARIANT(NAMEBASE, DIM) DIM_NAME(NAMEBASE, 1)
406 #elif DIM_MAX == 2
407 # define DIM_VARIANT(NAMEBASE, DIM) \
408 ((DIM) == 1 ? DIM_NAME(NAMEBASE, 1) : DIM_NAME(NAMEBASE, 2))
409 #elif DIM_MAX == 3
410 # define DIM_VARIANT(NAMEBASE, DIM) \
411 ((DIM) == 1 \
412 ? DIM_NAME(NAMEBASE, 1) \
413 : ((DIM) == 2 \
414 ? DIM_NAME(NAMEBASE, 2) \
415 : DIM_NAME(NAMEBASE, 3)))
416 #else
417 # error unsupported DIM_MAX
418 #endif
419
420 /*******************************************************************************
421 * stack data structure for non-recursive mesh traversal
422 ******************************************************************************/
423
424 /* should not this be a private data-structure??? */
425 struct traverse_stack
426 {
427 MESH *traverse_mesh;
428 int traverse_level;
429 FLAGS traverse_flags;
430 FLAGS fill_flag;
431
432 const MACRO_EL *traverse_mel;
433 int stack_size;
434 int stack_used;
435 EL_INFO *elinfo_stack;
436 U_CHAR *info_stack;
437
438 const MACRO_EL *save_traverse_mel;
439 EL_INFO *save_elinfo_stack;
440 U_CHAR *save_info_stack;
441 int save_stack_used;
442
443 int el_count;
444
445 struct {
446 int pos;
447 int traverse_level;
448 FLAGS traverse_flags;
449 /* fill-flags are in stack->elinfo_stack[marker.pos] */
450 } marker;
451
452 TRAVERSE_STACK *next;
453 };
454
455 void __AI_enlarge_traverse_stack(TRAVERSE_STACK *stack);
456
457 /*** several lookup-tables for 3d *******************************************/
458 #if DIM_MAX > 2
459
460 /* slave_numbering_3d[type=0,1][orientation][face][vertex]: Define the local
461 * slave vertex index given in all cases for master orientation and face.
462 * The rules are: slave numbering is counterclockwise when looking
463 * at the slave triangle from outside the master tetrahedron. The numbering
464 * is defined so that a master refinement edge always corresponds to a slave
465 * refinement edge. Unfortunately, the numbering must depend on the parent
466 * orientation to remain consistent.
467 */
468 static const int slave_numbering_3d[2][2][N_WALLS_3D][N_VERTICES_3D] =
469 {/* type 0 */ {/* orientation + */ {/* face 0 */ {-1, 1, 2, 0},
470 /* face 1 */ { 1,-1, 0, 2},
471 /* face 2 */ { 0, 1,-1, 2},
472 /* face 3 */ { 1, 0, 2,-1}},
473 /* orientation - */ {/* face 0 */ {-1, 0, 2, 1},
474 /* face 1 */ { 0,-1, 1, 2},
475 /* face 2 */ { 1, 0,-1, 2},
476 /* face 3 */ { 0, 1, 2,-1}}},
477 /* type 1 */ {/* orientation + */ {/* face 0 */ {-1, 0, 1, 2},
478 /* face 1 */ { 1,-1, 0, 2},
479 /* face 2 */ { 0, 1,-1, 2},
480 /* face 3 */ { 1, 0, 2,-1}},
481 /* orientation - */ {/* face 0 */ {-1, 1, 0, 2},
482 /* face 1 */ { 0,-1, 1, 2},
483 /* face 2 */ { 1, 0,-1, 2},
484 /* face 3 */ { 0, 1, 2,-1}}}};
485
486 /****************************************************************************
487 * master_edge_3d[orientation][master face][slave edge]: Gives the master
488 * edge corresponding to the slave edge.
489 ****************************************************************************/
490
491 static const int master_edge_3d[2][2][N_WALLS_3D][N_EDGES_2D] =
492 {/* type 0 */ {/* orientation + */ {/* face 0 */ { 3, 5, 4},
493 /* face 1 */ { 2, 5, 1},
494 /* face 2 */ { 4, 2, 0},
495 /* face 3 */ { 1, 3, 0}},
496 /* orientation - */ {/* face 0 */ { 5, 3, 4},
497 /* face 1 */ { 5, 2, 1},
498 /* face 2 */ { 2, 4, 0},
499 /* face 3 */ { 3, 1, 0}}},
500 /* type 1 */ {/* orientation + */ {/* face 0 */ { 4, 3, 5},
501 /* face 1 */ { 2, 5, 1},
502 /* face 2 */ { 4, 2, 0},
503 /* face 3 */ { 1, 3, 0}},
504 /* orientation - */ {/* face 0 */ { 4, 5, 3},
505 /* face 1 */ { 5, 2, 1},
506 /* face 2 */ { 2, 4, 0},
507 /* face 3 */ { 3, 1, 0}}}};
508
509 /* child_face_3d[type][parent_face][child_no] give the local face
510 * number on the child, when parent_face is the local face number of
511 * on the parent.
512 */
513 static const int child_face_3d[3][4][2] = {
514 /* type 0: */ { /* face 0: */ {-1, 3},
515 /* face 1: */ { 3,-1},
516 /* face 2: */ { 1, 2},
517 /* face 3: */ { 2, 1}},
518 /* type 1: */ { /* face 0: */ {-1, 3},
519 /* face 1: */ { 3,-1},
520 /* face 2: */ { 1, 1},
521 /* face 3: */ { 2, 2}},
522 /* type 2: */ { /* face 0: */ {-1, 3},
523 /* face 1: */ { 3,-1},
524 /* face 2: */ { 1, 1},
525 /* face 3: */ { 2, 2}}
526 };
527
528 /* The inverse lookup-array: given the child_face, the parent_type and
529 * the child number, "return" the parent's face number of -1.
530 *
531 * parent_face_3d[parent_type][child_face][ichild]
532 */
533 static const int parent_face_3d[3][N_FACES_3D][2] = {
534 /* type 0: */ { /* face 0: */ {-1,-1 },
535 /* face 1: */ { 2, 3 },
536 /* face 2: */ { 3, 2 },
537 /* face 3: */ { 1, 0 }},
538 /* type 1: */ { /* face 0: */ {-1,-1 },
539 /* face 1: */ { 2, 2 },
540 /* face 2: */ { 3, 3 },
541 /* face 3: */ { 1, 0 }},
542 /* type 2: */ { /* face 0: */ {-1,-1 },
543 /* face 1: */ { 2, 2 },
544 /* face 2: */ { 3, 3 },
545 /* face 3: */ { 1, 0 }},
546 };
547
548 /*--------------------------------------------------------------------------*/
549 /* child_vertex_3d[el_type][child][i] = */
550 /* father's local vertex index of new vertex i */
551 /* 4 stands for the newly generated vertex */
552 /*--------------------------------------------------------------------------*/
553 static const int child_vertex_3d[3][2][4] = {
554 {{0, 2, 3, 4}, {1, 3, 2, 4}},
555 {{0, 2, 3, 4}, {1, 2, 3, 4}},
556 {{0, 2, 3, 4}, {1, 2, 3, 4}}
557 };
558
559 /*--------------------------------------------------------------------------*/
560 /* child_edge_3d[el_type][child][i] = */
561 /* father's local edge index of new edge i */
562 /* new edge 2 is half of old edge 0, */
563 /* new edges 4, 5 are really new edges, and value is different: */
564 /* child_edge_3d[][][4, 5] = index of same edge in other child */
565 /*--------------------------------------------------------------------------*/
566 static const int child_edge_3d[3][2][6] = {
567 {{1, 2, 0, 5, 5, 4}, {4, 3, 0, 5, 5, 4}},
568 {{1, 2, 0, 5, 4, 5}, {3, 4, 0, 5, 4, 5}},
569 {{1, 2, 0, 5, 4, 5}, {3, 4, 0, 5, 4, 5}}
570 };
571
572 /*--------------------------------------------------------------------------*/
573 /* child_orientation[el_type][child] = */
574 /* +1 if orientation is not changed during refinement */
575 /* -1 if orientation is changed during refinement */
576 /*--------------------------------------------------------------------------*/
577 static const S_CHAR child_orientation_3d[3][2] = {{1, 1}, {1, -1}, {1, -1}};
578
579 /*--------------------------------------------------------------------------*/
580 /* n_child_edge_3d[el_type][ichild][dir] */
581 /* gives local index of new edge on child[ichild] part of face [2+dir] on */
582 /* the parent */
583 /*--------------------------------------------------------------------------*/
584 static const int n_child_edge_3d[3][2][2] = {
585 {{5, 4}, {4, 5}},
586 {{5, 4}, {5, 4}},
587 {{5, 4}, {5, 4}}
588 };
589
590 /*--------------------------------------------------------------------------*/
591 /* n_child_face_3d[el_type][ichild][dir] */
592 /* gives local index of sub-face on child[ichild] part of face [2+dir] on */
593 /* the parent */
594 /*--------------------------------------------------------------------------*/
595 static const int n_child_face_3d[3][2][2] = {
596 {{1, 2}, {2, 1}}, {{1, 2}, {1, 2}}, {{1, 2}, {1, 2}}
597 };
598
599 /*--------------------------------------------------------------------------*/
600 /* adjacent_child_3d[position][ichild] */
601 /* gives number of the adjacent child on a neighbour element */
602 /* position = 0 same position of element and neigh at refinement edge */
603 /* position = 1 different ... */
604 /*--------------------------------------------------------------------------*/
605 static const int adjacent_child_3d[2][2] = {{0, 1}, {1, 0}};
606
607 #endif
608
609 /*** file memory.c ******************************************************/
610 extern void *AI_get_leaf_data(MESH *mesh);
611 extern DOF_ADMIN *AI_get_dof_admin(MESH *mesh, const char *name,
612 const int n_dof[N_NODE_TYPES]);
613 extern void AI_free_leaf_data(void *leaf_data, MESH *mesh);
614 extern void AI_reactivate_dof(MESH *mesh, const EL *el,
615 DOF **edge_twins, DOF **face_twins);
616 extern DOF *AI_get_dof_memory(MESH *mesh, int position);
617 extern void AI_free_dof_memory(DOF *dof, MESH *mesh, int position);
618 extern void AI_get_dof_list(MESH *mesh, int i);
619 extern void AI_get_dof_ptr_list(MESH *mesh);
620 extern void AI_fill_missing_dofs(MESH *mesh);
621 extern void AI_advance_cookies_rec(MESH *mesh);
622 extern MESH *
623 _AI_get_mesh(int dim, const char *name,
624 const MACRO_DATA *macro_data,
625 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
626 AFF_TRAFO *(*init_wall_trafos)(MESH *, MACRO_EL *, int wall),
627 bool strict_periodic);
628 extern DOF *_AI_get_dof(MESH *mesh, int position, bool alloc_index);
629
630 /*** file macro.c *******************************************************/
631 extern S_CHAR AI_get_orientation_3d(MACRO_EL *mel);
632 extern void _AI_compute_element_wall_transformations(MACRO_DATA *data);
633 extern void
634 _AI_fill_bound_info(MESH *mesh, int *mel_vertices, int nv, int ne, bool count);
635 extern void
636 _AI_macro_data2mesh(MESH *mesh, const MACRO_DATA *data,
637 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
638 AFF_TRAFO *(*init_wall_trafos)(
639 MESH *, MACRO_EL *, int wall),
640 bool strict_periodic);
641
642 /*** file dof_admin.c ***************************************************/
643 extern bool _AI_check_matrix_types(MATENT_TYPE mat_type, MATENT_TYPE elm_type);
644 extern void _AI_allocate_n_dofs(DOF_ADMIN *admin, int used_count);
645
646 /*** file write_mesh.c ***************************************************/
647 extern
648 size_t AI_fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream);
649
650 /*** file read_mesh.c ***************************************************/
651 extern
652 size_t AI_fread(const void *ptr, size_t size, size_t nmemb, FILE *stream);
653 extern bool_t AI_xdr_REAL(XDR *xdr, void *rp);
654 extern bool_t AI_xdr_U_CHAR(XDR *xdr, void *ucp);
655 extern bool_t AI_xdr_S_CHAR(XDR *xdr, void *cp);
656 extern bool_t AI_xdr_DOF(XDR *xdr, void *dp);
657 extern XDR *AI_xdr_open_file(const char *fn, enum xdr_op mode);
658 extern bool AI_xdr_close_file(XDR *xdr);
659 extern XDR *AI_xdr_fopen(FILE *file, enum xdr_op mode);
660 extern bool AI_xdr_close(XDR *xdr);
661
662 extern void _AI_read_REAL(REAL *val);
663 extern void _AI_read_int(int *val);
664 extern void _AI_read_string(char *string, int strileng);
665 extern void _AI_read_var_string(char **string);
666 extern void _AI_read_vector(void *start, int n, size_t size, xdrproc_t xdrproc);
667 extern void _AI_read_U_CHAR(U_CHAR *val);
668 extern void _AI_read_S_CHAR(S_CHAR *val);
669
670 /*** file read_mesh_xdr_1.2.c *******************************************/
671 extern MESH *_AI_read_mesh_1_2(REAL *timeptr);
672 extern void _AI_match_node_types(int *node_vec);
673
674 /*** file refine.c ******************************************************/
675 extern void AI_refine_fct_1d(const EL_INFO *el_info, void *data);
676 extern void AI_coarse_fct_1d(const EL_INFO *el_info, void *data);
677 extern void _AI_refine_update_bbox(MESH *mesh, const REAL_D new_coord);
678
679 #if DIM_MAX > 1
680 extern void AI_bisect_element_2d(MESH *mesh, EL *el, DOF *dof[3]);
681 extern void AI_bisect_patch_2d(MESH *mesh, RC_LIST_EL ref_list[],
682 int n_neighs);
683 extern void AI_coarse_patch_2d(MESH *mesh,RC_LIST_EL coarse_list[],
684 int n_neigh);
685 #endif
686 extern DOF_VEC_LIST *AI_get_dof_vec_list(MESH *mesh);
687 extern DOF_VEC_LIST *AI_get_dof_vec_list_np(MESH *mesh);
688 extern void AI_free_dof_vec_list(MESH *mesh);
689 extern void AI_free_dof_vec_list_np(MESH *mesh);
690 extern void AI_set_neighs_on_patch_3d(RC_LIST_EL ref_list[], int n_neigh,
691 int bound);
692 extern int AI_split_rc_list_3d(RC_LIST_EL *orig,
693 RC_LIST_EL *copy, int n_left);
694 extern void AI_reverse_rc_list_3d(RC_LIST_EL ref_list[],
695 int n_neigh, DOF *edge[2]);
696 extern RC_LIST_EL *AI_rotate_rc_list_3d(RC_LIST_EL *first,
697 int n_neigh, DOF *edge[2]);
698
699 /*** file submesh.c *****************************************************/
700 extern void AI_check_slavery(MESH *master);
701
702 /*** file traverse_nr.c ***************************************************/
703 extern void AI_test_traverse_nr(MESH *mesh, int level, FLAGS fill_flag);
704 extern void AI_update_elinfo_stack_3d(TRAVERSE_STACK *stack);
705
706 /*** file traverse_r.c ***************************************************/
707 extern void AI_test_traverse(MESH *mesh, int level, FLAGS fill_flag);
708 extern void AI_update_elinfo_3d(EL_INFO *elinfo);
709
710 /*** file periodic.c ********************************************************/
711 extern int
712 _AI_wall_trafo_vertex_orbit(int dim,
713 int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2],
714 int nwt,
715 int v, int *orbit, int nv);
716 extern int
717 _AI_wall_trafo_vertex_orbits(int dim,
718 int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2],
719 int nwt,
720 int *orbit_map, int *nvp);
721 #if DIM_MAX > 2
722 extern int
723 _AI_wall_trafo_edge_orbit(int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2],
724 int nwt,
725 int e, int *orbit,
726 int (*edges)[2], int ne);
727 extern int
728 _AI_wall_trafo_edge_orbits(int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2],
729 int nwt,
730 int *orbit_map,
731 int (*edges)[2], int ne);
732 #endif
733 extern int
734 _AI_compute_macro_wall_trafos(MESH *mesh,
735 int (**wall_trafos_p)[N_VERTICES(DIM_MAX-1)][2]);
736
737 /*** file gauss-quad.c *******************************************************/
738 void _AI_gauss_quad(int kind, int n, REAL alpha,
739 REAL beta, int kpts, REAL endpts[2],
740 REAL t[/*n*/], REAL w[/*n*/]);
741
742 /*** file eval.c *************************************************************/
743
744
745 extern REAL _AI_inter_fct_loc_param(const EL_INFO *el_info,
746 const QUAD *quad, int iq,
747 void *ud);
748 extern REAL _AI_inter_fct_loc(const EL_INFO *el_info,
749 const QUAD *quad, int iq,
750 void *ud);
751 extern const REAL *_AI_inter_fct_loc_d(REAL_D result,
752 const EL_INFO *el_info,
753 const QUAD *quad, int iq,
754 void *ud);
755 extern const REAL *_AI_inter_fct_loc_d_param(REAL_D result,
756 const EL_INFO *el_info,
757 const QUAD *quad, int iq,
758 void *ud);
759
760 /*** file parametric.c *******************************************************/
761
762 extern bool _AI_is_lagrange_parametric(MESH *mesh);
763 extern PARAM_STRATEGY _AI_lagrange_strategy(MESH *mesh);
764
765 /*** file block_precon.c ****************************************************/
766
767 extern const PRECON *_AI_vget_block_diag_precon(const DOF_MATRIX *A,
768 const DOF_SCHAR_VEC *mask,
769 int info,
770 va_list ap);
771 extern const PRECON *_AI_get_block_diag_precon(const DOF_MATRIX *A,
772 const DOF_SCHAR_VEC *mask,
773 int info,
774 ...);
775 extern const PRECON *_AI_vget_block_SSOR_precon(const DOF_MATRIX *A,
776 const DOF_SCHAR_VEC *mask,
777 int info,
778 va_list ap);
779 extern const PRECON *_AI_get_block_precon(const DOF_MATRIX *A,
780 const DOF_SCHAR_VEC *mask,
781 int info,
782 const PRECON_TYPE *prec_type);
783
784 #ifdef __cplusplus
785 }
786 #endif
787
788 #endif /* !_ALBERTA_INTERN_H_ */
789