Class xlifepp::Constraints#
-
class Constraints#
-
Collaboration diagram for xlifepp::Constraints:
class handles the matrix representation of xlifepp::EssentialConditions
Public Functions
-
Constraints(const Constraints&)#
-
copy constructor (full copy)
-
Constraints(const EssentialCondition&)#
-
construct Constraints system from a EssentialCondition
-
Constraints(MatrixEntry* = nullptr, VectorEntry* = nullptr)#
-
default constructor
-
~Constraints()#
-
destructor (delete allocated pointers)
-
template<typename T>
void buildRhs(const Function*, const std::vector<Point>&, const std::vector<Vector<real_t>>&, const T&)#
-
tool to build rhs constraint
-
void clear()#
-
clear all (reset)
-
void concatenateMatrix(MatrixEntry&, std::vector<DofComponent>&, const MatrixEntry&, const std::vector<DofComponent>&)#
-
concatenate 2 constraints matrix
concatenate 2 constraint matrices mat1 and mat2 with no common column dofs mat1 and mat2 are scalar matrix stored in cs column, mat1 and cdofsc1 are updated
-
MatrixEntry *constraintsMatrix(const OperatorOnUnknown&, const GeomDomain*, Space*, const complex_t&, const std::vector<Point>&, const Function*, std::vector<DofComponent>&)#
-
create submatrix of a constraint
create constraints matrix from a term a*op(u)|nodes where nodes is a collection of nodes:
the matrix is a scalar one and always stored in column sparse (csc) built also the column component dofs vector[ ... a*op(w_i1)(Mk)... a*op(w_i2)(Mk) .] for node Mk
opu : operator on unknown dom : domain where opu is restricted spu : space or subspace used for interpolation (may be different from space of unknown) coeff : coefficient applied to operator nodes : list of points fmap : mapping function of nodes, may be 0 meaning id mapping
cdofsc : component dofs vector built by this routine
-
void copy(const Constraints&)#
-
copy function (full copy)
-
inline bool coupledUnknowns() const#
-
true if one essential condition couples at least two unknowns
-
void createArgyris(const EssentialCondition &ec)#
-
create u=0 constraints for Argyris
create constraints for condition u=0 on gamma in case of Argyris approximations u=0 on gamma -> u_i=0 but also null tangential derivatives at vertices: implies dy.u*nx-dx.u*ny = 0 -> dxu_i*nix-dyu_i*niy = 0 dyy.u*nx*nx+dxx.u*ny*ny-2dxy.u*nx*ny =0 -> dyyu_i*nix*nix+dxxu_i*niy*niy-2dxyu_i*nix*niy =0 with u_i,dxu_i,dyu_i, … dof on vertex i NOTE: because side element are not defined for Argyris, we use more complex algorithm that travels geometric side elements to retry dof located on side gamma and then create constraints on value and derivative dofs because tangential derivatives at vertex may differ from side elements sharing vertex, both tangential derivative constraints are added and the reduction process will average them.
This average process may induce non local constraints and thus matrix storage change
-
bool createDirichlet(const EssentialCondition&)#
-
create Dirichlet condition like
construct constraints matrix in case of Dirichlet condition
Lagrange scalar/vector unknown u: a*u = 0/g Raviart-Thomas (3D), Nedelec face (3D) scalar unknown: a*u.n = 0/g Nedelec (2D), Nedelec Edge (3D) scalar unknown: a*u^n = 0 Morlay scalar unknown: a*u=g or grad(u).n=g
build Id matrix and rhs vector 0/g
-
void createLf(const EssentialCondition&)#
-
create from linear form condition
create constraint data for a lf condition lf(u)=c (one row constraint), lf represented by a TermVector v the constraints matrix has to be a scalar matrix.
It means that in case of vector unknown the constraints are split
scalar unknown case C = | v1 ... vn | row matrix vector unknown(2D) case C = | [v11 v12] [v21 v22] ...| -> | v11 v12 v21 v22 ....|
-
void createNodal(const EssentialCondition&)#
-
create general condition (nodal)
General method to construct “nodal” constraint a1*op1(u1)|d1 + a2*op2(u2)|d2 + … = g where ai are complex scalars, opi are operators on unknowns ui and di are geometric domain.
only one ui,di : u|d = g Dirichlet condition
ui may be the same: u|d1 - u|d2 = g periodic condition
di may be the same: u1|d - u2|d = g transmission condition
ui,di are different: u1|d1 - u2|d2 = g generalized transmission condition When domains are different it is mandatory to give the maps mi relating d1 to di (see defineMap function), so the condition reads a1*op1(u1)|d1 + a2*op2(u2)|m2(d1) + … = g In any case, some nodes of d1 drive the computation: a1*op1(u1)(Mk) + a2*op2(u2)(m2(Mk)) + … = g(Mk) Mk in d1 The nodes Mk depends on the status of domain d1 : 1) if domain d1 belongs to the mesh supporting u1 and u1 is a Lagrange interpolation then nodes used are all the nodes of the FE interpolation belonging to domain d1 2) if domain d1 belongs to the mesh supporting u1 and u1 is a Nedelec or Raviart interpolation then nodes used are all the virtual coordinates of the dofs belonging to domain d1 3) else nodes defining geometric domain are used
Be cautious: u1|d1 + u2|d2 = g is not equivalent to u2|d2 + u1|d1 = g
The algorithm is the following
construct the list of nodes (Nk) of d1 regarding case
build matrix ai*opi(ui)|mi(di) for i=1,nbterms (see constraintsMatrix function)
concatenate them (see concatenateMatrix function)
construct right hand side
Note
the constraints matrix has to be a scalar matrix. It means that in case of vector unknown the constraints are split
| a 0 ...| Dirichlet and scalar unknown case C = | 0 ... 0 | diagonal matrix |... 0 a | | [a 0] [0 0] ... ...| | a 0 0 0 ... ...| | [0 a] [0 0] ... ...| | 0 a 0 0 ... ...| Dirichlet and vector unknown(2D) case C = | [0 0] [a 0] [0 0] ...| -> | 0 0 a 0 0 0 ...| | [0 0] [0 a] [0 0] ...| | 0 0 0 a 0 0 ...| | [0 0] [0 0] [a 0] ...| | 0 0 0 0 a 0 ...| | ... ... ... ...| | ... ... ... ...|
-
void extractCdofs(const std::vector<DofComponent>&, std::map<DofComponent, number_t>&, std::map<DofComponent, number_t>&, bool useDual = false) const#
-
extract eliminated and reduced cdofs from list of cdofs
extract eliminated and reduced cdofs from list of cdofs cdofs: list of cdofs melcdofs: list of elimited cdofs with their ranks in cdofs list mrecdofs: list of reduced cdofs with their ranks in cdofs list useDual : if true use either unknown or dual unknown to get cdofs
-
inline const MatrixEntry *matrixp() const#
-
accessor to matrix
-
Constraints &operator=(const Constraints&)#
-
assign operator (full copy)
-
void penalizationReduction(MatrixEntry*, std::vector<DofComponent>&, std::vector<DofComponent>&, const complex_t&)#
-
penalization reduction
penalization reduction add alpha(U_E + F*U_R) to the block (Ve,Ue+Ur) of the matrix mat, assuming same constraints on u and v works only for Dirichlet condition (no reduced unknown component)
-
void print(std::ostream&) const#
-
print utillity
-
void pseudoColReduction(MatrixEntry*, std::vector<DofComponent>&, std::vector<DofComponent>&, MatrixEntry*&)#
-
pseudo reduction of columns
perform pseudo reduction of reduced essential conditions in a matrix.
Essential conditions have the following form : U_E + F*U_R = H for column unknown U V_E + G*V_R = 0 for row test function V (generally related to unknown U) where E are the eliminated unknown/test function indices and R are the reduced unknown/test function indices other indices S correspond to dofs not connected to the essential conditions dofs The pseudo reduction of matrix consists in
modifying column A.j for j in R by adding -Fkj*A.k for k in E and replacing column A.k for k in E by a 0 column
modifying row Ai. for i in R by adding -Gki*Ak. for k in E and replacing row Ak. for k in E by a 0 row
if eliminated v unknowns are dual of eliminated u unknowns (G=F), the VE rows are replaced by the equation (or a part of) a*U_E + a*F*U_R = a*H where a is a given scalar to delay the right hand side modification, the (A.k) columns (k in E) are stored in a new structure
At the end of the process, the eliminated system looks like (C the correction matrix mxE get from reduction of the constaints)
If F=G=0 (Dirichlet condition) , R={} and the system readsU_E U_R U_S U RHS ---------------------------- ------- ------- | | | | | | | | V_E | (a+b)*Id | a*F | 0 | | U_E | | a*H | => (a+b)*U_E + a*F*U_R = a*H (mimics the constraints) | | | | | | | | ---------------------------- ------- ------- | | | | | | | | V_R | 0 | A_RR | A_RS | * | U_R | = | B_R | = B0_R - C_RE*H -Gt*B0_E (rhs correction on reduced dof) | | | | | | | | ---------------------------- ------- ------ | | | | | | | | V_S | 0 | A_SR | A_SS | | U_S | | B_S | = B0_S - C_SE*H (rhs correction on free dofs) | | | | | | | | ---------------------------- ------- -------
In some cases (F=G=0 or non dof coupling condition…) the storage is not modified but in other cases (transmission condition for instance) the storage is modifiedU_E U_S U RHS -------------------- ------- ------- | | | | | | | V_E | (a+b)*Id | 0 | | U_E | | a*H | => (a+b)*U_E = a*H (constraints) | | | | | | | ------------------- * ------- = ------- | | | | | | | V_S | 0 | A_SS | | U_S | | B_S | = B0_S - C_SE*H (rhs correction) | | | | | | | -------------------- ------- ------- for homogeneous dirichlet condition H=0 and B_S = B0_S!
choosing a=0,b!=0 allows to deal with eigen value problems by shifing the spectra corresponding to eliminated dof
mat : matrix to be eliminated cdofr: row component dofs of matrix cdofc: col component dofs of matrix rhsmat: right hand side matrix pointer
-
void pseudoRowReduction(MatrixEntry*, std::vector<DofComponent>&, std::vector<DofComponent>&, const complex_t&, const complex_t&, bool)#
-
pseudo reduction of rows
-
void reduceConstraints(real_t aszero = 10. * theEpsilon)#
-
reduce constraints matrix to a set of reduced constraints
reduce constraints system to an upper triangular system using stable QR factorization transform the constraints into a new constraints where matrix is an upper triangular matrix
original C = | … C1 … | … … | … | Cn | D= | . | constraints | … … … | … … | | … | | . |u1 u2 ... un -------------------------------------------- | ... ... ... | ... ... | | ... | | . |
| 1 x … | … … | | … …| | . | reduce rC = | 0 1 x | … … | … | … …| rD= | . | constraints | 0 0 … | … … | | … …| | . |rearrangement of u1 , u2, ..., un
when algorithm ‘fails’ it means that some constraints are redundant or conflicting, so we eliminate on the fly some line constraints
In a second step, the upper triangular system is inverted to produce residual rectangular matrix and modified right hand side
residual constraints F = | x … | … | x …| b = | . | | … … | | … …| | . |part of u1 , u2, ..., un -------------------------------- | x x | | x x | | . |
The first p column indices (of squared upper triangular matrix) corresponds to the eliminated dofs (say U_e) while the other column indices corresponds to the residual dofs (say U_r). Finally the constraints system reads
| U_e + F*U_r = b <=> U_e= b - F*U_r |
which is the useful expression to deal with essential conditions in bilinear form computation
this process may be permutes column cdofs and produces elcdofs_ : the map of eliminated cdofs, cdof -> k (rank in cdofsc_) recdofs_ : the map of reduced cdofs: cdof -> k (rank in cdofsc_)
- Parameters:
-
aszero – : a positive value near 0 used to round to 0 very small values, say |v|<aszero, default value is 10*theEpsilon
Note
F may be a null matrix (classic Dirichlet condition for instance, condition reads U_e = b as usual)
-
inline const VectorEntry *rhsp() const#
-
accessor to rhs
Public Members
-
std::vector<DofComponent> cdofsc_#
-
col DofComponent’s of matrix_p (unknown cdofs involved in constraints)
-
std::vector<DofComponent> cdofsr_#
-
row DofComponent’s of matrix_p (constraint dofs, virtual cdofs)
-
EssentialConditions conditions_#
-
list of Essential conditions (analytic representation)
-
std::map<DofComponent, number_t> elcdofs_#
-
map of eliminated cdofs: cdof -> k (rank in cdofsc_)
-
bool isId#
-
true if matrix is Id matrix
-
bool local#
-
true if all conditions are local: only interaction with dofs in same element
-
std::map<DofComponent, number_t> recdofs_#
-
map of reduced cdofs: cdof -> k (rank in cdofsc_)
-
bool reduced#
-
true means that the essential conditions have been reduced to a diagonal form
-
real_t redundancyTolerance_#
-
tolerance factor used in reduction process (elimination of redundant constraints)
-
bool symmetric#
-
true if all constraints keep symmetry
Friends
-
friend std::ostream &operator<<(std::ostream&, const Constraints&)#
-
print operator
-
Constraints(const Constraints&)#