00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __deal2__mesh_worker_assembler_h
00014 #define __deal2__mesh_worker_assembler_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/meshworker/simple.h>
00023 #include <deal.II/multigrid/mg_constrained_dofs.h>
00024
00025
00026 DEAL_II_NAMESPACE_OPEN
00027
00028 namespace MeshWorker
00029 {
00081 namespace Assembler
00082 {
00104 template <class VECTOR>
00105 class ResidualLocalBlocksToGlobalBlocks
00106 {
00107 public:
00113 void initialize(const BlockInfo* block_info,
00114 NamedData<VECTOR*>& residuals);
00118 void initialize(const ConstraintMatrix& constraints);
00132 template <class DOFINFO>
00133 void initialize_info(DOFINFO& info, bool face) const;
00134
00135
00140 template<class DOFINFO>
00141 void assemble(const DOFINFO& info);
00142
00147 template<class DOFINFO>
00148 void assemble(const DOFINFO& info1,
00149 const DOFINFO& info2);
00150 private:
00155 void assemble(VECTOR& global,
00156 const BlockVector<double>& local,
00157 const std::vector<unsigned int>& dof);
00158
00164 NamedData<SmartPointer<VECTOR,
00165 ResidualLocalBlocksToGlobalBlocks<VECTOR> > > residuals;
00166
00170 SmartPointer<const BlockInfo,
00171 ResidualLocalBlocksToGlobalBlocks<VECTOR> > block_info;
00175 SmartPointer<const ConstraintMatrix,
00176 ResidualLocalBlocksToGlobalBlocks<VECTOR> > constraints;
00177 };
00178
00179
00207 template <class MATRIX, typename number = double>
00208 class MatrixLocalBlocksToGlobalBlocks
00209 {
00210 public:
00218 MatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
00219
00226 void initialize(const BlockInfo* block_info,
00227 MatrixBlockVector<MATRIX>& matrices);
00228
00232 void initialize(const ConstraintMatrix& constraints);
00246 template <class DOFINFO>
00247 void initialize_info(DOFINFO& info, bool face) const;
00248
00249
00254 template<class DOFINFO>
00255 void assemble(const DOFINFO& info);
00256
00261 template<class DOFINFO>
00262 void assemble(const DOFINFO& info1,
00263 const DOFINFO& info2);
00264
00265 private:
00270 void assemble(
00271 MatrixBlock<MATRIX>& global,
00272 const FullMatrix<number>& local,
00273 unsigned int block_row,
00274 unsigned int block_col,
00275 const std::vector<unsigned int>& dof1,
00276 const std::vector<unsigned int>& dof2);
00277
00283 SmartPointer<MatrixBlockVector<MATRIX>,
00284 MatrixLocalBlocksToGlobalBlocks<MATRIX, number> > matrices;
00285
00289 SmartPointer<const BlockInfo,
00290 MatrixLocalBlocksToGlobalBlocks<MATRIX, number> > block_info;
00294 SmartPointer<const ConstraintMatrix,
00295 MatrixLocalBlocksToGlobalBlocks<MATRIX,number> > constraints;
00296
00306 const double threshold;
00307
00308 };
00309
00334 template <class MATRIX, typename number = double>
00335 class MGMatrixLocalBlocksToGlobalBlocks
00336 {
00337 public:
00338 typedef MGMatrixBlockVector<MATRIX> MatrixPtrVector;
00339 typedef SmartPointer<MatrixPtrVector, MGMatrixLocalBlocksToGlobalBlocks<MATRIX,number> >
00340 MatrixPtrVectorPtr;
00341
00349 MGMatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
00350
00357 void initialize(const BlockInfo* block_info,
00358 MatrixPtrVector& matrices);
00359
00364 void initialize(const MGConstrainedDoFs& mg_constrained_dofs);
00365
00377 void initialize_edge_flux(MatrixPtrVector& up, MatrixPtrVector& down);
00378
00390 void initialize_interfaces (MatrixPtrVector& interface_in, MatrixPtrVector& interface_out);
00404 template <class DOFINFO>
00405 void initialize_info(DOFINFO& info, bool face) const;
00406
00407
00412 template<class DOFINFO>
00413 void assemble(const DOFINFO& info);
00414
00419 template<class DOFINFO>
00420 void assemble(const DOFINFO& info1,
00421 const DOFINFO& info2);
00422
00423 private:
00428 void assemble(
00429 MATRIX& global,
00430 const FullMatrix<number>& local,
00431 unsigned int block_row,
00432 unsigned int block_col,
00433 const std::vector<unsigned int>& dof1,
00434 const std::vector<unsigned int>& dof2,
00435 unsigned int level1,
00436 unsigned int level2,
00437 bool transpose = false);
00438
00443 void assemble_fluxes(
00444 MATRIX& global,
00445 const FullMatrix<number>& local,
00446 unsigned int block_row,
00447 unsigned int block_col,
00448 const std::vector<unsigned int>& dof1,
00449 const std::vector<unsigned int>& dof2,
00450 unsigned int level1,
00451 unsigned int level2);
00452
00457 void assemble_up(
00458 MATRIX& global,
00459 const FullMatrix<number>& local,
00460 unsigned int block_row,
00461 unsigned int block_col,
00462 const std::vector<unsigned int>& dof1,
00463 const std::vector<unsigned int>& dof2,
00464 unsigned int level1,
00465 unsigned int level2);
00466
00471 void assemble_down(
00472 MATRIX& global,
00473 const FullMatrix<number>& local,
00474 unsigned int block_row,
00475 unsigned int block_col,
00476 const std::vector<unsigned int>& dof1,
00477 const std::vector<unsigned int>& dof2,
00478 unsigned int level1,
00479 unsigned int level2);
00480
00485 void assemble_in(
00486 MATRIX& global,
00487 const FullMatrix<number>& local,
00488 unsigned int block_row,
00489 unsigned int block_col,
00490 const std::vector<unsigned int>& dof1,
00491 const std::vector<unsigned int>& dof2,
00492 unsigned int level1,
00493 unsigned int level2);
00494
00499 void assemble_out(
00500 MATRIX& global,
00501 const FullMatrix<number>& local,
00502 unsigned int block_row,
00503 unsigned int block_col,
00504 const std::vector<unsigned int>& dof1,
00505 const std::vector<unsigned int>& dof2,
00506 unsigned int level1,
00507 unsigned int level2);
00508
00514 MatrixPtrVectorPtr matrices;
00515
00521 MatrixPtrVectorPtr flux_down;
00522
00528 MatrixPtrVectorPtr flux_up;
00529
00535 MatrixPtrVectorPtr interface_out;
00536
00542 MatrixPtrVectorPtr interface_in;
00543
00547 SmartPointer<const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number> > block_info;
00548
00552 SmartPointer<const MGConstrainedDoFs,MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number> > mg_constrained_dofs;
00553
00554
00564 const double threshold;
00565
00566 };
00567
00568
00569
00570 template <class VECTOR>
00571 inline void
00572 ResidualLocalBlocksToGlobalBlocks<VECTOR>::initialize(const BlockInfo* b,
00573 NamedData<VECTOR*>& m)
00574 {
00575 block_info = b;
00576 residuals = m;
00577 }
00578
00579 template <class VECTOR>
00580 inline void
00581 ResidualLocalBlocksToGlobalBlocks<VECTOR>::initialize(
00582 const ConstraintMatrix& c)
00583 {
00584 constraints = &c;
00585 }
00586
00587
00588 template <class VECTOR>
00589 template <class DOFINFO>
00590 inline void
00591 ResidualLocalBlocksToGlobalBlocks<VECTOR>::initialize_info(
00592 DOFINFO& info, bool) const
00593 {
00594 info.initialize_vectors(residuals.size());
00595 }
00596
00597 template <class VECTOR>
00598 inline void
00599 ResidualLocalBlocksToGlobalBlocks<VECTOR>::assemble(
00600 VECTOR& global,
00601 const BlockVector<double>& local,
00602 const std::vector<unsigned int>& dof)
00603 {
00604 if(constraints == 0)
00605 {
00606 for (unsigned int b=0;b<local.n_blocks();++b)
00607 for (unsigned int j=0;j<local.block(b).size();++j)
00608 {
00609
00610
00611
00612
00613
00614
00615
00616
00617 const unsigned int jcell = this->block_info->local().local_to_global(b, j);
00618 global(dof[jcell]) += local.block(b)(j);
00619 }
00620 }
00621 else
00622 constraints->distribute_local_to_global(local, dof, global);
00623 }
00624
00625
00626 template <class VECTOR>
00627 template <class DOFINFO>
00628 inline void
00629 ResidualLocalBlocksToGlobalBlocks<VECTOR>::assemble(
00630 const DOFINFO& info)
00631 {
00632 for (unsigned int i=0;i<residuals.size();++i)
00633 assemble(*residuals(i), info.vector(i), info.indices);
00634 }
00635
00636
00637 template <class VECTOR>
00638 template <class DOFINFO>
00639 inline void
00640 ResidualLocalBlocksToGlobalBlocks<VECTOR>::assemble(
00641 const DOFINFO& info1,
00642 const DOFINFO& info2)
00643 {
00644 for (unsigned int i=0;i<residuals.size();++i)
00645 {
00646 assemble(*residuals(i), info1.vector(i), info1.indices);
00647 assemble(*residuals(i), info2.vector(i), info2.indices);
00648 }
00649 }
00650
00651
00652
00653
00654 template <class MATRIX, typename number>
00655 inline
00656 MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::MatrixLocalBlocksToGlobalBlocks(
00657 double threshold)
00658 :
00659 threshold(threshold)
00660 {}
00661
00662
00663 template <class MATRIX, typename number>
00664 inline void
00665 MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(
00666 const BlockInfo* b,
00667 MatrixBlockVector<MATRIX>& m)
00668 {
00669 block_info = b;
00670 matrices = &m;
00671 }
00672
00673
00674
00675 template <class MATRIX, typename number>
00676 inline void
00677 MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(
00678 const ConstraintMatrix& c)
00679 {
00680 constraints = &c;
00681 }
00682
00683
00684
00685 template <class MATRIX ,typename number>
00686 template <class DOFINFO>
00687 inline void
00688 MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize_info(
00689 DOFINFO& info,
00690 bool face) const
00691 {
00692 info.initialize_matrices(*matrices, face);
00693 }
00694
00695
00696
00697 template <class MATRIX, typename number>
00698 inline void
00699 MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
00700 MatrixBlock<MATRIX>& global,
00701 const FullMatrix<number>& local,
00702 unsigned int block_row,
00703 unsigned int block_col,
00704 const std::vector<unsigned int>& dof1,
00705 const std::vector<unsigned int>& dof2)
00706 {
00707 if(constraints == 0)
00708 {
00709 for (unsigned int j=0;j<local.n_rows();++j)
00710 for (unsigned int k=0;k<local.n_cols();++k)
00711 if (std::fabs(local(j,k)) >= threshold)
00712 {
00713
00714
00715
00716
00717
00718
00719
00720
00721 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
00722 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
00723
00724 global.add(dof1[jcell], dof2[kcell], local(j,k));
00725 }
00726 }
00727 else
00728 {
00729 const BlockIndices &bi = this->block_info->local();
00730 std::vector<unsigned int> sliced_row_indices (bi.block_size(block_row));
00731 for(unsigned int i=0; i<sliced_row_indices.size(); ++i)
00732 sliced_row_indices[i] = dof1[bi.block_start(block_row)+i];
00733
00734 std::vector<unsigned int> sliced_col_indices (bi.block_size(block_col));
00735 for(unsigned int i=0; i<sliced_col_indices.size(); ++i)
00736 sliced_col_indices[i] = dof2[bi.block_start(block_col)+i];
00737
00738 constraints->distribute_local_to_global(local,
00739 sliced_row_indices, sliced_col_indices, global);
00740 }
00741 }
00742
00743
00744 template <class MATRIX, typename number>
00745 template <class DOFINFO>
00746 inline void
00747 MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
00748 const DOFINFO& info)
00749 {
00750 for (unsigned int i=0;i<matrices->size();++i)
00751 {
00752
00753
00754 const unsigned int row = matrices->block(i).row;
00755 const unsigned int col = matrices->block(i).column;
00756
00757 assemble(matrices->block(i), info.matrix(i,false).matrix, row, col, info.indices, info.indices);
00758 }
00759 }
00760
00761
00762 template <class MATRIX, typename number>
00763 template <class DOFINFO>
00764 inline void
00765 MatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
00766 const DOFINFO& info1,
00767 const DOFINFO& info2)
00768 {
00769 for (unsigned int i=0;i<matrices->size();++i)
00770 {
00771
00772
00773 const unsigned int row = matrices->block(i).row;
00774 const unsigned int col = matrices->block(i).column;
00775
00776 assemble(matrices->block(i), info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices);
00777 assemble(matrices->block(i), info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices);
00778 assemble(matrices->block(i), info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices);
00779 assemble(matrices->block(i), info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices);
00780 }
00781 }
00782
00783
00784
00785
00786 template <class MATRIX, typename number>
00787 inline
00788 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::MGMatrixLocalBlocksToGlobalBlocks(
00789 double threshold)
00790 :
00791 threshold(threshold)
00792 {}
00793
00794
00795 template <class MATRIX, typename number>
00796 inline void
00797 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(
00798 const BlockInfo* b,
00799 MatrixPtrVector& m)
00800 {
00801 block_info = b;
00802 AssertDimension(block_info->local().size(), block_info->global().size());
00803 matrices = &m;
00804 }
00805
00806
00807 template <class MATRIX, typename number>
00808 inline void
00809 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize(
00810 const MGConstrainedDoFs& mg_c)
00811 {
00812 mg_constrained_dofs = &mg_c;
00813 }
00814
00815
00816 template <class MATRIX ,typename number>
00817 template <class DOFINFO>
00818 inline void
00819 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize_info(
00820 DOFINFO& info,
00821 bool face) const
00822 {
00823 info.initialize_matrices(*matrices, face);
00824 }
00825
00826
00827
00828 template <class MATRIX, typename number>
00829 inline void
00830 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize_edge_flux(
00831 MatrixPtrVector& up,
00832 MatrixPtrVector& down)
00833 {
00834 flux_up = up;
00835 flux_down = down;
00836 }
00837
00838
00839 template <class MATRIX, typename number>
00840 inline void
00841 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::initialize_interfaces(
00842 MatrixPtrVector& in,
00843 MatrixPtrVector& out)
00844 {
00845 interface_in = in;
00846 interface_out = out;
00847 }
00848
00849
00850 template <class MATRIX, typename number>
00851 inline void
00852 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
00853 MATRIX& global,
00854 const FullMatrix<number>& local,
00855 unsigned int block_row,
00856 unsigned int block_col,
00857 const std::vector<unsigned int>& dof1,
00858 const std::vector<unsigned int>& dof2,
00859 unsigned int level1,
00860 unsigned int level2,
00861 bool transpose)
00862 {
00863 for (unsigned int j=0;j<local.n_rows();++j)
00864 for (unsigned int k=0;k<local.n_cols();++k)
00865 if (std::fabs(local(j,k)) >= threshold)
00866 {
00867
00868
00869
00870
00871
00872
00873
00874
00875 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
00876 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
00887 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
00888
00889 if(mg_constrained_dofs == 0)
00890 {
00891 if (transpose)
00892 global.add(kglobal, jglobal, local(j,k));
00893 else
00894 global.add(jglobal, kglobal, local(j,k));
00895 }
00896 else
00897 {
00898 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
00899 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
00900 {
00901 if (mg_constrained_dofs->set_boundary_values())
00902 {
00903 if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
00904 !mg_constrained_dofs->is_boundary_index(level2, kglobal))
00905 ||
00906 (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
00907 mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
00908 jglobal == kglobal))
00909 {
00910 if (transpose)
00911 global.add(kglobal, jglobal, local(j,k));
00912 else
00913 global.add(jglobal, kglobal, local(j,k));
00914 }
00915 }
00916 else
00917 {
00918 if (transpose)
00919 global.add(kglobal, jglobal, local(j,k));
00920 else
00921 global.add(jglobal, kglobal, local(j,k));
00922 }
00923 }
00924 }
00925 }
00926 }
00927
00928
00929 template <class MATRIX, typename number>
00930 inline void
00931 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble_fluxes(
00932 MATRIX& global,
00933 const FullMatrix<number>& local,
00934 unsigned int block_row,
00935 unsigned int block_col,
00936 const std::vector<unsigned int>& dof1,
00937 const std::vector<unsigned int>& dof2,
00938 unsigned int level1,
00939 unsigned int level2)
00940 {
00941 for (unsigned int j=0;j<local.n_rows();++j)
00942 for (unsigned int k=0;k<local.n_cols();++k)
00943 if (std::fabs(local(j,k)) >= threshold)
00944 {
00945
00946
00947
00948
00949
00950
00951
00952
00953 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
00954 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
00965 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
00966
00967 if(mg_constrained_dofs == 0)
00968 global.add(jglobal, kglobal, local(j,k));
00969 else
00970 {
00971 if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
00972 !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
00973 {
00974 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
00975 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
00976 global.add(jglobal, kglobal, local(j,k));
00977 }
00978 }
00979 }
00980 }
00981
00982 template <class MATRIX, typename number>
00983 inline void
00984 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble_up(
00985 MATRIX& global,
00986 const FullMatrix<number>& local,
00987 unsigned int block_row,
00988 unsigned int block_col,
00989 const std::vector<unsigned int>& dof1,
00990 const std::vector<unsigned int>& dof2,
00991 unsigned int level1,
00992 unsigned int level2)
00993 {
00994 for (unsigned int j=0;j<local.n_rows();++j)
00995 for (unsigned int k=0;k<local.n_cols();++k)
00996 if (std::fabs(local(j,k)) >= threshold)
00997 {
00998
00999
01000
01001
01002
01003
01004
01005
01006 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
01007 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
01018 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
01019
01020 if(mg_constrained_dofs == 0)
01021 global.add(jglobal, kglobal, local(j,k));
01022 else
01023 {
01024 if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
01025 !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
01026 {
01027 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
01028 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
01029 global.add(jglobal, kglobal, local(j,k));
01030 }
01031 }
01032 }
01033 }
01034
01035 template <class MATRIX, typename number>
01036 inline void
01037 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble_down(
01038 MATRIX& global,
01039 const FullMatrix<number>& local,
01040 unsigned int block_row,
01041 unsigned int block_col,
01042 const std::vector<unsigned int>& dof1,
01043 const std::vector<unsigned int>& dof2,
01044 unsigned int level1,
01045 unsigned int level2)
01046 {
01047 for (unsigned int j=0;j<local.n_rows();++j)
01048 for (unsigned int k=0;k<local.n_cols();++k)
01049 if (std::fabs(local(k,j)) >= threshold)
01050 {
01051
01052
01053
01054
01055
01056
01057
01058
01059 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
01060 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
01071 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
01072
01073 if(mg_constrained_dofs == 0)
01074 global.add(jglobal, kglobal, local(k,j));
01075 else
01076 {
01077 if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
01078 !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
01079 {
01080 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
01081 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
01082 global.add(jglobal, kglobal, local(k,j));
01083 }
01084 }
01085 }
01086 }
01087
01088 template <class MATRIX, typename number>
01089 inline void
01090 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble_in(
01091 MATRIX& global,
01092 const FullMatrix<number>& local,
01093 unsigned int block_row,
01094 unsigned int block_col,
01095 const std::vector<unsigned int>& dof1,
01096 const std::vector<unsigned int>& dof2,
01097 unsigned int level1,
01098 unsigned int level2)
01099 {
01100
01101
01102
01103 for (unsigned int j=0;j<local.n_rows();++j)
01104 for (unsigned int k=0;k<local.n_cols();++k)
01105 if (std::fabs(local(j,k)) >= threshold)
01106 {
01107
01108
01109
01110
01111
01112
01113
01114
01115 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
01116 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
01127 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
01128
01129 if(mg_constrained_dofs == 0)
01130 global.add(jglobal, kglobal, local(j,k));
01131 else
01132 {
01133 if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
01134 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
01135 {
01136 if (mg_constrained_dofs->set_boundary_values())
01137 {
01138 if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
01139 !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
01140 ||
01141 (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
01142 mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
01143 jglobal == kglobal))
01144 global.add(jglobal, kglobal, local(j,k));
01145 }
01146 else
01147 global.add(jglobal, kglobal, local(j,k));
01148 }
01149 }
01150 }
01151 }
01152
01153 template <class MATRIX, typename number>
01154 inline void
01155 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble_out(
01156 MATRIX& global,
01157 const FullMatrix<number>& local,
01158 unsigned int block_row,
01159 unsigned int block_col,
01160 const std::vector<unsigned int>& dof1,
01161 const std::vector<unsigned int>& dof2,
01162 unsigned int level1,
01163 unsigned int level2)
01164 {
01165
01166
01167
01168 for (unsigned int j=0;j<local.n_rows();++j)
01169 for (unsigned int k=0;k<local.n_cols();++k)
01170 if (std::fabs(local(k,j)) >= threshold)
01171 {
01172
01173
01174
01175
01176
01177
01178
01179
01180 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
01181 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
01192 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
01193
01194 if(mg_constrained_dofs == 0)
01195 global.add(jglobal, kglobal, local(k,j));
01196 else
01197 {
01198 if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
01199 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
01200 {
01201 if (mg_constrained_dofs->set_boundary_values())
01202 {
01203 if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
01204 !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
01205 ||
01206 (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
01207 mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
01208 jglobal == kglobal))
01209 global.add(jglobal, kglobal, local(k,j));
01210 }
01211 else
01212 global.add(jglobal, kglobal, local(k,j));
01213 }
01214 }
01215 }
01216 }
01217
01218
01219 template <class MATRIX, typename number>
01220 template <class DOFINFO>
01221 inline void
01222 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(const DOFINFO& info)
01223 {
01224 const unsigned int level = info.cell->level();
01225
01226 for (unsigned int i=0;i<matrices->size();++i)
01227 {
01228
01229
01230 const unsigned int row = matrices->block(i)[level].row;
01231 const unsigned int col = matrices->block(i)[level].column;
01232
01233 assemble(matrices->block(i)[level].matrix, info.matrix(i,false).matrix, row, col,
01234 info.indices, info.indices, level, level);
01235 if(mg_constrained_dofs != 0)
01236 {
01237 if(interface_in != 0)
01238 assemble_in(interface_in->block(i)[level], info.matrix(i,false).matrix, row, col,
01239 info.indices, info.indices, level, level);
01240 if(interface_out != 0)
01241 assemble_in(interface_out->block(i)[level], info.matrix(i,false).matrix, row, col,
01242 info.indices, info.indices, level, level);
01243
01244 assemble_in(matrices->block_in(i)[level], info.matrix(i,false).matrix, row, col,
01245 info.indices, info.indices, level, level);
01246 assemble_out(matrices->block_out(i)[level], info.matrix(i,false).matrix, row, col,
01247 info.indices, info.indices, level, level);
01248 }
01249 }
01250 }
01251
01252
01253 template <class MATRIX, typename number>
01254 template <class DOFINFO>
01255 inline void
01256 MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number>::assemble(
01257 const DOFINFO& info1,
01258 const DOFINFO& info2)
01259 {
01260 const unsigned int level1 = info1.cell->level();
01261 const unsigned int level2 = info2.cell->level();
01262
01263 for (unsigned int i=0;i<matrices->size();++i)
01264 {
01265 MGLevelObject<MatrixBlock<MATRIX> >& o = matrices->block(i);
01266
01267
01268
01269 const unsigned int row = o[level1].row;
01270 const unsigned int col = o[level1].column;
01271
01272 if (level1 == level2)
01273 {
01274 if(mg_constrained_dofs == 0)
01275 {
01276 assemble(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
01277 info1.indices, info1.indices, level1, level1);
01278 assemble(o[level1].matrix, info1.matrix(i,true).matrix, row, col,
01279 info1.indices, info2.indices, level1, level2);
01280 assemble(o[level1].matrix, info2.matrix(i,false).matrix, row, col,
01281 info2.indices, info2.indices, level2, level2);
01282 assemble(o[level1].matrix, info2.matrix(i,true).matrix, row, col,
01283 info2.indices, info1.indices, level2, level1);
01284 }
01285 else
01286 {
01287 assemble_fluxes(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
01288 info1.indices, info1.indices, level1, level1);
01289 assemble_fluxes(o[level1].matrix, info1.matrix(i,true).matrix, row, col,
01290 info1.indices, info2.indices, level1, level2);
01291 assemble_fluxes(o[level1].matrix, info2.matrix(i,false).matrix, row, col,
01292 info2.indices, info2.indices, level2, level2);
01293 assemble_fluxes(o[level1].matrix, info2.matrix(i,true).matrix, row, col,
01294 info2.indices, info1.indices, level2, level1);
01295 }
01296 }
01297 else
01298 {
01299 Assert(level1 > level2, ExcNotImplemented());
01300 if (flux_up->size() != 0)
01301 {
01302
01303
01304
01305 assemble_fluxes(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
01306 info1.indices, info1.indices, level1, level1);
01307 assemble_up(flux_up->block(i)[level1].matrix, info1.matrix(i,true).matrix, row, col,
01308 info1.indices, info2.indices, level1, level2);
01309 assemble_down(flux_down->block(i)[level1].matrix, info2.matrix(i,true).matrix, row, col,
01310 info2.indices, info1.indices, level2, level1);
01311 }
01312 }
01313 }
01314 }
01315 }
01316 }
01317
01318 DEAL_II_NAMESPACE_CLOSE
01319
01320 #endif
01321