]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Restructure documentation
authorMatthias Maier <tamiko@43-1.org>
Mon, 27 Apr 2015 10:02:29 +0000 (12:02 +0200)
committerMatthias Maier <tamiko@43-1.org>
Mon, 11 May 2015 21:05:20 +0000 (23:05 +0200)
include/deal.II/lac/computation.h
include/deal.II/lac/linear_operator.h

index bcadf0fd75f219e3747ac3a4b88be4e4d02fc17a..a583ef267da6e642e124ecaee5bff043ad581397 100644 (file)
@@ -265,7 +265,7 @@ public:
 //@{
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Addition of two Computation objects @p first_comp and @p second_comp given by
  * vector space addition of the corresponding results.
@@ -300,7 +300,7 @@ operator+(const Computation<Range> &first_comp,
 }
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Subtraction of two Computation objects @p first_comp and @p second_comp
  * given by vector space addition of the corresponding results.
@@ -338,7 +338,7 @@ operator-(const Computation<Range> &first_comp,
 }
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Scalar multiplication of a Computation objects @p comp with a scalar @p number
  * given by a scaling Computation result with @p number.
@@ -371,7 +371,7 @@ operator*(const Computation<Range> &comp,
 }
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Scalar multiplication of a Computation objects @p comp with a scalar @p number
  * given by a scaling Computation result with @p number.
@@ -387,7 +387,7 @@ operator*(typename Range::value_type number,
 }
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Add a constant @p offset (of the @p Range space) to the result of a
  * Computation.
@@ -402,7 +402,7 @@ Computation<Range> operator+(const Computation<Range> &comp,
 }
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Add a constant @p offset (of the @p Range space) to the result of a
  * Computation.
@@ -417,7 +417,7 @@ Computation<Range> operator+(const Range &offset,
 }
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Subtract a constant @p offset (of the @p Range space) from the result of a
  * Computation.
@@ -433,7 +433,7 @@ Computation<Range> operator-(const Computation<Range> &comp,
 
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Subtract a computational result from a constant @p offset (of the @p
  * Range space). The result is a Computation object that applies this
@@ -570,7 +570,7 @@ Computation<Range> operator-(const Range &u, const Range &v)
 
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Create a Computation object from a LinearOperator and a reference to a
  * vector @p u of the Domain space. The object stores the computation
@@ -612,7 +612,7 @@ operator*(const LinearOperator<Range, Domain> &op,
 
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Create a Computation object from a LinearOperator and a reference to a
  * vector @p u of the Range space. The object stores the computation
@@ -654,7 +654,7 @@ operator*(const Range &u,
 
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Composition of a Computation object with a LinearOperator.
  * The object stores the computation $\text{op} \,comp$ (in matrix notation).
@@ -704,7 +704,7 @@ operator*(const LinearOperator<Range, Domain> &op,
 
 
 /**
- * \relates Computation
+ * @relates Computation
  *
  * Composition of a Computation object with a LinearOperator.
  * The object stores the computation $\text{op}^T \,comp$ (in matrix notation).
index c238222d4239df63e9eec735cb6b5e4a922b7d92..71915eb777ec215e7467e17e36ce16ca9359fd59 100644 (file)
@@ -225,6 +225,10 @@ public:
    */
   std::function<void(Domain &v, bool fast)> reinit_domain_vector;
 
+  /**
+   * @name In-place vector space operations
+   */
+  //@{
 
   /**
    * Addition with a LinearOperator @p second_op with the same @p Domain
@@ -270,12 +274,17 @@ public:
     return *this;
   }
 
+  //@}
 };
 
 
+/**
+ * @name Vector space operations
+ */
+//@{
 
 /**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
  * Addition of two linear operators @p first_op and @p second_op given by
  * $(\text{first\_op}+\text{second\_op})x:=\text{first\_op}(x)+\text{second\_op}(x)$
@@ -324,7 +333,7 @@ operator+(const LinearOperator<Range, Domain> &first_op,
 
 
 /**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
  * Subtraction of two linear operators @p first_op and @p second_op given by
  * $(\text{first\_op}-\text{second\_op})x:=\text{first\_op}(x)-\text{second\_op}(x)$
@@ -342,76 +351,7 @@ operator-(const LinearOperator<Range, Domain> &first_op,
 
 
 /**
- * \relates LinearOperator
- *
- * Composition of two linear operators @p first_op and @p second_op given by
- * $(\text{first\_op}*\text{second\_op})x:=\text{first\_op}(\text{second\_op}(x))$
- *
- * @ingroup LAOperators
- */
-template <typename Range, typename Intermediate, typename Domain>
-LinearOperator<Range, Domain>
-operator*(const LinearOperator<Range, Intermediate> &first_op,
-          const LinearOperator<Intermediate, Domain> &second_op)
-{
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
-  return_op.reinit_range_vector = first_op.reinit_range_vector;
-
-  // ensure to have valid computation objects by catching first_op and
-  // second_op by value
-
-  return_op.vmult = [first_op, second_op](Range &v, const Domain &u)
-  {
-    static GrowingVectorMemory<Intermediate> vector_memory;
-
-    Intermediate *i = vector_memory.alloc();
-    second_op.reinit_range_vector(*i, /*bool fast =*/ true);
-    second_op.vmult(*i, u);
-    first_op.vmult(v, *i);
-    vector_memory.free(i);
-  };
-
-  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u)
-  {
-    static GrowingVectorMemory<Intermediate> vector_memory;
-
-    Intermediate *i = vector_memory.alloc();
-    second_op.reinit_range_vector(*i, /*bool fast =*/ true);
-    second_op.vmult(*i, u);
-    first_op.vmult_add(v, *i);
-    vector_memory.free(i);
-  };
-
-  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u)
-  {
-    static GrowingVectorMemory<Intermediate> vector_memory;
-
-    Intermediate *i = vector_memory.alloc();
-    first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
-    first_op.Tvmult(*i, u);
-    second_op.Tvmult(v, *i);
-    vector_memory.free(i);
-  };
-
-  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u)
-  {
-    static GrowingVectorMemory<Intermediate> vector_memory;
-
-    Intermediate *i = vector_memory.alloc();
-    first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
-    first_op.Tvmult(*i, u);
-    second_op.Tvmult_add(v, *i);
-    vector_memory.free(i);
-  };
-
-  return return_op;
-}
-
-
-/**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
  * Scalar multiplication of a ScalarOperator object @p op with @p
  * number from the left.
@@ -472,7 +412,7 @@ operator*(typename Range::value_type  number,
 
 
 /**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
  * Scalar multiplication of a ScalarOperator object from the right.
  *
@@ -498,80 +438,108 @@ operator*(const LinearOperator<Range, Domain> &op,
   return number * op;
 }
 
+//@}
+
 
 /**
- * \relates LinearOperator
- *
- * Returns the transpose linear operations of @ref op.
- *
- * @ingroup LAOperators
+ * @name Composition and manipulation of a LinearOperator
  */
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-transpose_operator(const LinearOperator<Range, Domain> &op)
-{
-  LinearOperator<Domain, Range> return_op;
-
-  return_op.reinit_range_vector = op.reinit_domain_vector;
-  return_op.reinit_domain_vector = op.reinit_range_vector;
-
-  return_op.vmult = op.Tvmult;
-  return_op.vmult_add = op.Tvmult_add;
-  return_op.Tvmult = op.vmult;
-  return_op.Tvmult_add = op.vmult_add;
-
-  return return_op;
-}
-
+//@{
 
 /**
- * \relates LinearOperator
- *
- * Returns a LinearOperator that is the identity of the vector space
- * @p Range.
+ * @relates LinearOperator
  *
- * The function takes an <code>std::function</code> object @ref
- * reinit_vector as an argument to initialize the
- * <code>reinit_range_vector</code> and <code>reinit_domain_vector</code>
- * objects of the LinearOperator object.
+ * Composition of two linear operators @p first_op and @p second_op given by
+ * $(\text{first\_op}*\text{second\_op})x:=\text{first\_op}(\text{second\_op}(x))$
  *
  * @ingroup LAOperators
  */
-template <typename Range>
-LinearOperator<Range, Range>
-identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
+template <typename Range, typename Intermediate, typename Domain>
+LinearOperator<Range, Domain>
+operator*(const LinearOperator<Range, Intermediate> &first_op,
+          const LinearOperator<Intermediate, Domain> &second_op)
 {
-  LinearOperator<Range, Range> return_op;
+  LinearOperator<Range, Domain> return_op;
 
-  return_op.reinit_range_vector = reinit_vector;
-  return_op.reinit_domain_vector = reinit_vector;
+  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
+  return_op.reinit_range_vector = first_op.reinit_range_vector;
 
-  return_op.vmult = [](Range &v, const Range &u)
+  // ensure to have valid computation objects by catching first_op and
+  // second_op by value
+
+  return_op.vmult = [first_op, second_op](Range &v, const Domain &u)
   {
-    v = u;
+    static GrowingVectorMemory<Intermediate> vector_memory;
+
+    Intermediate *i = vector_memory.alloc();
+    second_op.reinit_range_vector(*i, /*bool fast =*/ true);
+    second_op.vmult(*i, u);
+    first_op.vmult(v, *i);
+    vector_memory.free(i);
   };
 
-  return_op.vmult_add = [](Range &v, const Range &u)
+  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u)
   {
-    v += u;
+    static GrowingVectorMemory<Intermediate> vector_memory;
+
+    Intermediate *i = vector_memory.alloc();
+    second_op.reinit_range_vector(*i, /*bool fast =*/ true);
+    second_op.vmult(*i, u);
+    first_op.vmult_add(v, *i);
+    vector_memory.free(i);
   };
 
-  return_op.Tvmult = [](Range &v, const Range &u)
+  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u)
   {
-    v = u;
+    static GrowingVectorMemory<Intermediate> vector_memory;
+
+    Intermediate *i = vector_memory.alloc();
+    first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
+    first_op.Tvmult(*i, u);
+    second_op.Tvmult(v, *i);
+    vector_memory.free(i);
   };
 
-  return_op.Tvmult_add = [](Range &v, const Range &u)
+  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u)
   {
-    v += u;
+    static GrowingVectorMemory<Intermediate> vector_memory;
+
+    Intermediate *i = vector_memory.alloc();
+    first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
+    first_op.Tvmult(*i, u);
+    second_op.Tvmult_add(v, *i);
+    vector_memory.free(i);
   };
 
   return return_op;
 }
 
+/**
+ * @relates LinearOperator
+ *
+ * Returns the transpose linear operations of @ref op.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+transpose_operator(const LinearOperator<Range, Domain> &op)
+{
+  LinearOperator<Domain, Range> return_op;
+
+  return_op.reinit_range_vector = op.reinit_domain_vector;
+  return_op.reinit_domain_vector = op.reinit_range_vector;
+
+  return_op.vmult = op.Tvmult;
+  return_op.vmult_add = op.Tvmult_add;
+  return_op.Tvmult = op.vmult;
+  return_op.Tvmult_add = op.vmult_add;
+
+  return return_op;
+}
 
 /**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
  * Returns an object representing the inverse of the LinearOperator @p op.
  *
@@ -639,297 +607,60 @@ inverse_operator(const LinearOperator<typename Solver::vector_type, typename Sol
   return return_op;
 }
 
+//@}
+
+
+/**
+ * @name Creation of a LinearOperator
+ */
+//@{
 
 /**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
- * A function that encapsulates a given collection @p ops of
- * LinearOperators into a block structure. Hereby, it is assumed that Range
- * and Domain are blockvectors, i.e., derived from @ref BlockVectorBase.
- * The individual linear operators in @p ops must act on a the underlying
- * vector type of the block vectors, i.e., on Domain::BlockType yielding a
- * result in Range::BlockType.
+ * Returns a LinearOperator that is the identity of the vector space
+ * @p Range.
  *
- * The list @p ops is best passed as an initializer list. Consider for
- * example a linear operator block (acting on Vector<double>)
- * @code
- *  op_a00 | op_a01
- *         |
- *  ---------------
- *         |
- *  op_a10 | op_a11
- * @endcode
- * The coresponding block_operator invocation takes the form
- * @code
- * block_operator<2, 2, BlockVector<double>>({op_a00, op_a01, op_a10, op_a11});
- * @endcode
+ * The function takes an <code>std::function</code> object @ref
+ * reinit_vector as an argument to initialize the
+ * <code>reinit_range_vector</code> and <code>reinit_domain_vector</code>
+ * objects of the LinearOperator object.
  *
  * @ingroup LAOperators
  */
-template <unsigned int m, unsigned int n,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
-block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &ops)
+template <typename Range>
+LinearOperator<Range, Range>
+identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
 {
-  static_assert(m > 0 && n > 0,
-                "a blocked LinearOperator must consist of at least one block");
+  LinearOperator<Range, Range> return_op;
 
-  LinearOperator<Range, Domain> return_op;
+  return_op.reinit_range_vector = reinit_vector;
+  return_op.reinit_domain_vector = reinit_vector;
 
-  return_op.reinit_range_vector = [ops](Range &v, bool fast)
+  return_op.vmult = [](Range &v, const Range &u)
   {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
-
-    // And reinitialize every individual block with reinit_range_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i][0].reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
+    v = u;
   };
 
-  return_op.reinit_domain_vector = [ops](Domain &v, bool fast)
+  return_op.vmult_add = [](Range &v, const Range &u)
   {
-    // Reinitialize the block vector to n blocks:
-    v.reinit(n);
-
-    // And reinitialize every individual block with reinit_domain_vectors:
-    for (unsigned int i = 0; i < n; ++i)
-      ops[0][i].reinit_domain_vector(v.block(i), fast);
+    v += u;
+  };
 
-    v.collect_sizes();
+  return_op.Tvmult = [](Range &v, const Range &u)
+  {
+    v = u;
   };
 
-  return_op.vmult = [ops](Range &v, const Domain &u)
+  return_op.Tvmult_add = [](Range &v, const Range &u)
   {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
-
-    for (unsigned int i = 0; i < m; ++i)
-      {
-        ops[i][0].vmult(v.block(i), u.block(0));
-        for (unsigned int j = 1; j < n; ++j)
-          ops[i][j].vmult_add(v.block(i), u.block(j));
-      }
-  };
-
-  return_op.vmult_add = [ops](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
-
-    for (unsigned int i = 0; i < m; ++i)
-      for (unsigned int j = 0; j < n; ++j)
-        ops[i][j].vmult_add(v.block(i), u.block(j));
-  };
-
-  return_op.Tvmult = [ops](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < n; ++i)
-      {
-        ops[0][i].Tvmult(v.block(i), u.block(0));
-        for (unsigned int j = 1; j < m; ++j)
-          ops[j][i].Tvmult_add(v.block(i), u.block(j));
-      }
-  };
-
-  return_op.Tvmult_add = [ops](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < n; ++i)
-      for (unsigned int j = 0; j < m; ++j)
-        ops[j][i].Tvmult_add(v.block(i), u.block(j));
-  };
-
-  return return_op;
-}
-
-
-/**
- * \relates LinearOperator
- *
- * A variant of above function that builds up a block diagonal linear
- * operator from an array @p ops of diagonal elements (off-diagonal blocks
- * are assumed to be 0).
- *
- * The list @p ops is best passed as an initializer list. Consider for
- * example a linear operator block (acting on Vector<double>)
- * <code>diag(op_a0, op_a1, ..., op_am)</code>. The coresponding
- * block_operator invocation takes the form
- * @code
- * block_diagonal_operator<m, BlockVector<double>>({op_00, op_a1, ..., op_am});
- * @endcode
- *
- * @ingroup LAOperators
- */
-template <unsigned int m,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
-block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &ops)
-{
-  static_assert(m > 0,
-                "a blockdiagonal LinearOperator must consist of at least one block");
-
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_range_vector = [ops](Range &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
-
-    // And reinitialize every individual block with reinit_range_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.reinit_domain_vector = [ops](Domain &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
-
-    // And reinitialize every individual block with reinit_domain_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.vmult = [ops](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].vmult(v.block(i), u.block(i));
-  };
-
-  return_op.vmult_add = [ops](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].vmult_add(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult = [ops](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].Tvmult(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult_add = [ops](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].Tvmult_add(v.block(i), u.block(i));
-  };
+    v += u;
+  };
 
   return return_op;
 }
 
 
-/**
- * \relates LinearOperator
- *
- * A variant of above function that only takes a single LinearOperator
- * argument @p op and creates a blockdiagonal linear operator with @p m
- * copies of it.
- *
- * @ingroup LAOperators
- */
-template <unsigned int m,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
-block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op)
-{
-
-  static_assert(m > 0,
-                "a blockdiagonal LinearOperator must consist of at least "
-                "one block");
-
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_range_vector = [op](Range &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
-
-    // And reinitialize every individual block with reinit_range_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      op.reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.reinit_domain_vector = [op](Domain &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
-
-    // And reinitialize every individual block with reinit_domain_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      op.reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.vmult = [op](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      op.vmult(v.block(i), u.block(i));
-  };
-
-  return_op.vmult_add = [op](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      op.vmult_add(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult = [op](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      op.Tvmult(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult_add = [op](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      op.Tvmult_add(v.block(i), u.block(i));
-  };
-
-  return return_op;
-}
-
-
-
 namespace internal
 {
   namespace LinearOperator
@@ -1100,7 +831,7 @@ namespace
 
 
 /**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
  * A function that encapsulates generic @p matrix objects that act on a
  * compatible Vector type into a LinearOperator. The LinearOperator object
@@ -1166,7 +897,7 @@ LinearOperator<Range, Domain> linear_operator(const Matrix &matrix)
 
 
 /**
- * \relates LinearOperator
+ * @relates LinearOperator
  *
  * Variant of above function that takes an operator object @p
  * operator_exemplar as an additional reference. This object is used to
@@ -1215,6 +946,303 @@ linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
   return return_op;
 }
 
+//@}
+
+
+/**
+ * @name Creation of LinearOperator block structures
+ */
+//@{
+
+/**
+ * @relates LinearOperator
+ *
+ * A function that encapsulates a given collection @p ops of
+ * LinearOperators into a block structure. Hereby, it is assumed that Range
+ * and Domain are blockvectors, i.e., derived from @ref BlockVectorBase.
+ * The individual linear operators in @p ops must act on a the underlying
+ * vector type of the block vectors, i.e., on Domain::BlockType yielding a
+ * result in Range::BlockType.
+ *
+ * The list @p ops is best passed as an initializer list. Consider for
+ * example a linear operator block (acting on Vector<double>)
+ * @code
+ *  op_a00 | op_a01
+ *         |
+ *  ---------------
+ *         |
+ *  op_a10 | op_a11
+ * @endcode
+ * The coresponding block_operator invocation takes the form
+ * @code
+ * block_operator<2, 2, BlockVector<double>>({op_a00, op_a01, op_a10, op_a11});
+ * @endcode
+ *
+ * @ingroup LAOperators
+ */
+template <unsigned int m, unsigned int n,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+LinearOperator<Range, Domain>
+block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &ops)
+{
+  static_assert(m > 0 && n > 0,
+                "a blocked LinearOperator must consist of at least one block");
+
+  LinearOperator<Range, Domain> return_op;
+
+  return_op.reinit_range_vector = [ops](Range &v, bool fast)
+  {
+    // Reinitialize the block vector to m blocks:
+    v.reinit(m);
+
+    // And reinitialize every individual block with reinit_range_vectors:
+    for (unsigned int i = 0; i < m; ++i)
+      ops[i][0].reinit_range_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.reinit_domain_vector = [ops](Domain &v, bool fast)
+  {
+    // Reinitialize the block vector to n blocks:
+    v.reinit(n);
+
+    // And reinitialize every individual block with reinit_domain_vectors:
+    for (unsigned int i = 0; i < n; ++i)
+      ops[0][i].reinit_domain_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.vmult = [ops](Range &v, const Domain &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
+
+    for (unsigned int i = 0; i < m; ++i)
+      {
+        ops[i][0].vmult(v.block(i), u.block(0));
+        for (unsigned int j = 1; j < n; ++j)
+          ops[i][j].vmult_add(v.block(i), u.block(j));
+      }
+  };
+
+  return_op.vmult_add = [ops](Range &v, const Domain &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
+
+    for (unsigned int i = 0; i < m; ++i)
+      for (unsigned int j = 0; j < n; ++j)
+        ops[i][j].vmult_add(v.block(i), u.block(j));
+  };
+
+  return_op.Tvmult = [ops](Domain &v, const Range &u)
+  {
+    Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < n; ++i)
+      {
+        ops[0][i].Tvmult(v.block(i), u.block(0));
+        for (unsigned int j = 1; j < m; ++j)
+          ops[j][i].Tvmult_add(v.block(i), u.block(j));
+      }
+  };
+
+  return_op.Tvmult_add = [ops](Domain &v, const Range &u)
+  {
+    Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < n; ++i)
+      for (unsigned int j = 0; j < m; ++j)
+        ops[j][i].Tvmult_add(v.block(i), u.block(j));
+  };
+
+  return return_op;
+}
+
+
+/**
+ * @relates LinearOperator
+ *
+ * A variant of above function that builds up a block diagonal linear
+ * operator from an array @p ops of diagonal elements (off-diagonal blocks
+ * are assumed to be 0).
+ *
+ * The list @p ops is best passed as an initializer list. Consider for
+ * example a linear operator block (acting on Vector<double>)
+ * <code>diag(op_a0, op_a1, ..., op_am)</code>. The coresponding
+ * block_operator invocation takes the form
+ * @code
+ * block_diagonal_operator<m, BlockVector<double>>({op_00, op_a1, ..., op_am});
+ * @endcode
+ *
+ * @ingroup LAOperators
+ */
+template <unsigned int m,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+LinearOperator<Range, Domain>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &ops)
+{
+  static_assert(m > 0,
+                "a blockdiagonal LinearOperator must consist of at least one block");
+
+  LinearOperator<Range, Domain> return_op;
+
+  return_op.reinit_range_vector = [ops](Range &v, bool fast)
+  {
+    // Reinitialize the block vector to m blocks:
+    v.reinit(m);
+
+    // And reinitialize every individual block with reinit_range_vectors:
+    for (unsigned int i = 0; i < m; ++i)
+      ops[i].reinit_range_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.reinit_domain_vector = [ops](Domain &v, bool fast)
+  {
+    // Reinitialize the block vector to m blocks:
+    v.reinit(m);
+
+    // And reinitialize every individual block with reinit_domain_vectors:
+    for (unsigned int i = 0; i < m; ++i)
+      ops[i].reinit_domain_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.vmult = [ops](Range &v, const Domain &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      ops[i].vmult(v.block(i), u.block(i));
+  };
+
+  return_op.vmult_add = [ops](Range &v, const Domain &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      ops[i].vmult_add(v.block(i), u.block(i));
+  };
+
+  return_op.Tvmult = [ops](Domain &v, const Range &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      ops[i].Tvmult(v.block(i), u.block(i));
+  };
+
+  return_op.Tvmult_add = [ops](Domain &v, const Range &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      ops[i].Tvmult_add(v.block(i), u.block(i));
+  };
+
+  return return_op;
+}
+
+
+/**
+ * @relates LinearOperator
+ *
+ * A variant of above function that only takes a single LinearOperator
+ * argument @p op and creates a blockdiagonal linear operator with @p m
+ * copies of it.
+ *
+ * @ingroup LAOperators
+ */
+template <unsigned int m,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+LinearOperator<Range, Domain>
+block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op)
+{
+
+  static_assert(m > 0,
+                "a blockdiagonal LinearOperator must consist of at least "
+                "one block");
+
+  LinearOperator<Range, Domain> return_op;
+
+  return_op.reinit_range_vector = [op](Range &v, bool fast)
+  {
+    // Reinitialize the block vector to m blocks:
+    v.reinit(m);
+
+    // And reinitialize every individual block with reinit_range_vectors:
+    for (unsigned int i = 0; i < m; ++i)
+      op.reinit_range_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.reinit_domain_vector = [op](Domain &v, bool fast)
+  {
+    // Reinitialize the block vector to m blocks:
+    v.reinit(m);
+
+    // And reinitialize every individual block with reinit_domain_vectors:
+    for (unsigned int i = 0; i < m; ++i)
+      op.reinit_domain_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.vmult = [op](Range &v, const Domain &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      op.vmult(v.block(i), u.block(i));
+  };
+
+  return_op.vmult_add = [op](Range &v, const Domain &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      op.vmult_add(v.block(i), u.block(i));
+  };
+
+  return_op.Tvmult = [op](Domain &v, const Range &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      op.Tvmult(v.block(i), u.block(i));
+  };
+
+  return_op.Tvmult_add = [op](Domain &v, const Range &u)
+  {
+    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
+    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+
+    for (unsigned int i = 0; i < m; ++i)
+      op.Tvmult_add(v.block(i), u.block(i));
+  };
+
+  return return_op;
+}
+
+//@}
 
 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.