src/Common/heat.c

/*******************************************************************************
 *
 * ALBERTA:  an Adaptive multi Level finite element toolbox  using
 *           Bisectioning refinement and Error control by Residual
 *           Techniques
 *
 * www.alberta-fem.de
 *
 * file:     heat.c
 *
 * description:  solver for parabolic model problem
 *
 *                   u,t - \Delta u = f  in \Omega
 *                                u = g  on \partial \Omega
 *
 * This version demonstrates the use of EL_SYS_INFO_INSTAT together with
 * update_system_instat()
 *
 ******************************************************************************
 *
 *  authors:   Alfred Schmidt
 *             Zentrum fuer Technomathematik
 *             Fachbereich 3 Mathematik/Informatik
 *             Universitaet Bremen
 *             Bibliothekstr. 2
 *             D-28359 Bremen, Germany
 *
 *             Kunibert G. Siebert
 *             Institut fuer Mathematik
 *             Universitaet Augsburg
 *             Universitaetsstr. 14
 *             D-86159 Augsburg, Germany
 *
 *  (c) by A. Schmidt and K.G. Siebert (1996-2003)
 *
 ******************************************************************************/


#include <alberta/alberta.h>

#include "alberta-demo.h"


static bool do_graphics = true; /* global graphics "kill-switch" */

/*******************************************************************************
 * global variables: finite element space, discrete solution, discrete
 *                   solution from previous time step, load vector,
 *                   system matrix, and structure for adaptive procedure
 *
 * These variables are kept global because they are shared across build(),
 * solve() and estimate().
 ******************************************************************************/

static const FE_SPACE *fe_space;      /* initialized by main() */
static DOF_REAL_VEC   *u_h = NULL;    /* initialized by main() */
static DOF_REAL_VEC   *u_old = NULL;  /* initialized by main() */
static DOF_REAL_VEC   *f_h = NULL;    /* initialized by main() */
static DOF_MATRIX     *matrix = NULL; /* initialized by main() */
static ADAPT_INSTAT   *adapt_instat;  /* initialized by main() */
static BNDRY_FLAGS    dirichlet_mask; /* initialized by main() */

static REAL theta = 0.5;   /* parameter of the time discretization */
static REAL err_L2  = 0.0; /* spatial error in a single time step */


/*******************************************************************************
 * struct heat_leaf_data: structure for storing one REAL value on each
 *                          leaf element as LEAF_DATA
 * rw_el_est():  return a pointer to the memory for storing the element
 *               estimate (stored as LEAF_DATA), called by heat_est()
 * get_el_est(): return the value of the element estimates (from LEAF_DATA),
 *               called by adapt_method_stat() and graphics()
 ******************************************************************************/

struct heat_leaf_data
{
  REAL estimate;            /*  one real for the estimate                   */
  REAL est_c;               /*  one real for the coarsening estimate        */
};

static REAL *rw_el_est(EL *el)
{
  if (IS_LEAF_EL(el))
    return &((struct heat_leaf_data *)LEAF_DATA(el))->estimate;
  else
    return NULL;
}

static REAL get_el_est(EL *el)
{
  if (IS_LEAF_EL(el))
    return ((struct heat_leaf_data *)LEAF_DATA(el))->estimate;
  else
    return 0.0;
}

static REAL *rw_el_estc(EL *el)
{
  if (IS_LEAF_EL(el))
    return &((struct heat_leaf_data *)LEAF_DATA(el))->est_c;
  else
    return NULL;
}

static REAL get_el_estc(EL *el)
{
  if (IS_LEAF_EL(el))
    return ((struct heat_leaf_data *)LEAF_DATA(el))->est_c;
  else
    return 0.0;
}

static REAL time_est = 0.0;

static REAL get_time_est(MESH *mesh, ADAPT_INSTAT *adapt)
{
  return time_est;
}

/*******************************************************************************
 * For test purposes: exact solution and its gradient (optional)
 ******************************************************************************/

static REAL eval_time_u = 0.0;
static REAL u(const REAL_D x)
{
  return sin(M_PI*eval_time_u)*exp(-10.0*SCP_DOW(x,x));
}

static REAL eval_time_u0 = 0.0;
static REAL u0(const REAL_D x)
{
  eval_time_u = eval_time_u0;
  return u(x);
}

/*******************************************************************************
 * problem data: right hand side, boundary values
 ******************************************************************************/

static REAL eval_time_g = 0.0;
static REAL g(const REAL_D x)  /* boundary values, not optional */
{
  eval_time_u = eval_time_g;
  return u(x);
}

static REAL eval_time_f = 0.0;
static REAL f(const REAL_D x)              /* -Delta u, not optional        */
{
  REAL  r2 = SCP_DOW(x,x), ux  = sin(M_PI*eval_time_f)*exp(-10.0*r2);
  REAL  ut = M_PI*cos(M_PI*eval_time_f)*exp(-10.0*r2);
  return ut - (400.0*r2 - 20.0*DIM_OF_WORLD)*ux;
}


/*******************************************************************************
 * write error and estimator data to files
 ******************************************************************************/

static void write_statistics(const char *path, ADAPT_INSTAT *adapt, int n_dof,
			     REAL space_est, REAL time_est, REAL err_L2)
{
  static FILE *file_ndof = NULL, *file_tau = NULL;
  static FILE *file_space_est = NULL, *file_time_est = NULL, *file_L2_err = NULL;
  const char  *name = fe_space->bas_fcts->name;
  REAL        time = adapt->time;
  char        filename[1024];

  if (!file_ndof)
  {
    sprintf(filename, "%s/n_dof-%s.agr", path ? path : ".", name);
    file_ndof = fopen(filename, "w");
  }

  if (!file_tau)
  {
    sprintf(filename, "%s/tau-%s.agr", path ? path : ".", name);
    file_tau       = fopen(filename, "w");
  }

  if (!file_space_est)
  {
    sprintf(filename, "%s/space_est-%s.agr", path ? path : ".", name);
    file_space_est = fopen(filename, "w");
  }

  if (!file_time_est)
  {
    sprintf(filename, "%s/time_est-%s.agr", path ? path : ".", name);
    file_time_est  = fopen(filename, "w");
  }

  if (!file_L2_err)
  {
    sprintf(filename, "%s/L2_err-%s.agr", path ? path : ".", name);
    file_L2_err    = fopen(filename, "w");
  }

  if (file_ndof)
    fprintf(file_ndof, "%.6le %d\n", time, n_dof);

/*---  don't print zeros, zeros do not allow log display of estimate ---*/
  if (file_space_est)
    fprintf(file_space_est, "%.6le %.6le\n", time, MAX(space_est,1.e-20));

  if (time > adapt->start_time)
  {
    if (file_tau)
    {
      fprintf(file_tau, "%.6le %.6le\n", time, adapt->timestep);
    }
/*---  don't print zeros, zeros do not allow log display of estimate ---*/
    if (file_time_est)
      fprintf(file_time_est, "%.6le %.6le\n", time, MAX(time_est,1.e-20));
  }

/*---  don't print zeros, zeros do not allow log display of error    ---*/
  if (file_L2_err)
    fprintf(file_L2_err, "%.6le %.6le\n", time, MAX(err_L2,1.e-20));

  if (time >= adapt->end_time)
  {
    if (file_ndof)      fclose(file_ndof);
    if (file_tau)       fclose(file_tau);
    if (file_space_est) fclose(file_space_est);
    if (file_time_est)  fclose(file_time_est);
    if (file_L2_err)    fclose(file_L2_err);
  }
  else
  {
    fflush(NULL);
  }

  return;
}

/*******************************************************************************
 * interpolation is solve on the initial grid
 ******************************************************************************/

static void interpol_u0(MESH *mesh)
{
  dof_compress(mesh);
  interpol(u0, u_h);

  return;
}

static void init_timestep(MESH *mesh, ADAPT_INSTAT *adapt)
{
  FUNCNAME("init_timestep");

  INFO(adapt_instat->info,1,
    "---8<---------------------------------------------------\n");
  INFO(adapt_instat->info, 1,"starting new timestep\n");

  dof_copy(u_h, u_old);
  return;
}

static void set_time(MESH *mesh, ADAPT_INSTAT *adapt)
{
  FUNCNAME("set_time");

  INFO(adapt->info,1,
    "---8<---------------------------------------------------\n");
  if (adapt->time == adapt->start_time)
  {
    INFO(adapt->info, 1,"start time: %.4le\n", adapt->time);
  }
  else
  {
    INFO(adapt->info, 1,"timestep for (%.4le %.4le), tau = %.4le\n",
			 adapt->time-adapt->timestep, adapt->time,
			 adapt->timestep);
  }

  eval_time_f = adapt->time - (1 - theta)*adapt->timestep;
  eval_time_g = adapt->time;

  return;
}

static void close_timestep(MESH *mesh, ADAPT_INSTAT *adapt)
{
  FUNCNAME("close_timestep");
  static REAL err_max = 0.0;                     /* max space-time error    */
  static REAL est_max = 0.0;                     /* max space-time estimate */
  static int  write_fe_data = 0, write_stat_data = 0;
  static int  step = 0;
  static char path[256] = "./";

  REAL        space_est = adapt->adapt_space->err_sum;
  REAL        tolerance = adapt->rel_time_error*adapt->tolerance;

  err_max = MAX(err_max, err_L2);
  est_max = MAX(est_max, space_est + time_est);

  INFO(adapt->info,1,
    "---8<---------------------------------------------------\n");

  if (adapt->time == adapt->start_time) {
    tolerance = adapt->adapt_initial->tolerance;
    INFO(adapt->info,1,"start time: %.4le\n", adapt->time);
  } else {
    tolerance += adapt->adapt_space->tolerance;
    INFO(adapt->info,1,"timestep for (%.4le %.4le), tau = %.4le\n",
			adapt->time-adapt->timestep, adapt->time,
			adapt->timestep);
  }
  INFO(adapt->info,2,"max. est.  = %.4le, tolerance = %.4le\n",
		      est_max, tolerance);
  INFO(adapt->info,2,"max. error = %.4le, ratio = %.2lf\n",
		      err_max, err_max/MAX(est_max,1.0e-20));

  if (!step) {
    GET_PARAMETER(1, "write finite element data", "%d", &write_fe_data);
    GET_PARAMETER(1, "write statistical data", "%d", &write_stat_data);
    GET_PARAMETER(1, "data path", "%s", path);
  }

  /*****************************************************************************
   * write mesh and discrete solution to file for post-processing
   ****************************************************************************/

  if (write_fe_data) {
    const char *fn;

    fn= generate_filename(path, "mesh", step);
    write_mesh_xdr(mesh, fn, adapt->time);
    fn= generate_filename(path, "u_h", step);
    write_dof_real_vec_xdr(u_h, fn);
  }

  step++;

  /*****************************************************************************
   * write data about estimate, error, time step size, etc.
   ****************************************************************************/

  if (write_stat_data) {
    int n_dof = fe_space->admin->size_used;
    write_statistics(path, adapt, n_dof, space_est, time_est, err_L2);
  }

  if (do_graphics) {
    graphics(mesh, u_h, get_el_est, u, adapt->time);
  }

  return;
}


/*--------------------------------------------------------------------------*/
/* build(): assemblage of the linear system: matrix, load vector,           */
/*          boundary values, called by adapt_method_stat()                  */
/*          on the first call initialize u_h, f_h, matrix and information   */
/*          for assembling the system matrix                                */
/*                                                                          */
/* struct op_data: structure for passing information from init_element() to*/
/*                 LALt()                                                   */
/* init_element(): initialization on the element; calculates the            */
/*                 coordinates and |det DF_S| used by LALt; passes these    */
/*                 values to LALt via user_data,                            */
/*                 called on each element by update_matrix()                */
/* LALt():         implementation of -Lambda id Lambda^t for -Delta u,      */
/*                 called by update_matrix()                                */
/* c():            implementation of 1/tau*m(,)                             */
/*--------------------------------------------------------------------------*/

struct op_data /* application data (AKA "user_data") */
{
  REAL_BD Lambda; /*  the gradient of the barycentric coordinates */
  REAL    det;    /*  |det D F_S|                                 */

  REAL    tau_1;
};

static const REAL_B *LALt(const EL_INFO *el_info, const QUAD *quad,
			  int iq, void *ud)
{
  struct op_data *info = (struct op_data *)ud;
  int i, j, k, dim = el_info->mesh->dim;
  static REAL_BB LALt;

  for (i = 0; i <= dim; i++)
    for (j = i; j <= dim; j++)
    {
      for (LALt[i][j] = k = 0; k < DIM_OF_WORLD; k++)
	LALt[i][j] += info->Lambda[i][k]*info->Lambda[j][k];
      LALt[i][j] *= info->det;
      LALt[j][i] = LALt[i][j];
    }
  return (const REAL_B *)LALt;
}

static REAL c(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
{
  struct op_data *info = (struct op_data *)ud;

  return info->tau_1*info->det;
}

static void assemble(DOF_MATRIX *matrix, DOF_REAL_VEC *fh, DOF_REAL_VEC *uh,
		     const DOF_REAL_VEC *u_old,
		     REAL theta, REAL tau,
		     REAL (*f)(const REAL_D), REAL (*g)(const REAL_D),
		     const BNDRY_FLAGS dirichlet_mask)
{
  /* Some quantities remembered across calls. Think of this routine
   * like being a "library function" ... The stuff is re-initialized
   * whenever the finite element space changes. We use fe_space->admin
   * to check for changes in the finite element space because
   * DOF_ADMIN's are persisitent within ALBERTA, while fe-space are
   * not.
   */
  static EL_MATRIX_INFO stiff_elmi, mass_elmi;
  static const DOF_ADMIN *admin = NULL;
  static const QUAD *quad = NULL;
  static struct op_data op_data[1]; /* storage for det and Lambda */

  /* Remaining (non-static) variables. */
  const BAS_FCTS *bas_fcts;
  FLAGS          fill_flag;
  EL_SCHAR_VEC   *bound;

  /* Initialize persistent variables. */
  if (admin != uh->fe_space->admin) {
    OPERATOR_INFO stiff_opi = { NULL, }, mass_opi = { NULL, };

    admin    = uh->fe_space->admin;

    stiff_opi.row_fe_space   = uh->fe_space;
    stiff_opi.quad[2]        = quad;
    stiff_opi.LALt.real      = LALt;
    stiff_opi.LALt_pw_const  = true;
    stiff_opi.LALt_symmetric = true;
    stiff_opi.user_data      = op_data;

    fill_matrix_info(&stiff_opi, &stiff_elmi);

    mass_opi.row_fe_space = uh->fe_space;
    mass_opi.quad[0]      = quad;
    mass_opi.c.real       = c;
    mass_opi.c_pw_const   = true;
    mass_opi.user_data    = op_data;

    fill_matrix_info(&mass_opi, &mass_elmi);

    quad = get_quadrature(uh->fe_space->bas_fcts->dim,
			  2*uh->fe_space->bas_fcts->degree);
  }

  op_data->tau_1 = 1.0/tau;

  /* Assemble the matrix and the right hand side. The idea is to
   * assemble the local mass and stiffness matrices only once, and to
   * use it to update both, the system matrix and the contribution of
   * the time discretisation to the RHS.
   */
  clear_dof_matrix(matrix);
  dof_set(0.0, fh);

  bas_fcts = uh->fe_space->bas_fcts;

  bound = get_el_schar_vec(bas_fcts);

  fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_BOUND;
  TRAVERSE_FIRST(uh->fe_space->mesh, -1, fill_flag) {
    const EL_REAL_VEC  *u_old_loc;
    const EL_DOF_VEC   *dof;
    const EL_BNDRY_VEC *bndry_bits;
    const EL_MATRIX *stiff_loc, *mass_loc;

    /* Get the local coefficients of u_old, boundary info, dof-mapping */
    u_old_loc  = fill_el_real_vec(NULL, el_info->el, u_old);
    dof        = get_dof_indices(NULL, uh->fe_space, el_info->el);
    bndry_bits = get_bound(NULL, bas_fcts, el_info);

    /* Initialization of values used by LALt and c. It is not
     * necessary to introduce an extra "init_element()" hook for our
     * OPERATOR_INFO structures; the line below is just what would be
     * contained in that function (compare with ellipt.c).
     *
     * Beware to replace the "..._0cd()" for co-dimension 0 by its
     * proper ..._dim() variant if ever "porting" this stuff to
     * parametric meshes on manifolds.
     */
    op_data->det = el_grd_lambda_0cd(el_info, op_data->Lambda);

    /* Obtain the local (i.e. per-element) matrices. */
    stiff_loc = stiff_elmi.el_matrix_fct(el_info, stiff_elmi.fill_info);
    mass_loc  = mass_elmi.el_matrix_fct(el_info, mass_elmi.fill_info);

    /* Translate the geometric boundary classification into
     * Dirichlet/Neumann/Interior boundary condition
     * interpretation. Inside the loop over the mesh-elements we need
     * only to care about Dirichlet boundary conditions.
     */
    dirichlet_map(bound, bndry_bits, dirichlet_mask);

    /* add theta*a(psi_i,psi_j) + 1/tau*m(4*u^3*psi_i,psi_j) */
    if (theta) {
      add_element_matrix(matrix,
			 theta, stiff_loc, NoTranspose, dof, dof, bound);
    }
    add_element_matrix(matrix, 1.0, mass_loc, NoTranspose, dof, dof, bound);

      /* compute the contributions from the old time-step:
       *
       * f += -(1-theta)*a(u_old,psi_i) + 1/tau*m(u_old,psi_i)
       */
    bi_mat_el_vec(-(1.0 - theta), stiff_loc,
		  1.0, mass_loc, u_old_loc,
		  1.0, fh, dof, bound);

  } TRAVERSE_NEXT();

  free_el_schar_vec(bound);

  /* Indicate that the boundary conditions are built into the matrix,
   * needed e.g. by the hierarchical preconditioners.
   */
  BNDRY_FLAGS_CPY(matrix->dirichlet_bndry, dirichlet_mask);

  /* Add the "force-term" to the right hand side (L2scp_...() is additive) */
  L2scp_fct_bas(f, quad, fh);

  /* Close the system by imposing suitable boundary conditions. Have a
   * look at ellipt.c for how to impose more complicated stuff; here
   * we only use Dirichlet b.c.
   */
  dirichlet_bound(fh, uh, NULL, dirichlet_mask, g);
}

static void build(MESH *mesh, U_CHAR flag)
{
  FUNCNAME("build");

  dof_compress(mesh);

  INFO(adapt_instat->adapt_space->info, 2,
    "%d DOFs for %s\n", fe_space->admin->size_used, fe_space->name);

  assemble(matrix, f_h, u_h, u_old, theta, adapt_instat->timestep,
	   f, g, dirichlet_mask);
}

/*--------------------------------------------------------------------------*/
/* solve(): solve the linear system, called by adapt_method_stat()          */
/*--------------------------------------------------------------------------*/

static void solve(MESH *mesh)
{
  FUNCNAME("solve");
  static REAL tol = 1.e-8, ssor_omega = 1.0;
  static int  max_iter = 1000, info = 2, restart = 0, ssor_iter = 1, ilu_k = 8;
  static OEM_PRECON icon = DiagPrecon;
  static OEM_SOLVER solver = NoSolver;
  const PRECON *precon;

  if (solver == NoSolver) {
    GET_PARAMETER(1, "solver", "%d", &solver);
    GET_PARAMETER(1, "solver tolerance", "%f", &tol);
    GET_PARAMETER(1, "solver precon", "%d", &icon);
    GET_PARAMETER(1, "solver max iteration", "%d", &max_iter);
    GET_PARAMETER(1, "solver info", "%d", &info);
    if (icon == __SSORPrecon) {
    	GET_PARAMETER(1, "precon ssor omega", "%f", &ssor_omega);
    	GET_PARAMETER(1, "precon ssor iter", "%d", &ssor_iter);
    }
    if (icon == ILUkPrecon)
        GET_PARAMETER(1, "precon ilu(k)", "%d", &ilu_k);
    if (solver == GMRes) {
      GET_PARAMETER(1, "solver restart", "%d", &restart);
    }
  }

  if (icon == ILUkPrecon)
	  precon = init_oem_precon(matrix, NULL, info, ILUkPrecon, ilu_k);
  else
	  precon = init_oem_precon(matrix, NULL, info, icon, ssor_omega, ssor_iter);
  oem_solve_s(matrix, NULL, f_h, u_h,
	      solver, tol, precon, restart, max_iter, info);
}

/*--------------------------------------------------------------------------*/
/* Functions for error estimate:                                            */
/* estimate():   calculates error estimate via heat_est()                   */
/*               calculates exact error also (only for test purpose),       */
/*               called by adapt_method_stat()                              */
/*--------------------------------------------------------------------------*/

static REAL r(const EL_INFO *el_info,
	      const QUAD *quad, int iq,
	      REAL uh_at_qp, const REAL_D grd_uh_at_qp,
	      REAL t)
{
  REAL_D      x;

  coord_to_world(el_info, quad->lambda[iq], x);
  eval_time_f = t;

  return -f(x);
}


static REAL estimate(MESH *mesh, ADAPT_STAT *adapt)
{
  FUNCNAME("estimate");
  static REAL    C[4] = {-1.0, 1.0, 1.0, 1.0};
  static REAL_DD A = {{0.0}};
  FLAGS          r_flag = 0;  /* = (INIT_UH|INIT_GRD_UH), if needed by r()  */
  int            n;
  REAL           space_est;

  eval_time_u = adapt_instat->time;

  if (C[0] < 0) {
    C[0] = 1.0;
    GET_PARAMETER(1, "estimator C0", "%f", &C[0]);
    GET_PARAMETER(1, "estimator C1", "%f", &C[1]);
    GET_PARAMETER(1, "estimator C2", "%f", &C[2]);
    GET_PARAMETER(1, "estimator C3", "%f", &C[3]);

    for (n = 0; n < DIM_OF_WORLD; n++) {
      A[n][n] = 1.0;   /* set diogonal of A; all other elements are zero */
    }
  }

  time_est = heat_est(u_h, u_old, adapt_instat, rw_el_est, rw_el_estc,
		      -1 /* quad_degree */,
		      C, (const REAL_D *)A, dirichlet_mask,
		      r, r_flag, NULL /* gn() */, 0 /* gn_flag */);

  space_est = adapt_instat->adapt_space->err_sum;
  err_L2 = L2_err(u, u_h, NULL, false, false, NULL, NULL);

  INFO(adapt_instat->info,2,
    "---8<---------------------------------------------------\n");
  INFO(adapt_instat->info, 2,"time = %.4le with timestep = %.4le\n",
			      adapt_instat->time, adapt_instat->timestep);
  INFO(adapt_instat->info, 2,"estimate   = %.4le, max = %.4le\n", space_est,
			      sqrt(adapt_instat->adapt_space->err_max));
  INFO(adapt_instat->info, 2,"||u-uh||L2 = %.4le, ratio = %.2lf\n", err_L2,
			      err_L2/MAX(space_est,1.e-20));

  return adapt_instat->adapt_space->err_sum;
}

static REAL est_initial(MESH *mesh, ADAPT_STAT *adapt)
{
  err_L2 = adapt->err_sum =
    L2_err(u0, u_h, NULL, false, false, rw_el_est, &adapt->err_max);
  return adapt->err_sum;
}

/*--------------------------------------------------------------------------*/
/* main program                                                             */
/*--------------------------------------------------------------------------*/

int main(int argc, char **argv)
{
  FUNCNAME("main");
  MACRO_DATA     *data;
  MESH           *mesh;
  const BAS_FCTS *lagrange;
  int             n_refine = 0, p = 1, dim;
  char            filename[PATH_MAX];
  REAL            fac = 1.0;

  /*****************************************************************************
   * first of all, initialize the access to parameters of the init file
   ****************************************************************************/
  parse_parameters(argc, argv, "INIT/heat.dat");

  GET_PARAMETER(1, "mesh dimension", "%d", &dim);
  GET_PARAMETER(1, "macro file name", "%s", filename);
  GET_PARAMETER(1, "global refinements", "%d", &n_refine);
  GET_PARAMETER(1, "parameter theta", "%e", &theta);
  GET_PARAMETER(1, "polynomial degree", "%d", &p);
  GET_PARAMETER(1, "online graphics", "%d", &do_graphics);

  /*****************************************************************************
  *  get a mesh, and read the macro triangulation from file
  ****************************************************************************/
  data = read_macro(filename);
  mesh = GET_MESH(dim, "ALBERTA mesh", data,
		  NULL /* init_node_projection() */,
		  NULL /* init_wall_trafos() */);
  free_macro_data(data);

  init_leaf_data(mesh, sizeof(struct heat_leaf_data),
		 NULL /* refine_leaf_data() */,
		 NULL /* coarsen_leaf_data() */);

  /* Finite element spaces can be added at any time, but it is more
   * efficient to do so before refining the mesh a lot.
   */
  lagrange = get_lagrange(mesh->dim, p);
  TEST_EXIT(lagrange, "no lagrange BAS_FCTS\n");

  fe_space = get_fe_space(mesh, lagrange->name, lagrange, 1, ADM_FLAGS_DFLT);

  global_refine(mesh, n_refine * dim, FILL_NOTHING);

  /*****************************************************************************
   *  initialize the global variables shared across build(), solve()
   *  and estimate().
   ****************************************************************************/
  matrix = get_dof_matrix("A", fe_space, NULL);
  f_h    = get_dof_real_vec("f_h", fe_space);
  u_h    = get_dof_real_vec("u_h", fe_space);
  u_h->refine_interpol = fe_space->bas_fcts->real_refine_inter;
  u_h->coarse_restrict = fe_space->bas_fcts->real_coarse_inter;
  u_old  = get_dof_real_vec("u_old", fe_space);
  u_old->refine_interpol = fe_space->bas_fcts->real_refine_inter;
  u_old->coarse_restrict = fe_space->bas_fcts->real_coarse_inter;
  dof_set(0.0, u_h);      /*  initialize u_h  !                          */

  BNDRY_FLAGS_ALL(dirichlet_mask); /* Only Dirichlet b.c. supported here */

  /*****************************************************************************
   *  init adapt_instat structure
   ****************************************************************************/
  adapt_instat = get_adapt_instat(dim, "heat", "adapt", 2, adapt_instat);

  /* Some animation in between ... */
  if (do_graphics) {
    graphics(mesh, NULL, NULL, NULL, adapt_instat->start_time);
  }

  /*****************************************************************************
   *  adapt time step size to refinement level and polynomial degree,
   *  based on the known L2-error error estimates.
   ****************************************************************************/
  if (theta < 0.5) {
    WARNING("You are using the explicit Euler scheme\n");
    WARNING("Use a sufficiently small time step size!!!\n");
    fac = 1.0e-3;
  }

  if (theta == 0.5) {
    adapt_instat->timestep *= fac*pow(2, -(REAL)(p*(n_refine))/2.0);
  } else {
    adapt_instat->timestep *= fac*pow(2, -(REAL)(p*(n_refine)));
  }
  MSG("using initial timestep size = %.4le\n", adapt_instat->timestep);

  eval_time_u0 = adapt_instat->start_time;

  adapt_instat->adapt_initial->get_el_est = get_el_est;
  adapt_instat->adapt_initial->estimate = est_initial;
  adapt_instat->adapt_initial->solve = interpol_u0;

  adapt_instat->adapt_space->get_el_est   = get_el_est;
  adapt_instat->adapt_space->get_el_estc  = get_el_estc;
  adapt_instat->adapt_space->estimate = estimate;
  adapt_instat->adapt_space->build_after_coarsen = build;
  adapt_instat->adapt_space->solve = solve;

  adapt_instat->init_timestep  = init_timestep;
  adapt_instat->set_time       = set_time;
  adapt_instat->get_time_est   = get_time_est;
  adapt_instat->close_timestep = close_timestep;

  /*****************************************************************************
   * ... off we go ...
   ****************************************************************************/
  adapt_method_instat(mesh, adapt_instat);

  WAIT_REALLY;

  return 0;
}