vtf-logo

GFMBoundary.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, ralf@amroc.net
00005 
00006 #ifndef AMROC_GFMBOUNDARYBASE_H
00007 #define AMROC_GFMBOUNDARYBASE_H
00008 
00016 #define REFLECTION     2
00017 #define EXTRAPOLATION  1
00018 
00025 #include "GFMGeom.h"
00026 #include "Interpolation.h"
00027 
00028 template <class VectorType, int dim>
00029 class GFMBoundaryBase : public AMRBase<VectorType,dim>, 
00030                         public Geom<typename VectorType::InternalDataType,dim>,
00031                         public GFMGeom<typename VectorType::InternalDataType,dim> {
00032   typedef typename VectorType::InternalDataType DataType;
00033   
00034   typedef AMRBase<VectorType,dim> base;
00035   typedef Geom<DataType,dim> geom_base;
00036   typedef GFMGeom<DataType,dim> gfm_geom_base;
00037 public:
00038   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00039   typedef typename base::vec_grid_data_type vec_grid_data_type;
00040 
00041   typedef GridData<DataType,dim> grid_data_type;   
00042   typedef GridFunction<DataType,dim> grid_fct_type;  
00043   typedef GridData<bool,dim> bool_grid_data_type;
00044   typedef GridFunction<bool,dim> bool_grid_fct_type;  
00045   typedef Vector<DataType,dim> point_type;
00046 
00047   typedef Interpolation<VectorType,dim> interpolation_type;
00048 
00049   GFMBoundaryBase() : base(), geom_base(), _BoundaryWidth(1), _bndtype(2), _naux(0), _Verbose(1),
00050                       _Cutoff(0.0), _ErrorValue(-1.e37), _normal_factor(2.0) {
00051     _interpolation = new interpolation_type(); 
00052     _far_away = 1.e36;    
00053   }
00054 
00055   virtual ~GFMBoundaryBase() {
00056     if (_interpolation) delete _interpolation;
00057   }
00058 
00059   //******************************************************************************
00060   // Abstract class interface
00061   //******************************************************************************
00062   virtual void SetGrid(vec_grid_data_type& gdu, grid_data_type& gdphi, 
00063                        const BBox& bb, const int& Level, double t, 
00064                        const int& nc, const int* idx, const point_type* xc, 
00065                        DataType* distance, point_type* normal) = 0;  
00066   virtual void GetBndryCells(vec_grid_fct_type& u, grid_data_type& gdphi, bool_grid_data_type& gdbf, 
00067                              const BBox& bb, const int& Time, const int& Level, const int& c, 
00068                              int& nc, int*& idx, int& ni, int*& idx_wrg) = 0;
00069   virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi, const BBox& bb, 
00070                            const VectorType* u, DataType* aux, const int& Level, 
00071                            double t, const int& nc, const int* idx, 
00072                            const point_type* xc, const DataType* distance, 
00073                            const point_type* normal) {}
00074   //******************************************************************************
00075 
00076   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00077     base::LocCtrl = Ctrl.getSubDevice(prefix+"InternalBoundary");
00078     RegisterAt(base::LocCtrl,"BoundaryWidth",_BoundaryWidth);
00079     RegisterAt(base::LocCtrl,"Cutoff",_Cutoff);
00080     RegisterAt(base::LocCtrl,"ErrorValue",_ErrorValue);
00081     RegisterAt(base::LocCtrl,"BndType",_bndtype);
00082     RegisterAt(base::LocCtrl,"NAux",_naux);
00083     RegisterAt(base::LocCtrl,"Verbose",_Verbose);    
00084   } 
00085   virtual void register_at(ControlDevice& Ctrl) {
00086     register_at(Ctrl, "");
00087   }
00088   virtual void init() {}
00089   virtual void update() {}
00090   virtual void finish() {} 
00091   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00092     base::SetupData(gh, ghosts);
00093     if (_bndtype == EXTRAPOLATION) _normal_factor = 1.0;
00094     else if (_bndtype == REFLECTION) _normal_factor = 2.0;
00095     else _normal_factor = 0.0;
00096     _Cutoff -= 1.e-12;
00097   }
00098 
00099   void SetPhysbd(vec_grid_fct_type& u, grid_fct_type& phi, 
00100                  bool_grid_fct_type& bf, const int Time, 
00101                  const int Level, double t, bool extrapolate=false) { 
00102     START_WATCH
00103     forall (u,Time,Level,c)
00104       int nc, ni;
00105       int* idx = (int *) 0;
00106       int* idx_wrg = (int *) 0;
00107       BBox bb = u.interiorbbox(Time,Level,c);
00108       bb.grow(std::min(NBoundaryWidth(),base::NGhosts()-1));
00109 
00110       START_WATCH
00111         GetBndryCells(u,phi(Time,Level,c),bf(Time,Level,c),bb,Time,Level,c,
00112                       nc,idx,ni,idx_wrg);
00113       END_WATCH(GFM_FINDING_CELLS)
00114 
00115       if (_Verbose && (idx_wrg!=(int *) 0)) {
00116         for (register int cnt=0; cnt<ni; cnt++) {
00117           std::cout << "Unrecoverable ambigious ghost cell " 
00118                     << Coords(base::Dim(),&(idx_wrg[base::Dim()*cnt])) 
00119                     << " on Level " << Level << std::endl;
00120 #ifdef DEBUG_PRINT      
00121           ( comm_service::log() << "GFMBoundaryBase: " << bb 
00122             << " Unrecoverable ambigious ghost cell " 
00123             << Coords(base::Dim(),&(idx_wrg[base::Dim()*cnt])) 
00124             << " on Level " << Level << "\n" ).flush();
00125 #endif
00126         }
00127       }
00128 
00129       point_type* xc = new point_type[nc];
00130       point_type* normal = new point_type[nc];
00131       DataType* distance = new DataType[nc];
00132 
00133       START_WATCH
00134         geom_base::CellCoords(base::GH(),bb.stepsize(),nc,idx,xc);
00135         gfm_geom_base::CellNormals(base::GH(),phi(Time,Level,c),nc,idx,normal);
00136       END_WATCH(GFM_GEOMETRY)
00137 
00138       phi(Time,Level,c).plus(Cutoff());
00139 
00140       START_WATCH
00141         gfm_geom_base::CellDistances(base::GH(),phi(Time,Level,c),nc,idx,distance);
00142       END_WATCH(GFM_GEOMETRY)
00143 
00144       if (extrapolate) {
00145         DataType factor = 1.0;
00146         VectorType* uv = (VectorType*) 0;
00147         Extrapolation(u(Time,Level,c),phi(Time,Level,c),uv,Level,nc,idx,xc,
00148                       distance,normal,factor);
00149         for (register int cnt=0; cnt<nc; cnt++) 
00150           u(Time,Level,c)(Coords(dim,idx+dim*cnt)) = uv[cnt];
00151         if (uv) delete [] uv;
00152       }
00153       else
00154         SetGrid(u(Time,Level,c),phi(Time,Level,c),bb,Level,t,
00155                 nc,idx,xc,distance,normal);  
00156 
00157       phi(Time,Level,c).minus(Cutoff());
00158 
00159       if (idx) delete [] idx;
00160       if (idx_wrg) delete [] idx_wrg;
00161       if (xc) delete [] xc;
00162       if (normal) delete [] normal;
00163       if (distance) delete [] distance;
00164 
00165     end_forall
00166     END_WATCH(GFM_SETBNDRY_WHOLE)
00167 
00168     START_WATCH
00169       Sync(u,Time,Level);
00170     END_WATCH(BOUNDARIES_SYNC)
00171     START_WATCH
00172       ExternalBoundaryUpdate(u,Time,Level);
00173     END_WATCH(BOUNDARIES_EXTERNAL)
00174   }
00175 
00176   virtual void Extrapolation(vec_grid_data_type& gdu, grid_data_type& gdphi, 
00177                              VectorType*& u, const int& Level, 
00178                              const int& nc, const int* idx, 
00179                              const point_type* xc, DataType* distance, 
00180                              point_type* normal, const DataType& normal_factor) {
00181     
00182     START_WATCH
00183     if (! ValidNormalFactor(normal_factor)) {
00184       std::cout << base::NGhosts() << " ghost cells too few for interpolation. "
00185                 << "At least " << int(NBoundaryWidth()*(normal_factor+1)) 
00186                 << " cells recommended.\n";
00187 #ifdef DEBUG_PRINT
00188       ( comm_service::log() << "GFMBoundaryBase: " << base::NGhosts() 
00189         << " ghost cells too few for interpolation. " << "At least " 
00190         << int(NBoundaryWidth()*(normal_factor+1)) << " cells recommended.\n" ).flush();
00191 #endif
00192     }
00193 
00194     if (u) delete [] u;
00195     u = new VectorType[nc];
00196     point_type* pos = new point_type[nc];
00197 
00198     register int n;
00199     for (n=0; n<nc; n++) {
00200       pos[n] = xc[n]+normal_factor*distance[n]*normal[n];
00201       u[n] = _ErrorValue;
00202     }
00203 
00204     if (_interpolation)
00205       Interpolation_().Interpolate(base::GH(),gdu,gdphi,nc,pos,u,_ErrorValue);
00206 
00207     for (n = 0; n<nc; n++) 
00208       if (u[n](0) == _ErrorValue) {
00209         Coords co(dim,idx+dim*n);
00210         u[n] = gdu(co);
00211 #ifdef DEBUG_PRINT
00212         if (_Verbose) {
00213           ( comm_service::log() << "GFMBoundaryBase: " << gdu.bbox() 
00214             << " Interpolation error in " << co 
00215             << "=" << xc[n] << " from " << pos[n] << " on Level " 
00216             << Level << "\n" ).flush();
00217           ( comm_service::log() << "                 " << distance[n] 
00218             << " " << normal[n] << "\n" ).flush();
00219         }
00220 #endif
00221         normal[n] = static_cast<DataType>(0.0);
00222         distance[n] = static_cast<DataType>(0.0);
00223       }
00224     delete [] pos;
00225     END_WATCH(GFM_EXTRAPOLATION)
00226   }
00227 
00228   inline void SetInterpolation(interpolation_type* inter) { 
00229     if (_interpolation) delete _interpolation;
00230     _interpolation = inter; 
00231   }
00232   inline interpolation_type& Interpolation_() { return *_interpolation; }
00233   inline const interpolation_type& Interpolation_() const { return *_interpolation; }
00234 
00235   inline const int& NAux() const { return _naux; }
00236   inline void SetNAux(const int naux) { _naux = naux; }
00237 
00238   inline bool ValidNormalFactor(const DataType& nf) const 
00239   { return (NBoundaryWidth()*(nf+1) <= base::NGhosts()); }
00240   inline const DataType& NormalFactor() const { return _normal_factor; }
00241   inline const int& BndType() const { return _bndtype; }
00242   inline int isgn(const DataType val)
00243   { return (val>0 ? 1 : (val<0 ? -1 : 0)); }
00244   inline const int& NBoundaryWidth() const { return _BoundaryWidth; }
00245   inline const DataType& Cutoff() const { return _Cutoff; }
00246   inline const DataType& FarAway() const { return _far_away; }
00247   inline void SetVerbose(const int v) { _Verbose = v; }  
00248   inline const int& GetVerbose() const { return _Verbose; }
00249 
00250 protected:
00251   int _BoundaryWidth, _bndtype, _naux, _Verbose;
00252   DataType _Cutoff, _ErrorValue, _normal_factor, _far_away;
00253   interpolation_type* _interpolation;
00254 };
00255 
00256 template <class VectorType, int dim> class GFMBoundary;
00257 
00258 
00264 template <class VectorType>
00265 class GFMBoundary<VectorType,2> : public GFMBoundaryBase<VectorType,2> {
00266   typedef typename VectorType::InternalDataType DataType;
00267   typedef GFMBoundaryBase<VectorType,2> base; 
00268 
00269 public:
00270   typedef typename base::vec_grid_data_type vec_grid_data_type;
00271   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00272   typedef typename base::grid_data_type grid_data_type;
00273   typedef typename base::bool_grid_data_type bool_grid_data_type;
00274   typedef typename base::point_type point_type;
00275 
00276   GFMBoundary() : base() {}
00277 
00278   virtual ~GFMBoundary() {}
00279 
00280   virtual void GetBndryCells(vec_grid_fct_type& u, grid_data_type& gdphi, bool_grid_data_type& gdbf, 
00281                              const BBox& bb, const int& Time, const int& Level, const int& c,
00282                              int& nc, int*& idx, int& Nillegal, int*& idx_wrg) {
00283 
00284     GridData<char,2> flag(gdphi.bbox());
00285     flag.equals(0);
00286     bool_grid_data_type pf(gdphi.bbox());
00287     BeginFastIndex2(phi, gdphi.bbox(), gdphi.data(), DataType);
00288     BeginFastIndex2(flag, flag.bbox(), flag.data(), char);
00289     BeginFastIndex2(bflag, gdbf.bbox(), gdbf.data(), bool);
00290     BeginFastIndex2(pf, pf.bbox(), pf.data(), bool);
00291     
00292     Sgn(gdphi,pf,base::Cutoff());
00293 
00294     Coords gdss(gdphi.bbox().stepsize());
00295     DCoords dx = base::GH().worldStep(gdss);    
00296     double distance = 0.;
00297     for (register int n=0; n<base::Dim(); n++)
00298       distance += dx[n]*dx[n];
00299     distance = std::sqrt(distance)-base::Cutoff();
00300     for_2 (n, m, gdphi.bbox(), gdss)
00301       if (FastIndex2(phi,n,m)<-base::Cutoff() || 
00302           FastIndex2(phi,n,m)>distance) continue;
00303       BBox mark(2,n,m,n,m,gdss(0),gdss(1));
00304       mark.grow(base::NBoundaryWidth());
00305       flag.equals(1,mark);
00306     end_for
00307 
00308     Coords ss(bb.stepsize());
00309     nc = 0;   
00310     Nillegal = 0;   
00311     for_2 (n, m, bb, ss)
00312       if (FastIndex2(pf,n,m))
00313         FastIndex2(flag,n,m) = 0;
00314       else {
00315         if (FastIndex2(flag,n,m)==1) {
00316           if (FastIndex2(pf,n+ss(0),m) && FastIndex2(pf,n-ss(0),m) &&
00317               FastIndex2(pf,n,m+ss(1)) && FastIndex2(pf,n,m-ss(1))) {
00318             if (u.interiorbbox(Time,Level,c).inside(n,m)) {
00319               bool illegal_point = true;
00320               for (register int ci=0; ci<u.children(Time,Level,c); ci++) 
00321                 if (u.interiorbbox(CurrentTime(base::GH(),Level+1),Level+1,
00322                                    u.childidx(Time,Level,c,ci)).inside(n,m)) {
00323                   illegal_point = false;
00324                   break;
00325                 }      
00326               if (illegal_point) FastIndex2(flag,n,m) = -1;
00327             }
00328             else 
00329               FastIndex2(flag,n,m) = -2;
00330           }
00331           if (FastIndex2(flag,n,m)==-1)
00332             Nillegal++;
00333           else if (FastIndex2(flag,n,m)==1)
00334             nc++;
00335         }
00336         else 
00337           if (FastIndex2(phi,n,m)>=-base::FarAway() && 
00338               FastIndex2(phi,n,m)<-base::Cutoff())
00339           FastIndex2(flag,n,m) = 2;
00340       }
00341     end_for
00342       
00343     if (idx) delete [] idx;
00344     if (idx_wrg) delete [] idx_wrg;
00345     
00346     idx = new int[2*nc];
00347     if (base::_Verbose) 
00348       idx_wrg = new int[2*Nillegal];
00349     int cnti = 0, cntill = 0;
00350    
00351     for_2 (n, m, gdbf.bbox(), ss)
00352       if (bb.inside(n,m)) {
00353         if (FastIndex2(flag,n,m)>0)
00354           FastIndex2(bflag,n,m) = true; 
00355         if (FastIndex2(flag,n,m)==1) {
00356           idx[2*cnti] = n; idx[2*cnti+1] = m; cnti++;
00357         }
00358         else if (FastIndex2(flag,n,m)==-1 && base::_Verbose) { 
00359           idx_wrg[2*cntill] = n; idx_wrg[2*cntill+1] = m; cntill++;
00360         }
00361       }
00362       else
00363         FastIndex2(bflag,n,m) = true; 
00364     end_for
00365  
00366     EndFastIndex2(pf);
00367     EndFastIndex2(bflag);
00368     EndFastIndex2(flag);
00369     EndFastIndex2(phi);
00370   }
00371 
00372   virtual void Sgn(grid_data_type& gdphi, bool_grid_data_type& gdbf, const DataType& c) {
00373     BBox OpBox = gdphi.bbox();
00374     Coords& OpBox_stepsize = OpBox.stepsize();
00375     BeginFastIndex2(phi, OpBox, gdphi.data(), DataType);
00376     BeginFastIndex2(bf, OpBox, gdbf.data(), bool);
00377     for_2 (n, m, OpBox, OpBox_stepsize)
00378       FastIndex2(bf,n,m) = (FastIndex2(phi,n,m)>=-c); 
00379     end_for
00380     EndFastIndex2(bf);
00381     EndFastIndex2(phi);
00382   }
00383 };
00384 
00391 template <class VectorType>
00392 class GFMBoundary<VectorType,3> : public GFMBoundaryBase<VectorType,3> {
00393   typedef typename VectorType::InternalDataType DataType;
00394   typedef GFMBoundaryBase<VectorType,3> base; 
00395 
00396 public:
00397   typedef typename base::vec_grid_data_type vec_grid_data_type;
00398   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00399   typedef typename base::grid_data_type grid_data_type;
00400   typedef typename base::bool_grid_data_type bool_grid_data_type;
00401   typedef typename base::point_type point_type;
00402 
00403   GFMBoundary() : base() {}
00404 
00405   virtual ~GFMBoundary() {}
00406 
00407   virtual void GetBndryCells(vec_grid_fct_type& u, grid_data_type& gdphi, bool_grid_data_type& gdbf, 
00408                              const BBox& bb, const int& Time, const int& Level, const int& c,
00409                              int& nc, int*& idx, int& Nillegal, int*& idx_wrg) {
00410 
00411     GridData<char,3> flag(gdphi.bbox());
00412     flag.equals(0);
00413     bool_grid_data_type pf(gdphi.bbox());
00414     BeginFastIndex3(phi, gdphi.bbox(), gdphi.data(), DataType);
00415     BeginFastIndex3(flag, flag.bbox(), flag.data(), char);
00416     BeginFastIndex3(bflag, gdbf.bbox(), gdbf.data(), bool);
00417     BeginFastIndex3(pf, pf.bbox(), pf.data(), bool);
00418 
00419     Sgn(gdphi,pf,base::Cutoff());
00420 
00421     Coords gdss(gdphi.bbox().stepsize());
00422     DCoords dx = base::GH().worldStep(gdss);    
00423     double distance = 0.;
00424     for (register int n=0; n<base::Dim(); n++)
00425       distance += dx[n]*dx[n];
00426     distance = std::sqrt(distance)-base::Cutoff();
00427     for_3 (n, m, l, gdphi.bbox(), gdss)
00428       if (FastIndex3(phi,n,m,l)<-base::Cutoff() ||
00429           FastIndex3(phi,n,m,l)>distance) continue;
00430       BBox mark(3,n,m,l,n,m,l,gdss(0),gdss(1),gdss(2));
00431       mark.grow(base::NBoundaryWidth());
00432       flag.equals(1,mark);
00433     end_for
00434 
00435     Coords ss(bb.stepsize());
00436     nc = 0;   
00437     Nillegal = 0;   
00438     for_3 (n, m, l, bb, ss)
00439       if (FastIndex3(pf,n,m,l))
00440         FastIndex3(flag,n,m,l) = 0;
00441       else {
00442         if (FastIndex3(flag,n,m,l)==1) {
00443           if (FastIndex3(pf,n+ss(0),m,l) && FastIndex3(pf,n-ss(0),m,l) &&
00444               FastIndex3(pf,n,m+ss(1),l) && FastIndex3(pf,n,m-ss(1),l) &&
00445               FastIndex3(pf,n,m,l+ss(2)) && FastIndex3(pf,n,m,l-ss(2))) {
00446             if (u.interiorbbox(Time,Level,c).inside(n,m,l)) {
00447               bool illegal_point = true;
00448               for (register int ci=0; ci<u.children(Time,Level,c); ci++) 
00449                 if (u.interiorbbox(CurrentTime(base::GH(),Level+1),Level+1,
00450                                    u.childidx(Time,Level,c,ci)).inside(n,m,l)) {
00451                   illegal_point = false;
00452                   break;
00453                 }      
00454               if (illegal_point) FastIndex3(flag,n,m,l) = -1;
00455             }
00456             else 
00457               FastIndex3(flag,n,m,l) = -2;
00458           }
00459           if (FastIndex3(flag,n,m,l)==-1)
00460             Nillegal++;
00461           else if (FastIndex3(flag,n,m,l)==1)
00462             nc++;
00463         }
00464         else if (FastIndex3(phi,n,m,l)>=-base::FarAway() && 
00465                  FastIndex3(phi,n,m,l)<-base::Cutoff())
00466           FastIndex3(flag,n,m,l) = 2;
00467       }
00468     end_for
00469       
00470     if (idx) delete [] idx;
00471     if (idx_wrg) delete [] idx_wrg;
00472     
00473     idx = new int[3*nc];
00474     if (base::_Verbose)
00475       idx_wrg = new int[3*Nillegal];
00476     int cnti = 0, cntill = 0;
00477    
00478     for_3 (n, m, l, gdbf.bbox(), ss)
00479       if (bb.inside(n,m,l)) {
00480         if (FastIndex3(flag,n,m,l)>0)
00481           FastIndex3(bflag,n,m,l) = true; 
00482         if (FastIndex3(flag,n,m,l)==1) {
00483           idx[3*cnti] = n; idx[3*cnti+1] = m; idx[3*cnti+2] = l; cnti++;
00484         }
00485         else if (FastIndex3(flag,n,m,l)==-1 && base::_Verbose) { 
00486           idx_wrg[3*cntill] = n; idx_wrg[3*cntill+1] = m; 
00487           idx_wrg[3*cntill+2] = l; cntill++;
00488         }
00489       }
00490       else
00491         FastIndex3(bflag,n,m,l) = true; 
00492     end_for
00493 
00494     EndFastIndex3(pf);
00495     EndFastIndex3(bflag);
00496     EndFastIndex3(flag);
00497     EndFastIndex3(phi);
00498   }
00499 
00500   virtual void Sgn(grid_data_type& gdphi, bool_grid_data_type& gdbf, const DataType& c) {
00501     BBox OpBox = gdphi.bbox();
00502     Coords& OpBox_stepsize = OpBox.stepsize();
00503     BeginFastIndex3(phi, OpBox, gdphi.data(), DataType);
00504     BeginFastIndex3(bf, OpBox, gdbf.data(), bool);
00505     for_3 (n, m, l, OpBox, OpBox_stepsize)
00506       FastIndex3(bf,n,m,l) = (FastIndex3(phi,n,m,l)>=-c); 
00507     end_for
00508     EndFastIndex3(bf);
00509     EndFastIndex3(phi);
00510   }
00511 };
00512 
00513 #endif

Generated on Fri Aug 24 13:00:51 2007 for AMROC Fluid-solver Framework - by  doxygen 1.4.7