Detailed Description

template<typename T>
class gismo::gsPrimalSystem< T >

This class represents the primal system for a IETI-DP algorithm.

The class gsIetiSystem does not know anything about the primal problem. For that class the primal problem is just another subdomain.

The class at hand allows to handle primal degrees of freedom (dofs). It is purely algebraic and assumes that the caller provides the constraint matrix C (or, more precisely, its individual rows) and a mapping (as vector of indices) that tells the class to which primal dof the corresponding constraint belongs. So, for example, if a primal dof is a vertex value, then that vertex gets one identifying index identifying it as one single primal dof, while there is a constraint for each patch the corner belongs to.

It is assumed that the class receives the matrix \( A_k \) and the vectors \( c_k^{(n)} \) that make up the constraints matrix \( C_k \) such that the local saddle point formulation reads as follows:

\[ \tilde A_k = \begin{pmatrix} A_k & C_k^\top \\ C_k & \end{pmatrix} \]

Using setEliminatePointwiseConstraints, it is possible to instruct the class to eliminate pointwise constraints. A constraint is pointwise iff the vector \( c_k^{(n)} \) only has one non-zero entry accordning to the sparsity pattern of the sparse vector data structure. This is typically the case for vertex values.

Usually, it will be the best only to call handleConstraints, which modifies the local system matrix, the jump matrix and the right-hand side such that they can then handed over to the gsIetiSystem. That function simultainously collects the contributions for the primal problem, which is handed over to gsIetiSystem as last subspace.

The member handleConstraints is based on the static members incorporateConstraints, primalBasis and the member addContribution. incorporateConstraints sets up the matrix \( \tilde A_k \), \( R_c \) (right-hand-side for basis), \( P_c \) (embedding for basis) and \( P_0 \) (local embedding). The member primalBasis constructs the basis for the primal space by the principle of energy minimization by calculating \( \Psi = P_c^\top \tilde{A}_k^{-1} R_c \). The member addContribution adds the local contributions to the primal problem. Finally, the matrix \( P_0 \) is used to modify the jump matrix and the local rhs.

The member handleConstraints combines the member incorporateConstraints, the setup of a sparse LU solver for the local system, and the members primalBasis and addContribution into one single function.

After going through all patches, the class gsIetiSystem gets one more subdomain: the primal one (which is set up based on jumpMatrix, localMatrix and localRhs).

After solving, the member distributePrimalSolution distributes the solution obtained for the primal problem back to the individual patches.

Collaboration diagram for gsPrimalSystem< T >:

Public Member Functions
void	addContribution (const JumpMatrix &jumpMatrix, const SparseMatrix &localMatrix, const Matrix &localRhs, SparseMatrix primalBasis, OpPtr embedding=OpPtr())
	Adds contributions for a patch to the data hold in the class. More...

std::vector< Matrix >	distributePrimalSolution (std::vector< Matrix > sol)
	Distributes the given solution for K+1 subdomains to the K patches. More...

index_t	eliminatePointwiseConstraints () const
	Returns true iff handleConstraints will eliminate pointwise constraints (typically vertex values)

	gsPrimalSystem (index_t nPrimalDofs)
	Constructor. More...

void	handleConstraints (const std::vector< SparseVector > &primalConstraints, const std::vector< index_t > &primalDofIndices, JumpMatrix &jumpMatrix, SparseMatrix &localMatrix, Matrix &localRhs)
	Convenience function for handling the primal constraints. More...

JumpMatrix &	jumpMatrix ()
	Returns the jump matrix for the primal problem.

SparseMatrix &	localMatrix ()
	Returns the local stiffness matrix for the primal problem.

Matrix &	localRhs ()
	Returns the right-hand-side for the primal problem.

index_t	nPrimalDofs () const
	Returns the size of the primal problem (number of primal dofs)

void	setEliminatePointwiseConstraints (bool v)
	Iff true, handleConstraints will eliminate pointwise constraints (typically vertex values)

Static Public Member Functions
static void	incorporateConstraints (const std::vector< SparseVector > &primalConstraints, bool eliminatePointwiseConstraints, const SparseMatrix &localMatrix, SparseMatrix &modifiedLocalMatrix, SparseMatrix &localEmbedding, SparseMatrix &embeddingForBasis, Matrix &rhsForBasis)
	Incorporates the given constraints in the local system. More...

static gsSparseMatrix< T >	primalBasis (OpPtr localSaddlePointSolver, const SparseMatrix &embeddingForBasis, const Matrix &rhsForBasis, const std::vector< index_t > &primalDofIndices, index_t nPrimalDofs)
	Returns the matrix representation of the energy minimizing primal basis. More...

Private Types
typedef gsSparseMatrix< T, RowMajor >	JumpMatrix
	Sparse matrix type for jumps.

typedef gsMatrix< T >	Matrix
	Matrix type.

typedef gsLinearOperator< T >	Op
	Linear operator.

typedef memory::shared_ptr< Op >	OpPtr
	Shared pointer to linear operator.

typedef gsSparseMatrix< T >	SparseMatrix
	Sparse matrix type.

typedef gsSparseVector< T >	SparseVector
	Sparse vector type.

Private Attributes
bool	m_eliminatePointwiseConstraints
	handleConstraints will eliminate pointwise constraints

std::vector< OpPtr >	m_embeddings
	For each patch, the map \( \tilde u_k \) to \( u_k \).

JumpMatrix	m_jumpMatrix
	The jump matrix for the primal problem.

SparseMatrix	m_localMatrix
	The overall matrix for the primal problem.

Matrix	m_localRhs
	The right-hand side for the primal problem.

std::vector< SparseMatrix >	m_primalBases
	The bases for the primal dofs on the patches.

Constructor & Destructor Documentation

gsPrimalSystem ( index_t nPrimalDofs )

Constructor.

Parameters

nPrimalDofs Number of primal constraints in total

Member Function Documentation

void addContribution	(	const JumpMatrix &	jumpMatrix,
		const SparseMatrix &	localMatrix,
		const Matrix &	localRhs,
		SparseMatrix	primalBasis,
		OpPtr	embedding = `OpPtr()`
	)

Adds contributions for a patch to the data hold in the class.

Parameters

jumpMatrix	Jump matrix \( B_k \)
localMatrix	Local stiffness matrix \( A_k \)
localRhs	Local right-hand side \( f_k \)
primalBasis	Matrix representation of the primal basis
embedding	The map \( \tilde u_k \) to \( u_k \); if not provided, defaulted to (possibly rectangular) identity

std::vector< typename gsPrimalSystem< T >::Matrix > distributePrimalSolution ( std::vector< Matrix > sol )

Distributes the given solution for K+1 subdomains to the K patches.

Parameters

sol The solution, first for the K patches, followed by the contribution for the primal dofs. The solutions for the K patches is expected to have first the values for all patch-local degrees of freedom, possibly followed by degrees of freedom from Lagrange-mutlipliers used for enforcing primal constraints.

Returns: The solution for the K patches

void handleConstraints	(	const std::vector< SparseVector > &	primalConstraints,
		const std::vector< index_t > &	primalDofIndices,
		JumpMatrix &	jumpMatrix,
		SparseMatrix &	localMatrix,
		Matrix &	localRhs
	)

Convenience function for handling the primal constraints.

Parameters

[in]	primalConstraints	Primal constraints; the given vectors make up the matrix \( C_k \)
[in]	primalDofIndices	Vector, that contains for every primal constraint the index of the respective primal dof
[in,out]	jumpMatrix	Jump matrix \( B_k \)
[in,out]	localMatrix	Local stiffness matrix \( A_k \)
[in,out]	localRhs	Local right-hand side \( f_k \)

The implementation is basically:

gsSparseMatrix<T> modifiedLocalMatrix, localEmbedding, embeddingForBasis;
gsMatrix<T> rhsForBasis;
gsPrimalSystem<T>::incorporateConstraints(
    primalConstraints, this->eliminatePointwiseConstraints(), localMatrix,
    modifiedLocalMatrix, localEmbedding, embeddingForBasis, rhsForBasis);
this->addContribution(
    jumpMatrix, localMatrix, localRhs,
    gsPrimalSystem<T>::primalBasis(
        makeSparseLUSolver(modifiedLocalMatrix),
        embeddingForBasis, rhsForBasis, primalDofIndices, this->nPrimalDofs()
    )
);
localMatrix  = give(modifiedLocalMatrix);
localRhs     = localEmbedding * localRhs;
jumpMatrix   = jumpMatrix * localEmbedding.transpose();

void incorporateConstraints	(	const std::vector< SparseVector > &	primalConstraints,
		bool	eliminatePointwiseConstraints,
		const SparseMatrix &	localMatrix,
		SparseMatrix &	modifiedLocalMatrix,
		SparseMatrix &	localEmbedding,
		SparseMatrix &	embeddingForBasis,
		Matrix &	rhsForBasis
	)

static

Incorporates the given constraints in the local system.

Parameters

[in]	primalConstraints	Primal constraints; the given vectors make up the matrix \( C_k \)
[in]	eliminatePointwiseConstraints	True iff pointwise constraints (typically vertex values are eliminated (and not treated as local saddle point problems). A constraint is pointwise iff \( C_k \) has one non-zero entry.
[in]	localMatrix	Local matrix \( A_k \)
[out]	modifiedLocalMatrix	Local matrix \( \tilde A_k \)
[out]	localEmbedding	Embedding \( P_0 \). This matrix is used by to modify the right-hand side and the jump matrix (cf. handleConstraints)
[out]	embeddingForBasis	Embedding \( P_c \) for energy-minimizing basis (cf. primalBasis)
[out]	rhsForBasis	Right-hand side \( R_c \) for energy-minimizing basis (cf. primalBasis)

Iff the constraint is not handled by elimination, the output is:

\[ \tilde{A}_k = \begin{pmatrix} A_k & C_k^\top \\ C_k & 0 \end{pmatrix} \quad\text{and}\quad P_0 = P_c = \begin{pmatrix} I & 0 \end{pmatrix} \quad\text{and}\quad R_c = \begin{pmatrix} 0 \\ I \end{pmatrix}. \]

Iff the constraint is handled by elimination and if we assume that the first block component of

\[ A_k = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \]

is to be eliminated, then we have:

\[ \tilde{A}_k = \begin{pmatrix} I \\ & A_{22} \end{pmatrix} \quad\text{and}\quad P_0 = \begin{pmatrix} 0 \\ & I \end{pmatrix} \quad\text{and}\quad P_c = \begin{pmatrix} I \\ & I \end{pmatrix} \quad\text{and}\quad R_c = \begin{pmatrix} I \\ -A_{21} \end{pmatrix}. \]

gsPrimalSystem< T >::SparseMatrix primalBasis	(	OpPtr	localSaddlePointSolver,
		const SparseMatrix &	embeddingForBasis,
		const Matrix &	rhsForBasis,
		const std::vector< index_t > &	primalDofIndices,
		index_t	nPrimalDofs
	)

static

Returns the matrix representation of the energy minimizing primal basis.

Parameters

localSaddlePointSolver	Solver that realizes \( \tilde{A}_k^{-1} \)
embeddingForBasis	Embedding \( P_c \) as provided by incorporateConstraints
rhsForBasis	Embedding \( R_c \) as provided by incorporateConstraints
primalDofIndices	Vector, that contains for every primal constraint the index of the respective primal dof
nPrimalDofs	The total number of primal dofs (nPrimalDofs())

The local basis is given by \( \Psi = P_c^\top \tilde{A}_k^{-1} R_c \).

Returns: a version of \( \Psi \) where the column indices are changed according to primalDofIndices.

Detailed Description

template<typename T> class gismo::gsPrimalSystem< T >

Public Member Functions

Static Public Member Functions

Private Types

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

template<typename T>
class gismo::gsPrimalSystem< T >