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 /* */
7 /* file: parametric_2d.c */
8 /* */
9 /* */
10 /* description: Support for parametric elements in 2D */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Daniel Koester */
14 /* Institut fuer Mathematik */
15 /* Universitaet Augsburg */
16 /* Universitaetsstr. 14 */
17 /* D-86159 Augsburg, Germany */
18 /* */
19 /* Claus-Justus Heine */
20 /* Numerische Mathematik fuer Hoechstleistungsrechner */
21 /* IANS / Universitaet Stuttgart */
22 /* Pfaffenwaldring 57 */
23 /* 70569 Stuttgart, Germany */
24 /* */
25 /* http://www.alberta-fem.de */
26 /* */
27 /* (c) by D. Koester (2005), */
28 /* C.-J. Heine (2006-2012). */
29 /*--------------------------------------------------------------------------*/
30
31 #ifdef MESH_DIM
32 # undef MESH_DIM
33 #endif
34 #define MESH_DIM 2
35
36 #ifdef N_BAS_MAX
37 # undef N_BAS_MAX
38 #endif
39 #define N_BAS_MAX N_BAS_LAGRANGE(LAGRANGE_DEG_MAX, MESH_DIM)
40
41 static const REAL_B mid_lambda_2d = INIT_BARY_2D(0.5, 0.5, 0.0);
42 static const REAL_B bary2_2d[N_BAS_LAG_2D(2)] = {
43 INIT_BARY_2D(1.0, 0.0, 0.0),
44 INIT_BARY_2D(0.0, 1.0, 0.0),
45 INIT_BARY_2D(0.0, 0.0, 1.0),
46 INIT_BARY_2D(0.0, 0.5, 0.5),
47 INIT_BARY_2D(0.5, 0.0, 0.5),
48 INIT_BARY_2D(0.5, 0.5, 0.0)
49 };
50
51 /*--------------------------------------------------------------------------*/
52 /* Functions for affine elements as parametric elements. (suffix 1_2d) */
53 /*--------------------------------------------------------------------------*/
54
param_init_element1_2d(const EL_INFO * el_info,const PARAMETRIC * parametric)55 static bool param_init_element1_2d(const EL_INFO *el_info,
56 const PARAMETRIC *parametric)
57 {
58 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)parametric->data;
59 DOF_REAL_D_VEC *coords = data->coords;
60 int node_v, n0_v, i;
61 EL *el = el_info->el;
62 EL_INFO *mod_el_info = (EL_INFO *)el_info; /* modifyable pointer reference */
63
64 data->el = el_info->el;
65
66 node_v = el_info->mesh->node[VERTEX];
67 n0_v = coords->fe_space->admin->n0_dof[VERTEX];
68
69 if (!parametric->use_reference_mesh) {
70 data->local_coords = mod_el_info->coord;
71 mod_el_info->fill_flag |= FILL_COORDS;
72 } else {
73 data->local_coords = data->param_local_coords;
74 }
75 for (i = 0; i < N_VERTICES_2D; i++) {
76 COPY_DOW(coords->vec[el->dof[node_v+i][n0_v]], data->local_coords[i]);
77 }
78
79 return false;
80 }
81
det1_2d(const EL_INFO * el_info,const QUAD * quad,int N,const REAL_B lambda[],REAL dets[])82 static void det1_2d(const EL_INFO *el_info, const QUAD *quad, int N,
83 const REAL_B lambda[], REAL dets[])
84 {
85 REAL det;
86 int n;
87
88 det = el_det_2d(el_info);
89
90 if (quad) {
91 N = quad->n_points;
92 }
93
94 for (n = 0; n < N; n++) {
95 dets[n] = det;
96 }
97
98 return;
99 }
100
grd_lambda1_2d(const EL_INFO * el_info,const QUAD * quad,int N,const REAL_B lambda[],REAL_BD grd_lam[],REAL_BDD D2_lam[],REAL dets[])101 static void grd_lambda1_2d(const EL_INFO *el_info, const QUAD *quad,
102 int N, const REAL_B lambda[],
103 REAL_BD grd_lam[], REAL_BDD D2_lam[], REAL dets[])
104 {
105 int i, n;
106
107 dets[0] = el_grd_lambda_2d(el_info, grd_lam[0]);
108
109 if (quad) {
110 N = quad->n_points;
111 }
112
113 for (n = 1; n < N; n++) {
114 for (i = 0; i < N_LAMBDA_2D; i++) {
115 COPY_DOW(grd_lam[0][i], grd_lam[n][i]);
116 }
117 for (; i < N_LAMBDA_MAX; i++) {
118 SET_DOW(0.0, grd_lam[n][i]);
119 }
120 if (dets) {
121 dets[n] = dets[0];
122 }
123 }
124
125 if (D2_lam) {
126 for (n = 0; n < N; n++) {
127 for (i = 0; i < N_LAMBDA_MAX;i++) {
128 MSET_DOW(0.0, D2_lam[n][i]);
129 }
130 }
131 }
132
133 return;
134 }
135
136 static void
grd_world1_2d(const EL_INFO * el_info,const QUAD * quad,int N,const REAL_B lambda[],REAL_BD grd_Xtr[],REAL_BDB D2_Xtr[],REAL_BDBB D3_Xtr[])137 grd_world1_2d(const EL_INFO *el_info, const QUAD *quad,
138 int N, const REAL_B lambda[],
139 REAL_BD grd_Xtr[], REAL_BDB D2_Xtr[], REAL_BDBB D3_Xtr[])
140 {
141 int i;
142
143 if (quad) {
144 N = quad->n_points;
145 }
146
147 for (i = 0; i < N_LAMBDA_2D; i++) {
148 COPY_DOW(el_info->coord[i], grd_Xtr[0][i]);
149 }
150 for (; i < N_LAMBDA_MAX; i++) {
151 SET_DOW(0.0, grd_Xtr[0][i]);
152 }
153
154 memcpy(grd_Xtr+1, grd_Xtr[0], (N-1)*sizeof(REAL_BDB));
155
156 if (D2_Xtr) {
157 memset(D2_Xtr, 0, N*sizeof(REAL_BBD));
158 }
159
160 if (D3_Xtr) {
161 memset(D3_Xtr, 0, N*sizeof(REAL_BDBB));
162 }
163 }
164
165 /* Compute the co-normal for WALL at the barycentric coordinates
166 * specified by LAMBDA. LAMBDA are 2d barycentric coordinates,
167 * i.e. lambda[wall] == 0.0.
168 */
169 static void
wall_normal1_2d(const EL_INFO * el_info,int wall,const QUAD * quad,int n,const REAL_B lambda[],REAL_D normals[],REAL_DB grd_normals[],REAL_DBB D2_normals[],REAL dets[])170 wall_normal1_2d(const EL_INFO *el_info, int wall,
171 const QUAD *quad, int n, const REAL_B lambda[],
172 REAL_D normals[], REAL_DB grd_normals[], REAL_DBB D2_normals[],
173 REAL dets[])
174 {
175 int i;
176
177 if (quad) {
178 n = quad->n_points;
179 }
180
181 if (grd_normals) {
182 memset(grd_normals, 0, n*sizeof(REAL_DB));
183 }
184
185 if (D2_normals) {
186 memset(D2_normals, 0, n*sizeof(REAL_DBB));
187 }
188
189 if (normals) {
190 REAL detsbuffer[n];
191
192 if (!dets) {
193 dets = detsbuffer;
194 }
195 dets[0] = get_wall_normal_2d(el_info, wall, normals[0]);
196
197 for (i = 1; i < n; i++) {
198 dets[i] = dets[0];
199 COPY_DOW(normals[0], normals[i]);
200 }
201 } else {
202 dets[0] = get_wall_normal_2d(el_info, wall, NULL);
203 for (i = 1; i < n; i++) {
204 dets[i] = dets[0];
205 }
206 }
207 }
208
209 /****************************************************************************/
210 /* fill_coords1_2d(data): initialize the DOF_REAL_D_VEC coords containing */
211 /* the position data of the parametric elements (coordinates of vertices in */
212 /* this case). */
213 /****************************************************************************/
214
fill_coords1_2d(LAGRANGE_PARAM_DATA * data)215 static void fill_coords1_2d(LAGRANGE_PARAM_DATA *data)
216 {
217 DOF_REAL_D_VEC *coords = data->coords;
218 NODE_PROJECTION *n_proj = data->n_proj;
219 bool selective = n_proj != NULL;
220 const NODE_PROJECTION *act_proj;
221 FLAGS fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_PROJECTION;
222 MESH *mesh;
223 const DOF_ADMIN *admin;
224 const BAS_FCTS *bas_fcts;
225 int i;
226 DOF dof[N_VERTICES_2D];
227
228 mesh = coords->fe_space->mesh;
229 admin = coords->fe_space->admin;
230 bas_fcts = coords->fe_space->bas_fcts;
231
232 TRAVERSE_FIRST(mesh, -1, fill_flag|FILL_BOUND) {
233
234 GET_DOF_INDICES(bas_fcts, el_info->el, admin, dof);
235
236 for (i = 0; i < N_VERTICES_2D; i++) {
237 REAL *vec = coords->vec[dof[i]];
238
239 COPY_DOW(el_info->coord[i], vec);
240
241 /* Look for a projection function that applies to vertex[i]. */
242 /* Apply this projection if found. */
243
244 if (!selective || n_proj->func) {
245 act_proj = wall_proj(el_info, (i+1) % N_WALLS_2D);
246 if (act_proj == NULL) {
247 act_proj = wall_proj(el_info, (i+2) % N_WALLS_2D);
248 }
249 if (!act_proj) {
250 act_proj = wall_proj(el_info, -1);
251 }
252 if (act_proj && act_proj->func && (!selective || act_proj == n_proj)){
253 act_proj->func(vec, el_info, bary2_2d[i]);
254 }
255 }
256 }
257 } TRAVERSE_NEXT();
258 }
259
260
261 /****************************************************************************/
262 /* refine_interpol1_2d(drdv,list,n): update coords vector during refinement.*/
263 /****************************************************************************/
264
refine_interpol1_2d(DOF_REAL_D_VEC * drdv,RC_LIST_EL * list,int n)265 static void refine_interpol1_2d(DOF_REAL_D_VEC *drdv, RC_LIST_EL *list, int n)
266 {
267 /*FUNCNAME("refine_interpol1_2d");*/
268 MESH *mesh = drdv->fe_space->mesh;
269 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data;
270 REAL_D *vec = drdv->vec;
271 NODE_PROJECTION *n_proj = data->n_proj;
272 bool selective = n_proj != NULL;
273 const NODE_PROJECTION *act_proj;
274 EL *el;
275 DOF dof_new, dof0, dof1;
276 int n0, j;
277
278 n0 = drdv->fe_space->admin->n0_dof[VERTEX];
279 el = list->el_info.el;
280
281 dof0 = el->dof[0][n0]; /* 1st endpoint of refinement edge */
282 dof1 = el->dof[1][n0]; /* 2nd endpoint of refinement edge */
283 dof_new = el->child[0]->dof[2][n0]; /* vertex[2] is newest vertex */
284
285 for (j = 0; j < DIM_OF_WORLD; j++)
286 vec[dof_new][j] = 0.5*(vec[dof0][j] + vec[dof1][j]);
287
288 act_proj = list->el_info.active_projection;
289
290 if (act_proj && act_proj->func && (!selective || n_proj == act_proj)) {
291 act_proj->func(vec[dof_new], &list->el_info, mid_lambda_2d);
292 _AI_refine_update_bbox(mesh, vec[dof_new]);
293 }
294 }
295
vertex_coordsY_2d(EL_INFO * el_info)296 static void vertex_coordsY_2d(EL_INFO *el_info)
297 {
298 PARAMETRIC *parametric = el_info->mesh->parametric;
299 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)parametric->data;
300 DOF_REAL_D_VEC *coords = data->coords;
301 EL *el = el_info->el;
302 int node_v, n0_v, i;
303
304 node_v = el_info->mesh->node[VERTEX];
305 n0_v = coords->fe_space->admin->n0_dof[VERTEX];
306
307 el_info->fill_flag |= FILL_COORDS;
308 for (i = 0; i < N_VERTICES_2D; i++) {
309 COPY_DOW(coords->vec[el->dof[node_v+i][n0_v]], el_info->coord[i]);
310 }
311 }
312
313 /*--------------------------------------------------------------------------*/
314 /* Common functions for higher order elements. (suffix Y_2d) */
315 /*--------------------------------------------------------------------------*/
316
param_init_elementY_2d(const EL_INFO * el_info,const PARAMETRIC * parametric)317 static bool param_init_elementY_2d(const EL_INFO *el_info,
318 const PARAMETRIC *parametric)
319 {
320 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)parametric->data;
321 DOF_PTR_VEC *edge_pr = data->edge_projections;
322 DOF_REAL_D_VEC *coords = data->coords;
323 EL *el = el_info->el;
324 int i;
325
326 if (data->el != el) {
327
328 data->el = el;
329
330 if (data->strategy == PARAM_ALL) {
331 /* full-featured parametric element */
332 const BAS_FCTS *bas_fcts = data->coords->fe_space->bas_fcts;
333
334 bas_fcts->get_real_d_vec(data->local_coords, el, coords);
335
336 return true;
337 } else {
338 int node_e = el_info->mesh->node[EDGE];
339 int n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
340
341 /* An element is affine iff none of its edges has suffered a
342 * projection.
343 */
344 data->i_am_affine = true;
345 for (i = 0; i < N_EDGES_2D; i++) {
346 if (edge_pr->vec[el->dof[node_e+i][n0_edge_pr]] != NULL) {
347 /* full-featured parametric element */
348 const BAS_FCTS *bas_fcts = data->coords->fe_space->bas_fcts;
349
350 data->i_am_affine = false;
351 data->local_coords = data->param_local_coords;
352 bas_fcts->get_real_d_vec(data->local_coords, el, coords);
353
354 return true;
355 }
356 }
357 if (parametric->use_reference_mesh) {
358 const BAS_FCTS *bas_fcts = data->coords->fe_space->bas_fcts;
359 data->local_coords = data->param_local_coords;
360 bas_fcts->get_real_d_vec(data->local_coords, el, coords);
361 }
362 }
363 }
364
365 if (!parametric->use_reference_mesh) {
366 EL_INFO *mod_el_info = (EL_INFO *)el_info;
367 if (data->i_am_affine) {
368 /* Need only the vertices. */
369 int node_v, n0_v;
370
371 node_v = el_info->mesh->node[VERTEX];
372 n0_v = coords->fe_space->admin->n0_dof[VERTEX];
373
374 data->local_coords = (REAL_D *)el_info->coord;
375 mod_el_info->fill_flag |= FILL_COORDS;
376 for (i = 0; i < N_VERTICES_2D; i++) {
377 COPY_DOW(coords->vec[el->dof[node_v+i][n0_v]], data->local_coords[i]);
378 }
379 } else {
380 mod_el_info->fill_flag &= ~FILL_COORDS;
381 }
382 }
383
384 return !data->i_am_affine;
385 }
386
387 /* Much of the parametric stuff is actually dimension independent,
388 * therefore a bunch of functions shared between (at least)
389 * parametric_2d.c and parametric_3d.c is implemented in the following
390 * header-file:
391 */
392 #include "parametric_intern.h"
393
394 /* Compute the co-normal for WALL at a given quadrature point. Helper
395 * routine for wall_normalY_2d().
396 *
397 * On request we compute also the barycentric gradient of the normal,
398 * this is needed by vector-valued basis function which point into
399 * normal direction.
400 */
401 static inline
wall_normal_iq_2d(REAL_D * const coords,REAL DD[][MESH_DIM],REAL dDD[][MESH_DIM][MESH_DIM],REAL ddDD[][MESH_DIM][MESH_DIM][MESH_DIM],int n_bas,int wall,REAL_D normal,REAL_DB grd_normal,REAL_DBB D2_normal)402 REAL wall_normal_iq_2d(REAL_D *const coords,
403 REAL DD[][MESH_DIM],
404 REAL dDD[][MESH_DIM][MESH_DIM],
405 REAL ddDD[][MESH_DIM][MESH_DIM][MESH_DIM],
406 int n_bas, int wall,
407 REAL_D normal, REAL_DB grd_normal, REAL_DBB D2_normal)
408 {
409 REAL w[MESH_DIM], g[MESH_DIM][MESH_DIM], ncoeff[MESH_DIM], detg, det;
410 REAL_D e[MESH_DIM];
411 int alpha;
412
413 /* Make the stuff a little bit more similar to the 3d case */
414 detg = Dt_and_DtD_2d(coords, DD, n_bas, e, g);
415
416 for (alpha = 0; alpha < MESH_DIM; alpha++) {
417 w[alpha] = g[alpha][0] - g[alpha][1];
418 }
419
420 /* The coefficient vector of the normal must be orthogonal to the
421 * coefficient vector w of (e0 - e1) w.r.t. to g, so the following
422 * should be the outer normal (at least the orientation of the
423 * normal is determined, plese FIXME if it should be the inner
424 * normal ... at least it has a defined sign).
425 */
426 ncoeff[0] = -w[1];
427 ncoeff[1] = +w[0];
428
429 AXPBY_DOW(ncoeff[0], e[0], ncoeff[1], e[1], normal);
430
431 /* The length of "normal" is sqrt(det g)*|wall|, as in the 3d
432 * case. We correct this later. By construction "normal" is the
433 * outer normal to the simplex.
434 *
435 * We can now proceed as in the 3d case and compute the Jacobian of
436 * the normal field.
437 */
438 if (grd_normal || D2_normal) {
439 REAL_D De[MESH_DIM][MESH_DIM];
440 REAL Dg[MESH_DIM][MESH_DIM][MESH_DIM];
441 REAL Dncoeff[MESH_DIM][MESH_DIM];
442 REAL nuDnu[MESH_DIM];
443 REAL_D grd_tmp[MESH_DIM], grd_loc[MESH_DIM];
444 REAL inv_nrm2, inv_norm;
445 int i, beta, gamma;
446
447 dDt_and_dDtD_2d(coords, dDD, e, n_bas, De, Dg);
448
449 /* normal = (g_10 - g11)*e0 - (g00 -g01) * e1 */
450
451 /* First part: differentiate the basis vector of the tangent space */
452 for (alpha = 0; alpha < MESH_DIM; alpha++) {
453 AXPBY_DOW(ncoeff[0], De[alpha][0], ncoeff[1], De[alpha][1],
454 grd_tmp[alpha]);
455 }
456
457 /* second part: differentiate the coefficients. */
458 for (alpha = 0; alpha < MESH_DIM; alpha++) {
459 Dncoeff[alpha][0] = Dg[alpha][1][0]-Dg[alpha][1][1];
460 Dncoeff[alpha][1] = Dg[alpha][0][1]-Dg[alpha][0][0];
461 AXPY_DOW(Dncoeff[alpha][0], e[0], grd_tmp[alpha]);
462 AXPY_DOW(Dncoeff[alpha][1], e[1], grd_tmp[alpha]);
463 }
464
465 /* Now it remains to take the normalization into account, and to
466 * re-order the barycentric co-ordinates s.t. they match ordering
467 * of the bulk element.
468 */
469 inv_nrm2 =
470 grd_normalize_2d(grd_loc, nuDnu, (const REAL_D *)grd_tmp, normal);
471 inv_norm = sqrt(inv_nrm2);
472
473 if (grd_normal) {
474 /* Sort the stuff s.t. it uses the proper order of the
475 * barycentric co-ordinates
476 */
477 for (i = 0; i < DIM_OF_WORLD; i++) {
478 /* constant extension along \lambda_wall co-ordinate lines. */
479 grd_normal[i][wall] = 0.0;
480 for (alpha = 0; alpha < MESH_DIM; alpha++) {
481 grd_normal[i][(wall+alpha+1) % N_LAMBDA(MESH_DIM)] =
482 inv_norm * grd_loc[alpha][i];
483 }
484 }
485 }
486
487 if (D2_normal) {
488 REAL_D D2e[MESH_DIM][MESH_DIM][MESH_DIM];
489 REAL_DDDD_2D D2g;
490 REAL_D D2_tmp[MESH_DIM][MESH_DIM], D2_loc[MESH_DIM][MESH_DIM];
491
492 ddDt_and_ddDtD_2d(coords, ddDD, e, De, Dg, n_bas, D2e, D2g);
493
494 /* First part: take the second derivatives of the tangent vectors */
495 for (alpha = 0; alpha < MESH_DIM; alpha++) {
496 for (gamma = 0; gamma < MESH_DIM; gamma++) {
497 AXPY_DOW(ncoeff[gamma], D2e[alpha][alpha][gamma],
498 D2_tmp[alpha][alpha]);
499 }
500 for (beta = alpha+1; beta < MESH_DIM; beta++) {
501 for (gamma = 0; gamma < MESH_DIM; gamma++) {
502 AXPY_DOW(ncoeff[gamma], D2e[alpha][beta][gamma],
503 D2_tmp[alpha][beta]);
504 }
505 }
506 }
507
508 /* Second part: first derivatives of coefficients, with first
509 * derivatives of the tangent vectors.
510 */
511 for (alpha = 0; alpha < MESH_DIM; alpha++) {
512 for (gamma = 0; gamma < MESH_DIM; gamma++) {
513 AXPY_DOW(2.0*Dncoeff[alpha][gamma], De[alpha][gamma],
514 D2_tmp[alpha][beta]);
515 }
516 for (beta = alpha+1; beta < MESH_DIM; beta++) {
517 for (gamma = 0; gamma < MESH_DIM; gamma++) {
518 AXPY_DOW(Dncoeff[alpha][gamma], De[beta][gamma],
519 D2_tmp[alpha][beta]);
520 AXPY_DOW(Dncoeff[beta][gamma], De[alpha][gamma],
521 D2_tmp[alpha][beta]);
522 }
523 }
524 }
525
526 /* Third part: second derivatives of coefficients, with values
527 * of the tangent vectors.
528 */
529 for (alpha = 0; alpha < MESH_DIM; alpha++) {
530 for (beta = alpha+1; beta < MESH_DIM; beta++) {
531 AXPY_DOW(D2g[alpha][beta][1][0] - D2g[alpha][beta][1][1],
532 e[0], D2_tmp[alpha][beta]);
533 AXPY_DOW(D2g[alpha][beta][0][1] - D2g[alpha][beta][0][0],
534 e[1], D2_tmp[alpha][beta]);
535 }
536 }
537
538 /* Now we have computed the second derivatives of the non-unit
539 * normals. We have to take the normalization into account.
540 */
541 D2_normalize_2d(D2_loc, normal,
542 (const REAL_D *)grd_tmp,
543 (const REAL_D (*)[MESH_DIM])D2_tmp, inv_nrm2, nuDnu);
544
545 /* Sort the stuff s.t. it uses the proper order of the
546 * barycentric co-ordinates
547 */
548 for (i = 0; i < DIM_OF_WORLD; i++) {
549 D2_normal[i][wall][wall] = 0.0;
550 for (alpha = 0; alpha < MESH_DIM; alpha++) {
551 int a = (wall + alpha + 1) % N_LAMBDA(MESH_DIM);
552 D2_normal[i][a][a] = inv_norm * D2_loc[alpha][alpha][i];
553 D2_normal[i][wall][a] = D2_normal[i][a][wall] = 0.0;
554 for (beta = alpha + 1; beta < MESH_DIM; beta++) {
555 int b = (wall + beta + 1) % N_LAMBDA(MESH_DIM);
556
557 D2_normal[i][a][b] =
558 D2_normal[i][b][a] = inv_norm * D2_loc[alpha][beta][i];
559 }
560 }
561 }
562 }
563 }
564
565 /* Now scale the normal s.t. it reflects the magnitude of the line element */
566 SCAL_DOW(1.0/sqrt(detg), normal);
567
568 #if 1
569 /* In principle the orientation should be fixed ... */
570 if (SCP_DOW(e[1], normal) <= 0.0) {
571 if (ALBERTA_DEBUG) {
572 WARNING("Wrong orientation?\n");
573 }
574 SCAL_DOW(-1.0, normal);
575 }
576 #endif
577
578 det = NORM_DOW(normal);
579
580 /* cH: shouldn't the following test the relative size of det w.r.t.,
581 * e.g. the edge length?
582 */
583 TEST_EXIT(det > 1.e-30, "face det = 0 on face %d.\n", wall);
584
585 return det;
586 }
587
588 /*--------------------------------------------------------------------------*/
589 /* Functions for quadratic elements. (suffix 2_2d) */
590 /*--------------------------------------------------------------------------*/
591
refine_interpol2_2d(DOF_REAL_D_VEC * drdv,RC_LIST_EL * list,int n)592 static void refine_interpol2_2d(DOF_REAL_D_VEC *drdv, RC_LIST_EL *list, int n)
593 {
594 /*FUNCNAME("refine_interpol2_2d");*/
595 MESH *mesh = drdv->fe_space->mesh;
596 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data;
597 const DOF_ADMIN *admin = drdv->fe_space->admin;
598 const BAS_FCTS *bas_fcts = drdv->fe_space->bas_fcts;
599 DOF_PTR_VEC *edge_pr = data->edge_projections;
600 REAL_D *vec = drdv->vec;
601 const EL_INFO *el_info;
602 EL *el;
603 int node_v, node_e, n0_v, n0_e, i, j, i_child;
604 DOF cdof, cdof_e[N_EDGES_2D], cdof_v[N_VERTICES_2D], pdof;
605 REAL *x[3];
606 NODE_PROJECTION *n_proj;
607 bool selective;
608 const NODE_PROJECTION *act_proj;
609 int n0_edge_pr = -1;
610 DOF cdof_edge_pr[3] = { -1, };
611 static const REAL_B lambda[3] = {INIT_BARY_2D(0.25, 0.25, 0.5),
612 INIT_BARY_2D(0.75, 0.25, 0.0),
613 INIT_BARY_2D(0.25, 0.75, 0.0) };
614
615 el_info = &list->el_info;
616 el = el_info->el;
617 n_proj = data->n_proj;
618 selective = n_proj != NULL;
619
620 node_v = mesh->node[VERTEX];
621 node_e = mesh->node[EDGE];
622 n0_v = admin->n0_dof[VERTEX];
623 n0_e = admin->n0_dof[EDGE];
624
625 /****************************************************************************/
626 /* Step 1: Now we need to determine three points: The midpoint on the edge */
627 /* between the two children and the two midpoints of the children along the */
628 /* current refinement edge. These three points correspond to the three new */
629 /* DOFS that were created. */
630 /* If strategy is PARAM_ALL or PARAM_CURVED_CHILDS we use the (nonlinear) */
631 /* barycentric coords. For strategy PARAM_STRAIGHT_CHILD we use geometric */
632 /* midpoints. */
633 /****************************************************************************/
634 x[0] = vec[el->child[0]->dof[node_e+1][n0_e]];
635 x[1] = vec[el->child[0]->dof[node_e+0][n0_e]];
636 x[2] = vec[el->child[1]->dof[node_e+1][n0_e]];
637
638 if (data->strategy != PARAM_STRAIGHT_CHILDS) {
639 /* Just call the default interpolation routine */
640 bas_fcts->real_d_refine_inter(drdv, list, n);
641 } else {
642 cdof_e[2] = el->dof[node_e+2][n0_e];
643 for (i = 0; i < N_EDGES_2D; i++)
644 cdof_v[i] = el->dof[node_v+i][n0_v];
645
646 for (i = 0; i < DIM_OF_WORLD; i++) {
647 x[0][i] = 0.5 * (vec[cdof_v[2]][i] + vec[cdof_e[2]][i]);
648 x[1][i] = 0.5 * (vec[cdof_v[0]][i] + vec[cdof_e[2]][i]);
649 x[2][i] = 0.5 * (vec[cdof_v[1]][i] + vec[cdof_e[2]][i]);
650 }
651 }
652
653 /****************************************************************************/
654 /* Step 2: We check if any projections need to be done on these three */
655 /* points. While doing this, we keep track of projected coordinates. */
656 /****************************************************************************/
657 if (edge_pr) {
658 n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
659
660 cdof_edge_pr[0] = el->child[0]->dof[node_e+1][n0_edge_pr];
661 cdof_edge_pr[1] = el->child[0]->dof[node_e+0][n0_edge_pr];
662 cdof_edge_pr[2] = el->child[1]->dof[node_e+1][n0_edge_pr];
663
664 edge_pr->vec[cdof_edge_pr[0]] =
665 edge_pr->vec[cdof_edge_pr[1]] =
666 edge_pr->vec[cdof_edge_pr[2]] = NULL;
667
668 act_proj = wall_proj(el_info, -1);
669 if (act_proj && (!selective || act_proj == n_proj)) {
670 edge_pr->vec[cdof_edge_pr[0]] = (void *)act_proj;
671 }
672 act_proj = el_info->active_projection;
673 if (act_proj && (!selective || act_proj == n_proj)) {
674 edge_pr->vec[cdof_edge_pr[1]] = (void *)act_proj;
675 edge_pr->vec[cdof_edge_pr[2]] = (void *)act_proj;
676 }
677 }
678
679 act_proj = wall_proj(el_info, -1);
680 if (act_proj && act_proj->func && (!selective || act_proj == n_proj)) {
681 act_proj->func(x[0], el_info, lambda[0]);
682 _AI_refine_update_bbox(mesh, x[0]);
683 }
684 act_proj = el_info->active_projection;
685 if (act_proj && act_proj->func && (!selective || act_proj == n_proj)) {
686 act_proj->func(x[1], el_info, lambda[1]);
687 _AI_refine_update_bbox(mesh, x[1]);
688 act_proj->func(x[2], el_info, lambda[2]);
689 _AI_refine_update_bbox(mesh, x[2]);
690 }
691
692 /****************************************************************************/
693 /* Step 3: We hand down the data corresponding to the parent edge vertex */
694 /* in the refinement edge (its DOF could be deleted!). */
695 /****************************************************************************/
696 pdof = el->dof[node_e+2][n0_e];
697 cdof = el->child[0]->dof[node_v+2][n0_v];
698 COPY_DOW(vec[pdof], vec[cdof]);
699
700 /****************************************************************************/
701 /* Step 4: Correct the children, if not all elements are to be parametric. */
702 /* If the parent element was already affine, then the children will already */
703 /* have the correct coordinates. */
704 /****************************************************************************/
705 if (edge_pr) {
706 int parent_affine = true;
707
708 for (j = 0; j < N_EDGES_2D; j++) {
709 DOF edge_dof = el->dof[node_e+j][n0_edge_pr];
710
711 if (edge_pr->vec[edge_dof] != NULL) {
712 parent_affine = false;
713 break;
714 }
715 }
716
717 if (!parent_affine) for (i_child = 0; i_child < 2; i_child++) {
718 EL *child = el->child[i_child];
719 int child_affine = true;
720
721 for (j = 0; j < N_EDGES_2D; j++) {
722 DOF edge_dof = child->dof[node_e+j][n0_edge_pr];
723
724 if (edge_pr->vec[edge_dof] != NULL) {
725 child_affine = false;
726 break;
727 }
728 }
729 if (child_affine) {
730 for (j = 0; j < N_EDGES_2D; j++) {
731 cdof_e[j] = child->dof[node_e+j][n0_e];
732 cdof_v[j] = child->dof[node_v+j][n0_v];
733 }
734
735 for (j = 0; j < DIM_OF_WORLD; j++) {
736 vec[cdof_e[0]][j] = 0.5*(vec[cdof_v[1]][j] + vec[cdof_v[2]][j]);
737 vec[cdof_e[1]][j] = 0.5*(vec[cdof_v[0]][j] + vec[cdof_v[2]][j]);
738 vec[cdof_e[2]][j] = 0.5*(vec[cdof_v[0]][j] + vec[cdof_v[1]][j]);
739 }
740 }
741 }
742 }
743
744 /*--------------------------------------------------------------------------*/
745 /*--- Take care of the neighbour element. ---*/
746 /*--------------------------------------------------------------------------*/
747 if (n > 1) {
748 el_info = &list[1].el_info;
749 el = el_info->el;
750
751 /****************************************************************************/
752 /* Step 1: Determine the location of the midpoint on the edge between the */
753 /* two children. */
754 /****************************************************************************/
755 x[0] = vec[el->child[0]->dof[node_e+1][n0_e]];
756
757 if (data->strategy == PARAM_STRAIGHT_CHILDS) {
758 cdof_e[2] = el->dof[node_e+2][n0_e];
759 cdof_v[2] = el->dof[node_v+2][n0_v];
760
761 AXPBY_DOW(0.5, vec[cdof_v[2]], 0.5, vec[cdof_e[2]], x[0]);
762 }
763
764 /****************************************************************************/
765 /* Step 2: We check if any projections need to be done on this point. */
766 /* While doing this, we keep track of projected coordinates. */
767 /****************************************************************************/
768
769 if (edge_pr) {
770 cdof_edge_pr[0] = el->child[0]->dof[node_e+1][n0_edge_pr];
771 edge_pr->vec[cdof_edge_pr[0]] = NULL;
772
773 act_proj = wall_proj(el_info, -1);
774 if (act_proj && (!selective || act_proj == n_proj)) {
775 cdof_edge_pr[0] = el->child[0]->dof[node_e+1][n0_edge_pr];
776 edge_pr->vec[cdof_edge_pr[0]] = (void *)act_proj;
777 }
778 }
779
780 act_proj = wall_proj(el_info, -1);
781 if (act_proj && act_proj->func && (!selective || act_proj == n_proj)) {
782 act_proj->func(x[0], el_info, lambda[0]);
783 _AI_refine_update_bbox(mesh, x[0]);
784 }
785
786 /****************************************************************************/
787 /* Step 3: Correct the children. */
788 /***************************************************************************/
789 if (edge_pr) {
790 int parent_affine = true;
791
792 for (j = 0; j < N_EDGES_2D; j++) {
793 DOF edge_dof = el->dof[node_e+j][n0_edge_pr];
794
795 if (edge_pr->vec[edge_dof] != NULL) {
796 parent_affine = false;
797 break;
798 }
799 }
800 if (!parent_affine) for (i_child = 0; i_child < 2; i_child++) {
801 EL *child = el->child[i_child];
802 int child_affine = true;
803
804 for (j = 0; j < N_EDGES_2D; j++) {
805 DOF edge_dof = child->dof[node_e+j][n0_edge_pr];
806
807 if (edge_pr->vec[edge_dof] != NULL) {
808 child_affine = false;
809 break;
810 }
811 }
812 if (child_affine) {
813 for (j = 0; j < N_EDGES_2D; j++) {
814 cdof_e[j] = child->dof[node_e+j][n0_e];
815 cdof_v[j] = child->dof[node_v+j][n0_v];
816 }
817
818 for (j = 0; j < DIM_OF_WORLD; j++) {
819 vec[cdof_e[0]][j] = 0.5*(vec[cdof_v[1]][j] + vec[cdof_v[2]][j]);
820 vec[cdof_e[1]][j] = 0.5*(vec[cdof_v[0]][j] + vec[cdof_v[2]][j]);
821 vec[cdof_e[2]][j] = 0.5*(vec[cdof_v[0]][j] + vec[cdof_v[1]][j]);
822 }
823 }
824 }
825 }
826 }
827
828 return;
829 }
830
coarse_interpol2_2d(DOF_REAL_D_VEC * drdv,RC_LIST_EL * list,int n)831 static void coarse_interpol2_2d(DOF_REAL_D_VEC *drdv, RC_LIST_EL *list, int n)
832 {
833 FUNCNAME("coarse_interpol2_2d");
834 EL *el;
835 REAL_D *vec = NULL;
836 int node_v, node_e, n0_v, n0_e, j;
837 DOF cdof, pdof;
838 const MESH *mesh = drdv->fe_space->mesh;
839 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data;
840 DOF_PTR_VEC *edge_pr = data->edge_projections;
841
842 GET_DOF_VEC(vec, drdv);
843 el = list->el_info.el;
844
845 node_v = mesh->node[VERTEX];
846 node_e = mesh->node[EDGE];
847 n0_v = drdv->fe_space->admin->n0_dof[VERTEX];
848 n0_e = drdv->fe_space->admin->n0_dof[EDGE];
849
850 /****************************************************************************/
851 /* copy values at refinement vertex to the parent refinement edge. */
852 /****************************************************************************/
853
854 cdof = el->child[0]->dof[node_v+2][n0_v]; /* newest vertex is dim */
855 pdof = el->dof[node_e+2][n0_e];
856
857 for (j = 0; j < DIM_OF_WORLD; j++)
858 vec[pdof][j] = vec[cdof][j];
859
860 if (edge_pr) {
861 int n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
862 DOF pdof = el->dof[node_e+2][n0_edge_pr];
863 DOF cdof = el->child[0]->dof[node_e+0][n0_edge_pr];
864
865 edge_pr->vec[pdof] = edge_pr->vec[cdof];
866 }
867
868 return;
869 }
870
fill_coords2_2d(LAGRANGE_PARAM_DATA * data)871 static void fill_coords2_2d(LAGRANGE_PARAM_DATA *data)
872 {
873 DOF_REAL_D_VEC *coords = data->coords;
874 MESH *mesh = coords->fe_space->mesh;
875 const DOF_ADMIN *admin = coords->fe_space->admin;
876 const BAS_FCTS *bas_fcts = coords->fe_space->bas_fcts;
877 DOF_PTR_VEC *edge_pr = data->edge_projections;
878 NODE_PROJECTION *n_proj = data->n_proj;
879 bool selective = n_proj != NULL;
880 const NODE_PROJECTION *act_proj;
881 FLAGS fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_PROJECTION;
882 DOF dof[bas_fcts->n_bas_fcts];
883 int i, node_e = -1, n0_edge_pr = -1;
884 DOF dof_edge_pr = -1;
885
886 if (edge_pr) {
887 node_e = mesh->node[EDGE];
888 n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
889 }
890
891 TRAVERSE_FIRST(mesh, -1, fill_flag|FILL_BOUND) {
892 REAL *vec, *vec0, *vec1;
893
894 GET_DOF_INDICES(bas_fcts, el_info->el, admin, dof);
895
896 for (i = 0; i < N_VERTICES_2D; i++) {
897 vec = coords->vec[dof[i]];
898
899 COPY_DOW(el_info->coord[i], vec);
900
901 /* Look for a projection function that applies to vertex[i].
902 * Apply this projection if found.
903 */
904
905 if (!selective || n_proj->func) {
906 act_proj = wall_proj(el_info, (i+1) % N_WALLS_2D);
907
908 if (!act_proj) {
909 act_proj = wall_proj(el_info, (i+2) % N_WALLS_2D);
910 }
911 if (!act_proj) {
912 act_proj = wall_proj(el_info, -1);
913 }
914 if (act_proj && act_proj->func && (!selective || act_proj == n_proj)) {
915 act_proj->func(vec, el_info, bary2_2d[i]);
916 }
917 }
918 }
919
920 for (i = 0; i < N_EDGES_2D; i++) {
921 const int *voe = vertex_of_edge_2d[i];
922
923 vec = coords->vec[dof[N_VERTICES_2D+i]];
924 vec0 = coords->vec[dof[voe[0]]];
925 vec1 = coords->vec[dof[voe[1]]];
926
927 AXPBY_DOW(0.5, vec0, 0.5, vec1, vec);
928
929 /* Look for a projection function that applies to edge[i]. */
930 /* Apply this projection if found. */
931
932 if (edge_pr) {
933 dof_edge_pr = el_info->el->dof[node_e+i][n0_edge_pr];
934 edge_pr->vec[dof_edge_pr] = NULL;
935
936 act_proj = wall_proj(el_info, i);
937 if (!act_proj) {
938 act_proj = wall_proj(el_info, -1);
939 }
940 if (act_proj && (!selective || act_proj == n_proj)) {
941 edge_pr->vec[dof_edge_pr] = (void *)act_proj;
942 }
943 }
944
945 if (!selective || n_proj->func) {
946 act_proj = wall_proj(el_info, i);
947
948 if (!act_proj) {
949 act_proj = wall_proj(el_info, -1);
950 }
951 if (act_proj && act_proj->func && (!selective || act_proj == n_proj)) {
952 act_proj->func(vec, el_info, bary2_2d[N_VERTICES_2D + i]);
953 }
954 }
955 }
956 } TRAVERSE_NEXT();
957 }
958
959 /* In order to obtain optimal error estimates for higher-order
960 * boundary approximations we need to be careful with the choice of
961 * the interior interpolation points. We do this only for strategy !=
962 * PARAM_ALL.
963 *
964 * We follow the procedure described in: M. Lenoir, Optimal
965 * Isoparametric Finite Elements and Error Estimates for Domains
966 * Involving Curved Bounaries, SIAM JNA 23(3) 1986.
967 *
968 * The procedure probably only makes sense when one and only one edge
969 * is projected, but we try to be a little bit more general here.
970 *
971 * This function assumes that center nodes are ordered this way:
972 *
973 * - the nodes are ordered in rows, row 0 corresponds to the smallest
974 * value of lambda[2].
975 */
adjust_center_nodesY_2d(const DOF * dof,const REAL_B * nodes,REAL_D * vec,NODE_PROJECTION ** edge_pr_loc,int n_e,int n_c)976 static void adjust_center_nodesY_2d(const DOF *dof, const REAL_B *nodes,
977 REAL_D *vec, NODE_PROJECTION **edge_pr_loc,
978 int n_e, int n_c)
979 {
980 int i, j, n_edge_pr;
981 REAL scale;
982
983 n_edge_pr = 0;
984 for (i = 0; i < N_WALLS_2D; i++) {
985 if (edge_pr_loc[i] != NULL) { /* This edge is curved */
986 n_edge_pr++;
987 }
988 }
989 scale = 1.0/(REAL)n_edge_pr;
990
991 for (i = 0; i < N_WALLS_2D; i++) {
992 if (edge_pr_loc[i] != NULL) { /* This edge is curved */
993 REAL_B lambda;
994 REAL_D d;
995 int v1, v2;
996
997 lambda[i] = 0.0;
998 v1 = vertex_of_edge_2d[i][0];
999 v2 = vertex_of_edge_2d[i][1];
1000
1001 for (j = 0; j < n_c; j++) {
1002 int ldof = N_VERTICES_2D+N_EDGES_2D*n_e+j;
1003 int e_ldof;
1004
1005 lambda[v1] = 1.0 - nodes[ldof][v2];
1006 lambda[v2] = nodes[ldof][v2];
1007
1008 /* This is actually a Lagrange-node on the wall, determine its
1009 * number.
1010 */
1011 e_ldof = (int)(lambda[v2]*(REAL)(n_e + 1) + 0.5) - 1;
1012 e_ldof += N_VERTICES_2D + n_e*i;
1013 AXPBY_DOW(lambda[v1], vec[dof[v1]],
1014 lambda[v2], vec[dof[v2]], d);
1015 AXPY_DOW(-1.0, vec[dof[e_ldof]], d);
1016 AXPY_DOW(-scale*0.5*nodes[ldof][v1]/lambda[v1], d, vec[dof[ldof]]);
1017
1018 /* Symmetrize a little bit */
1019 lambda[v1] = nodes[ldof][v1];
1020 lambda[v2] = 1.0 - nodes[ldof][v1];
1021
1022 /* This is actually a Lagrange-node on the wall,
1023 * determine its number.
1024 */
1025 e_ldof = n_e - 1 - (int)(lambda[v1]*(REAL)(n_e + 1) + 0.5) + 1;
1026 e_ldof += N_VERTICES_2D + n_e*i;
1027 AXPBY_DOW(lambda[v1], vec[dof[v1]],
1028 lambda[v2], vec[dof[v2]], d);
1029 AXPY_DOW(-1.0, vec[dof[e_ldof]], d);
1030 AXPY_DOW(-scale*0.5*nodes[ldof][v2]/lambda[v2], d, vec[dof[ldof]]);
1031 }
1032 }
1033 }
1034 }
1035
child_nodes_2d(int degree)1036 static const REAL_B *const*child_nodes_2d(int degree)
1037 {
1038 static const REAL_B child_loc[2][N_VERTICES_2D] = {
1039 { INIT_BARY_2D(0.0, 0.0, 1.0),
1040 INIT_BARY_2D(1.0, 0.0, 0.0),
1041 INIT_BARY_2D(0.5, 0.5, 0.0) },
1042 { INIT_BARY_2D(0.0, 1.0, 0.0),
1043 INIT_BARY_2D(0.0, 0.0, 1.0),
1044 INIT_BARY_2D(0.5, 0.5, 0.0) },
1045 };
1046 static REAL_B *(*child_nodes)[2];
1047 static int degree_max;
1048 int i, j, k, d;
1049
1050 if (!child_nodes) {
1051 child_nodes = (REAL_B *(*)[2])MEM_ALLOC(2*(degree+1), REAL_B *);
1052 } else if (degree > degree_max) {
1053 child_nodes = (REAL_B *(*)[2])
1054 MEM_REALLOC(child_nodes, 3*(degree_max+1), 2*(degree+1), REAL_B *);
1055 }
1056
1057 /* we generate all child nodes up to degree_max */
1058 if (degree_max < degree) {
1059 for (d = MAX(degree_max, 1); d <= degree; d++) {
1060 const BAS_FCTS *bas_fcts = get_lagrange(MESH_DIM, d);
1061 const REAL_B *nodes = LAGRANGE_NODES(bas_fcts);
1062 int n_bas_fcts;
1063
1064 n_bas_fcts = (d+1)*(d+2)/2;
1065 child_nodes[d][0] = MEM_ALLOC(n_bas_fcts, REAL_B);
1066 child_nodes[d][1] = MEM_ALLOC(n_bas_fcts, REAL_B);
1067
1068 for (i = 0; i < n_bas_fcts; i++) {
1069 for (j = 0; j < 2; j++) {
1070 AXEY_BAR(MESH_DIM, nodes[i][0], child_loc[j][0],
1071 child_nodes[d][j][i]);
1072 for (k = 1; k < N_VERTICES_2D; k++) {
1073 AXPY_BAR(MESH_DIM, nodes[i][k], child_loc[j][k],
1074 child_nodes[d][j][i]);
1075 }
1076 }
1077 }
1078 }
1079 degree_max = degree;
1080 }
1081 return (const REAL_B *const*)child_nodes[degree];
1082 }
1083
refine_interpolY_2d(DOF_REAL_D_VEC * drdv,RC_LIST_EL * list,int n)1084 static void refine_interpolY_2d(DOF_REAL_D_VEC *drdv, RC_LIST_EL *list, int n)
1085 {
1086 const FE_SPACE *fe_space = drdv->fe_space;
1087 MESH *mesh = fe_space->mesh;
1088 const BAS_FCTS *bas_fcts = fe_space->bas_fcts;
1089 const DOF_ADMIN *admin = fe_space->admin;
1090 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data;
1091 PARAM_STRATEGY strategy = data->strategy;
1092 DOF_PTR_VEC *edge_pr = data->edge_projections;
1093 NODE_PROJECTION *n_proj = data->n_proj;
1094 bool selective = n_proj != NULL;
1095 const NODE_PROJECTION *act_proj;
1096 REAL_D *vec = drdv->vec;
1097 int node_v, node_e, n0_v, n_c, n_e;
1098 const EL_INFO *el_info;
1099 EL *el;
1100 const REAL_B *nodes = LAGRANGE_NODES(bas_fcts);
1101 const REAL_B *const*child_nodes;
1102 int ref_edge_untouched = false;
1103 int n0_edge_pr = -1;
1104 int i, j, elem, i_child;
1105
1106 TEST_EXIT(bas_fcts->degree <= DEGREE_MAX,
1107 "degree %d unsupported, only up to degree %d.\n",
1108 bas_fcts->degree, DEGREE_MAX);
1109
1110 child_nodes = child_nodes_2d(bas_fcts->degree);
1111
1112 node_e = mesh->node[EDGE];
1113 node_v = mesh->node[VERTEX];
1114 n0_v = admin->n0_dof[VERTEX];
1115 n_e = admin->n_dof[EDGE];
1116 n_c = admin->n_dof[CENTER];
1117
1118 if (edge_pr) {
1119 n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
1120 }
1121
1122 /* If strategy != PARAM_STRAIGHT_CHILDS first call the default
1123 * interpolation routine for the coordinates. Why should we
1124 * duplicate work. Otherwise we interpolate between the vertices
1125 * according to the Lagrange-nodes.
1126 */
1127 if (strategy != PARAM_STRAIGHT_CHILDS) {
1128 bas_fcts->real_d_refine_inter(drdv, list, n);
1129 }
1130
1131 for (elem = 0; elem < n; elem++) {
1132 el_info = &list[elem].el_info;
1133 el = el_info->el;
1134
1135 for (i_child = 0; i_child < 2; i_child++) {
1136 const REAL_B *cnodes;
1137 DOF dof[N_BAS_MAX];
1138
1139 cnodes = child_nodes[i_child];
1140
1141 /* Get the DOF-indices on the child */
1142 GET_DOF_INDICES(bas_fcts, el->child[i_child], admin, dof);
1143
1144 /* handle the refinement edge, once only */
1145 if (elem == 0) {
1146 /* <<< refinement edge stuff */
1147
1148 if (i_child == 0) {
1149
1150 /* Compute coordinates if not already done */
1151 if (strategy == PARAM_STRAIGHT_CHILDS) {
1152 /* vertex 2 coordinate */
1153 AXPBY_DOW(0.5, vec[el->dof[node_v+0][n0_v]],
1154 0.5, vec[el->dof[node_v+1][n0_v]], vec[dof[2]]);
1155 }
1156
1157 /* Update edge_projections and do the necessary
1158 * projections. The children just inherit the
1159 * edge_projections value from the parent.
1160 */
1161 if (edge_pr) {
1162 DOF pdof = el->dof[node_e+2][n0_edge_pr];
1163 ref_edge_untouched = edge_pr->vec[pdof] == NULL;
1164 }
1165
1166 act_proj = list->el_info.active_projection;
1167
1168 if (i_child == 0 &&
1169 !ref_edge_untouched &&
1170 act_proj && act_proj->func &&
1171 (!selective || act_proj == n_proj)) {
1172 /* Project the midpoint of the refinement edge */
1173 act_proj->func(vec[dof[2]], el_info, mid_lambda_2d);
1174 _AI_refine_update_bbox(mesh, vec[dof[2]]);
1175 }
1176 }
1177
1178 /* Compute coordinates if not already done. Do this after
1179 * (potentially) projecting the refinement edge.
1180 */
1181 if (strategy == PARAM_STRAIGHT_CHILDS) {
1182 for (i = 0; i < n_e; i++) {
1183 int local_dof = N_VERTICES_2D+i_child*n_e+i;
1184
1185 AXPBY_DOW(nodes[local_dof][1-i_child], vec[dof[1-i_child]],
1186 nodes[local_dof][2], vec[dof[2]],
1187 vec[dof[local_dof]]);
1188 }
1189 }
1190
1191 /* Apply the necessary projections */
1192 act_proj = list->el_info.active_projection;
1193
1194 if (!ref_edge_untouched &&
1195 act_proj && (!selective || n_proj == act_proj)) {
1196 if (act_proj->func) {
1197 /* Project all other edge DOFs */
1198 for (i = 0; i < n_e; i++) {
1199 int ldof = N_VERTICES_2D+i_child*n_e+i;
1200
1201 act_proj->func(vec[dof[ldof]], el_info, cnodes[ldof]);
1202 _AI_refine_update_bbox(mesh, vec[dof[ldof]]);
1203 }
1204 }
1205
1206 if (edge_pr) {
1207 DOF edge_pr_dof = el->child[i_child]->dof[node_e+i_child][n0_edge_pr];
1208 edge_pr->vec[edge_pr_dof] = (void *)act_proj;
1209 }
1210 } else if (strategy != PARAM_ALL) {
1211 DOF edge_pr_dof = el->child[i_child]->dof[node_e+i_child][n0_edge_pr];
1212 edge_pr->vec[edge_pr_dof] = NULL;
1213 }
1214 /* >>> */
1215 }
1216
1217 /* handle the common edge between the two children, do this only
1218 * once for each parent
1219 */
1220 if (i_child == 0 /* && n_e */) {
1221 /* <<< interior edge */
1222
1223 /* Compute coordinates if not already done */
1224 if (data->strategy == PARAM_STRAIGHT_CHILDS) {
1225 for (i = 0; i < n_e; i++) {
1226 int local_dof = N_VERTICES_2D+(1-i_child)*n_e+i;
1227
1228 AXPBY_DOW(nodes[local_dof][i_child], vec[dof[i_child]],
1229 nodes[local_dof][2], vec[dof[2]],
1230 vec[dof[local_dof]]);
1231 }
1232 }
1233
1234 act_proj = wall_proj(el_info, -1);
1235 if (!ref_edge_untouched &&
1236 act_proj && (!selective || n_proj == act_proj)) {
1237 if (act_proj->func) {
1238 for (i = 0; i < n_e; i++) {
1239 int local_dof = N_VERTICES_2D+(1-i_child)*n_e+i;
1240
1241 act_proj->func(vec[dof[local_dof]], el_info, cnodes[local_dof]);
1242 _AI_refine_update_bbox(mesh, vec[dof[local_dof]]);
1243 }
1244 }
1245
1246 if (edge_pr) {
1247 DOF edge_pr_dof = el->child[0]->dof[node_e+(1-i_child)][n0_edge_pr];
1248 edge_pr->vec[edge_pr_dof] = (void *)act_proj;
1249 }
1250 } else if (edge_pr) {
1251 DOF edge_pr_dof = el->child[0]->dof[node_e+(1-i_child)][n0_edge_pr];
1252 edge_pr->vec[edge_pr_dof] = NULL;
1253 }
1254
1255 /* >>> */
1256 }
1257
1258 /* then the interior DOFs, if any */
1259 if (n_c) {
1260 if (data->strategy == PARAM_STRAIGHT_CHILDS) {
1261 NODE_PROJECTION *edge_pr_loc[N_EDGES_2D];
1262
1263 /* First compute default coordinates if not already done */
1264 for (i = 0; i < n_c; i++) {
1265 int local_dof = N_VERTICES_2D+N_EDGES_2D*n_e+i;
1266
1267 AXEY_DOW(nodes[local_dof][0], vec[dof[0]], vec[dof[local_dof]]);
1268 for (j = 1; j < N_VERTICES_2D; j++) {
1269 AXPY_DOW(nodes[local_dof][j], vec[dof[j]], vec[dof[local_dof]]);
1270 }
1271 }
1272
1273 /* Then add a correction which guarantees that the higher
1274 * derivative of the transformation to the reference element
1275 * vanish fast enough.
1276 */
1277 for (i = 0; i < N_EDGES_2D; i++) {
1278 edge_pr_loc[i] =
1279 edge_pr->vec[el->child[i_child]->dof[node_e+i][n0_edge_pr]];
1280 }
1281
1282 adjust_center_nodesY_2d(dof, nodes, vec, edge_pr_loc, n_e, n_c);
1283 } /* !PARAM_ALL */
1284
1285 act_proj = wall_proj(el_info, -1);
1286
1287 if (!ref_edge_untouched &&
1288 act_proj && act_proj->func && (!selective || n_proj == act_proj)){
1289 for (i = 0; i < n_c; i++) {
1290 int local_dof = N_VERTICES_2D+N_EDGES_2D*n_e+i;
1291 act_proj->func(vec[dof[local_dof]], el_info, cnodes[local_dof]);
1292 }
1293 }
1294 }
1295 } /* child loop */
1296 } /* element loop */
1297
1298 /* Co-ordinates are already correct when the parent was affine
1299 * or with strategy == PARAM_STRAIGHT_CHILDS.
1300 *
1301 * We have to adjust all edge and center DOFs. Even though some of
1302 * those DOFs were not projected, they were (potentially)
1303 * interpolated by the default real_d_refine_inter(). The effects
1304 * are undone here.
1305 *
1306 * We have to do the correction after doing anything else.
1307 */
1308 if (ref_edge_untouched && data->strategy == PARAM_CURVED_CHILDS) {
1309
1310 for (elem = 0; elem < n; elem++) {
1311 int parent_affine = true;
1312
1313 el_info = &list[elem].el_info;
1314 el = el_info->el;
1315
1316 for (j = 0; j < N_EDGES_2D; j++) {
1317 DOF edge_dof = el->dof[node_e+j][n0_edge_pr];
1318
1319 if (edge_pr->vec[edge_dof]) {
1320 parent_affine = false;
1321 break;
1322 }
1323 }
1324
1325 if (!parent_affine) {
1326 for (i_child = 0; i_child < 2; i_child++) {
1327 EL *child = el->child[i_child];
1328 int child_affine = true;
1329
1330 /* The child is untouched iff all edges are untouched. */
1331 for (j = 0; j < N_EDGES_2D; j++) {
1332 DOF edge_dof = child->dof[node_e+j][n0_edge_pr];
1333
1334 if (edge_pr->vec[edge_dof]) {
1335 child_affine = false;
1336 break;
1337 }
1338 }
1339 if (child_affine) {
1340 DOF dof[N_BAS_MAX];
1341
1342 /* Get the DOF-indices on the child */
1343 GET_DOF_INDICES(bas_fcts, child, admin, dof);
1344
1345 /* Adjust all new co-ordinates. */
1346 if (n_e > 0) {
1347 for (i = 0; i < N_EDGES_2D; i++) {
1348 const int *voe = vertex_of_edge_2d[i];
1349 REAL *vec0, *vec1;
1350
1351 vec0 = vec[dof[voe[0]]];
1352 vec1 = vec[dof[voe[1]]];
1353 for (j = 0; j < n_e; j++) {
1354 int ldof = N_VERTICES_2D + i*n_e + j;
1355
1356 AXPBY_DOW(nodes[ldof][voe[0]], vec0,
1357 nodes[ldof][voe[1]], vec1,
1358 vec[dof[ldof]]);
1359 }
1360 }
1361 }
1362 for (i = 0; i < n_c; i++) {
1363 int ldof = N_VERTICES_2D + N_EDGES_2D*n_e + i;
1364
1365 AXEY_DOW(nodes[ldof][0], vec[dof[0]], vec[dof[ldof]]);
1366 for (j = 1; j < N_VERTICES_2D; j++) {
1367 AXPY_DOW(nodes[ldof][j], vec[dof[j]], vec[dof[ldof]]);
1368 }
1369 }
1370 }
1371 } /* child loop */
1372 } /* I am not affine */
1373 } /* element loop */
1374 } /* child correction necessary */
1375 }
1376
coarse_interpolY_2d(DOF_REAL_D_VEC * drdv,RC_LIST_EL * list,int n)1377 static void coarse_interpolY_2d(DOF_REAL_D_VEC *drdv, RC_LIST_EL *list, int n)
1378 {
1379 const FE_SPACE *fe_space = drdv->fe_space;
1380 const MESH *mesh = fe_space->mesh;
1381 const BAS_FCTS *bas_fcts = fe_space->bas_fcts;
1382 LAGRANGE_PARAM_DATA *data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data;
1383 DOF_PTR_VEC *edge_pr = data->edge_projections;
1384
1385 /* first call the default interpolation routine for the coordinates */
1386 bas_fcts->real_d_coarse_inter(drdv, list, n);
1387
1388 /* then take care of edge_projections */
1389 if (edge_pr) {
1390 int n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
1391 EL *el = list->el_info.el;
1392 int node_e = mesh->node[EDGE];
1393 DOF pdof = el->dof[node_e+2][n0_edge_pr];
1394 DOF cdof = el->child[0]->dof[node_e+0][n0_edge_pr];
1395
1396 edge_pr->vec[pdof] = edge_pr->vec[cdof];
1397
1398 /* BIG FAT TODO: center nodes have to be re-adjusted after
1399 * coarsening
1400 */
1401 }
1402
1403 return;
1404 }
1405
fill_coordsY_2d(LAGRANGE_PARAM_DATA * data)1406 static void fill_coordsY_2d(LAGRANGE_PARAM_DATA *data)
1407 {
1408 DOF_REAL_D_VEC *coords = data->coords;
1409 DOF_PTR_VEC *edge_pr = data->edge_projections;
1410 NODE_PROJECTION *n_proj = data->n_proj;
1411 bool selective = n_proj != NULL;
1412 const NODE_PROJECTION *act_proj;
1413 FLAGS fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_PROJECTION;
1414 MESH *mesh = coords->fe_space->mesh;
1415 const DOF_ADMIN *admin = coords->fe_space->admin;
1416 const BAS_FCTS *bas_fcts = coords->fe_space->bas_fcts;
1417 int i, j, k, n_e, n_c, node_e = -1, n0_edge_pr = -1;
1418 const REAL_B *nodes = LAGRANGE_NODES(bas_fcts);
1419
1420 n_e = admin->n_dof[EDGE];
1421 n_c = admin->n_dof[CENTER];
1422
1423 if (edge_pr) {
1424 node_e = mesh->node[EDGE];
1425 n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
1426 }
1427
1428 TRAVERSE_FIRST(mesh, -1, fill_flag|FILL_BOUND) {
1429 DOF dof[N_BAS_MAX];
1430 REAL *vec, *vec0, *vec1;
1431
1432 GET_DOF_INDICES(bas_fcts, el_info->el, admin, dof);
1433
1434 for (i = 0; i < N_VERTICES_2D; i++) {
1435 vec = coords->vec[dof[i]];
1436
1437 COPY_DOW(el_info->coord[i], vec);
1438
1439 /* Look for a projection function that applies to vertex[i].
1440 * Apply this projection if found.
1441 */
1442 if (!selective || n_proj->func) {
1443 act_proj = wall_proj(el_info, (i+1) % N_WALLS_2D);
1444
1445 if (!act_proj) {
1446 act_proj = wall_proj(el_info, (i+2) % N_WALLS_2D);
1447 }
1448 if (!act_proj) {
1449 act_proj = wall_proj(el_info, -1);
1450 }
1451
1452 if (act_proj && act_proj->func && (!selective || act_proj == n_proj)) {
1453 act_proj->func(vec, el_info, nodes[i]);
1454 }
1455 }
1456 }
1457
1458 if (true /* n_e */) {
1459 for (i = 0; i < N_EDGES_2D; i++) {
1460 const int *voe = vertex_of_edge_2d[i];
1461
1462 vec0 = coords->vec[dof[voe[0]]];
1463 vec1 = coords->vec[dof[voe[1]]];
1464
1465 act_proj = wall_proj(el_info, i);
1466 if (!act_proj) {
1467 act_proj = wall_proj(el_info, -1);
1468 }
1469
1470 if (edge_pr) {
1471 DOF edge_dof = el_info->el->dof[node_e+i][n0_edge_pr];
1472 edge_pr->vec[edge_dof] =
1473 act_proj && (!selective || act_proj == n_proj)
1474 ? (void *)act_proj : NULL;
1475 }
1476
1477 for (j = 0; j < n_e; j++) {
1478 int local_dof = N_VERTICES_2D + i*n_e + j;
1479
1480 vec = coords->vec[dof[local_dof]];
1481
1482 AXPBY_DOW(nodes[local_dof][voe[0]], vec0,
1483 nodes[local_dof][voe[1]], vec1,
1484 vec);
1485
1486 if (act_proj && act_proj->func &&
1487 (!selective || act_proj == n_proj)) {
1488 act_proj->func(vec, el_info, nodes[local_dof]);
1489 }
1490 }
1491 }
1492
1493 if (n_c) {
1494 act_proj = wall_proj(el_info, -1);
1495
1496 for (j = 0; j < n_c; j++) {
1497 int local_dof = N_VERTICES_2D + N_EDGES_2D*n_e + j;
1498
1499 vec = coords->vec[dof[local_dof]];
1500
1501 AXEY_DOW(nodes[local_dof][0], coords->vec[dof[0]], vec);
1502 for (k = 1; k < N_VERTICES_2D; k++) {
1503 AXPY_DOW(nodes[local_dof][k], coords->vec[dof[k]], vec);
1504 }
1505 }
1506
1507 if (data->strategy != PARAM_ALL) {
1508 /* Add a correction which guarantees that the 2nd
1509 * derivative of the transformation to the reference element
1510 * vanishes fast enough.
1511 */
1512 NODE_PROJECTION *edge_pr_loc[N_EDGES_2D];
1513
1514 for (i = 0; i < N_EDGES_2D; i++) {
1515 edge_pr_loc[i] =
1516 edge_pr->vec[el_info->el->dof[node_e+i][n0_edge_pr]];
1517 }
1518
1519 adjust_center_nodesY_2d(dof,
1520 nodes, coords->vec, edge_pr_loc, n_e, n_c);
1521 }
1522
1523 for (j = 0; j < n_c; j++) {
1524 int local_dof = N_VERTICES_2D + N_EDGES_2D*n_e + j;
1525
1526 vec = coords->vec[dof[local_dof]];
1527
1528 if (act_proj && act_proj->func &&
1529 (!selective || act_proj == n_proj)) {
1530 act_proj->func(vec, el_info, nodes[local_dof]);
1531 }
1532 }
1533 }
1534 }
1535 } TRAVERSE_NEXT();
1536
1537 return;
1538 }
1539
1540 #if DIM_MAX > 2
1541 static void
slave_refine_interpolY_2d(DOF_REAL_D_VEC * s_coords,RC_LIST_EL * list,int n)1542 slave_refine_interpolY_2d(DOF_REAL_D_VEC *s_coords, RC_LIST_EL *list, int n)
1543 {
1544 const FE_SPACE *fe_space = s_coords->fe_space;
1545 MESH *slave = fe_space->mesh;
1546 const BAS_FCTS *s_bfcts = fe_space->bas_fcts;
1547 const DOF_ADMIN *s_admin = fe_space->admin;
1548 PARAMETRIC *s_parametric = slave->parametric;
1549 LAGRANGE_PARAM_DATA *s_data = (LAGRANGE_PARAM_DATA *)s_parametric->data;
1550 MESH *master = get_master(slave);
1551 DOF_REAL_D_VEC *m_coords =
1552 ((LAGRANGE_PARAM_DATA *)master->parametric->data)->coords;
1553 const BAS_FCTS *m_bfcts = m_coords->fe_space->bas_fcts;
1554 const DOF_ADMIN *m_admin = m_coords->fe_space->admin;
1555 DOF_PTR_VEC *s_edge_pr;
1556 DOF_PTR_VEC *m_edge_pr = NULL;
1557 const EL_INFO *el_info;
1558 EL *m_el, *s_el, *s_child, *m_child;
1559 int wall, orientation, type;
1560 int p_wall, p_orientation, p_type;
1561 int node_e = -1;
1562 int n_c, n_e;
1563 int m_node_e = -1;
1564 int m_n_e = -1;
1565 int n0_edge_pr = -1, m_n0_edge_pr = -1;
1566 int s_edge, m_edge, s_dof, m_dof, s_dof_loc, m_dof_loc;
1567 int elem, i_child, i, m_i_child;
1568 DOF s_dofs[s_bfcts->n_bas_fcts];
1569 DOF m_dofs[m_bfcts->n_bas_fcts];
1570 const int *trace_map;
1571
1572 if ((s_edge_pr = s_data->edge_projections)) {
1573 node_e = slave->node[EDGE];
1574 m_node_e = master->node[EDGE];
1575 n0_edge_pr = s_edge_pr->fe_space->admin->n0_dof[EDGE];
1576 m_edge_pr =
1577 ((LAGRANGE_PARAM_DATA *)master->parametric->data)->edge_projections;
1578 m_n0_edge_pr = m_edge_pr->fe_space->admin->n0_dof[EDGE];
1579 m_n_e = m_admin->n_dof[EDGE];
1580 }
1581
1582 n_e = s_admin->n_dof[EDGE];
1583 n_c = s_admin->n_dof[CENTER];
1584
1585 for (elem = 0; elem < n; elem++) {
1586
1587 el_info = &list[elem].el_info;
1588 s_el = el_info->el;
1589
1590 /* slave -> master binding has not been set up yet, so duplicate
1591 * the logic from submesh_3d.c here, given the master element,
1592 * orientation and type
1593 */
1594 m_el = el_info->master.el;
1595 p_wall = el_info->master.opp_vertex;
1596 p_orientation = el_info->master.orientation;
1597 p_type = el_info->master.el_type;
1598
1599 for (i_child = 0; i_child < 2; i_child++) {
1600
1601 if (p_wall == 2) {
1602 m_i_child = p_orientation > 0 ? i_child : 1 - i_child;
1603 } else {
1604 m_i_child = p_orientation < 0 ? i_child : 1 - i_child;
1605 }
1606
1607 s_child = s_el->child[i_child];
1608 m_child = m_el->child[m_i_child];
1609
1610 wall = child_face_3d[p_type][p_wall][m_i_child];
1611 orientation = p_orientation*child_orientation_3d[p_type][m_i_child];
1612 type = (p_type + 1) % 3;
1613 trace_map = m_bfcts->trace_dof_map[type > 0][orientation < 0][wall];
1614
1615 GET_DOF_INDICES(s_bfcts, s_child, s_admin, s_dofs);
1616 GET_DOF_INDICES(m_bfcts, m_child, m_admin, m_dofs);
1617
1618 if (i_child == 0) {
1619 /* New vertex */
1620 s_dof = s_dofs[2];
1621 m_dof = m_dofs[trace_map[2]];
1622 COPY_DOW(m_coords->vec[m_dof], s_coords->vec[s_dof]);
1623
1624 /* common wall between the two children */
1625 s_edge = 1-i_child;
1626 for (i = 0; i < n_e; i++) {
1627 s_dof_loc = N_VERTICES_2D + s_edge*n_e + i;
1628 m_dof_loc = trace_map[s_dof_loc];
1629 s_dof = s_dofs[s_dof_loc];
1630 m_dof = m_dofs[m_dof_loc];
1631 COPY_DOW(m_coords->vec[m_dof], s_coords->vec[s_dof]);
1632
1633 if (s_edge_pr) {
1634 s_dof = s_child->dof[node_e + s_edge][n0_edge_pr];
1635 m_edge = (m_dof_loc - N_VERTICES_3D) / m_n_e;
1636 m_dof = m_child->dof[m_node_e + m_edge][m_n0_edge_pr];
1637 s_edge_pr->vec[s_dof] = m_edge_pr->vec[m_dof];
1638 }
1639 }
1640 }
1641
1642 if (elem == 0) {
1643 /* wall to neighbour (or boundary) */
1644 s_edge = i_child;
1645 for (i = 0; i < n_e; i++) {
1646 s_dof_loc = N_VERTICES_2D + s_edge*n_e + i;
1647 m_dof_loc = trace_map[s_dof_loc];
1648 s_dof = s_dofs[s_dof_loc];
1649 m_dof = m_dofs[m_dof_loc];
1650 COPY_DOW(m_coords->vec[m_dof], s_coords->vec[s_dof]);
1651
1652 if (s_edge_pr) {
1653 s_dof = s_child->dof[node_e + s_edge][n0_edge_pr];
1654 m_edge = (m_dof_loc - N_VERTICES_3D) / m_n_e;
1655 m_dof = m_child->dof[m_node_e + m_edge][m_n0_edge_pr];
1656 s_edge_pr->vec[s_dof] = m_edge_pr->vec[m_dof];
1657 }
1658 }
1659 }
1660
1661 /* interior DOFs (face-dofs on the master) */
1662 for (i = 0; i < n_c; i++) {
1663 s_dof_loc = N_VERTICES_2D + N_EDGES_2D*n_e + i;
1664 m_dof_loc = trace_map[s_dof_loc];
1665 s_dof = s_dofs[s_dof_loc];
1666 m_dof = m_dofs[m_dof_loc];
1667 COPY_DOW(m_coords->vec[m_dof], s_coords->vec[s_dof]);
1668 }
1669 }
1670 }
1671
1672 return;
1673 }
1674
slave_fill_coordsY_2d(LAGRANGE_PARAM_DATA * s_data)1675 static void slave_fill_coordsY_2d(LAGRANGE_PARAM_DATA *s_data)
1676 {
1677 DOF_REAL_D_VEC *s_coords = s_data->coords;
1678 DOF_PTR_VEC *s_edge_pr = s_data->edge_projections;
1679 MESH *slave = s_coords->fe_space->mesh;
1680 const BAS_FCTS *s_bfcts = s_coords->fe_space->bas_fcts;
1681 const DOF_ADMIN *s_admin = s_coords->fe_space->admin;
1682 MESH *master = get_master(slave);
1683 PARAMETRIC *m_parametric = master->parametric;
1684 LAGRANGE_PARAM_DATA *m_data = (LAGRANGE_PARAM_DATA *)m_parametric->data;
1685 DOF_REAL_D_VEC *m_coords = m_data->coords;
1686 DOF_PTR_VEC *m_edge_pr = m_data->edge_projections;
1687 const BAS_FCTS *m_bfcts = m_coords->fe_space->bas_fcts;
1688 const DOF_ADMIN *m_admin = m_coords->fe_space->admin;
1689 int n0_edge_pr =
1690 s_edge_pr ? s_edge_pr->fe_space->admin->n0_dof[EDGE] : -1;
1691 int m_n0_edge_pr =
1692 m_edge_pr ? m_edge_pr->fe_space->admin->n0_dof[EDGE] : -1;
1693 EL *m_el, *s_el;
1694 int wall, orientation, type;
1695 int node_e;
1696 int m_node_e;
1697 int m_n_e, n_e;
1698 int s_edge, m_edge, s_dof, m_dof, s_dof_loc, m_dof_loc, i;
1699 const int *trace_map;
1700 DOF m_dofs[m_bfcts->n_bas_fcts];
1701 DOF s_dofs[s_bfcts->n_bas_fcts];
1702
1703 n_e = s_admin->n_dof[EDGE];
1704 node_e = slave->node[EDGE];
1705
1706 m_n_e = m_admin->n_dof[EDGE];
1707 m_node_e = master->node[EDGE];
1708
1709 TRAVERSE_FIRST(slave, -1, CALL_LEAF_EL|FILL_MASTER_INFO) {
1710
1711 s_el = el_info->el;
1712 m_el = el_info->master.el;
1713 wall = el_info->master.opp_vertex;
1714 orientation = el_info->master.orientation;
1715 type = el_info->master.el_type;
1716
1717 trace_map = m_bfcts->trace_dof_map[type > 0][orientation < 0][wall];
1718
1719 GET_DOF_INDICES(s_bfcts, s_el, s_admin, s_dofs);
1720 GET_DOF_INDICES(m_bfcts, m_el, m_admin, m_dofs);
1721
1722 for (i = 0; i < s_data->n_local_coords; i++) {
1723 COPY_DOW(m_coords->vec[m_dofs[trace_map[i]]], s_coords->vec[s_dofs[i]]);
1724 }
1725
1726 if (s_edge_pr) {
1727 for (s_edge = 0; s_edge < N_EDGES_2D; s_edge++) {
1728 /* first determine the matching edge on the master */
1729 s_dof_loc = N_VERTICES_2D + s_edge*n_e;
1730 m_dof_loc = trace_map[s_dof_loc];
1731 m_edge = (m_dof_loc - N_VERTICES_3D) / m_n_e;
1732
1733 /* Copy over the projeted edges status */
1734 s_dof = s_el->dof[node_e+s_edge][n0_edge_pr];
1735 m_dof = m_el->dof[m_node_e+m_edge][m_n0_edge_pr];
1736 s_edge_pr->vec[s_dof] = m_edge_pr->vec[m_dof];
1737 }
1738 }
1739 } TRAVERSE_NEXT();
1740
1741 return;
1742 }
1743 #endif
1744
1745 static const PARAMETRIC lagrange_parametric1_2d = {
1746 "2D Lagrange parametric elements of degree 1",
1747 false, /* not_all */
1748 false, /* use_reference_mesh */
1749 param_init_element1_2d,
1750 vertex_coordsY_2d,
1751 param_coord_to_world,
1752 param_world_to_coord,
1753 det1_2d,
1754 grd_lambda1_2d,
1755 NULL,
1756 wall_normal1_2d,
1757 NULL /* data */
1758 };
1759
1760 static const PARAMETRIC lagrange_parametric2_2d = {
1761 "2D Lagrange parametric elements of degree 2",
1762 false, /* not_all */
1763 false, /* use_reference_mesh */
1764 param_init_elementY_2d,
1765 vertex_coordsY_2d,
1766 param_coord_to_world,
1767 param_world_to_coord,
1768 detY_2d,
1769 grd_lambdaY_2d,
1770 grd_worldY_2d,
1771 wall_normalY_2d,
1772 NULL /* data */
1773 };
1774
1775 static const PARAMETRIC lagrange_parametricY_2d = {
1776 "2D Lagrange parametric elements of generic degree",
1777 false, /* not_all */
1778 false, /* use_reference_mesh */
1779 param_init_elementY_2d,
1780 vertex_coordsY_2d,
1781 param_coord_to_world,
1782 param_world_to_coord,
1783 detY_2d,
1784 grd_lambdaY_2d,
1785 grd_worldY_2d,
1786 wall_normalY_2d,
1787 NULL /* data */
1788 };
1789