/[escript]/trunk/tools/libescriptreader/src/escriptreader/MPDataSet.cpp
ViewVC logotype

Annotation of /trunk/tools/libescriptreader/src/escriptreader/MPDataSet.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2183 - (hide annotations)
Fri Dec 19 03:52:50 2008 UTC (12 years, 7 months ago) by caltinay
File size: 11852 byte(s)
Added escriptreader library with tools. For more information refer to README
file or wait for user documentation.

1 caltinay 2183
2     /*******************************************************
3     *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     //
15     // MPDataSet.cpp
16     //
17     #include <escriptreader/MPDataSet.h>
18     #include <escriptreader/MeshWithElements.h>
19     #include <escriptreader/DataVar.h>
20     #include <cstring>
21     #include <netcdf.hh>
22     #if HAVE_SILO
23     #include <silo.h>
24     #endif
25    
26     using namespace std;
27    
28     //
29     // Constructor
30     //
31     MPDataSet::MPDataSet() : numParts(0), keepMesh(false)
32     {
33     }
34    
35     //
36     // Destructor
37     //
38     MPDataSet::~MPDataSet()
39     {
40     VarVector::iterator viIt;
41     for (viIt = variables.begin(); viIt != variables.end(); viIt++)
42     delete (*viIt).dataVar;
43    
44     for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++)
45     delete (*viIt).dataVar;
46    
47     if (!keepMesh) {
48     MeshBlocks::iterator meshIt;
49     for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++)
50     delete *meshIt;
51     }
52     }
53    
54     //
55     //
56     //
57     bool MPDataSet::load(const string meshFile, const StringVec& varFiles,
58     const StringVec& varNames, int nBlocks)
59     {
60     numParts = nBlocks;
61     meshFmt = meshFile;
62    
63     // initialize variables
64     StringVec::const_iterator fileIt = varFiles.begin();
65     StringVec::const_iterator nameIt = varNames.begin();
66     for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
67     VarInfo vi;
68     vi.fileName = *fileIt;
69     vi.varName = *nameIt;
70     vi.valid = true;
71     variables.push_back(vi);
72     }
73    
74     // Load the mesh files
75     if (!readMeshes()) {
76     cerr << "Reading the domain failed." << endl;
77     return false;
78     }
79    
80     // Retrieve the mesh variables
81     convertMeshVariables();
82    
83     // Load the variables
84     if (!readVariables()) {
85     cerr << "Reading the variables failed." << endl;
86     return false;
87     }
88    
89     return true;
90     }
91    
92     //
93     // Load only variables using provided mesh
94     //
95     bool MPDataSet::load(const MeshBlocks& m, const string siloFile,
96     const StringVec& varFiles, const StringVec& varNames)
97     {
98     siloMeshFile = siloFile;
99     meshBlocks = m;
100     numParts = meshBlocks.size();
101     keepMesh = true;
102    
103     // initialize variables
104     StringVec::const_iterator fileIt = varFiles.begin();
105     StringVec::const_iterator nameIt = varNames.begin();
106     for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
107     VarInfo vi;
108     vi.fileName = *fileIt;
109     vi.varName = *nameIt;
110     vi.valid = true;
111     variables.push_back(vi);
112     }
113    
114     // Load the variables
115     if (!readVariables()) {
116     cerr << "Reading the variables failed." << endl;
117     return false;
118     }
119    
120     return true;
121     }
122    
123    
124     //
125     //
126     //
127     bool MPDataSet::readMeshes()
128     {
129     bool ok = true;
130     char* str = new char[meshFmt.length()+10];
131     for (int idx=0; idx < numParts; idx++) {
132     MeshWithElements* meshPart = new MeshWithElements();
133     sprintf(str, meshFmt.c_str(), idx);
134     string meshfile = str;
135     if (meshPart->readFromNc(meshfile)) {
136     if (numParts > 1)
137     meshPart->handleGhostZones(idx);
138     meshBlocks.push_back(meshPart);
139     } else {
140     delete meshPart;
141     ok = false;
142     break;
143     }
144     }
145     delete[] str;
146     return ok;
147     }
148    
149     void MPDataSet::convertMeshVariables()
150     {
151     const StringVec& varNames = meshBlocks[0]->getVarNames();
152     StringVec::const_iterator it;
153     for (it = varNames.begin(); it != varNames.end(); it++) {
154     VarInfo vi;
155     vi.varName = *it;
156     vi.valid = true;
157     DataParts parts;
158     // get all parts of current variable
159     for (int idx=0; idx < numParts; idx++) {
160     const IntVec& data = meshBlocks[idx]->getVarDataByName(*it);
161     DataVar* var = new DataVar(*it, data, meshBlocks[idx]);
162     parts.push_back(var);
163     }
164     // if we have more than one part, assemble variable
165     if (numParts > 1) {
166     assembleVariable(vi, parts);
167     for (size_t i=0; i<parts.size(); i++)
168     delete parts[i];
169     } else
170     vi.dataVar = parts.back();
171    
172     meshVariables.push_back(vi);
173     }
174     }
175    
176     //
177     //
178     //
179     bool MPDataSet::readVariables()
180     {
181     if (variables.size() == 0)
182     return true;
183    
184     VarVector::iterator it;
185     bool ok = false;
186     for (it = variables.begin(); it != variables.end(); it++) {
187     bool curVarOk = true;
188     VarInfo& vi = (*it);
189     DataParts parts;
190     char* str = new char[vi.fileName.length()+10];
191     // read all parts of current variable
192     for (int idx=0; idx < numParts; idx++) {
193     sprintf(str, vi.fileName.c_str(), idx);
194     string dfile = str;
195     DataVar* var = new DataVar(vi.varName);
196     if (var->readFromNc(dfile))
197     parts.push_back(var);
198     else {
199     delete var;
200     curVarOk = false;
201     break;
202     }
203     }
204     delete[] str;
205     if (curVarOk) {
206     // at least one variable was read without problems
207     ok = true;
208    
209     // if we have more than one part, assemble variable
210     if (numParts > 1) {
211     assembleVariable(vi, parts);
212     for (size_t i=0; i<parts.size(); i++)
213     delete parts[i];
214     } else
215     vi.dataVar = parts.back();
216     } else {
217     vi.dataVar = NULL;
218     vi.valid = false;
219     }
220     }
221     return ok;
222     }
223    
224     //
225     //
226     //
227     void MPDataSet::assembleVariable(VarInfo& vi, const DataParts& parts)
228     {
229     DataVar* var = new DataVar(*parts[0]);
230     vi.dataVar = var;
231    
232     DataParts::const_iterator it;
233     for (it = parts.begin()+1; it != parts.end(); it++) {
234     var->append(*(*it));
235     }
236     }
237    
238     bool MPDataSet::saveAsSilo(string siloFile, bool useMultiMesh)
239     {
240     #if HAVE_SILO
241     if (numParts == 0)
242     return false;
243    
244     DBfile* dbfile;
245     dbfile = DBCreate(siloFile.c_str(), DB_CLOBBER, DB_LOCAL, "Data", DB_PDB);
246     if (!dbfile) {
247     cerr << "Could not create Silo file." << endl;
248     return false;
249     }
250    
251     MeshBlocks::iterator meshIt;
252     VarVector::iterator viIt;
253     int idx = 0;
254     for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, idx++) {
255     string siloPath("");
256     //if (numParts > 1) {
257     char str[64];
258     snprintf(str, 64, "/block%04d", idx);
259     siloPath = str;
260     DBMkdir(dbfile, siloPath.c_str());
261     //}
262     // write block of the mesh if we don't use an external mesh
263     if (!siloMeshFile.length()) {
264     if (! (*meshIt)->writeToSilo(dbfile, siloPath)) {
265     cerr << "Error writing mesh of block " << idx << " to Silo!\n";
266     break;
267     }
268     }
269    
270     // write variables for current mesh block
271     for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
272     // do not attempt to write this variable if previous steps failed
273     if (!(*viIt).valid) continue;
274     DataVar* var = (*viIt).dataVar;
275     if (!var->setMesh(*meshIt) || !var->writeToSilo(dbfile, siloPath)) {
276     cerr << "Error writing block " << idx << " of "
277     << var->getName() << " to Silo!" << endl;
278     (*viIt).valid = false;
279     }
280     }
281     }
282    
283     if (useMultiMesh) {
284     const StringVec& meshNames = meshBlocks[0]->getMeshNames();
285     StringVec::const_iterator it;
286     for (it = meshNames.begin(); it != meshNames.end(); it++)
287     putSiloMultiMesh(dbfile, *it);
288    
289     const StringVec& varNames = meshBlocks[0]->getVarNames();
290     for (it = varNames.begin(); it != varNames.end(); it++)
291     putSiloMultiVar(dbfile, *it, true);
292    
293     for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
294     if (!(*viIt).valid) continue;
295     DataVar* var = (*viIt).dataVar;
296     if (var->getRank() < 2)
297     putSiloMultiVar(dbfile, var->getName());
298     else
299     putSiloMultiTensor(dbfile, var);
300     }
301     }
302    
303     vector<char*> tensorNames;
304     vector<string> tensorDefStrings;
305     vector<char*> tensorDefs;
306    
307     // collect tensors for their Silo definitions
308     for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
309     if (!(*viIt).valid) continue;
310     DataVar* var = (*viIt).dataVar;
311     if (var->getRank() == 2) {
312     tensorDefStrings.push_back(var->getTensorDef());
313     tensorDefs.push_back((char*)tensorDefStrings.back().c_str());
314     tensorNames.push_back((char*)var->getName().c_str());
315     }
316     }
317    
318     if (tensorDefs.size()) {
319     vector<int> defTypes(tensorDefs.size(), DB_VARTYPE_TENSOR);
320     DBPutDefvars(dbfile, "tensors", tensorDefs.size(), &tensorNames[0],
321     &defTypes[0], &tensorDefs[0], NULL);
322     }
323    
324     DBClose(dbfile);
325     return true;
326    
327     #else // !HAVE_SILO
328     return false;
329     #endif
330     }
331    
332     void MPDataSet::putSiloMultiMesh(DBfile* dbfile, string meshName)
333     {
334     #if HAVE_SILO
335     vector<int> meshtypes(meshBlocks.size(), DB_UCDMESH);
336     vector<string> tempstrings;
337     vector<char*> meshnames;
338     vector<Mesh*>::iterator it;
339    
340     for (size_t idx = 0; idx < meshBlocks.size(); idx++) {
341     string siloPath = meshBlocks[idx]->getSiloPath();
342     tempstrings.push_back(siloPath + string("/") + meshName);
343     meshnames.push_back((char*)tempstrings.back().c_str());
344     }
345     DBPutMultimesh(dbfile, meshName.c_str(), meshBlocks.size(), &meshnames[0],
346     &meshtypes[0], NULL);
347     #endif
348     }
349    
350     void MPDataSet::putSiloMultiVar(DBfile* dbfile, string varName,
351     bool useMeshFile)
352     {
353     #if HAVE_SILO
354     vector<int> vartypes(meshBlocks.size(), DB_UCDVAR);
355     vector<string> tempstrings;
356     vector<char*> varnames;
357     for (size_t idx = 0; idx < meshBlocks.size(); idx++) {
358     string siloPath;
359     if (useMeshFile)
360     siloPath = meshBlocks[idx]->getSiloPath();
361     else {
362     char str[64];
363     snprintf(str, 64, "/block%04d", static_cast<int>(idx));
364     siloPath = str;
365     }
366     tempstrings.push_back(siloPath + string("/") + varName);
367     varnames.push_back((char*)tempstrings.back().c_str());
368     }
369     DBPutMultivar(dbfile, varName.c_str(), meshBlocks.size(), &varnames[0],
370     &vartypes[0], NULL);
371     #endif
372     }
373    
374     void MPDataSet::putSiloMultiTensor(DBfile* dbfile, const DataVar* var)
375     {
376     #if HAVE_SILO
377     string tensorDir = var->getName()+string("_comps/");
378     DBSetDir(dbfile, "/");
379     DBMkdir(dbfile, tensorDir.c_str());
380     int one = 1;
381     DBoptlist* optList = DBMakeOptlist(1);
382     DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
383     vector<int> vartypes(meshBlocks.size(), DB_UCDVAR);
384     const IntVec& shape = var->getShape();
385     for (int i=0; i<shape[1]; i++) {
386     for (int j=0; j<shape[0]; j++) {
387     vector<string> tempstrings;
388     vector<char*> varnames;
389     char comp[255];
390     snprintf(comp, 255, "%s_comps/a_%d%d", var->getName().c_str(), i,j);
391     for (size_t idx = 0; idx < meshBlocks.size(); idx++) {
392     string siloPath;
393     char str[64];
394     snprintf(str, 64, "/block%04d", static_cast<int>(idx));
395     siloPath = str;
396     tempstrings.push_back(siloPath + string("/") + string(comp));
397     varnames.push_back((char*)tempstrings.back().c_str());
398     }
399     DBPutMultivar(dbfile, comp, meshBlocks.size(), &varnames[0],
400     &vartypes[0], optList);
401     }
402     }
403     DBFreeOptlist(optList);
404     #endif
405     }
406    

  ViewVC Help
Powered by ViewVC 1.1.26