/[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 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 3 months ago) by jfenwick
Original Path: trunk/weipa/src/DataVar.cpp
File size: 19455 byte(s)
Merging dudley and scons updates from branches

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

  ViewVC Help
Powered by ViewVC 1.1.26