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: l2scp.c */
7 /* */
8 /* description: compute L2 scalar products of the right hand side and */
9 /* the basis functions of the fe-space */
10 /* sets Dirichlet boundary values */
11 /* */
12 /*--------------------------------------------------------------------------*/
13 /* */
14 /* authors: Alfred Schmidt */
15 /* Zentrum fuer Technomathematik */
16 /* Fachbereich 3 Mathematik/Informatik */
17 /* Universitaet Bremen */
18 /* Bibliothekstr. 2 */
19 /* D-28359 Bremen, Germany */
20 /* */
21 /* Kunibert G. Siebert */
22 /* Institut fuer Mathematik */
23 /* Universitaet Augsburg */
24 /* Universitaetsstr. 14 */
25 /* D-86159 Augsburg, Germany */
26 /* */
27 /* Claus-Justus Heine */
28 /* Abteilung fuer Angewandte Mathematik */
29 /* Universitaet Freiburg */
30 /* Hermann-Herder-Strasse 10 */
31 /* 79104 Freiburg, Germany */
32 /* */
33 /* http://www.alberta-fem.de */
34 /* */
35 /* (c) by A. Schmidt and K.G. Siebert (1996-2003), */
36 /* C.-J. Heine (2006-2009). */
37 /* */
38 /*--------------------------------------------------------------------------*/
39
40 #include "alberta_intern.h"
41 #include "alberta.h"
42
43 /*--------------------------------------------------------------------------*/
44 /* calculates the L2 scalar product of the function f with each basis */
45 /* function of the mesh belonging to the dof admin of fh; the scalar */
46 /* products are stored in fh; */
47 /* */
48 /* fh->vec[j] = (f, phi_j)_\Omega */
49 /* */
50 /* where phi_j are the basis functions belonging to fh->fe_space->admin */
51 /* */
52 /*--------------------------------------------------------------------------*/
53
_AI_L2scp_fct_bas(DOF_REAL_VEC * fh,LOC_FCT_AT_QP f_loc,void * fd,FCT_AT_X f,FLAGS fill_flag,const QUAD * quad)54 static void _AI_L2scp_fct_bas(DOF_REAL_VEC *fh,
55 LOC_FCT_AT_QP f_loc, void *fd, FCT_AT_X f,
56 FLAGS fill_flag,
57 const QUAD *quad)
58 {
59 FUNCNAME("L2scp_fct_bas");
60 MESH *mesh = NULL;
61 const QUAD_FAST *quad_fast;
62 const FE_SPACE *fe_space;
63 const REAL *w;
64 const REAL_B *lambda;
65 REAL val, det;
66 REAL_D x;
67 int iq, j, n_phi, dim;
68 const BAS_FCTS *bas_fcts;
69 INIT_EL_TAG tag = INIT_EL_TAG_DFLT;
70
71 TEST_EXIT(fh, "no DOF_REAL_VEC fh\n");
72 if (!f && !f_loc)
73 return;
74
75 TEST_EXIT(fh->fe_space, "no fe_space in DOF_REAL_VEC %s\n", NAME(fh));
76
77 GET_STRUCT(mesh, fh->fe_space);
78
79 dim = mesh->dim;
80
81 fe_space = fh->fe_space;
82 bas_fcts = fe_space->bas_fcts;
83
84 if (!quad) {
85 quad = get_quadrature(dim, 2*fh->fe_space->bas_fcts->degree-2);
86 }
87 quad_fast = get_quad_fast(fh->fe_space->bas_fcts, quad, INIT_PHI);
88
89 w = quad_fast->w;
90 lambda = quad_fast->quad->lambda;
91 n_phi = bas_fcts->n_bas_fcts;
92
93 fill_flag |= CALL_LEAF_EL|FILL_COORDS;
94 if (mesh->is_periodic && !(fh->fe_space->admin->flags & ADM_PERIODIC)) {
95 fill_flag |= FILL_NON_PERIODIC;
96 }
97 fill_flag |= quad_fast->fill_flags;
98
99 {
100 PARAMETRIC *parametric = mesh->parametric;
101 int is_parametric = false;
102 int dim = mesh->dim;
103 REAL det_p[quad->n_points_max];
104 REAL wdetf_qp[quad->n_points_max];
105 REAL_D x_p[quad->n_points_max];
106 const EL_DOF_VEC *dof;
107
108 TRAVERSE_FIRST(mesh, -1, fill_flag) {
109
110 INIT_ELEMENT_SETUP(el_info, quad_fast, tag, continue, {
111 w = quad_fast->w;
112 lambda = quad_fast->quad->lambda;
113 n_phi = bas_fcts->n_bas_fcts;
114 });
115
116 if (parametric) {
117 is_parametric = parametric->init_element(el_info, parametric);
118 }
119
120 if (is_parametric) {
121 parametric->det(el_info, quad, 0, NULL, det_p);
122
123 if (f) {
124 parametric->coord_to_world(el_info, quad, 0, NULL, x_p);
125 for (iq = 0; iq < quad->n_points; iq++) {
126 wdetf_qp[iq] = w[iq]*det_p[iq]*f(x_p[iq]);
127 }
128 } else {
129 for (iq = 0; iq < quad->n_points; iq++) {
130 wdetf_qp[iq] = w[iq]*det_p[iq]*f_loc(el_info, quad, iq, fd);
131 }
132 }
133 } else {
134 det = el_det_dim(dim, el_info);
135 if (f) {
136 for (iq = 0; iq < quad->n_points; iq++) {
137 coord_to_world(el_info, lambda[iq], x);
138 wdetf_qp[iq] = w[iq]*det*f(x);
139 }
140 } else {
141 for (iq = 0; iq < quad->n_points; iq++) {
142 wdetf_qp[iq] = w[iq]*det*f_loc(el_info, quad, iq, fd);
143 }
144 }
145 }
146
147 INIT_ELEMENT(el_info, quad_fast);
148
149 CHAIN_DO(quad_fast, const QUAD_FAST) {
150 const BAS_FCTS *bas_fcts;
151 const REAL *const*phi;
152
153 bas_fcts = quad_fast->bas_fcts;
154 phi = quad_fast->phi;
155 n_phi = bas_fcts->n_bas_fcts;
156
157 dof = GET_DOF_INDICES(bas_fcts, el_info->el, fh->fe_space->admin, NULL);
158
159 for (j = 0; j < n_phi; j++) {
160 val = 0.0;
161 for (iq = 0; iq < quad->n_points; iq++) {
162 val += phi[iq][j]*wdetf_qp[iq];
163 }
164 fh->vec[dof->vec[j]] += val;
165 }
166 fh = CHAIN_NEXT(fh, DOF_REAL_VEC);
167 } CHAIN_WHILE(quad_fast, const QUAD_FAST);
168 } TRAVERSE_NEXT();
169 }
170 }
171
172 __FLATTEN_ATTRIBUTE__
L2scp_fct_bas(FCT_AT_X f,const QUAD * quad,DOF_REAL_VEC * fh)173 void L2scp_fct_bas(FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh)
174 {
175 _AI_L2scp_fct_bas(fh, NULL, NULL, f, FILL_NOTHING, quad);
176 }
177
178 __FLATTEN_ATTRIBUTE__
L2scp_fct_bas_loc(DOF_REAL_VEC * fh,LOC_FCT_AT_QP f,void * fd,FLAGS fill_flag,const QUAD * quad)179 void L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
180 LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
181 const QUAD *quad)
182 {
183 _AI_L2scp_fct_bas(fh, f, fd, NULL, fill_flag, quad);
184 }
185
_AI_L2scp_fct_bas_dow(DOF_REAL_VEC_D * fh,FCT_D_AT_X f,LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flag,const QUAD * quad)186 static void _AI_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
187 FCT_D_AT_X f,
188 LOC_FCT_D_AT_QP f_loc, void *fd,
189 FLAGS fill_flag,
190 const QUAD *quad)
191 {
192 FUNCNAME("L2scp_fct_bas_dow");
193 MESH *mesh = NULL;
194 const QUAD_FAST *quad_fast;
195 const FE_SPACE *fe_space;
196 REAL_D x;
197 REAL det;
198 const REAL *w;
199 const REAL_B *lambda;
200 REAL_D val;
201 int iq, j, n_phi, dim, n_points;
202 const BAS_FCTS *bas_fcts;
203 INIT_EL_TAG qtag = INIT_EL_TAG_DFLT;
204
205 TEST_EXIT(fh, "no DOF_REAL_VEC fh\n");
206 if (!f && !f_loc)
207 return;
208
209 TEST_EXIT(fh->fe_space, "no fe_space in DOF_REAL_D_VEC \"%s\"\n", NAME(fh));
210
211 TEST_EXIT(fh->fe_space->rdim == DIM_OF_WORLD,
212 "Called for scalar finite element space \"%s\".\n",
213 NAME(fh->fe_space));
214
215 GET_STRUCT(mesh, fh->fe_space);
216
217 dim = mesh->dim;
218
219 fe_space = fh->fe_space;
220 bas_fcts = fe_space->bas_fcts;
221
222 if (!quad) {
223 quad = get_quadrature(dim, 2*bas_fcts->degree-2);
224 }
225 quad_fast = get_quad_fast(bas_fcts, quad, INIT_PHI);
226
227 w = quad->w;
228 lambda = quad->lambda;
229 n_points = quad->n_points;
230
231 fill_flag |= CALL_LEAF_EL|FILL_COORDS;
232 if (mesh->is_periodic && !(fh->fe_space->admin->flags & ADM_PERIODIC)) {
233 fill_flag |= FILL_NON_PERIODIC;
234 }
235 fill_flag |= quad_fast->fill_flags;
236
237 {
238 PARAMETRIC *parametric = mesh->parametric;
239 int is_parametric = false;
240 int dim = mesh->dim;
241 REAL det_p[quad->n_points_max];
242 REAL_D wdetf_qp[quad->n_points_max];
243 REAL_D x_p[quad->n_points_max];
244 const EL_DOF_VEC *dof;
245
246 TRAVERSE_FIRST(mesh, -1, fill_flag) {
247
248 INIT_ELEMENT_SETUP(el_info, quad, qtag, continue, {
249 w = quad->w;
250 n_points = quad->n_points;
251 lambda = quad->lambda;
252 });
253
254 if (parametric) {
255 is_parametric = parametric->init_element(el_info, parametric);
256 }
257
258 if (is_parametric) {
259 parametric->det(el_info, quad, -1, NULL, det_p);
260 if (f) {
261 parametric->coord_to_world(el_info, quad, -1, NULL, x_p);
262 for (iq = 0; iq < quad->n_points; iq++) {
263 AXEY_DOW(w[iq]*det_p[iq], f(x_p[iq], wdetf_qp[iq]), wdetf_qp[iq]);
264 }
265 } else {
266 for (iq = 0; iq < quad->n_points; iq++) {
267 AXEY_DOW(w[iq]*det_p[iq],
268 f_loc(wdetf_qp[iq], el_info, quad, iq, fd),
269 wdetf_qp[iq]);
270 }
271 }
272 } else {
273 det = el_det_dim(dim, el_info);
274
275 if (f) {
276 for (iq = 0; iq < quad->n_points; iq++) {
277 coord_to_world(el_info, lambda[iq], x);
278 AXEY_DOW(w[iq]*det, f(x, wdetf_qp[iq]), wdetf_qp[iq]);
279 }
280 } else {
281 for (iq = 0; iq < quad->n_points; iq++) {
282 AXEY_DOW(w[iq]*det,
283 f_loc(wdetf_qp[iq], el_info, quad, iq, fd),
284 wdetf_qp[iq]);
285 }
286 }
287 }
288
289 INIT_ELEMENT(el_info, quad_fast);
290
291 CHAIN_DO(quad_fast, const QUAD_FAST) {
292 const BAS_FCTS *bas_fcts;
293
294 bas_fcts = quad_fast->bas_fcts;
295 n_phi = bas_fcts->n_bas_fcts;
296
297 dof = GET_DOF_INDICES(bas_fcts, el_info->el, fh->fe_space->admin, NULL);
298
299 if (fh->stride != 1) {
300 /* Cartesian product space */
301 DOF_REAL_D_VEC *fhd = (DOF_REAL_D_VEC *)fh;
302
303 for (j = 0; j < n_phi; j++) {
304 REAL_D val = { 0.0, };
305 for (iq = 0; iq < quad->n_points; iq++) {
306 AXPY_DOW(quad_fast->phi[iq][j], wdetf_qp[iq], val);
307 }
308 AXPY_DOW(1.0, val, fhd->vec[dof->vec[j]]);
309 }
310 } else {
311 if (bas_fcts->dir_pw_const) {
312 for (j = 0; j < n_phi; j++) {
313 SET_DOW(0.0, val);
314 for (iq = 0; iq < n_points; iq++) {
315 AXPY_DOW(quad_fast->phi[iq][j], wdetf_qp[iq], val);
316 }
317 fh->vec[dof->vec[j]] +=
318 SCP_DOW(val, PHI_D(quad_fast->bas_fcts, j, NULL));
319 }
320 } else {
321 for (j = 0; j < n_phi; j++) {
322 REAL val = 0.0;
323 for (iq = 0; iq < n_points; iq++) {
324 val +=
325 quad_fast->phi[iq][j]
326 *
327 SCP_DOW(wdetf_qp[iq],
328 PHI_D(quad_fast->bas_fcts, j, lambda[iq]));
329 }
330 fh->vec[dof->vec[j]] += val;
331 }
332 }
333 }
334 fh = CHAIN_NEXT(fh, DOF_REAL_VEC_D);
335 } CHAIN_WHILE(quad_fast, const QUAD_FAST);
336
337 } TRAVERSE_NEXT();
338 }
339 }
340
341 __FLATTEN_ATTRIBUTE__
L2scp_fct_bas_dow(FCT_D_AT_X f,const QUAD * quad,DOF_REAL_VEC_D * fh)342 void L2scp_fct_bas_dow(FCT_D_AT_X f, const QUAD *quad, DOF_REAL_VEC_D *fh)
343 {
344 _AI_L2scp_fct_bas_dow(fh, f, NULL, NULL, FILL_NOTHING, quad);
345 }
346
347 __FLATTEN_ATTRIBUTE__
L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D * fh,LOC_FCT_D_AT_QP f,void * fd,FLAGS fill_flag,const QUAD * quad)348 void L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
349 LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
350 const QUAD *quad)
351 {
352 _AI_L2scp_fct_bas_dow(fh, NULL, f, fd, fill_flag, quad);
353 }
354
_AI_H1scp_fct_bas(DOF_REAL_VEC * fh,GRD_FCT_AT_X f,GRD_LOC_FCT_AT_QP f_loc,void * fd,FLAGS fill_flag,const QUAD * quad)355 static void _AI_H1scp_fct_bas(DOF_REAL_VEC *fh,
356 GRD_FCT_AT_X f,
357 GRD_LOC_FCT_AT_QP f_loc,
358 void *fd /* userdata */,
359 FLAGS fill_flag,
360 const QUAD *quad)
361 {
362 FUNCNAME("H1scp_fct_bas");
363 MESH *mesh = NULL;
364 const QUAD_FAST *quad_fast;
365 const FE_SPACE *fe_space;
366 const REAL *w;
367 int iq, j, alpha, n_phi, dim, n_points;
368 const BAS_FCTS *bas_fcts;
369 INIT_EL_TAG qtag = INIT_EL_TAG_DFLT;
370
371 TEST_EXIT(fh, "no DOF_REAL_VEC fh\n");
372 if (!f && !f_loc)
373 return;
374
375 TEST_EXIT(fh->fe_space, "no fe_space in DOF_REAL_VEC_D \"%s\"\n", NAME(fh));
376
377 TEST_EXIT(fh->fe_space->rdim == 1,
378 "Called for vector valued finite element space \"%s\".\n",
379 NAME(fh->fe_space));
380
381 GET_STRUCT(mesh, fh->fe_space);
382
383 dim = mesh->dim;
384
385 fe_space = fh->fe_space;
386 bas_fcts = fe_space->bas_fcts;
387
388 if (!quad) {
389 quad = get_quadrature(dim, 2*bas_fcts->degree-2);
390 }
391 quad_fast = get_quad_fast(bas_fcts, quad, INIT_GRD_PHI);
392
393 w = quad->w;
394 n_points = quad->n_points;
395
396 fill_flag |= CALL_LEAF_EL|FILL_COORDS;
397 if (mesh->is_periodic && !(fh->fe_space->admin->flags & ADM_PERIODIC)) {
398 fill_flag |= FILL_NON_PERIODIC;
399 }
400 fill_flag |= quad_fast->fill_flags;
401
402 {
403 PARAMETRIC *parametric = mesh->parametric;
404 int is_parametric = false;
405 int dim = mesh->dim;
406 REAL_B wdetgrdf_qp[quad->n_points_max];
407 const EL_DOF_VEC *dof;
408
409 TRAVERSE_FIRST(mesh, -1, fill_flag) {
410
411 INIT_ELEMENT_SETUP(el_info, quad, qtag, continue, {
412 w = quad->w;
413 n_points = quad->n_points;
414 });
415
416 if (parametric) {
417 is_parametric = parametric->init_element(el_info, parametric);
418 }
419
420 if (is_parametric) {
421 const QUAD_EL_CACHE *qelc;
422
423 qelc = fill_quad_el_cache(el_info, quad,
424 FILL_EL_QUAD_DET|FILL_EL_QUAD_LAMBDA);
425 if (f) {
426 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
427 for (iq = 0; iq < quad->n_points; iq++) {
428 REAL_D grd;
429
430 f(qelc->world[iq], grd);
431 SCAL_DOW(w[iq]*qelc->param.det[iq], grd);
432 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
433 wdetgrdf_qp[iq][alpha] =
434 SCP_DOW(grd, qelc->param.Lambda[iq][alpha]);
435 }
436 }
437 } else {
438 for (iq = 0; iq < quad->n_points; iq++) {
439 REAL_D grd;
440
441 f_loc(grd, el_info,
442 (const REAL_D *)qelc->param.Lambda[iq], quad, iq, fd);
443 SCAL_DOW(w[iq]*qelc->param.det[iq], grd);
444 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
445 wdetgrdf_qp[iq][alpha] =
446 SCP_DOW(grd, qelc->param.Lambda[iq][alpha]);
447 }
448 }
449 }
450 } else {
451 if (f) {
452 const QUAD_EL_CACHE *qelc;
453 const EL_GEOM_CACHE *elgc;
454
455 elgc = fill_el_geom_cache(el_info, FILL_EL_DET|FILL_EL_LAMBDA);
456 qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
457 for (iq = 0; iq < quad->n_points; iq++) {
458 REAL_D grd;
459
460 f(qelc->world[iq], grd);
461 SCAL_DOW(w[iq]*elgc->det, grd);
462 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
463 wdetgrdf_qp[iq][alpha] =
464 SCP_DOW(grd, elgc->Lambda[alpha]);
465 }
466 }
467 } else {
468 const EL_GEOM_CACHE *elgc;
469
470 elgc = fill_el_geom_cache(el_info, FILL_EL_DET|FILL_EL_LAMBDA);
471 for (iq = 0; iq < quad->n_points; iq++) {
472 REAL_D grd;
473 f_loc(grd, el_info, elgc->Lambda, quad, iq, fd);
474 SCAL_DOW(w[iq]*elgc->det, grd);
475 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
476 wdetgrdf_qp[iq][alpha] =
477 SCP_DOW(grd, elgc->Lambda[alpha]);
478 }
479 }
480 }
481 }
482
483 INIT_ELEMENT(el_info, quad_fast);
484
485 CHAIN_DO(quad_fast, const QUAD_FAST) {
486 const BAS_FCTS *bas_fcts;
487 const REAL_B *const*grd_phi = quad_fast->grd_phi;
488 REAL accu;
489
490 bas_fcts = quad_fast->bas_fcts;
491 n_phi = bas_fcts->n_bas_fcts;
492
493 dof = GET_DOF_INDICES(bas_fcts, el_info->el, fh->fe_space->admin, NULL);
494
495 for (j = 0; j < n_phi; j++) {
496 accu = 0.0;
497 for (iq = 0; iq < n_points; iq++) {
498 accu += SCP_BAR(dim, wdetgrdf_qp[iq], grd_phi[iq][j]);
499 }
500 fh->vec[dof->vec[j]] += accu;
501 }
502 fh = CHAIN_NEXT(fh, DOF_REAL_VEC);
503 } CHAIN_WHILE(quad_fast, const QUAD_FAST);
504
505 } TRAVERSE_NEXT();
506 }
507 }
508 __FLATTEN_ATTRIBUTE__
H1scp_fct_bas(GRD_FCT_AT_X f,const QUAD * quad,DOF_REAL_VEC * fh)509 void H1scp_fct_bas(GRD_FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh)
510 {
511 _AI_H1scp_fct_bas(fh, f, NULL, NULL,FILL_NOTHING, quad);
512 }
513
514 __FLATTEN_ATTRIBUTE__
H1scp_fct_bas_loc(DOF_REAL_VEC * fh,GRD_LOC_FCT_AT_QP f,void * fd,FLAGS fill_flag,const QUAD * quad)515 void H1scp_fct_bas_loc(DOF_REAL_VEC *fh,
516 GRD_LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
517 const QUAD *quad)
518 {
519 _AI_H1scp_fct_bas(fh, NULL, f, fd, fill_flag, quad);
520 }
521
_AI_H1scp_fct_bas_dow(DOF_REAL_VEC_D * fh,GRD_FCT_D_AT_X f,GRD_LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flag,const QUAD * quad)522 static void _AI_H1scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
523 GRD_FCT_D_AT_X f,
524 GRD_LOC_FCT_D_AT_QP f_loc,
525 void *fd /* userdata */,
526 FLAGS fill_flag,
527 const QUAD *quad)
528 {
529 FUNCNAME("H1scp_fct_bas_dow");
530 MESH *mesh = NULL;
531 const QUAD_FAST *quad_fast;
532 const FE_SPACE *fe_space;
533 const REAL *w;
534 int iq, i, j, alpha, n_phi, dim, n_points;
535 const BAS_FCTS *bas_fcts;
536 INIT_EL_TAG qtag = INIT_EL_TAG_DFLT;
537
538 TEST_EXIT(fh, "no DOF_REAL_VEC fh\n");
539 if (!f && !f_loc)
540 return;
541
542 TEST_EXIT(fh->fe_space, "no fe_space in DOF_REAL_VEC_D \"%s\"\n", NAME(fh));
543
544 TEST_EXIT(fh->fe_space->rdim == DIM_OF_WORLD,
545 "Called for scalar finite element space \"%s\".\n",
546 NAME(fh->fe_space));
547
548 GET_STRUCT(mesh, fh->fe_space);
549
550 dim = mesh->dim;
551
552 fe_space = fh->fe_space;
553 bas_fcts = fe_space->bas_fcts;
554
555 if (!quad) {
556 quad = get_quadrature(dim, 2*bas_fcts->degree-2);
557 }
558 quad_fast = get_quad_fast(bas_fcts, quad, INIT_GRD_PHI);
559
560 w = quad->w;
561 n_points = quad->n_points;
562
563 fill_flag |= CALL_LEAF_EL|FILL_COORDS;
564 if (mesh->is_periodic && !(fh->fe_space->admin->flags & ADM_PERIODIC)) {
565 fill_flag |= FILL_NON_PERIODIC;
566 }
567 fill_flag |= quad_fast->fill_flags;
568
569 {
570 PARAMETRIC *parametric = mesh->parametric;
571 int is_parametric = false;
572 int dim = mesh->dim;
573 REAL_DB wdetgrdf_qp[quad->n_points_max];
574 const EL_DOF_VEC *dof;
575
576 TRAVERSE_FIRST(mesh, -1, fill_flag) {
577
578 INIT_ELEMENT_SETUP(el_info, quad, qtag, continue, {
579 w = quad->w;
580 n_points = quad->n_points;
581 });
582
583 if (parametric) {
584 is_parametric = parametric->init_element(el_info, parametric);
585 }
586
587 if (is_parametric) {
588 const QUAD_EL_CACHE *qelc;
589
590 qelc = fill_quad_el_cache(el_info, quad,
591 FILL_EL_QUAD_DET|FILL_EL_QUAD_LAMBDA);
592 if (f) {
593 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
594 for (iq = 0; iq < quad->n_points; iq++) {
595 REAL_DD grd;
596
597 f(qelc->world[iq], grd);
598 MSCAL_DOW(w[iq]*qelc->param.det[iq], grd);
599 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
600 for (i = 0; i < DIM_OF_WORLD; i++) {
601 wdetgrdf_qp[iq][i][alpha] =
602 SCP_DOW(grd[i], qelc->param.Lambda[iq][alpha]);
603 }
604 }
605 }
606 } else {
607 for (iq = 0; iq < quad->n_points; iq++) {
608 REAL_DD grd;
609
610 f_loc(grd, el_info,
611 (const REAL_D *)qelc->param.Lambda[iq], quad, iq, fd);
612 MSCAL_DOW(w[iq]*qelc->param.det[iq], grd);
613 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
614 for (i = 0; i < DIM_OF_WORLD; i++) {
615 wdetgrdf_qp[iq][i][alpha] =
616 SCP_DOW(grd[i], qelc->param.Lambda[iq][alpha]);
617 }
618 }
619 }
620 }
621 } else {
622 if (f) {
623 const QUAD_EL_CACHE *qelc;
624 const EL_GEOM_CACHE *elgc;
625
626 elgc = fill_el_geom_cache(el_info, FILL_EL_DET|FILL_EL_LAMBDA);
627 qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
628 for (iq = 0; iq < quad->n_points; iq++) {
629 REAL_DD grd;
630
631 f(qelc->world[iq], grd);
632 MSCAL_DOW(w[iq]*elgc->det, grd);
633 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
634 for (i = 0; i < DIM_OF_WORLD; i++) {
635 wdetgrdf_qp[iq][i][alpha] =
636 SCP_DOW(grd[i], elgc->Lambda[alpha]);
637 }
638 }
639 }
640 } else {
641 const EL_GEOM_CACHE *elgc;
642
643 elgc = fill_el_geom_cache(el_info, FILL_EL_DET|FILL_EL_LAMBDA);
644 for (iq = 0; iq < quad->n_points; iq++) {
645 REAL_DD grd;
646 f_loc(grd, el_info, elgc->Lambda, quad, iq, fd);
647 MSCAL_DOW(w[iq]*elgc->det, grd);
648 for (alpha = 0; alpha < N_LAMBDA(dim); alpha++) {
649 for (i = 0; i < DIM_OF_WORLD; i++) {
650 wdetgrdf_qp[iq][i][alpha] =
651 SCP_DOW(grd[i], elgc->Lambda[alpha]);
652 }
653 }
654 }
655 }
656 }
657
658 INIT_ELEMENT(el_info, quad_fast);
659
660 CHAIN_DO(quad_fast, const QUAD_FAST) {
661 const BAS_FCTS *bas_fcts;
662
663 bas_fcts = quad_fast->bas_fcts;
664 n_phi = bas_fcts->n_bas_fcts;
665
666 dof = GET_DOF_INDICES(bas_fcts, el_info->el, fh->fe_space->admin, NULL);
667
668 if (fh->stride != 1) {
669 /* Cartesian product space */
670 DOF_REAL_D_VEC *fhd = (DOF_REAL_D_VEC *)fh;
671 const REAL_B *const*grd_phi = quad_fast->grd_phi;
672
673 for (j = 0; j < n_phi; j++) {
674 REAL_D val = { 0.0, };
675 for (iq = 0; iq < quad->n_points; iq++) {
676 for (i = 0; i < DIM_OF_WORLD; i++) {
677 val[i] += SCP_BAR(dim, grd_phi[iq][j], wdetgrdf_qp[iq][i]);
678 }
679 }
680 AXPY_DOW(1.0, val, fhd->vec[dof->vec[j]]);
681 }
682 } else {
683 const REAL_DB *const*grd_phi;
684 REAL accu;
685
686 grd_phi = get_quad_fast_grd_phi_dow(quad_fast);
687
688 for (j = 0; j < n_phi; j++) {
689 accu = 0.0;
690 for (iq = 0; iq < n_points; iq++) {
691 for (i = 0; i < DIM_OF_WORLD; i++) {
692 accu += SCP_BAR(dim, wdetgrdf_qp[iq][i], grd_phi[iq][j][i]);
693 }
694 }
695 fh->vec[dof->vec[j]] += accu;
696 }
697 }
698 fh = CHAIN_NEXT(fh, DOF_REAL_VEC_D);
699 } CHAIN_WHILE(quad_fast, const QUAD_FAST);
700
701 } TRAVERSE_NEXT();
702 }
703 }
704
705 __FLATTEN_ATTRIBUTE__
H1scp_fct_bas_dow(GRD_FCT_D_AT_X f,const QUAD * quad,DOF_REAL_VEC_D * fh)706 void H1scp_fct_bas_dow(GRD_FCT_D_AT_X f, const QUAD *quad, DOF_REAL_VEC_D *fh)
707 {
708 _AI_H1scp_fct_bas_dow(fh, f, NULL, NULL, FILL_NOTHING, quad);
709 }
710
711 __FLATTEN_ATTRIBUTE__
H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D * fh,GRD_LOC_FCT_D_AT_QP f,void * fd,FLAGS fill_flag,const QUAD * quad)712 void H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
713 GRD_LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
714 const QUAD *quad)
715 {
716 _AI_H1scp_fct_bas_dow(fh, NULL, f, fd, fill_flag, quad);
717 }
718
719 /* Compute the L2 scalar product over the boundary of the domain.
720 * Return "true" if not all boundary segments of the mesh belong to
721 * the segment specified by "bndry_seg". If "bndry_seg == NULL" then
722 * the scalar product is computed over the entire boundary (i.e. over
723 * all walls without neighbour).
724 */
_AI_bndry_L2scp_fct_bas(DOF_REAL_VEC * fh,LOC_FCT_AT_QP f_loc,void * fd,FLAGS fill_flag,REAL (* f)(const REAL_D x,const REAL_D normal),const BNDRY_FLAGS bndry_seg,const WALL_QUAD * wall_quad)725 static bool _AI_bndry_L2scp_fct_bas(DOF_REAL_VEC *fh,
726 LOC_FCT_AT_QP f_loc,
727 void *fd, FLAGS fill_flag,
728 REAL (*f)(const REAL_D x,
729 const REAL_D normal),
730 const BNDRY_FLAGS bndry_seg,
731 const WALL_QUAD *wall_quad)
732 {
733 MESH *mesh = NULL;
734 const FE_SPACE *fe_space;
735 const BAS_FCTS *bas_fcts;
736 const WALL_QUAD_FAST *wqfast;
737 const QUAD_FAST *qfast;
738 REAL val;
739 int iq, not_all = false;
740
741 if (!f && !f_loc) {
742 return not_all;
743 }
744
745 if (!fh) {
746 return not_all;
747 }
748
749 if (BNDRY_FLAGS_IS_INTERIOR(bndry_seg)) {
750 /* nothing to do */
751 return not_all;
752 }
753
754 fe_space = fh->fe_space;
755 bas_fcts = fe_space->bas_fcts;
756 mesh = fe_space->mesh;
757
758 if (wall_quad == NULL) {
759 wall_quad = get_wall_quad(mesh->dim, 2*bas_fcts->degree);
760 }
761 wqfast = get_wall_quad_fast(bas_fcts, wall_quad, INIT_PHI);
762
763 fill_flag |= CALL_LEAF_EL|FILL_COORDS|FILL_MACRO_WALLS;
764
765 if (mesh->is_periodic && !(fe_space->admin->flags & ADM_PERIODIC)) {
766 fill_flag |= FILL_NON_PERIODIC;
767 }
768 fill_flag |= wqfast->fill_flags;
769
770 {
771 PARAMETRIC *parametric = mesh->parametric;
772 int is_parametric = false;
773 int dim = mesh->dim;
774 REAL wdetf_qp[wall_quad->n_points_max];
775 const REAL *w[N_WALLS_MAX];
776 int n_points[N_WALLS_MAX];
777 INIT_EL_TAG wall_tag[N_WALLS_MAX];
778 EL_DOF_VEC *dof;
779 int wall;
780
781 for (wall = 0; wall < N_WALLS(dim); wall++) {
782 wall_tag[wall] = INIT_EL_TAG_DFLT;
783 w[wall] = wqfast->quad_fast[wall]->w;
784 n_points[wall] = wall_quad->quad[wall].n_points;
785 }
786
787 dof = get_el_dof_vec(bas_fcts);
788
789 TRAVERSE_FIRST(mesh, -1, fill_flag) {
790 BNDRY_TYPE wall_bnd;
791 int walls[N_WALLS_MAX];
792 int n_walls, wnr, wall;
793
794 n_walls = 0;
795 for (wnr = 0; wnr < N_WALLS(dim); wnr++) {
796 wall_bnd = wall_bound(el_info, wnr);
797 if (wall_bnd == INTERIOR) {
798 continue;
799 }
800 if (bndry_seg == NULL ||
801 BNDRY_FLAGS_IS_AT_BNDRY(bndry_seg, wall_bnd)) {
802 walls[n_walls++] = wnr;
803 } else {
804 not_all = true;
805 }
806 }
807 if (n_walls == 0) {
808 continue;
809 }
810
811 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
812 continue;
813 }
814
815 /* We assume that it is somewhat likely that we have to do some
816 * work here and fetch the local DOF-indices here.
817 */
818 get_dof_indices(dof, fe_space, el_info->el);
819
820 if (parametric) {
821 is_parametric = parametric->init_element(el_info, parametric);
822 }
823
824 for (wnr = 0; wnr < n_walls; wnr++) {
825
826 wall = walls[wnr];
827
828 INIT_ELEMENT_SETUP(el_info,
829 wqfast->quad_fast[wall], wall_tag[wall], continue, {
830 w[wall] = wqfast->quad_fast[wall]->w;
831 n_points[wall] = wall_quad->quad[wall].n_points;
832 });
833
834 if (is_parametric) {
835 const QUAD_EL_CACHE *qelc;
836 if (f) {
837 qelc = fill_quad_el_cache(el_info, &wall_quad->quad[wall],
838 FILL_EL_QUAD_WALL_NORMAL|
839 FILL_EL_QUAD_WALL_DET|
840 FILL_EL_QUAD_WORLD);
841 for (iq = 0; iq < n_points[wall]; iq++) {
842 wdetf_qp[iq] = (w[wall][iq]*qelc->param.wall_det[iq]
843 *
844 f(qelc->world[iq], qelc->param.wall_normal[iq]));
845 }
846 } else {
847 qelc = fill_quad_el_cache(el_info, &wall_quad->quad[wall],
848 FILL_EL_QUAD_WALL_NORMAL|
849 FILL_EL_QUAD_WALL_DET);
850 for (iq = 0; iq < n_points[wall]; iq++) {
851 wdetf_qp[iq] =
852 w[wall][iq] * qelc->param.wall_det[iq]
853 *
854 f_loc(el_info, &wall_quad->quad[wall], iq, fd);
855 }
856 }
857 } else {
858 const EL_GEOM_CACHE *elgc =
859 fill_el_geom_cache(el_info, FILL_EL_WALL_DET(wall));
860 if (f) {
861 const QUAD_EL_CACHE *qelc =
862 fill_quad_el_cache(el_info,
863 &wall_quad->quad[wall], FILL_EL_QUAD_WORLD);
864 for (iq = 0; iq < n_points[wall]; iq++) {
865 wdetf_qp[iq] =
866 w[wall][iq]*elgc->wall_det[wall]
867 *
868 f(qelc->world[iq], elgc->wall_normal[wall]);
869 }
870 } else {
871 for (iq = 0; iq < n_points[wall]; iq++) {
872 wdetf_qp[iq] =
873 w[wall][iq] * elgc->wall_det[wall]
874 *
875 f_loc(el_info, &wall_quad->quad[wall], iq, fd);
876 }
877 }
878 }
879
880 /*****************************/
881 qfast = wqfast->quad_fast[wall];
882 INIT_ELEMENT(el_info, qfast); /* must be: tag != INIT_EL_TAG_NULL */
883
884 CHAIN_DO(qfast, const QUAD_FAST) {
885 const REAL *const*phi;
886 int n_phi, j_tr;
887 const BAS_FCTS *bas_fcts;
888
889 bas_fcts = qfast->bas_fcts;
890 n_phi = bas_fcts->n_trace_bas_fcts[wall];
891 phi = qfast->phi;
892
893 for (j_tr = 0; j_tr < n_phi; j_tr++) {
894 int j = bas_fcts->trace_dof_map[0][0][wall][j_tr];
895 val = 0.0;
896 for (iq = 0; iq < n_points[wall]; iq++) {
897 val += phi[iq][j]*wdetf_qp[iq];
898 }
899 fh->vec[dof->vec[j]] += val;
900 }
901
902 fh = CHAIN_NEXT(fh, DOF_REAL_VEC);
903 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
904 } CHAIN_WHILE(qfast, const QUAD_FAST);
905 }
906 } TRAVERSE_NEXT();
907
908 free_el_dof_vec(dof);
909 }
910
911 return not_all;
912 }
913
bndry_L2scp_fct_bas(DOF_REAL_VEC * fh,REAL (* f)(const REAL_D x,const REAL_D normal),const BNDRY_FLAGS scp_segment,const WALL_QUAD * quad)914 bool bndry_L2scp_fct_bas(DOF_REAL_VEC *fh,
915 REAL (*f)(const REAL_D x, const REAL_D normal),
916 const BNDRY_FLAGS scp_segment,
917 const WALL_QUAD *quad)
918 {
919 return
920 _AI_bndry_L2scp_fct_bas(fh, NULL, NULL, FILL_NOTHING, f, scp_segment, quad);
921 }
922
bndry_L2scp_fct_bas_loc(DOF_REAL_VEC * fh,LOC_FCT_AT_QP f,void * fd,FLAGS fill_flag,const BNDRY_FLAGS scp_segment,const WALL_QUAD * quad)923 bool bndry_L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
924 LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
925 const BNDRY_FLAGS scp_segment,
926 const WALL_QUAD *quad)
927 {
928 return
929 _AI_bndry_L2scp_fct_bas(fh, f, fd, fill_flag, NULL, scp_segment, quad);
930 }
931
932 /* Compute the L2 scalar product over the boundary of the domain. This
933 * version is for vector-valued fe-spaces.
934 */
935 static bool
_AI_bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D * fh,LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flag,const REAL * (* f)(const REAL_D x,const REAL_D normal,REAL_D result),const BNDRY_FLAGS bndry_seg,const WALL_QUAD * wall_quad)936 _AI_bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
937 LOC_FCT_D_AT_QP f_loc, void *fd, FLAGS fill_flag,
938 const REAL *(*f)(const REAL_D x,
939 const REAL_D normal,
940 REAL_D result),
941 const BNDRY_FLAGS bndry_seg,
942 const WALL_QUAD *wall_quad)
943 {
944 MESH *mesh = NULL;
945 const FE_SPACE *fe_space;
946 const BAS_FCTS *bas_fcts;
947 const WALL_QUAD_FAST *wqfast;
948 const QUAD_FAST *qfast; /* for a single wall */
949 REAL_D x, normal;
950 REAL det;
951 int iq, not_all = false;
952
953 if (!f && !f_loc) {
954 return not_all;
955 }
956
957 if (!fh) {
958 return not_all;
959 }
960
961 if (BNDRY_FLAGS_IS_INTERIOR(bndry_seg)) {
962 /* nothing to do */
963 return not_all;
964 }
965
966 fe_space = fh->fe_space;
967 bas_fcts = fe_space->bas_fcts;
968 mesh = fe_space->mesh;
969
970 if (wall_quad == NULL) {
971 wall_quad = get_wall_quad(mesh->dim, 2*bas_fcts->degree);
972 }
973 wqfast = get_wall_quad_fast(bas_fcts, wall_quad, INIT_PHI);
974
975 fill_flag |= CALL_LEAF_EL|FILL_COORDS|FILL_MACRO_WALLS;
976
977 if (mesh->is_periodic && !(fe_space->admin->flags & ADM_PERIODIC)) {
978 fill_flag |= FILL_NON_PERIODIC;
979 }
980 fill_flag |= wqfast->fill_flags;
981
982 {
983 PARAMETRIC *parametric = mesh->parametric;
984 int is_parametric = false;
985 int dim = mesh->dim;
986 REAL_D normal_p[wall_quad->n_points_max];
987 REAL det_p[wall_quad->n_points_max];
988 REAL_D wdetf_qp[wall_quad->n_points_max];
989 REAL_D x_p[wall_quad->n_points_max];
990 const REAL *w[N_WALLS_MAX];
991 const REAL_B *lambda[N_WALLS_MAX];
992 int n_points[N_WALLS_MAX];
993 INIT_EL_TAG wall_tag[N_WALLS_MAX];
994 int wall;
995 EL_DOF_VEC *dof;
996
997 for (wall = 0; wall < N_WALLS(dim); wall++) {
998 wall_tag[wall] = INIT_EL_TAG_DFLT;
999 lambda[wall] = wall_quad->quad[wall].lambda;
1000 w[wall] = wqfast->quad_fast[wall]->w;
1001 n_points[wall] = wall_quad->quad[wall].n_points;
1002 }
1003
1004 dof = get_el_dof_vec(bas_fcts);
1005
1006 TRAVERSE_FIRST(mesh, -1, fill_flag) {
1007 BNDRY_TYPE wall_bnd;
1008 int walls[N_WALLS_MAX];
1009 int n_walls, wnr, wall;
1010
1011 n_walls = 0;
1012 for (wnr = 0; wnr < N_WALLS(dim); wnr++) {
1013 wall_bnd = wall_bound(el_info, wnr);
1014 if (wall_bnd == INTERIOR) {
1015 continue;
1016 }
1017 if (bndry_seg == NULL ||
1018 BNDRY_FLAGS_IS_AT_BNDRY(bndry_seg, wall_bnd)) {
1019 walls[n_walls++] = wnr;
1020 } else {
1021 not_all = true;
1022 }
1023 }
1024 if (n_walls == 0) {
1025 continue;
1026 }
1027
1028 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
1029 continue;
1030 }
1031
1032 /* We assume that it is somewhat likely that we have to do some
1033 * work here and fetch the local DOF-indices here.
1034 */
1035 get_dof_indices(dof, fe_space, el_info->el);
1036
1037 if (parametric) {
1038 is_parametric = parametric->init_element(el_info, parametric);
1039 }
1040
1041 for (wnr = 0; wnr < n_walls; wnr++) {
1042
1043 wall = walls[wnr];
1044
1045 INIT_ELEMENT_SETUP(el_info,
1046 &wall_quad->quad[wall], wall_tag[wall], continue, {
1047 lambda[wall] = wall_quad->quad[wall].lambda;
1048 w[wall] = wall_quad->quad[wall].w;
1049 n_points[wall] = wall_quad->quad[wall].n_points;
1050 });
1051
1052 if (is_parametric) {
1053 if (f) {
1054 parametric->wall_normal(
1055 el_info, wall, &wall_quad->quad[wall], -1, NULL,
1056 normal_p, NULL, NULL, det_p);
1057 parametric->coord_to_world(el_info,
1058 &wall_quad->quad[wall], -1, NULL, x_p);
1059 for (iq = 0; iq < n_points[wall]; iq++) {
1060 AXEY_DOW(w[wall][iq]*det_p[iq],
1061 f(x_p[iq], normal_p[iq], wdetf_qp[iq]), wdetf_qp[iq]);
1062 }
1063 } else {
1064 const QUAD_EL_CACHE *qelc =
1065 fill_quad_el_cache(el_info, &wall_quad->quad[wall],
1066 FILL_EL_QUAD_WALL_DET|
1067 FILL_EL_QUAD_WALL_NORMAL);
1068 for (iq = 0; iq < n_points[wall]; iq++) {
1069 f_loc(wdetf_qp[iq], el_info, &wall_quad->quad[wall], iq, fd);
1070 SCAL_DOW(w[wall][iq]*qelc->param.wall_det[iq], wdetf_qp[iq]);
1071 }
1072 }
1073 } else {
1074 if (f) {
1075 det = get_wall_normal_dim(dim, el_info, wall, normal);
1076 for (iq = 0; iq < n_points[wall]; iq++) {
1077 coord_to_world(el_info, lambda[wall][iq], x);
1078 AXEY_DOW(w[wall][iq]*det, f(x, normal, wdetf_qp[iq]),
1079 wdetf_qp[iq]);
1080 }
1081 } else {
1082 const EL_GEOM_CACHE *elgc =
1083 fill_el_geom_cache(el_info,
1084 FILL_EL_WALL_DET(wall)|
1085 FILL_EL_WALL_NORMAL(wall));
1086 for (iq = 0; iq < n_points[wall]; iq++) {
1087 f_loc(wdetf_qp[iq], el_info, &wall_quad->quad[wall], iq, fd);
1088 SCAL_DOW(w[wall][iq]*elgc->wall_det[wall], wdetf_qp[iq]);
1089 }
1090 }
1091 }
1092
1093 /*****************************/
1094 qfast = wqfast->quad_fast[wall];
1095 INIT_ELEMENT(el_info, qfast); /* must be: tag != INIT_EL_TAG_NULL */
1096
1097 CHAIN_DO(qfast, const QUAD_FAST) {
1098 const BAS_FCTS *bas_fcts;
1099 const REAL_D *const*phi_d;
1100 int n_phi, j_tr;
1101
1102 bas_fcts = qfast->bas_fcts;
1103 n_phi = bas_fcts->n_trace_bas_fcts[wall];
1104
1105 if (fh->stride != 1) {
1106 DOF_REAL_D_VEC *fhd = (DOF_REAL_D_VEC *)fh;
1107
1108 for (j_tr = 0; j_tr < n_phi; j_tr++) {
1109 int j = bas_fcts->trace_dof_map[0][0][wall][j_tr];
1110 REAL_D val = { 0.0, };
1111
1112 for (iq = 0; iq < n_points[wall]; iq++) {
1113 AXPY_DOW(qfast->phi[iq][j], wdetf_qp[iq], val);
1114 }
1115 AXPY_DOW(1.0, val, fhd->vec[dof->vec[j]]);
1116 }
1117 } else {
1118 phi_d = get_quad_fast_phi_dow(qfast);
1119 for (j_tr = 0; j_tr < n_phi; j_tr++) {
1120 int j = bas_fcts->trace_dof_map[0][0][wall][j_tr];
1121 REAL val = 0.0;
1122
1123 for (iq = 0; iq < n_points[wall]; iq++) {
1124 val += SCP_DOW(phi_d[iq][j], wdetf_qp[iq]);
1125 }
1126 fh->vec[dof->vec[j]] += val;
1127 }
1128 }
1129
1130 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
1131 fh = CHAIN_NEXT(fh, DOF_REAL_VEC_D);
1132
1133 } CHAIN_WHILE(qfast, const QUAD_FAST);
1134
1135 }
1136 } TRAVERSE_NEXT();
1137
1138 free_el_dof_vec(dof);
1139 }
1140
1141 return not_all;
1142 }
1143
bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D * fh,const REAL * (* f)(const REAL_D x,const REAL_D normal,REAL_D result),const BNDRY_FLAGS bndry_seg,const WALL_QUAD * quad)1144 bool bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
1145 const REAL *(*f)(const REAL_D x,
1146 const REAL_D normal,
1147 REAL_D result),
1148 const BNDRY_FLAGS bndry_seg,
1149 const WALL_QUAD *quad)
1150 {
1151 return _AI_bndry_L2scp_fct_bas_dow(fh, NULL, NULL, FILL_NOTHING,
1152 f, bndry_seg, quad);
1153 }
1154
bndry_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D * fh,const REAL * (* f_loc)(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int iq,void * fd),void * fd,FLAGS fill_flag,const BNDRY_FLAGS bndry_seg,const WALL_QUAD * quad)1155 bool bndry_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
1156 const REAL *(*f_loc)(REAL_D result,
1157 const EL_INFO *el_info,
1158 const QUAD *quad,
1159 int iq,
1160 void *fd),
1161 void *fd, FLAGS fill_flag,
1162 const BNDRY_FLAGS bndry_seg,
1163 const WALL_QUAD *quad)
1164 {
1165 return _AI_bndry_L2scp_fct_bas_dow(fh, f_loc, fd, fill_flag,
1166 NULL, bndry_seg, quad);
1167 }
1168
1169 /* Compute the tangential H1 scalar product over the boundary of the
1170 * domain. Return "true" if not all boundary segments of the mesh
1171 * belong to the segment specified by "bndry_seg". If "bndry_seg ==
1172 * NULL" then the scalar product is computed over the entire boundary
1173 * (i.e. over all walls without neighbour).
1174 */
_AI_bndry_H1scp_fct_bas(DOF_REAL_VEC * fh,GRD_LOC_FCT_AT_QP f_loc,void * fd,FLAGS fill_flag,const REAL * (* f)(REAL_D result,const REAL_D x,const REAL_D normal),const BNDRY_FLAGS bndry_seg,const WALL_QUAD * wall_quad)1175 static bool _AI_bndry_H1scp_fct_bas(DOF_REAL_VEC *fh,
1176 GRD_LOC_FCT_AT_QP f_loc,
1177 void *fd, FLAGS fill_flag,
1178 const REAL *(*f)(REAL_D result,
1179 const REAL_D x,
1180 const REAL_D normal),
1181 const BNDRY_FLAGS bndry_seg,
1182 const WALL_QUAD *wall_quad)
1183 {
1184 MESH *mesh = NULL;
1185 const FE_SPACE *fe_space;
1186 const BAS_FCTS *bas_fcts;
1187 const WALL_QUAD_FAST *wqfast;
1188 const QUAD_FAST *qfast;
1189 REAL val;
1190 int iq, not_all = false;
1191
1192 if (!f && !f_loc) {
1193 return not_all;
1194 }
1195
1196 if (!fh) {
1197 return not_all;
1198 }
1199
1200 if (BNDRY_FLAGS_IS_INTERIOR(bndry_seg)) {
1201 /* nothing to do */
1202 return not_all;
1203 }
1204
1205 fe_space = fh->fe_space;
1206 bas_fcts = fe_space->bas_fcts;
1207 mesh = fe_space->mesh;
1208
1209 if (wall_quad == NULL) {
1210 wall_quad = get_wall_quad(mesh->dim, 2*bas_fcts->degree);
1211 }
1212 wqfast =
1213 get_wall_quad_fast(bas_fcts, wall_quad, INIT_GRD_PHI|INIT_TANGENTIAL);
1214
1215 fill_flag |= CALL_LEAF_EL|FILL_COORDS|FILL_MACRO_WALLS;
1216
1217 if (mesh->is_periodic && !(fe_space->admin->flags & ADM_PERIODIC)) {
1218 fill_flag |= FILL_NON_PERIODIC;
1219 }
1220 fill_flag |= wqfast->fill_flags;
1221
1222 {
1223 PARAMETRIC *parametric = mesh->parametric;
1224 int is_parametric = false;
1225 int dim = mesh->dim;
1226 REAL_B wdetgrdf_qp[wall_quad->n_points_max];
1227 const REAL *w[N_WALLS_MAX];
1228 int n_points[N_WALLS_MAX];
1229 INIT_EL_TAG wall_tag[N_WALLS_MAX];
1230 EL_DOF_VEC *dof;
1231 int wall;
1232
1233 for (wall = 0; wall < N_WALLS(dim); wall++) {
1234 wall_tag[wall] = INIT_EL_TAG_DFLT;
1235 w[wall] = wqfast->quad_fast[wall]->w;
1236 n_points[wall] = wall_quad->quad[wall].n_points;
1237 }
1238
1239 dof = get_el_dof_vec(bas_fcts);
1240
1241 TRAVERSE_FIRST(mesh, -1, fill_flag) {
1242 BNDRY_TYPE wall_bnd;
1243 int walls[N_WALLS_MAX];
1244 int n_walls, wnr, wall;
1245
1246 n_walls = 0;
1247 for (wnr = 0; wnr < N_WALLS(dim); wnr++) {
1248 wall_bnd = wall_bound(el_info, wnr);
1249 if (wall_bnd == INTERIOR) {
1250 continue;
1251 }
1252 if (bndry_seg == NULL ||
1253 BNDRY_FLAGS_IS_AT_BNDRY(bndry_seg, wall_bnd)) {
1254 walls[n_walls++] = wnr;
1255 } else {
1256 not_all = true;
1257 }
1258 }
1259 if (n_walls == 0) {
1260 continue;
1261 }
1262
1263 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
1264 continue;
1265 }
1266
1267 /* We assume that it is somewhat likely that we have to do some
1268 * work here and fetch the local DOF-indices here.
1269 */
1270 get_dof_indices(dof, fe_space, el_info->el);
1271
1272 if (parametric) {
1273 is_parametric = parametric->init_element(el_info, parametric);
1274 }
1275
1276 for (wnr = 0; wnr < n_walls; wnr++) {
1277
1278 wall = walls[wnr];
1279
1280 INIT_ELEMENT_SETUP(el_info,
1281 wqfast->quad_fast[wall], wall_tag[wall], continue, {
1282 w[wall] = wqfast->quad_fast[wall]->w;
1283 n_points[wall] = wall_quad->quad[wall].n_points;
1284 });
1285
1286 if (is_parametric) {
1287 const QUAD_EL_CACHE *qelc;
1288 if (f) {
1289 qelc = fill_quad_el_cache(el_info, &wall_quad->quad[wall],
1290 FILL_EL_QUAD_LAMBDA|
1291 FILL_EL_QUAD_WALL_NORMAL|
1292 FILL_EL_QUAD_WALL_DET|
1293 FILL_EL_QUAD_WORLD);
1294 for (iq = 0; iq < n_points[wall]; iq++) {
1295 REAL_D grd;
1296 int alpha;
1297
1298 f(grd, qelc->world[iq], qelc->param.wall_normal[iq]);
1299 /* Project to the tangent space */
1300 AXPY_DOW(-SCP_DOW(grd, qelc->param.wall_normal[iq]),
1301 qelc->param.wall_normal[iq], grd);
1302 SCAL_DOW(w[wall][iq]*qelc->param.wall_det[iq], grd);
1303 for (alpha = 0; alpha < wall; alpha++) {
1304 wdetgrdf_qp[iq][alpha] =
1305 SCP_DOW(grd, qelc->param.Lambda[iq][alpha]);
1306 }
1307 wdetgrdf_qp[iq][alpha++] = 0.0;
1308 for (; alpha < N_LAMBDA(dim); alpha++) {
1309 wdetgrdf_qp[iq][alpha] =
1310 SCP_DOW(grd, qelc->param.Lambda[iq][alpha]);
1311 }
1312 }
1313 } else {
1314 qelc = fill_quad_el_cache(el_info, &wall_quad->quad[wall],
1315 FILL_EL_QUAD_LAMBDA|
1316 FILL_EL_QUAD_WALL_NORMAL|
1317 FILL_EL_QUAD_WALL_DET);
1318 for (iq = 0; iq < n_points[wall]; iq++) {
1319 REAL_D grd;
1320 int alpha;
1321
1322 f_loc(grd, el_info, (const REAL_D *)qelc->param.Lambda[iq],
1323 &wall_quad->quad[wall], iq, fd);
1324 /* Project to the tangent space */
1325 AXPY_DOW(-SCP_DOW(grd, qelc->param.wall_normal[iq]),
1326 qelc->param.wall_normal[iq], grd);
1327 SCAL_DOW(w[wall][iq] * qelc->param.wall_det[iq], grd);
1328 for (alpha = 0; alpha < wall; alpha++) {
1329 wdetgrdf_qp[iq][alpha] =
1330 SCP_DOW(grd, qelc->param.Lambda[iq][alpha]);
1331 }
1332 wdetgrdf_qp[iq][alpha++] = 0.0;
1333 for (; alpha < N_LAMBDA(dim); alpha++) {
1334 wdetgrdf_qp[iq][alpha] =
1335 SCP_DOW(grd, qelc->param.Lambda[iq][alpha]);
1336 }
1337 }
1338 }
1339 } else { /* non-parametric version */
1340 const EL_GEOM_CACHE *elgc =
1341 fill_el_geom_cache(el_info,
1342 FILL_EL_WALL_NORMAL(wall)|
1343 FILL_EL_WALL_DET(wall)|
1344 FILL_EL_LAMBDA);
1345 if (f) {
1346 const QUAD_EL_CACHE *qelc =
1347 fill_quad_el_cache(el_info,
1348 &wall_quad->quad[wall], FILL_EL_QUAD_WORLD);
1349 for (iq = 0; iq < n_points[wall]; iq++) {
1350 REAL_D grd;
1351 int alpha;
1352
1353 f(grd, qelc->world[iq], elgc->wall_normal[wall]);
1354 /* Project to the tangent space */
1355 AXPY_DOW(-SCP_DOW(grd, elgc->wall_normal[wall]),
1356 elgc->wall_normal[wall], grd);
1357 SCAL_DOW(w[wall][iq]*elgc->wall_det[wall], grd);
1358 for (alpha = 0; alpha < wall; alpha++) {
1359 wdetgrdf_qp[iq][alpha] = SCP_DOW(grd, elgc->Lambda[alpha]);
1360 }
1361 wdetgrdf_qp[iq][alpha++] = 0.0;
1362 for (; alpha < N_LAMBDA(dim); alpha++) {
1363 wdetgrdf_qp[iq][alpha] = SCP_DOW(grd, elgc->Lambda[alpha]);
1364 }
1365 }
1366 } else {
1367 for (iq = 0; iq < n_points[wall]; iq++) {
1368 REAL_D grd;
1369 int alpha;
1370
1371 f_loc(grd, el_info, elgc->Lambda, &wall_quad->quad[wall], iq, fd);
1372 /* Project to the tangent space */
1373 AXPY_DOW(-SCP_DOW(grd, elgc->wall_normal[wall]),
1374 elgc->wall_normal[wall], grd);
1375 SCAL_DOW(w[wall][iq]*elgc->wall_det[wall], grd);
1376 for (alpha = 0; alpha < wall; alpha++) {
1377 wdetgrdf_qp[iq][alpha] = SCP_DOW(grd, elgc->Lambda[alpha]);
1378 }
1379 wdetgrdf_qp[iq][alpha++] = 0.0;
1380 for (; alpha < N_LAMBDA(dim); alpha++) {
1381 wdetgrdf_qp[iq][alpha] = SCP_DOW(grd, elgc->Lambda[alpha]);
1382 }
1383 }
1384 }
1385 }
1386
1387 /*****************************/
1388 qfast = wqfast->quad_fast[wall];
1389 INIT_ELEMENT(el_info, qfast); /* must be: tag != INIT_EL_TAG_NULL */
1390
1391 CHAIN_DO(qfast, const QUAD_FAST) {
1392 int n_phi, j_tr;
1393 const BAS_FCTS *bas_fcts;
1394 const REAL_B *const*grd_phi = qfast->grd_phi;
1395
1396 bas_fcts = qfast->bas_fcts;
1397 n_phi = bas_fcts->n_trace_bas_fcts[wall];
1398
1399 for (j_tr = 0; j_tr < n_phi; j_tr++) {
1400 int j = bas_fcts->trace_dof_map[0][0][wall][j_tr];
1401 val = 0.0;
1402 for (iq = 0; iq < n_points[wall]; iq++) {
1403 val += SCP_BAR(dim, grd_phi[iq][j], wdetgrdf_qp[iq]);
1404 }
1405 fh->vec[dof->vec[j]] += val;
1406 }
1407
1408 fh = CHAIN_NEXT(fh, DOF_REAL_VEC);
1409 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
1410 } CHAIN_WHILE(qfast, const QUAD_FAST);
1411 }
1412 } TRAVERSE_NEXT();
1413
1414 free_el_dof_vec(dof);
1415 }
1416
1417 return not_all;
1418 }
1419
bndry_H1scp_fct_bas(DOF_REAL_VEC * fh,const REAL * (* f)(REAL_D result,const REAL_D x,const REAL_D normal),const BNDRY_FLAGS scp_segment,const WALL_QUAD * quad)1420 bool bndry_H1scp_fct_bas(DOF_REAL_VEC *fh,
1421 const REAL *(*f)(REAL_D result,
1422 const REAL_D x, const REAL_D normal),
1423 const BNDRY_FLAGS scp_segment,
1424 const WALL_QUAD *quad)
1425 {
1426 return
1427 _AI_bndry_H1scp_fct_bas(fh, NULL, NULL, FILL_NOTHING, f, scp_segment, quad);
1428 }
1429
bndry_H1scp_fct_bas_loc(DOF_REAL_VEC * fh,GRD_LOC_FCT_AT_QP f,void * fd,FLAGS fill_flag,const BNDRY_FLAGS scp_segment,const WALL_QUAD * quad)1430 bool bndry_H1scp_fct_bas_loc(DOF_REAL_VEC *fh,
1431 GRD_LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
1432 const BNDRY_FLAGS scp_segment,
1433 const WALL_QUAD *quad)
1434 {
1435 return
1436 _AI_bndry_H1scp_fct_bas(fh, f, fd, fill_flag, NULL, scp_segment, quad);
1437 }
1438
1439 /* Compute the L2 scalar product over the boundary of the domain. This
1440 * version is for vector-valued fe-spaces.
1441 */
1442 static bool
_AI_bndry_H1scp_fct_bas_dow(DOF_REAL_VEC_D * fh,GRD_LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flag,const REAL_D * (* f)(REAL_DD result,const REAL_D x,const REAL_D normal),const BNDRY_FLAGS bndry_seg,const WALL_QUAD * wall_quad)1443 _AI_bndry_H1scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
1444 GRD_LOC_FCT_D_AT_QP f_loc,
1445 void *fd, FLAGS fill_flag,
1446 const REAL_D *(*f)(REAL_DD result,
1447 const REAL_D x,
1448 const REAL_D normal),
1449 const BNDRY_FLAGS bndry_seg,
1450 const WALL_QUAD *wall_quad)
1451 {
1452 MESH *mesh = NULL;
1453 const FE_SPACE *fe_space;
1454 const BAS_FCTS *bas_fcts;
1455 const WALL_QUAD_FAST *wqfast;
1456 const QUAD_FAST *qfast; /* for a single wall */
1457 int iq, not_all = false;
1458
1459 if (!f && !f_loc) {
1460 return not_all;
1461 }
1462
1463 if (!fh) {
1464 return not_all;
1465 }
1466
1467 if (BNDRY_FLAGS_IS_INTERIOR(bndry_seg)) {
1468 /* nothing to do */
1469 return not_all;
1470 }
1471
1472 fe_space = fh->fe_space;
1473 bas_fcts = fe_space->bas_fcts;
1474 mesh = fe_space->mesh;
1475
1476 if (wall_quad == NULL) {
1477 wall_quad = get_wall_quad(mesh->dim, 2*bas_fcts->degree);
1478 }
1479 wqfast =
1480 get_wall_quad_fast(bas_fcts, wall_quad, INIT_GRD_PHI|INIT_TANGENTIAL);
1481
1482 fill_flag |= CALL_LEAF_EL|FILL_COORDS|FILL_MACRO_WALLS;
1483
1484 if (mesh->is_periodic && !(fe_space->admin->flags & ADM_PERIODIC)) {
1485 fill_flag |= FILL_NON_PERIODIC;
1486 }
1487 fill_flag |= wqfast->fill_flags;
1488
1489 {
1490 PARAMETRIC *parametric = mesh->parametric;
1491 int is_parametric = false;
1492 int dim = mesh->dim;
1493 REAL_D normal_p[wall_quad->n_points_max];
1494 REAL_DB wdetgrdf_qp[wall_quad->n_points_max];
1495 REAL_D x_p[wall_quad->n_points_max];
1496 const REAL *w[N_WALLS_MAX];
1497 int n_points[N_WALLS_MAX];
1498 INIT_EL_TAG wall_tag[N_WALLS_MAX];
1499 int wall;
1500 EL_DOF_VEC *dof;
1501
1502 for (wall = 0; wall < N_WALLS(dim); wall++) {
1503 wall_tag[wall] = INIT_EL_TAG_DFLT;
1504 w[wall] = wqfast->quad_fast[wall]->w;
1505 n_points[wall] = wall_quad->quad[wall].n_points;
1506 }
1507
1508 dof = get_el_dof_vec(bas_fcts);
1509
1510 TRAVERSE_FIRST(mesh, -1, fill_flag) {
1511 BNDRY_TYPE wall_bnd;
1512 int walls[N_WALLS_MAX];
1513 int n_walls, wnr, wall;
1514
1515 n_walls = 0;
1516 for (wnr = 0; wnr < N_WALLS(dim); wnr++) {
1517 wall_bnd = wall_bound(el_info, wnr);
1518 if (wall_bnd == INTERIOR) {
1519 continue;
1520 }
1521 if (bndry_seg == NULL ||
1522 BNDRY_FLAGS_IS_AT_BNDRY(bndry_seg, wall_bnd)) {
1523 walls[n_walls++] = wnr;
1524 } else {
1525 not_all = true;
1526 }
1527 }
1528 if (n_walls == 0) {
1529 continue;
1530 }
1531
1532 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
1533 continue;
1534 }
1535
1536 /* We assume that it is somewhat likely that we have to do some
1537 * work here and fetch the local DOF-indices here.
1538 */
1539 get_dof_indices(dof, fe_space, el_info->el);
1540
1541 if (parametric) {
1542 is_parametric = parametric->init_element(el_info, parametric);
1543 }
1544
1545 for (wnr = 0; wnr < n_walls; wnr++) {
1546
1547 wall = walls[wnr];
1548
1549 INIT_ELEMENT_SETUP(el_info,
1550 &wall_quad->quad[wall], wall_tag[wall], continue, {
1551 w[wall] = wall_quad->quad[wall].w;
1552 n_points[wall] = wall_quad->quad[wall].n_points;
1553 });
1554
1555 if (is_parametric) {
1556 const QUAD_EL_CACHE *qelc =
1557 fill_quad_el_cache(el_info, &wall_quad->quad[wall],
1558 FILL_EL_QUAD_LAMBDA|
1559 FILL_EL_QUAD_WALL_DET|
1560 FILL_EL_QUAD_WALL_NORMAL);
1561 if (f) {
1562 fill_quad_el_cache(el_info, &wall_quad->quad[wall],
1563 FILL_EL_QUAD_WORLD);
1564 for (iq = 0; iq < n_points[wall]; iq++) {
1565 REAL_DD grd;
1566 int alpha, i;
1567
1568 f(grd, x_p[iq], normal_p[iq]);
1569 /* Project to the tangent space */
1570 for (i = 0; i < DIM_OF_WORLD; i++) {
1571 AXPY_DOW(-SCP_DOW(grd[i], qelc->param.wall_normal[iq]),
1572 qelc->param.wall_normal[iq], grd[i]);
1573 }
1574 MSCAL_DOW(w[wall][iq]*qelc->param.wall_det[iq], grd);
1575 for (alpha = 0; alpha < wall; alpha++) {
1576 for (i = 0; i < DIM_OF_WORLD; i++) {
1577 wdetgrdf_qp[iq][i][alpha] =
1578 SCP_DOW(grd[i], qelc->param.Lambda[iq][alpha]);
1579 }
1580 }
1581 for (i = 0; i < DIM_OF_WORLD; i++) {
1582 wdetgrdf_qp[iq][i][alpha] = 0.0;
1583 }
1584 for (++alpha; alpha < N_LAMBDA(dim); alpha++) {
1585 for (i = 0; i < DIM_OF_WORLD; i++) {
1586 wdetgrdf_qp[iq][i][alpha] =
1587 SCP_DOW(grd[i], qelc->param.Lambda[iq][alpha]);
1588 }
1589 }
1590 }
1591 } else {
1592 for (iq = 0; iq < n_points[wall]; iq++) {
1593 REAL_DD grd;
1594 int alpha, i;
1595
1596 f_loc(grd, el_info, (const REAL_D *)qelc->param.Lambda[iq],
1597 &wall_quad->quad[wall], iq, fd);
1598 /* Project to the tangent space */
1599 for (i = 0; i < DIM_OF_WORLD; i++) {
1600 AXPY_DOW(-SCP_DOW(grd[i], qelc->param.wall_normal[iq]),
1601 qelc->param.wall_normal[iq], grd[i]);
1602 }
1603 MSCAL_DOW(w[wall][iq]*qelc->param.wall_det[iq], grd);
1604 for (alpha = 0; alpha < wall; alpha++) {
1605 for (i = 0; i < DIM_OF_WORLD; i++) {
1606 wdetgrdf_qp[iq][i][alpha] =
1607 SCP_DOW(grd[i], qelc->param.Lambda[iq][alpha]);
1608 }
1609 }
1610 for (i = 0; i < DIM_OF_WORLD; i++) {
1611 wdetgrdf_qp[iq][i][alpha] = 0.0;
1612 }
1613 for (++alpha; alpha < N_LAMBDA(dim); alpha++) {
1614 for (i = 0; i < DIM_OF_WORLD; i++) {
1615 wdetgrdf_qp[iq][i][alpha] =
1616 SCP_DOW(grd[i], qelc->param.Lambda[iq][alpha]);
1617 }
1618 }
1619 }
1620 }
1621 } else {
1622 const EL_GEOM_CACHE *elgc =
1623 fill_el_geom_cache(el_info,
1624 FILL_EL_LAMBDA|
1625 FILL_EL_WALL_DET(wall)|
1626 FILL_EL_WALL_NORMAL(wall));
1627 if (f) {
1628 const QUAD_EL_CACHE *qelc =
1629 fill_quad_el_cache(el_info, &wall_quad->quad[wall],
1630 FILL_EL_QUAD_WORLD);
1631 for (iq = 0; iq < n_points[wall]; iq++) {
1632 REAL_DD grd;
1633 int alpha, i;
1634
1635 f(grd, qelc->world[iq], elgc->wall_normal[wall]);
1636 /* Project to the tangent space */
1637 for (i = 0; i < DIM_OF_WORLD; i++) {
1638 AXPY_DOW(-SCP_DOW(grd[i], elgc->wall_normal[wall]),
1639 elgc->wall_normal[wall], grd[i]);
1640 }
1641 MSCAL_DOW(w[wall][iq]*elgc->wall_det[wall], grd);
1642 for (alpha = 0; alpha < wall; alpha++) {
1643 for (i = 0; i < DIM_OF_WORLD; i++) {
1644 wdetgrdf_qp[iq][i][alpha] =
1645 SCP_DOW(grd[i], elgc->Lambda[alpha]);
1646 }
1647 }
1648 for (i = 0; i < DIM_OF_WORLD; i++) {
1649 wdetgrdf_qp[iq][i][alpha] = 0.0;
1650 }
1651 for (++alpha; alpha < N_LAMBDA(dim); alpha++) {
1652 for (i = 0; i < DIM_OF_WORLD; i++) {
1653 wdetgrdf_qp[iq][i][alpha] =
1654 SCP_DOW(grd[i], elgc->Lambda[alpha]);
1655 }
1656 }
1657 }
1658 } else {
1659 for (iq = 0; iq < n_points[wall]; iq++) {
1660 REAL_DD grd;
1661 int alpha, i;
1662
1663 f_loc(grd, el_info, elgc->Lambda, &wall_quad->quad[wall], iq, fd);
1664 /* Project to the tangent space */
1665 for (i = 0; i < DIM_OF_WORLD; i++) {
1666 AXPY_DOW(-SCP_DOW(grd[i], elgc->wall_normal[wall]),
1667 elgc->wall_normal[wall], grd[i]);
1668 }
1669 MSCAL_DOW(w[wall][iq]*elgc->wall_det[wall], grd);
1670 for (alpha = 0; alpha < wall; alpha++) {
1671 for (i = 0; i < DIM_OF_WORLD; i++) {
1672 wdetgrdf_qp[iq][i][alpha] =
1673 SCP_DOW(grd[i], elgc->Lambda[alpha]);
1674 }
1675 }
1676 for (i = 0; i < DIM_OF_WORLD; i++) {
1677 wdetgrdf_qp[iq][i][alpha] = 0.0;
1678 }
1679 for (++alpha; alpha < N_LAMBDA(dim); alpha++) {
1680 for (i = 0; i < DIM_OF_WORLD; i++) {
1681 wdetgrdf_qp[iq][i][alpha] =
1682 SCP_DOW(grd[i], elgc->Lambda[alpha]);
1683 }
1684 }
1685 }
1686 }
1687 }
1688
1689 /*****************************/
1690 qfast = wqfast->quad_fast[wall];
1691 INIT_ELEMENT(el_info, qfast); /* must be: tag != INIT_EL_TAG_NULL */
1692
1693 CHAIN_DO(qfast, const QUAD_FAST) {
1694 const BAS_FCTS *bas_fcts;
1695 int n_phi, j_tr;
1696
1697 bas_fcts = qfast->bas_fcts;
1698 n_phi = bas_fcts->n_trace_bas_fcts[wall];
1699
1700 if (fh->stride != 1) {
1701 DOF_REAL_D_VEC *fhd = (DOF_REAL_D_VEC *)fh;
1702 const REAL_B *const*grd_phi = qfast->grd_phi;
1703
1704 for (j_tr = 0; j_tr < n_phi; j_tr++) {
1705 int j = bas_fcts->trace_dof_map[0][0][wall][j_tr];
1706 REAL_D val = { 0.0, };
1707 int i;
1708
1709 for (iq = 0; iq < n_points[wall]; iq++) {
1710 for (i = 0; i < DIM_OF_WORLD; i++) {
1711 val[i] += SCP_BAR(dim, grd_phi[iq][j], wdetgrdf_qp[iq][i]);
1712 }
1713 }
1714 AXPY_DOW(1.0, val, fhd->vec[dof->vec[j]]);
1715 }
1716 } else {
1717 const REAL_DB *const*grd_phi;
1718
1719 grd_phi = get_quad_fast_grd_phi_dow(qfast);
1720
1721 for (j_tr = 0; j_tr < n_phi; j_tr++) {
1722 int j = bas_fcts->trace_dof_map[0][0][wall][j_tr];
1723 REAL val = 0.0;
1724 int i;
1725
1726 for (iq = 0; iq < n_points[wall]; iq++) {
1727 for (i = 0; i < DIM_OF_WORLD; i++) {
1728 val += SCP_BAR(dim, wdetgrdf_qp[iq][i], grd_phi[iq][j][i]);
1729 }
1730 }
1731 fh->vec[dof->vec[j]] += val;
1732 }
1733 }
1734
1735 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
1736 fh = CHAIN_NEXT(fh, DOF_REAL_VEC_D);
1737
1738 } CHAIN_WHILE(qfast, const QUAD_FAST);
1739
1740 }
1741 } TRAVERSE_NEXT();
1742
1743 free_el_dof_vec(dof);
1744 }
1745
1746 return not_all;
1747 }
1748
bndry_H1scp_fct_bas_dow(DOF_REAL_VEC_D * fh,const REAL_D * (* f)(REAL_DD result,const REAL_D x,const REAL_D normal),const BNDRY_FLAGS bndry_seg,const WALL_QUAD * quad)1749 bool bndry_H1scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
1750 const REAL_D *(*f)(REAL_DD result,
1751 const REAL_D x,
1752 const REAL_D normal),
1753 const BNDRY_FLAGS bndry_seg,
1754 const WALL_QUAD *quad)
1755 {
1756 return _AI_bndry_H1scp_fct_bas_dow(fh, NULL, NULL, FILL_NOTHING,
1757 f, bndry_seg, quad);
1758 }
1759
bndry_H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D * fh,GRD_LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flag,const BNDRY_FLAGS bndry_seg,const WALL_QUAD * quad)1760 bool bndry_H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
1761 GRD_LOC_FCT_D_AT_QP f_loc,
1762 void *fd, FLAGS fill_flag,
1763 const BNDRY_FLAGS bndry_seg,
1764 const WALL_QUAD *quad)
1765 {
1766 return _AI_bndry_H1scp_fct_bas_dow(fh, f_loc, fd, fill_flag,
1767 NULL, bndry_seg, quad);
1768 }
1769
1770 /* Compute the L2-scalar-product over the given trace-mesh. fh must
1771 * belong to the master mesh. quad must belong to the trace mesh.
1772 */
_AI_trace_L2scp_fct_bas(DOF_REAL_VEC * fh,LOC_FCT_AT_QP f_loc,void * fd,FLAGS fill_flag,FCT_AT_X f,MESH * mesh,const QUAD * quad)1773 static void _AI_trace_L2scp_fct_bas(DOF_REAL_VEC *fh,
1774 LOC_FCT_AT_QP f_loc,
1775 void *fd, FLAGS fill_flag,
1776 FCT_AT_X f,
1777 MESH *mesh,
1778 const QUAD *quad)
1779 {
1780 FUNCNAME("trace_L2scp_fct_bas");
1781 const QUAD_FAST *quad_fast;
1782 const FE_SPACE *fe_space;
1783 const REAL *w;
1784 REAL val;
1785 int iq, j, n_phi, dim;
1786 const BAS_FCTS *bas_fcts;
1787 INIT_EL_TAG tag = INIT_EL_TAG_DFLT;
1788
1789 TEST_EXIT(fh, "no DOF_REAL_VEC fh\n");
1790 if (!f && !f_loc)
1791 return;
1792
1793 TEST_EXIT(fh->fe_space, "no fe_space in DOF_REAL_VEC %s\n", NAME(fh));
1794 TEST_EXIT(mesh != NULL && fh->fe_space->mesh == get_master(mesh),
1795 "mesh is not a trace mesh of fh->fe_space->mesh\n");
1796
1797 dim = mesh->dim;
1798
1799 fe_space = fh->fe_space;
1800 bas_fcts = fe_space->bas_fcts->trace_bas_fcts;
1801
1802 if (!quad) {
1803 quad = get_quadrature(dim, 2*bas_fcts->degree);
1804 }
1805 quad_fast = get_quad_fast(bas_fcts, quad, INIT_PHI);
1806
1807 w = quad_fast->w;
1808 n_phi = bas_fcts->n_bas_fcts;
1809
1810 fill_flag |= CALL_LEAF_EL|FILL_COORDS|FILL_MASTER_INFO;
1811 if (mesh->is_periodic && !(fh->fe_space->admin->flags & ADM_PERIODIC)) {
1812 fill_flag |= FILL_NON_PERIODIC;
1813 }
1814 fill_flag |= quad_fast->fill_flags;
1815
1816 {
1817 PARAMETRIC *parametric = mesh->parametric;
1818 int is_parametric = false;
1819 REAL wdetf_qp[quad->n_points_max];
1820 EL_DOF_VEC *dof = get_el_dof_vec(fe_space->bas_fcts);
1821
1822 TRAVERSE_FIRST(mesh, -1, fill_flag) {
1823
1824 INIT_ELEMENT_SETUP(el_info, quad_fast, tag, continue, {
1825 w = quad_fast->w;
1826 n_phi = bas_fcts->n_bas_fcts;
1827 });
1828
1829 if (parametric) {
1830 is_parametric = parametric->init_element(el_info, parametric);
1831 }
1832
1833 if (is_parametric) {
1834 const QUAD_EL_CACHE *qelc =
1835 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_DET);
1836
1837 if (f) {
1838 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
1839 for (iq = 0; iq < quad->n_points; iq++) {
1840 wdetf_qp[iq] = w[iq]*qelc->param.det[iq]*f(qelc->world[iq]);
1841 }
1842 } else {
1843 for (iq = 0; iq < quad->n_points; iq++) {
1844 wdetf_qp[iq] =
1845 w[iq]*qelc->param.det[iq]*f_loc(el_info, quad, iq, fd);
1846 }
1847 }
1848 } else {
1849 const EL_GEOM_CACHE *elgc = fill_el_geom_cache(el_info, FILL_EL_DET);
1850 if (f) {
1851 const QUAD_EL_CACHE *qelc =
1852 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
1853 for (iq = 0; iq < quad->n_points; iq++) {
1854 wdetf_qp[iq] = w[iq]*elgc->det*f(qelc->world[iq]);
1855 }
1856 } else {
1857 for (iq = 0; iq < quad->n_points; iq++) {
1858 wdetf_qp[iq] = w[iq]*elgc->det*f_loc(el_info, quad, iq, fd);
1859 }
1860 }
1861 }
1862
1863 INIT_ELEMENT(el_info, quad_fast);
1864
1865 /* The following line is the only actual difference to L2scp_fct_bas: */
1866 get_master_dof_indices(dof, el_info, fe_space);
1867
1868 CHAIN_DO(quad_fast, const QUAD_FAST) {
1869 const BAS_FCTS *bas_fcts;
1870 const REAL *const*phi;
1871
1872 bas_fcts = quad_fast->bas_fcts;
1873 phi = quad_fast->phi;
1874 n_phi = bas_fcts->n_bas_fcts;
1875
1876 for (j = 0; j < n_phi; j++) {
1877 val = 0.0;
1878 for (iq = 0; iq < quad->n_points; iq++) {
1879 val += phi[iq][j]*wdetf_qp[iq];
1880 }
1881 fh->vec[dof->vec[j]] += val;
1882 }
1883 fh = CHAIN_NEXT(fh, DOF_REAL_VEC);
1884 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
1885 } CHAIN_WHILE(quad_fast, const QUAD_FAST);
1886 } TRAVERSE_NEXT();
1887
1888 free_el_dof_vec(dof);
1889 }
1890 }
1891
1892 __FLATTEN_ATTRIBUTE__
trace_L2scp_fct_bas(DOF_REAL_VEC * fh,FCT_AT_X f,MESH * trace_mesh,const QUAD * quad)1893 void trace_L2scp_fct_bas(DOF_REAL_VEC *fh, FCT_AT_X f,
1894 MESH *trace_mesh, const QUAD *quad)
1895 {
1896 _AI_trace_L2scp_fct_bas(fh, NULL, NULL, FILL_NOTHING, f, trace_mesh, quad);
1897 }
1898
1899 __FLATTEN_ATTRIBUTE__
trace_L2scp_fct_bas_loc(DOF_REAL_VEC * fh,LOC_FCT_AT_QP f,void * fd,FLAGS fill_flag,MESH * trace_mesh,const QUAD * quad)1900 void trace_L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
1901 LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
1902 MESH *trace_mesh,
1903 const QUAD *quad)
1904 {
1905 _AI_trace_L2scp_fct_bas(fh, f, fd, fill_flag, NULL, trace_mesh, quad);
1906 }
1907
_AI_trace_L2scp_fct_bas_dow(DOF_REAL_VEC_D * fh,FCT_D_AT_X f,LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flag,MESH * mesh,const QUAD * quad)1908 static void _AI_trace_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
1909 FCT_D_AT_X f,
1910 LOC_FCT_D_AT_QP f_loc,
1911 void *fd, FLAGS fill_flag,
1912 MESH *mesh,
1913 const QUAD *quad)
1914 {
1915 FUNCNAME("trace_L2scp_fct_bas_dow");
1916 const QUAD_FAST *quad_fast;
1917 const FE_SPACE *fe_space;
1918 const REAL *w;
1919 const REAL_B *lambda;
1920 REAL_D val;
1921 int iq, j, n_phi, dim, n_points;
1922 const BAS_FCTS *bas_fcts;
1923 INIT_EL_TAG qtag = INIT_EL_TAG_DFLT;
1924
1925 TEST_EXIT(fh, "no DOF_REAL_VEC fh\n");
1926 if (!f && !f_loc)
1927 return;
1928
1929 TEST_EXIT(fh->fe_space, "no fe_space in DOF_REAL_D_VEC \"%s\"\n", NAME(fh));
1930 TEST_EXIT(mesh != NULL && fh->fe_space->mesh == get_master(mesh),
1931 "mesh is not a trace mesh of fh->fe_space->mesh\n");
1932 TEST_EXIT(fh->fe_space->rdim == DIM_OF_WORLD,
1933 "Called for scalar finite element space \"%s\".\n",
1934 NAME(fh->fe_space));
1935
1936 dim = mesh->dim;
1937
1938 fe_space = fh->fe_space;
1939 bas_fcts = fe_space->bas_fcts->trace_bas_fcts;
1940
1941 if (!quad) {
1942 quad = get_quadrature(dim, 2*bas_fcts->degree-2);
1943 }
1944 quad_fast = get_quad_fast(bas_fcts, quad, INIT_PHI);
1945
1946 w = quad->w;
1947 lambda = quad->lambda;
1948 n_points = quad->n_points;
1949
1950 fill_flag |= CALL_LEAF_EL|FILL_COORDS|FILL_MASTER_INFO;
1951 if (mesh->is_periodic && !(fh->fe_space->admin->flags & ADM_PERIODIC)) {
1952 fill_flag |= FILL_NON_PERIODIC;
1953 }
1954 fill_flag |= quad_fast->fill_flags;
1955
1956 {
1957 PARAMETRIC *parametric = mesh->parametric;
1958 int is_parametric = false;
1959 REAL_D wdetf_qp[quad->n_points_max];
1960 EL_DOF_VEC *dof = get_el_dof_vec(fe_space->bas_fcts);
1961
1962 TRAVERSE_FIRST(mesh, -1, fill_flag) {
1963
1964 INIT_ELEMENT_SETUP(el_info, quad, qtag, continue, {
1965 w = quad->w;
1966 n_points = quad->n_points;
1967 lambda = quad->lambda;
1968 });
1969
1970 if (parametric) {
1971 is_parametric = parametric->init_element(el_info, parametric);
1972 }
1973
1974 if (is_parametric) {
1975 const QUAD_EL_CACHE *qelc =
1976 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_DET);
1977
1978 if (f) {
1979 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
1980 for (iq = 0; iq < quad->n_points; iq++) {
1981 AXEY_DOW(w[iq]*qelc->param.det[iq],
1982 f(qelc->world[iq], wdetf_qp[iq]),
1983 wdetf_qp[iq]);
1984 }
1985 } else {
1986 for (iq = 0; iq < quad->n_points; iq++) {
1987 AXEY_DOW(w[iq]*qelc->param.det[iq],
1988 f_loc(wdetf_qp[iq], el_info, quad, iq, fd),
1989 wdetf_qp[iq]);
1990 }
1991 }
1992 } else {
1993 const EL_GEOM_CACHE *elgc = fill_el_geom_cache(el_info, FILL_EL_DET);
1994
1995 if (f) {
1996 const QUAD_EL_CACHE *qelc =
1997 fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
1998 for (iq = 0; iq < quad->n_points; iq++) {
1999 AXEY_DOW(w[iq]*elgc->det,
2000 f(qelc->world[iq], wdetf_qp[iq]),
2001 wdetf_qp[iq]);
2002 }
2003 } else {
2004 for (iq = 0; iq < quad->n_points; iq++) {
2005 AXEY_DOW(w[iq]*elgc->det,
2006 f_loc(wdetf_qp[iq], el_info, quad, iq, fd),
2007 wdetf_qp[iq]);
2008 }
2009 }
2010 }
2011
2012 INIT_ELEMENT(el_info, quad_fast);
2013
2014 /* The following line is the only actual difference to L2scp_fct_bas: */
2015 get_master_dof_indices(dof, el_info, fe_space);
2016
2017 CHAIN_DO(quad_fast, const QUAD_FAST) {
2018 const BAS_FCTS *bas_fcts;
2019
2020 bas_fcts = quad_fast->bas_fcts;
2021 n_phi = bas_fcts->n_bas_fcts;
2022
2023 if (fh->stride != 1) {
2024 /* Cartesian product space */
2025 DOF_REAL_D_VEC *fhd = (DOF_REAL_D_VEC *)fh;
2026
2027 for (j = 0; j < n_phi; j++) {
2028 REAL_D val = { 0.0, };
2029 for (iq = 0; iq < quad->n_points; iq++) {
2030 AXPY_DOW(quad_fast->phi[iq][j], wdetf_qp[iq], val);
2031 }
2032 AXPY_DOW(1.0, val, fhd->vec[dof->vec[j]]);
2033 }
2034 } else {
2035 if (bas_fcts->dir_pw_const) {
2036 for (j = 0; j < n_phi; j++) {
2037 SET_DOW(0.0, val);
2038 for (iq = 0; iq < n_points; iq++) {
2039 AXPY_DOW(quad_fast->phi[iq][j], wdetf_qp[iq], val);
2040 }
2041 fh->vec[dof->vec[j]] +=
2042 SCP_DOW(val, PHI_D(quad_fast->bas_fcts, j, NULL));
2043 }
2044 } else {
2045 for (j = 0; j < n_phi; j++) {
2046 REAL val = 0.0;
2047 for (iq = 0; iq < n_points; iq++) {
2048 val +=
2049 quad_fast->phi[iq][j]
2050 *
2051 SCP_DOW(wdetf_qp[iq],
2052 PHI_D(quad_fast->bas_fcts, j, lambda[iq]));
2053 }
2054 fh->vec[dof->vec[j]] += val;
2055 }
2056 }
2057 }
2058 fh = CHAIN_NEXT(fh, DOF_REAL_VEC_D);
2059 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
2060 } CHAIN_WHILE(quad_fast, const QUAD_FAST);
2061
2062 } TRAVERSE_NEXT();
2063
2064 free_el_dof_vec(dof);
2065 }
2066 }
2067
2068 __FLATTEN_ATTRIBUTE__
trace_L2scp_fct_bas_dow(DOF_REAL_VEC_D * fh,FCT_D_AT_X f,MESH * trace_mesh,const QUAD * quad)2069 void trace_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh, FCT_D_AT_X f,
2070 MESH *trace_mesh,
2071 const QUAD *quad)
2072 {
2073 _AI_trace_L2scp_fct_bas_dow(fh, f, NULL, NULL, FILL_NOTHING,
2074 trace_mesh, quad);
2075 }
2076
2077 __FLATTEN_ATTRIBUTE__
trace_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D * fh,LOC_FCT_D_AT_QP f,void * fd,FLAGS fill_flag,MESH * trace_mesh,const QUAD * quad)2078 void trace_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
2079 LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
2080 MESH *trace_mesh,
2081 const QUAD *quad)
2082 {
2083 _AI_trace_L2scp_fct_bas_dow(fh, NULL, f, fd, fill_flag, trace_mesh, quad);
2084 }
2085
2086 /*--------------------------------------------------------------------------*/
2087 /* sets dirichlet boundary values g() in the vectors fh and uh and fills */
2088 /* the vector boundary with the boundary type of the corresponding DOF */
2089 /* return true if any Dirichlet nodes were found, false otherwise */
2090 /*--------------------------------------------------------------------------*/
2091
dirichlet_bound_loc(DOF_REAL_VEC * fh,DOF_REAL_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,LOC_FCT_AT_QP g,void * g_data,FLAGS fill_flags)2092 bool dirichlet_bound_loc(DOF_REAL_VEC *fh, DOF_REAL_VEC *uh,
2093 DOF_SCHAR_VEC *bound,
2094 const BNDRY_FLAGS dirichlet_segment,
2095 LOC_FCT_AT_QP g, void *g_data, FLAGS fill_flags)
2096 {
2097 FUNCNAME("dirichlet_bound_loc");
2098 MESH *mesh = NULL;
2099 const PARAMETRIC *parametric = NULL;
2100 int j, dim;
2101 const FE_SPACE *fe_space;
2102 const BAS_FCTS *bas_fcts = NULL;
2103 bool have_dirichlet = false;
2104 EL_REAL_VEC *vec_loc;
2105 BNDRY_FLAGS bndry_all;
2106
2107 if (!fh && !uh && !bound) {
2108 return false;
2109 }
2110
2111 if (dirichlet_segment == NULL) {
2112 BNDRY_FLAGS_ALL(bndry_all);
2113 dirichlet_segment = bndry_all;
2114 }
2115
2116 if (fh && uh) {
2117 TEST_EXIT(FE_SPACE_EQ_P(fh->fe_space, uh->fe_space),
2118 "fe_spaces in fh and uh differ\n");
2119 } else if (!uh && fh) {
2120 uh = fh;
2121 fh = NULL;
2122 }
2123
2124 if (uh == NULL && bound != NULL) {
2125
2126 fe_space = bound->fe_space;
2127 mesh = fe_space->mesh;
2128 bas_fcts = fe_space->bas_fcts;
2129 dim = mesh->dim;
2130
2131 /* Don't forget this: */
2132 FOREACH_DOF(fe_space,
2133 bound->vec[dof] = INTERIOR,
2134 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC));
2135
2136 /* If we have a periodic mesh but a non-periodic admin make sure
2137 * that we use non-periodic boundary information.
2138 */
2139 if (mesh->is_periodic && !(bound->fe_space->admin->flags & ADM_PERIODIC)) {
2140 fill_flags |= FILL_NON_PERIODIC;
2141 }
2142 fill_flags |= bas_fcts->fill_flags;
2143
2144 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_MACRO_WALLS|fill_flags) {
2145 int t, o, wall, n_tr_phi;
2146 const int *trmap;
2147
2148 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
2149 continue;
2150 }
2151
2152 t = el_info->el_type > 0;
2153 o = el_info->orientation < 0;
2154
2155 CHAIN_DO(fe_space, const FE_SPACE) {
2156 DOF dof[bas_fcts->n_bas_fcts];
2157
2158 GET_DOF_INDICES(bas_fcts, el_info->el, fe_space->admin, dof);
2159
2160 for (wall = 0; wall < N_WALLS(dim); wall++) {
2161 BNDRY_TYPE btype = wall_bound(el_info, wall);
2162
2163 if (btype == INTERIOR) {
2164 continue;
2165 }
2166
2167 /* Boundary conditions are attached to boundary segments,
2168 * most elements will be interior elements, and most
2169 * boundary elements will only have a single boundary facet,
2170 * so in almost all cases we will reach this code path only
2171 * once for each element.
2172 */
2173 n_tr_phi = bas_fcts->n_trace_bas_fcts[wall];
2174 trmap = bas_fcts->trace_dof_map[t][o][wall];
2175
2176 if (BNDRY_FLAGS_IS_AT_BNDRY(dirichlet_segment, btype)) {
2177 /* Mark boundary segment as DIRICHLET */
2178 for (j = 0; j < n_tr_phi; j++) {
2179 bound->vec[dof[trmap[j]]] = DIRICHLET;
2180 }
2181 have_dirichlet = true;
2182 } else { /* NEUMANN if not already DIRICHLET */
2183 for (j = 0; j < n_tr_phi; j++) {
2184 if (bound->vec[dof[trmap[j]]] == INTERIOR) {
2185 bound->vec[dof[trmap[j]]] = NEUMANN;
2186 }
2187 }
2188 }
2189 }
2190 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC );
2191 bas_fcts = CHAIN_NEXT(bas_fcts, const BAS_FCTS);
2192 } CHAIN_WHILE(fe_space, const FE_SPACE);
2193 } TRAVERSE_NEXT();
2194
2195 return have_dirichlet;
2196 }
2197
2198 /* uh != NULL below this point */
2199
2200 fe_space = uh->fe_space;
2201 mesh = fe_space->mesh;
2202 bas_fcts = uh->fe_space->bas_fcts;
2203 dim = mesh->dim;
2204
2205 /* If we have a periodic mesh but a non-periodic admin make sure
2206 * that we use non-periodic boundary information.
2207 */
2208 if (mesh->is_periodic && !(uh->fe_space->admin->flags & ADM_PERIODIC)) {
2209 fill_flags |= FILL_NON_PERIODIC;
2210 }
2211
2212 if (bound) {
2213 TEST_EXIT(FE_SPACE_EQ_P(uh->fe_space, bound->fe_space),
2214 "uh->fe_space, bound->fe_space differ!\n");
2215 FOREACH_DOF(fe_space,
2216 bound->vec[dof] = INTERIOR,
2217 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC));
2218 }
2219
2220 parametric = mesh->parametric;
2221
2222 vec_loc = get_el_real_vec(bas_fcts);
2223 fill_flags |= bas_fcts->fill_flags;
2224
2225 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_MACRO_WALLS|fill_flags) {
2226 int t, o, wall, n_tr_phi;
2227 const int *trmap;
2228
2229 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
2230 continue;
2231 }
2232
2233 if (parametric) {
2234 parametric->init_element(el_info, parametric);
2235 }
2236
2237 for (wall = 0; wall < N_WALLS(dim); wall++) {
2238 BNDRY_TYPE btype = wall_bound(el_info, wall);
2239
2240 if (btype == INTERIOR) {
2241 continue;
2242 }
2243
2244 /* Boundary conditions are attached to boundary segments, most
2245 * elements will be interior elements, and most boundary
2246 * elements will only have a single boundary facet, so in
2247 * almost all cases we will reach this code path only once for
2248 * each element.
2249 */
2250
2251 t = el_info->el_type > 0;
2252 o = el_info->orientation < 0;
2253
2254 if (!BNDRY_FLAGS_IS_AT_BNDRY(dirichlet_segment, btype)) {
2255 if (bound) {
2256 /* Mark boundary segment as NEUMANN */
2257 CHAIN_DO(fe_space, const FE_SPACE) {
2258 DOF dofvec[bas_fcts->n_bas_fcts];
2259
2260 GET_DOF_INDICES(bas_fcts, el_info->el, fe_space->admin, dofvec);
2261 n_tr_phi = bas_fcts->n_trace_bas_fcts[wall];
2262 trmap = bas_fcts->trace_dof_map[t][o][wall];
2263 for (j = 0; j < n_tr_phi; j++) {
2264 DOF dof = dofvec[trmap[j]];
2265 if (bound->vec[dof] == INTERIOR) {
2266 bound->vec[dof] = NEUMANN;
2267 }
2268 }
2269 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC);
2270 bas_fcts = CHAIN_NEXT(bas_fcts, const BAS_FCTS);
2271 } CHAIN_WHILE(fe_space, const FE_SPACE);
2272 }
2273 continue; /* not on a Dirichlet boundary segment */
2274 }
2275
2276 have_dirichlet = true;
2277
2278 /* Otherwise set Dirichlet boundary conditions for this
2279 * boundary segment. The coefficients returned in "dirichlet"
2280 * are relative to the trace-dof-mapping for the respective
2281 * type and orientation of the given wall.
2282 */
2283
2284 /* Could wrap this in blas-like inline functions for
2285 * element-vector -> global-vector updates
2286 */
2287 CHAIN_DO(fe_space, const FE_SPACE) {
2288 DOF dofvec[bas_fcts->n_bas_fcts];
2289
2290 GET_DOF_INDICES(bas_fcts, el_info->el, fe_space->admin, dofvec);
2291 n_tr_phi = bas_fcts->n_trace_bas_fcts[wall];
2292 trmap = bas_fcts->trace_dof_map[t][o][wall];
2293
2294 if (g) {
2295 INTERPOL(bas_fcts, vec_loc, el_info, wall, -1, NULL, g, g_data);
2296 }
2297
2298 for (j = 0; j < n_tr_phi; j++) {
2299 int ib = trmap[j];
2300 DOF dof = dofvec[ib];
2301 if (bound) {
2302 bound->vec[dof] = DIRICHLET;
2303 }
2304 uh->vec[dof] = g ? vec_loc->vec[ib] : 0.0;
2305 if (fh) {
2306 fh->vec[dof] = (bound == NULL || true) ? uh->vec[dof] : 0.0;
2307 }
2308 }
2309 bas_fcts = CHAIN_NEXT(bas_fcts, const BAS_FCTS);
2310 uh = CHAIN_NEXT(uh, DOF_REAL_VEC);
2311 vec_loc = CHAIN_NEXT(vec_loc, EL_REAL_VEC);
2312 fh = fh ? CHAIN_NEXT(fh, DOF_REAL_VEC) : NULL;
2313 bound = bound ? CHAIN_NEXT(bound, DOF_SCHAR_VEC) : NULL;
2314 } CHAIN_WHILE(fe_space, const FE_SPACE);
2315 } /* end of loop over the "walls" of this element. */
2316 } TRAVERSE_NEXT();
2317
2318 free_el_real_vec(vec_loc);
2319
2320 return have_dirichlet;
2321 }
2322
dirichlet_bound(DOF_REAL_VEC * fh,DOF_REAL_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,FCT_AT_X g)2323 bool dirichlet_bound(DOF_REAL_VEC *fh,
2324 DOF_REAL_VEC *uh,
2325 DOF_SCHAR_VEC *bound,
2326 const BNDRY_FLAGS dirichlet_segment,
2327 FCT_AT_X g)
2328 {
2329 const PARAMETRIC *parametric = NULL;
2330
2331 if (fh) {
2332 parametric = fh->fe_space->mesh->parametric;
2333 } else if (uh) {
2334 parametric = uh->fe_space->mesh->parametric;
2335 } else if (bound) {
2336 parametric = bound->fe_space->mesh->parametric;
2337 }
2338 if (parametric) {
2339 return dirichlet_bound_loc(fh, uh, bound, dirichlet_segment,
2340 _AI_inter_fct_loc_param, &g, FILL_COORDS);
2341 } else {
2342 return dirichlet_bound_loc(fh, uh, bound, dirichlet_segment,
2343 _AI_inter_fct_loc, &g, FILL_COORDS);
2344 }
2345 }
2346
dirichlet_bound_loc_dow(DOF_REAL_VEC_D * fh,DOF_REAL_VEC_D * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,LOC_FCT_D_AT_QP g,void * g_data,FLAGS fill_flags)2347 bool dirichlet_bound_loc_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
2348 DOF_SCHAR_VEC *bound,
2349 const BNDRY_FLAGS dirichlet_segment,
2350 LOC_FCT_D_AT_QP g, void *g_data, FLAGS fill_flags)
2351 {
2352 FUNCNAME("dirichlet_bound_loc_dow");
2353 MESH *mesh = NULL;
2354 const PARAMETRIC *parametric = NULL;
2355 int j, dim;
2356 const FE_SPACE *fe_space;
2357 const BAS_FCTS *bas_fcts = NULL;
2358 bool have_dirichlet = false;
2359 BNDRY_FLAGS bndry_all;
2360 EL_REAL_VEC_D *vec_loc;
2361
2362 if (!fh && !uh && !bound) {
2363 return false;
2364 }
2365
2366 if (dirichlet_segment == NULL) {
2367 BNDRY_FLAGS_ALL(bndry_all);
2368 dirichlet_segment = bndry_all;
2369 }
2370
2371 if (fh && uh) {
2372 TEST_EXIT(FE_SPACE_EQ_P(fh->fe_space, uh->fe_space),
2373 "fe_spaces in fh and uh differ\n");
2374 } else if (!uh && fh) {
2375 uh = fh;
2376 fh = NULL;
2377 }
2378
2379 if (uh == NULL && bound != NULL) {
2380
2381 fe_space = bound->fe_space;
2382 mesh = fe_space->mesh;
2383 bas_fcts = fe_space->bas_fcts;
2384 dim = mesh->dim;
2385
2386 /* Don't forget this: */
2387 FOREACH_DOF(fe_space,
2388 bound->vec[dof] = INTERIOR,
2389 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC));
2390
2391 /* If we have a periodic mesh but a non-periodic admin make sure
2392 * that we use non-periodic boundary information.
2393 */
2394 if (mesh->is_periodic && !(bound->fe_space->admin->flags & ADM_PERIODIC)) {
2395 fill_flags |= FILL_NON_PERIODIC;
2396 }
2397
2398 fill_flags |= bas_fcts->fill_flags;
2399 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_MACRO_WALLS|fill_flags) {
2400 int t, o, wall, n_tr_phi;
2401 const int *trmap;
2402
2403 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
2404 continue;
2405 }
2406
2407 t = el_info->el_type > 0;
2408 o = el_info->orientation < 0;
2409 CHAIN_DO(fe_space, const FE_SPACE) {
2410 DOF dof[bas_fcts->n_bas_fcts];
2411 GET_DOF_INDICES(bas_fcts, el_info->el, fe_space->admin, dof);
2412
2413 for (wall = 0; wall < N_WALLS(dim); wall++) {
2414 BNDRY_TYPE btype = wall_bound(el_info, wall);
2415
2416 if (btype == INTERIOR) {
2417 continue;
2418 }
2419
2420 /* Boundary conditions are attached to boundary segments,
2421 * most elements will be interior elements, and most
2422 * boundary elements will only have a single boundary facet,
2423 * so in almost all cases we will reach this code path only
2424 * once for each element.
2425 */
2426 n_tr_phi = bas_fcts->n_trace_bas_fcts[wall];
2427 trmap = bas_fcts->trace_dof_map[t][o][wall];
2428
2429 if (BNDRY_FLAGS_IS_AT_BNDRY(dirichlet_segment, btype)) {
2430 /* Mark boundary segment as DIRICHLET */
2431 for (j = 0; j < n_tr_phi; j++) {
2432 bound->vec[dof[trmap[j]]] = DIRICHLET;
2433 }
2434 have_dirichlet = true;
2435 } else { /* NEUMANN if not already DIRICHLET */
2436 for (j = 0; j < n_tr_phi; j++) {
2437 if (bound->vec[dof[trmap[j]]] == INTERIOR) {
2438 bound->vec[dof[trmap[j]]] = NEUMANN;
2439 }
2440 }
2441 }
2442 }
2443 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC);
2444 bas_fcts = CHAIN_NEXT(bas_fcts, const BAS_FCTS);
2445 } CHAIN_WHILE(fe_space, const FE_SPACE);
2446 } TRAVERSE_NEXT();
2447
2448 return have_dirichlet;
2449 }
2450
2451 /* uh != NULL below this point */
2452
2453 fe_space = uh->fe_space;
2454 mesh = fe_space->mesh;
2455 bas_fcts = uh->fe_space->bas_fcts;
2456 dim = mesh->dim;
2457
2458 /* If we have a periodic mesh but a non-periodic admin make sure
2459 * that we use non-periodic boundary information.
2460 */
2461 if (mesh->is_periodic && !(uh->fe_space->admin->flags & ADM_PERIODIC)) {
2462 fill_flags |= FILL_NON_PERIODIC;
2463 }
2464
2465 if (bound) {
2466 TEST_EXIT(FE_SPACE_EQ_P(uh->fe_space, bound->fe_space),
2467 "uh->fe_space, bound->fe_space differ!\n");
2468 FOREACH_DOF(fe_space,
2469 bound->vec[dof] = INTERIOR,
2470 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC));
2471 }
2472
2473 parametric = mesh->parametric;
2474
2475 vec_loc = get_el_real_vec_d(bas_fcts);
2476
2477 fill_flags |= bas_fcts->fill_flags;
2478 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_MACRO_WALLS|fill_flags) {
2479 int t, o, wall, n_tr_phi;
2480 const int *trmap;
2481
2482 if (INIT_ELEMENT(el_info, bas_fcts) == INIT_EL_TAG_NULL) {
2483 continue;
2484 }
2485
2486 if (parametric) {
2487 parametric->init_element(el_info, parametric);
2488 }
2489
2490 for (wall = 0; wall < N_WALLS(dim); wall++) {
2491 BNDRY_TYPE btype = wall_bound(el_info, wall);
2492
2493 if (btype == INTERIOR) {
2494 continue;
2495 }
2496
2497 /* Boundary conditions are attached to boundary segments, most
2498 * elements will be interior elements, and most boundary
2499 * elements will only have a single boundary facet, so in
2500 * almost all cases we will reach this code path only once for
2501 * each element.
2502 */
2503
2504 t = el_info->el_type > 0;
2505 o = el_info->orientation < 0;
2506
2507 if (!BNDRY_FLAGS_IS_AT_BNDRY(dirichlet_segment, btype)) {
2508 if (bound) {
2509 /* Mark boundary segment as NEUMANN */
2510 CHAIN_DO(fe_space, const FE_SPACE) {
2511 DOF dofvec[bas_fcts->n_bas_fcts];
2512
2513 GET_DOF_INDICES(bas_fcts, el_info->el, fe_space->admin, dofvec);
2514 n_tr_phi = bas_fcts->n_trace_bas_fcts[wall];
2515 trmap = bas_fcts->trace_dof_map[t][o][wall];
2516 for (j = 0; j < n_tr_phi; j++) {
2517 DOF dof = dofvec[trmap[j]];
2518 if (bound->vec[dof] == INTERIOR) {
2519 bound->vec[dof] = NEUMANN;
2520 }
2521 }
2522 bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC);
2523 bas_fcts = CHAIN_NEXT(bas_fcts, const BAS_FCTS);
2524 } CHAIN_WHILE(fe_space, const FE_SPACE);
2525 }
2526 continue; /* not on a Dirichlet boundary segment */
2527 }
2528
2529 have_dirichlet = true;
2530
2531 /* Otherwise set Dirichlet boundary conditions for this
2532 * boundary segment. The coefficients returned in "dirichlet"
2533 * are relative to the trace-dof-mapping for the respective
2534 * type and orientation of the given wall.
2535 */
2536
2537 /* Could wrap this in blas-like inline functions for
2538 * element-vector -> global-vector updates
2539 *
2540 * This following few lines are the only place where this
2541 * functions differs from its scalar counterpart.
2542 */
2543 CHAIN_DO(fe_space, const FE_SPACE) {
2544 DOF dofvec[bas_fcts->n_bas_fcts];
2545
2546 GET_DOF_INDICES(bas_fcts, el_info->el, fe_space->admin, dofvec);
2547 n_tr_phi = bas_fcts->n_trace_bas_fcts[wall];
2548 trmap = bas_fcts->trace_dof_map[t][o][wall];
2549
2550 if (uh->stride != 1) {
2551 /* quasi-scalar Cartesian product case */
2552 DOF_REAL_D_VEC *fhd = (DOF_REAL_D_VEC *)fh;
2553 DOF_REAL_D_VEC *uhd = (DOF_REAL_D_VEC *)uh;
2554
2555 if (g) {
2556 INTERPOL_D(bas_fcts, (EL_REAL_D_VEC *)vec_loc,
2557 el_info, wall, -1, NULL, g, g_data);
2558 }
2559
2560 for (j = 0; j < n_tr_phi; j++) {
2561 int ib = trmap[j];
2562 DOF dof = dofvec[ib];
2563 if (bound) {
2564 bound->vec[dof] = DIRICHLET;
2565 }
2566 if (g) {
2567 COPY_DOW(((EL_REAL_D_VEC *)vec_loc)->vec[ib], uhd->vec[dof]);
2568 } else {
2569 SET_DOW(0.0, uhd->vec[dof]);
2570 }
2571 if (fh) {
2572 if (bound == NULL || true) {
2573 COPY_DOW(uhd->vec[dof], fhd->vec[dof]);
2574 } else {
2575 SET_DOW(0.0, uhd->vec[dof]);
2576 }
2577 }
2578 }
2579 bas_fcts = CHAIN_NEXT(bas_fcts, const BAS_FCTS);
2580 uh = CHAIN_NEXT(uh, DOF_REAL_VEC_D);
2581 fh = fh ? CHAIN_NEXT(fh, DOF_REAL_VEC_D) : NULL;
2582 bound = bound ? CHAIN_NEXT(bound, DOF_SCHAR_VEC) : NULL;
2583 } else {
2584
2585 if (g) {
2586 INTERPOL_DOW(bas_fcts, vec_loc, el_info, wall, -1, NULL, g, g_data);
2587 }
2588
2589 for (j = 0; j < n_tr_phi; j++) {
2590 int ib = trmap[j];
2591 DOF dof = dofvec[ib];
2592 if (bound) {
2593 bound->vec[dof] = DIRICHLET;
2594 }
2595 uh->vec[dof] = g ? vec_loc->vec[ib] : 0.0;
2596 if (fh) {
2597 fh->vec[dof] = (bound == NULL || true) ? uh->vec[dof] : 0.0;
2598 }
2599 }
2600 bas_fcts = CHAIN_NEXT(bas_fcts, const BAS_FCTS);
2601 uh = CHAIN_NEXT(uh, DOF_REAL_VEC_D);
2602 fh = fh ? CHAIN_NEXT(fh, DOF_REAL_VEC_D) : NULL;
2603 bound = bound ? CHAIN_NEXT(bound, DOF_SCHAR_VEC) : NULL;
2604 }
2605 vec_loc = CHAIN_NEXT(vec_loc, EL_REAL_VEC_D);
2606 } CHAIN_WHILE(fe_space, const FE_SPACE);
2607 } /* end of loop over the "walls" of this element. */
2608 } TRAVERSE_NEXT();
2609
2610 free_el_real_vec_d(vec_loc);
2611
2612 return have_dirichlet;
2613 }
2614
dirichlet_bound_dow(DOF_REAL_VEC_D * fh,DOF_REAL_VEC_D * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,FCT_D_AT_X g)2615 bool dirichlet_bound_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
2616 DOF_SCHAR_VEC *bound,
2617 const BNDRY_FLAGS dirichlet_segment,
2618 FCT_D_AT_X g)
2619 {
2620 const PARAMETRIC *parametric = NULL;
2621
2622 if (fh) {
2623 parametric = fh->fe_space->mesh->parametric;
2624 } else if (uh) {
2625 parametric = uh->fe_space->mesh->parametric;
2626 } else if (bound) {
2627 parametric = bound->fe_space->mesh->parametric;
2628 }
2629 if (parametric) {
2630 return dirichlet_bound_loc_dow(fh, uh, bound, dirichlet_segment,
2631 _AI_inter_fct_loc_d_param, &g, FILL_COORDS);
2632 } else {
2633 return dirichlet_bound_loc_dow(fh, uh, bound, dirichlet_segment,
2634 _AI_inter_fct_loc_d, &g, FILL_COORDS);
2635 }
2636 }
2637
2638 /* Compute the mean value of either a finite element function or a
2639 * non-discrete function. If both are given return the difference of
2640 * their mean values (f-fh).
2641 */
2642 static inline
_AI_mean_value(MESH * mesh,FCT_AT_X f,LOC_FCT_AT_QP f_loc,void * fd,FLAGS fill_flags,const DOF_REAL_VEC * fh,const QUAD * quad)2643 REAL _AI_mean_value(MESH *mesh,
2644 FCT_AT_X f,
2645 LOC_FCT_AT_QP f_loc, void *fd, FLAGS fill_flags,
2646 const DOF_REAL_VEC *fh,
2647 const QUAD *quad)
2648 {
2649 /*FUNCNAME("mean_value");*/
2650 const QUAD_FAST *quad_fast = NULL;
2651 REAL accu1, accu2, det;
2652 int iq, dim;
2653 const BAS_FCTS *bas_fcts = NULL;
2654 REAL mean;
2655 REAL meas;
2656
2657 dim = mesh->dim;
2658
2659 if (!quad) {
2660 quad = get_quadrature(dim, fh ? fh->fe_space->bas_fcts->degree : 1);
2661 }
2662
2663 if (fh) {
2664 quad_fast = get_quad_fast(fh->fe_space->bas_fcts, quad, INIT_PHI);
2665 bas_fcts = fh->fe_space->bas_fcts;
2666 }
2667
2668 mean = 0.0;
2669 meas = 0.0;
2670 {
2671 PARAMETRIC *parametric = mesh->parametric;
2672 REAL det_p[quad->n_points_max];
2673 REAL fh_vec[quad->n_points_max];
2674 REAL f_vec[quad->n_points_max];
2675 EL_REAL_VEC *fh_loc;
2676 int is_parametric = false;
2677 int dim = mesh->dim;
2678
2679 fh_loc = get_el_real_vec(bas_fcts);
2680
2681 fill_flags |= quad_fast->fill_flags;
2682 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_COORDS|fill_flags) {
2683
2684 if (fh) {
2685 if (INIT_ELEMENT(el_info, quad_fast) == INIT_EL_TAG_NULL) {
2686 continue;
2687 }
2688 fill_el_real_vec(fh_loc, el_info->el, fh);
2689 }
2690
2691 if (parametric) {
2692 is_parametric = parametric->init_element(el_info, parametric);
2693 }
2694
2695 if (f && fh) {
2696 uh_at_qp(fh_vec, quad_fast, fh_loc);
2697 fx_at_qp(f_vec, el_info, quad, f);
2698
2699 for (iq = 0; iq < quad->n_points; iq++) {
2700 f_vec[iq] -= fh_vec[iq];
2701 }
2702 } else if (f_loc && fh) {
2703 uh_at_qp(fh_vec, quad_fast, fh_loc);
2704 f_loc_at_qp(f_vec, el_info, quad, f_loc, fd);
2705
2706 for (iq = 0; iq < quad->n_points; iq++) {
2707 f_vec[iq] -= fh_vec[iq];
2708 }
2709 } else if (f) {
2710 fx_at_qp(f_vec, el_info, quad, f);
2711 } else if (f_loc) {
2712 f_loc_at_qp(f_vec, el_info, quad, f_loc, fd);
2713 } else if (fh) {
2714 uh_at_qp(f_vec, quad_fast, fh_loc);
2715 }
2716
2717 if (is_parametric) {
2718 parametric->det(el_info, quad, 0, NULL, det_p);
2719 for (iq = 0; iq < quad->n_points; iq++) {
2720 mean += quad->w[iq]*det_p[iq]*f_vec[iq];
2721 meas += quad->w[iq]*det_p[iq];
2722 }
2723 } else {
2724 det = el_det_dim(dim, el_info);
2725 accu1 = accu2 = 0.0;
2726 for (iq = 0; iq < quad->n_points; iq++) {
2727 accu1 += quad->w[iq]*f_vec[iq];
2728 accu2 += quad->w[iq];
2729 }
2730 mean += det*accu1;
2731 meas += det*accu2;
2732 }
2733 } TRAVERSE_NEXT();
2734 }
2735
2736 return mean/meas;
2737 }
2738
2739 __FLATTEN_ATTRIBUTE__
mean_value(MESH * mesh,REAL (* f)(const REAL_D),const DOF_REAL_VEC * fh,const QUAD * quad)2740 REAL mean_value(MESH *mesh,
2741 REAL (*f)(const REAL_D),
2742 const DOF_REAL_VEC *fh,
2743 const QUAD *quad)
2744 {
2745 return _AI_mean_value(mesh, f, NULL, NULL, FILL_NOTHING, fh, quad);
2746 }
2747
2748 __FLATTEN_ATTRIBUTE__
mean_value_loc(MESH * mesh,REAL (* f_loc)(const EL_INFO * el_info,const QUAD * quad,int iq,void * fd),void * fd,FLAGS fill_flags,const DOF_REAL_VEC * fh,const QUAD * quad)2749 REAL mean_value_loc(MESH *mesh,
2750 REAL (*f_loc)(const EL_INFO *el_info,
2751 const QUAD *quad, int iq,
2752 void *fd),
2753 void *fd, FLAGS fill_flags,
2754 const DOF_REAL_VEC *fh,
2755 const QUAD *quad)
2756 {
2757 return _AI_mean_value(mesh, NULL, f_loc, fd, fill_flags, fh, quad);
2758 }
2759
2760 static const REAL *
_AI_mean_value_dow(MESH * mesh,FCT_D_AT_X f,LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flags,const DOF_REAL_VEC_D * fh,const QUAD * quad,REAL_D mean)2761 _AI_mean_value_dow(MESH *mesh,
2762 FCT_D_AT_X f,
2763 LOC_FCT_D_AT_QP f_loc, void *fd, FLAGS fill_flags,
2764 const DOF_REAL_VEC_D *fh,
2765 const QUAD *quad,
2766 REAL_D mean)
2767 {
2768 /* FUNCNAME("mean_value_dow"); */
2769 static REAL_D _mean;
2770 const QUAD_FAST *quad_fast = NULL;
2771 REAL det;
2772 int iq, dim = mesh->dim;
2773 REAL meas;
2774
2775 if (!mean) {
2776 mean = _mean;
2777 }
2778
2779 if (!quad) {
2780 quad = get_quadrature(dim, fh ? fh->fe_space->bas_fcts->degree : 1);
2781 }
2782
2783 if (fh) {
2784 quad_fast = get_quad_fast(fh->fe_space->bas_fcts, quad, INIT_PHI);
2785 }
2786
2787 SET_DOW(0.0, mean);
2788 meas = 0.0;
2789 {
2790 PARAMETRIC *parametric = mesh->parametric;
2791 REAL_D accu1;
2792 REAL accu2;
2793 REAL det_p[quad->n_points_max];
2794 REAL_D fh_vec[quad->n_points_max];
2795 REAL_D f_vec[quad->n_points_max];
2796 int is_parametric = false;
2797 int dim = mesh->dim;
2798 const EL_REAL_VEC_D *fh_loc = NULL;
2799
2800 fill_flags |= quad_fast->fill_flags;
2801 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_COORDS|fill_flags) {
2802
2803 if (fh) {
2804 if (INIT_ELEMENT(el_info, quad_fast) == INIT_EL_TAG_NULL) {
2805 continue;
2806 }
2807 fh_loc = fill_el_real_vec_d(NULL, el_info->el, fh);
2808 }
2809
2810 if (parametric) {
2811 is_parametric = parametric->init_element(el_info, parametric);
2812 }
2813
2814 if (f && fh) {
2815 uh_dow_at_qp(fh_vec, quad_fast, fh_loc);
2816 fx_d_at_qp(f_vec, el_info, quad, f);
2817
2818 for (iq = 0; iq < quad->n_points; iq++) {
2819 AXPY_DOW(-1.0, fh_vec[iq], f_vec[iq]);
2820 }
2821 } else if (f_loc && fh) {
2822 uh_dow_at_qp(fh_vec, quad_fast, fh_loc);
2823 f_loc_d_at_qp(f_vec, el_info, quad, f_loc, fd);
2824
2825 for (iq = 0; iq < quad->n_points; iq++) {
2826 AXPY_DOW(-1.0, fh_vec[iq], f_vec[iq]);
2827 }
2828 } else if (f) {
2829 fx_d_at_qp(f_vec, el_info, quad, f);
2830 } else if (f_loc) {
2831 f_loc_d_at_qp(f_vec, el_info, quad, f_loc, fd);
2832 } else if (fh) {
2833 uh_dow_at_qp(f_vec, quad_fast, fh_loc);
2834 }
2835
2836 if (is_parametric) {
2837 parametric->det(el_info, quad, 0, NULL, det_p);
2838 for (iq = 0; iq < quad->n_points; iq++) {
2839 AXPY_DOW(quad->w[iq]*det_p[iq], f_vec[iq], mean);
2840 meas += quad->w[iq]*det_p[iq];
2841 }
2842 } else {
2843 det = el_det_dim(dim, el_info);
2844 SET_DOW(0.0, accu1);
2845 accu2 = 0.0;
2846 for (iq = 0; iq < quad->n_points; iq++) {
2847 AXPY_DOW(quad->w[iq], f_vec[iq], accu1);
2848 accu2 += quad->w[iq];
2849 }
2850 AXPY_DOW(det, accu1, mean);
2851 meas += det*accu2;
2852 }
2853 } TRAVERSE_NEXT();
2854 }
2855
2856 AX_DOW(1.0/meas, mean);
2857
2858 return mean;
2859 }
2860
2861 __FLATTEN_ATTRIBUTE__
mean_value_dow(MESH * mesh,FCT_D_AT_X f,const DOF_REAL_VEC_D * fh,const QUAD * quad,REAL_D mean)2862 const REAL *mean_value_dow(MESH *mesh,
2863 FCT_D_AT_X f, const DOF_REAL_VEC_D *fh,
2864 const QUAD *quad,
2865 REAL_D mean)
2866 {
2867 return _AI_mean_value_dow(mesh, f, NULL, NULL, FILL_NOTHING, fh, quad, mean);
2868 }
2869
2870 __FLATTEN_ATTRIBUTE__
mean_value_loc_dow(REAL_D mean,MESH * mesh,LOC_FCT_D_AT_QP f_loc,void * fd,FLAGS fill_flags,const DOF_REAL_VEC_D * fh,const QUAD * quad)2871 const REAL *mean_value_loc_dow(REAL_D mean,
2872 MESH *mesh,
2873 LOC_FCT_D_AT_QP f_loc, void *fd,
2874 FLAGS fill_flags,
2875 const DOF_REAL_VEC_D *fh,
2876 const QUAD *quad)
2877 {
2878 return _AI_mean_value_dow(mesh, NULL, f_loc, fd, fill_flags, fh, quad, mean);
2879 }
2880
2881 typedef struct robin_ass_data {
2882 struct robin_ass_data *next;
2883 const REAL *dets;
2884 BNDRY_FLAGS bndry_seg;
2885 REAL factor;
2886 REAL exponent;
2887 } ROBIN_ASS_DATA;
2888
h2_from_det(int dim,REAL det)2889 static inline REAL h2_from_det(int dim, REAL det)
2890 {
2891 FUNCNAME("h2_from_det");
2892
2893 switch(dim) {
2894 case 0:
2895 return 1.0;
2896 case 1:
2897 return det*det;
2898 case 2:
2899 return det;
2900 case 3:
2901 return pow(det, 2.0/3.0);
2902 default:
2903 ERROR_EXIT("Illegal dim!\n");
2904 return 0.0; /* shut up the compiler */
2905 }
2906 }
2907
2908 /* Implementation of Robin boundary conditions */
param_robin_bndry_init_element(const EL_INFO * el_info,int wall,const WALL_QUAD * quad[3],void * vd)2909 static bool param_robin_bndry_init_element(const EL_INFO *el_info, int wall,
2910 const WALL_QUAD *quad[3], void *vd)
2911 {
2912 ROBIN_ASS_DATA *data = (ROBIN_ASS_DATA *)vd;
2913 PARAMETRIC *parametric = el_info->mesh->parametric;
2914 int is_parametric = false;
2915 const EL_GEOM_CACHE *elgc;
2916 const QUAD_EL_CACHE *qelc;
2917 REAL area = 0.0;
2918
2919 is_parametric = parametric->init_element(el_info, parametric);
2920
2921 if (is_parametric) {
2922 int iq;
2923 INIT_ELEMENT(el_info, &quad[0]->quad[wall]);
2924 qelc = fill_quad_el_cache(
2925 el_info, &quad[0]->quad[wall], FILL_EL_QUAD_WALL_DET);
2926 data->dets = qelc->param.wall_det;
2927 if (data->exponent > 0.0) {
2928 area = 0.0;
2929 for (iq = 0; iq < quad[0]->quad[wall].n_points; ++iq) {
2930 area += quad[0]->quad[wall].w[iq] * data->dets[iq];
2931 }
2932 }
2933 } else {
2934 elgc = fill_el_geom_cache(el_info, FILL_EL_WALL_DET(wall));
2935 data->dets = &elgc->wall_det[wall];
2936 area = *data->dets;
2937 }
2938
2939 if (data->exponent > 0.0) {
2940 data->factor =
2941 pow(h2_from_det(el_info->mesh->dim - 1, area), -0.5 * data->exponent);
2942 } else {
2943 data->factor = 1.0;
2944 }
2945
2946 return is_parametric;
2947 }
2948
robin_bndry_init_element(const EL_INFO * el_info,int wall,const WALL_QUAD * quad[3],void * vd)2949 static bool robin_bndry_init_element(const EL_INFO *el_info, int wall,
2950 const WALL_QUAD *quad[3], void *vd)
2951 {
2952 ROBIN_ASS_DATA *data = (ROBIN_ASS_DATA *)vd;
2953 const EL_GEOM_CACHE *elgc;
2954
2955 elgc = fill_el_geom_cache(el_info, FILL_EL_WALL_DET(wall));
2956 data->dets = &elgc->wall_det[wall];
2957
2958 if (data->exponent > 0.0) {
2959 data->factor = pow(
2960 h2_from_det(el_info->mesh->dim - 1, *data->dets), -0.5 * data->exponent);
2961 } else {
2962 data->factor = 1.0;
2963 }
2964
2965 return false;
2966 }
2967
robin_bndry_c(const EL_INFO * el_info,const QUAD * quad,int iq,void * vd)2968 static REAL robin_bndry_c(const EL_INFO *el_info,
2969 const QUAD *quad, int iq, void *vd)
2970 {
2971 ROBIN_ASS_DATA *data = (ROBIN_ASS_DATA *)vd;
2972
2973 /* iq will never be != 0 for pw constant and non-param */
2974 return data->factor * data->dets[iq];
2975 }
2976
2977 /* Add Robin boundary condition to a DOF-matrix. */
robin_bound_matrix_info(EL_MATRIX_INFO * robin_info,const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space,const BNDRY_FLAGS robin_segment,REAL alpha_r,const WALL_QUAD * wall_quad,REAL exponent)2978 void robin_bound_matrix_info(EL_MATRIX_INFO *robin_info,
2979 const FE_SPACE *row_fe_space,
2980 const FE_SPACE *col_fe_space,
2981 const BNDRY_FLAGS robin_segment,
2982 REAL alpha_r,
2983 const WALL_QUAD *wall_quad,
2984 REAL exponent)
2985 {
2986 MESH *mesh;
2987 const FE_SPACE *fe_space;
2988 const BAS_FCTS *bas_fcts;
2989 BNDRY_OPERATOR_INFO oinfo = { NULL, };
2990 static ROBIN_ASS_DATA *first_data;
2991 ROBIN_ASS_DATA *data;
2992 BNDRY_FLAGS segment;
2993
2994 if (robin_segment == NULL) {
2995 BNDRY_FLAGS_ALL(segment);
2996 } else {
2997 BNDRY_FLAGS_CPY(segment, robin_segment);
2998 }
2999
3000 for (data = first_data; data != NULL; data = data->next) {
3001 if (data->factor == alpha_r &&
3002 data->exponent == exponent &&
3003 BNDRY_FLAGS_CMP(data->bndry_seg, segment)) {
3004 break;
3005 }
3006 }
3007 if (data == NULL) {
3008 data = MEM_CALLOC(1, ROBIN_ASS_DATA);
3009 data->next = first_data;
3010 first_data = data;
3011 BNDRY_FLAGS_CPY(data->bndry_seg, robin_segment);
3012 BNDRY_FLAGS_MARK_BNDRY(data->bndry_seg);
3013 data->factor = alpha_r;
3014 data->exponent = exponent;
3015 }
3016
3017 fe_space = row_fe_space;
3018 mesh = fe_space->mesh;
3019 bas_fcts = fe_space->bas_fcts;
3020
3021 if (wall_quad == NULL) {
3022 wall_quad = get_wall_quad(mesh->dim, 2*bas_fcts->degree);
3023 }
3024
3025 oinfo.row_fe_space = row_fe_space;
3026 oinfo.col_fe_space = col_fe_space;
3027 oinfo.quad[0] = wall_quad;
3028 if (mesh->parametric) {
3029 oinfo.init_element = param_robin_bndry_init_element;
3030 } else {
3031 oinfo.init_element = robin_bndry_init_element;
3032 }
3033 oinfo.c.real = robin_bndry_c;
3034 oinfo.c_pw_const = true;
3035 BNDRY_FLAGS_CPY(oinfo.bndry_type, data->bndry_seg);
3036 oinfo.user_data = data;
3037 oinfo.fill_flag = FILL_COORDS|CALL_LEAF_EL;
3038 fill_matrix_info_ext(robin_info, NULL, &oinfo, NULL);
3039 robin_info->factor = alpha_r;
3040 }
3041
3042 /* Add Robin boundary condition to a DOF-matrix. */
robin_bound(DOF_MATRIX * matrix,const BNDRY_FLAGS robin_segment,REAL alpha_r,const WALL_QUAD * wall_quad,REAL exponent)3043 void robin_bound(DOF_MATRIX *matrix,
3044 const BNDRY_FLAGS robin_segment,
3045 REAL alpha_r,
3046 const WALL_QUAD *wall_quad,
3047 REAL exponent)
3048 {
3049 EL_MATRIX_INFO robin_info;
3050
3051 if (alpha_r == 0.0 || matrix == NULL) {
3052 return;
3053 }
3054
3055 robin_bound_matrix_info(&robin_info,
3056 matrix->row_fe_space,
3057 matrix->col_fe_space
3058 ? matrix->col_fe_space : matrix->row_fe_space,
3059 robin_segment,
3060 alpha_r,
3061 wall_quad,
3062 exponent);
3063
3064 update_matrix(matrix, &robin_info, NoTranspose);
3065 }
3066
one(const REAL_D x)3067 static REAL one(const REAL_D x)
3068 {
3069 return 1.0;
3070 }
3071
one_d(const REAL_D x,REAL_D result)3072 static const REAL *one_d(const REAL_D x, REAL_D result)
3073 {
3074 static REAL_D space;
3075
3076 if (result == NULL) {
3077 result = space;
3078 }
3079
3080 SET_DOW(1.0, result);
3081
3082 return result;
3083 }
3084
3085 /* One more convenience function: implement Dirichlet and
3086 * Neumann/Robin boundary conditions. If we have pure Neumann boundary
3087 * conditions, we also perform a mean-value correction of the RHS;
3088 * otherwise the compatibility condition for the RHS may not be
3089 * satisfied (because of quadrature errors).
3090 */
_AI_boundary_conditions(DOF_MATRIX * matrix,DOF_REAL_VEC * fh,DOF_REAL_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,FCT_AT_X g,REAL (* gn)(const REAL_D x,const REAL_D normal),LOC_FCT_AT_QP g_loc,LOC_FCT_AT_QP gn_loc,void * fd,FLAGS fill_flag,REAL alpha_r,const WALL_QUAD * wall_quad)3091 static bool _AI_boundary_conditions(DOF_MATRIX *matrix,
3092 DOF_REAL_VEC *fh,
3093 DOF_REAL_VEC *uh,
3094 DOF_SCHAR_VEC *bound,
3095 const BNDRY_FLAGS dirichlet_segment,
3096 FCT_AT_X g,
3097 REAL (*gn)(const REAL_D x,
3098 const REAL_D normal),
3099 LOC_FCT_AT_QP g_loc,
3100 LOC_FCT_AT_QP gn_loc,
3101 void *fd, FLAGS fill_flag,
3102 REAL alpha_r,
3103 const WALL_QUAD *wall_quad)
3104 {
3105 FUNCNAME("boundary_conditions");
3106 bool not_all_neumann = false;
3107 BNDRY_FLAGS non_dirichlet, dirichlet;
3108
3109 BNDRY_FLAGS_ALL(non_dirichlet);
3110 if (dirichlet_segment) {
3111 BNDRY_FLAGS_XOR(non_dirichlet, dirichlet_segment);
3112 BNDRY_FLAGS_MARK_BNDRY(non_dirichlet);
3113 } else {
3114 BNDRY_FLAGS_INIT(dirichlet);
3115 dirichlet_segment = dirichlet;
3116 }
3117
3118 if (fh && (gn || gn_loc)) {
3119 not_all_neumann =
3120 _AI_bndry_L2scp_fct_bas(fh, gn_loc, fd, fill_flag, gn,
3121 non_dirichlet, wall_quad);
3122 }
3123 if (alpha_r > 0.0 && matrix) {
3124 robin_bound(matrix, non_dirichlet, alpha_r, wall_quad, 0.0);
3125 not_all_neumann = true;
3126 }
3127 if (((g || g_loc) && (fh || uh)) || bound) {
3128 if (g_loc) {
3129 not_all_neumann =
3130 dirichlet_bound_loc(fh, uh, bound, dirichlet_segment,
3131 g_loc, fd, fill_flag) || not_all_neumann;
3132 } else {
3133 not_all_neumann =
3134 dirichlet_bound(fh, uh, bound, dirichlet_segment, g) || not_all_neumann;
3135 }
3136 }
3137
3138 /* do a mean-value correction for the rhs if we have pure Neumann
3139 * boundary-conditions. We use the combination "bndry_type ==
3140 * NEUMANN && alpha_r == -1.0" to flag the request for automatic
3141 * mean-value adjustment.
3142 */
3143 if (!not_all_neumann && fh && alpha_r < 0.0) {
3144 const DOF_ADMIN *admin = fh->fe_space->admin;
3145 const BAS_FCTS *bas_fcts = fh->fe_space->bas_fcts;
3146 REAL mv;
3147 int dim = bas_fcts->unchained->dim;
3148 int degree = bas_fcts->unchained->degree;
3149
3150 if (bas_fcts->unchained == get_lagrange(dim, degree)
3151 ||
3152 bas_fcts->unchained == get_discontinuous_lagrange(dim, degree)) {
3153
3154 mv = 0.0;
3155 FOR_ALL_DOFS(admin, mv += fh->vec[dof]);
3156
3157 if (fabs(mv) > REAL_EPSILON) {
3158 MSG("Mean-value adjustment by %e\n", mv);
3159 mv /= (REAL)admin->used_count;
3160 FOR_ALL_DOFS(admin, fh->vec[dof] -= mv);
3161 }
3162 } else {
3163 /* More work to do. The code below works iff the interpolation
3164 * operator preserves polynomials and if the constant 1 function
3165 * belongs to the finite element space.
3166 */
3167 DOF_REAL_VEC *one_h = get_dof_real_vec("one", fh->fe_space);
3168 REAL mv, weight;
3169
3170 interpol(one, one_h);
3171 weight = mv = 0.0;
3172 FOR_ALL_DOFS(admin, if (fabs(one_h->vec[dof]) > 100*REAL_EPSILON) {
3173 weight += one_h->vec[dof];
3174 mv += fh->vec[dof] * one_h->vec[dof];
3175 } else {
3176 one_h->vec[dof] = 0.0;
3177 });
3178 if (fabs(mv) > REAL_EPSILON) {
3179 MSG("Mean-value adjustment by %e\n", mv);
3180 mv /= weight;
3181 FOR_ALL_DOFS(admin, if (one_h->vec[dof] != 0.0) {
3182 fh->vec[dof] -= mv;
3183 });
3184 }
3185 free_dof_real_vec(one_h);
3186 }
3187 }
3188
3189 return not_all_neumann;
3190 }
3191
boundary_conditions(DOF_MATRIX * matrix,DOF_REAL_VEC * fh,DOF_REAL_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,FCT_AT_X g,REAL (* gn)(const REAL_D x,const REAL_D normal),REAL alpha_r,const WALL_QUAD * wall_quad)3192 bool boundary_conditions(DOF_MATRIX *matrix,
3193 DOF_REAL_VEC *fh,
3194 DOF_REAL_VEC *uh,
3195 DOF_SCHAR_VEC *bound,
3196 const BNDRY_FLAGS dirichlet_segment,
3197 FCT_AT_X g,
3198 REAL (*gn)(const REAL_D x, const REAL_D normal),
3199 REAL alpha_r,
3200 const WALL_QUAD *wall_quad)
3201 {
3202 return _AI_boundary_conditions(matrix, fh, uh, bound,
3203 dirichlet_segment,
3204 g, gn,
3205 NULL, NULL, NULL, FILL_NOTHING,
3206 alpha_r, wall_quad);
3207 }
3208
boundary_conditions_loc(DOF_MATRIX * matrix,DOF_REAL_VEC * fh,DOF_REAL_VEC * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,LOC_FCT_AT_QP g_loc,LOC_FCT_AT_QP gn_loc,void * fd,FLAGS fill_flag,REAL alpha_r,const WALL_QUAD * wall_quad)3209 bool boundary_conditions_loc(DOF_MATRIX *matrix,
3210 DOF_REAL_VEC *fh,
3211 DOF_REAL_VEC *uh,
3212 DOF_SCHAR_VEC *bound,
3213 const BNDRY_FLAGS dirichlet_segment,
3214 LOC_FCT_AT_QP g_loc,
3215 LOC_FCT_AT_QP gn_loc,
3216 void *fd, FLAGS fill_flag,
3217 REAL alpha_r,
3218 const WALL_QUAD *wall_quad)
3219 {
3220 return _AI_boundary_conditions(matrix, fh, uh, bound,
3221 dirichlet_segment,
3222 NULL, NULL,
3223 g_loc, gn_loc, fd, fill_flag,
3224 alpha_r, wall_quad);
3225 }
3226
3227 static bool
_AI_boundary_conditions_dow(DOF_MATRIX * matrix,DOF_REAL_VEC_D * fh,DOF_REAL_VEC_D * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,FCT_D_AT_X g,const REAL * (* gn)(const REAL_D x,const REAL_D normal,REAL_D res),LOC_FCT_D_AT_QP g_loc,LOC_FCT_D_AT_QP gn_loc,void * fd,FLAGS fill_flag,REAL alpha_r,const WALL_QUAD * wall_quad)3228 _AI_boundary_conditions_dow(DOF_MATRIX *matrix,
3229 DOF_REAL_VEC_D *fh,
3230 DOF_REAL_VEC_D *uh,
3231 DOF_SCHAR_VEC *bound,
3232 const BNDRY_FLAGS dirichlet_segment,
3233 FCT_D_AT_X g,
3234 const REAL *(*gn)(const REAL_D x,
3235 const REAL_D normal,
3236 REAL_D res),
3237 LOC_FCT_D_AT_QP g_loc,
3238 LOC_FCT_D_AT_QP gn_loc, void *fd, FLAGS fill_flag,
3239 REAL alpha_r,
3240 const WALL_QUAD *wall_quad)
3241 {
3242 FUNCNAME("boundary_conditions_dow");
3243 bool not_all_neumann = false;
3244 BNDRY_FLAGS non_dirichlet, dirichlet;
3245
3246 BNDRY_FLAGS_ALL(non_dirichlet);
3247 if (dirichlet_segment) {
3248 BNDRY_FLAGS_XOR(non_dirichlet, dirichlet_segment);
3249 BNDRY_FLAGS_MARK_BNDRY(non_dirichlet);
3250 } else {
3251 BNDRY_FLAGS_INIT(dirichlet);
3252 dirichlet_segment = dirichlet;
3253 }
3254
3255 if (fh && (gn || gn_loc)) {
3256 not_all_neumann =
3257 _AI_bndry_L2scp_fct_bas_dow(fh, gn_loc, fd, fill_flag, gn,
3258 non_dirichlet, wall_quad);
3259 }
3260 if (alpha_r > 0.0 && matrix) {
3261 robin_bound(matrix, non_dirichlet, alpha_r, wall_quad, 0.0);
3262 not_all_neumann = true;
3263 }
3264 if (((g && g_loc) && (fh || uh)) || bound) {
3265 if (g_loc) {
3266 not_all_neumann =
3267 dirichlet_bound_loc_dow(fh, uh, bound, dirichlet_segment,
3268 g_loc, fd, fill_flag)
3269 ||
3270 not_all_neumann;
3271 } else {
3272 not_all_neumann =
3273 dirichlet_bound_dow(fh, uh, bound, dirichlet_segment, g)
3274 ||
3275 not_all_neumann;
3276 }
3277 }
3278 /* do a mean-value correction for the rhs if we have pure Neumann
3279 * boundary-conditions. We use the combination "bndry_type ==
3280 * NEUMANN && alpha_r == -1.0" to flag the request for automatic
3281 * mean-value adjustment.
3282 */
3283 if (!not_all_neumann && fh && alpha_r < 0.0) {
3284 /* do a mean-value correction for the rhs if we have pure Neumann
3285 * boundary-conditions. Hope that the first component of a direct
3286 * sum corresponds to the Lagrange basis functions.
3287 */
3288 const DOF_ADMIN *admin = fh->fe_space->admin;
3289 const BAS_FCTS *bas_fcts = fh->fe_space->bas_fcts;
3290 DOF_REAL_D_VEC *fhd = (DOF_REAL_D_VEC *)fh;
3291 REAL_D mv;
3292 int dim = bas_fcts->unchained->dim;
3293 int degree = bas_fcts->unchained->degree;
3294
3295 if (bas_fcts->unchained == get_lagrange(dim, degree)
3296 ||
3297 bas_fcts->unchained == get_discontinuous_lagrange(dim, degree)) {
3298
3299 SET_DOW(0.0, mv);
3300 FOR_ALL_DOFS(admin, AXPY_DOW(1.0, fhd->vec[dof], mv));
3301 if (NORM_DOW(mv) > REAL_EPSILON) {
3302 MSG("Mean-value adjustment by %e\n", NORM_DOW(mv));
3303 }
3304 SCAL_DOW(1.0/(REAL)admin->used_count, mv);
3305 FOR_ALL_DOFS(admin, AXPY_DOW(-1.0, mv, fhd->vec[dof]));
3306 } else if (fh->stride == DIM_OF_WORLD) {
3307 /* More work to do. The code below works iff the interpolation
3308 * operator preserves polynomials and if the constant 1 function
3309 * belongs to the finite element space.
3310 */
3311 DOF_REAL_D_VEC *one_h =
3312 (DOF_REAL_D_VEC *)get_dof_real_vec_d("one", fhd->fe_space);
3313 REAL_D weight;
3314 int i;
3315
3316 interpol_d(one_d, one_h);
3317 SET_DOW(0.0, weight);
3318 SET_DOW(0.0, mv);
3319 FOR_ALL_DOFS(admin, if (NORM8_DOW(one_h->vec[dof]) > 100*REAL_EPSILON) {
3320 AXPY_DOW(1.0, one_h->vec[dof], weight);
3321 for (i = 0; i < DIM_OF_WORLD; i++) {
3322 mv[i] += one_h->vec[dof][i] * fhd->vec[dof][i];
3323 }
3324 } else {
3325 SET_DOW(0.0, one_h->vec[dof]);
3326 });
3327 if (NORM8_DOW(mv) > REAL_EPSILON) {
3328 MSG("Mean-value adjustment by "FORMAT_DOW"\n", EXPAND_DOW(mv));
3329 for (i = 0; i < DIM_OF_WORLD; i++) {
3330 mv[i] /= weight[i];
3331 }
3332 FOR_ALL_DOFS(admin, if (one_h->vec[dof][0] != 0.0) {
3333 AXPY_DOW(-1.0, mv, fhd->vec[dof]);
3334 });
3335 }
3336 free_dof_real_d_vec(one_h);
3337 } else {
3338 ERROR_EXIT("Sorry, mean-value correction "
3339 "not implemented for exotic finite element spaces.\n");
3340 }
3341 }
3342
3343 return not_all_neumann;
3344 }
3345
boundary_conditions_dow(DOF_MATRIX * matrix,DOF_REAL_VEC_D * fh,DOF_REAL_VEC_D * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,FCT_D_AT_X g,const REAL * (* gn)(const REAL_D x,const REAL_D normal,REAL_D res),REAL alpha_r,const WALL_QUAD * wall_quad)3346 bool boundary_conditions_dow(DOF_MATRIX *matrix,
3347 DOF_REAL_VEC_D *fh,
3348 DOF_REAL_VEC_D *uh,
3349 DOF_SCHAR_VEC *bound,
3350 const BNDRY_FLAGS dirichlet_segment,
3351 FCT_D_AT_X g,
3352 const REAL *(*gn)(const REAL_D x,
3353 const REAL_D normal,
3354 REAL_D res),
3355 REAL alpha_r,
3356 const WALL_QUAD *wall_quad)
3357 {
3358 return _AI_boundary_conditions_dow(matrix, fh, uh, bound,
3359 dirichlet_segment,
3360 g, gn,
3361 NULL, NULL, NULL, FILL_NOTHING,
3362 alpha_r, wall_quad);
3363 }
3364
boundary_conditions_loc_dow(DOF_MATRIX * matrix,DOF_REAL_VEC_D * fh,DOF_REAL_VEC_D * uh,DOF_SCHAR_VEC * bound,const BNDRY_FLAGS dirichlet_segment,LOC_FCT_D_AT_QP g_loc,LOC_FCT_D_AT_QP gn_loc,void * fd,FLAGS fill_flag,REAL alpha_r,const WALL_QUAD * wall_quad)3365 bool boundary_conditions_loc_dow(DOF_MATRIX *matrix,
3366 DOF_REAL_VEC_D *fh,
3367 DOF_REAL_VEC_D *uh,
3368 DOF_SCHAR_VEC *bound,
3369 const BNDRY_FLAGS dirichlet_segment,
3370 LOC_FCT_D_AT_QP g_loc,
3371 LOC_FCT_D_AT_QP gn_loc,
3372 void *fd, FLAGS fill_flag,
3373 REAL alpha_r,
3374 const WALL_QUAD *wall_quad)
3375 {
3376 return _AI_boundary_conditions_dow(matrix, fh, uh, bound,
3377 dirichlet_segment,
3378 NULL, NULL,
3379 g_loc, gn_loc, fd, fill_flag,
3380 alpha_r, wall_quad);
3381 }
3382
3383
3384