include/deal.II/meshworker/assembler.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: assembler.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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               // The coordinates of
00610               // the current entry in
00611               // DoFHandler
00612               // numbering, which
00613               // differs from the
00614               // block-wise local
00615               // numbering we use in
00616               // our local vectors
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               // The coordinates of
00714               // the current entry in
00715               // DoFHandler
00716               // numbering, which
00717               // differs from the
00718               // block-wise local
00719               // numbering we use in
00720               // our local matrices
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                                            // Row and column index of
00753                                            // the block we are dealing with
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                                            // Row and column index of
00772                                            // the block we are dealing with
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                                                // The coordinates of
00868                                                // the current entry in
00869                                                // DoFHandler
00870                                                // numbering, which
00871                                                // differs from the
00872                                                // block-wise local
00873                                                // numbering we use in
00874                                                // our local matrices
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                                                // The global dof
00879                                                // indices to assemble
00880                                                // in. Since we may
00881                                                // have face matrices
00882                                                // coupling two
00883                                                // different cells, we
00884                                                // provide two sets of
00885                                                // dof indices.
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                                                // The coordinates of
00946                                                // the current entry in
00947                                                // DoFHandler
00948                                                // numbering, which
00949                                                // differs from the
00950                                                // block-wise local
00951                                                // numbering we use in
00952                                                // our local matrices
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                                                // The global dof
00957                                                // indices to assemble
00958                                                // in. Since we may
00959                                                // have face matrices
00960                                                // coupling two
00961                                                // different cells, we
00962                                                // provide two sets of
00963                                                // dof indices.
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                                                // The coordinates of
00999                                                // the current entry in
01000                                                // DoFHandler
01001                                                // numbering, which
01002                                                // differs from the
01003                                                // block-wise local
01004                                                // numbering we use in
01005                                                // our local matrices
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                                                // The global dof
01010                                                // indices to assemble
01011                                                // in. Since we may
01012                                                // have face matrices
01013                                                // coupling two
01014                                                // different cells, we
01015                                                // provide two sets of
01016                                                // dof indices.
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                                                // The coordinates of
01052                                                // the current entry in
01053                                                // DoFHandler
01054                                                // numbering, which
01055                                                // differs from the
01056                                                // block-wise local
01057                                                // numbering we use in
01058                                                // our local matrices
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                                                // The global dof
01063                                                // indices to assemble
01064                                                // in. Since we may
01065                                                // have face matrices
01066                                                // coupling two
01067                                                // different cells, we
01068                                                // provide two sets of
01069                                                // dof indices.
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 //      AssertDimension(local.n(), dof1.size());
01101 //      AssertDimension(local.m(), dof2.size());
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                                                // The coordinates of
01108                                                // the current entry in
01109                                                // DoFHandler
01110                                                // numbering, which
01111                                                // differs from the
01112                                                // block-wise local
01113                                                // numbering we use in
01114                                                // our local matrices
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                                                // The global dof
01119                                                // indices to assemble
01120                                                // in. Since we may
01121                                                // have face matrices
01122                                                // coupling two
01123                                                // different cells, we
01124                                                // provide two sets of
01125                                                // dof indices.
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 //      AssertDimension(local.n(), dof1.size());
01166 //      AssertDimension(local.m(), dof2.size());
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                                                // The coordinates of
01173                                                // the current entry in
01174                                                // DoFHandler
01175                                                // numbering, which
01176                                                // differs from the
01177                                                // block-wise local
01178                                                // numbering we use in
01179                                                // our local matrices
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                                                // The global dof
01184                                                // indices to assemble
01185                                                // in. Since we may
01186                                                // have face matrices
01187                                                // coupling two
01188                                                // different cells, we
01189                                                // provide two sets of
01190                                                // dof indices.
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                                            // Row and column index of
01229                                            // the block we are dealing with
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                                            // Row and column index of
01268                                            // the block we are dealing with
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                                                    // Do not add M22,
01303                                                    // which is done by
01304                                                    // the coarser cell
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:46 by doxygen 1.7.3