Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
filtered_matrix.h
1 // ---------------------------------------------------------------------
2 // @f$Id: filtered_matrix.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2001 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__filtered_matrix_h
18 #define __deal2__filtered_matrix_h
19 
20 
21 
22 #include <deal.II/base/config.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/base/thread_management.h>
25 #include <deal.II/base/memory_consumption.h>
26 #include <deal.II/lac/pointer_matrix.h>
27 #include <deal.II/lac/vector_memory.h>
28 #include <vector>
29 #include <algorithm>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 template <typename number> class Vector;
34 template <class VECTOR> class FilteredMatrixBlock;
35 
36 
208 template <class VECTOR>
210 {
211 public:
213 
218 
222  class Accessor
223  {
231  const size_type index);
232 
233  public:
239  size_type row() const;
240 
246  size_type column() const;
247 
252  double value() const;
253 
254  private:
258  void advance ();
259 
264 
269  /*
270  * Make enclosing class a
271  * friend.
272  */
273  friend class const_iterator;
274  };
275 
280  {
281  public:
286  const size_type index);
287 
292 
297 
301  const Accessor &operator* () const;
302 
306  const Accessor *operator-> () const;
307 
314  bool operator == (const const_iterator &) const;
318  bool operator != (const const_iterator &) const;
319 
327  bool operator < (const const_iterator &) const;
328 
334  bool operator > (const const_iterator &) const;
335 
336  private:
342  };
343 
350  typedef std::pair<size_type, double> IndexValuePair;
351 
362  FilteredMatrix ();
363 
370  FilteredMatrix (const FilteredMatrix &fm);
371 
383  template <class MATRIX>
384  FilteredMatrix (const MATRIX &matrix,
385  bool expect_constrained_source = false);
386 
393 
409  template <class MATRIX>
410  void initialize (const MATRIX &m,
411  bool expect_constrained_source = false);
412 
417  void clear ();
419 
429  void add_constraint (const size_type i, const double v);
430 
461  template <class ConstraintList>
462  void add_constraints (const ConstraintList &new_constraints);
463 
468  void clear_constraints ();
470 
486  void apply_constraints (VECTOR &v,
487  const bool matrix_is_symmetric) const;
488 
496  void vmult (VECTOR &dst,
497  const VECTOR &src) const;
498 
507  void Tvmult (VECTOR &dst,
508  const VECTOR &src) const;
509 
522  void vmult_add (VECTOR &dst,
523  const VECTOR &src) const;
524 
537  void Tvmult_add (VECTOR &dst,
538  const VECTOR &src) const;
540 
549  const_iterator begin () const;
553  const_iterator end () const;
555 
564  std::size_t memory_consumption () const;
565 
566 private:
587 
596  typedef typename std::vector<IndexValuePair>::const_iterator const_index_value_iterator;
597 
605  {
611  bool operator () (const IndexValuePair &i1,
612  const IndexValuePair &i2) const;
613  };
614 
620  std_cxx1x::shared_ptr<PointerMatrixBase<VECTOR> > matrix;
621 
628  std::vector<IndexValuePair> constraints;
629 
636  void pre_filter (VECTOR &v) const;
637 
647  void post_filter (const VECTOR &in,
648  VECTOR &out) const;
649 
650  friend class Accessor;
655  friend class FilteredMatrixBlock<VECTOR>;
656 };
657 
659 /*---------------------- Inline functions -----------------------------------*/
660 
661 
662 //--------------------------------Iterators--------------------------------------//
663 
664 template<class VECTOR>
665 inline
667  const FilteredMatrix<VECTOR> *matrix,
668  const size_type index)
669  :
670  matrix(matrix),
671  index(index)
672 {
673  Assert (index <= matrix->constraints.size(),
674  ExcIndexRange(index, 0, matrix->constraints.size()));
675 }
676 
677 
678 
679 template<class VECTOR>
680 inline
683 {
684  return matrix->constraints[index].first;
685 }
686 
687 
688 
689 template<class VECTOR>
690 inline
693 {
694  return matrix->constraints[index].first;
695 }
696 
697 
698 
699 template<class VECTOR>
700 inline
701 double
703 {
704  return matrix->constraints[index].second;
705 }
706 
707 
708 
709 template<class VECTOR>
710 inline
711 void
713 {
714  Assert (index < matrix->constraints.size(), ExcIteratorPastEnd());
715  ++index;
716 }
717 
718 
719 
720 
721 template<class VECTOR>
722 inline
725  const size_type index)
726  :
727  accessor(matrix, index)
728 {}
729 
730 
731 
732 template<class VECTOR>
733 inline
736 {
737  accessor.advance();
738  return *this;
739 }
740 
741 
742 template <typename number>
743 inline
744 const typename FilteredMatrix<number>::Accessor &
746 {
747  return accessor;
748 }
749 
750 
751 template <typename number>
752 inline
753 const typename FilteredMatrix<number>::Accessor *
755 {
756  return &accessor;
757 }
758 
759 
760 template <typename number>
761 inline
762 bool
764 operator == (const const_iterator &other) const
765 {
766  return (accessor.index == other.accessor.index
767  && accessor.matrix == other.accessor.matrix);
768 }
769 
770 
771 template <typename number>
772 inline
773 bool
775 operator != (const const_iterator &other) const
776 {
777  return ! (*this == other);
778 }
779 
780 
781 
782 //------------------------------- FilteredMatrix ---------------------------------------//
783 
784 template <typename number>
785 inline
788 {
789  return const_iterator(this, 0);
790 }
791 
792 
793 template <typename number>
794 inline
797 {
798  return const_iterator(this, constraints.size());
799 }
800 
801 
802 template <class VECTOR>
803 inline
804 bool
807  const IndexValuePair &i2) const
808 {
809  return (i1.first < i2.first);
810 }
811 
812 
813 
814 template <class VECTOR>
815 template <class MATRIX>
816 inline
817 void
819 {
820  matrix.reset (new_pointer_matrix_base(m, VECTOR()));
821 
823 }
824 
825 
826 
827 template <class VECTOR>
828 inline
830 {}
831 
832 
833 
834 template <class VECTOR>
835 inline
837  :
838  Subscriptor(),
840  matrix(fm.matrix),
842 {}
843 
844 
845 
846 template <class VECTOR>
847 template <class MATRIX>
848 inline
850 FilteredMatrix (const MATRIX &m, bool ecs)
851 {
852  initialize (m, ecs);
853 }
854 
855 
856 
857 template <class VECTOR>
858 inline
861 {
862  matrix = fm.matrix;
863  expect_constrained_source = fm.expect_constrained_source;
864  constraints = fm.constraints;
865  return *this;
866 }
867 
868 
869 
870 template <class VECTOR>
871 inline
872 void
873 FilteredMatrix<VECTOR>::add_constraint (const size_type index, const double value)
874 {
875  // add new constraint to end
876  constraints.push_back(IndexValuePair(index, value));
877 }
878 
879 
880 
881 template <class VECTOR>
882 template <class ConstraintList>
883 inline
884 void
885 FilteredMatrix<VECTOR>::add_constraints (const ConstraintList &new_constraints)
886 {
887  // add new constraints to end
888  const size_type old_size = constraints.size();
889  constraints.reserve (old_size + new_constraints.size());
890  constraints.insert (constraints.end(),
891  new_constraints.begin(),
892  new_constraints.end());
893  // then merge the two arrays to
894  // form one sorted one
895  std::inplace_merge (constraints.begin(),
896  constraints.begin()+old_size,
897  constraints.end(),
898  PairComparison());
899 }
900 
901 
902 
903 template <class VECTOR>
904 inline
905 void
907 {
908  // swap vectors to release memory
909  std::vector<IndexValuePair> empty;
910  constraints.swap (empty);
911 }
912 
913 
914 
915 template <class VECTOR>
916 inline
917 void
919 {
920  clear_constraints();
921  matrix.reset();
922 }
923 
924 
925 
926 template <class VECTOR>
927 inline
928 void
930  VECTOR &v,
931  const bool /* matrix_is_symmetric */) const
932 {
934  typename VectorMemory<VECTOR>::Pointer tmp_vector(mem);
935  tmp_vector->reinit(v);
936  const_index_value_iterator i = constraints.begin();
937  const const_index_value_iterator e = constraints.end();
938  for (; i!=e; ++i)
939  {
941  (*tmp_vector)(i->first) = -i->second;
942  }
943 
944  // This vmult is without bc, to get
945  // the rhs correction in a correct
946  // way.
947  matrix->vmult_add(v, *tmp_vector);
948  // finally set constrained
949  // entries themselves
950  for (i=constraints.begin(); i!=e; ++i)
951  {
953  v(i->first) = i->second;
954  }
955 }
956 
957 
958 
959 template <class VECTOR>
960 inline
961 void
963 {
964  // iterate over all constraints and
965  // zero out value
966  const_index_value_iterator i = constraints.begin();
967  const const_index_value_iterator e = constraints.end();
968  for (; i!=e; ++i)
969  v(i->first) = 0;
970 }
971 
972 
973 
974 template <class VECTOR>
975 inline
976 void
978  VECTOR &out) const
979 {
980  // iterate over all constraints and
981  // set value correctly
982  const_index_value_iterator i = constraints.begin();
983  const const_index_value_iterator e = constraints.end();
984  for (; i!=e; ++i)
985  {
986  Assert(numbers::is_finite(in(i->first)), ExcNumberNotFinite());
987  out(i->first) = in(i->first);
988  }
989 }
990 
991 
992 
993 template <class VECTOR>
994 inline
995 void
996 FilteredMatrix<VECTOR>::vmult (VECTOR &dst, const VECTOR &src) const
997 {
998  if (!expect_constrained_source)
999  {
1001  VECTOR *tmp_vector = mem.alloc();
1002  // first copy over src vector and
1003  // pre-filter
1004  tmp_vector->reinit(src, true);
1005  *tmp_vector = src;
1006  pre_filter (*tmp_vector);
1007  // then let matrix do its work
1008  matrix->vmult (dst, *tmp_vector);
1009  mem.free(tmp_vector);
1010  }
1011  else
1012  {
1013  matrix->vmult (dst, src);
1014  }
1015 
1016  // finally do post-filtering
1017  post_filter (src, dst);
1018 }
1019 
1020 
1021 
1022 template <class VECTOR>
1023 inline
1024 void
1025 FilteredMatrix<VECTOR>::Tvmult (VECTOR &dst, const VECTOR &src) const
1026 {
1027  if (!expect_constrained_source)
1028  {
1030  VECTOR *tmp_vector = mem.alloc();
1031  // first copy over src vector and
1032  // pre-filter
1033  tmp_vector->reinit(src, true);
1034  *tmp_vector = src;
1035  pre_filter (*tmp_vector);
1036  // then let matrix do its work
1037  matrix->Tvmult (dst, *tmp_vector);
1038  mem.free(tmp_vector);
1039  }
1040  else
1041  {
1042  matrix->Tvmult (dst, src);
1043  }
1044 
1045  // finally do post-filtering
1046  post_filter (src, dst);
1047 }
1048 
1049 
1050 
1051 template <class VECTOR>
1052 inline
1053 void
1054 FilteredMatrix<VECTOR>::vmult_add (VECTOR &dst, const VECTOR &src) const
1055 {
1056  if (!expect_constrained_source)
1057  {
1059  VECTOR *tmp_vector = mem.alloc();
1060  // first copy over src vector and
1061  // pre-filter
1062  tmp_vector->reinit(src, true);
1063  *tmp_vector = src;
1064  pre_filter (*tmp_vector);
1065  // then let matrix do its work
1066  matrix->vmult_add (dst, *tmp_vector);
1067  mem.free(tmp_vector);
1068  }
1069  else
1070  {
1071  matrix->vmult_add (dst, src);
1072  }
1073 
1074  // finally do post-filtering
1075  post_filter (src, dst);
1076 }
1077 
1078 
1079 
1080 template <class VECTOR>
1081 inline
1082 void
1083 FilteredMatrix<VECTOR>::Tvmult_add (VECTOR &dst, const VECTOR &src) const
1084 {
1085  if (!expect_constrained_source)
1086  {
1088  VECTOR *tmp_vector = mem.alloc();
1089  // first copy over src vector and
1090  // pre-filter
1091  tmp_vector->reinit(src, true);
1092  *tmp_vector = src;
1093  pre_filter (*tmp_vector);
1094  // then let matrix do its work
1095  matrix->Tvmult_add (dst, *tmp_vector);
1096  mem.free(tmp_vector);
1097  }
1098  else
1099  {
1100  matrix->Tvmult_add (dst, src);
1101  }
1102 
1103  // finally do post-filtering
1104  post_filter (src, dst);
1105 }
1106 
1107 
1108 
1109 template <class VECTOR>
1110 inline
1111 std::size_t
1113 {
1114  return (MemoryConsumption::memory_consumption (matrix) +
1116 }
1117 
1118 
1119 
1120 DEAL_II_NAMESPACE_CLOSE
1121 
1122 #endif
1123 /*---------------------------- filtered_matrix.h ---------------------------*/
void vmult(VECTOR &dst, const VECTOR &src) const
const_iterator(const FilteredMatrix< VECTOR > *matrix, const size_type index)
std::size_t memory_consumption() const
void add_constraints(const ConstraintList &new_constraints)
void clear_constraints()
bool operator>(const const_iterator &) const
void initialize(const MATRIX &m, bool expect_constrained_source=false)
virtual VECTOR * alloc()
std::vector< IndexValuePair > constraints
bool is_finite(const double x)
bool expect_constrained_source
FilteredMatrix & operator=(const FilteredMatrix &fm)
const_iterator end() const
unsigned int global_dof_index
Definition: types.h:100
size_type column() const
void Tvmult(VECTOR &dst, const VECTOR &src) const
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption(const T &t)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
bool operator()(const IndexValuePair &i1, const IndexValuePair &i2) const
void post_filter(const VECTOR &in, VECTOR &out) const
void apply_constraints(VECTOR &v, const bool matrix_is_symmetric) const
std::vector< IndexValuePair >::const_iterator const_index_value_iterator
const FilteredMatrix< VECTOR > * matrix
types::global_dof_index size_type
::ExceptionBase & ExcIteratorPastEnd()
std_cxx1x::shared_ptr< PointerMatrixBase< VECTOR > > matrix
::ExceptionBase & ExcNumberNotFinite()
bool operator<(const const_iterator &) const
void vmult_add(VECTOR &dst, const VECTOR &src) const
const_iterator begin() const
void pre_filter(VECTOR &v) const
void add_constraint(const size_type i, const double v)
virtual void free(const VECTOR *const)
bool operator==(const const_iterator &) const
void Tvmult_add(VECTOR &dst, const VECTOR &src) const
bool operator!=(const const_iterator &) const
const Accessor * operator->() const
size_type row() const
Accessor(const FilteredMatrix< VECTOR > *matrix, const size_type index)
std::pair< size_type, double > IndexValuePair
const Accessor & operator*() const