File List
amroc/converters/HDFToFile/FileVisualGrid.h

00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 //
00006 // Copyright (C) 2003-2007 California Institute of Technology
00007 // Ralf Deiterding, ralf@cacr.caltech.edu
00008 // Santiago Lombeyda (VTU, SILO)
00009 
00010 #ifndef AMROC_FILEVISUALGRID_H
00011 #define AMROC_FILEVISUALGRID_H
00012 
00020 #include "VisualGrid.h" 
00021 #include "HDFToFile/FileVisualGridBase.h"
00022 #include "EncodeToBase64.h"
00023 #include <strings.h>
00024 #include <cctype>
00025 #include <hlimits.h>
00026 
00027 #if defined(HAVE_LIBSILO) || defined(HAVE_LIBSILOH5)
00028 #include <silo.h>
00029 #endif
00030 
00031 #define PROBE_TOLERANCE 1.e-6
00032 
00033 typedef float points_list_type[][3];
00034 typedef int   cons_list_type[][8];
00035 
00042 template <class VectorType>
00043 class FileVisualGrid : public FileVisualGridBase, 
00044                        public VisualGrid<VectorType> {
00045   typedef typename VectorType::InternalDataType DataType;
00046   typedef VisualGrid<VectorType> base;
00047 public:
00048   typedef typename base::data_block_type data_block_type;
00049   typedef typename base::block_list_type block_list_type;
00050   typedef typename base::hierarchy_storage_type hierarchy_storage_type;
00051   typedef typename base::hdata_block_type hdata_block_type;
00052   typedef typename base::evaluator_type evaluator_type;
00053   typedef typename base::point_list_type point_list_type;
00054 
00055   typedef std::multimap<float, int> index_ordered_map_type;
00056   typedef std::pair<float, int> ordered_index_type;
00057 
00058 public:
00059   FileVisualGrid(evaluator_type* ev) : FileVisualGridBase(ev), 
00060                                        base(ev), OutputType(0), OutputStep(-1) {
00061     base::_BlockList = (block_list_type *) new block_list_type;
00062     nkeys = 0;
00063     keys = (int *) 0;
00064     OutputName = "-";
00065     for (int d=0; d<3; d++) {
00066       LinePoint[d]  = 0.0;   
00067       LineVector[d] = 0.0;   
00068       RegData[d] = 10;
00069     } 
00070     RegLim[0] = FLT_MIN;
00071     RegLim[1] = FLT_MAX;
00072     RegMax = 1;
00073     LevelCutOut = 0;
00074     PointOutput = 0;
00075     _Step = -1;
00076     FileMerge = 1;
00077     InputFile = "";
00078   }
00079 
00080   ~FileVisualGrid() {
00081     if (keys) delete [] keys;
00082     for (typename block_list_type::iterator bit=BlockList().begin();
00083          bit!=BlockList().end(); bit++)
00084       delete (*bit);
00085     BlockList().clear();      
00086     delete base::_BlockList;
00087   }
00088 
00089   virtual void update() {    
00090     std::string iname;
00091     if (UpdateName()) iname = *UpdateName();
00092     else iname = "display_file.in";
00093 
00094     // Generating display_file.in if file is missing; 
00095     std::ifstream infile;
00096     std::ofstream outfile;
00097     infile.open(iname.c_str(), std::ios::in);
00098     if (!infile) {
00099       if (!UpdateName()) {
00100         outfile.open(iname.c_str(), std::ios::out);
00101         outfile << "Type  6\n\nDisplayMinLevel    0\nDisplayMaxLevel    10\n\n";
00102         outfile << "FileName convert\nFileType 6\nKeys p,t,s,i\n"; 
00103         outfile.close();
00104       }
00105     }
00106     else{
00107       infile.close();
00108     }
00109 
00110     ControlDevice DisplayCtrl(GetFileControlDevice(iname.c_str(),""));
00111     RegisterAt(DisplayCtrl,"Type",base::_FKeys); 
00112     RegisterAt(DisplayCtrl,"FileName",OutputName);     
00113     RegisterAt(DisplayCtrl,"FileType",OutputType);   
00114     RegisterAt(DisplayCtrl,"FileStep",OutputStep);   
00115     RegisterAt(DisplayCtrl,"Keys",StrKeys);     
00116     RegisterAt(DisplayCtrl,"DisplayMinLevel",base::_MinLevel); 
00117     RegisterAt(DisplayCtrl,"DisplayMaxLevel",base::_MaxLevel);
00118     RegisterAt(DisplayCtrl,"CutOut",base::_Eliminate);
00119     RegisterAt(DisplayCtrl,"CutOutValue",base::_EliminateValue);
00120     RegisterAt(DisplayCtrl,"RegularMax",RegMax);
00121     RegisterAt(DisplayCtrl,"RegularLimit(1)",RegLim[0]);
00122     RegisterAt(DisplayCtrl,"RegularLimit(2)",RegLim[1]);
00123     RegisterAt(DisplayCtrl,"PointOutput",PointOutput);
00124     RegisterAt(DisplayCtrl,"LevelCutOut",LevelCutOut);
00125     RegisterAt(DisplayCtrl,"FileMerge",FileMerge);
00126     RegisterAt(DisplayCtrl,"InputFile",InputFile);
00127     for (int d=0; d<3; d++) {
00128       char VariableName[16];
00129       std::sprintf(VariableName,"ShowMin(%d)",d+1);
00130       RegisterAt(DisplayCtrl,VariableName,base::show[2*d]);
00131       std::sprintf(VariableName,"ShowMax(%d)",d+1);
00132       RegisterAt(DisplayCtrl,VariableName,base::show[2*d+1]);     
00133       std::sprintf(VariableName,"ShowCut(%d)",d+1);
00134       RegisterAt(DisplayCtrl,VariableName,base::cut[d]);     
00135       std::sprintf(VariableName,"Symmetry(%d)",d+1);
00136       RegisterAt(DisplayCtrl,VariableName,base::Symmetry[d]);     
00137       std::sprintf(VariableName,"Periodic(%d)",d+1);
00138       RegisterAt(DisplayCtrl,VariableName,base::Periodic[d]);     
00139       std::sprintf(VariableName,"LinePoint(%d)",d+1);
00140       RegisterAt(DisplayCtrl,VariableName,LinePoint[d]);     
00141       std::sprintf(VariableName,"LineVector(%d)",d+1);
00142       RegisterAt(DisplayCtrl,VariableName,LineVector[d]);          
00143       std::sprintf(VariableName,"RegularData(%d)",d+1);
00144       RegisterAt(DisplayCtrl,VariableName,RegData[d]);
00145     }   
00146 
00147     ControlDevice WhereCtrl, ToCtrl, FromCtrl;
00148     for (register int l=0; l<MAXLEV; l++) {
00149       char devno[8];
00150       std::sprintf(devno,"(%d)",l);
00151       WhereCtrl = DisplayCtrl.getSubDevice("Where"+std::string(devno));
00152       for (int d=0; d<3; d++) {
00153         char text[8];
00154         std::sprintf(text,"lb(%d)",d+1);
00155         RegisterAt(WhereCtrl,text,base::_where[l].lower(d));
00156         std::sprintf(text,"ub(%d)",d+1);
00157         RegisterAt(WhereCtrl,text,base::_where[l].upper(d));
00158         std::sprintf(text,"s(%d)",d+1);
00159         RegisterAt(WhereCtrl,text,base::_where[l].stepsize(d));
00160       }
00161       ToCtrl = DisplayCtrl.getSubDevice("To"+std::string(devno));
00162       for (int d=0; d<3; d++) {
00163         char text[8];
00164         std::sprintf(text,"lb(%d)",d+1);
00165         RegisterAt(ToCtrl,text,base::_to[l].lower(d));
00166         std::sprintf(text,"ub(%d)",d+1);
00167         RegisterAt(ToCtrl,text,base::_to[l].upper(d));
00168         std::sprintf(text,"s(%d)",d+1);
00169         RegisterAt(ToCtrl,text,base::_to[l].stepsize(d));
00170       }
00171       FromCtrl = DisplayCtrl.getSubDevice("From"+std::string(devno));
00172       for (int d=0; d<3; d++) {
00173         char text[8];
00174         std::sprintf(text,"lb(%d)",d+1);
00175         RegisterAt(FromCtrl,text,base::_from[l].lower(d));
00176         std::sprintf(text,"ub(%d)",d+1);
00177         RegisterAt(FromCtrl,text,base::_from[l].upper(d));
00178         std::sprintf(text,"s(%d)",d+1);
00179         RegisterAt(FromCtrl,text,base::_from[l].stepsize(d));
00180       }
00181     }
00182     RegisterAt(DisplayCtrl,"RestrictData",base::RestrictData);
00183 
00184     TheEvaluator().register_at(DisplayCtrl,"");
00185     DisplayCtrl.update();  
00186 
00187     TheEvaluator().update();
00188   }
00189 
00190   virtual int SetUp(int step, std::string** filenames, int Nfilenames) { 
00191     int res;
00192     if ((OutputType==0 || OutputType==1 || OutputType==20) && !InputFile.empty())
00193       base::_Eliminate = 0;
00194     if (OutputType==6 || OutputType==7) { 
00195       base::_CutOut = 0;
00196       _Step = step;
00197     }
00198     if (OutputType==15)
00199       if (LevelCutOut<=0 && FKeys()!=1)
00200         base::_CutOut = 0;
00201       
00202     if ((res=base::SetUp(step,filenames,Nfilenames)) >= 0) {
00203       char append_str[32], time[20];
00204       base::OutputFileTime(time, step);
00205       if (OutputName != "-") {
00206         if (OutputType>=0 && OutputType<=2) 
00207           std::sprintf(append_str,"_%s.txt",time);
00208         else if (OutputType>=3 && OutputType<=5) 
00209           std::sprintf(append_str,"_%s.dx",time);
00210         else if (OutputType==8) 
00211           std::sprintf(append_str,"_%s.dat",time);
00212         else if (OutputType==10 || OutputType==11) 
00213           std::sprintf(append_str,"_%s.vtk",time);
00214         else if (OutputType==12 || OutputType==13) 
00215           std::sprintf(append_str,"_%s.vtu",time);
00216         else if (OutputType==15) 
00217           std::sprintf(append_str,"_%s.silo",time);
00218         else if (OutputType==20) 
00219           std::sprintf(append_str,"_%s.tec",time);
00220         else if (OutputType==-1) 
00221           std::sprintf(append_str,"_%s.txt",time);
00222         if (OutputType<=5 || OutputType>=8) OutputName += append_str;
00223       }
00224 
00225       int* work = new int[StrKeys.length()];
00226       for (unsigned int sk=0; sk<StrKeys.length(); sk++) 
00227         for (int i=0; i<TheEvaluator().NKeys(); i++) 
00228           if (TheEvaluator().KeyboardKeys()[i] == StrKeys.c_str()[sk]) {
00229             work[nkeys] = i+1; nkeys++; 
00230           }
00231       
00232       if (keys) delete [] keys;
00233       keys = new int[nkeys];
00234       for (int k=0; k<nkeys; k++)
00235         keys[k] = work[k];
00236 
00237       delete work;
00238     }
00239     return res;
00240   }
00241 
00242 
00243   virtual void Write() {
00244     if (OutputType>=6 && OutputType<=7 && OutputName=="-")
00245       OutputName = "out";
00246     if (OutputName != "-") {
00247       if (OutputType <= 5 || ( OutputType >= 8 && OutputType < 15 ) || 
00248           OutputType == 20) {
00249         std::ofstream outfile;
00250         outfile.open(OutputName.c_str(), std::ios::out);
00251         std::ostream* out = new std::ostream(outfile.rdbuf());
00252         if (OutputType == 0) 
00253           WriteASCIITabular(0,0,*out);
00254         else if (OutputType == 1) 
00255           WriteASCIITabular(1,0,*out); 
00256         else if (OutputType == 2) 
00257           WriteRegularTabular(*out);
00258         else if (OutputType == 3) 
00259           WriteDXFile(*out);
00260         else if (OutputType == 4) 
00261           WriteDXASCIIFile(*out);
00262         else if (OutputType == 5) 
00263           WriteDXExternalFilter(*out);
00264         else if (OutputType == 8) 
00265           WriteTecplotASCII(*out);
00266         else if (OutputType == 10) 
00267           WriteVTKASCIIFile(*out);
00268         else if (OutputType == 11) 
00269           WriteVTKFile(*out);
00270         else if (OutputType == 12)
00271           WriteVTUASCIIFile(*out);
00272         else if (OutputType == 13)
00273           WriteVTUFile(*out);
00274         else if (OutputType == 20)
00275           WriteTecplotASCIIMeshFile(*out);
00276         else if (OutputType == -1)
00277           WriteMeshSweep(*out);
00278         outfile.close();
00279         delete out;
00280       }
00281       else if (OutputType == 6) 
00282         WriteHDFConvert(0);
00283       else if (OutputType == 7) 
00284         WriteHDFConvert(1);
00285       else if (OutputType == 15)
00286         WriteSiloFile((char *)OutputName.c_str());
00287     }
00288     else {
00289       if (OutputType == 0) 
00290         WriteASCIITabular(0,0,std::cout);
00291       else if (OutputType == 1) 
00292         WriteASCIITabular(1,0,std::cout); 
00293       else if (OutputType == 2) 
00294         WriteRegularTabular(std::cout);
00295       else if (OutputType == 3) 
00296         WriteDXFile(std::cout);
00297       else if (OutputType == 4) 
00298         WriteDXASCIIFile(std::cout);
00299       else if (OutputType == 5) 
00300         WriteDXExternalFilter(std::cout);
00301       else if (OutputType == 8) 
00302         WriteTecplotASCII(std::cout);
00303       else if (OutputType == 10)
00304         WriteVTKASCIIFile(std::cout);
00305       else if (OutputType == 11)
00306         WriteVTKFile(std::cout);
00307       else if (OutputType == 12) 
00308         WriteVTUASCIIFile(std::cout);
00309       else if (OutputType == 13) 
00310         WriteVTUFile(std::cout);
00311       else if (OutputType == 15)
00312         WriteSiloFile(0);
00313       else if (OutputType == 20)
00314         WriteTecplotASCIIMeshFile(std::cout);
00315       else if (OutputType == -1)
00316         WriteMeshSweep(std::cout);
00317     }
00318   }
00319 
00320   void WriteMeshSweep(std::ostream& out) {
00321     if (!keys) return;
00322     int length;
00323     if (FKeys() == 1) 
00324       length = NNodes();
00325     else
00326       length = NCells();
00327     
00328     float* dat = new float[length];
00329     for (int k=0; k<nkeys; k++) {
00330       for (register int i=0; i<length; i++) *(dat+i) = 0.;
00331       SetScal(&(keys[k]),dat);
00332       float min=FLT_MAX, max=FLT_MIN, sum=0.;
00333       for (register int i=0; i<length; i++) {
00334         min = Min(min,*(dat+i));
00335         max = Max(max,*(dat+i));
00336         sum += *(dat+i);
00337       }
00338       out << min << " " << max << " " << sum << std::endl;
00339     }
00340     delete [] dat;
00341   }
00342 
00343   void WriteRegularTabular(std::ostream& out) {
00344     if (!keys) return;
00345     int length;
00346     if (FKeys() == 1) 
00347       length = NNodes();
00348     else
00349       length = NCells();
00350     
00351     float* xyz = new float[3*length];
00352     SetGrid((float (*)[3]) xyz); 
00353     float* dat = new float[nkeys*length];
00354     float* val = new float[nkeys];
00355     for (int i=0; i<nkeys*length; i++) *(dat+i) = 0.;
00356     for (int k=0; k<nkeys; k++) 
00357       SetScal(&(keys[k]),&(dat[k*length]));
00358 
00359     float regdx[3], reggeom[6];
00360     int xd = -1, yd = -1, fd = -1;
00361     for (int d=0; d<3; d++) {
00362       if (base::DimUsed(d)==0 || base::show[2*d]>=base::show[2*d+1]) {
00363         reggeom[2*d] = base::geom[2*d]; reggeom[2*d+1] = base::geom[2*d+1]; 
00364       }
00365       else {
00366         reggeom[2*d] = Max(base::show[2*d],base::geom[2*d]); 
00367         reggeom[2*d+1] = Min(base::show[2*d+1],base::geom[2*d+1]); 
00368       }
00369       if (RegData[d]>0) regdx[d] = (reggeom[2*d+1] - reggeom[2*d])/RegData[d];
00370       else regdx[d] = 0.0;
00371       if (yd<0 && xd>=0 && regdx[d]>0.0) yd = d;
00372       if (xd<0 && regdx[d]>0.0) xd = d;
00373       if (fd<0 && regdx[d]==0.0) fd = d;
00374     }
00375     float xs = reggeom[2*xd] + (FKeys()==1 ? 0.0 : regdx[xd]/2.0);
00376     while (xs<=reggeom[2*xd+1] && regdx[xd]>0.0) {
00377       float ys = reggeom[2*yd] + (FKeys()==1 ? 0.0 : regdx[yd]/2.0);
00378       while (ys<=reggeom[2*yd+1] && regdx[yd]>0.0) {
00379         float mindist=FLT_MAX;
00380         float fddist=(RegMax ? FLT_MIN : FLT_MAX);
00381         int inearest = -1;
00382         for (register int i=0; i<length; i++) {
00383           float sum=std::sqrt((xs-xyz[3*i+xd])*(xs-xyz[3*i+xd])+
00384                               (ys-xyz[3*i+yd])*(ys-xyz[3*i+yd]));
00385           if (sum<mindist-PROBE_TOLERANCE*(mindist>0?mindist:1.0) ||
00386               inearest<0) {
00387             mindist = sum;
00388             inearest = i;
00389           }
00390           if (fd>=0 && sum<mindist+PROBE_TOLERANCE*(mindist>0?mindist:1.0)) {
00391             if (RegMax && xyz[3*i+fd]>fddist &&
00392                 dat[i]>=RegLim[0] && dat[i]<=RegLim[1]) {
00393               fddist = xyz[3*i+fd]; inearest = i;
00394             }
00395             if (!RegMax && xyz[3*i+fd]<fddist &&
00396                 dat[i]>=RegLim[0] && dat[i]<=RegLim[1]) {
00397               fddist = xyz[3*i+fd]; inearest = i;
00398             }
00399           }
00400         }
00401         out << xs << " " << ys;
00402         for (register int k=0; k<nkeys; k++) 
00403           out << " " << dat[k*length+inearest];
00404         out << std::endl;
00405         std::cerr << ".";
00406         ys += regdx[yd];
00407       }
00408       out << std::endl;
00409       std::cerr << std::endl;
00410       xs += regdx[xd];
00411     }
00412 
00413     delete [] val;
00414     delete [] dat;
00415     delete [] xyz;
00416   }
00417 
00418   void WriteASCIITabular(int ordered, int tecplot, std::ostream& out) {
00419     if (!keys) return;
00420 
00421     if (InputFile.empty()) {
00422       int length;
00423       if (FKeys() == 1)
00424         length = NNodes();
00425       else
00426         length = NCells();
00427       
00428       float* xyz = new float[3*length];
00429       SetGrid((float (*)[3]) xyz);  
00430       int nkl = nkeys*length;
00431       float* dat = new float[nkl];
00432       for (register int i=0; i<nkl; i++) *(dat+i) = 0.;
00433       register int k;
00434       for (k=0; k<nkeys; k++) 
00435         SetScal(&(keys[k]),&(dat[k*length]));
00436       
00437       if (!tecplot) {
00438         out << RealTime() << ";";
00439         if (!LineProbe()) {
00440           if (!ordered) {
00441             if (base::DimUsed(0)) out << "x ;";
00442             if (base::DimUsed(1)) out << "y ;";
00443             if (base::DimUsed(2)) out << "z ;";
00444           }
00445           else
00446             if (base::DimUsed(0)) out << "x ;";
00447             else if (base::DimUsed(1)) out << "y ;";
00448             else if (base::DimUsed(2)) out << "z ;";
00449         }
00450         else
00451           out << "r ;";
00452         for (k=0; k<nkeys; k++) 
00453           out << TheEvaluator().TitleKeys(keys[k]-1) << ";";
00454       }
00455       else {
00456         out << "title = \"t = " << RealTime() << "\"" << std::endl;
00457         out << "variables = ";
00458         if (!LineProbe()) out << "\"x\"";
00459         else out << "\"r\"";
00460         for (k=0; k<nkeys; k++) 
00461           out << ", \"" << TheEvaluator().TitleKeys(keys[k]-1) << "\"";
00462       }
00463       out << std::endl;
00464    
00465       if (!ordered) {
00466         for (register int i=0; i<length; i++) {
00467           if (!LineProbe()) {
00468             for (register int d=0; d<3; d++) 
00469               if (base::DimUsed(d))
00470                 out << xyz[3*i+d] << " ";        
00471             for (k=0; k<nkeys; k++) 
00472                 out << dat[k*length+i] << " ";
00473             out << std::endl;
00474           }
00475           else {
00476             int lp = 1, d;
00477             float di = 0.0, xs = 0.0, c = 0.0;
00478             for (d=0; d<3; d++) 
00479               if (base::DimUsed(d)) {
00480                 di += (xyz[3*i+d]-LinePoint[d])*(xyz[3*i+d]-LinePoint[d]);
00481                 if (LineVector[d] != 0.0) {
00482                   xs += (xyz[3*i+d]-LinePoint[d])/LineVector[d]; c += 1.0;
00483                 }
00484               }   
00485             c = xs/c;
00486             for (d=0; d<3; d++) 
00487               if (base::DimUsed(d) && 
00488                   std::fabs(c*LineVector[d]-xyz[3*i+d]+LinePoint[d])>
00489                   PROBE_TOLERANCE) {
00490                 lp = 0; break;
00491               }
00492             if (lp) {
00493               out << std::sqrt(di) << " ";        
00494               for (k=0; k<nkeys; k++) 
00495                 out << dat[k*length+i] << " ";
00496               out << std::endl;
00497             }
00498           }
00499         }
00500       }
00501       else {
00502         index_ordered_map_type imap;
00503         for (register int i=0; i<length; i++) 
00504           if (!LineProbe()) {
00505             for (register int d=0; d<3; d++)
00506                 if (base::DimUsed(d)) {
00507                   imap.insert(ordered_index_type(xyz[3*i+d],i));
00508                   break;
00509                 }
00510           }
00511           else {
00512             int lp = 1, d;
00513             float di = 0.0, xs = 0.0, c = 0.0;
00514             for (d=0; d<3; d++) 
00515               if (base::DimUsed(d)) {
00516                 di += (xyz[3*i+d]-LinePoint[d])*(xyz[3*i+d]-LinePoint[d]);
00517                 if (LineVector[d] != 0.0) {
00518                   xs += (xyz[3*i+d]-LinePoint[d])/LineVector[d]; c += 1.0;
00519                 }
00520               }
00521             c = xs/c;
00522             for (d=0; d<3; d++) 
00523               if (base::DimUsed(d) && 
00524                   std::fabs(c*LineVector[d]-xyz[3*i+d]+LinePoint[d])>
00525                   PROBE_TOLERANCE) {
00526                 lp = 0; break;
00527               }
00528             if (lp) 
00529               imap.insert(ordered_index_type(std::sqrt(di),i));
00530           }
00531         
00532         for (index_ordered_map_type::iterator it = imap.begin();
00533              it != imap.end(); ++it) {
00534           out << (*it).first << " ";
00535           for (k=0; k<nkeys; k++) 
00536             out << dat[k*length+(*it).second] << " ";
00537           out << std::endl;
00538         }
00539       }
00540     
00541       delete [] dat;
00542       delete [] xyz;
00543     }    
00544     else {
00545       std::ifstream in(InputFile.c_str(), std::ios::in);
00546       if (!in) {
00547         std::cerr << "Cannot open " << InputFile << " for reading!" << std::endl;
00548         return;
00549       }
00550       float SP[3];
00551       int d, k;
00552       point_list_type p;
00553       while (!in.eof()) {
00554         for (d=0; d<3; d++) 
00555           if (base::DimUsed(d))
00556             in >> SP[d];
00557           else SP[d]=0.0;
00558         if (in.eof()) continue;
00559         p.insert(p.end(), SP, SP+3);
00560       }
00561       in.close();
00562       if (!p.empty()) {
00563         unsigned int Np=p.size()/3;
00564         float *v=new float[Np*nkeys];
00565         for (k=0; k<nkeys; k++) 
00566           if (FKeys() == 1)
00567             base::SetPointsNodes(&(keys[k]),p,&(v[k*Np]));
00568           else
00569             base::SetPointsCells(&(keys[k]),p,&(v[k*Np]));
00570         
00571         for (register unsigned int i=0;i<Np;i++) {
00572           bool print=true;
00573           for (k=0; k<nkeys; k++) 
00574             if (v[k*Np+i]<=float(-1.e37) && !ordered) { print=false; break; }
00575           if (print) {
00576             for (d=0; d<3; d++) 
00577               if (base::DimUsed(d))
00578                 out << p[3*i+d] << " ";
00579             for (k=0; k<nkeys; k++) 
00580               out << v[k*Np+i] << " ";
00581             out << std::endl;
00582           }
00583           else {
00584             std::cerr << "Point ( ";
00585             for (d=0; d<3; d++) 
00586               if (base::DimUsed(d))
00587                 std::cerr << p[3*i+d] << " ";
00588             std::cerr << ") outside of domain. Skipping..." << std::endl; 
00589           }
00590         }
00591         delete [] v;
00592       }
00593     }
00594   }
00595   
00596   void DimShape(int& dimshape, int& conshape) {
00597     if (base::Make3d) {
00598       dimshape = 3;
00599       conshape = 8;
00600     }
00601     else {
00602       dimshape = 0;
00603       conshape = 1;
00604       for (int d=0; d<3; d++) {
00605         dimshape += base::DimUsed(d);
00606         conshape *= (base::DimUsed(d) ? 2 : 1);
00607       }
00608     }
00609   }
00610 
00611   //  
00612   // All nodes are dumped out. Equal nodes in different subgrids are NOT deleted, but
00613   // DX does not seem to be disturbed.
00614   //
00615   void WriteDXFile(std::ostream& out) {
00616     int dxfile_dataoffset=0;
00617     int dimshape, conshape;
00618     DimShape(dimshape,conshape);
00619 
00620     int nnodes = NNodes(), ncells = NCells();    
00621     float* xyz = new float[3*nnodes];
00622     out << "object 1 class array type float rank 1 shape " << dimshape;
00623     out << " items " << nnodes;
00624 #if defined DX_LSB
00625     out << " lsb";
00626 #else
00627     out << " msb";
00628 #endif   
00629     out << " binary data " << dxfile_dataoffset << std::endl;
00630     SetGridNodes((float (*)[3]) xyz);
00631     dxfile_dataoffset+=dimshape*nnodes*sizeof(float);
00632 
00633     int* con = new int[8*ncells];    
00634     out << "object 2 class array type int rank 1 shape " << conshape;
00635     out << " items " << ncells;
00636 #if defined DX_LSB
00637     out << " lsb";
00638 #else
00639     out << " msb";
00640 #endif   
00641     out << " binary data " << dxfile_dataoffset << std::endl;
00642     SetGridConns((int (*)[8]) con);  
00643     out << "attribute \"element type\" string ";
00644     if (dimshape==3) out << "\"cubes\"";
00645     else if (dimshape==2) out << "\"quads\"";
00646     else if (dimshape==1) out << "\"lines\"";
00647 
00648     out << std::endl << "attribute \"ref\" string \"positions\"" << std::endl;
00649     dxfile_dataoffset+=conshape*ncells*sizeof(int);
00650 
00651     int length;
00652     if (FKeys() == 1)
00653       length = nnodes;
00654     else
00655       length = ncells;
00656     int nkl = nkeys*length;
00657     float* dat = new float[nkl];
00658     for (register int i=0; i<nkl; i++) *(dat+i) = 0.;
00659     out << "object 3 class array type float";
00660     if (nkeys>1) out << " rank 1 shape " << nkeys;
00661     out << " items " << length;
00662 #if defined DX_LSB
00663     out << " lsb";
00664 #else
00665     out << " msb";
00666 #endif   
00667     out << " binary data " << dxfile_dataoffset << std::endl;
00668     register int k;
00669     for (k=0; k<nkeys; k++) 
00670       SetScal(&(keys[k]),&(dat[k*length]));
00671     if (FKeys() == 6)
00672       out << "attribute \"dep\" string \"connections\"" << std::endl;
00673     dxfile_dataoffset+= nkeys*length*sizeof(float);
00674 
00675     out << "object \"default\" class field" << std::endl;
00676     out << "component \"positions\" value 1" << std::endl;
00677     out << "component \"connections\" value 2" << std::endl;
00678     out << "component \"data\" value 3" << std::endl;
00679     out << "end" << std::endl;
00680     
00681     register int n,t;
00682     if (dimshape==1) {
00683       for (n=0,t=0; n<nnodes; n++)
00684         xyz[t++] = xyz[3*n];
00685     }
00686     else if (dimshape==2) {
00687       for (n=0,t=0; n<nnodes; n++) {
00688         xyz[t++] = xyz[3*n]; xyz[t++] = xyz[3*n+1];
00689       }
00690     }
00691     out.write((char*) xyz, dimshape*nnodes*sizeof(float));
00692 
00693     if (conshape==2) {
00694       for (n=0,t=0; n<ncells; n++) {
00695         con[t++] = con[8*n]; con[t++] = con[8*n+1];
00696       }
00697     }
00698     if (conshape==4) {
00699       for (n=0,t=0; n<ncells; n++) {
00700         con[t++] = con[8*n  ]; con[t++] = con[8*n+1]; 
00701         con[t++] = con[8*n+2]; con[t++] = con[8*n+3];
00702       }
00703     }
00704     out.write((char*) con, conshape*ncells*sizeof(int));
00705 
00706     if (nkeys>1) {
00707       float* datout = new float[nkeys*length];
00708       int nn = 0;
00709       for (n=0; n<length; n++)
00710         for (k=0; k<nkeys; k++) {
00711           datout[nn] = dat[k*length+n]; nn++;
00712         }
00713       out.write((char*) datout, sizeof(float)*nkeys*length);
00714       delete [] datout;
00715     }
00716     else 
00717       out.write((char*) dat, sizeof(float)*length);
00718       
00719     delete [] xyz;
00720     delete [] con;
00721     delete [] dat;
00722   }    
00723 
00724   void WriteDXASCIIFile(std::ostream& out) {
00725     int dimshape, conshape;
00726     DimShape(dimshape,conshape);
00727     
00728     int nnodes = NNodes(), ncells = NCells();
00729     float* xyz = new float[3*nnodes];
00730     out << "object 1 class array type float rank 1 shape " << dimshape;
00731     out << " items " << nnodes << " data follows" << std::endl;
00732     SetGridNodes((float (*)[3]) xyz);
00733     register int i;
00734     for (i=0; i<nnodes; i++) {
00735       for (int d=0; d<dimshape; d++) 
00736         out << xyz[3*i+d] << " ";
00737       out << std::endl;
00738     }
00739 
00740     int* con = new int[8*ncells];    
00741     out << "object 2 class array type int rank 1 shape " << conshape;
00742     out << " items " << ncells << " data follows" << std::endl;
00743     SetGridConns((int (*)[8]) con);  
00744     for (i=0; i<ncells; i++) {
00745       for (int d=0; d<conshape; d++) 
00746         out << con[8*i+d] << " ";
00747       out << std::endl;
00748     }
00749 
00750     out << "attribute \"element type\" string ";
00751     if (dimshape==3) out << "\"cubes\"";
00752     else if (dimshape==2) out << "\"quads\"";
00753     else if (dimshape==1) out << "\"lines\"";
00754 
00755     out << std::endl << "attribute \"ref\" string \"positions\"" << std::endl;
00756     int length;
00757     if (FKeys() == 1)
00758       length = nnodes;
00759     else
00760       length = ncells;
00761     float* dat = new float[nkeys*length];
00762     for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
00763     out << "object 3 class array type float";
00764     if (nkeys>1) out << " rank 1 shape " << nkeys;
00765     out << " items " << length << " data follows" << std::endl;
00766     register int k;
00767     for (k=0; k<nkeys; k++) 
00768       SetScal(&(keys[k]),&(dat[k*length]));
00769     for (i=0; i<length; i++) {
00770       for (k=0; k<nkeys; k++) 
00771         out << " " << dat[k*length+i];
00772       out << std::endl;
00773     }
00774     if (FKeys() == 6)
00775       out << "attribute \"dep\" string \"connections\"" << std::endl;
00776 
00777     out << "object \"default\" class field" << std::endl;
00778     out << "component \"positions\" value 1" << std::endl;
00779     out << "component \"data\" value 3" << std::endl;
00780     out << "component \"connections\" value 2" << std::endl;
00781     out << "end" << std::endl;
00782         
00783     delete [] xyz;
00784     delete [] con;
00785     delete [] dat;
00786   }    
00787 
00788   void WriteDXExternalFilter(std::ostream& out) {
00789     int dimshape, conshape;
00790     DimShape(dimshape,conshape);
00791     
00792     int nnodes = NNodes(), ncells = NCells();
00793     float* xyz = new float[3*nnodes];
00794     out << "object 1 class array type float rank 1 shape " << dimshape;
00795     out << " items " << nnodes;
00796 #if defined DX_LSB
00797     out << " lsb";
00798 #else
00799     out << " msb";
00800 #endif   
00801     out << " binary data follows" << std::endl;
00802     SetGridNodes((float (*)[3]) xyz);
00803     register int n,t;
00804     if (dimshape==1) {
00805       for (n=0,t=0; n<nnodes; n++)
00806         xyz[t++] = xyz[3*n];
00807     }
00808     else if (dimshape==2) {
00809       for (n=0,t=0; n<nnodes; n++) {
00810         xyz[t++] = xyz[3*n]; xyz[t++] = xyz[3*n+1];
00811       }
00812     }
00813     out.write((char*) xyz, dimshape*nnodes*sizeof(float));
00814     delete [] xyz;
00815 
00816     int* con = new int[8*ncells];    
00817     out << "object 2 class array type int rank 1 shape " << conshape;
00818     out << " items " << ncells;
00819 #if defined DX_LSB
00820     out << " lsb";
00821 #else
00822     out << " msb";
00823 #endif   
00824     out << " binary data follows" << std::endl;
00825     SetGridConns((int (*)[8]) con);  
00826     if (conshape==2) {
00827       for (n=0,t=0; n<ncells; n++) {
00828         con[t++] = con[8*n]; con[t++] = con[8*n+1];
00829       }
00830     }
00831     if (conshape==4) {
00832       for (n=0,t=0; n<ncells; n++) {
00833         con[t++] = con[8*n  ]; con[t++] = con[8*n+1]; 
00834         con[t++] = con[8*n+2]; con[t++] = con[8*n+3];
00835       }
00836     }
00837     out.write((char*) con, conshape*ncells*sizeof(int));
00838     delete [] con;
00839 
00840     out << "attribute \"element type\" string ";
00841     if (dimshape==3) out << "\"cubes\"";
00842     else if (dimshape==2) out << "\"quads\"";
00843     else if (dimshape==1) out << "\"lines\"";
00844 
00845     out << "attribute \"ref\" string \"positions\"" << std::endl;
00846 
00847     int length;
00848     if (FKeys() == 1)
00849       length = nnodes;
00850     else
00851       length = ncells;
00852     float* dat = new float[nkeys*length];
00853     for (int i=0; i<nkeys*length; i++) *(dat+i) = 0.;
00854     out << "object 3 class array type float ";
00855     if (nkeys>1) out << "rank 1 shape " << nkeys << " ";
00856     if (FKeys() == 6)
00857       out << "attribute \"dep\" string \"connections\"" << std::endl;
00858     out << "items " << length;
00859 #if defined DX_LSB
00860     out << " lsb";
00861 #else
00862     out << " msb";
00863 #endif   
00864     out << " binary data follows" << std::endl;
00865     register int k;
00866     for (k=0; k<nkeys; k++) 
00867       SetScal(&(keys[k]),&(dat[k*length]));
00868     if (nkeys>1) {
00869       float* datout = new float[nkeys*length];
00870       int nn = 0;
00871       for (n=0; n<length; n++)
00872         for (k=0; k<nkeys; k++) {
00873           datout[nn] = dat[k*length+n]; nn++;
00874         }
00875       out.write((char*) datout, sizeof(float)*nkeys*length);
00876       delete [] datout;
00877     }
00878     else 
00879       out.write((char*) dat, sizeof(float)*length);      
00880     delete [] dat;
00881 
00882     out << "object \"default\" class field" << std::endl;
00883     out << "component \"positions\" value 1" << std::endl;
00884     out << "component \"connections\" value 2" << std::endl;
00885     out << "component \"data\" value 3" << std::endl;
00886     out << "end" << std::endl;  
00887   }    
00888 
00889   void WriteHDFConvert(int DimReduce) {
00890     if (!keys || _Step<0 || OutputName == "-") return;
00891     float64 realtime = RealTime();
00892     int32 ori[3],dims[3],res[4];
00893     int32 proc,level,gridID;
00894     int32 step = _Step;
00895     int32 rank = 3;
00896     char target[256];
00897     char dname[16];
00898     
00899     if (BlockList().size()>MAX_NC_VARS && FileMerge>=0) {
00900       std::cerr << "Maximal block count in HDF4 file exceeded. Setting FileMerge to 0." << std::endl;
00901       FileMerge = 0;
00902     }
00903 
00904     if (OutputStep>=0) step = OutputStep;
00905 
00906     for (register int k=0; k<nkeys; k++) {
00907       std::string TargetName = OutputName;
00908       dname[0] = TheEvaluator().KeyboardKeys()[keys[k]-1]; dname[1] = 0;
00909       char time[20];
00910       base::OutputFileTime(time, step);
00911       if (nkeys>1)
00912         std::sprintf(target,"%s_%s_%s.hdf",OutputName.c_str(),dname,time);   
00913       else
00914         std::sprintf(target,"%s_%s.hdf",OutputName.c_str(),time);   
00915       
00916       std::cerr << "SetScalCells()";
00917       gridID=0;
00918       int last_proc = -1;
00919       for (typename block_list_type::reverse_iterator bit=BlockList().rbegin();
00920            bit!=BlockList().rend(); bit++, gridID++) {
00921 
00922         proc = (**bit).Processor();
00923         level = (**bit).Level();
00924         
00925         if (base::_parData==1 && FileMerge<=0) 
00926           if (proc!=last_proc) {
00927             if (last_proc>=0) {
00928               SDSflush(target);
00929               SDSclose(target);
00930             }
00931             if (nkeys>1)
00932               std::sprintf(target,"%s_%s_%s.%d.hdf",OutputName.c_str(),dname,time,(int)proc);   
00933             else
00934               std::sprintf(target,"%s_%s.%d.hdf",OutputName.c_str(),time,(int)proc); 
00935             last_proc = proc;
00936           }
00937 
00938         int offset=0;
00939         hdata_block_type tmpblock((**bit).DB().bbox());
00940         float* work = new DataType[(**bit).DB().bbox().size()];
00941 
00942         (**bit).SetScalCells(&(keys[k]),work,offset,base::CompRead);    
00943         BeginFastIndex3(tmp, tmpblock.bbox(), tmpblock.data(), DataType);
00944         int idx = 0;
00945         for_3 (n,m,l,tmpblock.bbox(),tmpblock.bbox().stepsize())
00946           FastIndex3(tmp,n,m,l) = work[idx];
00947           idx++;              
00948         end_for
00949         EndFastIndex3(tmp); 
00950         delete [] work;
00951 
00952         register int d;
00953         for (d=0; d<3; d++) {
00954           ori[d] = (int32) tmpblock.lower(d);
00955           if (!base::show_exact[d]) res[d] = (int32) tmpblock.stepsize(d);
00956           else res[d] = 1;
00957           dims[d] = (int32) tmpblock.extents(d);
00958         }
00959         res[3] = 1;
00960         if (!base::DimUsed(2)) {
00961           dims[2] = 1;
00962           ori[2] = 0;
00963           res[2] = 1;
00964           res[3] = 1;
00965           if (!base::DimUsed(1)) {
00966             dims[1] = 1;
00967             ori[1] = 0;
00968             res[1] = 1;
00969           }
00970         }
00971         if (DimReduce) {
00972           int dc = 0;
00973           for (d=0; d<3; d++) 
00974             if (base::show_exact[d]) {
00975               for (register int dd=d+1; dd<3; dd++) {
00976                 ori[dd-1] = ori[dd];
00977                 res[dd-1] = res[dd];
00978                 dims[dd-1] = dims[dd]; 
00979               }
00980               dc++;
00981               ori[3-dc] = 0; res[3-dc] = 1; dims[3-dc] = 1; 
00982             }
00983         }       
00984 
00985         if (sizeof(DataType) == SIZE_FLOAT32)
00986           SDSsetNT(target,DFNT_FLOAT32);    
00987         else if (sizeof(DataType) == SIZE_FLOAT64)
00988           SDSsetNT(target,DFNT_FLOAT64);    
00989         else if (sizeof(DataType) == SIZE_FLOAT128)
00990           SDSsetNT(target,DFNT_FLOAT128); 
00991         else
00992           std::cerr << "SDSsetNT() not called!" << std::endl;
00993         AMRwriteData(target,dname,proc,level,gridID,
00994                      step,realtime,rank,dims,ori,res,
00995                      tmpblock.data());
00996       }
00997       SDSflush(target);
00998       SDSclose(target);
00999       std::cerr << std::endl;
01000     }
01001   }    
01002   
01003   void WriteTecplotASCII(std::ostream& out) {
01004     int dimshape, conshape;
01005     DimShape(dimshape,conshape);
01006 
01007     if (dimshape==1) { 
01008       WriteASCIITabular(1,1,out); 
01009       return;
01010     }
01011     
01012     int length;
01013     if (FKeys() == 1)
01014       length = NNodes();
01015     else
01016       length = NCells();
01017     out << "title = \"t = " << RealTime() << "\"" << std::endl;
01018     float* xyz = new float[3*length];
01019     SetGrid((float (*)[3]) xyz);
01020 
01021     out << "variables = \"x\", \"y\"";
01022     if (dimshape>2) out << ", \"z\"";
01023     register int i,k;
01024     for (k=0; k<nkeys; k++) {
01025       char* title = new char[LEN_TKEYS];
01026       std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
01027       i=std::strlen(title)-1;
01028       while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }      
01029       out << ", \"" << title << "\"";
01030     }
01031     out << std::endl;
01032     float* dat = new float[nkeys*length];
01033     for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
01034     for (k=0; k<nkeys; k++) 
01035       SetScal(&(keys[k]),&(dat[k*length]));
01036 
01037     int cnt=0,cbit=0,blength;
01038     for (typename block_list_type::iterator bit=BlockList().begin();
01039          bit!=BlockList().end(); bit++, cbit++) {
01040 
01041       Coords ext = (*bit)->DB().extents();
01042       if (FKeys() == 1) {
01043         ext += base::DimUsed; 
01044         blength = (*bit)->NNodes();
01045       }
01046       else
01047         blength = (*bit)->NCells();
01048 
01049       out << "zone t=\"B-" << cbit << " L-" << (*bit)->Level() << " P-" 
01050           << (*bit)->Processor() << "\", i=" << ext(0) << ", j=" << ext(1);
01051       if (dimshape>2) out << ", k=" << ext(2);
01052 
01053       out << ", solutiontime=" << RealTime(); 
01054       if (PointOutput<=0) {
01055         out << ", datapacking=block, varlocation=([1-" << 2+(dimshape>2?1:0)+nkeys << "]=";
01056         if (FKeys() == 1) 
01057           out << "nodal)" << std::endl;
01058         else
01059           out << "cellcentered)" << std::endl;    
01060         for (register int d=0; d<dimshape; d++)
01061           for (i=0; i<blength; i++) {
01062             out << xyz[3*(cnt+i)+d] << " ";
01063             if ((i+1)%16==0) out << std::endl; 
01064           }
01065         out << std::endl << std::endl;
01066         for (k=0; k<nkeys; k++) {
01067           for (i=0; i<blength; i++) {
01068             out << dat[k*length+cnt+i] << " ";
01069             if ((i+1)%16==0) out << std::endl; 
01070           }
01071           out << std::endl << std::endl;
01072         }
01073       }
01074       else {
01075         out << ", datapacking=point, varlocation=([1-" << 2+(dimshape>2?1:0)+nkeys << "]=";
01076         if (FKeys() == 1) 
01077           out << "nodal)" << std::endl;
01078         else
01079           out << "cellcentered)" << std::endl;    
01080         for (i=0; i<blength; i++) {
01081           for (register int d=0; d<dimshape; d++) 
01082             out << xyz[3*(cnt+i)+d] << " ";
01083           for (k=0; k<nkeys; k++) 
01084             out << " " << dat[k*length+cnt+i];
01085           out << std::endl;
01086         }
01087         out << std::endl;
01088       }
01089       cnt += blength;
01090     }
01091         
01092     delete [] xyz;
01093     delete [] dat;
01094   }
01095 
01096 
01097   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01098   //                                            WRITE VTU ASCII FILE
01099   //
01100   // this function converts the amroc dumped hdf files (in memory)
01101   // into an xml based VTK unstructured grid format. notice that the
01102   // levels have been flattened. it is quite inefficient.
01103   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01104 # define VTK_USED_DIMS 3
01105   void WriteVTUASCIIFile(std::ostream& out)  {
01106     // VTI ASCI File writer. based of VTKASCII writer.
01107     // http://www.cacr.caltech.edu/~slombey/asci/vtk/
01108     int dimshape, conshape;
01109     DimShape(dimshape,conshape);
01110 
01111     out << "<?xml version=\"1.0\"?>" << std::endl;
01112     out << "<?meta name=\"source\" content=\"HDF data created with AMROC's hdf2file\"?>" 
01113         << std::endl;
01114     out << "<?meta name=\"timestamp\" content=\"" << RealTime() << "\"?>" << std::endl;
01115 
01116     // VTU XML HEADER
01117     int nnodes = NNodes(), ncells = NCells();
01118     float* xyz = new float[3*nnodes];
01119     out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" 
01120         << std::endl;
01121     out << "<UnstructuredGrid>" << std::endl;
01122     out << "<Piece NumberOfPoints=\"" << nnodes << "\" NumberOfCells=\"" << ncells << "\">" 
01123         << std::endl;
01124 
01125     // flatten points (i.e. remove duplicates) and print
01126     out << "<Points>" << std::endl;
01127     out << "<DataArray type=\"Float32\" NumberOfComponents=\""<< VTK_USED_DIMS 
01128         << "\" format=\"ascii\" Name=\"vertex\">" << std::endl;
01129     SetGridNodes((float (*)[3]) xyz);
01130     register int i;
01131     for (i=0; i<nnodes; i++) {
01132       for (register int d=0; d<VTK_USED_DIMS; d++) 
01133         out << xyz[3*i+d] << " ";
01134       out << std::endl;
01135     }
01136     out << "</DataArray>" << std::endl;
01137     out << "</Points>" << std::endl;
01138 
01139     // cells
01140     int* con = new int[8*ncells];    
01141     out << "<Cells>" << std::endl;
01142     out << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
01143     // print cells->connectivities
01144     SetGridConns((int (*)[8]) con);  
01145     for (i=0; i<ncells; i++) {
01146       for (register int d=0; d<conshape; d++) 
01147         out << con[8*i+d] << " ";
01148       out << std::endl;
01149     }
01150     out << "</DataArray>" << std::endl;
01151 
01152     // print cells->element offsets
01153     out << "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
01154     for (i=1; i<=ncells; i++)
01155       out << i*conshape << std::endl;
01156     out << "</DataArray>" << std::endl;
01157     
01158     // print cells->types
01159     int cell_type = 1;
01160     if (dimshape==3) cell_type = 11;
01161     else if (dimshape==2) cell_type = 8;
01162     else if (dimshape==1) cell_type = 3;
01163 
01164     out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">" << std::endl;
01165     for (i=0; i<ncells; i++) 
01166       out << cell_type << std::endl;
01167     out << "</DataArray>" << std::endl;
01168 
01169     out << "</Cells>" << std::endl;
01170     
01171     int length;
01172     if (FKeys() == 1) {
01173       length = nnodes;
01174       out << "<PointData>" << std::endl;
01175     }
01176     else {
01177       length = ncells;
01178       out << "<CellData>" << std::endl;
01179     }
01180     float* dat = new float[nkeys*length];
01181     for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
01182     register int k;
01183     for (k=0; k<nkeys; k++) 
01184       SetScal(&(keys[k]),&(dat[k*length]));
01185     for (k=0; k<nkeys; k++) {
01186       char* title = new char[LEN_TKEYS];
01187       // clean up title string
01188         std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
01189         i=std::strlen(title)-1;
01190         while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
01191         int n=std::strlen(title);
01192         for (i=0; i<n; i++) if (title[i]==' ') title[i]='_';
01193       out << "<DataArray type=\"Float32\" Name=\"" << title 
01194           << "\" format=\"ascii\" NumberOfComponents=\"1\">" << std::endl;
01195       for (i=0; i<length; i++) 
01196         out << dat[k*length+i] << std::endl; 
01197       out << "</DataArray>" << std::endl;
01198     }
01199 
01200     if (FKeys() == 1)
01201       out << "</PointData>" << std::endl;
01202     else 
01203       out << "</CellData>" << std::endl;
01204     out << "</Piece>" << std::endl;
01205     out << "</UnstructuredGrid>" << std::endl;
01206     out << "</VTKFile>" << std::endl;
01207  
01208     delete [] xyz;
01209     delete [] con;
01210     delete [] dat;
01211   }
01212 
01213   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01214   //                                                     WRITE VTU FILE
01215   //
01216   // this function converts the amroc dumped hdf files (in memory)
01217   // into an xml based VTK unstructured grid format. the data is dumped
01218   // in binary after being 64base converted. (see mime). an extra integer
01219   // is added at the beginning of each binary array with the length.
01220   // + some loops have been flattened for efficiency.
01221   // + the data is dumped to a buffer (2megs) to make more efficient
01222   //   writes. the size of this buffer can be altered (hardwired). see
01223   //   encodeToBase64.h
01224   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01225   void WriteVTUFile(std::ostream& out)  {
01226     // VTI ASCI File writer. based of VTKASCII writer.
01227     // http://www.cacr.caltech.edu/~slombey/asci/vtk/
01228     int dimshape, conshape;
01229     DimShape(dimshape,conshape);
01230     unsigned long appendedDataOffset=0;
01231 
01232      // hack or needed? not in vtk specs, but seems necessary
01233     bool addHeaderSizeAsInt32=true;
01234 
01235     out << "<?xml version=\"1.0\"?>" << std::endl;
01236     out << "<?meta name=\"source\" content=\"HDF data created with AMROC's hdf2file\"?>" 
01237         << std::endl;
01238     out << "<?meta name=\"timestamp\" content=\"" << RealTime() << "\"?>" << std::endl;
01239 
01240     // VTU XML HEADER
01241     int nnodes = NNodes(), ncells = NCells();
01242     float* xyz = new float[3*nnodes];
01243     out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=";
01244     if (isBigEndian())
01245         out << "\"BigEndian\">" << std::endl;
01246     else
01247         out << "\"LittleEndian\">" << std::endl;
01248     out << "<UnstructuredGrid>" << std::endl;
01249     out << "<Piece NumberOfPoints=\"" << nnodes << "\" NumberOfCells=\"" << ncells << "\">" 
01250         << std::endl;
01251 
01252     // flatten points (i.e. remove duplicates) and print
01253     out << "<Points>" << std::endl;
01254     out << "<DataArray type=";
01255     if (sizeof(float)==4)
01256       out << "\"Float32\"";
01257     else if (sizeof(float)==8)
01258       out << "\"Float64\"";
01259     else out << "\"Float" << sizeof(float)*8 <<"\""; // usupported!
01260     out << " NumberOfComponents=\""<< VTK_USED_DIMS 
01261         << "\" Name=\"vertex\" format=\"appended\" offset=\"" << appendedDataOffset 
01262         << "\" />" << std::endl;
01263 
01264     unsigned long deltaDataOffset= VTK_USED_DIMS*nnodes;
01265     deltaDataOffset*=sizeof(float);
01266     if (addHeaderSizeAsInt32) deltaDataOffset+=4;
01267     if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
01268     deltaDataOffset=deltaDataOffset/3*4;
01269     appendedDataOffset+=deltaDataOffset;
01270 
01271     out << "</Points>" << std::endl;
01272     out << "<Cells>" << std::endl;
01273     out << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" 
01274         << appendedDataOffset << "\" />" << std::endl;
01275 
01276     deltaDataOffset= conshape*ncells;
01277     deltaDataOffset*=sizeof(int);
01278     if (addHeaderSizeAsInt32) deltaDataOffset+=4;
01279     if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
01280     deltaDataOffset=deltaDataOffset/3*4;
01281     appendedDataOffset+=deltaDataOffset;
01282 
01283     out << "<DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" 
01284         << appendedDataOffset << "\" />" << std::endl;
01285     deltaDataOffset= ncells;
01286     deltaDataOffset*=sizeof(int);
01287     if (addHeaderSizeAsInt32) deltaDataOffset+=4;
01288     if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
01289     deltaDataOffset=deltaDataOffset/3*4;
01290     appendedDataOffset+=deltaDataOffset;
01291 
01292     out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\"" 
01293         << appendedDataOffset << "\" />" << std::endl;
01294     deltaDataOffset= ncells;
01295     deltaDataOffset*=sizeof(char);
01296     if (addHeaderSizeAsInt32) deltaDataOffset+=4;
01297     if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
01298     deltaDataOffset=deltaDataOffset/3*4;
01299     appendedDataOffset+=deltaDataOffset;
01300 
01301     out << "</Cells>" << std::endl;
01302     
01303     int length;
01304     if (FKeys() == 1) {
01305       length = nnodes;
01306       out << "<PointData>" << std::endl;
01307     } else {
01308       length = ncells;
01309       out << "<CellData>" << std::endl;
01310     }
01311 
01312     for (int k=0; k<nkeys; k++) {
01313       char* title = new char[LEN_TKEYS];
01314       // clean up title string
01315         std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
01316         int i=std::strlen(title)-1;
01317         while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
01318         int n=std::strlen(title);
01319         for (i=0; i<n; i++) if (title[i]==' ') title[i]='_';
01320 
01321       out << "<DataArray type=";
01322       if (sizeof(float)==4) out << "\"Float32\"";
01323       else if (sizeof(float)==8) out << "\"Float64\"";
01324       out << " Name=\"" << title << "\" NumberOfComponents=\"1\" format=\"appended\" offset=\"" 
01325           << appendedDataOffset << "\" />" << std::endl;
01326 
01327       deltaDataOffset=length;
01328       deltaDataOffset*=sizeof(float);
01329       if (addHeaderSizeAsInt32) deltaDataOffset+=4;
01330       if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
01331       deltaDataOffset=deltaDataOffset/3*4;
01332       appendedDataOffset+=deltaDataOffset;
01333     }
01334 
01335     if (FKeys() == 1)
01336       out << "</PointData>" << std::endl;
01337     else 
01338       out << "</CellData>" << std::endl;
01339     out << "</Piece>" << std::endl;
01340     out << "</UnstructuredGrid>" << std::endl;
01341 
01342     // APPENDED DATA 
01343 
01344     out << "<AppendedData encoding=\"base64\">" << std::endl;
01345     out << "_";
01346 
01347     SetGridNodes((float (*)[3]) xyz);
01348 
01349     // encode x,y,z (go in sets of 3, and encode leftovers at the end)
01350     register int i;
01351     int dataBlockSize;
01352     dataBlockSize=nnodes*VTK_USED_DIMS*sizeof(float);
01353     if (addHeaderSizeAsInt32)
01354       encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
01355     encodeToBase64AndWrite((unsigned char *)xyz,dataBlockSize,out,true);
01356 
01357     // cells
01358     int* con = new int[8*ncells];    
01359 
01360     // print cells->connectivities
01361     SetGridConns((int (*)[8]) con);  
01362 
01363     // encode conns(go in sets of 3, and encode leftovers at the end)
01364     dataBlockSize=ncells*conshape*sizeof(int);
01365     if (addHeaderSizeAsInt32)
01366       encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
01367     int deltamaxcon=8*sizeof(int);
01368     int deltacon   =conshape*sizeof(int);
01369     unsigned char *srccon=(unsigned char *)con;
01370     unsigned char *dstcon=(unsigned char *)con;
01371     if (deltamaxcon!=deltacon && conshape<8) {
01372       for (i=1; i<ncells; i++) {
01373          srccon+=deltamaxcon;
01374          dstcon+=deltacon;
01375          memcpy(dstcon,srccon,deltacon);
01376       }
01377     }
01378     encodeToBase64AndWrite((unsigned char *)con,dataBlockSize,out,true);
01379 
01380     // print cells->element offsets
01381     // encode element offsets (go in sets of 3, and encode leftovers at the end)
01382     dataBlockSize=ncells*sizeof(int);
01383     if (addHeaderSizeAsInt32)
01384       encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
01385     for (i=0; i<ncells; i++) con[i]= (i+1)*conshape;
01386     encodeToBase64AndWrite((unsigned char *)con,dataBlockSize,out,true);
01387 
01388 
01389     // print cells->types
01390     unsigned char cell_type = 1;
01391     if (dimshape==3) cell_type = 11;
01392     else if (dimshape==2) cell_type = 8;
01393     else if (dimshape==1) cell_type = 3;
01394 
01395 
01396     // encode element celltypes (go in sets of 3, and encode leftovers at the end)
01397     dataBlockSize=ncells;
01398     if (addHeaderSizeAsInt32) {
01399       encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
01400     }
01401     int four_cell_types;
01402     ((char *)(&four_cell_types))[0]=cell_type;
01403     ((char *)(&four_cell_types))[1]=cell_type;
01404     ((char *)(&four_cell_types))[2]=cell_type;
01405     ((char *)(&four_cell_types))[3]=cell_type;
01406     for (i=0; i<(ncells+3)/4; i++) con[i]= four_cell_types;
01407     encodeToBase64AndWrite((unsigned char *)con,dataBlockSize,out,true);
01408 
01409     float* dat = new float[nkeys*length];
01410     bzero(dat,nkeys*length*sizeof(float));
01411     register int k;
01412     for (k=0; k<nkeys; k++) 
01413       SetScal(&(keys[k]),&(dat[k*length]));
01414     for (k=0; k<nkeys; k++) {
01415       // encode each of the data keys [vars] (go in sets of 3, and encode leftovers at the end)
01416       dataBlockSize=length*sizeof(float);
01417       if (addHeaderSizeAsInt32) {
01418         encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
01419       }
01420       encodeToBase64AndWrite((unsigned char *)&(dat[k*length]),dataBlockSize,out,true);
01421     }
01422     out <<  std::endl;
01423 
01424     out << "<? meta name=\"datasize\" size=\"" << appendedDataOffset << "\" ?>" << std::endl;
01425 
01426     out << std::endl;
01427     out << "</AppendedData>" << std::endl;
01428     out << "</VTKFile>" << std::endl;
01429  
01430     delete [] xyz;
01431     delete [] con;
01432     delete [] dat;
01433   }
01434 
01435   void WriteVTKASCIIFile(std::ostream& out)  {
01436     int dimshape, conshape;
01437     DimShape(dimshape,conshape);
01438 
01439     int nnodes = NNodes(), ncells = NCells();
01440     float* xyz = new float[3*nnodes]; 
01441     out << "# vtk DataFile Version 3.0" << std::endl;
01442     out << "HDF data created with AMROC's hdf2file" << std::endl;
01443     out << "ASCII" << std::endl;
01444     out << "DATASET UNSTRUCTURED_GRID" << std::endl;
01445     out << "FIELD FieldData 1" << std::endl;
01446     out << "TIME 1 1 float" << std::endl;
01447     out << float(RealTime()) << std::endl; 
01448     out << "POINTS " << nnodes << " float" << std::endl;
01449     SetGridNodes((float (*)[3]) xyz);
01450     register int i;
01451     for (i=0; i<nnodes; i++) {
01452       for (int d=0; d<VTK_USED_DIMS; d++) 
01453         out << xyz[3*i+d] << " ";
01454       out << std::endl;
01455     }
01456 
01457     int* con = new int[8*ncells];    
01458     out << "CELLS " << ncells << " " << ncells*(conshape+1) << std::endl;
01459     SetGridConns((int (*)[8]) con);  
01460     for (i=0; i<ncells; i++) {
01461       out << conshape << " ";
01462       for (int d=0; d<conshape; d++) 
01463         out << con[8*i+d] << " ";
01464       out << std::endl;
01465     }
01466 
01467     int cell_type = 1;
01468     if (dimshape==3) cell_type = 11;
01469     else if (dimshape==2) cell_type = 8;
01470     else if (dimshape==1) cell_type = 3;
01471 
01472     out << "CELL_TYPES " << ncells << std::endl;
01473     for (i=0; i<ncells; i++) 
01474       out << cell_type << std::endl;
01475     
01476     int length;
01477     if (FKeys() == 1) {
01478       length = nnodes;
01479       out << "POINT_DATA ";
01480     }
01481     else {
01482       length = ncells;
01483       out << "CELL_DATA ";
01484     }
01485     out << length << std::endl;    
01486     float* dat = new float[nkeys*length];
01487     bzero(dat,nkeys*length*sizeof(float));
01488 
01489     register int k;
01490     for (k=0; k<nkeys; k++) 
01491       SetScal(&(keys[k]),&(dat[k*length]));
01492     out << "FIELD Data " << nkeys << std::endl;     
01493     for (k=0; k<nkeys; k++) {
01494       char* title = new char[LEN_TKEYS];
01495       std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
01496       i=std::strlen(title)-1;
01497       while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
01498       int n=std::strlen(title);
01499       for (i=0; i<n; i++) 
01500         if (title[i]==' ') title[i]='_';
01501       out << title << " 1 " << length << " float" << std::endl;
01502       for (i=0; i<length; i++) 
01503         out << dat[k*length+i] << std::endl; 
01504     }
01505     
01506     delete [] xyz;
01507     delete [] con;
01508     delete [] dat;
01509   }   
01510 
01511 
01512   void flipBytes(char *c, int n) {
01513      char *ptr=c;
01514      char x;
01515      n/=4;
01516      for (register int i=0; i<n; i++) {
01517        x=ptr[3]; ptr[3]=ptr[0]; ptr[0]=x;
01518        x=ptr[2]; ptr[2]=ptr[1]; ptr[1]=x;
01519        ptr+=4;
01520      }
01521   }
01522   void WriteVTKFile(std::ostream& out) {
01523     int dimshape, conshape;
01524     DimShape(dimshape,conshape);
01525     bool writeSizeAsHeader=false;       // equivalent to addHeaderSizeAsInt32 on binary VTU
01526     bool flipdatabytes=((char *)(&conshape))[3]==0;
01527 
01528     int nnodes = NNodes(), ncells = NCells();
01529     float* xyz = new float[3*nnodes];
01530     out << "# vtk DataFile Version 3.0" << std::endl;
01531     out << "HDF data created with AMROC's hdf2file" << std::endl;
01532     out << "BINARY" << std::endl;
01533     out << "DATASET UNSTRUCTURED_GRID" << std::endl;
01534     out << "FIELD FieldData 1" << std::endl;
01535     out << "TIME 1 1 float" << std::endl;
01536     float timeout = float(RealTime());
01537     if (flipdatabytes) flipBytes((char*) &timeout, sizeof(float));
01538     out.write((char*) &timeout, sizeof(float));
01539     out << std::endl; 
01540     out << "POINTS " << nnodes << " float" << std::endl;
01541     SetGridNodes((float (*)[3]) xyz);
01542     register int n;
01543     if (writeSizeAsHeader) {
01544       int binaryBlockSize= VTK_USED_DIMS*nnodes;
01545       out.write((char*)&binaryBlockSize,sizeof(int));
01546     }
01547     if (flipdatabytes) flipBytes((char*) xyz, VTK_USED_DIMS*nnodes*sizeof(float));
01548     out.write((char*) xyz, VTK_USED_DIMS*nnodes*sizeof(float));
01549     out << std::endl;
01550     delete [] xyz;
01551 
01552     int* con = new int[8*ncells];    
01553     int* con2 = new int[9*ncells];  
01554     int c2shape = conshape+1;
01555     out << "CELLS " << ncells << " " 
01556         << ncells*(conshape+1) << std::endl;
01557     SetGridConns((int (*)[8]) con);  
01558     for (n=0; n<ncells; n++) 
01559       con2[c2shape*n] = conshape;
01560     for (register int i=0; i<conshape; i++)
01561       for (n=0; n<ncells; n++) 
01562         con2[c2shape*n+i+1] = con[8*n+i];
01563     if (writeSizeAsHeader) {
01564       int binaryBlockSize= (conshape+1)*ncells;
01565       out.write((char*)&binaryBlockSize,sizeof(int));
01566     }
01567     if (flipdatabytes) flipBytes((char*) con2, (conshape+1)*ncells*sizeof(int));
01568     out.write((char*) con2, (conshape+1)*ncells*sizeof(int));
01569     out << std::endl;
01570     delete [] con; delete [] con2;
01571     
01572     int cell_type = 1;
01573     if (dimshape==3) cell_type = 11;
01574     else if (dimshape==2) cell_type = 8;
01575     else if (dimshape==1) cell_type = 3;
01576     int* type = new int[ncells];    
01577     for (n=0; n<ncells; n++)
01578       type[n] = cell_type;
01579     out << "CELL_TYPES " << ncells << std::endl;
01580     if (writeSizeAsHeader) {
01581       int binaryBlockSize=  ncells;
01582       out.write((char*)&binaryBlockSize,sizeof(int));
01583     }
01584     if (flipdatabytes) flipBytes((char*) type, ncells*sizeof(int));
01585     out.write((char*) type, ncells*sizeof(int));
01586     out << std::endl;
01587     delete [] type;
01588 
01589     int length;
01590     if (FKeys() == 1) {
01591       length = nnodes;
01592       out << "POINT_DATA ";
01593     }
01594     else {
01595       length = ncells;
01596       out << "CELL_DATA ";
01597     }
01598     out << length << std::endl;    
01599     float* dat = new float[nkeys*length];
01600     for (n=0; n<nkeys*length; n++) *(dat+n) = 0.;
01601     register int k;
01602     for (k=0; k<nkeys; k++) 
01603       SetScal(&(keys[k]),&(dat[k*length]));
01604     out << "FIELD Data " << nkeys << std::endl;  
01605     for (k=0; k<nkeys; k++) {
01606       char* title = new char[LEN_TKEYS];
01607       std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
01608       int i=std::strlen(title)-1;
01609       while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
01610       int n=std::strlen(title);
01611       for (i=0; i<n; i++) 
01612         if (title[i]==' ') title[i]='_';
01613       out << title << " 1 " << length << " float" << std::endl;
01614       if (writeSizeAsHeader) {
01615         int binaryBlockSize=  length;
01616         out.write((char*)&binaryBlockSize,sizeof(int));
01617       }
01618 
01619       if (flipdatabytes) flipBytes((char*) &(dat[k*length]), length*sizeof(float));
01620       out.write((char*) &(dat[k*length]), length*sizeof(float));
01621       out << std::endl;
01622     }
01623     delete [] dat;
01624   }    
01625 
01626   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01627   //                                                   WRITE SILO FILE
01628   //
01629   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01630 #if defined(HAVE_LIBSILO) || defined(HAVE_LIBSILOH5)
01631   void WriteSiloFile(char *silofilename)  {
01632     if (silofilename==0 || silofilename[0]==0) {
01633       std::cerr << "Abort: Cannot write SILO to stdout. "
01634                 << "Please provide FileName in your display_file.in" << std::endl;
01635       return;
01636     }
01637 #if defined(HAVE_LIBSILOH5)
01638     DBfile * silofptr= DBCreate(silofilename, DB_CLOBBER, DB_LOCAL, NULL, DB_HDF5);
01639 #else
01640     DBfile * silofptr= DBCreate(silofilename, DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);
01641 #endif
01642     if (!silofptr) {
01643       std::cerr << "Abort: could not create SILO file." << std::endl;
01644       return;
01645     }
01646 
01647     int dimshape, conshape;
01648     DimShape(dimshape,conshape);
01649 
01650     int nnodes = NNodes(), ncells = NCells();
01651     int vardatacentering, length;
01652     if (FKeys() == 1) {
01653       vardatacentering=DB_NODECENT;
01654       length = nnodes;
01655     } else {
01656       vardatacentering=DB_ZONECENT;
01657       length = ncells;
01658     }
01659 
01660     std::cerr << std::endl;
01661 
01662     float * xyz = new float[3*nnodes];
01663     char  * coordnames [3]= {"X","Y","Z"};
01664     float * coordinates[3];
01665 
01666     coordinates[0]=xyz;
01667     coordinates[1]=coordinates[0]+nnodes;
01668     coordinates[2]=coordinates[1]+nnodes;
01669 
01670 #   ifdef USE_PRESET_AMROC_BLOCK_COORDS
01671     SetGridNodes((float (*)[3]) xyz);
01672     register int n,t;
01673     if (dimshape==1) {
01674       for (n=0,t=0; n<nnodes; n++)
01675         xyz[t++] = xyz[3*n];
01676     }
01677     else if (dimshape==2) {
01678       for (n=0,t=0; n<nnodes; n++) {
01679         xyz[t++] = xyz[3*n]; xyz[t++] = xyz[3*n+1];
01680       }
01681     }
01682 #   endif
01683 
01684     float* dat = new float[nkeys*length];
01685     register int i,k;
01686     for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
01687     for (k=0; k<nkeys; k++) 
01688       SetScal(&(keys[k]),&(dat[k*length]));
01689 
01690     int    *blocklevel = new int[BlockList().size()];
01691     char  **vartitles  = new char* [nkeys];
01692     int     leveltitle = -1;
01693 
01694     // get variable name
01695     for (k=0; k<nkeys; k++) {
01696       char *title = new char[LEN_TKEYS];
01697       std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
01698       i=std::strlen(title)-1;
01699       while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
01700       int n=std::strlen(title);
01701       for (i=0; i<n; i++) {
01702         if (title[i]==' ' ||
01703             title[i]=='.' ||
01704             title[i]=='-' ||
01705             title[i]=='+' ||
01706             title[i]=='/'    ) title[i]='_';
01707       }
01708       vartitles[k] = new char[std::strlen(title)+1];
01709       for (i=0; i<n; i++) 
01710         vartitles[k][i]=std::toupper(title[i]);
01711       vartitles[k][i]='\0';
01712       delete [] title;
01713       if (!std::strcmp(vartitles[k],"LEVEL") || !std::strcmp(vartitles[k],"LEVELS"))
01714       leveltitle=k;
01715     }
01716 
01717     int offset=0, maxlevel=-1, nblock=0;
01718     char **meshptrs= new char *[BlockList().size()];
01719     char **varptrs=  new char *[BlockList().size()];
01720 
01721     for (typename block_list_type::iterator bit=BlockList().begin();
01722          bit!=BlockList().end(); bit++, nblock++) {
01723       float fmin[3], fmax[3];
01724       int idims[3], vardims[3];
01725       for (i=0; i<dimshape; i++) {
01726         fmin[i] = (**bit).DB().bbox().lower(i)*base::dx[i]+base::geom[i*2];
01727         fmax[i] = ((**bit).DB().bbox().upper(i)+
01728                    (**bit).DB().bbox().stepsize(i))*base::dx[i]+base::geom[i*2];
01729         idims[i]= (**bit).DB().bbox().extents(i)+1;
01730         vardims[i]=idims[i];
01731         if (vardatacentering==DB_ZONECENT) vardims[i]--;
01732         if (vardims[i]<=0) vardims[i]=1;
01733       }
01734 
01735       // creating a directory per block
01736       char *meshblockdirname= new char[LEN_TKEYS];
01737       meshblockdirname[0]='/';
01738       meshblockdirname++;
01739       sprintf(meshblockdirname,"BLOCK%04i",nblock);
01740       DBSetDir(silofptr,"/");
01741       DBMkDir(silofptr,meshblockdirname);
01742       DBSetDir(silofptr,meshblockdirname);
01743 
01744       char *meshblockname= meshblockdirname + std::strlen(meshblockdirname)+1;
01745       meshblockname[-1]='/';
01746       sprintf(meshblockname,"block%04i",nblock);
01747 
01748       meshptrs[nblock]=meshblockdirname;
01749       varptrs[nblock]= meshblockname;
01750 
01751 #     ifdef USE_PRESET_AMROC_BLOCK_COORDS
01752       std::cerr<<"using USE_PRESET_AMROC_BLOCK_COORDS" << std::endl;
01753       for (i=0; i<dimshape; i++)
01754         coordinates[i]=&(xyz[dimshape*offset+i]);
01755       DBPutQuadmesh(silofptr, meshblockname ,coordnames, coordinates, idims, dimshape, DB_FLOAT, DB_COLLINEAR, NULL);
01756 #     else
01757       for (i=0; i<dimshape; i++) {
01758         for (register int idx=0; idx<idims[i]; idx++)
01759           coordinates[i][idx]=fmin[i]+(float)idx*(fmax[i]-fmin[i])/(float)(idims[i]-1);
01760       }
01761       DBPutQuadmesh(silofptr, meshblockname ,coordnames, coordinates, idims, dimshape, DB_FLOAT, DB_COLLINEAR, NULL);
01762 #     endif
01763 
01764       for (k=0; k<nkeys; k++) {
01765         if (k==leveltitle) {
01766           blocklevel[nblock]=(**bit).Level();
01767           if (blocklevel[nblock]>maxlevel) maxlevel=blocklevel[nblock];
01768         }
01769         DBPutQuadvar1(silofptr, vartitles[k], meshblockname, &(dat[offset+k*length]),
01770                       vardims, dimshape, NULL, 0,DB_FLOAT,  vardatacentering, NULL);
01771       }
01772       if (FKeys()==1) offset += (**bit).NNodes();
01773       else offset += (**bit).NCells();
01774     }
01775 
01776     bool   dump_sets=true;
01777     if (dump_sets) {
01778     DBSetDir(silofptr,"/");
01779     char * levelmesh = new char[LEN_TKEYS];
01780     char * varmeshname = new char[LEN_TKEYS];
01781  
01782     int   * meshtypes=new int[BlockList().size()];
01783     char ** imeshptrs=new char *[nblock];
01784     for (int i=0; i<nblock; i++) {
01785       meshtypes[i]= DB_QUAD_RECT;
01786       imeshptrs[i]=meshptrs[nblock-i-1];
01787     }
01788     DBPutMultimesh(silofptr, "/ALL_LEVELS", nblock, meshptrs, meshtypes, NULL);
01789 
01790     DBSetDir(silofptr,"/");
01791     DBMkDir(silofptr,"LEVELS");
01792     for (int l=maxlevel; l>=0; l--) {
01793       int ninlevel=0;
01794       for (register int idx=0; idx<nblock; idx++) {
01795         if (blocklevel[idx]==l) {
01796           imeshptrs[ninlevel++]=meshptrs[idx]-1;
01797         }
01798       }  
01799       if (ninlevel>0) {
01800         sprintf(levelmesh,"/LEVELS/LEVEL_%03i",l);
01801         DBPutMultimesh(silofptr, levelmesh, ninlevel, imeshptrs, meshtypes, NULL);
01802       }
01803     }
01804 
01805     for (register int idx=0; idx<nblock; idx++) meshtypes[idx]= DB_QUADVAR;
01806     for (k=0; k<nkeys; k++) {
01807       sprintf(varmeshname,"VAR_%s",vartitles[k]);
01808       for (register int idx=0; idx<nblock; idx++) std::strcpy(varptrs[idx],vartitles[k]);
01809 
01810       DBSetDir(silofptr,"/");
01811       DBMkDir(silofptr,varmeshname);
01812       DBSetDir(silofptr,varmeshname);
01813       sprintf(varmeshname,"/VAR_%s/ALL_LEVELS",vartitles[k]);
01814       for (register int idx=0; idx<nblock; idx++) imeshptrs[idx]=meshptrs[idx]-1;
01815       DBPutMultivar(silofptr, varmeshname, nblock, imeshptrs, meshtypes, NULL);
01816       DBSetDir(silofptr,"/");
01817 
01818       for (int l=maxlevel; l>=0; l--) {
01819         int ninlevel=0;
01820         for (register int idx=0; idx<nblock; idx++)
01821           if (blocklevel[idx]==l)
01822             imeshptrs[ninlevel++]=meshptrs[idx]-1;
01823         if (ninlevel>0) {
01824           sprintf(varmeshname,"/VAR_%s/LEVEL_%03i",vartitles[k],l);
01825           DBPutMultivar(silofptr, varmeshname, ninlevel, imeshptrs, meshtypes, NULL);
01826         }
01827       }
01828     }
01829     delete [] levelmesh;
01830     delete [] varmeshname;
01831     }
01832     DBClose(silofptr);
01833 
01834 
01835     //for (i=0; i<nblock; i++) if (meshptrs[i]) delete meshptrs[i];
01836     delete [] meshptrs;
01837     delete [] varptrs;
01838     delete [] xyz;
01839     delete [] dat;
01840     delete [] blocklevel;
01841     delete [] vartitles;
01842 
01843   }
01844 #else
01845   void WriteSiloFile(char *silofilename)  {
01846      std::cerr << "hdf2file SILO package not compiled. "
01847                << "Please reconfigure with SILO_DIR defined and recompile." << std::endl;
01848   }  
01849 #endif
01850 
01851   void WriteTecplotASCIIMeshFile(std::ostream& out) {
01852     if (!keys) return;
01853 
01854     std::ifstream in(InputFile.c_str(), std::ios::in);
01855     if (!in) {
01856       std::cerr << "Cannot open " << InputFile << " for reading!" << std::endl;
01857       return;
01858     }
01859 
01860     std::string line;
01861     int d, k;
01862     do {
01863       std::getline(in, line);
01864       if (in.eof()) continue;
01865       out << line;
01866       std::transform(line.begin(), line.end(), line.begin(), toupper);
01867       if (line.find("VARIABLES")!=std::string::npos) break;
01868       out << std::endl;
01869     } while (!in.eof());
01870     
01871     register int i;
01872     for (k=0; k<nkeys; k++) {
01873       char* title = new char[LEN_TKEYS];
01874       std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
01875       i=std::strlen(title)-1;
01876       while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }      
01877       int n=std::strlen(title);
01878       for (i=0; i<n; i++) if (title[i]==' ') title[i]='_';
01879       out << " " << title;
01880     }
01881     out << std::endl;
01882     
01883     unsigned int num_vertices=0;
01884     do {
01885       std::getline(in, line);
01886       if (in.eof()) continue;
01887       out << line << std::endl;      
01888       std::transform(line.begin(), line.end(), line.begin(), toupper);
01889       if (line.find("ZONE")!=std::string::npos) {
01890         std::string::size_type idx = line.find("I=");
01891         std::string strnumber;
01892         strnumber.assign(line, idx+2, line.size());
01893         std::istringstream numberstr(strnumber);
01894         numberstr >> num_vertices;
01895         break;
01896       }
01897     } while (!in.eof());
01898 
01899     float SP[3];
01900     point_list_type p;
01901     std::vector<std::string> ds;
01902     do {
01903       std::getline(in, line);
01904       if (in.eof()) continue;
01905       std::istringstream linestr(line);
01906       ds.insert(ds.end(), line);
01907       for (d=0; d<3; d++) 
01908         if (base::DimUsed(d))
01909           linestr >> SP[d];
01910         else
01911           SP[d]=0.0;
01912       p.insert(p.end(), SP, SP+3);
01913     } while (!in.eof() && p.size()<3*num_vertices);
01914 
01915     if (!p.empty()) {
01916       unsigned int Np=p.size()/3;
01917       float *v=new float[Np*nkeys];
01918       for (k=0; k<nkeys; k++) 
01919         if (FKeys() == 1)
01920           base::SetPointsNodes(&(keys[k]),p,&(v[k*Np]));
01921         else
01922           base::SetPointsCells(&(keys[k]),p,&(v[k*Np]));
01923       
01924       for (register unsigned int i=0;i<Np;i++) {
01925         out << ds[i];
01926         for (k=0; k<nkeys; k++) 
01927           out << " " << v[k*Np+i];
01928         out << std::endl;
01929       }
01930       delete [] v;
01931     } 
01932 
01933     do {
01934       std::getline(in, line);
01935       if (in.eof()) continue;
01936       out << line << std::endl;
01937     } while (!in.eof());   
01938     in.close();
01939   }
01940 
01941   virtual void SetGrid(float xyz[][3]) {
01942     if (FKeys() == 1)
01943       SetGridNodes(xyz);
01944     else
01945       SetGridCells(xyz);
01946   }
01947 
01948   virtual void SetScal(int *key,float v[]) {
01949     if (FKeys() == 1)
01950       SetScalNodes(key, v);
01951     else
01952       SetScalCells(key, v);
01953   }
01954 
01955   virtual void Vect(int *key,float v[][3]) {
01956     if (FKeys() == 1)
01957       VectNodes(key, v);
01958     else
01959       VectCells(key, v);
01960   }
01961 
01962   virtual void BlockListAppend(hdata_block_type& workdata, BBox& bb, int comp, int& Level,
01963                                int& proc, float* dx, float* geom) {
01964     BlockList().push_back(new data_block_type(*this, workdata, bb, comp,
01965                                               TheEvaluator(), Level, proc, dx, geom));
01966   }
01967 
01968   block_list_type& BlockList() { return (block_list_type &)*base::_BlockList; }
01969   block_list_type& BlockList() const { return (block_list_type &)*base::_BlockList; }
01970 
01971   int LineProbe() const 
01972     { return (LinePoint[0]!=LineVector[0] || LinePoint[1]!=LineVector[1] || 
01973               LinePoint[2]!=LineVector[2]); }
01974  
01975   virtual void register_at(ControlDevice& Ctrl,const std::string& prefix) {
01976     base::register_at(Ctrl,prefix); }
01977   virtual void register_at(ControlDevice& Ctrl) {
01978     register_at(Ctrl, ""); }
01979   virtual int NCells() { return base::NCells(); }
01980   virtual int NNodes() { return base::NNodes(); }
01981   virtual void SetBlocks(int blocks[]) { base::SetBlocks(blocks); }
01982   virtual void SetGridCells(float xyz[][3]) {
01983     base::SetGridCells(xyz); }
01984   virtual void SetGridNodes(float xyz[][3]) {
01985     base::SetGridNodes(xyz); }
01986   virtual void SetGridConns(int con[][8]) {
01987     base::SetGridConns(con); }
01988   virtual void SetScalCells(int *key,float v[]) {
01989     base::SetScalCells(key,v); }
01990   virtual void SetScalNodes(int *key,float v[]) {
01991     base::SetScalNodes(key,v); }
01992   virtual void VectNodes(int *key,float v[][3]) {
01993     base::VectNodes(key,v); }
01994   virtual void VectCells(int *key,float v[][3]) {
01995     base::VectCells(key,v); }
01996   virtual int Size() const { return base::Size(); }
01997   virtual const int& MinLevel() const { return base::MinLevel(); }
01998   virtual const int& MaxLevel() const { return base::MaxLevel(); }
01999   virtual const int& NLevels() const { return base::NLevels(); }
02000   virtual const int& RefineFactor(const int l) const {
02001     return base::RefineFactor(l); }
02002   virtual const BBox& GlobalBBox() const {
02003     return base::GlobalBBox(); }
02004   virtual const int& FKeys() const { return base::FKeys(); }
02005   virtual evaluator_type& TheEvaluator() const {
02006     return base::TheEvaluator(); }
02007   virtual evaluator_type& TheEvaluator() {
02008     return base::TheEvaluator(); }
02009   virtual float64 RealTime() const { return base::RealTime(); }
02010   virtual void SetUpdateName(std::string* iname) {
02011     base::SetUpdateName(iname); }
02012   virtual std::string* UpdateName() { return base::UpdateName(); }
02013   virtual std::string* UpdateName() const {
02014     return base::UpdateName(); }
02015 
02016 private:
02017   std::string OutputName;
02018   int OutputType, OutputStep;   
02019   std::string StrKeys;   
02020   int* keys;
02021   int nkeys;
02022   float LinePoint[3], LineVector[3];
02023   int RegData[3];
02024   float RegLim[2];
02025   int RegMax;
02026   int LevelCutOut;
02027   int PointOutput;
02028   int _Step;
02029   int FileMerge;
02030   std::string InputFile;
02031 };
02032 
02033 #endif
Generated on Thu Jun 30 02:19:02 2016 for AMROC's HDF Tools by 1.6.3