diff --git a/howto-port-ellipt.txt b/howto-port-ellipt.txt
new file mode 100644
index 0000000000000000000000000000000000000000..14ac6d72079a9d2ca04b0ee8b79cfc3e1db98259
--- /dev/null
+++ b/howto-port-ellipt.txt
@@ -0,0 +1,813 @@
+Hey, Emacs, we're -*- text -*- mode.
+
+Porting the v1.2 `ellipt.c' demo-program to post-V2.0 ALBERTA
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The source code of the changed ellipt program is attached at the end
+of this file.
+
+First, one has to compile and install the new ALBERTA version, very
+briefly:
+
+ cd alberta-XYZ/
+ ./configure --prefix=PREFIX
+ make install
+
+PREFIX is the installation PREFIX and should be an absolute pathname
+(e.g. `${HOME}/mysoftware/'). The `configure' script searches first
+below `PREFIX/lib/' and `PREFIX/include/' for support libraries, then
+in the system default paths.. For example, the easiest way to tell
+`configure' where to find the gltools-package is to copy
+libgltools.{a,so} to `PREFIX/lib/' and the various header files to
+`PREFIX/include/'. `configure --help' will give a somewhat lengthy
+list of the many other options, also the file `README' in the ALBERTA
+package contains some information to fine tune the installation, if
+necessary.
+
+After installing ALBERTA the new demo-package is located below
+
+ PREFIX/share/alberta/alberta-XYZ-demo.tar.gz
+
+The layout of the demo-package has only changed very little, also the
+structure of the Makefile's has not changed much.
+
+What follows is a descriptions of the necessary changes which have to
+be made to port the old V1.2 ellipt.c to the new post-V2.0 version of
+ALBERTA such that it compiles when copied to the new demo-package:
+
+- Some functions shared between the demo programs have been moved into a
+ small convenience library, "libdemo.a". The prototypes for the
+ functions defined in the library can be found in the header-file
+ `src/Common/alberta-demo.h'. The new makefiles in the `src/Xd/'
+ subdirectories contain targets to build this library, and the
+ demo-programs are linked against this library.
+
+- The function for displaying mesh, discrete solution, and/or estimate
+ is different:
+
+ graphics(mesh, u_h, get_el_est, u_exact(), time);
+
+ where the last two are new arguments. (Remark: the new version of
+ the demo-program also has a global variable 'bool do_graphics =
+ true/false', which determines whether or not online-graphics is
+ disabled).
+
+- In the old (i.e. V1.2) version of ALBERTA, the function
+ `init_dof_admin()' was passed to GET_MESH() in order to initialise
+ the finite element spaces needed by the application at the proper
+ stage during forming of the initial mesh. The post-1.2 versions of
+ ALBERTA allow for adding finite element spaces at any time,
+ consequently it is possible to call
+
+ lagrange = get_lagrange(...);
+ fe_space = get_fe_space(...);
+
+ in the main-function _after_ the call to GET_MESH().
+
+- Formerly, an `init_leaf_data()' had to be passed to the
+ GET_MESH()-routine. The new style is to call the library function
+
+ init_leaf_data(mesh, size, refine_leaf_data(), coarsen_leaf_data());
+
+ after initialisation the grid by calling `GET_MESH()'.
+
+- The symbolic constant `DIM' does not exist anymore, one has to use
+ `mesh->dim'. Of course, for the ellipt demo-program we have
+ DIM_OF_WORLD == mesh->dim. DIM_OF_WORLD is still fixed; there are
+ distinct ALBERTA-libraries for each supported value of DIM_OF_WORLD.
+
+- There are some new type-definitions, e.g. 'REAL_BD Lambda' it is
+ equivalent to 'REAL_D Lambda[N_LAMBDA_MAX]'.
+
+- The function `init_element()' has to be a `bool'-function, instead
+ of `void'. It is supposed to return `true' if the current element is
+ a curved element, and `false' otherwise. For the standard case one
+ should simply return `false'.
+
+- As above mentioned there exists new types for barycentric coordinates,
+ so you may want to change, e.g., the function 'LALt' like this:
+
+const REAL_B *LALt(const EL_INFO *el_info, const QUAD *quad,
+ int iq, void *ud)
+{
+ struct op_info *info = (struct op_info *)ud;
+ int i, j, dim = el_info->mesh->dim;
+ static REAL_BB LALt;
+
+ 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);
+}
+
+ Of course, nothing substantial has changed here, it just makes
+ LALt() a little bit more readable.
+
+- The structure OPERATOR_INFO has also been modified, instead of
+
+ o_info.LALt = LALt;
+
+ you have to use
+
+ o_info.LALt.real = LALt;
+
+ There are different MATENT_TYPES in the structure for LALt
+ (e.g. real, real_d or real_dd). The reason for the change is that
+ the old block-matrix support has been folded into the OPERATOR_INFO,
+ DOF_MATRIX and EL_MATRIX_INFO structures. In effect this means: the
+ old DOF_MATRIX etc. structures have been deleted, and the old
+ DOWB_... stuff has been renamed to old non-block names (Note: in the
+ vanilla v1.2 version there was no block-matrix support).
+
+ Another change in the OPERATOR_INFO structure is caused by the
+ changed treatment of boundary conditions:
+
+ o_info.use_get_bound = true;
+
+ does not exist anymore. Instead there is a component
+
+ o_info.dirichlet_bndry;
+
+ which is a bit-mask flagging those boundary segments which should be
+ treated as Dirichlet boundaries. In the new demo-program the
+ bit-mask of Dirichlet boundary segments is stored in the global
+ variable `dirichlet_mask', the following can be used to copy
+ boundary bit-masks:
+
+ BNDRY_FLAGS_CPY(o_info.dirichlet_bndry, dirichlet_mask);
+
+- The following functions have now different/more/less arguments:
+
+get_dof_matrix("A", row_fe_space, col_fe_space);
+ `col_fe_space' is a new argument and may be different from
+ `row_fe_space'. Passing NULL for `col_fe_space' means to use
+ `row_fe_space' also for `col_fe_space'.
+
+update_matrix(matrix, matrix_info, Transpose/NoTranspose);
+ `Transpose'/`NoTranspose' is the new argument; `Transpose'
+ instructs update_matrix() to add the transpose of the element
+ matrix to the system matrix, so `NoTranspose' should be used
+ to trigger the old behaviour.
+
+dirichlet_bound(f_h, u_h, bound, dirichlet_mask, g);
+ `dirichlet_mask' is the new argument. In the new demo program
+ it is declared as a global variable
+
+ static BNDRY_FLAGS dirichlet_mask;
+
+ The rationale is to separate the classification of boundary
+ segments from the definition of boundary conditions. The
+ boundary types assigned to boundary segments in the
+ macro-triangulation are now mere numbers, the interpretation
+ of those numbers is done by the application. To impose
+ Dirichlet b.c. on a boundary segment with boundary "type" N,
+ then it has to set bit N in that BNDRY_FLAGS bit-mask and pass
+ that mask to `dirichlet_bound()'.
+
+oem_solve_s(matrix, bound, f_h, u_h, solver, tol, precon, restart, miter, info);
+ `(const DOF_SCHAR_VEC *)bound' and `(const PRECON *)precon'
+ are new arguments. `precon' is an opaque handle describing a
+ preconditioner. It is initialised by a call to
+
+ precon =
+ init_oem_precon(matrix, NULL, info, icon, ssor_omega, ssor_iter);
+
+ellipt_est(u_h, adapt, rw_el_est, NULL/*rw_est_c*/, degree, norm, C,
+ (const REAL_D *)A, dirichlet_mask, r, r_flag,
+ NULL/*gn*/, 0/*gn_flags*/);
+
+ `(BNDRY_FLAGS)dirichlet_mask', `gn' and `(FLAGS)gn_flags' are
+ new arguments. `dirichlet_mask' has the same meaning as the
+ parameter passed to `dirichlet_bound()'; if a boundary segment
+ marked with a number `N' in the macro-triangulation should be
+ treated as Dirichlet boundary, then bit number `N' has to be
+ set in `dirichet_mask'. The parameter
+
+ (REAL (*gn)(const EL_INFO *el_info, const QUAD *quad, int qp,
+ REAL uh_qp, const REAL_D normal))
+
+ is a function pointer which should point to a function
+ implementing inhomogeneous Neumann b.c. It may be NULL.
+
+L2_err(u, u_h, quad, 0, false/*mean-value adjust*/,
+ NULL /* rw_err_el()*/, NULL /* max_err_el2 */);
+ L2_err() now can optionally do an automatic mean-value
+ adjustment. This is meant, e.g., for pure Neumann problems or
+ to compute the error for the pressure after solving a Stokes
+ problem with pure Dirichlet conditions for the velocity.
+
+H1_err(grd_u, u_h, NULL /* quad */, false /* relative error */,
+ NULL /* rw_err_el()*/, NULL /* max_err_el2 */);
+ Should be unchanged.
+
+GET_MESH(dim, "ALBERTA mesh", data, NULL /* init_node_projection */,
+ NULL /* init_wall_trafo */);
+ The additional parameter `init_node_projection()' is meant for
+ parametric meshes (see, e.g., ellipt-isoparam.c in the demo
+ package for the new version) and periodic boundary conditions
+ (see ellipt-periodic.c)
+
+MACRO_DATA *macro_data = read_macro(filename)
+ This belongs to the changes around the `GET_MESH()' and
+ `init_dof_admin()' stuff. Formerly `read_macro()' had to be
+ called with a half-initialised mesh, this is no longer
+ necessary. The returned data-handle `macro_data' has to be
+ freed again by `free_macro_data()'.
+
+get_adapt_stat(mesh->dim, "ellipt", "adapt", 2,
+ NULL /* ADAPT_STAT storage area, optional */);
+ The first parameter, the mesh-dimension, is new.
+
+################################################################################
+#
+# v1.2 ellipt.c ported to the most recent ALBERTA follows below. By
+# setting the NEW_VERSION `#define' at the start of the program (after
+# the copyright header) it is possible to switch between the old and
+# the new version. The `#if NEW_VERSION'-blocks also make it possible
+# to compare the old and the new version side-by-side (well, one above
+# the other).
+#
+/*--------------------------------------------------------------------------*/
+/* ALBERTA: an Adaptive multi Level finite element toolbox using */
+/* Bisectioning refinement and Error control by Residual */
+/* Techniques for scientific Applications */
+/* */
+/*--------------------------------------------------------------------------*/
+/* */
+/* file: ellipt.c */
+/* */
+/* description: solver for an elliptic model problem */
+/* */
+/* -\Delta u = f in \Omega */
+/* u = g on \partial \Omega */
+/* */
+/*--------------------------------------------------------------------------*/
+/* */
+/* 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 */
+/* */
+/* http://www.alberta-fem.de/ */
+/* */
+/*--------------------------------------------------------------------------*/
+/* (c) by A. Schmidt and K.G. Siebert (1996-2004) */
+/* */
+/* This program is free software; you can redistribute it and/or modify */
+/* it under the terms of the GNU General Public License as published by */
+/* the Free Software Foundation; either version 2 of the License, or */
+/* any later version. */
+/* */
+/* This program is distributed in the hope that it will be useful, */
+/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
+/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
+/* GNU General Public License for more details. */
+/*--------------------------------------------------------------------------*/
+
+#define NEW_VERSION 1
+
+#if NEW_VERSION
+#include "alberta-demo.h" /* proto-types for some helper functions */
+#else /* old version */
+#include <alberta.h>
+#endif
+
+
+/*--------------------------------------------------------------------------*/
+/* function for displaying mesh, discrete solution, and/or estimate */
+/* defined in graphics.c */
+/*--------------------------------------------------------------------------*/
+#if NEW_VERSION
+static bool do_graphics = true; /* global graphics "kill-switch" */
+#else
+void graphics(MESH *mesh, DOF_REAL_VEC *u_h, REAL (*get_est)(EL *el));
+#endif
+
+/*--------------------------------------------------------------------------*/
+/* global variables: finite element space, discrete solution */
+/* load vector and system matrix */
+/*--------------------------------------------------------------------------*/
+
+static const FE_SPACE *fe_space; /* initialized by init_dof_admin() */
+static DOF_REAL_VEC *u_h = nil; /* initialized by build() */
+static DOF_REAL_VEC *f_h = nil; /* initialized by build() */
+static DOF_MATRIX *matrix = nil; /* initialized by build() */
+
+#if NEW_VERSION
+static BNDRY_FLAGS dirichlet_mask; /* bit-mask of Dirichlet segments */
+#endif
+
+/*--------------------------------------------------------------------------*/
+/* init_dof_admin(): init DOFs for Lagrange elements of order */
+/* <polynomial degree>, called by GET_MESH() */
+/*--------------------------------------------------------------------------*/
+#if NEW_VERSION
+#else /* old version */
+void init_dof_admin(MESH *mesh)
+{
+ FUNCNAME("init_dof_admin");
+ int degree = 1;
+ const BAS_FCTS *lagrange;
+
+ GET_PARAMETER(1, "polynomial degree", "%d", °ree);
+ lagrange = get_lagrange(degree);
+ TEST_EXIT(lagrange)("no lagrange BAS_FCTS\n");
+ fe_space = get_fe_space(mesh, lagrange->name, nil, lagrange);
+ return;
+}
+#endif
+
+/*--------------------------------------------------------------------------*/
+/* struct ellipt_leaf_data: structure for storing one REAL value on each */
+/* leaf element as LEAF_DATA */
+/* init_leaf_data(): initialization of leaf data, called by GET_MESH() */
+/* 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 */
+};
+
+#if NEW_VERSION
+#else /* old version */
+void init_leaf_data(LEAF_DATA_INFO *leaf_data_info)
+{
+ leaf_data_info->leaf_data_size = sizeof(struct ellipt_leaf_data);
+ leaf_data_info->coarsen_leaf_data = nil; /* no transformation */
+ leaf_data_info->refine_leaf_data = nil; /* no transformation */
+ return;
+}
+#endif
+
+static REAL *rw_el_est(EL *el)
+{
+ if (IS_LEAF_EL(el))
+ return(&((struct ellipt_leaf_data *)LEAF_DATA(el))->estimate);
+ else
+ return(nil);
+}
+
+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) */
+/*--------------------------------------------------------------------------*/
+
+static REAL u(const REAL_D x)
+{
+ return(exp(-10.0*SCP_DOW(x,x)));
+}
+
+static const REAL *grd_u(const REAL_D x, REAL_D grd)
+{
+ static REAL_D buffer;
+ REAL ux = exp(-10.0*SCP_DOW(x,x));
+ int n;
+
+ if(!grd) {
+ grd = buffer;
+ }
+
+ for (n = 0; n < DIM_OF_WORLD; n++)
+ grd[n] = -20.0*x[n]*ux;
+
+ return(grd);
+}
+
+/*--------------------------------------------------------------------------*/
+/* problem data: right hand side, boundary values */
+/*--------------------------------------------------------------------------*/
+
+static REAL g(const REAL_D x) /* boundary values, not optional */
+{
+ return(u(x));
+}
+
+static REAL f(const REAL_D x) /* -Delta u, not optional */
+{
+ REAL r2 = SCP_DOW(x,x), ux = exp(-10.0*r2);
+ return(-(400.0*r2 - 20.0*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_info: 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() after init_element() */
+/*--------------------------------------------------------------------------*/
+
+struct op_info
+{
+#if NEW_VERSION
+ REAL_BD Lambda; /* equivalent to REAL_D Lambda[N_LAMBDA_MAX]; */
+#else /* old version */
+ REAL_D Lambda[DIM+1]; /* the gradient of the barycentric coordinates */
+#endif
+ REAL det; /* |det D F_S| */
+};
+
+#if NEW_VERSION
+static bool init_element(const EL_INFO *el_info, const QUAD *quad[3], void *ud)
+{
+ struct op_info *info = ud;
+
+ info->det = el_grd_lambda(el_info, info->Lambda);
+ return false; /* not parametric */
+}
+#else /* old version */
+static void init_element(const EL_INFO *el_info, const QUAD *quad[3], void *ud)
+{
+ struct op_info *info = ud;
+
+ info->det = el_grd_lambda(el_info, info->Lambda);
+ return;
+}
+#endif
+
+#if NEW_VERSION
+const REAL_B *LALt(const EL_INFO *el_info, const QUAD *quad,
+ int iq, void *ud)
+{
+ struct op_info *info = (struct op_info *)ud;
+ int i, j, dim = el_info->mesh->dim;
+ static REAL_BB LALt;
+
+ 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);
+}
+#else /* old version */
+const REAL (*LALt(const EL_INFO *el_info, const QUAD *quad,
+ int iq, void *ud))[DIM+1]
+{
+ struct op_info *info = ud;
+ int i, j, k;
+ static REAL LALt[DIM+1][DIM+1];
+
+ 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 (*)[DIM+1]) LALt);
+}
+#endif
+
+static void build(MESH *mesh, U_CHAR flag)
+{
+ FUNCNAME("build");
+ static const EL_MATRIX_INFO *matrix_info = nil;
+ const QUAD *quad;
+
+ dof_compress(mesh);
+ MSG("%d DOFs for %s\n", fe_space->admin->size_used, fe_space->name);
+
+ if (!u_h) /* access matrix and vector for linear system */
+ {
+#if NEW_VERSION
+ matrix = get_dof_matrix("A", fe_space, NULL /* col_fe_space */);
+#else
+ matrix = get_dof_matrix("A", fe_space);
+#endif
+ 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 ! */
+ }
+
+ if (!matrix_info) /* information for matrix assembling */
+ {
+ OPERATOR_INFO o_info = {nil};
+
+ o_info.row_fe_space = o_info.col_fe_space = fe_space;
+ o_info.init_element = init_element;
+#if NEW_VERSION
+ o_info.LALt.real = LALt;
+#else
+ o_info.LALt = LALt;
+#endif
+ o_info.LALt_pw_const = true; /* pw const. assemblage is faster */
+ o_info.LALt_symmetric = true; /* symmetric assemblage is faster */
+#if NEW_VERSION
+ BNDRY_FLAGS_CPY(o_info.dirichlet_bndry,
+ dirichlet_mask); /* Dirichlet bndry conditions */
+#else
+ o_info.use_get_bound = true; /* Dirichlet boundary conditions! */
+#endif
+ o_info.user_data = MEM_ALLOC(1, struct op_info); /* user data! */
+ o_info.fill_flag = CALL_LEAF_EL|FILL_COORDS;
+
+ matrix_info = fill_matrix_info(&o_info, nil);
+ }
+
+ clear_dof_matrix(matrix); /* assembling of matrix */
+#if NEW_VERSION
+ update_matrix(matrix, matrix_info, NoTranspose);
+#else
+ update_matrix(matrix, matrix_info);
+#endif
+
+ dof_set(0.0, f_h); /* assembling of load vector */
+ quad = get_quadrature(mesh->dim, 2*fe_space->bas_fcts->degree - 2);
+ L2scp_fct_bas(f, quad, f_h);
+
+#if NEW_VERSION
+ dirichlet_bound(f_h, u_h, NULL /* bound */, dirichlet_mask, g);
+#else
+ dirichlet_bound(g, f_h, u_h, nil); /* boundary values */
+#endif
+
+ return;
+}
+
+
+/*--------------------------------------------------------------------------*/
+/* solve(): solve the linear system, called by adapt_method_stat() */
+/*--------------------------------------------------------------------------*/
+
+static void solve(MESH *mesh)
+{
+ FUNCNAME("solve");
+#if NEW_VERSION
+ static REAL tol = 1.e-8, ssor_omega = 1.0;;
+ static int miter = 1000, info = 2, icon = 1;
+ static int restart = 0,ssor_iter = 1;
+ const PRECON *precon;
+#else
+ static REAL tol = 1.e-8;
+ static int miter = 1000, info = 2, icon = 1, restart = 0;
+#endif
+ static OEM_SOLVER solver = NoSolver;
+
+ 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 NEW_VERSION
+ GET_PARAMETER(1, "precon ssor omega", "%f", &ssor_omega);
+ GET_PARAMETER(1, "precon ssor iter", "%d", &ssor_iter);
+#endif
+ if (solver == GMRes)
+ GET_PARAMETER(1, "solver restart", "%d", &restart);
+ }
+
+#if NEW_VERSION
+ precon = init_oem_precon(matrix, NULL, info, icon, ssor_omega, ssor_iter);
+ oem_solve_s(matrix, NULL /* bound */, 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 */);
+ }
+#else /* old version */
+ oem_solve_s(matrix, f_h, u_h, solver, tol, icon, restart, miter, info);
+
+ graphics(mesh, u_h, nil);
+#endif
+
+ 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_iq,
+ const REAL_D grd_uh_iq)
+{
+ REAL_D x;
+ coord_to_world(el_info, quad->lambda[iq], x);
+ return(-f(x));
+}
+
+#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 degree, norm = -1;
+ static REAL C[3] = {1.0, 1.0, 0.0};
+ static REAL est, est_old = -1.0, err, err_old = -1.0;
+ static FLAGS r_flag = 0; /* = (INIT_UH | INIT_GRD_UH), if needed by r() */
+ REAL_DD A = {{0.0}};
+ int n;
+ const QUAD *quad;
+
+ for (n = 0; n < DIM_OF_WORLD; n++)
+ A[n][n] = 1.0; /* set diogonal 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]);
+ }
+ degree = 2*u_h->fe_space->bas_fcts->degree;
+#if NEW_VERSION
+ est = ellipt_est(u_h, adapt, rw_el_est, nil, degree, norm, C,
+ (const REAL_D *) A, dirichlet_mask, r, r_flag,
+ NULL /* gn */, 0 /* gn_flags*/);
+#else
+ est = ellipt_est(u_h, adapt, rw_el_est, nil, degree, norm, C,
+ (const REAL_D *) A, r, r_flag);
+#endif
+
+ MSG("estimate = %.8le", est);
+ if (est_old >= 0)
+ print_msg(", EOC: %.2lf\n", EOC(est,est_old));
+ else
+ print_msg("\n");
+ est_old = est;
+
+ quad = get_quadrature(mesh->dim, degree);
+ if (norm == L2_NORM)
+#if NEW_VERSION
+ err = L2_err(u, u_h, quad, 0, false/* 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 */);
+
+#else
+ err = L2_err(u, u_h, quad, 0, nil, nil);
+ else
+ err = H1_err(grd_u, u_h, quad, 0, nil, nil);
+#endif
+
+ 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 NEW_VERSION
+ if (do_graphics) {
+ MSG("Displaying the estimate.\n");
+ graphics(mesh, NULL /* u_h */, get_el_est, NULL /* u_exact() */,
+ HUGE_VAL /* time */);
+ }
+#else
+ graphics(mesh, nil, get_el_est);
+#endif
+
+ return(adapt->err_sum);
+}
+
+/*--------------------------------------------------------------------------*/
+/* main program */
+/*--------------------------------------------------------------------------*/
+
+int main(int argc, char **argv)
+{
+ FUNCNAME("main");
+ MESH *mesh;
+#if NEW_VERSION
+ MACRO_DATA *data;
+ const BAS_FCTS *lagrange;
+ int degree = 1, dirichlet_bit = 1, dim;
+#endif
+ int n_refine = 0;
+ static ADAPT_STAT *adapt;
+ char filename[100];
+
+/*--------------------------------------------------------------------------*/
+/* first of all, init parameters of the init file */
+/*--------------------------------------------------------------------------*/
+
+ init_parameters(0, "INIT/ellipt.dat");
+
+/*--------------------------------------------------------------------------*/
+/* get a mesh, and read the macro triangulation from file */
+/*--------------------------------------------------------------------------*/
+
+#if NEW_VERSION
+ GET_PARAMETER(1, "mesh dimension", "%d", &dim);/* neu eingefuegt */
+ GET_PARAMETER(1, "polynomial degree", "%d", °ree);/* neu eingefuegt */
+ GET_PARAMETER(1, "macro file name", "%s", filename);
+ GET_PARAMETER(1, "global refinements", "%d", &n_refine);
+ GET_PARAMETER(1, "dirichlet boundary", "%d", &dirichlet_bit);
+
+ data = read_macro(filename);
+ mesh = GET_MESH(dim, "ALBERTA mesh", data, NULL /* init_node_projection */,
+ NULL /* init_wall_trafo */);
+ 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);
+
+ if (do_graphics) {
+ MSG("Displaying the mesh.\n");
+ graphics(mesh, NULL /* u_h */, NULL /* get_est()*/ , NULL /* u_exact() */,
+ HUGE_VAL /* time */);
+ }
+
+#else /* old version */
+ mesh = GET_MESH("ALBERTA mesh", init_dof_admin, init_leaf_data);
+ GET_PARAMETER(1, "macro file name", "%s", filename);
+ read_macro(mesh, filename, nil);
+ GET_PARAMETER(1, "global refinements", "%d", &n_refine);
+ global_refine(mesh, n_refine*DIM);
+ graphics(mesh, nil, nil);
+
+#endif
+
+
+#if NEW_VERSION
+ lagrange = get_lagrange(mesh->dim, degree);
+ TEST_EXIT(lagrange, "no lagrange BAS_FCTS\n");
+ fe_space = get_fe_space(mesh, lagrange->name, lagrange, 1 /* rdim */,
+ ADM_FLAGS_DFLT);
+
+ if (dirichlet_bit > 0) {
+ BNDRY_FLAGS_SET(dirichlet_mask, dirichlet_bit);
+ }
+#endif /* analogue to the function 'init_dof_admin()', see above */
+
+/*--------------------------------------------------------------------------*/
+/* init adapt structure and start adaptive method */
+/*--------------------------------------------------------------------------*/
+
+#if NEW_VERSION
+ adapt = get_adapt_stat(mesh->dim, "ellipt", "adapt", 2,
+ NULL /* ADAPT_STAT storage area, optional */);
+#else
+ adapt = get_adapt_stat("ellipt", "adapt", 2, nil);
+#endif
+ adapt->estimate = estimate;
+ adapt->get_el_est = get_el_est;
+ adapt->build_after_coarsen = build;
+ adapt->solve = solve;
+
+ adapt_method_stat(mesh, adapt);
+
+ WAIT_REALLY;
+ return(0);
+}
+################################### END #######################################
+
+ LocalWords: ellipt gltools libgltools libdemo barycentric