/[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 4763 - (hide annotations)
Tue Mar 18 05:20:23 2014 UTC (4 years, 7 months ago) by jfenwick
File size: 20527 byte(s)
Replaced references to getSampleDataRO(0) with getDataRO()
This was to deal with the case where processes with no data where trying to compute the offset of sample 0 when all that was really needed was a null.

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

  ViewVC Help
Powered by ViewVC 1.1.26