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: bas_fct.c */
7 /* */
8 /* description: collecting information of all Lagrange elements */
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 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
26 /* */
27 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
28 /* */
29 /*--------------------------------------------------------------------------*/
30
31 #include "alberta_intern.h"
32 #include "alberta.h"
33
34 /* Macro for a general scalar interpolation function for Lagrange basis
35 * functions.
36 *
37 * This new version only allows "local" functions. The complexity of
38 * parametric->init_element() etc. is moved to the high level code.
39 *
40 * If "wall >= 0" then the numbering of the basis function in the
41 * returned coefficient array is according to BAS_FCTS::trace_dof_map.
42 */
43 #define GENERATE_INTERPOL(PREFIX, DEGREE, DIM, N_BAS_FCT) \
44 static void PREFIX##interpol##DEGREE##_##DIM##d( \
45 EL_REAL_VEC *el_vec, \
46 const EL_INFO *el_info, int wall, \
47 int no, const int *b_no, \
48 LOC_FCT_AT_QP f, void *f_data, \
49 const BAS_FCTS *thisptr) \
50 { \
51 FUNCNAME(#PREFIX"interpol"#DEGREE"_"#DIM"d"); \
52 LAGRANGE_DATA *ld = &lag_##DEGREE##_##DIM##d_data; \
53 REAL *rvec = el_vec->vec; \
54 const QUAD *lq; \
55 const int *trace_map; \
56 int i; \
57 \
58 DEBUG_TEST_EXIT(ld->lumping_quad != NULL, \
59 "called for uninitialized Lagrange basis functions\n"); \
60 \
61 if (wall < 0) { \
62 lq = ld->lumping_quad; \
63 trace_map = NULL; \
64 } else { \
65 int type = el_info->el_type > 0; \
66 int orient = el_info->orientation < 0; \
67 lq = &ld->trace_lumping_quad[type][orient][wall]; \
68 trace_map = thisptr->trace_dof_map[type][orient][wall]; \
69 } \
70 \
71 DEBUG_TEST_EXIT(!b_no || (no >= 0 && no <= lq->n_points), \
72 "not for %d points\n", no); \
73 \
74 el_vec->n_components = thisptr->n_bas_fcts; \
75 \
76 if (b_no) { \
77 for (i = 0; i < no; i++) { \
78 int ib = wall < 0 ? b_no[i] : trace_map[b_no[i]]; \
79 rvec[ib] = f(el_info, lq, b_no[i], f_data); \
80 } \
81 } else { \
82 for (i = 0; i < lq->n_points; i++) { \
83 int ib = wall < 0 ? i : trace_map[i]; \
84 rvec[ib] = f(el_info, lq, i, f_data); \
85 } \
86 } \
87 } \
88 struct _AI_semicolon_dummy
89
90 /* Macro for a general vector interpolation function (Lagrange basis fcts) */
91
92 #define GENERATE_INTERPOL_D(PREFIX, DEGREE, DIM, N_BAS_FCT) \
93 static void PREFIX##interpol_d_##DEGREE##_##DIM##d( \
94 EL_REAL_D_VEC *el_vec, \
95 const EL_INFO *el_info, \
96 int wall, \
97 int no, const int *b_no, \
98 LOC_FCT_D_AT_QP f, void *f_data, \
99 const BAS_FCTS *thisptr) \
100 { \
101 FUNCNAME("interpol_d"#DEGREE"_"#DIM"d"); \
102 LAGRANGE_DATA *ld = &lag_##DEGREE##_##DIM##d_data; \
103 const QUAD *lq; \
104 REAL_D *rvec = el_vec->vec; \
105 const int *trace_map; \
106 int i; \
107 \
108 DEBUG_TEST_EXIT(ld->lumping_quad != NULL, \
109 "called for uninitialized Lagrange basis functions\n"); \
110 \
111 if (wall < 0) { \
112 lq = ld->lumping_quad; \
113 trace_map = NULL; \
114 } else { \
115 int type = el_info->el_type > 0; \
116 int orient = el_info->orientation < 0; \
117 lq = &ld->trace_lumping_quad[type][orient][wall]; \
118 trace_map = thisptr->trace_dof_map[type][orient][wall]; \
119 } \
120 \
121 DEBUG_TEST_EXIT(!b_no || (no >= 0 && no <= lq->n_points), \
122 "not for %d points\n", no); \
123 \
124 el_vec->n_components = thisptr->n_bas_fcts; \
125 \
126 if (b_no) { \
127 for (i = 0; i < no; i++) { \
128 int ib = wall < 0 ? b_no[i] : trace_map[b_no[i]]; \
129 f(rvec[ib], el_info, lq, b_no[i], f_data); \
130 } \
131 } else { \
132 for (i = 0; i < lq->n_points; i++) { \
133 int ib = wall < 0 ? i : trace_map[i]; \
134 f(rvec[ib], el_info, lq, i, f_data); \
135 } \
136 } \
137 } \
138 struct _AI_semicolon_dummy
139
140 /* Macro for a more general vector interpolation function (Lagrange
141 * basis fcts _ONLY_).
142 */
143
144 #define GENERATE_INTERPOL_DOW(PREFIX, DEGREE, DIM, N_BAS_FCT) \
145 static void PREFIX##interpol_dow_##DEGREE##_##DIM##d( \
146 EL_REAL_VEC_D *vec, \
147 const EL_INFO *el_info, \
148 int wall, \
149 int no, const int *b_no, \
150 LOC_FCT_D_AT_QP f, void *f_data, \
151 const BAS_FCTS *thisptr) \
152 { \
153 PREFIX##interpol_d_##DEGREE##_##DIM##d( \
154 (EL_REAL_D_VEC *)vec, el_info, wall, no, b_no, f, f_data, thisptr); \
155 } \
156 struct _AI_semicolon_dummy
157
158 typedef struct lagrange_data {
159 const REAL_B *nodes; /* Lagrange nodes */
160 const QUAD *lumping_quad; /* Lumping quadrature on element ... */
161 QUAD trace_lumping_quad[2][2][N_WALLS_MAX]; /* ... and on walls */
162 } LAGRANGE_DATA;
163
164 typedef struct ortho_data {
165 const QUAD_FAST *qfast; /* Quadrature for interpolation */
166 } ORTHO_DATA;
167
168 typedef const EL_REAL_VEC_D *
169 (*GET_REAL_VEC_D_TYPE)(REAL *, const EL *, const DOF_REAL_VEC_D *);
170
171 #define INIT_BFCTS_CHAIN(bstruct) CHAIN_INITIALIZER(bstruct), &bstruct
172
173 #define LAGRANGE_NAME(deg) "lagrange"CPP_STRINGIZE(deg)
174 #define LAG_NAME_LEN(deg) (sizeof(LAGRANGE_NAME(deg))-1)
175
176 #define DISC_LAGRANGE_NAME(deg) "disc_"LAGRANGE_NAME(deg)
177 #define DISC_LAG_NAME_LEN(deg) (sizeof(DISC_LAGRANGE_NAME(deg))-1)
178
179 #define DISC_ORTHO_NAME(deg) "disc_ortho_poly"CPP_STRINGIZE(deg)
180 #define DISC_ORTHO_NAME_LEN(deg) (sizeof(DISC_ORTHO_NAME(deg))-1)
181
182 #include "bas_fct_0d.c"
183 #include "bas_fct_1d.c"
184 #if DIM_MAX > 1
185 #include "bas_fct_2d.c"
186 #if DIM_MAX > 2
187 #include "bas_fct_3d.c"
188 #endif
189 #endif
190
191 #undef GENERATE_INTERPOL
192 #undef GENERATE_INTERPOL_D
193
194 typedef struct bfcts_node BFCTS_NODE;
195 struct bfcts_node
196 {
197 const BAS_FCTS *bas_fcts;
198 size_t name_len;/* The length of the _relevant_ part of the name */
199 BFCTS_NODE *next;
200 };
201
lagrange_lumping_quadrature(const BAS_FCTS * bfcts)202 const QUAD *lagrange_lumping_quadrature(const BAS_FCTS *bfcts)
203 {
204 const QUAD *wq = get_quadrature(bfcts->dim, bfcts->degree);
205 QUAD *lq;
206 int ib, iq;
207 REAL *lw; /* lumping weights */
208 char *name;
209
210 lq = MEM_CALLOC(1, QUAD);
211 lq->w = lw = MEM_CALLOC(bfcts->n_bas_fcts, REAL);
212
213 lq->name = name = MEM_ALLOC(sizeof("Lagrange X Xd Lumping Quadrature"), char);
214 sprintf(name, "Lagrange %d %dd Lumping Quadrature",
215 bfcts->degree, bfcts->dim);
216
217 lq->degree = bfcts->degree;
218 lq->dim = bfcts->dim;
219 lq->codim = 0;
220 lq->subsplx = -1;
221 lq->n_points =
222 lq->n_points_max = bfcts->n_bas_fcts;
223 lq->lambda = LAGRANGE_NODES(bfcts);
224 lq->w = lw;
225
226 /* Now compute the weights using one sufficiently good quadrature
227 * formula s.t. the integrals are exact for all basis functions (and
228 * thus for the entire space of polynomials up to
229 * bfcts->degree). This means that the weights are simply the
230 * integral of the basis functions taken over the standard simplex.
231 */
232 for (ib = 0; ib < bfcts->n_bas_fcts; ib++) {
233 for (iq = 0; iq < wq->n_points; iq++) {
234 lw[ib] += wq->w[iq]*PHI(bfcts, ib, wq->lambda[iq]);
235 }
236 }
237
238 /* Allocate the meta-data structure, but do not register the
239 * quadrature for general use.
240 */
241 register_quadrature(lq);
242
243 return lq;
244 }
245
246 /*******************************************************************************
247 * List of all discontinuous orthogonal basis functions.
248 ******************************************************************************/
249 #define MAX_DEG 2
250 static BFCTS_NODE all_disc_ortho_poly[DIM_MAX+1][MAX_DEG] = {
251 { { NULL, 0, NULL } },
252 #if DIM_MAX > 0 /* always */
253 { { &disc_ortho1_1d, DISC_ORTHO_NAME_LEN(1), &all_disc_ortho_poly[1][1] },
254 { &disc_ortho2_1d, DISC_ORTHO_NAME_LEN(2), NULL }, },
255 #endif
256 #if DIM_MAX > 1
257 { { &disc_ortho1_2d, DISC_ORTHO_NAME_LEN(1), &all_disc_ortho_poly[2][1] },
258 { &disc_ortho2_2d, DISC_ORTHO_NAME_LEN(2), NULL }, },
259 #endif
260 #if DIM_MAX > 2
261 { { &disc_ortho1_3d, DISC_ORTHO_NAME_LEN(1), &all_disc_ortho_poly[3][1] },
262 { &disc_ortho2_3d, DISC_ORTHO_NAME_LEN(2), NULL }, },
263 #endif
264 };
265
get_disc_ortho_poly(int dim,int degree)266 const BAS_FCTS *get_disc_ortho_poly(int dim, int degree)
267 {
268 FUNCNAME("get_disc_ortho_poly");
269 const BAS_FCTS *bfcts;
270 ORTHO_DATA *od;
271
272 if ((dim < 0 || dim > DIM_MAX)) {
273 ERROR("Discontinuous orthogonal basis functions of dimension %d are "
274 "not available for DIM_MAX == %d!\n", dim, DIM_MAX);
275 return NULL;
276 }
277
278 if (degree < 0 || degree > MAX_DEG) {
279 ERROR("Discontinuous orthogonal basis functions of degree %d are "
280 "not available\n", degree);
281 return NULL;
282 }
283
284 if (dim == 0) {
285 bfcts = get_lagrange(0, 4);
286 } else if (degree == 0) {
287 bfcts = get_lagrange(dim, degree);
288 } else {
289 bfcts = all_disc_ortho_poly[dim][degree-1].bas_fcts;
290 }
291
292 od = (ORTHO_DATA *)bfcts->ext_data;
293
294 if (od->qfast == NULL) {
295 od->qfast = get_quad_fast(bfcts, get_quadrature(dim, 2*degree), INIT_PHI);
296 }
297
298 return bfcts;
299 }
300
301 /*******************************************************************************
302 * List of all discontinuous Lagrange basisfunctions.
303 ******************************************************************************/
304 #define MAX_DEG 2
305
306 static BFCTS_NODE all_disc_lagrange[DIM_MAX+1][MAX_DEG+1] = {
307 { { NULL, 0, NULL } },
308 #if DIM_MAX > 0 /* always */
309 { { &disc_lagrange0_1d, DISC_LAG_NAME_LEN(0), &all_disc_lagrange[1][1] },
310 { &disc_lagrange1_1d, DISC_LAG_NAME_LEN(1), &all_disc_lagrange[1][2] },
311 { &disc_lagrange2_1d, DISC_LAG_NAME_LEN(2), &all_disc_ortho_poly[1][0] }, },
312 #endif
313 #if DIM_MAX > 1
314 { { &disc_lagrange0_2d, DISC_LAG_NAME_LEN(0), &all_disc_lagrange[2][1] },
315 { &disc_lagrange1_2d, DISC_LAG_NAME_LEN(1), &all_disc_lagrange[2][2] },
316 { &disc_lagrange2_2d, DISC_LAG_NAME_LEN(2), &all_disc_ortho_poly[2][0] }, },
317 #endif
318 #if DIM_MAX > 2
319 { { &disc_lagrange0_3d, DISC_LAG_NAME_LEN(0), &all_disc_lagrange[3][1] },
320 { &disc_lagrange1_3d, DISC_LAG_NAME_LEN(1), &all_disc_lagrange[3][2] },
321 { &disc_lagrange2_3d, DISC_LAG_NAME_LEN(2), &all_disc_ortho_poly[3][0] }, },
322 #endif
323 };
324
get_discontinuous_lagrange(int dim,int degree)325 const BAS_FCTS *get_discontinuous_lagrange(int dim, int degree)
326 {
327 FUNCNAME("get_discontinuous_lagrange");
328 const BAS_FCTS *bfcts;
329 LAGRANGE_DATA *ld, *tr_ld;
330
331 if (dim < 0 || dim > DIM_MAX) {
332 ERROR("Discontinuous Lagrange basis functions of dimension %d are "
333 "not available for DIM_MAX == %d!\n", dim, DIM_MAX);
334 return NULL;
335 }
336
337 if (degree < 0 || degree > MAX_DEG) {
338 ERROR("Discontinuous Lagrange basis functions of degree %d are "
339 "not available\n", degree);
340 return NULL;
341 }
342
343 if (dim == 0) {
344 bfcts = get_lagrange(0, 4);
345 } else {
346 bfcts = all_disc_lagrange[dim][degree].bas_fcts;
347 }
348
349 ld = (LAGRANGE_DATA *)bfcts->ext_data;
350 if (ld->lumping_quad == NULL) {
351 ld->lumping_quad = lagrange_lumping_quadrature(bfcts);
352 if (dim >= 1) {
353 REAL_B *lambda;
354 QUAD *lq, *tlq;
355 int o, t, w;
356 /* recursively also initialize the trace-space quadratures. */
357 get_discontinuous_lagrange(dim-1, degree);
358 tr_ld = (LAGRANGE_DATA *)bfcts->trace_bas_fcts->ext_data;
359 lq = (QUAD *)ld->lumping_quad;
360 for (t = 0; t <= (dim > 2); t++) {
361 for (o = 0; o <= (dim > 2); o++) {
362 for (w = 0; w < N_WALLS(dim); w++) {
363 tlq = &ld->trace_lumping_quad[t][o][w];
364 *tlq = *tr_ld->lumping_quad;
365 tlq->codim = 1;
366 tlq->subsplx = w;
367 tlq->lambda = (const REAL_B *)
368 (lambda = MEM_CALLOC(tr_ld->lumping_quad->n_points, REAL_B));
369 if (degree == 0) {
370 /* Rather choose the mid-point on the respective wall;
371 * degree 0 is special in that there is -- in principle
372 * -- no well defined Lagrange-node.
373 */
374 int i;
375 for (i = 0; i < w; i++) {
376 lambda[0][i] = 1.0/(REAL)N_LAMBDA(dim);
377 }
378 for (++i; i < N_LAMBDA(dim); i++) {
379 lambda[0][i] = 1.0/(REAL)N_LAMBDA(dim);
380 }
381 } else {
382 int b;
383 for (b = 0; b < N_BAS_LAGRANGE(degree, dim-1); b++) {
384 COPY_BAR(dim,
385 lq->lambda[bfcts->trace_dof_map[t][o][w][b]],
386 lambda[b]);
387 }
388 }
389 }
390 }
391 }
392 }
393 }
394
395 return bfcts;
396 }
397
398 #undef MAX_DEG
399
400 /*--------------------------------------------------------------------------*/
401 /* linked list of all used basis functions: Lagrange basisfunctions are */
402 /* always members of the list */
403 /*--------------------------------------------------------------------------*/
404 #define MAX_DEG 4
405
406 static BFCTS_NODE all_lagrange[DIM_MAX+1][MAX_DEG+1] = {
407 #if DIM_MAX >= 0 /* ;) */
408 { { &lagrange_0d, sizeof("lagrange")-1, NULL }, },
409 #endif
410 #if DIM_MAX >= 1
411 { { &lagrange1_1d, LAG_NAME_LEN(1), &all_lagrange[1][1] },
412 { &lagrange2_1d, LAG_NAME_LEN(2), &all_lagrange[1][2] },
413 { &lagrange3_1d, LAG_NAME_LEN(3), &all_lagrange[1][3] },
414 { &lagrange4_1d, LAG_NAME_LEN(4), &all_disc_lagrange[1][0] }, },
415 #endif
416 #if DIM_MAX >= 2
417 { { &lagrange1_2d, LAG_NAME_LEN(1), &all_lagrange[2][1] },
418 { &lagrange2_2d, LAG_NAME_LEN(2), &all_lagrange[2][2] },
419 { &lagrange3_2d, LAG_NAME_LEN(3), &all_lagrange[2][3] },
420 { &lagrange4_2d, LAG_NAME_LEN(4), &all_disc_lagrange[2][0] }, },
421 #endif
422 #if DIM_MAX >= 3
423 { { &lagrange1_3d, LAG_NAME_LEN(1), &all_lagrange[3][1] },
424 { &lagrange2_3d, LAG_NAME_LEN(2), &all_lagrange[3][2] },
425 { &lagrange3_3d, LAG_NAME_LEN(3), &all_lagrange[3][3] },
426 { &lagrange4_3d, LAG_NAME_LEN(4), &all_disc_lagrange[3][0] }, },
427 #endif
428 };
429
430 int n_bas_fcts_max[] = {
431 N_BAS_LAG_0D,
432 N_BAS_LAG_4_1D,
433 #if DIM_MAX > 1
434 N_BAS_LAG_4_2D,
435 # if DIM_MAX > 2
436 N_BAS_LAG_4_3D,
437 # if DIM_MAX > 3
438 # error FIXME!
439 # endif
440 # endif
441 #endif
442 };
443
get_lagrange(int dim,int degree)444 const BAS_FCTS *get_lagrange(int dim, int degree)
445 {
446 FUNCNAME("get_lagrange");
447 const BAS_FCTS *bfcts;
448 LAGRANGE_DATA *ld, *tr_ld;
449
450 if (degree == 0) {
451 return get_discontinuous_lagrange(dim, degree);
452 }
453
454 if (dim < 0 || dim > DIM_MAX) {
455 #if ALBERTA_DEBUG
456 WARNING("Lagrange basis functions of dimension %d are "
457 "not available for DIM_MAX == %d!\n", dim, DIM_MAX);
458 #endif
459 return NULL;
460 }
461
462 if (degree < 1 || degree > MAX_DEG) {
463 #if ALBERTA_DEBUG
464 WARNING("no lagrangian basis functions of degree %d\n", degree);
465 #endif
466 return NULL;
467 }
468
469 if (dim == 0) {
470 degree = 1;
471 }
472
473 bfcts = all_lagrange[dim][degree-1].bas_fcts;
474
475 ld = (LAGRANGE_DATA *)bfcts->ext_data;
476 if (ld->lumping_quad == NULL) {
477 ld->lumping_quad = lagrange_lumping_quadrature(bfcts);
478 if (dim >= 1) {
479 REAL_B *lambda;
480 QUAD *lq, *tlq;
481 int o, t, w;
482 /* recursively also initialize the trace-space quadratures. */
483 get_lagrange(dim-1, degree);
484 tr_ld = (LAGRANGE_DATA *)bfcts->trace_bas_fcts->ext_data;
485 lq = (QUAD *)ld->lumping_quad;
486 for (t = 0; t <= (dim > 2); t++) {
487 for (o = 0; o <= (dim > 2); o++) {
488 for (w = 0; w < N_WALLS(dim); w++) {
489 tlq = &ld->trace_lumping_quad[t][o][w];
490 *tlq = *tr_ld->lumping_quad;
491 tlq->codim = 1;
492 tlq->subsplx = w;
493 tlq->lambda = (const REAL_B *)
494 (lambda = MEM_CALLOC(tr_ld->lumping_quad->n_points, REAL_B));
495 if (degree == 0) {
496 /* Rather choose the mid-point on the respective wall;
497 * degree 0 is special in that there is -- in principle
498 * -- no well defined Lagrange-node.
499 */
500 int i;
501 for (i = 0; i < w; i++) {
502 lambda[0][i] = 1.0/(REAL)N_LAMBDA(dim);
503 }
504 for (++i; i < N_LAMBDA(dim); i++) {
505 lambda[0][i] = 1.0/(REAL)N_LAMBDA(dim);
506 }
507 } else {
508 int b;
509 for (b = 0; b < N_BAS_LAGRANGE(degree, dim-1); b++) {
510 COPY_BAR(dim,
511 lq->lambda[bfcts->trace_dof_map[t][o][w][b]],
512 lambda[b]);
513 }
514 }
515 }
516 }
517 }
518 }
519 }
520
521 return bfcts;
522 }
523
524 /* Default get_FOOBAR_vec() implementation */
525 #define EQ_COPY(src, dst) (dst) = (src)
526 #define DEFUN_GET_VEC(VECNAME, COPY) \
527 const EL_##VECNAME##_VEC * \
528 CPP_CONCAT3(default_get_, VECNAME##_name, _vec)( \
529 VECNAME##_VEC_TYPE *vec, \
530 const EL *el, \
531 const DOF_##VECNAME##_VEC *dof_vec) \
532 { \
533 FUNCNAME("get_"CPP_STRINGIZE(VECNAME##_name)"_vec"); \
534 VECNAME##_VEC_TYPE *rvec = vec ? vec : dof_vec->vec_loc->vec; \
535 const BAS_FCTS *bas_fcts = dof_vec->fe_space->bas_fcts; \
536 int n_bas_fcts = bas_fcts->n_bas_fcts; \
537 DOF index[n_bas_fcts]; \
538 int i; \
539 \
540 GET_DOF_INDICES(dof_vec->fe_space->bas_fcts, \
541 el, dof_vec->fe_space->admin, index); \
542 \
543 for (i = 0; i < n_bas_fcts; i++) { \
544 COPY(dof_vec->vec[index[i]], rvec[i]); \
545 } \
546 \
547 return vec ? NULL : dof_vec->vec_loc; \
548 } \
549 struct _AI_semicolon_dummy
550
551 DEFUN_GET_VEC(INT, EQ_COPY);
552 DEFUN_GET_VEC(REAL, EQ_COPY);
553 DEFUN_GET_VEC(REAL_D, COPY_DOW);
554 DEFUN_GET_VEC(REAL_DD, _AI_MCOPY_DOW);
555 DEFUN_GET_VEC(UCHAR, EQ_COPY);
556 DEFUN_GET_VEC(SCHAR, EQ_COPY);
557 DEFUN_GET_VEC(PTR, EQ_COPY);
558
559 const EL_REAL_VEC_D *
default_get_real_vec_d(REAL vec[],const EL * el,const DOF_REAL_VEC_D * dof_vec)560 default_get_real_vec_d(REAL vec[],
561 const EL *el,
562 const DOF_REAL_VEC_D *dof_vec)
563 {
564 if (dof_vec->stride != 1) {
565 return (EL_REAL_VEC_D *)default_get_real_d_vec(
566 (REAL_D *)vec, el, (const DOF_REAL_D_VEC *)dof_vec);
567 } else {
568 return (EL_REAL_VEC_D *)default_get_real_vec(
569 vec, el, (const DOF_REAL_VEC *)dof_vec);
570 }
571 }
572
573 /*--------------------------------------------------------------------------*/
574 /* add a set of new basis functions to the list; return true, if possible */
575 /* else false */
576 /*--------------------------------------------------------------------------*/
577
578 static struct bfcts_node *bas_fcts_list[DIM_MAX+1] = {
579 #if DIM_MAX >= 0
580 &all_lagrange[0][0],
581 #endif
582 #if DIM_MAX >= 1
583 &all_lagrange[1][0],
584 #endif
585 #if DIM_MAX >= 2
586 &all_lagrange[2][0],
587 #endif
588 #if DIM_MAX >= 3
589 &all_lagrange[3][0],
590 #endif
591 };
592
593 /* The return value is the old definition, if basis-functions of the
594 * same name already exist. The functions ignores an options _Xd
595 * suffix (e.g.: "lagrange_3_2d", only the "lagrange_3" part of the
596 * name is taken into account).
597 */
new_bas_fcts(const BAS_FCTS * bas_fcts)598 const BAS_FCTS *new_bas_fcts(const BAS_FCTS * bas_fcts)
599 {
600 FUNCNAME("new_bas_fcts");
601 BFCTS_NODE *bf_node;
602 char suffix[3] = { '_', 'X', 'd' };
603 size_t len;
604 int dim;
605
606 if (!bas_fcts) {
607 ERROR("no basis functions specified; bas_fcts pointer to NULL\n");
608 return(0);
609 }
610
611 TEST_EXIT(bas_fcts->name,
612 "new basis functions must have name; bas_fcts->name pointer to NULL\n");
613
614 TEST_EXIT(strlen(bas_fcts->name),
615 "new basis functions must have a non empty name\n");
616 TEST_EXIT(bas_fcts->dim >= 0 && bas_fcts->dim <= DIM_MAX,
617 "new basis functions must have a dimension between 1 and %d\n",
618 DIM_MAX);
619 dim = bas_fcts->dim;
620 if (dim >= 1) {
621 TEST_EXIT(bas_fcts->trace_bas_fcts != NULL,
622 "new basis functions must define their trace-space.\n");
623 new_bas_fcts(bas_fcts->trace_bas_fcts);
624 }
625 TEST_EXIT(bas_fcts->rdim == 1 || bas_fcts->rdim == DIM_OF_WORLD,
626 "Rand dimension must be either 1 or DIM_OF_WORLD.\n");
627 TEST_EXIT(bas_fcts->degree >= 0,
628 "new basis functions must have a non negative quadrature degree\n");
629 if (bas_fcts->n_bas_fcts > 0) { /* allow for the NULL space */
630 TEST_EXIT(bas_fcts->phi,
631 "new basis functions: phi not set\n");
632 TEST_EXIT(bas_fcts->grd_phi,
633 "new basis functions: grd_phi not set\n");
634 TEST_EXIT(bas_fcts->rdim == 1 || bas_fcts->phi_d != NULL,
635 "new basis functions: rdim == DIM_OF_WORLD, "
636 "but phi_d == NULL.\n");
637 TEST(bas_fcts->D2_phi,
638 "Warning: new basis functions: D2_phi not set\n");
639 }
640 TEST_EXIT(bas_fcts->get_dof_indices,
641 "new basis functions: get_dof_indices not set\n");
642 TEST_EXIT(bas_fcts->get_bound,
643 "new basis functions: get_bound not set\n");
644 if (bas_fcts->rdim == 1) {
645 TEST(bas_fcts->interpol,
646 "Warning: new basis functions \"%s\": interpol not set\n",
647 bas_fcts->name);
648 TEST(bas_fcts->interpol_d,
649 "Warning: new basis functions \"%s\": interpol_d not set\n",
650 bas_fcts->name);
651 }
652 TEST(bas_fcts->interpol_dow,
653 "Warning: new basis functions: interpol_dow not set\n");
654
655 TEST_EXIT(bas_fcts->n_bas_fcts <= bas_fcts->n_bas_fcts_max,
656 "Error: n_bas_fcts must be < n_bas_fcts_max.\n");
657
658 dim = bas_fcts->dim;
659 suffix[1] = '0' + dim;
660 len = strlen(bas_fcts->name);
661 if (memcmp(bas_fcts->name+len-3, suffix, 3) == 0) {
662 len -= 3; /* ignore the dimension part */
663 }
664 for (bf_node = bas_fcts_list[dim]; bf_node; bf_node = bf_node->next) {
665 if (len == bf_node->name_len &&
666 !strncmp(bas_fcts->name, bf_node->bas_fcts->name, len)) {
667 const BAS_FCTS *old_bf = bf_node->bas_fcts;
668 if (bas_fcts != old_bf) {
669 WARNING("pointer to new and existing basis functions differ %p!=%p\n",
670 bas_fcts, old_bf);
671 WARNING("overriding old definition.\n");
672 }
673 bf_node->bas_fcts = bas_fcts;
674 return old_bf;
675 }
676 }
677
678 bf_node = MEM_ALLOC(1, struct bfcts_node);
679 bf_node->bas_fcts = bas_fcts;
680 bf_node->name_len = len; /* length discarding the "_Xd" suffix */
681 bf_node->next = bas_fcts_list[dim];
682 bas_fcts_list[dim] = bf_node;
683
684 n_bas_fcts_max[dim] = MAX(n_bas_fcts_max[dim], bas_fcts->n_bas_fcts_max);
685
686 return NULL; /* no old definition */
687 }
688
689 /*--------------------------------------------------------------------------*/
690 /* get a pointer to a set of basis functions from the list; identifier is */
691 /* the name of the basis functions; */
692 /* returns a pointer to the BAS_FCTS structure if a corresponding set was */
693 /* found in the list, else pointer to NULL */
694 /*--------------------------------------------------------------------------*/
695
696 #if HAVE_LIBLTDL && HAVE_LTDL_H
697 # include <ltdl.h>
698 #endif
699
700 #if 0
701 typedef struct bfcts_metadata
702 {
703 EL_DOF_VEC *dof;
704 EL_BNDRY_VEC *bndry;
705 EL_REAL_VEC *coeff;
706 EL_REAL_D_VEC *coeff_d;
707 EL_INT_VEC *int_vec;
708 EL_REAL_VEC *real_vec;
709 EL_REAL_D_VEC *real_d_vec;
710 EL_UCHAR_VEC *uchar_vec;
711 EL_SCHAR_VEC *schar_vec;
712 EL_PTR_VEC *ptr_vec;
713 } BFCTS_METADATA;
714
715 /* Allocate a suitable meta-data object (storage for results etc.) */
716 const BFCTS_METADATA *get_bfcts_meta_data(const BAS_FCTS *bfcts)
717 {
718 BFCTS_METADATA *md;
719
720 md = MEM_ALLOC(1, BFCTS_METADATA);
721 md->dof = get_el_dof_vec(bfcts);
722 md->bndry = get_el_bndry_vec(bfcts);
723 md->coeff = get_el_real_vec(bfcts);
724 md->coeff_d = get_el_real_d_vec(bfcts);
725 md->int_vec = get_el_int_vec(bfcts);
726 md->real_vec = get_el_real_vec(bfcts);
727 md->real_d_vec = get_el_real_d_vec(bfcts);
728 md->uchar_vec = get_el_uchar_vec(bfcts);
729 md->schar_vec = get_el_schar_vec(bfcts);
730 md->ptr_vec = get_el_ptr_vec(bfcts);
731 }
732 #endif
733
734 #define BAS_FCTS_INIT_NAME "bas_fcts_init"
735
736 typedef struct bas_fcts_plugin BFCTS_PLUGIN;
737 struct bas_fcts_plugin
738 {
739 BAS_FCTS_INIT_FCT bas_fcts_init;
740 BFCTS_PLUGIN *next;
741 };
742
743 static BFCTS_PLUGIN *plugin_list;
744
add_bas_fcts_plugin(BAS_FCTS_INIT_FCT init_fct)745 void add_bas_fcts_plugin(BAS_FCTS_INIT_FCT init_fct)
746 {
747 BFCTS_PLUGIN *plugin = MEM_ALLOC(1, BFCTS_PLUGIN);
748
749 plugin->bas_fcts_init = init_fct;
750 plugin->next = plugin_list;
751 plugin_list = plugin;
752 }
753
plugin_from_module(const char * dl_bas_fcts)754 static void plugin_from_module(const char *dl_bas_fcts)
755 {
756 FUNCNAME("plugin_from_module");
757 #if HAVE_LIBLTDL && HAVE_LTDL_H
758 static bool init_done;
759 lt_dlhandle dl_handle;
760 BAS_FCTS_INIT_FCT init_fct;
761
762 if (!init_done) {
763 init_done = true;
764
765 TEST_EXIT(lt_dlinit() == 0,
766 "Could not initialize libltdl (%s).\n", lt_dlerror);
767 }
768 if (dl_bas_fcts != NULL) {
769 TEST_EXIT((dl_handle = lt_dlopenext(dl_bas_fcts)) != NULL,
770 "Could not dlopen \"%s\" (%s)\n", dl_bas_fcts, lt_dlerror());
771 TEST_EXIT((init_fct = (BAS_FCTS_INIT_FCT)(unsigned long)
772 lt_dlsym(dl_handle, BAS_FCTS_INIT_NAME)) != NULL,
773 "Could not resolve \"%s\" (%s)\n",
774 BAS_FCTS_INIT_NAME, lt_dlerror());
775 (void)lt_dlmakeresident(dl_handle);
776 add_bas_fcts_plugin(init_fct);
777 } else {
778 if ((dl_handle = lt_dlopenext(dl_bas_fcts)) == NULL) {
779 return;
780 }
781 if ((init_fct = (BAS_FCTS_INIT_FCT)(unsigned long)
782 lt_dlsym(dl_handle, BAS_FCTS_INIT_NAME)) == NULL) {
783 return;
784 }
785 (void)lt_dlmakeresident(dl_handle);
786 add_bas_fcts_plugin(init_fct);
787 }
788 #else
789 WARNING("No dynamic loading facilities available.\n");
790 #endif
791 }
792
793 /* Once-only initialization. We load any environment specified plugins
794 * first, then we load the linked-in plugin (if defined).
795 *
796 * We search for plugins specified by the ALBERTA_BAS_FCTS_LIB_XD
797 * environment variable, where X stands for DIM_OF_WORLD.
798 */
plugin_init(void)799 static void plugin_init(void)
800 {
801 static bool done;
802 char *env_plugin;
803
804 if (done) {
805 return;
806 }
807
808 done = true;
809
810 #if ALBERTA_DEBUG
811 env_plugin =
812 getenv("ALBERTA_BAS_FCTS_LIB_"CPP_STRINGIZE(DIM_OF_WORLD)"D_DEBUG");
813 MSG("Trying to load \"%s\"\n",
814 "ALBERTA_BAS_FCTS_LIB_"CPP_STRINGIZE(DIM_OF_WORLD)"D_DEBUG");
815 #elif ALBERTA_PROFILE
816 env_plugin =
817 getenv("ALBERTA_BAS_FCTS_LIB_"CPP_STRINGIZE(DIM_OF_WORLD)"D_PROFILE");
818 MSG("Trying to load \"%s\"\n",
819 "ALBERTA_BAS_FCTS_LIB_"CPP_STRINGIZE(DIM_OF_WORLD)"D_PROFILE");
820 #else
821 env_plugin = getenv("ALBERTA_BAS_FCTS_LIB_"CPP_STRINGIZE(DIM_OF_WORLD)"D");
822 MSG("Trying to load \"%s\"\n",
823 "ALBERTA_BAS_FCTS_LIB_"CPP_STRINGIZE(DIM_OF_WORLD)"D");
824 #endif
825 if (env_plugin) {
826 plugin_from_module(env_plugin);
827 }
828 plugin_from_module(NULL); /* Search for linked-in module. */
829 }
830
get_bas_fcts(int dim,const char * name)831 const BAS_FCTS *get_bas_fcts(int dim, const char *name)
832 {
833 FUNCNAME("get_bas_fcts");
834 static bool initialized;
835 const BAS_FCTS *bas_fcts;
836 BFCTS_NODE *bfcts_node;
837 BFCTS_PLUGIN *plug;
838 size_t len;
839
840 if (!initialized) {
841 int dim, deg;
842
843 for (dim = 0; dim <= DIM_MAX; dim++) {
844 for (deg = 0; deg <= 4; deg++) {
845 get_lagrange(dim, deg);
846 }
847 for (deg = 0; deg <= 2; deg++) {
848 get_discontinuous_lagrange(dim, deg);
849 }
850 for (deg = 1; deg <= 2; deg++) {
851 get_disc_ortho_poly(dim, deg);
852 }
853 }
854 initialized = true;
855 }
856
857 if (!name) {
858 ERROR("no name specified; cannot return pointer to basis functions\n");
859 return NULL;
860 }
861 if ((len = strlen(name)) == 0) {
862 ERROR("empty name; cannot return pointer to basis functions\n");
863 return NULL;
864 }
865 /* discard the "_Xd" suffix */
866 if (name[len-3] == '_' && name[len-2] == '0'+dim && name[len-1] == 'd') {
867 len -= 3;
868 }
869
870 /* Special compatibility hack: in the course of proceeding to
871 * ALBERTA-2 the p.w. constant Lagrange basis functions were renamed
872 * to "disc_lagrange0".
873 */
874 if (strncmp(name, "lagrange0", len) == 0) {
875 name = "disc_lagrange0";
876 len += 5;
877 }
878
879 for (bfcts_node = bas_fcts_list[dim];
880 bfcts_node; bfcts_node = bfcts_node->next) {
881 if (len == bfcts_node->name_len &&
882 !strncmp(bfcts_node->bas_fcts->name, name, len)) {
883 return bfcts_node->bas_fcts;
884 }
885 }
886
887 /* Not found: traverse the list of plugin-modules. */
888 plugin_init();
889
890 for (plug = plugin_list; plug != NULL; plug = plug->next) {
891 if ((bas_fcts = plug->bas_fcts_init(dim, DIM_OF_WORLD, name)) != NULL) {
892 new_bas_fcts(bas_fcts);
893 return bas_fcts;
894 }
895 }
896
897 ERROR("basis functions with name %s not found in list of all functions\n", name);
898 return NULL;
899 }
900
901 /* Clone the set of basis-functions specified as "head", if tail !=
902 * NULL, then add the copy of "head" as head to the list specified by
903 * tail.
904 *
905 * Example:
906 *
907 * lagrange = get_lagrange(dim, degree);
908 * bubble = get_bas_fcts(dim, "Bubble");
909 * bas_fcts = chain_bas_fcts(lagrange, chain_bas_fcts(bubble, NULL));
910 * old_fcts = new_bas_fcts(bas_fcts);
911 *
912 * The chained basis functions will be added to the list of known
913 * basis functions using the name
914 *
915 * "OWN_NAME#NEXT_NAME#NEXT_NEXT_NAME"
916 *
917 * During the construction of the name any "_Xd" suffixes are
918 * discarded in order not to make the name too complicated. At the end
919 * of the name-chain an appropriate "_Xd" suffix is added. The _Xd
920 * suffixes are only for debugging, they are ignore everywhere else.
921 *
922 * A basis function chain is cyclic. The trace-spaces of the given
923 * sets of basis functions are chained together accordingly.
924 *
925 * NOTE: this function does NOT call new_bas_fcts(); the caller has to
926 * do so after constructing the desired chain.
927 */
928 static
929 INIT_EL_TAG bfcts_chain_init_element(const EL_INFO *el_info, void *thisptr);
930
chain_bas_fcts(const BAS_FCTS * head,BAS_FCTS * tail)931 BAS_FCTS *chain_bas_fcts(const BAS_FCTS *head, BAS_FCTS *tail)
932 {
933 FUNCNAME("chain_bas_fcts");
934 BAS_FCTS *bfcts, *pos;
935 FLAGS fill_flags = head->fill_flags;
936 char *name;
937 size_t head_len = 0, tail_len = 0, len;
938 bool init_element_needed = false;
939 int dim = head->dim;
940
941 bfcts = MEM_ALLOC(1, BAS_FCTS);
942
943 *bfcts = *head;
944 DBL_LIST_INIT(&bfcts->chain);
945 bfcts->unchained = head; /* remember the unchained instance */
946
947 if ((tail && tail->init_element) || bfcts->init_element) {
948 init_element_needed = true;
949 if (tail != NULL) {
950 fill_flags |= tail->fill_flags;
951 }
952 /* we need to install a proxy-initializer */
953 INIT_ELEMENT_DEFUN(bfcts, bfcts_chain_init_element, fill_flags);
954 }
955
956 head_len = strlen(head->name);
957 if (head->name[head_len-3] == '_' &&
958 head->name[head_len-2] == '0'+dim &&
959 head->name[head_len-1] == 'd') {
960 head_len -= 3;
961 }
962 if (tail) {
963 TEST_EXIT(dim == tail->dim,
964 "Trying to chain basis function with different dimensions.\n");
965 tail_len = strlen(tail->name);
966 if (tail->name[tail_len-3] == '_' &&
967 tail->name[tail_len-2] == '0'+dim &&
968 tail->name[tail_len-1] == 'd') {
969 tail_len -= 3;
970 }
971 }
972
973 len = head_len + (tail ? strlen("#") + tail_len : 0) + sizeof("_Xd");
974
975 /* Use plain malloc() because of strdup() stuff */
976 bfcts->name = name = (char *)malloc(len);
977 sprintf(name, "%.*s%s%.*s_%dd",
978 (int)head_len, head->name,
979 tail ? "#" : "", (int)tail_len, tail ? tail->name : "",
980 dim);
981
982 if (dim > 0) {
983 TEST_EXIT(head->trace_bas_fcts != NULL &&
984 (tail == NULL || tail->trace_bas_fcts != NULL),
985 "Missing trace basis functions.\n");
986 bfcts->trace_bas_fcts =
987 chain_bas_fcts(head->trace_bas_fcts,
988 tail ? (BAS_FCTS *)tail->trace_bas_fcts : NULL);
989 }
990
991 if (tail) {
992 /* Adjust names in cyclic list and add a copy of head to tail */
993 CHAIN_ADD_TAIL(tail, bfcts);
994 bfcts->degree = MAX(bfcts->degree, tail->degree);
995 CHAIN_FOREACH(pos, bfcts, BAS_FCTS) {
996 name = (char *)malloc(len);
997 sprintf(name, "%.*s#%.*s_%dd",
998 (int)tail_len, pos->name, (int)head_len, head->name, dim);
999 free((char *)pos->name);
1000 pos->name = name;
1001 if (init_element_needed) {
1002 /* we need to install a proxy-initializer */
1003 INIT_ELEMENT_DEFUN(pos, bfcts_chain_init_element, fill_flags);
1004 }
1005 }
1006 }
1007
1008 return bfcts;
1009 }
1010
1011 /* Chain initializer, walk the entire chain and initialize each
1012 * element in turn. This should be (algorithmically) the same as
1013 * qfast_chain_init_element().
1014 */
1015 static INIT_EL_TAG
bfcts_chain_init_element(const EL_INFO * el_info,void * thisptr)1016 bfcts_chain_init_element(const EL_INFO *el_info, void *thisptr)
1017 {
1018 BAS_FCTS *bas_fcts = (BAS_FCTS *)thisptr;
1019 INIT_EL_TAG chain_tag = INIT_EL_TAG_NONE;
1020 INIT_EL_TAG old_tag;
1021 bool new_tag = false;
1022
1023 CHAIN_DO(bas_fcts, BAS_FCTS) {
1024 if (INIT_ELEMENT_NEEDED(bas_fcts->unchained)) {
1025 old_tag = INIT_EL_TAG_CTX_TAG(&bas_fcts->tag_ctx);
1026 chain_tag |= INIT_ELEMENT_SINGLE(el_info, bas_fcts);
1027 if (old_tag != INIT_EL_TAG_CTX_TAG(&bas_fcts->tag_ctx)) {
1028 new_tag = true;
1029 }
1030 } else {
1031 chain_tag |= INIT_EL_TAG_DFLT;
1032 }
1033 } CHAIN_WHILE(bas_fcts, BAS_FCTS);
1034
1035 if (chain_tag == INIT_EL_TAG_NONE) {
1036 chain_tag = INIT_EL_TAG_DFLT;
1037 }
1038 if (chain_tag == INIT_EL_TAG_DFLT || chain_tag == INIT_EL_TAG_NULL) {
1039 return chain_tag;
1040 } else {
1041 if (new_tag) {
1042 INIT_EL_TAG_CTX_UNIQ(&bas_fcts->tag_ctx);
1043 }
1044 return INIT_EL_TAG_CTX_TAG(&bas_fcts->tag_ctx);
1045 }
1046 }
1047
1048