/[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 3185 - (show annotations)
Thu Sep 16 00:30:09 2010 UTC (8 years, 4 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 19423 byte(s)
Weipa now supports initialization from DataConstant instances.

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

  ViewVC Help
Powered by ViewVC 1.1.26