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 /* */
8 /* file: heat.c */
9 /* */
10 /* description: solver for parabolic model problem */
11 /* */
12 /* u,t - \Delta u = f in \Omega */
13 /* u = g on \partial \Omega */
14 /* */
15 /*--------------------------------------------------------------------------*/
16 /* */
17 /* authors: Alfred Schmidt */
18 /* Zentrum fuer Technomathematik */
19 /* Fachbereich 3 Mathematik/Informatik */
20 /* Universitaet Bremen */
21 /* Bibliothekstr. 2 */
22 /* D-28359 Bremen, Germany */
23 /* */
24 /* Kunibert G. Siebert */
25 /* Institut fuer Mathematik */
26 /* Universitaet Augsburg */
27 /* Universitaetsstr. 14 */
28 /* D-86159 Augsburg, Germany */
29 /* */
30 /* http://www.alberta-fem.de/ */
31 /* */
32 /*--------------------------------------------------------------------------*/
33 /* (c) by A. Schmidt and K.G. Siebert (1996-2004) */
34 /* */
35 /* This program is free software; you can redistribute it and/or modify */
36 /* it under the terms of the GNU General Public License as published by */
37 /* the Free Software Foundation; either version 2 of the License, or */
38 /* any later version. */
39 /* */
40 /* This program is distributed in the hope that it will be useful, */
41 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
42 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
43 /* GNU General Public License for more details. */
44 /*--------------------------------------------------------------------------*/
45
46 #include <alberta.h>
47
48 /*--------------------------------------------------------------------------*/
49 /* function for displaying mesh, discrete solution, and/or estimate */
50 /* defined in graphics.c */
51 /*--------------------------------------------------------------------------*/
52 /*--------------------------------------------------------------------------*/
53 /* global variables: finite element space, discrete solution, discrete */
54 /* solution from previous time step, load vector, */
55 /* system matrix, and structure for adaptive procedure */
56 /*--------------------------------------------------------------------------*/
57 #include "demoheat.h"
58 #include "dxexport.hh"
59 static REAL time_est;/*muss hier hin!*/
60 #include "demographics.h"
61 #include "problem.cc"
62
63 /*--------------------------------------------------------------------------*/
64 /* init_dof_admin(): init DOFs for Lagrange elements of order */
65 /* <polynomial degree>, once called by GET_MESH(), now by main */
66 /*--------------------------------------------------------------------------*/
67
init_dof_admin(MESH * mesh)68 static void init_dof_admin(MESH *mesh)
69 {
70 FUNCNAME("init_dof_admin");
71 int degree = 1;
72 const BAS_FCTS *lagrange, *lagrange0;
73
74 GET_PARAMETER(1, "polynomial degree", "%d", °ree);
75 lagrange = get_lagrange(mesh->dim, degree);
76 TEST_EXIT(lagrange, "no lagrange BAS_FCTS\n");
77
78 fe_space = get_fe_space(mesh, lagrange->name, NULL, lagrange, false);
79 lagrange0 = get_lagrange(mesh->dim,0);
80 const_fe_space = get_fe_space(mesh, lagrange0->name,NULL, lagrange0, false);
81 #if USE_NCDF
82 init_export_dof_admin(mesh);
83 #endif
84 return;
85 }
86
87
88 /*--------------------------------------------------------------------------*/
89 /* struct heat_leaf_data: structure for storing one REAL value on each */
90 /* leaf element as LEAF_DATA */
91 /* init_leaf_data(): initialization of leaf data, called by GET_MESH() */
92 /* rw_el_est(): return a pointer to the memory for storing the element */
93 /* estimate (stored as LEAF_DATA), called by heat_est() */
94 /* get_el_est(): return the value of the element estimates (from LEAF_DATA),*/
95 /* called by adapt_method_stat() and graphics() */
96 /*--------------------------------------------------------------------------*/
97
98 struct heat_leaf_data
99 {
100 REAL estimate; /* one real for the estimate */
101 REAL est_c; /* one real for the coarsening estimate */
102 };
103 #if 0
104 static void init_leaf_data(LEAF_DATA_INFO *leaf_data_info)
105 {
106 leaf_data_info->leaf_data_size = sizeof(struct heat_leaf_data);
107 leaf_data_info->coarsen_leaf_data = nil; /* no transformation */
108 leaf_data_info->refine_leaf_data = nil; /* no transformation */
109 return;
110 #endif
111 #define USE_DOF_ERR 1
112 #if USE_DOF_ERR
113 /*Kopie aus est.cc*/
get_el_est(EL * el)114 static REAL get_el_est(EL *el)
115 {REAL err_est;
116 if (IS_LEAF_EL(el))
117 {
118 (const_fe_space->bas_fcts->get_real_vec)(el, error_estimate, &err_est);
119 return(err_est);
120 }
121 else
122 return(0.0);
123 }
124
rw_el_est(EL * el)125 static REAL * rw_el_est(EL *el)
126 {
127 static const DOF *(*get_dof)(const EL *, const DOF_ADMIN *, DOF *)
128 = const_fe_space->bas_fcts->get_dof_indices;
129 if (IS_LEAF_EL(el))
130 {
131 const DOF *dof = (*get_dof)(el, const_fe_space->admin, nil);
132 return(&error_estimate->vec[*dof]);
133 }
134 else
135 return(nil);
136 }
137 #else
rw_el_est(EL * el)138 static REAL *rw_el_est(EL *el)
139 {
140 if (IS_LEAF_EL(el))
141 return(&((struct heat_leaf_data *)LEAF_DATA(el))->estimate);
142 else
143 return(nil);
144 }
get_el_est(EL * el)145 static REAL get_el_est(EL *el)
146 {
147 if (IS_LEAF_EL(el))
148 return(((struct heat_leaf_data *)LEAF_DATA(el))->estimate);
149 else
150 return(0.0);
151 }
152 #endif
153
rw_el_estc(EL * el)154 static REAL *rw_el_estc(EL *el)
155 {
156 if (IS_LEAF_EL(el))
157 return(&((struct heat_leaf_data *)LEAF_DATA(el))->est_c);
158 else
159 return(nil);
160 }
161
get_el_estc(EL * el)162 static REAL get_el_estc(EL *el)
163 {
164 if (IS_LEAF_EL(el))
165 return(((struct heat_leaf_data *)LEAF_DATA(el))->est_c);
166 else
167 return(0.0);
168 }
169
170 //static REAL time_est = 0.0;
171
get_time_est(MESH * mesh,ADAPT_INSTAT * adapt)172 static REAL get_time_est(MESH *mesh, ADAPT_INSTAT *adapt)
173 {
174 return(time_est);
175 }
176
177 /*--------------------------------------------------------------------------*/
178 /* For test purposes: exact solution and its gradient (optional) */
179 /*--------------------------------------------------------------------------*/
180
181
182
183
184 /*---8<---------------------------------------------------------------------*/
185 /*--- write error and estimator data to files ---*/
186 /*--------------------------------------------------------------------->8---*/
187
188
189 /*---8<---------------------------------------------------------------------*/
190 /*--- interpolation is solve on the initial grid ---*/
191 /*--------------------------------------------------------------------->8---*/
192
interpol_u0(MESH * mesh)193 static void interpol_u0(MESH *mesh)
194 {
195 dof_compress(mesh);
196 interpol(u0, u_h);
197
198 return;
199 }
200
init_timestep(MESH * mesh,ADAPT_INSTAT * adapt)201 static void init_timestep(MESH *mesh, ADAPT_INSTAT *adapt)
202 {
203 FUNCNAME("init_timestep");
204
205 INFO(adapt_instat->info,1,
206 "---8<---------------------------------------------------\n");
207 INFO(adapt_instat->info, 1,"starting new timestep\n");
208
209 dof_copy(u_h, u_old);
210 return;
211 }
212
set_time(MESH * mesh,ADAPT_INSTAT * adapt)213 static void set_time(MESH *mesh, ADAPT_INSTAT *adapt)
214 {
215 FUNCNAME("set_time");
216
217 INFO(adapt->info,1,
218 "---8<---------------------------------------------------\n");
219 if (adapt->time == adapt->start_time)
220 {
221 INFO(adapt->info, 1,"start time: %.4le\n", adapt->time);
222 }
223 else
224 {
225 INFO(adapt->info, 1,"timestep for (%.4le %.4le), tau = %.4le\n",
226 adapt->time-adapt->timestep, adapt->time,
227 adapt->timestep);
228 }
229
230 eval_time_f = adapt->time - (1 - theta)*adapt->timestep;
231 eval_time_g = adapt->time;
232
233 return;
234 }
235
236
237
238 /*--------------------------------------------------------------------------*/
239 /* build(): assemblage of the linear system: matrix, load vector, */
240 /* boundary values, called by adapt_method_stat() */
241 /* on the first call initialize u_h, f_h, matrix and information */
242 /* for assembling the system matrix */
243 /* */
244 /* struct op_info: structure for passing information from init_element() to */
245 /* LALt() */
246 /* init_element(): initialization on the element; calculates the */
247 /* coordinates and |det DF_S| used by LALt; passes these */
248 /* values to LALt via user_data, */
249 /* called on each element by update_matrix() */
250 /* LALt(): implementation of -Lambda id Lambda^t for -Delta u, */
251 /* called by update_matrix() */
252 /* c(): implementation of 1/tau*m(,) */
253 /*--------------------------------------------------------------------------*/
254
255 struct op_info
256 {
257 REAL_D Lambda[DIM+1]; /* the gradient of the barycentric coordinates */
258 REAL det; /* |det D F_S| */
259
260 REAL tau_1;
261 };
262
LALt(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)263 static const REAL (*LALt(const EL_INFO *el_info, const QUAD *quad,
264 int iq, void *ud))[DIM+1]
265 {
266 struct op_info *info =(struct op_info *) ud;
267 int i, j, k;
268 static REAL LALt[DIM+1][DIM+1];
269
270 for (i = 0; i <= DIM; i++)
271 for (j = i; j <= DIM; j++)
272 {
273 for (LALt[i][j] = k = 0; k < DIM_OF_WORLD; k++)
274 LALt[i][j] += info->Lambda[i][k]*info->Lambda[j][k];
275 LALt[i][j] *= info->det;
276 LALt[j][i] = LALt[i][j];
277 }
278 return((const REAL (*)[DIM+1]) LALt);
279 }
280
c(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)281 static REAL c(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
282 {
283 struct op_info *info = (struct op_info *)ud;
284
285 return(info->tau_1*info->det);
286 }
287
assemble(DOF_REAL_VEC * u_old,DOF_MATRIX * matrix,DOF_REAL_VEC * fh,DOF_REAL_VEC * u_h,REAL theta,REAL tau,REAL (* f)(const REAL_D),REAL (* g)(const REAL_D))288 static void assemble(DOF_REAL_VEC *u_old, DOF_MATRIX *matrix, DOF_REAL_VEC *fh,
289 DOF_REAL_VEC *u_h, REAL theta, REAL tau,
290 REAL (*f)(const REAL_D), REAL (*g)(const REAL_D))
291 {
292 FUNCNAME("assemble");
293 static struct op_info *op_info = nil;
294 // static const REAL **(*fill_a)(const EL_INFO *, void *) = nil;
295 static const REAL *const*(*fill_a)(const EL_INFO *, void *) = NULL;
296 static void *a_info = nil;
297 // static const REAL **(*fill_c)(const EL_INFO *, void *) = nil;
298 static const REAL *const*(*fill_c)(const EL_INFO *, void *) = NULL;
299 static void *c_info = nil;
300
301 static const DOF_ADMIN *admin = nil;
302 static int n;
303 static const REAL *(*get_u_loc)(const EL *, const DOF_REAL_VEC *, REAL *);
304 static const DOF *(*get_dof)(const EL *, const DOF_ADMIN *, DOF *);
305 static const S_CHAR *(*get_bound)(const EL_INFO *, S_CHAR *);
306
307 TRAVERSE_STACK *stack = get_traverse_stack();
308 const EL_INFO *el_info;
309 FLAGS fill_flag;
310 // const REAL **a_mat, **c_mat;
311 const REAL *const*a_mat, *const*c_mat;
312 REAL *f_vec;
313 const QUAD *quad;
314 int i, j;
315
316 quad = get_quadrature(DIM, 2*u_h->fe_space->bas_fcts->degree);
317
318 /*--------------------------------------------------------------------------*/
319 /* init functions for matrix assembling */
320 /*--------------------------------------------------------------------------*/
321
322 if (admin != u_h->fe_space->admin)
323 {
324 OPERATOR_INFO o_info2 = {nil}, o_info0 = {nil};
325 const EL_MATRIX_INFO *matrix_info;
326 const BAS_FCTS *bas_fcts = u_h->fe_space->bas_fcts;
327
328 admin = u_h->fe_space->admin;
329
330 n = bas_fcts->n_bas_fcts;
331 get_dof = bas_fcts->get_dof_indices;
332 get_bound = bas_fcts->get_bound;
333 get_u_loc = bas_fcts->get_real_vec;
334
335 if (!op_info)
336 op_info = MEM_ALLOC(1, struct op_info);
337
338 o_info2.row_fe_space = o_info2.col_fe_space = u_h->fe_space;
339
340 o_info2.quad[2] = quad;
341 o_info2.LALt = LALt;
342 o_info2.LALt_pw_const = true;
343 o_info2.LALt_symmetric = true;
344 o_info2.user_data = op_info;
345
346 matrix_info = fill_matrix_info(&o_info2, nil);
347 fill_a = matrix_info->el_matrix_fct;
348 a_info = matrix_info->fill_info;
349
350 o_info0.row_fe_space = o_info0.col_fe_space = u_h->fe_space;
351
352 o_info0.quad[0] = quad;
353 o_info0.c = c;
354 o_info0.c_pw_const = true;
355 o_info0.user_data = op_info;
356
357 matrix_info = fill_matrix_info(&o_info0, nil);
358 fill_c = matrix_info->el_matrix_fct;
359 c_info = matrix_info->fill_info;
360 }
361
362 op_info->tau_1 = 1.0/tau;
363
364 /*--------------------------------------------------------------------------*/
365 /* and now assemble the matrix and right hand side */
366 /*--------------------------------------------------------------------------*/
367
368 clear_dof_matrix(matrix);
369 dof_set(0.0, fh);
370 f_vec = fh->vec;
371
372 fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_BOUND;
373 el_info = traverse_first(stack, u_h->fe_space->mesh, -1, fill_flag);
374 while (el_info)
375 {
376 const REAL *u_old_loc = (*get_u_loc)(el_info->el, u_old, nil);
377 const DOF *dof = (*get_dof)(el_info->el, admin, nil);
378 const S_CHAR *bound = (*get_bound)(el_info, nil);
379
380 /*--------------------------------------------------------------------------*/
381 /* initialization of values used by LALt and c */
382 /*--------------------------------------------------------------------------*/
383 op_info->det = el_grd_lambda(el_info, op_info->Lambda);
384
385 a_mat = fill_a(el_info, a_info);
386 c_mat = fill_c(el_info, c_info);
387
388 /*--------------------------------------------------------------------------*/
389 /* add theta*a(psi_i,psi_j) + 1/tau*m(4*u^3*psi_i,psi_j) */
390 /*--------------------------------------------------------------------------*/
391
392 if (theta)
393 {
394 add_element_matrix(matrix, theta, n, n, dof, dof, a_mat, bound);
395 }
396 add_element_matrix(matrix, 1.0, n, n, dof, dof, c_mat, bound);
397
398 /*--------------------------------------------------------------------------*/
399 /* f += -(1-theta)*a(u_old,psi_i) + 1/tau*m(u_old,psi_i) */
400 /*--------------------------------------------------------------------------*/
401
402 if (1.0 - theta)
403 {
404 REAL theta1 = 1.0 - theta;
405 for (i = 0; i < n; i++)
406 {
407 if (bound[i] < DIRICHLET)
408 {
409 REAL val = 0.0;
410 for (j = 0; j < n; j++)
411 val += (-theta1*a_mat[i][j] + c_mat[i][j])*u_old_loc[j];
412 f_vec[dof[i]] += val;
413 }
414 }
415 }
416 else
417 {
418 for (i = 0; i < n; i++)
419 {
420 if (bound[i] < DIRICHLET)
421 {
422 REAL val = 0.0;
423 for (j = 0; j < n; j++)
424 val += c_mat[i][j]*u_old_loc[j];
425 f_vec[dof[i]] += val;
426 }
427 }
428 }
429 el_info = traverse_next(stack, el_info);
430 }
431
432 free_traverse_stack(stack);
433
434 L2scp_fct_bas(f, quad, fh);
435 dirichlet_bound(g, fh, u_h, nil);
436
437 return;
438 }
439
build(MESH * mesh,U_CHAR flag)440 static void build(MESH *mesh, U_CHAR flag)
441 {
442 FUNCNAME("build");
443
444 dof_compress(mesh);
445
446 INFO(adapt_instat->adapt_space->info, 2 ,"%d DOFs for %s\n", fe_space->admin->size_used, fe_space->name);
447
448 assemble(u_old, matrix, f_h, u_h, theta, adapt_instat->timestep, f, g);
449 if (adapt_instat->adapt_space->info >= 8)
450 {
451 print_dof_real_vec(f_h);
452 if (adapt_instat->adapt_space->info >= 9)
453 print_dof_matrix(matrix);
454 }
455 return;
456 }
457
458 /*--------------------------------------------------------------------------*/
459 /* solve(): solve the linear system, called by adapt_method_stat() */
460 /*--------------------------------------------------------------------------*/
461
solve(MESH * mesh)462 static void solve(MESH *mesh)
463 {
464 FUNCNAME("solve");
465 static REAL tol = 1.e-8;
466 static int max_iter = 1000, info = 2, icon = 1, restart = 0;
467 static OEM_SOLVER solver = NoSolver;
468
469 if (solver == NoSolver)
470 {
471 tol = 1.e-8;
472 GET_PARAMETER(1, "heat solver", "%d", &solver);
473 GET_PARAMETER(1, "heat solver tolerance", "%f", &tol);
474 GET_PARAMETER(1, "heat solver precon", "%d", &icon);
475 GET_PARAMETER(1, "heat solver max iteration", "%d", &max_iter);
476 GET_PARAMETER(1, "heat solver info", "%d", &info);
477 if (solver == GMRes)
478 GET_PARAMETER(1, "solver restart", "%d", &restart);
479 }
480 //print_dof_real_vec(u_h);
481 //print_dof_real_vec(f_h);
482 oem_solve_s(matrix, NULL/*bound*/,f_h, u_h, solver, tol, (OEM_PRECON)icon, restart, max_iter, info);
483
484 return;
485 }
486
487 /*--------------------------------------------------------------------------*/
488 /* Functions for error estimate: */
489 /* estimate(): calculates error estimate via heat_est() */
490 /* calculates exact error also (only for test purpose), */
491 /* called by adapt_method_stat() */
492 /*--------------------------------------------------------------------------*/
493
494 /*static REAL r(const EL_INFO *el_info, const QUAD *quad, int iq, REAL t,
495 REAL uh_iq, const REAL_D grd_uh_iq)
r(const EL_INFO * el_info,const QUAD * quad,int iq,REAL uh_at_qp,const REAL_D grd_uh_at_qp,REAL t)496 */static REAL r(const EL_INFO *el_info,
497 const QUAD *quad, int iq,
498 REAL uh_at_qp, const REAL_D grd_uh_at_qp,
499 REAL t)
500 {
501 REAL_D x;
502 coord_to_world(el_info, quad->lambda[iq], x);
503 eval_time_f = t;
504 return(-f(x));
505 }
506
507
estimate(MESH * mesh,ADAPT_STAT * adapt)508 static REAL estimate(MESH *mesh, ADAPT_STAT *adapt)
509 {
510 FUNCNAME("estimate");
511 static int degree;
512 static REAL C[4] = {-1.0, 1.0, 1.0, 1.0};
513 REAL_DD A = {{0.0}};
514 FLAGS r_flag = 0; /* = (INIT_UH|INIT_GRD_UH), if needed by r() */
515 int n;
516 REAL space_est;
517
518 for (n = 0; n < DIM_OF_WORLD; n++)
519 A[n][n] = 1.0; /* set diogonal of A; all other elements are zero */
520
521 eval_time_u = adapt_instat->time;
522
523 if (C[0] < 0)
524 {
525 C[0] = 1.0;
526 GET_PARAMETER(1, "th estimator C0", "%f", &C[0]);
527 GET_PARAMETER(1, "th estimator C1", "%f", &C[1]);
528 GET_PARAMETER(1, "th estimator C2", "%f", &C[2]);
529 GET_PARAMETER(1, "th estimator C3", "%f", &C[3]);
530 }
531
532 degree = 2*u_h->fe_space->bas_fcts->degree;
533 time_est = heat_est(u_h, adapt_instat, rw_el_est, rw_el_estc,
534 degree, C, u_old, (const REAL_D *) A, r, r_flag
535 , NULL /* gn() */, 0 /* gn_flag */);
536 #if USE_DOF_ERR
537 // print_dof_real_vec(error_estimate);
538 #endif
539 space_est = adapt_instat->adapt_space->err_sum;
540 // err_L2 = L2_err(u, u_h, nil, 0, nil, nil);
541 err_L2 = L2_err(u, u_h, NULL, false, false, NULL, NULL);
542
543 INFO(adapt_instat->info,2
544 ,"---8<---------------------------------------------------\n");
545 INFO(adapt_instat->info, 6,"time_est = %.4lg\n", time_est);
546 INFO(adapt_instat->info, 2,"time = %.4le with timestep = %.4le\n",
547 adapt_instat->time, adapt_instat->timestep);
548 INFO(adapt_instat->info, 2,"||u-uh||L2 = %.4le, ratio = %.2lf\n", err_L2,
549 err_L2/MAX(space_est,1.e-20));
550 INFO(adapt_instat->info, 2,"estimate = %.4le, max = %.4le\n", space_est,
551 sqrt(adapt_instat->adapt_space->err_max));
552
553 return(adapt_instat->adapt_space->err_sum);
554 }
555
est_initial(MESH * mesh,ADAPT_STAT * adapt)556 static REAL est_initial(MESH *mesh, ADAPT_STAT *adapt)
557 {
558 FUNCNAME("est_initial");
559 err_L2 = adapt->err_sum = L2_err(u0, u_h, nil, false, false ,rw_el_est, &adapt->err_max);
560 MSG("L2 error = %.3le\n", adapt->err_sum);
561 return(adapt->err_sum);
562 }
563
564 /*--------------------------------------------------------------------------*/
565 /* main program */
566 /*--------------------------------------------------------------------------*/
567
main(int argc,char ** argv)568 int main(int argc, char **argv)
569 {
570 FUNCNAME("main");
571 MACRO_DATA *data;
572 MESH *mesh;
573 int n_refine = 0, k, p = 1, dim=2;
574 char filename[128];
575 REAL fac = 1.0;
576
577 /*--------------------------------------------------------------------------*/
578 /* first of all, init parameters of the init file */
579 /*--------------------------------------------------------------------------*/
580
581 init_parameters(0, "INIT/heat.dat");
582 for (k = 1; k+1 < argc; k += 2)
583 ADD_PARAMETER(0, argv[k], argv[k+1]);
584
585 /*--------------------------------------------------------------------------*/
586 /* get a mesh, and read the macro triangulation from file */
587 /*--------------------------------------------------------------------------*/
588 //von 2.0
589 GET_PARAMETER(1, "mesh dimension", "%d", &dim);
590 GET_PARAMETER(1, "macro file name", "%s", filename);
591 GET_PARAMETER(1, "global refinements", "%d", &n_refine);
592 GET_PARAMETER(1, "polynomial degree", "%d", &p);
593 data = read_macro(filename);
594
595 mesh = GET_MESH(dim, "ALBERTA mesh", data, NULL, NULL);
596 free_macro_data(data);
597 init_dof_admin( mesh);
598 init_leaf_data(mesh, sizeof(struct heat_leaf_data), NULL, NULL);
599 /*init_leaf_data(MESH *mesh, size_t size,
600 void (*refine_leaf_data)(EL *parent,
601 EL *child[2]),
602 void (*coarsen_leaf_data)(EL *parent,
603 EL *child[2]));
604 */
605
606 /* mesh = GET_MESH("ALBERTA mesh", init_dof_admin, init_leaf_data);
607 read_macro(mesh, filename, nil);
608 */ global_refine(mesh, n_refine*DIM);
609 graphics(mesh, nil, nil, nil, 0.0);
610
611 GET_PARAMETER(1, "parameter theta", "%e", &theta);
612 if (theta < 0.5)
613 {
614 WARNING("You are using the explicit Euler scheme\n");
615 WARNING("Use a sufficiently small time step size!!!\n");
616 fac = 1.0e-3;
617 }
618
619 matrix = get_dof_matrix("A", fe_space, NULL);
620 f_h = get_dof_real_vec("f_h", fe_space);
621 u_h = get_dof_real_vec("u_h", fe_space);
622 u_h->refine_interpol = fe_space->bas_fcts->real_refine_inter;
623 u_h->coarse_restrict = fe_space->bas_fcts->real_coarse_inter;
624 u_old = get_dof_real_vec("u_old", fe_space);
625 u_old->refine_interpol = fe_space->bas_fcts->real_refine_inter;
626 u_old->coarse_restrict = fe_space->bas_fcts->real_coarse_inter;
627 dof_set(0.0, u_h); /* initialize u_h ! */
628 /*in order to be able to print vector of error estimate. Switch between element data and
629 dof_real_vec mode on top of get_el_est definition.*/
630 error_estimate = get_dof_real_vec("error_estimate", const_fe_space);
631 error_estimate->refine_interpol = nil;
632 error_estimate->coarse_restrict = nil;
633 dof_set(0.0, error_estimate); /* initialize u_h ! */
634 /*--------------------------------------------------------------------------*/
635 /* init adapt structure and start adaptive method */
636 /*--------------------------------------------------------------------------*/
637
638 adapt_instat = get_adapt_instat(dim, "heat", "adapt", 2, adapt_instat);
639
640 /*--------------------------------------------------------------------------*/
641 /* adapt time step size to refinement level and polynomial degree */
642 /*--------------------------------------------------------------------------*/
643
644
645 // if (theta == 0.5)
646 // adapt_instat->timestep *= fac*pow(2, -(REAL)(p*(n_refine))/2.0);
647 // else
648 // adapt_instat->timestep *= fac*pow(2, -(REAL)(p*(n_refine)));
649 MSG("using initial timestep size = %.4le\n", adapt_instat->timestep);
650
651 eval_time_u0 = adapt_instat->start_time;
652
653 adapt_instat->adapt_initial->get_el_est = get_el_est;
654 adapt_instat->adapt_initial->estimate = est_initial;
655 adapt_instat->adapt_initial->solve = interpol_u0;
656
657 adapt_instat->adapt_space->get_el_est = get_el_est;
658 adapt_instat->adapt_space->get_el_estc = get_el_estc;
659 adapt_instat->adapt_space->estimate = estimate;
660 adapt_instat->adapt_space->build_after_coarsen = build;
661 adapt_instat->adapt_space->solve = solve;
662
663 adapt_instat->init_timestep = init_timestep;
664 adapt_instat->set_time = set_time;
665 adapt_instat->get_time_est = get_time_est;
666 adapt_instat->close_timestep = close_timestep;
667
668 adapt_method_instat(mesh, adapt_instat);
669
670 WAIT_REALLY;
671 return(0);
672 }
673