/[escript]/trunk/weipa/src/DataVar.cpp
ViewVC logotype

Contents of /trunk/weipa/src/DataVar.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2910 - (show annotations)
Wed Feb 3 03:22:31 2010 UTC (9 years, 2 months ago) by caltinay
Original Path: trunk/dataexporter/src/DataVar.cpp
File size: 19622 byte(s)
Added a libescriptreader.a target which does not depend on escript or finley
at runtime. This will be used for the VisIt plugin so is not built by default.

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

  ViewVC Help
Powered by ViewVC 1.1.26