/[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 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 3 months ago) by jfenwick
Original Path: trunk/weipa/src/DataVar.cpp
File size: 19455 byte(s)
Merging dudley and scons updates from branches

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

  ViewVC Help
Powered by ViewVC 1.1.26