/[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 2190 - (show annotations)
Wed Dec 24 02:04:16 2008 UTC (12 years, 6 months ago) by caltinay
File size: 12114 byte(s)
escriptreader: write special mesh variables into a "subdirectory" in the Silo
file. That way, "real" variables are better visible.


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 const char* MESH_VARS = "mesh_vars/";
31
32 //
33 // Constructor
34 //
35 MPDataSet::MPDataSet() : numParts(0), keepMesh(false)
36 {
37 }
38
39 //
40 // Destructor
41 //
42 MPDataSet::~MPDataSet()
43 {
44 VarVector::iterator viIt;
45 for (viIt = variables.begin(); viIt != variables.end(); viIt++)
46 delete (*viIt).dataVar;
47
48 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++)
49 delete (*viIt).dataVar;
50
51 if (!keepMesh) {
52 MeshBlocks::iterator meshIt;
53 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++)
54 delete *meshIt;
55 }
56 }
57
58 //
59 //
60 //
61 bool MPDataSet::load(const string meshFile, const StringVec& varFiles,
62 const StringVec& varNames, int nBlocks)
63 {
64 numParts = nBlocks;
65 meshFmt = meshFile;
66
67 // initialize variables
68 StringVec::const_iterator fileIt = varFiles.begin();
69 StringVec::const_iterator nameIt = varNames.begin();
70 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
71 VarInfo vi;
72 vi.fileName = *fileIt;
73 vi.varName = *nameIt;
74 vi.valid = true;
75 variables.push_back(vi);
76 }
77
78 // Load the mesh files
79 if (!readMeshes()) {
80 cerr << "Reading the domain failed." << endl;
81 return false;
82 }
83
84 // Retrieve the mesh variables
85 convertMeshVariables();
86
87 // Load the variables
88 if (!readVariables()) {
89 cerr << "Reading the variables failed." << endl;
90 return false;
91 }
92
93 return true;
94 }
95
96 //
97 // Load only variables using provided mesh
98 //
99 bool MPDataSet::load(const MeshBlocks& m, const string siloFile,
100 const StringVec& varFiles, const StringVec& varNames)
101 {
102 siloMeshFile = siloFile;
103 meshBlocks = m;
104 numParts = meshBlocks.size();
105 keepMesh = true;
106
107 // initialize variables
108 StringVec::const_iterator fileIt = varFiles.begin();
109 StringVec::const_iterator nameIt = varNames.begin();
110 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
111 VarInfo vi;
112 vi.fileName = *fileIt;
113 vi.varName = *nameIt;
114 vi.valid = true;
115 variables.push_back(vi);
116 }
117
118 // Load the variables
119 if (!readVariables()) {
120 cerr << "Reading the variables failed." << endl;
121 return false;
122 }
123
124 return true;
125 }
126
127
128 //
129 //
130 //
131 bool MPDataSet::readMeshes()
132 {
133 bool ok = true;
134 char* str = new char[meshFmt.length()+10];
135 for (int idx=0; idx < numParts; idx++) {
136 MeshWithElements* meshPart = new MeshWithElements();
137 sprintf(str, meshFmt.c_str(), idx);
138 string meshfile = str;
139 if (meshPart->readFromNc(meshfile)) {
140 if (numParts > 1)
141 meshPart->handleGhostZones(idx);
142 meshBlocks.push_back(meshPart);
143 } else {
144 delete meshPart;
145 ok = false;
146 break;
147 }
148 }
149 delete[] str;
150 return ok;
151 }
152
153 void MPDataSet::convertMeshVariables()
154 {
155 const StringVec& varNames = meshBlocks[0]->getVarNames();
156 StringVec::const_iterator it;
157 for (it = varNames.begin(); it != varNames.end(); it++) {
158 VarInfo vi;
159 vi.varName = *it;
160 vi.valid = true;
161 DataParts parts;
162 // get all parts of current variable
163 for (int idx=0; idx < numParts; idx++) {
164 const IntVec& data = meshBlocks[idx]->getVarDataByName(*it);
165 DataVar* var = new DataVar(*it, data, meshBlocks[idx]);
166 parts.push_back(var);
167 }
168 // if we have more than one part, assemble variable
169 if (numParts > 1) {
170 assembleVariable(vi, parts);
171 for (size_t i=0; i<parts.size(); i++)
172 delete parts[i];
173 } else
174 vi.dataVar = parts.back();
175
176 meshVariables.push_back(vi);
177 }
178 }
179
180 //
181 //
182 //
183 bool MPDataSet::readVariables()
184 {
185 if (variables.size() == 0)
186 return true;
187
188 VarVector::iterator it;
189 bool ok = false;
190 for (it = variables.begin(); it != variables.end(); it++) {
191 bool curVarOk = true;
192 VarInfo& vi = (*it);
193 DataParts parts;
194 char* str = new char[vi.fileName.length()+10];
195 // read all parts of current variable
196 for (int idx=0; idx < numParts; idx++) {
197 sprintf(str, vi.fileName.c_str(), idx);
198 string dfile = str;
199 DataVar* var = new DataVar(vi.varName);
200 if (var->readFromNc(dfile))
201 parts.push_back(var);
202 else {
203 delete var;
204 curVarOk = false;
205 break;
206 }
207 }
208 delete[] str;
209 if (curVarOk) {
210 // at least one variable was read without problems
211 ok = true;
212
213 // if we have more than one part, assemble variable
214 if (numParts > 1) {
215 assembleVariable(vi, parts);
216 for (size_t i=0; i<parts.size(); i++)
217 delete parts[i];
218 } else
219 vi.dataVar = parts.back();
220 } else {
221 vi.dataVar = NULL;
222 vi.valid = false;
223 }
224 }
225 return ok;
226 }
227
228 //
229 //
230 //
231 void MPDataSet::assembleVariable(VarInfo& vi, const DataParts& parts)
232 {
233 DataVar* var = new DataVar(*parts[0]);
234 vi.dataVar = var;
235
236 DataParts::const_iterator it;
237 for (it = parts.begin()+1; it != parts.end(); it++) {
238 var->append(*(*it));
239 }
240 }
241
242 bool MPDataSet::saveAsSilo(string siloFile, bool useMultiMesh)
243 {
244 #if HAVE_SILO
245 if (numParts == 0)
246 return false;
247
248 DBfile* dbfile;
249 dbfile = DBCreate(siloFile.c_str(), DB_CLOBBER, DB_LOCAL, "Data", DB_PDB);
250 if (!dbfile) {
251 cerr << "Could not create Silo file." << endl;
252 return false;
253 }
254 DBMkdir(dbfile, MESH_VARS);
255
256 MeshBlocks::iterator meshIt;
257 VarVector::iterator viIt;
258 int idx = 0;
259 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, idx++) {
260 char str[64];
261 snprintf(str, 64, "/block%04d", idx);
262 string siloPath(str);
263 DBMkdir(dbfile, siloPath.c_str());
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 if (useMeshFile) {
372 string vpath = string(MESH_VARS)+varName;
373 DBPutMultivar(dbfile, vpath.c_str(), meshBlocks.size(), &varnames[0],
374 &vartypes[0], NULL);
375 } else {
376 DBPutMultivar(dbfile, varName.c_str(), meshBlocks.size(), &varnames[0],
377 &vartypes[0], NULL);
378 }
379 #endif
380 }
381
382 void MPDataSet::putSiloMultiTensor(DBfile* dbfile, const DataVar* var)
383 {
384 #if HAVE_SILO
385 string tensorDir = var->getName()+string("_comps/");
386 DBSetDir(dbfile, "/");
387 DBMkdir(dbfile, tensorDir.c_str());
388 int one = 1;
389 DBoptlist* optList = DBMakeOptlist(1);
390 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
391 vector<int> vartypes(meshBlocks.size(), DB_UCDVAR);
392 const IntVec& shape = var->getShape();
393 for (int i=0; i<shape[1]; i++) {
394 for (int j=0; j<shape[0]; j++) {
395 vector<string> tempstrings;
396 vector<char*> varnames;
397 char comp[255];
398 snprintf(comp, 255, "%s_comps/a_%d%d", var->getName().c_str(), i,j);
399 for (size_t idx = 0; idx < meshBlocks.size(); idx++) {
400 string siloPath;
401 char str[64];
402 snprintf(str, 64, "/block%04d", static_cast<int>(idx));
403 siloPath = str;
404 tempstrings.push_back(siloPath + string("/") + string(comp));
405 varnames.push_back((char*)tempstrings.back().c_str());
406 }
407 DBPutMultivar(dbfile, comp, meshBlocks.size(), &varnames[0],
408 &vartypes[0], optList);
409 }
410 }
411 DBFreeOptlist(optList);
412 #endif
413 }
414
415 } // namespace EscriptReader
416

  ViewVC Help
Powered by ViewVC 1.1.26