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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3398 - (hide annotations)
Mon Dec 6 00:02:33 2010 UTC (8 years, 4 months ago) by caltinay
File size: 20201 byte(s)
Fixed the quad mask in a very specific case and changed saveVTK to write a
specific value for unused nodes instead of the first one which breaks unit
tests depending on the number of ranks used.
Still need to rebaseline some VTK files so the tests pass.

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

  ViewVC Help
Powered by ViewVC 1.1.26