/[escript]/trunk/weipa/src/DataVar.cpp
ViewVC logotype

Annotation of /trunk/weipa/src/DataVar.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (12 months, 2 weeks ago) by jfenwick
File size: 23346 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26