/[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 2183 - (show annotations)
Fri Dec 19 03:52:50 2008 UTC (11 years, 11 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
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