/[escript]/branches/split/weipa/src/DataVar.cpp
ViewVC logotype

Diff of /branches/split/weipa/src/DataVar.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/tools/libescriptreader/src/escriptreader/DataVar.cpp revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC trunk/weipa/src/DataVar.cpp revision 4657 by jfenwick, Thu Feb 6 06:12:20 2014 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14    *
15    *****************************************************************************/
16    
17  //  #include <weipa/DataVar.h>
18  // DataVar.cpp  #include <weipa/DomainChunk.h>
19  //  #include <weipa/ElementData.h>
20  #include <escriptreader/DataVar.h>  #include <weipa/NodeData.h>
21  #include <escriptreader/ElementData.h>  #ifndef VISIT_PLUGIN
22  #include <escriptreader/MeshWithElements.h>  #include <escript/Data.h>
23  #include <netcdf.hh>  #endif
24  #if HAVE_SILO  
25    #if USE_NETCDF
26    #include <netcdfcpp.h>
27    #endif
28    
29    #if USE_SILO
30  #include <silo.h>  #include <silo.h>
31  #endif  #endif
32    
33    #include <numeric> // for accumulate
34    #include <iostream> // for cerr
35    #include <stdio.h>
36    
37  using namespace std;  using namespace std;
38    
39  namespace EscriptReader {  namespace weipa {
40            
 enum {  
     NODE_CENTERED = 1,  
     ZONE_CENTERED = 2  
 };  
   
41  //  //
42  // Constructor  // Constructor
43  //  //
44  DataVar::DataVar(const string& name) :  DataVar::DataVar(const string& name) :
45      varName(name), numSamples(0), rank(0), ptsPerSample(0), centering(0),      initialized(false), varName(name),
46      reorderedNumSamples(0), fullMesh(NULL)      numSamples(0), rank(0), ptsPerSample(0)
47  {  {
48  }  }
49    
50  //  //
 // Destructor  
 //  
 DataVar::~DataVar()  
 {  
     CoordArray::iterator it;  
     for (it = reorderedData.begin(); it != reorderedData.end(); it++)  
         delete[] *it;  
     for (it = rawData.begin(); it != rawData.end(); it++)  
         delete[] *it;  
 }  
   
 //  
51  // Copy constructor  // Copy constructor
52  //  //
53  DataVar::DataVar(const DataVar& d) :  DataVar::DataVar(const DataVar& d) :
54        initialized(d.initialized), domain(d.domain),
55      varName(d.varName), numSamples(d.numSamples),      varName(d.varName), numSamples(d.numSamples),
56      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),      rank(d.rank), ptsPerSample(d.ptsPerSample), funcSpace(d.funcSpace),
57      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID),      centering(d.centering), shape(d.shape), sampleID(d.sampleID)
     reorderedNumSamples(d.reorderedNumSamples), fullMesh(d.fullMesh)  
58  {  {
59      CoordArray::const_iterator it;      if (numSamples > 0) {
60      for (it = d.rawData.begin(); it != d.rawData.end(); it++) {          CoordArray::const_iterator it;
61          float* c = new float[numSamples];          for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
62          copy(*it, (*it)+numSamples, c);              float* c = new float[numSamples];
63          rawData.push_back(c);              copy(*it, (*it)+numSamples, c);
64      }              dataArray.push_back(c);
65      for (it = d.reorderedData.begin(); it != d.reorderedData.end(); it++) {          }
         float* c = new float[reorderedNumSamples];  
         copy(*it, (*it)+reorderedNumSamples, c);  
         reorderedData.push_back(c);  
66      }      }
67  }  }
68    
69  //  //
70  // Special constructor for mesh data  // Destructor
71  //  //
72  DataVar::DataVar(const string& name, const IntVec& data,  DataVar::~DataVar()
                  MeshWithElements* mesh) :  
     varName(name)  
73  {  {
74      numSamples = data.size();      cleanup();
   
     float* c = new float[numSamples];  
     rawData.push_back(c);  
     IntVec::const_iterator it;  
     for (it=data.begin(); it != data.end(); it++)  
         *c++ = static_cast<float>(*it);  
   
     rank = 0;  
     ptsPerSample = 1;  
     if (name.compare(0, 6, "Nodes_") == 0) {  
         funcSpace = FINLEY_NODES;  
         centering = NODE_CENTERED;  
         sampleID.insert(sampleID.end(), mesh->getNodeIDs().begin(),  
                 mesh->getNodeIDs().end());  
     } else if (name.compare(0, 9, "Elements_") == 0) {  
         funcSpace = FINLEY_ELEMENTS;  
         centering = ZONE_CENTERED;  
         sampleID.insert(sampleID.end(), mesh->getElements()->getIDs().begin(),  
                 mesh->getElements()->getIDs().end());  
     } else if (name.compare(0, 13, "FaceElements_") == 0) {  
         funcSpace = FINLEY_FACE_ELEMENTS;  
         centering = ZONE_CENTERED;  
         sampleID.insert(sampleID.end(),  
                 mesh->getFaceElements()->getIDs().begin(),  
                 mesh->getFaceElements()->getIDs().end());  
     } else if (name.compare(0, 16, "ContactElements_") == 0) {  
         funcSpace = FINLEY_CONTACT_ELEMENTS_1;  
         centering = ZONE_CENTERED;  
         sampleID.insert(sampleID.end(),  
                 mesh->getContactElements()->getIDs().begin(),  
                 mesh->getContactElements()->getIDs().end());  
     } else if (name.compare(0, 7, "Points_") == 0) {  
         funcSpace = FINLEY_POINTS;  
         centering = NODE_CENTERED;  
         sampleID.insert(sampleID.end(), mesh->getPoints()->getIDs().begin(),  
                 mesh->getPoints()->getIDs().end());  
     }  
   
     shape.clear();  
     reorderedNumSamples = 0;  
75  }  }
76    
77  //  //
 // Appends raw data including IDs from rhs.  
78  //  //
79  bool DataVar::append(const DataVar& rhs)  //
80    void DataVar::cleanup()
81  {  {
     // check if variables are compatible  
     if (varName != rhs.varName || ptsPerSample != rhs.ptsPerSample ||  
             rank != rhs.rank || shape.size() != rhs.shape.size() ||  
             centering != rhs.centering)  
         return false;  
   
     for (size_t i=0; i<shape.size(); i++)  
         if (shape[i] != rhs.shape[i])  
             return false;  
   
     sampleID.insert(sampleID.end(), rhs.sampleID.begin(), rhs.sampleID.end());  
     for (size_t i=0; i<rawData.size(); i++) {  
         float* c = new float[numSamples+rhs.numSamples];  
         copy(rawData[i], rawData[i]+numSamples, c);  
         copy(rhs.rawData[i], rhs.rawData[i]+rhs.numSamples, c+numSamples);  
         delete[] rawData[i];  
         rawData[i] = c;  
     }  
     numSamples += rhs.numSamples;  
   
     // invalidate previously reordered data  
82      CoordArray::iterator it;      CoordArray::iterator it;
83      for (it = reorderedData.begin(); it != reorderedData.end(); it++)      for (it = dataArray.begin(); it != dataArray.end(); it++)
84          delete[] *it;          delete[] *it;
85      reorderedData.clear();      dataArray.clear();
86      reorderedNumSamples = 0;      shape.clear();
87            sampleID.clear();
88      return true;      numSamples = 0;
89        initialized = false;
90  }  }
91    
92  //  //
 // Returns a subset of the src array according to stride parameter.  
 // If samples consist of multiple values they are averaged beforehand.  
 // Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and  
 // (z0,z1,...)  
93  //  //
94  float* DataVar::averageData(const float* src, size_t stride)  //
95    bool DataVar::initFromEscript(escript::Data& escriptData, const_DomainChunk_ptr dom)
96  {  {
97      float* res = new float[numSamples];  #ifndef VISIT_PLUGIN
98        cleanup();
99    
100      if (ptsPerSample == 1) {      if (!escriptData.isConstant() && !escriptData.actsExpanded()) {
101          float* dest = res;          cerr << "WARNING: Weipa only supports constant & expanded data, "
102          for (int i=0; i<numSamples; i++, src+=stride)              << "not initializing " << varName << endl;
103              *dest++ = *src;          return false;
     } else {  
         float* dest = res;  
         for (int i=0; i<numSamples; i++) {  
             double tmpVal = 0.0;  
             for (int j=0; j<ptsPerSample; j++, src+=stride)  
                 tmpVal += *src;  
             *dest++ = (float)(tmpVal / ptsPerSample);  
         }  
104      }      }
     return res;  
 }  
105    
106  //      domain = dom;
107  // Reads scalar data (rank=0) from NetCDF file and stores the values      rank = escriptData.getDataPointRank();
108  // after averaging.      ptsPerSample = escriptData.getNumDataPointsPerSample();
109  //      shape = escriptData.getDataPointShape();
110  void DataVar::readRank0Data(NcFile* ncfile)      funcSpace = escriptData.getFunctionSpace().getTypeCode();
111  {      numSamples = escriptData.getNumSamples();
112      shape.clear();      centering = domain->getCenteringForFunctionSpace(funcSpace);
     float* tempData = new float[ptsPerSample*numSamples];  
     NcVar* var = ncfile->get_var("data");  
     var->get(tempData, ptsPerSample, numSamples);  
113    
114      float* c = averageData(tempData, 1);  #ifdef _DEBUG
115      rawData.push_back(c);      cout << varName << ":\t" << numSamples << " samples,  "
116            << ptsPerSample << " pts/s,  rank: " << rank << endl;
117    #endif
118    
119      delete[] tempData;      NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
120  }      if (nodes == NULL)
121            return false;
122    
123  //      meshName = nodes->getName();
124  // Reads vector data (rank=1) from NetCDF file and stores the components      siloMeshName = nodes->getFullSiloName();
125  // separately after averaging.      initialized = true;
126  //  
127  void DataVar::readRank1Data(NcFile* ncfile)      // no samples? Nothing more to do.
128  {      if (numSamples == 0)
129      shape.clear();          return true;
130      NcDim* dim = ncfile->get_dim("d0");  
131      shape.push_back(dim->size());      const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
132        sampleID.insert(sampleID.end(), numSamples, 0);
133        copy(iPtr, iPtr+numSamples, sampleID.begin());
134    
135      float* tempData = new float[shape[0]*ptsPerSample*numSamples];      size_t dimSize = 1;
136      NcVar* var = ncfile->get_var("data");      if (rank > 0)
137      var->get(tempData, shape[0], ptsPerSample, numSamples);          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    
145        // 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        if (initialized) {
152            size_t dataSize = dimSize * ptsPerSample;
153            float* tempData = new float[dataSize*numSamples];
154            float* destPtr = tempData;
155            if (escriptData.isConstant()) {
156                const escript::DataAbstract::ValueType::value_type* values =
157                    escriptData.getSampleDataRO(0);
158                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            }
170    
171      for (int i=0; i<shape[0]; i++) {          const float* srcPtr = tempData;
172          const float* src = tempData;          for (int i=0; i < dimSize; i++, srcPtr++) {
173          src+=i;              float* c = averageData(srcPtr, dimSize);
174          float* c = averageData(src, shape[0]);              dataArray.push_back(c);
175          rawData.push_back(c);          }
176            delete[] tempData;
177    
178            initialized = reorderSamples();
179      }      }
180      delete[] tempData;  
181        return initialized;
182    
183    #else // VISIT_PLUGIN
184        return false;
185    #endif
186  }  }
187    
188  //  //
189  // Like readRank1Data() but for tensor data (rank=2).  // Initialise with domain data
190  //  //
191  void DataVar::readRank2Data(NcFile* ncfile)  bool DataVar::initFromMeshData(const_DomainChunk_ptr dom, const IntVec& data,
192            int fsCode, Centering c, NodeData_ptr nodes, const IntVec& id)
193  {  {
194      shape.clear();      cleanup();
195      NcDim* dim = ncfile->get_dim("d0");      
196      shape.push_back(dim->size());      domain = dom;
197      dim = ncfile->get_dim("d1");      rank = 0;
198      shape.push_back(dim->size());      ptsPerSample = 1;
199        centering = c;
200      float* tempData = new float[shape[0]*shape[1]*ptsPerSample*numSamples];      sampleID = id;
201      NcVar* var = ncfile->get_var("data");      meshName = nodes->getName();
202      var->get(tempData, shape[0], shape[1], ptsPerSample, numSamples);      siloMeshName = nodes->getFullSiloName();
203        numSamples = data.size();
204    
205      for (int i=0; i<shape[1]; i++) {      if (numSamples > 0) {
206          for (int j=0; j<shape[0]; j++) {          float* c = new float[numSamples];
207              const float* src = tempData;          dataArray.push_back(c);
208              src+=i*shape[0]+j;          IntVec::const_iterator it;
209              float* c = averageData(src, shape[0]*shape[1]);          for (it=data.begin(); it != data.end(); it++)
210              rawData.push_back(c);              *c++ = static_cast<float>(*it);
         }  
211      }      }
212      delete[] tempData;      initialized = true;
213    
214        return initialized;
215  }  }
216    
217  //  //
218  // Reads a NetCDF file in escript/finley format.  // Reads variable data from dump file
219  //  //
220  bool DataVar::readFromNc(const string& filename)  bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
221  {  {
222        cleanup();
223        
224    #if USE_NETCDF
225      NcError ncerr(NcError::silent_nonfatal);          NcError ncerr(NcError::silent_nonfatal);    
226      NcFile* input = new NcFile(filename.c_str());      NcFile* input = new NcFile(filename.c_str());
227      if (!input->is_valid()) {      if (!input->is_valid()) {
# Line 267  bool DataVar::readFromNc(const string& f Line 233  bool DataVar::readFromNc(const string& f
233      NcDim* dim;      NcDim* dim;
234      NcAtt* att;      NcAtt* att;
235    
     dim = input->get_dim("num_samples");  
     numSamples = dim->size();  
   
     // if there are no data samples bail out  
     if (numSamples == 0) {  
         delete input;  
         return false;  
     }  
   
236      att = input->get_att("type_id");      att = input->get_att("type_id");
237      int typeID = att->as_int(0);      int typeID = att->as_int(0);
238      if (typeID != 2) {      if (typeID != 2) {
239          cerr << "WARNING: Only expanded data supported at the moment!" << endl;          cerr << "WARNING: Only expanded data supported!" << endl;
240          delete input;          delete input;
241          return false;          return false;
242      }      }
# Line 293  bool DataVar::readFromNc(const string& f Line 250  bool DataVar::readFromNc(const string& f
250      att = input->get_att("function_space_type");      att = input->get_att("function_space_type");
251      funcSpace = att->as_int(0);      funcSpace = att->as_int(0);
252    
253        centering = dom->getCenteringForFunctionSpace(funcSpace);
254    
255        dim = input->get_dim("num_samples");
256        numSamples = dim->size();
257    
258  #ifdef _DEBUG  #ifdef _DEBUG
259      cout << varName << ":\t" << numSamples << " samples,  "      cout << varName << ":\t" << numSamples << " samples,  "
260          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
261  #endif  #endif
262    
263      sampleID.clear();      domain = dom;
264      sampleID.insert(sampleID.end(), numSamples, 0);      NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
265      NcVar* var = input->get_var("id");      if (nodes == NULL) {
266      var->get(&sampleID[0], numSamples);          delete input;
267            return false;
     switch (rank) {  
         case 0:  
             readRank0Data(input);  
             break;  
         case 1:  
             readRank1Data(input);  
             break;  
         case 2:  
             readRank2Data(input);  
             break;  
         default:  
             cerr << "WARNING: Rank " << rank << " data is not supported!\n";  
             delete input;  
             return false;  
268      }      }
269    
270      delete input;      meshName = nodes->getName();
271      return true;      siloMeshName = nodes->getFullSiloName();
272  }      initialized = true;
273    
274  //      size_t dimSize = 1;
275  // Returns one of the mesh names provided by mainMesh that matches the      vector<long> counts;
276  // data variable's function space type and reduced/unreduced state.  
277  //      if (rank > 0) {
278  string DataVar::getMeshName(MeshWithElements* mainMesh) const          dim = input->get_dim("d0");
279  {          int d = dim->size();
280      string name;          shape.push_back(d);
281            counts.push_back(d);
282      switch (funcSpace) {          dimSize *= d;
283          case FINLEY_REDUCED_NODES:      }
284          case FINLEY_REDUCED_DEGREES_OF_FREEDOM:      if (rank > 1) {
285          case FINLEY_REDUCED_ELEMENTS:          dim = input->get_dim("d1");
286          case FINLEY_ELEMENTS:          int d = dim->size();
287              if (mainMesh->getElements()->reducedCount > 0) {          shape.push_back(d);
288                  name = mainMesh->getElements()->reducedMesh->getName();          counts.push_back(d);
289              } else {          dimSize *= d;
290                  name = mainMesh->getElements()->fullMesh->getName();      }
291              }      if (rank > 2) {
292              break;          cerr << "WARNING: Rank " << rank << " data is not supported!\n";
293            initialized = false;
294        }
295    
296        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            size_t dataSize = dimSize*numSamples*ptsPerSample;
302            counts.push_back(ptsPerSample);
303            counts.push_back(numSamples);
304            float* tempData = new float[dataSize];
305            var = input->get_var("data");
306            var->get(tempData, &counts[0]);
307    
308            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    
315          case FINLEY_REDUCED_FACE_ELEMENTS:          initialized = reorderSamples();
316          case FINLEY_FACE_ELEMENTS:      }
             if (mainMesh->getFaceElements()->reducedCount > 0) {  
                 name = mainMesh->getFaceElements()->reducedMesh->getName();  
             } else {  
                 name = mainMesh->getFaceElements()->fullMesh->getName();  
             }  
             break;  
317    
318          case FINLEY_REDUCED_CONTACT_ELEMENTS_1:      delete input;
319          case FINLEY_REDUCED_CONTACT_ELEMENTS_2:  #endif // USE_NETCDF
         case FINLEY_CONTACT_ELEMENTS_1:  
         case FINLEY_CONTACT_ELEMENTS_2:  
             if (mainMesh->getContactElements()->reducedCount > 0) {  
                 name = mainMesh->getContactElements()->reducedMesh->getName();  
             } else {  
                 name = mainMesh->getContactElements()->fullMesh->getName();  
             }  
             break;  
320    
321          case FINLEY_NODES:      return initialized;
         case FINLEY_DEGREES_OF_FREEDOM:  
             name = mainMesh->getElements()->fullMesh->getName();  
             break;  
   
         case FINLEY_POINTS:  
             name = mainMesh->getPoints()->fullMesh->getName();  
             break;  
     }  
     return name;  
322  }  }
323    
324  //  //
# Line 380  string DataVar::getMeshName(MeshWithElem Line 326  string DataVar::getMeshName(MeshWithElem
326  //  //
327  bool DataVar::isNodeCentered() const  bool DataVar::isNodeCentered() const
328  {  {
329      return (funcSpace == FINLEY_REDUCED_NODES ||      return (centering == NODE_CENTERED);
             funcSpace == FINLEY_REDUCED_DEGREES_OF_FREEDOM ||  
             funcSpace == FINLEY_NODES ||  
             funcSpace == FINLEY_DEGREES_OF_FREEDOM ||  
             funcSpace == FINLEY_POINTS);  
330  }  }
331    
332  //  //
333  // Filters and reorders the raw sample values according to the IDs provided  // Returns a subset of the src array according to stride parameter.
334  // in 'requiredIDs'. This is used to have data arrays ordered according to  // If samples consist of multiple values they are averaged beforehand.
335  // the underlying mesh (i.e. DataID[i]==MeshNodeID[i])  // Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
336    // (z0,z1,...)
337  //  //
338  void DataVar::reorderSamples(const IndexMap& id2idxMap,  float* DataVar::averageData(const float* src, size_t stride)
                              const IntVec& requiredIDs)  
339  {  {
340      CoordArray::iterator it;      float* res;
     for (it = reorderedData.begin(); it != reorderedData.end(); it++)  
         delete[] *it;  
     reorderedData.clear();  
341    
342      buildIndexMap();      if (ptsPerSample == 1) {
343      for (size_t i=0; i < rawData.size(); i++) {          res = new float[numSamples];
344          float* c = new float[reorderedNumSamples];          float* dest = res;
345          reorderedData.push_back(c);          for (int i=0; i<numSamples; i++, src+=stride)
346          const float* src = rawData[i];              *dest++ = *src;
347          IntVec::const_iterator idIt = requiredIDs.begin();      } else {
348          for (; idIt != requiredIDs.end(); idIt++) {          ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
349              size_t srcIdx = sampleID2idx.find(*idIt)->second;          int cellFactor = cells->getElementFactor();
350              size_t destIdx = id2idxMap.find(*idIt)->second;          res = new float[cellFactor * numSamples];
351              c[destIdx] = src[srcIdx];          float* dest = res;
352            QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
353            if (!qmi.mask.empty()) {
354                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          }          }
378      }      }
379        return res;
380  }  }
381    
382  //  //
383  // For zonal data this method reorders the values according to the indices  // Filters and reorders the raw sample values according to the node/element
384  // given in reorderArray. This is used to move ghost zones to the end of  // IDs. This is used to have data arrays ordered according to the underlying
385  // the arrays which conforms to Silo's expected format.  // mesh (i.e. DataID[i]==MeshNodeID[i])
 // Nodal data is not changed by this method.  
386  //  //
387  void DataVar::handleGhostZones(const IntVec& reorderArray)  bool DataVar::reorderSamples()
388  {  {
389      if (centering == NODE_CENTERED)      if (numSamples == 0)
390          return;          return true;
391    
392      vector<float*>::iterator it;      const IntVec* requiredIDs = NULL;
393      for (it = reorderedData.begin(); it!=reorderedData.end(); it++) {      int requiredNumSamples = 0;
394          float* original = *it;      int cellFactor = 1;
395          float* reordered = new float[reorderedNumSamples];  
396          float* arrIt = reordered;      if (centering == NODE_CENTERED) {
397          IntVec::const_iterator idxIt;          NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
398          for (idxIt=reorderArray.begin(); idxIt!=reorderArray.end(); idxIt++)          requiredIDs = &nodes->getNodeIDs();
399              *arrIt++ = original[*idxIt];          requiredNumSamples = nodes->getNumNodes();
400        } else {
401            ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
402            if (cells == NULL)
403                return false;
404    
405          delete[] *it;          requiredIDs = &cells->getIDs();
406          *it = reordered;          requiredNumSamples = cells->getNumElements();
407            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        }
420    
421        if (requiredNumSamples > numSamples) {
422            cerr << "ERROR: " << varName << " has " << numSamples
423                << " instead of " << requiredNumSamples << " samples!" << endl;
424            return false;
425        }
426    
427        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            size_t destIdx = 0;
436            for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
437                size_t srcIdx = sampleID2idx.find(*idIt)->second;
438                copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
439            }
440            delete[] dataArray[i];
441            dataArray[i] = c;
442      }      }
443    
444        // sample IDs now = mesh node/element IDs
445        sampleID = *requiredIDs;
446    
447        return true;
448  }  }
449    
450  //  //
 // Makes the top-level mesh known to this data variable. The mesh is used  
 // to reorder and filter the samples and inquire the number of ghost zones.  
451  //  //
452  bool DataVar::setMesh(MeshWithElements* mesh)  //
453    int DataVar::getNumberOfComponents() const
454  {  {
455      if (fullMesh == mesh)      return (rank == 0 ? 1 : accumulate(shape.begin(), shape.end(), 0));
456          return true;  }
   
     const IndexMap* id2idxMap;  
     const IntVec* reqIDs;  
     const IntVec* reorderArray = NULL;  
   
     switch (funcSpace) {  
         case FINLEY_REDUCED_NODES:  
         case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
             {  
                 centering = NODE_CENTERED;  
                 ElementData* cells = mesh->getElements();  
                 if (cells->reducedCount > 0) {  
                     if (cells->getReducedGhostCount())  
                         reorderArray = &cells->reducedIndexArray;  
                     siloMeshName = cells->reducedMesh->getFullSiloName();  
                     id2idxMap = &cells->reducedMesh->getIndexMap();  
                     reqIDs = &cells->reducedMesh->getNodeIDs();  
                     reorderedNumSamples = cells->reducedMesh->getNumNodes();  
                 } else {  
                     if (cells->getGhostCount())  
                         reorderArray = &cells->indexArray;  
                     siloMeshName = cells->fullMesh->getFullSiloName();  
                     id2idxMap = &cells->fullMesh->getIndexMap();  
                     reqIDs = &cells->fullMesh->getNodeIDs();  
                     reorderedNumSamples = cells->fullMesh->getNumNodes();  
                 }  
             }  
             break;  
   
         case FINLEY_NODES:  
         case FINLEY_DEGREES_OF_FREEDOM:  
             {  
                 centering = NODE_CENTERED;  
                 ElementData* cells = mesh->getElements();  
                 if (cells->getGhostCount())  
                     reorderArray = &cells->indexArray;  
                 siloMeshName = cells->fullMesh->getFullSiloName();  
                 id2idxMap = &cells->fullMesh->getIndexMap();  
                 reqIDs = &cells->fullMesh->getNodeIDs();  
                 reorderedNumSamples = cells->fullMesh->getNumNodes();  
             }  
             break;  
   
         case FINLEY_REDUCED_ELEMENTS:  
         case FINLEY_ELEMENTS:  
             {  
                 centering = ZONE_CENTERED;  
                 ElementData* cells = mesh->getElements();  
                 id2idxMap = &cells->ID2idx;  
                 reqIDs = &cells->getIDs();  
                 if (cells->reducedCount > 0) {  
                     if (cells->getReducedGhostCount())  
                         reorderArray = &cells->reducedIndexArray;  
                     reorderedNumSamples = cells->reducedCount;  
                     siloMeshName = cells->reducedMesh->getFullSiloName();  
                 } else {  
                     if (cells->getGhostCount())  
                         reorderArray = &cells->indexArray;  
                     reorderedNumSamples = cells->count;  
                     siloMeshName = cells->fullMesh->getFullSiloName();  
                 }  
             }  
             break;  
457    
458          case FINLEY_REDUCED_FACE_ELEMENTS:  //
459          case FINLEY_FACE_ELEMENTS:  //
460              {  //
461                  centering = ZONE_CENTERED;  float* DataVar::getDataFlat() const
462                  ElementData* cells = mesh->getFaceElements();  {
463                  id2idxMap = &cells->ID2idx;      int totalSize = numSamples * getNumberOfComponents();
464                  reqIDs = &cells->getIDs();      float* res = new float[totalSize];
465                  if (cells->reducedCount > 0) {      if (rank == 0) {
466                      if (cells->getReducedGhostCount())          copy(dataArray[0], dataArray[0]+numSamples, res);
467                          reorderArray = &cells->reducedIndexArray;      } else if (rank == 1) {
468                      reorderedNumSamples = cells->reducedCount;          float *dest = res;
469                      siloMeshName = cells->reducedMesh->getFullSiloName();          for (size_t c=0; c<numSamples; c++) {
470                  } else {              for (size_t i=0; i<shape[0]; i++) {
471                      if (cells->getGhostCount())                  *dest++ = dataArray[i][c];
                         reorderArray = &cells->indexArray;  
                     reorderedNumSamples = cells->count;  
                     siloMeshName = cells->fullMesh->getFullSiloName();  
                 }  
472              }              }
473              break;          }
474        } else if (rank == 2) {
475          case FINLEY_REDUCED_CONTACT_ELEMENTS_1:          float *dest = res;
476          case FINLEY_REDUCED_CONTACT_ELEMENTS_2:          for (size_t c=0; c<numSamples; c++) {
477          case FINLEY_CONTACT_ELEMENTS_1:              for (int i=0; i<shape[1]; i++) {
478          case FINLEY_CONTACT_ELEMENTS_2:                  for (int j=0; j<shape[0]; j++) {
479              {                      *dest++ = dataArray[i*shape[0]+j][c];
                 centering = ZONE_CENTERED;  
                 ElementData* cells = mesh->getContactElements();  
                 id2idxMap = &cells->ID2idx;  
                 reqIDs = &cells->getIDs();  
                 if (cells->reducedCount > 0) {  
                     if (cells->getReducedGhostCount())  
                         reorderArray = &cells->reducedIndexArray;  
                     reorderedNumSamples = cells->reducedCount;  
                     siloMeshName = cells->reducedMesh->getFullSiloName();  
                 } else {  
                     if (cells->getGhostCount())  
                         reorderArray = &cells->indexArray;  
                     reorderedNumSamples = cells->count;  
                     siloMeshName = cells->fullMesh->getFullSiloName();  
480                  }                  }
481              }              }
482              break;          }
   
         case FINLEY_POINTS:  
             {  
                 centering = NODE_CENTERED;  
                 ElementData* cells = mesh->getPoints();  
                 if (cells->getGhostCount())  
                     reorderArray = &cells->indexArray;  
                 siloMeshName = cells->fullMesh->getFullSiloName();  
                 id2idxMap = &cells->ID2idx;  
                 reqIDs = &cells->getIDs();  
                 reorderedNumSamples = cells->count;  
             }  
             break;  
   
         default:  
             cerr << "Unknown function space type " << funcSpace << "!\n";  
             return false;  
483      }      }
484    
485      if (reorderedNumSamples > numSamples) {      return res;
486          cerr << "WARNING: " << varName << " has " << numSamples  }
487              << " instead of " << reorderedNumSamples << " samples!" << endl;  
488          return false;  //
489    //
490    //
491    void DataVar::sampleToStream(ostream& os, int index)
492    {
493        // index is -1 for dummy samples, i.e. if writing the full mesh but
494        // only a reduced number of samples is required
495        if (rank == 0) {
496            if (index < 0) {
497                os << 0.;
498            } else {
499                os << dataArray[0][index];
500            }
501        } else if (rank == 1) {
502            if (index < 0) {
503                os << 0. << " " << 0.  << " " << 0.;
504            } else if (shape[0] < 3) {
505                os << dataArray[0][index] << " " << dataArray[1][index]
506                    << " " << 0.;
507            } else {
508                os << dataArray[0][index] << " " << dataArray[1][index]
509                    << " " << dataArray[2][index];
510            }
511        } else if (rank == 2) {
512            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                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      fullMesh = mesh;  //
535    //
536    //
537    void DataVar::writeToVTK(ostream& os, int ownIndex)
538    {
539        if (numSamples == 0)
540            return;
541    
542      reorderSamples(*id2idxMap, *reqIDs);      if (isNodeCentered()) {
543      if (reorderArray)          // data was reordered in reorderSamples() but for VTK we write the
544          handleGhostZones(*reorderArray);          // original node mesh and thus need the original ordering.
545      return true;          // 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            const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
548            const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
549            const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
550            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                    IndexMap::const_iterator it = sampleID2idx.find(requiredIDs[i]);
556                    int idx = (it==sampleID2idx.end() ? -1 : (int)it->second);
557                    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            int toWrite = domain->getElementsByName(meshName)->getNumElements();
564            for (int i=0; i<toWrite; i++) {
565                sampleToStream(os, i);
566            }
567        }
568  }  }
569    
570  ///////////////////////////////  ///////////////////////////////
# Line 599  bool DataVar::setMesh(MeshWithElements* Line 579  bool DataVar::setMesh(MeshWithElements*
579  //  //
580  string DataVar::getTensorDef() const  string DataVar::getTensorDef() const
581  {  {
582      if (rank < 2)      if (rank < 2 || !initialized)
583          return string();          return string();
584            
585      /// Format string for Silo 2x2 tensor      /// Format string for Silo 2x2 tensor
# Line 637  string DataVar::getTensorDef() const Line 617  string DataVar::getTensorDef() const
617    
618  //  //
619  // Writes the data to given Silo file under the virtual path provided.  // Writes the data to given Silo file under the virtual path provided.
620  // The corresponding mesh must have been written already and made known  // The corresponding mesh must have been written already.
 // to this variable by a call to setMesh().  
621  //  //
622  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
623                              const string& units)
624  {  {
625  #if HAVE_SILO  #if USE_SILO
626        if (!initialized)
627            return false;
628    
629      if (numSamples == 0)      if (numSamples == 0)
630          return true;          return true;
631    
     // have to set mesh first  
     if (!fullMesh)  
         return false;  
   
632      int ret;      int ret;
633    
634      if (siloPath != "") {      if (siloPath != "") {
# Line 657  bool DataVar::writeToSilo(DBfile* dbfile Line 636  bool DataVar::writeToSilo(DBfile* dbfile
636          if (ret != 0)          if (ret != 0)
637              return false;              return false;
638      }      }
639        
640      char* meshName = (char*)siloMeshName.c_str();      char* siloMesh = const_cast<char*>(siloMeshName.c_str());
641      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
642        DBoptlist* optList = DBMakeOptlist(2);
643        if (units.length()>0) {
644            DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
645        }
646    
647      if (rank == 0) {      if (rank == 0) {
648          ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName, reorderedData[0],          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
649                  reorderedNumSamples, NULL, 0, DB_FLOAT, dcenter, NULL);                  numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
650      }      }
651      else if (rank == 1) {      else if (rank == 1) {
652          const string comps[3] = {          const string comps[3] = {
# Line 673  bool DataVar::writeToSilo(DBfile* dbfile Line 656  bool DataVar::writeToSilo(DBfile* dbfile
656              comps[0].c_str(), comps[1].c_str(), comps[2].c_str()              comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
657          };          };
658    
659          ret = DBPutUcdvar(dbfile, varName.c_str(), meshName, shape[0],          ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
660                  (char**)varnames, &reorderedData[0], reorderedNumSamples, NULL,                  (char**)varnames, &dataArray[0], numSamples, NULL,
661                  0, DB_FLOAT, dcenter, NULL);                  0, DB_FLOAT, dcenter, optList);
662      }      }
663      else {      else {
664          string tensorDir = varName+string("_comps/");          string tensorDir = varName+string("_comps/");
665          ret = DBMkdir(dbfile, tensorDir.c_str());          ret = DBMkdir(dbfile, tensorDir.c_str());
666          if (ret == 0) {          if (ret == 0) {
667              int one = 1;              int one = 1;
             DBoptlist* optList = DBMakeOptlist(1);  
668              DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);              DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
669    
670              for (int i=0; i<shape[1]; i++) {              for (int i=0; i<shape[1]; i++) {
671                  for (int j=0; j<shape[0]; j++) {                  for (int j=0; j<shape[0]; j++) {
672                      ostringstream varname;                      ostringstream varname;
673                      varname << tensorDir << "a_" << i << j;                      varname << tensorDir << "a_" << i << j;
674                      ret = DBPutUcdvar1(dbfile, varname.str().c_str(), meshName,                      ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
675                              reorderedData[i*shape[0]+j], reorderedNumSamples,                              dataArray[i*shape[0]+j], numSamples,
676                              NULL, 0, DB_FLOAT, dcenter, optList);                              NULL, 0, DB_FLOAT, dcenter, optList);
677                      if (ret != 0) break;                      if (ret != 0) break;
678                  }                  }
679                  if (ret != 0) break;                  if (ret != 0) break;
680              }              }
             DBFreeOptlist(optList);  
681          } // ret==0          } // ret==0
682      } // rank      } // rank
683    
684        DBFreeOptlist(optList);
685      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
686      return (ret == 0);      return (ret == 0);
687    
688  #else // !HAVE_SILO  #else // !USE_SILO
689      return false;      return false;
690  #endif  #endif
691  }  }
692    
693  } // namespace EscriptReader  } // namespace weipa
694    
Legend:
Removed from v.2548  
changed lines
  Added in v.4657
  ViewVC Help
Powered by ViewVC 1.1.26