NGSolve
4.9
|
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