/[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 2187 - (hide annotations)
Tue Dec 23 04:13:15 2008 UTC (12 years, 7 months ago) by caltinay
File size: 11909 byte(s)
Moved escriptreader related classes into own namespace.

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

  ViewVC Help
Powered by ViewVC 1.1.26