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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3357 - (show annotations)
Wed Nov 17 06:21:37 2010 UTC (8 years, 6 months ago) by caltinay
File size: 19471 byte(s)
Fixed a crash when loading data from netCDF files.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 <weipa/DataVar.h>
15 #include <weipa/DomainChunk.h>
16 #include <weipa/ElementData.h>
17 #include <weipa/NodeData.h>
18 #ifndef VISIT_PLUGIN
19 #include <escript/Data.h>
20 #endif
21
22 #if USE_NETCDF
23 #include <netcdfcpp.h>
24 #endif
25
26 #if USE_SILO
27 #include <silo.h>
28 #endif
29
30 #include <numeric> // for accumulate
31 #include <iostream> // for cerr
32 #include <stdio.h>
33
34 using namespace std;
35
36 namespace weipa {
37
38 //
39 // Constructor
40 //
41 DataVar::DataVar(const string& name) :
42 initialized(false), varName(name),
43 numSamples(0), rank(0), ptsPerSample(0)
44 {
45 }
46
47 //
48 // Copy constructor
49 //
50 DataVar::DataVar(const DataVar& d) :
51 initialized(d.initialized), domain(d.domain),
52 varName(d.varName), numSamples(d.numSamples),
53 rank(d.rank), ptsPerSample(d.ptsPerSample), funcSpace(d.funcSpace),
54 centering(d.centering), shape(d.shape), sampleID(d.sampleID)
55 {
56 if (numSamples > 0) {
57 CoordArray::const_iterator it;
58 for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
59 float* c = new float[numSamples];
60 copy(*it, (*it)+numSamples, c);
61 dataArray.push_back(c);
62 }
63 }
64 }
65
66 //
67 // Destructor
68 //
69 DataVar::~DataVar()
70 {
71 cleanup();
72 }
73
74 //
75 //
76 //
77 void DataVar::cleanup()
78 {
79 CoordArray::iterator it;
80 for (it = dataArray.begin(); it != dataArray.end(); it++)
81 delete[] *it;
82 dataArray.clear();
83 shape.clear();
84 sampleID.clear();
85 numSamples = 0;
86 initialized = false;
87 }
88
89 //
90 //
91 //
92 bool DataVar::initFromEscript(escript::Data& escriptData, const_DomainChunk_ptr dom)
93 {
94 #ifndef VISIT_PLUGIN
95 cleanup();
96
97 if (!escriptData.isConstant() && !escriptData.actsExpanded()) {
98 cerr << "WARNING: Weipa only supports constant & expanded data, "
99 << "not initializing " << varName << endl;
100 return false;
101 }
102
103 domain = dom;
104 rank = escriptData.getDataPointRank();
105 ptsPerSample = escriptData.getNumDataPointsPerSample();
106 shape = escriptData.getDataPointShape();
107 funcSpace = escriptData.getFunctionSpace().getTypeCode();
108 numSamples = escriptData.getNumSamples();
109 centering = domain->getCenteringForFunctionSpace(funcSpace);
110
111 #ifdef _DEBUG
112 cout << varName << ":\t" << numSamples << " samples, "
113 << ptsPerSample << " pts/s, rank: " << rank << endl;
114 #endif
115
116 NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
117 if (nodes == NULL)
118 return false;
119
120 meshName = nodes->getName();
121 siloMeshName = nodes->getFullSiloName();
122 initialized = true;
123
124 // no samples? Nothing more to do.
125 if (numSamples == 0)
126 return true;
127
128 const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
129 sampleID.insert(sampleID.end(), numSamples, 0);
130 copy(iPtr, iPtr+numSamples, sampleID.begin());
131
132 size_t dimSize = 1;
133 if (rank > 0)
134 dimSize *= shape[0];
135 if (rank > 1)
136 dimSize *= shape[1];
137 if (rank > 2) {
138 cerr << "WARNING: Rank " << rank << " data is not supported!\n";
139 initialized = false;
140 }
141
142 if (initialized) {
143 size_t dataSize = dimSize * ptsPerSample;
144 float* tempData = new float[dataSize*numSamples];
145 float* destPtr = tempData;
146 if (escriptData.isConstant()) {
147 const escript::DataAbstract::ValueType::value_type* values =
148 escriptData.getSampleDataRO(0);
149 for (int pointNo=0; pointNo<numSamples*ptsPerSample; pointNo++) {
150 copy(values, values+dimSize, destPtr);
151 destPtr += dimSize;
152 }
153 } else {
154 for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
155 const escript::DataAbstract::ValueType::value_type* values =
156 escriptData.getSampleDataRO(sampleNo);
157 copy(values, values+dataSize, destPtr);
158 destPtr += dataSize;
159 }
160 }
161
162 const float* srcPtr = tempData;
163 for (int i=0; i < dimSize; i++, srcPtr++) {
164 float* c = averageData(srcPtr, dimSize);
165 dataArray.push_back(c);
166 }
167 delete[] tempData;
168
169 initialized = reorderSamples();
170 }
171
172 return initialized;
173
174 #else // VISIT_PLUGIN
175 return false;
176 #endif
177 }
178
179 //
180 // Initialise with domain data
181 //
182 bool DataVar::initFromMeshData(const_DomainChunk_ptr dom, const IntVec& data,
183 int fsCode, Centering c, NodeData_ptr nodes, const IntVec& id)
184 {
185 cleanup();
186
187 domain = dom;
188 rank = 0;
189 ptsPerSample = 1;
190 centering = c;
191 sampleID = id;
192 meshName = nodes->getName();
193 siloMeshName = nodes->getFullSiloName();
194 numSamples = data.size();
195
196 if (numSamples > 0) {
197 float* c = new float[numSamples];
198 dataArray.push_back(c);
199 IntVec::const_iterator it;
200 for (it=data.begin(); it != data.end(); it++)
201 *c++ = static_cast<float>(*it);
202 }
203 initialized = true;
204
205 return initialized;
206 }
207
208 //
209 // Reads variable data from dump file
210 //
211 bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
212 {
213 cleanup();
214
215 #if USE_NETCDF
216 NcError ncerr(NcError::silent_nonfatal);
217 NcFile* input = new NcFile(filename.c_str());
218 if (!input->is_valid()) {
219 cerr << "Could not open input file " << filename << "." << endl;
220 delete input;
221 return false;
222 }
223
224 NcDim* dim;
225 NcAtt* att;
226
227 att = input->get_att("type_id");
228 int typeID = att->as_int(0);
229 if (typeID != 2) {
230 cerr << "WARNING: Only expanded data supported!" << endl;
231 delete input;
232 return false;
233 }
234
235 att = input->get_att("rank");
236 rank = att->as_int(0);
237
238 dim = input->get_dim("num_data_points_per_sample");
239 ptsPerSample = dim->size();
240
241 att = input->get_att("function_space_type");
242 funcSpace = att->as_int(0);
243
244 centering = dom->getCenteringForFunctionSpace(funcSpace);
245
246 dim = input->get_dim("num_samples");
247 numSamples = dim->size();
248
249 #ifdef _DEBUG
250 cout << varName << ":\t" << numSamples << " samples, "
251 << ptsPerSample << " pts/s, rank: " << rank << endl;
252 #endif
253
254 domain = dom;
255 NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
256 if (nodes == NULL) {
257 delete input;
258 return false;
259 }
260
261 meshName = nodes->getName();
262 siloMeshName = nodes->getFullSiloName();
263 initialized = true;
264
265 size_t dimSize = 1;
266 vector<long> counts;
267
268 if (rank > 0) {
269 dim = input->get_dim("d0");
270 int d = dim->size();
271 shape.push_back(d);
272 counts.push_back(d);
273 dimSize *= d;
274 }
275 if (rank > 1) {
276 dim = input->get_dim("d1");
277 int d = dim->size();
278 shape.push_back(d);
279 counts.push_back(d);
280 dimSize *= d;
281 }
282 if (rank > 2) {
283 cerr << "WARNING: Rank " << rank << " data is not supported!\n";
284 initialized = false;
285 }
286
287 if (initialized && numSamples > 0) {
288 sampleID.insert(sampleID.end(), numSamples, 0);
289 NcVar* var = input->get_var("id");
290 var->get(&sampleID[0], numSamples);
291
292 size_t dataSize = dimSize*numSamples*ptsPerSample;
293 counts.push_back(ptsPerSample);
294 counts.push_back(numSamples);
295 float* tempData = new float[dataSize];
296 var = input->get_var("data");
297 var->get(tempData, &counts[0]);
298
299 const float* srcPtr = tempData;
300 for (int i=0; i < dimSize; i++, srcPtr++) {
301 float* c = averageData(srcPtr, dimSize);
302 dataArray.push_back(c);
303 }
304 delete[] tempData;
305
306 initialized = reorderSamples();
307 }
308
309 delete input;
310 #endif // USE_NETCDF
311
312 return initialized;
313 }
314
315 //
316 // Returns true if the data values are nodal, false if they are zonal.
317 //
318 bool DataVar::isNodeCentered() const
319 {
320 return (centering == NODE_CENTERED);
321 }
322
323 //
324 // Returns a subset of the src array according to stride parameter.
325 // If samples consist of multiple values they are averaged beforehand.
326 // Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
327 // (z0,z1,...)
328 //
329 float* DataVar::averageData(const float* src, size_t stride)
330 {
331 float* res;
332
333 if (ptsPerSample == 1) {
334 res = new float[numSamples];
335 float* dest = res;
336 for (int i=0; i<numSamples; i++, src+=stride)
337 *dest++ = *src;
338 } else {
339 ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
340 int cellFactor = cells->getElementFactor();
341 res = new float[cellFactor * numSamples];
342 float* dest = res;
343 QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
344 if (!qmi.mask.empty()) {
345 const float* tmpSrc = src;
346 for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
347 for (int l=0; l<cellFactor; l++) {
348 double tmpVal = 0.0;
349 for (int j=0; j<ptsPerSample; j++) {
350 if (qmi.mask[l][j] != 0) {
351 tmpVal += *(tmpSrc+stride*j);
352 }
353 }
354 *dest++ = (float)(tmpVal / qmi.factor[l]);
355 }
356 }
357 } else {
358 for (int i=0; i<numSamples; i++) {
359 double tmpVal = 0.0;
360 for (int j=0; j<ptsPerSample; j++, src+=stride) {
361 tmpVal += *src;
362 }
363 tmpVal /= ptsPerSample;
364 for (int l=0; l<cellFactor; l++) {
365 *dest++ = static_cast<float>(tmpVal);
366 }
367 }
368 }
369 }
370 return res;
371 }
372
373 //
374 // Filters and reorders the raw sample values according to the node/element
375 // IDs. This is used to have data arrays ordered according to the underlying
376 // mesh (i.e. DataID[i]==MeshNodeID[i])
377 //
378 bool DataVar::reorderSamples()
379 {
380 if (numSamples == 0)
381 return true;
382
383 const IntVec* requiredIDs = NULL;
384 int requiredNumSamples = 0;
385 int cellFactor = 1;
386
387 if (centering == NODE_CENTERED) {
388 NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
389 requiredIDs = &nodes->getNodeIDs();
390 requiredNumSamples = nodes->getNumNodes();
391 } else {
392 ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
393 if (cells == NULL)
394 return false;
395
396 requiredIDs = &cells->getIDs();
397 requiredNumSamples = cells->getNumElements();
398 cellFactor = cells->getElementFactor();
399 if (cellFactor > 1) {
400 numSamples *= cellFactor;
401 // update sample IDs
402 IntVec newSampleID(numSamples);
403 IntVec::const_iterator idIt = sampleID.begin();
404 IntVec::iterator newIDit = newSampleID.begin();
405 for (; idIt != sampleID.end(); idIt++, newIDit+=cellFactor) {
406 fill(newIDit, newIDit+cellFactor, *idIt);
407 }
408 sampleID.swap(newSampleID);
409 }
410 }
411
412 if (requiredNumSamples > numSamples) {
413 cerr << "ERROR: " << varName << " has " << numSamples
414 << " instead of " << requiredNumSamples << " samples!" << endl;
415 return false;
416 }
417
418 IndexMap sampleID2idx = buildIndexMap();
419 numSamples = requiredNumSamples;
420
421 // now filter the data
422 for (size_t i=0; i < dataArray.size(); i++) {
423 float* c = new float[numSamples];
424 const float* src = dataArray[i];
425 IntVec::const_iterator idIt = requiredIDs->begin();
426 size_t destIdx = 0;
427 for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
428 size_t srcIdx = sampleID2idx.find(*idIt)->second;
429 copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
430 }
431 delete[] dataArray[i];
432 dataArray[i] = c;
433 }
434
435 // sample IDs now = mesh node/element IDs
436 sampleID = *requiredIDs;
437
438 return true;
439 }
440
441 //
442 //
443 //
444 int DataVar::getNumberOfComponents() const
445 {
446 return (rank == 0 ? 1 : accumulate(shape.begin(), shape.end(), 0));
447 }
448
449 //
450 //
451 //
452 float* DataVar::getDataFlat() const
453 {
454 int totalSize = numSamples * getNumberOfComponents();
455 float* res = new float[totalSize];
456 if (rank == 0) {
457 copy(dataArray[0], dataArray[0]+numSamples, res);
458 } else if (rank == 1) {
459 float *dest = res;
460 for (size_t c=0; c<numSamples; c++) {
461 for (size_t i=0; i<shape[0]; i++) {
462 *dest++ = dataArray[i][c];
463 }
464 }
465 } else if (rank == 2) {
466 float *dest = res;
467 for (size_t c=0; c<numSamples; c++) {
468 for (int i=0; i<shape[1]; i++) {
469 for (int j=0; j<shape[0]; j++) {
470 *dest++ = dataArray[i*shape[0]+j][c];
471 }
472 }
473 }
474 }
475
476 return res;
477 }
478
479 //
480 //
481 //
482 void DataVar::sampleToStream(ostream& os, int index)
483 {
484 if (rank == 0) {
485 os << dataArray[0][index];
486 } else if (rank == 1) {
487 if (shape[0] < 3)
488 os << dataArray[0][index] << " " << dataArray[1][index]
489 << " " << 0.;
490 else
491 os << dataArray[0][index] << " " << dataArray[1][index]
492 << " " << dataArray[2][index];
493 } else if (rank == 2) {
494 if (shape[1] < 3) {
495 os << dataArray[0][index] << " " << dataArray[1][index]
496 << " " << 0. << " ";
497 os << dataArray[2][index] << " " << dataArray[3][index]
498 << " " << 0. << " ";
499 os << 0. << " " << 0. << " " << 0.;
500 } else {
501 os << dataArray[0][index] << " " << dataArray[1][index]
502 << " " << dataArray[2][index] << " ";
503 os << dataArray[3][index] << " " << dataArray[4][index]
504 << " " << dataArray[5][index] << " ";
505 os << dataArray[6][index] << " " << dataArray[7][index]
506 << " " << dataArray[8][index];
507 }
508 }
509 os << endl;
510 }
511
512 //
513 //
514 //
515 void DataVar::writeToVTK(ostream& os, int ownIndex)
516 {
517 if (numSamples == 0)
518 return;
519
520 if (isNodeCentered()) {
521 // data was reordered in reorderSamples() but for VTK we write the
522 // original node mesh and thus need the original ordering...
523 const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
524 const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
525 const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
526 int firstId = nodeDist[ownIndex];
527 int lastId = nodeDist[ownIndex+1];
528 IndexMap sampleID2idx = buildIndexMap();
529 for (int i=0; i<nodeGNI.size(); i++) {
530 if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
531 int idx = sampleID2idx[requiredIDs[i]];
532 sampleToStream(os, idx);
533 }
534 }
535 } else {
536 // cell data: ghost cells have been removed so do not write ghost
537 // samples (which are the last elements in the arrays)
538 int toWrite = domain->getElementsByName(meshName)->getNumElements();
539 for (int i=0; i<toWrite; i++) {
540 sampleToStream(os, i);
541 }
542 }
543 }
544
545 ///////////////////////////////
546 // SILO related methods follow
547 ///////////////////////////////
548
549 //
550 // If the data is tensor data then the components of the tensor are stored
551 // separately in the Silo file. This method then returns a string that
552 // contains the proper Silo expression to put the tensor together again.
553 // For non-tensor data this method returns an empty string.
554 //
555 string DataVar::getTensorDef() const
556 {
557 if (rank < 2 || !initialized)
558 return string();
559
560 /// Format string for Silo 2x2 tensor
561 const string tensor2DefFmt =
562 "{{ <%sa_00>, <%sa_01> },"
563 " { <%sa_10>, <%sa_11> }}";
564
565 /// Format string for Silo 3x3 tensor
566 const string tensor3DefFmt =
567 "{{ <%sa_00>, <%sa_01>, <%sa_02> },"
568 " { <%sa_10>, <%sa_11>, <%sa_12> },"
569 " { <%sa_20>, <%sa_21>, <%sa_22> }}";
570
571 string tensorDef;
572 string tensorDir = varName+string("_comps/");
573 if (shape[1] == 3) {
574 char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
575 sprintf(tDef, tensor3DefFmt.c_str(),
576 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
577 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
578 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
579 tensorDef = tDef;
580 delete[] tDef;
581 } else {
582 char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
583 sprintf(tDef, tensor2DefFmt.c_str(),
584 tensorDir.c_str(), tensorDir.c_str(),
585 tensorDir.c_str(), tensorDir.c_str(),
586 tensorDir.c_str(), tensorDir.c_str());
587 tensorDef = tDef;
588 delete[] tDef;
589 }
590 return tensorDef;
591 }
592
593 //
594 // Writes the data to given Silo file under the virtual path provided.
595 // The corresponding mesh must have been written already.
596 //
597 bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
598 const string& units)
599 {
600 #if USE_SILO
601 if (!initialized)
602 return false;
603
604 if (numSamples == 0)
605 return true;
606
607 int ret;
608
609 if (siloPath != "") {
610 ret = DBSetDir(dbfile, siloPath.c_str());
611 if (ret != 0)
612 return false;
613 }
614
615 char* siloMesh = const_cast<char*>(siloMeshName.c_str());
616 int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
617 DBoptlist* optList = DBMakeOptlist(2);
618 if (units.length()>0) {
619 DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
620 }
621
622 if (rank == 0) {
623 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
624 numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
625 }
626 else if (rank == 1) {
627 const string comps[3] = {
628 varName+string("_x"), varName+string("_y"), varName+string("_z")
629 };
630 const char* varnames[3] = {
631 comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
632 };
633
634 ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
635 (char**)varnames, &dataArray[0], numSamples, NULL,
636 0, DB_FLOAT, dcenter, optList);
637 }
638 else {
639 string tensorDir = varName+string("_comps/");
640 ret = DBMkdir(dbfile, tensorDir.c_str());
641 if (ret == 0) {
642 int one = 1;
643 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
644
645 for (int i=0; i<shape[1]; i++) {
646 for (int j=0; j<shape[0]; j++) {
647 ostringstream varname;
648 varname << tensorDir << "a_" << i << j;
649 ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
650 dataArray[i*shape[0]+j], numSamples,
651 NULL, 0, DB_FLOAT, dcenter, optList);
652 if (ret != 0) break;
653 }
654 if (ret != 0) break;
655 }
656 } // ret==0
657 } // rank
658
659 DBFreeOptlist(optList);
660 DBSetDir(dbfile, "/");
661 return (ret == 0);
662
663 #else // !USE_SILO
664 return false;
665 #endif
666 }
667
668 } // namespace weipa
669

  ViewVC Help
Powered by ViewVC 1.1.26