00001
00002
00003
00004
00005
00006 # ifndef strequals
00007 # define strequals(src,dst) (!strcmp(src,dst))
00008 # endif
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 void WriteSiloFile(char *silofilename) {
00023 # if DEBUG
00024 # ifdef OBNOXIOUS_DEBUG
00025 bool verbose=true;
00026 # else
00027 bool verbose=false;
00028 # endif
00029 # else
00030 bool verbose=false;
00031 # endif
00032
00033 using std::malloc;
00034 using std::toupper;
00035
00036 bool matchCoordAndVarDims=true;
00037 int ndims=3;
00038 int nvardims;
00039 int nodespercell;
00040
00041 DimShape(nvardims,nodespercell);
00042
00043 if (matchCoordAndVarDims) ndims=nvardims;
00044
00045 if (silofilename==NULL || silofilename[0]==0) {
00046 cerr << "Abort: currentl cannot write SILO to stdout. "
00047 "Please provide FileName in your display_file.in" << endl;
00048 return;
00049 }
00050
00051 if (verbose) {
00052 cerr << endl << LEN_TKEYS << endl;
00053 cerr << "Starting convesion to Silo." << endl;
00054 }
00055
00056
00057 SiloCreator * silofptr= new SiloCreator(silofilename);
00058
00059 if (silofptr==NULL) {
00060 cerr << "Abort: could not create SILO file." << endl;
00061 return;
00062 }
00063
00064 if (verbose) {
00065 cerr << "Volume bounds " << geom[0] << ","
00066 << geom[2] << ","
00067 << geom[4] << " - "
00068 << geom[1] << ","
00069 << geom[3] << ","
00070 << geom[5] << " :: "
00071 << dx[0] << ","
00072 << dx[1] << ","
00073 << dx[2] << endl;
00074 }
00075
00076 char * meshblockname;
00077 char * meshblockdirname;
00078 char * title;
00079 double min[3],max[3];
00080 int dims[3]={0,0,0};
00081 float fmin[3];
00082 float fmax[3];
00083 double spacing[3];
00084 int idims[3];
00085 int vardims[3];
00086 int nchunk=0;
00087 int nblock=0;
00088 int maxlevel=-1;
00089
00090 float **blockptrs= (float **)malloc(nkeys*BlockList().size()*sizeof(float *));
00091 int *blocklevel=(int *)malloc(BlockList().size()*sizeof(int));
00092 char **vartitles= (char **)malloc(sizeof(char *)*nkeys);
00093
00094 SiloMesh **meshptrs= (SiloMesh **)malloc(BlockList().size()*sizeof(SiloMesh *));
00095 SiloVariable **varptrs= (SiloVariable **)malloc(BlockList().size()*sizeof(SiloVariable *)*nkeys);
00096 int vardatacentering=DB_NODECENT;
00097
00098 for (typename block_list_type::iterator bit=BlockList().begin();
00099 bit!=BlockList().end();
00100 bit++, nblock++) {
00101 int length;
00102
00103 if (FKeys() == 1) {
00104 vardatacentering=DB_NODECENT;
00105 length = (**bit).NNodes();
00106 } else {
00107 vardatacentering=DB_ZONECENT;
00108 length = (**bit).NCells();
00109 }
00110
00111 float* dat = new float[length*nkeys];
00112
00113 if (dat==NULL) {
00114 cerr<<"Abort: could not allocate memory for block no." << nblock << endl;
00115 return;
00116 }
00117
00118 if (verbose) {
00119 cerr << "Encoding block #" << nblock << " of size " << length << endl;
00120 }
00121
00122 for (int i=0; i<ndims; i++) {
00123
00124 min[i]=(**bit).DB().bbox().lower(i);
00125 max[i]=(**bit).DB().bbox().upper(i);
00126
00127
00128 dims[i]=(**bit).DB().bbox().extents(i);
00129 dims[i]++;
00130
00131
00132 spacing[i]=(**bit).DB().bbox().stepsize(i);
00133 max[i]+=spacing[i];
00134
00135
00136 if (verbose) {
00137 if (max[i]>48) {
00138 cerr << "for axis " << i << " got max " << max[i] << endl;
00139 }
00140 }
00141
00142 fmin[i]=(float)(min[i]*(double)dx[i]+(double)geom[i*2]);
00143 fmax[i]=(float)(max[i]*(double)dx[i]+(double)geom[i*2]);
00144 idims[i]=(int)dims[i];
00145 vardims[i]=idims[i];
00146 if (vardatacentering==DB_ZONECENT) vardims[i]--;
00147 if (vardims[i]<=0) vardims[i]=1;
00148
00149
00150 if (verbose) {
00151 if (fmax[i]>geom[i*2+1]) {
00152 cerr << "for axis " << i << " got fmax " << fmax[i] << endl;
00153 }
00154 }
00155 }
00156
00157 if (verbose) {
00158 cerr << " [" << min[0] << "," <<
00159 min[1] << "," <<
00160 min[2] << "] - [" <<
00161 max[0] << "," <<
00162 max[1] << "," <<
00163 max[2] << "] steps of (" <<
00164 spacing[0] << "," <<
00165 spacing[1] << "," <<
00166 spacing[2] << ")" << endl;
00167 }
00168
00169
00170 meshblockdirname= new char[LEN_TKEYS];
00171
00172 sprintf(meshblockdirname,"BLOCK%04i",nblock);
00173 silofptr->SetToRoot();
00174 silofptr->MakeAndSetDirectory(meshblockdirname);
00175
00176 meshblockname= new char[LEN_TKEYS];
00177 sprintf(meshblockname,"block%04i",nblock);
00178
00179 if (verbose) {
00180 cerr << "* ";
00181 cerr << meshblockname << ":" << idims[0] << "x" << idims[1] << "x" << idims[2] << endl;
00182 cerr << " ";
00183 cerr << fmin[0] << "," << fmin[1] << "," << fmin[2] << " - ";
00184 cerr << fmax[0] << "," << fmax[1] << "," << fmax[2] << endl;
00185 }
00186
00187 # ifdef USE_PRESET_AMROC_BLOCK_COORDS
00188
00189
00190 float* xyz = new float[3*(**bit).DB().bbox().size()];
00191 SetGridNodes((float (*)[3]) xyz);
00192 if (ndims<3)
00193 for (int idx=0,cidx=0; cidx<NNodes();cidx++)
00194 for (int icidx=0; icidx<ndims; icidx++)
00195 xyz[idx++]=xyz[cidx*3+icidx];
00196 meshptrs[nblock]= new SiloMesh((char *)meshblockname,ndims,(int *)idims,(float *)xyz);
00197 # else
00198 meshptrs[nblock]= new SiloMesh((char *)meshblockname,ndims,(int *)idims,(float *)fmin,(float *)fmax);
00199 # endif
00200 silofptr->SetMesh(meshptrs[nblock]);
00201 silofptr->DrawMesh();
00202
00203 int offset=0;
00204 for (register int k=0; k<nkeys; k++) {
00205
00206 blockptrs[nchunk]=dat+offset;
00207
00208 int poffset=offset;
00209 if (FKeys() == 1) {
00210
00211 (**bit).SetScalNodes(&(keys[k]),dat,offset,base::CompRead);
00212 } else {
00213
00214 (**bit).SetScalCells(&(keys[k]),dat,offset,base::CompRead);
00215 }
00216 if (offset==poffset) {
00217 std::cerr << "Warning: SetScal::CompRead of key #" << k << " created no data offset. Variable may not be defined." << std::endl;
00218 }
00219
00220
00221
00222 title = new char[LEN_TKEYS];
00223 strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
00224 int i=strlen(title)-1;
00225 while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
00226 int n=strlen(title);
00227 for (i=0; i<n; i++) {
00228 if (title[i]==' ') title[i]='_';
00229 if (title[i]=='-') title[i]='_';
00230 if (title[i]=='+') title[i]='_';
00231 if (title[i]=='/') title[i]='_';
00232 }
00233
00234
00235
00236 char* varTITLE = (char*) malloc(strlen(title)+1);
00237 strcpy(varTITLE,title);
00238 for (unsigned int i=0; i<strlen(title); i++) varTITLE[i]=toupper(title[i]);
00239
00240 if (nblock==0) {
00241
00242
00243 vartitles[nchunk] = (char*) malloc(strlen(varTITLE)+1);
00244 strcpy(vartitles[nchunk],varTITLE);
00245 }
00246
00247 if (strequals(varTITLE,"LEVEL") || strequals(varTITLE,"LEVELS")) {
00248 blocklevel[nblock]=(**bit).Level();
00249 if (blocklevel[nblock]>maxlevel) maxlevel=blocklevel[nblock];
00250 }
00251 # ifdef DEBUG_REPLACE_VALUES
00252 if (nvardims==3) {
00253 int idx=0;
00254 for (int idx_k=0; idx_k<idims[2]; idx_k++)
00255 for (int idx_j=0; idx_j<idims[1]; idx_j++)
00256 for (int idx_i=0; idx_i<idims[0]; idx_i++,idx++)
00257 blockptrs[nchunk][idx]=(((k%2)==1)?k:-k)*(min[0]+min[1]+min[2]+idx_i+idx_j+idx_k);
00258
00259 }
00260 # endif
00261 varptrs[nchunk]=new SiloVariable(varTITLE,nvardims,vardims,meshptrs[nblock],1,blockptrs+nchunk,vardatacentering);
00262 silofptr->SetVariable(varptrs[nchunk]);
00263 silofptr->DrawVariable();
00264 if (verbose) cerr << " +" << varTITLE << endl;
00265 nchunk++;
00266 }
00267
00268 }
00269
00270 # ifndef SILO_NO_GROUPING
00271 silofptr->SetToRoot();
00272 silofptr->CreateMultiMesh("/ALL_LEVELS");
00273 SiloMesh **imeshptrs= (SiloMesh **)malloc(BlockList().size()*sizeof(SiloMesh *));
00274 for (int i=0; i<nblock; i++) imeshptrs[i]=meshptrs[nblock-i-1];
00275 silofptr->SetMultiMesh(nblock,imeshptrs);
00276 silofptr->WriteMultiMesh();
00277
00278 silofptr->SetToRoot();
00279 silofptr->MakeAndSetDirectory("LEVELS");
00280
00281 if (verbose) cerr << "Creating per level collections." << endl;
00282 for (int i=maxlevel; i>=0; i--) {
00283 int ninlevel=0;
00284 for (int idx=0; idx<nblock; idx++) if (blocklevel[idx]==i) ninlevel++;
00285 if (ninlevel>0) {
00286 if (verbose) cerr << "Creating mesh collection for level " << i << " [" << ninlevel << "] ";
00287 SiloMesh ** blocksinlevel= (SiloMesh **)malloc(sizeof(SiloMesh *)*ninlevel);
00288 for (int idx=0, ni=0; idx<nblock; idx++)
00289 if (blocklevel[idx]==i)
00290 blocksinlevel[ni++]=meshptrs[idx];
00291 char * levelmesh = new char[LEN_TKEYS];
00292 sprintf(levelmesh,"/LEVELS/LEVEL_%03i",i);
00293 if (verbose) cerr << levelmesh << endl;
00294 silofptr->CreateMultiMesh(levelmesh);
00295 silofptr->SetMultiMesh(ninlevel,blocksinlevel);
00296 silofptr->WriteMultiMesh();
00297 }
00298 }
00299
00300 int nvars= nkeys;
00301 for (int i=0; i<nvars; i++) {
00302 SiloVariable ** blocksinvar= (SiloVariable **)malloc(sizeof(SiloVariable *)*nblock);
00303
00304 for (int ni=nblock, idx=i; idx<nchunk; idx+=nvars) blocksinvar[--ni]=varptrs[idx];
00305
00306 char * varmeshname = new char[LEN_TKEYS];
00307 silofptr->SetToRoot();
00308 sprintf(varmeshname,"VAR_%s",vartitles[i]);
00309 silofptr->MakeAndSetDirectory(varmeshname);
00310
00311 sprintf(varmeshname,"/VAR_%s/ALL_LEVELS",vartitles[i]);
00312
00313 silofptr->CreateMultiVar(varmeshname);
00314 silofptr->SetMultiVar(nblock,blocksinvar);
00315 silofptr->WriteMultiVar();
00316
00317 for (int l=maxlevel; l>=0; l--) {
00318 int ninlevel=0;
00319 for (int idx=0; idx<nblock; idx++)
00320 if (blocklevel[idx]==l) ninlevel++;
00321 if (ninlevel>0) {
00322 SiloVariable ** blocksinlevel= (SiloVariable **)malloc(sizeof(SiloVariable *)*ninlevel);
00323 for (int idx=0, ni=0; idx<nblock; idx++)
00324 if (blocklevel[idx]==l)
00325 blocksinlevel[ni++]=varptrs[idx*nvars+i];
00326 sprintf(varmeshname,"/VAR_%s/LEVEL_%03i",vartitles[i],l);
00327 if (verbose)
00328 cerr << "Writing collection " << varmeshname << " " << ninlevel << "/" << nblock << endl;
00329 silofptr->CreateMultiVar(varmeshname);
00330 silofptr->SetMultiVar(ninlevel,blocksinlevel);
00331 silofptr->WriteMultiVar();
00332 }
00333 }
00334 }
00335 # endif
00336
00337
00338
00339 if (verbose) cerr << "Closing silo file " << silofilename << endl ;
00340 silofptr->Close();
00341
00342 if (verbose) cerr << "DONE" << endl;
00343
00344
00345
00346 cerr << endl;
00347
00348 }