/[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 3398 - (show annotations)
Mon Dec 6 00:02:33 2010 UTC (8 years, 1 month ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 20201 byte(s)
Fixed the quad mask in a very specific case and changed saveVTK to write a
specific value for unused nodes instead of the first one which breaks unit
tests depending on the number of ranks used.
Still need to rebaseline some VTK files so the tests pass.

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 // index is -1 for dummy samples, i.e. if writing the full mesh but
485 // only a reduced number of samples is required
486 if (rank == 0) {
487 if (index < 0) {
488 os << 0.;
489 } else {
490 os << dataArray[0][index];
491 }
492 } else if (rank == 1) {
493 if (index < 0) {
494 os << 0. << " " << 0. << " " << 0.;
495 } else if (shape[0] < 3) {
496 os << dataArray[0][index] << " " << dataArray[1][index]
497 << " " << 0.;
498 } else {
499 os << dataArray[0][index] << " " << dataArray[1][index]
500 << " " << dataArray[2][index];
501 }
502 } else if (rank == 2) {
503 if (index < 0) {
504 os << 0. << " " << 0. << " " << 0. << " ";
505 os << 0. << " " << 0. << " " << 0. << " ";
506 os << 0. << " " << 0. << " " << 0.;
507 } else if (shape[1] < 3) {
508 os << dataArray[0][index] << " " << dataArray[1][index]
509 << " " << 0. << " ";
510 os << dataArray[2][index] << " " << dataArray[3][index]
511 << " " << 0. << " ";
512 os << 0. << " " << 0. << " " << 0.;
513 } else {
514 os << dataArray[0][index] << " " << dataArray[1][index]
515 << " " << dataArray[2][index] << " ";
516 os << dataArray[3][index] << " " << dataArray[4][index]
517 << " " << dataArray[5][index] << " ";
518 os << dataArray[6][index] << " " << dataArray[7][index]
519 << " " << dataArray[8][index];
520 }
521 }
522 os << endl;
523 }
524
525 //
526 //
527 //
528 void DataVar::writeToVTK(ostream& os, int ownIndex)
529 {
530 if (numSamples == 0)
531 return;
532
533 if (isNodeCentered()) {
534 // data was reordered in reorderSamples() but for VTK we write the
535 // original node mesh and thus need the original ordering.
536 // Note, that this also means we may not have samples for all nodes
537 // in which case we set idx to -1 and write a dummy sample
538 const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
539 const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
540 const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
541 int firstId = nodeDist[ownIndex];
542 int lastId = nodeDist[ownIndex+1];
543 IndexMap sampleID2idx = buildIndexMap();
544 for (int i=0; i<nodeGNI.size(); i++) {
545 if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
546 IndexMap::const_iterator it = sampleID2idx.find(requiredIDs[i]);
547 int idx = (it==sampleID2idx.end() ? -1 : (int)it->second);
548 sampleToStream(os, idx);
549 }
550 }
551 } else {
552 // cell data: ghost cells have been removed so do not write ghost
553 // samples (which are the last elements in the arrays)
554 int toWrite = domain->getElementsByName(meshName)->getNumElements();
555 for (int i=0; i<toWrite; i++) {
556 sampleToStream(os, i);
557 }
558 }
559 }
560
561 ///////////////////////////////
562 // SILO related methods follow
563 ///////////////////////////////
564
565 //
566 // If the data is tensor data then the components of the tensor are stored
567 // separately in the Silo file. This method then returns a string that
568 // contains the proper Silo expression to put the tensor together again.
569 // For non-tensor data this method returns an empty string.
570 //
571 string DataVar::getTensorDef() const
572 {
573 if (rank < 2 || !initialized)
574 return string();
575
576 /// Format string for Silo 2x2 tensor
577 const string tensor2DefFmt =
578 "{{ <%sa_00>, <%sa_01> },"
579 " { <%sa_10>, <%sa_11> }}";
580
581 /// Format string for Silo 3x3 tensor
582 const string tensor3DefFmt =
583 "{{ <%sa_00>, <%sa_01>, <%sa_02> },"
584 " { <%sa_10>, <%sa_11>, <%sa_12> },"
585 " { <%sa_20>, <%sa_21>, <%sa_22> }}";
586
587 string tensorDef;
588 string tensorDir = varName+string("_comps/");
589 if (shape[1] == 3) {
590 char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
591 sprintf(tDef, tensor3DefFmt.c_str(),
592 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
593 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
594 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
595 tensorDef = tDef;
596 delete[] tDef;
597 } else {
598 char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
599 sprintf(tDef, tensor2DefFmt.c_str(),
600 tensorDir.c_str(), tensorDir.c_str(),
601 tensorDir.c_str(), tensorDir.c_str(),
602 tensorDir.c_str(), tensorDir.c_str());
603 tensorDef = tDef;
604 delete[] tDef;
605 }
606 return tensorDef;
607 }
608
609 //
610 // Writes the data to given Silo file under the virtual path provided.
611 // The corresponding mesh must have been written already.
612 //
613 bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
614 const string& units)
615 {
616 #if USE_SILO
617 if (!initialized)
618 return false;
619
620 if (numSamples == 0)
621 return true;
622
623 int ret;
624
625 if (siloPath != "") {
626 ret = DBSetDir(dbfile, siloPath.c_str());
627 if (ret != 0)
628 return false;
629 }
630
631 char* siloMesh = const_cast<char*>(siloMeshName.c_str());
632 int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
633 DBoptlist* optList = DBMakeOptlist(2);
634 if (units.length()>0) {
635 DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
636 }
637
638 if (rank == 0) {
639 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
640 numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
641 }
642 else if (rank == 1) {
643 const string comps[3] = {
644 varName+string("_x"), varName+string("_y"), varName+string("_z")
645 };
646 const char* varnames[3] = {
647 comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
648 };
649
650 ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
651 (char**)varnames, &dataArray[0], numSamples, NULL,
652 0, DB_FLOAT, dcenter, optList);
653 }
654 else {
655 string tensorDir = varName+string("_comps/");
656 ret = DBMkdir(dbfile, tensorDir.c_str());
657 if (ret == 0) {
658 int one = 1;
659 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
660
661 for (int i=0; i<shape[1]; i++) {
662 for (int j=0; j<shape[0]; j++) {
663 ostringstream varname;
664 varname << tensorDir << "a_" << i << j;
665 ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
666 dataArray[i*shape[0]+j], numSamples,
667 NULL, 0, DB_FLOAT, dcenter, optList);
668 if (ret != 0) break;
669 }
670 if (ret != 0) break;
671 }
672 } // ret==0
673 } // rank
674
675 DBFreeOptlist(optList);
676 DBSetDir(dbfile, "/");
677 return (ret == 0);
678
679 #else // !USE_SILO
680 return false;
681 #endif
682 }
683
684 } // namespace weipa
685

  ViewVC Help
Powered by ViewVC 1.1.26