xref: /dports/math/alberta/alberta3-920315ae1bbfd1b1fb6672d916619ac37a411e95/demo/src/Common/ellipt-dg.c

/*******************************************************************************
 *
 * ALBERTA:  an Adaptive multi Level finite element toolbox  using
 *           Bisectioning refinement and Error control by Residual
 *           Techniques
 *
 * www.alberta-fem.de
 *
 *
 * file:     ellipt-dg.c
 *
 * description:  solver for an elliptic model problem
 *
 *                        -\Delta u = f  in \Omega
 *                          D_\nu u = g  on \partial \Omega
 *
 * This version of the Poisson solver is just a quick and dirty hack
 * to test the basic functionality of the crude DG support built into
 * ALBERTA. We use the following scheme:
 *
 * (\nabla u, \nabla v) + \alpha \int_E [u][v] = (f, v) + b.c. in \Omega.
 *
 * alpha is chosen to be 1/h^q for some q. The discretisation is of
 * course not state of the art, however, this is really only meant as
 * test program.
 *
 ******************************************************************************
 *
 * author(s): Claus-Justus Heine
 *            Abteilung fuer Angewandte Mathematik
 *            Albert-Ludwigs-Universitaet Freiburg
 *            Hermann-Herder-Str. 10
 *            79104 Freiburg
 *            Germany
 *            Claus.Heine@Mathematik.Uni-Freiburg.DE
 *
 * (c) by C.-J. Heine (2008)
 *
 * Based on the famous ellipt.c demo-program by K.G. Siebert and A. Schmidt.
 *
 ******************************************************************************/

#include "alberta-demo.h" /* proto-types for some helper functions */

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

/*******************************************************************************
 * global variables: finite element space, discrete solution
 *                   load vector and system matrix
 *
 * 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;       /* initialized by main() */
static DOF_REAL_VEC   *f_h;       /* initialized by main() */
static DOF_MATRIX     *matrix;    /* initialized by main() */

static REAL robin_alpha = -1.0;    /* GET_PARAMTER() in main() */
static REAL jmp_alpha   = 1.0;     /* penalty for jump over walls */
static REAL jmp_exp     = 1.0;     /* jmp_alpha / pow(h, jmp_exp) */

/*******************************************************************************
 * struct ellipt_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 ellipt_est()
 * get_el_est(): return the value of the element estimates (from LEAF_DATA),
 *               called by adapt_method_stat() and graphics()
 ******************************************************************************/

struct ellipt_leaf_data
{
  REAL estimate; /*  one real for the estimate */
};

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

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

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

#define GAUSS_SCALE 10.0

static REAL u(const REAL_D x)
{
  return exp(-GAUSS_SCALE*SCP_DOW(x,x));
}

static const REAL *grd_u(const REAL_D x, REAL_D grd)
{
  static REAL_D buffer;
  REAL          ux = exp(-GAUSS_SCALE*SCP_DOW(x,x));
  int           n;

  if (!grd) {
    grd = buffer;
  }

  for (n = 0;  n < DIM_OF_WORLD; n++)
    grd[n] = -2.0*GAUSS_SCALE*x[n]*ux;

  return grd;
}

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

#if 0
static REAL g(const REAL_D x) /* boundary values, not optional */
{
  return u(x);
}
#endif

static REAL gn(const REAL_D x, const REAL_D normal) /* Neumann b.c. */
{
  return robin_alpha > 0.0
    ? SCP_DOW(grd_u(x, NULL), normal) + robin_alpha * u(x)
    : SCP_DOW(grd_u(x, NULL), normal);
}

static REAL f(const REAL_D x) /* -Delta u, not optional        */
{
  REAL  r2 = SCP_DOW(x,x), ux  = exp(-GAUSS_SCALE*r2);
  return -(4.0*SQR(GAUSS_SCALE)*r2 - 2.0*GAUSS_SCALE*DIM_OF_WORLD)*ux;
}

/*******************************************************************************
 * 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() ("operator data").
 * 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() after init_element()
 ******************************************************************************/

struct op_data
{
  REAL_BD Lambda; /*  the gradient of the barycentric coordinates */
  REAL    det;    /*  |det D F_S| */
  REAL    wall_det[N_WALLS_MAX];
  REAL    c_factor;
};

static bool init_element(const EL_INFO *el_info, const QUAD *quad[3], void *ud)
{
  struct op_data *info = (struct op_data *)ud;

  /* ..._0cd: co-dimension 0 version of el_grd_lambda(dim, ...) */
  info->det = el_grd_lambda_0cd(el_info, info->Lambda);

  return false; /* not parametric */
}

static bool bndry_init_element(const EL_INFO *el_info, int wall,
			       const WALL_QUAD *quad[3], void *ud)
{
  struct op_data *info = (struct op_data *)ud;

  info->wall_det[wall] = get_wall_normal(el_info, wall, NULL);

  if (fe_space->bas_fcts->degree > 0) {
    REAL det = info->wall_det[wall];
    info->wall_det[wall] = info->c_factor =
      (det *= jmp_alpha * pow(det, -(jmp_exp / (REAL)(DIM_OF_WORLD-1))));
  } else {
    /* Use a very simplistic finite-volume method */
#define VAL (1.0/(REAL)N_LAMBDA_MAX)
    const REAL_B c_b = INIT_BARY_MAX(VAL, VAL, VAL, VAL);
#undef VAL
    EL_INFO neigh_info[1];
    const EL_GEOM_CACHE *elgc =
      fill_el_geom_cache(el_info, FILL_EL_WALL_REL_ORIENTATION(wall));
    REAL_D c, cn;
    REAL delta;

    fill_neigh_el_info(neigh_info, el_info, wall, elgc->rel_orientation[wall]);

    coord_to_world(el_info, c_b, c);
    coord_to_world(neigh_info, c_b, cn);
    delta = DIST_DOW(c, cn);

    info->c_factor = (info->wall_det[wall] /= delta);
  }

  return false;
}

static bool neigh_init_element(const EL_INFO *el_info, int wall,
			       const WALL_QUAD *quad[3], void *ud)
{
  struct op_data *info = (struct op_data *)ud;

  info->c_factor = -info->wall_det[wall];

  return false;
}

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

  for (i = 0; i < N_VERTICES(dim); i++) {
    LALt[i][i] = info->det*SCP_DOW(info->Lambda[i], info->Lambda[i]);
    for (j = i+1; j < N_VERTICES(dim); j++) {
      LALt[i][j] = SCP_DOW(info->Lambda[i], info->Lambda[j]);
      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->c_factor;
}

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

  dof_compress(mesh);
  MSG("%d DOFs for %s\n", fe_space->admin->size_used, fe_space->name);

  if (!matrix_info) {
    /* information for matrix assembling (only once) */
    OPERATOR_INFO       o_info = { NULL, }, *oi_ptr;
    BNDRY_OPERATOR_INFO el_bop_info = { NULL, };
    BNDRY_OPERATOR_INFO neigh_bop_info = { NULL, };
    static struct op_data user_data; /* storage for det and Lambda,
				      * must be static.
				      */

    o_info.row_fe_space    = fe_space;
    o_info.init_element    = init_element;
    o_info.LALt.real       = LALt;
    o_info.LALt_pw_const   = true;        /* pw const. assemblage is faster */
    o_info.LALt_symmetric  = true;        /* symmetric assemblage is faster */
    o_info.user_data = (void *)&user_data; /* application data */

    /* only FILL_BOUND is added automatically. */
    o_info.fill_flag = CALL_LEAF_EL|FILL_COORDS;

    /* ALBERTA requires that the boundary contributions are split into
     * one part which is assembled using only the current's element
     * local basis functions, and a second part which pairs the given
     * element with its neighbour. In our case both contributions are
     * simple mass-matrices
     */
    el_bop_info.row_fe_space = fe_space;
    el_bop_info.init_element = bndry_init_element;
    el_bop_info.c.real       = c;
    el_bop_info.c_pw_const   = true;
    el_bop_info.user_data    = (void *)&user_data;
    BNDRY_FLAGS_INIT(el_bop_info.bndry_type); /* all interior edges */

    el_bop_info.fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_OPP_COORDS;

    /* The contribution of the pairing with the neighbouring basis
     * functions is almost the same:
     */
    neigh_bop_info = el_bop_info;
    neigh_bop_info.init_element = neigh_init_element;
    neigh_bop_info.discontinuous = true;

    oi_ptr = fe_space->bas_fcts->degree > 0 ? &o_info : NULL;
    matrix_info =
      fill_matrix_info_ext(NULL, oi_ptr, &el_bop_info, &neigh_bop_info, NULL);
  }

  /* assembling of matrix */
  clear_dof_matrix(matrix);
  update_matrix(matrix, matrix_info, NoTranspose);

  /* assembling of load vector */
  dof_set(0.0, f_h);
  L2scp_fct_bas(f, get_quadrature(DIM_OF_WORLD, 4) /* quadrature */, f_h);

  /* Boundary values, the combination alpha_r < 0.0 flags automatic
   * mean-value correction iff f_h has non-zero mean-value and no
   * non-Neumann boundary conditions were detected during mesh
   * traversal.
   */
  boundary_conditions(matrix, f_h, u_h, NULL /* bound */,
		      NULL /* Dirichlet_mask */,
		      NULL /* g */, gn,
		      robin_alpha, /* < 0: automatic mean-value correction */
		      NULL /* wall_quad, use default */);
}

/*******************************************************************************
 * 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        miter = 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)
  {
    tol = 1.e-8;
    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", &miter);
    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, miter, info);

  if (do_graphics) {
    MSG("Displaying u_h, u and (u_h-u).\n");
    graphics(mesh, u_h, NULL /* get_el_est */, u, HUGE_VAL /* time */);
  }

  return;
}

/*******************************************************************************
 * Functions for error estimate:
 * estimate():   calculates error estimate via ellipt_est()
 *               calculates exact error also (only for test purpose),
 *               called by adapt_method_stat()
 * r():          calculates the lower order terms of the element residual
 *               on each element at the quadrature node iq of quad
 *               argument to ellipt_est() and called by ellipt_est()
 ******************************************************************************/

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_D x;

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

  return -f(x);
}

static REAL est_gn(const EL_INFO *el_info,
		   const QUAD *quad,
		   int qp,
		   REAL uh_at_qp,
		   const REAL_D normal)
{
  /* we simply return gn(), expoiting the fact that the geometry cache
   * of the quadrature already contains the world-coordinates of the
   * quadrature points.
   */
  const QUAD_EL_CACHE *qelc =
    fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);

  if (robin_alpha > 0.0) {
    return gn(qelc->world[qp], normal) - robin_alpha * uh_at_qp;
  } else {
    return gn(qelc->world[qp], normal);
  }
}

#define EOC(e,eo) log(eo/MAX(e,1.0e-15))/M_LN2

static REAL estimate(MESH *mesh, ADAPT_STAT *adapt)
{
  FUNCNAME("estimate");
  static int   norm = -1;
  static REAL  C[3] = {1.0, 1.0, 0.0};
  static REAL  est_old = -1.0, err_old = -1.0;
  REAL         est, err;
  REAL_DD      A = {{0.0}};
  int          n;

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

  if (norm < 0) {
    norm = H1_NORM;
    GET_PARAMETER(1, "error norm", "%d", &norm);
    GET_PARAMETER(1, "estimator C0", "%f", &C[0]);
    GET_PARAMETER(1, "estimator C1", "%f", &C[1]);
    GET_PARAMETER(1, "estimator C2", "%f", &C[2]);
  }

  est = ellipt_est(u_h, adapt, rw_el_est, NULL /* rw_est_c() */,
		   -1 /* quad_degree */,
		   norm, C,
		   (const REAL_D *) A,
		   NULL,
		   r, 0 /* (INIT_UH | INIT_GRD_UH), if needed by r() */,
		   est_gn, robin_alpha > 0.0 ? INIT_UH : 0);

  MSG("estimate   = %.8le", est);
  if (est_old >= 0)
    print_msg(", EOC: %.2lf\n", EOC(est,est_old));
  else
    print_msg("\n");
  est_old = est;

  if (norm == L2_NORM) {
    err = L2_err(u, u_h, NULL /* quad */,
		 false /* relative error*/,
		 true /* mean-value adjust */,
		 NULL /* rw_err_el()*/, NULL /* max_err_el2 */);
  } else {
    err = H1_err(grd_u, u_h, NULL /* quad */,
		 false /* relative error */,
		 NULL /* rw_err_el()*/, NULL /* max_err_el2 */);
  }

  MSG("||u-uh||%s = %.8le", norm == L2_NORM ? "L2" : "H1", err);
  if (err_old >= 0)
    print_msg(", EOC: %.2lf\n", EOC(err,err_old));
  else
    print_msg("\n");
  err_old = err;
  MSG("||u-uh||%s/estimate = %.2lf\n", norm == L2_NORM ? "L2" : "H1",
      err/MAX(est,1.e-15));

  if (do_graphics) {
    MSG("Displaying the estimate.\n");
    graphics(mesh, NULL /* u_h */, get_el_est, NULL /* u_exact() */,
	     HUGE_VAL /* time */);
  }

  return adapt->err_sum;
}

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

int main(int argc, char **argv)
{
  FUNCNAME("main");
  MACRO_DATA     *data;
  MESH           *mesh;
  int            n_refine = 0, dim, degree = 1;
  const BAS_FCTS *bas_fcts;
  ADAPT_STAT     *adapt;
  char           filename[PATH_MAX];
  char           bfcts_name[PATH_MAX] = "disc_lagrange";

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

  GET_PARAMETER(1, "mesh dimension", "%d", &dim);
  GET_PARAMETER(1, "macro file name", "%s", filename);
  GET_PARAMETER(1, "polynomial degree", "%d", &degree);
  GET_PARAMETER(1, "global refinements", "%d", &n_refine);
  GET_PARAMETER(1, "online graphics", "%d", &do_graphics);
  GET_PARAMETER(1, "robin factor", "%f", &robin_alpha);
  GET_PARAMETER(1, "jump penalty", "%f", &jmp_alpha);
  GET_PARAMETER(1, "penalty exponent", "%f", &jmp_exp);
  GET_PARAMETER(1, "basis functions", "%s", bfcts_name);

  /*****************************************************************************
   *  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 ellipt_leaf_data),
		 NULL /* refine_leaf_data() */,
		 NULL /* coarsen_leaf_data() */);

  global_refine(mesh, n_refine * mesh->dim, FILL_NOTHING);

  if (do_graphics) {
    MSG("Displaying the mesh.\n");
    graphics(mesh, NULL /* u_h */, NULL /* get_est()*/ , NULL /* u_exact() */,
	     HUGE_VAL /* time */);
  }

  /*****************************************************************************
   *  initialize the global variables shared across build(), solve()
   *  and estimate().
   ****************************************************************************/
  sprintf(bfcts_name + strlen(bfcts_name), "%d", degree);
  bas_fcts = get_bas_fcts(mesh->dim, bfcts_name);

#if 0
  {
    const QUAD *quad = get_quadrature(bas_fcts->dim, 2 * bas_fcts->degree);
    const QUAD_FAST *qfast = get_quad_fast(bas_fcts, quad, INIT_PHI);
    REAL mass_matrix[bas_fcts->n_bas_fcts][bas_fcts->n_bas_fcts];
    int i, j, iq;

    for (i = 0; i < bas_fcts->n_bas_fcts; i++) {
      for (j = 0; j < bas_fcts->n_bas_fcts; j++) {
	mass_matrix[i][j] = 0.0;
	for (iq = 0; iq < quad->n_points; iq++) {
	  mass_matrix[i][j] +=
	    qfast->w[iq] * qfast->phi[iq][i] *  qfast->phi[iq][j];
	}
	MSG("m[%d,%d] = %e\n", i, j, mass_matrix[i][j]);
      }
    }

  }
#endif

  TEST_EXIT(bas_fcts, "no BAS_FCTS of name \"%s\"\n", bfcts_name);
  fe_space =
    get_fe_space(mesh, bas_fcts->name, bas_fcts, 1 /* rdim */, ADM_FLAGS_DFLT);

  matrix = get_dof_matrix("A", fe_space, NULL /* col_fe_space */);
  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;
  dof_set(0.0, u_h);      /*  initialize u_h */

  /*****************************************************************************
   *  init adapt structure and start adaptive method
   ****************************************************************************/
  adapt = get_adapt_stat(mesh->dim, "ellipt", "adapt", 2,
			 NULL /* ADAPT_STAT storage area, optional */);
  adapt->estimate = estimate;
  adapt->get_el_est = get_el_est;
  adapt->build_after_coarsen = build;
  adapt->solve = solve;

#if 0
  {
    const QUAD *quad = get_quadrature(mesh->dim, 10);
    interpol(u, u_h);
    MSG("L2-interpolation error before global refine: %e\n",
	L2_err(u, u_h, quad,
	       false /* relative error*/,
	       false /* mean-value adjust */,
	       NULL /* rw_err_el()*/, NULL /* max_err_el2 */));
    if (do_graphics) {
      graphics(mesh, u_h, NULL, u, HUGE_VAL /* time */);
    }
    global_refine(mesh, mesh->dim);
    MSG("L2-interpolation error after global refine: %e\n",
	L2_err(u, u_h, quad,
	       false /* relative error*/,
	       false /* mean-value adjust */,
	       NULL /* rw_err_el()*/, NULL /* max_err_el2 */));
    if (do_graphics) {
      graphics(mesh, u_h, NULL, u, HUGE_VAL /* time */);
    }
    global_coarsen(mesh, -mesh->dim);
    MSG("L2-interpolation error after global coarsen: %e\n",
	L2_err(u, u_h, quad,
	       false /* relative error*/,
	       false /* mean-value adjust */,
	       NULL /* rw_err_el()*/, NULL /* max_err_el2 */));
    if (do_graphics) {
      graphics(mesh, u_h, NULL, u, HUGE_VAL /* time */);
    }
  }
#endif

  adapt_method_stat(mesh, adapt);

  if (do_graphics) {
    MSG("Displaying u_h, u, (u_h-u) and the final estimate.\n");
    graphics(mesh, u_h, get_el_est, u, HUGE_VAL /* time */);
  }
  WAIT_REALLY;

  return 0;
}