NAME
ElmItgRules - numerical quadrature over a finite element
INCLUDE
include "ElmItgRules.h"
SYNTAX
class ElmItgRules : public HandleId
{
private:
NumItgPoints pt_tp; // type of points (Gauss or nodal pt)
NumItgDomain domainshape; // shape of domain, found from an ElmDef
ElementType el_tp; // identifier for element
Ptv(int) dir_order; // absolute order in the space dir.
Ptv(int) relative_order; // order = ElmDef's order + relative_order
int nsd; // number of space dimensions (physical space)
int nsd_itg_domain; // nsd for integration domain (=2 for 3D surf.)
VecSimple(Ptv(real)) points;
VecSimple(real) weights;
// help arrays for combining 1D rules in, e.g., rectangles/boxes:
VecSimplest(VecSimple(Ptv(real))) dimension_dep_points;
VecSimplest(VecSimple(real)) dimension_dep_weights;
VecSimple(Ptv(real)) points_side; // help array for integr. on element sides
int k; // local loop-counter for efficiency in member func.
int rel_order; // used to remember menu item or constructor arg
bool ignore_refill;
bool need_refill; // this variable is set true if status changed
public:
ElmItgRules (NumItgPoints point_type = GAUSS_POINTS, int relative_order = 0);
~ElmItgRules ();
bool ok () const;
void print (Os os) const; // print this object
void print (Os os, // print any points and weights data structure
const VecSimple(Ptv(real))& pt,
const VecSimple(real)& w) const;
void scan (Is is);
static void defineStatic (MenuSystem& menu, int level = MAIN);
void define (MenuSystem& menu, int level = MAIN) { defineStatic(menu,level);}
void scan (MenuSystem& menu);
Ptv(int) getRelativeOrder () const { return relative_order; }
void setRelativeOrder (int rel_order);
void setRelativeOrder (const Ptv(int)& directionwise_rel_order);
NumItgPoints getPointType () const { return pt_tp; }
void setPointType (NumItgPoints pt_tp);
// if the user wants to give the integration points and weights
// explicitly and avoid initialization actions for every element by
// class FiniteElement, these functions are useful:
void ignoreRefill (); // refill calls have no effect
void acceptRefill (); // refill calls have their normal effect
void setPointsAndWeights (const VecSimple(Ptv(real))& points,
const VecSimple(real)& weights);
bool operator == (const ElmItgRules& r) const;
bool operator != (const ElmItgRules& r) const
{ return notbool(operator==(r)); }
// functions for defining the standard textbook rules for intervals
// (-1,1), triangles and tetrahedras:
static void rulesInterval (VecSimple(Ptv(real))& pt,
VecSimple(real)& w,
NumItgPoints pt_type,
int order);
static void rulesTriangle (VecSimple(Ptv(real))& pt,
VecSimple(real) & w,
NumItgPoints pt_type,
int order);
static void rulesTetrahedra (VecSimple(Ptv(real))& pt,
VecSimple(real) & w,
NumItgPoints pt_type,
int order);
// function for reducing the length of a standard local (-1,1) interval
// to a (0,1) interval. Applied by, e.g., side/surface integration of
// prism elements.
static void reduceInterval (VecSimple(Ptv(real))& pt,
VecSimple(real)& w);
static void combine
(
const VecSimplest(VecSimple(Ptv(real)))& pt,
const VecSimplest(VecSimple(real))& w,
VecSimple(Ptv(real))& combined_points,
VecSimple(real)& combined_weights
);
// function for initializing the rule explicitly:
void refill (int order, int nsd = 1, NumItgDomain domainshape = BOX);
// functions for initializing the rule via a finite element
// (these two refill functions are used by class FiniteElement):
bool refill (ElmDef& elmdef);
bool refill (ElmDef& elmdef, int side); // side integration
// functions for reading the rule:
int getNoPoints () const { return points.size(); }
int getNoSpaceDim () const { return nsd; }
Ptv(real) point (int pt_no);
void point (Ptv(real)& pt_coor, int pt_no);
real weight (int pt_no) const;
// for integration along sides:
static void normalVectorAndDsInBox
(
const Mat(real)& jacobi, // Jacobi matrix
int side, // side of element
real& ds, // computed surface element
Ptv(real)& normal_vector // computed normal vector to surface
);
static void normalVectorAndDsInTriangle
(const Mat(real)& jacobi,int side,real& ds,Ptv(real)& normal_vector);
static void normalVectorAndDsInPrism
(const Mat(real)& jacobi,int side,real& ds,Ptv(real)& normal_vector);
static void normalVectorAndDs
(NumItgDomain domainshape, const Mat(real)& jacobi, int side,
real& ds, Ptv(real)& normal_vector);
void normalVectorAndDs
(const Mat(real)& jacobi, int side, real& ds, Ptv(real)& normal_vector);
private:
void makeRules (VecSimple(Ptv(real))& pt, int side); // called from refill
void redimWhenNsdIsKnown (int nsd); // redim relative_order and dir_order
static void sideVarsInBox
(int side, int nsd, int& i, int& j, int& k, real& icoor);
static void refSide2SideBox
(int side, const VecSimple(Ptv(real))& pt_side, VecSimple(Ptv(real))& pt);
void refSide2SideTriangle
(int side, const VecSimple(Ptv(real))& pt_side, VecSimple(Ptv(real))& pt);
void refSide2SideTetrahedra
(int side, const VecSimple(Ptv(real))& pt_side, VecSimple(Ptv(real))& pt);
void refSide2SidePrism
(int side, const VecSimple(Ptv(real))& pt_side, VecSimple(Ptv(real))& pt);
static void side1Dto2D
(real pt1D, const Ptv(real)& pa, const Ptv(real)& pb, Ptv(real)& pt);
static void redim (VecSimple(Ptv(real))& pt,
int npt,
int ndim);
static void redim (VecSimple(Ptv(real))& pt,
VecSimple(real) & w,
int npt,
int ndim);
};
KEYWORDS
finite element, numerical integration, quadrature
DESCRIPTION
The class performs numerical integration over a finite element,
either the interior or the sides of the element. Usually the
class is accessed internally in a FiniteElement object. The pro
grammer must, however, declare (and, if the default behavior is
not sufficient, also initialize) an ElmItgRules object that can
be used in the FiniteElement::refill function. When using class
FEM, there is an object (ElmItgRules) itg_rules that can be used
- see, e.g., FEM::calcElmMatVec. (An ElmItgRules object is ini
tialized for an element by the refill function, but it is not
necessary for the programmer to refill the ElmItgRules object
prior to each call to FiniteElement::refill, the latter function
will perform a call to ElmItgRules::refill.)
The standard use of the object is hence to declare it, and maybe
use the scan function to adjust the integration point type (Gauss
or nodal points) and the relative order of the rule, and there
after feed the object into FiniteElement::refill.
The type of integration point is given by the enum NumItgPoints,
see the enum man page. Recall that Gauss rules are more accurate
than nodal point rules, but the latter is advantageous when
reproducing finite difference schemes in a finite element con
text. The relative order is a term that allows the user to
define the order of the rule. When used in conjunction with class
FiniteElement, the ElmItgRules object is updated for a new ele
ment by a call to its refill function. That function will check
the current finite element and the order of the basis function
polynomials. That order (plus 1) is the default order of the
integration rule. To this default order, the relative order is
added. Hence, to use default accuracy, dictated by the basis
functions, the relative order should be zero. A value of -1 for
the relative order gives an under integrated element, e.g., the
well known selective reduced integration technique from penalty
methods. A value of +1 or larger for the relative order may be
necessary if the user has products of several basis functions in
the integrands. This may be the case in nonlinear problems or
problems with strongly variable coefficients. Why use relative
instead of absolute order? If the element mesh consists of a mix
ture of elements with polynomials of different degrees, the abso
lute order of the integration rule must vary. By giving a rela
tive order to adjust the default order, one can easily get all
elements under integrated, regardless of whether they may be lin
ear (absolute order 1) or quadratic (absolute order 2).
After declaration, and a possible use of scan or setPointType and
setRelativeOrder, a FiniteElement object performs the element-
wise refilling and computing.
Although the ElmItgRules can almost be hidden for a finite ele
ment programmer, there are several features that are convenient
for the programmer who wants to compose her own rule. This may be
necessary when constructing new elements or when speeding up the
code (a higher order integration rule is for example only
required in some space directions). In this case, one must (1)
avoid that FiniteElement performs an automatic update of the
ElmItgRules object, (2) initialize the ElmItgRules completely
with points and weights. Problem (1) is solved by calling the
ignoreRefill function. Any call to ElmItgRules::refill will then
result in an immediate return such that the composition of the
integration rule is not dictated by the statements in refill.
Recall that FiniteElement::refill calls ElmItgRules::refill.
Action (2) consists in two steps, first the programmer must com
pose a rule, e.g., by applying several of the ElmItgRules fea
tures, and then the points and weights must be "submitted" to the
ElmItgRules object by calling setPointsAndWeights.
The composition of numerical integration rules is normally based
on combining some standard rules for intervals, triangles or
tetrahedras. The ElmItgRules class has functions for computing
the points and weights of such standard rules, see the functions
rulesInterval, rulesTriangle, rulesTetrahedra. These functions
are static and can hence be used anywhere too compute points and
weights of standard integration rules.
Having some standard rules, the function combine combines them
into a composite rule. For example, when integrating over a 3D
box, one can define standard interval rules on (-1,1) for each
space direction and then use combine to combine the three rules
into one rule. The programmer has of course full freedom to use
different number of integration points in the different space
directions if this is desired. As an example, consider a prism
element in 3D with linear triangle basis functions combined with
second order basis functions on (-1,1) in the third direction.
The appropriate integration rule will be a combination of a tri
angle rule and an interval rule. By calling rulesTriangle and
rulesInterval and thereafter combining these two rules using com
bine, the result will be an integration rule for the prism ele
ment.
With the flexibility described above, class ElmItgRules can be
used for integration functions in general, see the Diffpack
report A B-Spline Package in C++ for an example where class
ElmItgRules is applied for integrating B-Spline functions.
CONSTRUCTORS AND INITIALIZATION
There is one constructor and two basic parameters that the user
must initialize. The parameters are the type of integration point
(Gauss points or nodal points) and the relative order of the
rule. These two quantities are described above. For an example
of reduced integration (relative order equal to -1), see the
report The Diffpack Implementation of a Navier-Stokes Solver
Based on the Penalty Function Method.
The (default) values given to the constructor can later be
adjusted by calling scan (using the menu system) or by calling
the functions setPointType and setRelativeOrder.
MEMBER FUNCTIONS
setPointsAndWeights - enables the user to supply her own rule.
See reference to a demonstration program under the examples sec
tion below.
getNoPoints - returns the total number of integration points.
getNoSpaceDim - returns the number of space dimensions of the
points. For example, for surface integrals in 3D, the integra
tion domain is two dimensional but the points are in 3D and hence
the routine returns 3.
point - return an integration point. Recall that the integration
point is given with respect to the local element coordinate sys
tem. The point number is given as argument.
weight - returns the weight associated with an integration point.
The point number is given as argument.
normalVectorAndDs - these functions are aimed at integration over
a side in a finite element. The routines compute the normal vec
tor of the side (at the current integration point) and the
infinitesimal surface area (which corresponds to the Jacobi
determinant in volume integrals).
ignoreRefill - ignores future calls to refill. This is convenient
if the programmer applies a standard finite element assembly
algorithm from class FEM (for example), but want to set the
numerical integration points explicitly. Since the calls to
refill are ignored, refill will not be able to distroy the user
given integration points and weights.
acceptRefill - accepts calls to refill. Can be used after a call
to ignoreRefill (for example).
refill - initializes the integration rules over a domain. One
overloaded version takes the order of the rule and the shape of
the domain as arguments. This function is quite general and can
be used for any domain of hypercube shape or 2D triangles or 3D
tetrahedras. (For triangles and tetrahedras it is equivalent to
calling rulesTriangle or rulesTetrahedra, but for box elements it
computes 1D rules and combines them). Another version takes an
ElmDef object as argument and gains information from this element
object to set the appropriate points and weights for this ele
ment. An overloaded version has the side number as additional
argument. It is aimed at numerical integration over a side in the
element, while the former version of refill was aimed at integra
tion over the interior of the element. refill returns a true
value if the rules in the object have changed (this indicates
recomputation in other classes, e.g., class FiniteElement).
Under the examples section below there is a reference to a demon
stration program where combine is used to create a rule for a
prism element based on a triangle rule and an interval rule.
sideVarsInBox - used for determining the fixed and variable coor
dinates on a side side in a box-shaped element. At return "i is
the coordinate that is fixed (at the value icoor equal to -1 or
1, usually), j and k varies over the side. In 3D problems k=0. In
2D problems, j=0.
rulesInterval - returns the points and weights for Guass-Legendre
or nodal points (approx. Gauss-Lobatto) rules on the interval
(-1,1). Gauss rules up to order 10 are implemented, while nodal
points rules are limited to order 3 (Simpsons rule, which is
applicable to second order elements).
rulesTriangle - returns the points and weights for Gauss or nodal
points rules in a triangle with vertices in (0,0), (1,0) and
(0,1). Gauss rules are supported up to degree 5 (7 points),
while nodal points rules are limited to 2nd order (3 points).
rulesTetrahedra - returns the points and weights for Gauss or
nodal point rules in a tetrahedra with vertices (0,0,0), (1,0,0),
(0,1,0) and (0,0,1). Nodal points rules up to order 2 (4 points)
are supported, while Gauss rules are implemented up to order 3 (5
points).
reduceInterval - reduces the extension of a standard local 1D
(-1,1) interval to an interval defined by (0,1). Integration
points and weights are modified due to the reduction in the
lenght of the interval. Remark that the rulesInterval function
has to be called prior to the present function (because the pre
sent function only modifies already introduced integration points
and weights).
combine - given several rules (with points and weights), this
functions combines the rules into one single rule. For example,
multi-dimensional quadrature based on (-1,1) interval rules in
each space direction can be computed by this routine. The com
bined rule can later be "submitted" to the ElmItgRules object by
calling setPointsAndWeights. Note that when constructing rules
using combine, the points and weights of the subrules must be
computed on beforehand. Neither "combine, nor the rulesInterval,
rulesTriangle, rulesTetrahedra and setPointsAndWeights functions
are influence by the integration point type and relative order
given to the constructor (or set with setPointType and setRela
tiveOrder). These parameters have only application when the
refill(ElmDef& functions are used.
EXAMPLES
Usage of ElmItgRules from a FiniteElement object is best demon
strated by reading the FiniteElement.cpp file. How to create
user defined rules is demonstrated in the program
$DPR/src/app/class-verify/category1/ElmItgRules2/composeIt
gRule.cpp
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO. Some additions by Wen Shen,
SINTEF/UiO. Revised, extended and optimized by Klas Samuelsson,
UiO.