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

Diff of /trunk/dataexporter/src/DataVar.cpp

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

trunk/tools/libescriptreader/src/escriptreader/DataVar.cpp revision 2806 by jfenwick, Mon Jul 20 06:20:06 2009 UTC trunk/dataexporter/src/DataVar.cpp revision 2810 by caltinay, Mon Dec 7 04:13:49 2009 UTC
# 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    #if USE_NETCDF
21  #include <netcdf.hh>  #include <netcdf.hh>
22  #if HAVE_SILO  #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 EscriptReader {  namespace escriptexport {
31            
32  enum {  enum {
33      NODE_CENTERED = 1,      NODE_CENTERED = 1,
# Line 35  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
48  //  //
49  DataVar::DataVar(const DataVar& d) :  DataVar::DataVar(const DataVar& d) :
50      varName(d.varName), numSamples(d.numSamples),      varName(d.varName), numSamples(d.numSamples),
51      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),
52      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID),      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)
     reorderedNumSamples(d.reorderedNumSamples), fullMesh(d.fullMesh)  
53  {  {
54      CoordArray::const_iterator it;      if (numSamples > 0) {
55      for (it = d.rawData.begin(); it != d.rawData.end(); it++) {          CoordArray::const_iterator it;
56          float* c = new float[numSamples];          for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
57          copy(*it, (*it)+numSamples, c);              float* c = new float[numSamples];
58          rawData.push_back(c);              copy(*it, (*it)+numSamples, c);
59      }              dataArray.push_back(c);
60      for (it = d.reorderedData.begin(); it != d.reorderedData.end(); it++) {          }
         float* c = new float[reorderedNumSamples];  
         copy(*it, (*it)+reorderedNumSamples, c);  
         reorderedData.push_back(c);  
61      }      }
62        initialized = d.initialized;
63  }  }
64    
65  //  //
66  // Special constructor for mesh data  // Destructor
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.  
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        rank = escriptData.getDataPointRank();
101        ptsPerSample = escriptData.getNumDataPointsPerSample();
102        shape = escriptData.getDataPointShape();
103        funcSpace = escriptData.getFunctionSpace().getTypeCode();
104        numSamples = escriptData.getNumSamples();
105    
106        if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
107            centering = NODE_CENTERED;
108      } else {      } else {
109          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);  
         }  
110      }      }
     return res;  
 }  
111    
112  //  #ifdef _DEBUG
113  // Reads scalar data (rank=0) from NetCDF file and stores the values      cout << varName << ":\t" << numSamples << " samples,  "
114  // after averaging.          << ptsPerSample << " pts/s,  rank: " << rank << endl;
115  //  #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);  
116    
117      float* c = averageData(tempData, 1);      initialized = true;
     rawData.push_back(c);  
118    
119      delete[] tempData;      if (numSamples == 0)
120  }          return true;
121    
122  //      const int* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
123  // Reads vector data (rank=1) from NetCDF file and stores the components      sampleID.insert(sampleID.end(), numSamples, 0);
124  // separately after averaging.      copy(iPtr, iPtr+numSamples, sampleID.begin());
125  //  
126  void DataVar::readRank1Data(NcFile* ncfile)      size_t dimSize = 1;
127  {      if (rank > 0)
128      shape.clear();          dimSize *= shape[0];
129      NcDim* dim = ncfile->get_dim("d0");      if (rank > 1)
130      shape.push_back(dim->size());          dimSize *= shape[1];
131        if (rank > 2) {
132            cerr << "WARNING: Rank " << rank << " data is not supported!\n";
133            initialized = false;
134        }
135    
136        if (initialized) {
137            size_t dataSize = dimSize * ptsPerSample;
138            float* tempData = new float[dataSize*numSamples];
139            float* destPtr = tempData;
140            for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
141                const escript::DataAbstract::ValueType::value_type* values =
142                    escriptData.getSampleDataRO(sampleNo);
143                copy(values, values+dataSize, destPtr);
144                destPtr += dataSize;
145            }
146    
147      float* tempData = new float[shape[0]*ptsPerSample*numSamples];          const float* srcPtr = tempData;
148      NcVar* var = ncfile->get_var("data");          for (int i=0; i < dimSize; i++, srcPtr++) {
149      var->get(tempData, shape[0], ptsPerSample, numSamples);              float* c = averageData(srcPtr, dimSize);
150                dataArray.push_back(c);
151            }
152            delete[] tempData;
153    
154      for (int i=0; i<shape[0]; i++) {          initialized = filterSamples(mesh);
         const float* src = tempData;  
         src+=i;  
         float* c = averageData(src, shape[0]);  
         rawData.push_back(c);  
155      }      }
156      delete[] tempData;  
157        return initialized;
158  }  }
159    
160  //  //
161  // Like readRank1Data() but for tensor data (rank=2).  // Initialise with mesh data
162  //  //
163  void DataVar::readRank2Data(NcFile* ncfile)  bool DataVar::initFromMesh(FinleyMesh_ptr mesh)
164  {  {
165      shape.clear();      cleanup();
166      NcDim* dim = ncfile->get_dim("d0");      
167      shape.push_back(dim->size());      const IntVec& data = mesh->getVarDataByName(varName);
168      dim = ncfile->get_dim("d1");      rank = 0;
169      shape.push_back(dim->size());      ptsPerSample = 1;
170        numSamples = data.size();
171      float* tempData = new float[shape[0]*shape[1]*ptsPerSample*numSamples];  
172      NcVar* var = ncfile->get_var("data");      if (numSamples > 0) {
173      var->get(tempData, shape[0], shape[1], ptsPerSample, numSamples);          float* c = new float[numSamples];
174            dataArray.push_back(c);
175      for (int i=0; i<shape[1]; i++) {          IntVec::const_iterator it;
176          for (int j=0; j<shape[0]; j++) {          for (it=data.begin(); it != data.end(); it++)
177              const float* src = tempData;              *c++ = static_cast<float>(*it);
178              src+=i*shape[0]+j;  
179              float* c = averageData(src, shape[0]*shape[1]);          if (varName.compare(0, 6, "Nodes_") == 0) {
180              rawData.push_back(c);              funcSpace = FINLEY_NODES;
181                centering = NODE_CENTERED;
182                sampleID.insert(sampleID.end(),
183                        mesh->getNodes()->getNodeIDs().begin(),
184                        mesh->getNodes()->getNodeIDs().end());
185            } else if (varName.compare(0, 9, "Elements_") == 0) {
186                funcSpace = FINLEY_ELEMENTS;
187                centering = ZONE_CENTERED;
188                sampleID.insert(sampleID.end(),
189                        mesh->getElements()->getIDs().begin(),
190                        mesh->getElements()->getIDs().end());
191            } else if (varName.compare(0, 13, "FaceElements_") == 0) {
192                funcSpace = FINLEY_FACE_ELEMENTS;
193                centering = ZONE_CENTERED;
194                sampleID.insert(sampleID.end(),
195                        mesh->getFaceElements()->getIDs().begin(),
196                        mesh->getFaceElements()->getIDs().end());
197            } else if (varName.compare(0, 16, "ContactElements_") == 0) {
198                funcSpace = FINLEY_CONTACT_ELEMENTS_1;
199                centering = ZONE_CENTERED;
200                sampleID.insert(sampleID.end(),
201                        mesh->getContactElements()->getIDs().begin(),
202                        mesh->getContactElements()->getIDs().end());
203            } else {
204                return false;
205          }          }
206    
207            NodeData_ptr nodes = mesh->getMeshForFinleyFS(funcSpace);
208            meshName = nodes->getName();
209            siloMeshName = nodes->getFullSiloName();
210      }      }
211      delete[] tempData;      initialized = true;
212    
213        return initialized;
214  }  }
215    
216  //  //
217  // Reads a NetCDF file in escript/finley format.  // Reads variable data from NetCDF file
218  //  //
219  bool DataVar::readFromNc(const string& filename)  bool DataVar::initFromNetCDF(const string& filename, FinleyMesh_ptr mesh)
220  {  {
221        cleanup();
222        
223    #if USE_NETCDF
224      NcError ncerr(NcError::silent_nonfatal);          NcError ncerr(NcError::silent_nonfatal);    
225      NcFile* input = new NcFile(filename.c_str());      NcFile* input = new NcFile(filename.c_str());
226      if (!input->is_valid()) {      if (!input->is_valid()) {
# Line 267  bool DataVar::readFromNc(const string& f Line 232  bool DataVar::readFromNc(const string& f
232      NcDim* dim;      NcDim* dim;
233      NcAtt* att;      NcAtt* att;
234    
     dim = input->get_dim("num_samples");  
     numSamples = dim->size();  
   
     // if there are no data samples bail out  
     if (numSamples == 0) {  
         delete input;  
         return false;  
     }  
   
235      att = input->get_att("type_id");      att = input->get_att("type_id");
236      int typeID = att->as_int(0);      int typeID = att->as_int(0);
237      if (typeID != 2) {      if (typeID != 2) {
238          cerr << "WARNING: Only expanded data supported at the moment!" << endl;          cerr << "WARNING: Only expanded data supported!" << endl;
239          delete input;          delete input;
240          return false;          return false;
241      }      }
# Line 293  bool DataVar::readFromNc(const string& f Line 249  bool DataVar::readFromNc(const string& f
249      att = input->get_att("function_space_type");      att = input->get_att("function_space_type");
250      funcSpace = att->as_int(0);      funcSpace = att->as_int(0);
251    
252        if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {
253            centering = NODE_CENTERED;
254        } else {
255            centering = ZONE_CENTERED;
256        }
257    
258        dim = input->get_dim("num_samples");
259        numSamples = dim->size();
260    
261  #ifdef _DEBUG  #ifdef _DEBUG
262      cout << varName << ":\t" << numSamples << " samples,  "      cout << varName << ":\t" << numSamples << " samples,  "
263          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
264  #endif  #endif
265    
266      sampleID.clear();      initialized = true;
     sampleID.insert(sampleID.end(), numSamples, 0);  
     NcVar* var = input->get_var("id");  
     var->get(&sampleID[0], numSamples);  
267    
268      switch (rank) {      // if there are no data samples we're done
269          case 0:      if (numSamples == 0) {
270              readRank0Data(input);          delete input;
271              break;          return true;
         case 1:  
             readRank1Data(input);  
             break;  
         case 2:  
             readRank2Data(input);  
             break;  
         default:  
             cerr << "WARNING: Rank " << rank << " data is not supported!\n";  
             delete input;  
             return false;  
272      }      }
273    
274      delete input;      sampleID.insert(sampleID.end(), numSamples, 0);
275      return true;      NcVar* var = input->get_var("id");
276  }      var->get(&sampleID[0], numSamples);
277    
278  //      size_t dimSize = 1;
279  // Returns one of the mesh names provided by mainMesh that matches the      vector<long> counts;
 // data variable's function space type and reduced/unreduced state.  
 //  
 string DataVar::getMeshName(MeshWithElements* mainMesh) const  
 {  
     string name;  
280    
281      switch (funcSpace) {      if (rank > 0) {
282          case FINLEY_REDUCED_NODES:          dim = input->get_dim("d0");
283          case FINLEY_REDUCED_DEGREES_OF_FREEDOM:          int d = dim->size();
284          case FINLEY_REDUCED_ELEMENTS:          shape.push_back(d);
285          case FINLEY_ELEMENTS:          counts.push_back(d);
286              if (mainMesh->getElements()->reducedCount > 0) {          dimSize *= d;
287                  name = mainMesh->getElements()->reducedMesh->getName();      }
288              } else {      if (rank > 1) {
289                  name = mainMesh->getElements()->fullMesh->getName();          dim = input->get_dim("d1");
290              }          int d = dim->size();
291              break;          shape.push_back(d);
292            counts.push_back(d);
293            dimSize *= d;
294        }
295        if (rank > 2) {
296            cerr << "WARNING: Rank " << rank << " data is not supported!\n";
297            initialized = false;
298        }
299    
300        if (initialized) {
301            size_t dataSize = dimSize*numSamples*ptsPerSample;
302            counts.push_back(ptsPerSample);
303            counts.push_back(numSamples);
304            float* tempData = new float[dataSize];
305            NcVar* 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 = filterSamples(mesh);
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 = new float[numSamples];
     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++) {          float* dest = res;
344          float* c = new float[reorderedNumSamples];          for (int i=0; i<numSamples; i++, src+=stride)
345          reorderedData.push_back(c);              *dest++ = *src;
346          const float* src = rawData[i];      } else {
347          IntVec::const_iterator idIt = requiredIDs.begin();          float* dest = res;
348          for (; idIt != requiredIDs.end(); idIt++) {          for (int i=0; i<numSamples; i++) {
349              size_t srcIdx = sampleID2idx.find(*idIt)->second;              double tmpVal = 0.0;
350              size_t destIdx = id2idxMap.find(*idIt)->second;              for (int j=0; j<ptsPerSample; j++, src+=stride)
351              c[destIdx] = src[srcIdx];                  tmpVal += *src;
352                *dest++ = (float)(tmpVal / ptsPerSample);
353          }          }
354      }      }
355        return res;
356  }  }
357    
358  //  //
359  // For zonal data this method reorders the values according to the indices  // Filters and reorders the raw sample values according to the IDs provided
360  // given in reorderArray. This is used to move ghost zones to the end of  // in 'requiredIDs'. This is used to have data arrays ordered according to
361  // the arrays which conforms to Silo's expected format.  // the underlying 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.  
362  //  //
363  bool DataVar::setMesh(MeshWithElements* mesh)  bool DataVar::filterSamples(FinleyMesh_ptr finleyMesh)
364  {  {
365      if (fullMesh == mesh)      if (numSamples == 0)
366          return true;          return true;
367    
368      const IndexMap* id2idxMap;      IndexMap id2idxMap;
369      const IntVec* reqIDs;      const IntVec* requiredIDs = NULL;
     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;  
370    
371          case FINLEY_REDUCED_ELEMENTS:      NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
372          case FINLEY_ELEMENTS:      if (nodes == NULL)
373              {          return false;
                 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;  
   
         case FINLEY_REDUCED_FACE_ELEMENTS:  
         case FINLEY_FACE_ELEMENTS:  
             {  
                 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;  
   
         case FINLEY_REDUCED_CONTACT_ELEMENTS_1:  
         case FINLEY_REDUCED_CONTACT_ELEMENTS_2:  
         case FINLEY_CONTACT_ELEMENTS_1:  
         case FINLEY_CONTACT_ELEMENTS_2:  
             {  
                 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();  
                 }  
             }  
             break;  
374    
375          case FINLEY_POINTS:      int requiredNumSamples = 0;
             {  
                 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;  
376    
377          default:      if (centering == NODE_CENTERED) {
378              cerr << "Unknown function space type " << funcSpace << "!\n";          id2idxMap = nodes->getIndexMap();
379            requiredIDs = &nodes->getNodeIDs();
380            requiredNumSamples = nodes->getNumNodes();
381        } else {
382            ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
383            if (cells == NULL)
384              return false;              return false;
385    
386            id2idxMap = cells->getIndexMap();
387            requiredIDs = &cells->getIDs();
388            if (cells->getReducedNumElements() > 0) {
389                requiredNumSamples = cells->getReducedNumElements();
390            } else {
391                requiredNumSamples = cells->getNumElements();
392            }
393      }      }
394    
395      if (reorderedNumSamples > numSamples) {      if (requiredNumSamples > numSamples) {
396          cerr << "WARNING: " << varName << " has " << numSamples          cerr << "ERROR: " << varName << " has " << numSamples
397              << " instead of " << reorderedNumSamples << " samples!" << endl;              << " instead of " << requiredNumSamples << " samples!" << endl;
398          return false;          return false;
399      }      }
400    
401      fullMesh = mesh;      meshName = nodes->getName();
402        siloMeshName = nodes->getFullSiloName();
403    
404      reorderSamples(*id2idxMap, *reqIDs);      IndexMap sampleID2idx = buildIndexMap();
405      if (reorderArray)      numSamples = requiredNumSamples;
406          handleGhostZones(*reorderArray);  
407        // now filter the data
408        for (size_t i=0; i < dataArray.size(); i++) {
409            float* c = new float[numSamples];
410            const float* src = dataArray[i];
411            IntVec::const_iterator idIt = requiredIDs->begin();
412            for (; idIt != requiredIDs->end(); idIt++) {
413                size_t srcIdx = sampleID2idx.find(*idIt)->second;
414                size_t destIdx = id2idxMap.find(*idIt)->second;
415                c[destIdx] = src[srcIdx];
416            }
417            delete[] dataArray[i];
418            dataArray[i] = c;
419        }
420      return true;      return true;
421  }  }
422    
# Line 599  bool DataVar::setMesh(MeshWithElements* Line 432  bool DataVar::setMesh(MeshWithElements*
432  //  //
433  string DataVar::getTensorDef() const  string DataVar::getTensorDef() const
434  {  {
435      if (rank < 2)      if (rank < 2 || !initialized)
436          return string();          return string();
437            
438      /// Format string for Silo 2x2 tensor      /// Format string for Silo 2x2 tensor
# Line 642  string DataVar::getTensorDef() const Line 475  string DataVar::getTensorDef() const
475  //  //
476  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)
477  {  {
478  #if HAVE_SILO  #if USE_SILO
479        if (!initialized)
480            return false;
481    
482      if (numSamples == 0)      if (numSamples == 0)
483          return true;          return true;
484    
     // have to set mesh first  
     if (!fullMesh)  
         return false;  
   
485      int ret;      int ret;
486    
487      if (siloPath != "") {      if (siloPath != "") {
# Line 657  bool DataVar::writeToSilo(DBfile* dbfile Line 489  bool DataVar::writeToSilo(DBfile* dbfile
489          if (ret != 0)          if (ret != 0)
490              return false;              return false;
491      }      }
492        
493      char* meshName = (char*)siloMeshName.c_str();      char* siloMesh = const_cast<char*>(siloMeshName.c_str());
494      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
495    
496      if (rank == 0) {      if (rank == 0) {
497          ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName, reorderedData[0],          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
498                  reorderedNumSamples, NULL, 0, DB_FLOAT, dcenter, NULL);                  numSamples, NULL, 0, DB_FLOAT, dcenter, NULL);
499      }      }
500      else if (rank == 1) {      else if (rank == 1) {
501          const string comps[3] = {          const string comps[3] = {
# Line 673  bool DataVar::writeToSilo(DBfile* dbfile Line 505  bool DataVar::writeToSilo(DBfile* dbfile
505              comps[0].c_str(), comps[1].c_str(), comps[2].c_str()              comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
506          };          };
507    
508          ret = DBPutUcdvar(dbfile, varName.c_str(), meshName, shape[0],          ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
509                  (char**)varnames, &reorderedData[0], reorderedNumSamples, NULL,                  (char**)varnames, &dataArray[0], numSamples, NULL,
510                  0, DB_FLOAT, dcenter, NULL);                  0, DB_FLOAT, dcenter, NULL);
511      }      }
512      else {      else {
# Line 689  bool DataVar::writeToSilo(DBfile* dbfile Line 521  bool DataVar::writeToSilo(DBfile* dbfile
521                  for (int j=0; j<shape[0]; j++) {                  for (int j=0; j<shape[0]; j++) {
522                      ostringstream varname;                      ostringstream varname;
523                      varname << tensorDir << "a_" << i << j;                      varname << tensorDir << "a_" << i << j;
524                      ret = DBPutUcdvar1(dbfile, varname.str().c_str(), meshName,                      ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
525                              reorderedData[i*shape[0]+j], reorderedNumSamples,                              dataArray[i*shape[0]+j], numSamples,
526                              NULL, 0, DB_FLOAT, dcenter, optList);                              NULL, 0, DB_FLOAT, dcenter, optList);
527                      if (ret != 0) break;                      if (ret != 0) break;
528                  }                  }
# Line 703  bool DataVar::writeToSilo(DBfile* dbfile Line 535  bool DataVar::writeToSilo(DBfile* dbfile
535      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
536      return (ret == 0);      return (ret == 0);
537    
538  #else // !HAVE_SILO  #else // !USE_SILO
539      return false;      return false;
540  #endif  #endif
541  }  }
542    
543  } // namespace EscriptReader  } // namespace escriptexport
544    
Legend:
Removed from v.2806  
changed lines
  Added in v.2810
  ViewVC Help
Powered by ViewVC 1.1.26