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

  ViewVC Help
Powered by ViewVC 1.1.26