/[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 3143 - (show annotations)
Fri Sep 3 00:31:55 2010 UTC (8 years, 4 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 18917 byte(s)
-Moved finley specifics into weipa subclasses
-EscriptDataset is now the only weipa class exported in Windows
-Some weipa code cleanup


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

  ViewVC Help
Powered by ViewVC 1.1.26