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

  ViewVC Help
Powered by ViewVC 1.1.26