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: adapt.c */
7 /* */
8 /* description: adaptive procedures for stationary and instationary */
9 /* problems */
10 /* */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Alfred Schmidt */
14 /* Zentrum fuer Technomathematik */
15 /* Fachbereich 3 Mathematik/Informatik */
16 /* Universitaet Bremen */
17 /* Bibliothekstr. 2 */
18 /* D-28359 Bremen, Germany */
19 /* */
20 /* Kunibert G. Siebert */
21 /* Institut fuer Mathematik */
22 /* Universitaet Augsburg */
23 /* Universitaetsstr. 14 */
24 /* D-86159 Augsburg, Germany */
25 /* */
26 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
27 /* */
28 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
29 /* */
30 /*--------------------------------------------------------------------------*/
31
32 #ifdef HAVE_CONFIG_H
33 # include "config.h"
34 #endif
35
36 #include "alberta.h"
37 #include <time.h>
38
39 #define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC)
40
41 /*--------------------------------------------------------------------------*/
42 /* adapt_mesh: call marking(), refine(), and coarsen() */
43 /*--------------------------------------------------------------------------*/
44
adapt_mesh(MESH * mesh,ADAPT_STAT * adapt)45 static U_CHAR adapt_mesh(MESH *mesh, ADAPT_STAT *adapt)
46 {
47 FUNCNAME("adapt_mesh");
48 U_CHAR flag = 0;
49 U_CHAR mark_flag;
50 int n_elements, iadmin;
51 clock_t first = clock();
52
53 TEST_EXIT(adapt, "no ADAPT_STAT\n");
54
55 if (adapt->marking)
56 mark_flag = adapt->marking(mesh, adapt);
57 else
58 mark_flag = marking(mesh, adapt);
59
60 if ((!adapt->coarsen_allowed))
61 mark_flag &= MESH_REFINED; /* use refine mark only */
62
63 if (adapt->build_before_refine)
64 adapt->build_before_refine(mesh, mark_flag);
65
66 n_elements = mesh->n_elements;
67
68 if (mark_flag & MESH_REFINED)
69 flag = refine(mesh, adapt->adaptation_fill_flags);
70
71 if (flag & MESH_REFINED)
72 {
73 n_elements = mesh->n_elements - n_elements;
74 INFO(adapt->info,8,
75 "%d element%s refined, giving %d element%s\n",
76 n_elements, n_elements > 1 ? "s" : "",
77 mesh->n_elements, mesh->n_elements > 1 ? "s" : "");
78 for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++)
79 INFO(adapt->info,7,"%d DOFs of admin <%s>\n",
80 mesh->dof_admin[iadmin]->used_count,
81 NAME(mesh->dof_admin[iadmin]));
82 }
83 else
84 INFO(adapt->info,8,"no element refined\n");
85
86
87 if (adapt->build_before_coarsen)
88 adapt->build_before_coarsen(mesh, mark_flag);
89
90 n_elements = mesh->n_elements;
91
92 if (mark_flag & MESH_COARSENED)
93 flag |= coarsen(mesh, adapt->adaptation_fill_flags);
94
95 if (flag & MESH_COARSENED)
96 {
97 n_elements -= mesh->n_elements;
98 INFO(adapt->info,8,
99 "%d element%s coarsened, giving %d element%s\n",
100 n_elements, n_elements > 1 ? "s" : "",
101 mesh->n_elements, mesh->n_elements > 1 ? "s" : "");
102 for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++)
103 INFO(adapt->info,7,"%d DOFs of dof_admin <%s>\n",
104 mesh->dof_admin[iadmin]->used_count,
105 NAME(mesh->dof_admin[iadmin]));
106 }
107 else
108 INFO(adapt->info,8,"no element coarsened\n");
109
110
111 if (adapt->build_after_coarsen)
112 adapt->build_after_coarsen(mesh, flag);
113
114 INFO(adapt->info,6,"adapting mesh and build needed %.5lg seconds\n",
115 TIME_USED(first,clock()));
116
117 return(flag);
118 }
119
120 /*--------------------------------------------------------------------------*/
121 /* adapt_method_stat: adapt for stationary problems */
122 /*--------------------------------------------------------------------------*/
123
124
adapt_method_stat(MESH * mesh,ADAPT_STAT * adapt)125 void adapt_method_stat(MESH *mesh, ADAPT_STAT *adapt)
126 {
127 FUNCNAME("adapt_method_stat");
128 int iter;
129 REAL est;
130 clock_t first;
131
132 TEST_EXIT(mesh, "no MESH\n");
133 TEST_EXIT(adapt, "no ADAPT_STAT\n");
134
135 /* get solution on initial mesh */
136 if (adapt->build_before_refine) adapt->build_before_refine(mesh, 0);
137 if (adapt->build_before_coarsen) adapt->build_before_coarsen(mesh, 0);
138 if (adapt->build_after_coarsen) adapt->build_after_coarsen(mesh, 0);
139 if (adapt->solve)
140 {
141 first = clock();
142 adapt->solve(mesh);
143 INFO(adapt->info,8,"solution of discrete system needed %.5lg seconds\n",
144 TIME_USED(first,clock()));
145 }
146
147 first = clock();
148 est = adapt->estimate ? adapt->estimate(mesh, adapt) : 0.0;
149 INFO(adapt->info,8,"estimation of the error needed %.5lg seconds\n",
150 TIME_USED(first,clock()));
151
152 for (iter = 0;
153 (est > adapt->tolerance) &&
154 ((adapt->max_iteration <= 0) || (iter < adapt->max_iteration));
155 iter++)
156 {
157 if (adapt_mesh(mesh, adapt))
158 {
159 if (adapt->solve)
160 {
161 first = clock();
162 adapt->solve(mesh);
163 INFO(adapt->info,8,
164 "solution of discrete system needed %.5lg seconds\n",
165 TIME_USED(first,clock()));
166 }
167 first = clock();
168 est = adapt->estimate ? adapt->estimate(mesh, adapt) : 0.0;
169 INFO(adapt->info,8,"estimation of the error needed %.5lg seconds\n",
170 TIME_USED(first,clock()));
171 }
172 else
173 {
174 ERROR("no mesh adaption, but estimate above tolerance ???\n");
175 break;
176 }
177 INFO(adapt->info, 4,"iter: %d", iter);
178 PRINT_INFO(adapt->info, 4,", tol = %.4le", adapt->tolerance);
179 PRINT_INFO(adapt->info, 4,", estimate = %.4le\n", est);
180 }
181
182 if (est > adapt->tolerance)
183 {
184 MSG("max_iterations REACHED: %d\n", adapt->max_iteration);
185 MSG("prescribed tolerance %le\n", adapt->tolerance);
186 MSG("finished with estimate %le\n", est);
187 }
188 else
189 {
190 INFO(adapt->info, 2,"no of iterations: %d\n", iter);
191 INFO(adapt->info, 2,"prescribed tolerance %.4le\n", adapt->tolerance);
192 INFO(adapt->info, 2,"finished with estimate %.4le\n", est);
193 }
194 return;
195 }
196
197 /*--------------------------------------------------------------------------*/
198 /* adapt_method_instat(): adapt for timedependent problems */
199 /*--------------------------------------------------------------------------*/
200
explicit_time_strategy(MESH * mesh,ADAPT_INSTAT * adapt)201 static void explicit_time_strategy(MESH *mesh, ADAPT_INSTAT *adapt)
202 {
203 FUNCNAME("explicit_time_strategy");
204 ADAPT_STAT *adapt_s = adapt->adapt_space;
205
206 /* first estimate must be computed before adapt_mesh() is called */
207 if (adapt->time <= adapt->start_time) {
208 if (adapt_s->estimate)
209 adapt_s->estimate(mesh, adapt_s);
210 }
211
212 adapt->time += adapt->timestep;
213 if (adapt->set_time)
214 adapt->set_time(mesh, adapt);
215
216 INFO(adapt->info,6,"time = %.4le, timestep = %.4le\n", adapt->time, adapt->timestep);
217
218 adapt_mesh(mesh, adapt_s);
219
220 if (adapt_s->solve)
221 adapt_s->solve(mesh);
222
223 if (adapt_s->estimate)
224 adapt_s->estimate(mesh, adapt_s); /* return value is not used */
225
226 return;
227 }
228
229
implicit_time_strategy(MESH * mesh,ADAPT_INSTAT * adapt)230 static void implicit_time_strategy(MESH *mesh, ADAPT_INSTAT *adapt)
231 {
232 FUNCNAME("implicit_time_strategy");
233 int iter, iter_s;
234 ADAPT_STAT *adapt_s = adapt->adapt_space;
235 REAL err_space, err_time, space_err_limit, time_err_limit, time_err_low;
236
237 space_err_limit = adapt->tolerance * adapt->rel_space_error;
238 time_err_limit = adapt->tolerance * adapt->rel_time_error;
239 time_err_low = adapt->tolerance * adapt->rel_time_error
240 * adapt->time_theta_2;
241 iter = iter_s = 0;
242 err_time = 0.0;
243
244 do
245 {
246 adapt->time += adapt->timestep;
247 if (adapt->set_time)
248 adapt->set_time(mesh, adapt);
249
250 INFO(adapt->info,6,"time = %.4le, try timestep = %.4le\n",
251 adapt->time, adapt->timestep);
252
253 if (adapt_s->build_before_refine) adapt_s->build_before_refine(mesh, 0);
254 if (adapt_s->build_before_coarsen) adapt_s->build_before_coarsen(mesh, 0);
255 if (adapt_s->build_after_coarsen) adapt_s->build_after_coarsen(mesh, 0);
256 if (adapt_s->solve) adapt_s->solve(mesh);
257 err_space = adapt_s->estimate ? adapt_s->estimate(mesh, adapt_s) : 0.0;
258 if (adapt->get_time_est)
259 err_time = adapt->get_time_est(mesh, adapt);
260
261 iter++;
262 if (iter > adapt->max_iteration) break;
263
264 if (err_time > time_err_limit)
265 {
266 adapt->time -= adapt->timestep;
267 adapt->timestep *= adapt->time_delta_1;
268 continue;
269 }
270
271 do
272 {
273 if (adapt_mesh(mesh, adapt_s))
274 {
275 adapt_s->solve(mesh);
276 err_space = adapt_s->estimate ? adapt_s->estimate(mesh, adapt_s) : 0.0;
277 if (adapt->get_time_est) {
278 err_time = adapt->get_time_est(mesh, adapt);
279 if (err_time > time_err_limit)
280 {
281 adapt->time -= adapt->timestep;
282 adapt->timestep *= adapt->time_delta_1;
283 break;
284 }
285 }
286 }
287
288 iter_s++;
289 if (iter_s > adapt_s->max_iteration) break;
290
291 } while(err_space > space_err_limit);
292
293 } while(err_time > time_err_limit);
294
295 if (adapt->get_time_est)
296 if (err_time <= time_err_low)
297 {
298 adapt->timestep *= adapt->time_delta_2;
299 }
300
301 return;
302 }
303
304
one_timestep(MESH * mesh,ADAPT_INSTAT * adapt)305 static void one_timestep(MESH *mesh, ADAPT_INSTAT *adapt)
306 {
307 FUNCNAME("one_timestep");
308
309 switch(adapt->strategy)
310 {
311 case 0:
312 explicit_time_strategy(mesh, adapt);
313 break;
314
315 case 1:
316 implicit_time_strategy(mesh, adapt);
317 break;
318
319 default:
320 MSG("unknown adapt->strategy = %d; use explicit strategy\n");
321 explicit_time_strategy(mesh, adapt);
322 }
323
324 return;
325 }
326
327 /*--------------------------------------------------------------------------*/
328
adapt_method_instat(MESH * mesh,ADAPT_INSTAT * adapt)329 void adapt_method_instat(MESH *mesh, ADAPT_INSTAT *adapt)
330 {
331 FUNCNAME("adapt_method_instat");
332
333 TEST_EXIT(adapt, "no ADAPT_INSTAT\n");
334
335 /*--------------------------------------------------------------------------*/
336 /* adaptation of the initial grid: done by adapt_method_stat */
337 /*--------------------------------------------------------------------------*/
338
339 adapt->time = adapt->start_time;
340 if (adapt->set_time) adapt->set_time(mesh, adapt);
341
342 adapt->adapt_initial->tolerance
343 = adapt->tolerance * adapt->rel_initial_error;
344 adapt->adapt_space->tolerance
345 = adapt->tolerance * adapt->rel_space_error;
346
347 adapt_method_stat(mesh, adapt->adapt_initial);
348 if (adapt->close_timestep)
349 adapt->close_timestep(mesh, adapt);
350
351 /*--------------------------------------------------------------------------*/
352 /* adaptation of timestepsize and mesh: done by one_timestep */
353 /*--------------------------------------------------------------------------*/
354
355 while (adapt->time < adapt->end_time)
356 {
357 if (adapt->init_timestep)
358 adapt->init_timestep(mesh, adapt);
359
360 if (adapt->one_timestep)
361 adapt->one_timestep(mesh, adapt);
362 else
363 one_timestep(mesh, adapt);
364
365 if (adapt->close_timestep)
366 adapt->close_timestep(mesh, adapt);
367 }
368 }
369
370 /*--------------------------------------------------------------------------*/
371 /*--------------------------------------------------------------------------*/
372 /*--------------------------------------------------------------------------*/
373 /* marking functions!!! */
374 /*--------------------------------------------------------------------------*/
375 /*--------------------------------------------------------------------------*/
376 /*--------------------------------------------------------------------------*/
377
378 typedef struct adapt_traverse_data {
379 REAL (*get_el_est)(EL *el);
380 REAL (*get_el_estc)(EL *el);
381 int el_mark, el_mark_c;
382 S_CHAR g_mark_refine;
383 S_CHAR g_mark_coarse;
384 REAL err_max, err_sum;
385 int mark_flag;
386 REAL mark_r_limit, mark_c_limit;
387 REAL GERS_sum;
388 } ADAPT_TRAVERSE_DATA;
389
390 /*--------------------------------------------------------------------------*/
391 /*--------------------------------------------------------------------------*/
392 /**** marking functions for several adaptive strategies ****/
393 /*--------------------------------------------------------------------------*/
394 /*--------------------------------------------------------------------------*/
395
marking_fct_1(const EL_INFO * el_info,void * data)396 static void marking_fct_1(const EL_INFO *el_info, void *data)
397 {
398 el_info->el->mark = ((ADAPT_TRAVERSE_DATA *)data)->g_mark_refine;
399 ((ADAPT_TRAVERSE_DATA *)data)->mark_flag = 1;
400 ((ADAPT_TRAVERSE_DATA *)data)->el_mark++;
401
402 return;
403 }
404
405
marking_fct_2(const EL_INFO * el_info,void * data)406 static void marking_fct_2(const EL_INFO *el_info, void *data)
407 {
408 ADAPT_TRAVERSE_DATA *ud = (ADAPT_TRAVERSE_DATA *)data;
409 REAL error = ud->get_el_est(el_info->el);
410
411 if (error > ud->mark_r_limit)
412 {
413 el_info->el->mark = ud->g_mark_refine;
414 ud->mark_flag = 1;
415 ud->el_mark++;
416 }
417 else if (error <= ud->mark_c_limit)
418 {
419 if (!ud->get_el_estc ||
420 (error + ud->get_el_estc(el_info->el) <= ud->mark_c_limit))
421 {
422 el_info->el->mark = ud->g_mark_coarse;
423 ud->mark_flag = 1;
424 ud->el_mark_c++;
425 }
426 }
427 return;
428 }
429
430
marking_fct_3(const EL_INFO * el_info,void * data)431 static void marking_fct_3(const EL_INFO *el_info, void *data)
432 {
433 ADAPT_TRAVERSE_DATA *ud = (ADAPT_TRAVERSE_DATA *)data;
434 REAL error = ud->get_el_est(el_info->el);
435
436 if (error > ud->mark_r_limit)
437 {
438 el_info->el->mark = ud->g_mark_refine;
439 ud->mark_flag = 1;
440 ud->el_mark++;
441 }
442 else if (error <= ud->mark_c_limit)
443 {
444 if (!ud->get_el_estc ||
445 (error + ud->get_el_estc(el_info->el) <= ud->mark_c_limit))
446 {
447 el_info->el->mark = ud->g_mark_coarse;
448 ud->mark_flag = 1;
449 ud->el_mark_c++;
450 }
451 }
452 return;
453 }
454
marking_fct_4(const EL_INFO * el_info,void * data)455 static void marking_fct_4(const EL_INFO *el_info, void *data)
456 {
457 ADAPT_TRAVERSE_DATA *ud = (ADAPT_TRAVERSE_DATA *)data;
458 REAL error = ud->get_el_est(el_info->el);
459
460 if (error > ud->mark_r_limit)
461 {
462 ud->GERS_sum += error;
463 el_info->el->mark = ud->g_mark_refine;
464 ud->mark_flag = 1;
465 ud->el_mark++;
466 }
467 return;
468 }
469
marking_fct_4c(const EL_INFO * el_info,void * data)470 static void marking_fct_4c(const EL_INFO *el_info, void *data)
471 {
472 ADAPT_TRAVERSE_DATA *ud = (ADAPT_TRAVERSE_DATA *)data;
473 REAL error;
474
475 if (el_info->el->mark <= 0)
476 {
477 error = ud->get_el_est(el_info->el);
478 if (ud->get_el_estc) error += ud->get_el_estc(el_info->el);
479
480 if (error <= ud->mark_c_limit)
481 {
482 ud->GERS_sum += error;
483 el_info->el->mark = ud->g_mark_coarse;
484 ud->mark_flag = 1;
485 ud->el_mark_c++;
486 }
487 else
488 el_info->el->mark = 0;
489 }
490
491 return;
492 }
493
494 /*--------------------------------------------------------------------------*/
495
marking(MESH * mesh,ADAPT_STAT * adapt)496 int marking(MESH *mesh, ADAPT_STAT *adapt)
497 {
498 FUNCNAME("marking");
499 static REAL old_err_sum = 0.0;
500 static ADAPT_TRAVERSE_DATA td[1];
501 REAL improv, Ltheta, redfac, wanted;
502 REAL eps_p, ES_theta_p,ES_theta_c_p;
503 REAL MS_gamma_p, MS_gamma_c_p, GERS_gamma;
504
505 TEST_EXIT(adapt, "no adapt_stat\n");
506 if(adapt->strategy > 1)
507 TEST_EXIT(td->get_el_est = adapt->get_el_est, "no adapt->get_el_est\n");
508 TEST_EXIT(adapt->p >= 1.0, "ADAPT_STAT->p < 1\n");
509
510 td->get_el_estc = adapt->get_el_estc;
511 td->g_mark_refine = adapt->refine_bisections;
512 td->g_mark_coarse = -adapt->coarse_bisections;
513
514 td->mark_flag = 0;
515 td->el_mark = td->el_mark_c = 0;
516
517 eps_p = pow(adapt->tolerance, adapt->p);
518 td->err_sum = pow(adapt->err_sum, adapt->p);
519
520 td->err_max = adapt->err_max;
521
522 switch(adapt->strategy) {
523 case GR:
524 if (adapt->err_sum > adapt->tolerance) {
525 mesh_traverse(mesh, -1, CALL_LEAF_EL, marking_fct_1, td);
526 }
527 break;
528
529 case MS:
530 MS_gamma_p = pow(adapt->MS_gamma, adapt->p);
531 if (adapt->coarsen_allowed)
532 MS_gamma_c_p = pow(adapt->MS_gamma_c, adapt->p);
533 else
534 MS_gamma_c_p = -1.0;
535
536 td->mark_r_limit = MS_gamma_p * td->err_max;
537 if (adapt->coarsen_allowed)
538 td->mark_c_limit = MS_gamma_c_p * td->err_max;
539
540 INFO(adapt->info, 4,
541 "start mark_limits: %.3le %.3le err_max = %.3le\n",
542 td->mark_r_limit, td->mark_c_limit, td->err_max);
543
544 mesh_traverse(mesh, -1, CALL_LEAF_EL, marking_fct_2, td);
545 break;
546
547 case ES:
548 ES_theta_p = pow(adapt->ES_theta, adapt->p);
549 td->mark_r_limit = ES_theta_p * eps_p / mesh->n_elements;
550 if (adapt->coarsen_allowed)
551 {
552 ES_theta_c_p = pow(adapt->ES_theta_c, adapt->p);
553 td->mark_c_limit = ES_theta_c_p * eps_p / mesh->n_elements;
554 }
555 else
556 td->mark_c_limit = -1.0;
557
558 INFO(adapt->info, 4,
559 "start mark_limits: %.3le %.3le n_elements = %d\n",
560 td->mark_r_limit, td->mark_c_limit, mesh->n_elements);
561
562 mesh_traverse(mesh, -1, CALL_LEAF_EL, marking_fct_3, td);
563 break;
564
565 case GERS:
566 /*--------------------------------------------------------------------------*/
567 /* Graranteed error reduction strategy */
568 /* using extrapolation for theta */
569 /*--------------------------------------------------------------------------*/
570
571 Ltheta= pow(1.0 - adapt->GERS_theta_star, adapt->p);
572
573 if (td->err_sum < old_err_sum)
574 {
575 improv = td->err_sum/old_err_sum;
576 wanted = 0.8 * eps_p /td->err_sum;
577 redfac = MIN((1.0-wanted)/(1.0-improv),1.0);
578 redfac = MAX(redfac, 0.0);
579
580 if (redfac < 1.0)
581 {
582 Ltheta = Ltheta*redfac;
583 INFO(adapt->info, 2,
584 "GERS: use extrapolated theta_star = %.3lf\n",
585 pow(Ltheta, 1.0/adapt->p));
586 }
587 }
588 old_err_sum = td->err_sum;
589
590 GERS_gamma = 1.0;
591
592 if (Ltheta > 0)
593 {
594 do
595 {
596 td->GERS_sum = 0.0;
597 GERS_gamma -= adapt->GERS_nu;
598 td->mark_r_limit = GERS_gamma * td->err_max;
599 mesh_traverse(mesh, -1, CALL_LEAF_EL, marking_fct_4, td);
600 }
601 while ((GERS_gamma > 0) && (td->GERS_sum < Ltheta * td->err_sum));
602 }
603
604 INFO(adapt->info, 4,
605 "GERS refinement with gamma = %.3lf\n", GERS_gamma);
606
607 if (adapt->coarsen_allowed)
608 {
609 GERS_gamma = 0.3;
610 Ltheta = adapt->GERS_theta_c * eps_p;
611
612 do
613 {
614 td->GERS_sum = 0.0;
615 GERS_gamma -= adapt->GERS_nu;
616 td->mark_c_limit = GERS_gamma * td->err_max;
617 mesh_traverse(mesh, -1, CALL_LEAF_EL, marking_fct_4c, td);
618
619 INFO(adapt->info, 6,
620 "coarse loop: gamma = %.3e, sum = %.3e, limit = %.3e\n",
621 GERS_gamma, td->GERS_sum, Ltheta);
622 }
623 while (td->GERS_sum > Ltheta);
624
625 INFO(adapt->info, 4,
626 "GERS coarsening with gamma = %.3lf\n", GERS_gamma);
627 }
628 break;
629 default:
630 break;
631 }
632
633 INFO(adapt->info, 4, "%d elements marked for refinement\n", td->el_mark);
634 INFO(adapt->info, 4, "%d elements marked for coarsening\n", td->el_mark_c);
635
636 td->mark_flag = 0;
637 if (td->el_mark) td->mark_flag = 1;
638 if (td->el_mark_c) td->mark_flag |= 2;
639
640 return(td->mark_flag);
641 }
642
643 /*--------------------------------------------------------------------------*/
644
645
646 /*--------------------------------------------------------------------------*/
647 /* get_adapt_stat(dim, name, prefix, info, adapt_stat): */
648 /* */
649 /* if adapt_stat is non NULL: */
650 /* allocates a new adpat_stat structure adapt and predefines the */
651 /* parameter entries of this structure; */
652 /* if name is non NULL adpat->name = name */
653 /* */
654 /* Now, either for adapt_stat or the newly allocated adapt_stat the */
655 /* parameter entries of the structures are initialized by GET_PARAMETER() */
656 /* if prefix is non NULL */
657 /* the info parameter is the the first argument of GET_PARAMETER() */
658 /* */
659 /* the return value is a pointer to the initialized adapt_stat structure */
660 /* */
661 /* for the initialization of entry we use the call */
662 /* GET_PARAMETER(info, "prefix->entry", "%{d,f}", &adapt->entry */
663 /* !! You do not have to define all parameters in your init file !! */
664 /* entries that are not defined in the init file have so standard value */
665 /* */
666 /* !!! before a call of get_adapt_stat() with non NULL second argument !!! */
667 /* !!! you have to initialize your parameters via init_parameters() !!! */
668 /*--------------------------------------------------------------------------*/
669
init_strategy(const char * funcName,const char * prefix,int info,ADAPT_STAT * adapt)670 static void init_strategy(const char *funcName, const char *prefix, int info,
671 ADAPT_STAT *adapt)
672 {
673 char key[1024];
674
675 sprintf(key, "%s->strategy", prefix);
676 GET_PARAMETER(info, key, "%d", &adapt->strategy);
677
678 switch (adapt->strategy)
679 {
680 case MS:
681 sprintf(key, "%s->MS_gamma", prefix);
682 GET_PARAMETER(info, key, "%f", &adapt->MS_gamma);
683 if (adapt->coarsen_allowed)
684 {
685 sprintf(key, "%s->MS_gamma_c", prefix);
686 GET_PARAMETER(info, key, "%f", &adapt->MS_gamma_c);
687 }
688 return;
689 case ES:
690 sprintf(key, "%s->ES_theta", prefix);
691 GET_PARAMETER(info, key, "%f", &adapt->ES_theta);
692 if (adapt->coarsen_allowed)
693 {
694 sprintf(key, "%s->ES_theta_c", prefix);
695 GET_PARAMETER(info-1, key, "%f", &adapt->ES_theta_c);
696 }
697 return;
698 case GERS:
699 sprintf(key, "%s->GERS_theta_star", prefix);
700 GET_PARAMETER(info, key, "%f", &adapt->GERS_theta_star);
701 sprintf(key, "%s->GERS_nu", prefix);
702 GET_PARAMETER(info, key, "%f", &adapt->GERS_nu);
703 if (adapt->coarsen_allowed)
704 {
705 sprintf(key, "%s->GERS_theta_c", prefix);
706 GET_PARAMETER(info, key, "%f", &adapt->GERS_theta_c);
707 }
708 return;
709 default:
710 break;
711 }
712 return;
713 }
714
get_adapt_stat(int dim,const char * name,const char * prefix,int info,ADAPT_STAT * adapt_stat)715 ADAPT_STAT *get_adapt_stat(int dim, const char *name,
716 const char *prefix, int info,
717 ADAPT_STAT *adapt_stat)
718 {
719 FUNCNAME("get_adapt_stat");
720 ADAPT_STAT adapt_stand = {NULL, 1.0, 2, 30, 2, NULL, NULL, NULL, NULL,
721 NULL, 0.0, 0.0, NULL, NULL, NULL, NULL,
722 -1 /* refine_bisections */,
723 false /* coarsen_allowed*/,
724 -1 /* coarse_bisections */,
725 FILL_NOTHING /* adaptation_fill_flags */,
726 GR /* strategy */,
727 0.5, 0.1, 0.9, 0.2,
728 0.6, 0.1, 0.1 };
729 char key[1024];
730 ADAPT_STAT *adapt;
731
732 if(dim == 0) {
733 WARNING("Adaption does not make sense for dim == 0!\n");
734 return NULL;
735 }
736
737 adapt_stand.refine_bisections =
738 adapt_stand.coarse_bisections = dim;
739
740 if (adapt_stat)
741 adapt = adapt_stat;
742 else {
743 adapt = MEM_ALLOC(1, ADAPT_STAT);
744 *adapt = adapt_stand;
745 if (name)
746 adapt->name = strdup(name);
747 if (!adapt->name && prefix)
748 adapt->name = strdup(prefix);
749 }
750
751 if (!prefix)
752 return(adapt);
753
754 sprintf(key, "%s->tolerance", prefix);
755 GET_PARAMETER(info-1, key, "%f", &adapt->tolerance);
756 sprintf(key, "%s->p", prefix);
757 GET_PARAMETER(info-2, key, "%f", &adapt->p);
758 sprintf(key, "%s->max_iteration", prefix);
759 GET_PARAMETER(info-1, key, "%d", &adapt->max_iteration);
760 sprintf(key, "%s->info", prefix);
761 GET_PARAMETER(info-1, key, "%d", &adapt->info);
762
763 sprintf(key, "%s->refine_bisections", prefix);
764 GET_PARAMETER(info-2, key, "%d", &adapt->refine_bisections);
765 sprintf(key, "%s->coarsen_allowed", prefix);
766 GET_PARAMETER(info-2, key, "%B", &adapt->coarsen_allowed);
767 if (adapt->coarsen_allowed) {
768 sprintf(key, "%s->coarse_bisections", prefix);
769 GET_PARAMETER(info-2, key, "%d", &adapt->coarse_bisections);
770 }
771 sprintf(key, "%s->adaptation fill flags", prefix);
772 GET_PARAMETER(info-2, key, "%i", &adapt->adaptation_fill_flags);
773
774 init_strategy(funcName, prefix, info-1, adapt);
775
776 return(adapt);
777 }
778
779 /*--------------------------------------------------------------------------*/
780 /* get_adapt_instat(dim, name, prefix, info, adapt_instat): */
781 /* */
782 /* if adapt_instat is non NULL: */
783 /* allocates a new adpat_instat structure adapt and predefines the */
784 /* parameter entries of this structure; */
785 /* if name is non NULL adapt->name = name */
786 /* */
787 /* Now, either for adapt_instat or the newly allocated adapt_instat the */
788 /* parameter entries of the structures are initialized by GET_PARAMETER() */
789 /* if prefix is non NULL */
790 /* the info parameter is the the first argument of GET_PARAMETER() */
791 /* */
792 /* the return value is a pointer to the initialized adapt_stat structure */
793 /* */
794 /* for the initialization of entry we use the call */
795 /* GET_PARAMETER(info, "prefix->entry", "%{d,f}", &adapt->entry) */
796 /* sub-structures adapt_initial/adapt_space are initialized with prefix */
797 /* prefix->initial resp. prefix->space */
798 /* */
799 /* !! You do not have to define all parameters in your init file !! */
800 /* entries that are not defined in the init file have so standard value */
801 /* */
802 /* !!! before a call of get_adapt_stat() with non NULL second argument !!! */
803 /* !!! you have to initialize your parameters via init_parameters() !!! */
804 /*--------------------------------------------------------------------------*/
805
get_adapt_instat(int dim,const char * name,const char * prefix,int info,ADAPT_INSTAT * adapt_instat)806 ADAPT_INSTAT *get_adapt_instat(int dim, const char *name,
807 const char *prefix, int info,
808 ADAPT_INSTAT *adapt_instat)
809 {
810 FUNCNAME("get_adapt_instat");
811 ADAPT_INSTAT adapt_stand = {NULL,
812 {{NULL,
813 1.0, 2.0, 1, -1, NULL, NULL, NULL, NULL,
814 NULL, 0.0, 0.0, NULL, NULL, NULL, NULL,
815 -1, 0, -1, 2, 0.5, 0.1, 0.9, 0.2,
816 0.6, 0.1, 0.1}},
817 {{NULL,
818 1.0, 2.0, 1, -1, NULL, NULL, NULL, NULL,
819 NULL, 0.0, 0.0, NULL, NULL, NULL, NULL,
820 -1, 1, -1, 2, 0.5, 0.1, 0.9, 0.2,
821 0.6, 0.1, 0.1}},
822 0.0, 0.0, 1.0, 0.01, NULL, NULL, NULL, NULL, NULL,
823 0, 0, 1.0, 0.1, 0.4, 0.4, 1.0, 0.3,
824 0.7071, 1.4142, 8};
825 char key[1024];
826 ADAPT_INSTAT *adapt;
827
828 if(dim == 0) {
829 WARNING("Adaption does not make sense for dim == 0!\n");
830 return NULL;
831 }
832
833 adapt_stand.adapt_initial->refine_bisections =
834 adapt_stand.adapt_initial->coarse_bisections =
835 adapt_stand.adapt_space->refine_bisections =
836 adapt_stand.adapt_space->coarse_bisections = dim;
837
838 if (adapt_instat)
839 adapt = adapt_instat;
840 else
841 {
842 adapt = MEM_ALLOC(1, ADAPT_INSTAT);
843 *adapt = adapt_stand;
844 if (name)
845 adapt->name = strdup(name);
846 if (!adapt->name && prefix)
847 adapt->name = strdup(prefix);
848 }
849
850 if (!prefix)
851 return(adapt);
852
853 sprintf(key, "%s initial", adapt->name);
854 adapt->adapt_initial->name = strdup(key);
855 sprintf(key, "%s space", adapt->name);
856 adapt->adapt_space->name = strdup(key);
857
858 /*---8<---------------------------------------------------------------------*/
859 /*--- and now, all other entries ---*/
860 /*--------------------------------------------------------------------->8---*/
861
862 sprintf(key, "%s->start_time", prefix);
863 GET_PARAMETER(info-1, key, "%f", &adapt->start_time);
864 adapt->time = adapt->start_time;
865 sprintf(key, "%s->end_time", prefix);
866 GET_PARAMETER(info-1, key, "%f", &adapt->end_time);
867 sprintf(key, "%s->timestep", prefix);
868 GET_PARAMETER(info-1, key, "%f", &adapt->timestep);
869 sprintf(key, "%s->strategy", prefix);
870 GET_PARAMETER(info-1, key, "%d", &adapt->strategy);
871 sprintf(key, "%s->max_iteration", prefix);
872 GET_PARAMETER(info-1, key, "%d", &adapt->max_iteration);
873 sprintf(key, "%s->tolerance", prefix);
874 GET_PARAMETER(info-1, key, "%f", &adapt->tolerance);
875 sprintf(key, "%s->rel_initial_error", prefix);
876 GET_PARAMETER(info-1, key, "%f", &adapt->rel_initial_error);
877 sprintf(key, "%s->rel_space_error", prefix);
878 GET_PARAMETER(info-1, key, "%f", &adapt->rel_space_error);
879 sprintf(key, "%s->rel_time_error", prefix);
880 GET_PARAMETER(info-1, key, "%f", &adapt->rel_time_error);
881 sprintf(key, "%s->time_theta_1", prefix);
882 GET_PARAMETER(info-2, key, "%f", &adapt->time_theta_1);
883 sprintf(key, "%s->time_theta_2", prefix);
884 GET_PARAMETER(info-2, key, "%f", &adapt->time_theta_2);
885 sprintf(key, "%s->time_delta_1", prefix);
886 GET_PARAMETER(info-2, key, "%f", &adapt->time_delta_1);
887 sprintf(key, "%s->time_delta_2", prefix);
888 GET_PARAMETER(info-2, key, "%f", &adapt->time_delta_2);
889 sprintf(key, "%s->info", prefix);
890 GET_PARAMETER(info-1, key, "%d", &adapt->info);
891
892 /*---8<---------------------------------------------------------------------*/
893 /*--- initialization of the adapt_stat for the initial grid ---*/
894 /*--------------------------------------------------------------------->8---*/
895
896 /*--- tolerance does not have to be initialized, is set! ---*/
897 adapt->adapt_initial->tolerance = adapt->tolerance*adapt->rel_initial_error;
898 sprintf(key, "%s->initial->p", prefix);
899 GET_PARAMETER(info-2, key, "%f", &adapt->adapt_initial->p);
900 sprintf(key, "%s->initial->max_iteration", prefix);
901 GET_PARAMETER(info-1, key, "%d", &adapt->adapt_initial->max_iteration);
902 sprintf(key, "%s->initial->info", prefix);
903 GET_PARAMETER(info-2, key, "%d", &adapt->adapt_initial->info);
904 if (adapt->adapt_initial->info < 0)
905 adapt->adapt_initial->info = adapt->info-2;
906
907 sprintf(key, "%s->initial->refine_bisections", prefix);
908 GET_PARAMETER(info-2, key, "%d", &adapt->adapt_initial->refine_bisections);
909 sprintf(key, "%s->initial->coarsen_allowed", prefix);
910 GET_PARAMETER(info-2, key, "%B", &adapt->adapt_initial->coarsen_allowed);
911 if (adapt->adapt_initial->coarsen_allowed) {
912 sprintf(key, "%s->initial->coarse_bisections", prefix);
913 GET_PARAMETER(info-2, key, "%d", &adapt->adapt_initial->coarse_bisections);
914 }
915 sprintf(key, "%s->initial", prefix);
916 init_strategy(funcName, key, info-1, adapt->adapt_initial);
917
918
919 /*---8<---------------------------------------------------------------------*/
920 /*--- initialization of the adapt_stat for the time-step iteration ---*/
921 /*--------------------------------------------------------------------->8---*/
922
923 /*--- tolerance does not have to be initialized, is set! ---*/
924 adapt->adapt_space->tolerance = adapt->tolerance*adapt->rel_space_error;
925 sprintf(key, "%s->space->p", prefix);
926 GET_PARAMETER(info-2, key, "%f", &adapt->adapt_space->p);
927 sprintf(key, "%s->space->max_iteration", prefix);
928 GET_PARAMETER(info-1, key, "%d", &adapt->adapt_space->max_iteration);
929 sprintf(key, "%s->space->info", prefix);
930 GET_PARAMETER(info-2, key, "%d", &adapt->adapt_space->info);
931 if (adapt->adapt_space->info < 0)
932 adapt->adapt_space->info = adapt->info-2;
933
934 sprintf(key, "%s->space->refine_bisections", prefix);
935 GET_PARAMETER(info-2, key, "%d", &adapt->adapt_space->refine_bisections);
936 sprintf(key, "%s->space->coarsen_allowed", prefix);
937 GET_PARAMETER(info-2, key, "%B", &adapt->adapt_space->coarsen_allowed);
938 if (adapt->adapt_space->coarsen_allowed) {
939 sprintf(key, "%s->space->coarse_bisections", prefix);
940 GET_PARAMETER(info-2, key, "%d", &adapt->adapt_space->coarse_bisections);
941 }
942 sprintf(key, "%s->space", prefix);
943 init_strategy(funcName, key, info-1, adapt->adapt_space);
944
945 return(adapt);
946 }
947