Class xlifepp::Constraints#

class 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)

void buildFrom(const EssentialCondition &ec)#: create constraints from an essential condition

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

std::vector<std::vector<std::pair<DofComponent, complex_t>>> constraintRows() const#

rows of C=[I F] as sparse lists (cdof, coefficient)

rows of the reduced constraints matrix C=[I F] (constraints read U_E + F*U_R = g, see reduceConstraints): row i (i=1..nbe, rank of the eliminated cdof e_i) -> [(e_i,1), (r_j,F(i,j)) for non zero F(i,j)]

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:

 [ ... a*op(w_i1)(Mk)... a*op(w_i2)(Mk) .]  for node Mk

the matrix is a scalar one and always stored in column sparse (csc) built also the column component dofs vector

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)

bool coupledUnknowns() const#: true if one essential condition couples at least two (parent) 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

inline number_t numberOfConstraints() const#: number of constraints

Constraints &operator=(const Constraints&)#: assign operator (full copy)

void penalizationReduction(MatrixEntry*, std::vector<DofComponent>&, std::vector<DofComponent>&, const complex_t&)#

penalization reduction

symmetric penalization reduction, assuming same constraints on u and v the reduced constraints read U_E + F*U_R = g, say C*U = g with C=[I F]; the penalized problem reads (A + alpha C’C) U = B + alpha C’g (C’ = adjoint of C) this function adds alpha C’C to the matrix mat, i.e.

alpha I to (E,E), alpha F to (E,R), alpha F’ to (R,E) and alpha F’F to (R,R) (for a Dirichlet condition F=0: alpha on the diagonal of the eliminated dofs); the right hand side is corrected by appliedRhsCorrectorTo; the storage has to be extended before (extendStorage with doPenal=true). With a symmetric storage, only the lower part is updated.

void print(std::ostream&) const#: print utillity

void pseudoColReduction(MatrixEntry*, std::vector<DofComponent>&, std::vector<DofComponent>&, MatrixEntry*&)#: pseudo reduction of columns

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

                   u1           u2         ...      un
             --------------------------------------------
             | ... ... ... | ... ... |           |  ... |          | . |

original C = | … C1 … | … … | … | Cn | D= | . | constraints | … … … | … … | | … | | . |

  rearrangement of u1 , u2, ..., un

| 1 x … | … … | | … …| | . | reduce rC = | 0 1 x | … … | … | … …| rD= | . | constraints | 0 0 … | … … | | … …| | . |

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

                              part of  u1 , u2, ..., un
                          --------------------------------
                          |  x   x  |           |  x   x |          | . |

residual constraints F = | x … | … | x …| b = | . | | … … | | … …| | . |

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

void toParentCdofs()#

cdofs of components of unknown -> cdofs of the parent unknown

cdofs (u_k, dof, 1) of a component u_k=u[k] of a vector unknown u -> cdofs (u, dof, k) of the parent unknown, so that the constraints match the numbering of the matrices (ComponentsOfUnknown u_K, K={k1,k2,…}: (u_K, dof, j) -> (u, dof, kj))

const Unknown *unknown() const#: return (parent) unknown of constraint, 0 if more than one unknown

std::set<const Unknown*> unknowns() const#: return (parent) unknowns of constraint

ValueType valueType() const#: return value type of set of constraints (i.e involved real or complex)

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

const Unknown *lambda_p = nullptr#: Lagrange multiplier (dual reduction): one component by reduced constraint (not owned)

Space *lambdaSpace_p = nullptr#: virtual space of the Lagrange multiplier (dual reduction, not owned)

bool local#: true if all conditions are local: only interaction with dofs in same element

const Unknown *mu_p = nullptr#: test function associated to the Lagrange multiplier (dual reduction, not owned)

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