NAME
FiniteElement - all info on a finite element for a programmer
INCLUDE
include "FiniteElement.h"
SYNTAX
class FiniteElement : public HandleId
{
protected:
int elm_no; // current element number
int current_side; // set in refill
int nsd;
BasisFuncAtPt* basis; // basis func, Jacobi at current itg.pt.
BasisFuncAtPt basis_at_one_itg_pt;
VecSimplest(BasisFuncAtPt) basis_at_itg_points; // fixed internal itg pts
int eval_derivatives; // order of deriv. in basis func.
Handle(BasisFuncGrid) bfnodes; // finite element grid for the unknown
HandleElmDefs elmdef; // element definition
Handle(ElmItgRules) itg_rules;// numerical integration points/weights
Mat(real) elmcoor; // coordinates of element nodes
int present_numitg_pt_no;
int no_of_itg_pts;
Ptv(real) present_eval_pt;
Ptv(real) global_eval_pt;
bool global_eval_pt_is_computed;
bool identical_elements; // true: all elms are equal
bool report_detJ_lt0; // error msg when det J <= 0
bool elm_tp_changed; // used in refill/setLocalEvalPts
int basis_pts_init4e; // avoids extra calcJacobi calls
void preprocess (); // if all elements are identical => increased effic.
void init (); // used by constructors
void initNewPoint (bool itg_rules_given_pt); // compute basis object
bool basisOk (const char* function_name) const;
bool redim (int elm_no, FEM* solver = NULL); // called by refill(int,FEM*)
int ncalls2print_since_last_refill; // help variable for the print function
public:
FiniteElement ();
FiniteElement
(
const GridFE& grid,
bool isoparametric_elements = true // only true is legal
);
FiniteElement (const BasisFuncGrid& bfnodes);
~FiniteElement ();
void attach (const GridFE& grid);
void attach (const BasisFuncGrid& bfnodes);
void detach (); // detach all handles
bool refill (int elm_no, FEM* solver = NULL);
bool setLocalEvalPts (ElmItgRules& elmitg,
bool automatic_rule_adjustment = true);
bool refill (int elm_no_, ElmItgRules& elmitg,
bool automatic_rule_adjustment = true)
{
bool b = refill(elm_no_);
setLocalEvalPts(elmitg,automatic_rule_adjustment);
return b;
}
void refill4side (int side, // surface itg, called after refill
bool automatic_rule_adjustment = true);
bool ok () const; // to be called after refill
const GridFE& grid () const { return bfnodes->grid(); }
GridFE& grid () { return bfnodes->grid(); }
const BasisFuncGrid& grid4basisFunc() const { return *bfnodes; }
BasisFuncGrid& grid4basisFunc() { return *bfnodes; }
const ElmDef& getElmDef () const { return *elmdef; }
int getNoBasisFunc () const // length of N/dN
{ return elmdef->getNoBasisFunc(); }
int getNoSpaceDim () const { return nsd; }
int getNoPhysicalSpaceDim () const
{ return bfnodes->getNoPhysicalSpaceDim(); }
int getElmNo () const { return elm_no; }
int getNoSides () const { return elmdef->getNoSides(); }
int getCurrentSide () const { return current_side; } // =0: volume int.
void setLocalEvalPt (const Ptv(real)& loc_pt);
void getLocalEvalPt (Ptv(real)& loc_pt) const;
void getGlobalEvalPt (Ptv(real)& glob_pt) const;
const Ptv(real)& getGlobalEvalPt () const;
int initNumItg (); // init and get no of points
bool moreItgPoints () const; // for while-loop over itg points
void update4nextItgPt (); // for while-loop over itg points
void update4nextItgPt (int itg_pt_no); // for for-loop over itg points
int getNoItgPts () const { return no_of_itg_pts; }
int getCurrentItgPt () const { return present_numitg_pt_no; }
real N (int i) const;
real dN (int i, int dir) const;
real d2N (int i, int dir1, int dir2) const;
const Vec(real)& N () const;
const Mat(real)& dN () const;
real dNloc (int i, int dir) const; // derivatives in local coordinates
real detJxW () const; // detJ times num.itg. weight (at current pt)
real detSideJxW () const; // detSideJ times the integration weight
bool identicalElements () const;
void optimize (bool onoff);
void errorMsg4detJlt0 (bool onoff);
const BasisFuncAtPt* getDataAtNumItgPt (int itg_pt_no);
static real massMat (bool lumped, int i, int j, real N_j);
NUMT massVec (bool lumped, FieldFE& field, int i,
NUMT field_valuePt) const;
bool boSide (int side, int indno) const; // is side subj. to the bo.ind.?
void getNormalVectorOnSide (Ptv(real)& vec) const;
const Mat(real)& JacobiMatrix () const { return basis->JacobiMatrix(); }
const Mat(real)& invJacobiMatrix () const { return basis->invJacobiMatrix();}
const Mat(real)& geomtNodeCoor () const { return elmcoor; }
void print (Os os, int degree_of_output = 1);
// seldom used (only if a special weight is to be applied):
real detJ () const; // Jacobi matrix determinant for volumes
real detSideJ () const; // infinitesimal surfaces element (sides)
real weight (int itg_pt_no);
void evalDerivatives (int i); // compute derivatives of order i
// not used anymore (just for backward compatibility), gives error :
int loc2glob (int loc_node) const;
};
KEYWORDS
finite element, jacobian, isoparametric mapping, global element
DESCRIPTION
The class implements a global finite element consisting of the
definition of the element in its local coordinate system (pro
vided by an object from the ElmDef hierarchy), the basis func
tions and the geometric mapping at an evaluation point (provided
by a BasisFuncAtPt object), the relation between local and global
node numbers (provided by a BasisFuncGrid object) and the nodal
coordinates of the element (given by a GridFE object).
Using an object of class FiniteElement enables calculation of
basis functions, and their derivatives, at a given point in the
element. The main application of the class is hence when evalu
ating integrands of the weak form in finite element problems.
The local evaluation point can be given in two ways: Either
explicitly by the user (in terms of the setLocalEvalPts function)
or automatically according to information on numerical integra
tion rules for the element. The information on the numerical
integration points is provided by an ElmItgRules object.
One updates the FiniteElement object for a new element by calling
the refill function that take an ElmItgRules argument. This
ElmItgRules object must be properly initialized (the integration
point type and the relative order of the rule must be set), but
refilling of that object for the current element is carried out
by FiniteElement::setLocalEvalPts (which is called from the
refill function that takes an ElmItgRules). The integration loop
is most conveniently written as a while loop, see documentation
of class FEM or the report Details of finite element programming
in Diffpack for examples. The initNumItg function is called to
initialize the numerical integration just before the loop starts,
and the while condition tests on the return value of moreNumItg
Points. Inside the loop one calls update4nextItgPt to set the
present evaluation point of FiniteElement to the relevant point
in the quadrature rule and perform the necessary initialization
of the FiniteElement object (such as evaluating the basis func
tions, Jacobi determinant etc).
Instead of having a while loop over the integration points, one
can use a for loop. The initNumItg function returns an integer
which equals the number of integration points in the current ele
ment. This will be the limit of the for loop. Inside the loop,
the update of the integration point is accomplished by calling
update4nextItgPt(k), where k is the current integration point
number. Examples of for and while loops for the numerical inte
gration is given below in the examples section.
Observe that the user supplies an ElmItgRules object with some
rough information on the class of integration rules (Gauss versus
nodal points) and general order of the rule (standard according
to the basis function order or deviation from this). The object
itself will provide more details about the numerical integration
rules when knowing the current element type inside the FiniteEle
ment member functions. In this way, the user can control the
basic principles of the numerical integration rules, like the
rule type, without having to adjust the information when the ele
ment type changes. For example, if the ElmItgRules object says
that Gauss-Legendre rules are to be applied and the relative
order of the rule is -1, the ElmItgRules is able to pick out the
right points even when the element type changes from quadrilat
eral to triangular.
When the grid has elements that differ in size and/or type, the
basis functions, their (global) derivatives and the Jacobian are
recomputed for every integration point. However, if all the ele
ments are of the same type and size, the basis function data will
only vary with the integration point, not the element. One can
hence avoid recomputation of the data in every element. Class
FiniteElement has a feature for doing this. In the general case,
the basis function data is contained in a "BasisFuncAtPt object,
but when the elements are equal, FiniteElement employs a vector
of BasisFuncAtPt objects, where entry no. k contains the basis
function data for integration point no. k. This array of Basis
FuncAtPt objects is computed in "refill(int,ElmItgRules&,int) or
refill(ElmItgRules&,int) and is accessible from the function get
DataAtNumItgPt. See more comments under the documentation of the
identicalElements function.
To initialize a FiniteElement object for a new element, one calls
either refill(int) plus setLocalEvalPt or setLocalEvalPts, or
refill(int,ElmItgRules&,int).
CONSTRUCTORS AND INITIALIZATION
There are three constructors. One requires a GridFE reference
and a bool indicating whether all elements are isoparametric or
not. The only legal value of this bool argument is true, indicat
ing that this constructor can only be used for grids consisting
of isoparametric elements. To initialize the object with data
from an element, the refill member function is used. By default,
the object is initialized with data from element number one.
Another constructor takes a BasisFuncGrid object that should be
filled on beforehand. This constructor can be used for meshes
containing non-isoparametric elements (not yet available).
Finally there is a constructor without arguments. It requires an
additional call to attach in order to attach a GridFE or Basis
FuncGrid object. As usual in Diffpack, calling the constructor
without arguments and thereafter the attach function is com
pletely equivalent to calling the constructor that has the same
argument as the attach function.
To fully initialize a FiniteElement object, one must call refill
for the element in question (see below) and setLocalEvalPt or
update4nextItgPt. The refill function just redimensions internal
data structures and makes some initialization for a new element,
and the setLocalEvalPt or update4nextItgPt sets the internal
evaluation point and evaluates the basis functions, their deriva
tives, the Jacobian etc. at this point.
OPTIMIZATION
When all the elements are of the same type and size, the computa
tions of the basis functions, their derivatives, the Jacobian
etc. can be made particularly efficient by storing the computed
basis function data for each integration point in a table. By
default, class FiniteElement recomputes the basis function infor
mation at every new integration point, and this is a safe, but
not always a very efficient, procedure. The command line argu
ment -bfeopt ON can be used to turn on the optimization where
basis function information is gained from a table look-up. Alter
natively, the programmer can call the function optimize(ON).
Note that the command line option -bfeopt will overrule any pro
gram calls to optimize, but if the option is missing, the argu
ment to optimize determines if the optimization should be
allowed. The default state is to not allow any optimization.
With -bfeopt ON or an optimize(ON) call, class FiniteElement will
- when the grid is attached - check if all the elements are of
the same type and size and apply efficient procedures for comput
ing information related to the basis functions. Be aware of pos
sible computational errors when -bfeopt ON or optimize(ON) is
used and the grid coordinates are manipulated by the user. The
manipulation may turn a uniform mesh into a non-uniform mesh, and
if no the function GridFE::setNonUniformMesh is called, the grid
may still think it is uniform and pass this information to class
FiniteElement. When deforming or moving the grid one should
always use the GridFE::move function (and not GridFE::putCoor)
which ensures safety with respect to updating the grid's status.
MEMBER FUNCTIONS
Several member functions concern quantities related to the "Jaco
bian" of the mapping between local and global coordinates. Our
term "Jacobi matrix" refers actually to the transposed Jacobi
matrix as it is defined in standard mathematical literature. How
ever, in finite element literature it is common to use the term
"Jacobi matrix" as we do it here.
refill - initializes the object by using data from the element
number given as argument to this function. The refill function
that takes one argument, the element number, must be followed by
a call to setLocalEvalPt or setLocalEvalPts in order to enforce
evaluation of the basis functions, the Jacobian etc. at a point.
The refill functions that also take an ElmItgRules reference as
argument enables the FiniteElement object to fill its data struc
ture with information on the integration (evaluation) points.
Instead of calling setLocalEvalPt, the FiniteElement object then
applies the ElmItgRules object for setting the local points
through the function setLocaEvalPts. To enforce calculation of
basis functions and the associated quantities, one usually calls
update4nextItgPt inside an integration loop. Of course, one can
also call setLocalEvalPt (the supplied ElmItgRules object is then
overrid). The user can just provide an ElmItgRules object (e.g.
from class FEM), the FiniteElement::refill function will update
the ElmItgRules object when necessary if the third argument auto
matic_rule_adjustment is true (the ElmItgRules object must only
be initialized). The updating of the ElmItgRules object is based
on information about the element shape, order of the basis func
tions and the so called relative order of the rule (set in ElmIt
gRules). A false value of third argument implies that the ElmIt
gRules object is already completely initialized by a special rule
and this rule should be used in the computations. Normally, one
sticks to the default true value for the third argument, but in
cases where the user want to have full control of the integration
rule, the third argument gives the necessary flexibility. The
refill(int, ElmItgRules&, int) function is actually just a call
to refill(int) and setLocalEvalPts. The reason for splitting the
initialization up like this, is that some special cases need cus
tomization of, e.g., the ElmDef subclass object after refill(e),
but before the array of basis functions etc. are computed in set
LocalEvalPts.
The refill4side function can be called after refill(ElmItgRules&)
and each update of the evaluation points (update4nextItgPt) then
loads a numerical integration point on the side of the element.
The side number is specified as an argument to refill4side. The
other argument is explained for the other refill function.
getNoBasisFunc - returns the number of basis functions in the
elements. For isoparametric finite elements this number coin
cides with the number of nodes in the element. Observe that the
function is not inline. In addition, the function performs a test
for safety reasons that the FiniteElement object is properly ini
tialized. For these reasons, one should never call getNoBasisFunc
inside loops (or as a part of a for-statement); extract the value
and store it in an int outside the loop.
getNoSpaceDim - returns the number of space dimensions of the
element.
getElmNo - returns the current element number for which this
finite element is initialized.
getNoSides - returns the number of sides in the current element.
getCurrentSide - returns the current side number. The number is
set in the all to refill4side. If a zero value is returned, it
means that no side is set, that is, the finite element object is
not evaluated at a numerical integration point on a side (usually
it means that the finite element object is performing a volume
integration over the interior of the element).
grid - returns access to the underlying grid.
grid4basisFunc - returns access to the BasisFuncGrid object. For
isoparametric elements one can instead use grid (which will be
slightly more efficient since the functions in BasisFuncGrid in
the isoparametric case first make a test and then call the simi
lar GridFE function).
getElmDef - returns access to the element definition (as a base
class reference).
setLocalEvalPt - sets the coordinates, in the local element coor
dinate system, of the point at which basis functions are to be
evaluated.
getLocalEvalPt - extracts the current evaluation point in the
element, expressed in the local element coordinate system. The
current evaluation point is usually automatically set as a part
of the numerical integration facilities, but can also be explic
itly set using setLocalEvalPt.
getGlobalEvalPt - finds the corresponding global coordinates of
the current (integration) point. In other words, getGlobalEvalPt
maps the result of getLocalEvalPt to global/physical coordinates.
initNumItg, moreNumItgPoints, update4nextItgPt - these functions
are explained above.
N - returns the value of basis function no. i at the present
evaluation point.
dN - returns the value of the derivative, in global (physical)
coordinates, in dir direction of basis function no. i.
dNloc - returns the value of the derivative, in local coordinates
(reference element coordinates), in dir direction of basis func
tion no. i.
d2N - returns the value of the 2nd order derivative in dir1 and
"dir2 directions of basis function no. i. The functionality is
only available if the evalDerivatives(i) is called with i equal
to 2. This must be done in the beginning of the program (the
evalDerivatives function is static so no FiniteElement object is
needed).
detJxW - returns the Jacobi determinant times the numerical inte
gration weight at the current evaluation point (as computed by a
call to update4nextItgPt or explicitly set by setLocalEvalPt).
This function is not as efficient as detJ (because it performs
some consistency checks) and it should be called once for each
evaluation point and the value should be stored in a local vari
able in the user's program to maximize efficiency. If one wants
to compute the detJxW at an integration point with a prescribed
number k (typically the case if one codes the element-loop as the
innermost loop), one should call
fe.getDataAtNumItgPt(k).detJ()*fe.weight(k)
This requires that the getDataAtNumItgPt function really works,
that is, that the grid consists of elements of the same type and
size. When compiled in non-optimized mode, the user will get a
run-time message if getDataAtNumItgPt does not work (in optimized
mode one will dereference a NULL pointer which probably causes a
core dump).
detSideJxW - the surface integral factor multiplied by the weight
for surface integration. Take a copy of the value and use that
copy for maximal efficiency when evaluating the various contribu
tions in the element matrix/vector.
identicalElements - returns true if the FiniteElement object con
siders all elements to be of the same type and shape (the prepro
cess function called at construction time is used to determine if
all the elements are identical by asking the GridFE functions
hasUniformMesh and allElmsOfSameType). In the case where all the
elements are of the same type and shape, one can compute the
basis functions, their (global) derivatives and the Jacobian only
once at an integration point. In other words, the BasisFuncAtPt
information is independent of the element number, and one can
have an array of BasisFuncAtPt objects, with length equal to the
number of (volume) integration points in the element, that holds
all the necessary information on basis functions which is needed
in the finite element algorithms. The FiniteElement object will
automatically .. If the weak form does not contain complicated
expressions, Only initial computation of the basis functions, the
Jacobian etc, will speed up the assembly process (provided that
evaluation of the weak form in integrands routines does not over
shadow the work associated with computation of a BasisFuncAtPt
object for a new integration point).
getDataAtNumItgPt - returns a BasisFuncAtPt object with informa
tion about basis functions, their derivatives, the Jacobian etc.
at an internal integration point in the element. The function is
only useful when all elements are of the same type and size (if
not, an error message will be issued). The user can then access
the basis functions at the integration points in any order. The
standard use of FiniteElement is sequential in the sense that one
assumes that the "integration loop" is outside the "basis func
tion loop". With getDataAtNumItgPt one can alter the sequence of
the loops and e.g. have the loop over the elements as the inner
most loop (this is advantageous for vectorization). The function
is not inline (and performs a test for safety reasons) so one
must avoid calling the function inside long loops. Instead,
declare a const BasisFuncAtPt& reference and store the returned
value in this local reference outside the loop.
massMat - utility function for computing the mass matrix term in
the integrand of a weak form. The boolean argument can be set to
achieve a lumped (true value) or consistent mass matrix. Note
that the function is static and can hence be used outside the
standard Diffpack finite element framework. The final argument
N_j should be the value of basis function no. j, typically
fe.N(j) if fe is a FiniteElement object. The N_j argument is only
used in case of consistent mass so if the aim is to compute a
lumped mass matrix, N_j can be any value.
massVec - companion function to massVec. massVec computes the
contribution of the mass matrix on the right hand side. Again
both the lumped and consistent case can be handled by specifying
the boolean argument. The form of the term is Mu, where M is the
mass matrix and u is some field, represented as a FieldFE object.
The FieldFE object is only used in case of a consistent mass
matrix. The final argument field_valuePt is the evaluation of u
at the current integration point (see the function for comments
on the need for this). This argument is only used in case of a
lumped mass matrix.
print - prints the contents of the object, mostly needed for edu
cational or debugging purposes.
JacobiMatrix - returns the Jacobi matrix (note that this is the
transposed Jacobi matrix according to the definition of the
Jacobi matrix in mathematical literature).
invJacobiMatrix - returns access to the inverse Jacobi matrix.
geomtNodeCoor - returns access to all the coordinates of the
nodes defining the geometry of the element. This function may be
particularly useful when finite element data structures from
Diffpack are to be sent to C or Fortran functions (simply extract
the underlying array by using getPtr0 or getData).
detJ - returns the determinant of the Jacobian of the isoparamet
ric mapping. It is intended for applications where the user wants
to adjust the weights in an integration rule and hence does not
want the determinant mixed with the weight as in detJxW.
detSideJ - the "Jacobi determinant" on a side, or in other words,
the surface integral factor for surface integration. See the
comments in the detJ function documentation.
loc2glob - returns the global node number corresponding to
the given local node number in the present element (note that the
present element is changed by the refill function). The node
numbers correspond to those used for the basis functions.
getNormalVectorOnSide - returns the unit outward normal vector at
the evaluation point. The function is only meaningful when the
local evaluation is at a side. The computation of the normal vec
tor is performed by the refill4side function. Hence the getNor
malVectorOnSide function is restricted to the cases where the
evaluation point is automatically computed as an integration
point.
BoSide - returns a true value if the given side is subjected to
the given boundary indicator. A side is subjected to an indicator
if all nodes on the side are subjected to the indicator.
evalDerivatives - used to set the order of the derivatives that
should be computed. Order 1 is default, order 0 has not been
tested, order 2 implies additional computations in class Basis
FuncAtPt for second order derivatives, and this option should
only be used when strictly necessary.
errorMsg4NegativeJacobian - enables setting of a flag for con
trolling the writing of error messages when the determinant of
the Jacobian of the mapping between local and global coordinates
is less than or equal to zero. By default, the update4nextItgPt
function tests if the determinant of the Jacobian is non-positive
and reports an error. In problems with moving grids, it can be
advantagous to turn off these error messages and instead check
the determinant in the simulator. If the determinant is non-pos
itive, a new re-meshing procedure should be invoked. By calling
errorMsg4detJlt0(OFF) the FiniteElement object will avoid error
messages due to a non-positive determinant of the Jacobian. (The
BasisFuncAtPt class also tests on the determinant of the Jacobian
when the parameter to the command line option +verbose is 1 or
greater.)
EXAMPLES
Given a FiniteElement object fe and a ElmItgRules object rules,
here are the typical steps in a numerical integration over an
element e. First a while loop.
fe.refill (e, rules);
fe.initNumItg ();
while (fe.moreNumItgPoints()) {
fe.update4nextItgPt();
// fe.N(i), fe.dN(i,j), fe.detJxW etc are ready for use!
// fe.print (s_o, 3); // check information!
}
This can also be coded as a for loop:
fe.refill (e, rules);
int npt = fe.initNumItg ();
for (int k = 1; k <= npt; k++) {
fe.update4nextItgPt(k);
// fe.N(i), fe.dN(i,j), fe.detJxW etc are ready for use!
// fe.print (s_o, 3); // check information!
}
It is not necessary to call initNumItg prior to every such loop.
If you know the number of integration points, only update4nextIt
gPt is critical for initialization of the fe object at a point.
Here is an example of how to make a table of all the basis func
tions and their derivatives in all the integration points in an
element.
// given an ElmItgRules itg_rules
FiniteElement fe (grid, true);
! fe.optimize (ON); // allow optimization, i.e., a table of all
// basis functions at all itg. points in an
elm.
fe.refill (1, itg_rules); // the table is automatically made
use getDataAtNumItgPt(k) to recover basis function data at inte
gration point no. k, check the man page of BasisFuncAtPt for how
to extract the various basis function information.
SEE ALSO
class BasisFuncAtPt
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO