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