/[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 4763 - (show annotations)
Tue Mar 18 05:20:23 2014 UTC (3 years, 9 months ago) by jfenwick
File size: 20527 byte(s)
Replaced references to getSampleDataRO(0) with getDataRO()
This was to deal with the case where processes with no data where trying to compute the offset of sample 0 when all that was really needed was a null.

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

  ViewVC Help
Powered by ViewVC 1.1.26