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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2880 - (show annotations)
Thu Jan 28 01:21:30 2010 UTC (9 years, 2 months ago) by caltinay
Original Path: trunk/dataexporter/src/DataVar.cpp
File size: 17184 byte(s)
Better resolution for Hex27 & Rec9 data in dataexporter (mantis #486).
Saved Silo data resolution is now equal to finley.SaveVTK generated data in
all cases.

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

  ViewVC Help
Powered by ViewVC 1.1.26