Reference documentation for deal.II version 9.1.1
\(\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 template <class VectorType>
36 class FilteredMatrixBlock;
37 
198 template <typename VectorType>
199 class DEAL_II_DEPRECATED FilteredMatrix : public Subscriptor
200 {
201 public:
202  class const_iterator;
203 
208 
212  class Accessor
213  {
218  Accessor(const FilteredMatrix<VectorType> *matrix, const size_type index);
219 
220  public:
224  size_type
225  row() const;
226 
230  size_type
231  column() const;
232 
236  double
237  value() const;
238 
239  private:
243  void
244  advance();
245 
250 
255  /*
256  * Make enclosing class a
257  * friend.
258  */
259  friend class const_iterator;
260  };
261 
266  {
267  public:
272  const size_type index);
273 
278  operator++();
279 
284  operator++(int);
285 
289  const Accessor &operator*() const;
290 
294  const Accessor *operator->() const;
295 
299  bool
300  operator==(const const_iterator &) const;
304  bool
305  operator!=(const const_iterator &) const;
306 
311  bool
312  operator<(const const_iterator &) const;
313 
318  bool
319  operator>(const const_iterator &) const;
320 
321  private:
326  };
327 
332  using IndexValuePair = std::pair<size_type, double>;
333 
342  FilteredMatrix();
343 
348  FilteredMatrix(const FilteredMatrix &fm);
349 
358  template <typename MatrixType>
359  FilteredMatrix(const MatrixType &matrix,
360  const bool expect_constrained_source = false);
361 
366  operator=(const FilteredMatrix &fm);
367 
377  template <typename MatrixType>
378  void
379  initialize(const MatrixType &m, const bool expect_constrained_source = false);
380 
384  void
385  clear();
387 
395  void
396  add_constraint(const size_type i, const double v);
397 
413  template <class ConstraintList>
414  void
415  add_constraints(const ConstraintList &new_constraints);
416 
420  void
421  clear_constraints();
423 
434  DEAL_II_DEPRECATED
435  void
436  apply_constraints(VectorType &v, const bool matrix_is_symmetric) const;
441  void
442  apply_constraints(VectorType &v) const;
443 
448  void
449  vmult(VectorType &dst, const VectorType &src) const;
450 
456  void
457  Tvmult(VectorType &dst, const VectorType &src) const;
458 
466  void
467  vmult_add(VectorType &dst, const VectorType &src) const;
468 
476  void
477  Tvmult_add(VectorType &dst, const VectorType &src) const;
479 
488  begin() const;
493  end() const;
495 
501  std::size_t
502  memory_consumption() const;
503 
504 private:
517 
524  typename std::vector<IndexValuePair>::const_iterator;
525 
531  {
535  bool
536  operator()(const IndexValuePair &i1, const IndexValuePair &i2) const;
537  };
538 
542  std::shared_ptr<PointerMatrixBase<VectorType>> matrix;
543 
548  std::vector<IndexValuePair> constraints;
549 
554  void
555  pre_filter(VectorType &v) const;
556 
562  void
563  post_filter(const VectorType &in, VectorType &out) const;
564 
565  friend class Accessor;
569  friend class FilteredMatrixBlock<VectorType>;
570 };
571 
573 /*---------------------- Inline functions -----------------------------------*/
574 
575 
576 //--------------------------------Iterators--------------------------------------//
577 
578 template <typename VectorType>
580  const FilteredMatrix<VectorType> *matrix,
581  const size_type index)
582  : matrix(matrix)
583  , index(index)
584 {
585  Assert(index <= matrix->constraints.size(),
586  ExcIndexRange(index, 0, matrix->constraints.size()));
587 }
588 
589 
590 
591 template <typename VectorType>
594 {
595  return matrix->constraints[index].first;
596 }
597 
598 
599 
600 template <typename VectorType>
603 {
604  return matrix->constraints[index].first;
605 }
606 
607 
608 
609 template <typename VectorType>
610 inline double
612 {
613  return matrix->constraints[index].second;
614 }
615 
616 
617 
618 template <typename VectorType>
619 inline void
621 {
622  Assert(index < matrix->constraints.size(), ExcIteratorPastEnd());
623  ++index;
624 }
625 
626 
627 
628 template <typename VectorType>
631  const size_type index)
632  : accessor(matrix, index)
633 {}
634 
635 
636 
637 template <typename VectorType>
640 {
641  accessor.advance();
642  return *this;
643 }
644 
645 
646 template <typename number>
647 inline const typename FilteredMatrix<number>::Accessor &
649 {
650  return accessor;
651 }
652 
653 
654 template <typename number>
655 inline const typename FilteredMatrix<number>::Accessor *
657 {
658  return &accessor;
659 }
660 
661 
662 template <typename number>
663 inline bool
665 operator==(const const_iterator &other) const
666 {
667  return (accessor.index == other.accessor.index &&
668  accessor.matrix == other.accessor.matrix);
669 }
670 
671 
672 template <typename number>
673 inline bool
675 operator!=(const const_iterator &other) const
676 {
677  return !(*this == other);
678 }
679 
680 
681 
682 //------------------------- FilteredMatrix ----------------------------------//
683 
684 template <typename number>
687 {
688  return const_iterator(this, 0);
689 }
690 
691 
692 template <typename number>
695 {
696  return const_iterator(this, constraints.size());
697 }
698 
699 
700 template <typename VectorType>
701 inline bool
703 operator()(const IndexValuePair &i1, const IndexValuePair &i2) const
704 {
705  return (i1.first < i2.first);
706 }
707 
708 
709 
710 template <typename VectorType>
711 template <typename MatrixType>
712 inline void
713 FilteredMatrix<VectorType>::initialize(const MatrixType &m, bool ecs)
714 {
715  matrix.reset(new_pointer_matrix_base(m, VectorType()));
716 
718 }
719 
720 
721 
722 template <typename VectorType>
725 {}
726 
727 
728 
729 template <typename VectorType>
731  : Subscriptor()
732  , expect_constrained_source(fm.expect_constrained_source)
733  , matrix(fm.matrix)
734  , constraints(fm.constraints)
735 {}
736 
737 
738 
739 template <typename VectorType>
740 template <typename MatrixType>
742  const bool ecs)
743  : expect_constrained_source(false)
744 {
745  initialize(m, ecs);
746 }
747 
748 
749 
750 template <typename VectorType>
753 {
754  matrix = fm.matrix;
755  expect_constrained_source = fm.expect_constrained_source;
756  constraints = fm.constraints;
757  return *this;
758 }
759 
760 
761 
762 template <typename VectorType>
763 inline void
765  const double value)
766 {
767  // add new constraint to end
768  constraints.push_back(IndexValuePair(index, value));
769 }
770 
771 
772 
773 template <typename VectorType>
774 template <class ConstraintList>
775 inline void
777  const ConstraintList &new_constraints)
778 {
779  // add new constraints to end
780  const size_type old_size = constraints.size();
781  constraints.reserve(old_size + new_constraints.size());
782  constraints.insert(constraints.end(),
783  new_constraints.begin(),
784  new_constraints.end());
785  // then merge the two arrays to
786  // form one sorted one
787  std::inplace_merge(constraints.begin(),
788  constraints.begin() + old_size,
789  constraints.end(),
790  PairComparison());
791 }
792 
793 
794 
795 template <typename VectorType>
796 inline void
798 {
799  // swap vectors to release memory
800  std::vector<IndexValuePair> empty;
801  constraints.swap(empty);
802 }
803 
804 
805 
806 template <typename VectorType>
807 inline void
809 {
810  clear_constraints();
811  matrix.reset();
812 }
813 
814 
815 
816 template <typename VectorType>
817 inline void
819  VectorType &v,
820  const bool /* matrix_is_symmetric */) const
821 {
822  apply_constraints(v);
823 }
824 
825 
826 template <typename VectorType>
827 inline void
829 {
831  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
832  tmp_vector->reinit(v);
833  const_index_value_iterator i = constraints.begin();
834  const const_index_value_iterator e = constraints.end();
835  for (; i != e; ++i)
836  {
837  AssertIsFinite(i->second);
838  (*tmp_vector)(i->first) = -i->second;
839  }
840 
841  // This vmult is without bc, to get
842  // the rhs correction in a correct
843  // way.
844  matrix->vmult_add(v, *tmp_vector);
845  // finally set constrained
846  // entries themselves
847  for (i = constraints.begin(); i != e; ++i)
848  {
849  AssertIsFinite(i->second);
850  v(i->first) = i->second;
851  }
852 }
853 
854 
855 template <typename VectorType>
856 inline void
858 {
859  // iterate over all constraints and
860  // zero out value
861  const_index_value_iterator i = constraints.begin();
862  const const_index_value_iterator e = constraints.end();
863  for (; i != e; ++i)
864  v(i->first) = 0;
865 }
866 
867 
868 
869 template <typename VectorType>
870 inline void
872  VectorType & out) const
873 {
874  // iterate over all constraints and
875  // set value correctly
876  const_index_value_iterator i = constraints.begin();
877  const const_index_value_iterator e = constraints.end();
878  for (; i != e; ++i)
879  {
880  AssertIsFinite(in(i->first));
881  out(i->first) = in(i->first);
882  }
883 }
884 
885 
886 
887 template <typename VectorType>
888 inline void
889 FilteredMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
890 {
891  if (!expect_constrained_source)
892  {
894  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
895  // first copy over src vector and
896  // pre-filter
897  tmp_vector->reinit(src, true);
898  *tmp_vector = src;
899  pre_filter(*tmp_vector);
900  // then let matrix do its work
901  matrix->vmult(dst, *tmp_vector);
902  }
903  else
904  {
905  matrix->vmult(dst, src);
906  }
907 
908  // finally do post-filtering
909  post_filter(src, dst);
910 }
911 
912 
913 
914 template <typename VectorType>
915 inline void
916 FilteredMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
917 {
918  if (!expect_constrained_source)
919  {
921  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
922  // first copy over src vector and
923  // pre-filter
924  tmp_vector->reinit(src, true);
925  *tmp_vector = src;
926  pre_filter(*tmp_vector);
927  // then let matrix do its work
928  matrix->Tvmult(dst, *tmp_vector);
929  }
930  else
931  {
932  matrix->Tvmult(dst, src);
933  }
934 
935  // finally do post-filtering
936  post_filter(src, dst);
937 }
938 
939 
940 
941 template <typename VectorType>
942 inline void
944  const VectorType &src) const
945 {
946  if (!expect_constrained_source)
947  {
949  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
950  // first copy over src vector and
951  // pre-filter
952  tmp_vector->reinit(src, true);
953  *tmp_vector = src;
954  pre_filter(*tmp_vector);
955  // then let matrix do its work
956  matrix->vmult_add(dst, *tmp_vector);
957  }
958  else
959  {
960  matrix->vmult_add(dst, src);
961  }
962 
963  // finally do post-filtering
964  post_filter(src, dst);
965 }
966 
967 
968 
969 template <typename VectorType>
970 inline void
972  const VectorType &src) const
973 {
974  if (!expect_constrained_source)
975  {
977  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
978  // first copy over src vector and
979  // pre-filter
980  tmp_vector->reinit(src, true);
981  *tmp_vector = src;
982  pre_filter(*tmp_vector);
983  // then let matrix do its work
984  matrix->Tvmult_add(dst, *tmp_vector);
985  }
986  else
987  {
988  matrix->Tvmult_add(dst, src);
989  }
990 
991  // finally do post-filtering
992  post_filter(src, dst);
993 }
994 
995 
996 
997 template <typename VectorType>
998 inline std::size_t
1000 {
1001  return (MemoryConsumption::memory_consumption(matrix) +
1003 }
1004 
1005 
1006 
1007 DEAL_II_NAMESPACE_CLOSE
1008 
1009 #endif
1010 /*---------------------------- 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
Subscriptor & operator=(const Subscriptor &)
Definition: subscriptor.h:283
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)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
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