/[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/dataexporter/src/DataVar.cpp revision 2910 by caltinay, Wed Feb 3 03:22:31 2010 UTC trunk/weipa/src/DataVar.cpp revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  #include <escriptexport/DataVar.h>  #include <weipa/DataVar.h>
17  #include <escriptexport/ElementData.h>  #include <weipa/DomainChunk.h>
18  #include <escriptexport/FinleyMesh.h>  #include <weipa/ElementData.h>
19  #include <escriptexport/NodeData.h>  #include <weipa/NodeData.h>
20  #ifndef VISIT_PLUGIN  #ifndef VISIT_PLUGIN
21  #include <escript/Data.h>  #include <escript/Data.h>
22  #endif  #endif
# Line 27  Line 29 
29  #include <silo.h>  #include <silo.h>
30  #endif  #endif
31    
32    #include <numeric> // for accumulate
33    #include <iostream> // for cerr
34    #include <stdio.h>
35    
36  using namespace std;  using namespace std;
37    
38  namespace escriptexport {  namespace weipa {
39            
 enum {  
     NODE_CENTERED = 1,  
     ZONE_CENTERED = 2  
 };  
   
40  //  //
41  // Constructor  // Constructor
42  //  //
43  DataVar::DataVar(const string& name) :  DataVar::DataVar(const string& name) :
44      initialized(false), varName(name),      initialized(false), varName(name),
45      numSamples(0), rank(0), ptsPerSample(0), centering(0)      numSamples(0), rank(0), ptsPerSample(0)
46  {  {
47  }  }
48    
# Line 49  DataVar::DataVar(const string& name) : Line 50  DataVar::DataVar(const string& name) :
50  // Copy constructor  // Copy constructor
51  //  //
52  DataVar::DataVar(const DataVar& d) :  DataVar::DataVar(const DataVar& d) :
53      initialized(d.initialized), finleyMesh(d.finleyMesh),      initialized(d.initialized), domain(d.domain),
54      varName(d.varName), numSamples(d.numSamples),      varName(d.varName), numSamples(d.numSamples),
55      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),      rank(d.rank), ptsPerSample(d.ptsPerSample), funcSpace(d.funcSpace),
56      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)      centering(d.centering), shape(d.shape), sampleID(d.sampleID)
57  {  {
58      if (numSamples > 0) {      if (numSamples > 0) {
59          CoordArray::const_iterator it;          CoordArray::const_iterator it;
# Line 90  void DataVar::cleanup() Line 91  void DataVar::cleanup()
91  //  //
92  //  //
93  //  //
94  bool DataVar::initFromEscript(escript::Data& escriptData, FinleyMesh_ptr mesh)  bool DataVar::initFromEscript(escript::Data& escriptData, const_DomainChunk_ptr dom)
95  {  {
96  #ifndef VISIT_PLUGIN  #ifndef VISIT_PLUGIN
97      cleanup();      cleanup();
98    
99      if (!escriptData.actsExpanded()) {      if (!escriptData.isConstant() && !escriptData.actsExpanded()) {
100          cerr << "WARNING: Only expanded data supported!" << endl;          cerr << "WARNING: Weipa only supports constant & expanded data, "
101                << "not initializing " << varName << endl;
102          return false;          return false;
103      }      }
104    
105      finleyMesh = mesh;      domain = dom;
106      rank = escriptData.getDataPointRank();      rank = escriptData.getDataPointRank();
107      ptsPerSample = escriptData.getNumDataPointsPerSample();      ptsPerSample = escriptData.getNumDataPointsPerSample();
108      shape = escriptData.getDataPointShape();      shape = escriptData.getDataPointShape();
109      funcSpace = escriptData.getFunctionSpace().getTypeCode();      funcSpace = escriptData.getFunctionSpace().getTypeCode();
110      numSamples = escriptData.getNumSamples();      numSamples = escriptData.getNumSamples();
111        centering = domain->getCenteringForFunctionSpace(funcSpace);
     if (funcSpace == FINLEY_REDUCED_NODES || funcSpace == FINLEY_NODES) {  
         centering = NODE_CENTERED;  
     } else {  
         centering = ZONE_CENTERED;  
     }  
112    
113  #ifdef _DEBUG  #ifdef _DEBUG
114      cout << varName << ":\t" << numSamples << " samples,  "      cout << varName << ":\t" << numSamples << " samples,  "
115          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
116  #endif  #endif
117    
118      NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);      NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
119      if (nodes == NULL)      if (nodes == NULL)
120          return false;          return false;
121    
# Line 144  bool DataVar::initFromEscript(escript::D Line 141  bool DataVar::initFromEscript(escript::D
141          initialized = false;          initialized = false;
142      }      }
143    
144        // special case: shape=(1,) or shape=(1,1) -> convert to scalar
145        if (dimSize==1 && rank>0) {
146            rank=0;
147            shape.clear();
148        }
149    
150      if (initialized) {      if (initialized) {
151          size_t dataSize = dimSize * ptsPerSample;          size_t dataSize = dimSize * ptsPerSample;
152          float* tempData = new float[dataSize*numSamples];          float* tempData = new float[dataSize*numSamples];
153          float* destPtr = tempData;          float* destPtr = tempData;
154          for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {          if (escriptData.isConstant()) {
155              const escript::DataAbstract::ValueType::value_type* values =              const escript::DataAbstract::ValueType::value_type* values =
156                  escriptData.getSampleDataRO(sampleNo);                  escriptData.getSampleDataRO(0);
157              copy(values, values+dataSize, destPtr);              for (int pointNo=0; pointNo<numSamples*ptsPerSample; pointNo++) {
158              destPtr += dataSize;                  copy(values, values+dimSize, destPtr);
159                    destPtr += dimSize;
160                }
161            } else {
162                for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
163                    const escript::DataAbstract::ValueType::value_type* values =
164                        escriptData.getSampleDataRO(sampleNo);
165                    copy(values, values+dataSize, destPtr);
166                    destPtr += dataSize;
167                }
168          }          }
169    
170          const float* srcPtr = tempData;          const float* srcPtr = tempData;
# Line 173  bool DataVar::initFromEscript(escript::D Line 185  bool DataVar::initFromEscript(escript::D
185  }  }
186    
187  //  //
188  // Initialise with mesh data  // Initialise with domain data
189  //  //
190  bool DataVar::initFromMesh(FinleyMesh_ptr mesh)  bool DataVar::initFromMeshData(const_DomainChunk_ptr dom, const IntVec& data,
191            int fsCode, Centering c, NodeData_ptr nodes, const IntVec& id)
192  {  {
193      cleanup();      cleanup();
194            
195      finleyMesh = mesh;      domain = dom;
196      rank = 0;      rank = 0;
197      ptsPerSample = 1;      ptsPerSample = 1;
198      NodeData_ptr nodes;      centering = c;
199        sampleID = id;
     if (varName.find("ContactElements_") != varName.npos) {  
         funcSpace = FINLEY_CONTACT_ELEMENTS_1;  
         centering = ZONE_CENTERED;  
         string elementName = varName.substr(0, varName.find('_'));  
         ElementData_ptr elements = mesh->getElementsByName(elementName);  
         nodes = elements->getNodeMesh();  
         sampleID = elements->getIDs();  
     } else if (varName.find("FaceElements_") != varName.npos) {  
         funcSpace = FINLEY_FACE_ELEMENTS;  
         centering = ZONE_CENTERED;  
         string elementName = varName.substr(0, varName.find('_'));  
         ElementData_ptr elements = mesh->getElementsByName(elementName);  
         nodes = elements->getNodeMesh();  
         sampleID = elements->getIDs();  
     } else if (varName.find("Elements_") != varName.npos) {  
         funcSpace = FINLEY_ELEMENTS;  
         centering = ZONE_CENTERED;  
         string elementName = varName.substr(0, varName.find('_'));  
         ElementData_ptr elements = mesh->getElementsByName(elementName);  
         nodes = elements->getNodeMesh();  
         sampleID = elements->getIDs();  
     } else if (varName.find("Nodes_") != varName.npos) {  
         funcSpace = FINLEY_NODES;  
         centering = NODE_CENTERED;  
         nodes = mesh->getNodes();  
         sampleID = nodes->getNodeIDs();  
     } else {  
         cerr << "WARNING: Unrecognized mesh variable '" << varName << "'\n";  
         return false;  
     }  
   
200      meshName = nodes->getName();      meshName = nodes->getName();
201      siloMeshName = nodes->getFullSiloName();      siloMeshName = nodes->getFullSiloName();
   
     const IntVec& data = mesh->getVarDataByName(varName);  
202      numSamples = data.size();      numSamples = data.size();
203    
204      if (numSamples > 0) {      if (numSamples > 0) {
# Line 234  bool DataVar::initFromMesh(FinleyMesh_pt Line 214  bool DataVar::initFromMesh(FinleyMesh_pt
214  }  }
215    
216  //  //
217  // Reads variable data from NetCDF file  // Reads variable data from dump file
218  //  //
219  bool DataVar::initFromNetCDF(const string& filename, FinleyMesh_ptr mesh)  bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
220  {  {
221      cleanup();      cleanup();
222            
# Line 269  bool DataVar::initFromNetCDF(const strin Line 249  bool DataVar::initFromNetCDF(const strin
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) {      centering = dom->getCenteringForFunctionSpace(funcSpace);
         centering = NODE_CENTERED;  
     } else {  
         centering = ZONE_CENTERED;  
     }  
253    
254      dim = input->get_dim("num_samples");      dim = input->get_dim("num_samples");
255      numSamples = dim->size();      numSamples = dim->size();
# Line 283  bool DataVar::initFromNetCDF(const strin Line 259  bool DataVar::initFromNetCDF(const strin
259          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
260  #endif  #endif
261    
262      finleyMesh = mesh;      domain = dom;
263      NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);      NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
264      if (nodes == NULL) {      if (nodes == NULL) {
265          delete input;          delete input;
266          return false;          return false;
# Line 368  float* DataVar::averageData(const float* Line 344  float* DataVar::averageData(const float*
344          for (int i=0; i<numSamples; i++, src+=stride)          for (int i=0; i<numSamples; i++, src+=stride)
345              *dest++ = *src;              *dest++ = *src;
346      } else {      } else {
347          ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);          ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
348          int cellFactor = cells->getElementFactor();          int cellFactor = cells->getElementFactor();
349          res = new float[cellFactor * numSamples];          res = new float[cellFactor * numSamples];
350          float* dest = res;          float* dest = res;
351          QuadMaskInfo qmi = cells->getQuadMask(funcSpace);          QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
352          if (qmi.mask.size() > 0) {          if (!qmi.mask.empty()) {
353              const float* tmpSrc = src;              const float* tmpSrc = src;
354              for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {              for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
355                  for (int l=0; l<cellFactor; l++) {                  for (int l=0; l<cellFactor; l++) {
# Line 417  bool DataVar::reorderSamples() Line 393  bool DataVar::reorderSamples()
393      int cellFactor = 1;      int cellFactor = 1;
394    
395      if (centering == NODE_CENTERED) {      if (centering == NODE_CENTERED) {
396          NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);          NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
397          requiredIDs = &nodes->getNodeIDs();          requiredIDs = &nodes->getNodeIDs();
398          requiredNumSamples = nodes->getNumNodes();          requiredNumSamples = nodes->getNumNodes();
399      } else {      } else {
400          ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);          ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
401          if (cells == NULL)          if (cells == NULL)
402              return false;              return false;
403    
# Line 473  bool DataVar::reorderSamples() Line 449  bool DataVar::reorderSamples()
449  //  //
450  //  //
451  //  //
452    int DataVar::getNumberOfComponents() const
453    {
454        return (rank == 0 ? 1 : accumulate(shape.begin(), shape.end(), 0));
455    }
456    
457    //
458    //
459    //
460    float* DataVar::getDataFlat() const
461    {
462        int totalSize = numSamples * getNumberOfComponents();
463        float* res = new float[totalSize];
464        if (rank == 0) {
465            copy(dataArray[0], dataArray[0]+numSamples, res);
466        } else if (rank == 1) {
467            float *dest = res;
468            for (size_t c=0; c<numSamples; c++) {
469                for (size_t i=0; i<shape[0]; i++) {
470                    *dest++ = dataArray[i][c];
471                }
472            }
473        } else if (rank == 2) {
474            float *dest = res;
475            for (size_t c=0; c<numSamples; c++) {
476                for (int i=0; i<shape[1]; i++) {
477                    for (int j=0; j<shape[0]; j++) {
478                        *dest++ = dataArray[i*shape[0]+j][c];
479                    }
480                }
481            }
482        }
483    
484        return res;
485    }
486    
487    //
488    //
489    //
490  void DataVar::sampleToStream(ostream& os, int index)  void DataVar::sampleToStream(ostream& os, int index)
491  {  {
492        // index is -1 for dummy samples, i.e. if writing the full mesh but
493        // only a reduced number of samples is required
494      if (rank == 0) {      if (rank == 0) {
495          os << dataArray[0][index];          if (index < 0) {
496                os << 0.;
497            } else {
498                os << dataArray[0][index];
499            }
500      } else if (rank == 1) {      } else if (rank == 1) {
501          if (shape[0] < 3)          if (index < 0) {
502                os << 0. << " " << 0.  << " " << 0.;
503            } else if (shape[0] < 3) {
504              os << dataArray[0][index] << " " << dataArray[1][index]              os << dataArray[0][index] << " " << dataArray[1][index]
505                  << " " << 0.;                  << " " << 0.;
506          else          } else {
507              os << dataArray[0][index] << " " << dataArray[1][index]              os << dataArray[0][index] << " " << dataArray[1][index]
508                  << " " << dataArray[2][index];                  << " " << dataArray[2][index];
509            }
510      } else if (rank == 2) {      } else if (rank == 2) {
511          if (shape[1] < 3) {          if (index < 0) {
512                os << 0. << " " << 0. << " " << 0. << " ";
513                os << 0. << " " << 0. << " " << 0. << " ";
514                os << 0. << " " << 0. << " " << 0.;
515            } else if (shape[1] < 3) {
516              os << dataArray[0][index] << " " << dataArray[1][index]              os << dataArray[0][index] << " " << dataArray[1][index]
517                  << " " << 0. << " ";                  << " " << 0. << " ";
518              os << dataArray[2][index] << " " << dataArray[3][index]              os << dataArray[2][index] << " " << dataArray[3][index]
# Line 513  void DataVar::writeToVTK(ostream& os, in Line 540  void DataVar::writeToVTK(ostream& os, in
540    
541      if (isNodeCentered()) {      if (isNodeCentered()) {
542          // data was reordered in reorderSamples() but for VTK we write the          // data was reordered in reorderSamples() but for VTK we write the
543          // original node mesh and thus need the original ordering...          // original node mesh and thus need the original ordering.
544          const IntVec& requiredIDs = finleyMesh->getNodes()->getNodeIDs();          // Note, that this also means we may not have samples for all nodes
545          const IntVec& nodeGNI = finleyMesh->getNodes()->getGlobalNodeIndices();          // in which case we set idx to -1 and write a dummy sample
546          const IntVec& nodeDist = finleyMesh->getNodes()->getNodeDistribution();          const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
547            const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
548            const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
549          int firstId = nodeDist[ownIndex];          int firstId = nodeDist[ownIndex];
550          int lastId = nodeDist[ownIndex+1];          int lastId = nodeDist[ownIndex+1];
551          IndexMap sampleID2idx = buildIndexMap();          IndexMap sampleID2idx = buildIndexMap();
552          for (int i=0; i<nodeGNI.size(); i++) {          for (int i=0; i<nodeGNI.size(); i++) {
553              if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {              if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
554                  int idx = sampleID2idx[requiredIDs[i]];                  IndexMap::const_iterator it = sampleID2idx.find(requiredIDs[i]);
555                    int idx = (it==sampleID2idx.end() ? -1 : (int)it->second);
556                  sampleToStream(os, idx);                  sampleToStream(os, idx);
557              }              }
558          }          }
559      } else {      } else {
560          // cell data: ghost cells have been removed so do not write ghost          // cell data: ghost cells have been removed so do not write ghost
561          // samples (which are the last elements in the arrays)          // samples (which are the last elements in the arrays)
562          int toWrite =          int toWrite = domain->getElementsByName(meshName)->getNumElements();
             finleyMesh->getElementsByName(meshName)->getNumElements();  
563          for (int i=0; i<toWrite; i++) {          for (int i=0; i<toWrite; i++) {
564              sampleToStream(os, i);              sampleToStream(os, i);
565          }          }
# Line 587  string DataVar::getTensorDef() const Line 616  string DataVar::getTensorDef() const
616    
617  //  //
618  // Writes the data to given Silo file under the virtual path provided.  // Writes the data to given Silo file under the virtual path provided.
619  // The corresponding mesh must have been written already and made known  // The corresponding mesh must have been written already.
 // to this variable by a call to setMesh().  
620  //  //
621  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath)  bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
622                              const string& units)
623  {  {
624  #if USE_SILO  #if USE_SILO
625      if (!initialized)      if (!initialized)
# Line 609  bool DataVar::writeToSilo(DBfile* dbfile Line 638  bool DataVar::writeToSilo(DBfile* dbfile
638    
639      char* siloMesh = const_cast<char*>(siloMeshName.c_str());      char* siloMesh = const_cast<char*>(siloMeshName.c_str());
640      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);      int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
641        DBoptlist* optList = DBMakeOptlist(2);
642        if (units.length()>0) {
643            DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
644        }
645    
646      if (rank == 0) {      if (rank == 0) {
647          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
648                  numSamples, NULL, 0, DB_FLOAT, dcenter, NULL);                  numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
649      }      }
650      else if (rank == 1) {      else if (rank == 1) {
651          const string comps[3] = {          const string comps[3] = {
# Line 624  bool DataVar::writeToSilo(DBfile* dbfile Line 657  bool DataVar::writeToSilo(DBfile* dbfile
657    
658          ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],          ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
659                  (char**)varnames, &dataArray[0], numSamples, NULL,                  (char**)varnames, &dataArray[0], numSamples, NULL,
660                  0, DB_FLOAT, dcenter, NULL);                  0, DB_FLOAT, dcenter, optList);
661      }      }
662      else {      else {
663          string tensorDir = varName+string("_comps/");          string tensorDir = varName+string("_comps/");
664          ret = DBMkdir(dbfile, tensorDir.c_str());          ret = DBMkdir(dbfile, tensorDir.c_str());
665          if (ret == 0) {          if (ret == 0) {
666              int one = 1;              int one = 1;
             DBoptlist* optList = DBMakeOptlist(1);  
667              DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);              DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
668    
669              for (int i=0; i<shape[1]; i++) {              for (int i=0; i<shape[1]; i++) {
# Line 645  bool DataVar::writeToSilo(DBfile* dbfile Line 677  bool DataVar::writeToSilo(DBfile* dbfile
677                  }                  }
678                  if (ret != 0) break;                  if (ret != 0) break;
679              }              }
             DBFreeOptlist(optList);  
680          } // ret==0          } // ret==0
681      } // rank      } // rank
682    
683        DBFreeOptlist(optList);
684      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
685      return (ret == 0);      return (ret == 0);
686    
# Line 657  bool DataVar::writeToSilo(DBfile* dbfile Line 689  bool DataVar::writeToSilo(DBfile* dbfile
689  #endif  #endif
690  }  }
691    
692  } // namespace escriptexport  } // namespace weipa
693    

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

  ViewVC Help
Powered by ViewVC 1.1.26