00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __deal2__constraint_matrix_h
00014 #define __deal2__constraint_matrix_h
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/index_set.h>
00019 #include <deal.II/base/subscriptor.h>
00020 #include <deal.II/base/template_constraints.h>
00021
00022 #include <deal.II/lac/vector.h>
00023
00024 #include <vector>
00025 #include <map>
00026 #include <set>
00027 #include <utility>
00028 #include <complex>
00029
00030 #include <boost/scoped_ptr.hpp>
00031
00032
00033 DEAL_II_NAMESPACE_OPEN
00034
00035 template<int dim, class T> class Table;
00036 template <typename> class FullMatrix;
00037 class SparsityPattern;
00038 class CompressedSparsityPattern;
00039 class CompressedSetSparsityPattern;
00040 class CompressedSimpleSparsityPattern;
00041 class BlockSparsityPattern;
00042 class BlockCompressedSparsityPattern;
00043 class BlockCompressedSetSparsityPattern;
00044 class BlockCompressedSimpleSparsityPattern;
00045 template <typename number> class SparseMatrix;
00046 template <typename number> class BlockSparseMatrix;
00047 class BlockIndices;
00048 template <typename MatrixType> struct IsBlockMatrix;
00049
00050 namespace internals
00051 {
00052 class GlobalRowsFromLocal;
00053 }
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00123 class ConstraintMatrix : public Subscriptor
00124 {
00125 public:
00134 enum MergeConflictBehavior
00135 {
00142 no_conflicts_allowed,
00143
00153 left_object_wins,
00154
00164 right_object_wins
00165 };
00166
00180 ConstraintMatrix (const IndexSet & local_constraints = IndexSet());
00181
00185 ConstraintMatrix (const ConstraintMatrix &constraint_matrix);
00186
00197 void reinit (const IndexSet & local_constraints = IndexSet());
00198
00208 bool can_store_line (const unsigned int line_index) const;
00209
00245 void add_selected_constraints (const ConstraintMatrix &constraints_in,
00246 const IndexSet &filter);
00247
00258 void add_line (const unsigned int line);
00259
00280 void add_lines (const std::vector<bool> &lines);
00281
00302 void add_lines (const std::set<unsigned int> &lines);
00303
00324 void add_lines (const IndexSet &lines);
00325
00347 void add_entry (const unsigned int line,
00348 const unsigned int column,
00349 const double value);
00350
00360 void add_entries (const unsigned int line,
00361 const std::vector<std::pair<unsigned int,double> > &col_val_pairs);
00362
00371 void set_inhomogeneity (const unsigned int line,
00372 const double value);
00373
00407 void close ();
00408
00436 void merge (const ConstraintMatrix &other_constraints,
00437 const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed);
00438
00457 void shift (const unsigned int offset);
00458
00469 void clear ();
00470
00485 unsigned int n_constraints () const;
00486
00501 bool is_constrained (const unsigned int index) const;
00502
00522 bool is_identity_constrained (const unsigned int index) const;
00523
00541 unsigned int max_constraint_indirections () const;
00542
00549 bool is_inhomogeneously_constrained (const unsigned int index) const;
00550
00558 bool has_inhomogeneities () const;
00559
00566 const std::vector<std::pair<unsigned int,double> >*
00567 get_constraint_entries (const unsigned int line) const;
00568
00575 double get_inhomogeneity (const unsigned int line) const;
00576
00589 void print (std::ostream &) const;
00590
00610 void write_dot (std::ostream &) const;
00611
00617 std::size_t memory_consumption () const;
00618
00649 void condense (const SparsityPattern &uncondensed,
00650 SparsityPattern &condensed) const;
00651
00652
00671 void condense (SparsityPattern &sparsity) const;
00672
00678 void condense (BlockSparsityPattern &sparsity) const;
00679
00705 void condense (CompressedSparsityPattern &sparsity) const;
00706
00713 void condense (CompressedSetSparsityPattern &sparsity) const;
00714
00721 void condense (CompressedSimpleSparsityPattern &sparsity) const;
00722
00746 void condense (BlockCompressedSparsityPattern &sparsity) const;
00747
00753 void condense (BlockCompressedSetSparsityPattern &sparsity) const;
00754
00760 void condense (BlockCompressedSimpleSparsityPattern &sparsity) const;
00761
00762
00774 template<typename number>
00775 void condense (const SparseMatrix<number> &uncondensed,
00776 SparseMatrix<number> &condensed) const;
00777
00785 template<typename number>
00786 void condense (SparseMatrix<number> &matrix) const;
00787
00793 template <typename number>
00794 void condense (BlockSparseMatrix<number> &matrix) const;
00795
00816 template <class VectorType>
00817 void condense (const VectorType &uncondensed,
00818 VectorType &condensed) const;
00819
00835 template <class VectorType>
00836 void condense (VectorType &vec) const;
00837
00852 template<typename number, class VectorType>
00853 void condense (const SparseMatrix<number> &uncondensed_matrix,
00854 const VectorType &uncondensed_vector,
00855 SparseMatrix<number> &condensed_matrix,
00856 VectorType &condensed_vector) const;
00857
00866 template<typename number, class VectorType>
00867 void condense (SparseMatrix<number> &matrix,
00868 VectorType &vector) const;
00869
00875 template <typename number, class BlockVectorType>
00876 void condense (BlockSparseMatrix<number> &matrix,
00877 BlockVectorType &vector) const;
00878
00891 template <class VectorType>
00892 void set_zero (VectorType &vec) const;
00893
00958 template <class InVector, class OutVector>
00959 void
00960 distribute_local_to_global (const InVector &local_vector,
00961 const std::vector<unsigned int> &local_dof_indices,
00962 OutVector &global_vector) const;
00963
01035 template <typename VectorType>
01036 void
01037 distribute_local_to_global (const Vector<double> &local_vector,
01038 const std::vector<unsigned int> &local_dof_indices,
01039 VectorType &global_vector,
01040 const FullMatrix<double> &local_matrix) const;
01041
01042
01084 template <typename ForwardIteratorVec, typename ForwardIteratorInd,
01085 class VectorType>
01086 void
01087 distribute_local_to_global (ForwardIteratorVec local_vector_begin,
01088 ForwardIteratorVec local_vector_end,
01089 ForwardIteratorInd local_indices_begin,
01090 VectorType &global_vector) const;
01091
01159 template <typename MatrixType>
01160 void
01161 distribute_local_to_global (const FullMatrix<double> &local_matrix,
01162 const std::vector<unsigned int> &local_dof_indices,
01163 MatrixType &global_matrix) const;
01164
01170 template <typename MatrixType>
01171 void
01172 distribute_local_to_global (const FullMatrix<double> &local_matrix,
01173 const std::vector<unsigned int> &row_indices,
01174 const std::vector<unsigned int> &col_indices,
01175 MatrixType &global_matrix) const;
01176
01195 template <typename MatrixType, typename VectorType>
01196 void
01197 distribute_local_to_global (const FullMatrix<double> &local_matrix,
01198 const Vector<double> &local_vector,
01199 const std::vector<unsigned int> &local_dof_indices,
01200 MatrixType &global_matrix,
01201 VectorType &global_vector,
01202 bool use_inhomogeneities_for_rhs = false) const;
01203
01283 template <typename SparsityType>
01284 void
01285 add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
01286 SparsityType &sparsity_pattern,
01287 const bool keep_constrained_entries = true,
01288 const Table<2,bool> &dof_mask = default_empty_table) const;
01289
01296 template <typename SparsityType>
01297 void
01298 add_entries_local_to_global (const std::vector<unsigned int> &row_indices,
01299 const std::vector<unsigned int> &col_indices,
01300 SparsityType &sparsity_pattern,
01301 const bool keep_constrained_entries = true,
01302 const Table<2,bool> &dof_mask = default_empty_table) const;
01303
01338 template <typename ForwardIteratorVec, typename ForwardIteratorInd,
01339 class VectorType>
01340 void
01341 get_dof_values (const VectorType &global_vector,
01342 ForwardIteratorInd local_indices_begin,
01343 ForwardIteratorVec local_vector_begin,
01344 ForwardIteratorVec local_vector_end) const;
01345
01374 template <class VectorType>
01375 void distribute (const VectorType &condensed,
01376 VectorType &uncondensed) const;
01377
01392 template <class VectorType>
01393 void distribute (VectorType &vec) const;
01394
01404 DeclException0 (ExcMatrixIsClosed);
01410 DeclException0 (ExcMatrixNotClosed);
01416 DeclException1 (ExcLineInexistant,
01417 unsigned int,
01418 << "The specified line " << arg1
01419 << " does not exist.");
01425 DeclException4 (ExcEntryAlreadyExists,
01426 int, int, double, double,
01427 << "The entry for the indices " << arg1 << " and "
01428 << arg2 << " already exists, but the values "
01429 << arg3 << " (old) and " << arg4 << " (new) differ "
01430 << "by " << (arg4-arg3) << ".");
01436 DeclException2 (ExcDoFConstrainedToConstrainedDoF,
01437 int, int,
01438 << "You tried to constrain DoF " << arg1
01439 << " to DoF " << arg2
01440 << ", but that one is also constrained. This is not allowed!");
01446 DeclException1 (ExcDoFIsConstrainedFromBothObjects,
01447 int,
01448 << "Degree of freedom " << arg1
01449 << " is constrained from both object in a merge operation.");
01455 DeclException1 (ExcDoFIsConstrainedToConstrainedDoF,
01456 int,
01457 << "In the given argument a degree of freedom is constrained "
01458 << "to another DoF with number " << arg1
01459 << ", which however is constrained by this object. This is not"
01460 << " allowed.");
01466 DeclException1 (ExcRowNotStoredHere,
01467 int,
01468 << "The index set given to this constraint matrix indicates "
01469 << "constraints for degree of freedom " << arg1
01470 << " should not be stored by this object, but a constraint "
01471 << "is being added.");
01472
01473 private:
01474
01479 struct ConstraintLine
01480 {
01486 typedef std::vector<std::pair<unsigned int,double> > Entries;
01487
01495 unsigned int line;
01496
01507 Entries entries;
01508
01512 double inhomogeneity;
01513
01526 bool operator < (const ConstraintLine &) const;
01527
01536 bool operator == (const ConstraintLine &) const;
01537
01543 std::size_t memory_consumption () const;
01544 };
01545
01564 std::vector<ConstraintLine> lines;
01565
01623 std::vector<unsigned int> lines_cache;
01624
01633 IndexSet local_lines;
01634
01639 bool sorted;
01640
01646 unsigned int calculate_line_index (const unsigned int line) const;
01647
01655 static bool check_zero_weight (const std::pair<unsigned int, double> &p);
01656
01662 static const Table<2,bool> default_empty_table;
01663
01669 template <typename MatrixType, typename VectorType>
01670 void
01671 distribute_local_to_global (const FullMatrix<double> &local_matrix,
01672 const Vector<double> &local_vector,
01673 const std::vector<unsigned int> &local_dof_indices,
01674 MatrixType &global_matrix,
01675 VectorType &global_vector,
01676 bool use_inhomogeneities_for_rhs,
01677 internal::bool2type<false>) const;
01678
01684 template <typename MatrixType, typename VectorType>
01685 void
01686 distribute_local_to_global (const FullMatrix<double> &local_matrix,
01687 const Vector<double> &local_vector,
01688 const std::vector<unsigned int> &local_dof_indices,
01689 MatrixType &global_matrix,
01690 VectorType &global_vector,
01691 bool use_inhomogeneities_for_rhs,
01692 internal::bool2type<true>) const;
01693
01699 template <typename SparsityType>
01700 void
01701 add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
01702 SparsityType &sparsity_pattern,
01703 const bool keep_constrained_entries,
01704 const Table<2,bool> &dof_mask,
01705 internal::bool2type<false>) const;
01706
01712 template <typename SparsityType>
01713 void
01714 add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
01715 SparsityType &sparsity_pattern,
01716 const bool keep_constrained_entries,
01717 const Table<2,bool> &dof_mask,
01718 internal::bool2type<true>) const;
01719
01730 void
01731 make_sorted_row_list (const std::vector<unsigned int> &local_dof_indices,
01732 internals::GlobalRowsFromLocal &global_rows) const;
01733
01744 void
01745 make_sorted_row_list (const std::vector<unsigned int> &local_dof_indices,
01746 std::vector<unsigned int> &active_dofs) const;
01747
01752 double
01753 resolve_vector_entry (const unsigned int i,
01754 const internals::GlobalRowsFromLocal &global_rows,
01755 const Vector<double> &local_vector,
01756 const std::vector<unsigned int> &local_dof_indices,
01757 const FullMatrix<double> &local_matrix) const;
01758 };
01759
01760
01761
01762
01763
01764 inline
01765 ConstraintMatrix::ConstraintMatrix (const IndexSet &local_constraints)
01766 :
01767 lines (),
01768 local_lines (local_constraints),
01769 sorted (false)
01770 {}
01771
01772
01773
01774 inline
01775 ConstraintMatrix::ConstraintMatrix (const ConstraintMatrix &constraint_matrix)
01776 :
01777 Subscriptor (),
01778 lines (constraint_matrix.lines),
01779 lines_cache (constraint_matrix.lines_cache),
01780 local_lines (constraint_matrix.local_lines),
01781 sorted (constraint_matrix.sorted)
01782 {}
01783
01784
01785 inline
01786 void
01787 ConstraintMatrix::add_line (const unsigned int line)
01788 {
01789 Assert (sorted==false, ExcMatrixIsClosed());
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801 Assert (line != numbers::invalid_unsigned_int,
01802 ExcInternalError());
01803 const unsigned int line_index = calculate_line_index (line);
01804
01805
01806
01807 if (is_constrained(line))
01808 return;
01809
01810
01811
01812 if (line_index >= lines_cache.size())
01813 lines_cache.resize (std::max(2*static_cast<unsigned int>(lines_cache.size()),
01814 line_index+1),
01815 numbers::invalid_unsigned_int);
01816
01817
01818
01819 lines.push_back (ConstraintLine());
01820 lines.back().line = line;
01821 lines.back().inhomogeneity = 0.;
01822 lines_cache[line_index] = lines.size()-1;
01823 }
01824
01825
01826
01827 inline
01828 void
01829 ConstraintMatrix::add_entry (const unsigned int line,
01830 const unsigned int column,
01831 const double value)
01832 {
01833 Assert (sorted==false, ExcMatrixIsClosed());
01834 Assert (line != column,
01835 ExcMessage ("Can't constrain a degree of freedom to itself"));
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845 Assert (lines_cache[calculate_line_index(line)] != numbers::invalid_unsigned_int,
01846 ExcInternalError());
01847 ConstraintLine* line_ptr = &lines[lines_cache[calculate_line_index(line)]];
01848 Assert (line_ptr->line == line, ExcInternalError());
01849 for (ConstraintLine::Entries::const_iterator
01850 p=line_ptr->entries.begin();
01851 p != line_ptr->entries.end(); ++p)
01852 if (p->first == column)
01853 {
01854 Assert (std::fabs(p->second - value) < 1.e-14,
01855 ExcEntryAlreadyExists(line, column, p->second, value));
01856 return;
01857 }
01858
01859 line_ptr->entries.push_back (std::make_pair(column,value));
01860 }
01861
01862
01863
01864 inline
01865 void
01866 ConstraintMatrix::set_inhomogeneity (const unsigned int line,
01867 const double value)
01868 {
01869 const unsigned int line_index = calculate_line_index(line);
01870 Assert( line_index < lines_cache.size() &&
01871 lines_cache[line_index] != numbers::invalid_unsigned_int,
01872 ExcMessage("call add_line() before calling set_inhomogeneity()"));
01873 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
01874 ConstraintLine* line_ptr = &lines[lines_cache[line_index]];
01875 line_ptr->inhomogeneity = value;
01876 }
01877
01878
01879
01880 inline
01881 bool
01882 ConstraintMatrix::is_constrained (const unsigned int index) const
01883 {
01884 const unsigned int line_index = calculate_line_index(index);
01885 return ((line_index < lines_cache.size())
01886 &&
01887 (lines_cache[line_index] != numbers::invalid_unsigned_int));
01888 }
01889
01890
01891
01892 inline
01893 bool
01894 ConstraintMatrix::is_inhomogeneously_constrained (const unsigned int index) const
01895 {
01896
01897
01898
01899 const unsigned int line_index = calculate_line_index(index);
01900 if (line_index >= lines_cache.size() ||
01901 lines_cache[line_index] == numbers::invalid_unsigned_int)
01902 return false;
01903 else
01904 {
01905 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
01906 return !(lines[lines_cache[line_index]].inhomogeneity == 0);
01907 }
01908 }
01909
01910
01911
01912 inline
01913 const std::vector<std::pair<unsigned int,double> >*
01914 ConstraintMatrix::get_constraint_entries (unsigned int line) const
01915 {
01916
01917
01918
01919 const unsigned int line_index = calculate_line_index(line);
01920 if (line_index >= lines_cache.size() ||
01921 lines_cache[line_index] == numbers::invalid_unsigned_int)
01922 return 0;
01923 else
01924 return &lines[lines_cache[line_index]].entries;
01925 }
01926
01927
01928
01929 inline
01930 double
01931 ConstraintMatrix::get_inhomogeneity (const unsigned int line) const
01932 {
01933
01934
01935
01936 const unsigned int line_index = calculate_line_index(line);
01937 if (line_index >= lines_cache.size() ||
01938 lines_cache[line_index] == numbers::invalid_unsigned_int)
01939 return 0;
01940 else
01941 return lines[lines_cache[line_index]].inhomogeneity;
01942 }
01943
01944
01945
01946 inline unsigned int
01947 ConstraintMatrix::calculate_line_index (const unsigned int line) const
01948 {
01949
01950 if (!local_lines.size())
01951 return line;
01952
01953 Assert(local_lines.is_element(line),
01954 ExcRowNotStoredHere(line));
01955
01956 return local_lines.index_within_set(line);
01957 }
01958
01959
01960
01961 inline bool
01962 ConstraintMatrix::can_store_line (unsigned int line_index) const
01963 {
01964 return !local_lines.size() || local_lines.is_element(line_index);
01965 }
01966
01967
01968
01969 template <class InVector, class OutVector>
01970 inline
01971 void
01972 ConstraintMatrix::
01973 distribute_local_to_global (const InVector &local_vector,
01974 const std::vector<unsigned int> &local_dof_indices,
01975 OutVector &global_vector) const
01976 {
01977 Assert (local_vector.size() == local_dof_indices.size(),
01978 ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
01979 distribute_local_to_global (local_vector.begin(), local_vector.end(),
01980 local_dof_indices.begin(), global_vector);
01981 }
01982
01983
01984
01985 template <typename ForwardIteratorVec, typename ForwardIteratorInd,
01986 class VectorType>
01987 inline
01988 void ConstraintMatrix::
01989 distribute_local_to_global (ForwardIteratorVec local_vector_begin,
01990 ForwardIteratorVec local_vector_end,
01991 ForwardIteratorInd local_indices_begin,
01992 VectorType &global_vector) const
01993 {
01994 Assert (sorted == true, ExcMatrixNotClosed());
01995 for ( ; local_vector_begin != local_vector_end;
01996 ++local_vector_begin, ++local_indices_begin)
01997 {
01998 if (is_constrained(*local_indices_begin) == false)
01999 global_vector(*local_indices_begin) += *local_vector_begin;
02000 else
02001 {
02002 const ConstraintLine& position =
02003 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
02004 for (unsigned int j=0; j<position.entries.size(); ++j)
02005 global_vector(position.entries[j].first)
02006 += *local_vector_begin * position.entries[j].second;
02007 }
02008 }
02009 }
02010
02011
02012
02013 template <typename ForwardIteratorVec, typename ForwardIteratorInd,
02014 class VectorType>
02015 inline
02016 void ConstraintMatrix::get_dof_values (const VectorType &global_vector,
02017 ForwardIteratorInd local_indices_begin,
02018 ForwardIteratorVec local_vector_begin,
02019 ForwardIteratorVec local_vector_end) const
02020 {
02021 Assert (sorted == true, ExcMatrixNotClosed());
02022 for ( ; local_vector_begin != local_vector_end;
02023 ++local_vector_begin, ++local_indices_begin)
02024 {
02025 if (is_constrained(*local_indices_begin) == false)
02026 *local_vector_begin = global_vector(*local_indices_begin);
02027 else
02028 {
02029 const ConstraintLine & position =
02030 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
02031 typename VectorType::value_type value = position.inhomogeneity;
02032 for (unsigned int j=0; j<position.entries.size(); ++j)
02033 value += (global_vector(position.entries[j].first) *
02034 position.entries[j].second);
02035 *local_vector_begin = value;
02036 }
02037 }
02038 }
02039
02040
02041
02042 template <typename MatrixType>
02043 inline
02044 void
02045 ConstraintMatrix::
02046 distribute_local_to_global (const FullMatrix<double> &local_matrix,
02047 const std::vector<unsigned int> &local_dof_indices,
02048 MatrixType &global_matrix) const
02049 {
02050
02051
02052
02053 Vector<double> dummy(0);
02054 distribute_local_to_global (local_matrix, dummy, local_dof_indices,
02055 global_matrix, dummy, false,
02056 internal::bool2type<IsBlockMatrix<MatrixType>::value>());
02057 }
02058
02059
02060
02061 template <typename MatrixType, typename VectorType>
02062 inline
02063 void
02064 ConstraintMatrix::
02065 distribute_local_to_global (const FullMatrix<double> &local_matrix,
02066 const Vector<double> &local_vector,
02067 const std::vector<unsigned int> &local_dof_indices,
02068 MatrixType &global_matrix,
02069 VectorType &global_vector,
02070 bool use_inhomogeneities_for_rhs) const
02071 {
02072
02073
02074
02075
02076 distribute_local_to_global (local_matrix, local_vector, local_dof_indices,
02077 global_matrix, global_vector, use_inhomogeneities_for_rhs,
02078 internal::bool2type<IsBlockMatrix<MatrixType>::value>());
02079 }
02080
02081
02082
02083 template <typename SparsityType>
02084 inline
02085 void
02086 ConstraintMatrix::
02087 add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
02088 SparsityType &sparsity_pattern,
02089 const bool keep_constrained_entries,
02090 const Table<2,bool> &dof_mask) const
02091 {
02092
02093
02094
02095
02096 add_entries_local_to_global (local_dof_indices, sparsity_pattern,
02097 keep_constrained_entries, dof_mask,
02098 internal::bool2type<IsBlockMatrix<SparsityType>::value>());
02099 }
02100
02101
02102 DEAL_II_NAMESPACE_CLOSE
02103
02104 #endif
02105
02106