/[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 3357 - (hide annotations)
Wed Nov 17 06:21:37 2010 UTC (8 years, 4 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 19471 byte(s)
Fixed a crash when loading data from netCDF files.

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     if (rank == 0) {
485     os << dataArray[0][index];
486     } else if (rank == 1) {
487     if (shape[0] < 3)
488     os << dataArray[0][index] << " " << dataArray[1][index]
489     << " " << 0.;
490     else
491     os << dataArray[0][index] << " " << dataArray[1][index]
492     << " " << dataArray[2][index];
493     } else if (rank == 2) {
494     if (shape[1] < 3) {
495     os << dataArray[0][index] << " " << dataArray[1][index]
496     << " " << 0. << " ";
497     os << dataArray[2][index] << " " << dataArray[3][index]
498     << " " << 0. << " ";
499     os << 0. << " " << 0. << " " << 0.;
500     } else {
501     os << dataArray[0][index] << " " << dataArray[1][index]
502     << " " << dataArray[2][index] << " ";
503     os << dataArray[3][index] << " " << dataArray[4][index]
504     << " " << dataArray[5][index] << " ";
505     os << dataArray[6][index] << " " << dataArray[7][index]
506     << " " << dataArray[8][index];
507     }
508     }
509     os << endl;
510     }
511    
512     //
513     //
514     //
515     void DataVar::writeToVTK(ostream& os, int ownIndex)
516     {
517     if (numSamples == 0)
518     return;
519    
520     if (isNodeCentered()) {
521     // data was reordered in reorderSamples() but for VTK we write the
522     // original node mesh and thus need the original ordering...
523 caltinay 3143 const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
524     const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
525     const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
526 caltinay 2886 int firstId = nodeDist[ownIndex];
527     int lastId = nodeDist[ownIndex+1];
528     IndexMap sampleID2idx = buildIndexMap();
529     for (int i=0; i<nodeGNI.size(); i++) {
530     if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
531     int idx = sampleID2idx[requiredIDs[i]];
532     sampleToStream(os, idx);
533     }
534     }
535     } else {
536     // cell data: ghost cells have been removed so do not write ghost
537     // samples (which are the last elements in the arrays)
538 caltinay 3143 int toWrite = domain->getElementsByName(meshName)->getNumElements();
539 caltinay 2886 for (int i=0; i<toWrite; i++) {
540     sampleToStream(os, i);
541     }
542     }
543     }
544    
545 caltinay 2183 ///////////////////////////////
546     // SILO related methods follow
547     ///////////////////////////////
548    
549     //
550     // If the data is tensor data then the components of the tensor are stored
551     // separately in the Silo file. This method then returns a string that
552     // contains the proper Silo expression to put the tensor together again.
553     // For non-tensor data this method returns an empty string.
554     //
555     string DataVar::getTensorDef() const
556     {
557 caltinay 2810 if (rank < 2 || !initialized)
558 caltinay 2183 return string();
559    
560     /// Format string for Silo 2x2 tensor
561     const string tensor2DefFmt =
562     "{{ <%sa_00>, <%sa_01> },"
563     " { <%sa_10>, <%sa_11> }}";
564    
565     /// Format string for Silo 3x3 tensor
566     const string tensor3DefFmt =
567     "{{ <%sa_00>, <%sa_01>, <%sa_02> },"
568     " { <%sa_10>, <%sa_11>, <%sa_12> },"
569     " { <%sa_20>, <%sa_21>, <%sa_22> }}";
570    
571     string tensorDef;
572     string tensorDir = varName+string("_comps/");
573     if (shape[1] == 3) {
574     char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
575     sprintf(tDef, tensor3DefFmt.c_str(),
576     tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
577     tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
578     tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
579     tensorDef = tDef;
580     delete[] tDef;
581     } else {
582     char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
583     sprintf(tDef, tensor2DefFmt.c_str(),
584     tensorDir.c_str(), tensorDir.c_str(),
585     tensorDir.c_str(), tensorDir.c_str(),
586     tensorDir.c_str(), tensorDir.c_str());
587     tensorDef = tDef;
588     delete[] tDef;
589     }
590     return tensorDef;
591     }
592    
593     //
594     // Writes the data to given Silo file under the virtual path provided.
595 caltinay 3143 // The corresponding mesh must have been written already.
596 caltinay 2183 //
597 caltinay 3128 bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
598     const string& units)
599 caltinay 2183 {
600 caltinay 2810 #if USE_SILO
601     if (!initialized)
602     return false;
603    
604 caltinay 2183 if (numSamples == 0)
605     return true;
606    
607     int ret;
608    
609     if (siloPath != "") {
610     ret = DBSetDir(dbfile, siloPath.c_str());
611     if (ret != 0)
612     return false;
613     }
614 caltinay 2810
615     char* siloMesh = const_cast<char*>(siloMeshName.c_str());
616 caltinay 2183 int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
617 caltinay 3128 DBoptlist* optList = DBMakeOptlist(2);
618     if (units.length()>0) {
619     DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
620     }
621 caltinay 2183
622     if (rank == 0) {
623 caltinay 2810 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
624 caltinay 3128 numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
625 caltinay 2183 }
626     else if (rank == 1) {
627     const string comps[3] = {
628     varName+string("_x"), varName+string("_y"), varName+string("_z")
629     };
630     const char* varnames[3] = {
631     comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
632     };
633    
634 caltinay 2810 ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
635     (char**)varnames, &dataArray[0], numSamples, NULL,
636 caltinay 3128 0, DB_FLOAT, dcenter, optList);
637 caltinay 2183 }
638     else {
639     string tensorDir = varName+string("_comps/");
640     ret = DBMkdir(dbfile, tensorDir.c_str());
641     if (ret == 0) {
642     int one = 1;
643     DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
644    
645     for (int i=0; i<shape[1]; i++) {
646     for (int j=0; j<shape[0]; j++) {
647     ostringstream varname;
648     varname << tensorDir << "a_" << i << j;
649 caltinay 2810 ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
650     dataArray[i*shape[0]+j], numSamples,
651 caltinay 2183 NULL, 0, DB_FLOAT, dcenter, optList);
652     if (ret != 0) break;
653     }
654     if (ret != 0) break;
655     }
656     } // ret==0
657     } // rank
658    
659 caltinay 3128 DBFreeOptlist(optList);
660 caltinay 2183 DBSetDir(dbfile, "/");
661     return (ret == 0);
662    
663 caltinay 2810 #else // !USE_SILO
664 caltinay 2183 return false;
665     #endif
666     }
667    
668 caltinay 3037 } // namespace weipa
669 caltinay 2187

  ViewVC Help
Powered by ViewVC 1.1.26