NAME
DegFreeFD - degree of freedom manager for finite difference lin
ear systems
INCLUDE
include "DegFreeFD.h"
SYNTAX
class DegFreeFD : public HandleId
{
friend class DegFreeFDLattice;
friend class DegFreeFDUnstruct;
public:
// public for use in makeSparsityPattern functions:
Handle(GridWithPts) mesh; // basic info for the numbering
Handle(FieldsWithPtValues) fields; // basic info for the numbering
Handle(IndexSet) indexset;
Handle(IndexSetCollection) indexsetcoll;
Handle(AlgebraicContribs) alg_contribs; // extra algebraic eqns/dofs
private:
int npdedof; // total no of dof in the system
int npdeeqs; // total no of equations in the system
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; // essential boundary conditions
bool no_ess_bc; // true if no conditions are set
bool interpret_ess_bc_as_zero; // needed in e.g. Newton iter
// indicators for various implementations of essential boundary conditions:
bool symm_mod_due2essBC;
bool modify_vec_due2essBC;
bool modify_mat_due2essBC;
void init ();
protected:
virtual int calcPDEHalfBandWidth (StencilCollection& stencils) =0;
public:
DegFreeFD () {}
DegFreeFD (const GridWithPts& grid, int ndof_per_node = 1);
DegFreeFD (const GridWithPts& grid, const IndexSet& indexset,
int ndof_per_node = 1);
DegFreeFD (const GridWithPts& grid,
const IndexSetCollection& indexsetcoll,
int ndof_per_node = 1);
DegFreeFD (FieldsWithPtValues& collection);
DegFreeFD (FieldsWithPtValues& collection,
const IndexSet& indexset);
DegFreeFD (FieldsWithPtValues& collection,
const IndexSetCollection& indexsetcoll);
~DegFreeFD ();
bool redim (const GridWithPts& grid, int ndof_per_node = 1);
bool redim (const GridWithPts& grid, const IndexSet& indexset,
int ndof_per_node = 1);
bool redim (const GridWithPts& grid,
const IndexSetCollection& indexsetcoll,
int ndof_per_node = 1);
bool redim (FieldsWithPtValues& collection);
bool redim (FieldsWithPtValues& collection,
const IndexSet& indexset);
bool redim (FieldsWithPtValues& collection,
const IndexSetCollection& indexsetcoll);
void printNumbering (Os os) const;
virtual bool ok () const;
bool update (); // update/redim due to a change in grid size
GridWithPts& grid ();
const GridWithPts& grid () const;
// For treatment of algebraic eqns/dofs
void connect (const AlgebraicContribs& alg_contribs_);
void disconnect ();
void getBaseIndex (int& eqn_base, int& dof_base);
int getAlgNoDof () const;
int getAlgNoEqs () const;
int getPDENoDof () const;
int getPDENoEqs () const;
int getTotalNoDof () const;
int getTotalNoEqs () const;
int getNoDofPerNode () const { return ndof_per_node; }
virtual int loc2glob (int node_no, int node_dof) const;
virtual void loc2glob (VecSimple(int)& l2g, int node_no) const;
void initEssBC ();
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 getEssBC (Vec(NUMT)& v);
void printEssBC (Os os, int output_level = 1);
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 fillEssBC2zero ();
void unfillEssBC2zero ();
bool zeroInterpretationOfEssBC () const
{ return interpret_ess_bc_as_zero; }
// linear system numbering <-> field numbering:
void vec2field (const Vec(NUMT)& x, FieldWithPtValues& f);
void field2vec (const FieldWithPtValues& f, Vec(NUMT)& x,
bool allocate_alg_dof = true);
void vec2field (const Vec(NUMT)& x, FieldsWithPtValues& f);
void field2vec (const FieldsWithPtValues& f, Vec(NUMT)& x,
bool allocate_alg_dof = true);
int fields2dof (int field_dof, int field_no) const;
void dof2fields (int global_dof, int& field_dof, int& field_no) const;
int getHalfBandWidth (StencilCollection& stencilcoll);
// this vector is used only if(!modifyVecDue2essBC()):
Vec(NUMT) b_mod; // separate storage of modifications to rhs due to ess.b.c.
void initAssemble (); // init b_mod
void insertEssBC (VecSimple(NUMT)& ebc);
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; }
// Handling of algebraic eqns and dofs
Handle(VecSimpleH(VecSimple(int))) alg_pattern;
void initAlgPattern ();
// Utility functions for computation of sparsity patterns
virtual void findSparsePatternFromIndexSet
(
IndexSet& indexset,
StencilCollection& stencils,
VecSimple(VecSimple(int))& col_idx,
int indexset_no
) =0;
virtual void findSparsePatternFromGrid
(
StencilCollection& stencils,
VecSimple(VecSimple(int))& col_idx
) =0;
virtual void findOffsetPatternFromIndexSet
(
IndexSet& indexset,
StencilCollection& stencils,
ArrayGenSimple(int)& offsetmarks,
int indexset_no
) =0;
virtual void findOffsetPatternFromGrid
(
StencilCollection& stencils,
ArrayGenSimple(int)& offsetmarks
) =0;
COPY_CONSTRUCTOR(DegFreeFD);
ASSIGNMENT_OPERATOR(DegFreeFD);
CLASS_INFO
// virtual cast functions, one for each subclass:
DEF_VIRTUAL_CAST(DegFreeFDLattice)
DEF_VIRTUAL_CAST(DegFreeFDUnstruct)
};
// general FDM grid coupling (MatSparse):
//-----------------------------------------------------------------------------
extern void makeSparsityPattern (SparseDS& ds,
StencilCollection& stencils,
DegFreeFD& dof);
//-----------------------------------------------------------------------------
// general FDM grid coupling (MatStructSparse):
//-----------------------------------------------------------------------------
extern void makeSparsityPattern (VecSort(int)& offset, int& ndiags,
StencilCollection& stencils,
DegFreeFD& dof);
//-----------------------------------------------------------------------------
KEYWORDS
finite difference, 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 differencet fields. By degrees of
freedom we mean the unknowns in a linear system associated with
some finite difference discretization of a system of partial dif
ferential equations. For example, when solving a scalar partial
differential equation (e.g. Poisson's equation) by the finite
difference method, there is usually one degree of freedom associ
ated with each node in the finite difference mesh. When solving
a vector equation in d dimensions there are d component fields.
The relation between, for example, component no j at node k and
the corresponding global degree of freedom in the linear system
is computed by a DegFreeFD object.
The DegFreeFD class also mark the degrees of freedom that are
subjected to essential boundary conditions. Class MatVecContribFD
uses DegFreeFD for this purpose. Moreover, the bandwidth of the
coefficient 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 DegFreeFD object. That is, if a vector equation is solved
sequentially, component by component, a DegFreeFD object with one
degree of freedom per node should be used for each component
equation.
CONSTRUCTORS AND INITIALIZATION
The most common constructors to use are the ones that take a
GridWithPts objects and an IndexSetCollection or a FieldsWithPt
Values and an IndexSetCollection. In the constructors that take a
grid, one must specify the number og degrees of freedom at every
node if it differs from 1. The constructors that take a
FieldsWithPtValues sets the number of degrees of freedom per node
to the number of components in the FieldsWithPtValues object.
Optionally one can use the constructors that take a GridLattice
object and the number of degrees of freedom per node, or the sim
ilar constructor that take an IndexSet. Correspondingly, we have
a constructor taking a FieldsWithPtValues objects and one taking
a FieldsWithPtValues and an IndexSet.
There is also a constructor without arguments. It performs no
initialization, but one can at a later time call one of the redim
functions. These correspond exactly to the other 6 constructors
(in fact, the constructors simply call the relevant redim func
tion - this is a standard in Diffpack).
MEMBER FUNCTIONS
printNumbering - prints the numbering of the global degrees of
freedom as they are computed by the DegFreeFD class. This func
tion may be of great help to gain a complete understanding of the
numbering techniques. Make a simple grid with few nodes, study
the output of printNumbering and make a sketch of the grid and
the degrees of freedom on a piece of paper.
grid - gives access to the associated grid.
connect - attach the supplied AlgebraicContribs object to the
"DegFreeFD
disconnect - disconnect the AlgebraicContribs object from the
DegFreeFD
getBaseIndex - returns the base index for the algebraic degrees
of freedom and equations. The base index for the algebraic
degrees of freedom equals getPDENoDof and the base index for the
algebraic equations equals getPDENoEqs.
getAlgNoDof - returns the number of degrees of freedom from the
algebraic contributions.
getAlgNoEqn - returns the number of equations from the algebraic
contributions.
getPDENoDof - returns the total number of degrees of freedom from
the partial differential equation part of the system. This is
number of nodes in the grid times the number of degrees of free
dom per node.
getPDENoEqn - returns the total number of equations from the par
tial differential equation part of the system. This is number of
nodes in the grid times the number of degrees of freedom per
node.
getTotalNoDof - returns the total number of degrees of freedom in
the global system. This is number of nodes in the grid times the
number of degrees of freedom per node, plus the number of alge
braic degrees of freedom.
getTotalNoEqs - returns the total number of equations in the
global system. This is number of nodes in the grid times the
number of degrees of freedom per node, plus the number of alge
braic equations.
loc2glob - returns the global degree of freedom number corre
sponding to local degree of freedom node_dof in node number
node_no. Note that node_dof runs from 1 to the total number of
degrees of freedom in a node (as given by getNoDofInNode). The
loc2glob function does not cover local2global mapping of the
algebraic contributions, this is handled by the class Algebraic
Contribs.
initEssBC - a function that erases all essential boundary condi
tions (must be called prior to fillEssBC).
noEssBC - returns true if one ore more boundary condition has
been set, otherwise false.
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. Another
overloaded version fills a vector with the essential boundary
condition if this has been set. Otherwise the entry is left as it
is. This is for instance used when one should fill the vector
holding the non-linear solution with the essential boundary
conditions.
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 DegFreeFD 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.
fillEssBC - set an essential boundary condition. The entry is
either specified by the global degree of freedom or the grid
point and the local degree of freedom in the node.
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 MatVecContribFD:: enforceEssBC function calls the zeroInter
pretationOfEssBC of DegFreeFD for determining whether the user
has indicated to the DegFreeFD 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 DegFreeFD 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 DegFreeFD 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 FieldWithPtValues or FieldsWithPtValues
representations. The field2vec and vec2field functions are very
efficient since they use either very simple formulas or pre-com
puted internal tables for switching between vector and field for
mats. If the supplied vector is not initialized, it will be redi
mensioned to the correct size. If the vector is already initial
ized but with the wrong size, an error message will be given. If
the optional third argument to this function is true, the sup
plied vector will be set to have dimension getTotalNoDofs, while
if this argument is false, no space is allocated for the alge
braic degrees of freedom, and the vector will then have size
equal to getPDENoDof.
vec2field - given a Vec from a linear system (e.g., the solu
tion), the function loads the vector into a field. The vector
must have size equal to the number of nodes in the field or total
number of degrees of freedoms in the system. In the latter case,
the only the PDE part of the vector is transferred to the field,
the algebraic degrees of freedom are ignored. 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_n |i-j|+1,
where n denotes a node and i and j are degree of freedom numbers
in that node (i.e., i and j are unknowns in the node).
fields2dof - returns the global degree of freedom number associ
ated with a coefficient in a component of the FieldsWithPtValues
structure that holds all the coupled, unknown scalar fields. 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.
dof2fields - this is the inverse function of fields2dof. The
input arguments to fields2dof are the return values of dof2fields
and vice versa.
b_mod - data structure for storing right hand side modifications
due to essential boundary conditions. See modifyVecDue2essBC, and
insertEssBC. The updating of b_mod is usually done by MatVecCon
tribFD:: enforceEssBC.
initAssemble - initializes the DegFreeFD object prior to an
assembly process. This means that b_mod is reset to zero, if
these vectors are used (the vectors are used only if modifyVec
Due2essBC() 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 function.
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 an internal array (b_mod) for later use. The
insertEssBC function utilizes the stored information 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 coefficient matrix
due to essential boundary are accomplished at the node level
(since that is much more straightforward for, e.g., sparse matri
ces). The updating of b_mod is usually done by MtVecContribFD::
enforceEssBC.
modifyMatDue2essBC - indicator for modification of the coeffi
cient matrix in the linear system due to essential boundary con
ditions. The function is used by MatVecContribFD:: enforceEssBC
to check if the contribution 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 difference assembly pro
cess, the call must be performed prior to the call to FDM::
makeSystem. The other version of the indicator function returns
a bool, and this function is used by MatVecContribFD::
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 possible symmetry of the coefficient matrix. If
symmModDue2essBC is ON, zeroes will also be introduced in column
number i in the coefficient matrix. This additional symmetric
modification 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 difference tools
(class MatVecContribFD) 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.
initAlgPattern - initialize tha pattern for the algebraic contri
butions.
SEE ALSO
class FieldWithPtValues, class FieldsWithPtValues, class Grid
WithPts, class LinEqAdmFD
DEVELOPED BY
Numerical Objects AS, Oslo, Norway
AUTHOR
Adapted from DegFreeFE by Elizabeth Acklam, NO