/[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 3185 - (hide annotations)
Thu Sep 16 00:30:09 2010 UTC (8 years, 7 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 19423 byte(s)
Weipa now supports initialization from DataConstant instances.

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

  ViewVC Help
Powered by ViewVC 1.1.26