/[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 2183 by caltinay, Fri Dec 19 03:52:50 2008 UTC trunk/dataexporter/src/DataVar.cpp revision 2888 by caltinay, Fri Jan 29 00:07:00 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
14  //  #include <escriptexport/DataVar.h>
15  // DataVar.cpp  #include <escriptexport/ElementData.h>
16  //  #include <escriptexport/FinleyMesh.h>
17  #include <escriptreader/DataVar.h>  #include <escriptexport/NodeData.h>
18  #include <escriptreader/ElementData.h>  #include <escript/Data.h>
19  #include <escriptreader/MeshWithElements.h>  
20  #include <netcdf.hh>  #if USE_NETCDF
21  #if HAVE_SILO  #include <netcdfcpp.h>
22    #endif
23    
24    #if USE_SILO
25  #include <silo.h>  #include <silo.h>
26  #endif  #endif
27    
28  using namespace std;  using namespace std;
29    
30    namespace escriptexport {
31        
32  enum {  enum {
33      NODE_CENTERED = 1,      NODE_CENTERED = 1,
34      ZONE_CENTERED = 2      ZONE_CENTERED = 2
# Line 33  enum { Line 38  enum {
38  // Constructor  // Constructor
39  //  //
40  DataVar::DataVar(const string& name) :  DataVar::DataVar(const string& name) :
41      varName(name), numSamples(0), rank(0), ptsPerSample(0), centering(0),      initialized(false), varName(name),
42      reorderedNumSamples(0), fullMesh(NULL)      numSamples(0), rank(0), ptsPerSample(0), centering(0)
43  {  {
44  }  }
45    
46  //  //
 // 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;  
 }  
   
 //  
47  // Copy constructor  // Copy constructor
 // Performs a deep copy of the data values  
48  //  //
49  DataVar::DataVar(const DataVar& d) :  DataVar::DataVar(const DataVar& d) :
50        initialized(d.initialized), finleyMesh(d.finleyMesh),
51      varName(d.varName), numSamples(d.numSamples),      varName(d.varName), numSamples(d.numSamples),
52      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),
53      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID),      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)
     reorderedNumSamples(d.reorderedNumSamples), fullMesh(d.fullMesh)  
54  {  {
55      CoordArray::const_iterator it;      if (numSamples > 0) {
56      for (it = d.rawData.begin(); it != d.rawData.end(); it++) {          CoordArray::const_iterator it;
57          float* c = new float[numSamples];          for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
58          copy(*it, (*it)+numSamples, c);              float* c = new float[numSamples];
59          rawData.push_back(c);              copy(*it, (*it)+numSamples, c);
60      }              dataArray.push_back(c);
61      for (it = d.reorderedData.begin(); it != d.reorderedData.end(); it++) {          }
         float* c = new float[reorderedNumSamples];  
         copy(*it, (*it)+reorderedNumSamples, c);  
         reorderedData.push_back(c);  
62      }      }
63  }  }
64    
65  //  //
66  // Special constructor for mesh data  // Destructor
 // Used to store "special" integral mesh variables such as IDs or tags  
67  //  //
68  DataVar::DataVar(const string& name, const IntVec& data,  DataVar::~DataVar()
                  MeshWithElements* mesh) :  
     varName(name)  
69  {  {
70      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;  
71  }  }
72    
73  //  //
 // Appends raw data including IDs from rhs. Reordered data becomes invalid.  
74  //  //
75  bool DataVar::append(const DataVar& rhs)  //
76    void DataVar::cleanup()
77  {  {
     // 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  
78      CoordArray::iterator it;      CoordArray::iterator it;
79      for (it = reorderedData.begin(); it != reorderedData.end(); it++)      for (it = dataArray.begin(); it != dataArray.end(); it++)
80          delete[] *it;          delete[] *it;
81      reorderedData.clear();      dataArray.clear();
82      reorderedNumSamples = 0;      shape.clear();
83            sampleID.clear();
84      return true;      numSamples = 0;
85        initialized = false;
86  }  }
87    
88  //  //
 // 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,...)  
89  //  //
90  float* DataVar::averageData(const float* src, size_t stride)  //
91    bool DataVar::initFromEscript(escript::Data& escriptData, FinleyMesh_ptr mesh)
92  {  {
93      float* res = new float[numSamples];      cleanup();
94    
95      if (ptsPerSample == 1) {      if (!escriptData.actsExpanded()) {
96          float* dest = res;          cerr << "WARNING: Only expanded data supported!" << endl;
97          for (int i=0; i<numSamples; i++, src+=stride)          return false;
98              *dest++ = *src;      }
99    
100        finleyMesh = mesh;
101        rank = escriptData.getDataPointRank();
102        ptsPerSample = escriptData.getNumDataPointsPerSample();
103        shape = escriptData.getDataPointShape();
104        funcSpace = escriptData.getFunctionSpace().getTypeCode();
105        numSamples = escriptData.getNumSamples();
106    
107        if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
108            centering = NODE_CENTERED;
109      } else {      } else {
110          float* dest = res;          centering = ZONE_CENTERED;
         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);  
         }  
111      }      }
     return res;  
 }  
112    
113  //  #ifdef _DEBUG
114  // Reads scalar data (rank=0) from NetCDF file and stores the values      cout << varName << ":\t" << numSamples << " samples,  "
115  // after averaging.          << ptsPerSample << " pts/s,  rank: " << rank << endl;
116  //  #endif
 void DataVar::readRank0Data(NcFile* ncfile)  
 {  
     shape.clear();  
     float* tempData = new float[ptsPerSample*numSamples];  
     NcVar* var = ncfile->get_var("data");  
     var->get(tempData, ptsPerSample, numSamples);  
117    
118      float* c = averageData(tempData, 1);      NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
119      rawData.push_back(c);      if (nodes == NULL)
120            return false;
121    
122      delete[] tempData;      meshName = nodes->getName();
123  }      siloMeshName = nodes->getFullSiloName();
124        initialized = true;
125    
126  //      // no samples? Nothing more to do.
127  // Reads vector data (rank=1) from NetCDF file and stores the components      if (numSamples == 0)
128  // separately after averaging.          return true;
129  //  
130  void DataVar::readRank1Data(NcFile* ncfile)      const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
131  {      sampleID.insert(sampleID.end(), numSamples, 0);
132      shape.clear();      copy(iPtr, iPtr+numSamples, sampleID.begin());
     NcDim* dim = ncfile->get_dim("d0");  
     shape.push_back(dim->size());  
133    
134      float* tempData = new float[shape[0]*ptsPerSample*numSamples];      size_t dimSize = 1;
135      NcVar* var = ncfile->get_var("data");      if (rank > 0)
136      var->get(tempData, shape[0], ptsPerSample, numSamples);          dimSize *= shape[0];
137        if (rank > 1)
138            dimSize *= shape[1];
139        if (rank > 2) {
140            cerr << "WARNING: Rank " << rank << " data is not supported!\n";
141            initialized = false;
142        }
143    
144        if (initialized) {
145            size_t dataSize = dimSize * ptsPerSample;
146            float* tempData = new float[dataSize*numSamples];
147            float* destPtr = tempData;
148            for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
149                const escript::DataAbstract::ValueType::value_type* values =
150                    escriptData.getSampleDataRO(sampleNo);
151                copy(values, values+dataSize, destPtr);
152                destPtr += dataSize;
153            }
154    
155            const float* srcPtr = tempData;
156            for (int i=0; i < dimSize; i++, srcPtr++) {
157                float* c = averageData(srcPtr, dimSize);
158                dataArray.push_back(c);
159            }
160            delete[] tempData;
161    
162      for (int i=0; i<shape[0]; i++) {          initialized = reorderSamples();
         const float* src = tempData;  
         src+=i;  
         float* c = averageData(src, shape[0]);  
         rawData.push_back(c);  
163      }      }
164      delete[] tempData;  
165        return initialized;
166  }  }
167    
168  //  //
169  // Like readRank1Data() but for tensor data (rank=2).  // Initialise with mesh data
170  //  //
171  void DataVar::readRank2Data(NcFile* ncfile)  bool DataVar::initFromMesh(FinleyMesh_ptr mesh)
172  {  {
173      shape.clear();      cleanup();
174      NcDim* dim = ncfile->get_dim("d0");      
175      shape.push_back(dim->size());      finleyMesh = mesh;
176      dim = ncfile->get_dim("d1");      rank = 0;
177      shape.push_back(dim->size());      ptsPerSample = 1;
178        NodeData_ptr nodes;
179      float* tempData = new float[shape[0]*shape[1]*ptsPerSample*numSamples];  
180      NcVar* var = ncfile->get_var("data");      if (varName.find("ContactElements_") != varName.npos) {
181      var->get(tempData, shape[0], shape[1], ptsPerSample, numSamples);          funcSpace = FINLEY_CONTACT_ELEMENTS_1;
182            centering = ZONE_CENTERED;
183      for (int i=0; i<shape[1]; i++) {          string elementName = varName.substr(0, varName.find('_'));
184          for (int j=0; j<shape[0]; j++) {          ElementData_ptr elements = mesh->getElementsByName(elementName);
185              const float* src = tempData;          nodes = elements->getNodeMesh();
186              src+=i*shape[0]+j;          sampleID = elements->getIDs();
187              float* c = averageData(src, shape[0]*shape[1]);      } else if (varName.find("FaceElements_") != varName.npos) {
188              rawData.push_back(c);          funcSpace = FINLEY_FACE_ELEMENTS;
189          }          centering = ZONE_CENTERED;
190            string elementName = varName.substr(0, varName.find('_'));
191            ElementData_ptr elements = mesh->getElementsByName(elementName);
192            nodes = elements->getNodeMesh();
193            sampleID = elements->getIDs();
194        } else if (varName.find("Elements_") != varName.npos) {
195            funcSpace = FINLEY_ELEMENTS;
196            centering = ZONE_CENTERED;
197            string elementName = varName.substr(0, varName.find('_'));
198            ElementData_ptr elements = mesh->getElementsByName(elementName);
199            nodes = elements->getNodeMesh();
200            sampleID = elements->getIDs();
201        } else if (varName.find("Nodes_") != varName.npos) {
202            funcSpace = FINLEY_NODES;
203            centering = NODE_CENTERED;
204            nodes = mesh->getNodes();
205            sampleID = nodes->getNodeIDs();
206        } else {
207            cerr << "WARNING: Unrecognized mesh variable '" << varName << "'\n";
208            return false;
209        }
210    
211        meshName = nodes->getName();
212        siloMeshName = nodes->getFullSiloName();
213    
214        const IntVec& data = mesh->getVarDataByName(varName);
215        numSamples = data.size();
216    
217        if (numSamples > 0) {
218            float* c = new float[numSamples];
219            dataArray.push_back(c);
220            IntVec::const_iterator it;
221            for (it=data.begin(); it != data.end(); it++)
222                *c++ = static_cast<float>(*it);
223      }      }
224      delete[] tempData;      initialized = true;
225    
226        return initialized;
227  }  }
228    
229  //  //
230  // Reads a NetCDF file in escript/finley format.  // Reads variable data from NetCDF file
231  //  //
232  bool DataVar::readFromNc(const string& ncFile)  bool DataVar::initFromNetCDF(const string& filename, FinleyMesh_ptr mesh)
233  {  {
234        cleanup();
235        
236    #if USE_NETCDF
237      NcError ncerr(NcError::silent_nonfatal);          NcError ncerr(NcError::silent_nonfatal);    
238      NcFile* input = new NcFile(ncFile.c_str());      NcFile* input = new NcFile(filename.c_str());
239      if (!input->is_valid()) {      if (!input->is_valid()) {
240          cerr << "Could not open input file " << ncFile.c_str() << "." << endl;          cerr << "Could not open input file " << filename << "." << endl;
241          delete input;          delete input;
242          return false;          return false;
243      }      }
# Line 267  bool DataVar::readFromNc(const string& n Line 245  bool DataVar::readFromNc(const string& n
245      NcDim* dim;      NcDim* dim;
246      NcAtt* att;      NcAtt* att;
247    
     dim = input->get_dim("num_samples");  
     numSamples = dim->size();  
   
     // if there are no data samples bail out  
     if (numSamples == 0) {  
         delete input;  
         return false;  
     }  
   
248      att = input->get_att("type_id");      att = input->get_att("type_id");
249      int typeID = att->as_int(0);      int typeID = att->as_int(0);
250      if (typeID != 2) {      if (typeID != 2) {
251          cerr << "WARNING: Only expanded data supported at the moment!" << endl;          cerr << "WARNING: Only expanded data supported!" << endl;
252          delete input;          delete input;
253          return false;          return false;
254      }      }
# Line 293  bool DataVar::readFromNc(const string& n Line 262  bool DataVar::readFromNc(const string& n
262      att = input->get_att("function_space_type");      att = input->get_att("function_space_type");
263      funcSpace = att->as_int(0);      funcSpace = att->as_int(0);
264    
265        if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
266            centering = NODE_CENTERED;
267        } else {
268            centering = ZONE_CENTERED;
269        }
270    
271        dim = input->get_dim("num_samples");
272        numSamples = dim->size();
273    
274  #ifdef _DEBUG  #ifdef _DEBUG
275      cout << varName << ":\t" << numSamples << " samples,  "      cout << varName << ":\t" << numSamples << " samples,  "
276          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
277  #endif  #endif
278    
279      sampleID.clear();      finleyMesh = mesh;
280      sampleID.insert(sampleID.end(), numSamples, 0);      NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
281      NcVar* var = input->get_var("id");      if (nodes == NULL) {
282      var->get(&sampleID[0], numSamples);          delete input;
283            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;  
284      }      }
285    
286      delete input;      meshName = nodes->getName();
287      return true;      siloMeshName = nodes->getFullSiloName();
288  }      initialized = true;
289    
290  //      size_t dimSize = 1;
291  // Returns one of the mesh names provided by mainMesh that matches the      vector<long> counts;
292  // data variable's function space type and reduced/unreduced state.  
293  //      if (rank > 0) {
294  string DataVar::getMeshName(MeshWithElements* mainMesh) const          dim = input->get_dim("d0");
295  {          int d = dim->size();
296      string name;          shape.push_back(d);
297            counts.push_back(d);
298      switch (funcSpace) {          dimSize *= d;
299          case FINLEY_REDUCED_NODES:      }
300          case FINLEY_REDUCED_DEGREES_OF_FREEDOM:      if (rank > 1) {
301          case FINLEY_REDUCED_ELEMENTS:          dim = input->get_dim("d1");
302          case FINLEY_ELEMENTS:          int d = dim->size();
303              if (mainMesh->getElements()->reducedCount > 0) {          shape.push_back(d);
304                  name = mainMesh->getElements()->reducedMesh->getName();          counts.push_back(d);
305              } else {          dimSize *= d;
306                  name = mainMesh->getElements()->fullMesh->getName();      }
307              }      if (rank > 2) {
308              break;          cerr << "WARNING: Rank " << rank << " data is not supported!\n";
309            initialized = false;
310        }
311    
312        if (initialized && numSamples > 0) {
313            sampleID.insert(sampleID.end(), numSamples, 0);
314            NcVar* var = input->get_var("id");
315            var->get(&sampleID[0], numSamples);
316    
317            size_t dataSize = dimSize*numSamples*ptsPerSample;
318            counts.push_back(ptsPerSample);
319            counts.push_back(numSamples);
320            float* tempData = new float[dataSize];
321            var = input->get_var("data");
322            var->get(tempData, &counts[0]);
323    
324            const float* srcPtr = tempData;
325            for (int i=0; i < dimSize; i++, srcPtr++) {
326                float* c = averageData(srcPtr, dimSize);
327                dataArray.push_back(c);
328            }
329            delete[] tempData;
330    
331          case FINLEY_REDUCED_FACE_ELEMENTS:          initialized = reorderSamples();
332          case FINLEY_FACE_ELEMENTS:      }
             if (mainMesh->getFaceElements()->reducedCount > 0) {  
                 name = mainMesh->getFaceElements()->reducedMesh->getName();  
             } else {  
                 name = mainMesh->getFaceElements()->fullMesh->getName();  
             }  
             break;  
333    
334          case FINLEY_REDUCED_CONTACT_ELEMENTS_1:      delete input;
335          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;  
336    
337          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;  
338  }  }
339    
340  //  //
# Line 380  string DataVar::getMeshName(MeshWithElem Line 342  string DataVar::getMeshName(MeshWithElem
342  //  //
343  bool DataVar::isNodeCentered() const  bool DataVar::isNodeCentered() const
344  {  {
345      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);  
346  }  }
347    
348  //  //
349  // Filters and reorders the raw sample values according to the IDs provided  // Returns a subset of the src array according to stride parameter.
350  // in 'requiredIDs'. This is used to have data arrays ordered according to  // If samples consist of multiple values they are averaged beforehand.
351  // the underlying mesh (i.e. DataID[i]==MeshNodeID[i])  // Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
352    // (z0,z1,...)
353  //  //
354  void DataVar::reorderSamples(const IndexMap& id2idxMap,  float* DataVar::averageData(const float* src, size_t stride)
                              const IntVec& requiredIDs)  
355  {  {
356      CoordArray::iterator it;      float* res;
     for (it = reorderedData.begin(); it != reorderedData.end(); it++)  
         delete[] *it;  
     reorderedData.clear();  
357    
358      buildIndexMap();      if (ptsPerSample == 1) {
359      for (size_t i=0; i < rawData.size(); i++) {          res = new float[numSamples];
360          float* c = new float[reorderedNumSamples];          float* dest = res;
361          reorderedData.push_back(c);          for (int i=0; i<numSamples; i++, src+=stride)
362          const float* src = rawData[i];              *dest++ = *src;
363          IntVec::const_iterator idIt = requiredIDs.begin();      } else {
364          for (; idIt != requiredIDs.end(); idIt++) {          ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
365              size_t srcIdx = sampleID2idx.find(*idIt)->second;          int cellFactor = cells->getElementFactor();
366              size_t destIdx = id2idxMap.find(*idIt)->second;          res = new float[cellFactor * numSamples];
367              c[destIdx] = src[srcIdx];          float* dest = res;
368            QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
369            if (qmi.mask.size() > 0) {
370                const float* tmpSrc = src;
371                for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
372                    for (int l=0; l<cellFactor; l++) {
373                        double tmpVal = 0.0;
374                        for (int j=0; j<ptsPerSample; j++) {
375                            if (qmi.mask[l][j] != 0) {
376                                tmpVal += *(tmpSrc+stride*j);
377                            }
378                        }
379                        *dest++ = (float)(tmpVal / qmi.factor[l]);
380                    }
381                }
382            } else {
383                for (int i=0; i<numSamples; i++) {
384                    double tmpVal = 0.0;
385                    for (int j=0; j<ptsPerSample; j++, src+=stride) {
386                        tmpVal += *src;
387                    }
388                    tmpVal /= ptsPerSample;
389                    for (int l=0; l<cellFactor; l++) {
390                        *dest++ = static_cast<float>(tmpVal);
391                    }
392                }
393          }          }
394      }      }
395        return res;
396  }  }
397    
398  //  //
399  // For zonal data this method reorders the values according to the indices  // Filters and reorders the raw sample values according to the node/element
400  // 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
401  // the arrays which conforms to Silo's expected format.  // mesh (i.e. DataID[i]==MeshNodeID[i])
 // Nodal data is not changed by this method.  
 //  
 void DataVar::handleGhostZones(const IntVec& reorderArray)  
 {  
     if (centering == NODE_CENTERED)  
         return;  
   
     vector<float*>::iterator it;  
     for (it = reorderedData.begin(); it!=reorderedData.end(); it++) {  
         float* original = *it;  
         float* reordered = new float[reorderedNumSamples];  
         float* arrIt = reordered;  
         IntVec::const_iterator idxIt;  
         for (idxIt=reorderArray.begin(); idxIt!=reorderArray.end(); idxIt++)  
             *arrIt++ = original[*idxIt];  
   
         delete[] *it;  
         *it = reordered;  
     }  
 }  
   
 //  
 // 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.  
402  //  //
403  bool DataVar::setMesh(MeshWithElements* mesh)  bool DataVar::reorderSamples()
404  {  {
405      if (fullMesh == mesh)      if (numSamples == 0)
406          return true;          return true;
407    
408      const IndexMap* id2idxMap;      const IntVec* requiredIDs = NULL;
409      const IntVec* reqIDs;      int requiredNumSamples = 0;
410      const IntVec* reorderArray = NULL;      int cellFactor = 1;
411    
412      switch (funcSpace) {      if (centering == NODE_CENTERED) {
413          case FINLEY_REDUCED_NODES:          NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
414          case FINLEY_REDUCED_DEGREES_OF_FREEDOM:          requiredIDs = &nodes->getNodeIDs();
415              {          requiredNumSamples = nodes->getNumNodes();
416                  centering = NODE_CENTERED;      } else {
417                  ElementData* cells = mesh->getElements();          ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
418                  if (cells->reducedCount > 0) {          if (cells == NULL)
419                      if (cells->getReducedGhostCount())              return false;
                         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;  
420    
421          case FINLEY_NODES:          requiredIDs = &cells->getIDs();
422          case FINLEY_DEGREES_OF_FREEDOM:          requiredNumSamples = cells->getNumElements();
423              {          cellFactor = cells->getElementFactor();
424                  centering = NODE_CENTERED;          if (cellFactor > 1) {
425                  ElementData* cells = mesh->getElements();              numSamples *= cellFactor;
426                  if (cells->getGhostCount())              // update sample IDs
427                      reorderArray = &cells->indexArray;              IntVec newSampleID(numSamples);
428                  siloMeshName = cells->fullMesh->getFullSiloName();              IntVec::const_iterator idIt = sampleID.begin();
429                  id2idxMap = &cells->fullMesh->getIndexMap();              IntVec::iterator newIDit = newSampleID.begin();
430                  reqIDs = &cells->fullMesh->getNodeIDs();              for (; idIt != sampleID.end(); idIt++, newIDit+=cellFactor) {
431                  reorderedNumSamples = cells->fullMesh->getNumNodes();                  fill(newIDit, newIDit+cellFactor, *idIt);
432              }              }
433              break;              sampleID.swap(newSampleID);
434            }
435        }
436    
437          case FINLEY_REDUCED_ELEMENTS:      if (requiredNumSamples > numSamples) {
438          case FINLEY_ELEMENTS:          cerr << "ERROR: " << varName << " has " << numSamples
439              {              << " instead of " << requiredNumSamples << " samples!" << endl;
440                  centering = ZONE_CENTERED;          return false;
441                  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;  
442    
443          case FINLEY_REDUCED_FACE_ELEMENTS:      IndexMap sampleID2idx = buildIndexMap();
444          case FINLEY_FACE_ELEMENTS:      numSamples = requiredNumSamples;
             {  
                 centering = ZONE_CENTERED;  
                 ElementData* cells = mesh->getFaceElements();  
                 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;  
445    
446          case FINLEY_REDUCED_CONTACT_ELEMENTS_1:      // now filter the data
447          case FINLEY_REDUCED_CONTACT_ELEMENTS_2:      for (size_t i=0; i < dataArray.size(); i++) {
448          case FINLEY_CONTACT_ELEMENTS_1:          float* c = new float[numSamples];
449          case FINLEY_CONTACT_ELEMENTS_2:          const float* src = dataArray[i];
450              {          IntVec::const_iterator idIt = requiredIDs->begin();
451                  centering = ZONE_CENTERED;          size_t destIdx = 0;
452                  ElementData* cells = mesh->getContactElements();          for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
453                  id2idxMap = &cells->ID2idx;              size_t srcIdx = sampleID2idx.find(*idIt)->second;
454                  reqIDs = &cells->getIDs();              copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
455                  if (cells->reducedCount > 0) {          }
456                      if (cells->getReducedGhostCount())          delete[] dataArray[i];
457                          reorderArray = &cells->reducedIndexArray;          dataArray[i] = c;
458                      reorderedNumSamples = cells->reducedCount;      }
                     siloMeshName = cells->reducedMesh->getFullSiloName();  
                 } else {  
                     if (cells->getGhostCount())  
                         reorderArray = &cells->indexArray;  
                     reorderedNumSamples = cells->count;  
                     siloMeshName = cells->fullMesh->getFullSiloName();  
                 }  
             }  
             break;  
459    
460          case FINLEY_POINTS:      // sample IDs now = mesh node/element IDs
461              {      sampleID = *requiredIDs;
                 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;  
462    
463          default:      return true;
464              cerr << "Unknown function space type " << funcSpace << "!\n";  }
             return false;  
     }  
465    
466      if (reorderedNumSamples > numSamples) {  //
467          cerr << "WARNING: " << varName << " has " << numSamples  //
468              << " instead of " << reorderedNumSamples << " samples!" << endl;  //
469          return false;  void DataVar::sampleToStream(ostream& os, int index)
470    {
471        if (rank == 0) {
472            os << dataArray[0][index];
473        } else if (rank == 1) {
474            if (shape[0] < 3)
475                os << dataArray[0][index] << " " << dataArray[1][index]
476                    << " " << 0.;
477            else
478                os << dataArray[0][index] << " " << dataArray[1][index]
479                    << " " << dataArray[2][index];
480        } else if (rank == 2) {
481            if (shape[1] < 3) {
482                os << dataArray[0][index] << " " << dataArray[1][index]
483                    << " " << 0. << " ";
484                os << dataArray[2][index] << " " << dataArray[3][index]
485                    << " " << 0. << " ";
486                os << 0. << " " << 0. << " " << 0.;
487            } else {
488                os << dataArray[0][index] << " " << dataArray[1][index]
489                    << " " << dataArray[2][index] << " ";
490                os << dataArray[3][index] << " " << dataArray[4][index]
491                    << " " << dataArray[5][index] << " ";
492                os << dataArray[6][index] << " " << dataArray[7][index]
493                    << " " << dataArray[8][index];
494            }
495      }      }
496        os << endl;
497    }
498    
499      fullMesh = mesh;  //
500    //
501    //
502    void DataVar::writeToVTK(ostream& os, int ownIndex)
503    {
504        if (numSamples == 0)
505            return;
506    
507      reorderSamples(*id2idxMap, *reqIDs);      if (isNodeCentered()) {
508      if (reorderArray)          // data was reordered in reorderSamples() but for VTK we write the
509          handleGhostZones(*reorderArray);          // original node mesh and thus need the original ordering...
510      return true;          const IntVec& requiredIDs = finleyMesh->getNodes()->getNodeIDs();
511            const IntVec& nodeGNI = finleyMesh->getNodes()->getGlobalNodeIndices();
512            const IntVec& nodeDist = finleyMesh->getNodes()->getNodeDistribution();
513            int firstId = nodeDist[ownIndex];
514            int lastId = nodeDist[ownIndex+1];
515            IndexMap sampleID2idx = buildIndexMap();
516            for (int i=0; i<nodeGNI.size(); i++) {
517                if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
518                    int idx = sampleID2idx[requiredIDs[i]];
519                    sampleToStream(os, idx);
520                }
521            }
522        } else {
523            // cell data: ghost cells have been removed so do not write ghost
524            // samples (which are the last elements in the arrays)
525            int toWrite =
526                finleyMesh->getElementsByName(meshName)->getNumElements();
527            for (int i=0; i<toWrite; i++) {
528                sampleToStream(os, i);
529            }
530        }
531  }  }
532    
533  ///////////////////////////////  ///////////////////////////////
# Line 599  bool DataVar::setMesh(MeshWithElements* Line 542  bool DataVar::setMesh(MeshWithElements*
542  //  //
543  string DataVar::getTensorDef() const  string DataVar::getTensorDef() const
544  {  {
545      if (rank < 2)      if (rank < 2 || !initialized)
546          return string();          return string();
547            
548      /// Format string for Silo 2x2 tensor      /// Format string for Silo 2x2 tensor
# Line 642  string DataVar::getTensorDef() const Line 585  string DataVar::getTensorDef() const
585  //  //
586  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)
587  {  {
588  #if HAVE_SILO  #if USE_SILO
589        if (!initialized)
590            return false;
591    
592      if (numSamples == 0)      if (numSamples == 0)
593          return true;          return true;
594    
     // have to set mesh first  
     if (!fullMesh)  
         return false;  
   
595      int ret;      int ret;
596    
597      if (siloPath != "") {      if (siloPath != "") {
# Line 657  bool DataVar::writeToSilo(DBfile* dbfile Line 599  bool DataVar::writeToSilo(DBfile* dbfile
599          if (ret != 0)          if (ret != 0)
600              return false;              return false;
601      }      }
602        
603      char* meshName = (char*)siloMeshName.c_str();      char* siloMesh = const_cast<char*>(siloMeshName.c_str());
604      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
605    
606      if (rank == 0) {      if (rank == 0) {
607          ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName, reorderedData[0],          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
608                  reorderedNumSamples, NULL, 0, DB_FLOAT, dcenter, NULL);                  numSamples, NULL, 0, DB_FLOAT, dcenter, NULL);
609      }      }
610      else if (rank == 1) {      else if (rank == 1) {
611          const string comps[3] = {          const string comps[3] = {
# Line 673  bool DataVar::writeToSilo(DBfile* dbfile Line 615  bool DataVar::writeToSilo(DBfile* dbfile
615              comps[0].c_str(), comps[1].c_str(), comps[2].c_str()              comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
616          };          };
617    
618          ret = DBPutUcdvar(dbfile, varName.c_str(), meshName, shape[0],          ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
619                  (char**)varnames, &reorderedData[0], reorderedNumSamples, NULL,                  (char**)varnames, &dataArray[0], numSamples, NULL,
620                  0, DB_FLOAT, dcenter, NULL);                  0, DB_FLOAT, dcenter, NULL);
621      }      }
622      else {      else {
# Line 689  bool DataVar::writeToSilo(DBfile* dbfile Line 631  bool DataVar::writeToSilo(DBfile* dbfile
631                  for (int j=0; j<shape[0]; j++) {                  for (int j=0; j<shape[0]; j++) {
632                      ostringstream varname;                      ostringstream varname;
633                      varname << tensorDir << "a_" << i << j;                      varname << tensorDir << "a_" << i << j;
634                      ret = DBPutUcdvar1(dbfile, varname.str().c_str(), meshName,                      ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
635                              reorderedData[i*shape[0]+j], reorderedNumSamples,                              dataArray[i*shape[0]+j], numSamples,
636                              NULL, 0, DB_FLOAT, dcenter, optList);                              NULL, 0, DB_FLOAT, dcenter, optList);
637                      if (ret != 0) break;                      if (ret != 0) break;
638                  }                  }
# Line 703  bool DataVar::writeToSilo(DBfile* dbfile Line 645  bool DataVar::writeToSilo(DBfile* dbfile
645      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
646      return (ret == 0);      return (ret == 0);
647    
648  #else // !HAVE_SILO  #else // !USE_SILO
649      return false;      return false;
650  #endif  #endif
651  }  }
652    
653    } // namespace escriptexport
654    

Legend:
Removed from v.2183  
changed lines
  Added in v.2888

  ViewVC Help
Powered by ViewVC 1.1.26