/[escript]/branches/split/weipa/src/DataVar.cpp
ViewVC logotype

Annotation of /branches/split/weipa/src/DataVar.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2880 - (hide annotations)
Thu Jan 28 01:21:30 2010 UTC (9 years, 1 month 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 caltinay 2183
2     /*******************************************************
3     *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 caltinay 2183 * 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 caltinay 2810 #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 caltinay 2183 #include <netcdf.hh>
22 caltinay 2810 #endif
23    
24     #if USE_SILO
25 caltinay 2183 #include <silo.h>
26     #endif
27    
28     using namespace std;
29    
30 caltinay 2810 namespace escriptexport {
31 caltinay 2187
32 caltinay 2183 enum {
33     NODE_CENTERED = 1,
34     ZONE_CENTERED = 2
35     };
36    
37     //
38     // Constructor
39     //
40     DataVar::DataVar(const string& name) :
41 caltinay 2810 initialized(false), varName(name),
42     numSamples(0), rank(0), ptsPerSample(0), centering(0)
43 caltinay 2183 {
44     }
45    
46     //
47     // Copy constructor
48     //
49     DataVar::DataVar(const DataVar& d) :
50 caltinay 2880 initialized(d.initialized), finleyMesh(d.finleyMesh),
51 caltinay 2183 varName(d.varName), numSamples(d.numSamples),
52     rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),
53 caltinay 2810 funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)
54 caltinay 2183 {
55 caltinay 2810 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 caltinay 2183 }
63     }
64    
65     //
66 caltinay 2810 // Destructor
67 caltinay 2183 //
68 caltinay 2810 DataVar::~DataVar()
69 caltinay 2183 {
70 caltinay 2810 cleanup();
71 caltinay 2183 }
72    
73     //
74     //
75 caltinay 2810 //
76     void DataVar::cleanup()
77 caltinay 2183 {
78     CoordArray::iterator it;
79 caltinay 2810 for (it = dataArray.begin(); it != dataArray.end(); it++)
80 caltinay 2183 delete[] *it;
81 caltinay 2810 dataArray.clear();
82     shape.clear();
83     sampleID.clear();
84     numSamples = 0;
85     initialized = false;
86 caltinay 2183 }
87    
88     //
89     //
90 caltinay 2810 //
91     bool DataVar::initFromEscript(escript::Data& escriptData, FinleyMesh_ptr mesh)
92 caltinay 2183 {
93 caltinay 2810 cleanup();
94 caltinay 2183
95 caltinay 2810 if (!escriptData.actsExpanded()) {
96     cerr << "WARNING: Only expanded data supported!" << endl;
97     return false;
98     }
99    
100 caltinay 2880 finleyMesh = mesh;
101 caltinay 2810 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 caltinay 2183 } else {
110 caltinay 2810 centering = ZONE_CENTERED;
111 caltinay 2183 }
112    
113 caltinay 2810 #ifdef _DEBUG
114     cout << varName << ":\t" << numSamples << " samples, "
115     << ptsPerSample << " pts/s, rank: " << rank << endl;
116     #endif
117 caltinay 2183
118 caltinay 2880 NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
119     if (nodes == NULL)
120     return false;
121    
122     meshName = nodes->getName();
123     siloMeshName = nodes->getFullSiloName();
124 caltinay 2810 initialized = true;
125 caltinay 2183
126 caltinay 2880 // no samples? Nothing more to do.
127 caltinay 2810 if (numSamples == 0)
128     return true;
129 caltinay 2183
130 caltinay 2810 const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
131     sampleID.insert(sampleID.end(), numSamples, 0);
132     copy(iPtr, iPtr+numSamples, sampleID.begin());
133 caltinay 2183
134 caltinay 2810 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 caltinay 2183
144 caltinay 2810 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 caltinay 2880 initialized = reorderSamples();
163 caltinay 2183 }
164 caltinay 2810
165     return initialized;
166 caltinay 2183 }
167    
168     //
169 caltinay 2810 // Initialise with mesh data
170 caltinay 2183 //
171 caltinay 2810 bool DataVar::initFromMesh(FinleyMesh_ptr mesh)
172 caltinay 2183 {
173 caltinay 2810 cleanup();
174    
175 caltinay 2880 finleyMesh = mesh;
176 caltinay 2810 rank = 0;
177     ptsPerSample = 1;
178 caltinay 2880 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 caltinay 2810 numSamples = data.size();
216 caltinay 2183
217 caltinay 2810 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 caltinay 2183 }
224 caltinay 2810 initialized = true;
225    
226     return initialized;
227 caltinay 2183 }
228    
229     //
230 caltinay 2810 // Reads variable data from NetCDF file
231 caltinay 2183 //
232 caltinay 2810 bool DataVar::initFromNetCDF(const string& filename, FinleyMesh_ptr mesh)
233 caltinay 2183 {
234 caltinay 2810 cleanup();
235    
236     #if USE_NETCDF
237 caltinay 2183 NcError ncerr(NcError::silent_nonfatal);
238 caltinay 2196 NcFile* input = new NcFile(filename.c_str());
239 caltinay 2183 if (!input->is_valid()) {
240 caltinay 2196 cerr << "Could not open input file " << filename << "." << endl;
241 caltinay 2183 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 caltinay 2810 cerr << "WARNING: Only expanded data supported!" << endl;
252 caltinay 2183 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 caltinay 2810 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 caltinay 2183 #ifdef _DEBUG
275     cout << varName << ":\t" << numSamples << " samples, "
276     << ptsPerSample << " pts/s, rank: " << rank << endl;
277     #endif
278    
279 caltinay 2880 finleyMesh = mesh;
280     NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
281     if (nodes == NULL) {
282 caltinay 2810 delete input;
283 caltinay 2880 return false;
284 caltinay 2810 }
285    
286 caltinay 2880 meshName = nodes->getName();
287     siloMeshName = nodes->getFullSiloName();
288     initialized = true;
289 caltinay 2183
290 caltinay 2810 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 caltinay 2183 }
300 caltinay 2810 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 caltinay 2880 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 caltinay 2810 size_t dataSize = dimSize*numSamples*ptsPerSample;
318     counts.push_back(ptsPerSample);
319     counts.push_back(numSamples);
320     float* tempData = new float[dataSize];
321 caltinay 2880 var = input->get_var("data");
322 caltinay 2810 var->get(tempData, &counts[0]);
323 caltinay 2183
324 caltinay 2810 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 caltinay 2183
331 caltinay 2880 initialized = reorderSamples();
332 caltinay 2810 }
333 caltinay 2183
334 caltinay 2810 delete input;
335     #endif // USE_NETCDF
336 caltinay 2183
337 caltinay 2810 return initialized;
338 caltinay 2183 }
339    
340     //
341     // Returns true if the data values are nodal, false if they are zonal.
342     //
343     bool DataVar::isNodeCentered() const
344     {
345 caltinay 2810 return (centering == NODE_CENTERED);
346 caltinay 2183 }
347    
348     //
349 caltinay 2810 // 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 caltinay 2183 //
354 caltinay 2810 float* DataVar::averageData(const float* src, size_t stride)
355 caltinay 2183 {
356 caltinay 2880 float* res;
357 caltinay 2183
358 caltinay 2810 if (ptsPerSample == 1) {
359 caltinay 2880 res = new float[numSamples];
360 caltinay 2810 float* dest = res;
361     for (int i=0; i<numSamples; i++, src+=stride)
362     *dest++ = *src;
363     } else {
364 caltinay 2880 ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
365     int cellFactor = cells->getElementFactor();
366     res = new float[cellFactor * numSamples];
367 caltinay 2810 float* dest = res;
368 caltinay 2880 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 caltinay 2183 }
394     }
395 caltinay 2810 return res;
396 caltinay 2183 }
397    
398     //
399 caltinay 2880 // 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 caltinay 2183 //
403 caltinay 2880 bool DataVar::reorderSamples()
404 caltinay 2183 {
405 caltinay 2810 if (numSamples == 0)
406 caltinay 2183 return true;
407    
408 caltinay 2810 const IntVec* requiredIDs = NULL;
409     int requiredNumSamples = 0;
410 caltinay 2880 int cellFactor = 1;
411 caltinay 2183
412 caltinay 2810 if (centering == NODE_CENTERED) {
413 caltinay 2880 NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
414 caltinay 2810 requiredIDs = &nodes->getNodeIDs();
415     requiredNumSamples = nodes->getNumNodes();
416     } else {
417     ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
418     if (cells == NULL)
419     return false;
420 caltinay 2183
421 caltinay 2810 requiredIDs = &cells->getIDs();
422 caltinay 2834 requiredNumSamples = cells->getNumElements();
423 caltinay 2880 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 caltinay 2183 }
436    
437 caltinay 2810 if (requiredNumSamples > numSamples) {
438     cerr << "ERROR: " << varName << " has " << numSamples
439     << " instead of " << requiredNumSamples << " samples!" << endl;
440 caltinay 2183 return false;
441     }
442    
443 caltinay 2810 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 caltinay 2880 size_t destIdx = 0;
452     for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
453 caltinay 2810 size_t srcIdx = sampleID2idx.find(*idIt)->second;
454 caltinay 2880 copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
455 caltinay 2810 }
456     delete[] dataArray[i];
457     dataArray[i] = c;
458     }
459 caltinay 2880
460     // sample IDs now = mesh node/element IDs
461     sampleID = *requiredIDs;
462    
463 caltinay 2183 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 caltinay 2810 if (rank < 2 || !initialized)
479 caltinay 2183 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 caltinay 2810 #if USE_SILO
522     if (!initialized)
523     return false;
524    
525 caltinay 2183 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 caltinay 2810
536     char* siloMesh = const_cast<char*>(siloMeshName.c_str());
537 caltinay 2183 int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
538    
539     if (rank == 0) {
540 caltinay 2810 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
541     numSamples, NULL, 0, DB_FLOAT, dcenter, NULL);
542 caltinay 2183 }
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 caltinay 2810 ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
552     (char**)varnames, &dataArray[0], numSamples, NULL,
553 caltinay 2183 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 caltinay 2810 ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
568     dataArray[i*shape[0]+j], numSamples,
569 caltinay 2183 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 caltinay 2810 #else // !USE_SILO
582 caltinay 2183 return false;
583     #endif
584     }
585    
586 caltinay 2810 } // namespace escriptexport
587 caltinay 2187

  ViewVC Help
Powered by ViewVC 1.1.26