inputs/dev/EssentialConditions.tex

%%%%%%%%%%%%%%%%%%%
% XLiFE++ is an extended library of finite elements written in C++
%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%

To take into account linear essential conditions (conditions acting on unknown or test function spaces), \xlifepp proposes a general approach based on an abstract algebraic representation of essential conditions, say
$$\displaystyle \mathbb{C}U=G$$
where $U$ represents the vector of components of a function $u$ in space ${\cal U}=\left\{\sum_{i=1,n}U_i\,w_i\right\}$, $(w_i)_{i=1,n}$ the basis of discrete space,  $\mathbb{C}$ a $p\times n$ matrix and $G$ a $p$ vector. In other words, we consider the discrete constrained space :
\begin{equation} \label{cons_Uspace}
\displaystyle \widetilde{\cal U}=\left\{\sum_{i=1,n}U_i\,w_i,\ \mathbb{C}U=G\right\}.
\end{equation}
Such representation may be also used for constraints on test function :
\begin{equation} \label{cons_Vspace}
\displaystyle \widetilde{\cal V}=\left\{\sum_{i=1,m}V_i\tau_i,\ \mathbb{D}V=0\right\}.
\end{equation}
Both may exist in the same time and may be different. Nevertheless, constraints on test functions are always homogeneous ones! In usual cases, the test function constraints are the linear part of unknown constraints ($\mathbb{C}=\mathbb{D}$).

\section{A general process}

The first step consists in writing essential conditions given in a symbolic form (see documentation of {\class EssentialConditions} class)  as a constraints system, details being exposed in documentation of {\class Constraints} class). As constraints may be redundant, the constraints system has to be reduced to a minimal system, eliminating redundant conditions or identifying conflicts in constraints. This is achieved using QR factorization. Once the constraints systems are reduced, the main work consists to use them in matrix problems we deal with. As we will see, it is possible to consider different techniques : \emph{real reduction} (some rows and columns of system are deleted), \emph{pseudo reduction} (some rows and columns of system are replaced), \emph{penalization} (some rows and columns of system are penalized) and \emph{duality} (constraints system and new dual unknowns are introduced in the problem).

\subsection*{Reinterpretation of a linear system with constraints}

We consider a general discrete variational problem :
\begin{equation} \label{fv1}
\displaystyle \text{find } u\in \widetilde{\cal U} \text{ such that } a(u,v)=l(v)\ \ \forall v\in \widetilde{\cal V},
\end{equation}
where the spaces $\widetilde{\cal U}$ and $\widetilde{\cal V}$ are given by (\ref{cons_Uspace}-\ref{cons_Vspace}), $a(.,.)$ a bilinear form on ${\cal U}\times {\cal V}$ and $l(.)$ a linear form on ${\cal V}$.

We assume in the following, that the constraints systems $\mathbb{C}U=G$ and $\mathbb{D}V=0$ are reduced, so the $m\times n$ matrix $\mathbb{C}$ and the $p\times q$ matrix $\mathbb{D}$ have the following representation :
\begin{equation} \label{fv2}
\displaystyle  \mathbb{C} = [\mathbb{E}\  \mathbb{R}] \text{ and }  \mathbb{D} = [\mathbb{F}\ \mathbb{S}]
\end{equation}
where $\mathbb{E}$ is a  $m\times m$ upper triangular matrix (invertible), $\mathbb{R}$ a  $m\times n-m$ matrix (may be null, $n\geq m$), $\mathbb{F}$ a $p\times p$ upper triangular matrix (invertible)  and $\mathbb{S}$ a  $p\times q-p$ matrix (may be null, $q\geq p$). The column indices ${\cal E}$ (resp.${\cal F}$) of $\mathbb{E}$ (resp. $\mathbb{F}$) correspond to the eliminated unknowns (resp. test functions) and the column indices ${\cal R}$ (resp. ${\cal S}$) of $\mathbb{R}$ (resp. $\mathbb{S}$) correspond to reduced unknowns (resp. test functions).

Introducing   $\mathbb{X}=\mathbb{E}^{-1}\,\mathbb{R}$ and $\mathbb{Y}=\mathbb{F}^{-1}\,\mathbb{S}$, the constraints read :
\begin{equation} \label{cons_Ured}
\displaystyle U_{\cal E} + \mathbb{X}\,U_{\cal R} = \mathbb{E}^{-1}G \text{  and  } V_{\cal F} + \mathbb{Y}\,V_{\cal S} = 0
\end{equation}
where $U_{\cal E}$  (resp. $U_{\cal R}$) denotes the unknown components restricted to indices ${\cal E}$ (resp. ${\cal R}$) and $V_{\cal F}$  (resp. $V_{\cal S}$) denotes the test function components restricted to indices ${\cal F}$ (resp. ${\cal S}$).

The space $\widetilde{\cal U}$ may be reinterpreted by eliminating the components $U_{\cal E}$ :
$$
\begin{array}{ll}
\displaystyle \sum_{i=1,n}U_i\,w_i & \displaystyle=\sum_{i\in \cal E}U_i\,w_i+\sum_{i\notin {\cal E}}U_i\,w_i\\
&\displaystyle =-\sum_{i\in \cal E}\sum_{j\in {\cal R}}\mathbb{X}_{ij}U_j\,w_i+\sum_{j\notin {\cal E}}U_j\,w_j +\sum_{i\in \cal E}(\mathbb{E}^{-1}G)_iw_i\\
&\displaystyle =\sum_{j\notin {\cal E}}u_j\left(w_j-\sum_{i\in \cal E}\mathbb{X}_{ij}w_i\right)+\sum_{i\in \cal E}(\mathbb{E}^{-1}G)_iw_i
\ \ (\text{setting } \mathbb{X}_{ij}=0\ \forall i,\ \forall j\notin \cal E). \\
\end{array}
$$
Considering the functions
\begin{equation} \label{fun_cu}
\widetilde{w}_j=w_j-\sum_{i\in \cal E}\mathbb{X}_{ij}w_i\ \ j\notin \cal E,
\end{equation}
the vector space $\widetilde{\cal U}_0$ associated to $\widetilde{\cal U}$ is generated by the functions $\widetilde{w}_j$ :
$$\displaystyle \widetilde{\cal U}_0=\left\{\sum_{i\notin \cal E}U_i\widetilde{w}_i\right\}.$$
In a same way, the vector space $\widetilde{\cal V}$ is generated by the functions :
\begin{equation} \label{fun_cv}
\widetilde{\tau}_j=\tau_j-\sum_{i\in \cal F}\mathbb{Y}_{ij}w_i\ \ j\notin \cal F,
\end{equation}
and
$$\displaystyle \widetilde{\cal V}=\left\{\sum_{i\notin \cal F}U_i\widetilde{\tau}_i\right\}.$$
Setting $g=\sum_{i\in \cal E}(\mathbb{E}^{-1}G)_iw_i$, the problem (\ref{fv1}) is equivalent to
\begin{equation} \label{fv0}
\displaystyle \text{find } \widetilde{u}\in \widetilde{\cal U}_0 \text{ such that } a(\widetilde{u},\overline{\widetilde{v}})=
l(\overline{\widetilde{v}})-a(g,\overline{\widetilde{v}})\ \ \forall v\in \widetilde{\cal V}
\end{equation}
equivalent to
\begin{equation} \label{sys0}
\displaystyle \text{find } \widetilde{u}=\sum_{j\notin \cal E}\widetilde{U}_j\widetilde{w}_j \text{ such that }
\sum_{j\notin \cal E}\widetilde{U}_j a(\widetilde{w}_j,\overline{\widetilde{\tau}}_i)=l(\overline{\widetilde{\tau}}_i)-a(g,\overline{\widetilde{\tau}}_i)\ \ \forall i \notin \cal F.
\end{equation}

The system (\ref{sys0}) may be explicited by using (assuming $\tau_i$ is a real function) :
\begin{equation} \label{mat_el}
a(\widetilde{w}_j,\overline{\widetilde{\tau}}_i)=a(w_j,\tau_i)
-\sum_{k\in {\cal E}}\mathbb{X}_{kj}a(w_k,\tau_i)
-\sum_{k\in {\cal F}}\overline{\mathbb{Y}}_{ki}a(w_j,\tau_k)
+\sum_{k\in {\cal E}}\sum_{l\in {\cal F}}\mathbb{X}_{kj}\overline{\mathbb{Y}}_{li}a(w_k,\tau_l)
\end{equation}
and
\begin{equation} \label{rhs_el}
l(\overline{\widetilde{\tau}}_i)=l(\tau_i)-\sum_{k\in {\cal F}}\overline{\mathbb{Y}}_{ki}l(\tau_k).
\end{equation}

Finally, from the knowledge of all $a(w_j,\tau_i)$ and $l(\tau_i)$ corresponding to the computation without essential conditions, it is possible to construct the system  (\ref{sys0}). In terms of matrix, the different operations are mainly row or column combinations. \xlifepp implements this general scheme. In a simple case, the Dirichlet problem for instance, the previous expression are trivial because $\mathbb{X}=\mathbb{Y}=0$, $\widetilde{w}_j=w_j,\ \forall j\notin {\cal E}$ and $\tau_j=w_j$.

When condition is not homogeneous ($g\neq 0$), there is a contribution from $g$ and a part of matrix representing $a(.,.)$ (columns with index in $\cal E$). In order to avoid additional reduction process when an other data $g$ is considered, it is possible to keep in memory this part of matrix.

\subsection*{Overview of classes involved}
To describe essential conditions, \xlifepp uses :
\begin{itemize}
\item some classes used to describe (bi)linearform : {\class Unknown DifferentialOperator OperatorOnUnknown}
\item the class {\class LcOperatorOnUnknown} which handles linear combination of OperatorOnUnknown and allow to associate geometric domains ({\class GeomDomain}). The general syntaxes look like ($\diamond$ represent any algebraic operation consistant with objects involved) :
$$a_1*op_1(f_1\diamond u_1\diamond g_1)|dom_1 + a_2*op_2(f_2\diamond u_2\diamond g_2)|dom_2 + ...$$
$$(a_1*op_1(f_1\diamond u_1\diamond g_1) + a_2*op_2(f_2\diamond u_2\diamond g_2)+ ...)|dom$$
This class supports a lot of expressions but just a few ones could be used as essential condition ! This class belongs to the {\lib Operator} library.
\item To transform a previous expression as an essential condition, it is required to give a right hand side :
\begin{center}
LcOperatorOnUnknown = constant or LcOperatorOnUnknown = function
\end{center}
it produces an object of {\class EssentialCondition} class.
\item Concatanating {\class EssentialCondition} objects using the operator \& produces a set of essential conditions managed by {\class EssentialConditions} class.
\end{itemize}

All these classes are user classes!\\

{\class Constraints} is the core class implementing most of the algorithms : transformation of symbolic conditions into constraints system, reduction of constraints system, reduction or pseudo-reduction of matrix, right hand side correction.

\section{The {\classtitle EssentialCondition} class}

The {\class EssentialCondition} class handles the symbolic representation of an essential condition in the form $expression(u_1|dom_1,u_2|dom_2,...) = f$ where $expression(...)$ is a {\class LcOperatorOnUnknown} object and $f$ is a {\class Function} object.
\vspace{.3cm}
\begin{lstlisting}[]{}
class EssentialCondition
{
  protected :
  LcOperatorOnUnknown ecTerms_;  // combination of operator on unknown
  Function *fun_p;               // function defining data of EC
                                 // (=0 when EC is homogeneous)
  EcType type_;                  // type of essential condition
}
\end{lstlisting}
\vspace{.3cm}
The type of essential condition are listed in the \emph{EcType} enumeration:
\begin{lstlisting}[]{}
enum EcType {_undefEcType,      //undefined type
             _DirichletEc,      //one unknown, one domain
             _transmissionEc,   //two unknowns, one domain
             _periodicEc,       //one unknown, two domains
             _crackEc,          //two unknowns, two domains
             _meanEc            //intg_D(opu)=f
            };
\end{lstlisting}
\vspace{.3cm}
The class proposes basic constructors which do full copy of  {\class LcOperatorOnUnknown} and {\class Function} objects:
\begin{lstlisting}[]{}
EssentialCondition()
EssentialCondition(const LcOperatorOnUnknown&);
EssentialCondition(const LcOperatorOnUnknown&, const Function&);
EssentialCondition(Domain&, const LcOperatorOnUnknown&);
EssentialCondition(Domain&, const LcOperatorOnUnknown&, const Function&);
EssentialCondition(const EssentialCondition&);
EssentialCondition& operator=(const EssentialCondition&);
\end{lstlisting}
\vspace{.3cm}
some accessors or shotrcuts :
\begin{lstlisting}[]{}
const std::vector<OpuValPair>& bcTerms() const;
const Function& fun() const;
Function*& funp();
it_opuval begin();
cit_opuval begin() const;
it_opuval end();
cit_opuval end() const;
Number size() const;
EcType type() const;
\end{lstlisting}
\vspace{.3cm}
some function giving properties
\begin{lstlisting}[]{}
bool isSingleUnknown() const;
Dimen nbOfUnknowns() const;
const Unknown* unknown() const;
std::set<const Unknown*> unknowns() const;
bool isSingleDomain() const;
Domain* domain() const;
std::set<Domain*> domains() const;
Dimen nbOfDomains() const;
Number nbTerms() const;
Complex coefficient() const;
std::vector<Complex> coefficients() const;
std::set<DiffOpType> diffOperators() const;
bool nbOfDifOp() const;
DiffOpType diffOperator() const;
bool isHomogeneous() const;
String name() const;
\end{lstlisting}
\vspace{.3cm}
some functions to set type and domain:
\begin{lstlisting}[]{}
EcType setType();              //set the type of the Essential condition
void setDomain(Domain& dom)    //affect domain to EC
EssentialCondition& operator |(Domain&)    //affect domain of EC
\end{lstlisting}
\vspace{.3cm}
When affect the domain using \emph{setDomain} function or | operator, all the domains defined in EcTerms are reset to the new domain. This domain affectation is also done by the non member function:
\begin{lstlisting}[]{}
EssentialCondition on(const Domain&, const EssentialCondition&);
\end{lstlisting}
\vspace{.3cm}
and printing facilities:
\begin{lstlisting}[]{}
void print(std::ostream&) const;
friend std::ostream& operator<<(std::ostream&, const EssentialCondition&);
\end{lstlisting}
\vspace{.3cm}
To construct essential condition on the fly, the = operator of classes {\class Unknown}, {\class OperatorOnUnknown} and {\class LcOperatorOnUnknown} are overloaded:
\begin{lstlisting}[]{}
EssentialCondition LcOperatorOnUnknown::operator = (const Real &):
EssentialCondition LcOperatorOnUnknown::operator = (const Complex &);
EssentialCondition LcOperatorOnUnknown::operator = (const Function &);
EssentialCondition OperatorOnUnknown::operator = (const Real &);
EssentialCondition OperatorOnUnknown::operator = (const Complex &);
EssentialCondition OperatorOnUnknown::operator = (const Function &);
EssentialCondition Unknown::operator = (const Real &);
EssentialCondition Unknown::operator = (const Complex &);
EssentialCondition Unknown::operator = (const Function &);
\end{lstlisting}
\vspace{.3cm}
All these member functions are implemented in the file \emph{EssentialCondition.cpp}.\\

\displayInfos{
library=essentialCondition,
header=EssentialCondition.hpp,
implementation=EssentialCondition.cpp,
test={test\_EssentialCondition.cpp},
header dep={Operator.h, config.h, utils.h}
}

\section{The {\classtitle EssentialConditions} class}

The {\class EssentialConditions} class manages a list of {\class EssentialCondition}. It inherits from {\class list<EssentialConditions>} with no additional member data:
\vspace{.3cm}
\begin{lstlisting}[deletekeywords={[3] list}]
class EssentialConditions : public std::list<EssentialCondition>
\end{lstlisting}
\vspace{.3cm}
As it is a user class, it mainly interfaces function of {\class list} class
\begin{lstlisting}[]{}
//constructors and assign
EssentialConditions() {}
EssentialConditions(const EssentialCondition&);
EssentialConditions(const EssentialCondition&, const EssentialCondition&);
EssentialConditions(const EssentialConditions&);
EssentialConditions& operator = (const EssentialConditions&);
//utilities
bool coupledUnknowns() const;
std::set<const Unknown*> unknowns() const;
const Unknown* unknown() const;
void print(std::ostream&) const;
friend std::ostream& operator << (std::ostream&, const EssentialConditions &);
\end{lstlisting}
\vspace{.3cm}
It overloads the \& operator to concatenate {\class EssentialCondition}:
\begin{lstlisting}[]{}
EssentialConditions operator & (const EssentialCondition&, const EssentialCondition&);
EssentialConditions operator & (const EssentialCondition&, const EssentialConditions&);
EssentialConditions operator & (const EssentialConditions&, const EssentialCondition&);
EssentialConditions operator & (const EssentialConditions&, const EssentialConditions&);
\end{lstlisting}
\displayInfos{
library=essentialCondition,
header=EssentialConditions.hpp,
implementation=EssentialConditions.cpp,
test={test\_EssentialCondition.cpp},
header dep={EssentialCondition.hpp, Operator.h, config.h, utils.h}
}

\section{The {\classtitle Constraints} class}

The {\class Constraints} class handles the algebraic representation $\mathbb{C}U=G$ of a set of essential conditions ({\class EssentialConditions}) :

\begin{lstlisting}[deletekeywords={[3] map}]
class Constraints
{
 protected :
 MatrixEntry* matrix_p;                   //constraints matrix
 VectorEntry* rhs_p;                      //constraints right hand side
 public :
 EssentialConditions conditions_;         //list of essential conditions
 std::vector<DofComponent> cdofsr_;       //row DofComponent's of matrix_p
 std::vector<DofComponent> cdofsc_;       //col DofComponent's of matrix_p
 std::map<DofComponent, Number> elcdofs_; //map of eliminated cdofs
 std::map<DofComponent, Number> recdofs_; //map of reduced cdofs
 bool reduced;                            //true if reduced system
 bool local;                              //true if local conditions
 bool symmetric;                          //true if keep symmetry
 bool isId;                               //true if Id matrix
}
\end{lstlisting}
\vspace{.3cm}
Even if unknowns involved in constraints system are vector unknowns, the representation uses scalar unknown numbering (\emph{cdofsc\_} vector). The \emph{cdofsr\_} handles row numbering, each row corresponding to one constraint equation linked to an "artificial" test function. The \emph{elcdofs\_} map carries the eliminated index set ${\cal E}$ and \emph{recdofs\_} map carries the reduced index set ${\cal R}$ . {\class Constraints} is processed in two steps :
\begin{itemize}
\item construction of the constraint system from the list of essential conditions : each condition are translated in its matrix representation, then all the matrix representations are merged in a large one
\item reduction of the large constraints system using QR algorithm with removing the redundant constraints and identifying the conflicting constraints
\end{itemize}
A constraint is local (\emph{local = true}) if it involves only one dof. Dirichlet conditions are local and transmission or periodic conditions are not.
The class proposes mainly the constructor from a list of essential conditions and usual copy and delete stuff (full copy):
\begin{lstlisting}[]{}
Constraints(MatrixEntry* =0, VectorEntry* =0);
Constraints(const EssentialCondition&);
Constraints(const Constraints&);
Constraints& operator=(const Constraints&);
void copy(const Constraints&);
~Constraints();
void clear();
\end{lstlisting}
\vspace{.3cm}
some accessors and property functions :
\begin{lstlisting}[deletekeywords={[3] set}]
const MatrixEntry* matrixp() const;
const VectorEntry* rhsp() const;
bool coupledUnknowns() const;
const Unknown* unknown() const;
std::set<const Unknown*> unknowns() const;
\end{lstlisting}
\vspace{.3cm}
The fundamental functions are the building functions
\begin{lstlisting}[]{}
bool createDirichlet(const EssentialCondition&);//create from Dirichlet cond.
void createLf(const EssentialCondition&);        //create from linear form cond.
void createNodal(const EssentialCondition&);     //create general condition
MatrixEntry* constraintsMatrix(const OperatorOnUnknown&, const GeomDomain*,
       Space*, const complex_t&, vector<Point>&, const Function*,
       vector<DofComponent>&);             //create submatrix of a constraint
void concatenateMatrix(MatrixEntry&, std::vector<DofComponent>&, const MatrixEntry&,
       const std::vector<DofComponent>&);  //concatenate 2 constraints matrices
void reduceConstraints(real_t aszero);     //reduce constraints matrix
template <typename T>
void buildRhs(const Function*, vector<Point>&, const T&);  //tool building rhs
\end{lstlisting}
\vspace{.3cm}
Functions \emph{createXXX(const EssentialCondition\&)} build the matrix representation for \emph{XXX} condition. Now, there exist Dirichlet, transmission and periodic type. It is quite easy to add new translation functions to take into account new ones. The \emph{reduceConstraints} function processes the reduction of any constraints sytem using QR algorithm. The \emph{asZero} variable is used to round to zero very small values induced by rounding errors in QR algorithm. This function modifies the current matrix and right hand side vector. So the original representation is lost after reduction!\\

To build the matrix representation of a set of essential conditions, two cases have to be considered:
\begin{itemize}
\item no condition that couples two unknowns (independent conditions); in that case one constraints system is built for each unknown
\item there exists one coupling condition (coupled conditions); in that case a unique constraints system for all unknowns is built
\end{itemize}
For instance, the conditions ($u_1$, $u_2$ two unknowns):
$$ u_1=0|\Sigma\ \  u_1=0|\Gamma\ \  u_2=0|\Gamma$$
gives two constraints system, one merging $u_1$ conditions and an other one for $u_2$. While the conditions :
$$ u_1=0|\Sigma\ \  u_1=0|\Gamma\ \  u_2-u_1=g|\Gamma$$
gives a global constraints system merging all conditions.\\

The merging and building process is done by an external functions that return a list of {\class Constraints} indexed by unknown ({\class Unknown} pointer). In case of coupled conditions, the list contains only one {\class Constraints} indexed by 0.
\begin{lstlisting}[deletekeywords={[3] map}]
map<const Unknown*, Constraints*> mergeConstraints(vector<Constraints*>&);
map<const Unknown*, Constraints*> buildConstraints(const EssentialConditions&);
\end{lstlisting}
\vspace{.3cm}
The \emph{buildConstraints} function looks like (iterator declaration omitted):
\begin{lstlisting}[deletekeywords={[3] map}]
map<const Unknown*, Constraints*> buildConstraints(const EssentialConditions& ecs){
map<const Unknown*, Constraints*> mconstraints
vector<Constraints *> constraints(ecs.size());
itc=constraints.begin();
//translate conditions into constraints
for(ite=ecs.begin(); ite!=ecs.end(); ite++, itc++)
	  *itc=new Constraints(*ite);
//merge constraints
mconstraints =  mergeConstraints(constraints);
//reduce constraints
for(itm=mconstraints.begin(); itm!=mconstraints.end(); itm++)
	  itm->second->reduceConstraints();
return mconstraints;}
\end{lstlisting}
\vspace{.3cm}
The {\class Constraints} class provides important tools that take into account essential conditions in linear systems:
\begin{lstlisting}[deletekeywords={[3] map}]
void pseudoColReduction(MatrixEntry*, vector<DofComponent>&,
                        vector<DofComponent>&, MatrixEntry*&);
void pseudoRowReduction(MatrixEntry*, vector<DofComponent>&,
                        vector<DofComponent>&);
friend void extendStorage(MatrixEntry*, vector<DofComponent>&,
            vector<DofComponent>&, const Constraints*, const Constraints*);
void extractCdofs(const vector<DofComponent>&, map<DofComponent, Number>&,
                  map<DofComponent, Number>&, bool useDual=false) const;
friend void appliedRhsCorrectorTo(VectorEntry*, vector<DofComponent>&,
             MatrixEntry*, const Constraints*, const Constraints*,
             ReductionMethod);
\end{lstlisting}
\vspace{.3cm}
The functions \emph{pseudoColReduction} and \emph{pseudoRowReduction} performs the reduction of system matrix according to relation \ref{mat_el}. The \emph{pseudoColReduction} also stores the part of matrix system to be memorized to process right hand side correction. If case of non homogeneous constraints, the function \emph{appliedRhsCorrectorTo} modifies the system right hand side according to relation \ref{rhs_el}. Regarding unknown and test function constraints, \emph{extendStorage} function extends the storage of a matrix. This extension occurs when constraints are not local. \emph{extractCdofs} function is a facility to localized eliminated or reduced rows/columns in matrix. All of these functions are called by compute functions of {\class TermMatrix} class.\\

In order to deal with constraints acting on different matrix blocs, the {\class SetOfConstraints} class is provided. It manages mainly a map relating unknown and constraint:
\begin{lstlisting}[deletekeywords={[3] map}]
class SetOfConstraints : public map<const Unknown*, Constraints*>
{
public:
SetOfConstraints(const Unknown*, Constraints*);  //from a constraints pointer
SetOfConstraints(const map<const Unknown*, Constraints*>& mc);   //from pointers
SetOfConstraints(const EssentialConditions&);    //from EssentialConditions
SetOfConstraints(const SetOfConstraints&); //copy constructor with full copy
SetOfConstraints& operator = (const SetOfConstraints&);//assignment with copy
~SetOfConstraints();    //destructor deallocating constraints pointers
void clear();           //deallocate constraints pointer and reset the map

Constraints* operator()(const Unknown*) const; //acces operator from unknown
bool isGlobal() const;              //true if a unique global constraint
bool isReduced() const;             //true if  all constraints have been reduced
void print(std::ostream&) const;    //print utility
friend std::ostream& operator<<(std::ostream&, const SetOfConstraints&);
};
\end{lstlisting}
\vspace{.2cm}
The {\class SetOfConstraints} objects are very important because they are used by the {\class TermMatrix} objects when they have to deal with essential conditions!

\vspace{.1cm}
\displayInfos{
library=essentialCondition,
header=Constraints.hpp,
implementation=Constraints.cpp,
test=test\_Constraints.cpp,
header dep={EssentialCondition.hpp, EssentialConditions.hpp, LargeMatrix.h, config.h, utils.h}
}