/[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 3128 - (show annotations)
Wed Sep 1 03:54:09 2010 UTC (8 years, 4 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 20635 byte(s)
-Added mesh axis label/unit and variable unit support to weipa's Silo writer
-Created EscriptDataset python wrapper
-Some code cleanup and a ghost zone fix for the VisIt interface

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

  ViewVC Help
Powered by ViewVC 1.1.26