/[escript]/trunk/weipa/src/EscriptDataset.cpp
ViewVC logotype

Annotation of /trunk/weipa/src/EscriptDataset.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6144 - (hide annotations)
Wed Apr 6 05:25:13 2016 UTC (2 years, 8 months ago) by caltinay
File size: 37320 byte(s)
last round of namespacing defines.

1 caltinay 2183
2 jfenwick 3981 /*****************************************************************************
3 caltinay 2183 *
4 jfenwick 5863 * Copyright (c) 2003-2016 by The University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 caltinay 2183 *
7     * Primary Business: Queensland, Australia
8 jfenwick 6112 * Licensed under the Apache License, version 2.0
9     * http://www.apache.org/licenses/LICENSE-2.0
10 caltinay 2183 *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
13     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 *
15     *****************************************************************************/
16 caltinay 2183
17 caltinay 3037 #include <weipa/EscriptDataset.h>
18     #include <weipa/DataVar.h>
19     #include <weipa/ElementData.h>
20 caltinay 5184 #include <weipa/NodeData.h>
21     #if defined USE_FINLEY || defined USE_DUDLEY
22 caltinay 3143 #include <weipa/FinleyDomain.h>
23 caltinay 5184 #endif
24     #ifdef USE_RIPLEY
25 caltinay 3792 #include <weipa/RipleyDomain.h>
26 caltinay 5184 #endif
27     #ifdef USE_SPECKLEY
28 sshaw 5123 #include <weipa/SpeckleyDomain.h>
29 caltinay 5184 #endif
30 caltinay 2910
31     #ifndef VISIT_PLUGIN
32 caltinay 2810 #include <escript/Data.h>
33 caltinay 5968 #include <escript/FileWriter.h>
34 caltinay 5184 #ifdef USE_DUDLEY
35 jfenwick 3259 #include <dudley/CppAdapter/MeshAdapter.h>
36 caltinay 5184 #endif
37     #ifdef USE_FINLEY
38 caltinay 3143 #include <finley/CppAdapter/MeshAdapter.h>
39 caltinay 5184 #endif
40     #ifdef USE_RIPLEY
41 caltinay 3792 #include <ripley/RipleyDomain.h>
42 caltinay 5184 #endif
43     #ifdef USE_SPECKLEY
44 sshaw 5123 #include <speckley/SpeckleyDomain.h>
45 caltinay 5184 #endif
46 caltinay 4320
47 caltinay 5968 using escript::FileWriter;
48 caltinay 2910 #endif
49 caltinay 2810
50 caltinay 3355 #include <cstring>
51 jfenwick 3259 #include <iostream>
52 caltinay 2886 #include <numeric> // for std::accumulate
53 caltinay 3191 #include <sstream> // for std::ostringstream
54 caltinay 2810
55 caltinay 6144 #if ESYS_HAVE_SILO
56 caltinay 2183 #include <silo.h>
57 caltinay 2810
58 caltinay 6119 #if WEIPA_HAVE_MPI
59 caltinay 2810 #include <pmpio.h>
60 caltinay 4920 const int NUM_SILO_FILES = 1; // number of Silo files to produce per save
61 caltinay 2183 #endif
62    
63 caltinay 6144 #endif // ESYS_HAVE_SILO
64 caltinay 2810
65 caltinay 2183 using namespace std;
66    
67 caltinay 3037 namespace weipa {
68 caltinay 2187
69 caltinay 2190 const char* MESH_VARS = "mesh_vars/";
70 jfenwick 4911
71 caltinay 2183 //
72 caltinay 2821 // Default constructor
73 caltinay 2183 //
74 caltinay 2810 EscriptDataset::EscriptDataset() :
75     cycle(0),
76     time(0.),
77 caltinay 3191 externalDomain(false),
78 caltinay 3623 wantsMeshVars(false),
79 caltinay 2810 mpiRank(0),
80     mpiSize(1)
81 caltinay 2183 {
82 caltinay 2821 }
83    
84     //
85     // Constructor with communicator
86     //
87 caltinay 6119 #if WEIPA_HAVE_MPI
88 caltinay 2821 EscriptDataset::EscriptDataset(MPI_Comm comm) :
89     cycle(0),
90     time(0.),
91 caltinay 3191 externalDomain(false),
92 caltinay 3623 wantsMeshVars(false),
93 caltinay 2821 mpiComm(comm)
94     {
95 caltinay 2810 MPI_Comm_rank(mpiComm, &mpiRank);
96     MPI_Comm_size(mpiComm, &mpiSize);
97 caltinay 2821 }
98 caltinay 2810 #endif
99 caltinay 2183
100     //
101     // Destructor
102     //
103 caltinay 2810 EscriptDataset::~EscriptDataset()
104 caltinay 2183 {
105 caltinay 2810 }
106 caltinay 2183
107 caltinay 2810 //
108 caltinay 3128 // Sets the domain using an escript domain instance.
109 caltinay 2810 //
110 caltinay 3128 bool EscriptDataset::setDomain(const escript::AbstractDomain* domain)
111 caltinay 2810 {
112 caltinay 3095 #ifndef VISIT_PLUGIN
113 caltinay 3128 int myError = 0, gError;
114    
115     // fail if the domain has already been set
116 caltinay 3191 if (domainChunks.size() > 0) {
117 caltinay 3128 cerr << "Domain has already been set!" << endl;
118     myError = 1;
119     } else if (!domain) {
120     cerr << "Domain is NULL!" << endl;
121     myError = 1;
122     } else {
123 caltinay 6119 #if WEIPA_HAVE_MPI
124 caltinay 3128 mpiComm = domain->getMPIComm();
125     mpiRank = domain->getMPIRank();
126     mpiSize = domain->getMPISize();
127     #endif
128 caltinay 5184
129     // this is here to allow using 'else if' for all selectively compiled
130     // paths
131     if (0) {
132     }
133     #if USE_FINLEY
134     else if (dynamic_cast<const finley::MeshAdapter*>(domain)) {
135 caltinay 3143 DomainChunk_ptr dom(new FinleyDomain());
136     if (dom->initFromEscript(domain)) {
137     if (mpiSize > 1)
138     dom->reorderGhostZones(mpiRank);
139 caltinay 3191 domainChunks.push_back(dom);
140 caltinay 3143 } else {
141     cerr << "Error initializing domain!" << endl;
142     myError = 2;
143     }
144 caltinay 5184 }
145     #endif
146     #if USE_DUDLEY
147     else if (dynamic_cast<const dudley::MeshAdapter*>(domain)) {
148     DomainChunk_ptr dom(new FinleyDomain());
149     if (dom->initFromEscript(domain)) {
150     if (mpiSize > 1)
151     dom->reorderGhostZones(mpiRank);
152     domainChunks.push_back(dom);
153     } else {
154     cerr << "Error initializing domain!" << endl;
155     myError = 2;
156     }
157     }
158     #endif
159     #if USE_RIPLEY
160     else if (dynamic_cast<const ripley::RipleyDomain*>(domain)) {
161 caltinay 3792 DomainChunk_ptr dom(new RipleyDomain());
162     if (dom->initFromEscript(domain)) {
163     if (mpiSize > 1)
164     dom->reorderGhostZones(mpiRank);
165     domainChunks.push_back(dom);
166     } else {
167     cerr << "Error initializing domain!" << endl;
168     myError = 2;
169     }
170 caltinay 5184 }
171     #endif
172     #if USE_SPECKLEY
173     else if (dynamic_cast<const speckley::SpeckleyDomain*>(domain)) {
174 sshaw 5123 DomainChunk_ptr dom(new SpeckleyDomain());
175     if (dom->initFromEscript(domain)) {
176     if (mpiSize > 1)
177     dom->reorderGhostZones(mpiRank);
178     domainChunks.push_back(dom);
179     } else {
180     cerr << "Error initializing domain!" << endl;
181     myError = 2;
182     }
183 caltinay 5184 }
184     #endif
185     else {
186 caltinay 3143 cerr << "Unsupported domain type!" << endl;
187 caltinay 3128 myError = 2;
188     }
189 caltinay 3107 }
190    
191 caltinay 3128 if (mpiSize > 1) {
192 caltinay 6119 #if WEIPA_HAVE_MPI
193 caltinay 3128 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
194     #else
195     gError = myError;
196     #endif
197     } else {
198     gError = myError;
199 caltinay 2810 }
200 caltinay 2183
201 caltinay 3128 if (gError>1) {
202 caltinay 3191 domainChunks.clear();
203 caltinay 3128 } else if (gError==0) {
204     // Convert mesh data to variables
205     convertMeshVariables();
206     }
207     return (gError==0);
208    
209     #else // VISIT_PLUGIN
210     return false;
211     #endif
212     }
213    
214     //
215     //
216     //
217     bool EscriptDataset::addData(escript::Data& data, const string name,
218     const string units)
219     {
220     #ifndef VISIT_PLUGIN
221     bool success = true;
222    
223     // fail if no domain has been set
224 caltinay 3191 if (domainChunks.size() == 0) {
225 caltinay 3128 success = false;
226     } else {
227     // initialize variable
228     VarInfo vi;
229     vi.varName = name;
230     vi.units = units;
231    
232     DataVar_ptr var(new DataVar(vi.varName));
233 caltinay 3191 if (var->initFromEscript(data, domainChunks[0])) {
234     vi.dataChunks.push_back(var);
235 caltinay 3128 updateSampleDistribution(vi);
236     vi.valid = true;
237 caltinay 3107 } else {
238 caltinay 3128 var.reset();
239     vi.valid = false;
240 caltinay 3107 }
241 caltinay 3128 variables.push_back(vi);
242 caltinay 2183 }
243 caltinay 3128 return success;
244 caltinay 2810
245 caltinay 3128 #else // VISIT_PLUGIN
246 caltinay 3095 return false;
247     #endif
248 caltinay 2183 }
249    
250     //
251     //
252     //
253 caltinay 3191 bool EscriptDataset::loadNetCDF(const string domainFilePattern,
254 caltinay 2810 const StringVec& varFiles,
255 caltinay 3191 const StringVec& varNames, int nChunks)
256 caltinay 2183 {
257 caltinay 3107 // sanity check
258     if (varFiles.size() != varNames.size()) {
259     return false;
260     }
261    
262     // load the domain files
263 caltinay 3191 if (!loadDomain(domainFilePattern, nChunks)) {
264 caltinay 2822 return false;
265     }
266    
267 caltinay 3107 // load the variables
268 caltinay 2847 StringVec::const_iterator fileIt = varFiles.begin();
269     StringVec::const_iterator nameIt = varNames.begin();
270     for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
271 caltinay 3128 loadData(*fileIt, *nameIt, "");
272 caltinay 2183 }
273    
274     return true;
275     }
276    
277     //
278 caltinay 3091 // Load only variables using provided domain
279 caltinay 2183 //
280 caltinay 3191 bool EscriptDataset::loadNetCDF(const DomainChunks& domain,
281 caltinay 2810 const StringVec& varFiles,
282     const StringVec& varNames)
283 caltinay 2183 {
284 caltinay 3107 // sanity check
285     if (varFiles.size() != varNames.size()) {
286     return false;
287     }
288    
289     // set the domain
290 caltinay 3191 if (!setExternalDomain(domain)) {
291 caltinay 2822 return false;
292     }
293    
294 caltinay 3107 // load the variables
295 caltinay 2183 StringVec::const_iterator fileIt = varFiles.begin();
296     StringVec::const_iterator nameIt = varNames.begin();
297     for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
298 caltinay 3128 loadData(*fileIt, *nameIt, "");
299 caltinay 2183 }
300    
301     return true;
302     }
303    
304 caltinay 2810 //
305     //
306     //
307     bool EscriptDataset::saveSilo(string fileName, bool useMultiMesh)
308     {
309 caltinay 6144 #if ESYS_HAVE_SILO
310 caltinay 3191 if (domainChunks.size() == 0)
311 caltinay 2810 return false;
312 caltinay 2183
313 caltinay 2810 const char* blockDirFmt = "/block%04d";
314     string siloPath;
315     DBfile* dbfile = NULL;
316 caltinay 3116 int driver = DB_HDF5; // prefer HDF5 if available
317     //FIXME: Silo's HDF5 driver and NetCDF 4 are currently incompatible because
318     //NetCDF calls H5close() when all its files are closed.
319     //Unidata has been contacted, Ticket ID: YTC-894489.
320     //When this issue is resolved, remove the following line.
321     driver = DB_PDB;
322 caltinay 6119 #if WEIPA_HAVE_MPI
323 caltinay 2810 PMPIO_baton_t* baton = NULL;
324     #endif
325    
326 caltinay 3397 if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
327 caltinay 3192 fileName+=".silo";
328     }
329    
330 caltinay 2810 if (mpiSize > 1) {
331 caltinay 6119 #if WEIPA_HAVE_MPI
332 caltinay 2810 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE,
333     mpiComm, 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
334 caltinay 3116 PMPIO_DefaultClose, (void*)&driver);
335     // try the fallback driver in case of error
336     if (!baton && driver != DB_PDB) {
337     driver = DB_PDB;
338     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE,
339     mpiComm, 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
340     PMPIO_DefaultClose, (void*)&driver);
341     }
342 caltinay 2810 if (baton) {
343     char str[64];
344     snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, mpiRank));
345     siloPath = str;
346     dbfile = (DBfile*) PMPIO_WaitForBaton(
347     baton, fileName.c_str(), siloPath.c_str());
348     }
349     #endif
350     } else {
351     dbfile = DBCreate(fileName.c_str(), DB_CLOBBER, DB_LOCAL,
352 caltinay 3116 "escriptData", driver);
353     // try the fallback driver in case of error
354     if (!dbfile && driver != DB_PDB) {
355     driver = DB_PDB;
356     dbfile = DBCreate(fileName.c_str(), DB_CLOBBER, DB_LOCAL,
357     "escriptData", driver);
358     }
359 caltinay 2810 }
360    
361     if (!dbfile) {
362     cerr << "Could not create Silo file." << endl;
363 caltinay 3116 if (mpiSize > 1) {
364 caltinay 6119 #if WEIPA_HAVE_MPI
365 caltinay 3116 PMPIO_HandOffBaton(baton, dbfile);
366     PMPIO_Finish(baton);
367     #endif
368     }
369 caltinay 2810 return false;
370     }
371    
372 caltinay 3116 if (driver==DB_HDF5) {
373     // gzip level 1 already provides good compression with minimal
374     // performance penalty. Some tests showed that gzip levels >3 performed
375     // rather badly on escript data both in terms of time and space
376     DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
377     }
378    
379 caltinay 3191 DomainChunks::iterator domIt;
380 caltinay 2810 VarVector::iterator viIt;
381     int idx = 0;
382 caltinay 3191 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, idx++) {
383 caltinay 2810 if (mpiSize == 1) {
384     char str[64];
385     snprintf(str, 64, blockDirFmt, idx);
386     siloPath = str;
387     DBMkdir(dbfile, siloPath.c_str());
388     }
389 caltinay 3191 // write block of the mesh if we don't use an external domain
390     if (!externalDomain) {
391     if (! (*domIt)->writeToSilo(
392 caltinay 3623 dbfile, siloPath, meshLabels, meshUnits, wantsMeshVars)) {
393 caltinay 3128 cerr << "Error writing block " << idx
394     << " of mesh to Silo file!" << endl;
395 caltinay 2810 break;
396     }
397     }
398    
399 caltinay 3191 // write variables for current domain chunk
400 caltinay 2810 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
401     // do not attempt to write this variable if previous steps failed
402 caltinay 3128 if (!viIt->valid) continue;
403 caltinay 3191 DataVar_ptr var = viIt->dataChunks[idx];
404 caltinay 3128 if (!var->writeToSilo(dbfile, siloPath, viIt->units)) {
405 caltinay 2810 cerr << "Error writing block " << idx << " of '"
406     << var->getName() << "' to Silo file!" << endl;
407 caltinay 3128 viIt->valid = false;
408 caltinay 2810 }
409     }
410     }
411    
412     // rank 0 writes additional data that describe how the parts fit together
413     if (mpiRank == 0) {
414     if (useMultiMesh) {
415 caltinay 3191 const StringVec& meshNames = domainChunks[0]->getMeshNames();
416 caltinay 2810 StringVec::const_iterator it;
417     for (it = meshNames.begin(); it != meshNames.end(); it++)
418     putSiloMultiMesh(dbfile, *it);
419    
420 caltinay 3623 if (wantsMeshVars) {
421     DBMkdir(dbfile, MESH_VARS);
422     for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++)
423     putSiloMultiVar(dbfile, *viIt, true);
424     }
425 caltinay 2810
426     for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
427 caltinay 3128 if (!viIt->valid) continue;
428 caltinay 3191 DataVar_ptr var = viIt->dataChunks[0];
429 caltinay 2810 if (var->getRank() < 2)
430     putSiloMultiVar(dbfile, *viIt);
431     else
432     putSiloMultiTensor(dbfile, *viIt);
433     }
434     }
435    
436     vector<char*> tensorNames;
437     vector<string> tensorDefStrings;
438     vector<char*> tensorDefs;
439    
440     // collect tensors for their Silo definitions
441     for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
442 caltinay 3128 if (!viIt->valid) continue;
443 caltinay 3191 DataVar_ptr var = viIt->dataChunks[0];
444 caltinay 2810 if (var->getRank() == 2) {
445     tensorDefStrings.push_back(var->getTensorDef());
446     tensorNames.push_back((char*)var->getName().c_str());
447     }
448     }
449 caltinay 5794 for (vector<string>::const_iterator sit=tensorDefStrings.begin();
450     sit != tensorDefStrings.end(); sit++)
451     tensorDefs.push_back((char*)sit->c_str());
452 caltinay 2810
453     if (tensorDefs.size()) {
454 caltinay 2877 DBSetDir(dbfile, "/");
455 caltinay 2810 DBoptlist* optList = DBMakeOptlist(2);
456     DBAddOption(optList, DBOPT_CYCLE, &cycle);
457     DBAddOption(optList, DBOPT_DTIME, &time);
458     vector<DBoptlist*> defOpts(tensorDefs.size(), optList);
459     vector<int> defTypes(tensorDefs.size(), DB_VARTYPE_TENSOR);
460     DBPutDefvars(dbfile, "tensors", tensorDefs.size(), &tensorNames[0],
461     &defTypes[0], &tensorDefs[0], &defOpts[0]);
462     DBFreeOptlist(optList);
463     }
464     }
465    
466     if (mpiSize > 1) {
467 caltinay 6119 #if WEIPA_HAVE_MPI
468 caltinay 2810 PMPIO_HandOffBaton(baton, dbfile);
469     PMPIO_Finish(baton);
470     #endif
471     } else {
472     DBClose(dbfile);
473     }
474    
475     return true;
476    
477 caltinay 6144 #else // !ESYS_HAVE_SILO
478 sshaw 4695 std::cerr << "WARNING: saving to silo file requested but escript was not built"
479     " with silo support";
480 caltinay 2810 return false;
481     #endif
482     }
483    
484 caltinay 2183 //
485     //
486     //
487 caltinay 2810 bool EscriptDataset::saveVTK(string fileName)
488 caltinay 2183 {
489 caltinay 3191 if (domainChunks.size() == 0)
490 caltinay 2810 return false;
491    
492 caltinay 3192 map<string,VarVector> varsPerMesh;
493 caltinay 2810
494 caltinay 3192 // determine meshes and variables to write
495 caltinay 2886 VarVector::iterator viIt;
496     for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
497     // skip empty variable
498     int numSamples = accumulate(viIt->sampleDistribution.begin(),
499     viIt->sampleDistribution.end(), 0);
500     if (numSamples == 0 || !viIt->valid) {
501     continue;
502     }
503 caltinay 3192 string meshName = viIt->dataChunks[0]->getMeshName();
504 caltinay 3347
505     // ranks with 0 samples can't know the correct mesh name.
506     // We assume rank 0 always has samples, if this turns out to be a
507     // wrong assumption then a bit more work is needed to get the correct
508     // mesh name to all ranks.
509 caltinay 6119 #if WEIPA_HAVE_MPI
510 caltinay 3347 if (mpiSize > 1) {
511     char name[100];
512     if (mpiRank == 0) {
513     strncpy(&name[0], meshName.c_str(), 100);
514     }
515    
516     MPI_Bcast(name, 100, MPI_CHAR, 0, mpiComm);
517     meshName = name;
518     }
519     #endif
520    
521 caltinay 3192 map<string,VarVector>::iterator it = varsPerMesh.find(meshName);
522     if (it != varsPerMesh.end()) {
523     it->second.push_back(*viIt);
524 caltinay 2886 } else {
525 caltinay 3192 VarVector v;
526     v.push_back(*viIt);
527     varsPerMesh[meshName] = v;
528 caltinay 2886 }
529 caltinay 3192 }
530 caltinay 2886
531 caltinay 3397 if (fileName.length() < 5 || fileName.compare(fileName.length()-4, 4, ".vtu") != 0) {
532 caltinay 3192 fileName+=".vtu";
533     }
534    
535     bool ret = true;
536     if (varsPerMesh.empty()) {
537     // no valid variables so just write default mesh
538     ret = saveVTKsingle(fileName, "Elements", VarVector());
539     } else {
540     // write one file per required mesh
541     string newName(fileName);
542     string filePrefix(fileName.substr(0, fileName.length()-4));
543     bool prependMeshName=(varsPerMesh.size()>1);
544     map<string,VarVector>::const_iterator vpmIt;
545     for (vpmIt=varsPerMesh.begin(); vpmIt!=varsPerMesh.end(); vpmIt++) {
546     if (prependMeshName) {
547     newName=filePrefix+"_"+vpmIt->first+".vtu";
548     }
549     // attempt to write all files even if one fails
550     if (!saveVTKsingle(newName, vpmIt->first, vpmIt->second)) {
551     ret = false;
552     }
553     }
554     }
555     return ret;
556     }
557    
558     //
559     //
560     //
561     void EscriptDataset::setMeshLabels(const string x, const string y, const string z)
562     {
563     meshLabels.clear();
564     meshLabels.push_back(x);
565     meshLabels.push_back(y);
566     if (z.length()>0)
567     meshLabels.push_back(z);
568     }
569    
570     //
571     //
572     //
573     void EscriptDataset::setMeshUnits(const string x, const string y, const string z)
574     {
575     meshUnits.clear();
576     meshUnits.push_back(x);
577     meshUnits.push_back(y);
578     if (z.length()>0)
579     meshUnits.push_back(z);
580     }
581    
582     //
583     //
584     //
585     bool EscriptDataset::saveVTKsingle(const string& fileName,
586     const string& meshName,
587     const VarVector& vars)
588     {
589 caltinay 4320 #ifndef VISIT_PLUGIN
590 caltinay 3192 VarVector nodalVars, cellVars;
591     VarVector::const_iterator viIt;
592     for (viIt = vars.begin(); viIt != vars.end(); viIt++) {
593     const DataChunks& varChunks = viIt->dataChunks;
594 caltinay 3191 if (varChunks[0]->isNodeCentered()) {
595 caltinay 2886 nodalVars.push_back(*viIt);
596     } else {
597     cellVars.push_back(*viIt);
598     }
599     }
600    
601 caltinay 3623 if (wantsMeshVars) {
602     // add mesh variables if requested
603     for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
604     DataVar_ptr var = viIt->dataChunks[0];
605     if (meshName == var->getMeshName()) {
606     VarInfo vi = *viIt;
607     vi.varName = string(MESH_VARS)+vi.varName;
608     if (var->isNodeCentered()) {
609     nodalVars.push_back(vi);
610     } else {
611     cellVars.push_back(vi);
612     }
613 caltinay 2886 }
614     }
615     }
616    
617 caltinay 3191 DomainChunks::iterator domIt;
618 caltinay 2886 int gNumPoints;
619     int gNumCells = 0;
620     int gCellSizeAndType[2] = { 0, 0 };
621    
622 caltinay 4482 boost::scoped_ptr<FileWriter> fw(NULL);
623 caltinay 2886
624 caltinay 3091 if (mpiSize > 1) {
625 caltinay 6119 #if WEIPA_HAVE_MPI
626 caltinay 4482 fw.reset(new FileWriter(mpiComm));
627 caltinay 3191 domainChunks[0]->removeGhostZones(mpiRank);
628     ElementData_ptr elements = domainChunks[0]->getElementsByName(meshName);
629 caltinay 3347 int myNumCells = 0;
630     if (elements) {
631     myNumCells = elements->getNumElements();
632     }
633 caltinay 3091 MPI_Reduce(&myNumCells, &gNumCells, 1, MPI_INT, MPI_SUM, 0, mpiComm);
634    
635 caltinay 3347 int myCellSizeAndType[2] = { 0, 0 };
636     if (elements) {
637     myCellSizeAndType[0] = elements->getNodesPerElement();
638     myCellSizeAndType[1] = elements->getType();
639     }
640 caltinay 3091
641     // rank 0 needs to know element type and size but it's possible that
642     // this information is only available on other ranks (value=0) so
643     // retrieve it
644     MPI_Reduce(&myCellSizeAndType, &gCellSizeAndType, 2, MPI_INT, MPI_MAX,
645     0, mpiComm);
646 caltinay 2886 #endif
647 caltinay 3091 } else {
648 caltinay 4482 fw.reset(new FileWriter());
649 caltinay 3091 int idx = 0;
650 caltinay 3191 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, idx++) {
651     if (domainChunks.size() > 1)
652     (*domIt)->removeGhostZones(idx);
653     ElementData_ptr elements = (*domIt)->getElementsByName(meshName);
654 caltinay 3091 gNumCells += elements->getNumElements();
655     if (gCellSizeAndType[0] == 0)
656     gCellSizeAndType[0] = elements->getNodesPerElement();
657     if (gCellSizeAndType[1] == 0)
658     gCellSizeAndType[1] = elements->getType();
659     }
660     }
661 caltinay 2886
662 caltinay 3191 gNumPoints = domainChunks[0]->getNodes()->getGlobalNumNodes();
663 caltinay 2886
664     ostringstream oss;
665     oss.setf(ios_base::scientific, ios_base::floatfield);
666    
667 caltinay 3107 if (!fw->openFile(fileName)) {
668     return false;
669     }
670    
671 caltinay 2886 if (mpiRank == 0) {
672     #ifdef _DEBUG
673     cout << meshName << ": pts=" << gNumPoints << ", cells=" << gNumCells
674     << ", cellsize=" << gCellSizeAndType[0] << ", type="
675     << gCellSizeAndType[1] << ", numNodalVars=" << nodalVars.size()
676     << ", numCellVars=" << cellVars.size()
677     << endl;
678     #endif
679    
680     oss << "<?xml version=\"1.0\"?>" << endl;
681 caltinay 3101 oss << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
682     if (mdSchema.length()>0) {
683     oss << " " << mdSchema;
684     }
685     oss << ">" << endl;
686     if (mdString.length()>0) {
687 caltinay 3147 oss << "<MetaData>" << endl << mdString << endl
688     << "</MetaData>" << endl;
689 caltinay 3101 }
690 caltinay 2886 oss << "<UnstructuredGrid>" << endl;
691    
692     // write time and cycle values
693     oss << "<FieldData>" << endl;
694     oss << "<DataArray Name=\"TIME\" type=\"Float64\" format=\"ascii\" NumberOfTuples=\"1\">" << endl;
695     oss << time << endl;
696     oss << "</DataArray>" << endl << "</FieldData>" << endl;
697     oss << "<FieldData>" << endl;
698     oss << "<DataArray Name=\"CYCLE\" type=\"Int32\" format=\"ascii\" NumberOfTuples=\"1\">" << endl;
699     oss << cycle << endl;
700     oss << "</DataArray>" << endl << "</FieldData>" << endl;
701    
702     // write mesh header
703     oss << "<Piece NumberOfPoints=\"" << gNumPoints
704     << "\" NumberOfCells=\"" << gNumCells << "\">" << endl;
705     oss << "<Points>" << endl;
706     oss << "<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">" << endl;
707     }
708    
709     //
710     // coordinates (note that we are using the original nodes)
711 caltinay 3091 //
712     int blockNum = (mpiSize>1 ? mpiRank : 0);
713 caltinay 3191 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, blockNum++) {
714     (*domIt)->getNodes()->writeCoordinatesVTK(oss, blockNum);
715 caltinay 2810 }
716 caltinay 2886
717 caltinay 3107 if (!fw->writeOrdered(oss)) {
718     cerr << "Warning, ignoring file write error!" << endl;
719     }
720 caltinay 2886
721     if (mpiRank == 0) {
722     oss << "</DataArray>" << endl << "</Points>" << endl;
723     oss << "<Cells>" << endl;
724     oss << "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">" << endl;
725     }
726    
727     //
728     // connectivity
729     //
730 caltinay 3191 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++) {
731     ElementData_ptr el = (*domIt)->getElementsByName(meshName);
732 caltinay 3347 if (el) {
733     el->writeConnectivityVTK(oss);
734     }
735 caltinay 2886 }
736    
737 caltinay 3107 if (!fw->writeOrdered(oss)) {
738     cerr << "Warning, ignoring file write error!" << endl;
739     }
740 caltinay 2886
741     //
742     // offsets & types
743     //
744     if (mpiRank == 0) {
745     oss << "</DataArray>" << endl;
746     oss << "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">" << endl;
747     for (int i=1; i < gNumCells+1; i++) {
748     oss << i*gCellSizeAndType[0] << endl;
749     }
750     oss << "</DataArray>" << endl;
751     oss << "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">" << endl;
752     for (int i=1; i < gNumCells+1; i++) {
753     oss << gCellSizeAndType[1] << endl;
754     }
755     oss << "</DataArray>" << endl << "</Cells>" << endl;
756 caltinay 3107 if (!fw->writeShared(oss)) {
757     cerr << "Warning, ignoring file write error!" << endl;
758     }
759 caltinay 2886 }
760    
761     // now write all variables - first the nodal data, then cell data
762    
763     // write nodal data if any
764 caltinay 2940 if (!nodalVars.empty()) {
765 caltinay 2886 if (mpiRank == 0)
766     oss << "<PointData>" << endl;
767     for (viIt = nodalVars.begin(); viIt != nodalVars.end(); viIt++) {
768     writeVarToVTK(*viIt, oss);
769 caltinay 3107 if (!fw->writeOrdered(oss)) {
770     cerr << "Warning, ignoring file write error!" << endl;
771     }
772 caltinay 2886 if (mpiRank == 0)
773     oss << "</DataArray>" << endl;
774     }
775     if (mpiRank == 0)
776     oss << "</PointData>" << endl;
777     }
778    
779     // write cell data if any
780 caltinay 2940 if (!cellVars.empty()) {
781 caltinay 2886 if (mpiRank == 0)
782     oss << "<CellData>" << endl;
783     for (viIt = cellVars.begin(); viIt != cellVars.end(); viIt++) {
784     writeVarToVTK(*viIt, oss);
785 caltinay 3107 if (!fw->writeOrdered(oss)) {
786     cerr << "Warning, ignoring file write error!" << endl;
787     }
788 caltinay 2886 if (mpiRank == 0)
789     oss << "</DataArray>" << endl;
790     }
791     if (mpiRank == 0)
792     oss << "</CellData>" << endl;
793     }
794    
795     if (mpiRank == 0) {
796     oss << "</Piece>" << endl << "</UnstructuredGrid>" << endl
797     << "</VTKFile>" << endl;
798 caltinay 3107 if (!fw->writeShared(oss)) {
799     cerr << "Warning, ignoring file write error!" << endl;
800     }
801 caltinay 2886 }
802    
803 caltinay 3107 fw->close();
804 caltinay 2886 return true;
805 caltinay 4320 #else // VISIT_PLUGIN
806     return false;
807     #endif
808 caltinay 2810 }
809    
810     //
811     //
812     //
813 caltinay 2886 void EscriptDataset::writeVarToVTK(const VarInfo& varInfo, ostream& os)
814     {
815 caltinay 3191 const DataChunks& varChunks = varInfo.dataChunks;
816     int rank = varChunks[0]->getRank();
817 caltinay 2886 int numComps = 1;
818     if (rank > 0)
819     numComps *= 3;
820     if (rank > 1)
821     numComps *= 3;
822    
823     if (mpiRank == 0) {
824     os << "<DataArray Name=\"" << varInfo.varName
825     << "\" type=\"Float64\" NumberOfComponents=\"" << numComps
826     << "\" format=\"ascii\">" << endl;
827     }
828    
829 caltinay 3091 // this is required in case we read a dataset with more than one chunk on
830     // one rank
831     int blockNum = (mpiSize>1 ? mpiRank : 0);
832 caltinay 3191 DataChunks::const_iterator blockIt;
833     for (blockIt = varChunks.begin(); blockIt != varChunks.end(); blockIt++, blockNum++) {
834 caltinay 3091 (*blockIt)->writeToVTK(os, blockNum);
835 caltinay 2886 }
836     }
837    
838     //
839 caltinay 3091 // Sets the domain from dump files.
840     //
841 caltinay 3191 bool EscriptDataset::loadDomain(const string filePattern, int nChunks)
842 caltinay 3091 {
843     int myError = 0, gError;
844    
845 caltinay 3191 if (mpiSize > 1 && nChunks != mpiSize) {
846     cerr << "Cannot load " << nChunks << " chunks on " << mpiSize
847 caltinay 3091 << " MPI ranks!" << endl;
848     myError = 1;
849    
850 caltinay 3191 } else if (domainChunks.size() > 0) {
851 caltinay 3091 cerr << "Domain has already been set!" << endl;
852     myError = 1;
853    
854     } else {
855 caltinay 5184 #ifdef USE_FINLEY
856 caltinay 3091 char* str = new char[filePattern.length()+10];
857 caltinay 3143 // FIXME: This assumes a finley domain!
858 caltinay 3091 if (mpiSize > 1) {
859 caltinay 3143 DomainChunk_ptr chunk(new FinleyDomain());
860 caltinay 3091 sprintf(str, filePattern.c_str(), mpiRank);
861 caltinay 3143 string domainFile = str;
862     if (chunk->initFromFile(domainFile)) {
863     chunk->reorderGhostZones(mpiRank);
864 caltinay 3191 domainChunks.push_back(chunk);
865 caltinay 2822 } else {
866 caltinay 3091 cerr << "Error initializing domain!" << endl;
867 caltinay 3147 myError = 2;
868 caltinay 2822 }
869 caltinay 3091 } else {
870 caltinay 3191 for (int idx=0; idx < nChunks; idx++) {
871 caltinay 3143 DomainChunk_ptr chunk(new FinleyDomain());
872 caltinay 3091 sprintf(str, filePattern.c_str(), idx);
873 caltinay 3143 string domainFile = str;
874     if (chunk->initFromFile(domainFile)) {
875 caltinay 3191 if (nChunks > 1)
876 caltinay 3143 chunk->reorderGhostZones(idx);
877 caltinay 3191 domainChunks.push_back(chunk);
878 caltinay 3091 } else {
879     cerr << "Error initializing domain block " << idx << endl;
880 caltinay 3147 myError = 2;
881 caltinay 3091 break;
882     }
883     }
884 caltinay 2822 }
885 caltinay 3091 delete[] str;
886 caltinay 5184 #else // USE_FINLEY
887     cerr << "weipa not compiled with finley support." << endl;
888     myError = 3;
889     #endif
890 caltinay 2183 }
891 caltinay 3091
892     if (mpiSize > 1) {
893 caltinay 6119 #if WEIPA_HAVE_MPI
894 caltinay 3091 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
895 gross 3093 #else
896     gError = myError;
897 caltinay 3091 #endif
898     } else {
899     gError = myError;
900     }
901    
902 caltinay 3147 if (gError>1) {
903 caltinay 3191 domainChunks.clear();
904 caltinay 3147 } else if (gError==0) {
905 caltinay 3091 // Convert mesh data to variables
906     convertMeshVariables();
907     }
908 caltinay 3147 return (gError==0);
909 caltinay 2183 }
910    
911 caltinay 2810 //
912 caltinay 3091 // Sets an already converted domain.
913 caltinay 2810 //
914 caltinay 3191 bool EscriptDataset::setExternalDomain(const DomainChunks& domain)
915 caltinay 3091 {
916     int myError = 0, gError;
917    
918     if (mpiSize > 1 && domain.size() > 1) {
919     cerr << "Can only add one domain block per rank when using MPI!"
920     << endl;
921     myError = 1;
922    
923 caltinay 3191 } else if (domainChunks.size() > 0) {
924 caltinay 3091 cerr << "Domain has already been set!" << endl;
925     myError = 1;
926    
927     }
928    
929     if (mpiSize > 1) {
930 caltinay 6119 #if WEIPA_HAVE_MPI
931 caltinay 3091 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
932 gross 3093 #else
933     gError = myError;
934 caltinay 3091 #endif
935     } else {
936     gError = myError;
937     }
938    
939     if (!gError) {
940 caltinay 3191 externalDomain = true;
941     domainChunks = domain;
942 caltinay 3091 }
943    
944     return !gError;
945     }
946    
947 caltinay 2810 //
948 caltinay 3091 //
949     //
950 caltinay 3128 bool EscriptDataset::loadData(const string filePattern, const string name,
951     const string units)
952 caltinay 2183 {
953 caltinay 3091 int myError = 0, gError;
954    
955     // fail if no domain has been set
956 caltinay 3191 if (domainChunks.size() == 0) {
957 caltinay 3091 gError = 1;
958    
959     } else {
960     // initialize variable
961     VarInfo vi;
962     vi.varName = name;
963 caltinay 3128 vi.units = units;
964 caltinay 3091 vi.valid = true;
965     char* str = new char[filePattern.length()+10];
966    
967     // read all parts of the variable
968 caltinay 3191 DomainChunks::iterator domIt;
969 caltinay 3091 int idx = (mpiSize > 1) ? mpiRank : 0;
970 caltinay 3191 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, idx++) {
971 caltinay 3091 sprintf(str, filePattern.c_str(), idx);
972     string dfile = str;
973     DataVar_ptr var(new DataVar(name));
974 caltinay 3191 if (var->initFromFile(dfile, *domIt))
975     vi.dataChunks.push_back(var);
976 caltinay 3091 else {
977     cerr << "Error reading " << dfile << endl;
978     myError = 1;
979     break;
980     }
981     }
982     delete[] str;
983    
984     if (mpiSize > 1) {
985 caltinay 6119 #if WEIPA_HAVE_MPI
986 caltinay 3091 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
987 gross 3093 #else
988     gError = myError;
989 caltinay 2847 #endif
990 caltinay 3091 } else {
991     gError = myError;
992     }
993 caltinay 2847
994 caltinay 3091 if (!gError) {
995     // only add variable if all chunks have been read without error
996     updateSampleDistribution(vi);
997     variables.push_back(vi);
998     }
999 caltinay 2847 }
1000    
1001 caltinay 3091 return !gError;
1002 caltinay 2847 }
1003    
1004 caltinay 2886 //
1005     //
1006     //
1007     void EscriptDataset::convertMeshVariables()
1008     {
1009 caltinay 3191 const StringVec& varNames = domainChunks[0]->getVarNames();
1010 caltinay 2886 StringVec::const_iterator it;
1011     for (it = varNames.begin(); it != varNames.end(); it++) {
1012     VarInfo vi;
1013     vi.varName = *it;
1014     vi.valid = true;
1015     // get all parts of current variable
1016 caltinay 3191 DomainChunks::iterator domIt;
1017     for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++) {
1018     DataVar_ptr var = (*domIt)->getDataVarByName(*it);
1019 caltinay 3143 if (var != NULL) {
1020 caltinay 3191 vi.dataChunks.push_back(var);
1021 caltinay 2886 } else {
1022     cerr << "Error converting mesh variable " << *it << endl;
1023     vi.valid = false;
1024     break;
1025     }
1026     }
1027     updateSampleDistribution(vi);
1028     meshVariables.push_back(vi);
1029     }
1030     }
1031    
1032     // retrieves the number of samples at each block - used to determine which
1033 caltinay 2847 // blocks contribute to given variable if any.
1034     void EscriptDataset::updateSampleDistribution(VarInfo& vi)
1035     {
1036     IntVec sampleDist;
1037 caltinay 3191 const DataChunks& varChunks = vi.dataChunks;
1038 caltinay 2847
1039     if (mpiSize > 1) {
1040 caltinay 6119 #if WEIPA_HAVE_MPI
1041 caltinay 3191 int myNumSamples = varChunks[0]->getNumberOfSamples();
1042 caltinay 2877 sampleDist.insert(sampleDist.end(), mpiSize, 0);
1043     MPI_Allgather(
1044     &myNumSamples, 1, MPI_INT, &sampleDist[0], 1, MPI_INT, mpiComm);
1045 caltinay 2847 #endif
1046     } else {
1047 caltinay 3191 DataChunks::const_iterator it;
1048     for (it = varChunks.begin(); it != varChunks.end(); it++) {
1049 caltinay 2847 sampleDist.push_back((*it)->getNumberOfSamples());
1050     }
1051 caltinay 2183 }
1052 caltinay 2847 vi.sampleDistribution = sampleDist;
1053 caltinay 2183 }
1054    
1055     //
1056     //
1057     //
1058 caltinay 2810 void EscriptDataset::putSiloMultiMesh(DBfile* dbfile, const string& meshName)
1059 caltinay 2183 {
1060 caltinay 6144 #if ESYS_HAVE_SILO
1061 caltinay 2877 vector<int> meshtypes;
1062 caltinay 2183 vector<string> tempstrings;
1063     vector<char*> meshnames;
1064 caltinay 2810 string pathPrefix;
1065 caltinay 2877
1066 caltinay 3191 int ppIndex = domainChunks[0]->getSiloPath().find(':');
1067 caltinay 2810 if (ppIndex != string::npos) {
1068 caltinay 3191 pathPrefix = domainChunks[0]->getSiloPath().substr(0, ppIndex+1);
1069 caltinay 2810 }
1070 caltinay 2183
1071 caltinay 3191 // find a variable that is defined on this mesh to get the sample
1072 caltinay 2877 // distribution (which tells us which ranks contribute to this mesh).
1073     // Try mesh variables first, then regular ones.
1074     VarVector::const_iterator viIt;
1075     for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
1076 caltinay 3191 if (meshName == viIt->dataChunks[0]->getMeshName())
1077 caltinay 2877 break;
1078 caltinay 2183 }
1079 caltinay 2847
1080 caltinay 2877 if (viIt == meshVariables.end()) {
1081     for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
1082 caltinay 3191 if (meshName == viIt->dataChunks[0]->getMeshName())
1083 caltinay 2877 break;
1084     }
1085 caltinay 5097 // this probably means that the mesh is empty
1086     if (viIt == variables.end()) {
1087     return;
1088     }
1089 caltinay 2877 }
1090    
1091     for (size_t idx = 0; idx < viIt->sampleDistribution.size(); idx++) {
1092     if (viIt->sampleDistribution[idx] > 0) {
1093     stringstream siloPath;
1094     siloPath << pathPrefix << "/block";
1095     int prevWidth = siloPath.width(4);
1096     char prevFill = siloPath.fill('0');
1097     siloPath << right << idx;
1098     siloPath.width(prevWidth);
1099     siloPath.fill(prevFill);
1100     siloPath << "/";
1101     siloPath << meshName;
1102     tempstrings.push_back(siloPath.str());
1103     meshtypes.push_back(DB_UCDMESH);
1104     }
1105     }
1106 caltinay 5794 for (vector<string>::const_iterator sit=tempstrings.begin();
1107     sit != tempstrings.end(); sit++)
1108     meshnames.push_back((char*)sit->c_str());
1109 caltinay 2877
1110 caltinay 2847 // ignore empty mesh
1111 caltinay 2940 if (!meshnames.empty()) {
1112 caltinay 2877 DBSetDir(dbfile, "/");
1113 caltinay 2847 DBoptlist* optList = DBMakeOptlist(2);
1114     DBAddOption(optList, DBOPT_CYCLE, &cycle);
1115     DBAddOption(optList, DBOPT_DTIME, &time);
1116 caltinay 2877 DBPutMultimesh(dbfile, meshName.c_str(), meshnames.size(),
1117     &meshnames[0], &meshtypes[0], optList);
1118 caltinay 2847 DBFreeOptlist(optList);
1119     }
1120 caltinay 2183 #endif
1121     }
1122    
1123 caltinay 2810 //
1124     //
1125     //
1126     void EscriptDataset::putSiloMultiVar(DBfile* dbfile, const VarInfo& vi,
1127     bool useMeshFile)
1128 caltinay 2183 {
1129 caltinay 6144 #if ESYS_HAVE_SILO
1130 caltinay 2847 vector<int> vartypes;
1131 caltinay 2183 vector<string> tempstrings;
1132     vector<char*> varnames;
1133 caltinay 2810 string pathPrefix;
1134     if (useMeshFile) {
1135 caltinay 3191 int ppIndex = domainChunks[0]->getSiloPath().find(':');
1136 caltinay 2810 if (ppIndex != string::npos) {
1137 caltinay 3191 pathPrefix = domainChunks[0]->getSiloPath().substr(0, ppIndex+1);
1138 caltinay 2183 }
1139 caltinay 2810 }
1140    
1141 caltinay 2847 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
1142     if (vi.sampleDistribution[idx] > 0) {
1143     stringstream siloPath;
1144     siloPath << pathPrefix << "/block";
1145     int prevWidth = siloPath.width(4);
1146     char prevFill = siloPath.fill('0');
1147     siloPath << right << idx;
1148     siloPath.width(prevWidth);
1149     siloPath.fill(prevFill);
1150     siloPath << "/";
1151     siloPath << vi.varName;
1152     tempstrings.push_back(siloPath.str());
1153     vartypes.push_back(DB_UCDVAR);
1154     }
1155 caltinay 2183 }
1156 caltinay 5794 for (vector<string>::const_iterator sit=tempstrings.begin();
1157     sit != tempstrings.end(); sit++)
1158     varnames.push_back((char*)sit->c_str());
1159 caltinay 2847
1160     // ignore empty variables
1161 caltinay 2940 if (!varnames.empty()) {
1162 caltinay 2877 DBSetDir(dbfile, "/");
1163 caltinay 2847 DBoptlist* optList = DBMakeOptlist(2);
1164     DBAddOption(optList, DBOPT_CYCLE, &cycle);
1165     DBAddOption(optList, DBOPT_DTIME, &time);
1166     if (useMeshFile) {
1167     string vpath = string(MESH_VARS)+vi.varName;
1168     DBPutMultivar(dbfile, vpath.c_str(), varnames.size(),
1169     &varnames[0], &vartypes[0], optList);
1170     } else {
1171     DBPutMultivar(dbfile, vi.varName.c_str(), varnames.size(),
1172     &varnames[0], &vartypes[0], optList);
1173     }
1174     DBFreeOptlist(optList);
1175 caltinay 2190 }
1176 caltinay 2183 #endif
1177     }
1178    
1179 caltinay 2810 //
1180     //
1181     //
1182     void EscriptDataset::putSiloMultiTensor(DBfile* dbfile, const VarInfo& vi)
1183 caltinay 2183 {
1184 caltinay 6144 #if ESYS_HAVE_SILO
1185 caltinay 2810 string tensorDir = vi.varName+string("_comps/");
1186 caltinay 2183 DBSetDir(dbfile, "/");
1187     DBMkdir(dbfile, tensorDir.c_str());
1188     int one = 1;
1189 caltinay 2810 DBoptlist* optList = DBMakeOptlist(3);
1190     DBAddOption(optList, DBOPT_CYCLE, &cycle);
1191     DBAddOption(optList, DBOPT_DTIME, &time);
1192 caltinay 2183 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
1193 caltinay 3191 const IntVec& shape = vi.dataChunks[0]->getShape();
1194 caltinay 2877
1195 caltinay 2183 for (int i=0; i<shape[1]; i++) {
1196     for (int j=0; j<shape[0]; j++) {
1197     vector<string> tempstrings;
1198     vector<char*> varnames;
1199 caltinay 2877 vector<int> vartypes;
1200 caltinay 2810 stringstream comp;
1201     comp << vi.varName << "_comps/a_";
1202     comp << i;
1203     comp << j;
1204 caltinay 2877 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
1205     if (vi.sampleDistribution[idx] > 0) {
1206     stringstream siloPath;
1207     siloPath << "/block";
1208     int prevWidth = siloPath.width(4);
1209     char prevFill = siloPath.fill('0');
1210     siloPath << right << idx;
1211     siloPath.width(prevWidth);
1212     siloPath.fill(prevFill);
1213     siloPath << "/" << comp.str();
1214     tempstrings.push_back(siloPath.str());
1215     varnames.push_back((char*)tempstrings.back().c_str());
1216     vartypes.push_back(DB_UCDVAR);
1217     }
1218 caltinay 2183 }
1219 caltinay 2940 if (!varnames.empty()) {
1220 caltinay 2877 DBPutMultivar(dbfile, comp.str().c_str(), varnames.size(),
1221     &varnames[0], &vartypes[0], optList);
1222     }
1223 caltinay 2183 }
1224     }
1225     DBFreeOptlist(optList);
1226     #endif
1227     }
1228    
1229 caltinay 3037 } // namespace weipa
1230 caltinay 2187

  ViewVC Help
Powered by ViewVC 1.1.26