NGSolve  4.9
multigrid/prolongation.hpp
00001 #ifndef FILE_PROLONGATION
00002 #define FILE_PROLONGATION
00003 
00004 /*********************************************************************/
00005 /* File:   prolongation.hh                                           */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   20. Apr. 2000                                             */
00008 /*********************************************************************/
00009 
00010 namespace ngmg
00011 {
00012 
00016   class NGS_DLL_HEADER Prolongation
00017   {
00018   public:
00020     Prolongation();
00022     virtual ~Prolongation();
00023   
00025     virtual void Update () = 0;
00026 
00028     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const = 0;
00030     virtual void ProlongateInline (int finelevel, BaseVector & v) const = 0;
00032     virtual void RestrictInline (int finelevel, BaseVector & v) const = 0;
00033 
00034     virtual BitArray * GetInnerDofs () const { return 0; }
00035   };
00036 
00037 
00042   // template <class TV>
00043   class LinearProlongation : public Prolongation
00044   {
00046     const MeshAccess & ma;
00048     const FESpace & space;
00050     Array<int> nvlevel;
00051   public:
00053     LinearProlongation(const MeshAccess & ama,
00054                        const FESpace & aspace)
00055       : ma(ama), space(aspace) { ; }
00057     LinearProlongation(const FESpace & aspace)
00058       : ma(aspace.GetMeshAccess()), space(aspace) { ; }
00059 
00061     virtual ~LinearProlongation() { ; }
00062   
00064     virtual void Update () 
00065     { 
00066       if (ma.GetNLevels() > nvlevel.Size())
00067         nvlevel.Append (ma.GetNV());
00068     }
00069 
00070 
00072     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const;
00073 
00074 
00076     virtual void ProlongateInline (int finelevel, BaseVector & v) const
00077     {
00078       int parents[2];
00079 
00080       int nc = nvlevel[finelevel-1];
00081       int nf = nvlevel[finelevel];
00082     
00083       //    FlatVector<TV> & fv = 
00084       //      dynamic_cast<VFlatVector<TV> &> (v).FV();
00085       //    FlatVector<TV> fv = 
00086       //      dynamic_cast<T_BaseVector<TV> &> (v).FV();
00087 
00088       FlatSysVector<> fv (v.Size(), v.EntrySize(), static_cast<double*>(v.Memory()));
00089 
00090 
00091       int i;
00092       for (i = nf; i < fv.Size(); i++)
00093         fv(i) = 0;
00094 
00095       for (i = nc; i < nf; i++)
00096         {
00097           ma.GetParentNodes (i, parents);
00098           fv(i) = 0.5 * (fv(parents[0]) + fv(parents[1]));
00099         }
00100     }
00101 
00102 
00104     virtual void RestrictInline (int finelevel, BaseVector & v) const
00105     {
00106       int parents[2];
00107       // int nc = space.GetNDofLevel (finelevel-1);
00108       // int nf = space.GetNDofLevel (finelevel);
00109       int nc = nvlevel[finelevel-1];
00110       int nf = nvlevel[finelevel];
00111 
00112       /*
00113         cout << "rest, h1, typeid(v) = " << typeid(v).name() << endl;
00114         cout << "nvlevel = " << nvlevel << ", level = " << finelevel << endl;
00115         cout << "nc = " << nc << ", nf = " << nf << endl;
00116         cout << "v.size = " << v.Size() << ", entrysize = " << v.EntrySize() << endl;
00117       */
00118 
00119       // FlatVector<TV> fv = dynamic_cast<T_BaseVector<TV> &> (v).FV();
00120 
00121       FlatSysVector<> fv (v.Size(), v.EntrySize(), static_cast<double*>(v.Memory()));
00122 
00123       int i;
00124       for (i = nf-1; i >= nc; i--)
00125         {
00126           ma.GetParentNodes (i, parents);
00127           fv(parents[0]) += 0.5 * fv(i);
00128           fv(parents[1]) += 0.5 * fv(i);
00129         }
00130 
00131       for (i = nf; i < fv.Size(); i++)
00132         fv(i) = 0;  
00133     }
00134   };
00135 
00136 
00137   /*
00139   class NonConformingProlongation : public Prolongation
00140   {
00142     const MeshAccess & ma;
00144     const NonConformingFESpace & space;
00145   public:
00147     NonConformingProlongation(const MeshAccess & ama,
00148                               const NonConformingFESpace & aspace);
00150     virtual ~NonConformingProlongation();
00151   
00153     virtual void Update ();
00154 
00156     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const
00157     { return NULL; }
00159     virtual void ProlongateInline (int finelevel, BaseVector & v) const;
00161     virtual void RestrictInline (int finelevel, BaseVector & v) const;
00162   };
00163   */
00164 
00165 
00167   class ElementProlongation : public Prolongation
00168   {
00170     const MeshAccess & ma;
00172     const ElementFESpace & space;
00173   public:
00175     ElementProlongation(const ElementFESpace & aspace)
00176       : ma(aspace.GetMeshAccess()), space(aspace) { ; }
00178     virtual ~ElementProlongation()
00179     { ; }
00180   
00182     virtual void Update ()
00183     { ; }
00184 
00186     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const
00187     { return NULL; }
00188 
00190     virtual void ProlongateInline (int finelevel, BaseVector & v) const
00191     {
00192       // FlatVector<TV> fv = dynamic_cast<T_BaseVector<TV> &> (v).FV();    
00193 
00194       FlatSysVector<> fv (v.Size(), v.EntrySize(), static_cast<double*>(v.Memory()));
00195     
00196       int nc = space.GetNDofLevel (finelevel-1);
00197       int nf = space.GetNDofLevel (finelevel);
00198 
00199       for (int i = nc; i < nf; i++)
00200         {
00201           int parent = ma.GetParentElement (i);
00202           fv(i) = fv(parent);
00203         }
00204     
00205       for (int i = nf; i < fv.Size(); i++)
00206         fv(i) = 0;
00207     }
00208 
00210     virtual void RestrictInline (int finelevel, BaseVector & v) const
00211     {
00212       //    FlatVector<TV> fv = dynamic_cast<T_BaseVector<TV> &> (v).FV();    
00213 
00214       FlatSysVector<> fv (v.Size(), v.EntrySize(), static_cast<double*>(v.Memory()));
00215     
00216       int nc = space.GetNDofLevel (finelevel-1);
00217       int nf = space.GetNDofLevel (finelevel);
00218 
00219       for (int i = nf-1; i >= nc; i--)
00220         {
00221           int parent = ma.GetParentElement (i);
00222           fv(parent) += fv(i);
00223           fv(i) = 0;
00224         }
00225     }
00226   };
00227 
00228 
00229 
00231   class SurfaceElementProlongation : public Prolongation
00232   {
00234     const MeshAccess & ma;
00236     const SurfaceElementFESpace & space;
00237   public:
00239     SurfaceElementProlongation(const MeshAccess & ama,
00240                                const SurfaceElementFESpace & aspace);
00242     virtual ~SurfaceElementProlongation();
00243   
00245     virtual void Update ();
00246 
00248     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const
00249     { return NULL; }
00251     virtual void ProlongateInline (int finelevel, BaseVector & v) const;
00253     virtual void RestrictInline (int finelevel, BaseVector & v) const;
00254   };
00255 
00256 
00257 
00259                                      // template <class TV>
00260   class EdgeProlongation : public Prolongation
00261   {
00263     const MeshAccess & ma;
00265     const NedelecFESpace & space;
00266   public:
00268     EdgeProlongation(const NedelecFESpace & aspace)
00269       : ma(aspace.GetMeshAccess()), space(aspace) { ; }
00271     virtual ~EdgeProlongation() { ; }
00272   
00274     virtual void Update () { ; }
00275 
00277     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const
00278     { return NULL; }
00279 
00281     virtual void ProlongateInline (int finelevel, BaseVector & v) const
00282     {
00283       int nc = space.GetNDofLevel (finelevel-1);
00284       int nf = space.GetNDofLevel (finelevel);
00285       /*    
00286             FlatVector<TV> & fv = 
00287             dynamic_cast<VFlatVector<TV> &> (v).FV();
00288       */
00289       //    FlatVector<TV> fv = dynamic_cast<T_BaseVector<TV> &> (v).FV();
00290       FlatSysVector<> fv (v.Size(), v.EntrySize(), static_cast<double*>(v.Memory()));
00291 
00292       int i, k;
00293 
00294       for (i = nf; i < fv.Size(); i++)
00295         fv(i) = 0;
00296 
00297       for (k = 1; k <= 10; k++)
00298         for (i = nc; i < nf; i++)
00299           {
00300             int pa1 = space.ParentEdge1 (i);
00301             int pa2 = space.ParentEdge2 (i);
00302           
00303             fv(i) = 0;
00304             if (pa1 != -1)
00305               {
00306                 if (pa1 & 1)
00307                   fv(i) += 0.5 * fv(pa1/2);
00308                 else
00309                   fv(i) -= 0.5 * fv(pa1/2);
00310               }
00311             if (pa2 != -1)
00312               {
00313                 if (pa2 & 1)
00314                   fv(i) += 0.5 * fv(pa2/2);
00315                 else
00316                   fv(i) -= 0.5 * fv(pa2/2);
00317               }
00318           }
00319 
00320       for (i = 0; i < nf; i++)
00321         if (space.FineLevelOfEdge(i) < finelevel)
00322           fv(i) = 0;
00323     }
00324 
00325 
00327     virtual void RestrictInline (int finelevel, BaseVector & v) const
00328     {
00329       int nc = space.GetNDofLevel (finelevel-1);
00330       int nf = space.GetNDofLevel (finelevel);
00331 
00332       //    FlatVector<TV> & fv = 
00333       //      dynamic_cast<VFlatVector<TV> &> (v).FV();
00334       //    FlatVector<TV> fv = dynamic_cast<T_BaseVector<TV> &> (v).FV();
00335 
00336       FlatSysVector<> fv (v.Size(), v.EntrySize(), static_cast<double*>(v.Memory()));
00337 
00338       for (int i = 0; i < nf; i++)
00339         if (space.FineLevelOfEdge(i) < finelevel)
00340           fv(i) = 0;
00341         
00342       for (int k = 1; k <= 10; k++)
00343         for (int i = nf-1; i >= nc; i--)
00344           {
00345             int pa1 = space.ParentEdge1 (i);
00346             int pa2 = space.ParentEdge2 (i);
00347           
00348             if (pa1 != -1)
00349               {
00350                 if (pa1 & 1)
00351                   fv(pa1/2) += 0.5 * fv(i);
00352                 else
00353                   fv(pa1/2) -= 0.5 * fv(i);
00354               }
00355             if (pa2 != -1)
00356               {
00357                 if (pa2 & 1)
00358                   fv(pa2/2) += 0.5 * fv(i);
00359                 else
00360                   fv(pa2/2) -= 0.5 * fv(i);
00361               }
00362             fv(i) = 0;
00363           }
00364 
00365       for (int i = nf; i < fv.Size(); i++)
00366         fv(i) = 0;  
00367     }
00368 
00370     void ApplyGradient (int level, const BaseVector & pot, BaseVector & grad) const
00371     {
00372       cout << "apply grad" << endl;
00373     }
00374   };
00375 
00376 
00377 
00379   class L2HoProlongation : public Prolongation
00380   {
00382     const MeshAccess & ma;
00384     const Array<int> & first_dofs;
00385   public:
00387     L2HoProlongation(const MeshAccess & ama, const Array<int> & afirst_dofs);
00389     virtual ~L2HoProlongation()
00390     { ; }
00392     virtual void Update ()
00393         { ; }
00394 
00396     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const
00397     { return NULL; }
00398 
00400     virtual void ProlongateInline (int finelevel, BaseVector & v) const;
00401 
00403     virtual void RestrictInline (int finelevel, BaseVector & v) const
00404     {
00405                 cout << "RestrictInline not implemented for L2HoProlongation" << endl;
00406     }
00407  
00408   };
00409 
00410 
00412   class CompoundProlongation : public Prolongation
00413   {
00414   protected:
00416     const CompoundFESpace * space;
00418     Array<const Prolongation*> prols;
00419   public:
00421     CompoundProlongation(const CompoundFESpace * aspace);
00423     CompoundProlongation(const CompoundFESpace * aspace,
00424                          Array<const Prolongation*> & aprols);
00426     virtual ~CompoundProlongation();
00427     // { ; }
00428   
00430     virtual void Update ();
00431 
00432     void AddProlongation (const Prolongation * prol)
00433     {
00434       prols.Append (prol);
00435     }
00436 
00438     virtual SparseMatrix< double >* CreateProlongationMatrix( int finelevel ) const
00439     { return NULL; }
00440 
00441 
00443     virtual void ProlongateInline (int finelevel, BaseVector & v) const;
00444 
00446     virtual void RestrictInline (int finelevel, BaseVector & v) const;
00447   };
00448 }
00449 
00450 
00451 #endif