00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __deal2__mesh_worker_simple_h
00014 #define __deal2__mesh_worker_simple_h
00015
00016 #include <deal.II/base/named_data.h>
00017 #include <deal.II/base/smartpointer.h>
00018 #include <deal.II/base/mg_level_object.h>
00019 #include <deal.II/lac/block_vector.h>
00020 #include <deal.II/meshworker/dof_info.h>
00021 #include <deal.II/meshworker/functional.h>
00022 #include <deal.II/multigrid/mg_constrained_dofs.h>
00023
00035 DEAL_II_NAMESPACE_OPEN
00036
00037 namespace MeshWorker
00038 {
00039 namespace Assembler
00040 {
00054 template <class VECTOR>
00055 class ResidualSimple
00056 {
00057 public:
00058 void initialize(NamedData<VECTOR*>& results);
00062 void initialize(const ConstraintMatrix& constraints);
00074 template <class DOFINFO>
00075 void initialize_info(DOFINFO& info, bool face) const;
00076
00087 template <class DOFINFO>
00088 void assemble(const DOFINFO& info);
00089
00094 template<class DOFINFO>
00095 void assemble(const DOFINFO& info1,
00096 const DOFINFO& info2);
00097 private:
00102 NamedData<SmartPointer<VECTOR,ResidualSimple<VECTOR> > > residuals;
00106 SmartPointer<const ConstraintMatrix,ResidualSimple<VECTOR> > constraints;
00107 };
00108
00137 template <class MATRIX>
00138 class MatrixSimple
00139 {
00140 public:
00148 MatrixSimple(double threshold = 1.e-12);
00149
00154 void initialize(MATRIX& m);
00168 void initialize(const ConstraintMatrix& constraints);
00169
00201 void initialize_local_blocks(const BlockIndices& local_indices);
00202
00214 template <class DOFINFO>
00215 void initialize_info(DOFINFO& info, bool face) const;
00216
00222 template<class DOFINFO>
00223 void assemble(const DOFINFO& info);
00224
00231 template<class DOFINFO>
00232 void assemble(const DOFINFO& info1,
00233 const DOFINFO& info2);
00234 private:
00239 void assemble(const FullMatrix<double>& M,
00240 const std::vector<unsigned int>& i1,
00241 const std::vector<unsigned int>& i2);
00242
00247 SmartPointer<MATRIX,MatrixSimple<MATRIX> > matrix;
00252 SmartPointer<const ConstraintMatrix,MatrixSimple<MATRIX> > constraints;
00253
00262 BlockIndices local_indices;
00263
00273 const double threshold;
00274
00275 };
00276
00277
00288 template <class MATRIX>
00289 class MGMatrixSimple
00290 {
00291 public:
00299 MGMatrixSimple(double threshold = 1.e-12);
00300
00305 void initialize(MGLevelObject<MATRIX>& m);
00306
00311 void initialize(const MGConstrainedDoFs& mg_constrained_dofs);
00312
00343 void initialize_local_blocks(const BlockIndices& local_indices);
00344
00352 void initialize_fluxes(MGLevelObject<MATRIX>& flux_up,
00353 MGLevelObject<MATRIX>& flux_down);
00354
00363 void initialize_interfaces(MGLevelObject<MATRIX>& interface_in,
00364 MGLevelObject<MATRIX>& interface_out);
00378 template <class DOFINFO>
00379 void initialize_info(DOFINFO& info, bool face) const;
00380
00386 template<class DOFINFO>
00387 void assemble(const DOFINFO& info);
00388
00395 template<class DOFINFO>
00396 void assemble(const DOFINFO& info1,
00397 const DOFINFO& info2);
00398 private:
00403 void assemble(MATRIX& G,
00404 const FullMatrix<double>& M,
00405 const std::vector<unsigned int>& i1,
00406 const std::vector<unsigned int>& i2);
00407
00412 void assemble(MATRIX& G,
00413 const FullMatrix<double>& M,
00414 const std::vector<unsigned int>& i1,
00415 const std::vector<unsigned int>& i2,
00416 const unsigned int level);
00417
00423 void assemble_up(MATRIX& G,
00424 const FullMatrix<double>& M,
00425 const std::vector<unsigned int>& i1,
00426 const std::vector<unsigned int>& i2,
00427 const unsigned int level = numbers::invalid_unsigned_int);
00433 void assemble_down(MATRIX& G,
00434 const FullMatrix<double>& M,
00435 const std::vector<unsigned int>& i1,
00436 const std::vector<unsigned int>& i2,
00437 const unsigned int level = numbers::invalid_unsigned_int);
00438
00444 void assemble_in(MATRIX& G,
00445 const FullMatrix<double>& M,
00446 const std::vector<unsigned int>& i1,
00447 const std::vector<unsigned int>& i2,
00448 const unsigned int level = numbers::invalid_unsigned_int);
00449
00455 void assemble_out(MATRIX& G,
00456 const FullMatrix<double>& M,
00457 const std::vector<unsigned int>& i1,
00458 const std::vector<unsigned int>& i2,
00459 const unsigned int level = numbers::invalid_unsigned_int);
00460
00465 SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > matrix;
00466
00473 SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > flux_up;
00474
00481 SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > flux_down;
00482
00490 SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_in;
00491
00499 SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_out;
00503 SmartPointer<const MGConstrainedDoFs,MGMatrixSimple<MATRIX> > mg_constrained_dofs;
00504
00513 BlockIndices local_indices;
00514
00524 const double threshold;
00525
00526 };
00527
00528
00539 template <class MATRIX, class VECTOR>
00540 class SystemSimple :
00541 private MatrixSimple<MATRIX>,
00542 private ResidualSimple<VECTOR>
00543 {
00544 public:
00550 SystemSimple(double threshold = 1.e-12);
00551
00556 void initialize(MATRIX& m, VECTOR& rhs);
00557
00571 template <class DOFINFO>
00572 void initialize_info(DOFINFO& info, bool face) const;
00573
00579 template<class DOFINFO>
00580 void assemble(const DOFINFO& info);
00581
00588 template<class DOFINFO>
00589 void assemble(const DOFINFO& info1,
00590 const DOFINFO& info2);
00591 };
00592
00593
00594
00595
00596 template <class VECTOR>
00597 inline void
00598 ResidualSimple<VECTOR>::initialize(NamedData<VECTOR*>& results)
00599 {
00600 residuals = results;
00601 }
00602
00603 template <class VECTOR>
00604 inline void
00605 ResidualSimple<VECTOR>::initialize(const ConstraintMatrix& c)
00606 {
00607 constraints = &c;
00608 }
00609
00610
00611 template <class VECTOR>
00612 template <class DOFINFO>
00613 inline void
00614 ResidualSimple<VECTOR>::initialize_info(DOFINFO& info, bool) const
00615 {
00616 info.initialize_vectors(residuals.size());
00617 }
00618
00619
00620 template <class VECTOR>
00621 template <class DOFINFO>
00622 inline void
00623 ResidualSimple<VECTOR>::assemble(const DOFINFO& info)
00624 {
00625 for (unsigned int k=0;k<residuals.size();++k)
00626 {
00627 if(constraints == 0)
00628 {
00629 for (unsigned int i=0;i<info.vector(k).block(0).size();++i)
00630 (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
00631 }
00632 else
00633 constraints->distribute_local_to_global(
00634 info.vector(k).block(0), info.indices, (*residuals(k)));
00635 }
00636 }
00637
00638
00639 template <class VECTOR>
00640 template <class DOFINFO>
00641 inline void
00642 ResidualSimple<VECTOR>::assemble(const DOFINFO& info1,
00643 const DOFINFO& info2)
00644 {
00645 for (unsigned int k=0;k<residuals.size();++k)
00646 {
00647 if(constraints == 0)
00648 {
00649 for (unsigned int i=0;i<info1.vector(k).block(0).size();++i)
00650 (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
00651 for (unsigned int i=0;i<info2.vector(k).block(0).size();++i)
00652 (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
00653 }
00654 else
00655 {
00656 constraints->distribute_local_to_global(
00657 info1.vector(k).block(0), info1.indices, (*residuals(k)));
00658 constraints->distribute_local_to_global(
00659 info2.vector(k).block(0), info2.indices, (*residuals(k)));
00660 }
00661 }
00662 }
00663
00664
00665
00666
00667 template <class MATRIX>
00668 inline
00669 MatrixSimple<MATRIX>::MatrixSimple(double threshold)
00670 :
00671 threshold(threshold)
00672 {}
00673
00674
00675 template <class MATRIX>
00676 inline void
00677 MatrixSimple<MATRIX>::initialize(MATRIX& m)
00678 {
00679 matrix = &m;
00680 }
00681
00682
00683 template <class MATRIX>
00684 inline void
00685 MatrixSimple<MATRIX>::initialize(const ConstraintMatrix& c)
00686 {
00687 constraints = &c;
00688 }
00689
00690
00691 template <class MATRIX>
00692 inline void
00693 MatrixSimple<MATRIX>::initialize_local_blocks(const BlockIndices& b)
00694 {
00695 local_indices = b;
00696 }
00697
00698
00699 template <class MATRIX >
00700 template <class DOFINFO>
00701 inline void
00702 MatrixSimple<MATRIX>::initialize_info(DOFINFO& info, bool face) const
00703 {
00704 const unsigned int n = local_indices.size();
00705
00706 if (n == 0)
00707 info.initialize_matrices(1, face);
00708 else
00709 {
00710 info.initialize_matrices(n*n, face);
00711 unsigned int k=0;
00712 for (unsigned int i=0;i<n;++i)
00713 for (unsigned int j=0;j<n;++j,++k)
00714 {
00715 info.matrix(k,false).row = i;
00716 info.matrix(k,false).column = j;
00717 if (face)
00718 {
00719 info.matrix(k,true).row = i;
00720 info.matrix(k,true).column = j;
00721 }
00722 }
00723 }
00724 }
00725
00726
00727
00728 template <class MATRIX>
00729 inline void
00730 MatrixSimple<MATRIX>::assemble(const FullMatrix<double>& M,
00731 const std::vector<unsigned int>& i1,
00732 const std::vector<unsigned int>& i2)
00733 {
00734 AssertDimension(M.m(), i1.size());
00735 AssertDimension(M.n(), i2.size());
00736
00737 if(constraints == 0)
00738 {
00739 for (unsigned int j=0; j<i1.size(); ++j)
00740 for (unsigned int k=0; k<i2.size(); ++k)
00741 if (std::fabs(M(j,k)) >= threshold)
00742 matrix->add(i1[j], i2[k], M(j,k));
00743 }
00744 else
00745 constraints->distribute_local_to_global(M, i1, i2, *matrix);
00746 }
00747
00748
00749 template <class MATRIX>
00750 template <class DOFINFO>
00751 inline void
00752 MatrixSimple<MATRIX>::assemble(const DOFINFO& info)
00753 {
00754 if (local_indices.size() == 0)
00755 assemble(info.matrix(0,false).matrix, info.indices, info.indices);
00756 else
00757 {
00758 for (unsigned int k=0;k<info.n_matrices();++k)
00759 {
00760 assemble(info.matrix(k,false).matrix,
00761 info.indices_by_block[info.matrix(k,false).row],
00762 info.indices_by_block[info.matrix(k,false).column]);
00763 }
00764 }
00765 }
00766
00767
00768 template <class MATRIX>
00769 template <class DOFINFO>
00770 inline void
00771 MatrixSimple<MATRIX>::assemble(const DOFINFO& info1, const DOFINFO& info2)
00772 {
00773 if (local_indices.size() == 0)
00774 {
00775 assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices);
00776 assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices);
00777 assemble(info2.matrix(0,false).matrix, info2.indices, info2.indices);
00778 assemble(info2.matrix(0,true).matrix, info2.indices, info1.indices);
00779 }
00780 else
00781 for (unsigned int k=0;k<info1.n_matrices();++k)
00782 {
00783 const unsigned int row = info1.matrix(k,false).row;
00784 const unsigned int column = info1.matrix(k,false).column;
00785
00786 assemble(info1.matrix(k,false).matrix,
00787 info1.indices_by_block[row], info1.indices_by_block[column]);
00788 assemble(info1.matrix(k,true).matrix,
00789 info1.indices_by_block[row], info2.indices_by_block[column]);
00790 assemble(info2.matrix(k,false).matrix,
00791 info2.indices_by_block[row], info2.indices_by_block[column]);
00792 assemble(info2.matrix(k,true).matrix,
00793 info2.indices_by_block[row], info1.indices_by_block[column]);
00794 }
00795 }
00796
00797
00798
00799
00800 template <class MATRIX>
00801 inline
00802 MGMatrixSimple<MATRIX>::MGMatrixSimple(double threshold)
00803 :
00804 threshold(threshold)
00805 {}
00806
00807
00808 template <class MATRIX>
00809 inline void
00810 MGMatrixSimple<MATRIX>::initialize(MGLevelObject<MATRIX>& m)
00811 {
00812 matrix = &m;
00813 }
00814
00815 template <class MATRIX>
00816 inline void
00817 MGMatrixSimple<MATRIX>::initialize(const MGConstrainedDoFs& c)
00818 {
00819 mg_constrained_dofs = &c;
00820 }
00821
00822 template <class MATRIX>
00823 inline void
00824 MGMatrixSimple<MATRIX>::initialize_local_blocks(const BlockIndices& b)
00825 {
00826 local_indices = b;
00827 }
00828
00829
00830 template <class MATRIX>
00831 inline void
00832 MGMatrixSimple<MATRIX>::initialize_fluxes(
00833 MGLevelObject<MATRIX>& up, MGLevelObject<MATRIX>& down)
00834 {
00835 flux_up = &up;
00836 flux_down = &down;
00837 }
00838
00839
00840 template <class MATRIX>
00841 inline void
00842 MGMatrixSimple<MATRIX>::initialize_interfaces(
00843 MGLevelObject<MATRIX>& in, MGLevelObject<MATRIX>& out)
00844 {
00845 interface_in = ∈
00846 interface_out = &out;
00847 }
00848
00849
00850 template <class MATRIX >
00851 template <class DOFINFO>
00852 inline void
00853 MGMatrixSimple<MATRIX>::initialize_info(DOFINFO& info, bool face) const
00854 {
00855 const unsigned int n = local_indices.size();
00856
00857 if (n == 0)
00858 info.initialize_matrices(1, face);
00859 else
00860 {
00861 info.initialize_matrices(n*n, face);
00862 unsigned int k=0;
00863 for (unsigned int i=0;i<n;++i)
00864 for (unsigned int j=0;j<n;++j,++k)
00865 {
00866 info.matrix(k,false).row = i;
00867 info.matrix(k,false).column = j;
00868 if (face)
00869 {
00870 info.matrix(k,true).row = i;
00871 info.matrix(k,true).column = j;
00872 }
00873 }
00874 }
00875 }
00876
00877
00878 template <class MATRIX>
00879 inline void
00880 MGMatrixSimple<MATRIX>::assemble(
00881 MATRIX& G,
00882 const FullMatrix<double>& M,
00883 const std::vector<unsigned int>& i1,
00884 const std::vector<unsigned int>& i2)
00885 {
00886 AssertDimension(M.m(), i1.size());
00887 AssertDimension(M.n(), i2.size());
00888
00889 if(mg_constrained_dofs == 0)
00890 {
00891 for (unsigned int j=0; j<i1.size(); ++j)
00892 for (unsigned int k=0; k<i2.size(); ++k)
00893 if (std::fabs(M(j,k)) >= threshold)
00894 G.add(i1[j], i2[k], M(j,k));
00895 }
00896 else
00897 {
00898 for (unsigned int j=0; j<i1.size(); ++j)
00899 for (unsigned int k=0; k<i2.size(); ++k)
00900 if (std::fabs(M(j,k)) >= threshold)
00901 {
00902 if(!mg_constrained_dofs->continuity_across_refinement_edges())
00903 G.add(i1[j], i2[k], M(j,k));
00904 }
00905 }
00906 }
00907
00908
00909 template <class MATRIX>
00910 inline void
00911 MGMatrixSimple<MATRIX>::assemble(
00912 MATRIX& G,
00913 const FullMatrix<double>& M,
00914 const std::vector<unsigned int>& i1,
00915 const std::vector<unsigned int>& i2,
00916 const unsigned int level)
00917 {
00918 AssertDimension(M.m(), i1.size());
00919 AssertDimension(M.n(), i2.size());
00920
00921 if(mg_constrained_dofs == 0)
00922 {
00923 for (unsigned int j=0; j<i1.size(); ++j)
00924 for (unsigned int k=0; k<i2.size(); ++k)
00925 if (std::fabs(M(j,k)) >= threshold)
00926 G.add(i1[j], i2[k], M(j,k));
00927 }
00928 else
00929 {
00930 for (unsigned int j=0; j<i1.size(); ++j)
00931 for (unsigned int k=0; k<i2.size(); ++k)
00932 if (std::fabs(M(j,k)) >= threshold)
00933 if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
00934 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
00935 {
00936 if (mg_constrained_dofs->set_boundary_values())
00937 {
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947 if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
00948 !mg_constrained_dofs->is_boundary_index(level, i2[k]))
00949 ||
00950 (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
00951 mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
00952 i1[j] == i2[k]))
00953 G.add(i1[j], i2[k], M(j,k));
00954 }
00955 else
00956 G.add(i1[j], i2[k], M(j,k));
00957 }
00958 }
00959 }
00960
00961
00962 template <class MATRIX>
00963 inline void
00964 MGMatrixSimple<MATRIX>::assemble_up(
00965 MATRIX& G,
00966 const FullMatrix<double>& M,
00967 const std::vector<unsigned int>& i1,
00968 const std::vector<unsigned int>& i2,
00969 const unsigned int level)
00970 {
00971 AssertDimension(M.n(), i1.size());
00972 AssertDimension(M.m(), i2.size());
00973
00974 if(mg_constrained_dofs == 0)
00975 {
00976 for (unsigned int j=0; j<i1.size(); ++j)
00977 for (unsigned int k=0; k<i2.size(); ++k)
00978 if (std::fabs(M(k,j)) >= threshold)
00979 G.add(i1[j], i2[k], M(k,j));
00980 }
00981 else
00982 {
00983 for (unsigned int j=0; j<i1.size(); ++j)
00984 for (unsigned int k=0; k<i2.size(); ++k)
00985 if (std::fabs(M(k,j)) >= threshold)
00986 if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
00987 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
00988 G.add(i1[j], i2[k], M(k,j));
00989 }
00990 }
00991
00992 template <class MATRIX>
00993 inline void
00994 MGMatrixSimple<MATRIX>::assemble_down(
00995 MATRIX& G,
00996 const FullMatrix<double>& M,
00997 const std::vector<unsigned int>& i1,
00998 const std::vector<unsigned int>& i2,
00999 const unsigned int level)
01000 {
01001 AssertDimension(M.m(), i1.size());
01002 AssertDimension(M.n(), i2.size());
01003
01004 if(mg_constrained_dofs == 0)
01005 {
01006 for (unsigned int j=0; j<i1.size(); ++j)
01007 for (unsigned int k=0; k<i2.size(); ++k)
01008 if (std::fabs(M(j,k)) >= threshold)
01009 G.add(i1[j], i2[k], M(j,k));
01010 }
01011 else
01012 {
01013 for (unsigned int j=0; j<i1.size(); ++j)
01014 for (unsigned int k=0; k<i2.size(); ++k)
01015 if (std::fabs(M(j,k)) >= threshold)
01016 if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
01017 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
01018 G.add(i1[j], i2[k], M(j,k));
01019 }
01020 }
01021
01022 template <class MATRIX>
01023 inline void
01024 MGMatrixSimple<MATRIX>::assemble_in(
01025 MATRIX& G,
01026 const FullMatrix<double>& M,
01027 const std::vector<unsigned int>& i1,
01028 const std::vector<unsigned int>& i2,
01029 const unsigned int level)
01030 {
01031 AssertDimension(M.m(), i1.size());
01032 AssertDimension(M.n(), i2.size());
01033
01034 if(mg_constrained_dofs == 0)
01035 {
01036 for (unsigned int j=0; j<i1.size(); ++j)
01037 for (unsigned int k=0; k<i2.size(); ++k)
01038 if (std::fabs(M(j,k)) >= threshold)
01039 G.add(i1[j], i2[k], M(j,k));
01040 }
01041 else
01042 {
01043 for (unsigned int j=0; j<i1.size(); ++j)
01044 for (unsigned int k=0; k<i2.size(); ++k)
01045 if (std::fabs(M(j,k)) >= threshold)
01046 if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
01047 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
01048 {
01049 if (mg_constrained_dofs->set_boundary_values())
01050 {
01051 if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
01052 !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
01053 ||
01054 (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
01055 mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
01056 i1[j] == i2[k]))
01057 G.add(i1[j], i2[k], M(j,k));
01058 }
01059 else
01060 G.add(i1[j], i2[k], M(j,k));
01061 }
01062 }
01063 }
01064
01065 template <class MATRIX>
01066 inline void
01067 MGMatrixSimple<MATRIX>::assemble_out(
01068 MATRIX& G,
01069 const FullMatrix<double>& M,
01070 const std::vector<unsigned int>& i1,
01071 const std::vector<unsigned int>& i2,
01072 const unsigned int level)
01073 {
01074 AssertDimension(M.n(), i1.size());
01075 AssertDimension(M.m(), i2.size());
01076
01077 if(mg_constrained_dofs == 0)
01078 {
01079 for (unsigned int j=0; j<i1.size(); ++j)
01080 for (unsigned int k=0; k<i2.size(); ++k)
01081 if (std::fabs(M(k,j)) >= threshold)
01082 G.add(i1[j], i2[k], M(k,j));
01083 }
01084 else
01085 {
01086 for (unsigned int j=0; j<i1.size(); ++j)
01087 for (unsigned int k=0; k<i2.size(); ++k)
01088 if (std::fabs(M(k,j)) >= threshold)
01089 if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
01090 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
01091 {
01092 if (mg_constrained_dofs->set_boundary_values())
01093 {
01094 if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
01095 !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
01096 ||
01097 (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
01098 mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
01099 i1[j] == i2[k]))
01100 G.add(i1[j], i2[k], M(k,j));
01101 }
01102 else
01103 G.add(i1[j], i2[k], M(k,j));
01104 }
01105 }
01106 }
01107
01108
01109 template <class MATRIX>
01110 template <class DOFINFO>
01111 inline void
01112 MGMatrixSimple<MATRIX>::assemble(const DOFINFO& info)
01113 {
01114 const unsigned int level = info.cell->level();
01115
01116 if (local_indices.size() == 0)
01117 {
01118 assemble((*matrix)[level], info.matrix(0,false).matrix,
01119 info.indices, info.indices, level);
01120 if(mg_constrained_dofs != 0)
01121 {
01122 assemble_in((*interface_in)[level], info.matrix(0,false).matrix,
01123 info.indices, info.indices, level);
01124 assemble_out((*interface_out)[level],info.matrix(0,false).matrix,
01125 info.indices, info.indices, level);
01126 }
01127 }
01128 else
01129 for (unsigned int k=0;k<info.n_matrices();++k)
01130 {
01131 const unsigned int row = info.matrix(k,false).row;
01132 const unsigned int column = info.matrix(k,false).column;
01133
01134 assemble((*matrix)[level], info.matrix(k,false).matrix,
01135 info.indices_by_block[row], info.indices_by_block[column], level);
01136
01137 if(mg_constrained_dofs != 0)
01138 {
01139 assemble_in((*interface_in)[level], info.matrix(k,false).matrix,
01140 info.indices_by_block[row], info.indices_by_block[column], level);
01141 assemble_out((*interface_out)[level],info.matrix(k,false).matrix,
01142 info.indices_by_block[column], info.indices_by_block[row], level);
01143 }
01144 }
01145 }
01146
01147
01148 template <class MATRIX>
01149 template <class DOFINFO>
01150 inline void
01151 MGMatrixSimple<MATRIX>::assemble(const DOFINFO& info1,
01152 const DOFINFO& info2)
01153 {
01154 const unsigned int level1 = info1.cell->level();
01155 const unsigned int level2 = info2.cell->level();
01156
01157 if (local_indices.size() == 0)
01158 {
01159 if (level1 == level2)
01160 {
01161 if(mg_constrained_dofs == 0)
01162 {
01163 assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
01164 assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices);
01165 assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
01166 assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
01167 }
01168 else
01169 {
01170 assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
01171 assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
01172 assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
01173 assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
01174 }
01175 }
01176 else
01177 {
01178 Assert(level1 > level2, ExcInternalError());
01179
01180
01181
01182 assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
01183 if(level1>0)
01184 {
01185 assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
01186 assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
01187 }
01188 }
01189 }
01190 else
01191 for (unsigned int k=0;k<info1.n_matrices();++k)
01192 {
01193 const unsigned int row = info1.matrix(k,false).row;
01194 const unsigned int column = info1.matrix(k,false).column;
01195
01196 if (level1 == level2)
01197 {
01198 if(mg_constrained_dofs == 0)
01199 {
01200 assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
01201 assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column]);
01202 assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column]);
01203 assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column]);
01204 }
01205 else
01206 {
01207 assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
01208 assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
01209 assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
01210 assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
01211 }
01212 }
01213 else
01214 {
01215 Assert(level1 > level2, ExcInternalError());
01216
01217
01218
01219 assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
01220 if(level1>0)
01221 {
01222 assemble_up((*flux_up)[level1],info1.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
01223 assemble_down((*flux_down)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
01224 }
01225 }
01226 }
01227 }
01228
01229
01230
01231 template <class MATRIX, class VECTOR>
01232 SystemSimple<MATRIX,VECTOR>::SystemSimple(double t)
01233 :
01234 MatrixSimple<MATRIX>(t)
01235 {}
01236
01237
01238 template <class MATRIX, class VECTOR>
01239 inline void
01240 SystemSimple<MATRIX,VECTOR>::initialize(MATRIX& m, VECTOR& rhs)
01241 {
01242 NamedData<VECTOR*> data;
01243 VECTOR* p = &rhs;
01244 data.add(p, "right hand side");
01245
01246 MatrixSimple<MATRIX>::initialize(m);
01247 ResidualSimple<VECTOR>::initialize(data);
01248 }
01249
01250
01251 template <class MATRIX, class VECTOR>
01252 template <class DOFINFO>
01253 inline void
01254 SystemSimple<MATRIX,VECTOR>::initialize_info(DOFINFO& info,
01255 bool face) const
01256 {
01257 MatrixSimple<MATRIX>::initialize_info(info, face);
01258 ResidualSimple<VECTOR>::initialize_info(info, face);
01259 }
01260
01261
01262 template <class MATRIX, class VECTOR>
01263 template <class DOFINFO>
01264 inline void
01265 SystemSimple<MATRIX,VECTOR>::assemble(const DOFINFO& info)
01266 {
01267 MatrixSimple<MATRIX>::assemble(info);
01268 ResidualSimple<VECTOR>::assemble(info);
01269 }
01270
01271
01272 template <class MATRIX, class VECTOR>
01273 template <class DOFINFO>
01274 inline void
01275 SystemSimple<MATRIX,VECTOR>::assemble(const DOFINFO& info1,
01276 const DOFINFO& info2)
01277 {
01278 MatrixSimple<MATRIX>::assemble(info1, info2);
01279 ResidualSimple<VECTOR>::assemble(info1, info2);
01280 }
01281 }
01282 }
01283
01284 DEAL_II_NAMESPACE_CLOSE
01285
01286 #endif
01287