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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2187 - (show 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
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 namespace EscriptReader {
29
30 //
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 } // namespace EscriptReader
410

  ViewVC Help
Powered by ViewVC 1.1.26