From d223077dcac606c4837855ab880a34322576a312 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sat, 24 Oct 2015 12:24:54 -0500 Subject: [PATCH] Refactor into constraint_linear_operator header file This commit refactors everything into its own header file and addresses a number of issues: - solve (Ct * A C + Id_c) x = Ct (b - Ak) - make code compile and refactor implementation - rewrite documentation --- .../deal.II/lac/constraint_linear_operator.h | 291 ++++++++++++++++++ include/deal.II/lac/linear_operator.h | 157 ---------- include/deal.II/lac/packaged_operation.h | 71 ----- 3 files changed, 291 insertions(+), 228 deletions(-) create mode 100644 include/deal.II/lac/constraint_linear_operator.h diff --git a/include/deal.II/lac/constraint_linear_operator.h b/include/deal.II/lac/constraint_linear_operator.h new file mode 100644 index 0000000000..347d845da1 --- /dev/null +++ b/include/deal.II/lac/constraint_linear_operator.h @@ -0,0 +1,291 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__constraint_linear_operator_h +#define dealii__constraint_linear_operator_h + +#include +#include +#include + +#ifdef DEAL_II_WITH_CXX11 + +DEAL_II_NAMESPACE_OPEN + + +/** + * @name Indirectly applying constraints to LinearOperator + */ +//@{ + + +/** + * This function takes a ConstraintMatrix @p cm and an operator exemplar @p + * exemplar (this exemplar is usually a linear operator that describes the + * system matrix) and returns a LinearOperator object associated with the + * "homogeneous action" of the underlying ConstraintMatrix object: + * + * Applying the LinearOperator object on a vector u results in + * a vector v that stores the result of calling + * ConstraintMatrix::distribute() on u - with one important + * difference; inhomogeneities are note applied, but always treated as 0 + * instead. + * + * The LinearOperator object created by this function is primarily used + * internally in constrained_linear_operator() to build up a modified + * system of linear equations. How to solve a linear system of equations + * with this approach is explained in detail in the @ref constraints + * module. + * + * @author Mauro Bardelloni, Matthias Maier, 2015 + * + * @relates LinearOperator + * @ingroup constraints + */ +template +LinearOperator distribute_constraints_linear_operator( + const ConstraintMatrix &constraint_matrix, + const LinearOperator &exemplar) +{ + LinearOperator return_op = exemplar; + + return_op.vmult_add = [&constraint_matrix](Range &v, const Domain &u) + { + for (auto i : v.locally_owned_elements()) + { + if (constraint_matrix.is_constrained(i)) + { + const auto *entries = constraint_matrix.get_constraint_entries (i); + for (types::global_dof_index j = 0; j < entries->size(); ++j) + { + const auto pos = (*entries)[j].first; + v(i) += u(pos) * (*entries)[j].second; + } + } + else + v(i) += u(i); + } + }; + + return_op.Tvmult_add = [&constraint_matrix](Domain &v, const Range &u) + { + for (auto i : v.locally_owned_elements()) + { + if (constraint_matrix.is_constrained(i)) + { + const auto *entries = constraint_matrix.get_constraint_entries(i); + for (types::global_dof_index j = 0; j < entries->size(); ++j) + { + const auto pos = (*entries)[j].first; + v(pos) += u(i) * (*entries)[j].second; + } + } + else + v(i)+=u(i); + } + }; + + // lambda capture expressions are a C++14 feature... + const auto vmult_add = return_op.vmult_add; + return_op.vmult = [vmult_add](Range &v, const Domain &u) + { + v = 0.; + vmult_add(v, u); + }; + + // lambda capture expressions are a C++14 feature... + const auto Tvmult_add = return_op.Tvmult_add; + return_op.Tvmult = [Tvmult_add](Domain &v, const Range &u) + { + v = 0.; + Tvmult_add(v, u); + }; + + return return_op; +} + + +/** + * Given a ConstraintMatrix @p constraint_matrix and an operator exemplar + * @p exemplar, return an LinearOperator that is the projection to the + * subspace of constrained degrees of freedom. + * + * @author Mauro Bardelloni, Matthias Maier, 2015 + * + * @relates LinearOperator + * @ingroup constraints + */ +template +LinearOperator project_to_constrained_linear_operator( + const ConstraintMatrix &constraint_matrix, + const LinearOperator &exemplar) +{ + LinearOperator return_op = exemplar; + + return_op.vmult_add = [&constraint_matrix](Range &v, const Domain &u) + { + for (auto i : v.locally_owned_elements()) + if (constraint_matrix.is_constrained(i)) + v(i) += u(i); + }; + + return_op.Tvmult_add = [&constraint_matrix](Domain &v, const Range &u) + { + for (auto i : v.locally_owned_elements()) + if (constraint_matrix.is_constrained(i)) + v(i) += u(i); + }; + + return_op.vmult = [&constraint_matrix](Range &v, const Domain &u) + { + v = 0.; + for (auto i : v.locally_owned_elements()) + if (constraint_matrix.is_constrained(i)) + v(i) += u(i); + }; + + return_op.Tvmult = [&constraint_matrix](Domain &v, const Range &u) + { + v = 0.; + for (auto i : v.locally_owned_elements()) + if (constraint_matrix.is_constrained(i)) + v(i) += u(i); + }; + + return return_op; +} + + +/** + * Given a ConstraintMatrix object @p constraint_matrix and a + * LinearOperator @p linop, this function creates a LinearOperator object + * consisting of the composition of three operations and a regularization: + * @code + * Ct * linop * C + Id_c; + * @endcode + * with + * @code + * C = distribute_constraints_linear_operator(constraint_matrix, linop); + * Ct = transpose_operator(C); + * Id_c = project_to_constrained_linear_operator(constraint_matrix, linop); + * @endcode + * and Id_c is the projection to the subspace consisting of + * all vector entries associated with constrained degrees of freedoms. + * + * This LinearOperator object is used together with + * constrained_right_hand_side() to build up the following modified system + * of linear equations: + * @f[ + * (C^T \cdot A \cdot C + Id_c) x = C^T (b - A\,k) + * @f] + * with a given (unconstrained) system matrix $A$, right hand side $b$, and + * linear constraints $C$ with inhomogeneities $k$. + * + * A detailed explanation of this approach is given in the @ref constraints + * module. + * + * @author Mauro Bardelloni, Matthias Maier, 2015 + * + * @relates LinearOperator + * @ingroup constraints + */ +template +LinearOperator +constrained_linear_operator(const ConstraintMatrix &constraint_matrix, + const LinearOperator &linop) +{ + const auto C = + distribute_constraints_linear_operator(constraint_matrix, linop); + const auto Ct = transpose_operator(C); + const auto Id_c = + project_to_constrained_linear_operator(constraint_matrix, linop); + return Ct * linop * C + Id_c; +} + + +/** + * Given a ConstraintMatrix object @p constraint_matrix, a LinearOperator + * @p linop and a right-hand side @p right_hand_side, this function creates + * a PackagedOperation that stores the following computation: + * @code + * Ct * (right_hand_side - linop * k) + * @endcode + * with + * @code + * C = distribute_constraints_linear_operator(constraint_matrix, linop); + * Ct = transpose_operator(C); + * @endcode + * + * This LinearOperator object is used together with + * constrained_right_hand_side() to build up the following modified system + * of linear equations: + * @f[ + * (C^T \cdot A \cdot C + Id_c) x = C^T (b - A\,k) + * @f] + * with a given (unconstrained) system matrix $A$, right hand side $b$, and + * linear constraints $C$ with inhomogeneities $k$. + * + * A detailed explanation of this approach is given in the @ref constraints + * module. + * + * @author Mauro Bardelloni, Matthias Maier, 2015 + * + * @relates LinearOperator + * @ingroup constraints + */ +template +PackagedOperation +constrained_right_hand_side(const ConstraintMatrix &constraint_matrix, + const LinearOperator &linop, + const Range &right_hand_side) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = linop.reinit_range_vector; + + return_comp.apply_add = + [&constraint_matrix, &linop, &right_hand_side](Range &v) + { + const auto C = + distribute_constraints_linear_operator(constraint_matrix, linop); + const auto Ct = transpose_operator(C); + + static GrowingVectorMemory vector_memory; + Domain *k = vector_memory.alloc(); + linop.reinit_domain_vector(*k, /*bool fast=*/ false); + constraint_matrix.distribute(*k); + + v += Ct * (right_hand_side - linop **k); + + vector_memory.free(k); + }; + + // lambda capture expressions are a C++14 feature... + const auto apply_add = return_comp.apply_add; + return_comp.apply = [apply_add](Range &v) + { + v = 0.; + apply_add(v); + }; + + return return_comp; +} + +//@} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_CXX11 +#endif diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index d63f733b81..eccf331ab1 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -1116,163 +1116,6 @@ linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix) return return_op; } -template , - typename Domain = Range, - typename Matrix> -LinearOperator -constraints_linear_operator(const ConstraintMatrix &, const Matrix &); - -/** - * @relates LinearOperator - * - * constraints_linear_operator takes a constraint matrix @p cm and an exemplar - * matrix @p m (this matrix could be the system matrix). - * The result of this function is the LinearOperator associated to the constraint - * matrix. - * In detail, a constrained problem can be expressed in this form: - * Ct * S * C * x = Ct * ( rhs + S * k) - * Where: - * - C has the role of constraint matrix - * - k is the vector representing constraints inhomogeneities. - * (a generic contraint can be expressed as x = Cx + k) - * - S is the system matrix - * - rhs is the original right-hand-side - * This function returns a LinearOperator representing the matrix C. - * - * @note Suppose we have n dof and m constraints. W.l.o.g. we can assume that it - * is possible to express n-m variables in terms of remainder variables - * (constrained variables). - * Therefore, $ x_i = C_{i,j} x_j + k_i $ for $j = 1, ..., n-m$ - * and $i = 1, ..., n$. - * Notice that Ct * S * C is a problem in ${\mathbb R}^{m-n}$, remainder - * variables are treated solving x = 0 in order to have a well-posed problem - * on ${\mathbb R}^n$. - * At the end a solution of the problem holding for constrained variables - * can be found applying the constraint equations. (This is equivalent to - * cm.distribute(x)). - * - * @see M. S. Shephard: Linear multipoint constraints applied via - * transformation as part of a direct stiffness assembly process, 1985. - * For more details. - * - * @relates ConstraintMatrix - * @ingroup LAOperators - */ -template < typename Range, - typename Domain, - typename Matrix > -LinearOperator -constraints_linear_operator(const ConstraintMatrix &cm, const Matrix &m) -{ - LinearOperator return_op = linear_operator(m); - - return_op.vmult_add = [&cm](Range &v, const Domain &u) - { - for (auto i : v.locally_owned_elements()) - { - if (cm.is_constrained(i)) - { - const std::vector< std::pair < types::global_dof_index, double > > - *entries = cm.get_constraint_entries (i); - for (types::global_dof_index j=0; j < entries->size(); ++j) - { - types::global_dof_index pos = (*entries)[j].first; - v(i) += u(pos) * (*entries)[j].second; - } - } - else - v(i) += u(i); - } - }; - - return_op.Tvmult_add = [&cm](Range &v, const Domain &u) - { - for (auto i : u.locally_owned_elements()) - { - if (cm.is_constrained(i)) - { - const std::vector< std::pair < types::global_dof_index, double > > - *entries = cm.get_constraint_entries (i); - for (types::global_dof_index j=0; j < entries->size(); ++j) - { - types::global_dof_index pos = (*entries)[j].first; - v(pos) += u(i) * (*entries)[j].second; - } - } - else - v(i)+=u(i); - - } - }; - - return_op.vmult = [&cm](Range &v, const Domain &u) - { - v = 0.; - vmult_add(v, u); - }; - - return_op.Tvmult = [&cm](Range &v, const Domain &u) - { - v = 0.; - Tvmult_add(v, u); - }; - - return return_op; -} - -template , - typename Domain = Range, - typename Matrix> -LinearOperator -constrained_linear_operator(const ConstraintMatrix &, const Matrix &); - -/** - * @relates LinearOperator - * - * Let @p m be a system matrix and assume we are interested in this problem - * constrained with a constraint matrix @p cm. - * This function returns the LinearOperator that replace the system matrix and - * is associated to this problem. - * In detail, a constrained problem can be expressed in this form: - * Ct * S * C * x = Ct * ( rhs + S * k) - * Where: - * - C has the role of constraint matrix - * - k is the vector representing constraints inhomogeneities. - * (a generic constraint can be expressed as x = Cx + k) - * - S is the system matrix - * - rhs is the original right-hand-side - * This function returns a LinearOperator representing the matrix Ct * S * C. - * - * @note Suppose to express the form of $C$, let us consider that it corresponds - * to a ConstraintMatrix object on $n$ degrees of freedom, of which $m\le n$ are - * constrained. Without loss of generality, let us assume that the first $m$ are - * the ones that are constrained, and have the form $x_i = \sum_{j=m+1}^n c_{ij} - * x_j + b_j$ for $i=1,\ldots, m$. Then the rest of the degrees of freedom could - * be written as $x_i=x_i$ for $i=m+1,\ldots, n$. The matrix $C$ is then that - * matrix that has entries $C_{ij}=c_{ij}$ for $1\le i \le m, m+1\le j\le n$, - * $c_{ii}=1$ for $m+1\le i \le n$, and $C_{ij}=0$ for all other elements. - * - * @see M. S. Shephard: Linear multipoint constraints applied via - * transformation as part of a direct stiffness assembly process, 1985. - * For more details. - * - * @relates ConstraintMatrix - * @ingroup LAOperators - */ -template < typename Range, - typename Domain, - typename Matrix > -LinearOperator -constrained_linear_operator(const ConstraintMatrix &cm, const Matrix &m) -{ - auto C = constraints_linear_operator(cm, m); - auto Ct = transpose_operator(C); - auto A = linear_operator(m); - - return Ct * A * C; -} - - //@} DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/lac/packaged_operation.h b/include/deal.II/lac/packaged_operation.h index 6b097cbe53..76d3144ffe 100644 --- a/include/deal.II/lac/packaged_operation.h +++ b/include/deal.II/lac/packaged_operation.h @@ -827,77 +827,6 @@ operator*(const PackagedOperation &comp, return return_comp; } -template , - typename Domain = Range, - typename Matrix> -PackagedOperation -constrained_rhs(const ConstraintMatrix &, const Matrix &, const Range &); - -/** - * @relates LinearOperator - * - * Let @p m be a system matrix, @p rhs a vector, and @p cm a constraint matrix. - * - * constrained_rhs returns the rhs associated to @p rhs of the constrained system - * induced by @p m and @p cm. - * In detail, a constrained problem can be expressed in this form: - * Ct * S * C * x = Ct * ( rhs + S * k) - * Where: - * - C has the role of constraint matrix - * - k is the vector representing constraints inhomogeneities. - * (a generic contraint can be expressed as x = Cx + k) - * - S is the system matrix - * - rhs is the original right-hand-side - * This function returns a LinearOperator representing the vector - * Ct * ( rhs + S * k). - * - * @note Suppose we have n dof and m constraints. W.l.o.g. we can assume that it - * is possible to express n-m variables in terms of remainder variables - * (constrained variables). - * Therefore, $ x_i = C_{i,j} x_j + k_i $ for $j = 1, ..., n-m$ - * and $i = 1, ..., n$. - * Notice that Ct * S * C is a problem in ${\mathbb R}^{m-n}$, remainder - * variables are treated solving x = 0 in order to have a well-posed problem - * on ${\mathbb R}^n$. - * At the end a solution of the problem holding for constrained variables - * can be found applying the constraint equations. (This is equivalent to - * cm.distribute(x)). - * - * @see M. S. Shephard: Linear multipoint constraints applied via - * transformation as part of a direct stiffness assembly process, 1985. - * For more details. - * - * @relates ConstraintMatrix - * @ingroup LAOperators - */ -template -PackagedOperation -constrained_rhs(const ConstraintMatrix &cm, const Matrix &m, const Range &rhs) -{ - PackagedOperation return_comp; - - return_comp.reinit_vector = [&m](Range &v, bool fast) - { - internal::LinearOperator::ReinitHelper::reinit_range_vector(m, v, fast); - }; - - return_comp.apply = [&cm, &m, &rhs](Range &v) - { - const auto M = linear_operator(m); - const auto C = constraints_linear_operator(cm, m); - const auto Ct = transpose_operator(C); - - for (auto i : v.locally_owned_elements()) - v[i] = -1 * cm.get_inhomogeneity(i); - - v = Ct * ( rhs + M * v); - }; - - return return_comp; -} - //@} DEAL_II_NAMESPACE_CLOSE -- 2.39.5