--- /dev/null
+Removed: The deprecated FilteredMatrix class has been removed.
+Use the LinearOperator class instead.
+<br>
+(Daniel Arndt, 2020/04/01)
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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 ---------------------------*/
* 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.
*
*
#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>
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,
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
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]));
#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>
* 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;
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
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();
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)
// 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
#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;
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)
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())
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];
%%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
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
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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);
-}
+++ /dev/null
-
-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
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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);
-}
+++ /dev/null
-
-DEAL::vmult
-DEAL::1.0000 28.0000 43.0000
-DEAL::vmult
-DEAL::1.0000 32.0000 50.0000
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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);
-}
+++ /dev/null
-
-DEAL::Tvmult
-DEAL::1.0000 34.0000 39.0000
-DEAL::Tvmult
-DEAL::1.0000 36.0000 42.0000
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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);
-}
+++ /dev/null
-
-DEAL::vmult_add
-DEAL::1.0000 29.0000 44.0000
-DEAL::vmult_add
-DEAL::1.0000 33.0000 51.0000
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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);
-}
+++ /dev/null
-
-DEAL::Tvmult_add
-DEAL::1.0000 35.0000 40.0000
-DEAL::Tvmult_add
-DEAL::1.0000 37.0000 43.0000
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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);
-}
+++ /dev/null
-
-DEAL::apply_constraints
-DEAL::1.0000 -2.0000 -4.0000
-DEAL::apply_constraints
-DEAL::1.0000 -2.0000 -4.0000
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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();
-}
+++ /dev/null
-
-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
#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>
#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>
#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>
#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>
#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>
#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>
#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>
#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>
projects / dealii.git / commitdiff