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: parametric.c
7 *
8 *
9 * description: Support for parametric elements.
10 *
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 * Numerik fuer Hoechstleistungsrechner
21 * Universit�t 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
32 #include "alberta_intern.h"
33 #include "alberta.h"
34
35 #define LAGRANGE_MAGIC { 'L', 'P', 'A', 'R' }
36 #define MAGIC_STRING "LPAR"
37
38 #define DEGREE_GEN 3 /* minimal polynomial degree with generic implementation */
39 #define DEGREE_MAX 4 /* maximal polynomial degree supported */
40
41 typedef struct lagrange_param_data LAGRANGE_PARAM_DATA;
42 struct lagrange_param_data
43 {
44 char magic[4];
45 int degree;
46 PARAM_STRATEGY strategy;
47 NODE_PROJECTION *n_proj;
48 DOF_REAL_D_VEC *coords;
49
50 DOF_PTR_VEC *edge_projections; /* The projections for each edge
51 * or NULL if this edge is not
52 * subject to a projection.
53 */
54
55 REAL_D *local_coords; /* points to either
56 * LAGRANGE_PARAM_DATA::param_local_coords
57 * or to the current EL_INFO::coord
58 * space if the current element is
59 * affine and after
60 * PARAMETRIC::init_element() has
61 * been called.
62 */
63
64 REAL_D *param_local_coords; /* allocated space for the
65 * local coordinate vector,
66 * filled by
67 * PARAMETRIC::init_element().
68 */
69
70 int n_local_coords; /* == n_bas_fcts */
71 int i_am_affine; /* cached return value of
72 * PARAMETRIC::init_element()
73 * (actually its negation).
74 */
75
76 EL *el; /* pointer to the element from the last
77 * PARAMETRIC::init_element() invocation.
78 */
79
80 int max_iter; /* maximum number of iterations and */
81 REAL newton_tolerance; /* tolerance for newton's method in */
82 /* param_world_to_coord */
83 REAL lambda_tolerance; /* specify when lambda[k] is considered < 0 */
84 };
85
86 static void param_coord_to_world(const EL_INFO *el_info, const QUAD *quad,
87 int N, const REAL_B lambda[],
88 REAL_D *world);
89
90 static void param_world_to_coord(const EL_INFO *el_info, int N,
91 const REAL_D world[],
92 REAL_B lambda[], int k[]);
93
94 /* Interface to ALBERTA's Newton iterators. */
95 static void nl_update(void *data,
96 int dim, const REAL *lambda_k, bool upd_DF, REAL *Fx);
97 static int nl_solve(void *data, int dim, const REAL *rhs, REAL *d);
98 typedef struct isp_data
99 {
100 REAL_D x; /* world */
101 REAL_D *coords; /* coords of the quad points */
102 REAL Df[DIM_OF_WORLD+1][DIM_OF_WORLD+1];
103 BAS_FCTS *bfcts;
104 const EL_INFO *el_info;
105 } ISP_DATA;
106
107 #include "../0d/parametric_0d.c"
108 #if DIM_MAX > 0
109 #include "../1d/parametric_1d.c"
110 #if DIM_MAX > 1
111 #include "../2d/parametric_2d.c"
112 #if DIM_MAX > 2
113 #include "../3d/parametric_3d.c"
114 #endif
115 #endif
116 #endif
117
118 static const PARAMETRIC *all_lagrange[DIM_MAX+1][DEGREE_GEN] = {
119 { &lagrange_parametric_0d,
120 &lagrange_parametric_0d,
121 &lagrange_parametric_0d },
122 { &lagrange_parametric1_1d,
123 &lagrange_parametric2_1d,
124 &lagrange_parametricY_1d }
125 #if DIM_MAX > 1
126 ,
127 { &lagrange_parametric1_2d,
128 &lagrange_parametric2_2d,
129 &lagrange_parametricY_2d }
130 #if DIM_MAX > 2
131 ,
132 { &lagrange_parametric1_3d,
133 &lagrange_parametric2_3d,
134 &lagrange_parametricY_3d }
135 #endif
136 #endif
137 };
138
139 static void (*all_refine_interpol[DIM_MAX+1][DEGREE_GEN])
140 (DOF_REAL_D_VEC *, RC_LIST_EL *, int) = {
141 { NULL,
142 NULL,
143 NULL },
144 { refine_interpol1_1d,
145 refine_interpol2_1d,
146 refine_interpolY_1d }
147 #if DIM_MAX > 1
148 ,
149 { refine_interpol1_2d,
150 refine_interpol2_2d,
151 refine_interpolY_2d }
152 #if DIM_MAX > 2
153 ,
154 { refine_interpol1_3d,
155 refine_interpol2_3d,
156 refine_interpolY_3d }
157 #endif
158 #endif
159 };
160
161 static void (*all_coarse_interpol[DIM_MAX+1][DEGREE_GEN])
162 (DOF_REAL_D_VEC *, RC_LIST_EL *, int) = {
163 { NULL,
164 NULL,
165 NULL },
166 { NULL,
167 coarse_interpol2_1d,
168 coarse_interpolY_1d }
169 #if DIM_MAX > 1
170 ,
171 { NULL,
172 coarse_interpol2_2d,
173 coarse_interpolY_2d }
174 #if DIM_MAX > 2
175 ,
176 { NULL,
177 coarse_interpol2_3d,
178 coarse_interpolY_3d }
179 #endif
180 #endif
181 };
182
183 static void (*all_fill_coords[DIM_MAX+1][DEGREE_GEN])
184 (LAGRANGE_PARAM_DATA *data) = {
185 { fill_coords_0d,
186 fill_coords_0d,
187 fill_coords_0d },
188 { fill_coords1_1d,
189 fill_coords2_1d,
190 fill_coordsY_1d }
191 #if DIM_MAX > 1
192 ,
193 { fill_coords1_2d,
194 fill_coords2_2d,
195 fill_coordsY_2d }
196 #if DIM_MAX > 2
197 ,
198 { fill_coords1_3d,
199 fill_coords2_3d,
200 fill_coordsY_3d }
201 #endif
202 #endif
203 };
204
205 static void (*alloc_quad_metadata[DIM_MAX+1])(const QUAD *quad,
206 const BAS_FCTS *bas_fcts) = {
207 NULL,
208 alloc_param_quad_metadata_1d,
209 #if DIM_OF_WORLD > 1
210 alloc_param_quad_metadata_2d,
211 # if DIM_OF_WORLD > 2
212 alloc_param_quad_metadata_3d,
213 # endif
214 #endif
215 };
216
217 typedef struct param_quad_metadata
218 {
219 const QUAD_FAST *quad_fast;
220 } PARAM_QUAD_METADATA;
221
param_coord_to_world(const EL_INFO * el_info,const QUAD * quad,int N,const REAL_B lambda[],REAL_D * world)222 static void param_coord_to_world(const EL_INFO *el_info, const QUAD *quad,
223 int N, const REAL_B lambda[],
224 REAL_D *world)
225 {
226 int i, iq, dim;
227 LAGRANGE_PARAM_DATA *data =
228 (LAGRANGE_PARAM_DATA *)el_info->mesh->parametric->data;
229 REAL_D *local_coords = data->local_coords;
230 const BAS_FCTS *bas_fcts = data->coords->fe_space->bas_fcts;
231
232 if (quad) {
233 if (!data->i_am_affine) {
234 QUAD_METADATA *qmd = (QUAD_METADATA *)quad->metadata;
235 PARAM_QUAD_METADATA *metadata;
236 const QUAD_FAST *quad_fast;
237 const REAL *const*phi;
238
239 metadata = (PARAM_QUAD_METADATA *)qmd->param_data[bas_fcts->degree];
240 if (metadata == NULL) {
241 alloc_quad_metadata[bas_fcts->dim](quad, bas_fcts);
242 metadata = (PARAM_QUAD_METADATA *)qmd->param_data[bas_fcts->degree];
243 }
244
245 quad_fast = metadata->quad_fast;
246
247 (void)INIT_ELEMENT(el_info, quad_fast);
248
249 phi = quad_fast->phi;
250
251 for (iq = 0; iq < quad->n_points; iq++) {
252 SET_DOW(0.0, world[iq]);
253
254 for (i = 0; i < quad_fast->n_bas_fcts; i++) {
255 AXPY_DOW(phi[iq][i], local_coords[i], world[iq]);
256 }
257 }
258 } else {
259
260 (void)INIT_ELEMENT(el_info, quad);
261
262 dim = quad->dim;
263
264 for (iq = 0; iq < quad->n_points; iq++) {
265 SET_DOW(0.0, world[iq]);
266
267 for (i = 0; i < N_VERTICES(dim); i++) {
268 AXPY_DOW(quad->lambda[iq][i], local_coords[i], world[iq]);
269 }
270 }
271 }
272 } else {
273 REAL phi_lambda;
274
275 if (!data->i_am_affine) {
276 for (iq = 0; iq < N; iq++) {
277 SET_DOW(0.0, world[iq]);
278
279 for (i = 0; i < bas_fcts->n_bas_fcts; i++) {
280 phi_lambda = PHI(bas_fcts, i, lambda[iq]);
281
282 AXPY_DOW(phi_lambda, local_coords[i], world[iq]);
283 }
284 }
285 } else {
286 dim = el_info->mesh->dim;
287
288 for (iq = 0; iq < N; iq++) {
289 SET_DOW(0.0, world[iq]);
290
291 for (i = 0; i < N_VERTICES(dim); i++) {
292 AXPY_DOW(lambda[iq][i], local_coords[i], world[iq]);
293 }
294 }
295 }
296 }
297
298 return;
299 }
300
301 /* return: -2 if newton did not converge,
302 * -1 if inside,
303 * otherwise index of lambda <0
304 */
param_world_to_coord(const EL_INFO * el_info,int N,const REAL_D world[],REAL_B lambda[],int k[])305 static void param_world_to_coord(const EL_INFO *el_info, int N,
306 const REAL_D world[],
307 REAL_B lambda[], int k[])
308 {
309 FUNCNAME("param_world_to_coord");
310
311 /* stuff for newton */
312 static ISP_DATA isp_data;
313 static NLS_DATA nls_data = {
314 nl_update,
315 &isp_data,
316 nl_solve,
317 &isp_data,
318 NULL, /* norm */
319 NULL, /* norm-data */
320 NULL, /* work-space */
321 1e-15, /* tolerance */
322 1, /* restart */
323 1000, /* max_iter */
324 1, /* info */
325 -1.0, /* initial residual */
326 -1.0 /* residual */
327 };
328 REAL_B lambda0_const, lambda0;
329 int iter, max_iter;
330 REAL tolerance, lambda_tolerance;
331 int dim_of_mesh = MIN(DIM_MAX, el_info->mesh->dim);
332 int fail_count = 0;
333 int i, iq/* , dim */;
334 int k_dummy[N];
335
336 LAGRANGE_PARAM_DATA *data =
337 (LAGRANGE_PARAM_DATA *)el_info->mesh->parametric->data;
338 REAL_D *local_coords = data->local_coords;
339 const BAS_FCTS *bas_fcts = data->coords->fe_space->bas_fcts;
340
341 REAL_D coord_backup[N_VERTICES(dim_of_mesh)];
342
343 nls_data.max_iter = max_iter = data->max_iter;
344 nls_data.tolerance = tolerance = data->newton_tolerance;
345 lambda_tolerance = data->lambda_tolerance;
346
347 if (!k) k = k_dummy;
348
349 if (DIM_OF_WORLD != dim_of_mesh) {
350 ERROR_EXIT("DIM_OF_WORLD = %d != %d = dim_of_mesh.",
351 DIM_OF_WORLD, dim_of_mesh);
352 }
353 DEBUG_TEST_EXIT((el_info->fill_flag & FILL_COORDS) ||
354 el_info->mesh->parametric->use_reference_mesh,
355 "You must enable the use_reference_mesh entry in the "
356 "PARAMETRIC structure to use this function.\n");
357
358 if (!data->i_am_affine) {
359 for (iq = 0; iq < N; iq++) {
360 /* fill data */
361 COPY_DOW(world[iq], isp_data.x);
362 isp_data.el_info = el_info;
363 isp_data.bfcts = (BAS_FCTS*)bas_fcts;
364 isp_data.coords = local_coords;
365
366 /* For planar and not too heavily deformed meshes the standard
367 * world_to_coord should provide a reasonable initial value for
368 * points inside the simplex. If world lies outside the simplex
369 * we can neither guarantee the convergence of newton's method
370 * nor be sure that our initial lambda ist near the solution.
371 */
372 /* As global_refine does not fill the coords of projected nodes
373 * into el_info_coord when use_reference_mesh == true (and not
374 * at all if use_reference_mesh == false), we update el_info
375 * with the correct information.
376 */
377 for (i = 0; i < N_VERTICES(dim_of_mesh); i++) {
378 COPY_DOW(el_info->coord[i], coord_backup[i]);
379 COPY_DOW(local_coords[i], (REAL *)el_info->coord[i]);
380 }
381
382 world_to_coord(el_info, world[iq], lambda0_const);
383
384 /* Undo the hack above. */
385 for (i=0; i < N_VERTICES(dim_of_mesh); i++) {
386 COPY_DOW(coord_backup[i], (REAL *)el_info->coord[i]);
387 }
388
389 /* start newton, and retry if not successful */
390 for (fail_count = 0; fail_count < N_LAMBDA(dim_of_mesh)+1;) {
391 REAL sum = 0.0;
392
393 COPY_BAR(DIM_MAX, lambda0_const, lambda0);
394
395 if (fail_count > 0) {
396 #if 0
397 WARNING("Newton failed! current number of fails: %d.\n",
398 fail_count);
399 MSG("Retry with different lambda0...\n");
400 #endif
401 lambda0[fail_count-1] += 0.1;
402
403 for (i = 0; i < N_LAMBDA(dim_of_mesh); i++) sum += lambda0[i];
404 for (i = 0; i < N_LAMBDA(dim_of_mesh); i++) lambda0[i] /= sum;
405 }
406
407 iter = nls_newton(&nls_data, N_LAMBDA(dim_of_mesh), (REAL *)&lambda0);
408
409 if (iter > max_iter) {
410 fail_count++;
411 } else {
412 sum = 0;
413 for (i = 0; i < N_LAMBDA(dim_of_mesh); i++) {
414 lambda[iq][i] = lambda0[i];
415 sum += lambda0[i];
416 }
417 if (ABS(sum-1) > 10.0*REAL_EPSILON) {
418 #if 0
419 MSG("l_r[0]=%e, l_r[1]=%e, l_r[2]=%e\n",
420 lambda0[0], lambda0[1], lambda0[2]);
421 #endif
422 ERROR_EXIT("%e = sum{lambda[i]} != 1\n", ABS(sum-1));
423 }
424
425 #if 0
426 MSG("Success! Newton needed %d iterations and failed %d times.\n",
427 iter, fail_count);
428 #endif
429 break;
430 }
431 }
432 if (fail_count < N_LAMBDA(dim_of_mesh)+1) {
433 REAL lmin = 0.0;
434 k[iq] = -1;
435 for (i = 0; i < N_LAMBDA(dim_of_mesh); i++) {
436 if (lambda[iq][i] < lambda_tolerance) {
437 if (lambda[iq][i] < lmin) {
438 k[iq] = i;
439 lmin = lambda[iq][i];
440 }
441 }
442 }
443 } else {
444 WARNING("\a Newton failed totally. "
445 "Maybe you have to lift the tolerance a little, "
446 "or there is no solution.\n");
447 k[iq] = -2;
448 }
449 }
450 return;
451 } else {
452 for (iq = 0; iq < N; iq++) {
453 k[iq] = world_to_coord(el_info, world[iq], lambda[iq]);
454 }
455 return;
456 }
457 }
458
459 /* Return true if MESH carries our Lagrange parametric structure. */
_AI_is_lagrange_parametric(MESH * mesh)460 bool _AI_is_lagrange_parametric(MESH *mesh)
461 {
462 LAGRANGE_PARAM_DATA *data;
463
464 return (mesh->parametric != NULL &&
465 (data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data) != NULL &&
466 memcmp(data->magic, MAGIC_STRING, sizeof(data->magic)) == 0);
467 }
468
_AI_lagrange_strategy(MESH * mesh)469 PARAM_STRATEGY _AI_lagrange_strategy(MESH *mesh)
470 {
471 if (!_AI_is_lagrange_parametric(mesh)) {
472 return (PARAM_STRATEGY)-1;
473 }
474 return ((LAGRANGE_PARAM_DATA *)mesh->parametric->data)->strategy;
475 }
476
477 static void _AI_use_lagrange_parametric(MESH *mesh, int degree,
478 NODE_PROJECTION *n_proj,
479 PARAM_STRATEGY strategy,
480 FLAGS adm_flags,
481 MESH *master);
482
inherit_lagrange_parametric(MESH * slave)483 static void inherit_lagrange_parametric(MESH *slave)
484 {
485 FUNCNAME("inherit_lagrange_parametric");
486 MESH *master;
487 PARAMETRIC *m_parametric;
488 LAGRANGE_PARAM_DATA *m_data;
489
490 TEST_EXIT(slave, "No slave mesh given!\n");
491
492 master = ((MESH_MEM_INFO *)slave->mem_info)->master;
493 TEST_EXIT(master, "'%s' is not a slave mesh!\n", NAME(slave));
494
495 m_parametric = master->parametric;
496 TEST_EXIT(m_parametric, "'%s' is not a parametric mesh!\n", NAME(master));
497
498 m_data = (LAGRANGE_PARAM_DATA *)m_parametric->data;
499
500 /* Turn the slave mesh into a parametric mesh. This also transfers
501 * coordinates and "edge_projections".
502 */
503
504 _AI_use_lagrange_parametric(slave,
505 m_data->degree,
506 m_data->n_proj,
507 m_data->strategy,
508 m_data->coords->fe_space->admin->flags,
509 master);
510
511 return;
512 }
513
unchain_lagrange_parametric(MESH * slave)514 static void unchain_lagrange_parametric(MESH *slave)
515 {
516 PARAMETRIC *s_parametric = slave->parametric;
517 LAGRANGE_PARAM_DATA *s_data = (LAGRANGE_PARAM_DATA *)s_parametric->data;
518 int dim, degree, deg_idx;
519
520 dim = slave->dim;
521 degree = s_data->degree;
522
523 deg_idx = degree >= DEGREE_GEN ? DEGREE_GEN-1 : degree-1;
524 s_data->coords->refine_interpol = all_refine_interpol[dim][deg_idx];
525
526 return;
527 }
528
529 /*******************************************************************************
530 * use_lagrange_parametric(mesh, degree, n_proj, strategy):
531 * Transforms the mesh to use parametric simplices of degree "degree".
532 * All mesh elements which are affected by the node projection "n_proj"
533 * will be deformed. The parameter "strategy" determines how elements will
534 * become parametric. The following choices are available at the moment:
535 *
536 * strategy == PARAM_ALL:
537 * all elements of the mesh will be treated as parametric
538 * elements, implying that determinants and Jacobeans will
539 * be calculated at all quadrature points during assembly.
540 * This is useful e.g. for triangulations of embedded
541 * curved manifolds. Please note that during refinement a
542 * parent element will be split along the surface defined
543 * by the following equation
544 *
545 * lambda_0 = lambda_1
546 *
547 * defined in the nonlinear barycentric coords lambda_i.
548 * For some meshes, this will create better shaped
549 * simplices.
550 *
551 * strategy == PARAM_CURVED_CHILDS:
552 * only those elements of the mesh affected by "n_proj"
553 * will be treated as parametric elements. This enables
554 * some optimization during assembly. All children of the
555 * original curved simplices will be determined as above,
556 * implying that they too will be parametric.
557 *
558 * strategy == PARAM_STRAIGHT_CHILDS:
559 * only those elements of the mesh affected by "n_proj"
560 * will be treated as parametric elements. The children of
561 * parametric elements will only be parametric elements
562 * if they too are affected by "n_proj". Parent elements
563 * are split along straight lines or planes, which could
564 * lead to degenerate elements. The purpose of this
565 * algorithm is to keep the number of curved elements as
566 * small as possible. This is useful if one only wishes to
567 * use parametric elements to approximate a curved boundary
568 * or interface manifold.
569 *
570 * If n_proj == NULL, then all elements which have just some
571 * projection attached to it will become parametric. If n_proj !=
572 * NULL, then only those elements are changed into curved elements
573 * which have are subject to the same projection.
574 *
575 *****************************************************************************/
576
_AI_use_lagrange_parametric(MESH * mesh,int degree,NODE_PROJECTION * n_proj,PARAM_STRATEGY strategy,FLAGS adm_flags,MESH * master)577 static void _AI_use_lagrange_parametric(MESH *mesh, int degree,
578 NODE_PROJECTION *n_proj,
579 PARAM_STRATEGY strategy,
580 FLAGS adm_flags,
581 MESH *master)
582 {
583 FUNCNAME("use_lagrange_parametric");
584 LAGRANGE_PARAM_DATA *data;
585 PARAMETRIC *parametric;
586 DOF_REAL_D_VEC *coords;
587 DOF_PTR_VEC *edge_pr = NULL;
588 const BAS_FCTS *lagrange = NULL;
589 int dim, i, deg_idx;
590 MESH *slave;
591 bool selective = n_proj != NULL;
592 const FE_SPACE *fe_space;
593
594 TEST_EXIT(mesh,"No fe_space given!\n");
595
596 if (mesh->parametric != NULL) {
597 WARNING("There is already a parametric structure defined on this mesh!\n");
598 }
599
600 dim = mesh->dim;
601
602 TEST_EXIT((dim >= 0) && (dim <= DIM_MAX),
603 "Parametric elements of dimension %d "
604 "are not available for DIM_MAX == %d!\n", dim, DIM_MAX);
605
606 TEST_EXIT((degree > 0) && (degree <= DEGREE_MAX),
607 "Only implemented for 1 <= degree <= %d.\n", DEGREE_MAX);
608
609 deg_idx = degree >= DEGREE_GEN ? DEGREE_GEN-1 : degree-1;
610
611 TEST_EXIT(strategy >= PARAM_ALL && strategy <= PARAM_STRAIGHT_CHILDS,
612 "Only strategy 0, 1, 2 are implemented!\n");
613
614 /* Promote the strategy to PARAM_ALL if all elements have a
615 * projection attached; do this only for mesh-hierarchies.
616 */
617 if (master) {
618 MACRO_EL *mel;
619 int not_all;
620
621 not_all = false;
622 if (strategy != PARAM_ALL) {
623 for (mel = mesh->macro_els;
624 mel < mesh->macro_els + mesh->n_macro_el;
625 mel++) {
626 if (mel->projection[0] == NULL ||
627 (selective && n_proj != mel->projection[0])) {
628 not_all = true;
629 break;
630 }
631 }
632 if (!not_all) {
633 strategy = PARAM_ALL;
634 }
635 }
636 }
637
638 if (dim > 0 && degree > 1 && strategy != PARAM_ALL) {
639 int n_dof[N_NODE_TYPES] = { 0, };
640 const FE_SPACE *edge_fe_space;
641
642 if (dim > 1) {
643 n_dof[EDGE] = 1;
644 } else {
645 n_dof[CENTER] = 1;
646 }
647
648 /* "edge_projections" belongs to a periodic DOF_ADMIN on periodic
649 * meshes.
650 */
651 edge_fe_space = get_dof_space(mesh, "Edge dof fe_space", n_dof,
652 ADM_FLAGS_DFLT|ADM_PERIODIC);
653
654 edge_pr = get_dof_ptr_vec("Edge projections", edge_fe_space);
655 FOR_ALL_DOFS(edge_fe_space->admin, edge_pr->vec[dof] = NULL);
656 free_fe_space(edge_fe_space);
657 }
658
659 lagrange = get_lagrange(dim, degree);
660 fe_space = get_fe_space(mesh,
661 lagrange->name, lagrange, DIM_OF_WORLD, adm_flags);
662 coords = get_dof_real_d_vec("Lagrange parametric coordinates", fe_space);
663 coords->refine_interpol = all_refine_interpol[dim][deg_idx];
664 coords->coarse_restrict = all_coarse_interpol[dim][deg_idx];
665
666 data = MEM_CALLOC(1, LAGRANGE_PARAM_DATA);
667 data->degree = degree;
668 data->strategy = strategy;
669 data->n_proj = n_proj;
670 data->coords = coords;
671 data->edge_projections = edge_pr;
672 data->n_local_coords = fe_space->bas_fcts->n_bas_fcts;
673
674 data->max_iter = 1000;
675 data->newton_tolerance = 5e-14;
676 data->lambda_tolerance = -75*REAL_EPSILON;
677
678 GET_PARAMETER(0, "parametric->lagrange->param_world_to_coord"
679 "->newton->max_iter" , "%d", &data->max_iter);
680 GET_PARAMETER(0, "parametric->lagrange->param_world_to_coord"
681 "->newton->tolerance", "%f", &data->newton_tolerance);
682 GET_PARAMETER(0, "parametric->lagrange->param_world_to_coord"
683 "->lambda->tolerance", "%f", &data->lambda_tolerance);
684
685 if (dim > 0 && degree > 1) {
686 data->param_local_coords = MEM_CALLOC(lagrange->n_bas_fcts, REAL_D);
687 data->i_am_affine = false;
688 } else {
689 data->param_local_coords = NULL;
690 data->i_am_affine = true;
691 }
692 data->local_coords = strategy == PARAM_ALL ? data->param_local_coords : NULL;
693 memcpy(data->magic, MAGIC_STRING, sizeof(data->magic));
694
695 /* Inherit the data from the master to assure consistent coordinate
696 * information for the entire mesh hierarchy.
697 */
698
699 if (master) {
700 #if DIM_MAX > 0
701 if (dim == 0) {
702 slave_fill_coords_0d(data);
703 coords->refine_interpol = NULL;
704 # if DIM_MAX > 1
705 } else if (dim == 1) {
706 slave_fill_coordsY_1d(data);
707 coords->refine_interpol = slave_refine_interpolY_1d;
708 # if DIM_MAX > 2
709 } else if (dim == 2) {
710 slave_fill_coordsY_2d(data);
711 coords->refine_interpol = slave_refine_interpolY_2d;
712 # endif
713 # endif
714 } else {
715 ERROR_EXIT("Strange combinations of dimensions: %d / %d\n",
716 dim, master->dim);
717 }
718 #endif
719 } else {
720 all_fill_coords[dim][deg_idx](data);
721 }
722
723 /* update mesh->bbox and mesh->diam */
724 SET_DOW(REAL_MAX, mesh->bbox[0]);
725 SET_DOW(REAL_MIN, mesh->bbox[1]);
726 FOR_ALL_DOFS(coords->fe_space->admin,
727 REAL *value = coords->vec[dof];
728 for (i = 0; i < DIM_OF_WORLD; i++) {
729 mesh->bbox[0][i] = MIN(mesh->bbox[0][i], value[i]);
730 mesh->bbox[1][i] = MAX(mesh->bbox[1][i], value[i]);
731 });
732 AXPBY_DOW(1.0, mesh->bbox[1], -1.0, mesh->bbox[0], mesh->diam);
733
734 parametric = MEM_CALLOC(1, PARAMETRIC);
735 *parametric = *all_lagrange[dim][deg_idx];
736 parametric->data = data;
737 mesh->parametric = parametric;
738 parametric->not_all = degree == 1 || data->strategy != PARAM_ALL;
739
740 /* Set inheritance function. */
741 parametric->inherit_parametric = inherit_lagrange_parametric;
742 parametric->unchain_parametric = unchain_lagrange_parametric;
743
744 /* Now handle submeshes. Pass ourselves as "master" to the sub-meshes */
745
746 if(mesh->dim > 0) {
747 for (i = 0; i < ((MESH_MEM_INFO *)mesh->mem_info)->n_slaves; i++) {
748 slave = ((MESH_MEM_INFO *)mesh->mem_info)->slaves[i];
749
750 _AI_use_lagrange_parametric(slave, degree, n_proj,
751 strategy, adm_flags, mesh);
752 }
753 }
754 }
755
use_lagrange_parametric(MESH * mesh,int degree,NODE_PROJECTION * n_proj,FLAGS flags)756 void use_lagrange_parametric(MESH *mesh, int degree,
757 NODE_PROJECTION *n_proj,
758 FLAGS flags)
759 {
760 FUNCNAME("use_lagrange_parametric");
761 PARAM_STRATEGY strategy;
762 FLAGS adm_flags;
763
764 if (_AI_is_lagrange_parametric(mesh)) {
765 WARNING("The mesh already has a parametric structure! "
766 "A second call will likely corrupt your mesh. "
767 "Returning without change.\n");
768 return;
769 }
770
771 TEST_EXIT(get_master(mesh) == NULL,
772 "ERROR: Parametric structures must be added on the "
773 "top-most master mesh of a sub-mesh hierarchy.\n");
774
775 strategy = (PARAM_STRATEGY)(flags & ~PARAM_PERIODIC_COORDS);
776 adm_flags = (flags & PARAM_PERIODIC_COORDS) ? ADM_PERIODIC : ADM_FLAGS_DFLT;
777 _AI_use_lagrange_parametric(mesh, degree, n_proj, strategy, adm_flags, NULL);
778 }
779
780 /*******************************************************************************
781 * get_lagrange_coords(mesh): Return the DOF_REAL_D_VEC coordinate vector
782 * used to defined the Lagrange type parametric elements implemented above.
783 * USE WITH CARE, ESPECIALLY IF YOU ARE CHANGING VALUES! Problems may arise
784 * because elements are not recognized as curved elements unless the
785 * internally used "edge_projections" vector is also set.
786 ******************************************************************************/
get_lagrange_coords(MESH * mesh)787 DOF_REAL_D_VEC *get_lagrange_coords(MESH *mesh)
788 {
789 FUNCNAME("get_lagrange_coords");
790
791 TEST_EXIT(mesh, "No mesh given!\n");
792
793 if (!_AI_is_lagrange_parametric(mesh)) {
794 /* Return NULL rather than bailing out. Let the caller take the
795 * proper actions.
796 */
797 return NULL;
798 } else {
799 return ((LAGRANGE_PARAM_DATA *)mesh->parametric->data)->coords;
800 }
801 }
802
803
get_lagrange_edge_projections(MESH * mesh)804 DOF_PTR_VEC *get_lagrange_edge_projections(MESH *mesh)
805 {
806 FUNCNAME("get_lagrange_edge_projections");
807
808 TEST_EXIT(mesh, "No mesh given!\n");
809
810 if (!_AI_is_lagrange_parametric(mesh)) {
811 /* Return NULL rather than bailing out. Let the caller take the
812 * proper actions.
813 */
814 return NULL;
815 } else {
816 return ((LAGRANGE_PARAM_DATA *)mesh->parametric->data)->edge_projections;
817 }
818 }
819
820 /** A more secure interface to the coordinate-vector: only the @b
821 * values are copied, the application does not get access to the
822 * DOF_REAL_D_VEC itself. This function makes sure that affine
823 * elements remain affine. We also handle the case when we do not have
824 * a real parameterisation but just affine linear elements with
825 * quasi-parametric coordinates stored in the EL->new_coord. The
826 * degree of the "coords" vector has to match the degree of the
827 * iso-parametric parameterisation. We also recompute mesh->diam here.
828 *
829 * @param[in] mesh The mesh to change the coordinates of.
830 *
831 * @param[in,out] coords The coordinate vector to install or where to
832 * store the mesh's current coordinates in.
833 *
834 * @param[in] dir Direction of the copy action: dir == 0 means to
835 * store the mesh's current coordinates in coords. dir == 1 means
836 * to install the coordinates specified by coords in the mesh, either
837 * using EL->new_coord or by copying the values into the internal
838 * coordinate vector.
839 */
copy_lagrange_coords(MESH * mesh,DOF_REAL_D_VEC * coords,bool tomesh)840 void copy_lagrange_coords(MESH *mesh, DOF_REAL_D_VEC *coords, bool tomesh)
841 {
842 PARAMETRIC *parametric = mesh->parametric;
843 int i, j;
844 int dim = mesh->dim;
845
846 if (tomesh) {
847 SET_DOW(REAL_MAX, mesh->bbox[0]);
848 SET_DOW(REAL_MIN, mesh->bbox[1]);
849 FOR_ALL_DOFS(coords->fe_space->admin,
850 REAL *value = coords->vec[dof];
851 for (i = 0; i < DIM_OF_WORLD; i++) {
852 mesh->bbox[0][i] = MIN(mesh->bbox[0][i], value[i]);
853 mesh->bbox[1][i] = MAX(mesh->bbox[1][i], value[i]);
854 });
855 AXPBY_DOW(1.0, mesh->bbox[1], -1.0, mesh->bbox[0], mesh->diam);
856 }
857
858 if (parametric) {
859 LAGRANGE_PARAM_DATA *data;
860
861 TEST_EXIT(_AI_is_lagrange_parametric(mesh),
862 "Parametric data has not type LAGRANGE_PARAM_DATA.\n");
863
864 data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data;
865
866 TEST_EXIT(coords->fe_space->bas_fcts == data->coords->fe_space->bas_fcts,
867 "basis function mismatch.\n");
868
869 data = (LAGRANGE_PARAM_DATA *)mesh->parametric->data;
870
871 if (data->strategy == PARAM_ALL || tomesh == false) { /* the easy case */
872 if (tomesh) {
873 dof_copy_d(coords, data->coords);
874 } else {
875 dof_copy_d(data->coords, coords);
876 }
877 } else {
878 const BAS_FCTS *bas_fcts = coords->fe_space->bas_fcts;
879 const DOF_ADMIN *admin = coords->fe_space->admin;
880 const REAL_B *nodes = LAGRANGE_NODES(bas_fcts);
881 NODE_PROJECTION **edge_pr = (NODE_PROJECTION **)data->edge_projections->vec;
882 int n0_edge_pr, node_e;
883 DOF dof[bas_fcts->n_bas_fcts];
884
885 dof_copy_d(coords, data->coords);
886
887 node_e = mesh->node[EDGE];
888 n0_edge_pr = data->edge_projections->fe_space->admin->n0_dof[EDGE];
889
890 /* We need a mesh-traversal to make sure that affine elements
891 * remain affine.
892 *
893 * inefficient like hell, but so what
894 */
895
896 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) {
897 int is_affine = true;
898
899 for (i = 0; i < N_EDGES(dim); i++) {
900 if (edge_pr[el_info->el->dof[node_e+i][n0_edge_pr]] != NULL) {
901 is_affine = false;
902 break;
903 }
904 }
905
906 if (is_affine) {
907 GET_DOF_INDICES(bas_fcts, el_info->el, admin, dof);
908 for (i = N_VERTICES(dim); i < bas_fcts->n_bas_fcts; i++) {
909 AXEY_DOW(nodes[i][0], data->coords->vec[dof[0]],
910 data->coords->vec[dof[i]]);
911 for (j = 1; j < N_VERTICES(dim); j++) {
912 AXPY_DOW(nodes[i][j], data->coords->vec[dof[j]],
913 data->coords->vec[dof[i]]);
914 }
915 }
916 }
917 } TRAVERSE_NEXT();
918 }
919 } else {
920 /* quasi-parametric case, "parametric" coords are stored in el->new_coord
921 *
922 * We do not use FILL_COORDS but fill the coordinates ourselves to
923 * avoid too many copy instructions.
924 */
925
926 if (tomesh) { /* coords -> mesh */
927 const DOF_ADMIN *admin = coords->fe_space->admin;
928 int node_v, n0_v, v;
929
930 node_v = mesh->node[VERTEX];
931 n0_v = admin->n0_dof[VERTEX];
932
933 TRAVERSE_FIRST(mesh, -1, CALL_EVERY_EL_PREORDER|FILL_NEIGH) {
934 if (el_info->level == 0) {
935 for (v = 0; v < N_VERTICES(dim); v++) {
936 COPY_DOW(coords->vec[el_info->el->dof[node_v+v][n0_v]],
937 *el_info->macro_el->coord[v]);
938 }
939 }
940 if (el_info->el->child[0]) {
941 DOF cdof = el_info->el->child[0]->dof[node_v+dim][n0_v];
942
943 TEST_EXIT(el_info->el->new_coord != NULL,
944 "el_info->el->new_coord == NULL");
945
946 COPY_DOW(coords->vec[cdof], el_info->el->new_coord);
947
948 }
949 } TRAVERSE_NEXT();
950
951 } else { /* mesh -> coords */
952 const DOF_ADMIN *admin = coords->fe_space->admin;
953 int node_v, n0_v, v, dim = mesh->dim;
954
955 node_v = mesh->node[VERTEX];
956 n0_v = admin->n0_dof[VERTEX];
957
958 TRAVERSE_FIRST(mesh, -1, CALL_EVERY_EL_PREORDER) {
959 if (el_info->level == 0) {
960 for (v = 0; v < N_VERTICES(dim); v++) {
961 COPY_DOW(*el_info->macro_el->coord[v],
962 coords->vec[el_info->el->dof[node_v+v][n0_v]]);
963 }
964 }
965 if (el_info->el->child[0]) {
966 DOF cdof = el_info->el->child[0]->dof[node_v+dim][n0_v];
967 if (el_info->el->new_coord) {
968 COPY_DOW(el_info->el->new_coord, coords->vec[cdof]);
969 } else {
970 DOF v0 = el_info->el->dof[node_v+0][n0_v];
971 DOF v1 = el_info->el->dof[node_v+1][n0_v];
972
973 AXPBY_DOW(0.5, coords->vec[v0], 0.5, coords->vec[v1],
974 coords->vec[cdof]);
975 }
976 }
977 } TRAVERSE_NEXT();
978 }
979 }
980 }
981
982 /* update(update data, dim, current iterate , update matrix true/false, rhs) */
nl_update(void * vdata,int dim,const REAL * lambda_n,bool upd_DF,REAL * f)983 static void nl_update(void *vdata, int dim, const REAL *lambda_n,
984 bool upd_DF, REAL *f)
985 {
986 int i, j, k;
987 ISP_DATA *data = (ISP_DATA *)vdata;
988 REAL_D *coords = data->coords;
989 BAS_FCTS *bfcts = data->bfcts;
990 REAL_D x_n;
991 const REAL *grd_phi;
992
993 /* Fs: IR^(n+1) -> IR^(n), Fs(lambda) = sum_i{a_i \phi_i(lambda)},
994 * a_i: world coords of the edges, lambda: barycentric coords,
995 * phi_i: i-th lagrange basis fct.
996 * f: IR^(n+1) -> IR^(n+1)
997 * f(lambda):= / Fs(lambda) - x \
998 * \ sum(lambda[i]) - 1 /
999 */
1000
1001 if (f) {
1002 param_coord_to_world(data->el_info, NULL, 1,
1003 (const REAL_B *)lambda_n, &x_n);
1004 AXPBY_DOW(1.0, x_n, -1.0, data->x, f);
1005 f[DIM_OF_WORLD] = -1;
1006 for (i=0; i<DIM_OF_WORLD+1; i++) {
1007 f[DIM_OF_WORLD] += lambda_n[i];
1008 }
1009 }
1010
1011 /* Df(lambda_n) */
1012 if (upd_DF) {
1013 for (j=0; j < DIM_OF_WORLD+1; j++) {
1014 for (i=0; i < DIM_OF_WORLD; i++) {
1015 data->Df[i][j] = 0;
1016 for (k=0; k < data->bfcts->n_bas_fcts; k++) {
1017 grd_phi = GRD_PHI(bfcts, k, lambda_n);
1018 data->Df[i][j] += coords[k][i] * grd_phi[j];
1019 }
1020 }
1021 data->Df[DIM_OF_WORLD][j] = 1;
1022 }
1023 }
1024 return;
1025 }
1026
1027 /* solve(solve data, dim, F, d) */
nl_solve(void * vdata,int dim,const REAL * rhs,REAL * d)1028 static int nl_solve(void *vdata, int dim, const REAL *rhs, REAL *d)
1029 {
1030 ISP_DATA *data = (ISP_DATA *)vdata;
1031 REAL Df[DIM_OF_WORLD+1][DIM_OF_WORLD+1];
1032 REAL rhs_cp[DIM_OF_WORLD+1];
1033
1034 memcpy(Df, data->Df, sizeof(Df));
1035
1036 memcpy(rhs_cp, rhs, sizeof(rhs_cp));
1037
1038 square_gauss((REAL *)Df, rhs_cp, d, DIM_OF_WORLD+1, 1);
1039
1040 return 0;
1041 }
1042