Skip to content
Snippets Groups Projects
Commit 620ea575 authored by Claus-Justus Heine's avatar Claus-Justus Heine
Browse files

V1.2 -> post-2.0 porting guide.

parent 92125c4d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 813 additions and 0 deletions
howto-port-ellipt.txt 0 → 100644
+ 813
− 0
View file @ 620ea575
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", &degree);
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", &degree);/* 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment