/[escript]/branches/split/weipa/src/DataVar.cpp
ViewVC logotype

Contents of /branches/split/weipa/src/DataVar.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26