Reference documentation for deal.II version Git b84270a1d4 2020-01-17 16:21:02 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
filtered_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_filtered_matrix_h
17 # define dealii_filtered_matrix_h
18 
19 
20 
21 # include <deal.II/base/config.h>
22 
23 # include <deal.II/base/memory_consumption.h>
24 # include <deal.II/base/smartpointer.h>
25 # include <deal.II/base/thread_management.h>
26 
27 # include <deal.II/lac/pointer_matrix.h>
28 # include <deal.II/lac/vector_memory.h>
29 
30 # include <algorithm>
31 # include <vector>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 // Forward declaration
36 # ifndef DOXYGEN
37 template <class VectorType>
38 class FilteredMatrixBlock;
39 # endif
40 
201 template <typename VectorType>
202 class DEAL_II_DEPRECATED FilteredMatrix : public Subscriptor
203 {
204 public:
205  class const_iterator;
206 
211 
215  class Accessor
216  {
221  Accessor(const FilteredMatrix<VectorType> *matrix, const size_type index);
222 
223  public:
227  size_type
228  row() const;
229 
233  size_type
234  column() const;
235 
239  double
240  value() const;
241 
242  private:
246  void
247  advance();
248 
253 
258 
259  // Make enclosing class a friend.
260  friend class const_iterator;
261  };
262 
267  {
268  public:
273  const size_type index);
274 
279  operator++();
280 
285  operator++(int);
286 
290  const Accessor &operator*() const;
291 
295  const Accessor *operator->() const;
296 
300  bool
301  operator==(const const_iterator &) const;
305  bool
306  operator!=(const const_iterator &) const;
307 
312  bool
313  operator<(const const_iterator &) const;
314 
319  bool
320  operator>(const const_iterator &) const;
321 
322  private:
327  };
328 
333  using IndexValuePair = std::pair<size_type, double>;
334 
343  FilteredMatrix();
344 
349  FilteredMatrix(const FilteredMatrix &fm);
350 
359  template <typename MatrixType>
360  FilteredMatrix(const MatrixType &matrix,
361  const bool expect_constrained_source = false);
362 
367  operator=(const FilteredMatrix &fm);
368 
378  template <typename MatrixType>
379  void
380  initialize(const MatrixType &m, const bool expect_constrained_source = false);
381 
385  void
386  clear();
388 
396  void
397  add_constraint(const size_type i, const double v);
398 
414  template <class ConstraintList>
415  void
416  add_constraints(const ConstraintList &new_constraints);
417 
421  void
422  clear_constraints();
424 
435  DEAL_II_DEPRECATED
436  void
437  apply_constraints(VectorType &v, const bool matrix_is_symmetric) const;
442  void
443  apply_constraints(VectorType &v) const;
444 
449  void
450  vmult(VectorType &dst, const VectorType &src) const;
451 
457  void
458  Tvmult(VectorType &dst, const VectorType &src) const;
459 
467  void
468  vmult_add(VectorType &dst, const VectorType &src) const;
469 
477  void
478  Tvmult_add(VectorType &dst, const VectorType &src) const;
480 
489  begin() const;
494  end() const;
496 
502  std::size_t
503  memory_consumption() const;
504 
505 private:
518 
525  typename std::vector<IndexValuePair>::const_iterator;
526 
532  {
536  bool
537  operator()(const IndexValuePair &i1, const IndexValuePair &i2) const;
538  };
539 
543  std::shared_ptr<PointerMatrixBase<VectorType>> matrix;
544 
549  std::vector<IndexValuePair> constraints;
550 
555  void
556  pre_filter(VectorType &v) const;
557 
563  void
564  post_filter(const VectorType &in, VectorType &out) const;
565 
566  friend class Accessor;
567  // FilteredMatrixBlock accesses pre_filter() and post_filter().
568  friend class FilteredMatrixBlock<VectorType>;
569 };
570 
572 /*---------------------- Inline functions -----------------------------------*/
573 
574 
575 //--------------------------------Iterators--------------------------------------//
576 
577 template <typename VectorType>
579  const FilteredMatrix<VectorType> *matrix,
580  const size_type index)
581  : matrix(matrix)
582  , index(index)
583 {
584  AssertIndexRange(index, matrix->constraints.size() + 1);
585 }
586 
587 
588 
589 template <typename VectorType>
592 {
593  return matrix->constraints[index].first;
594 }
595 
596 
597 
598 template <typename VectorType>
601 {
602  return matrix->constraints[index].first;
603 }
604 
605 
606 
607 template <typename VectorType>
608 inline double
610 {
611  return matrix->constraints[index].second;
612 }
613 
614 
615 
616 template <typename VectorType>
617 inline void
619 {
620  Assert(index < matrix->constraints.size(), ExcIteratorPastEnd());
621  ++index;
622 }
623 
624 
625 
626 template <typename VectorType>
629  const size_type index)
630  : accessor(matrix, index)
631 {}
632 
633 
634 
635 template <typename VectorType>
638 {
639  accessor.advance();
640  return *this;
641 }
642 
643 
644 template <typename number>
645 inline const typename FilteredMatrix<number>::Accessor &
647 {
648  return accessor;
649 }
650 
651 
652 template <typename number>
653 inline const typename FilteredMatrix<number>::Accessor *
655 {
656  return &accessor;
657 }
658 
659 
660 template <typename number>
661 inline bool
663 operator==(const const_iterator &other) const
664 {
665  return (accessor.index == other.accessor.index &&
666  accessor.matrix == other.accessor.matrix);
667 }
668 
669 
670 template <typename number>
671 inline bool
673 operator!=(const const_iterator &other) const
674 {
675  return !(*this == other);
676 }
677 
678 
679 
680 //------------------------- FilteredMatrix ----------------------------------//
681 
682 template <typename number>
685 {
686  return const_iterator(this, 0);
687 }
688 
689 
690 template <typename number>
693 {
694  return const_iterator(this, constraints.size());
695 }
696 
697 
698 template <typename VectorType>
699 inline bool
701 operator()(const IndexValuePair &i1, const IndexValuePair &i2) const
702 {
703  return (i1.first < i2.first);
704 }
705 
706 
707 
708 template <typename VectorType>
709 template <typename MatrixType>
710 inline void
711 FilteredMatrix<VectorType>::initialize(const MatrixType &m, bool ecs)
712 {
713  matrix.reset(new_pointer_matrix_base(m, VectorType()));
714 
716 }
717 
718 
719 
720 template <typename VectorType>
723 {}
724 
725 
726 
727 template <typename VectorType>
729  : Subscriptor()
731  , matrix(fm.matrix)
733 {}
734 
735 
736 
737 template <typename VectorType>
738 template <typename MatrixType>
740  const bool ecs)
742 {
743  initialize(m, ecs);
744 }
745 
746 
747 
748 template <typename VectorType>
751 {
752  matrix = fm.matrix;
755  return *this;
756 }
757 
758 
759 
760 template <typename VectorType>
761 inline void
763  const double value)
764 {
765  // add new constraint to end
766  constraints.push_back(IndexValuePair(index, value));
767 }
768 
769 
770 
771 template <typename VectorType>
772 template <class ConstraintList>
773 inline void
775  const ConstraintList &new_constraints)
776 {
777  // add new constraints to end
778  const size_type old_size = constraints.size();
779  constraints.reserve(old_size + new_constraints.size());
780  constraints.insert(constraints.end(),
781  new_constraints.begin(),
782  new_constraints.end());
783  // then merge the two arrays to
784  // form one sorted one
785  std::inplace_merge(constraints.begin(),
786  constraints.begin() + old_size,
787  constraints.end(),
788  PairComparison());
789 }
790 
791 
792 
793 template <typename VectorType>
794 inline void
796 {
797  // swap vectors to release memory
798  std::vector<IndexValuePair> empty;
799  constraints.swap(empty);
800 }
801 
802 
803 
804 template <typename VectorType>
805 inline void
807 {
809  matrix.reset();
810 }
811 
812 
813 
814 template <typename VectorType>
815 inline void
817  VectorType &v,
818  const bool /* matrix_is_symmetric */) const
819 {
821 }
822 
823 
824 template <typename VectorType>
825 inline void
827 {
829  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
830  tmp_vector->reinit(v);
832  const const_index_value_iterator e = constraints.end();
833  for (; i != e; ++i)
834  {
835  AssertIsFinite(i->second);
836  (*tmp_vector)(i->first) = -i->second;
837  }
838 
839  // This vmult is without bc, to get
840  // the rhs correction in a correct
841  // way.
842  matrix->vmult_add(v, *tmp_vector);
843  // finally set constrained
844  // entries themselves
845  for (i = constraints.begin(); i != e; ++i)
846  {
847  AssertIsFinite(i->second);
848  v(i->first) = i->second;
849  }
850 }
851 
852 
853 template <typename VectorType>
854 inline void
856 {
857  // iterate over all constraints and
858  // zero out value
860  const const_index_value_iterator e = constraints.end();
861  for (; i != e; ++i)
862  v(i->first) = 0;
863 }
864 
865 
866 
867 template <typename VectorType>
868 inline void
870  VectorType & out) const
871 {
872  // iterate over all constraints and
873  // set value correctly
875  const const_index_value_iterator e = constraints.end();
876  for (; i != e; ++i)
877  {
878  AssertIsFinite(in(i->first));
879  out(i->first) = in(i->first);
880  }
881 }
882 
883 
884 
885 template <typename VectorType>
886 inline void
887 FilteredMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
888 {
890  {
892  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
893  // first copy over src vector and
894  // pre-filter
895  tmp_vector->reinit(src, true);
896  *tmp_vector = src;
897  pre_filter(*tmp_vector);
898  // then let matrix do its work
899  matrix->vmult(dst, *tmp_vector);
900  }
901  else
902  {
903  matrix->vmult(dst, src);
904  }
905 
906  // finally do post-filtering
907  post_filter(src, dst);
908 }
909 
910 
911 
912 template <typename VectorType>
913 inline void
914 FilteredMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
915 {
917  {
919  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
920  // first copy over src vector and
921  // pre-filter
922  tmp_vector->reinit(src, true);
923  *tmp_vector = src;
924  pre_filter(*tmp_vector);
925  // then let matrix do its work
926  matrix->Tvmult(dst, *tmp_vector);
927  }
928  else
929  {
930  matrix->Tvmult(dst, src);
931  }
932 
933  // finally do post-filtering
934  post_filter(src, dst);
935 }
936 
937 
938 
939 template <typename VectorType>
940 inline void
942  const VectorType &src) const
943 {
945  {
947  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
948  // first copy over src vector and
949  // pre-filter
950  tmp_vector->reinit(src, true);
951  *tmp_vector = src;
952  pre_filter(*tmp_vector);
953  // then let matrix do its work
954  matrix->vmult_add(dst, *tmp_vector);
955  }
956  else
957  {
958  matrix->vmult_add(dst, src);
959  }
960 
961  // finally do post-filtering
962  post_filter(src, dst);
963 }
964 
965 
966 
967 template <typename VectorType>
968 inline void
970  const VectorType &src) const
971 {
973  {
975  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
976  // first copy over src vector and
977  // pre-filter
978  tmp_vector->reinit(src, true);
979  *tmp_vector = src;
980  pre_filter(*tmp_vector);
981  // then let matrix do its work
982  matrix->Tvmult_add(dst, *tmp_vector);
983  }
984  else
985  {
986  matrix->Tvmult_add(dst, src);
987  }
988 
989  // finally do post-filtering
990  post_filter(src, dst);
991 }
992 
993 
994 
995 template <typename VectorType>
996 inline std::size_t
998 {
1001 }
1002 
1003 
1004 
1005 DEAL_II_NAMESPACE_CLOSE
1006 
1007 #endif
1008 /*---------------------------- filtered_matrix.h ---------------------------*/
void initialize(const MatrixType &m, const bool expect_constrained_source=false)
const FilteredMatrix< VectorType > * matrix
const Accessor * operator->() const
const_iterator begin() const
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
void vmult_add(VectorType &dst, const VectorType &src) const
size_type column() const
void pre_filter(VectorType &v) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
std::shared_ptr< PointerMatrixBase< VectorType > > matrix
void apply_constraints(VectorType &v, const bool matrix_is_symmetric) const
void post_filter(const VectorType &in, VectorType &out) const
std::vector< IndexValuePair > constraints
FilteredMatrix & operator=(const FilteredMatrix &fm)
Accessor(const FilteredMatrix< VectorType > *matrix, const size_type index)
void add_constraint(const size_type i, const double v)
PointerMatrixBase< VectorType > * new_pointer_matrix_base(MatrixType &matrix, const VectorType &, const char *name="PointerMatrixAux")
const Accessor & operator*() const
void vmult(VectorType &dst, const VectorType &src) const
void Tvmult(VectorType &dst, const VectorType &src) const
#define Assert(cond, exc)
Definition: exceptions.h:1411
std::pair< size_type, double > IndexValuePair
bool expect_constrained_source
bool operator()(const IndexValuePair &i1, const IndexValuePair &i2) const
bool operator!=(const const_iterator &) const
bool operator==(const const_iterator &) const
void Tvmult_add(VectorType &dst, const VectorType &src) const
typename std::vector< IndexValuePair >::const_iterator const_index_value_iterator
void add_constraints(const ConstraintList &new_constraints)
void clear_constraints()
static ::ExceptionBase & ExcIteratorPastEnd()
unsigned int global_dof_index
Definition: types.h:89
const_iterator(const FilteredMatrix< VectorType > *matrix, const size_type index)
#define AssertIsFinite(number)
Definition: exceptions.h:1673
std::size_t memory_consumption() const
types::global_dof_index size_type
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
const_iterator end() const