/[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

revision 2810 by caltinay, Mon Dec 7 04:13:49 2009 UTC revision 2910 by caltinay, Wed Feb 3 03:22:31 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2009 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 15  Line 15 
15  #include <escriptexport/ElementData.h>  #include <escriptexport/ElementData.h>
16  #include <escriptexport/FinleyMesh.h>  #include <escriptexport/FinleyMesh.h>
17  #include <escriptexport/NodeData.h>  #include <escriptexport/NodeData.h>
18    #ifndef VISIT_PLUGIN
19  #include <escript/Data.h>  #include <escript/Data.h>
20    #endif
21    
22  #if USE_NETCDF  #if USE_NETCDF
23  #include <netcdf.hh>  #include <netcdfcpp.h>
24  #endif  #endif
25    
26  #if USE_SILO  #if USE_SILO
# Line 47  DataVar::DataVar(const string& name) : Line 49  DataVar::DataVar(const string& name) :
49  // Copy constructor  // Copy constructor
50  //  //
51  DataVar::DataVar(const DataVar& d) :  DataVar::DataVar(const DataVar& d) :
52        initialized(d.initialized), finleyMesh(d.finleyMesh),
53      varName(d.varName), numSamples(d.numSamples),      varName(d.varName), numSamples(d.numSamples),
54      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),
55      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)
# Line 59  DataVar::DataVar(const DataVar& d) : Line 62  DataVar::DataVar(const DataVar& d) :
62              dataArray.push_back(c);              dataArray.push_back(c);
63          }          }
64      }      }
     initialized = d.initialized;  
65  }  }
66    
67  //  //
# Line 90  void DataVar::cleanup() Line 92  void DataVar::cleanup()
92  //  //
93  bool DataVar::initFromEscript(escript::Data& escriptData, FinleyMesh_ptr mesh)  bool DataVar::initFromEscript(escript::Data& escriptData, FinleyMesh_ptr mesh)
94  {  {
95    #ifndef VISIT_PLUGIN
96      cleanup();      cleanup();
97    
98      if (!escriptData.actsExpanded()) {      if (!escriptData.actsExpanded()) {
# Line 97  bool DataVar::initFromEscript(escript::D Line 100  bool DataVar::initFromEscript(escript::D
100          return false;          return false;
101      }      }
102    
103        finleyMesh = mesh;
104      rank = escriptData.getDataPointRank();      rank = escriptData.getDataPointRank();
105      ptsPerSample = escriptData.getNumDataPointsPerSample();      ptsPerSample = escriptData.getNumDataPointsPerSample();
106      shape = escriptData.getDataPointShape();      shape = escriptData.getDataPointShape();
# Line 114  bool DataVar::initFromEscript(escript::D Line 118  bool DataVar::initFromEscript(escript::D
118          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
119  #endif  #endif
120    
121        NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
122        if (nodes == NULL)
123            return false;
124    
125        meshName = nodes->getName();
126        siloMeshName = nodes->getFullSiloName();
127      initialized = true;      initialized = true;
128    
129        // no samples? Nothing more to do.
130      if (numSamples == 0)      if (numSamples == 0)
131          return true;          return true;
132    
# Line 151  bool DataVar::initFromEscript(escript::D Line 162  bool DataVar::initFromEscript(escript::D
162          }          }
163          delete[] tempData;          delete[] tempData;
164    
165          initialized = filterSamples(mesh);          initialized = reorderSamples();
166      }      }
167    
168      return initialized;      return initialized;
169    
170    #else // VISIT_PLUGIN
171        return false;
172    #endif
173  }  }
174    
175  //  //
# Line 164  bool DataVar::initFromMesh(FinleyMesh_pt Line 179  bool DataVar::initFromMesh(FinleyMesh_pt
179  {  {
180      cleanup();      cleanup();
181            
182      const IntVec& data = mesh->getVarDataByName(varName);      finleyMesh = mesh;
183      rank = 0;      rank = 0;
184      ptsPerSample = 1;      ptsPerSample = 1;
185        NodeData_ptr nodes;
186    
187        if (varName.find("ContactElements_") != varName.npos) {
188            funcSpace = FINLEY_CONTACT_ELEMENTS_1;
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("FaceElements_") != varName.npos) {
195            funcSpace = FINLEY_FACE_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("Elements_") != varName.npos) {
202            funcSpace = FINLEY_ELEMENTS;
203            centering = ZONE_CENTERED;
204            string elementName = varName.substr(0, varName.find('_'));
205            ElementData_ptr elements = mesh->getElementsByName(elementName);
206            nodes = elements->getNodeMesh();
207            sampleID = elements->getIDs();
208        } else if (varName.find("Nodes_") != varName.npos) {
209            funcSpace = FINLEY_NODES;
210            centering = NODE_CENTERED;
211            nodes = mesh->getNodes();
212            sampleID = nodes->getNodeIDs();
213        } else {
214            cerr << "WARNING: Unrecognized mesh variable '" << varName << "'\n";
215            return false;
216        }
217    
218        meshName = nodes->getName();
219        siloMeshName = nodes->getFullSiloName();
220    
221        const IntVec& data = mesh->getVarDataByName(varName);
222      numSamples = data.size();      numSamples = data.size();
223    
224      if (numSamples > 0) {      if (numSamples > 0) {
# Line 175  bool DataVar::initFromMesh(FinleyMesh_pt Line 227  bool DataVar::initFromMesh(FinleyMesh_pt
227          IntVec::const_iterator it;          IntVec::const_iterator it;
228          for (it=data.begin(); it != data.end(); it++)          for (it=data.begin(); it != data.end(); it++)
229              *c++ = static_cast<float>(*it);              *c++ = static_cast<float>(*it);
   
         if (varName.compare(0, 6, "Nodes_") == 0) {  
             funcSpace = FINLEY_NODES;  
             centering = NODE_CENTERED;  
             sampleID.insert(sampleID.end(),  
                     mesh->getNodes()->getNodeIDs().begin(),  
                     mesh->getNodes()->getNodeIDs().end());  
         } else if (varName.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 (varName.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 (varName.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 {  
             return false;  
         }  
   
         NodeData_ptr nodes = mesh->getMeshForFinleyFS(funcSpace);  
         meshName = nodes->getName();  
         siloMeshName = nodes->getFullSiloName();  
230      }      }
231      initialized = true;      initialized = true;
232    
# Line 263  bool DataVar::initFromNetCDF(const strin Line 283  bool DataVar::initFromNetCDF(const strin
283          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
284  #endif  #endif
285    
286      initialized = true;      finleyMesh = mesh;
287        NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
288      // if there are no data samples we're done      if (nodes == NULL) {
     if (numSamples == 0) {  
289          delete input;          delete input;
290          return true;          return false;
291      }      }
292    
293      sampleID.insert(sampleID.end(), numSamples, 0);      meshName = nodes->getName();
294      NcVar* var = input->get_var("id");      siloMeshName = nodes->getFullSiloName();
295      var->get(&sampleID[0], numSamples);      initialized = true;
296    
297      size_t dimSize = 1;      size_t dimSize = 1;
298      vector<long> counts;      vector<long> counts;
# Line 297  bool DataVar::initFromNetCDF(const strin Line 316  bool DataVar::initFromNetCDF(const strin
316          initialized = false;          initialized = false;
317      }      }
318    
319      if (initialized) {      if (initialized && numSamples > 0) {
320            sampleID.insert(sampleID.end(), numSamples, 0);
321            NcVar* var = input->get_var("id");
322            var->get(&sampleID[0], numSamples);
323    
324          size_t dataSize = dimSize*numSamples*ptsPerSample;          size_t dataSize = dimSize*numSamples*ptsPerSample;
325          counts.push_back(ptsPerSample);          counts.push_back(ptsPerSample);
326          counts.push_back(numSamples);          counts.push_back(numSamples);
327          float* tempData = new float[dataSize];          float* tempData = new float[dataSize];
328          NcVar* var = input->get_var("data");          var = input->get_var("data");
329          var->get(tempData, &counts[0]);          var->get(tempData, &counts[0]);
330    
331          const float* srcPtr = tempData;          const float* srcPtr = tempData;
# Line 312  bool DataVar::initFromNetCDF(const strin Line 335  bool DataVar::initFromNetCDF(const strin
335          }          }
336          delete[] tempData;          delete[] tempData;
337    
338          initialized = filterSamples(mesh);          initialized = reorderSamples();
339      }      }
340    
341      delete input;      delete input;
# Line 337  bool DataVar::isNodeCentered() const Line 360  bool DataVar::isNodeCentered() const
360  //  //
361  float* DataVar::averageData(const float* src, size_t stride)  float* DataVar::averageData(const float* src, size_t stride)
362  {  {
363      float* res = new float[numSamples];      float* res;
364    
365      if (ptsPerSample == 1) {      if (ptsPerSample == 1) {
366            res = new float[numSamples];
367          float* dest = res;          float* dest = res;
368          for (int i=0; i<numSamples; i++, src+=stride)          for (int i=0; i<numSamples; i++, src+=stride)
369              *dest++ = *src;              *dest++ = *src;
370      } else {      } else {
371            ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
372            int cellFactor = cells->getElementFactor();
373            res = new float[cellFactor * numSamples];
374          float* dest = res;          float* dest = res;
375          for (int i=0; i<numSamples; i++) {          QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
376              double tmpVal = 0.0;          if (qmi.mask.size() > 0) {
377              for (int j=0; j<ptsPerSample; j++, src+=stride)              const float* tmpSrc = src;
378                  tmpVal += *src;              for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
379              *dest++ = (float)(tmpVal / ptsPerSample);                  for (int l=0; l<cellFactor; l++) {
380                        double tmpVal = 0.0;
381                        for (int j=0; j<ptsPerSample; j++) {
382                            if (qmi.mask[l][j] != 0) {
383                                tmpVal += *(tmpSrc+stride*j);
384                            }
385                        }
386                        *dest++ = (float)(tmpVal / qmi.factor[l]);
387                    }
388                }
389            } else {
390                for (int i=0; i<numSamples; i++) {
391                    double tmpVal = 0.0;
392                    for (int j=0; j<ptsPerSample; j++, src+=stride) {
393                        tmpVal += *src;
394                    }
395                    tmpVal /= ptsPerSample;
396                    for (int l=0; l<cellFactor; l++) {
397                        *dest++ = static_cast<float>(tmpVal);
398                    }
399                }
400          }          }
401      }      }
402      return res;      return res;
403  }  }
404    
405  //  //
406  // Filters and reorders the raw sample values according to the IDs provided  // Filters and reorders the raw sample values according to the node/element
407  // in 'requiredIDs'. This is used to have data arrays ordered according to  // IDs. This is used to have data arrays ordered according to the underlying
408  // the underlying mesh (i.e. DataID[i]==MeshNodeID[i])  // mesh (i.e. DataID[i]==MeshNodeID[i])
409  //  //
410  bool DataVar::filterSamples(FinleyMesh_ptr finleyMesh)  bool DataVar::reorderSamples()
411  {  {
412      if (numSamples == 0)      if (numSamples == 0)
413          return true;          return true;
414    
     IndexMap id2idxMap;  
415      const IntVec* requiredIDs = NULL;      const IntVec* requiredIDs = NULL;
   
     NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);  
     if (nodes == NULL)  
         return false;  
   
416      int requiredNumSamples = 0;      int requiredNumSamples = 0;
417        int cellFactor = 1;
418    
419      if (centering == NODE_CENTERED) {      if (centering == NODE_CENTERED) {
420          id2idxMap = nodes->getIndexMap();          NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
421          requiredIDs = &nodes->getNodeIDs();          requiredIDs = &nodes->getNodeIDs();
422          requiredNumSamples = nodes->getNumNodes();          requiredNumSamples = nodes->getNumNodes();
423      } else {      } else {
# Line 383  bool DataVar::filterSamples(FinleyMesh_p Line 425  bool DataVar::filterSamples(FinleyMesh_p
425          if (cells == NULL)          if (cells == NULL)
426              return false;              return false;
427    
         id2idxMap = cells->getIndexMap();  
428          requiredIDs = &cells->getIDs();          requiredIDs = &cells->getIDs();
429          if (cells->getReducedNumElements() > 0) {          requiredNumSamples = cells->getNumElements();
430              requiredNumSamples = cells->getReducedNumElements();          cellFactor = cells->getElementFactor();
431          } else {          if (cellFactor > 1) {
432              requiredNumSamples = cells->getNumElements();              numSamples *= cellFactor;
433                // update sample IDs
434                IntVec newSampleID(numSamples);
435                IntVec::const_iterator idIt = sampleID.begin();
436                IntVec::iterator newIDit = newSampleID.begin();
437                for (; idIt != sampleID.end(); idIt++, newIDit+=cellFactor) {
438                    fill(newIDit, newIDit+cellFactor, *idIt);
439                }
440                sampleID.swap(newSampleID);
441          }          }
442      }      }
443    
# Line 398  bool DataVar::filterSamples(FinleyMesh_p Line 447  bool DataVar::filterSamples(FinleyMesh_p
447          return false;          return false;
448      }      }
449    
     meshName = nodes->getName();  
     siloMeshName = nodes->getFullSiloName();  
   
450      IndexMap sampleID2idx = buildIndexMap();      IndexMap sampleID2idx = buildIndexMap();
451      numSamples = requiredNumSamples;      numSamples = requiredNumSamples;
452    
# Line 409  bool DataVar::filterSamples(FinleyMesh_p Line 455  bool DataVar::filterSamples(FinleyMesh_p
455          float* c = new float[numSamples];          float* c = new float[numSamples];
456          const float* src = dataArray[i];          const float* src = dataArray[i];
457          IntVec::const_iterator idIt = requiredIDs->begin();          IntVec::const_iterator idIt = requiredIDs->begin();
458          for (; idIt != requiredIDs->end(); idIt++) {          size_t destIdx = 0;
459            for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
460              size_t srcIdx = sampleID2idx.find(*idIt)->second;              size_t srcIdx = sampleID2idx.find(*idIt)->second;
461              size_t destIdx = id2idxMap.find(*idIt)->second;              copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
             c[destIdx] = src[srcIdx];  
462          }          }
463          delete[] dataArray[i];          delete[] dataArray[i];
464          dataArray[i] = c;          dataArray[i] = c;
465      }      }
466    
467        // sample IDs now = mesh node/element IDs
468        sampleID = *requiredIDs;
469    
470      return true;      return true;
471  }  }
472    
473    //
474    //
475    //
476    void DataVar::sampleToStream(ostream& os, int index)
477    {
478        if (rank == 0) {
479            os << dataArray[0][index];
480        } else if (rank == 1) {
481            if (shape[0] < 3)
482                os << dataArray[0][index] << " " << dataArray[1][index]
483                    << " " << 0.;
484            else
485                os << dataArray[0][index] << " " << dataArray[1][index]
486                    << " " << dataArray[2][index];
487        } else if (rank == 2) {
488            if (shape[1] < 3) {
489                os << dataArray[0][index] << " " << dataArray[1][index]
490                    << " " << 0. << " ";
491                os << dataArray[2][index] << " " << dataArray[3][index]
492                    << " " << 0. << " ";
493                os << 0. << " " << 0. << " " << 0.;
494            } else {
495                os << dataArray[0][index] << " " << dataArray[1][index]
496                    << " " << dataArray[2][index] << " ";
497                os << dataArray[3][index] << " " << dataArray[4][index]
498                    << " " << dataArray[5][index] << " ";
499                os << dataArray[6][index] << " " << dataArray[7][index]
500                    << " " << dataArray[8][index];
501            }
502        }
503        os << endl;
504    }
505    
506    //
507    //
508    //
509    void DataVar::writeToVTK(ostream& os, int ownIndex)
510    {
511        if (numSamples == 0)
512            return;
513    
514        if (isNodeCentered()) {
515            // data was reordered in reorderSamples() but for VTK we write the
516            // original node mesh and thus need the original ordering...
517            const IntVec& requiredIDs = finleyMesh->getNodes()->getNodeIDs();
518            const IntVec& nodeGNI = finleyMesh->getNodes()->getGlobalNodeIndices();
519            const IntVec& nodeDist = finleyMesh->getNodes()->getNodeDistribution();
520            int firstId = nodeDist[ownIndex];
521            int lastId = nodeDist[ownIndex+1];
522            IndexMap sampleID2idx = buildIndexMap();
523            for (int i=0; i<nodeGNI.size(); i++) {
524                if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
525                    int idx = sampleID2idx[requiredIDs[i]];
526                    sampleToStream(os, idx);
527                }
528            }
529        } else {
530            // cell data: ghost cells have been removed so do not write ghost
531            // samples (which are the last elements in the arrays)
532            int toWrite =
533                finleyMesh->getElementsByName(meshName)->getNumElements();
534            for (int i=0; i<toWrite; i++) {
535                sampleToStream(os, i);
536            }
537        }
538    }
539    
540  ///////////////////////////////  ///////////////////////////////
541  // SILO related methods follow  // SILO related methods follow
542  ///////////////////////////////  ///////////////////////////////

Legend:
Removed from v.2810  
changed lines
  Added in v.2910

  ViewVC Help
Powered by ViewVC 1.1.26