NAME
ElmDef - abstract base class for definition of finite elements
(local system)
INCLUDE
include "ElmDef.h"
SYNTAX
enum RectanglePartition // rectangle vs. its two (upper/lower) triangles
{ // should be local. global for C++ v2.1/3.0 mix
RECTANGLE, // rectangle
LOWER_TRIANGLE, // lower triangle
UPPER_TRIANGLE // upper triangle
};
class ElmDef : public virtual HandleId
{
protected:
int nsd; // no of space dimensions
int nne_geomt; // no nodes defining element geometry functions
int nne_basis; // no nodes defining basis (trial) functions
int nmodes; // no of internal modes of basis functions (hierarchical)
int nbf; // no of basis functions (!=nne_basis for hierarchical func.)
int nsides; // no of sides in the element
Ptv(int) order; // order of interpolation polynomial (in each space dir)
bool isoparametric_elm; // true if it is an isoparametric element
NumItgDomain domainshape; // elm shape: BOX/TRIANGLE (for ElmItgRules)
ElementType elmtype; // enum variable for identification of elements
MatSimple(real) geomt_coor; // coordinates of the local geomt. nodes
MatSimple(real) basis_coor; // coordinates of the local basis func. nodes
Ptv(real) loc_centroid; // the centroid in local coordinates
Vec(real) workgeomt; // work array, used in mapping
Vec(real) workbasis; // work array, used in findLocPt
MatSimple(int) geomt_nodes_sides; // geomt. nodes on a side
MatSimple(int) basis_nodes_sides; // basis func. nodes on a side
public:
int ngno_on_a_side; // no of geometry nodes on a side (max if variable)
int nbno_on_a_side; // no of basis func. nodes on a side (max if variable)
protected:
ElmDef // the most general constructor
(
int nne_geomt,
int nne_basis,
int nmodes,
int nbf,
int nsd,
int nsides,
int max_ngno_on_a_side, // for geometry nodes
int max_nbno_on_a_side, // for basis func. nodes
const Ptv(int)& order, // order of basis func polynomial
NumItgDomain domainshape,
ElementType elmtype
);
ElmDef // constructor for standard isoparametric elements
(
int nne,
int nsd,
int nsides,
int nno_on_a_side,
int order,
NumItgDomain domainshape,
ElementType elmtype
);
public:
virtual ~ElmDef();
public:
int getNoSpaceDim () const { return nsd; }
int getNoGeomtFunc () const { return nne_geomt; }
int getNoBasisFunc () const { return nbf; }
int getNoBasisNodes () const { return nne_basis; }
// (it is assumed that no of geomt func equals no of geomt nodes)
int getNoSides () const { return nsides; }
Ptv(int) getOrder () const { return order; }
NumItgDomain getDomainShape () const { return domainshape; }
bool isoparametric () const { return isoparametric_elm; }
ElementType getElementType () const { return elmtype; }
virtual NumItgDomain getSideShape (int side) const =0;
// geomtFunc: functions for the geometric mapping of the element.
// basisFunc: functions for interpolating unknowns over the element.
// For isoparametric elements geomtFunc and basisFunc are equal
virtual void geomtFunc (Vec(real)& N, const Ptv(real)& loc_pt) const =0;
virtual void basisFunc (Vec(real)& N, const Ptv(real)& loc_pt) const =0;
// derivatives of N with respect to the local element coordinates:
virtual void dLocGeomtFunc (Mat(real)& dNloc, const Ptv(real)& loc_pt)
const=0;
virtual void dLocBasisFunc (Mat(real)& dNloc, const Ptv(real)& loc_pt)
const=0;
// the default version of this return zero 2nd order derivatives:
virtual void d2LocBasisFunc (ArrayGen(real)& d2Nloc, const Ptv(real)& loc_pt)
const;
void mapping
(
Ptv(real)& glob_pt,
const Ptv(real)& loc_pt,
const Mat(real)& coor
)const;
void mapping
(
Ptv(real)& glob_pt,
const Mat(real)& coor,
const Vec(real)& N
)const;
Ptv(real) mapping (const Ptv(real)& loc_pt, const Mat(real)& coor) const;
void consistencyCheck (); // is the element ok?
virtual bool isInside (const Ptv(real)& locpt) const;
virtual bool isStrictlyInside (const Ptv(real)& locpt) const;
virtual bool canBeInside
(const Ptv(real)& globpt, const Mat(real)& coor) const;
bool isOnElmBoundary (const Ptv(real)& locpt) const;
bool canBeInsideBox // version for curved element sides
(const Ptv(real)& globpt, const Mat(real)& coor, bool curved) const;
protected:
virtual Ptv(real) setCentroid () const; // only used internally
public:
Ptv(real) centroid () const { return loc_centroid; }
virtual bool findLocPt // find local pt corresponding to global pt
(
const Ptv(real)& globpt, // global pt
Ptv(real)& locpt, // local pt (output)
int& niter, // no of Newton iterations in algorithm
real& error, // error from the Newton iteration
const Mat(real)& coor // global coordinates of the nodes
);
virtual int getNoGeomtNodesOnSide (int /*side*/) const
{ return ngno_on_a_side;
/* NOTE: redefine this function if nno_on_a_side varies! */ }
virtual int getNoBasisNodesOnSide (int /*side*/) const
{ return nbno_on_a_side; }
int getNodeOnSide (int side, int local_node_on_side) const
{ return geomt_nodes_sides (side, local_node_on_side); }
int getBasisNodeOnSide (int side, int local_node_on_side) const
{ return basis_nodes_sides (side, local_node_on_side); }
virtual int getNoDof4scalarFunction (int /*local_node*/) const
{ return 1; }
void getNodeOnSide (int side, VecSimple(int)& nodes) const; // inefficient!
// some preprocessors need the number of nodes along an edge in order
// to compute the number of elements from the number of divisions:
virtual int getNoNodesOnEdge () const; // not defined in 1D and 2D
// used to indicate that the ElmDef subclass object must communicate with
// the application (PDE):
virtual bool applicationDepPrms () { return false; }
// --- drawing routines ---
virtual void drawElement
(
const Mat(real)& coor, // global, nodal coordinates
GraphBasics& plotter, // graphics: data, functions, dvi-code
int resolution = 1 // no of line segments along an edge
)const =0;
virtual void drawSide
(
int side, // side number
const Mat(real)& coor, // global, nodal coordinates
GraphBasics& plotter, // graphics: data, functions, dvi-code
int resolution = 1 // no of line segments along an edge
)const =0;
// draw edge (3D elements) between two sides, side1 and side2:
virtual void drawEdge
(
int side1, // first side
int side2, // second side
const Mat(real)& coor, // global, nodal coordinates
GraphBasics& plotter, // graphics: data, functions, dvi-code
int resolution = 1 // no of line segments along an edge
)const;
void setResolution2DGrid (int nintervals)
{ resol2Dgrid = nintervals; }
// the default virtual functions here assume that the 2D grid is
// defined on a square shaped domain (-1,1)x(-1,1).
virtual Ptv(real) getVertexPt2DGrid (int i, int j) const // 2D elements
{ return getPt2DGridSqr (double(i), double(j), resol2Dgrid+1); }
virtual Ptv(real) getCenterPt2DGrid
(int i, int j, RectanglePartition ind) const // 2D elements
{ return getCenterPt2DGridSqr (i, j, ind); }
virtual Ptv(real) getVertexPt2DGrid (int i, int j, int side) const // 3D
{ return getPt2DGridSqr (double(i), double(j), resol2Dgrid+1, side); }
virtual Ptv(real) getCenterPt2DGrid
(int i, int j, RectanglePartition ind, int side) const // 3D
{ return getCenterPt2DGridSqr (i, j, ind, side); }
virtual void calcMidSideNodes (Mat(real)& /*coor*/)
{ errorFP("ElmDef::calcMidSideNodes",
"not implemented for class"); }
// the below procedure is to be used from the virtual void calcMidSideNodes.
// however, in the virtual functions a call to the below routine will
// not be resolved by overloading rules - the compiler treats the call
// as a recursive call and reports a type mismatch error in the arguments.
// the below routine should be used in the virtual functions as
// ElmDef::calcMidSideNodes (...)
void calcMidSideNodes
(
Mat(real)& coor, // nodal coordinates
ElmDef* interpol_elm, // interpolation element
const Mat(real)& coor_interpol_elm // nodal coord. of interp.elm.
);
void print (Os os) const;
protected:
// these check functions are used by geomtFunc, basisFunc, dLocGeomtFunc etc:
bool checkNandPt (Vec(real)& N , const Ptv(real)& pt,bool geomt)const;
bool checkdNandPt (Mat(real)& dN, const Ptv(real)& pt,bool geomt)const;
void drawBoxElmEdge
(
int side1, // first side
int side2, // second side
const Mat(real)& coor, // global, nodal coordinates
GraphBasics& plotter, // graphics: data, functions, dvi-code
int resolution = 1 // no of line segments along an edge
)const;
void drawCurvedLine
(
int vc, // variable coordinate
const Ptv(real)& from, // local start coordinate
const Ptv(real)& to, // local stop coordinate
const Mat(real)& coor, // global nodal coordinates
int resolution, // no of straight line segments
GraphBasics& plotter // graphics, output tools
)const;
void drawBoxElmSide
(
int side, // side number
const Mat(real)& coor, // global, nodal coordinates
GraphBasics& plotter, // graphics, output tools
int resolution // no of line segments along an edge
)const;
void drawAllSides
(
const Mat(real)& coor, // global, nodal coordinates
GraphBasics& plotter, // graphics, output tools
int resolution // no of line segments along an edge
)const;
// for side = 1,...,nsides: drawSide
// ---------------- routines for drawing contour and colour plots:
int resol2Dgrid; // no of divisions of an element for contour/colour plot
// square grid (on (-1,1)x(-1,1):
Ptv(real) getPt2DGridSqr (double i, double j, int limit) const;
Ptv(real) getPt2DGridSqr (double i, double j, int limit, int side) const;
// triangle shaped domain, grid consists of squares and triangles
// for 0<x<1, 0<y<1 and y<x.
Ptv(real) getPt2DGridTri (double i, double j, int limit) const;
Ptv(real) getPt2DGridTri (double i, double j, int limit, int side) const;
// getCenterPt2DGrid depends on whether the domain/side is of
// square or trangle shape:
Ptv(real) getCenterPt2DGridSqr
(int i, int j, RectanglePartition ind) const;
Ptv(real) getCenterPt2DGridSqr
(int i, int j, RectanglePartition ind, int side) const;
Ptv(real) getCenterPt2DGridTri
(int i, int j, RectanglePartition ind) const;
Ptv(real) getCenterPt2DGridTri
(int i, int j, RectanglePartition ind, int side) const;
public:
Ptv(real) getLocalGeomtCoor (int local_node) const
{ Ptv(real) lc(nsd); getRow(local_node,lc,geomt_coor); return lc; }
Ptv(real) getLocalBasisCoor (int local_node) const
{ Ptv(real) lc(nsd); getRow(local_node,lc,basis_coor); return lc; }
private:
void compatibilityCheck (const Mat(real)& coor) const;
// do not define virtual casts for this hierarchy - the standard virtual
// functions and functions in the base should be sufficient for
// operations of all elements.
// CLASS_INFO is neither necessary in the derived classes.
public:
virtual bool basisFuncNodes(MatSimple(real)& loc_coor);
CLASS_INFO
};
KEYWORDS
finite elements, base class
DESCRIPTION
The class contains data structures and functions for describing a
finite element in its local element coordinate system. Particu
lar elements are implemented as derived classes of class ElmDef.
A programmer should normally not operate on element subclasses
using an ElmDef handle (or pointer/reference). Instead the pro
grammer should apply the HandleElmDefs handle which has addi
tional features for easy updating the current element given its
name and the number of space dimensions.
Each finite element is of a particular type, for example bilinear
element, triquadratic element and so on. In general the geometric
shape of the element is determined by the mapping from a local
coordinate system onto the physical coordinate system. The func
tions that are used in that mapping are here called element
geometry functions. In addition to these functions one must
specify the trial functions, that is, the basis functions that
are used to approximate the unknown functions in the problem over
the element. The basis and geometry functions may in principle be
different (that is the case in mixed finite elements methods).
However, in many cases they are identical and we refer to the
element as isoparametric.
Given a grid, it may be useful to approximate different unknowns
using different basis functions over an element with fixed shape.
For example, in fluid flow applications the elements may be
quadrilaterals with bilinear interpolation for velocities and
constants for the pressure. In that case, the geometry functions
are the standard bilinear functions while the basis functions are
a constant (for pressure) and the bilinear functions (for veloci
ties). Using ElmDef objects, the velocities would be non-isopara
metric elements with the basis function nodes on the middle of
each side, while the pressure requires a special element where
the geometry functions are bilinear and the basis function is
constant. Such elements can be found in the ElmDef hierarchy
(see, e.g., ElmB4gn1bn2D for an example on piecewise constant
basis functions and bilinear geometry).
Objects in the ElmDef hierarchy are used several places in finite
element algorithms (for example, in FiniteElement and FieldFE).
The typical use is when, e.g., basis functions are required for a
new element. One checks the element type in the grid, and sends
the type (as a String or ElementType) to the Han
dleElmDefs::refill function which then creates a subclass object
in the ElmDef hierarchy.
NOTE: Some comments in the documentation of this class refer to
mixed finite element methods. Support for such methods is only to
a very limited extent provided by the present version of the
class.
CONSTRUCTORS AND INITIALIZATION
This base class is initialized completely by its only constructor
(which is protected and only called from the subclass objects).
Derived classes must in addition initialize the arrays geomt_coor
and basis_coor containing the coordinates of the local geometry
function nodes and the local basis function nodes (in the local
element coordinate system). Similar arrays (geomt_nodes_sides and
basis_nodes_sides) containing the local node numbers that belong
to the different sides of the element must also be initialized.
This type of initialization is performed inside the constructor
of the derived classes and should be of no concern for a user of
already programmed classes.
MEMBER FUNCTIONS
Only the public member functions are described below. For most
users, the geometry and basis (trial) functions are equal
(isoparametric elements) and the description of geomtFunc and
basisFunc is unnecessarily complicated. It then suffices to know
that these two routines calculate the local "shape functions" on
the elements.
geomtFunc - calculates the geometry functions on the element (see
the description above). These functions are used for the isopara
metric mapping associated with the element (these functions
determine the geometric shape of the element in the global, phys
ical coordinate system). Recall that the geometry functions are
defined with respect to the local element coordinate system and
the local node numbering of the element. The return array N can
be empty (not ok) at input, the geomtFunc will redimension the
array properly. The N array contains the values of the geometry
functions at the evaluation point. The length of the N array
equals the number of (geometry) nodes in the element, and the
ordering of the entries in N follow the ordering of the local
nodes (cf. concrete implementations and conventions in the sub
classes).
basisFunc - calculates the basis (trial) functions over the ele
ment. In other words, the routine computes the interpolation
functions for functions defined over the element. Functions
entering differential equations etc should be interpolated
according to these basis functions. The return array N can be
empty (not ok) at input, the basisFunc will redimension the array
properly. The length of N equals the number of basis function
nodes in the element. For isoparametric elements basisFunc and
geomtFunc are identical. However, for mixed finite elements, one
may have a bilinear geometry function (four nodes) and a constant
basis function. The basis function has then only one node (for
example in the centroid of the element) and the length of N will
be 1. Confer subclasses for information on conventions regarding
basis function nodes and the ordering. It should be mentioned
that class BasisFuncGrid is used to assign basis nodes to a grid
(class GridFE is meant to contain the geometry of the grid,
including geometry nodes).
dLocGeomtFunc - the derivatives in the local element coordinates
of the functions defined by geomtFunc. As for the geomtFunc and
basisFunc functions the return array can be empty or have wrong
dimension at input. The function will redimension the array prop
erly.
dLocBasisFunc - the derivatives in the local element coordinates
of the functions defined by basisFunc. As for the geomtFunc and
basisFunc functions the return array can be empty or have wrong
dimension at input. The function will redimension the array prop
erly.
d2LocBasisFunc - the second order derivatives in the local ele
ment coordinates of the functions defined by basisFunc.
mapping - maps a point in the local element coordinate system to
its corresponding point in the global, physical coordinate sys
tem. The function makes use of the virtual geomtFunc member
function.
consistencyCheck - checks if the element is in an ok state. No
value is returned. If abnormalities are detected, error messages
are issued. This function is mainly used for checking the imple
mentations of new elements.
getNoSpaceDim - returns the number of physical space dimensions
of the element.
getNoGeomtFunc - returns the number of geometry functions.
getNoBasisFunc - returns the number of basis (trial) functions on
the element.
getNoSides - returns the number of sides of the element.
getOrder - returns the polynomial order of the basis functions on
the element.
getDomainShape - returns an enum that reflects the shape of the
element, BOX versus TRIANGLE. This is necessary for, among other
things, numerical integration over the element.
getSideShape - returns the shape of the sides of the element.
This is the counterpart to getDomainShape for example for numeri
cal integration over the sides.
isoparametric - returns a true value if the element is of
isoparametric type.
getElementType - returns the variable ElementType (int) which
informs about the element type the object ElmDef is filled with.
All possible values of ElementType is listed in ElmDef_prm.h.
Each implemented version of ElmDef has a unique value of the Ele
mentType.
canBeInside - returns a true value if a point in the global
(physical) coordinate system can be inside the element. The algo
rithm for detecting whether the point is inside or outside is
only approximative. To achieve an exact algorithm, one can use
findLocPt and then test the local point by using the function
isInside.
isInside - returns a true value if the local coordinates of a
point implies that the point is inside the element. The function
is mainly used to check if findLocPt has found a point that is
really inside the element.
isStrictlyInside - returns a true value if the local coordinates
of a point implies that the point is -strictly- inside the ele
ment, i.e. inside but not on the boundary.
isOnElmBoundary - returns true if the local coordinates is inside
but not strictly inside. The function uses isInside and isStrict
lyInside.
centroid - returns the centroid of the element in local coordi
nates.
getNoDof4scalarFunction - returns the number of degrees of free
dom at a nodal point that the element needs to represent a scalar
finite element function. For all elements where the nodal values
are used in the definition of the function, there is one degree
of freedom per node. However, special elements, like the 1D Her
mite cubic element, has two degrees of freedom per node (the
function value and the derivative). Elements in the p version of
the finite element method typically have many sidemodes, which in
the Diffpack terminology translates to many degrees of freedom
per side node (and correspondingly many degrees of freedom per
internal node - one internal node is associated with all the bub
ble functions). For an ElmDef subclass object one should have the
following relation: The sum of getNoDof4scalarFunction over all
nne_basis nodes must equal nbf. This can be used to implicitly
define the concept "basis function nodes" for a given element.
findLocPt - given a global point, the function solves the inverse
mapping equations to find the corresponding local coordinates.
Newton-Raphson iteration is used for the solution. The function
is mainly used in the findElmAndLocPt functions in class GridFE.
The next three functions draw parts of an element (the volume, a
side, an edge) and are used in drawing routines in class DrawFE.
drawElement - function for drawing an element in the global,
physical coordinate system. The drawing consists of geometric
primitives (line segments) in the GraphBasics format. Note that
this function does not flush or close the GraphBasics file so the
physical file may be incomplete after a call to this function.
An explicit call to \mp{GraphBasics::close} will complete the
file and make it ready for reading.
drawSide - function for drawing a specified side of an element in
the global, physical coordinate system. See drawElement for
potential problems with incomplete output files.
drawEdge - function for drawing the intersection of two specified
sides of an element in the global, physical coordinate system.
See drawElement for potential problems with incomplete output
files.
The next set of routines are used for dividing sides of the ele
ments (or the whole element in 2D) into triangles or quadrilater
als. This 2D grid on the element is used for contour-plot rou
tines or colour-fill routines in class DrawFE.
setResolution2DGrid - sets the resolution of grids for drawing
options. Over a 2D element or over a side in a 3D element one
can define a uniformly partitioned grid that can be used for var
ious drawing procedures. This function sets the number of grid
intervals in each direction of the grid. Note that the grid is
always two-dimensional. The grid can either be boxes or trian
gles. For accessing the grid points or the centers in the
boxes/triangles, see the functions "getVertexPt2DGrid and getCen
terPt2DGrid.
getVertexPt2DGrid - returns a Ptv(real) containing the coordi
nates of of vertex point (i,j) in the 2D grid over a 2D element
or over a side in a 3D element (see function setResolu
tion2DGrid). A vertex point is a corner of the boxes/triangles of
the 2D grid. The lower, leftmost vertex has coordinates (1,1).
getCenterPt2DGrid - as getVertexPt2DGrid but the center point of
the box or triangle with index (i,j) is returned. The lower,
leftmost box has index (1,1). If triangles are desired, the user
must indicate whether the center of the lower or the upper trian
gle is desired. The enum "RectanglePartition takes the names
LOWER, UPPER or RECTANGLE indicating lower triangle, upper trian
gle or four-sided box, respectively.
The next set of routines is related to nodes on sides in ele
ments. Applications are in surface integrals, drawing routines
and preprocessor classes. The user will seldom call these rou
tines directly, but it is necessary to supply the routines when
implementing new elements in the ElmDef hierarchy.
calcMidSideNodes - if mid-side nodes are missing, calculate their
coordinates by an interpolation scheme that is suited for the
element.
getNoGeomtNodesOnSide - returns the number of geometry nodes on a
specified side of the element.
getNoBasisNodesOnSide - returns the number of basis (trial) func
tion nodes on a specified side of the element.
getNodeOnSide - given a side and a local node on this side (the
local node is numbered from 1 to the number of nodes on a side)
the function returns the cooresponding node number in the local
element system. For example, ElmB4n2D::getNodeOnSide(2,1) equals
3. An overloaded version takes a VecSimple(int) argument and
fills it with all the node numbers that belong to the specified
side. This version takes a redim of the vector and may hence
require memory allocation.
getBasisNodeOnSide - returns the number of basis function nodes
on a side in the element.
getNoNodesOnEdge - returns the number of nodes on an edge (an
edge is the intersection between two adjacent sides in a 3D ele
ment). The function has no meaning for 1D and 2D elements.
print - prints the contents of the element (number of nodes,
coordinates of local nodes etc).
getLocalGeomtCoor - returns the local coordinates of the geometry
nodes.
getLocalBasisCoor - returns the local coordinates of the basis
function nodes.
basisFuncNodes - returns an array of the basis function nodal
coordinates.
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO. Some modifications by Wen
Shen, SINTEF/UiO, Hilde N. Lund, UiO, and extensive modifications
by Klas Samuelsson, UiO.