NAME
DegFreeFE - degree of freedom manager for finite element linear
systems
INCLUDE
include "DegFreeFE.h"
SYNTAX
class DegFreeFE : public HandleId
{
public:
// public for use in makeSparsityPattern functions:
Handle(FieldsFE) fields; // basic info for the general numbering
Handle(GridFE) mesh; // basic info for the special numbering
private:
// used for eqs numbering, if NULL then # eqs = # dof:
Handle(BasisFuncGrid) test_func;
int ntotdof; // total no of dof in the system
int ntoteqs; // total no of equations in the system
int nelmdof; // no of total dof per element
int nelmeqs; // no of total equations per element
int nno_grid; // used to detect if no of unknowns changes
int ndof_per_node; // no of dofs per node (if mesh.ok)
VecSimple(NUMT) ebc; // hold ess. b.c. for each degree of freedom
// these arrays are used only for the general numbering functionality:
MatSimple(int) loc2glob_arr; // loc->glob numbering (general numb.alg.)
MatSimple(int) u2x; // help array for field2vec
MatSimple(int) x2u; // help array for vec2field
bool interpret_ess_bc_as_zero; // needed in e.g. Newton iter
bool no_ess_bc; // true if no conditions are set
// indicators for various implementations of essential boundary conditions:
bool modify_mat_due2essBC;
bool modify_vec_due2essBC;
bool symm_mod_due2essBC;
bool unit_mat_diagonal_due2essBC;
bool keep_dof_due2essBC; // not impl. yet, must be true
bool special_numbering;
// this bool is set in the initilization
// special_numbering is ON for special numbering and OFF for the
// general numbering. Note that both can not be ON at the same time
bool test_func_is_ok; // set in attachTestFuncDescription
void initSpecialNumbering ();
void initGeneralNumbering ();
public:
DegFreeFE (const GridFE& grid, int ndof_per_node = 1); // isoparam. elements
DegFreeFE (FieldsFE& collection);
DegFreeFE (FieldFE& f); // see redim(FieldFE&) below
DegFreeFE () {}
~DegFreeFE ();
bool redim (const GridFE& grid, int ndof_per_node = 1);
bool redim (FieldsFE& collection);
bool redim (FieldFE& f); // make a FieldsFE from f and call redim(FieldsFE&)
bool generalNumbering () const { return notbool(special_numbering); }
bool specialNumbering () const { return special_numbering; }
void printNumbering (Os os) const;
bool ok () const;
bool update (); // update/redim due to a change in grid size
GridFE& grid ();
const GridFE& grid () const;
int getNoDofInElm (int e) const
{ return specialNumbering() ?
mesh->getNoNodesInElm(e)*ndof_per_node : nelmdof; }
int getNoEqsInElm (int e) const
{ return specialNumbering() ?
mesh->getNoNodesInElm(e)*ndof_per_node : nelmeqs; }
int getTotalNoDof () const { return ntotdof; }
int getTotalNoEqs () const { return ntoteqs; }
int loc2glob (int e, int element_dof) const;
int loc2globDof (int e, int element_dof) const
{ return loc2glob (e, element_dof); }
int loc2globEqs (int e, int element_dof) const
{ return test_func_is_ok ?
test_func->loc2glob (e, element_dof) : loc2glob (e, element_dof);}
void loc2glob (VecSimple(int)& l2g, int e) const;
void loc2globDof (VecSimple(int)& l2g, int e) const
{ loc2glob (l2g, e); }
void loc2globEqs (VecSimple(int)& l2g, int e) const
{ test_func_is_ok ? test_func->loc2glob(l2g,e) : loc2glob(l2g,e); }
int getNoDofPerNode () const { return ndof_per_node; }
int globalDegFree (int field_value_no, int field_no) const
{ return fields2dof (field_value_no, field_no); } // for backward comp.
void attachTestFuncDescription (const BasisFuncGrid& test_functions);
bool noEssBC () const { return getbool(!ebc.ok() || no_ess_bc); }
bool getEssBC (int iglobdof, NUMT& ebc_value) const;
void getEssBC (VecSimple(bool)& essential_dof) const;
void printEssBC (Os os, int output_level = 1);
void initEssBC ();
void fillEssBC (int globdof, NUMT ebc_value);
// this one is convenient for system of PDEs having a common linear system:
void fillEssBC (int field_dof, int field_no, NUMT ebc_value);
void insertEssBC // fill v vector with ess.b.c.
(VecSimple(NUMT)& u, bool account4multiplicity = false);
// fill in correct dimensions of a Matrix_prm object, based on
// current grid and element types:
void updateMatrixPrm (Matrix_prm(NUMT)& mat_prm);
// modified interpretation of EssBC in case of e.g. Newton-Raphson's method:
void fillEssBC2zero ();
void unfillEssBC2zero ();
bool zeroInterpretationOfEssBC () const
{ return interpret_ess_bc_as_zero; }
// linear system numbering <-> field numbering:
void field2vec (const FieldFE& f, Vec(NUMT)& x);
void vec2field (const Vec(NUMT)& x, FieldFE& f);
void field2vec (const FieldsFE& f, Vec(NUMT)& x);
void vec2field (const Vec(NUMT)& x, FieldsFE& f);
int fields2dof (int field_dof, int field_no) const;
void dof2fields (int global_dof, int& field_dof, int& field_no) const;
int computeHalfBandwidth (); // force computation of bandwidth
int getHalfBandwidth (); // avoid comp. if present bandwidth is ok
// these two vectors are used only if(!modifyVecDue2essBC()):
Vec(NUMT) b_mod; // separate storage of modifications to rhs due to ess.b.c.
VecSimple(int) nelm_contr; // no of element contribution for each d.o.f.
void initAssemble (); // init b_mod and nelm_contr
void modifyMatDue2essBC (bool onoff);
bool modifyMatDue2essBC () const { return modify_mat_due2essBC; }
void modifyVecDue2essBC (bool onoff);
bool modifyVecDue2essBC () const { return modify_vec_due2essBC; }
void symmModDue2essBC (bool onoff);
bool symmModDue2essBC () const { return symm_mod_due2essBC; }
void unitMatDiagonalDue2essBC (bool onoff);
bool unitMatDiagonalDue2essBC () const
{ return unit_mat_diagonal_due2essBC; }
void keepDofDue2essBC (bool onoff);
bool keepDofDue2essBC () const { return keep_dof_due2essBC; }
COPY_CONSTRUCTOR(DegFreeFE);
ASSIGNMENT_OPERATOR(DegFreeFE);
private:
int halfbw;
int specHalfBandwidth () const;
int genHalfBandwidth () const;
void init ();
};
class SparseDS;
// FDM/FEM coupling:
//-----------------------------------------------------------------------------
extern void makeSparsityPattern
(
SparseDS& ds, // sparsity pattern data structure
const Ptv(int)& divisions, // divisions in grid
int npt_per_div, // no of points in each division
int ndof_per_pt // no of unknowns per grid point
);
//-----------------------------------------------------------------------------
// general FEM grid coupling (MatSparse):
//-----------------------------------------------------------------------------
extern void makeSparsityPattern (SparseDS& ds, const DegFreeFE& dof);
//-----------------------------------------------------------------------------
// general FEM grid coupling (MatStructSparse):
//-----------------------------------------------------------------------------
extern void makeSparsityPattern (VecSort(int)& offset, int& ndiags,
const DegFreeFE& dof);
//-----------------------------------------------------------------------------
KEYWORDS
finite element, degree of freedom, grid, mesh, linear system
DESCRIPTION
The class relates degrees of freedom in a linear system to the
nodal values of discrete finite element fields. By degrees of
freedom we mean the unknowns in a linear system associated with
some finite element discretization of a system of partial differ
ential equations. For example, when solving a scalar partial
differential equation (e.g. Poisson's equation) by the finite
element method, there is usually one degree of freedom associated
with each node in the finite element mesh. When solving a vector
equation in d dimensions there are d component fields. The rela
tion between, for example, component no j at node k and the cor
responding global degree of freedom in the linear system is com
puted by a DegFreeFE object. The class is simple to use for both
isoparametric and non-isoparametric element meshes. We refer to
the document A Brief Documentation of Finite Element Fields over
Non-isoparametric Elements in Diffpack for the definition of
basic concepts used in and the general ideas behind class
DegFreeFE. This document defines the important concepts general
and special numbering used below.
Sometimes one creates stiffness matrices and corresponding vec
tors where the number of unknowns differ from the number of equa
tions. This is a frequent situation with mixed finite element
methods in combination with sophisticated linear solvers (utiliz
ing a block matrix structure) or in combination with explicit
time integration. Normally, class DegFreeFE assumes that the
number of equations equals the number of degrees of freedoms
(unknowns), but the interface shows that there are two different
functions for these quantities. In the case where the trial func
tions are of another type than the test functions, one can attach
a BasisFuncGrid, reflecting the test functions, to the DegFreeFE
object. This BasisFuncGrid will be used for obtaining information
about the number of equations (and the getNoTotalEqs and getNoE
qsInElm functions will base their information on quantities found
by the BasisFuncGrid description of the test functions). This
way of computing the number of equations is only legal for scalar
element fields and the general numbering, that is, the DegFreeFE
object must be initialized by the constructor or redim function
that takes a FieldsFE object with only one component. Note that
the mapping from local to global degrees of freedom, as provided
by the loc2glob function, is limited to the case where the number
of unknowns equals the number of equations. In the more general
case, the functions loc2globDof and loc2globEqs must be used.
There are overloaded cases of loc2glob, loc2globDof and
loc2globEqs to fill all nodes for an element at once.
The DegFreeFE class may also mark the degrees of freedom that are
subjected to essential boundary conditions. Class ElmMatVec uses
DegFreeFE for this purpose. Moreover, the bandwidth of the coef
ficient matrix in the associated system of linear equations can
be calculated. Note that it is assumed that the size of this
system equals the total number of degrees of freedom in the asso
ciated DegFreeFE object. That is, if a vector equation is solved
sequentially, component by component, a DegFreeFE object with one
degree of freedom per node should be used for each component
equation.
CONSTRUCTORS AND INITIALIZATION
The mostly used constructor is aimed at isoparametric elements
and takes a GridFE object as argument. The number of degrees of
freedom at every node can optionally be given, by default there
is one unknown per node. Note that this constructor implies the
so called special numbering (documented in the report referred to
above).
The other constructor is aimed at the general case where one can
have a vector field with components defined over different types
of non-isoparametric elements. The vector field is represented in
terms of a FieldsFE object, holding all the scalar fields that
are involved in the numbering of degrees of freedom in the linear
system. This FieldsFE object is the argument in the constructor.
Note that this constructor implies the the so called general num
bering when the degrees of freedom are assigned numbers in this
DegFreeFE object. Of course, this FieldsFE object can be very
simple, e.g., holding n FieldFE objects over isoparametric ele
ments. In this case one can also use the constructor that takes
the GridFE object as argument. In more complicated situations,
the FieldsFE object can hold FieldFE objects that have different
BasisFuncGrid. For example, when solving Navier-Stokes-like equa
tions, the FieldsFE may hold a velocity field over isoparametric
elements and a pressure field over non-isoparametric elements.
There is a constructor with no arguments. It performs no initial
ization, but one can at a later time call one of the redim func
tions. These correspond exactly to the other two constructors (in
fact, the constructors simply call the relevant redim function -
this is a standard in Diffpack).
MEMBER FUNCTIONS
redim - redefines the object (performs the same operations as one
of the constructors).
grid - gives access to the associated grid.
getNoDofInElm - returns the total number of degrees of freedom in
a specified element. (Actually, the present version of DegFreeFE
utilize the basic assumption in the present implementation of
BasisFuncGrid, namely that all elements must be of the same type.
Then the number of degrees of freedom in an element is indepen
dent of the e argument).
getTotalNoDof - returns the total number of degrees of freedom in
the (total) coupled global system of algebraic equations. This
will normally be the total number of unknowns in the system. For
a problem involving a scalar function over isoparametric ele
ments, the return value of getTotalNoDof is simply the number of
nodes.
fields2dof - returns the global degree of freedom number associ
ated with a coefficient in a component of the FieldsFE structure
that holds all the coupled, unknown scalar fields. In case of the
special numbering the parameter field_dof corresponds to the node
number and field_no is the local degree of freedom number (field
component number) at this node. In case of the general number
ing, field_dof is the global degree of freedom number for the
scalar component field field_no. The fields2dof function is
related to loc2glob, but the parameters to loc2glob are different
(element number of local degree of freedom number in the ele
ment). Note that the name of fields2dof was globalDegFree in
Diffpack versions earlier than 1.4. The name globalDegFree can
still be used (see below).
globalDegFree - this is a synonym for fields2dof and included for
backward compatibility.
loc2glob - returns the global degree of freedom number corre
sponding to local degree of freedom element_dof in element number
e. Note that element_dof runs from 1 to the total number of
degrees of freedom in the element (as given by getNoDofInElm).
The loc2glob function assumes that the number of unknowns and the
number of equations are the same. In the case where rectangular
matrices are involved and these two numbers differ (that is, a
BasisFuncGrid describing the test functions has been attached to
the DegFreeFE object), one must use loc2globDof and loc2globEqs
for the degree of freedom and equation numbering, respectively.
(The former is just a call to loc2glob; in other words, the
loc2glob function represents the local-global degree of freedom
mapping that be used for the equations in the special case where
the number of equations equals the number of unknowns). The
loc2glob, loc2globDof and loc2globEqs functions are mostly used
by class ElmMatVec and possibly by general matrix sparsity pat
tern computation routines.
"loc2globEqs, loc2globDof - see the description of the loc2glob
function. Overloaded versions of loc2glob, loc2globEqs and
loc2globDof redim and fill the whole local-global vector at once.
getNoDofPerNode - returns the number of degrees of freedom at a
specified node. If positive, the special numbering of the
unknowns is used and the return value equals the constant number
of degrees of freedom per node. A return value equal to zero
implies the general numbering and the function has no application
in this case (except for testing that the special numbering fea
tures do not apply).
dof2fields - this is the inverse function of fields2dof. The
input arguments to fields2dof are the return values of dof2fields
and vice versa.
getEssBC - returns true if a global (system) degree of freedom
number is subjected to an essential boundary condition. The
boundary value is available as an argument to the function.
There is an overloaded version that fills a vector of booleans
according to whether there is an essential dof or not.
printEssBC - prints the essential conditions. Aimed at debugging.
fillEssBC - sets an essential boundary condition at a degree of
freedom.
initEssBC - a function that erases all essential boundary condi
tions (usually called prior to fillEssBC).
printEssBC - prints the degrees of freedom where essential bound
ary conditions are assigned. This function is particularly useful
when debugging simulators as one can check that the programmer's
assignment of essential boundary conditions in a DegFreeFE object
is correct. This information can then be compared to the solution
to check that the boundary conditions propagate through the whole
simulation. There is an optional argument to printEssBC that
decides the level of output (or verbosity). A value of 0 makes
the output as compact as possible, the value 1 is standard, and 2
will print the coordinates of the point at which a degree of
freedom is associated.
printNumbering - prints the numbering of the global degrees of
freedom as they are computed by the DegFreeFE class. The output
is written element by element, local node by local node etc. This
function may be of great help to gain a complete understanding of
the numbering techniques. Make a simple grid with few elements,
study the output of printNumbering and make a sketch of the grid
and the degrees of freedom locations on a piece of paper.
fillEssBC2zero - forces all essential boundary condition values
to be interpreted as zero. If the system of linear equations have
an unknown which is a correction vector, special care must be
taken. This is the case when e.g. Newton's method is used for
solving systems of nonlinear equations. In each linear problem it
is assumed that the total solution has its boundary conditions
inserted and that the correction vector must have zeroes for all
entries where the value is known. That is, in the modification of
element matrices one must set the essential condition to be zero.
The ElmMatVec::enforceEssBC function calls the zeroInterpreta
tionOfEssBC of DegFreeFE for determining whether the user has
indicated to the DegFreeFE object that all non-zero essential
boundary condition values are to be interpreted as if they were
zero. When Newton's method is used, it is important that the pro
grammer calls the fillEssBC2zero function of the DegFreeFE class
object prior to any enforceEssBC.
unfillEssBC2zero - essential boundary condition values are inter
preted as their original values.
zeroInterpretationOfEssBC - returns true if essential boundary
conditions are to be treated as zero values (in e.g. Newton's
method). In other words, the function returns true if fil
lEssBC2zero has been called and not yet been cleared by unfil
lEssBC2zero. Note that only the last call to fillEssBC2zero or
unfillEssBC2zero is remembered.
field2vec - given a field, the function loads the field into a
Vec that can be used in a linear system. Recall that the main
purpose of the DegFreeFE class is to administer the numbering of
field values versus unknowns in associated linear systems, and
field2vec (and vec2field) contains basic functionality for
switching between Vec and FieldFE or FieldsFE representations.
The numbering of the components in a vector and a field represen
tation has a simple relation if the special numbering technique
is used. For the general numbering technique, the details of the
numbering strategy is documented elsewhere. Note that if the gen
eral numbering is used, the field2vec and vec2field functions
that take a FieldFE argument cannot be used; one must always com
municate with the DegFreeFE object through a FieldsFE object
(e.g. the object that was given as the collection argument when
initializing the DegFreeFE object). (The vec2field and field2vec
functions that take a FieldFE argument will result in a fatal
error if the general numbering is chosen). The field2vec and
vec2field functions are very efficient since they use either very
simple formulas or pre-computed internal tables for switching
between vector and field formats.
vec2field - given a Vec from a linear system (e.g., the solu
tion), the function loads the vector into a field. See field2vec.
getHalfBandwidth - computes the bandwidth of the coefficient
matrix in a linear system associated with the grid and the
degrees of freedom. The number to be returned is max_e |i-j|+1,
where e denotes an element and i and j are degree of freedom num
bers in that element (i.e., i and j are unknowns in the element).
The function returns a local variable containing the bandwidth if
the number of degrees of freedom in the object has not changed.
Hence, the function is quite efficient.
computeHalfBandwidth - as getHalfBandwidth but the computation of
the bandwidth is enforced (that is, the local variable containing
the bandwidth is not used).
updateMatrixPrm - fill in correct dimensions of a Matrix_prm
object, based on the current grid information.
modifyMatDue2essBC, modifyVecDue2essBC, symmModDue2essBC, unit
MatDiagonalDue2essBC, keepDofDue2essBC - these functions control
the implementation of essential boundary conditions in the linear
system. By default, all these modes are ON, that is, the func
tions return true boolean values. The programmer can set the
modes ON or OFF to control the modifications of the linear system
in detail. Especially in problems with time independent coeffi
cients in the PDEs, where the finite element assembly process can
be avoided at every time step, it is necessary to control the
modifications of various components of the linear system. This is
documented elsewhere (see the Users Guide or the book "Numerical
Solution of Partial Differential Equations; Models, Algorithms,
and Software").
modifyMatDue2essBC - indicator for modification of the coeffi
cient matrix in the linear system due to essential boundary con
ditions. The function is used by ElmMatVec::enforceEssBC to
check if the element matrix should be manipulated or not. The
function that takes the bool onoff argument is used to select
modification (ON) or no modification (OFF). It is usually called
from the programmers simulation code. Since the use of these
indicator functions are in the finite element assembly process,
the call must be performed prior to the call to FEM::makeSystem.
The other version of the indicator function returns a bool, and
this function is used by ElmMatVec::enforceEssBC.
modifyVecDue2essBC - indicator for modification of the right hand
side in the linear system due to essential boundary conditions.
See the documentation of modifyMatDue2essBC for more information.
The modifyVecDue2essBC automatically assembles all the contribu
tions to the right hand side corresponding to subtractions of
columns of the coefficient matrix. The contributions are avail
able as the Vec(NUMT) vector b_mod. Moreover, the values of the
unknowns are available from the function getEssBC. With b_mod and
getEssBC it is easy to enforce the essential boundary conditions
in the right hand side data structure at any stage in the compu
tations.
symmModDue2essBC - indicator for symmetric modification of the
linear system due to essential boundary conditions. If degree of
freedom number i is known, equation number i in the linear system
is replaced by the boundary condition. However, such a modifica
tion destroys the symmetry of the coefficient matrix. If symmMod
Due2essBC is ON, zeroes will also be introduced in column number
i in the coefficient matrix. This additional symmetric modifica
tion also requires the right hand side to be modified (column no
i in the coefficient matrix times the known solution value must
be subtracted from the right hand side).
Note that the modify functions tell the finite element tools
(class ElmMatVec) whether it is allowed at all to modify (touch)
the arrays in the linear system. The symmModDue2essBC function
tells the tools if the implementation of the boundary conditions
are to be done in a symmetric or non-symmetric way.
unitMatDiagonalDue2essBC - indicates if 1.0 is set on the diago
nal of the coefficient matrix in the position where there is an
essential boundary condition (this is the normal choice). Other
wise, 0.0 is set on the diagonal (this is convenient if the coef
ficient matrix is on the form A=M+a*K, where K and M are computed
separately - M may have 1.0 on the diagonal and K may have 0.0
such that the sum M+a*K implements the essential boundary condi
tion correctly).
keepDofDue2essBC - essential boundary conditions can be removed
from the linear system to reduce the system size. This function
indicates whether degrees of freedom that are known (due to
essential boundary conditions) should be kept (default) or elimi
nated. NOTE: At present, elimination is not yet implemented!
b_mod - data structure for storing right hand side modifications
due to essential boundary conditions. See modifyVecDue2essBC,
nelm_contr and insertEssBC. The updating of b_mod (and
nelm_contr) is usually done by ElmMatVec::enforceEssBC.
nelm_contr - vector that counts the number of elements contribut
ing to a degree of freedom. Useful when the right hand side of a
linear system is modified, due to essential boundary conditions,
at the global level while the coefficient matrix has been modi
fied at the element level. When inserting the essential boundary
conditions values in the right hand side, one must match the mul
tiplicity of that degree of freedom in the matrix, that is, the
value must be multiplied by the number of element contributions,
which is available nelm_contr. The updating of nelm_contr is
usually done by ElmMatVec::enforceEssBC. See also the
insertEssBC function.
initAssemble - initializes the DegFreeFE object prior to an
assembly process. This means that b_mod and nelm_contr are reset
to zero, if these vectors are used (the vectors are used only if
modifyVecDue2essBC() is false - in that case the vectors hold
information about the modifications of the right hand side due to
essential boundary conditions). See also the insertEssBC func
tion.
insertEssBC - inserts essential boundary conditions in a vector.
This is useful for start vectors for linear or nonlinear solvers
and when incorporating essential boundary conditions after a lin
ear system is computed. Here are some comments regarding the lat
ter application. When modifyVecDue2essBC() is false, the modifi
cations of the right hand side due to essential boundary condi
tions are not performed at the element level. Instead the modifi
cations is stored in internal arrays (b_mod and nelm_contr) for
later use. The insertEssBC function utilizes the stored informa
tion and inserts the correct boundary values in the rhs argument.
The argument account4multiplicity is by default true, which means
that the inserted value and its position are found by getEssBC
and then multiplied by the number of elements contributing to
this degree of freedom. Usually, the modifications of the coeffi
cient matrix due to essential boundary are accomplished at the
element level (since that is much more straightforward for, e.g.,
sparse matrices) and the "unit" values on the main diagonal are
then not 1, but n, where n is the number of elements contributing
to the degree of freedom in question. The modifications of the
right hand side must reflect this and the necessary information
is provided by nelm_contr. The updating of b_mod and nelm_contr
is usually done by ElmMatVec::enforceEssBC.
EXAMPLES
Here is an example where one sets the boundary conditions of a
velocity field. At all nodes where indicator 1 is on, the veloc
ity is to be set to zero. The DegFreeFE object (handle called
dof) is initialized in this way:
dof.rebind (new DegFreeFE(grid(),nsd));
where grid is a Handle(GridFE) and nsd is the number of space
dimensions, that is, the number of vector components.
dof->initEssBC ();
const int nno = grid->getNoNodes();
const int nsd = grid->getNoSpaceDim();
int i,k;
for (i = 1; i <= nno; i++)
if (grid->BoNode (i, 1))
for (k = 1; k <= nsd; k++)
dof->fillEssBC (dof->fields2dof(i, k, 0.0));
See class NavierStokes or class Elasticity for simulators that
works with vector fields and DegFreeFE.
SEE ALSO
class FieldFE, class FieldsFE, class GridFE, class LinEqAdm
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO. Optimizations and modifica
tions by Klas Samuelsson, UiO.