1 /*******************************************************************************
2 *
3 * ALBERTA: an Adaptive multi Level finite element toolbox using
4 * Bisectioning refinement and Error control by Residual
5 * Techniques for scientific Applications
6 *
7 * www.alberta-fem.de
8 *
9 * file: assemble_wall.c
10 *
11 * description: assembling of contributions across interior walls, meant as
12 * support for DG methods.
13 *
14 *
15 *******************************************************************************
16 *
17 * This file's authors: Claus-Justus Heine
18 * Abteilung f�r Angewandte Mathematik
19 * Universit�t Freiburg
20 * Hermann-Herder-Stra�e 10
21 * 79104 Freiburg, Germany
22 *
23 * (c) by C.-J. Heine (2009-2009)
24 *
25 ******************************************************************************/
26
27 #include "alberta_intern.h"
28 #include "alberta.h"
29 #include "assemble_bndry.h"
30
31 #if !HAVE_ROW_FCTS_V_TYPE && !HAVE_COL_FCTS_V_TYPE
32 # define VEC_PFX SS
33 # define EMIT_SS_VERSIONS 1
34 #elif HAVE_ROW_FCTS_V_TYPE && HAVE_COL_FCTS_V_TYPE
35 # define VEC_PFX VV
36 # define EMIT_VV_VERSIONS 1
37 #elif HAVE_ROW_FCTS_V_TYPE && HAVE_COL_FCTS_C_TYPE
38 # define VEC_PFX VC
39 # define EMIT_VC_VERSIONS 1
40 #elif HAVE_ROW_FCTS_C_TYPE && HAVE_COL_FCTS_V_TYPE
41 # define VEC_PFX CV
42 # define EMIT_CV_VERSIONS 1
43 #elif !HAVE_M_DST_TYPE && HAVE_ROW_FCTS_V_TYPE && !HAVE_COL_FCTS_V_TYPE
44 # define VEC_PFX VS
45 # define EMIT_VS_VERSIONS 1
46 #elif !HAVE_M_DST_TYPE && !HAVE_ROW_FCTS_V_TYPE && HAVE_COL_FCTS_V_TYPE
47 # define VEC_PFX SV
48 # define EMIT_SV_VERSIONS 1
49 #else
50 # error impossible block-matrix combinations
51 #endif
52
53 #define HAVE_VSSV (EMIT_SV_VERSIONS || EMIT_VS_VERSIONS)
54
55 #define VNAME(base) _AI_CONCAT(VEC_PFX, _AI_CONCAT(_, base))
56
57 /* <<< el_wall_fcts lookup table */
58
59 #define N_TYPE_DEFUNS 6 /* MM, MDM, MSCM, DMDM, DMSCM, SCMSCM */
60
61 #if !HAVE_VSSV
62 extern AI_EL_WALL_FCT VNAME(MM_el_wall_fcts)[DIM_MAX+1][N_WALLS_MAX][5][16];
63 extern AI_EL_WALL_FCT VNAME(MDM_el_wall_fcts)[DIM_MAX+1][N_WALLS_MAX][5][16];
64 extern AI_EL_WALL_FCT VNAME(MSCM_el_wall_fcts)[DIM_MAX+1][N_WALLS_MAX][5][16];
65 #endif
66 extern AI_EL_WALL_FCT VNAME(DMDM_el_wall_fcts)[DIM_MAX+1][N_WALLS_MAX][5][16];
67 extern AI_EL_WALL_FCT VNAME(DMSCM_el_wall_fcts)[DIM_MAX+1][N_WALLS_MAX][5][16];
68 extern AI_EL_WALL_FCT VNAME(SCMSCM_el_wall_fcts)[DIM_MAX+1][N_WALLS_MAX][5][16];
69
70 /* This function is called with srctype <= dsttype, it computes the
71 * index into the el_wall_fcts[] table.
72 */
type_index(MATENT_TYPE dsttype,MATENT_TYPE srctype)73 static inline unsigned int type_index(MATENT_TYPE dsttype,
74 MATENT_TYPE srctype)
75 {
76 unsigned int dstnr, srcnr;
77
78 switch (dsttype) {
79 case MATENT_REAL_DD: dstnr = 0; break;
80 case MATENT_REAL_D: dstnr = 1; break;
81 case MATENT_REAL: dstnr = 2; break;
82 default: dstnr = ~0; break;
83 }
84 switch (srctype) {
85 case MATENT_REAL_DD: srcnr = 0; break;
86 case MATENT_REAL_D: srcnr = 1; break;
87 case MATENT_REAL: srcnr = 2; break;
88 default: srcnr = ~0; break;
89 }
90
91 #define N_TYPES 3
92 /* just the difference between two arithmetic sums */
93 return ((2*N_TYPES+1)*dstnr - SQR(dstnr))/2 + srcnr - dstnr;
94 #undef N_TYPES
95 }
96
97 /* >>> */
98
99 /* <<< element_matrix() definitions */
100
101 /* Poor C-programmers explicit template specialization. Further below
102 * specific versions of the element matrices are emitted where those
103 * flags are constant s.t. the compiler can optimize the code.
104 *
105 * Note that "MIXED" -- meaning differing ansatz- and test-space -- is
106 * in some sense always true because we pair the test space on the
107 * current element with the ansatz space on the neighbour element.
108 */
109
110 #define ZERO_ORDER (1 << 0)
111 #define FRST_ORDER (1 << 1)
112 #define SCND_ORDER (1 << 2)
113 #define INIT_EL (1 << 3)
114 #define OI_INIT_EL (1 << 4)
115 #define MIXED (1 << 5)
116 #define TANGENTIAL (1 << 6)
117 #define PARTPARAM (1 << 7) /* only used in get_fill_data() */
118
119 /* NOTE: these function just adds something to an existing
120 * element-matrix, there is no BNDRY_EL_MATRIX_INFO().
121 */
122
123
124 /* <<< init_element helper */
125
init_objects(BNDRY_FILL_INFO * info,int wall,U_CHAR mask)126 static inline void init_objects(BNDRY_FILL_INFO *info, int wall, U_CHAR mask)
127 {
128 if (mask & INIT_EL) {
129 if (mask & SCND_ORDER) {
130 INIT_OBJECT(info->row_wquad_fast[2]);
131 }
132 if (mask & FRST_ORDER) {
133 INIT_OBJECT(info->row_wquad_fast[1]);
134 }
135 if (mask & ZERO_ORDER) {
136 INIT_OBJECT(info->row_wquad_fast[0]);
137 }
138 if (mask & MIXED) {
139 if (mask & SCND_ORDER) {
140 INIT_OBJECT(info->col_wquad_fast[2]);
141 }
142 if (mask & FRST_ORDER) {
143 INIT_OBJECT(info->col_wquad_fast[1]);
144 }
145 if (mask & ZERO_ORDER) {
146 INIT_OBJECT(info->col_wquad_fast[0]);
147 }
148 }
149
150 /* Reallocate scl_el_mat if necessary */
151 ROW_CHAIN_DO(info, BNDRY_FILL_INFO) {
152 COL_CHAIN_DO(info, BNDRY_FILL_INFO) {
153 if (mask & TANGENTIAL) {
154 info->row_fcts_trace_map[wall] =
155 info->op_info.row_fe_space->bas_fcts->trace_dof_map[0][0][wall];
156 info->n_trace_row_fcts[wall] =
157 info->op_info.row_fe_space->bas_fcts->n_trace_bas_fcts[wall];
158 }
159
160 if (info->scl_el_mat != NULL) {
161 int n_row_max =
162 info->op_info.row_fe_space->bas_fcts->n_bas_fcts_max;
163 int n_col_max =
164 info->op_info.col_fe_space->bas_fcts->n_bas_fcts_max;
165 if (n_row_max > info->n_row_max || n_col_max > info->n_col_max) {
166 #undef MAT_BODY
167 #define MAT_BODY(F, CC, C, S, TYPE) \
168 MAT_FREE(info->el_mat, info->n_row_max, info->n_col_max, TYPE); \
169 info->scl_el_mat = \
170 (void *)MAT_ALLOC(n_row_max, n_col_max, TYPE)
171 MAT_EMIT_BODY_SWITCH(info->krn_blk_type);
172 info->n_row_max = n_row_max;
173 info->n_col_max = n_col_max;
174 }
175 }
176 } COL_CHAIN_WHILE(info, BNDRY_FILL_INFO);
177 } ROW_CHAIN_WHILE(info, BNDRY_FILL_INFO);
178 }
179
180 info->cur_el = NULL;
181 info->cur_el_info = NULL;
182 }
183
184 static inline INIT_EL_TAG
init_element(const EL_INFO * el_info,BNDRY_FILL_INFO * info,U_CHAR mask)185 init_element(const EL_INFO *el_info, BNDRY_FILL_INFO *info, U_CHAR mask)
186 {
187 INIT_EL_TAG init_tag = INIT_EL_TAG_NONE;
188
189 if (info->cur_el != el_info->el || info->cur_el_info != el_info) {
190 #if HAVE_ROW_FCTS_V_TYPE
191 int iw;
192 #endif
193 int dim = el_info->mesh->dim;
194
195 ROW_CHAIN_DO(info, BNDRY_FILL_INFO) {
196 COL_CHAIN_DO(info, BNDRY_FILL_INFO) {
197 INIT_ELEMENT_SINGLE(el_info, info->op_info.row_fe_space->bas_fcts);
198 info->el_mat->n_row = info->op_info.row_fe_space->bas_fcts->n_bas_fcts;
199
200 if (mask & TANGENTIAL) {
201 int wall;
202 for (wall = 0; wall < N_WALLS(dim); ++wall) {
203 info->row_fcts_trace_map[wall] =
204 info->op_info.row_fe_space->bas_fcts->trace_dof_map[0][0][wall];
205 info->n_trace_row_fcts[wall] =
206 info->op_info.row_fe_space->bas_fcts->n_trace_bas_fcts[wall];
207 }
208 }
209 } COL_CHAIN_WHILE(info, BNDRY_FILL_INFO);
210 } ROW_CHAIN_WHILE(info, BNDRY_FILL_INFO);
211
212 if (mask & SCND_ORDER) {
213 init_tag |= INIT_ELEMENT(el_info, info->row_wquad_fast[2]);
214 #if HAVE_ROW_FCTS_V_TYPE
215 /* force the computation of the gradients */
216 if (init_tag != INIT_EL_TAG_NULL) {
217 for (iw = 0; iw < N_WALLS(dim); iw++) {
218 const QUAD_FAST *cache = info->row_wquad_fast[2]->quad_fast[iw];
219 CHAIN_DO(cache, const QUAD_FAST) {
220 if (!cache->bas_fcts->dir_pw_const) {
221 get_quad_fast_grd_phi_dow(cache);
222 }
223 } CHAIN_WHILE(cache, const QUAD_FAST);
224 }
225 }
226 #endif
227 }
228 if (mask & FRST_ORDER) {
229 init_tag |= INIT_ELEMENT(el_info, info->row_wquad_fast[1]);
230 #if HAVE_ROW_FCTS_V_TYPE
231 /* force the computation of the gradients and/or values */
232 if (init_tag != INIT_EL_TAG_NULL) {
233 for (iw = 0; iw < N_WALLS(dim); iw++) {
234 const QUAD_FAST *cache = info->row_wquad_fast[1]->quad_fast[iw];
235 CHAIN_DO(cache, const QUAD_FAST) {
236 if (!cache->bas_fcts->dir_pw_const) {
237 if (cache->init_flag & INIT_GRD_PHI) {
238 get_quad_fast_grd_phi_dow(cache);
239 }
240 if (cache->init_flag & INIT_PHI) {
241 get_quad_fast_phi_dow(cache);
242 }
243 }
244 } CHAIN_WHILE(cache, const QUAD_FAST);
245 }
246 }
247 #endif
248 }
249 if (mask & ZERO_ORDER) {
250 init_tag |= INIT_ELEMENT(el_info, info->row_wquad_fast[0]);
251 #if HAVE_ROW_FCTS_V_TYPE
252 /* force the computation of the gradients and/or values */
253 if (init_tag != INIT_EL_TAG_NULL) {
254 for (iw = 0; iw < N_WALLS(dim); iw++) {
255 const QUAD_FAST *cache = info->row_wquad_fast[0]->quad_fast[iw];
256 CHAIN_DO(cache, const QUAD_FAST) {
257 get_quad_fast_phi_dow(cache);
258 } CHAIN_WHILE(cache, const QUAD_FAST);
259 }
260 }
261 #endif
262 }
263 info->cur_el = el_info->el;
264 info->cur_el_info = el_info;
265 if (init_tag == INIT_EL_TAG_NULL) {
266 /* any none-NULL tag has set another bit in init_tag */
267 return init_tag;
268 }
269 }
270 return init_tag;
271 }
272
273 /* >>> */
274
275 /* <<< element_matrix_default */
276
277 static inline
VNAME(element_matrix_default)278 const EL_MATRIX *VNAME(element_matrix_default)(const EL_INFO *el_info,
279 int wall,
280 const void *fill_info,
281 U_CHAR mask)
282 {
283 BNDRY_FILL_INFO *info = (BNDRY_FILL_INFO *)fill_info;
284 EL_INFO neigh_info[1];
285
286 if (el_info == NULL) {
287 init_objects(info, wall, mask);
288 return NULL;
289 }
290
291 if (el_info->neigh[wall] == NULL) {
292 return NULL;
293 }
294
295 if (mask & INIT_EL) {
296 const EL_GEOM_CACHE *elgc;
297 if (init_element(el_info, info, mask) == INIT_EL_TAG_NULL) {
298 return NULL;
299 }
300 elgc = fill_el_geom_cache(el_info, FILL_EL_WALL_REL_ORIENTATION(wall));
301 fill_neigh_el_info(neigh_info, el_info, wall, elgc->rel_orientation[wall]);
302 INIT_ELEMENT(neigh_info, info->op_info.col_fe_space->bas_fcts);
303 }
304
305 ROW_CHAIN_DO(info, BNDRY_FILL_INFO) {
306 COL_CHAIN_DO(info, BNDRY_FILL_INFO) {
307 void **velmat = (void **)info->el_mat->data.real;
308 int i, j;
309
310 if ((mask & OI_INIT_EL)) {
311 info->op_info.init_element(
312 el_info, wall, info->op_info.quad, info->op_info.user_data);
313 }
314
315 if (mask & INIT_EL) {
316 info->el_mat->n_col = info->op_info.col_fe_space->bas_fcts->n_bas_fcts;
317 }
318
319 #undef MAT_BODY
320 #define MAT_BODY(F, CC, C, S, TYPE) \
321 for(i = 0; i < info->el_mat->n_row; i++) \
322 for (j = 0; j < info->el_mat->n_col; j++) \
323 F##SET_DOW(0.0, info->el_mat->data.S[i][j]);
324 MAT_EMIT_BODY_SWITCH(info->el_mat->type);
325
326 if (mask & SCND_ORDER) {
327 info->col_quad_fast[2] =
328 get_neigh_quad_fast(el_info, info->col_wquad_fast[2], wall);
329 if (mask & INIT_EL) {
330 if (info->col_quad_fast[2]) {
331 INIT_ELEMENT_SINGLE(neigh_info, info->col_quad_fast[2]);
332 info->dflt_second_order[wall](el_info, info, velmat);
333 }
334 } else {
335 info->dflt_second_order[wall](el_info, info, velmat);
336 }
337 }
338 if (mask & FRST_ORDER) {
339 info->col_quad_fast[1] =
340 get_neigh_quad_fast(el_info, info->col_wquad_fast[1], wall);
341 if (mask & INIT_EL) {
342 if (info->col_quad_fast[1]) {
343 INIT_ELEMENT_SINGLE(neigh_info, info->col_quad_fast[1]);
344 info->dflt_first_order[wall](el_info, info, velmat);
345 }
346 } else {
347 info->dflt_first_order[wall](el_info, info, velmat);
348 }
349 }
350 if (mask & ZERO_ORDER) {
351 info->col_quad_fast[0] =
352 get_neigh_quad_fast(el_info, info->col_wquad_fast[0], wall);
353 if (mask & INIT_EL) {
354 if (info->col_quad_fast[0]) {
355 INIT_ELEMENT_SINGLE(neigh_info, info->col_quad_fast[0]);
356 info->dflt_zero_order[wall](el_info, info, velmat);
357 }
358 } else {
359 info->dflt_zero_order[wall](el_info, info, velmat);
360 }
361 }
362 } COL_CHAIN_WHILE(info, BNDRY_FILL_INFO);
363 } ROW_CHAIN_WHILE(info, BNDRY_FILL_INFO);
364
365 return info->el_mat;
366 }
367
368 /* >>> */
369
370 /* <<< element_matrix_partparam */
371
372 static inline
VNAME(element_matrix_partparam)373 const EL_MATRIX *VNAME(element_matrix_partparam)(const EL_INFO *el_info,
374 int wall,
375 const void *fill_info,
376 U_CHAR mask)
377 {
378 BNDRY_FILL_INFO *info = (BNDRY_FILL_INFO *)fill_info;
379 bool is_parametric = true;
380 EL_INFO neigh_info[1];
381
382 if (el_info == NULL) {
383 init_objects(info, wall, mask);
384 return NULL;
385 }
386 if (el_info->neigh[wall] == NULL) {
387 return NULL;
388 }
389
390 if (info->cur_el != el_info->el || info->cur_el_info != el_info) {
391 is_parametric = info->parametric->init_element(el_info, info->parametric);
392 }
393
394 if (mask & INIT_EL) {
395 const EL_GEOM_CACHE *elgc;
396 if (init_element(el_info, info, mask) == INIT_EL_TAG_NULL) {
397 return NULL;
398 }
399 elgc = fill_el_geom_cache(el_info, FILL_EL_WALL_REL_ORIENTATION(wall));
400 fill_neigh_el_info(neigh_info, el_info, wall, elgc->rel_orientation[wall]);
401 INIT_ELEMENT(neigh_info, info->op_info.col_fe_space->bas_fcts);
402 }
403
404 ROW_CHAIN_DO(info, BNDRY_FILL_INFO) {
405 COL_CHAIN_DO(info, BNDRY_FILL_INFO) {
406 void **velmat = (void **)info->el_mat->data.real;
407 int i, j;
408
409 if (mask & OI_INIT_EL) {
410 info->op_info.init_element(
411 el_info, wall, info->op_info.quad, info->op_info.user_data);
412 }
413
414 if (mask & INIT_EL) {
415 info->el_mat->n_col = info->op_info.col_fe_space->bas_fcts->n_bas_fcts;
416 }
417
418 #undef MAT_BODY
419 #define MAT_BODY(F, CC, C, S, TYPE) \
420 for(i = 0; i < info->el_mat->n_row; i++) \
421 for (j = 0; j < info->el_mat->n_col; j++) \
422 F##SET_DOW(0.0, info->el_mat->data.S[i][j]);
423 MAT_EMIT_BODY_SWITCH(info->el_mat->type);
424
425 if (mask & SCND_ORDER) {
426 info->col_quad_fast[2] =
427 get_neigh_quad_fast(el_info, info->col_wquad_fast[2], wall);
428 if (mask & INIT_EL) {
429 if (info->col_quad_fast[2]) {
430 INIT_ELEMENT_SINGLE(neigh_info, info->col_quad_fast[2]);
431 if (!is_parametric && info->fast_second_order[wall]) {
432 info->fast_second_order[wall](el_info, info, velmat);
433 } else {
434 info->dflt_second_order[wall](el_info, info, velmat);
435 }
436 }
437 } else {
438 if (!is_parametric && info->fast_second_order[wall]) {
439 info->fast_second_order[wall](el_info, info, velmat);
440 } else {
441 info->dflt_second_order[wall](el_info, info, velmat);
442 }
443 }
444 }
445 if (mask & FRST_ORDER) {
446 info->col_quad_fast[1] =
447 get_neigh_quad_fast(el_info, info->col_wquad_fast[1], wall);
448 if (mask & INIT_EL) {
449 if (info->col_quad_fast[1]) {
450 INIT_ELEMENT_SINGLE(neigh_info, info->col_quad_fast[1]);
451 if (!is_parametric && info->fast_first_order[wall]) {
452 info->fast_first_order[wall](el_info, info, velmat);
453 } else {
454 info->dflt_first_order[wall](el_info, info, velmat);
455 }
456 }
457 } else {
458 if (!is_parametric && info->fast_first_order[wall]) {
459 info->fast_first_order[wall](el_info, info, velmat);
460 } else {
461 info->dflt_first_order[wall](el_info, info, velmat);
462 }
463 }
464 }
465 if (mask & ZERO_ORDER) {
466 info->col_quad_fast[0] =
467 get_neigh_quad_fast(el_info, info->col_wquad_fast[0], wall);
468 if (mask & INIT_EL) {
469 if (info->col_quad_fast[0]) {
470 INIT_ELEMENT_SINGLE(neigh_info, info->col_quad_fast[0]);
471 if (!is_parametric && info->fast_first_order[wall]) {
472 info->fast_zero_order[wall](el_info, info, velmat);
473 } else {
474 info->dflt_zero_order[wall](el_info, info, velmat);
475 }
476 }
477 } else {
478 if (!is_parametric && info->fast_first_order[wall]) {
479 info->fast_zero_order[wall](el_info, info, velmat);
480 } else {
481 info->dflt_zero_order[wall](el_info, info, velmat);
482 }
483 }
484 }
485 } COL_CHAIN_WHILE(info, BNDRY_FILL_INFO);
486 } ROW_CHAIN_WHILE(info, BNDRY_FILL_INFO);
487
488 return info->el_mat;
489 }
490
491 /* >>> */
492
493 /* <<< element matrix lookup-table */
494
495 /* We instantiate various versions of the element matrices where the
496 * flags bits defined above are kept constant s.t. the compiler can
497 * optimize the code for the various cases:
498 *
499 * #define ZERO_ORDER (1 << 0)
500 * #define FRST_ORDER (1 << 1)
501 * #define SCND_ORDER (1 << 2)
502 * #define INIT_EL (1 << 3)
503 * #define OI_INIT_EL (1 << 4)
504 * #define MIXED (1 << 5)
505 * #define TANGENTIAL (1 << 6)
506 * #define PARTPARAM (1 << 7)
507 *
508 */
509
510 #define ELEMENT_MATRIX_DEFAULT_FUN(wall, n2, n1) \
511 VNAME(element_matrix_default_##wall##_##n2##_##n1)
512 #define DEFUN_ELEMENT_MATRIX_DEFAULT(wall, n2, n1) \
513 static const EL_MATRIX * FLATTEN_ATTR \
514 ELEMENT_MATRIX_DEFAULT_FUN(wall, n2, n1)(const EL_INFO *el_info, \
515 void *fill_info) \
516 { \
517 return VNAME(element_matrix_default)(el_info, wall, fill_info, n2|n1); \
518 } \
519 struct _AI_semicolon_dummy
520 #if DIM_MAX == 1
521 # define DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, n1) \
522 DEFUN_ELEMENT_MATRIX_DEFAULT(0, n2, n1); \
523 DEFUN_ELEMENT_MATRIX_DEFAULT(1, n2, n1)
524 #elif DIM_MAX == 2
525 # define DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, n1) \
526 DEFUN_ELEMENT_MATRIX_DEFAULT(0, n2, n1); \
527 DEFUN_ELEMENT_MATRIX_DEFAULT(1, n2, n1); \
528 DEFUN_ELEMENT_MATRIX_DEFAULT(2, n2, n1)
529 #elif DIM_MAX == 3
530 # define DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, n1) \
531 DEFUN_ELEMENT_MATRIX_DEFAULT(0, n2, n1); \
532 DEFUN_ELEMENT_MATRIX_DEFAULT(1, n2, n1); \
533 DEFUN_ELEMENT_MATRIX_DEFAULT(2, n2, n1); \
534 DEFUN_ELEMENT_MATRIX_DEFAULT(3, n2, n1)
535 #else
536 # error unsupported DIM_MAX
537 #endif
538
539 #define ELEMENT_MATRIX_PARTPARAM_FUN(wall, n2, n1) \
540 VNAME(element_matrix_partparm_##wall##_##n2##_##n1)
541 #define DEFUN_ELEMENT_MATRIX_PARTPARAM(wall, n2, n1) \
542 static const EL_MATRIX * FLATTEN_ATTR \
543 ELEMENT_MATRIX_PARTPARAM_FUN(wall, n2, n1)(const EL_INFO *el_info, \
544 void *fill_info) \
545 { \
546 return VNAME(element_matrix_partparam)(el_info, wall, fill_info, n2|n1); \
547 } \
548 struct _AI_semicolon_dummy
549 #if DIM_MAX == 1
550 # define DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, n1) \
551 DEFUN_ELEMENT_MATRIX_PARTPARAM(0, n2, n1); \
552 DEFUN_ELEMENT_MATRIX_PARTPARAM(1, n2, n1)
553 #elif DIM_MAX == 2
554 # define DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, n1) \
555 DEFUN_ELEMENT_MATRIX_PARTPARAM(0, n2, n1); \
556 DEFUN_ELEMENT_MATRIX_PARTPARAM(1, n2, n1); \
557 DEFUN_ELEMENT_MATRIX_PARTPARAM(2, n2, n1)
558 #elif DIM_MAX >= 3
559 # define DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, n1) \
560 DEFUN_ELEMENT_MATRIX_PARTPARAM(0, n2, n1); \
561 DEFUN_ELEMENT_MATRIX_PARTPARAM(1, n2, n1); \
562 DEFUN_ELEMENT_MATRIX_PARTPARAM(2, n2, n1); \
563 DEFUN_ELEMENT_MATRIX_PARTPARAM(3, n2, n1)
564 #else
565 # error unsupported DIM_MAX
566 #endif
567
568 /* Stuff for non-parametric meshes, or entirely parametric meshes, at
569 * least one of the three least significant bits must be set, coding
570 * for the presence of a zero, first, second order part in the
571 * BNDRY_OPERATOR_INFO structure.
572 */
573 #define DEFUNS_16_EL_MAT_DEFAULT(n2) \
574 /* 0x0 */ \
575 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x1); \
576 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x2); \
577 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x3); \
578 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x4); \
579 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x5); \
580 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x6); \
581 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x7); \
582 /* 0x8 */ \
583 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0x9); \
584 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0xA); \
585 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0xB); \
586 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0xC); \
587 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0xD); \
588 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0xE); \
589 DEFUNS_ELEMENT_MATRIX_DEFAULT(n2, 0xF)
590
591 /* Instantiate the templates for each of the three least significant
592 * bits of the upper nibble, partparam is implemented using a
593 * different function template.
594 */
595 DEFUNS_16_EL_MAT_DEFAULT(0x00);
596 DEFUNS_16_EL_MAT_DEFAULT(0x10);
597 DEFUNS_16_EL_MAT_DEFAULT(0x20);
598 DEFUNS_16_EL_MAT_DEFAULT(0x30);
599 DEFUNS_16_EL_MAT_DEFAULT(0x40);
600 DEFUNS_16_EL_MAT_DEFAULT(0x50);
601 DEFUNS_16_EL_MAT_DEFAULT(0x60);
602 DEFUNS_16_EL_MAT_DEFAULT(0x70);
603
604 #if DIM_MAX == 1
605 # define NULL_N_WALLS { NULL, NULL }
606 # define ELEMENT_MATRIX_DEFAULT_FUNS(n2, n1) \
607 { ELEMENT_MATRIX_DEFAULT_FUN(0, n2, n1), \
608 ELEMENT_MATRIX_DEFAULT_FUN(1, n2, n1) }
609 # define ELEMENT_MATRIX_PARTPARAM_FUNS(n2, n1) \
610 { ELEMENT_MATRIX_PARTPARAM_FUN(0, n2, n1), \
611 ELEMENT_MATRIX_PARTPARAM_FUN(1, n2, n1) }
612 #elif DIM_MAX == 2
613 # define NULL_N_WALLS { NULL, NULL, NULL }
614 # define ELEMENT_MATRIX_DEFAULT_FUNS(n2, n1) \
615 { ELEMENT_MATRIX_DEFAULT_FUN(0, n2, n1), \
616 ELEMENT_MATRIX_DEFAULT_FUN(1, n2, n1), \
617 ELEMENT_MATRIX_DEFAULT_FUN(2, n2, n1) }
618 # define ELEMENT_MATRIX_PARTPARAM_FUNS(n2, n1) \
619 { ELEMENT_MATRIX_PARTPARAM_FUN(0, n2, n1), \
620 ELEMENT_MATRIX_PARTPARAM_FUN(1, n2, n1), \
621 ELEMENT_MATRIX_PARTPARAM_FUN(2, n2, n1) }
622 #elif DIM_MAX == 3
623 # define NULL_N_WALLS { NULL, NULL, NULL, NULL }
624 # define ELEMENT_MATRIX_DEFAULT_FUNS(n2, n1) \
625 { ELEMENT_MATRIX_DEFAULT_FUN(0, n2, n1), \
626 ELEMENT_MATRIX_DEFAULT_FUN(1, n2, n1), \
627 ELEMENT_MATRIX_DEFAULT_FUN(2, n2, n1), \
628 ELEMENT_MATRIX_DEFAULT_FUN(3, n2, n1) }
629 # define ELEMENT_MATRIX_PARTPARAM_FUNS(n2, n1) \
630 { ELEMENT_MATRIX_PARTPARAM_FUN(0, n2, n1), \
631 ELEMENT_MATRIX_PARTPARAM_FUN(1, n2, n1), \
632 ELEMENT_MATRIX_PARTPARAM_FUN(2, n2, n1), \
633 ELEMENT_MATRIX_PARTPARAM_FUN(3, n2, n1) }
634 #else
635 # error unsupported DIM_MAX
636 #endif
637
638 /* Initializer for the look-up table. "n2" must loop form 0 -- 7. */
639 #define EL_MAT_16_DEFAULT_FUNS(n2) \
640 NULL_N_WALLS, \
641 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x1), \
642 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x2), \
643 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x3), \
644 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x4), \
645 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x5), \
646 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x6), \
647 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x7), \
648 NULL_N_WALLS, \
649 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0x9), \
650 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0xA), \
651 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0xB), \
652 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0xC), \
653 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0xD), \
654 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0xE), \
655 ELEMENT_MATRIX_DEFAULT_FUNS(n2, 0xF)
656
657 /* Stuff for partly parametric meshes, at least one of the three least
658 * significant bits must be set, coding for the presence of a zero,
659 * first, second order part in the BNDRY_OPERATOR_INFO structure. The
660 * OI_INIT_EL bit (bit 4) must be set.
661 */
662 #define DEFUNS_16_EL_MAT_PARTPARAM(n2) \
663 /* 0x0 */ \
664 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x1); \
665 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x2); \
666 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x3); \
667 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x4); \
668 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x5); \
669 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x6); \
670 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x7); \
671 /* 0x8 */ \
672 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0x9); \
673 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0xA); \
674 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0xB); \
675 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0xC); \
676 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0xD); \
677 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0xE); \
678 DEFUNS_ELEMENT_MATRIX_PARTPARAM(n2, 0xF)
679
680 DEFUNS_16_EL_MAT_PARTPARAM(0x00);
681 DEFUNS_16_EL_MAT_PARTPARAM(0x10);
682 DEFUNS_16_EL_MAT_PARTPARAM(0x20);
683 DEFUNS_16_EL_MAT_PARTPARAM(0x30);
684 DEFUNS_16_EL_MAT_PARTPARAM(0x40);
685 DEFUNS_16_EL_MAT_PARTPARAM(0x50);
686 DEFUNS_16_EL_MAT_PARTPARAM(0x60);
687 DEFUNS_16_EL_MAT_PARTPARAM(0x70);
688
689 /* initializer for the lookup-table */
690 #define EL_MAT_16_PARTPARAM_FUNS(n2) \
691 NULL_N_WALLS, \
692 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x1), \
693 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x2), \
694 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x3), \
695 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x4), \
696 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x5), \
697 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x6), \
698 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x7), \
699 NULL_N_WALLS, \
700 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0x9), \
701 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0xA), \
702 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0xB), \
703 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0xC), \
704 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0xD), \
705 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0xE), \
706 ELEMENT_MATRIX_PARTPARAM_FUNS(n2, 0xF)
707
708 EL_MATRIX_FCT VNAME(neigh_el_mat_fct_table)[256][N_WALLS_MAX] =
709 {
710 /* Non-parametric or entirely parametric meshes. */
711 EL_MAT_16_DEFAULT_FUNS(0x00),
712 EL_MAT_16_DEFAULT_FUNS(0x10),
713 EL_MAT_16_DEFAULT_FUNS(0x20),
714 EL_MAT_16_DEFAULT_FUNS(0x30),
715 EL_MAT_16_DEFAULT_FUNS(0x40),
716 EL_MAT_16_DEFAULT_FUNS(0x50),
717 EL_MAT_16_DEFAULT_FUNS(0x60),
718 EL_MAT_16_DEFAULT_FUNS(0x70),
719 /* Partly parametric meshes. */
720 EL_MAT_16_PARTPARAM_FUNS(0x00),
721 EL_MAT_16_PARTPARAM_FUNS(0x10),
722 EL_MAT_16_PARTPARAM_FUNS(0x20),
723 EL_MAT_16_PARTPARAM_FUNS(0x30),
724 EL_MAT_16_PARTPARAM_FUNS(0x40),
725 EL_MAT_16_PARTPARAM_FUNS(0x50),
726 EL_MAT_16_PARTPARAM_FUNS(0x60),
727 EL_MAT_16_PARTPARAM_FUNS(0x70),
728 };
729
730 /* >>> */
731
732 /* >>> */
733
734 #if EMIT_SS_VERSIONS
735
736 /* lookup-tables for element function */
737 static EL_MATRIX_FCT
738 (*neigh_el_mat_fct_table[N_OP_BLOCK_TYPES])[N_WALLS_MAX] = {
739 /* [OP_TYPE_SS] = */ SS_neigh_el_mat_fct_table,
740 #if VECTOR_BASIS_FUNCTIONS
741 /* [OP_TYPE_SV] = */ SV_neigh_el_mat_fct_table,
742 /* [OP_TYPE_VS] = */ VS_neigh_el_mat_fct_table,
743 /* [OP_TYPE_CV] = */ CV_neigh_el_mat_fct_table,
744 /* [OP_TYPE_VC] = */ VC_neigh_el_mat_fct_table,
745 /* [OP_TYPE_VV] = */ VV_neigh_el_mat_fct_table
746 #endif
747 };
748
unify_bop_info(BNDRY_OPERATOR_INFO * oinfo,const BNDRY_OPERATOR_INFO * oi_orig,const WALL_QUAD_TENSOR * quad_tensor[3],const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space,MATENT_TYPE blk_type)749 static bool unify_bop_info(BNDRY_OPERATOR_INFO *oinfo,
750 const BNDRY_OPERATOR_INFO *oi_orig,
751 const WALL_QUAD_TENSOR *quad_tensor[3],
752 const FE_SPACE *row_fe_space,
753 const FE_SPACE *col_fe_space,
754 MATENT_TYPE blk_type)
755 {
756 const BAS_FCTS *row_fcts, *col_fcts;
757 PARAMETRIC *parametric = NULL;
758 int dim;
759 int row_fcts_deg, col_fcts_deg, max_deg, i;
760
761 *oinfo = *oi_orig;
762 for (i = 0; i < 3; i++) {
763 oinfo->quad_tensor[i] = quad_tensor[i];
764 if (quad_tensor[i]) {
765 oinfo->quad[i] = quad_tensor[i]->quad;
766 }
767 }
768
769 oinfo->row_fe_space = row_fe_space;
770 oinfo->col_fe_space = col_fe_space;
771
772 /* <<< consistency checks etc. */
773
774 row_fcts = oinfo->row_fe_space->bas_fcts;
775 col_fcts = oinfo->col_fe_space->bas_fcts;
776
777 if(col_fcts->dim != row_fcts->dim) {
778 ERROR("Support dimensions of col_fcts and row_fcts do not match!\n");
779 ERROR("cannot initialize EL_MATRIX_INFO; returning NULL\n");
780 return false;
781 }
782
783 dim = col_fcts->dim;
784
785 row_fcts_deg = row_fcts->unchained->degree;
786 col_fcts_deg = col_fcts->unchained->degree;
787
788 parametric = oinfo->row_fe_space->mesh->parametric;
789
790 if (!oinfo->c.real_dd &&
791 !oinfo->Lb0.real_dd && !oinfo->Lb1.real_dd &&
792 !oinfo->LALt.real_dd) {
793 ERROR("no function for 2nd, 1st, and 0 order term;\n");
794 ERROR("can not initialize EL_MATRIX_INFO; returning NULL\n");
795 return false;
796 }
797
798 if (oinfo->LALt.real_dd == NULL) {
799 /* reset all LALt related settings to zero */
800 oinfo->LALt_type = MATENT_REAL;
801 oinfo->LALt_pw_const = false;
802 oinfo->LALt_symmetric = false;
803 oinfo->LALt_degree = 0;
804 oinfo->quad[2] = NULL;
805 oinfo->quad_tensor[2] = NULL;
806 }
807 if (oinfo->Lb0.real_dd == NULL) {
808 oinfo->Lb0_pw_const = false;
809 }
810 if (oinfo->Lb1.real_dd == NULL) {
811 oinfo->Lb1_pw_const = false;
812 }
813 if (oinfo->Lb0.real_dd == NULL && oinfo->Lb1.real_dd == NULL) {
814 oinfo->Lb_type = MATENT_REAL;
815 oinfo->Lb0_Lb1_anti_symmetric = false;
816 oinfo->Lb_degree = 0;
817 oinfo->advection_field = NULL; /* unused, though */
818 oinfo->adv_fe_space = NULL; /* unused, though */
819 oinfo->quad[1] = NULL;
820 oinfo->quad_tensor[1] = NULL;
821 }
822 if (oinfo->c.real_dd == NULL) {
823 oinfo->c_pw_const = false;
824 oinfo->c_type = MATENT_REAL;
825 oinfo->c_degree = 0;
826 oinfo->quad[0] = NULL;
827 oinfo->quad_tensor[0] = NULL;
828 }
829
830 if(parametric) {
831 if(!oinfo->quad[0] && !oinfo->quad[1] && !oinfo->quad[2]) {
832 ERROR("User is responsible for providing at least one quadrature\n");
833 ERROR("when using a parametric mesh!\n");
834 ERROR("can not initialize EL_MATRIX_INFO; returning NULL\n");
835 return false;
836 }
837 }
838
839 /* the element matrices are always non-symmetric because we pair the
840 * test-space of the current element with the ansatz space on the
841 * neighbour element.
842 */
843 oinfo->LALt_symmetric =
844 oinfo->Lb0_Lb1_anti_symmetric = false;
845
846 max_deg = 0;
847
848 if (oinfo->c.real_dd && oinfo->quad[0] == NULL) {
849 if (oinfo->c_pw_const) {
850 oinfo->c_degree = 0;
851 }
852 max_deg = MAX(max_deg, row_fcts_deg + col_fcts_deg + oinfo->c_degree);
853 }
854
855 if ((oinfo->Lb0.real_dd || oinfo->Lb1.real_dd) && oinfo->quad[1] == NULL) {
856 if (oinfo->Lb0_pw_const && oinfo->Lb1_pw_const) {
857 oinfo->Lb_degree = 0;
858 }
859 max_deg = MAX(max_deg, row_fcts_deg + col_fcts_deg - 1 + oinfo->Lb_degree);
860 }
861
862 if (oinfo->LALt.real_dd && oinfo->quad[2] == NULL) {
863 if (oinfo->LALt_pw_const) {
864 oinfo->LALt_degree = 0;
865 }
866 max_deg =
867 MAX(max_deg, row_fcts_deg + col_fcts_deg - 2 + oinfo->LALt_degree);
868 }
869
870 #if 1
871 if (oinfo->LALt.real_dd && !oinfo->quad[2]) {
872 oinfo->quad[2] = get_wall_quad(dim, max_deg);
873 } else if (!oinfo->LALt.real_dd) {
874 oinfo->LALt_degree = 0;
875 oinfo->quad[2] = NULL;
876 }
877
878 if ((oinfo->Lb0.real_dd || oinfo->Lb1.real_dd) && !oinfo->quad[1]) {
879 if ((!oinfo->Lb0_pw_const || !oinfo->Lb1_pw_const) && oinfo->quad[2]) {
880 oinfo->quad[1] = oinfo->quad[2];
881 } else {
882 oinfo->quad[1] = get_wall_quad(dim, row_fcts_deg + col_fcts_deg - 1);
883 }
884 } else if (!oinfo->Lb0.real_dd && !oinfo->Lb1.real_dd) {
885 oinfo->Lb_degree = 0;
886 oinfo->quad[1] = NULL;
887 }
888
889 if (oinfo->c.real_dd && !oinfo->quad[0]) {
890 if (!oinfo->c_pw_const && oinfo->quad[2]) {
891 oinfo->quad[0] = oinfo->quad[2];
892 } else if (!oinfo->c_pw_const && oinfo->quad[1]) {
893 oinfo->quad[0] = oinfo->quad[1];
894 } else {
895 oinfo->quad[0] = get_wall_quad(dim, row_fcts_deg + col_fcts_deg);
896 }
897 } else if (!oinfo->c.real_dd) {
898 oinfo->c_degree = 0;
899 oinfo->quad[0] = NULL;
900 }
901 #else
902 if (oinfo->LALt.real_dd && !oinfo->quad[2]) {
903 oinfo->quad[2] = get_wall_quad(dim, row_fcts_deg + col_fcts_deg - 2);
904 } else if (!oinfo->LALt.real_dd) {
905 oinfo->quad[2] = NULL;
906 }
907
908 if ((oinfo->Lb0.real_dd || oinfo->Lb1.real_dd) && !oinfo->quad[1]) {
909 if ((!oinfo->Lb0_pw_const || !oinfo->Lb1_pw_const) && oinfo->quad[2]) {
910 oinfo->quad[1] = oinfo->quad[2];
911 } else {
912 oinfo->quad[1] = get_wall_quad(dim, row_fcts_deg + col_fcts_deg - 1);
913 }
914 } else if (!oinfo->Lb0.real_dd && !oinfo->Lb1.real_dd) {
915 oinfo->quad[1] = NULL;
916 }
917
918 if (oinfo->c.real_dd && !oinfo->quad[0]) {
919 if (!oinfo->c_pw_const && oinfo->quad[2]) {
920 oinfo->quad[0] = oinfo->quad[2];
921 } else if (!oinfo->c_pw_const && oinfo->quad[1]) {
922 oinfo->quad[0] = oinfo->quad[1];
923 } else {
924 oinfo->quad[0] = get_wall_quad(dim, row_fcts_deg + col_fcts_deg);
925 }
926 } else if (!oinfo->c.real_dd) {
927 oinfo->quad[0] = NULL;
928 }
929 #endif
930
931 /* >>> */
932
933 return true;
934 }
935
936 static BNDRY_FILL_INFO *first_fill_info = NULL;
937
__get_neigh_fill_info(const BNDRY_OPERATOR_INFO * oinfo,MATENT_TYPE krn_blk_type)938 static BNDRY_FILL_INFO *__get_neigh_fill_info(const BNDRY_OPERATOR_INFO *oinfo,
939 MATENT_TYPE krn_blk_type)
940 {
941 FUNCNAME("__get_neigh_fill_info");
942 BNDRY_FILL_INFO *fill_info;
943 const BAS_FCTS *row_fcts, *col_fcts;
944 PARAMETRIC *parametric = NULL;
945 int dim, not_all_param = false;
946 unsigned int dflt_flags[3];
947 unsigned int fast_flags[3];
948 U_CHAR init_row_fcts[3];
949 U_CHAR init_col_fcts[3];
950 int i, wall, wall_fcts_first;
951 unsigned int el_mat_idx, typeidx = 0;
952 OP_BLOCK_TYPE op_type;
953 const WALL_QUAD_FAST **row_fcts_fast, **col_fcts_fast;
954
955 op_type = operator_type(oinfo->row_fe_space, oinfo->col_fe_space);
956
957 /* <<< generate a new fill_info object */
958
959 fill_info = MEM_CALLOC(1, BNDRY_FILL_INFO);
960
961 ROW_CHAIN_INIT(fill_info);
962 COL_CHAIN_INIT(fill_info);
963
964 fill_info->next = first_fill_info;
965 first_fill_info = fill_info;
966
967 fill_info->krn_blk_type = krn_blk_type;
968
969 /* First clone the OPERATOR_INFO structure again */
970 fill_info->op_info = *oinfo;
971
972 row_fcts_fast = fill_info->row_wquad_fast;
973 col_fcts_fast = fill_info->col_wquad_fast;
974
975 row_fcts = oinfo->row_fe_space->bas_fcts;
976 col_fcts = oinfo->col_fe_space->bas_fcts;
977
978 dim = row_fcts->dim;
979
980 if (op_type != OP_TYPE_SS) {
981 #undef MAT_BODY
982 #define MAT_BODY(F, CC, C, S, TYPE) \
983 fill_info->scl_el_mat = \
984 (void *)MAT_ALLOC(row_fcts->n_bas_fcts_max, \
985 col_fcts->n_bas_fcts_max, \
986 TYPE)
987 MAT_EMIT_BODY_SWITCH(fill_info->krn_blk_type);
988 fill_info->n_row_max = row_fcts->n_bas_fcts_max;
989 fill_info->n_col_max = col_fcts->n_bas_fcts_max;
990 }
991
992 fill_info->parametric = parametric;
993
994 el_mat_idx = 0;
995
996 if (fill_info->op_info.tangential)
997 el_mat_idx |= TANGENTIAL;
998 if (fill_info->op_info.init_element)
999 el_mat_idx |= OI_INIT_EL;
1000 if (parametric && not_all_param)
1001 el_mat_idx |= PARTPARAM;
1002 if (row_fcts != col_fcts)
1003 el_mat_idx |= MIXED;
1004 if (INIT_ELEMENT_NEEDED(col_fcts) || INIT_ELEMENT_NEEDED(row_fcts))
1005 el_mat_idx |= INIT_EL;
1006
1007 /* <<< 2nd/1st/0th order parsing */
1008
1009 row_fcts_fast[2] = col_fcts_fast[2] = NULL;
1010 init_row_fcts[2] = init_col_fcts[2] = 0x00;
1011 dflt_flags[2] = fast_flags[2] = 0x00;
1012
1013 if (fill_info->op_info.LALt.real_dd) {
1014 /* <<< LALt stuff */
1015
1016 el_mat_idx |= SCND_ORDER;
1017
1018 /* always the "mixed" case: one of the local Ansatz spaces belongs
1019 * to the neighbour element.
1020 */
1021 dflt_flags[2] |= WALL_FCT_MIX;
1022
1023 if (fill_info->op_info.tangential) {
1024 dflt_flags[2] |= WALL_FCT_TAN;
1025 init_row_fcts[2] |= INIT_TANGENTIAL;
1026 init_col_fcts[2] |= INIT_TANGENTIAL;
1027 }
1028
1029 init_row_fcts[2] = init_col_fcts[2] = INIT_GRD_PHI;
1030
1031 if (INIT_ELEMENT_NEEDED(fill_info->op_info.quad[2])) {
1032 el_mat_idx |= INIT_EL;
1033 }
1034
1035 if (fill_info->op_info.LALt_pw_const) {
1036 dflt_flags[2] |= WALL_FCT_PWC;
1037
1038 if(parametric && !not_all_param) {
1039 WARNING("You have selected piecewise constant LALt but seem to\n");
1040 WARNING("have a parametric mesh without affine elements!\n");
1041 }
1042 }
1043
1044 if(!oinfo->LALt_pw_const || parametric) {
1045 fast_flags[2] = dflt_flags[2];
1046 dflt_flags[2] &= ~WALL_FCT_PWC;
1047 }
1048
1049 /* >>> */
1050 }
1051
1052 row_fcts_fast[1] = col_fcts_fast[1] = NULL;
1053 init_row_fcts[1] = init_col_fcts[1] = 0x00;
1054 dflt_flags[1] = fast_flags[1] = 0x00;
1055
1056 if (fill_info->op_info.Lb0.real_dd) {
1057 /* <<< Lb0 stuff */
1058
1059 el_mat_idx |= FRST_ORDER|TANGENTIAL;
1060
1061 dflt_flags[1] |= WALL_FCT_MIX;
1062
1063 if (fill_info->op_info.tangential) {
1064 dflt_flags[1] |= WALL_FCT_TAN;
1065 init_row_fcts[1] |= INIT_TANGENTIAL;
1066 init_col_fcts[1] |= INIT_TANGENTIAL;
1067 }
1068
1069 init_row_fcts[1] |= INIT_PHI;
1070 init_col_fcts[1] |= INIT_GRD_PHI;
1071
1072 if (INIT_ELEMENT_NEEDED(fill_info->op_info.quad[1])) {
1073 el_mat_idx |= INIT_EL;
1074 }
1075
1076 if (oinfo->Lb0_pw_const) {
1077
1078 dflt_flags[1] |= WALL_FCT_PWC;
1079
1080 if(parametric && !not_all_param) {
1081 WARNING("You have selected piecewise constant Lb0 but seem to\n");
1082 WARNING("have a parametric mesh without affine elements!\n");
1083 }
1084 }
1085
1086 if(!oinfo->Lb0_pw_const || parametric) {
1087
1088 /* We will need the slow routines on some elements at least. */
1089 fast_flags[1] = dflt_flags[1];
1090 dflt_flags[1] &= ~WALL_FCT_PWC;
1091 }
1092
1093 /* >>> */
1094 }
1095
1096 if (fill_info->op_info.Lb1.real_dd) {
1097 /* <<< Lb1 stuff */
1098
1099 el_mat_idx |= FRST_ORDER|TANGENTIAL;
1100
1101 dflt_flags[1] |= WALL_FCT_MIX;
1102
1103 if (fill_info->op_info.tangential) {
1104 dflt_flags[1] |= WALL_FCT_TAN;
1105 init_row_fcts[1] |= INIT_TANGENTIAL;
1106 init_col_fcts[1] |= INIT_TANGENTIAL;
1107 }
1108
1109 init_row_fcts[1] |= INIT_GRD_PHI;
1110 init_col_fcts[1] |= INIT_PHI;
1111
1112 if (INIT_ELEMENT_NEEDED(fill_info->op_info.quad[1])) {
1113 el_mat_idx |= INIT_EL;
1114 }
1115
1116 if (oinfo->Lb1_pw_const) {
1117
1118 dflt_flags[1] |= WALL_FCT_PWC;
1119
1120 if(parametric && !not_all_param) {
1121 WARNING("You have selected piecewise constant Lb1 but seem to\n");
1122 WARNING("have a parametric mesh without affine elements!\n");
1123 }
1124 }
1125
1126 if(!oinfo->Lb1_pw_const || parametric) {
1127 /* We will need the slow routines on some elements at least. */
1128 fast_flags[1] = dflt_flags[1];
1129 dflt_flags[1] &= ~WALL_FCT_PWC;
1130 }
1131
1132 /* >>> */
1133 }
1134
1135 row_fcts_fast[0] = col_fcts_fast[0] = NULL;
1136 init_row_fcts[0] = init_col_fcts[0] = 0x00;
1137 dflt_flags[0] = fast_flags[0] = 0x00;
1138
1139 if (fill_info->op_info.c.real_dd) {
1140 /* <<< c stuff */
1141
1142 el_mat_idx |= ZERO_ORDER;
1143
1144 dflt_flags[0] |= WALL_FCT_MIX;
1145
1146 if (true || fill_info->op_info.tangential) {
1147 /* The zero-order term is always "tangential", i.e. only the
1148 * basis functions with non-vanishing trace on the respective
1149 * wall need to be assembled.
1150 */
1151 dflt_flags[0] |= WALL_FCT_TAN;
1152 init_row_fcts[0] |= INIT_TANGENTIAL;
1153 init_col_fcts[0] |= INIT_TANGENTIAL;
1154 el_mat_idx |= TANGENTIAL;
1155 }
1156
1157 init_row_fcts[0] = INIT_PHI;
1158 init_col_fcts[0] = INIT_PHI;
1159
1160 if (INIT_ELEMENT_NEEDED(fill_info->op_info.quad[0])) {
1161 el_mat_idx |= INIT_EL;
1162 }
1163
1164 if (oinfo->c_pw_const) {
1165
1166 dflt_flags[0] |= WALL_FCT_PWC;
1167
1168 if(parametric && !not_all_param) {
1169 WARNING("You have selected piecewise constant c but seem to\n");
1170 WARNING("have a parametric mesh without affine elements!\n");
1171 }
1172 }
1173
1174 if(!oinfo->c_pw_const || parametric) {
1175
1176 /* We will need the slow routines on some elements at least. */
1177 fast_flags[0] = dflt_flags[0];
1178 dflt_flags[0] &= ~WALL_FCT_PWC;
1179 }
1180
1181 /* >>> */
1182 }
1183
1184 /* >>> */
1185
1186 if (el_mat_idx & TANGENTIAL) {
1187 for (wall = 0; wall < N_WALLS(dim); ++wall) {
1188 fill_info->row_fcts_trace_map[wall] = row_fcts->trace_dof_map[0][0][wall];
1189 fill_info->n_trace_row_fcts[wall] = row_fcts->n_trace_bas_fcts[wall];
1190 }
1191 }
1192
1193 /* <<< aquire quad_fast structs */
1194
1195 if (fill_info->op_info.quad[0] == fill_info->op_info.quad[1]) {
1196 init_row_fcts[1] |= init_row_fcts[0];
1197 init_row_fcts[0] = 0x00;
1198 init_col_fcts[1] |= init_col_fcts[0];
1199 init_col_fcts[0] = 0x00;
1200 }
1201 if (fill_info->op_info.quad[1] == fill_info->op_info.quad[2]) {
1202 init_row_fcts[2] |= init_row_fcts[1];
1203 init_row_fcts[1] = 0x00;
1204 init_col_fcts[2] |= init_col_fcts[1];
1205 init_col_fcts[1] = 0x00;
1206 }
1207
1208 if (row_fcts == col_fcts) {
1209 for (i = 0; i < 3; i++) {
1210 init_row_fcts[i] |= init_col_fcts[i];
1211 if (init_row_fcts[i]) {
1212 fill_info->row_wquad_fast[i] = fill_info->col_wquad_fast[i] =
1213 get_wall_quad_fast(row_fcts,
1214 fill_info->op_info.quad[i], init_row_fcts[i]);
1215 }
1216 }
1217 } else {
1218 for (i = 0; i < 3; i++) {
1219 if (init_row_fcts[i]) {
1220 fill_info->row_wquad_fast[i] =
1221 get_wall_quad_fast(row_fcts,
1222 fill_info->op_info.quad[i], init_row_fcts[i]);
1223 }
1224 if (init_col_fcts[i]) {
1225 fill_info->col_wquad_fast[i] =
1226 get_wall_quad_fast(col_fcts,
1227 fill_info->op_info.quad[i], init_col_fcts[i]);
1228 }
1229 }
1230 }
1231
1232 /* >>> */
1233
1234 /* <<< fast/dflt cleanup */
1235
1236 if (dflt_flags[2]) {
1237 typeidx = type_index(fill_info->krn_blk_type, oinfo->LALt_type);
1238 for (wall = 0; wall < N_WALLS(dim); wall++) {
1239 fill_info->dflt_second_order[wall] =
1240 _AI_el_wall_fcts[op_type][typeidx][dim][wall][4][dflt_flags[2]];
1241 }
1242 }
1243 if (fast_flags[2]) {
1244 for (wall = 0; wall < N_WALLS(dim); wall++) {
1245 fill_info->fast_second_order[wall] =
1246 _AI_el_wall_fcts[op_type][typeidx][dim][wall][4][
1247 fast_flags[2]*N_TYPE_DEFUNS];
1248 }
1249 }
1250
1251 if (fill_info->op_info.Lb0.real_dd && fill_info->op_info.Lb1.real_dd) {
1252 wall_fcts_first = 3;
1253 } else if (fill_info->op_info.Lb1.real_dd) {
1254 wall_fcts_first = 2;
1255 } else if (fill_info->op_info.Lb1.real_dd) {
1256 wall_fcts_first = 1;
1257 } else {
1258 wall_fcts_first = -1;
1259 }
1260
1261 if (dflt_flags[1]) {
1262 typeidx = type_index(fill_info->krn_blk_type, oinfo->Lb_type);
1263 for (wall = 0; wall < N_WALLS(dim); wall++) {
1264 fill_info->dflt_first_order[wall] =
1265 _AI_el_wall_fcts[op_type][
1266 typeidx][dim][wall][wall_fcts_first][dflt_flags[1]];
1267 }
1268 }
1269 if (fast_flags[1]) {
1270 for (wall = 0; wall < N_WALLS(dim); wall++) {
1271 fill_info->fast_first_order[wall] =
1272 _AI_el_wall_fcts[op_type][
1273 typeidx][dim][wall][wall_fcts_first][fast_flags[1]];
1274 }
1275 }
1276
1277 if (dflt_flags[0]) {
1278 typeidx = type_index(fill_info->krn_blk_type, oinfo->c_type);
1279 for (wall = 0; wall < N_WALLS(dim); wall++) {
1280 fill_info->dflt_zero_order[wall] =
1281 _AI_el_wall_fcts[op_type][typeidx][dim][wall][0][dflt_flags[0]];
1282 }
1283 }
1284 if (fast_flags[1]) {
1285 for (wall = 0; wall < N_WALLS(dim); wall++) {
1286 fill_info->fast_zero_order[wall] =
1287 _AI_el_wall_fcts[op_type][typeidx][dim][wall][0][fast_flags[0]];
1288 }
1289 }
1290 /* >>> */
1291
1292 fill_info->neigh_el_mat_fcts = neigh_el_mat_fct_table[op_type][el_mat_idx];
1293
1294 TEST_EXIT(fill_info->neigh_el_mat_fcts != NULL,
1295 "Bogus choice for element matrix.\n");
1296
1297 /* >>> */
1298
1299 return fill_info;
1300 }
1301
1302 /* <<< AI_get_bndry_el_mat_fct() */
1303
AI_get_neigh_fill_info(const BNDRY_OPERATOR_INFO * oi_orig,MATENT_TYPE blk_type)1304 BNDRY_FILL_INFO *AI_get_neigh_fill_info(const BNDRY_OPERATOR_INFO *oi_orig,
1305 MATENT_TYPE blk_type)
1306 {
1307 /* FUNCNAME("AI_get_neigh_fill_info"); */
1308 BNDRY_OPERATOR_INFO oinfo[1];
1309 BNDRY_FILL_INFO *fill_info;
1310 const FE_SPACE *row_fe_space, *col_fe_space;
1311 PARAMETRIC *parametric = NULL;
1312 const WALL_QUAD_TENSOR *quad_tensor[3];
1313 int i;
1314
1315 row_fe_space = oi_orig->row_fe_space;
1316 col_fe_space = oi_orig->col_fe_space;
1317 if (col_fe_space == NULL) {
1318 col_fe_space = row_fe_space;
1319 }
1320
1321 for (i = 0; i < 3; i++) {
1322 quad_tensor[i] = oi_orig->quad_tensor[i];
1323 }
1324
1325 if (!unify_bop_info(
1326 oinfo, oi_orig, quad_tensor, row_fe_space, col_fe_space, blk_type)) {
1327 return NULL;
1328 }
1329
1330 /* <<< look for an existing fill_info object */
1331
1332 for (fill_info = first_fill_info; fill_info; fill_info = fill_info->next) {
1333
1334 /* standard tests */
1335 if (!OP_INFO_EQ_P(&fill_info->op_info, oinfo))
1336 continue;
1337 if (fill_info->parametric != parametric)
1338 continue;
1339
1340 /* Check whether the block-types match */
1341 if (fill_info->krn_blk_type != blk_type)
1342 continue;
1343
1344 #if 0
1345 /* Check for special bndry-operator stuff */
1346 if (BNDRY_FLAGS_CMP(fill_info->op_info.bndry_type, oinfo->bndry_type) != 0)
1347 continue;
1348 #endif
1349 if (fill_info->op_info.tangential != oinfo->tangential)
1350 continue;
1351
1352 /* All tests passed, but still the application-data pointer could
1353 * be different. We remember the first matching operator info with
1354 * possibly non-matching application-data pointer and clone that
1355 * if we do not find a fully matching fill-info structure.
1356 */
1357 if (fill_info->op_info.user_data != oinfo->user_data) {
1358 continue;
1359 }
1360
1361 break;
1362 }
1363
1364 if (fill_info != NULL) {
1365 return fill_info;
1366 }
1367
1368 /* >>> */
1369
1370 oinfo->row_fe_space = row_fe_space = copy_fe_space(row_fe_space);
1371 oinfo->col_fe_space = col_fe_space = copy_fe_space(col_fe_space);
1372
1373 fill_info = __get_neigh_fill_info(oinfo, blk_type);
1374
1375 fill_info->el_mat = get_el_matrix(row_fe_space, col_fe_space, blk_type);
1376
1377 if (!CHAIN_SINGLE(row_fe_space) || !CHAIN_SINGLE(col_fe_space)) {
1378 const FE_SPACE *col_fesp, *row_fesp;
1379 BNDRY_FILL_INFO *row_chain, *col_chain, *chain_info;
1380 EL_MATRIX *elm;
1381
1382 /* See whether we need to allocate a block-operator */
1383 elm = fill_info->el_mat;
1384 row_fesp = row_fe_space;
1385 col_chain = fill_info;
1386 CHAIN_FOREACH(col_fesp, col_fe_space, const FE_SPACE) {
1387 elm = ROW_CHAIN_NEXT(elm, EL_MATRIX);
1388 for (i = 0; i < 3; i++) {
1389 if (quad_tensor[i]) {
1390 quad_tensor[i] =
1391 ROW_CHAIN_NEXT(quad_tensor[i], const WALL_QUAD_TENSOR);
1392 }
1393 }
1394 unify_bop_info(oinfo, oi_orig, quad_tensor, row_fesp, col_fesp, blk_type);
1395 row_chain = __get_neigh_fill_info(oinfo, blk_type);
1396 ROW_CHAIN_ADD_TAIL(col_chain, row_chain);
1397 row_chain->el_mat = elm;
1398 }
1399 /* reset to list-head */
1400 elm = fill_info->el_mat;
1401 col_fesp = col_fe_space;
1402 row_chain = fill_info;
1403 CHAIN_FOREACH(row_fesp, row_fe_space, const FE_SPACE) {
1404 elm = COL_CHAIN_NEXT(elm, EL_MATRIX);
1405 for (i = 0; i < 3; i++) {
1406 if (quad_tensor[i]) {
1407 quad_tensor[i] =
1408 COL_CHAIN_NEXT(quad_tensor[i], const WALL_QUAD_TENSOR);
1409 }
1410 }
1411 unify_bop_info(oinfo, oi_orig, quad_tensor, row_fesp, col_fesp, blk_type);
1412 col_chain = __get_neigh_fill_info(oinfo, blk_type);
1413 COL_CHAIN_ADD_TAIL(row_chain, col_chain);
1414 col_chain->el_mat = elm;
1415 CHAIN_FOREACH(col_fesp, col_fe_space, const FE_SPACE) {
1416 elm = ROW_CHAIN_NEXT(elm, EL_MATRIX);
1417 row_chain = ROW_CHAIN_NEXT(row_chain, BNDRY_FILL_INFO);
1418 for (i = 0; i < 3; i++) {
1419 if (quad_tensor[i]) {
1420 quad_tensor[i] =
1421 ROW_CHAIN_NEXT(quad_tensor[i], const WALL_QUAD_TENSOR);
1422 }
1423 }
1424 unify_bop_info(
1425 oinfo, oi_orig, quad_tensor, row_fesp, col_fesp, blk_type);
1426 chain_info = __get_neigh_fill_info(oinfo, blk_type);
1427 ROW_CHAIN_ADD_TAIL(col_chain, chain_info);
1428 COL_CHAIN_ADD_TAIL(row_chain, chain_info);
1429 }
1430 /* roll over to the list head */
1431 elm = ROW_CHAIN_NEXT(elm, EL_MATRIX);
1432 row_chain = ROW_CHAIN_NEXT(row_chain, BNDRY_FILL_INFO);
1433 for (i = 0; i < 3; i++) {
1434 if (quad_tensor[i]) {
1435 quad_tensor[i] =
1436 ROW_CHAIN_NEXT(quad_tensor[i], const WALL_QUAD_TENSOR);
1437 }
1438 }
1439 }
1440 }
1441
1442 return fill_info;
1443 }
1444
1445 /* >>> */
1446
1447 #endif /* EMIT_SS_VERSIONS */
1448