NAME
FEM - base class and methods for finite element programming
INCLUDE
include "FEM.h"
SYNTAX
class FEM : public SimCase
{
public:
static bool optimize; // true: allow optimizations in Dp-FEM
// for indicating arithmetic or harmonic mean in smoothing functions:
enum FieldAverage_type { ARITHMETIC = 0, HARMONIC = 1 };
// for smoothing method:
enum Smoothing_type { GLOBAL_LS = 0, MOVING_LS = 1 };
protected:
ElmItgRules itg_rules; // used in makeSystem
FiniteElement finite_elm; // used in makeSystem
ElmMatVec elm_matvec; // used in makeSystem
Handle(MassMatrix) mass_mat; // used for smoothing derivatives
//Handle(DegFreeFE) dof4massmat; // used in smoothFields
MovingLS MLS; // used for smoothing derivatives
FieldsFEatItgPt field_at_pts; // member for efficiency in makeFlux etc
FieldAverage_type average_tp; // used in smoothing proc., set on menu
// measure the CPU time of the assembly process:
CPUclock cpu;
double cpu_time_makeSystem;
void reportCPUtime (const char* algorithm, const char* comment = NULL);
void checkMassMatrix (GridFE& grid);
public:
FEM ();
virtual ~FEM ();
virtual void makeSystem
(
DegFreeFE& dof,
LinEqAdmFE& lineq,
bool compute_A = true,
bool compute_RHS = true,
bool only_safe_nopt = true
);
// the default version is sufficient for volume integration (not sides):
virtual void calcElmMatVec (int e, ElmMatVec& elmat, FiniteElement& fe);
// using FEM-derived integrand definitions:
virtual void numItgOverElm (ElmMatVec& elmat, FiniteElement& fe);
// called from numItgOverElm:
virtual void integrands (ElmMatVec& elmat, const FiniteElement& fe);
virtual void numItgOverSide
(
int side,
int boind,
ElmMatVec& elmat,
FiniteElement& fe
);
virtual void integrands4side
(
int side, int boind, ElmMatVec& elmat, const FiniteElement& fe
);
// --- some utility functions ---
void detach ();
virtual void derivedQuantitiesAtItgPt
(VecSimple(NUMT)& quantities, const FiniteElement& fe);
void makeMassMatrix (GridFE& grid, Matrix(NUMT)& mm, bool lumped);
void makeSystem
(
DegFreeFE& dof,
ElmMatVecCalc& emv,
IntegrandCalc* integrand, // if NULL: emv supplies the integrand
LinEqAdmFE& lineq,
bool compute_A = true,
bool compute_RHS = true
);
void makeSystem
(
DegFreeFE& dof,
ElmMatVecCalc& emv,
IntegrandCalc* integrand, // if NULL: emv supplies the integrand
Matrix(NUMT)& mt
);
void makeSystem
(
DegFreeFE& dof,
Matrix(NUMT)& mt // calls subclass of FEM integrands
);
void makeSystem
(
DegFreeFE& dof,
Matrix(NUMT)& mt, // calls subclass of FEM integrands
Vector(NUMT)& vc // calls subclass of FEM integrands
);
void makeSystem
(
DegFreeFE& dof,
ElmMatVecCalc& emv,
IntegrandCalc* integrand, // if NULL: emv supplies the integrand
Vector(NUMT)& vc
);
void makeSystem
(
DegFreeFE& dof,
Vector(NUMT)& vc
);
virtual void numItgOverElm
(
ElmMatVec& elmat,
FiniteElement& fe,
IntegrandCalc& integrand // element matrix/vector at an itg.pt.
);
virtual void numItgOverSide
(
int side,
int boind,
ElmMatVec& elmat,
FiniteElement& fe,
IntegrandCalc& integrand // element matrix/vector at an itg.pt.
);
// -------------- derivatives and fluxes (with/without smoothing) ----------
void attachMassMatrix (MassMatrix& mass_matrix)
{ mass_mat.rebind (mass_matrix); }
void makeFlux
(
FieldsFE& flux, // flux = smooth[c*k*grad(f)]
const FieldFE& f,
Field& k,
real c = -1.0,
enum Smoothing_type method = GLOBAL_LS,
real t = DUMMY
);
void makeFlux // requires virtual function k
(
FieldsFE& flux, // flux = smooth[-k*grad(f)]
const FieldFE& f,
enum Smoothing_type method = GLOBAL_LS,
real t = DUMMY
);
// used by the makeFlux function above (to be impl. in subclass of FEM):
virtual NUMT k (const FiniteElement& fe, real t = DUMMY);
void smoothFields (FieldsFE& smooth_field,
const FieldsFEatItgPt& some_field,
enum Smoothing_type method = GLOBAL_LS);
void smoothFields (FieldsFE& smooth_field,
const FieldsPiWisConst& piecewise_constant_field);
void smoothField (FieldFE& smooth_field,
const Field& some_field);
void smoothField
(
FieldFE& smooth_field,
IntegrandCalc& rhs_integrand, // defines field to be smoothed
int relative_itg_order = -1 // (-1: reduced Gauss points)
);
void smoothDerivative (FieldFE& df_dir, int dir, const FieldFE& f);
void smoothGradient (FieldsFE& gradf, const FieldFE& f);
void makeGradient (FieldsFEatItgPt& gradf, const FieldFE& f)
{ gradf.gradient (f, GAUSS_POINTS, 0); }
// computes grad(f) as the average of optimal grad(f) values in each elm:
void makeGradient (FieldsPiWisConst& gradf, const FieldFE& f);
// ||grad(f)|| as the average of optimal ||grad(f)|| values in each elm:
void makeGradient (FieldPiWisConst& gradf_norm, const FieldFE& f);
void smoothMultipleFields
(FieldsFE& collection, Matrix(NUMT)& mass_matrix);
virtual void integrand4smoothing
(FiniteElement& fe, MatSimple(NUMT)& contributions);
// ------- backward compatible functionality for smoothing/derivatives -----
// (these functions issue warning message and instructs the user to
// reimplement the functionality using the functions above)
void smoothField
(
FieldFE& smooth_field,
Matrix(NUMT)& mass_matrix,
Field& some_field // field to be smoothed
);
void smoothFields
(
FieldsFE& smooth_field,
Matrix(NUMT)& mass_matrix,
FieldsPiWisConst& piecewise_constant_field
);
void smoothFields
(
FieldsFE& smooth_field,
FieldsPiWisConst& piecewise_constant_field // field to be smoothed
); // makes an internal mass matrix
void smoothField
(
FieldFE& smooth_field,
Matrix(NUMT)& mass_matrix,
IntegrandCalc& rhs_integrand, // def. of field to be smoothed
int relative_itg_order = -1 // (-1: reduced Gauss points)
);
void smoothDerivative
(
FieldFE& df,
int dir,
Matrix(NUMT)& mass_matrix,
FieldFE& f
);
void smoothGradient
(
FieldsFE& df,
Matrix(NUMT)& mass_matrix,
FieldFE& f
);
static void makeGradientStatic // computes grad(f) at the centroids
(
FieldsPiWisConst& gradient,
FieldFE& f
);
static void makeGradientStatic // computes ||grad(f)|| at the centroids
(
FieldPiWisConst& gradient,
FieldFE& f
);
void makeFlux
(
FieldsPiWisConst& v, // piecewise constant flux: v = c*k*grad(u)
FieldFE& u,
real c = 1.0,
Field* k = NULL,
real t = DUMMY // used when evaluating k
);
void makeFlux
(
FieldsFE& v, // smooth flux: v = c*k*grad(u)
FieldFE& u,
int method, // 1: v=c*k*smooth[grad u], 2: v=smooth[c*k*grad u]
Matrix(NUMT)* mass_matrix, // if NULL: make an internal mass matrix
real c = 1.0,
Field* k = NULL,
real t = DUMMY
);
// --------------------------------------------------------------------------
static void defineStatic (MenuSystem& menu, int level = MAIN);
void define (MenuSystem& menu, int level = MAIN) { defineStatic(menu,level);}
void scan (MenuSystem& menu);
double getCPUtime4makeSystem () const { return cpu_time_makeSystem; }
// ------- simulator interface for adaptive grid and other functionality ----
virtual real integrands4energyErrorNorm (const FiniteElement& fe);
virtual void evalOwnRefInd (FieldPiWisConst& refinement_field,
real& evaluated_error);
// used to set solver dependent parameters in element definitions:
virtual void customizeElmDef (ElmDef& elmdef, const String& elm_name);
};
KEYWORDS
finite element method, base class, problem dependent data
DESCRIPTION
Class FEM offers the standard finite element algorithms for the
assembly loop over the elements, including numerical integration
over elements. In addition, class FEM offers computation and
smoothing of (discontinuous) derivatives of finite element
fields. Usually, FEM is used in conjuction with a simulator
class and serves as base class for the simulator. All the stan
dard finite element algorithms are therehy inherited in the simu
lator class. These standard algorithms set up the loops and
peform all actions that are common to "all" problems. The prob
lem dependent parts are defined in virtual functions that are
called from these algorithms. Some virtual functions must be
implemented in the simulator class (typically the one that
defines the integrands in the weak form - otherwise the algo
rithms do not know which PDE to solve!).
One can apply class FEM without knowing the details of the imple
mentation. However, it is strongly recommended that serious
finite element programmers study the bodies of the class FEM
functions to understand how the building blocks of the finite
element toolbox in Diffpack work. Start with makeSystem and fol
low the program flow until the element loop is finished. The
understanding of this class requires the reader to be familiar
with finite element programming, C++ and object-oriented program
ming.
Class FEM can perform some algorithmic optimizations if the
static variable FEM::optimize is true (or ON). This can be turned
on by the statement FEM::optimize=ON in the code, preferably in
main, or by the command line argument --FEM::optimize ON. The
value of this optimization indicator is written to the logfile
(casename.dp). The optimizations assumes that standard Diffpack
conventions are followed. If you program in a different way, or
work with a problem where you have your own special construc
tions, the optimizations may lead to wrong answers! Therefore,
the optimizations are off by default and the user should only
turn them on after having verified the code thoroughly.
CONSTRUCTORS AND INITIALIZATION
A constructor without arguments is available. No further ini
tialization is necessary. The define and scan functions makes it
possible to alter the numerical integration rules using a menu in
the standard way. There is also a menu item related to average
vs harmonic mean when calculating fluxes. The default menu
answers will suffice in many cases.
MEMBER FUNCTIONS
makeSystem - computes the linear system arising from a finite
element method. The function is completely general and calls vir
tual functions for defining the element matrix and vector. A
default version of this virtual function is provided. The default
version simply calls a function for numerical integration over
the element which again calls a virtual function integrands that
the programmer must implement in the derived simulator class. If
only_safe_opt is false, special optimization might be turned on.
The code is then not necessarily valid in all cases. What type of
optimizations that are available and their validity will be docu
mented here. Typically, all elements must be of the same type
and no surface integrals must appear in the weak form. (We can
test that all elements are of the same type, but not if the simu
lator developer also has a surface integral! The only_safe_opt
flag must therefore be used with care.)
There are several overloaded versions of the makeSystem function.
Some functions are aimed at computing the coefficient matrix only
or the right hand side only. The programmer is encouraged to read
the body of the makeSystem functions since this will provide all
details about the different versions. Some of the versions
require particular initializations of the involved data struc
tures and this is explained in comments in the source code. If
the user needs to form linear systems corresponding to different
weak forms, one can create special "integrands functors" upon
which special makeSystem functions can base the computation. If
none of the makeSystem functions are appropriate, one can easily
make one's own makeSystem in a subclass of FEM. This might be
required to construct an efficient assembly process, taking
advantage of special features in the problem.
The makeSystem functions measures the CPU time spent in the func
tion. This CPU time is written to the logfile (casename.dp) if
more than half a second is spent in the makeSystem function.
Similar CPU time information is provided by LinEqAdmFE::solve and
written on the same line in the logfile. The CPU time is avail
able in an internal (protected) variable that the programmer can
use in a derived simulation class if desired.
calcElmMatVec - calculates the element matrix and vector. When
the programmer has analytical expressions for the element vector
and matrix available, these expressions can simply be inserted in
a subclass (simulation class) implementation of calcElmMatVec.
If numerical integration is required, class FEM provides a
default implementation of calcElmMatVec that performs necessary
initializations and a call to numItgOverElement. The latter
administers the numerical integration and samples the integrand
by calling a function integrands that must be implemented in the
simulator class.
In cases where a surface integral enters the weak form in addi
tion to the volume integral over the element's interior, the pro
grammer must implement a local version of calcElmMatVec. Such an
implementation will typically consist of a call to FEM::calcElm
MatVec for performing the integration over the interior of the
element and thereafter a loop over the sides of the element is
coded. In the loop, one tests if the relevant boundary indicators
are set (these indicators mark parts of the boundaries where the
surface integral is to be evaluated), and thereafter the "inte
grands4side function is called. This is the counterpart to inte
grands and must of course be implemented in the subclass (simula
tor). The complete weak form is hence defined by integrands and
integrands4side.
If the programmer needs to compute the element matrix and vector
for different weak forms, the calcElmMatVec function can be rep
resented as a functor, using class ElmMatVecCalc. Particularly,
a default functor corresponding to the FEM::calcElmMatVec imple
mentation exists in class ElmMatVecCalcStd.
"numItgOverElm - numerical integration over the element. The
FiniteElement object must be initialized by calling its refill
function prior to the call to numItgOverElm. The function evalu
ates the integrand of the integral by calling the virtual func
tion integrands. If the integrand is defined by a functor
(derived from IntegrandCalc) rather than a class FEM virtual
function integrands, an overloaded version of numItgOverElm can
handle that case.
numItgOverSide - numerical integration over the side (surface) of
an element. Contrary to numItgOverElm the FiniteElement object
is automatically refilled for the current side (refill4side)
inside numItgOverSide. (However, refill must be called, but that
is usually performed prior to a numItgOverElm call which normally
preceeds the computation of the side integral). An overloaded
version can handle the case where the integrand in the side (sur
face) integral is defined in terms of an IntegrandCalc functor.
integrands4energyErrorNorm - a function similar to integrands,
but the purpose is to evaluate the contribution (from an integra
tion point in an element) to the energy norm of the error. Hence,
one must in this function have access to the primary unknown of
the problem and some reference solution (e.g. an exact solution
or a fine grid simulation).
attachMassMatrix - this function allows the user to explicitly
attach his own mass matrix. If the function is not called, the
various other functions for computing smooth derivatives and
fluxes will use an internal mass matrix, which is automatically
re-used in later computations (unless the grid changes during a
simulation). In practice, many simulators will never have or
attach a mass matrix, but rely on FEM and its internal mass
matrix computations. If the mass matrix is to be used directly
in the solver, it is available as the mass_mat data member.
The next functions concerns computations of derivatives of finite
element fields. As these derivatives are discontinuous, various
smoothing algorithms are offered. All the following functions
were revised in Diffpack v2.8. The old functions are available
for backward compatibility, but give warning messages that some
thing new and better is available.
smoothField, makeGradient, makeFlux, smoothDerivative, "smooth
Gradient - enable computation and smoothing of derivatives of
finite element fields. makeFlux is an all-round function that is
sufficient for most users. See the description of makeFlux
below.
The most fundamental smoothField function runs the least squares
method on a system where a IntegrandCalc functor defines the
right-hand side, i.e. the field to be smoothed. An assembly pro
cess over the elements is set up, the mass matrix is computed if
it is not computed before or attached by the user, and the
resulting system is solved.
Other versions of smoothField offer more user-friendly interfaces
in that they take ordinary Diffpack field objects as input and
define the corresponding functor for the right-hand side in the
least squares system. Of particular interest is the function
that takes a FieldsFEatItgPt object, usually containing deriva
tives at reduced Gauss points in the elements, and runs the least
squares process on this collection of discrete values.
The smoothField functions are usually used for smoothing deriva
tives. Direct computation of derivatives are offered by the
functions "smoothDerivative (smooth derivative in one direction),
makeGradient (all derivatives, but no smoothing) and smoothGradi
ent (computes all directional derivatives and smooths the
result). We also refer to class FieldsFEatItgPt for functional
ity concerning computation of derivatives (this functionality is
much used in the FEM functions as well).
Two makeGradient functions represent the gradient of a field in
terms of piecewise constant fields. When using multilinear or
linear elements this is a natural representation of derivatives.
For higher order elements these functions first computes the
derivatives at the reduced Gauss points and then assign the aver
age of these point values to the element value in a
Field(s)PiWisConst field. The FieldsPiWisConst field can be
plotted directly by the SimRes2mtv::plotmtvVector function or it
can be smoothed by FEM::smoothFields and then stored as a general
finite element vector field and later plotted by any visualiza
tion system that handles finite element vector fields.
The main function of all functions for smoothing and derivative
computation in class FEM is makeFlux. This is a flexible function
for computing c*k*grad(f), where c is scalar parameter (e.g.
minus one in many cases, k is a coefficient in the expression for
the flux, and f is finite element field. There are two basic rep
resentations of k, either as a subclass of Field or as a virtual
function k in the solver class. Simple solvers usually apply the
virtual function representation (see e.g. class Elliptic2 in the
basic Diffpack documentation), while the field representation
might be preferable in more complicated simulators. If one calls
the version of makeFlux where there is no explicit k argument,
the virtual k function must be implemented in the solver class
(otherwise an instructive error message is issued).
The makeFlux function first computes the derivatives at the opti
mal points in each element (the reduced Gauss points), using
class FieldsFEatItgPt functionality. Then these derivatives are
multiplied by the k values, and then calls a smoothField function
for smoothing the result. If the k field is not given as a param
eter, it is taken as unity. Finally, the smoothed field is mul
tiplied by c.
There are two smoothing methods in makeFlux; global least squares
and moving least squares. The former is most efficient, but the
latter has a higher convergence rate and is recommended when high
accuracy is important.
There is a data member average_tp (in clas FEM) which indicates
the type of averaging procedures that should be used in smoothing
procedures (smoothField, makeFlux and similar routines). By
default, the method gives arithmetic averages, but we have a
developed a special method that mimics harmonic means. You can
set the arithmetic or harmonic mean option on the menu, or set
the average_tp variable explicitly in the simulator class (which
is normally a subclass of FEM if finite element methods are used
in the simulator).
smoothMultipleFields - runs the least squares smoothing algorithm
for several fields simultaneously. The functions smoothFields
invokes a smoothing procedure for each component. In other words,
the assembly of the right hand side of the resulting linear sys
tem will be performed independently for each field (vector compo
nent) that is to be smoothed. The smoothMultipleFields function
has an assembly loop that computes several right hand sides
simultaneously. Of course, this requires a special "integrands"
function where the programmer can fill several element vectors at
the same time (each element vector corresponds to the elmat.b
quantity in the standard integrands routine). Instad of sending
an ElmMatVec object to the "integrands" routine, one sends
instead a Mat object, where each column is an element vector. The
function that defines the integrand has the name inte
grands4smoothing. The smoothMultipleFields function can give
significant speed up if the smoothing operations are a time con
suming part of a simulation. The programmer should consult the
source code before the function is applied such that a better
understanding of the use of the function arguments is obtained.
The function is not well tested (15/10-94 - no bug reports pr
1/11-97).
modifyField4Constraints - In case of a GridFEAdB refined grid and
a FieldFE object that has been computed by a standard smoothing
technique, the field values will usually not fulfill the con
straints. This function enforces the constraints explicitly on
the vector in the field. The modifyField4Constraint function is
automatically called from all functions in class FEM that perform
smoothing of fields. Hence, this function is only needed for
programmers that perform the smoothing locally in the simulator.
EXAMPLES
The member functions of class FEM are good examples on how both
other member functions and the general Diffpack building blocks
for finite element computations are used.
Here is an example on calling the functions that can compute the
derivatives of a field and smooth them:
// given FieldFE u, and GridFE grid
FieldsPiWisConst pwc (grid(), true, "pwc");
MatDiag(real) mass_matrix (grid->getNoNodes());
makeMassMatrix (grid(), mass_matrix, true);
mass_matrix.factLU(); // more efficient than factChol for
diagonal matrix!
FieldsFE gradient (grid(), "smooth_gradient_1");
makeGradient (pwc, u());
smoothFields (gradient, mass_matrix, pwc);
smoothFields (gradient, pwc); // FEM-internal
mass matrix
smoothGradient (gradient, mass_matrix, u());
FieldFE magnitude; // FieldsFE::magnitude can initialize the
field
gradient.magnitude (magnitude);
See also the report Details of Finite Element Programming in
Diffpack.
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO