Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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  Assert(index <= matrix->constraints.size(),
585  ExcIndexRange(index, 0, matrix->constraints.size()));
586 }
587 
588 
589 
590 template <typename VectorType>
593 {
594  return matrix->constraints[index].first;
595 }
596 
597 
598 
599 template <typename VectorType>
602 {
603  return matrix->constraints[index].first;
604 }
605 
606 
607 
608 template <typename VectorType>
609 inline double
611 {
612  return matrix->constraints[index].second;
613 }
614 
615 
616 
617 template <typename VectorType>
618 inline void
620 {
621  Assert(index < matrix->constraints.size(), ExcIteratorPastEnd());
622  ++index;
623 }
624 
625 
626 
627 template <typename VectorType>
630  const size_type index)
631  : accessor(matrix, index)
632 {}
633 
634 
635 
636 template <typename VectorType>
639 {
640  accessor.advance();
641  return *this;
642 }
643 
644 
645 template <typename number>
646 inline const typename FilteredMatrix<number>::Accessor &
648 {
649  return accessor;
650 }
651 
652 
653 template <typename number>
654 inline const typename FilteredMatrix<number>::Accessor *
656 {
657  return &accessor;
658 }
659 
660 
661 template <typename number>
662 inline bool
664 operator==(const const_iterator &other) const
665 {
666  return (accessor.index == other.accessor.index &&
667  accessor.matrix == other.accessor.matrix);
668 }
669 
670 
671 template <typename number>
672 inline bool
674 operator!=(const const_iterator &other) const
675 {
676  return !(*this == other);
677 }
678 
679 
680 
681 //------------------------- FilteredMatrix ----------------------------------//
682 
683 template <typename number>
686 {
687  return const_iterator(this, 0);
688 }
689 
690 
691 template <typename number>
694 {
695  return const_iterator(this, constraints.size());
696 }
697 
698 
699 template <typename VectorType>
700 inline bool
702 operator()(const IndexValuePair &i1, const IndexValuePair &i2) const
703 {
704  return (i1.first < i2.first);
705 }
706 
707 
708 
709 template <typename VectorType>
710 template <typename MatrixType>
711 inline void
712 FilteredMatrix<VectorType>::initialize(const MatrixType &m, bool ecs)
713 {
714  matrix.reset(new_pointer_matrix_base(m, VectorType()));
715 
717 }
718 
719 
720 
721 template <typename VectorType>
724 {}
725 
726 
727 
728 template <typename VectorType>
730  : Subscriptor()
732  , matrix(fm.matrix)
734 {}
735 
736 
737 
738 template <typename VectorType>
739 template <typename MatrixType>
741  const bool ecs)
743 {
744  initialize(m, ecs);
745 }
746 
747 
748 
749 template <typename VectorType>
752 {
753  matrix = fm.matrix;
756  return *this;
757 }
758 
759 
760 
761 template <typename VectorType>
762 inline void
764  const double value)
765 {
766  // add new constraint to end
767  constraints.push_back(IndexValuePair(index, value));
768 }
769 
770 
771 
772 template <typename VectorType>
773 template <class ConstraintList>
774 inline void
776  const ConstraintList &new_constraints)
777 {
778  // add new constraints to end
779  const size_type old_size = constraints.size();
780  constraints.reserve(old_size + new_constraints.size());
781  constraints.insert(constraints.end(),
782  new_constraints.begin(),
783  new_constraints.end());
784  // then merge the two arrays to
785  // form one sorted one
786  std::inplace_merge(constraints.begin(),
787  constraints.begin() + old_size,
788  constraints.end(),
789  PairComparison());
790 }
791 
792 
793 
794 template <typename VectorType>
795 inline void
797 {
798  // swap vectors to release memory
799  std::vector<IndexValuePair> empty;
800  constraints.swap(empty);
801 }
802 
803 
804 
805 template <typename VectorType>
806 inline void
808 {
810  matrix.reset();
811 }
812 
813 
814 
815 template <typename VectorType>
816 inline void
818  VectorType &v,
819  const bool /* matrix_is_symmetric */) const
820 {
822 }
823 
824 
825 template <typename VectorType>
826 inline void
828 {
830  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
831  tmp_vector->reinit(v);
833  const const_index_value_iterator e = constraints.end();
834  for (; i != e; ++i)
835  {
836  AssertIsFinite(i->second);
837  (*tmp_vector)(i->first) = -i->second;
838  }
839 
840  // This vmult is without bc, to get
841  // the rhs correction in a correct
842  // way.
843  matrix->vmult_add(v, *tmp_vector);
844  // finally set constrained
845  // entries themselves
846  for (i = constraints.begin(); i != e; ++i)
847  {
848  AssertIsFinite(i->second);
849  v(i->first) = i->second;
850  }
851 }
852 
853 
854 template <typename VectorType>
855 inline void
857 {
858  // iterate over all constraints and
859  // zero out value
861  const const_index_value_iterator e = constraints.end();
862  for (; i != e; ++i)
863  v(i->first) = 0;
864 }
865 
866 
867 
868 template <typename VectorType>
869 inline void
871  VectorType & out) const
872 {
873  // iterate over all constraints and
874  // set value correctly
876  const const_index_value_iterator e = constraints.end();
877  for (; i != e; ++i)
878  {
879  AssertIsFinite(in(i->first));
880  out(i->first) = in(i->first);
881  }
882 }
883 
884 
885 
886 template <typename VectorType>
887 inline void
888 FilteredMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
889 {
891  {
893  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
894  // first copy over src vector and
895  // pre-filter
896  tmp_vector->reinit(src, true);
897  *tmp_vector = src;
898  pre_filter(*tmp_vector);
899  // then let matrix do its work
900  matrix->vmult(dst, *tmp_vector);
901  }
902  else
903  {
904  matrix->vmult(dst, src);
905  }
906 
907  // finally do post-filtering
908  post_filter(src, dst);
909 }
910 
911 
912 
913 template <typename VectorType>
914 inline void
915 FilteredMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
916 {
918  {
920  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
921  // first copy over src vector and
922  // pre-filter
923  tmp_vector->reinit(src, true);
924  *tmp_vector = src;
925  pre_filter(*tmp_vector);
926  // then let matrix do its work
927  matrix->Tvmult(dst, *tmp_vector);
928  }
929  else
930  {
931  matrix->Tvmult(dst, src);
932  }
933 
934  // finally do post-filtering
935  post_filter(src, dst);
936 }
937 
938 
939 
940 template <typename VectorType>
941 inline void
943  const VectorType &src) const
944 {
946  {
948  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
949  // first copy over src vector and
950  // pre-filter
951  tmp_vector->reinit(src, true);
952  *tmp_vector = src;
953  pre_filter(*tmp_vector);
954  // then let matrix do its work
955  matrix->vmult_add(dst, *tmp_vector);
956  }
957  else
958  {
959  matrix->vmult_add(dst, src);
960  }
961 
962  // finally do post-filtering
963  post_filter(src, dst);
964 }
965 
966 
967 
968 template <typename VectorType>
969 inline void
971  const VectorType &src) const
972 {
974  {
976  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
977  // first copy over src vector and
978  // pre-filter
979  tmp_vector->reinit(src, true);
980  *tmp_vector = src;
981  pre_filter(*tmp_vector);
982  // then let matrix do its work
983  matrix->Tvmult_add(dst, *tmp_vector);
984  }
985  else
986  {
987  matrix->Tvmult_add(dst, src);
988  }
989 
990  // finally do post-filtering
991  post_filter(src, dst);
992 }
993 
994 
995 
996 template <typename VectorType>
997 inline std::size_t
999 {
1002 }
1003 
1004 
1005 
1006 DEAL_II_NAMESPACE_CLOSE
1007 
1008 #endif
1009 /*---------------------------- 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
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
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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:1407
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:1669
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