/[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 3092 - (hide annotations)
Fri Aug 13 07:55:03 2010 UTC (8 years, 5 months ago) by caltinay
Original Path: trunk/weipa/src/DataVar.cpp
File size: 20499 byte(s)
Added getDataFlat() which will be used by the VisIt simulation interface

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

  ViewVC Help
Powered by ViewVC 1.1.26