/[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 3867 - (hide annotations)
Thu Mar 15 05:45:54 2012 UTC (6 years, 10 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 20347 byte(s)
weipa treats shapes (1,) and (1,1) as scalars. Fixes bogus output in VTK writer.

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

  ViewVC Help
Powered by ViewVC 1.1.26