]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactor into constraint_linear_operator header file
authorMatthias Maier <tamiko@43-1.org>
Sat, 24 Oct 2015 17:24:54 +0000 (12:24 -0500)
committerMatthias Maier <tamiko@43-1.org>
Mon, 28 Dec 2015 11:36:24 +0000 (12:36 +0100)
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

include/deal.II/lac/constraint_linear_operator.h [new file with mode: 0644]
include/deal.II/lac/linear_operator.h
include/deal.II/lac/packaged_operation.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 (file)
index 0000000..347d845
--- /dev/null
@@ -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 <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#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 <code>u</code> results in
+ * a vector <code>v</code> that stores the result of calling
+ * ConstraintMatrix::distribute() on <code>u</code> - 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 <typename Range, typename Domain>
+LinearOperator<Range, Domain> distribute_constraints_linear_operator(
+  const ConstraintMatrix &constraint_matrix,
+  const LinearOperator<Range, Domain> &exemplar)
+{
+  LinearOperator<Range, Domain> 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 <typename Range, typename Domain>
+LinearOperator<Range, Domain> project_to_constrained_linear_operator(
+  const ConstraintMatrix &constraint_matrix,
+  const LinearOperator<Range, Domain> &exemplar)
+{
+  LinearOperator<Range, Domain> 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 <code>Id_c</code> 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 <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+constrained_linear_operator(const ConstraintMatrix &constraint_matrix,
+                            const LinearOperator<Range, Domain> &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 <typename Range, typename Domain>
+PackagedOperation<Range>
+constrained_right_hand_side(const ConstraintMatrix &constraint_matrix,
+                            const LinearOperator<Range, Domain> &linop,
+                            const Range &right_hand_side)
+{
+  PackagedOperation<Range> 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<Domain> 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
index d63f733b81a5b0c600a65080575b04ff1bf9da33..eccf331ab16fe8ad8588717cbbe97ea41f94b00d 100644 (file)
@@ -1116,163 +1116,6 @@ linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
   return return_op;
 }
 
-template <typename Range = Vector<double>,
-          typename Domain = Range,
-          typename Matrix>
-LinearOperator<Range, Domain>
-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<Range, Domain>
-constraints_linear_operator(const ConstraintMatrix &cm, const Matrix &m)
-{
-  LinearOperator<Range, Domain> return_op = linear_operator<Range, Domain>(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 Range = Vector<double>,
-          typename Domain = Range,
-          typename Matrix>
-LinearOperator<Range, Domain>
-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<Range, Domain>
-constrained_linear_operator(const ConstraintMatrix &cm, const Matrix &m)
-{
-  auto C  = constraints_linear_operator<Range, Domain, Matrix>(cm, m);
-  auto Ct = transpose_operator<Domain, Range>(C);
-  auto A  = linear_operator<Range, Domain>(m);
-
-  return Ct * A * C;
-}
-
-
 //@}
 
 DEAL_II_NAMESPACE_CLOSE
index 6b097cbe53ac4e5cc9b42f741c7b3f7e312bab17..76d3144ffe1284a091af55f47ce39db4082a75c8 100644 (file)
@@ -827,77 +827,6 @@ operator*(const PackagedOperation<Range> &comp,
   return return_comp;
 }
 
-template <typename Range = Vector<double>,
-          typename Domain = Range,
-          typename Matrix>
-PackagedOperation<Range>
-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 <typename Range,
-          typename Domain,
-          typename Matrix>
-PackagedOperation<Range>
-constrained_rhs(const ConstraintMatrix &cm, const Matrix &m, const Range &rhs)
-{
-  PackagedOperation<Range> return_comp;
-
-  return_comp.reinit_vector = [&m](Range &v, bool fast)
-  {
-    internal::LinearOperator::ReinitHelper<Range>::reinit_range_vector(m, v, fast);
-  };
-
-  return_comp.apply = [&cm, &m, &rhs](Range &v)
-  {
-    const auto M = linear_operator<Range, Domain>(m);
-    const auto C  = constraints_linear_operator<Range, Domain, Matrix>(cm, m);
-    const auto Ct = transpose_operator<Domain, Range>(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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.