/[escript]/trunk/weipa/src/DataVar.cpp
ViewVC logotype

Contents of /trunk/weipa/src/DataVar.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2834 - (show annotations)
Thu Jan 7 06:06:56 2010 UTC (9 years, 3 months ago) by caltinay
Original Path: trunk/dataexporter/src/DataVar.cpp
File size: 15353 byte(s)
Reduced elements are now handled separate from the main elements within the
data exporter. This simplifies usage considerably.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 #include <escriptexport/DataVar.h>
15 #include <escriptexport/ElementData.h>
16 #include <escriptexport/FinleyMesh.h>
17 #include <escriptexport/NodeData.h>
18 #include <escript/Data.h>
19
20 #if USE_NETCDF
21 #include <netcdf.hh>
22 #endif
23
24 #if USE_SILO
25 #include <silo.h>
26 #endif
27
28 using namespace std;
29
30 namespace escriptexport {
31
32 enum {
33 NODE_CENTERED = 1,
34 ZONE_CENTERED = 2
35 };
36
37 //
38 // Constructor
39 //
40 DataVar::DataVar(const string& name) :
41 initialized(false), varName(name),
42 numSamples(0), rank(0), ptsPerSample(0), centering(0)
43 {
44 }
45
46 //
47 // Copy constructor
48 //
49 DataVar::DataVar(const DataVar& d) :
50 varName(d.varName), numSamples(d.numSamples),
51 rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),
52 funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)
53 {
54 if (numSamples > 0) {
55 CoordArray::const_iterator it;
56 for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
57 float* c = new float[numSamples];
58 copy(*it, (*it)+numSamples, c);
59 dataArray.push_back(c);
60 }
61 }
62 initialized = d.initialized;
63 }
64
65 //
66 // Destructor
67 //
68 DataVar::~DataVar()
69 {
70 cleanup();
71 }
72
73 //
74 //
75 //
76 void DataVar::cleanup()
77 {
78 CoordArray::iterator it;
79 for (it = dataArray.begin(); it != dataArray.end(); it++)
80 delete[] *it;
81 dataArray.clear();
82 shape.clear();
83 sampleID.clear();
84 numSamples = 0;
85 initialized = false;
86 }
87
88 //
89 //
90 //
91 bool DataVar::initFromEscript(escript::Data& escriptData, FinleyMesh_ptr mesh)
92 {
93 cleanup();
94
95 if (!escriptData.actsExpanded()) {
96 cerr << "WARNING: Only expanded data supported!" << endl;
97 return false;
98 }
99
100 rank = escriptData.getDataPointRank();
101 ptsPerSample = escriptData.getNumDataPointsPerSample();
102 shape = escriptData.getDataPointShape();
103 funcSpace = escriptData.getFunctionSpace().getTypeCode();
104 numSamples = escriptData.getNumSamples();
105
106 if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
107 centering = NODE_CENTERED;
108 } else {
109 centering = ZONE_CENTERED;
110 }
111
112 #ifdef _DEBUG
113 cout << varName << ":\t" << numSamples << " samples, "
114 << ptsPerSample << " pts/s, rank: " << rank << endl;
115 #endif
116
117 initialized = true;
118
119 if (numSamples == 0)
120 return true;
121
122 const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
123 sampleID.insert(sampleID.end(), numSamples, 0);
124 copy(iPtr, iPtr+numSamples, sampleID.begin());
125
126 size_t dimSize = 1;
127 if (rank > 0)
128 dimSize *= shape[0];
129 if (rank > 1)
130 dimSize *= shape[1];
131 if (rank > 2) {
132 cerr << "WARNING: Rank " << rank << " data is not supported!\n";
133 initialized = false;
134 }
135
136 if (initialized) {
137 size_t dataSize = dimSize * ptsPerSample;
138 float* tempData = new float[dataSize*numSamples];
139 float* destPtr = tempData;
140 for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
141 const escript::DataAbstract::ValueType::value_type* values =
142 escriptData.getSampleDataRO(sampleNo);
143 copy(values, values+dataSize, destPtr);
144 destPtr += dataSize;
145 }
146
147 const float* srcPtr = tempData;
148 for (int i=0; i < dimSize; i++, srcPtr++) {
149 float* c = averageData(srcPtr, dimSize);
150 dataArray.push_back(c);
151 }
152 delete[] tempData;
153
154 initialized = filterSamples(mesh);
155 }
156
157 return initialized;
158 }
159
160 //
161 // Initialise with mesh data
162 //
163 bool DataVar::initFromMesh(FinleyMesh_ptr mesh)
164 {
165 cleanup();
166
167 const IntVec& data = mesh->getVarDataByName(varName);
168 rank = 0;
169 ptsPerSample = 1;
170 numSamples = data.size();
171
172 if (numSamples > 0) {
173 float* c = new float[numSamples];
174 dataArray.push_back(c);
175 IntVec::const_iterator it;
176 for (it=data.begin(); it != data.end(); it++)
177 *c++ = static_cast<float>(*it);
178
179 if (varName.find("ContactElements_") != varName.npos) {
180 funcSpace = FINLEY_CONTACT_ELEMENTS_1;
181 centering = ZONE_CENTERED;
182 sampleID.insert(sampleID.end(),
183 mesh->getContactElements()->getIDs().begin(),
184 mesh->getContactElements()->getIDs().end());
185 } else if (varName.find("FaceElements_") != varName.npos) {
186 funcSpace = FINLEY_FACE_ELEMENTS;
187 centering = ZONE_CENTERED;
188 sampleID.insert(sampleID.end(),
189 mesh->getFaceElements()->getIDs().begin(),
190 mesh->getFaceElements()->getIDs().end());
191 } else if (varName.find("Elements_") != varName.npos) {
192 funcSpace = FINLEY_ELEMENTS;
193 centering = ZONE_CENTERED;
194 sampleID.insert(sampleID.end(),
195 mesh->getElements()->getIDs().begin(),
196 mesh->getElements()->getIDs().end());
197 } else if (varName.find("Nodes_") != varName.npos) {
198 funcSpace = FINLEY_NODES;
199 centering = NODE_CENTERED;
200 sampleID.insert(sampleID.end(),
201 mesh->getNodes()->getNodeIDs().begin(),
202 mesh->getNodes()->getNodeIDs().end());
203 } else {
204 return false;
205 }
206
207 NodeData_ptr nodes = mesh->getMeshForFinleyFS(funcSpace);
208 meshName = nodes->getName();
209 siloMeshName = nodes->getFullSiloName();
210 }
211 initialized = true;
212
213 return initialized;
214 }
215
216 //
217 // Reads variable data from NetCDF file
218 //
219 bool DataVar::initFromNetCDF(const string& filename, FinleyMesh_ptr mesh)
220 {
221 cleanup();
222
223 #if USE_NETCDF
224 NcError ncerr(NcError::silent_nonfatal);
225 NcFile* input = new NcFile(filename.c_str());
226 if (!input->is_valid()) {
227 cerr << "Could not open input file " << filename << "." << endl;
228 delete input;
229 return false;
230 }
231
232 NcDim* dim;
233 NcAtt* att;
234
235 att = input->get_att("type_id");
236 int typeID = att->as_int(0);
237 if (typeID != 2) {
238 cerr << "WARNING: Only expanded data supported!" << endl;
239 delete input;
240 return false;
241 }
242
243 att = input->get_att("rank");
244 rank = att->as_int(0);
245
246 dim = input->get_dim("num_data_points_per_sample");
247 ptsPerSample = dim->size();
248
249 att = input->get_att("function_space_type");
250 funcSpace = att->as_int(0);
251
252 if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
253 centering = NODE_CENTERED;
254 } else {
255 centering = ZONE_CENTERED;
256 }
257
258 dim = input->get_dim("num_samples");
259 numSamples = dim->size();
260
261 #ifdef _DEBUG
262 cout << varName << ":\t" << numSamples << " samples, "
263 << ptsPerSample << " pts/s, rank: " << rank << endl;
264 #endif
265
266 initialized = true;
267
268 // if there are no data samples we're done
269 if (numSamples == 0) {
270 delete input;
271 return true;
272 }
273
274 sampleID.insert(sampleID.end(), numSamples, 0);
275 NcVar* var = input->get_var("id");
276 var->get(&sampleID[0], numSamples);
277
278 size_t dimSize = 1;
279 vector<long> counts;
280
281 if (rank > 0) {
282 dim = input->get_dim("d0");
283 int d = dim->size();
284 shape.push_back(d);
285 counts.push_back(d);
286 dimSize *= d;
287 }
288 if (rank > 1) {
289 dim = input->get_dim("d1");
290 int d = dim->size();
291 shape.push_back(d);
292 counts.push_back(d);
293 dimSize *= d;
294 }
295 if (rank > 2) {
296 cerr << "WARNING: Rank " << rank << " data is not supported!\n";
297 initialized = false;
298 }
299
300 if (initialized) {
301 size_t dataSize = dimSize*numSamples*ptsPerSample;
302 counts.push_back(ptsPerSample);
303 counts.push_back(numSamples);
304 float* tempData = new float[dataSize];
305 NcVar* var = input->get_var("data");
306 var->get(tempData, &counts[0]);
307
308 const float* srcPtr = tempData;
309 for (int i=0; i < dimSize; i++, srcPtr++) {
310 float* c = averageData(srcPtr, dimSize);
311 dataArray.push_back(c);
312 }
313 delete[] tempData;
314
315 initialized = filterSamples(mesh);
316 }
317
318 delete input;
319 #endif // USE_NETCDF
320
321 return initialized;
322 }
323
324 //
325 // Returns true if the data values are nodal, false if they are zonal.
326 //
327 bool DataVar::isNodeCentered() const
328 {
329 return (centering == NODE_CENTERED);
330 }
331
332 //
333 // Returns a subset of the src array according to stride parameter.
334 // If samples consist of multiple values they are averaged beforehand.
335 // Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
336 // (z0,z1,...)
337 //
338 float* DataVar::averageData(const float* src, size_t stride)
339 {
340 float* res = new float[numSamples];
341
342 if (ptsPerSample == 1) {
343 float* dest = res;
344 for (int i=0; i<numSamples; i++, src+=stride)
345 *dest++ = *src;
346 } else {
347 float* dest = res;
348 for (int i=0; i<numSamples; i++) {
349 double tmpVal = 0.0;
350 for (int j=0; j<ptsPerSample; j++, src+=stride)
351 tmpVal += *src;
352 *dest++ = (float)(tmpVal / ptsPerSample);
353 }
354 }
355 return res;
356 }
357
358 //
359 // Filters and reorders the raw sample values according to the IDs provided
360 // in 'requiredIDs'. This is used to have data arrays ordered according to
361 // the underlying mesh (i.e. DataID[i]==MeshNodeID[i])
362 //
363 bool DataVar::filterSamples(FinleyMesh_ptr finleyMesh)
364 {
365 if (numSamples == 0)
366 return true;
367
368 IndexMap id2idxMap;
369 const IntVec* requiredIDs = NULL;
370
371 NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
372 if (nodes == NULL)
373 return false;
374
375 int requiredNumSamples = 0;
376
377 if (centering == NODE_CENTERED) {
378 id2idxMap = nodes->getIndexMap();
379 requiredIDs = &nodes->getNodeIDs();
380 requiredNumSamples = nodes->getNumNodes();
381 } else {
382 ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
383 if (cells == NULL)
384 return false;
385
386 id2idxMap = cells->getIndexMap();
387 requiredIDs = &cells->getIDs();
388 requiredNumSamples = cells->getNumElements();
389 }
390
391 if (requiredNumSamples > numSamples) {
392 cerr << "ERROR: " << varName << " has " << numSamples
393 << " instead of " << requiredNumSamples << " samples!" << endl;
394 return false;
395 }
396
397 meshName = nodes->getName();
398 siloMeshName = nodes->getFullSiloName();
399
400 IndexMap sampleID2idx = buildIndexMap();
401 numSamples = requiredNumSamples;
402
403 // now filter the data
404 for (size_t i=0; i < dataArray.size(); i++) {
405 float* c = new float[numSamples];
406 const float* src = dataArray[i];
407 IntVec::const_iterator idIt = requiredIDs->begin();
408 for (; idIt != requiredIDs->end(); idIt++) {
409 size_t srcIdx = sampleID2idx.find(*idIt)->second;
410 size_t destIdx = id2idxMap.find(*idIt)->second;
411 c[destIdx] = src[srcIdx];
412 }
413 delete[] dataArray[i];
414 dataArray[i] = c;
415 }
416 return true;
417 }
418
419 ///////////////////////////////
420 // SILO related methods follow
421 ///////////////////////////////
422
423 //
424 // If the data is tensor data then the components of the tensor are stored
425 // separately in the Silo file. This method then returns a string that
426 // contains the proper Silo expression to put the tensor together again.
427 // For non-tensor data this method returns an empty string.
428 //
429 string DataVar::getTensorDef() const
430 {
431 if (rank < 2 || !initialized)
432 return string();
433
434 /// Format string for Silo 2x2 tensor
435 const string tensor2DefFmt =
436 "{{ <%sa_00>, <%sa_01> },"
437 " { <%sa_10>, <%sa_11> }}";
438
439 /// Format string for Silo 3x3 tensor
440 const string tensor3DefFmt =
441 "{{ <%sa_00>, <%sa_01>, <%sa_02> },"
442 " { <%sa_10>, <%sa_11>, <%sa_12> },"
443 " { <%sa_20>, <%sa_21>, <%sa_22> }}";
444
445 string tensorDef;
446 string tensorDir = varName+string("_comps/");
447 if (shape[1] == 3) {
448 char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
449 sprintf(tDef, tensor3DefFmt.c_str(),
450 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
451 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
452 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
453 tensorDef = tDef;
454 delete[] tDef;
455 } else {
456 char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
457 sprintf(tDef, tensor2DefFmt.c_str(),
458 tensorDir.c_str(), tensorDir.c_str(),
459 tensorDir.c_str(), tensorDir.c_str(),
460 tensorDir.c_str(), tensorDir.c_str());
461 tensorDef = tDef;
462 delete[] tDef;
463 }
464 return tensorDef;
465 }
466
467 //
468 // Writes the data to given Silo file under the virtual path provided.
469 // The corresponding mesh must have been written already and made known
470 // to this variable by a call to setMesh().
471 //
472 bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)
473 {
474 #if USE_SILO
475 if (!initialized)
476 return false;
477
478 if (numSamples == 0)
479 return true;
480
481 int ret;
482
483 if (siloPath != "") {
484 ret = DBSetDir(dbfile, siloPath.c_str());
485 if (ret != 0)
486 return false;
487 }
488
489 char* siloMesh = const_cast<char*>(siloMeshName.c_str());
490 int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
491
492 if (rank == 0) {
493 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
494 numSamples, NULL, 0, DB_FLOAT, dcenter, NULL);
495 }
496 else if (rank == 1) {
497 const string comps[3] = {
498 varName+string("_x"), varName+string("_y"), varName+string("_z")
499 };
500 const char* varnames[3] = {
501 comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
502 };
503
504 ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
505 (char**)varnames, &dataArray[0], numSamples, NULL,
506 0, DB_FLOAT, dcenter, NULL);
507 }
508 else {
509 string tensorDir = varName+string("_comps/");
510 ret = DBMkdir(dbfile, tensorDir.c_str());
511 if (ret == 0) {
512 int one = 1;
513 DBoptlist* optList = DBMakeOptlist(1);
514 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
515
516 for (int i=0; i<shape[1]; i++) {
517 for (int j=0; j<shape[0]; j++) {
518 ostringstream varname;
519 varname << tensorDir << "a_" << i << j;
520 ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
521 dataArray[i*shape[0]+j], numSamples,
522 NULL, 0, DB_FLOAT, dcenter, optList);
523 if (ret != 0) break;
524 }
525 if (ret != 0) break;
526 }
527 DBFreeOptlist(optList);
528 } // ret==0
529 } // rank
530
531 DBSetDir(dbfile, "/");
532 return (ret == 0);
533
534 #else // !USE_SILO
535 return false;
536 #endif
537 }
538
539 } // namespace escriptexport
540

  ViewVC Help
Powered by ViewVC 1.1.26