]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove deprecated class FilteredMatrix 9786/head
authorDaniel Arndt <arndtd@ornl.gov>
Tue, 31 Mar 2020 00:20:57 +0000 (20:20 -0400)
committerDaniel Arndt <arndtd@ornl.gov>
Wed, 1 Apr 2020 15:27:13 +0000 (11:27 -0400)
29 files changed:
doc/news/changes/incompatibilities/20200401DanielArndt [new file with mode: 0644]
include/deal.II/lac/filtered_matrix.h [deleted file]
include/deal.II/numerics/matrix_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/grid/grid_tools.cc
tests/data_out/data_out_curved_cells.cc
tests/grid/grid_transform_02.output
tests/lac/filtered_matrix.cc [deleted file]
tests/lac/filtered_matrix.output [deleted file]
tests/lac/filtered_matrix_01.cc [deleted file]
tests/lac/filtered_matrix_01.output [deleted file]
tests/lac/filtered_matrix_02.cc [deleted file]
tests/lac/filtered_matrix_02.output [deleted file]
tests/lac/filtered_matrix_03.cc [deleted file]
tests/lac/filtered_matrix_03.output [deleted file]
tests/lac/filtered_matrix_04.cc [deleted file]
tests/lac/filtered_matrix_04.output [deleted file]
tests/lac/filtered_matrix_05.cc [deleted file]
tests/lac/filtered_matrix_05.output [deleted file]
tests/lac/filtered_matrix_06.cc [deleted file]
tests/lac/filtered_matrix_06.output [deleted file]
tests/mappings/mapping_q1_eulerian_01.cc
tests/mappings/mapping_q_eulerian_01.cc
tests/mappings/mapping_q_eulerian_02.cc
tests/mappings/mapping_q_eulerian_03.cc
tests/mappings/mapping_q_eulerian_07.cc
tests/mappings/mapping_q_eulerian_08.cc
tests/mappings/mapping_q_eulerian_11.cc
tests/matrix_free/interpolate_to_mg.cc

diff --git a/doc/news/changes/incompatibilities/20200401DanielArndt b/doc/news/changes/incompatibilities/20200401DanielArndt
new file mode 100644 (file)
index 0000000..4ddea85
--- /dev/null
@@ -0,0 +1,4 @@
+Removed: The deprecated FilteredMatrix class has been removed.
+Use the LinearOperator class instead.
+<br>
+(Daniel Arndt, 2020/04/01)
diff --git a/include/deal.II/lac/filtered_matrix.h b/include/deal.II/lac/filtered_matrix.h
deleted file mode 100644 (file)
index 0fe20a7..0000000
+++ /dev/null
@@ -1,1008 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2001 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii_filtered_matrix_h
-#  define dealii_filtered_matrix_h
-
-
-
-#  include <deal.II/base/config.h>
-
-#  include <deal.II/base/memory_consumption.h>
-#  include <deal.II/base/smartpointer.h>
-#  include <deal.II/base/thread_management.h>
-
-#  include <deal.II/lac/pointer_matrix.h>
-#  include <deal.II/lac/vector_memory.h>
-
-#  include <algorithm>
-#  include <vector>
-
-DEAL_II_NAMESPACE_OPEN
-
-// Forward declaration
-#  ifndef DOXYGEN
-template <class VectorType>
-class FilteredMatrixBlock;
-#  endif
-
-/*! @addtogroup Matrix2
- *@{
- */
-
-
-/**
- * This class is a wrapper for linear systems of equations with simple
- * equality constraints fixing individual degrees of freedom to a certain
- * value such as when using Dirichlet boundary values.
- *
- * In order to accomplish this, the vmult(), Tvmult(), vmult_add() and
- * Tvmult_add functions modify the same function of the original matrix such
- * as if all constrained entries of the source vector were zero. Additionally,
- * all constrained entries of the destination vector are set to zero.
- *
- * <h3>Usage</h3>
- *
- * Usage is simple: create an object of this type, point it to a matrix that
- * shall be used for $A$ above (either through the constructor, the copy
- * constructor, or the set_referenced_matrix() function), specify the list of
- * boundary values or other constraints (through the add_constraints()
- * function), and then for each required solution modify the right hand side
- * vector (through apply_constraints()) and use this object as matrix object
- * in a linear solver. As linear solvers should only use vmult() and
- * residual() functions of a matrix class, this class should be as good a
- * matrix as any other for that purpose.
- *
- * Furthermore, also the precondition_Jacobi() function is provided (since the
- * computation of diagonal elements of the filtered matrix $A_X$ is simple),
- * so you can use this as a preconditioner. Some other functions useful for
- * matrices are also available.
- *
- * A typical code snippet showing the above steps is as follows:
- * @code
- *   // set up sparse matrix A and right hand side b somehow
- *   ...
- *
- *   // initialize filtered matrix with matrix and boundary value constraints
- *   FilteredMatrix<Vector<double> > filtered_A (A);
- *   filtered_A.add_constraints (boundary_values);
- *
- *   // set up a linear solver
- *   SolverControl control (1000, 1.e-10, false, false);
- *   GrowingVectorMemory<Vector<double> > mem;
- *   SolverCG<Vector<double> > solver (control, mem);
- *
- *   // set up a preconditioner object
- *   PreconditionJacobi<SparseMatrix<double> > prec;
- *   prec.initialize (A, 1.2);
- *   FilteredMatrix<Vector<double> > filtered_prec (prec);
- *   filtered_prec.add_constraints (boundary_values);
- *
- *   // compute modification of right hand side
- *   filtered_A.apply_constraints (b, true);
- *
- *   // solve for solution vector x
- *   solver.solve (filtered_A, x, b, filtered_prec);
- * @endcode
- *
- *
- * <h3>Connection to other classes</h3>
- *
- * The function MatrixTools::apply_boundary_values() does exactly the same
- * that this class does, except for the fact that that function actually
- * modifies the matrix. Consequently, it is only possible to solve with a
- * matrix to which MatrixTools::apply_boundary_values() was applied for one
- * right hand side and one set of boundary values since the modification of
- * the right hand side depends on the original matrix.
- *
- * While this is a feasible method in cases where only one solution of the
- * linear system is required, for example in solving linear stationary
- * systems, one would often like to have the ability to solve multiple times
- * with the same matrix in nonlinear problems (where one often does not want
- * to update the Hessian between Newton steps, despite having different right
- * hand sides in subsequent steps) or time dependent problems, without having
- * to re-assemble the matrix or copy it to temporary matrices with which one
- * then can work. For these cases, this class is meant.
- *
- *
- * <h3>Some background</h3> Mathematically speaking, it is used to represent a
- * system of linear equations $Ax=b$ with the constraint that $B_D x = g_D$,
- * where $B_D$ is a rectangular matrix with exactly one $1$ in each row, and
- * these $1$s in those columns representing constrained degrees of freedom
- * (e.g. for Dirichlet boundary nodes, thus the index $D$) and zeroes for all
- * other diagonal entries, and $g_D$ having the requested nodal values for
- * these constrained nodes. Thus, the underdetermined equation $B_D x = g_D$
- * fixes only the constrained nodes and does not impose any condition on the
- * others. We note that $B_D B_D^T = 1_D$, where $1_D$ is the identity matrix
- * with dimension as large as the number of constrained degrees of freedom.
- * Likewise, $B_D^T B_D$ is the diagonal matrix with diagonal entries $0$ or
- * $1$ that, when applied to a vector, leaves all constrained nodes untouched
- * and deletes all unconstrained ones.
- *
- * For solving such a system of equations, we first write down the Lagrangian
- * $L=1/2 x^T A x - x^T b + l^T B_D x$, where $l$ is a Lagrange multiplier for
- * the constraints. The stationarity condition then reads
- * @code
- * [ A   B_D^T ] [x] = [b  ]
- * [ B_D 0     ] [l] = [g_D]
- * @endcode
- *
- * The first equation then reads $B_D^T l = b-Ax$. On the other hand, if we
- * left-multiply the first equation by $B_D^T B_D$, we obtain $B_D^T B_D A x +
- * B_D^T l = B_D^T B_D b$ after equating $B_D B_D^T$ to the identity matrix.
- * Inserting the previous equality, this yields $(A - B_D^T B_D A) x = (1 -
- * B_D^T B_D)b$. Since $x=(1 - B_D^T B_D) x + B_D^T B_D x = (1 - B_D^T B_D) x
- * + B_D^T g_D$, we can restate the linear system: $A_D x = (1 - B_D^T B_D)b -
- * (1 - B_D^T B_D) A B^T g_D$, where $A_D = (1 - B_D^T B_D) A (1 - B_D^T B_D)$
- * is the matrix where all rows and columns corresponding to constrained nodes
- * have been deleted.
- *
- * The last system of equation only defines the value of the unconstrained
- * nodes, while the constrained ones are determined by the equation $B_D x =
- * g_D$. We can combine these two linear systems by using the zeroed out rows
- * of $A_D$: if we set the diagonal to $1$ and the corresponding zeroed out
- * element of the right hand side to that of $g_D$, then this fixes the
- * constrained elements as well. We can write this as follows: $A_X x = (1 -
- * B_D^T B_D)b - (1 - B_D^T B_D) A B^T g_D + B_D^T g_D$, where $A_X = A_D +
- * B_D^T B_D$. Note that the two parts of the latter matrix operate on
- * disjoint subspaces (the first on the unconstrained nodes, the latter on the
- * constrained ones).
- *
- * In iterative solvers, it is not actually necessary to compute $A_X$
- * explicitly, since only matrix-vector operations need to be performed. This
- * can be done in a three-step procedure that first clears all elements in the
- * incoming vector that belong to constrained nodes, then performs the product
- * with the matrix $A$, then clears again. This class is a wrapper to this
- * procedure, it takes a pointer to a matrix with which to perform matrix-
- * vector products, and does the cleaning of constrained elements itself. This
- * class therefore implements an overloaded @p vmult function that does the
- * matrix-vector product, as well as @p Tvmult for transpose matrix-vector
- * multiplication and @p residual for residual computation, and can thus be
- * used as a matrix replacement in linear solvers.
- *
- * It also has the ability to generate the modification of the right hand
- * side, through the apply_constraints() function.
- *
- *
- *
- * <h3>Template arguments</h3>
- *
- * This class takes as template arguments a matrix and a vector class. The
- * former must provide @p vmult, @p vmult_add,  @p Tvmult, and @p residual
- * member function that operate on the vector type (the second template
- * argument). The latter template parameter must provide access to individual
- * elements through <tt>operator()</tt>, assignment through
- * <tt>operator=</tt>.
- *
- *
- * <h3>Thread-safety</h3>
- *
- * The functions that operate as a matrix and do not change the internal state
- * of this object are synchronized and thus threadsafe. Consequently, you do
- * not need to serialize calls to @p vmult or @p residual .
- *
- * @author Wolfgang Bangerth 2001, Luca Heltai 2006, Guido Kanschat 2007, 2008
- *
- * @deprecated Use LinearOperator instead. See the documentation of
- * constrained_linear_operator().
- */
-template <typename VectorType>
-class DEAL_II_DEPRECATED FilteredMatrix : public Subscriptor
-{
-public:
-  class const_iterator;
-
-  /**
-   * Declare the type of container size.
-   */
-  using size_type = types::global_dof_index;
-
-  /**
-   * Accessor class for iterators
-   */
-  class Accessor
-  {
-    /**
-     * Constructor. Since we use accessors only for read access, a const
-     * matrix pointer is sufficient.
-     */
-    Accessor(const FilteredMatrix<VectorType> *matrix, const size_type index);
-
-  public:
-    /**
-     * Row number of the element represented by this object.
-     */
-    size_type
-    row() const;
-
-    /**
-     * Column number of the element represented by this object.
-     */
-    size_type
-    column() const;
-
-    /**
-     * Value of the right hand side for this row.
-     */
-    double
-    value() const;
-
-  private:
-    /**
-     * Advance to next entry
-     */
-    void
-    advance();
-
-    /**
-     * The matrix accessed.
-     */
-    const FilteredMatrix<VectorType> *matrix;
-
-    /**
-     * Current row number.
-     */
-    size_type index;
-
-    // Make enclosing class a friend.
-    friend class const_iterator;
-  };
-
-  /**
-   * Standard-conforming iterator.
-   */
-  class const_iterator
-  {
-  public:
-    /**
-     * Constructor.
-     */
-    const_iterator(const FilteredMatrix<VectorType> *matrix,
-                   const size_type                   index);
-
-    /**
-     * Prefix increment.
-     */
-    const_iterator &
-    operator++();
-
-    /**
-     * Postfix increment.
-     */
-    const_iterator &
-    operator++(int);
-
-    /**
-     * Dereferencing operator.
-     */
-    const Accessor &operator*() const;
-
-    /**
-     * Dereferencing operator.
-     */
-    const Accessor *operator->() const;
-
-    /**
-     * Comparison. True, if both iterators point to the same matrix position.
-     */
-    bool
-    operator==(const const_iterator &) const;
-    /**
-     * Inverse of <tt>==</tt>.
-     */
-    bool
-    operator!=(const const_iterator &) const;
-
-    /**
-     * Comparison operator. Result is true if either the first row number is
-     * smaller or if the row numbers are equal and the first index is smaller.
-     */
-    bool
-    operator<(const const_iterator &) const;
-
-    /**
-     * Comparison operator. Compares just the other way around than the
-     * operator above.
-     */
-    bool
-    operator>(const const_iterator &) const;
-
-  private:
-    /**
-     * Store an object of the accessor class.
-     */
-    Accessor accessor;
-  };
-
-  /**
-   * Typedef defining a type that represents a pair of degree of freedom index
-   * and the value it shall have.
-   */
-  using IndexValuePair = std::pair<size_type, double>;
-
-  /**
-   * @name Constructors and initialization
-   */
-  //@{
-  /**
-   * Default constructor. You will have to set the matrix to be used later
-   * using initialize().
-   */
-  FilteredMatrix();
-
-  /**
-   * Copy constructor. Use the matrix and the constraints set in the given
-   * object for the present one as well.
-   */
-  FilteredMatrix(const FilteredMatrix &fm);
-
-  /**
-   * Constructor. Use the given matrix for future operations.
-   *
-   * @arg @p m: The matrix being used in multiplications.
-   *
-   * @arg @p expect_constrained_source: See documentation of
-   * #expect_constrained_source.
-   */
-  template <typename MatrixType>
-  FilteredMatrix(const MatrixType &matrix,
-                 const bool        expect_constrained_source = false);
-
-  /**
-   * Copy operator. Take over matrix and constraints from the other object.
-   */
-  FilteredMatrix &
-  operator=(const FilteredMatrix &fm);
-
-  /**
-   * Set the matrix to be used further on. You will probably also want to call
-   * the clear_constraints() function if constraints were previously added.
-   *
-   * @arg @p m: The matrix being used in multiplications.
-   *
-   * @arg @p expect_constrained_source: See documentation of
-   * #expect_constrained_source.
-   */
-  template <typename MatrixType>
-  void
-  initialize(const MatrixType &m, const bool expect_constrained_source = false);
-
-  /**
-   * Delete all constraints and the matrix pointer.
-   */
-  void
-  clear();
-  //@}
-  /**
-   * @name Managing constraints
-   */
-  //@{
-  /**
-   * Add the constraint that the value with index <tt>i</tt> should have the
-   * value <tt>v</tt>.
-   */
-  void
-  add_constraint(const size_type i, const double v);
-
-  /**
-   * Add a list of constraints to the ones already managed by this object. The
-   * actual data type of this list must be so that dereferenced iterators are
-   * pairs of indices and the corresponding values to be enforced on the
-   * respective solution vector's entry. Thus, the data type might be, for
-   * example, a @p std::list or @p std::vector of IndexValuePair objects, but
-   * also a <tt>std::map<unsigned, double></tt>.
-   *
-   * The second component of these pairs will only be used in
-   * apply_constraints(). The first is used to set values to zero in matrix
-   * vector multiplications.
-   *
-   * It is an error if the argument contains an entry for a degree of freedom
-   * that has already been constrained previously.
-   */
-  template <class ConstraintList>
-  void
-  add_constraints(const ConstraintList &new_constraints);
-
-  /**
-   * Delete the list of constraints presently in use.
-   */
-  void
-  clear_constraints();
-  //@}
-  /**
-   * Vector operations
-   */
-  //@{
-  /**
-   * Apply the constraints to a right hand side vector. This needs to be done
-   * before starting to solve with the filtered matrix. If the matrix is
-   * symmetric (i.e. the matrix itself, not only its sparsity pattern), set
-   * the second parameter to @p true to use a faster algorithm. Note: This
-   * method is deprecated as matrix_is_symmetric parameter is no longer used.
-   */
-  DEAL_II_DEPRECATED
-  void
-  apply_constraints(VectorType &v, const bool matrix_is_symmetric) const;
-  /**
-   * Apply the constraints to a right hand side vector. This needs to be done
-   * before starting to solve with the filtered matrix.
-   */
-  void
-  apply_constraints(VectorType &v) const;
-
-  /**
-   * Matrix-vector multiplication: this operation performs pre_filter(),
-   * multiplication with the stored matrix and post_filter() in that order.
-   */
-  void
-  vmult(VectorType &dst, const VectorType &src) const;
-
-  /**
-   * Matrix-vector multiplication: this operation performs pre_filter(),
-   * transposed multiplication with the stored matrix and post_filter() in
-   * that order.
-   */
-  void
-  Tvmult(VectorType &dst, const VectorType &src) const;
-
-  /**
-   * Adding matrix-vector multiplication.
-   *
-   * @note The result vector of this multiplication will have the constraint
-   * entries set to zero, independent of the previous value of <tt>dst</tt>.
-   * We expect that in most cases this is the required behavior.
-   */
-  void
-  vmult_add(VectorType &dst, const VectorType &src) const;
-
-  /**
-   * Adding transpose matrix-vector multiplication:
-   *
-   * @note The result vector of this multiplication will have the constraint
-   * entries set to zero, independent of the previous value of <tt>dst</tt>.
-   * We expect that in most cases this is the required behavior.
-   */
-  void
-  Tvmult_add(VectorType &dst, const VectorType &src) const;
-  //@}
-
-  /**
-   * @name Iterators
-   */
-  //@{
-  /**
-   * Iterator to the first constraint.
-   */
-  const_iterator
-  begin() const;
-  /**
-   * Final iterator.
-   */
-  const_iterator
-  end() const;
-  //@}
-
-  /**
-   * Determine an estimate for the memory consumption (in bytes) of this
-   * object. Since we are not the owner of the matrix referenced, its memory
-   * consumption is not included.
-   */
-  std::size_t
-  memory_consumption() const;
-
-private:
-  /**
-   * Determine, whether multiplications can expect that the source vector has
-   * all constrained entries set to zero.
-   *
-   * If so, the auxiliary vector can be avoided and memory as well as time can
-   * be saved.
-   *
-   * We expect this for instance in Newton's method, where the residual
-   * already should be zero on constrained nodes. This is, because there is no
-   * test function in these nodes.
-   */
-  bool expect_constrained_source;
-
-  /**
-   * Declare an abbreviation for an iterator into the array constraint pairs,
-   * since that data type is so often used and is rather awkward to write out
-   * each time.
-   */
-  using const_index_value_iterator =
-    typename std::vector<IndexValuePair>::const_iterator;
-
-  /**
-   * Helper class used to sort pairs of indices and values. Only the index is
-   * considered as sort key.
-   */
-  struct PairComparison
-  {
-    /**
-     * Function comparing the pairs @p i1 and @p i2 for their keys.
-     */
-    bool
-    operator()(const IndexValuePair &i1, const IndexValuePair &i2) const;
-  };
-
-  /**
-   * Pointer to the sparsity pattern used for this matrix.
-   */
-  std::shared_ptr<PointerMatrixBase<VectorType>> matrix;
-
-  /**
-   * Sorted list of pairs denoting the index of the variable and the value to
-   * which it shall be fixed.
-   */
-  std::vector<IndexValuePair> constraints;
-
-  /**
-   * Do the pre-filtering step, i.e. zero out those components that belong to
-   * constrained degrees of freedom.
-   */
-  void
-  pre_filter(VectorType &v) const;
-
-  /**
-   * Do the postfiltering step, i.e. set constrained degrees of freedom to the
-   * value of the input vector, as the matrix contains only ones on the
-   * diagonal for these degrees of freedom.
-   */
-  void
-  post_filter(const VectorType &in, VectorType &out) const;
-
-  friend class Accessor;
-  // FilteredMatrixBlock accesses pre_filter() and post_filter().
-  friend class FilteredMatrixBlock<VectorType>;
-};
-
-/*@}*/
-/*---------------------- Inline functions -----------------------------------*/
-
-
-//--------------------------------Iterators--------------------------------------//
-
-template <typename VectorType>
-inline FilteredMatrix<VectorType>::Accessor::Accessor(
-  const FilteredMatrix<VectorType> *matrix,
-  const size_type                   index)
-  : matrix(matrix)
-  , index(index)
-{
-  AssertIndexRange(index, matrix->constraints.size() + 1);
-}
-
-
-
-template <typename VectorType>
-inline types::global_dof_index
-FilteredMatrix<VectorType>::Accessor::row() const
-{
-  return matrix->constraints[index].first;
-}
-
-
-
-template <typename VectorType>
-inline types::global_dof_index
-FilteredMatrix<VectorType>::Accessor::column() const
-{
-  return matrix->constraints[index].first;
-}
-
-
-
-template <typename VectorType>
-inline double
-FilteredMatrix<VectorType>::Accessor::value() const
-{
-  return matrix->constraints[index].second;
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::Accessor::advance()
-{
-  Assert(index < matrix->constraints.size(), ExcIteratorPastEnd());
-  ++index;
-}
-
-
-
-template <typename VectorType>
-inline FilteredMatrix<VectorType>::const_iterator::const_iterator(
-  const FilteredMatrix<VectorType> *matrix,
-  const size_type                   index)
-  : accessor(matrix, index)
-{}
-
-
-
-template <typename VectorType>
-inline typename FilteredMatrix<VectorType>::const_iterator &
-FilteredMatrix<VectorType>::const_iterator::operator++()
-{
-  accessor.advance();
-  return *this;
-}
-
-
-template <typename number>
-inline const typename FilteredMatrix<number>::Accessor &
-  FilteredMatrix<number>::const_iterator::operator*() const
-{
-  return accessor;
-}
-
-
-template <typename number>
-inline const typename FilteredMatrix<number>::Accessor *
-  FilteredMatrix<number>::const_iterator::operator->() const
-{
-  return &accessor;
-}
-
-
-template <typename number>
-inline bool
-FilteredMatrix<number>::const_iterator::
-operator==(const const_iterator &other) const
-{
-  return (accessor.index == other.accessor.index &&
-          accessor.matrix == other.accessor.matrix);
-}
-
-
-template <typename number>
-inline bool
-FilteredMatrix<number>::const_iterator::
-operator!=(const const_iterator &other) const
-{
-  return !(*this == other);
-}
-
-
-
-//------------------------- FilteredMatrix ----------------------------------//
-
-template <typename number>
-inline typename FilteredMatrix<number>::const_iterator
-FilteredMatrix<number>::begin() const
-{
-  return const_iterator(this, 0);
-}
-
-
-template <typename number>
-inline typename FilteredMatrix<number>::const_iterator
-FilteredMatrix<number>::end() const
-{
-  return const_iterator(this, constraints.size());
-}
-
-
-template <typename VectorType>
-inline bool
-FilteredMatrix<VectorType>::PairComparison::
-operator()(const IndexValuePair &i1, const IndexValuePair &i2) const
-{
-  return (i1.first < i2.first);
-}
-
-
-
-template <typename VectorType>
-template <typename MatrixType>
-inline void
-FilteredMatrix<VectorType>::initialize(const MatrixType &m, bool ecs)
-{
-  matrix.reset(new_pointer_matrix_base(m, VectorType()));
-
-  expect_constrained_source = ecs;
-}
-
-
-
-template <typename VectorType>
-inline FilteredMatrix<VectorType>::FilteredMatrix()
-  : expect_constrained_source(false)
-{}
-
-
-
-template <typename VectorType>
-inline FilteredMatrix<VectorType>::FilteredMatrix(const FilteredMatrix &fm)
-  : Subscriptor()
-  , expect_constrained_source(fm.expect_constrained_source)
-  , matrix(fm.matrix)
-  , constraints(fm.constraints)
-{}
-
-
-
-template <typename VectorType>
-template <typename MatrixType>
-inline FilteredMatrix<VectorType>::FilteredMatrix(const MatrixType &m,
-                                                  const bool        ecs)
-  : expect_constrained_source(false)
-{
-  initialize(m, ecs);
-}
-
-
-
-template <typename VectorType>
-inline FilteredMatrix<VectorType> &
-FilteredMatrix<VectorType>::operator=(const FilteredMatrix &fm)
-{
-  matrix                    = fm.matrix;
-  expect_constrained_source = fm.expect_constrained_source;
-  constraints               = fm.constraints;
-  return *this;
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::add_constraint(const size_type index,
-                                           const double    value)
-{
-  // add new constraint to end
-  constraints.push_back(IndexValuePair(index, value));
-}
-
-
-
-template <typename VectorType>
-template <class ConstraintList>
-inline void
-FilteredMatrix<VectorType>::add_constraints(
-  const ConstraintList &new_constraints)
-{
-  // add new constraints to end
-  const size_type old_size = constraints.size();
-  constraints.reserve(old_size + new_constraints.size());
-  constraints.insert(constraints.end(),
-                     new_constraints.begin(),
-                     new_constraints.end());
-  // then merge the two arrays to
-  // form one sorted one
-  std::inplace_merge(constraints.begin(),
-                     constraints.begin() + old_size,
-                     constraints.end(),
-                     PairComparison());
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::clear_constraints()
-{
-  // swap vectors to release memory
-  std::vector<IndexValuePair> empty;
-  constraints.swap(empty);
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::clear()
-{
-  clear_constraints();
-  matrix.reset();
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::apply_constraints(
-  VectorType &v,
-  const bool /* matrix_is_symmetric */) const
-{
-  apply_constraints(v);
-}
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::apply_constraints(VectorType &v) const
-{
-  GrowingVectorMemory<VectorType>            mem;
-  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
-  tmp_vector->reinit(v);
-  const_index_value_iterator       i = constraints.begin();
-  const const_index_value_iterator e = constraints.end();
-  for (; i != e; ++i)
-    {
-      AssertIsFinite(i->second);
-      (*tmp_vector)(i->first) = -i->second;
-    }
-
-  // This vmult is without bc, to get
-  // the rhs correction in a correct
-  // way.
-  matrix->vmult_add(v, *tmp_vector);
-  // finally set constrained
-  // entries themselves
-  for (i = constraints.begin(); i != e; ++i)
-    {
-      AssertIsFinite(i->second);
-      v(i->first) = i->second;
-    }
-}
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::pre_filter(VectorType &v) const
-{
-  // iterate over all constraints and
-  // zero out value
-  const_index_value_iterator       i = constraints.begin();
-  const const_index_value_iterator e = constraints.end();
-  for (; i != e; ++i)
-    v(i->first) = 0;
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::post_filter(const VectorType &in,
-                                        VectorType &      out) const
-{
-  // iterate over all constraints and
-  // set value correctly
-  const_index_value_iterator       i = constraints.begin();
-  const const_index_value_iterator e = constraints.end();
-  for (; i != e; ++i)
-    {
-      AssertIsFinite(in(i->first));
-      out(i->first) = in(i->first);
-    }
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
-{
-  if (!expect_constrained_source)
-    {
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
-      // first copy over src vector and
-      // pre-filter
-      tmp_vector->reinit(src, true);
-      *tmp_vector = src;
-      pre_filter(*tmp_vector);
-      // then let matrix do its work
-      matrix->vmult(dst, *tmp_vector);
-    }
-  else
-    {
-      matrix->vmult(dst, src);
-    }
-
-  // finally do post-filtering
-  post_filter(src, dst);
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
-{
-  if (!expect_constrained_source)
-    {
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
-      // first copy over src vector and
-      // pre-filter
-      tmp_vector->reinit(src, true);
-      *tmp_vector = src;
-      pre_filter(*tmp_vector);
-      // then let matrix do its work
-      matrix->Tvmult(dst, *tmp_vector);
-    }
-  else
-    {
-      matrix->Tvmult(dst, src);
-    }
-
-  // finally do post-filtering
-  post_filter(src, dst);
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::vmult_add(VectorType &      dst,
-                                      const VectorType &src) const
-{
-  if (!expect_constrained_source)
-    {
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
-      // first copy over src vector and
-      // pre-filter
-      tmp_vector->reinit(src, true);
-      *tmp_vector = src;
-      pre_filter(*tmp_vector);
-      // then let matrix do its work
-      matrix->vmult_add(dst, *tmp_vector);
-    }
-  else
-    {
-      matrix->vmult_add(dst, src);
-    }
-
-  // finally do post-filtering
-  post_filter(src, dst);
-}
-
-
-
-template <typename VectorType>
-inline void
-FilteredMatrix<VectorType>::Tvmult_add(VectorType &      dst,
-                                       const VectorType &src) const
-{
-  if (!expect_constrained_source)
-    {
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
-      // first copy over src vector and
-      // pre-filter
-      tmp_vector->reinit(src, true);
-      *tmp_vector = src;
-      pre_filter(*tmp_vector);
-      // then let matrix do its work
-      matrix->Tvmult_add(dst, *tmp_vector);
-    }
-  else
-    {
-      matrix->Tvmult_add(dst, src);
-    }
-
-  // finally do post-filtering
-  post_filter(src, dst);
-}
-
-
-
-template <typename VectorType>
-inline std::size_t
-FilteredMatrix<VectorType>::memory_consumption() const
-{
-  return (MemoryConsumption::memory_consumption(matrix) +
-          MemoryConsumption::memory_consumption(constraints));
-}
-
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
-/*----------------------------   filtered_matrix.h ---------------------------*/
index 588d5962da06a6b3954ffbafc89b491b4502756a..14a8bc559951ed2608e623980041dddb6aa9d8dc 100644 (file)
@@ -735,8 +735,8 @@ namespace MatrixCreator
  * depends on the original matrix, this is not possible without storing the
  * original matrix somewhere and applying the @p apply_boundary_conditions
  * function to a copy of it each time we want to solve. In that case, you can
- * use the FilteredMatrix class in the @p LAC sublibrary. There you can also
- * find a formal (mathematical) description of the process of modifying the
+ * use the constrained_linear_operator() function. In its documentation you can
+ * also find a formal (mathematical) description of the process of modifying the
  * matrix and right hand side vectors for boundary values.
  *
  *
index 366263d6a007ec501504030fd09f09b5d93b9495..8bce666efa878d4e9d86e3676ea5cd7e844429ae 100644 (file)
@@ -60,7 +60,6 @@
 #include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/block_vector.h>
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/la_parallel_block_vector.h>
 #include <deal.II/lac/la_parallel_vector.h>
 #include <deal.II/lac/la_vector.h>
@@ -3367,77 +3366,6 @@ namespace VectorTools
       return std::numeric_limits<number>::min();
     }
 
-
-
-    template <typename number>
-    void
-    invert_mass_matrix(
-      const SparseMatrix<number> &          mass_matrix,
-      const FilteredMatrix<Vector<number>> &filtered_mass_matrix,
-      FilteredMatrix<Vector<number>> &      filtered_preconditioner,
-      const Vector<number> &                rhs,
-      Vector<number> &                      boundary_projection)
-    {
-      // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
-      // steps may not be sufficient, since roundoff errors may accumulate for
-      // badly conditioned matrices. This behavior can be observed, e.g. for
-      // FE_Q_Hierarchical for degree higher than three.
-      ReductionControl control(5 * rhs.size(), 0., 1.e-12, false, false);
-      GrowingVectorMemory<Vector<number>> memory;
-      SolverCG<Vector<number>>            cg(control, memory);
-
-      PreconditionSSOR<SparseMatrix<number>> prec;
-      prec.initialize(mass_matrix, 1.2);
-      filtered_preconditioner.initialize(prec, true);
-      // solve
-      cg.solve(filtered_mass_matrix,
-               boundary_projection,
-               rhs,
-               filtered_preconditioner);
-      filtered_preconditioner.apply_constraints(boundary_projection, true);
-      filtered_preconditioner.clear();
-    }
-
-
-
-    template <typename number>
-    void
-    invert_mass_matrix(
-      const SparseMatrix<number> &          mass_matrix,
-      const FilteredMatrix<Vector<number>> &filtered_mass_matrix,
-      FilteredMatrix<Vector<number>> &      filtered_preconditioner,
-      const Vector<std::complex<number>> &  rhs,
-      Vector<std::complex<number>> &        boundary_projection)
-    {
-      auto solve_for_one_component = [&](const bool real_part) {
-        // copy the real or imaginary part out of the rhs vector
-        Vector<number> rhs_part(rhs.size());
-        for (unsigned int i = 0; i < rhs.size(); ++i)
-          rhs_part(i) = (real_part ? rhs(i).real() : rhs(i).imag());
-
-        // then solve the linear system for this part
-        Vector<number> boundary_projection_part(boundary_projection.size());
-        invert_mass_matrix(mass_matrix,
-                           filtered_mass_matrix,
-                           filtered_preconditioner,
-                           rhs_part,
-                           boundary_projection_part);
-
-        // finally copy the real or imaginary part of the
-        // solution back into the global solution vector
-        for (unsigned int i = 0; i < boundary_projection.size(); ++i)
-          if (real_part == true)
-            boundary_projection(i).real(boundary_projection_part(i));
-          else
-            boundary_projection(i).imag(boundary_projection_part(i));
-      };
-
-      // solve for real and imaginary parts of the solution separately
-      solve_for_one_component(true);
-      solve_for_one_component(false);
-    }
-
-
     template <int dim,
               int spacedim,
               template <int, int> class DoFHandlerType,
@@ -3572,42 +3500,6 @@ namespace VectorTools
         static_cast<const Function<spacedim, number> *>(nullptr),
         component_mapping);
 
-      // For certain weird elements,
-      // there might be degrees of
-      // freedom on the boundary, but
-      // their shape functions do not
-      // have support there. Let's
-      // eliminate them here.
-
-      // The Bogner-Fox-Schmidt element
-      // is an example for those.
-
-      // TODO: Maybe we should figure out if the element really needs this
-
-      FilteredMatrix<Vector<number>> filtered_mass_matrix(mass_matrix, true);
-      FilteredMatrix<Vector<number>> filtered_precondition;
-      std::vector<bool>              excluded_dofs(mass_matrix.m(), false);
-
-      // we assemble mass matrix with unit weight,
-      // thus it will be real-valued irrespectively of the underlying algebra
-      // with positive elements on diagonal.
-      // Thus in order to extend this filtering to complex-algebra simply take
-      // the real-part of element.
-      number max_element = 0.;
-      for (unsigned int i = 0; i < mass_matrix.m(); ++i)
-        if (real_part_bigger_than(mass_matrix.diag_element(i), max_element))
-          max_element = mass_matrix.diag_element(i);
-
-      for (unsigned int i = 0; i < mass_matrix.m(); ++i)
-        if (real_part_bigger_than(1.e-8 * max_element,
-                                  mass_matrix.diag_element(i)))
-          {
-            filtered_mass_matrix.add_constraint(i, 0.);
-            filtered_precondition.add_constraint(i, 0.);
-            mass_matrix.diag_element(i) = 1.;
-            excluded_dofs[i]            = true;
-          }
-
       Vector<number> boundary_projection(rhs.size());
 
       // cannot reduce residual in a useful way if we are close to the square
@@ -3616,16 +3508,11 @@ namespace VectorTools
         boundary_projection = 0;
       else
         {
-          invert_mass_matrix(mass_matrix,
-                             filtered_mass_matrix,
-                             filtered_precondition,
-                             rhs,
-                             boundary_projection);
+          invert_mass_matrix(mass_matrix, rhs, boundary_projection);
         }
       // fill in boundary values
       for (unsigned int i = 0; i < dof_to_boundary_mapping.size(); ++i)
-        if (dof_to_boundary_mapping[i] != numbers::invalid_dof_index &&
-            !excluded_dofs[dof_to_boundary_mapping[i]])
+        if (dof_to_boundary_mapping[i] != numbers::invalid_dof_index)
           {
             AssertIsFinite(boundary_projection(dof_to_boundary_mapping[i]));
 
index 4a4cee06befefb9d40b007b12a059e15780629d0..2fd1c5b98cb9ff47574b98e84cfa6502c936b8a7 100644 (file)
@@ -41,8 +41,8 @@
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
 
+#include <deal.II/lac/constrained_linear_operator.h>
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/precondition.h>
 #include <deal.II/lac/solver_cg.h>
 #include <deal.II/lac/sparse_matrix.h>
@@ -1104,15 +1104,15 @@ namespace GridTools
      * in order to allow parallel execution.
      */
     void
-    laplace_solve(const SparseMatrix<double> &                     S,
-                  const std::map<types::global_dof_index, double> &fixed_dofs,
-                  Vector<double> &                                 u)
+    laplace_solve(const SparseMatrix<double> &     S,
+                  const AffineConstraints<double> &constraints,
+                  Vector<double> &                 u)
     {
-      const unsigned int                       n_dofs = S.n();
-      FilteredMatrix<Vector<double>>           SF(S);
+      const unsigned int n_dofs = S.n();
+      const auto         op     = linear_operator(S);
+      const auto         SF     = constrained_linear_operator(constraints, op);
       PreconditionJacobi<SparseMatrix<double>> prec;
       prec.initialize(S, 1.2);
-      FilteredMatrix<Vector<double>> PF(prec);
 
       SolverControl                       control(n_dofs, 1.e-10, false, false);
       GrowingVectorMemory<Vector<double>> mem;
@@ -1120,9 +1120,11 @@ namespace GridTools
 
       Vector<double> f(n_dofs);
 
-      SF.add_constraints(fixed_dofs);
-      SF.apply_constraints(f, true);
-      solver.solve(SF, u, f, PF);
+      const auto constrained_rhs =
+        constrained_right_hand_side(constraints, op, f);
+      solver.solve(SF, u, constrained_rhs, prec);
+
+      constraints.distribute(u);
     }
   } // namespace
 
@@ -1171,7 +1173,7 @@ namespace GridTools
       StaticMappingQ1<dim>::mapping, dof_handler, quadrature, S, coefficient);
 
     // set up the boundary values for the laplace problem
-    std::map<types::global_dof_index, double>                   fixed_dofs[dim];
+    std::array<AffineConstraints<double>, dim>                  constraints;
     typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
       new_points.end();
 
@@ -1190,14 +1192,20 @@ namespace GridTools
 
             if (map_iter != map_end)
               for (unsigned int i = 0; i < dim; ++i)
-                fixed_dofs[i].insert(std::pair<types::global_dof_index, double>(
-                  cell->vertex_dof_index(vertex_no, 0),
-                  (solve_for_absolute_positions ?
-                     map_iter->second(i) :
-                     map_iter->second(i) - vertex_point[i])));
+                {
+                  constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
+                  constraints[i].set_inhomogeneity(
+                    cell->vertex_dof_index(vertex_no, 0),
+                    (solve_for_absolute_positions ?
+                       map_iter->second(i) :
+                       map_iter->second(i) - vertex_point[i]));
+                }
           }
       }
 
+    for (unsigned int i = 0; i < dim; ++i)
+      constraints[i].close();
+
     // solve the dim problems with different right hand sides.
     Vector<double> us[dim];
     for (unsigned int i = 0; i < dim; ++i)
@@ -1206,7 +1214,7 @@ namespace GridTools
     // solve linear systems in parallel
     Threads::TaskGroup<> tasks;
     for (unsigned int i = 0; i < dim; ++i)
-      tasks += Threads::new_task(&laplace_solve, S, fixed_dofs[i], us[i]);
+      tasks += Threads::new_task(&laplace_solve, S, constraints[i], us[i]);
     tasks.join_all();
 
     // change the coordinates of the points of the triangulation
index 1e590aad0922e869b1306a70520d2cf3c7aa5d4e..f40bf99020684352aef9101a901605679dc799cf 100644 (file)
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
+#include <deal.II/lac/constrained_linear_operator.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
 #include <deal.II/lac/precondition.h>
 #include <deal.II/lac/solver_cg.h>
 #include <deal.II/lac/sparse_matrix.h>
 #include "../tests.h"
 
 
-// this is copied from GridGenerator
 void
-laplace_solve(const SparseMatrix<double> &          S,
-              const std::map<unsigned int, double> &m,
-              Vector<double> &                      u)
+laplace_solve(const SparseMatrix<double> &     S,
+              const AffineConstraints<double> &constraints,
+              Vector<double> &                 u)
 {
-  const unsigned int                       n_dofs = S.n();
-  FilteredMatrix<Vector<double>>           SF(S);
-  PreconditionJacobi<SparseMatrix<double>> prec;
-  prec.initialize(S, 1.2);
-  FilteredMatrix<Vector<double>> PF(prec);
+  const unsigned int n_dofs = S.n();
+  const auto         op     = linear_operator(S);
+  const auto constrained_op = constrained_linear_operator(constraints, op);
 
   SolverControl                       control(10000, 1.e-10, false, false);
   GrowingVectorMemory<Vector<double>> mem;
@@ -74,13 +73,13 @@ laplace_solve(const SparseMatrix<double> &          S,
 
   Vector<double> f(n_dofs);
 
-  SF.add_constraints(m);
-  SF.apply_constraints(f, true);
-  solver.solve(SF, u, f, PF);
+  const auto constrained_rhs = constrained_right_hand_side(constraints, op, f);
+  solver.solve(constrained_op, u, constrained_rhs, PreconditionIdentity());
+  constraints.distribute(u);
 }
 
 
-// create a rinf grid and compute a MappingQEuler to represent the inner
+// create a ring grid and compute a MappingQEuler to represent the inner
 // boundary
 void
 curved_grid(std::ostream &out)
@@ -151,27 +150,25 @@ curved_grid(std::ostream &out)
   MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
   // set up the boundary values for
   // the laplace problem
-  std::vector<std::map<unsigned int, double>> m(2);
-  // fill these maps: on the inner boundary try
+  std::vector<AffineConstraints<double>> m(2);
+  // fill these constraints: on the inner boundary try
   // to approximate a circle, on the outer
   // boundary use straight lines
-  DoFHandler<2>::cell_iterator cell = dof_handler.begin_active(),
-                               endc = dof_handler.end();
   DoFHandler<2>::face_iterator face;
-  for (; cell != endc; ++cell)
+  for (const auto &cell : dof_handler.active_cell_iterators())
     {
       // fix all vertices
       for (const unsigned int vertex_no : GeometryInfo<2>::vertex_indices())
         {
           for (unsigned int i = 0; i < 2; ++i)
-            m[i][cell->vertex_dof_index(vertex_no, 0)] = 0.;
+            m[i].add_line(cell->vertex_dof_index(vertex_no, 0));
         }
 
       if (cell->at_boundary())
         for (const unsigned int face_no : GeometryInfo<2>::face_indices())
           {
             face = cell->face(face_no);
-            // insert a modifiued value for
+            // insert a modified value for
             // the middle point of boundary
             // faces
             if (face->at_boundary())
@@ -179,16 +176,24 @@ curved_grid(std::ostream &out)
                 const double eps = 1e-4;
                 if (std::fabs(face->vertex(1).norm() - r_i) < eps)
                   for (unsigned int i = 0; i < 2; ++i)
-                    m[i][face->dof_index(0)] =
-                      (face->center() * (r_i / face->center().norm() - 1))(i);
+                    {
+                      m[i].add_line(face->dof_index(0));
+                      m[i].set_inhomogeneity(face->dof_index(0),
+                                             (face->center() *
+                                              (r_i / face->center().norm() -
+                                               1))(i));
+                    }
                 else if (std::fabs(face->vertex(1).norm() - r_a) < eps)
                   for (unsigned int i = 0; i < 2; ++i)
-                    m[i][face->dof_index(0)] = 0.;
+                    m[i].add_line(face->dof_index(0));
                 else
                   Assert(false, ExcInternalError());
               }
           }
     }
+  m[0].close();
+  m[1].close();
+
   // solve the 2 problems with
   // different right hand sides.
   Vector<double> us[2];
index b35bbeabbdd29fd6324300d318b7d695de39d3fe..88d3c9a87e307991c5ae648fafbc143b508d6a74 100644 (file)
@@ -2,7 +2,7 @@
 %%Title: deal.II Output
 %%Creator: the deal.II library
 
-%%BoundingBox: 0 0 301 300
+%%BoundingBox: 0 0 301 301
 /m {moveto} bind def
 /x {lineto stroke} bind def
 /b {0 0 0 setrgbcolor} def
 %%EndProlog
 
 0.5 setlinewidth
-b 1.86517e-13 0 m 0 37.5 x
-b 37.5 2.66454e-14 m 38.4292 37.5 x
-b 1.86517e-13 0 m 37.5 2.66454e-14 x
+b 0 0 m 0 37.5 x
+b 37.5 0 m 38.4292 37.5 x
+b 0 0 m 37.5 0 x
 b 0 37.5 m 38.4292 37.5 x
-b 37.5 2.66454e-14 m 38.4292 37.5 x
-b 75 2.66454e-14 m 76.8146 37.5 x
-b 37.5 2.66454e-14 m 75 2.66454e-14 x
+b 37.5 0 m 38.4292 37.5 x
+b 75 0 m 76.8146 37.5 x
+b 37.5 0 m 75 0 x
 b 38.4292 37.5 m 76.8146 37.5 x
-b 75 2.66454e-14 m 76.8146 37.5 x
-b 112.5 2.66454e-14 m 115.008 37.5 x
-b 75 2.66454e-14 m 112.5 2.66454e-14 x
+b 75 0 m 76.8146 37.5 x
+b 112.5 0 m 115.008 37.5 x
+b 75 0 m 112.5 0 x
 b 76.8146 37.5 m 115.008 37.5 x
-b 112.5 2.66454e-14 m 115.008 37.5 x
-b 150 2.66454e-14 m 152.77 37.5 x
-b 112.5 2.66454e-14 m 150 2.66454e-14 x
+b 112.5 0 m 115.008 37.5 x
+b 150 0 m 152.77 37.5 x
+b 112.5 0 m 150 0 x
 b 115.008 37.5 m 152.77 37.5 x
-b 150 2.66454e-14 m 152.77 37.5 x
-b 187.5 2.66454e-14 m 190.008 37.5 x
-b 150 2.66454e-14 m 187.5 2.66454e-14 x
+b 150 0 m 152.77 37.5 x
+b 187.5 0 m 190.008 37.5 x
+b 150 0 m 187.5 0 x
 b 152.77 37.5 m 190.008 37.5 x
-b 187.5 2.66454e-14 m 190.008 37.5 x
-b 225 2.66454e-14 m 226.815 37.5 x
-b 187.5 2.66454e-14 m 225 2.66454e-14 x
+b 187.5 0 m 190.008 37.5 x
+b 225 0 m 226.815 37.5 x
+b 187.5 0 m 225 0 x
 b 190.008 37.5 m 226.815 37.5 x
-b 225 2.66454e-14 m 226.815 37.5 x
-b 262.5 2.66454e-14 m 263.429 37.5 x
-b 225 2.66454e-14 m 262.5 2.66454e-14 x
+b 225 0 m 226.815 37.5 x
+b 262.5 0 m 263.429 37.5 x
+b 225 0 m 262.5 0 x
 b 226.815 37.5 m 263.429 37.5 x
-b 262.5 2.66454e-14 m 263.429 37.5 x
+b 262.5 0 m 263.429 37.5 x
 b 300 0 m 300 37.5 x
-b 262.5 2.66454e-14 m 300 0 x
+b 262.5 0 m 300 0 x
 b 263.429 37.5 m 300 37.5 x
 b 0 37.5 m 0 75 x
 b 38.4292 37.5 m 39.3146 75 x
@@ -234,10 +234,10 @@ b 264.315 225 m 263.429 262.5 x
 b 300 225 m 300 262.5 x
 b 264.315 225 m 300 225 x
 b 263.429 262.5 m 300 262.5 x
-b 0 262.5 m 1.59872e-13 300 x
+b 0 262.5 m 0 300 x
 b 38.4292 262.5 m 37.5 300 x
 b 0 262.5 m 38.4292 262.5 x
-b 1.59872e-13 300 m 37.5 300 x
+b 0 300 m 37.5 300 x
 b 38.4292 262.5 m 37.5 300 x
 b 76.8146 262.5 m 75 300 x
 b 38.4292 262.5 m 76.8146 262.5 x
diff --git a/tests/lac/filtered_matrix.cc b/tests/lac/filtered_matrix.cc
deleted file mode 100644 (file)
index 1b91f02..0000000
+++ /dev/null
@@ -1,99 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2007 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-
-// Test properties of FilteredMatrix and iterators
-
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/linear_operator.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/vector.h>
-
-#include "../tests.h"
-
-
-template <typename VectorType>
-void
-test(const FilteredMatrix<VectorType> &M)
-{
-  deallog << "Iterator";
-
-  unsigned int max = 0;
-  for (typename FilteredMatrix<VectorType>::const_iterator i = M.begin();
-       i != M.end();
-       ++i)
-    {
-      Assert(i->row() == i->column(), ExcInternalError());
-      deallog << ' ' << i->row() << ':' << i->value();
-      max = i->row();
-    }
-  VectorType v(max + 1);
-  VectorType w(max + 1);
-
-  for (unsigned int i = 0; i < v.size(); ++i)
-    v(i) = 31 + i;
-
-  deallog << std::endl << "vmult ";
-
-  w = 0.;
-  M.vmult(w, v);
-  for (unsigned int i = 0; i < v.size(); ++i)
-    deallog << ' ' << w(i);
-
-  deallog << std::endl << "Tvmult";
-
-  w = 0.;
-  M.Tvmult(w, v);
-  for (unsigned int i = 0; i < v.size(); ++i)
-    deallog << ' ' << w(i);
-
-  deallog << std::endl << "vmult_add";
-
-  M.vmult_add(w, v);
-  for (unsigned int i = 0; i < v.size(); ++i)
-    deallog << ' ' << w(i);
-
-  deallog << std::endl << "boundary";
-
-  M.apply_constraints(w, true);
-  for (unsigned int i = 0; i < v.size(); ++i)
-    deallog << ' ' << w(i);
-
-  deallog << std::endl << "boundary";
-
-  M.apply_constraints(w, false);
-  for (unsigned int i = 0; i < v.size(); ++i)
-    deallog << ' ' << w(i);
-
-  deallog << std::endl;
-}
-
-int
-main()
-{
-  initlog();
-
-  PreconditionIdentity           identity;
-  FilteredMatrix<Vector<double>> f;
-  f.initialize(identity, false);
-  test(f);
-
-  for (unsigned int i = 0; i < 5; ++i)
-    f.add_constraint(i * i, i / 2.);
-  test(f);
-
-  f.initialize(identity, true);
-  test(f);
-}
diff --git a/tests/lac/filtered_matrix.output b/tests/lac/filtered_matrix.output
deleted file mode 100644 (file)
index c1fc45f..0000000
+++ /dev/null
@@ -1,19 +0,0 @@
-
-DEAL::Iterator
-DEAL::vmult  31.0000
-DEAL::Tvmult 31.0000
-DEAL::vmult_add 62.0000
-DEAL::boundary 62.0000
-DEAL::boundary 62.0000
-DEAL::Iterator 0:0.00000 1:0.500000 4:1.00000 9:1.50000 16:2.00000
-DEAL::vmult  31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 41.0000 42.0000 43.0000 44.0000 45.0000 46.0000 47.0000
-DEAL::Tvmult 31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 41.0000 42.0000 43.0000 44.0000 45.0000 46.0000 47.0000
-DEAL::vmult_add 31.0000 32.0000 66.0000 68.0000 35.0000 72.0000 74.0000 76.0000 78.0000 40.0000 82.0000 84.0000 86.0000 88.0000 90.0000 92.0000 47.0000
-DEAL::boundary 0.00000 0.500000 66.0000 68.0000 1.00000 72.0000 74.0000 76.0000 78.0000 1.50000 82.0000 84.0000 86.0000 88.0000 90.0000 92.0000 2.00000
-DEAL::boundary 0.00000 0.500000 66.0000 68.0000 1.00000 72.0000 74.0000 76.0000 78.0000 1.50000 82.0000 84.0000 86.0000 88.0000 90.0000 92.0000 2.00000
-DEAL::Iterator 0:0.00000 1:0.500000 4:1.00000 9:1.50000 16:2.00000
-DEAL::vmult  31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 41.0000 42.0000 43.0000 44.0000 45.0000 46.0000 47.0000
-DEAL::Tvmult 31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 41.0000 42.0000 43.0000 44.0000 45.0000 46.0000 47.0000
-DEAL::vmult_add 31.0000 32.0000 66.0000 68.0000 35.0000 72.0000 74.0000 76.0000 78.0000 40.0000 82.0000 84.0000 86.0000 88.0000 90.0000 92.0000 47.0000
-DEAL::boundary 0.00000 0.500000 66.0000 68.0000 1.00000 72.0000 74.0000 76.0000 78.0000 1.50000 82.0000 84.0000 86.0000 88.0000 90.0000 92.0000 2.00000
-DEAL::boundary 0.00000 0.500000 66.0000 68.0000 1.00000 72.0000 74.0000 76.0000 78.0000 1.50000 82.0000 84.0000 86.0000 88.0000 90.0000 92.0000 2.00000
diff --git a/tests/lac/filtered_matrix_01.cc b/tests/lac/filtered_matrix_01.cc
deleted file mode 100644 (file)
index b3958ee..0000000
+++ /dev/null
@@ -1,66 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2007 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#include "../tests.h"
-
-template <typename number>
-void
-checkVmult(FullMatrix<number> &A,
-           Vector<number> &    V,
-           bool                expect_constrained_source = false)
-{
-  deallog << "vmult" << std::endl;
-
-  FilteredMatrix<Vector<double>> F;
-  F.initialize(A, expect_constrained_source);
-  F.add_constraint(0, 1);
-
-  Vector<number> O(A.m());
-
-  F.vmult(O, V);
-
-  for (unsigned int i = 0; i < O.size(); ++i)
-    deallog << O(i) << '\t';
-  deallog << std::endl;
-}
-
-int
-main()
-{
-  std::ofstream logfile("output");
-  deallog << std::fixed;
-  deallog << std::setprecision(4);
-  deallog.attach(logfile);
-
-  const double Adata[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
-
-  FullMatrix<double> A(3, 3);
-
-  A.fill(Adata);
-
-  Vector<double> V1(3);
-  Vector<double> V2(3);
-
-  V1(0) = V2(0) = 1;
-  V1(1) = V2(1) = 2;
-  V1(2) = V2(2) = 3;
-
-  checkVmult<double>(A, V1, false);
-  checkVmult<double>(A, V2, true);
-}
diff --git a/tests/lac/filtered_matrix_01.output b/tests/lac/filtered_matrix_01.output
deleted file mode 100644 (file)
index e893ffe..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL::vmult
-DEAL::1.0000   28.0000 43.0000 
-DEAL::vmult
-DEAL::1.0000   32.0000 50.0000 
diff --git a/tests/lac/filtered_matrix_02.cc b/tests/lac/filtered_matrix_02.cc
deleted file mode 100644 (file)
index 86e6869..0000000
+++ /dev/null
@@ -1,66 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2007 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#include "../tests.h"
-
-template <typename number>
-void
-checkTvmult(FullMatrix<number> &A,
-            Vector<number> &    V,
-            bool                expect_constrained_source = false)
-{
-  deallog << "Tvmult" << std::endl;
-
-  FilteredMatrix<Vector<double>> F;
-  F.initialize(A, expect_constrained_source);
-  F.add_constraint(0, 1);
-
-  Vector<number> O(A.n());
-
-  F.Tvmult(O, V);
-
-  for (unsigned int i = 0; i < O.size(); ++i)
-    deallog << O(i) << '\t';
-  deallog << std::endl;
-}
-
-int
-main()
-{
-  std::ofstream logfile("output");
-  deallog << std::fixed;
-  deallog << std::setprecision(4);
-  deallog.attach(logfile);
-
-  const double Adata[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
-
-  FullMatrix<double> A(3, 3);
-
-  A.fill(Adata);
-
-  Vector<double> V1(3);
-  Vector<double> V2(3);
-
-  V1(0) = V2(0) = 1;
-  V1(1) = V2(1) = 2;
-  V1(2) = V2(2) = 3;
-
-  checkTvmult<double>(A, V1, false);
-  checkTvmult<double>(A, V2, true);
-}
diff --git a/tests/lac/filtered_matrix_02.output b/tests/lac/filtered_matrix_02.output
deleted file mode 100644 (file)
index 5c3f2be..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL::Tvmult
-DEAL::1.0000   34.0000 39.0000 
-DEAL::Tvmult
-DEAL::1.0000   36.0000 42.0000 
diff --git a/tests/lac/filtered_matrix_03.cc b/tests/lac/filtered_matrix_03.cc
deleted file mode 100644 (file)
index 1c334e1..0000000
+++ /dev/null
@@ -1,70 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2007 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#include "../tests.h"
-
-template <typename number>
-void
-checkVmult_Add(FullMatrix<number> &A,
-               Vector<number> &    V,
-               bool                expect_constrained_source = false)
-{
-  deallog << "vmult_add" << std::endl;
-
-  FilteredMatrix<Vector<double>> F;
-  F.initialize(A, expect_constrained_source);
-  F.add_constraint(0, 1);
-
-  Vector<number> O(A.m());
-  for (unsigned int i = 0; i < O.size(); ++i)
-    {
-      O(i) = 1;
-    }
-
-  F.vmult_add(O, V);
-
-  for (unsigned int i = 0; i < O.size(); ++i)
-    deallog << O(i) << '\t';
-  deallog << std::endl;
-}
-
-int
-main()
-{
-  std::ofstream logfile("output");
-  deallog << std::fixed;
-  deallog << std::setprecision(4);
-  deallog.attach(logfile);
-
-  const double Adata[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
-
-  FullMatrix<double> A(3, 3);
-
-  A.fill(Adata);
-
-  Vector<double> V1(3);
-  Vector<double> V2(3);
-
-  V1(0) = V2(0) = 1;
-  V1(1) = V2(1) = 2;
-  V1(2) = V2(2) = 3;
-
-  checkVmult_Add<double>(A, V1, false);
-  checkVmult_Add<double>(A, V2, true);
-}
diff --git a/tests/lac/filtered_matrix_03.output b/tests/lac/filtered_matrix_03.output
deleted file mode 100644 (file)
index 9ea8103..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL::vmult_add
-DEAL::1.0000   29.0000 44.0000 
-DEAL::vmult_add
-DEAL::1.0000   33.0000 51.0000 
diff --git a/tests/lac/filtered_matrix_04.cc b/tests/lac/filtered_matrix_04.cc
deleted file mode 100644 (file)
index 6ad73c9..0000000
+++ /dev/null
@@ -1,70 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2007 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#include "../tests.h"
-
-template <typename number>
-void
-checkTvmult_Add(FullMatrix<number> &A,
-                Vector<number> &    V,
-                bool                expect_constrained_source = false)
-{
-  deallog << "Tvmult_add" << std::endl;
-
-  FilteredMatrix<Vector<double>> F;
-  F.initialize(A, expect_constrained_source);
-  F.add_constraint(0, 1);
-
-  Vector<number> O(A.n());
-  for (unsigned int i = 0; i < O.size(); ++i)
-    {
-      O(i) = 1;
-    }
-
-  F.Tvmult_add(O, V);
-
-  for (unsigned int i = 0; i < O.size(); ++i)
-    deallog << O(i) << '\t';
-  deallog << std::endl;
-}
-
-int
-main()
-{
-  std::ofstream logfile("output");
-  deallog << std::fixed;
-  deallog << std::setprecision(4);
-  deallog.attach(logfile);
-
-  const double Adata[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
-
-  FullMatrix<double> A(3, 3);
-
-  A.fill(Adata);
-
-  Vector<double> V1(3);
-  Vector<double> V2(3);
-
-  V1(0) = V2(0) = 1;
-  V1(1) = V2(1) = 2;
-  V1(2) = V2(2) = 3;
-
-  checkTvmult_Add<double>(A, V1, false);
-  checkTvmult_Add<double>(A, V2, true);
-}
diff --git a/tests/lac/filtered_matrix_04.output b/tests/lac/filtered_matrix_04.output
deleted file mode 100644 (file)
index b455ab3..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL::Tvmult_add
-DEAL::1.0000   35.0000 40.0000 
-DEAL::Tvmult_add
-DEAL::1.0000   37.0000 43.0000 
diff --git a/tests/lac/filtered_matrix_05.cc b/tests/lac/filtered_matrix_05.cc
deleted file mode 100644 (file)
index 4e99667..0000000
+++ /dev/null
@@ -1,64 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2007 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#include "../tests.h"
-
-template <typename number>
-void
-checkApply_Constraints(FullMatrix<number> &A,
-                       Vector<number> &    V,
-                       bool                matrix_is_symmetric = false)
-{
-  deallog << "apply_constraints" << std::endl;
-
-  FilteredMatrix<Vector<double>> F;
-  F.initialize(A);
-  F.add_constraint(0, 1);
-
-  F.apply_constraints(V, matrix_is_symmetric);
-
-  for (unsigned int i = 0; i < V.size(); ++i)
-    deallog << V(i) << '\t';
-  deallog << std::endl;
-}
-
-int
-main()
-{
-  std::ofstream logfile("output");
-  deallog << std::fixed;
-  deallog << std::setprecision(4);
-  deallog.attach(logfile);
-
-  const double Adata[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
-
-  FullMatrix<double> A(3, 3);
-
-  A.fill(Adata);
-
-  Vector<double> V1(3);
-  Vector<double> V2(3);
-
-  V1(0) = V2(0) = 1;
-  V1(1) = V2(1) = 2;
-  V1(2) = V2(2) = 3;
-
-  checkApply_Constraints<double>(A, V1, false);
-  checkApply_Constraints<double>(A, V2, true);
-}
diff --git a/tests/lac/filtered_matrix_05.output b/tests/lac/filtered_matrix_05.output
deleted file mode 100644 (file)
index 664c0a3..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL::apply_constraints
-DEAL::1.0000   -2.0000 -4.0000 
-DEAL::apply_constraints
-DEAL::1.0000   -2.0000 -4.0000 
diff --git a/tests/lac/filtered_matrix_06.cc b/tests/lac/filtered_matrix_06.cc
deleted file mode 100644 (file)
index 4df8737..0000000
+++ /dev/null
@@ -1,224 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2013 - 2018 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.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-
-
-#include <deal.II/base/function_lib.h>
-#include <deal.II/base/quadrature_lib.h>
-
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping_q.h>
-
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/vector_memory.h>
-
-#include <deal.II/numerics/matrix_tools.h>
-#include <deal.II/numerics/vector_tools.h>
-
-#include <string>
-#include <vector>
-
-#include "../tests.h"
-
-
-
-void
-solve_filtered(std::map<types::global_dof_index, double> &bv,
-               SparseMatrix<double> &                     A,
-               Vector<double> &                           u,
-               Vector<double> &                           f)
-{
-  FilteredMatrix<Vector<double>> A1(A);
-  A1.add_constraints(bv);
-
-  SolverControl                            control(1000, 1.e-10, false, false);
-  PrimitiveVectorMemory<Vector<double>>    mem;
-  SolverCG<Vector<double>>                 solver(control, mem);
-  PreconditionJacobi<SparseMatrix<double>> prec;
-  FilteredMatrix<Vector<double>>           fprec;
-  prec.initialize(A, 1.2);
-  fprec.initialize(prec);
-
-  Vector<double> f1(f.size());
-  f1 = f;
-  A1.apply_constraints(f1, true);
-
-  solver.solve(A1, u, f1, fprec);
-
-  for (std::map<types::global_dof_index, double>::const_iterator i = bv.begin();
-       i != bv.end();
-       ++i)
-    AssertThrow(std::fabs(u(i->first) - i->second) < 1e-8, ExcInternalError());
-}
-
-
-
-template <int dim>
-void
-solve_eliminated(std::map<types::global_dof_index, double> &bv,
-                 SparseMatrix<double> &                     A,
-                 Vector<double> &                           u,
-                 Vector<double> &                           f)
-{
-  MatrixTools::apply_boundary_values(bv, A, u, f);
-
-  SolverControl                         control(1000, 1.e-10, false, false);
-  PrimitiveVectorMemory<Vector<double>> mem;
-  SolverCG<Vector<double>>              solver(control, mem);
-  PreconditionJacobi<>                  prec;
-  prec.initialize(A, 1.2);
-
-  solver.solve(A, u, f, prec);
-}
-
-
-
-template <int dim>
-void
-check()
-{
-  Triangulation<dim> tr;
-
-  Functions::CosineFunction<dim> cosine;
-
-  if (dim == 2)
-    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
-  else
-    GridGenerator::hyper_cube(tr, -1, 1);
-  tr.reset_manifold(0);
-
-  tr.refine_global(5 - dim);
-
-  MappingQ<dim> mapping(2);
-  FE_Q<dim>     element(1);
-  QGauss<dim>   quadrature(4);
-
-  DoFHandler<dim> dof(tr);
-  dof.distribute_dofs(element);
-
-  FEValues<dim> fe(mapping,
-                   element,
-                   quadrature,
-                   update_values | update_gradients | update_quadrature_points |
-                     update_JxW_values);
-
-  std::vector<types::global_dof_index> global_dofs(element.dofs_per_cell);
-  std::vector<double>                  function(quadrature.size());
-
-  Vector<double> f(dof.n_dofs());
-
-  SparsityPattern A_pattern(dof.n_dofs(),
-                            dof.n_dofs(),
-                            dof.max_couplings_between_dofs());
-  DoFTools::make_sparsity_pattern(dof, A_pattern);
-  A_pattern.compress();
-
-  SparseMatrix<double> A(A_pattern);
-
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
-  const typename DoFHandler<dim>::cell_iterator  end  = dof.end();
-
-  for (; cell != end; ++cell)
-    {
-      fe.reinit(cell);
-      cell->get_dof_indices(global_dofs);
-      cosine.value_list(fe.get_quadrature_points(), function);
-
-      for (unsigned int k = 0; k < quadrature.size(); ++k)
-        {
-          double dx = fe.JxW(k);
-
-          for (unsigned int i = 0; i < element.dofs_per_cell; ++i)
-            {
-              const double         v      = fe.shape_value(i, k);
-              const Tensor<1, dim> grad_v = fe.shape_grad(i, k);
-
-              double rhs = dx * v * (function[k]);
-
-              f(global_dofs[i]) += rhs;
-              for (unsigned int j = 0; j < element.dofs_per_cell; ++j)
-                {
-                  const Tensor<1, dim> grad_u = fe.shape_grad(j, k);
-                  double               el     = dx * (grad_u * grad_v);
-                  A.add(global_dofs[i], global_dofs[j], el);
-                }
-            }
-        }
-    }
-
-  // interpolate boundary values
-  std::map<types::global_dof_index, double> bv;
-  VectorTools::interpolate_boundary_values(
-    mapping, dof, 0, cosine, bv, std::vector<bool>());
-  // the cosine has too many zero
-  // values on the boundary of the
-  // domain, so reset the elements to
-  // some other value
-  for (typename std::map<types::global_dof_index, double>::iterator i =
-         bv.begin();
-       i != bv.end();
-       ++i)
-    i->second = std::sin(i->second + 0.5) + 1.0;
-
-  // first solve filtered. this does
-  // not change the matrix
-  Vector<double> u_filtered(dof.n_dofs());
-  solve_filtered(bv, A, u_filtered, f);
-
-  // then solve by eliminating in the
-  // matrix. since this changes the
-  // matrix, this call must come
-  // second
-  Vector<double> u_eliminated(dof.n_dofs());
-  solve_eliminated<dim>(bv, A, u_eliminated, f);
-
-  // output and check
-  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
-    {
-      deallog << u_filtered(i) << std::endl;
-      Assert(std::fabs(u_filtered(i) - u_eliminated(i)) < 1e-8,
-             ExcInternalError());
-    };
-}
-
-
-int
-main()
-{
-  initlog();
-  deallog << std::setprecision(2) << std::fixed;
-
-  deallog.push("1d");
-  check<1>();
-  deallog.pop();
-  deallog.push("2d");
-  check<2>();
-  deallog.pop();
-  deallog.push("3d");
-  check<3>();
-  deallog.pop();
-}
diff --git a/tests/lac/filtered_matrix_06.output b/tests/lac/filtered_matrix_06.output
deleted file mode 100644 (file)
index 20dfc60..0000000
+++ /dev/null
@@ -1,480 +0,0 @@
-
-DEAL:1d::1.48
-DEAL:1d::1.64
-DEAL:1d::1.79
-DEAL:1d::1.94
-DEAL:1d::2.08
-DEAL:1d::2.21
-DEAL:1d::2.33
-DEAL:1d::2.43
-DEAL:1d::2.52
-DEAL:1d::2.59
-DEAL:1d::2.65
-DEAL:1d::2.69
-DEAL:1d::2.72
-DEAL:1d::2.74
-DEAL:1d::2.75
-DEAL:1d::2.75
-DEAL:1d::2.75
-DEAL:2d::1.64
-DEAL:2d::1.72
-DEAL:2d::1.68
-DEAL:2d::1.74
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.72
-DEAL:2d::1.77
-DEAL:2d::1.80
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.82
-DEAL:2d::1.83
-DEAL:2d::1.75
-DEAL:2d::1.79
-DEAL:2d::1.81
-DEAL:2d::1.78
-DEAL:2d::1.81
-DEAL:2d::1.83
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.82
-DEAL:2d::1.80
-DEAL:2d::1.72
-DEAL:2d::1.74
-DEAL:2d::1.64
-DEAL:2d::1.68
-DEAL:2d::1.77
-DEAL:2d::1.72
-DEAL:2d::1.83
-DEAL:2d::1.81
-DEAL:2d::1.84
-DEAL:2d::1.83
-DEAL:2d::1.79
-DEAL:2d::1.75
-DEAL:2d::1.81
-DEAL:2d::1.78
-DEAL:2d::1.80
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.82
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.84
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.88
-DEAL:2d::1.88
-DEAL:2d::1.85
-DEAL:2d::1.84
-DEAL:2d::1.86
-DEAL:2d::1.85
-DEAL:2d::1.83
-DEAL:2d::1.80
-DEAL:2d::1.84
-DEAL:2d::1.82
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.86
-DEAL:2d::1.84
-DEAL:2d::1.87
-DEAL:2d::1.86
-DEAL:2d::1.72
-DEAL:2d::1.74
-DEAL:2d::1.77
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.80
-DEAL:2d::1.79
-DEAL:2d::1.81
-DEAL:2d::1.81
-DEAL:2d::1.83
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.83
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.88
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.80
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.81
-DEAL:2d::1.83
-DEAL:2d::1.72
-DEAL:2d::1.74
-DEAL:2d::1.77
-DEAL:2d::1.64
-DEAL:2d::1.68
-DEAL:2d::1.72
-DEAL:2d::1.79
-DEAL:2d::1.81
-DEAL:2d::1.75
-DEAL:2d::1.78
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.87
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.80
-DEAL:2d::1.82
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.84
-DEAL:2d::1.86
-DEAL:2d::1.88
-DEAL:2d::1.88
-DEAL:2d::1.88
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.90
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.90
-DEAL:2d::1.90
-DEAL:2d::1.90
-DEAL:2d::1.90
-DEAL:2d::1.90
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.90
-DEAL:2d::1.89
-DEAL:2d::1.90
-DEAL:2d::1.90
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.89
-DEAL:2d::1.90
-DEAL:2d::1.90
-DEAL:2d::1.89
-DEAL:2d::1.90
-DEAL:2d::1.88
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.88
-DEAL:2d::1.90
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.89
-DEAL:2d::1.88
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.88
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.86
-DEAL:2d::1.72
-DEAL:2d::1.74
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.77
-DEAL:2d::1.80
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.82
-DEAL:2d::1.83
-DEAL:2d::1.79
-DEAL:2d::1.81
-DEAL:2d::1.81
-DEAL:2d::1.83
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.82
-DEAL:2d::1.80
-DEAL:2d::1.72
-DEAL:2d::1.74
-DEAL:2d::1.64
-DEAL:2d::1.68
-DEAL:2d::1.77
-DEAL:2d::1.72
-DEAL:2d::1.83
-DEAL:2d::1.81
-DEAL:2d::1.84
-DEAL:2d::1.83
-DEAL:2d::1.79
-DEAL:2d::1.75
-DEAL:2d::1.81
-DEAL:2d::1.78
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.85
-DEAL:2d::1.84
-DEAL:2d::1.86
-DEAL:2d::1.85
-DEAL:2d::1.83
-DEAL:2d::1.80
-DEAL:2d::1.84
-DEAL:2d::1.82
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.86
-DEAL:2d::1.84
-DEAL:2d::1.72
-DEAL:2d::1.74
-DEAL:2d::1.77
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.80
-DEAL:2d::1.79
-DEAL:2d::1.81
-DEAL:2d::1.81
-DEAL:2d::1.83
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.83
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.86
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.80
-DEAL:2d::1.81
-DEAL:2d::1.82
-DEAL:2d::1.77
-DEAL:2d::1.78
-DEAL:2d::1.80
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.81
-DEAL:2d::1.83
-DEAL:2d::1.72
-DEAL:2d::1.74
-DEAL:2d::1.77
-DEAL:2d::1.79
-DEAL:2d::1.81
-DEAL:2d::1.85
-DEAL:2d::1.86
-DEAL:2d::1.84
-DEAL:2d::1.85
-DEAL:2d::1.87
-DEAL:2d::1.87
-DEAL:2d::1.83
-DEAL:2d::1.84
-DEAL:2d::1.86
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.59
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.59
-DEAL:3d::1.48
-DEAL:3d::1.59
-DEAL:3d::1.63
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.59
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.59
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.59
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.55
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.53
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
-DEAL:3d::1.48
index 31571942ff36c35bc49ba1493fd9b04478de98d4..d06bb86dfc0b3bdc773e38dbac8feee29a12a66a 100644 (file)
@@ -41,7 +41,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/precondition.h>
index a171851b92c3ecafd9f2b1fe75b04f060054d4ce..4728db96afb430511bb29b7e2275c7721832ab5c 100644 (file)
@@ -39,7 +39,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/precondition.h>
index 9002e1fa69609abeef6d63b5590d26382def3245..41362d5819cf58e44d4ee60f33fd7ccd5d1585ca 100644 (file)
@@ -41,7 +41,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/precondition.h>
index 953cb01b119a4a6442ddfeb32cd453080f37b447..95dff816942da0c9aaf47b85b866feacf7022aa9 100644 (file)
@@ -42,7 +42,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/precondition.h>
index 598478ea1eea5ffa901c0c412f0f595dfbd75075..539b18259196a319db9b51e77ebe7986568e9359 100644 (file)
@@ -42,7 +42,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/petsc_vector.h>
index c5a2cd74d2737f4afc4cfa48d8b54a2116383a73..870be2c220215c583de178416e4dcc532a17595d 100644 (file)
@@ -44,7 +44,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/petsc_vector.h>
index 78dd5191f9c172c89ac40111ab6010aaf19eb2d4..5c1a6e7497354fb2cb0c501ff3915309673f221d 100644 (file)
@@ -41,7 +41,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/precondition.h>
index a844db8c0aade93eefb88e59f3d0f32fad051088..6ddcccc3eeee62c80b4c0b3a21c9475f8b1eb220 100644 (file)
@@ -42,7 +42,6 @@
 #include <deal.II/grid/tria_iterator.h>
 
 #include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/petsc_vector.h>

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.