/[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 3146 - (show annotations)
Fri Sep 3 01:29:21 2010 UTC (8 years, 4 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 18936 byte(s)
Properly initialise data centering variable.

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

  ViewVC Help
Powered by ViewVC 1.1.26