1 /**@file
2 *
3 * Face-Bubble basis functions, meant to be chained with some standard
4 * Lagrange space, optionally with divergence preserving interpolation
5 * routine. This variant is a fancy space with DOFs defined on an
6 * admin belonging to a trace mesh. We assume that the mapping between
7 * bulk-elements and trace-elements is 1-1, i.e. at most one face of a
8 * bulk element belongs to the trace mesh.
9 *
10 * author: Claus-Justus Heine
11 * Abteilung fuer Angewandte Mathematik
12 * Albert-Ludwigs-Universitaet Freiburg
13 * Hermann-Herder-Str. 10
14 * 79104 Freiburg i.Br.
15 * Germany
16 * claus@mathematik.uni-freiburg.de
17 *
18 * Copyright (C) by Claus-Justus Heine (2008-2009).
19 *
20 */
21
22 #if HAVE_CONFIG_H
23 # include "config.h"
24 #endif
25
26 #include <alberta.h>
27
28 #include "albas.h"
29
30 /* The weights are chosen s.t. the integral over the faces of the
31 * reference simplex is 1.
32 */
33 #define WEIGHT_1D 1.0
34 #define WEIGHT_2D 6.0
35 #define WEIGHT_3D 120.0
36
37 typedef struct btb_data
38 {
39 const EL_INFO *cur_el_info;
40 const EL *cur_el;
41 const EL *cur_trace_el[N_WALLS_MAX];
42 int cur_wall[N_WALLS_MAX];
43 MESH *trace_mesh;
44 int trace_id;
45 REAL_D wall_normal[N_WALLS_MAX];
46 BAS_FCT phi[N_WALLS_MAX];
47 GRD_BAS_FCT grd_phi[N_WALLS_MAX];
48 D2_BAS_FCT D2_phi[N_WALLS_MAX];
49 BAS_FCT_D phi_d[N_WALLS_MAX];
50 int trace_dof_map[N_WALLS_MAX];
51 const WALL_QUAD *wquad;
52 const WALL_QUAD_FAST *wqfast;
53 int inter_deg;
54 } BTB_DATA;
55
56 /* scalar factors */
57
58 #if DIM_OF_WORLD >= 1
59 /* <<< 1d */
60
btb_phi_w0_1d(const REAL_B lambda,const BAS_FCTS * bfcts)61 static REAL btb_phi_w0_1d(const REAL_B lambda, const BAS_FCTS *bfcts)
62 {
63 return WEIGHT_1D*lambda[1];
64 }
65
btb_phi_w1_1d(const REAL_B lambda,const BAS_FCTS * bfcts)66 static REAL btb_phi_w1_1d(const REAL_B lambda, const BAS_FCTS *bfcts)
67 {
68 return WEIGHT_1D*lambda[0];
69 }
70
btb_grd_phi_w0_1d(const REAL_B lambda,const BAS_FCTS * bfcts)71 static const REAL *btb_grd_phi_w0_1d(const REAL_B lambda, const BAS_FCTS *bfcts)
72 {
73 static const REAL_B grd = { 0.0, WEIGHT_1D, };
74 return grd;
75 }
76
btb_grd_phi_w1_1d(const REAL_B lambda,const BAS_FCTS * bfcts)77 static const REAL *btb_grd_phi_w1_1d(const REAL_B lambda, const BAS_FCTS *bfcts)
78 {
79 static const REAL_B grd = { WEIGHT_1D, };
80 return grd;
81 }
82
btb_D2_phi_w0_1d(const REAL_B lambda,const BAS_FCTS * bfcts)83 static const REAL_B *btb_D2_phi_w0_1d(const REAL_B lambda, const BAS_FCTS *bfcts)
84 {
85 static const REAL_BB D2;
86 return D2;
87 }
88
btb_D2_phi_w1_1d(const REAL_B lambda,const BAS_FCTS * bfcts)89 static const REAL_B *btb_D2_phi_w1_1d(const REAL_B lambda, const BAS_FCTS *bfcts)
90 {
91 static const REAL_BB D2;
92 return D2;
93 }
94
95 /* >>> */
96 #endif
97
98 #if DIM_OF_WORLD >= 2
99 /* <<< 2d */
100
btb_phi_w0_2d(const REAL_B lambda,const BAS_FCTS * bfcts)101 static REAL btb_phi_w0_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
102 {
103 return WEIGHT_2D*lambda[1]*lambda[2];
104 }
105
btb_phi_w1_2d(const REAL_B lambda,const BAS_FCTS * bfcts)106 static REAL btb_phi_w1_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
107 {
108 return WEIGHT_2D*lambda[0]*lambda[2];
109 }
110
btb_phi_w2_2d(const REAL_B lambda,const BAS_FCTS * bfcts)111 static REAL btb_phi_w2_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
112 {
113 return WEIGHT_2D*lambda[0]*lambda[1];
114 }
115
btb_grd_phi_w0_2d(const REAL_B lambda,const BAS_FCTS * bfcts)116 static const REAL *btb_grd_phi_w0_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
117 {
118 static REAL_B grd;
119
120 grd[1] = WEIGHT_2D*lambda[2];
121 grd[2] = WEIGHT_2D*lambda[1];
122
123 return grd;
124 }
125
btb_grd_phi_w1_2d(const REAL_B lambda,const BAS_FCTS * bfcts)126 static const REAL *btb_grd_phi_w1_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
127 {
128 static REAL_B grd;
129
130 grd[0] = WEIGHT_2D*lambda[2];
131 grd[2] = WEIGHT_2D*lambda[0];
132
133 return grd;
134 }
135
btb_grd_phi_w2_2d(const REAL_B lambda,const BAS_FCTS * bfcts)136 static const REAL *btb_grd_phi_w2_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
137 {
138 static REAL_B grd;
139
140 grd[0] = WEIGHT_2D*lambda[1];
141 grd[1] = WEIGHT_2D*lambda[0];
142
143 return grd;
144 }
145
btb_D2_phi_w0_2d(const REAL_B lambda,const BAS_FCTS * bfcts)146 static const REAL_B *btb_D2_phi_w0_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
147 {
148 static const REAL_BB D2 = {
149 { 0.0, },
150 { 0.0, 0.0, WEIGHT_2D, },
151 { 0.0, WEIGHT_2D, 0.0, },
152 };
153
154 return D2;
155 }
156
btb_D2_phi_w1_2d(const REAL_B lambda,const BAS_FCTS * bfcts)157 static const REAL_B *btb_D2_phi_w1_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
158 {
159 static const REAL_BB D2 = {
160 { 0.0, 0.0, WEIGHT_2D, },
161 { 0.0, },
162 { WEIGHT_2D, },
163 };
164
165 return D2;
166 }
167
btb_D2_phi_w2_2d(const REAL_B lambda,const BAS_FCTS * bfcts)168 static const REAL_B *btb_D2_phi_w2_2d(const REAL_B lambda, const BAS_FCTS *bfcts)
169 {
170 static const REAL_BB D2 = {
171 { 0.0, WEIGHT_2D, },
172 { WEIGHT_2D, },
173 };
174
175 return D2;
176 }
177
178 /* >>> */
179 #endif
180
181 #if DIM_OF_WORLD >= 3
182 /* <<< 3d */
183
btb_phi_w0_3d(const REAL_B lambda,const BAS_FCTS * bfcts)184 static REAL btb_phi_w0_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
185 {
186 return WEIGHT_3D*lambda[1]*lambda[2]*lambda[3];
187 }
188
btb_phi_w1_3d(const REAL_B lambda,const BAS_FCTS * bfcts)189 static REAL btb_phi_w1_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
190 {
191 return WEIGHT_3D*lambda[0]*lambda[2]*lambda[3];
192 }
193
btb_phi_w2_3d(const REAL_B lambda,const BAS_FCTS * bfcts)194 static REAL btb_phi_w2_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
195 {
196 return WEIGHT_3D*lambda[0]*lambda[1]*lambda[3];
197 }
198
btb_phi_w3_3d(const REAL_B lambda,const BAS_FCTS * bfcts)199 static REAL btb_phi_w3_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
200 {
201 return WEIGHT_3D*lambda[0]*lambda[1]*lambda[2];
202 }
203
btb_grd_phi_w0_3d(const REAL_B lambda,const BAS_FCTS * bfcts)204 static const REAL *btb_grd_phi_w0_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
205 {
206 static REAL_B grd;
207
208 grd[1] = WEIGHT_3D*lambda[2]*lambda[3];
209 grd[2] = WEIGHT_3D*lambda[1]*lambda[3];
210 grd[3] = WEIGHT_3D*lambda[1]*lambda[2];
211
212 return grd;
213 }
214
btb_grd_phi_w1_3d(const REAL_B lambda,const BAS_FCTS * bfcts)215 static const REAL *btb_grd_phi_w1_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
216 {
217 static REAL_B grd;
218
219 grd[0] = WEIGHT_3D*lambda[2]*lambda[3];
220 grd[2] = WEIGHT_3D*lambda[0]*lambda[3];
221 grd[3] = WEIGHT_3D*lambda[0]*lambda[2];
222
223 return grd;
224 }
225
btb_grd_phi_w2_3d(const REAL_B lambda,const BAS_FCTS * bfcts)226 static const REAL *btb_grd_phi_w2_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
227 {
228 static REAL_B grd;
229
230 grd[0] = WEIGHT_3D*lambda[1]*lambda[3];
231 grd[1] = WEIGHT_3D*lambda[0]*lambda[3];
232 grd[3] = WEIGHT_3D*lambda[0]*lambda[1];
233
234 return grd;
235 }
236
btb_grd_phi_w3_3d(const REAL_B lambda,const BAS_FCTS * bfcts)237 static const REAL *btb_grd_phi_w3_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
238 {
239 static REAL_B grd;
240
241 grd[0] = WEIGHT_3D*lambda[1]*lambda[2];
242 grd[1] = WEIGHT_3D*lambda[0]*lambda[2];
243 grd[2] = WEIGHT_3D*lambda[0]*lambda[1];
244
245 return grd;
246 }
247
btb_D2_phi_w0_3d(const REAL_B lambda,const BAS_FCTS * bfcts)248 static const REAL_B *btb_D2_phi_w0_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
249 {
250 static REAL_BB D2;
251
252 D2[2][1] = D2[1][2] = WEIGHT_3D*lambda[3];
253 D2[3][1] = D2[1][3] = WEIGHT_3D*lambda[2];
254 D2[2][3] = D2[3][2] = WEIGHT_3D*lambda[1];
255
256 return (const REAL_B *)D2;
257 }
258
btb_D2_phi_w1_3d(const REAL_B lambda,const BAS_FCTS * bfcts)259 static const REAL_B *btb_D2_phi_w1_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
260 {
261 static REAL_BB D2;
262
263 D2[2][0] = D2[0][2] = WEIGHT_3D*lambda[3];
264 D2[3][0] = D2[0][3] = WEIGHT_3D*lambda[2];
265 D2[2][3] = D2[3][2] = WEIGHT_3D*lambda[0];
266
267 return (const REAL_B *)D2;
268 }
269
btb_D2_phi_w2_3d(const REAL_B lambda,const BAS_FCTS * bfcts)270 static const REAL_B *btb_D2_phi_w2_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
271 {
272 static REAL_BB D2;
273
274 D2[1][0] = D2[0][1] = WEIGHT_3D*lambda[3];
275 D2[3][0] = D2[0][3] = WEIGHT_3D*lambda[1];
276 D2[1][3] = D2[3][1] = WEIGHT_3D*lambda[0];
277
278 return (const REAL_B *)D2;
279 }
280
btb_D2_phi_w3_3d(const REAL_B lambda,const BAS_FCTS * bfcts)281 static const REAL_B *btb_D2_phi_w3_3d(const REAL_B lambda, const BAS_FCTS *bfcts)
282 {
283 static REAL_BB D2;
284
285 D2[1][0] = D2[0][1] = WEIGHT_3D*lambda[2];
286 D2[2][0] = D2[0][2] = WEIGHT_3D*lambda[1];
287 D2[1][2] = D2[2][1] = WEIGHT_3D*lambda[0];
288
289 return (const REAL_B *)D2;
290 }
291
292 /* >>> */
293 #endif
294
295 static const REAL *
btb_phi_d_0(const REAL_B lambda,const BAS_FCTS * bfcts)296 btb_phi_d_0(const REAL_B lambda, const BAS_FCTS *bfcts)
297 {
298 BTB_DATA *data = (BTB_DATA *)bfcts->ext_data;
299
300 return data->wall_normal[0];
301 }
302
303 static const REAL *
btb_phi_d_1(const REAL_B lambda,const BAS_FCTS * bfcts)304 btb_phi_d_1(const REAL_B lambda, const BAS_FCTS *bfcts)
305 {
306 BTB_DATA *data = (BTB_DATA *)bfcts->ext_data;
307
308 return data->wall_normal[1];
309 }
310
311 #if DIM_MAX >= 2
312 static const REAL *
btb_phi_d_2(const REAL_B lambda,const BAS_FCTS * bfcts)313 btb_phi_d_2(const REAL_B lambda, const BAS_FCTS *bfcts)
314 {
315 BTB_DATA *data = (BTB_DATA *)bfcts->ext_data;
316
317 return data->wall_normal[2];
318 }
319 #endif
320
321 #if DIM_MAX >= 3
322 static const REAL *
btb_phi_d_3(const REAL_B lambda,const BAS_FCTS * bfcts)323 btb_phi_d_3(const REAL_B lambda, const BAS_FCTS *bfcts)
324 {
325 BTB_DATA *data = (BTB_DATA *)bfcts->ext_data;
326
327 return data->wall_normal[3];
328 }
329 #endif
330
331 static const BAS_FCT btb_phi[DIM_MAX+1][N_WALLS_MAX] = {
332 { NULL, /* 0 dim */ },
333 #if DIM_MAX >= 1
334 { btb_phi_w0_1d, btb_phi_w1_1d, },
335 #endif
336 #if DIM_MAX >= 2
337 { btb_phi_w0_2d, btb_phi_w1_2d, btb_phi_w2_2d, },
338 #endif
339 #if DIM_MAX >= 3
340 { btb_phi_w0_3d, btb_phi_w1_3d, btb_phi_w2_3d, btb_phi_w3_3d, },
341 #endif
342 };
343
344 static const GRD_BAS_FCT btb_grd_phi[DIM_MAX+1][N_WALLS_MAX] = {
345 { NULL, /* 0 dim */ },
346 #if DIM_MAX >= 1
347 { btb_grd_phi_w0_1d, btb_grd_phi_w1_1d, },
348 #endif
349 #if DIM_MAX >= 2
350 { btb_grd_phi_w0_2d, btb_grd_phi_w1_2d, btb_grd_phi_w2_2d, },
351 #endif
352 #if DIM_MAX >= 3
353 { btb_grd_phi_w0_3d, btb_grd_phi_w1_3d, btb_grd_phi_w2_3d, btb_grd_phi_w3_3d, },
354 #endif
355 };
356
357 static const D2_BAS_FCT btb_D2_phi[DIM_MAX+1][N_WALLS_MAX] = {
358 { NULL, /* 0 dim */ },
359 #if DIM_MAX >= 1
360 { btb_D2_phi_w0_1d, btb_D2_phi_w1_1d, },
361 #endif
362 #if DIM_MAX >= 2
363 { btb_D2_phi_w0_2d, btb_D2_phi_w1_2d, btb_D2_phi_w2_2d, },
364 #endif
365 #if DIM_MAX >= 3
366 { btb_D2_phi_w0_3d, btb_D2_phi_w1_3d, btb_D2_phi_w2_3d, btb_D2_phi_w3_3d, },
367 #endif
368 };
369
370 static const BAS_FCT_D btb_phi_d[] = {
371 btb_phi_d_0, btb_phi_d_1,
372 #if DIM_MAX >= 2
373 btb_phi_d_2,
374 # if DIM_MAX >= 3
375 btb_phi_d_3
376 # endif
377 #endif
378 };
379
380 /* init_element() routines for vector-valued basis-functions: ALBERTA
381 * assumes that the directional part is never constant, so the return
382 * value should _ONLY_ reflect the status of the scalar factor.
383 *
384 * This version need FILL_NEIGH all the time. Maybe we should store
385 * the wall normals globally, a suitable
386 * refine_interpol()/coarse_restrict() routine could then be used to
387 * update the global normals during mesh adaption. The normals would
388 * be automatically oriented.
389 */
390 static
btb_init_element(const EL_INFO * el_info,void * vself)391 INIT_EL_TAG btb_init_element(const EL_INFO *el_info, void *vself)
392 {
393 FUNCNAME("btb_init_element");
394 BAS_FCTS *self = (BAS_FCTS *)vself;
395 BTB_DATA *data = (BTB_DATA *)self->ext_data;
396 MESH *mesh;
397 int w, dim, n_bas;
398 const EL_GEOM_CACHE *elgc;
399 static bool coord_warning = false;
400
401 if (el_info == NULL) {
402 data->cur_el = NULL;
403 data->cur_el_info = NULL;
404
405 self->dir_pw_const = true;
406 self->n_bas_fcts = 0;
407 for (w = 0; w < N_WALLS_MAX; w++) {
408 self->n_trace_bas_fcts[w] =
409 ((BAS_FCTS *)self->unchained)->n_trace_bas_fcts[w] = 0;
410 }
411 memset(data->cur_wall, -1, sizeof(data->cur_wall));
412 memset(data->cur_trace_el, 0, sizeof(data->cur_trace_el));
413 self->n_bas_fcts =
414 ((BAS_FCTS *)self->unchained)->n_bas_fcts = 0;
415
416 INIT_EL_TAG_CTX_DFLT(&self->tag_ctx);
417 return INIT_EL_TAG_CTX_TAG(&self->tag_ctx);
418 }
419
420 if (data->cur_el == el_info->el && data->cur_el_info == el_info) {
421 return INIT_EL_TAG_CTX_TAG(&self->tag_ctx);
422 }
423 data->cur_el = el_info->el;
424 data->cur_el_info = el_info;
425
426 mesh = el_info->mesh;
427 dim = mesh->dim;
428
429 if (mesh->parametric) {
430 ERROR_EXIT("Not yet implemented for parametric meshes.\n");
431 }
432
433 if (data->trace_mesh == NULL) {
434 data->trace_mesh = lookup_submesh_by_id(mesh, data->trace_id);
435 TEST_EXIT(data->trace_mesh != NULL,
436 "No trace-mesh with id %d\n", data->trace_id);
437 }
438
439 if ((el_info->fill_flag & FILL_COORDS) == 0) {
440 if (!coord_warning) {
441 WARNING("FILL_COORDS not set, doing nothing.\n");
442 coord_warning = true;
443 }
444 return INIT_EL_TAG_CTX_TAG(&self->tag_ctx);
445 }
446
447 /* Determine whether one of the faces belongs to our dedicated
448 * trace-mesh.
449 */
450 for (n_bas = w = 0; w < N_WALLS(dim); w++) {
451 const EL *tr_el;
452
453 if ((tr_el = get_slave_el(el_info->el, w, data->trace_mesh)) == NULL) {
454 self->n_trace_bas_fcts[w] =
455 ((BAS_FCTS *)self->unchained)->n_trace_bas_fcts[w] = 0;
456 data->cur_wall[n_bas] = -1;
457 continue;
458 }
459 data->cur_trace_el[n_bas] = tr_el;
460 data->cur_wall[n_bas] = w;
461 elgc = fill_el_geom_cache(el_info, FILL_EL_WALL_NORMAL(w));
462 COPY_DOW(elgc->wall_normal[w], data->wall_normal[n_bas]);
463
464 data->phi[n_bas] = btb_phi[dim][w];
465 data->grd_phi[n_bas] = btb_grd_phi[dim][w];
466 data->D2_phi[n_bas] = btb_D2_phi[dim][w];
467 data->phi_d[n_bas] = btb_phi_d[n_bas];
468
469 self->n_trace_bas_fcts[w] =
470 ((BAS_FCTS *)self->unchained)->n_trace_bas_fcts[w] = 1;
471
472 data->trace_dof_map[w] = n_bas;
473
474 ++n_bas;
475 }
476
477 if (n_bas == 0) {
478 /* The default (not NULL) case, actually */
479 if (INIT_EL_TAG_CTX_TAG(&self->tag_ctx) != INIT_EL_TAG_DFLT) {
480 for (w = 0; w < N_WALLS(dim); w++) {
481 self->n_trace_bas_fcts[w] =
482 ((BAS_FCTS *)self->unchained)->n_trace_bas_fcts[w] = 0;
483 }
484 memset(data->cur_wall, -1, sizeof(data->cur_wall));
485 memset(data->cur_trace_el, 0, sizeof(data->cur_trace_el));
486 self->n_bas_fcts =
487 ((BAS_FCTS *)self->unchained)->n_bas_fcts = 0;
488 }
489 INIT_EL_TAG_CTX_DFLT(&self->tag_ctx);
490 } else {
491 self->n_bas_fcts =
492 ((BAS_FCTS *)self->unchained)->n_bas_fcts = n_bas;
493
494 INIT_EL_TAG_CTX_UNIQ(&self->tag_ctx);
495 }
496
497 return INIT_EL_TAG_CTX_TAG(&self->tag_ctx);
498 }
499
500 #undef DEF_EL_VEC_BTB
501 #define DEF_EL_VEC_BTB(type, name) \
502 DEF_EL_VEC_CONST(type, name, 1, N_WALLS_MAX)
503
504 #undef DEFUN_GET_EL_VEC_BTB
505 #define DEFUN_GET_EL_VEC_BTB(dim, name, type, admin, self, body, ...) \
506 static const EL_##type##_VEC * \
507 btb_get_##name(type##_VEC_TYPE *vec, const EL *el, __VA_ARGS__) \
508 { \
509 FUNCNAME("bulk_trace_bubble_get_"#name); \
510 static DEF_EL_VEC_BTB(type, rvec_space); \
511 type##_VEC_TYPE *rvec = vec ? vec : rvec_space->vec; \
512 BTB_DATA *data = (BTB_DATA *)self->ext_data; \
513 int n0, node, i; \
514 \
515 DEBUG_TEST_EXIT(true, ""); \
516 \
517 node = (admin)->mesh->node[CENTER]; \
518 n0 = (admin)->n0_dof[CENTER]; \
519 for (i = 0; i < self->n_bas_fcts; i++) { \
520 DOF **dofptr = data->cur_trace_el[i]->dof, dof; \
521 \
522 dof = dofptr[node][n0]; \
523 body; \
524 } \
525 \
526 if (vec) { \
527 return NULL; \
528 } else { \
529 rvec_space->n_components = self->n_bas_fcts; \
530 return rvec_space; \
531 } \
532 } \
533 struct _AI_semicolon_dummy
534
535 #undef DEFUN_GET_EL_DOF_VEC_BTB
536 #define DEFUN_GET_EL_DOF_VEC_BTB(name, type, ASSIGN) \
537 DEFUN_GET_EL_VEC_BTB(dv->fe_space->admin->mesh->dim, \
538 _##name##_vec, type, dv->fe_space->admin, \
539 dv->fe_space->bas_fcts, \
540 ASSIGN(dv->vec[dof], rvec[i]), \
541 const DOF_##type##_VEC *dv); \
542 static const EL_##type##_VEC * \
543 btb_get_##name##_vec(type##_VEC_TYPE *vec, const EL *el, \
544 const DOF_##type##_VEC *dv) \
545 { \
546 EL_##type##_VEC *vec_loc = dv->vec_loc; \
547 \
548 if (vec != NULL || vec_loc == NULL) { \
549 return btb_get__##name##_vec(vec, el, dv); \
550 } else { \
551 btb_get__##name##_vec(vec_loc->vec, el, dv); \
552 return vec_loc; \
553 } \
554 } \
555 struct _AI_semicolon_dummy
556
557
558 #undef COPY_EQ
559 #define COPY_EQ(a, b) (b) = (a)
560
561 /*******************************************************************************
562 * functions for combining basisfunctions with coefficients
563 ******************************************************************************/
564
565 DEFUN_GET_EL_VEC_BTB(self->dim, dof_indices, DOF, admin, self,
566 rvec[i] = dof,
567 const DOF_ADMIN *admin, const BAS_FCTS *self);
568
569 static const EL_BNDRY_VEC *
btb_get_bound_1d(BNDRY_FLAGS * vec,const EL_INFO * el_info,const BAS_FCTS * self)570 btb_get_bound_1d(BNDRY_FLAGS *vec,
571 const EL_INFO *el_info, const BAS_FCTS *self)
572 {
573 static DEF_EL_VEC_BTB(BNDRY, rvec_space);
574 BTB_DATA *data = (BTB_DATA *)self->ext_data;
575 BNDRY_FLAGS *rvec = vec ? vec : rvec_space->vec;
576 int i;
577
578 for (i = 0; i < self->n_bas_fcts; i++) {
579 int wall = data->cur_wall[i];
580 BNDRY_FLAGS_CPY(rvec[i], el_info->vertex_bound[1-wall]);
581 }
582
583 return vec ? NULL : rvec_space;
584 }
585
586 #if DIM_OF_WORLD > 1
587 static const EL_BNDRY_VEC *
btb_get_bound_2d(BNDRY_FLAGS * vec,const EL_INFO * el_info,const BAS_FCTS * self)588 btb_get_bound_2d(BNDRY_FLAGS *vec,
589 const EL_INFO *el_info, const BAS_FCTS *self)
590 {
591 static DEF_EL_VEC_BTB(BNDRY, rvec_space);
592 BTB_DATA *data = (BTB_DATA *)self->ext_data;
593 BNDRY_FLAGS *rvec = vec ? vec : rvec_space->vec;
594 int i;
595
596 for (i = 0; i < self->n_bas_fcts; i++) {
597 int wall = data->cur_wall[i];
598 BNDRY_FLAGS_CPY(rvec[i], el_info->edge_bound[wall]);
599 }
600
601 return vec ? NULL : rvec_space;
602 }
603 #endif
604
605 #if DIM_OF_WORLD > 2
606 static const EL_BNDRY_VEC *
btb_get_bound_3d(BNDRY_FLAGS * vec,const EL_INFO * el_info,const BAS_FCTS * self)607 btb_get_bound_3d(BNDRY_FLAGS *vec,
608 const EL_INFO *el_info, const BAS_FCTS *self)
609 {
610 static DEF_EL_VEC_BTB(BNDRY, rvec_space);
611 BTB_DATA *data = (BTB_DATA *)self->ext_data;
612 BNDRY_FLAGS *rvec = vec ? vec : rvec_space->vec;
613 int i;
614
615 for (i = 0; i < self->n_bas_fcts; i++) {
616 int wall = data->cur_wall[i];
617 BNDRY_FLAGS_INIT(rvec[i]);
618 BNDRY_FLAGS_SET(rvec[i], el_info->face_bound[wall]);
619 }
620
621 return vec ? NULL : rvec_space;
622 }
623 #endif
624
625 /*******************************************************************************
626 * function for accessing a local DOF_INT_VEC
627 ******************************************************************************/
628
629 DEFUN_GET_EL_DOF_VEC_BTB(int, INT, COPY_EQ);
630
631 /*******************************************************************************
632 * function for accessing a local DOF_REAL_VEC
633 ******************************************************************************/
634
635 DEFUN_GET_EL_DOF_VEC_BTB(real, REAL, COPY_EQ);
636
637 /*******************************************************************************
638 * function for accessing a local DOF_REAL_D_VEC
639 ******************************************************************************/
640
641 DEFUN_GET_EL_DOF_VEC_BTB(real_d, REAL_D, COPY_DOW);
642
643 /*******************************************************************************
644 * function for accessing a local DOF_REAL_VEC_D
645 ******************************************************************************/
646
647 static const EL_REAL_VEC_D *
btb_get_real_vec_d(REAL result[],const EL * el,const DOF_REAL_VEC_D * dv)648 btb_get_real_vec_d(REAL result[],
649 const EL *el, const DOF_REAL_VEC_D *dv)
650 {
651 return (const EL_REAL_VEC_D *)
652 btb_get_real_vec(result, el, (const DOF_REAL_VEC *)dv);
653 }
654
655 /*******************************************************************************
656 * function for accessing a local DOF_SCHAR_VEC
657 ******************************************************************************/
658
659 DEFUN_GET_EL_DOF_VEC_BTB(schar, SCHAR, COPY_EQ);
660
661 /*******************************************************************************
662 * function for accessing a local DOF_UCHAR_VEC
663 ******************************************************************************/
664
665 DEFUN_GET_EL_DOF_VEC_BTB(uchar, UCHAR, COPY_EQ);
666
667 /*******************************************************************************
668 * function for accessing a local DOF_PTR_VEC
669 ******************************************************************************/
670
671 DEFUN_GET_EL_DOF_VEC_BTB(ptr, PTR, COPY_EQ);
672
673 /*******************************************************************************
674 * function for accessing a local DOF_REAL_DD_VEC
675 ******************************************************************************/
676
677 #define MYMCOPY(a, b) MCOPY_DOW((const REAL_D *)(a), (b))
678 DEFUN_GET_EL_DOF_VEC_BTB(real_dd, REAL_DD, MYMCOPY);
679
680 /*******************************************************************************
681 * functions for interpolation:
682 ******************************************************************************/
683
btb_interpol_dow(EL_REAL_VEC_D * vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_D_AT_QP f,void * f_data,const BAS_FCTS * self)684 static void btb_interpol_dow(EL_REAL_VEC_D *vec,
685 const EL_INFO *el_info, int wall,
686 int no, const int *b_no,
687 LOC_FCT_D_AT_QP f, void *f_data,
688 const BAS_FCTS *self)
689 {
690 int i, iq;
691 BTB_DATA *data = (BTB_DATA *)self->ext_data;
692 const WALL_QUAD_FAST *wquad_fast;
693 const QUAD_FAST *qfast;
694 REAL_D f_res;
695 REAL_D uh_res;
696 REAL_D diff;
697 REAL accu;
698
699 vec->n_components = self->n_bas_fcts;
700
701 if (data->wqfast->bas_fcts != self) {
702 data->wqfast = get_wall_quad_fast(self, data->wquad, INIT_PHI);
703 INIT_ELEMENT(el_info, self);
704 }
705 wquad_fast = data->wqfast;
706
707 if (wall >= 0) {
708 for (i = 0; i < self->n_bas_fcts; i++) {
709 if (wall != data->cur_wall[i]) {
710 continue;
711 }
712 if (b_no && *b_no != i) {
713 continue;
714 }
715 vec->vec[i] = 0.0;
716 qfast = wquad_fast->quad_fast[wall];
717 INIT_ELEMENT(el_info, qfast);
718 accu = 0.0;
719 for (iq = 0; iq < qfast->n_points; iq++) {
720 eval_uh_dow_fast(uh_res, vec, qfast, iq);
721 f(f_res, el_info, qfast->quad, iq, f_data);
722 AXPBY_DOW(1.0, f_res, -1.0, uh_res, diff);
723 accu += qfast->w[iq] * SCP_DOW(diff, data->wall_normal[i]);
724 }
725 vec->vec[i] = accu;
726 }
727 } else if (b_no) {
728 for (i = 0; i < no; i++) {
729 int b = b_no[i];
730
731 DEBUG_TEST_EXIT(b < self->n_bas_fcts,
732 "not so many basis functions (%d), only %d\n",
733 b, self->n_bas_fcts);
734
735 wall = data->cur_wall[b];
736 vec->vec[b] = 0.0;
737 qfast = wquad_fast->quad_fast[wall];
738 INIT_ELEMENT(el_info, qfast);
739 accu = 0.0;
740 for (iq = 0; iq < qfast->n_points; iq++) {
741 eval_uh_dow_fast(uh_res, vec, qfast, iq);
742 f(f_res, el_info, qfast->quad, iq, f_data);
743 /*
744 * (int_ei(v*nu_i)-(I_h v)*nu_i)/int_ei(b)
745 */
746 AXPBY_DOW(1.0, f_res, -1.0, uh_res, diff);
747 accu += qfast->w[iq] * SCP_DOW(diff, data->wall_normal[b]);
748 }
749 vec->vec[b] = accu;
750 }
751 } else {
752 for (i = 0; i < self->n_bas_fcts; i++) {
753 wall = data->cur_wall[i];
754 vec->vec[i] = 0.0;
755 qfast = wquad_fast->quad_fast[wall];
756 INIT_ELEMENT(el_info, qfast);
757 accu = 0.0;
758 for (iq = 0; iq < qfast->n_points; iq++) {
759 eval_uh_dow_fast(uh_res, vec, qfast, iq);
760 f(f_res, el_info, qfast->quad, iq, f_data);
761 /*
762 * (int_ei(v*nu_i)-(I_h v)*nu_i)/int_ei(b)
763 */
764 AXPBY_DOW(1.0, f_res, -1.0, uh_res, diff);
765 accu += qfast->w[iq] * SCP_DOW(diff, data->wall_normal[i]);
766 }
767 vec->vec[i] = accu;
768 }
769 }
770 }
771
btb_interpol(EL_REAL_VEC * vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_AT_QP f,void * f_data,const BAS_FCTS * self)772 static void btb_interpol(EL_REAL_VEC *vec,
773 const EL_INFO *el_info, int wall,
774 int no, const int *b_no,
775 LOC_FCT_AT_QP f, void *f_data,
776 const BAS_FCTS *self)
777 {
778 int i, iq;
779 BTB_DATA *data = (BTB_DATA *)self->ext_data;
780 const WALL_QUAD_FAST *wquad_fast;
781 const QUAD_FAST *qfast;
782 REAL f_res;
783 REAL uh_res;
784 REAL accu;
785
786 vec->n_components = self->n_bas_fcts;
787
788 if (data->wqfast->bas_fcts != self) {
789 data->wqfast = get_wall_quad_fast(self, data->wquad, INIT_PHI);
790 INIT_ELEMENT(el_info, self);
791 }
792 wquad_fast = data->wqfast;
793
794 if (wall >= 0) {
795 for (i = 0; i < self->n_bas_fcts; i++) {
796 if (wall != data->cur_wall[i]) {
797 continue;
798 }
799 if (b_no && *b_no != i) {
800 continue;
801 }
802 vec->vec[i] = 0.0;
803 qfast = wquad_fast->quad_fast[wall];
804 INIT_ELEMENT(el_info, qfast);
805 accu = 0.0;
806 for (iq = 0; iq < qfast->n_points; iq++) {
807 uh_res = eval_uh_fast(vec, qfast, iq);
808 f_res = f(el_info, qfast->quad, iq, f_data);
809 accu += qfast->w[iq] * (f_res - uh_res);
810 }
811 vec->vec[i] = accu;
812 }
813 } else if (b_no) {
814 for (i = 0; i < no; i++) {
815 int b = b_no[i];
816
817 DEBUG_TEST_EXIT(b < self->n_bas_fcts,
818 "not so many basis functions (%d), only %d\n",
819 b, self->n_bas_fcts);
820
821 wall = data->cur_wall[b];
822 vec->vec[b] = 0.0;
823 qfast = wquad_fast->quad_fast[wall];
824 INIT_ELEMENT(el_info, qfast);
825 accu = 0.0;
826 for (iq = 0; iq < qfast->n_points; iq++) {
827 uh_res = eval_uh_fast(vec, qfast, iq);
828 f_res = f(el_info, qfast->quad, iq, f_data);
829 /*
830 * (int_ei(v*nu_i)-(I_h v)*nu_i)/int_ei(b)
831 */
832 accu += qfast->w[iq] * (f_res - uh_res);
833 }
834 vec->vec[b] = accu;
835 }
836 } else {
837 for (i = 0; i < self->n_bas_fcts; i++) {
838 wall = data->cur_wall[i];
839 vec->vec[i] = 0.0;
840 qfast = wquad_fast->quad_fast[wall];
841 INIT_ELEMENT(el_info, qfast);
842 accu = 0.0;
843 for (iq = 0; iq < qfast->n_points; iq++) {
844 uh_res = eval_uh_fast(vec, qfast, iq);
845 f_res = f(el_info, qfast->quad, iq, f_data);
846 /*
847 * (int_ei(v*nu_i)-(I_h v)*nu_i)/int_ei(b)
848 */
849 accu += qfast->w[iq] * (f_res - uh_res);
850 }
851 vec->vec[i] = accu;
852 }
853 }
854 }
855
856 /* Called on the trace mesh */
btb_real_refine_inter_d(DOF_REAL_VEC_D * drv,RC_LIST_EL * rclist,int n)857 static void btb_real_refine_inter_d(DOF_REAL_VEC_D *drv,
858 RC_LIST_EL *rclist,
859 int n)
860 {
861 int node, n0;
862 DOF cdof, pdof;
863 EL *el;
864 int i;
865
866 node = drv->fe_space->admin->mesh->node[CENTER];
867 n0 = drv->fe_space->admin->n0_dof[CENTER];
868
869 for (i = 0; i < n; i++) {
870 el = rclist[i].el_info.el;
871 pdof = el->dof[node][n0];
872 cdof = el->child[0]->dof[node][n0];
873 drv->vec[cdof] = 0.5 * drv->vec[pdof];
874 cdof = el->child[1]->dof[node][n0];
875 drv->vec[cdof] = 0.5 * drv->vec[pdof];
876 }
877 }
878
879 /* Called on the trace mesh */
btb_real_coarse_inter_d(DOF_REAL_VEC_D * drv,RC_LIST_EL * rclist,int n)880 static void btb_real_coarse_inter_d(DOF_REAL_VEC_D *drv,
881 RC_LIST_EL *rclist,
882 int n)
883 {
884 int node, n0;
885 DOF cdof, pdof;
886 EL *el;
887 REAL value;
888 int i;
889
890 node = drv->fe_space->admin->mesh->node[CENTER];
891 n0 = drv->fe_space->admin->n0_dof[CENTER];
892
893 for (i = 0; i < n; i++) {
894 el = rclist[i].el_info.el;
895 pdof = el->dof[node][n0];
896 cdof = el->child[0]->dof[node][n0];
897 value = drv->vec[cdof];
898 cdof = el->child[1]->dof[node][n0];
899 value += drv->vec[cdof];
900 drv->vec[pdof] = value;
901 }
902 }
903
904 /* Called on the trace mesh */
btb_real_coarse_restr_d(DOF_REAL_VEC_D * drv,RC_LIST_EL * rclist,int n)905 static void btb_real_coarse_restr_d(DOF_REAL_VEC_D *drv,
906 RC_LIST_EL *rclist,
907 int n)
908 {
909 FUNCNAME("btb_real_coarse_restr_d");
910
911 WARNING("Not implemented.\n");
912 }
913
914 /*******************************************************************************
915 * constructor for the beast:
916 ******************************************************************************/
917
918
919 #define MAX_BUBBLE_INTER_DEG 20 /* Overkill */
920
get_bulk_trace_bubble(unsigned int dim,unsigned int inter_deg,int trace_id)921 const BAS_FCTS *get_bulk_trace_bubble(unsigned int dim,
922 unsigned int inter_deg,
923 int trace_id)
924 {
925 FUNCNAME("get_bulk_trace_bubble");
926 static BAS_FCTS *btb_bfcts[DIM_MAX+1][MAX_BUBBLE_INTER_DEG+1];
927 BTB_DATA *data;
928
929 TEST_EXIT(dim <= DIM_MAX, "dim = %d > DIM_MAX = %d.\n", dim, DIM_MAX);
930
931 if (inter_deg > MAX_BUBBLE_INTER_DEG) {
932 WARNING("Truncating quad-degree from %d to %d.\n",
933 inter_deg, MAX_BUBBLE_INTER_DEG);
934 inter_deg = MAX_BUBBLE_INTER_DEG;
935 }
936
937 if (btb_bfcts[inter_deg][dim] == NULL) {
938 char name[sizeof("BulkTraceBubble@XX_IXX_Xd")];
939
940 sprintf(name, "BulkTraceBubble@%02d_I%02d_%dd",
941 trace_id, inter_deg, dim);
942
943 btb_bfcts[dim][inter_deg] = MEM_CALLOC(1, BAS_FCTS);
944
945 data = btb_bfcts[dim][inter_deg]->ext_data = MEM_CALLOC(1, BTB_DATA);
946
947 btb_bfcts[dim][inter_deg]->name = strdup(name);
948 btb_bfcts[dim][inter_deg]->dim = dim;
949 btb_bfcts[dim][inter_deg]->rdim = DIM_OF_WORLD;
950 btb_bfcts[dim][inter_deg]->degree = N_LAMBDA(dim-1);
951 btb_bfcts[dim][inter_deg]->n_bas_fcts = 0;
952 btb_bfcts[dim][inter_deg]->n_bas_fcts_max = N_WALLS(dim);
953 btb_bfcts[dim][inter_deg]->n_dof[CENTER] = 1; /* trace mesh admin */
954 btb_bfcts[dim][inter_deg]->trace_admin = trace_id;
955 CHAIN_INIT(btb_bfcts[dim][inter_deg]);
956 btb_bfcts[dim][inter_deg]->unchained = btb_bfcts[dim][inter_deg];
957 btb_bfcts[dim][inter_deg]->phi = data->phi;
958 btb_bfcts[dim][inter_deg]->grd_phi = data->grd_phi;
959 btb_bfcts[dim][inter_deg]->D2_phi = data->D2_phi;
960 btb_bfcts[dim][inter_deg]->phi_d = data->phi_d;
961 if (dim > 0) {
962 int o, t, w;
963 btb_bfcts[dim][inter_deg]->trace_bas_fcts =
964 get_trace_bubble(dim-1, inter_deg);
965 for (w = 0; w < N_WALLS(dim); w++) {
966 btb_bfcts[dim][inter_deg]->n_trace_bas_fcts[w] = 1;
967 for (t = 0; t < 2; t++) {
968 for (o = 0; o < 2; o++) {
969 btb_bfcts[dim][inter_deg]->trace_dof_map[t][o][w] =
970 data->trace_dof_map + w;
971 }
972 }
973 }
974 } else {
975 btb_bfcts[dim][inter_deg]->trace_bas_fcts = get_null_bfcts(0);
976 }
977 btb_bfcts[dim][inter_deg]->get_dof_indices = btb_get_dof_indices;
978 switch (dim) {
979 case 1:
980 btb_bfcts[dim][inter_deg]->get_bound = btb_get_bound_1d;
981 break;
982 #if DIM_MAX > 1
983 case 2:
984 btb_bfcts[dim][inter_deg]->get_bound = btb_get_bound_2d;
985 break;
986 #endif
987 #if DIM_MAX > 2
988 case 3:
989 btb_bfcts[dim][inter_deg]->get_bound = btb_get_bound_3d;
990 break;
991 #endif
992 }
993 btb_bfcts[dim][inter_deg]->interpol = btb_interpol;
994 btb_bfcts[dim][inter_deg]->interpol_d = NULL;
995 btb_bfcts[dim][inter_deg]->interpol_dow = btb_interpol_dow;
996 btb_bfcts[dim][inter_deg]->dir_pw_const = true;
997
998 btb_bfcts[dim][inter_deg]->get_int_vec = btb_get_int_vec;
999 btb_bfcts[dim][inter_deg]->get_real_vec = btb_get_real_vec;
1000 btb_bfcts[dim][inter_deg]->get_real_d_vec = btb_get_real_d_vec;
1001 btb_bfcts[dim][inter_deg]->get_real_dd_vec = btb_get_real_dd_vec;
1002 btb_bfcts[dim][inter_deg]->get_real_vec_d = btb_get_real_vec_d;
1003 btb_bfcts[dim][inter_deg]->get_uchar_vec = btb_get_uchar_vec;
1004 btb_bfcts[dim][inter_deg]->get_schar_vec = btb_get_schar_vec;
1005 btb_bfcts[dim][inter_deg]->get_ptr_vec = btb_get_ptr_vec;
1006
1007 btb_bfcts[dim][inter_deg]->real_refine_inter =
1008 (void (*)(DOF_REAL_VEC *, RC_LIST_EL *, int))
1009 btb_real_refine_inter_d;
1010 btb_bfcts[dim][inter_deg]->real_coarse_inter =
1011 (void (*)(DOF_REAL_VEC *, RC_LIST_EL *, int))
1012 btb_real_coarse_inter_d;
1013 btb_bfcts[dim][inter_deg]->real_coarse_restr =
1014 (void (*)(DOF_REAL_VEC *, RC_LIST_EL *, int))
1015 btb_real_coarse_restr_d;
1016
1017 btb_bfcts[dim][inter_deg]->real_refine_inter_d =
1018 btb_real_refine_inter_d;
1019 btb_bfcts[dim][inter_deg]->real_coarse_inter_d =
1020 btb_real_coarse_inter_d;
1021 btb_bfcts[dim][inter_deg]->real_coarse_restr_d =
1022 btb_real_coarse_restr_d;
1023
1024 INIT_ELEMENT_DEFUN(btb_bfcts[dim][inter_deg], btb_init_element,
1025 FILL_NEIGH|FILL_COORDS);
1026 INIT_OBJECT(btb_bfcts[dim][inter_deg]);
1027
1028 data->trace_mesh = NULL;
1029 data->trace_id = trace_id;
1030 data->wquad = get_wall_quad(dim, inter_deg);
1031 data->inter_deg = inter_deg;
1032 data->wqfast =
1033 get_wall_quad_fast(btb_bfcts[dim][inter_deg], data->wquad, INIT_PHI);
1034 }
1035
1036 return btb_bfcts[dim][inter_deg];
1037 }
1038
1039