From 6bba6c7360cc4cc58135ca52b4f58d4dcaf633fc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@math.fau.de>
Date: Sun, 16 Feb 2025 11:48:22 +0100
Subject: [PATCH] [constraints] Improve doc of AffineConstraints
Document what exactly `AffineConstraints` is supposed
to do and how it is intended to be used and rename
`constrainMatrixPattern()` to extendMatrixPattern()`
since this in fact enlarges the pattern and never
reduces it.
---
dune/fufem/constraints/affineconstraints.hh | 206 ++++++++++++--------
examples/utilities.hh | 4 +-
2 files changed, 126 insertions(+), 84 deletions(-)
diff --git a/dune/fufem/constraints/affineconstraints.hh b/dune/fufem/constraints/affineconstraints.hh
index 6dcf45b3..4310053e 100644
--- a/dune/fufem/constraints/affineconstraints.hh
+++ b/dune/fufem/constraints/affineconstraints.hh
@@ -29,27 +29,68 @@ namespace Dune::Fufem {
*
* This class handles affine constraints on a subset of DOFs of the form
*
- * \f[ x_i = \sum_{j} c_{ij} x_j + b_i\f]
+ * \f[ x_i = \sum_{j} L_{ij} x_j + c_i\f]
*
- * or in short form
+ * or in short \f$ x = L x + c\f$ for a square matrix
+ * \f$ L \in \mathbb{R}^{n \times n} \f$ and a vector \f$c \in \mathbb{R}^n\f$.
+ * The main purpose is to simplify the handling
+ * of problems in the constrained affine subspace
*
- * \f[ x = Cx + b\f]
+ * \f[ V_{L,c} = \{ x\,|\, x = L x+c\}.\f]
*
- * for a square matrix \f$C\f$ and a vector \f$b\f$.
* We assume that the indices \f$ i \f$ can be classified
- * into non-constrained ones, where \f$ c_{ij} = \delta_{ij} \forall j \f$
- * and \f$b_i=0\f$ and constrained ones, where \f$c_{ii}=c_{ij}=0 \f$
- * if \f$x_j\f$ is also constrained. Then any constrained \f$x_i\f$
- * can be interpolated from the non-constrained \f$ x_j, j\neq i\f$
- * and \f$b_i\f$.
+ * into non-constrained ones, where (using the Kronecker delta)
+ * \f$ L_{ij} = \delta_{ij} \forall j \f$
+ * and \f$c_i=0\f$ and constrained ones, where \f$L_{ii}=L_{ij}=0 \f$
+ * if \f$j\f$ is another constrained index. Then any constrained \f$x_i\f$
+ * can be interpolated from \f$c_i\f$ and \f$ x_j \f$
+ * for the non-constraints indices \f$j\neq i\f$.
+ * Under these assumptions we have \f$L^2 = L\f$ and \f$L c=0\f$
+ * such that the affine map \f$ \Psi:\mathbb{R}^n \to V_{L,c}\f$
+ * given by \f$v \mapsto \Psi(v) = L v +c\f$ is surjective projection
+ * into \f$V_{L,c}\f$.
*
- * Under these assumptions we furthermore have \f$C^2 = C\f$ and \f$Cb=0\f$
- * such that the affine subspace
+ * Consider a linear variational problem
*
- * \f[ V_{C,b} = \{ x\,|\, x = Cx+b\}\f]
+ * \f[ u \in V_{L,c}: \qquad \langle A u, v \rangle = \langle b, v \rangle \qquad \forall v \in V_{L,0} \f]
*
- * of all \f$x\f$ satisfying the affine constraints is given by the range
- * of the affine map \f$v \mapsto Cv+b\f$.
+ * in the affine subspace \f$ V_{L,c} = \{ L x +c \,|\, x\in \mathbb{R}^n\}\f$.
+ * Using of the decomposition \f$u= u_0 + c\f$ with \f$ u_0 = L u_0 = L u \in V_{L,0} \f$
+ * this can then be written as
+ *
+ * \f[ u_0 \in V_{L,c}: \qquad \langle A u_0, v \rangle = \langle b - A c, v \rangle \qquad \forall v \in V_{L,0}. \f]
+ *
+ * which, because \f$L:V_{L,0} \to V_{L,0}\f$ is a linear projection, is equivalent to
+ *
+ * \f[ u_0 \in \mathbb{R}^n: \qquad \langle L^T A L u_0, v \rangle = \langle L^T(b - A c), v \rangle \qquad \forall v \in \mathbb{R}^n \f]
+ *
+ * The matrix \f$ L^T A L\f$ is obviously singular, since it is
+ * zero on \f$ker(L) = (V_{L,0})^\perp\f$. Hence this
+ * system does not have a unique solution.
+ * Assuming that \f$A\f$ is regular on \f$V_{L,0}\f$, i.e.
+ * \f$ker (A) \subset ker(L)\f$ we know
+ * that the solution is at least unique in \f$V_{L,0}\f$.
+ * To regularize the problem we consider the reduced identity matrix
+ * \f$ \hat{I} \in \mathbb{R}^{n \times n}\f$ which is zero for all non-constrained indices.
+ * Since \f$LI = L\f$ it is clear that the system
+ *
+ * \f[ \hat{u} \in \mathbb{R}^n: \qquad \langle \hat{A}, v \rangle = \langle \hat{b}, v \rangle \qquad \forall v \in \mathbb{R}^n \f]
+ *
+ * with \f$ \hat{A} = L^T A L + \hat{I} \f$ and \f$\hat{b} = L^T(b-A c) + c\f$
+ * is equivalent in the sense that \f$u_0 = L \hat{u}\f$. Hence we can solve
+ * the constrained problem by solving
+ *
+ * \f[ \hat{A}\hat{u} = \hat{b} \f]
+ *
+ * and then setting \f$u = \Psi(\hat{u}) = L\hat{u} + c = L u_0 + c\f$.
+ *
+ * The present class supports this approach as follows:
+ * The pair \f$(A,b)\f$ can be transformed to \f$(\hat{A}, \hat{b})\f$
+ * using the method `constrainLinearSystem(A,b)`. Since this is an inplace operation
+ * `extendMatrixPattern(...)` can be used before, to ensure that the pattern of the
+ * matrix is large enough to also store \f$\hat{A}\f$. Once the system is solved,
+ * \f$\hat{u}\f$ can be transformed to \f$u=\Psi(\hat{u})\f$ using
+ * `interpolate(u)` which is again an inplace operation.
*
* \warning This class is still experimental and the interface may change.
*/
@@ -132,8 +173,8 @@ public:
/**
* \brief Const access to linear part of affine map
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this result represents the sparse matrix \f$C\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * the result represents the sparse matrix \f$L\f$.
*/
const LinearTerm& linearTerm() const
{
@@ -143,8 +184,8 @@ public:
/**
* \brief Mutable access to linear part of affine map
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this result represents the sparse matrix \f$C\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * the result represents the sparse matrix \f$L\f$.
*/
LinearTerm& linearTerm()
{
@@ -154,8 +195,8 @@ public:
/**
* \brief Const access to i-th row of linear part of affine map
*
- * Assuming that the affine constraints is written as \f$ x = Cx + b\f$
- * this result represents the i-th row of the sparse matrix \f$C\f$.
+ * Assuming that the affine constraints is written as \f$ x = L x + c\f$
+ * the result represents the i-th row of the sparse matrix \f$L\f$.
*
* This method should only be used if the i-th DOF is constrained.
* For non-constrained DOFs this returns an empty range instead
@@ -174,8 +215,8 @@ public:
/**
* \brief Const access to constant part of affine map
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this result represents the vector \f$b\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * the result represents the vector \f$c\f$.
*/
const ConstantTerm& constantTerm() const
{
@@ -185,8 +226,8 @@ public:
/**
* \brief Mutable access to constant part of affine map
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this result represents the vector \f$b\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * the result represents the vector \f$c\f$.
*/
ConstantTerm& constantTerm()
{
@@ -196,8 +237,8 @@ public:
/**
* \brief Const access to i-th row of constant part of affine map
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this result represents the i-th row of the vector \f$b\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * the result represents the i-th row of the vector \f$c\f$.
*/
template<class MultiIndexType>
const auto& constantTerm(const MultiIndexType& i) const
@@ -208,8 +249,8 @@ public:
/**
* \brief Mutable access to i-th row of constant part of affine map
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this result represents the i-th row of the vector \f$b\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * the result represents the i-th row of the vector \f$c\f$.
*/
template<class MultiIndexType>
auto& constantTerm(const MultiIndexType& i)
@@ -220,8 +261,8 @@ public:
/**
* \brief Interpolate vector according to affine constraints
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this overwrites the given vector \f$v\f$ by \f$Cv+b\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * this overwrites the given vector \f$v\f$ by \f$L v+c\f$.
*
* This computed the values of constrained DOFs using the
* interpolation weights while non-constrained DOFs remain
@@ -237,8 +278,8 @@ public:
if (isConstrained(i))
{
vector_i = constantTerm(i);
- for(const auto& [j, c_ij] : linearTerm(i))
- vector_i += c_ij * vectorBackend[j];
+ for(const auto& [j, L_ij] : linearTerm(i))
+ vector_i += L_ij * vectorBackend[j];
}
});
}
@@ -246,12 +287,12 @@ public:
/**
* \brief Extend matrix pattern for a constrained basis
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
* and the given pattern is able to store a sparse matrix \f$A\f$
- * this enlarges the pattern for storing \f$C^TAC\f$.
+ * this enlarges the pattern for storing \f$L^T A L\f$.
*/
template<class MatrixPatternBackend, class Basis>
- void constrainMatrixPattern(MatrixPatternBackend& patternBackend, const Basis& basis) const
+ void extendMatrixPattern(MatrixPatternBackend& patternBackend, const Basis& basis) const
{
auto localView = basis.localView();
for(const auto& element : elements(basis.gridView()))
@@ -266,16 +307,16 @@ public:
if (isConstrained(i))
{
if (isConstrained(j))
- for(const auto& [k, c_ik] : linearTerm(i))
- for(const auto& [l, c_jl] : linearTerm(j))
+ for(const auto& [k, L_ik] : linearTerm(i))
+ for(const auto& [l, L_jl] : linearTerm(j))
patternBackend.insertEntry(k, l);
else
- for(const auto& [k, c_ik] : linearTerm(i))
+ for(const auto& [k, L_ik] : linearTerm(i))
patternBackend.insertEntry(k, j);
}
else if (isConstrained(j))
{
- for(const auto& [l, c_jl] : linearTerm(j))
+ for(const auto& [l, L_jl] : linearTerm(j))
patternBackend.insertEntry(i, l);
}
}
@@ -286,9 +327,9 @@ public:
/**
* \brief Constrain matrix to a constrained basis
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this overwrites the given matrix \f$A\f$ by \f$C^TAC+\hat{E}\f$ where
- * \f$E\f$ is the identity matrix constrained to the kernel of \f$C\f$
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * this overwrites the given matrix \f$A\f$ by \f$L^T A L+\hat{I}\f$ where
+ * \f$\hat{I}\f$ is the identity matrix constrained to the kernel of \f$L\f$
* given by the constrained rows.
*/
template<class Matrix>
@@ -306,18 +347,18 @@ public:
if (isConstrained(i))
{
if (isConstrained(j))
- for(const auto& [k, c_ik] : linearTerm(i))
- for(const auto& [l, c_jl] : linearTerm(j))
- matrixBackend(k, l) += c_ik*c_jl * matrix_ij;
+ for(const auto& [k, L_ik] : linearTerm(i))
+ for(const auto& [l, L_jl] : linearTerm(j))
+ matrixBackend(k, l) += L_ik*L_jl * matrix_ij;
else
- for(const auto& [k, c_ik] : linearTerm(i))
- matrixBackend(k, j) += c_ik * matrix_ij;
+ for(const auto& [k, L_ik] : linearTerm(i))
+ matrixBackend(k, j) += L_ik * matrix_ij;
matrix_ij = equals(i, j) ? 1 : 0;
}
else if (isConstrained(j))
{
- for(const auto& [l, c_jl] : linearTerm(j))
- matrixBackend(i, l) += c_jl * matrix_ij;
+ for(const auto& [l, L_jl] : linearTerm(j))
+ matrixBackend(i, l) += L_jl * matrix_ij;
matrix_ij = equals(i, j) ? 1 : 0;
}
});
@@ -326,19 +367,19 @@ public:
/**
* \brief Constrain vector to a constrained basis
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this overwrites the given vector \f$r\f$ by \f$C^Tr + b\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * this overwrites the given vector \f$b\f$ by \f$L^T b + c\f$.
*/
template<class Vector>
- void constrainVector(Vector& r) const
+ void constrainVector(Vector& b) const
{
- auto vectorBackend = Dune::Functions::istlVectorBackend(r);
- Dune::Fufem::recursiveForEachVectorEntry(r, [&](auto&& vector_i, const auto& i)
+ auto vectorBackend = Dune::Functions::istlVectorBackend(b);
+ Dune::Fufem::recursiveForEachVectorEntry(b, [&](auto&& vector_i, const auto& i)
{
if (isConstrained(i))
{
- for(const auto& [k, c_ik] : linearTerm(i))
- vectorBackend[k] += c_ik * vector_i;
+ for(const auto& [k, L_ik] : linearTerm(i))
+ vectorBackend[k] += L_ik * vector_i;
vector_i = constantTerm(i);
}
});
@@ -347,16 +388,17 @@ public:
/**
* \brief Constrain linear system affine subspace satisfying the constraints
*
- * Assuming that the affine constraints are written as \f$ x = Cx + b\f$
- * this overwrites the given vector \f$r\f$ by \f$C^T(r-Ab) + b\f$
- * and the matrix \f$A\f$ by \f$C^TAC\f$.
+ * Assuming that the affine constraints are written as \f$ x = L x + c\f$
+ * this overwrites the given vector \f$b\f$ by \f$L^T(b-A c) + c\f$
+ * and the matrix \f$A\f$ by \f$L^T A L + \hat{I}\f$ where
+ * \f$\hat{I}\f$ is the identity matrix constrained to the kernel of \f$L\f$.
*/
template<class Matrix, class Vector>
- void constrainLinearSystem(Matrix& A, Vector& r) const
+ void constrainLinearSystem(Matrix& A, Vector& b) const
{
- A.mmv(constantTerm_, r);
+ A.mmv(constantTerm_, b);
constrainMatrix(A);
- constrainVector(r);
+ constrainVector(b);
}
/**
@@ -364,11 +406,11 @@ public:
*/
void report() const
{
- for(const auto& [i, c_i] : interpolation_)
+ for(const auto& [i, L_i] : interpolation_)
{
std::cout << "x_" << i << " = ";
- for(const auto& [j, c_ij] : c_i)
- std::cout << c_ij << " * x_" << j << " + " ;
+ for(const auto& [j, L_ij] : L_i)
+ std::cout << L_ij << " * x_" << j << " + " ;
std::cout << constantTerm(i) << std::endl;
}
std::cout << "interpolated dofs : " << interpolation_.size() << std::endl;
@@ -382,9 +424,9 @@ public:
bool check() const
{
bool pass = true;
- for(const auto& [i, c_i] : interpolation_)
+ for(const auto& [i, L_i] : interpolation_)
{
- for(const auto& [j, c_ij] : c_i)
+ for(const auto& [j, L_ij] : L_i)
{
if (isConstrained(j))
{
@@ -412,39 +454,39 @@ public:
auto dependenciesResolved = isConstrained_;
auto dependenciesResolvedBackend = Dune::Functions::istlVectorBackend(dependenciesResolved);
dependenciesResolvedBackend = false;
- for(auto& [i, c_i] : interpolation_)
+ for(auto& [i, L_i] : interpolation_)
if (not(dependenciesResolvedBackend[i]))
- resolveDOF(dependenciesResolved, i, c_i);
+ resolveDOF(dependenciesResolved, i, L_i);
}
private:
- void resolveDOF(BitVector& dependenciesResolved, const MultiIndex& i, LinearTermRow& c_i)
+ void resolveDOF(BitVector& dependenciesResolved, const MultiIndex& i, LinearTermRow& L_i)
{
auto dependenciesResolvedBackend = Dune::Functions::istlVectorBackend(dependenciesResolved);
- auto b = Dune::Functions::istlVectorBackend(constantTerm_);
- auto it = c_i.begin();
- while (it != c_i.end())
+ auto c = Dune::Functions::istlVectorBackend(constantTerm_);
+ auto it = L_i.begin();
+ while (it != L_i.end())
{
- auto& [j, c_ij] = *it;
+ auto& [j, L_ij] = *it;
// If interpolation weight is itself constraint,
// first resolve it recursively (if needed)
// and then insert all indirect dependencies directly
if (isConstrained(j))
{
- b[i] += c_ij * b[j];
- auto c_j_it = interpolation_.find(j);
- if (c_j_it != interpolation_.end())
+ c[i] += L_ij * c[j];
+ auto L_j_it = interpolation_.find(j);
+ if (L_j_it != interpolation_.end())
{
- auto& c_j = c_j_it->second;
+ auto& L_j = L_j_it->second;
if (not(dependenciesResolvedBackend[j]))
- resolveDOF(dependenciesResolved, j, c_j);
- for (const auto& [k, c_jk] : c_j)
- c_i[k] += c_ij * c_jk;
+ resolveDOF(dependenciesResolved, j, L_j);
+ for (const auto& [k, L_jk] : L_j)
+ L_i[k] += L_ij * L_jk;
}
else
dependenciesResolvedBackend[j] = true;
- it = c_i.erase(it);
+ it = L_i.erase(it);
}
else
++it;
diff --git a/examples/utilities.hh b/examples/utilities.hh
index e4340143..6f94b854 100644
--- a/examples/utilities.hh
+++ b/examples/utilities.hh
@@ -75,7 +75,7 @@ public:
operatorAssembler.assembleBulkPattern(patternBackend, Dune::resolveRef(args)...);
- constraints.constrainMatrixPattern(patternBackend, testBasis_);
+ constraints.extendMatrixPattern(patternBackend, testBasis_);
patternBackend.setupMatrix();
@@ -126,7 +126,7 @@ public:
operatorAssembler.assembleBulkPattern(patternBackend, Dune::resolveRef(args)...);
- constraints.constrainMatrixPattern(patternBackend, testBasis_);
+ constraints.extendMatrixPattern(patternBackend, testBasis_);
patternBackend.setupMatrix();
--
GitLab