/[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 2834 - (hide annotations)
Thu Jan 7 06:06:56 2010 UTC (9 years, 2 months ago) by caltinay
Original Path: trunk/dataexporter/src/DataVar.cpp
File size: 15353 byte(s)
Reduced elements are now handled separate from the main elements within the
data exporter. This simplifies usage considerably.

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     varName(d.varName), numSamples(d.numSamples),
51     rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),
52 caltinay 2810 funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)
53 caltinay 2183 {
54 caltinay 2810 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 caltinay 2183 }
62 caltinay 2810 initialized = d.initialized;
63 caltinay 2183 }
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     rank = escriptData.getDataPointRank();
101     ptsPerSample = escriptData.getNumDataPointsPerSample();
102     shape = escriptData.getDataPointShape();
103     funcSpace = escriptData.getFunctionSpace().getTypeCode();
104     numSamples = escriptData.getNumSamples();
105    
106     if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
107     centering = NODE_CENTERED;
108 caltinay 2183 } else {
109 caltinay 2810 centering = ZONE_CENTERED;
110 caltinay 2183 }
111    
112 caltinay 2810 #ifdef _DEBUG
113     cout << varName << ":\t" << numSamples << " samples, "
114     << ptsPerSample << " pts/s, rank: " << rank << endl;
115     #endif
116 caltinay 2183
117 caltinay 2810 initialized = true;
118 caltinay 2183
119 caltinay 2810 if (numSamples == 0)
120     return true;
121 caltinay 2183
122 caltinay 2810 const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
123     sampleID.insert(sampleID.end(), numSamples, 0);
124     copy(iPtr, iPtr+numSamples, sampleID.begin());
125 caltinay 2183
126 caltinay 2810 size_t dimSize = 1;
127     if (rank > 0)
128     dimSize *= shape[0];
129     if (rank > 1)
130     dimSize *= shape[1];
131     if (rank > 2) {
132     cerr << "WARNING: Rank " << rank << " data is not supported!\n";
133     initialized = false;
134     }
135 caltinay 2183
136 caltinay 2810 if (initialized) {
137     size_t dataSize = dimSize * ptsPerSample;
138     float* tempData = new float[dataSize*numSamples];
139     float* destPtr = tempData;
140     for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
141     const escript::DataAbstract::ValueType::value_type* values =
142     escriptData.getSampleDataRO(sampleNo);
143     copy(values, values+dataSize, destPtr);
144     destPtr += dataSize;
145     }
146    
147     const float* srcPtr = tempData;
148     for (int i=0; i < dimSize; i++, srcPtr++) {
149     float* c = averageData(srcPtr, dimSize);
150     dataArray.push_back(c);
151     }
152     delete[] tempData;
153    
154     initialized = filterSamples(mesh);
155 caltinay 2183 }
156 caltinay 2810
157     return initialized;
158 caltinay 2183 }
159    
160     //
161 caltinay 2810 // Initialise with mesh data
162 caltinay 2183 //
163 caltinay 2810 bool DataVar::initFromMesh(FinleyMesh_ptr mesh)
164 caltinay 2183 {
165 caltinay 2810 cleanup();
166    
167     const IntVec& data = mesh->getVarDataByName(varName);
168     rank = 0;
169     ptsPerSample = 1;
170     numSamples = data.size();
171 caltinay 2183
172 caltinay 2810 if (numSamples > 0) {
173     float* c = new float[numSamples];
174     dataArray.push_back(c);
175     IntVec::const_iterator it;
176     for (it=data.begin(); it != data.end(); it++)
177     *c++ = static_cast<float>(*it);
178 caltinay 2183
179 caltinay 2834 if (varName.find("ContactElements_") != varName.npos) {
180     funcSpace = FINLEY_CONTACT_ELEMENTS_1;
181 caltinay 2810 centering = ZONE_CENTERED;
182     sampleID.insert(sampleID.end(),
183 caltinay 2834 mesh->getContactElements()->getIDs().begin(),
184     mesh->getContactElements()->getIDs().end());
185     } else if (varName.find("FaceElements_") != varName.npos) {
186 caltinay 2810 funcSpace = FINLEY_FACE_ELEMENTS;
187     centering = ZONE_CENTERED;
188     sampleID.insert(sampleID.end(),
189     mesh->getFaceElements()->getIDs().begin(),
190     mesh->getFaceElements()->getIDs().end());
191 caltinay 2834 } else if (varName.find("Elements_") != varName.npos) {
192     funcSpace = FINLEY_ELEMENTS;
193 caltinay 2810 centering = ZONE_CENTERED;
194     sampleID.insert(sampleID.end(),
195 caltinay 2834 mesh->getElements()->getIDs().begin(),
196     mesh->getElements()->getIDs().end());
197     } else if (varName.find("Nodes_") != varName.npos) {
198     funcSpace = FINLEY_NODES;
199     centering = NODE_CENTERED;
200     sampleID.insert(sampleID.end(),
201     mesh->getNodes()->getNodeIDs().begin(),
202     mesh->getNodes()->getNodeIDs().end());
203 caltinay 2810 } else {
204     return false;
205 caltinay 2183 }
206 caltinay 2810
207     NodeData_ptr nodes = mesh->getMeshForFinleyFS(funcSpace);
208     meshName = nodes->getName();
209     siloMeshName = nodes->getFullSiloName();
210 caltinay 2183 }
211 caltinay 2810 initialized = true;
212    
213     return initialized;
214 caltinay 2183 }
215    
216     //
217 caltinay 2810 // Reads variable data from NetCDF file
218 caltinay 2183 //
219 caltinay 2810 bool DataVar::initFromNetCDF(const string& filename, FinleyMesh_ptr mesh)
220 caltinay 2183 {
221 caltinay 2810 cleanup();
222    
223     #if USE_NETCDF
224 caltinay 2183 NcError ncerr(NcError::silent_nonfatal);
225 caltinay 2196 NcFile* input = new NcFile(filename.c_str());
226 caltinay 2183 if (!input->is_valid()) {
227 caltinay 2196 cerr << "Could not open input file " << filename << "." << endl;
228 caltinay 2183 delete input;
229     return false;
230     }
231    
232     NcDim* dim;
233     NcAtt* att;
234    
235     att = input->get_att("type_id");
236     int typeID = att->as_int(0);
237     if (typeID != 2) {
238 caltinay 2810 cerr << "WARNING: Only expanded data supported!" << endl;
239 caltinay 2183 delete input;
240     return false;
241     }
242    
243     att = input->get_att("rank");
244     rank = att->as_int(0);
245    
246     dim = input->get_dim("num_data_points_per_sample");
247     ptsPerSample = dim->size();
248    
249     att = input->get_att("function_space_type");
250     funcSpace = att->as_int(0);
251    
252 caltinay 2810 if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
253     centering = NODE_CENTERED;
254     } else {
255     centering = ZONE_CENTERED;
256     }
257    
258     dim = input->get_dim("num_samples");
259     numSamples = dim->size();
260    
261 caltinay 2183 #ifdef _DEBUG
262     cout << varName << ":\t" << numSamples << " samples, "
263     << ptsPerSample << " pts/s, rank: " << rank << endl;
264     #endif
265    
266 caltinay 2810 initialized = true;
267    
268     // if there are no data samples we're done
269     if (numSamples == 0) {
270     delete input;
271     return true;
272     }
273    
274 caltinay 2183 sampleID.insert(sampleID.end(), numSamples, 0);
275     NcVar* var = input->get_var("id");
276     var->get(&sampleID[0], numSamples);
277    
278 caltinay 2810 size_t dimSize = 1;
279     vector<long> counts;
280    
281     if (rank > 0) {
282     dim = input->get_dim("d0");
283     int d = dim->size();
284     shape.push_back(d);
285     counts.push_back(d);
286     dimSize *= d;
287 caltinay 2183 }
288 caltinay 2810 if (rank > 1) {
289     dim = input->get_dim("d1");
290     int d = dim->size();
291     shape.push_back(d);
292     counts.push_back(d);
293     dimSize *= d;
294     }
295     if (rank > 2) {
296     cerr << "WARNING: Rank " << rank << " data is not supported!\n";
297     initialized = false;
298     }
299    
300     if (initialized) {
301     size_t dataSize = dimSize*numSamples*ptsPerSample;
302     counts.push_back(ptsPerSample);
303     counts.push_back(numSamples);
304     float* tempData = new float[dataSize];
305     NcVar* var = input->get_var("data");
306     var->get(tempData, &counts[0]);
307 caltinay 2183
308 caltinay 2810 const float* srcPtr = tempData;
309     for (int i=0; i < dimSize; i++, srcPtr++) {
310     float* c = averageData(srcPtr, dimSize);
311     dataArray.push_back(c);
312     }
313     delete[] tempData;
314 caltinay 2183
315 caltinay 2810 initialized = filterSamples(mesh);
316     }
317 caltinay 2183
318 caltinay 2810 delete input;
319     #endif // USE_NETCDF
320 caltinay 2183
321 caltinay 2810 return initialized;
322 caltinay 2183 }
323    
324     //
325     // Returns true if the data values are nodal, false if they are zonal.
326     //
327     bool DataVar::isNodeCentered() const
328     {
329 caltinay 2810 return (centering == NODE_CENTERED);
330 caltinay 2183 }
331    
332     //
333 caltinay 2810 // Returns a subset of the src array according to stride parameter.
334     // If samples consist of multiple values they are averaged beforehand.
335     // Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
336     // (z0,z1,...)
337 caltinay 2183 //
338 caltinay 2810 float* DataVar::averageData(const float* src, size_t stride)
339 caltinay 2183 {
340 caltinay 2810 float* res = new float[numSamples];
341 caltinay 2183
342 caltinay 2810 if (ptsPerSample == 1) {
343     float* dest = res;
344     for (int i=0; i<numSamples; i++, src+=stride)
345     *dest++ = *src;
346     } else {
347     float* dest = res;
348     for (int i=0; i<numSamples; i++) {
349     double tmpVal = 0.0;
350     for (int j=0; j<ptsPerSample; j++, src+=stride)
351     tmpVal += *src;
352     *dest++ = (float)(tmpVal / ptsPerSample);
353 caltinay 2183 }
354     }
355 caltinay 2810 return res;
356 caltinay 2183 }
357    
358     //
359 caltinay 2810 // Filters and reorders the raw sample values according to the IDs provided
360     // in 'requiredIDs'. This is used to have data arrays ordered according to
361     // the underlying mesh (i.e. DataID[i]==MeshNodeID[i])
362 caltinay 2183 //
363 caltinay 2810 bool DataVar::filterSamples(FinleyMesh_ptr finleyMesh)
364 caltinay 2183 {
365 caltinay 2810 if (numSamples == 0)
366 caltinay 2183 return true;
367    
368 caltinay 2810 IndexMap id2idxMap;
369     const IntVec* requiredIDs = NULL;
370 caltinay 2183
371 caltinay 2810 NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
372     if (nodes == NULL)
373     return false;
374 caltinay 2183
375 caltinay 2810 int requiredNumSamples = 0;
376 caltinay 2183
377 caltinay 2810 if (centering == NODE_CENTERED) {
378     id2idxMap = nodes->getIndexMap();
379     requiredIDs = &nodes->getNodeIDs();
380     requiredNumSamples = nodes->getNumNodes();
381     } else {
382     ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
383     if (cells == NULL)
384     return false;
385 caltinay 2183
386 caltinay 2810 id2idxMap = cells->getIndexMap();
387     requiredIDs = &cells->getIDs();
388 caltinay 2834 requiredNumSamples = cells->getNumElements();
389 caltinay 2183 }
390    
391 caltinay 2810 if (requiredNumSamples > numSamples) {
392     cerr << "ERROR: " << varName << " has " << numSamples
393     << " instead of " << requiredNumSamples << " samples!" << endl;
394 caltinay 2183 return false;
395     }
396    
397 caltinay 2810 meshName = nodes->getName();
398     siloMeshName = nodes->getFullSiloName();
399 caltinay 2183
400 caltinay 2810 IndexMap sampleID2idx = buildIndexMap();
401     numSamples = requiredNumSamples;
402    
403     // now filter the data
404     for (size_t i=0; i < dataArray.size(); i++) {
405     float* c = new float[numSamples];
406     const float* src = dataArray[i];
407     IntVec::const_iterator idIt = requiredIDs->begin();
408     for (; idIt != requiredIDs->end(); idIt++) {
409     size_t srcIdx = sampleID2idx.find(*idIt)->second;
410     size_t destIdx = id2idxMap.find(*idIt)->second;
411     c[destIdx] = src[srcIdx];
412     }
413     delete[] dataArray[i];
414     dataArray[i] = c;
415     }
416 caltinay 2183 return true;
417     }
418    
419     ///////////////////////////////
420     // SILO related methods follow
421     ///////////////////////////////
422    
423     //
424     // If the data is tensor data then the components of the tensor are stored
425     // separately in the Silo file. This method then returns a string that
426     // contains the proper Silo expression to put the tensor together again.
427     // For non-tensor data this method returns an empty string.
428     //
429     string DataVar::getTensorDef() const
430     {
431 caltinay 2810 if (rank < 2 || !initialized)
432 caltinay 2183 return string();
433    
434     /// Format string for Silo 2x2 tensor
435     const string tensor2DefFmt =
436     "{{ <%sa_00>, <%sa_01> },"
437     " { <%sa_10>, <%sa_11> }}";
438    
439     /// Format string for Silo 3x3 tensor
440     const string tensor3DefFmt =
441     "{{ <%sa_00>, <%sa_01>, <%sa_02> },"
442     " { <%sa_10>, <%sa_11>, <%sa_12> },"
443     " { <%sa_20>, <%sa_21>, <%sa_22> }}";
444    
445     string tensorDef;
446     string tensorDir = varName+string("_comps/");
447     if (shape[1] == 3) {
448     char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
449     sprintf(tDef, tensor3DefFmt.c_str(),
450     tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
451     tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
452     tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
453     tensorDef = tDef;
454     delete[] tDef;
455     } else {
456     char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
457     sprintf(tDef, tensor2DefFmt.c_str(),
458     tensorDir.c_str(), tensorDir.c_str(),
459     tensorDir.c_str(), tensorDir.c_str(),
460     tensorDir.c_str(), tensorDir.c_str());
461     tensorDef = tDef;
462     delete[] tDef;
463     }
464     return tensorDef;
465     }
466    
467     //
468     // Writes the data to given Silo file under the virtual path provided.
469     // The corresponding mesh must have been written already and made known
470     // to this variable by a call to setMesh().
471     //
472     bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)
473     {
474 caltinay 2810 #if USE_SILO
475     if (!initialized)
476     return false;
477    
478 caltinay 2183 if (numSamples == 0)
479     return true;
480    
481     int ret;
482    
483     if (siloPath != "") {
484     ret = DBSetDir(dbfile, siloPath.c_str());
485     if (ret != 0)
486     return false;
487     }
488 caltinay 2810
489     char* siloMesh = const_cast<char*>(siloMeshName.c_str());
490 caltinay 2183 int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
491    
492     if (rank == 0) {
493 caltinay 2810 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
494     numSamples, NULL, 0, DB_FLOAT, dcenter, NULL);
495 caltinay 2183 }
496     else if (rank == 1) {
497     const string comps[3] = {
498     varName+string("_x"), varName+string("_y"), varName+string("_z")
499     };
500     const char* varnames[3] = {
501     comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
502     };
503    
504 caltinay 2810 ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
505     (char**)varnames, &dataArray[0], numSamples, NULL,
506 caltinay 2183 0, DB_FLOAT, dcenter, NULL);
507     }
508     else {
509     string tensorDir = varName+string("_comps/");
510     ret = DBMkdir(dbfile, tensorDir.c_str());
511     if (ret == 0) {
512     int one = 1;
513     DBoptlist* optList = DBMakeOptlist(1);
514     DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
515    
516     for (int i=0; i<shape[1]; i++) {
517     for (int j=0; j<shape[0]; j++) {
518     ostringstream varname;
519     varname << tensorDir << "a_" << i << j;
520 caltinay 2810 ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
521     dataArray[i*shape[0]+j], numSamples,
522 caltinay 2183 NULL, 0, DB_FLOAT, dcenter, optList);
523     if (ret != 0) break;
524     }
525     if (ret != 0) break;
526     }
527     DBFreeOptlist(optList);
528     } // ret==0
529     } // rank
530    
531     DBSetDir(dbfile, "/");
532     return (ret == 0);
533    
534 caltinay 2810 #else // !USE_SILO
535 caltinay 2183 return false;
536     #endif
537     }
538    
539 caltinay 2810 } // namespace escriptexport
540 caltinay 2187

  ViewVC Help
Powered by ViewVC 1.1.26