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: disc-ortho-poly.c
7 *
8 * description: Orthogonal discontinuous polynomials of degree 1 and 2.
9 *
10 *******************************************************************************
11 *
12 * authors: Claus-Justus Heine
13 * Abteilung fuer Angewandte Mathematik
14 * Albert-Ludwigs-Universitaet Freiburg
15 * Hermann-Herder-Str. 10
16 * D-79104 Freiburg im Breisgau, Germany
17 *
18 * http://www.alberta-fem.de
19 *
20 * (c) by C.-J. Heine (2008)
21 *
22 ******************************************************************************/
23
24 #include "alberta.h"
25
26 #define COMPILED_CONSTANTS 0
27
28 #if COMPILED_CONSTANTS
29 # define WEIGHT_V (2.0*sqrt(6.0))
30 # define WEIGHT_E (60.0/sqrt(35.0-10.0*sqrt(10.0)))
31 #else /* computed */
32 # define WEIGHT_V 4.898979485566356196394568149411782783932L
33 # define WEIGHT_E 32.64911064067351732799557417773087413491L
34 #endif
35
36 static ORTHO_DATA d_ortho_1_2d_data;
37 static ORTHO_DATA d_ortho_2_2d_data;
38
39 /* p.w. linear orthogonal polynomials */
d_ortho_phi_v_2d(const REAL_B lambda,int v)40 static inline REAL d_ortho_phi_v_2d(const REAL_B lambda, int v)
41 {
42 return WEIGHT_V * (lambda[v] - 0.5);
43 }
44
45 static inline const REAL *
d_ortho_grd_phi_v_2d(REAL_B result,const REAL_B lambda,int v,bool doinit)46 d_ortho_grd_phi_v_2d(REAL_B result, const REAL_B lambda, int v, bool doinit)
47 {
48 if (doinit) {
49 SET_BAR(2 /* dim */, 0.0, result);
50 }
51
52 result[v] = WEIGHT_V;
53
54 return result;
55 }
56
d_ortho_phi_v0_2d(const REAL_B lambda,const BAS_FCTS * thisptr)57 static REAL d_ortho_phi_v0_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
58 {
59 return d_ortho_phi_v_2d(lambda, 0);
60 }
61
d_ortho_phi_v1_2d(const REAL_B lambda,const BAS_FCTS * thisptr)62 static REAL d_ortho_phi_v1_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
63 {
64 return d_ortho_phi_v_2d(lambda, 1);
65 }
66
d_ortho_phi_v2_2d(const REAL_B lambda,const BAS_FCTS * thisptr)67 static REAL d_ortho_phi_v2_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
68 {
69 return d_ortho_phi_v_2d(lambda, 2);
70 }
71
72 static const REAL *
d_ortho_grd_phi_v0_2d(const REAL_B lambda,const BAS_FCTS * thisptr)73 d_ortho_grd_phi_v0_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
74 {
75 static REAL_B grd;
76
77 return d_ortho_grd_phi_v_2d(grd, lambda, 0, false);
78 }
79
80 static const REAL *
d_ortho_grd_phi_v1_2d(const REAL_B lambda,const BAS_FCTS * thisptr)81 d_ortho_grd_phi_v1_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
82 {
83 static REAL_B grd;
84
85 return d_ortho_grd_phi_v_2d(grd, lambda, 1, false);
86 }
87
88 static const REAL *
d_ortho_grd_phi_v2_2d(const REAL_B lambda,const BAS_FCTS * thisptr)89 d_ortho_grd_phi_v2_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
90 {
91 static REAL_B grd;
92
93 return d_ortho_grd_phi_v_2d(grd, lambda, 2, false);
94 }
95
d_ortho_D2_phi_v_2d(REAL_BB result,const REAL_B lambda,int v,bool doinit)96 static inline const REAL_B *d_ortho_D2_phi_v_2d(REAL_BB result,
97 const REAL_B lambda, int v,
98 bool doinit)
99 {
100 if (doinit) {
101 MSET_BAR(2 /* dim */, 0.0, result);
102 }
103
104 return (const REAL_B *)result;
105 }
106
107 static const REAL_B *
d_ortho_D2_phi_v0_2d(const REAL_B lambda,const BAS_FCTS * thisptr)108 d_ortho_D2_phi_v0_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
109 {
110 static REAL_BB D2;
111
112 return d_ortho_D2_phi_v_2d(D2, lambda, 0, false);
113 }
114
115 static const REAL_B *
d_ortho_D2_phi_v1_2d(const REAL_B lambda,const BAS_FCTS * thisptr)116 d_ortho_D2_phi_v1_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
117 {
118 static REAL_BB D2;
119
120 return d_ortho_D2_phi_v_2d(D2, lambda, 1, false);
121 }
122
123 static const REAL_B *
d_ortho_D2_phi_v2_2d(const REAL_B lambda,const BAS_FCTS * thisptr)124 d_ortho_D2_phi_v2_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
125 {
126 static REAL_BB D2;
127
128 return d_ortho_D2_phi_v_2d(D2, lambda, 2, false);
129 }
130
131 /*******************************************************************************
132 *
133 * p.w. quadratic orthogonal add-on
134 *
135 ******************************************************************************/
136 #if COMPILED_CONSTANTS
137 # define B (-3.0/2.0+0.5*sqrt(10.0)) /* (sqrt(46.0) / 120.0) */
138 # define D (sqrt(15.0) / 60.0)
139 # define E (-sqrt(6.0) / 24.0 + sqrt(15.0) / 30.0)
140 #else /* computed */
141 # define B 0.081138830084189665999446772216359266860L
142 # define D 0.06454972243679028141965442332970666018056L
143 # define E 0.0270373722576148087477553435466678456958L
144 #endif
145
d_ortho_phi_e_2d(const REAL_B lambda,int e)146 static inline REAL d_ortho_phi_e_2d(const REAL_B lambda, int e)
147 {
148 int i1 = (e + 1) % N_VERTICES_2D, i2 = (e + 2) % N_VERTICES_2D;
149
150 return
151 WEIGHT_E
152 *
153 (lambda[i1]*lambda[i2]
154 +
155 B*lambda[e]*(lambda[i1]+lambda[i2])
156 +
157 D*d_ortho_phi_v_2d(lambda, e)
158 +
159 E*(d_ortho_phi_v_2d(lambda, i1) + d_ortho_phi_v_2d(lambda, i2)));
160 }
161
162 static inline const REAL *
d_ortho_grd_phi_e_2d(REAL_B grd,const REAL_B lambda,int e)163 d_ortho_grd_phi_e_2d(REAL_B grd, const REAL_B lambda, int e)
164 {
165 REAL_B tmp;
166 int i1 = (e + 1) % N_VERTICES_2D, i2 = (e + 2) % N_VERTICES_2D;
167
168 grd[e] = B * (lambda[i1] + lambda[i2]);
169 grd[i1] = lambda[i2] + B * lambda[e];
170 grd[i2] = lambda[i1] + B * lambda[e];
171
172 AXPY_BAR(2 /* dim */, D, d_ortho_grd_phi_v_2d(tmp, lambda, e, true), grd);
173 AXPY_BAR(2 /* dim */, E, d_ortho_grd_phi_v_2d(tmp, lambda, i1, true), grd);
174 AXPY_BAR(2 /* dim */, E, d_ortho_grd_phi_v_2d(tmp, lambda, i2, true), grd);
175
176 SCAL_BAR(2 /* dim */, WEIGHT_E, grd);
177
178 return grd;
179 }
180
181 static inline const REAL_B *
d_ortho_D2_phi_e_2d(REAL_BB D2,const REAL_B lambda,int e)182 d_ortho_D2_phi_e_2d(REAL_BB D2, const REAL_B lambda, int e)
183 {
184 REAL_BB tmp = { { 0.0, } };
185 int i1 = (e + 1) % N_VERTICES_2D, i2 = (e + 2) % N_VERTICES_2D;
186
187 D2[e][i1] = D2[i1][e] = B;
188 D2[e][i2] = D2[i2][e] = B;
189 D2[i1][i2] = 1.0;
190 D2[i2][i1] = 1.0;
191
192 MAXPY_BAR(2 /* dim */, D, d_ortho_D2_phi_v_2d(tmp, lambda, e, true), D2);
193 MAXPY_BAR(2 /* dim */, E, d_ortho_D2_phi_v_2d(tmp, lambda, i1, true), D2);
194 MAXPY_BAR(2 /* dim */, E, d_ortho_D2_phi_v_2d(tmp, lambda, i2, true), D2);
195
196 MSCAL_BAR(2 /* dim */, WEIGHT_E, D2);
197
198 return (const REAL_B *)D2;
199 }
200
201 #undef B
202 #undef D
203 #undef E
204 #undef WEIGHT_E
205 #undef WEIGHT_V
206
d_ortho_phi_e0_2d(const REAL_B lambda,const BAS_FCTS * thisptr)207 static REAL d_ortho_phi_e0_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
208 {
209 return d_ortho_phi_e_2d(lambda, 0);
210 }
211
d_ortho_phi_e1_2d(const REAL_B lambda,const BAS_FCTS * thisptr)212 static REAL d_ortho_phi_e1_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
213 {
214 return d_ortho_phi_e_2d(lambda, 1);
215 }
216
d_ortho_phi_e2_2d(const REAL_B lambda,const BAS_FCTS * thisptr)217 static REAL d_ortho_phi_e2_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
218 {
219 return d_ortho_phi_e_2d(lambda, 2);
220 }
221
222 static const REAL *
d_ortho_grd_phi_e0_2d(const REAL_B lambda,const BAS_FCTS * thisptr)223 d_ortho_grd_phi_e0_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
224 {
225 static REAL_B grd;
226
227 return d_ortho_grd_phi_e_2d(grd, lambda, 0);
228 }
229
230 static const REAL *
d_ortho_grd_phi_e1_2d(const REAL_B lambda,const BAS_FCTS * thisptr)231 d_ortho_grd_phi_e1_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
232 {
233 static REAL_B grd;
234
235 return d_ortho_grd_phi_e_2d(grd, lambda, 1);
236 }
237
238 static const REAL *
d_ortho_grd_phi_e2_2d(const REAL_B lambda,const BAS_FCTS * thisptr)239 d_ortho_grd_phi_e2_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
240 {
241 static REAL_B grd;
242
243 return d_ortho_grd_phi_e_2d(grd, lambda, 2);
244 }
245
246 static const REAL_B *
d_ortho_D2_phi_e0_2d(const REAL_B lambda,const BAS_FCTS * thisptr)247 d_ortho_D2_phi_e0_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
248 {
249 static REAL_BB D2;
250
251 return d_ortho_D2_phi_e_2d(D2, lambda, 0);
252 }
253
254 static const REAL_B *
d_ortho_D2_phi_e1_2d(const REAL_B lambda,const BAS_FCTS * thisptr)255 d_ortho_D2_phi_e1_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
256 {
257 static REAL_BB D2;
258
259 return d_ortho_D2_phi_e_2d(D2, lambda, 1);
260 }
261
262 static const REAL_B *
d_ortho_D2_phi_e2_2d(const REAL_B lambda,const BAS_FCTS * thisptr)263 d_ortho_D2_phi_e2_2d(const REAL_B lambda, const BAS_FCTS *thisptr)
264 {
265 static REAL_BB D2;
266
267 return d_ortho_D2_phi_e_2d(D2, lambda, 2);
268 }
269
270 static const BAS_FCT d_ortho_phi_2d[N_BAS_LAG_2_2D] = {
271 d_ortho_phi_v0_2d, d_ortho_phi_v1_2d, d_ortho_phi_v2_2d,
272 d_ortho_phi_e0_2d, d_ortho_phi_e1_2d, d_ortho_phi_e2_2d
273 };
274
275 static const GRD_BAS_FCT d_ortho_grd_phi_2d[N_BAS_LAG_2_2D] = {
276 d_ortho_grd_phi_v0_2d, d_ortho_grd_phi_v1_2d, d_ortho_grd_phi_v2_2d,
277 d_ortho_grd_phi_e0_2d, d_ortho_grd_phi_e1_2d, d_ortho_grd_phi_e2_2d
278 };
279
280 static const D2_BAS_FCT d_ortho_D2_phi_2d[N_BAS_LAG_2_2D] = {
281 d_ortho_D2_phi_v0_2d, d_ortho_D2_phi_v1_2d, d_ortho_D2_phi_v2_2d,
282 d_ortho_D2_phi_e0_2d, d_ortho_D2_phi_e1_2d, d_ortho_D2_phi_e2_2d
283 };
284
285 /******************************************************************************
286 *
287 * Degree 1
288 *
289 ******************************************************************************/
290
291 #undef DEF_EL_VEC_D_ORTHO_1_2D
292 #define DEF_EL_VEC_D_ORTHO_1_2D(type, name) \
293 DEF_EL_VEC_CONST(type, name, N_BAS_LAG_1_2D, N_BAS_LAG_1_2D)
294
295 #undef DEFUN_GET_EL_VEC_D_ORTHO_1_2D
296 #define DEFUN_GET_EL_VEC_D_ORTHO_1_2D(name, type, admin, body, ...) \
297 static const EL_##type##_VEC * \
298 d_ortho_get_##name##_1_2d(type##_VEC_TYPE *vec, const EL *el, __VA_ARGS__) \
299 { \
300 FUNCNAME("d_get_"#name"d_1_2d"); \
301 static DEF_EL_VEC_D_ORTHO_1_2D(type, rvec_space); \
302 type##_VEC_TYPE *rvec = vec ? vec : rvec_space->vec; \
303 int n0, node, ibas; \
304 DOF **dofptr = el->dof, dof; \
305 \
306 DEBUG_TEST_EXIT(true, ""); \
307 \
308 node = (admin)->mesh->node[CENTER]; \
309 n0 = (admin)->n0_dof[CENTER]; \
310 for (ibas = 0; ibas < N_BAS_LAG_1_2D; ibas++) { \
311 dof = dofptr[node][n0+ibas]; \
312 body; \
313 } \
314 \
315 return vec ? NULL : rvec_space; \
316 } \
317 struct _AI_semicolon_dummy
318
319 #undef DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D
320 #define DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(name, type, ASSIGN) \
321 DEFUN_GET_EL_VEC_D_ORTHO_1_2D(_##name##_vec, type, dv->fe_space->admin, \
322 ASSIGN(dv->vec[dof], rvec[ibas]), \
323 const DOF_##type##_VEC *dv); \
324 static const EL_##type##_VEC * \
325 d_ortho_get_##name##_vec_1_2d(type##_VEC_TYPE *vec, const EL *el, \
326 const DOF_##type##_VEC *dv) \
327 { \
328 EL_##type##_VEC *vec_loc = dv->vec_loc; \
329 \
330 if (vec != NULL || vec_loc == NULL) { \
331 return d_ortho_get__##name##_vec_1_2d(vec, el, dv); \
332 } else { \
333 d_ortho_get__##name##_vec_1_2d(vec_loc->vec, el, dv); \
334 return vec_loc; \
335 } \
336 } \
337 struct _AI_semicolon_dummy
338
339 #undef COPY_EQ
340 #define COPY_EQ(a, b) (b) = (a)
341
342 /*--------------------------------------------------------------------------*/
343 /* functions for combining basisfunctions with coefficients */
344 /*--------------------------------------------------------------------------*/
345
346 DEFUN_GET_EL_VEC_D_ORTHO_1_2D(dof_indices, DOF, admin,
347 rvec[ibas] = dof,
348 const DOF_ADMIN *admin, const BAS_FCTS *thisptr);
349
d_ortho_get_bound_1_2d(BNDRY_FLAGS * vec,const EL_INFO * el_info,const BAS_FCTS * thisptr)350 static const EL_BNDRY_VEC *d_ortho_get_bound_1_2d(BNDRY_FLAGS *vec,
351 const EL_INFO *el_info,
352 const BAS_FCTS *thisptr)
353 {
354 FUNCNAME("d_get_bound2_2d");
355 static DEF_EL_VEC_D_ORTHO_1_2D(BNDRY, rvec_space);
356 BNDRY_FLAGS *rvec = vec ? vec : rvec_space->vec;
357 int i;
358
359 TEST_FLAG(FILL_BOUND, el_info);
360
361 for (i = 0; i < N_BAS_LAG_1_2D; i++) {
362 BNDRY_FLAGS_INIT(rvec[i]);
363 BNDRY_FLAGS_SET(rvec[i], el_info->face_bound[0]);
364 }
365
366 return vec ? NULL : rvec_space;
367 }
368
369 /*--------------------------------------------------------------------*/
370 /*--- function for accessing a local DOF_INT_VEC ---*/
371 /*--------------------------------------------------------------------*/
372
373 DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(int, INT, COPY_EQ);
374
375 /*--------------------------------------------------------------------*/
376 /*--- function for accessing a local DOF_REAL_VEC ---*/
377 /*--------------------------------------------------------------------*/
378
379 DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(real, REAL, COPY_EQ);
380
381 /*--------------------------------------------------------------------*/
382 /*--- function for accessing a local DOF_REAL_D_VEC ---*/
383 /*--------------------------------------------------------------------*/
384
385 DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(real_d, REAL_D, COPY_DOW);
386
387 /*--------------------------------------------------------------------*/
388 /*--- function for accessing a local DOF_SCHAR_VEC ---*/
389 /*--------------------------------------------------------------------*/
390
391 DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(schar, SCHAR, COPY_EQ);
392
393 /*--------------------------------------------------------------------*/
394 /*--- function for accessing a local DOF_UCHAR_VEC ---*/
395 /*--------------------------------------------------------------------*/
396
397 DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(uchar, UCHAR, COPY_EQ);
398
399 /*--------------------------------------------------------------------*/
400 /*--- function for accessing a local DOF_PTR_VEC ---*/
401 /*--------------------------------------------------------------------*/
402
403 DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(ptr, PTR, COPY_EQ);
404
405 /*--------------------------------------------------------------------*/
406 /*--- function for accessing a local DOF_REAL_DD_VEC ---*/
407 /*--------------------------------------------------------------------*/
408
409 DEFUN_GET_EL_DOF_VEC_D_ORTHO_1_2D(real_dd, REAL_DD, _AI_MCOPY_DOW);
410
411 /****************************************************************************/
412 /* functions for interpolation/ restriction during refinement/coarsening */
413 /****************************************************************************/
414
415 static void
d_ortho_real_refine_inter_1_2d(DOF_REAL_VEC * drv,RC_LIST_EL * list,int n)416 d_ortho_real_refine_inter_1_2d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
417 {
418 DOF pdof[N_VERTICES_2D];
419 DOF cdof;
420 EL *child;
421 int i, j, n0, node;
422
423 n0 = drv->fe_space->admin->n0_dof[CENTER];
424 node = drv->fe_space->admin->mesh->node[CENTER];
425
426 for (i = 0; i < n; i++) {
427 EL *el = list[i].el_info.el;
428 for (j = 0; j < N_VERTICES_2D; j++) {
429 pdof[j] = el->dof[node][n0+j];
430 }
431
432 /* child0
433 * phi0 = 1/2*phi1_c - 1/2*phi0_c;
434 * phi1 = phi2_c + 1/2*(phi0_c + phi1_c);
435 * phi2 = phi0_c;
436 */
437 child = list[i].el_info.el->child[0];
438
439 cdof = child->dof[node][n0+0];
440 drv->vec[cdof] =
441 -0.5*drv->vec[pdof[0]] + 0.5*drv->vec[pdof[1]] + drv->vec[pdof[2]];
442
443 cdof = child->dof[node][n0+1];
444 drv->vec[cdof] = 0.5*drv->vec[pdof[0]] + 0.5*drv->vec[pdof[1]];
445
446 cdof = child->dof[node][n0+2];
447 drv->vec[cdof] = drv->vec[pdof[1]];
448
449 /* child1
450 * phi0 = phi2_c + 1/2*(phi0_c + phi1_c);
451 * phi1 = 1/2*phi0_c -1/2*phi1_c;
452 * phi2 = phi1_c;
453 */
454 child = list[i].el_info.el->child[1];
455
456 cdof = child->dof[node][n0+0];
457 drv->vec[cdof] = 0.5*drv->vec[pdof[0]] + 0.5*drv->vec[pdof[1]];
458
459 cdof = child->dof[node][n0+1];
460 drv->vec[cdof] =
461 0.5*drv->vec[pdof[0]] - 0.5*drv->vec[pdof[1]] + drv->vec[pdof[2]];
462
463 cdof = child->dof[node][n0+2];
464 drv->vec[cdof] = drv->vec[pdof[0]];
465 }
466 }
467
468 static void
d_ortho_real_coarse_inter_1_2d(DOF_REAL_VEC * drv,RC_LIST_EL * list,int n)469 d_ortho_real_coarse_inter_1_2d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
470 {
471 DOF pdof[N_VERTICES_2D];
472 DOF cdof;
473 EL *child;
474 int i, j, n0, node;
475
476 n0 = drv->fe_space->admin->n0_dof[CENTER];
477 node = drv->fe_space->admin->mesh->node[CENTER];
478
479 for (i = 0; i < n; i++) {
480 EL *el = list[i].el_info.el;
481 for (j = 0; j < N_VERTICES_2D; j++) {
482 pdof[j] = el->dof[node][n0+j];
483 drv->vec[pdof[j]] = 0.0;
484 }
485
486 /* child0:
487 * phi0.phi0_c = 1/2 (phi1_c - phi0_c).phi0_c = -1/2
488 * phi0.phi1_c = ... = +1/2
489 * phi0.phi2_c = = 0
490 *
491 * phi1.phi0_c = (phi2_c + 1/2*(phi0_c + phi1_c)).phi0_c = 1/2
492 * phi1.phi1_c = (phi2_c + 1/2*(phi0_c + phi1_c)).phi1_c = 1/2
493 * phi1.phi2_c = (phi2_c + 1/2*(phi0_c + phi1_c)).phi2_c = 1
494 *
495 * phi2.phi0_c = phi0_c.phi0_c = 1
496 * phi2.phi1_c = 0
497 * phi2.phi2_c = 0
498 *
499 * child1:
500 * phi0.phi0_c = (phi2_c + 1/2*(phi0_c + phi1_c)).phi0_c = 1/2
501 * phi0.phi1_c = (phi2_c + 1/2*(phi0_c + phi1_c)).phi1_c = 1/2
502 * phi0.phi2_c = (phi2_c + 1/2*(phi0_c + phi1_c)).phi2_c = 1
503 *
504 * phi1.phi0_c = (1/2*phi0_c -1/2*phi1_c).phi0_c = +1/2
505 * phi1.phi1_c = (1/2*phi0_c -1/2*phi1_c).phi1_c = -1/2
506 * phi1.phi2_c = (1/2*phi0_c -1/2*phi1_c).phi2_c = 0
507 *
508 * phi2.phi0_c = phi1_c.phi0_c = 0
509 * phi2.phi1_c = phi1_c.phi1_c = 1
510 * phi2.phi2_c = phi1_c.phi2_c = 0
511 */
512
513 child = list[i].el_info.el->child[0];
514
515 cdof = child->dof[node][n0+0];
516 drv->vec[pdof[0]] += 0.5 * drv->vec[cdof] * (-0.5);
517 drv->vec[pdof[1]] += 0.5 * drv->vec[cdof] * (+0.5);
518 drv->vec[pdof[2]] += 0.5 * drv->vec[cdof] * (+1.0);
519 cdof = child->dof[node][n0+1];
520 drv->vec[pdof[0]] += 0.5 * drv->vec[cdof] * (+0.5);
521 drv->vec[pdof[1]] += 0.5 * drv->vec[cdof] * (+0.5);
522 cdof = child->dof[node][n0+2];
523 drv->vec[pdof[1]] += 0.5 * drv->vec[cdof] * (+1.0);
524
525 child = list[i].el_info.el->child[1];
526 cdof = child->dof[node][n0+0];
527 drv->vec[pdof[0]] += 0.5 * drv->vec[cdof] * (+0.5);
528 drv->vec[pdof[1]] += 0.5 * drv->vec[cdof] * (+0.5);
529 cdof = child->dof[node][n0+1];
530 drv->vec[pdof[0]] += 0.5 * drv->vec[cdof] * (+0.5);
531 drv->vec[pdof[1]] += 0.5 * drv->vec[cdof] * (-0.5);
532 drv->vec[pdof[2]] += 0.5 * drv->vec[cdof] * (+1.0);
533 cdof = child->dof[node][n0+2];
534 drv->vec[pdof[0]] += 0.5 * drv->vec[cdof] * (+1.0);
535 }
536 }
537
538 /* Interpolation: we always use L2-interpolation over the entire
539 * element, even if WALL >= 0. The quadrature degree is such that we
540 * would reproduce the polynomials of the given degree.
541 */
542 static void
d_ortho_interpol_1_2d(EL_REAL_VEC * el_vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_AT_QP f,void * f_data,const BAS_FCTS * thisptr)543 d_ortho_interpol_1_2d(EL_REAL_VEC *el_vec,
544 const EL_INFO *el_info, int wall,
545 int no, const int *b_no,
546 LOC_FCT_AT_QP f, void *f_data,
547 const BAS_FCTS *thisptr)
548 {
549 const QUAD_FAST *qfast = ((ORTHO_DATA *)thisptr->ext_data)->qfast;
550
551 if (b_no) {
552 int i, iq;
553
554 for (i = 0; i < no; i++) {
555 el_vec->vec[b_no[i]] = 0.0;
556 }
557 for (iq = 0; iq < qfast->n_points; iq++) {
558 REAL value = qfast->w[iq] * f(el_info, qfast->quad, iq, f_data);
559 for (i = 0; i < no; i++) {
560 int ib = b_no[i];
561 el_vec->vec[ib] += qfast->phi[iq][ib] * value;
562 }
563 }
564 } else {
565 int iq, ib;
566
567 for (ib = 0; ib < N_BAS_LAG_1_2D; ib++) {
568 el_vec->vec[ib] = 0.0;
569 }
570 for (iq = 0; iq < qfast->n_points; iq++) {
571 REAL value = qfast->w[iq] * f(el_info, qfast->quad, iq, f_data);
572 for (ib = 0; ib < N_BAS_LAG_1_2D; ib++) {
573 el_vec->vec[ib] += qfast->phi[iq][ib] * value;
574 }
575 }
576 }
577 }
578
579 static void
d_ortho_interpol_d_1_2d(EL_REAL_D_VEC * el_vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_D_AT_QP f,void * f_data,const BAS_FCTS * thisptr)580 d_ortho_interpol_d_1_2d(EL_REAL_D_VEC *el_vec,
581 const EL_INFO *el_info, int wall,
582 int no, const int *b_no,
583 LOC_FCT_D_AT_QP f, void *f_data,
584 const BAS_FCTS *thisptr)
585 {
586 const QUAD_FAST *qfast = ((ORTHO_DATA *)thisptr->ext_data)->qfast;
587
588 if (b_no) {
589 int i, iq;
590
591 for (i = 0; i < no; i++) {
592 SET_DOW(0.0, el_vec->vec[b_no[i]]);
593 }
594 for (iq = 0; iq < qfast->n_points; iq++) {
595 REAL_D value;
596
597 f(value, el_info, qfast->quad, iq, f_data);
598 SCAL_DOW(qfast->w[iq], value);
599 for (i = 0; i < no; i++) {
600 int ib = b_no[i];
601 AXPY_DOW( qfast->phi[iq][ib], value, el_vec->vec[ib]);
602 }
603 }
604 } else {
605 int iq, ib;
606
607 for (ib = 0; ib < N_BAS_LAG_1_2D; ib++) {
608 SET_DOW(0.0, el_vec->vec[ib]);
609 }
610 for (iq = 0; iq < qfast->n_points; iq++) {
611 REAL_D value;
612
613 f(value, el_info, qfast->quad, iq, f_data);
614 SCAL_DOW(qfast->w[iq], value);
615 for (ib = 0; ib < N_BAS_LAG_1_2D; ib++) {
616 AXPY_DOW(qfast->phi[iq][ib], value, el_vec->vec[ib]);
617 }
618 }
619 }
620 }
621
622 static void
d_ortho_interpol_dow_1_2d(EL_REAL_VEC_D * el_vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_D_AT_QP f,void * f_data,const BAS_FCTS * thisptr)623 d_ortho_interpol_dow_1_2d(EL_REAL_VEC_D *el_vec,
624 const EL_INFO *el_info, int wall,
625 int no, const int *b_no,
626 LOC_FCT_D_AT_QP f, void *f_data,
627 const BAS_FCTS *thisptr)
628 {
629 d_ortho_interpol_d_1_2d((EL_REAL_D_VEC *)el_vec,
630 el_info, wall, no, b_no, f, f_data, thisptr);
631 }
632
633 /* The trace-space is degenerated: all basis functions have non-zero
634 * trace on all walls. If permute the stuff cyclic s.t. we would
635 * eventually have the chance to really define the trace-space as
636 * degenerated set of basis functions.
637 */
638 static const int trace_mapping_d_ortho_1_2d[N_WALLS_2D][N_BAS_LAG_1_2D] = {
639 { 1, 2, 0 }, { 2, 0, 1 }, { 0, 1, 2 }
640 };
641
642 static const BAS_FCTS disc_ortho1_2d = {
643 DISC_ORTHO_NAME(1)"_2d", 2, 1, N_BAS_LAG_1_2D, N_BAS_LAG_1_2D, 1,
644 {0, N_BAS_LAG_1_2D, 0, 0},
645 -1, /* trace_admin */
646 INIT_BFCTS_CHAIN(disc_ortho1_2d),
647 INIT_ELEMENT_INITIALIZER(NULL, FILL_NOTHING), /* init_element + fill-flags */
648 d_ortho_phi_2d, d_ortho_grd_phi_2d, d_ortho_D2_phi_2d,
649 NULL, NULL, /* third and fourth derivatives */
650 NULL, NULL, NULL, false, /* phi_d etc. */
651 /********************/
652 NULL, /* trace space */
653 { { { trace_mapping_d_ortho_1_2d[0],
654 trace_mapping_d_ortho_1_2d[1],
655 trace_mapping_d_ortho_1_2d[2], }, }, }, /* trace mapping */
656 { N_BAS_LAG_1_2D,
657 N_BAS_LAG_1_2D,
658 N_BAS_LAG_1_2D, }, /* n_trace_bas_fcts */
659 d_ortho_get_dof_indices_1_2d,
660 d_ortho_get_bound_1_2d,
661 d_ortho_interpol_1_2d,
662 d_ortho_interpol_d_1_2d,
663 d_ortho_interpol_dow_1_2d,
664 d_ortho_get_int_vec_1_2d,
665 d_ortho_get_real_vec_1_2d,
666 d_ortho_get_real_d_vec_1_2d,
667 (GET_REAL_VEC_D_TYPE)d_ortho_get_real_d_vec_1_2d,
668 d_ortho_get_uchar_vec_1_2d,
669 d_ortho_get_schar_vec_1_2d,
670 d_ortho_get_ptr_vec_1_2d,
671 d_ortho_get_real_dd_vec_1_2d,
672 d_ortho_real_refine_inter_1_2d,
673 d_ortho_real_coarse_inter_1_2d,
674 NULL,
675 NULL,
676 NULL,
677 NULL,
678 NULL,
679 NULL,
680 NULL,
681 (void *)&d_ortho_1_2d_data
682 };
683
684 /******************************************************************************
685 *
686 * Degree 2
687 *
688 ******************************************************************************/
689
690 #undef DEF_EL_VEC_D_ORTHO_2_2D
691 #define DEF_EL_VEC_D_ORTHO_2_2D(type, name) \
692 DEF_EL_VEC_CONST(type, name, N_BAS_LAG_2_2D, N_BAS_LAG_2_2D)
693
694 #undef DEFUN_GET_EL_VEC_D_ORTHO_2_2D
695 #define DEFUN_GET_EL_VEC_D_ORTHO_2_2D(name, type, admin, body, ...) \
696 static const EL_##type##_VEC * \
697 d_ortho_get_##name##_2_2d(type##_VEC_TYPE *vec, const EL *el, __VA_ARGS__) \
698 { \
699 FUNCNAME("d_get_"#name"d_2_2d"); \
700 static DEF_EL_VEC_D_ORTHO_2_2D(type, rvec_space); \
701 type##_VEC_TYPE *rvec = vec ? vec : rvec_space->vec; \
702 int n0, node, ibas; \
703 DOF **dofptr = el->dof, dof; \
704 \
705 DEBUG_TEST_EXIT(true, ""); \
706 \
707 node = (admin)->mesh->node[CENTER]; \
708 n0 = (admin)->n0_dof[CENTER]; \
709 for (ibas = 0; ibas < N_BAS_LAG_2_2D; ibas++) { \
710 dof = dofptr[node][n0+ibas]; \
711 body; \
712 } \
713 \
714 return vec ? NULL : rvec_space; \
715 } \
716 struct _AI_semicolon_dummy
717
718 #undef DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D
719 #define DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(name, type, ASSIGN) \
720 DEFUN_GET_EL_VEC_D_ORTHO_2_2D(_##name##_vec, type, dv->fe_space->admin, \
721 ASSIGN(dv->vec[dof], rvec[ibas]), \
722 const DOF_##type##_VEC *dv); \
723 static const EL_##type##_VEC * \
724 d_ortho_get_##name##_vec_2_2d(type##_VEC_TYPE *vec, const EL *el, \
725 const DOF_##type##_VEC *dv) \
726 { \
727 EL_##type##_VEC *vec_loc = dv->vec_loc; \
728 \
729 if (vec != NULL || vec_loc == NULL) { \
730 return d_ortho_get__##name##_vec_2_2d(vec, el, dv); \
731 } else { \
732 d_ortho_get__##name##_vec_2_2d(vec_loc->vec, el, dv); \
733 return vec_loc; \
734 } \
735 } \
736 struct _AI_semicolon_dummy
737
738 #undef COPY_EQ
739 #define COPY_EQ(a, b) (b) = (a)
740
741 /*--------------------------------------------------------------------------*/
742 /* functions for combining basisfunctions with coefficients */
743 /*--------------------------------------------------------------------------*/
744
745 DEFUN_GET_EL_VEC_D_ORTHO_2_2D(dof_indices, DOF, admin,
746 rvec[ibas] = dof,
747 const DOF_ADMIN *admin, const BAS_FCTS *thisptr);
748
d_ortho_get_bound_2_2d(BNDRY_FLAGS * vec,const EL_INFO * el_info,const BAS_FCTS * thisptr)749 static const EL_BNDRY_VEC *d_ortho_get_bound_2_2d(BNDRY_FLAGS *vec,
750 const EL_INFO *el_info,
751 const BAS_FCTS *thisptr)
752 {
753 FUNCNAME("d_get_bound2_2d");
754 static DEF_EL_VEC_D_ORTHO_2_2D(BNDRY, rvec_space);
755 BNDRY_FLAGS *rvec = vec ? vec : rvec_space->vec;
756 int i;
757
758 TEST_FLAG(FILL_BOUND, el_info);
759
760 for (i = 0; i < N_BAS_LAG_2_2D; i++) {
761 BNDRY_FLAGS_INIT(rvec[i]);
762 BNDRY_FLAGS_SET(rvec[i], el_info->face_bound[0]);
763 }
764
765 return vec ? NULL : rvec_space;
766 }
767
768 /*--------------------------------------------------------------------*/
769 /*--- function for accessing a local DOF_INT_VEC ---*/
770 /*--------------------------------------------------------------------*/
771
772 DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(int, INT, COPY_EQ);
773
774 /*--------------------------------------------------------------------*/
775 /*--- function for accessing a local DOF_REAL_VEC ---*/
776 /*--------------------------------------------------------------------*/
777
778 DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(real, REAL, COPY_EQ);
779
780 /*--------------------------------------------------------------------*/
781 /*--- function for accessing a local DOF_REAL_D_VEC ---*/
782 /*--------------------------------------------------------------------*/
783
784 DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(real_d, REAL_D, COPY_DOW);
785
786 /*--------------------------------------------------------------------*/
787 /*--- function for accessing a local DOF_SCHAR_VEC ---*/
788 /*--------------------------------------------------------------------*/
789
790 DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(schar, SCHAR, COPY_EQ);
791
792 /*--------------------------------------------------------------------*/
793 /*--- function for accessing a local DOF_UCHAR_VEC ---*/
794 /*--------------------------------------------------------------------*/
795
796 DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(uchar, UCHAR, COPY_EQ);
797
798 /*--------------------------------------------------------------------*/
799 /*--- function for accessing a local DOF_PTR_VEC ---*/
800 /*--------------------------------------------------------------------*/
801
802 DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(ptr, PTR, COPY_EQ);
803
804 /*--------------------------------------------------------------------*/
805 /*--- function for accessing a local DOF_REAL_DD_VEC ---*/
806 /*--------------------------------------------------------------------*/
807
808 DEFUN_GET_EL_DOF_VEC_D_ORTHO_2_2D(real_dd, REAL_DD, _AI_MCOPY_DOW);
809
810 /****************************************************************************/
811 /* functions for interpolation/ restriction during refinement/coarsening */
812 /****************************************************************************/
813
814 static void
d_ortho_real_refine_inter_2_2d(DOF_REAL_VEC * drv,RC_LIST_EL * list,int n)815 d_ortho_real_refine_inter_2_2d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
816 {
817 REAL pcoeff[N_BAS_LAG_2_2D];
818 EL *ch0, *ch1;
819 DOF cdof;
820 int i, j, n0, node;
821
822 n0 = drv->fe_space->admin->n0_dof[CENTER];
823 node = drv->fe_space->admin->mesh->node[CENTER];
824
825 for (i = 0; i < n; i++) {
826 EL *el = list[i].el_info.el;
827 for (j = 0; j < N_BAS_LAG_2_2D; j++) {
828 pcoeff[j] = drv->vec[el->dof[node][n0+j]];
829 }
830
831 /* Define some "magic" numbers ... */
832 #define SQRT2 sqrt(2.0)
833 #define SQRT6 sqrt(6.0)
834 #define SQRT10 sqrt(10.0)
835 #define SQRT15 sqrt(15.0)
836
837 #define DENOM sqrt(-10.0*SQRT10+35.0)
838
839 #define COEFF11 (1.0/8.0*(5.0*SQRT6-2.0*SQRT15) / DENOM)
840 #define COEFF12 (11.0*SQRT6-6.0*SQRT15) / (8.0*DENOM)
841 #define COEFF13 (SQRT6-2.0*SQRT15)/(8.0*DENOM)
842 #define COEFF14 (SQRT6/(2.0*DENOM))
843
844 /* linear sub-space */
845 ch0 = list[i].el_info.el->child[0];
846 ch1 = list[i].el_info.el->child[1];
847
848 cdof = ch0->dof[node][n0+0];
849 drv->vec[cdof] =
850 0.5*(-pcoeff[0] + pcoeff[1]) + pcoeff[2] + COEFF11*(-pcoeff[3]+pcoeff[4]);
851 cdof = ch1->dof[node][n0+1];
852 drv->vec[cdof] =
853 0.5*(+pcoeff[0] - pcoeff[1]) + pcoeff[2] + COEFF11*(+pcoeff[3]-pcoeff[4]);
854
855 cdof = ch0->dof[node][n0+1];
856 drv->vec[cdof] =
857 0.5*(pcoeff[0]+pcoeff[1])
858 +
859 COEFF12*pcoeff[3] + COEFF13*pcoeff[4] - COEFF14*pcoeff[5];
860 cdof = ch1->dof[node][n0+0];
861 drv->vec[cdof] =
862 0.5*(pcoeff[0]+pcoeff[1])
863 +
864 COEFF13*pcoeff[3] + COEFF12*pcoeff[4] - COEFF14*pcoeff[5];
865
866 cdof = ch0->dof[node][n0+2];
867 drv->vec[cdof] =
868 pcoeff[1]
869 - COEFF13*pcoeff[3] - COEFF12*pcoeff[4] + COEFF14*pcoeff[5];
870 cdof = ch1->dof[node][n0+2];
871 drv->vec[cdof] =
872 pcoeff[0]
873 - COEFF12*pcoeff[3] - COEFF13*pcoeff[4] + COEFF14*pcoeff[5];
874
875 /* quadratic add-on */
876 #define DENOM2 (2.0*SQRT10-7.0)
877 cdof = ch0->dof[node][n0+3];
878 drv->vec[cdof] =
879 0.125*(-22.0+7.0*SQRT10)/DENOM2*pcoeff[3]
880 +
881 0.125*(-3.0*SQRT10+10.0)/DENOM2*pcoeff[4]
882 +
883 0.125*(-52.0+16.0*SQRT10)/DENOM2*pcoeff[5];
884 cdof = ch1->dof[node][n0+4];
885 drv->vec[cdof] =
886 0.125*(-22.0+7.0*SQRT10)/DENOM2*pcoeff[4]
887 +
888 0.125*(-3.0*SQRT10+10.0)/DENOM2*pcoeff[3]
889 +
890 0.125*(-52.0+16.0*SQRT10)/DENOM2*pcoeff[5];
891
892 cdof = ch0->dof[node][n0+4];
893 drv->vec[cdof] =
894 1.0/40.0*(43.0*SQRT10-150.0)/DENOM2*pcoeff[3]
895 +
896 1.0/40.0*(-7.0*SQRT10+10.0)/DENOM2*pcoeff[4]
897 +
898 1.0/40.0*(-20.0+8.0*SQRT10)/DENOM2*pcoeff[5];
899 cdof = ch1->dof[node][n0+3];
900 drv->vec[cdof] =
901 1.0/40.0*(43.0*SQRT10-150.0)/DENOM2*pcoeff[4]
902 +
903 1.0/40.0*(-7.0*SQRT10+10.0)/DENOM2*pcoeff[3]
904 +
905 1.0/40.0*(-20.0+8.0*SQRT10)/DENOM2*pcoeff[5];
906
907 cdof = ch0->dof[node][n0+5];
908 drv->vec[cdof] =
909 -1.0/40.0*(13.0*SQRT10-40.0)/DENOM2*pcoeff[3]
910 -1.0/40.0*(-17.0*SQRT10+80.0)/DENOM2*pcoeff[4]
911 -1.0/40.0*(-12.0*SQRT10+40.0)/DENOM2*pcoeff[5];
912 cdof = ch1->dof[node][n0+5];
913 drv->vec[cdof] =
914 -1.0/40.0*(13.0*SQRT10-40.0)/DENOM2*pcoeff[4]
915 -1.0/40.0*(-17.0*SQRT10+80.0)/DENOM2*pcoeff[3]
916 -1.0/40.0*(-12.0*SQRT10+40.0)/DENOM2*pcoeff[5];
917 }
918 }
919
920 static void
d_ortho_real_coarse_inter_2_2d(DOF_REAL_VEC * drv,RC_LIST_EL * list,int n)921 d_ortho_real_coarse_inter_2_2d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
922 {
923 DOF pdof[N_EDGES_2D];
924 DOF cdof;
925 EL *child;
926 int i, j, n0, node;
927
928 /* The linear sub-space restricts to the coarser grid just like in
929 * the linear case.
930 */
931 d_ortho_real_coarse_inter_1_2d(drv, list, n);
932
933 n0 = drv->fe_space->admin->n0_dof[CENTER];
934 node = drv->fe_space->admin->mesh->node[CENTER];
935
936 for (i = 0; i < n; i++) {
937 EL *el = list[i].el_info.el;
938 for (j = 0; j < N_EDGES_2D; j++) {
939 pdof[j] = el->dof[node][n0+j+N_VERTICES_2D];
940 drv->vec[pdof[j]] = 0.0;
941 }
942
943 child = list[i].el_info.el->child[0];
944
945 cdof = child->dof[node][n0+0];
946 drv->vec[pdof[0]] += -COEFF11 * drv->vec[cdof];
947 drv->vec[pdof[1]] += +COEFF11 * drv->vec[cdof];
948
949 cdof = child->dof[node][n0+1];
950 drv->vec[pdof[0]] += +COEFF12 * drv->vec[cdof];
951 drv->vec[pdof[1]] += +COEFF13 * drv->vec[cdof];
952 drv->vec[pdof[2]] += -COEFF14 * drv->vec[cdof];
953
954 cdof = child->dof[node][n0+2];
955 drv->vec[pdof[0]] += -COEFF13 * drv->vec[cdof];
956 drv->vec[pdof[1]] += -COEFF12 * drv->vec[cdof];
957 drv->vec[pdof[2]] += +COEFF14 * drv->vec[cdof];
958
959 cdof = child->dof[node][n0+3];
960 drv->vec[pdof[0]] += 0.125*(-22.0+7.0*SQRT10)/DENOM2 * drv->vec[cdof];
961 drv->vec[pdof[1]] += 0.125*(-3.0*SQRT10+10.0)/DENOM2 * drv->vec[cdof];
962 drv->vec[pdof[2]] += 0.125*(-52.0+16.0*SQRT10)/DENOM2 * drv->vec[cdof];
963
964 cdof = child->dof[node][n0+4];
965 drv->vec[pdof[0]] += 1.0/40.0*(43.0*SQRT10-150.0)/DENOM2 * drv->vec[cdof];
966 drv->vec[pdof[1]] += 1.0/40.0*(-7.0*SQRT10+10.0)/DENOM2 * drv->vec[cdof];
967 drv->vec[pdof[2]] += 1.0/40.0*(-20.0+8.0*SQRT10)/DENOM2 * drv->vec[cdof];
968
969 cdof = child->dof[node][n0+5];
970 drv->vec[pdof[0]] += -1.0/40.0*(13.0*SQRT10-40.0)/DENOM2 * drv->vec[cdof];
971 drv->vec[pdof[1]] += -1.0/40.0*(-17.0*SQRT10+80.0)/DENOM2 * drv->vec[cdof];
972 drv->vec[pdof[2]] += -1.0/40.0*(-12.0*SQRT10+40.0)/DENOM2 * drv->vec[cdof];
973
974 /* Second child*/
975 child = list[i].el_info.el->child[1];
976
977 cdof = child->dof[node][n0+1];
978 drv->vec[pdof[0]] += +COEFF11 * drv->vec[cdof];
979 drv->vec[pdof[1]] += -COEFF11 * drv->vec[cdof];
980
981 cdof = child->dof[node][n0+0];
982 drv->vec[pdof[0]] += +COEFF13 * drv->vec[cdof];
983 drv->vec[pdof[1]] += +COEFF12 * drv->vec[cdof];
984 drv->vec[pdof[2]] += -COEFF14 * drv->vec[cdof];
985
986 cdof = child->dof[node][n0+2];
987 drv->vec[pdof[0]] += -COEFF12 * drv->vec[cdof];
988 drv->vec[pdof[1]] += -COEFF13 * drv->vec[cdof];
989 drv->vec[pdof[2]] += +COEFF14 * drv->vec[cdof];
990
991 cdof = child->dof[node][n0+4];
992 drv->vec[pdof[1]] += 0.125*(-22.0+7.0*SQRT10)/DENOM2 * drv->vec[cdof];
993 drv->vec[pdof[0]] += 0.125*(-3.0*SQRT10+10.0)/DENOM2 * drv->vec[cdof];
994 drv->vec[pdof[2]] += 0.125*(-52.0+16.0*SQRT10)/DENOM2 * drv->vec[cdof];
995
996 cdof = child->dof[node][n0+3];
997 drv->vec[pdof[1]] += 1.0/40.0*(43.0*SQRT10-150.0)/DENOM2 * drv->vec[cdof];
998 drv->vec[pdof[0]] += 1.0/40.0*(-7.0*SQRT10+10.0)/DENOM2 * drv->vec[cdof];
999 drv->vec[pdof[2]] += 1.0/40.0*(-20.0+8.0*SQRT10)/DENOM2 * drv->vec[cdof];
1000
1001 cdof = child->dof[node][n0+5];
1002 drv->vec[pdof[1]] += -1.0/40.0*(13.0*SQRT10-40.0)/DENOM2 * drv->vec[cdof];
1003 drv->vec[pdof[0]] += -1.0/40.0*(-17.0*SQRT10+80.0)/DENOM2 * drv->vec[cdof];
1004 drv->vec[pdof[2]] += -1.0/40.0*(-12.0*SQRT10+40.0)/DENOM2 * drv->vec[cdof];
1005
1006 for (j = 0; j < N_EDGES_2D; j++) {
1007 pdof[j] = el->dof[node][n0+j+N_VERTICES_2D];
1008 drv->vec[pdof[j]] *= 0.5;
1009 }
1010 }
1011 }
1012
1013 #undef SQRT2
1014 #undef SQRT6
1015 #undef SQRT10
1016 #undef SQRT15
1017 #undef DENOM
1018 #undef DENOM2
1019 #undef COEFF11
1020 #undef COEFF12
1021 #undef COEFF13
1022 #undef COEFF14
1023
1024 /* Interpolation: we always use L2-interpolation over the entire
1025 * element, even if WALL >= 0. The quadrature degree is such that we
1026 * would reproduce the polynomials of the given degree.
1027 */
1028 static void
d_ortho_interpol_2_2d(EL_REAL_VEC * el_vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_AT_QP f,void * f_data,const BAS_FCTS * thisptr)1029 d_ortho_interpol_2_2d(EL_REAL_VEC *el_vec,
1030 const EL_INFO *el_info, int wall,
1031 int no, const int *b_no,
1032 LOC_FCT_AT_QP f, void *f_data,
1033 const BAS_FCTS *thisptr)
1034 {
1035 const QUAD_FAST *qfast = ((ORTHO_DATA *)thisptr->ext_data)->qfast;
1036
1037 if (b_no) {
1038 int i, iq;
1039
1040 for (i = 0; i < no; i++) {
1041 el_vec->vec[b_no[i]] = 0.0;
1042 }
1043 for (iq = 0; iq < qfast->n_points; iq++) {
1044 REAL value = qfast->w[iq] * f(el_info, qfast->quad, iq, f_data);
1045 for (i = 0; i < no; i++) {
1046 int ib = b_no[i];
1047 el_vec->vec[ib] += qfast->phi[iq][ib] * value;
1048 }
1049 }
1050 } else {
1051 int iq, ib;
1052
1053 for (ib = 0; ib < N_BAS_LAG_2_2D; ib++) {
1054 el_vec->vec[ib] = 0.0;
1055 }
1056 for (iq = 0; iq < qfast->n_points; iq++) {
1057 REAL value = qfast->w[iq] * f(el_info, qfast->quad, iq, f_data);
1058 for (ib = 0; ib < N_BAS_LAG_2_2D; ib++) {
1059 el_vec->vec[ib] += qfast->phi[iq][ib] * value;
1060 }
1061 }
1062 }
1063 }
1064
1065 static void
d_ortho_interpol_d_2_2d(EL_REAL_D_VEC * el_vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_D_AT_QP f,void * f_data,const BAS_FCTS * thisptr)1066 d_ortho_interpol_d_2_2d(EL_REAL_D_VEC *el_vec,
1067 const EL_INFO *el_info, int wall,
1068 int no, const int *b_no,
1069 LOC_FCT_D_AT_QP f, void *f_data,
1070 const BAS_FCTS *thisptr)
1071 {
1072 const QUAD_FAST *qfast = ((ORTHO_DATA *)thisptr->ext_data)->qfast;
1073
1074 if (b_no) {
1075 int i, iq;
1076
1077 for (i = 0; i < no; i++) {
1078 SET_DOW(0.0, el_vec->vec[b_no[i]]);
1079 }
1080 for (iq = 0; iq < qfast->n_points; iq++) {
1081 REAL_D value;
1082
1083 f(value, el_info, qfast->quad, iq, f_data);
1084 SCAL_DOW(qfast->w[iq], value);
1085 for (i = 0; i < no; i++) {
1086 int ib = b_no[i];
1087 AXPY_DOW( qfast->phi[iq][ib], value, el_vec->vec[ib]);
1088 }
1089 }
1090 } else {
1091 int iq, ib;
1092
1093 for (ib = 0; ib < N_BAS_LAG_2_2D; ib++) {
1094 SET_DOW(0.0, el_vec->vec[ib]);
1095 }
1096 for (iq = 0; iq < qfast->n_points; iq++) {
1097 REAL_D value;
1098
1099 f(value, el_info, qfast->quad, iq, f_data);
1100 SCAL_DOW(qfast->w[iq], value);
1101 for (ib = 0; ib < N_BAS_LAG_2_2D; ib++) {
1102 AXPY_DOW(qfast->phi[iq][ib], value, el_vec->vec[ib]);
1103 }
1104 }
1105 }
1106 }
1107
1108 static void
d_ortho_interpol_dow_2_2d(EL_REAL_VEC_D * el_vec,const EL_INFO * el_info,int wall,int no,const int * b_no,LOC_FCT_D_AT_QP f,void * f_data,const BAS_FCTS * thisptr)1109 d_ortho_interpol_dow_2_2d(EL_REAL_VEC_D *el_vec,
1110 const EL_INFO *el_info, int wall,
1111 int no, const int *b_no,
1112 LOC_FCT_D_AT_QP f, void *f_data,
1113 const BAS_FCTS *thisptr)
1114 {
1115 d_ortho_interpol_d_1_2d((EL_REAL_D_VEC *)el_vec,
1116 el_info, wall, no, b_no, f, f_data, thisptr);
1117 }
1118
1119 /* The trace-space is degenerated: all basis functions have non-zero
1120 * trace on all walls. If permute the stuff cyclic s.t. we would
1121 * eventually have the chance to really define the trace-space as
1122 * degenerated set of basis functions.
1123 */
1124 static const int trace_mapping_d_ortho_2_2d[N_WALLS_2D][N_BAS_LAG_2_2D] = {
1125 { 1, 2, 3, 0, 4, 5 }, { 2, 0, 4, 1, 5, 3 }, { 0, 1, 5, 2, 3, 4 }
1126 };
1127
1128
1129 static const BAS_FCTS disc_ortho2_2d = {
1130 DISC_ORTHO_NAME(2)"_2d", 2, 1, N_BAS_LAG_2_2D, N_BAS_LAG_2_2D, 2,
1131 {0, N_BAS_LAG_2_2D, 0, 0},
1132 -1, /* trace_admin */
1133 INIT_BFCTS_CHAIN(disc_ortho2_2d),
1134 INIT_ELEMENT_INITIALIZER(NULL, FILL_NOTHING), /* init_element + fill-flags */
1135 d_ortho_phi_2d, d_ortho_grd_phi_2d, d_ortho_D2_phi_2d,
1136 NULL, NULL, /* third and fourth derivatives */
1137 NULL, NULL, NULL, false, /* phi_d etc. */
1138 /********************/
1139 NULL /* &disc_lagrange2_1d */, /* trace space */
1140 { { { trace_mapping_d_ortho_2_2d[0],
1141 trace_mapping_d_ortho_2_2d[1],
1142 trace_mapping_d_ortho_2_2d[2], }, }, }, /* trace mapping */
1143 { N_BAS_LAG_2_2D,
1144 N_BAS_LAG_2_2D,
1145 N_BAS_LAG_2_2D, }, /* n_trace_bas_fcts */
1146 d_ortho_get_dof_indices_2_2d,
1147 d_ortho_get_bound_2_2d,
1148 d_ortho_interpol_2_2d,
1149 d_ortho_interpol_d_2_2d,
1150 d_ortho_interpol_dow_2_2d,
1151 d_ortho_get_int_vec_2_2d,
1152 d_ortho_get_real_vec_2_2d,
1153 d_ortho_get_real_d_vec_2_2d,
1154 (GET_REAL_VEC_D_TYPE)d_ortho_get_real_d_vec_2_2d,
1155 d_ortho_get_uchar_vec_2_2d,
1156 d_ortho_get_schar_vec_2_2d,
1157 d_ortho_get_ptr_vec_2_2d,
1158 d_ortho_get_real_dd_vec_2_2d,
1159 d_ortho_real_refine_inter_2_2d,
1160 d_ortho_real_coarse_inter_2_2d,
1161 NULL,
1162 NULL,
1163 NULL,
1164 NULL,
1165 NULL,
1166 NULL,
1167 NULL,
1168 (void *)&d_ortho_2_2d_data
1169 };
1170
1171