00001 
00002 
00003 
00004 
00005 
00006 #ifndef AMROC_FILE_INPUT_H
00007 #define AMROC_FILE_INPUT_H
00008 
00016 #include "DAGHIO.h"
00017 
00018 #define MAXLEV 11
00019 
00026 template <class VectorType, int dim>
00027 class FileInput : public AMRBase<VectorType,dim> {
00028   typedef AMRBase<VectorType,dim> base;
00029   typedef typename VectorType::InternalDataType DataType;
00030 
00031 public:
00032   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00033   typedef typename base::vec_grid_data_type vec_grid_data_type;
00034   typedef GridFunction<DataType,dim> grid_fct_type;
00035   typedef GridData<DataType,dim> grid_data_type;
00036 
00037   typedef LevelTransfer<DataType,dim> leveltransfer_type;
00038   typedef GFLevelTransferSpecificFunc<leveltransfer_type,DataType,dim> leveltransfer_functor_type;
00039 
00040   FileInput() : base(), _Ncnt(0), _InputType(DAGHIO_HDF_NCSA) {
00041     CompName = (std::string*) 0;
00042     _LevelTransfer = (leveltransfer_type *) 0;
00043     _ProlongFunc = (leveltransfer_functor_type *) 0;
00044     _RestrictFunc = (leveltransfer_functor_type *) 0;
00045     _where = new BBox[MAXLEV];
00046     _to = new BBox[MAXLEV];
00047     _from = new BBox[MAXLEV];
00048     for (register int l=0; l<MAXLEV; l++) {
00049       _where[l] = BBox(dim,1);
00050       _to[l] = BBox(dim,1);
00051       _from[l] = BBox(dim,1);
00052     }
00053   } 
00054 
00055   virtual ~FileInput() {
00056     if (CompName) delete [] CompName;
00057     if (_ProlongFunc) delete _ProlongFunc;
00058     if (_RestrictFunc) delete _RestrictFunc;
00059     delete [] _where;
00060     delete [] _to;
00061     delete [] _from;
00062   }
00063 
00064   virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
00065   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00066     LocCtrl = Ctrl.getSubDevice(prefix+"InputFiles");    
00067     RegisterAt(LocCtrl,"Type",_InputType);
00068     if (CompName) delete [] CompName;
00069     CompName = new std::string[MAXCOMPONENTS];
00070     for (int d=0; d<MAXCOMPONENTS; d++) {
00071       char Name[16];
00072       std::sprintf(Name,"-");
00073       CompName[d] = Name;
00074       std::sprintf(Name,"InputName(%d)",d+1);
00075       RegisterAt(LocCtrl,Name,CompName[d]);
00076     }    
00077     for (register int l=0; l<MAXLEV; l++) {
00078       char devno[8];
00079       std::sprintf(devno,"(%d)",l);
00080       WhereCtrl = LocCtrl.getSubDevice(prefix+"Where"+std::string(devno));
00081       for (int d=0; d<dim; d++) {
00082         char text[8];
00083         std::sprintf(text,"lb(%d)",d+1);
00084         RegisterAt(WhereCtrl,text,_where[l].lower(d));
00085         std::sprintf(text,"ub(%d)",d+1);
00086         RegisterAt(WhereCtrl,text,_where[l].upper(d));
00087         std::sprintf(text,"s(%d)",d+1);
00088         RegisterAt(WhereCtrl,text,_where[l].stepsize(d));
00089       }
00090       ToCtrl = LocCtrl.getSubDevice(prefix+"To"+std::string(devno));
00091       for (int d=0; d<dim; d++) {
00092         char text[8];
00093         std::sprintf(text,"lb(%d)",d+1);
00094         RegisterAt(ToCtrl,text,_to[l].lower(d));
00095         std::sprintf(text,"ub(%d)",d+1);
00096         RegisterAt(ToCtrl,text,_to[l].upper(d));
00097         std::sprintf(text,"s(%d)",d+1);
00098         RegisterAt(ToCtrl,text,_to[l].stepsize(d));
00099       }
00100       FromCtrl = LocCtrl.getSubDevice(prefix+"From"+std::string(devno));
00101       for (int d=0; d<dim; d++) {
00102         char text[8];
00103         std::sprintf(text,"lb(%d)",d+1);
00104         RegisterAt(FromCtrl,text,_from[l].lower(d));
00105         std::sprintf(text,"ub(%d)",d+1);
00106         RegisterAt(FromCtrl,text,_from[l].upper(d));
00107         std::sprintf(text,"s(%d)",d+1);
00108         RegisterAt(FromCtrl,text,_from[l].stepsize(d));
00109       }
00110     }
00111     if (_LevelTransfer) _LevelTransfer->register_at(Ctrl, ""); 
00112   }
00113   virtual void init() { 
00114     DAGHIOInit(); 
00115     if (_LevelTransfer) _LevelTransfer->init();
00116   }
00117   virtual void update() {
00118     int d=MAXCOMPONENTS;
00119     for (; d>0; d--) 
00120       if (CompName[d-1].c_str()[0] != '-') 
00121         break;
00122     if (d>0) SetNcnt(d);
00123     if (_LevelTransfer) _LevelTransfer->update();
00124   }
00125   virtual void finish() {
00126     if (base::_Hierarchy) CloseIO();
00127     if (base::_Hierarchy) DAGHIOEnd(base::GH());  
00128     if (CompName) {
00129       delete [] CompName;
00130       CompName = (std::string*) NULL;
00131     }
00132     if (_LevelTransfer) _LevelTransfer->finish();
00133   }
00134 
00135   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00136     base::SetupData(gh, ghosts);
00137     DAGHSetIOType(base::GH(), _InputType);   
00138     if (_LevelTransfer) _LevelTransfer->SetupData(gh, ghosts);
00139     if (_LevelTransfer) {
00140       _ProlongFunc = new leveltransfer_functor_type(_LevelTransfer, &leveltransfer_type::Prolong);
00141       _RestrictFunc = new leveltransfer_functor_type(_LevelTransfer, &leveltransfer_type::Restrict);
00142     }
00143     for (register int l=0; l<=MaxLevel(base::GH()); l++) {
00144       if (_to[l].empty()) {
00145         _to[l] = base::GH().wholebbox();
00146         _to[l].refine(base::GH().refinedby(l));
00147       }
00148     }
00149   }
00150 
00151   
00152   
00153   
00154   virtual void ReadIn(vec_grid_fct_type& u, grid_fct_type& IOfunc) {
00155     for (int cnt=0; cnt<Ncnt(); cnt++) { 
00156       if (CompName[cnt].c_str()[0] == '-') 
00157         continue;
00158       int lev;
00159       for (lev=0; lev<=FineLevel(base::GH()); lev++) {
00160         int Time = CurrentTime(base::GH(),lev); 
00161         forall (u,Time,lev,c)           
00162           equals_from(IOfunc(Time,lev,c), u(Time,lev,c), cnt);
00163         end_forall      
00164       }
00165       ReadIn(IOfunc,CompName[cnt].c_str());
00166       for (lev=0; lev<=FineLevel(base::GH()); lev++) {
00167         int Time = CurrentTime(base::GH(),lev); 
00168         forall (u,Time,lev,c)           
00169           equals_to(u(Time,lev,c), IOfunc(Time,lev,c), cnt);
00170         end_forall      
00171       }
00172     } 
00173   }   
00174 
00175   
00176   
00177   
00178   virtual void ReadIn(grid_fct_type& IOfunc, const char* name) {
00179     for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00180       int Time = CurrentTime(base::GH(),lev); 
00181       ReadIn(IOfunc,name,Time,lev);
00182     }
00183   }   
00184 
00185   
00186   
00187   
00188   virtual void ReadIn(grid_fct_type& IOfunc, const char* name, 
00189                       const int& Time, const int& Level) {
00190 
00191     int old_prolong_flag = DAGHTrue;
00192     if (!IOfunc.prolong()) old_prolong_flag = DAGHFalse;
00193     GFLevelTransferFunc<DataType,dim>* old_ProlongFunc = GetProlongFunction(IOfunc);
00194     int old_restrict_flag = DAGHTrue;
00195     if (!IOfunc.Restrict()) old_restrict_flag = DAGHFalse;
00196     GFLevelTransferFunc<DataType,dim>* old_RestrictFunc = GetRestrictFunction(IOfunc);
00197 
00198     if (_LevelTransfer && _ProlongFunc) {
00199       SetProlongFunction(IOfunc, _ProlongFunc);
00200       IOfunc.GF_SetProlongFlag(DAGHTrue);
00201       if (Level>0) {
00202         int TimeCoarse = CurrentTime(base::GH(),Level-1); 
00203         Prolong(IOfunc, TimeCoarse, Level-1, Time, Level);
00204       }
00205     }
00206     if (_LevelTransfer && _RestrictFunc) {
00207       SetRestrictFunction(IOfunc, _RestrictFunc);
00208       IOfunc.GF_SetRestrictFlag(DAGHTrue);
00209     }
00210 
00211     if (_where[Level].empty()) {
00212       if (_to[Level].empty() || _from[Level].empty())
00213         ReadPlain(IOfunc,name,Time,Level);
00214       else {
00215         BBox to(_to[Level]), from(_from[Level]);
00216         Coords ss(StepSize(IOfunc,Level));
00217         for (int i=0;i<dim;i++) {
00218           if (to.stepsize(i)>ss(i))
00219             to.refine(i,to.stepsize(i)/ss(i),0);
00220           else to.coarsen(i,ss(i)/to.stepsize(i));
00221         }
00222         ReadPlain(IOfunc,name,Time,Level,to,from);
00223       }
00224     }
00225     else {
00226       BBox where(_where[Level]);
00227       Coords ss(StepSize(IOfunc,Level));
00228       for (int i=0;i<dim;i++) {
00229         if (where.stepsize(i)>ss(i))
00230           where.refine(i,where.stepsize(i)/ss(i),0);
00231         else where.coarsen(i,ss(i)/where.stepsize(i));
00232       }
00233       ReadPlain(IOfunc,name,Time,Level,where);
00234     }
00235 
00236     SetProlongFunction(IOfunc, old_ProlongFunc);
00237     IOfunc.GF_SetProlongFlag(old_prolong_flag);
00238     SetRestrictFunction(IOfunc, old_RestrictFunc);
00239     IOfunc.GF_SetRestrictFlag(old_restrict_flag);
00240   }   
00241 
00242   void ReadPlain(grid_fct_type& IOfunc, const char* name, 
00243                  const int& Time, const int& Level,
00244                  const BBox& where=BBox::_empty_bbox) {
00245     int me = MY_PROC;     
00246     char ioname[DAGHBktGFNameWidth+16];
00247     std::sprintf(ioname,"%s",name);
00248     double t_old = GetPhysicalTime(IOfunc,Time,Level);
00249     if (me == VizServer) 
00250       ( std::cout << " *** Reading " << ioname << " " ).flush();
00251     if (where.empty()) 
00252       Read(IOfunc,Time,Level,DAGH_Double,ioname);      
00253     else
00254       Read(IOfunc,Time,Level,where,DAGH_Double,ioname);      
00255     double t = GetPhysicalTime(IOfunc,Time,Level);
00256     if (me == VizServer) 
00257       std::cout << "at t = " << t << " to level " << Level 
00258                 << " time index " << Time << std::endl;
00259     SetPhysicalTime(IOfunc,Time,Level,t_old);
00260 #ifndef DAGH_NO_MPI
00261     MPI_Barrier(comm_service::comm()); 
00262 #endif
00263   }   
00264 
00265   void ReadPlain(grid_fct_type& IOfunc, const char* name, 
00266                  const int& Time, const int& Level,
00267                  const BBox& to, const BBox& from) {
00268     int me = MY_PROC;     
00269     char ioname[DAGHBktGFNameWidth+16];
00270     std::sprintf(ioname,"%s",name);
00271     double t_old = GetPhysicalTime(IOfunc,Time,Level);
00272     if (me == VizServer) 
00273       ( std::cout << " *** Reading " << ioname << " " ).flush();
00274     if (to.empty() || from.empty()) 
00275       Read(IOfunc,Time,Level,DAGH_Double,ioname);      
00276     else
00277       Read(IOfunc,Time,Level,to,from,DAGH_Double,ioname);      
00278     double t = GetPhysicalTime(IOfunc,Time,Level);
00279     if (me == VizServer) 
00280       std::cout << "at t = " << t << " to level " << Level 
00281                 << " time index " << Time << std::endl;
00282     SetPhysicalTime(IOfunc,Time,Level,t_old);
00283 #ifndef DAGH_NO_MPI
00284     MPI_Barrier(comm_service::comm()); 
00285 #endif
00286   }   
00287 
00288   
00289   
00290   
00291   virtual void ReadIn(grid_data_type& IOdata, const char* name, 
00292                       const int& Time, const int& Level) {
00293     ReadPlain(IOdata,name,Time,Level);
00294   }   
00295 
00296   virtual void ReadPlain(grid_data_type& IOdata, const char* name, 
00297                          const int& Time, const int& Level) {
00298     int me = MY_PROC;     
00299     char ioname[DAGHBktGFNameWidth+16];
00300     std::sprintf(ioname,"%s_%d",name,Time); 
00301 
00302     gdhdr head;
00303     head.type = DAGHSingle;
00304     head.owner = me;
00305     head.gfid = 10000;
00306     head.gfdatatype = DAGH_Double;
00307     head.gfstaggertype = DAGHCellCentered;
00308     std::strncpy(head.gfname,ioname,DAGHBktGFNameWidth-1);
00309     head.bbox = IOdata.bbox();
00310     head.time = 0;
00311     head.time_value = Time;
00312     double Current_Time = 0.;
00313     std::memcpy(head.physical_time,(char *)&Current_Time,sizeof(double));
00314     head.level = Level;
00315     head.index = 0;
00316     if (me == VizServer) 
00317       ( std::cout << " *** Reading " << ioname << " " ).flush();
00318 
00319     (base::GH().DAGH_IORead())(base::GH(), &head, IOdata.data());
00320 
00321     std::memcpy((char *)&Current_Time,head.physical_time,sizeof(double));
00322     if (me == VizServer) 
00323       std::cout << "at t = " << Current_Time << std::endl;
00324   }  
00325 
00326   virtual void CloseIO() { DAGHIOClose(base::GH()); }
00327   inline void SetNcnt(const int& cnt) { _Ncnt = cnt; }
00328   inline const int& Ncnt() const { return _Ncnt; }
00329   inline const int& InputType() const { return _InputType; }
00330   inline const char* FileName(int d) const { return CompName[d].c_str(); }
00331 
00332   inline void SetLevelTransfer(leveltransfer_type* _leveltransfer) 
00333   { _LevelTransfer = _leveltransfer; }
00334   inline leveltransfer_type& LevelTransfer_() { return *_LevelTransfer; }
00335   inline const leveltransfer_type& LevelTransfer_() const { return *_LevelTransfer; }
00336 
00337 protected:
00338   int _Ncnt, _InputType;
00339   std::string* CompName;
00340   ControlDevice LocCtrl, WhereCtrl, ToCtrl, FromCtrl;
00341   leveltransfer_type* _LevelTransfer;
00342   leveltransfer_functor_type *_ProlongFunc, *_RestrictFunc;  
00343   BBox *_where, *_to, *_from;
00344 };
00345 
00346 
00347 #endif