include/deal.II/meshworker/simple.h

Go to the documentation of this file.
00001 //---------------------------------------------------------------------------
00002 //    @f$Id: simple.h 25370 2012-04-02 20:53:57Z kanschat @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_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 = &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                                                          // At the
00939                                                          // boundary,
00940                                                          // only enter
00941                                                          // the term
00942                                                          // on the
00943                                                          // diagonal,
00944                                                          // but not
00945                                                          // the
00946                                                          // coupling terms
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                                                // Do not add info2.M1,
01180                                                // which is done by
01181                                                // the coarser cell
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                                                // Do not add info2.M1,
01217                                                // which is done by
01218                                                // the coarser cell
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:42 by doxygen 1.7.3