/[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 2194 - (show annotations)
Tue Jan 6 01:23:37 2009 UTC (13 years, 7 months ago) by caltinay
File size: 12141 byte(s)
escript2silo: Allow different meshes for different timesteps.
escriptreader: Fixed segfaults for some error conditions.


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

  ViewVC Help
Powered by ViewVC 1.1.26