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

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

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

trunk/tools/libescriptreader/src/escriptreader/ElementData.cpp revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC trunk/dataexporter/src/ElementData.cpp revision 2810 by caltinay, Mon Dec 7 04:13:49 2009 UTC
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
14  //  #include <escriptexport/ElementData.h>
15  // ElementData.cpp  #include <escriptexport/NodeData.h>
16  //  
17  #include <escriptreader/ElementData.h>  extern "C" {
18    #include <finley/ElementFile.h>
19    }
20    
21    #if USE_NETCDF
22  #include <netcdf.hh>  #include <netcdf.hh>
23  #if HAVE_SILO  #endif
24    
25    #if USE_SILO
26  #include <silo.h>  #include <silo.h>
27  #endif  #endif
28    
# Line 40  static const size_t rec8indices[6*3] = { Line 46  static const size_t rec8indices[6*3] = {
46      4, 1, 5,      4, 1, 5,
47      5, 2, 6,      5, 2, 6,
48      6, 3, 7,      6, 3, 7,
49      4, 6, 7,      7, 5, 6,
50      4, 5, 6      7, 4, 5
51  };  };
52  static const size_t rec9indices[4*4] = {  static const size_t rec9indices[4*4] = {
53      0, 4, 8, 7,      0, 4, 8, 7,
# Line 84  static const size_t hex27indices[8*8] = Line 90  static const size_t hex27indices[8*8] =
90      26, 22, 14, 23, 25, 17,  6, 18      26, 22, 14, 23, 25, 17,  6, 18
91  };  };
92    
93  namespace EscriptReader {  namespace escriptexport {
94            
95  //  //
96  // Constructor  // Constructor
97  //  //
98  ElementData::ElementData(const string& elementName, const Mesh* mainMesh)  ElementData::ElementData(const string& elementName, NodeData_ptr nodeData)
99      : name(elementName), count(0), reducedCount(0), numGhostElements(0),      : name(elementName), numElements(0), reducedNumElements(0),
100      numReducedGhostElements(0), fullMesh(NULL), reducedMesh(NULL),      numGhostElements(0), numReducedGhostElements(0), originalNodes(nodeData),
101      originalMesh(mainMesh), fullMeshIsOriginalMesh(false)      fullMeshIsOriginalMesh(false)
102  {  {
103  }  }
104    
# Line 102  ElementData::ElementData(const string& e Line 108  ElementData::ElementData(const string& e
108  ElementData::ElementData(const ElementData& e)  ElementData::ElementData(const ElementData& e)
109  {  {
110      name = e.name;      name = e.name;
111      count = e.count;      numElements = e.numElements;
112      reducedCount = e.reducedCount;      reducedNumElements = e.reducedNumElements;
113      numGhostElements = e.numGhostElements;      numGhostElements = e.numGhostElements;
114      numReducedGhostElements = e.numReducedGhostElements;      numReducedGhostElements = e.numReducedGhostElements;
115      numDims = e.numDims;      numDims = e.numDims;
# Line 112  ElementData::ElementData(const ElementDa Line 118  ElementData::ElementData(const ElementDa
118      nodesPerElement = e.nodesPerElement;      nodesPerElement = e.nodesPerElement;
119      reducedNodesPerElement = e.reducedNodesPerElement;      reducedNodesPerElement = e.reducedNodesPerElement;
120    
121      originalMesh = e.originalMesh;      originalNodes = e.originalNodes;
122      fullMeshIsOriginalMesh = e.fullMeshIsOriginalMesh;      fullMeshIsOriginalMesh = e.fullMeshIsOriginalMesh;
123      if (fullMeshIsOriginalMesh)      if (fullMeshIsOriginalMesh)
124          fullMesh = const_cast<Mesh*>(originalMesh);          fullMesh = originalNodes;
125      else if (e.fullMesh)      else if (e.fullMesh)
126          fullMesh = new Mesh(*e.fullMesh);          fullMesh = NodeData_ptr(new NodeData(*e.fullMesh));
     else  
         fullMesh = NULL;  
127    
128      if (e.reducedMesh)      if (e.reducedMesh)
129          reducedMesh = new Mesh(*e.reducedMesh);          reducedMesh = NodeData_ptr(new NodeData(*e.reducedMesh));
     else  
         reducedMesh = NULL;  
130    
131      nodes = e.nodes;      nodes = e.nodes;
132      reducedNodes = e.reducedNodes;      reducedNodes = e.reducedNodes;
# Line 133  ElementData::ElementData(const ElementDa Line 135  ElementData::ElementData(const ElementDa
135      tag = e.tag;      tag = e.tag;
136      owner = e.owner;      owner = e.owner;
137      reducedOwner = e.reducedOwner;      reducedOwner = e.reducedOwner;
     indexArray = e.indexArray;  
     reducedIndexArray = e.reducedIndexArray;  
     ID2idx = e.ID2idx;  
138  }  }
139    
140  //  //
# Line 143  ElementData::ElementData(const ElementDa Line 142  ElementData::ElementData(const ElementDa
142  //  //
143  ElementData::~ElementData()  ElementData::~ElementData()
144  {  {
145      if (reducedMesh)  }
146          delete reducedMesh;  
147      if (fullMesh && !fullMeshIsOriginalMesh)  bool ElementData::initFromFinley(const Finley_ElementFile* finleyFile)
148          delete fullMesh;  {
149        numElements = finleyFile->numElements;
150    
151        if (numElements > 0) {
152            nodesPerElement = finleyFile->numNodes;
153    
154            int* iPtr;
155      
156            iPtr = finleyFile->Nodes;
157            nodes.clear();
158            nodes.insert(nodes.end(), numElements*nodesPerElement, 0);
159            copy(iPtr, iPtr+numElements*nodesPerElement, nodes.begin());
160    
161            iPtr = finleyFile->Color;
162            color.clear();
163            color.insert(color.end(), numElements, 0);
164            copy(iPtr, iPtr+numElements, color.begin());
165    
166            iPtr = finleyFile->Id;
167            ID.clear();
168            ID.insert(ID.end(), numElements, 0);
169            copy(iPtr, iPtr+numElements, ID.begin());
170    
171            iPtr = finleyFile->Owner;
172            owner.clear();
173            owner.insert(owner.end(), numElements, 0);
174            copy(iPtr, iPtr+numElements, owner.begin());
175    
176            iPtr = finleyFile->Tag;
177            tag.clear();
178            tag.insert(tag.end(), numElements, 0);
179            copy(iPtr, iPtr+numElements, tag.begin());
180    
181            ElementTypeId tid = finleyFile->referenceElementSet->
182                referenceElement->Type->TypeId;
183            FinleyElementInfo f = getFinleyTypeInfo(tid);
184            type = f.elementType;
185            reducedType = f.reducedElementType;
186            if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement)
187                buildReducedElements(f);
188    
189            buildMeshes();
190        }
191        return true;
192  }  }
193    
194  StringVec ElementData::getMeshNames() const  StringVec ElementData::getMeshNames() const
# Line 162  StringVec ElementData::getMeshNames() co Line 204  StringVec ElementData::getMeshNames() co
204  StringVec ElementData::getVarNames() const  StringVec ElementData::getVarNames() const
205  {  {
206      StringVec res;      StringVec res;
207      if (count > 0) {      if (numElements > 0) {
208          res.push_back(name + string("_Color"));          res.push_back(name + string("_Color"));
209          res.push_back(name + string("_Id"));          res.push_back(name + string("_Id"));
210          res.push_back(name + string("_Owner"));          res.push_back(name + string("_Owner"));
# Line 178  const IntVec& ElementData::getVarDataByN Line 220  const IntVec& ElementData::getVarDataByN
220      else if (varName == name+string("_Id"))      else if (varName == name+string("_Id"))
221          return ID;          return ID;
222      else if (varName == name+string("_Owner")) {      else if (varName == name+string("_Owner")) {
223          if (reducedCount > 0)          if (reducedNumElements > 0)
224              return reducedOwner;              return reducedOwner;
225          else          else
226              return owner;              return owner;
# Line 193  const IntVec& ElementData::getVarDataByN Line 235  const IntVec& ElementData::getVarDataByN
235  //  //
236  bool ElementData::readFromNc(NcFile* ncfile)  bool ElementData::readFromNc(NcFile* ncfile)
237  {  {
238      NcAtt* att;  #if USE_NETCDF
     NcVar* var;  
   
239      string num_str("num_");      string num_str("num_");
240      num_str += name;      num_str += name;
241    
242      att = ncfile->get_att(num_str.c_str());      NcAtt* att = ncfile->get_att(num_str.c_str());
243      count = att->as_int(0);      numElements = att->as_int(0);
244    
245      // Only attempt to read data if there are any elements.      // Only attempt to read data if there are any elements.
246      // Having no elements is not an error.      // Having no elements is not an error.
247      if (count == 0)      if (numElements > 0) {
248          return true;          att = ncfile->get_att((num_str + string("_numNodes")).c_str());
249            nodesPerElement = att->as_int(0);
250    
251            nodes.insert(nodes.end(), numElements*nodesPerElement, 0);
252            NcVar* var = ncfile->get_var((name + string("_Nodes")).c_str());
253            var->get(&nodes[0], numElements, nodesPerElement);
254    
255            color.insert(color.end(), numElements, 0);
256            var = ncfile->get_var((name + string("_Color")).c_str());
257            var->get(&color[0], numElements);
258    
259            ID.insert(ID.end(), numElements, 0);
260            var = ncfile->get_var((name + string("_Id")).c_str());
261            var->get(&ID[0], numElements);
262    
263            owner.insert(owner.end(), numElements, 0);
264            var = ncfile->get_var((name + string("_Owner")).c_str());
265            var->get(&owner[0], numElements);
266    
267            tag.insert(tag.end(), numElements, 0);
268            var = ncfile->get_var((name + string("_Tag")).c_str());
269            var->get(&tag[0], numElements);
270    
271            att = ncfile->get_att((name + string("_TypeId")).c_str());
272            FinleyElementInfo f = getFinleyTypeInfo((ElementTypeId)att->as_int(0));
273            type = f.elementType;
274            reducedType = f.reducedElementType;
275    
276      att = ncfile->get_att((num_str + string("_numNodes")).c_str());          if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement)
277      nodesPerElement = att->as_int(0);              buildReducedElements(f);
278    
279      nodes.insert(nodes.end(), count*nodesPerElement, 0);          buildMeshes();
280      var = ncfile->get_var((name + string("_Nodes")).c_str());      }
     var->get(&nodes[0], count, nodesPerElement);  
   
     color.insert(color.end(), count, 0);  
     var = ncfile->get_var((name + string("_Color")).c_str());  
     var->get(&color[0], count);  
   
     ID.insert(ID.end(), count, 0);  
     var = ncfile->get_var((name + string("_Id")).c_str());  
     var->get(&ID[0], count);  
   
     owner.insert(owner.end(), count, 0);  
     var = ncfile->get_var((name + string("_Owner")).c_str());  
     var->get(&owner[0], count);  
   
     tag.insert(tag.end(), count, 0);  
     var = ncfile->get_var((name + string("_Tag")).c_str());  
     var->get(&tag[0], count);  
   
     att = ncfile->get_att((name + string("_TypeId")).c_str());  
     FinleyElementInfo f = getFinleyTypeInfo((ElementTypeId)att->as_int(0));  
     type = f.elementType;  
     reducedType = f.reducedElementType;  
   
     // build elementID->index map  
     buildIndexMap();  
       
     if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement)  
         buildReducedElements(f);  
   
     buildMeshes();  
281    
282      return true;      return true;
283    #else // !USE_NETCDF
284        return false;
285    #endif
286  }  }
287    
288  //  //
# Line 276  void ElementData::buildReducedElements(c Line 315  void ElementData::buildReducedElements(c
315  {  {
316      reducedNodes.clear();      reducedNodes.clear();
317      reducedNodes.insert(reducedNodes.end(),      reducedNodes.insert(reducedNodes.end(),
318              f.reducedElementSize*count, 0);              f.reducedElementSize*numElements, 0);
319      IntVec::iterator reducedIt = reducedNodes.begin();      IntVec::iterator reducedIt = reducedNodes.begin();
320      IntVec::const_iterator origIt;      IntVec::const_iterator origIt;
321      for (origIt=nodes.begin(); origIt!=nodes.end();      for (origIt=nodes.begin(); origIt!=nodes.end();
# Line 290  void ElementData::buildReducedElements(c Line 329  void ElementData::buildReducedElements(c
329      // elements except for the main mesh      // elements except for the main mesh
330      if (f.elementFactor > 1 /*&& name == "Elements"*/) {      if (f.elementFactor > 1 /*&& name == "Elements"*/) {
331          // replace each element by multiple smaller ones          // replace each element by multiple smaller ones
332          IntVec fullNodes(f.elementSize*f.elementFactor*count);          IntVec fullNodes(f.elementSize*f.elementFactor*numElements);
333          IntVec::iterator cellIt = fullNodes.begin();          IntVec::iterator cellIt = fullNodes.begin();
334    
335          // copy owner data          // copy owner data
336          owner.swap(reducedOwner);          owner.swap(reducedOwner);
337          owner.clear();          owner.clear();
338          for (int i=0; i < count; i++) {          for (int i=0; i < numElements; i++) {
339              owner.insert(owner.end(), f.elementFactor, reducedOwner[i]);              owner.insert(owner.end(), f.elementFactor, reducedOwner[i]);
340              for (int j=0; j < f.elementFactor*f.elementSize; j++)              for (int j=0; j < f.elementFactor*f.elementSize; j++)
341                  *cellIt++ = nodes[                  *cellIt++ = nodes[
342                      nodesPerElement*i+f.multiCellIndices[j]];                      nodesPerElement*i+f.multiCellIndices[j]];
343          }          }
344          nodes.swap(fullNodes);          nodes.swap(fullNodes);
345          reducedCount = count;          reducedNumElements = numElements;
346          reducedNodesPerElement = f.reducedElementSize;          reducedNodesPerElement = f.reducedElementSize;
347          nodesPerElement = f.elementSize;          nodesPerElement = f.elementSize;
348          count *= f.elementFactor;          numElements *= f.elementFactor;
349      } else {      } else {
350          // we only keep the reduced elements but treat them as regular          // we only keep the reduced elements but treat them as regular
351          // ones, so replace the data accordingly          // ones, so replace the data accordingly
# Line 320  void ElementData::buildReducedElements(c Line 359  void ElementData::buildReducedElements(c
359  //  //
360  //  //
361  //  //
362  void ElementData::prepareGhostIndices(int ownIndex)  void ElementData::prepareGhostIndices(int ownIndex, IntVec& indexArray,
363                                          IntVec& reducedIndexArray)
364  {  {
365      indexArray.clear();      indexArray.clear();
366      reducedIndexArray.clear();      reducedIndexArray.clear();
# Line 329  void ElementData::prepareGhostIndices(in Line 369  void ElementData::prepareGhostIndices(in
369            
370      // move indices of "ghost zones" to the end to be able to reorder      // move indices of "ghost zones" to the end to be able to reorder
371      // data accordingly      // data accordingly
372      for (int i=0; i<count; i++)      for (int i=0; i<numElements; i++)
373          if (owner[i] == ownIndex)          if (owner[i] == ownIndex)
374              indexArray.push_back(i);              indexArray.push_back(i);
375    
376      for (int i=0; i<count; i++)      for (int i=0; i<numElements; i++)
377          if (owner[i] != ownIndex) {          if (owner[i] != ownIndex) {
378              numGhostElements++;              numGhostElements++;
379              indexArray.push_back(i);              indexArray.push_back(i);
380          }          }
381    
382      for (int i=0; i<reducedCount; i++)      for (int i=0; i<reducedNumElements; i++)
383          if (reducedOwner[i] == ownIndex)          if (reducedOwner[i] == ownIndex)
384              reducedIndexArray.push_back(i);              reducedIndexArray.push_back(i);
385    
386      for (int i=0; i<reducedCount; i++)      for (int i=0; i<reducedNumElements; i++)
387          if (reducedOwner[i] != ownIndex) {          if (reducedOwner[i] != ownIndex) {
388              numReducedGhostElements++;              numReducedGhostElements++;
389              reducedIndexArray.push_back(i);              reducedIndexArray.push_back(i);
# Line 353  void ElementData::prepareGhostIndices(in Line 393  void ElementData::prepareGhostIndices(in
393  //  //
394  //  //
395  //  //
396  void ElementData::handleGhostZones(int ownIndex)  void ElementData::reorderGhostZones(int ownIndex)
397  {  {
398      prepareGhostIndices(ownIndex);      IntVec indexArray, reducedIndexArray;
399        prepareGhostIndices(ownIndex, indexArray, reducedIndexArray);
400    
401      // move "ghost data" to the end of the arrays      // move "ghost data" to the end of the arrays
402      if (numGhostElements > 0) {      if (numGhostElements > 0) {
403          reorderArray(nodes, indexArray, nodesPerElement);          reorderArray(nodes, indexArray, nodesPerElement);
404          reorderArray(owner, indexArray, 1);          reorderArray(owner, indexArray, 1);
405          if (reducedCount == 0) {          if (reducedNumElements == 0) {
406              reorderArray(color, indexArray, 1);              reorderArray(color, indexArray, 1);
407              reorderArray(ID, indexArray, 1);              reorderArray(ID, indexArray, 1);
408              reorderArray(tag, indexArray, 1);              reorderArray(tag, indexArray, 1);
# Line 380  void ElementData::handleGhostZones(int o Line 421  void ElementData::handleGhostZones(int o
421  //  //
422  //  //
423  //  //
424  void ElementData::removeGhostZones()  void ElementData::removeGhostZones(int ownIndex)
425  {  {
426        reorderGhostZones(ownIndex);
427    
428      if (numGhostElements > 0) {      if (numGhostElements > 0) {
429          count -= numGhostElements;          numElements -= numGhostElements;
430          nodes.resize(count*nodesPerElement);          nodes.resize(numElements*nodesPerElement);
431          owner.resize(count);          owner.resize(numElements);
432          if (reducedCount == 0) {          if (reducedNumElements == 0) {
433              color.resize(count);              color.resize(numElements);
434              ID.resize(count);              ID.resize(numElements);
435              tag.resize(count);              tag.resize(numElements);
436          }          }
437          numGhostElements = 0;          numGhostElements = 0;
438      }      }
439    
440      if (numReducedGhostElements > 0) {      if (numReducedGhostElements > 0) {
441          reducedCount -= numReducedGhostElements;          reducedNumElements -= numReducedGhostElements;
442          reducedNodes.resize(reducedCount*reducedNodesPerElement);          reducedNodes.resize(reducedNumElements*reducedNodesPerElement);
443          reducedOwner.resize(reducedCount);          reducedOwner.resize(reducedNumElements);
444          color.resize(reducedCount);          color.resize(reducedNumElements);
445          ID.resize(reducedCount);          ID.resize(reducedNumElements);
446          tag.resize(reducedCount);          tag.resize(reducedNumElements);
447          numReducedGhostElements = 0;          numReducedGhostElements = 0;
448      }      }
449        buildMeshes();
450        if (numElements > 0)
451            fullMesh->removeGhostNodes(ownIndex);
452        if (reducedNumElements > 0)
453            reducedMesh->removeGhostNodes(ownIndex);
454  }  }
455    
456  //  //
# Line 410  void ElementData::removeGhostZones() Line 458  void ElementData::removeGhostZones()
458  //  //
459  void ElementData::buildMeshes()  void ElementData::buildMeshes()
460  {  {
461      if (count == 0)      // build a new mesh containing only the required nodes
462          return;      if (numElements > 0) {
463            fullMesh = NodeData_ptr(new NodeData(originalNodes, nodes, name));
     // use existing original mesh for Elements but build a new mesh  
     // containing only the required nodes for other types  
     if (name == "Elements") {  
         fullMesh = const_cast<Mesh*>(originalMesh);  
         fullMeshIsOriginalMesh = true;  
     } else {  
         // first: build a map of required IDs while translating  
         // the original node IDs at the same time  
         IntVec::iterator it;  
         IndexMap nodeID2idx;  
         IntVec nodeIDs;  
         size_t newIdx = 0;  
   
         for (it = nodes.begin(); it != nodes.end(); it++) {  
             IndexMap::iterator res = nodeID2idx.find(*it);  
             if (res == nodeID2idx.end()) {  
                 nodeIDs.push_back(*it);  
                 nodeID2idx[*it] = newIdx;  
                 *it = newIdx++;  
             } else {  
                 *it = res->second;  
             }  
         }  
   
         // second: use map to fill coordinates for the nodes  
         // newIdx contains the number of nodes  
         CoordArray coords;  
         const CoordArray& origCoords = originalMesh->getCoords();  
         for (int dim=0; dim < originalMesh->getNumDims(); dim++) {  
             float* c = new float[newIdx];  
             coords.push_back(c);  
             IndexMap::const_iterator mIt;  
             for (mIt = nodeID2idx.begin(); mIt != nodeID2idx.end(); mIt++)  
                 c[mIt->second] = origCoords[dim][mIt->first];  
         }  
   
         if (fullMesh)  
             delete fullMesh;  
   
         fullMesh = new Mesh(coords, originalMesh->getNumDims(), newIdx);  
         fullMesh->setName(name);  
         fullMesh->setIndexMap(nodeID2idx);  
         fullMesh->setNodeIDs(nodeIDs);  
     }  
464    
465  #ifdef _DEBUG  #ifdef _DEBUG
466      cout << fullMesh->getName() << " has " << fullMesh->getNumNodes()          cout << fullMesh->getName() << " has " << fullMesh->getNumNodes()
467          << " nodes, " << count << " elements" << endl;              << " nodes, " << numElements << " elements" << endl;
468  #endif  #endif
469        }
470    
471      // build a reduced mesh if necessary      // build a reduced mesh if necessary
472      if (reducedCount > 0) {      if (reducedNumElements > 0) {
473          IntVec::iterator it;          reducedMesh = NodeData_ptr(new NodeData(
474          IndexMap nodeID2idx;                      originalNodes, reducedNodes, "Reduced"+name));
         IntVec reducedNodeIDs;  
         size_t newIdx = 0;  
   
         for (it = reducedNodes.begin(); it != reducedNodes.end(); it++) {  
             int id = originalMesh->getNodeIDs()[*it];  
             IndexMap::iterator res = nodeID2idx.find(id);  
             if (res == nodeID2idx.end()) {  
                 reducedNodeIDs.push_back(id);  
                 nodeID2idx[id] = newIdx;  
                 *it = newIdx++;  
             } else {  
                 *it = res->second;  
             }  
         }  
   
         CoordArray coords;  
         const CoordArray& origCoords = originalMesh->getCoords();  
         for (int dim=0; dim < originalMesh->getNumDims(); dim++) {  
             float* c = new float[newIdx];  
             coords.push_back(c);  
             IndexMap::const_iterator mIt;  
             for (mIt = nodeID2idx.begin(); mIt != nodeID2idx.end(); mIt++) {  
                 IndexMap::const_iterator idx;  
                 idx = originalMesh->getIndexMap().find(mIt->first);  
                 c[mIt->second] = origCoords[dim][idx->second];  
             }  
         }  
         if (reducedMesh)  
             delete reducedMesh;  
           
         reducedMesh = new Mesh(coords, originalMesh->getNumDims(), newIdx);  
         reducedMesh->setName(string("Reduced") + name);  
         reducedMesh->setIndexMap(nodeID2idx);  
         reducedMesh->setNodeIDs(reducedNodeIDs);  
   
475  #ifdef _DEBUG  #ifdef _DEBUG
476          cout << reducedMesh->getName() << " has " << newIdx          cout << reducedMesh->getName() << " has " << newIdx
477              << " nodes, " << reducedCount << " elements" << endl;              << " nodes, " << reducedNumElements << " elements" << endl;
478  #endif  #endif
479      }      }
480  }  }
481    
482  #if HAVE_SILO  #if USE_SILO
483  //  //
484  //  //
485  //  //
# Line 532  inline int toSiloElementType(int type) Line 502  inline int toSiloElementType(int type)
502  //  //
503  bool ElementData::writeToSilo(DBfile* dbfile, const string& siloPath)  bool ElementData::writeToSilo(DBfile* dbfile, const string& siloPath)
504  {  {
505  #if HAVE_SILO  #if USE_SILO
506      if (count == 0)      if (numElements == 0)
507          return true;          return true;
508    
     const char* meshName;  
509      int numCells;      int numCells;
510      string varName, siloMeshName;      string varName, siloMeshNameStr;
511      int ret;      int ret;
512    
513      if (siloPath != "") {      if (siloPath != "") {
# Line 548  bool ElementData::writeToSilo(DBfile* db Line 517  bool ElementData::writeToSilo(DBfile* db
517      }      }
518    
519      // write out the full mesh in any case      // write out the full mesh in any case
520      if (siloPath == "/")      fullMesh->setSiloPath(siloPath);
521          siloMeshName = siloPath + fullMesh->getName();      siloMeshNameStr = fullMesh->getFullSiloName();
522      else      const char* siloMeshName = siloMeshNameStr.c_str();
523          siloMeshName = siloPath + string("/") + fullMesh->getName();      int arraylen = numElements * nodesPerElement;
   
     int arraylen = count * nodesPerElement;  
524      int eltype = toSiloElementType(type);      int eltype = toSiloElementType(type);
525    
526      varName = name + string("_zones");      varName = name + string("_zones");
527      ret = DBPutZonelist2(dbfile, varName.c_str(), count,      ret = DBPutZonelist2(dbfile, varName.c_str(), numElements,
528              originalMesh->getNumDims(), &nodes[0], arraylen, 0, 0,              fullMesh->getNumDims(), &nodes[0], arraylen, 0, 0,
529              numGhostElements, &eltype, &nodesPerElement, &count, 1, NULL);              numGhostElements, &eltype, &nodesPerElement, &numElements, 1, NULL);
530      if (ret == 0) {      if (ret == 0) {
531          CoordArray& coordbase = const_cast<CoordArray&>(fullMesh->getCoords());          CoordArray& coordbase = const_cast<CoordArray&>(fullMesh->getCoords());
532          ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(),          ret = DBPutUcdmesh(dbfile, siloMeshName,
533                  originalMesh->getNumDims(), NULL, &coordbase[0],                  fullMesh->getNumDims(), NULL, &coordbase[0],
534                  fullMesh->getNumNodes(), count, varName.c_str(),                  fullMesh->getNumNodes(), numElements, varName.c_str(),
535                  /*"facelist"*/NULL, DB_FLOAT, NULL);                  /*"facelist"*/NULL, DB_FLOAT, NULL);
536      }      }
537            
# Line 571  bool ElementData::writeToSilo(DBfile* db Line 539  bool ElementData::writeToSilo(DBfile* db
539      if (0) {      if (0) {
540          CoordArray& coordbase = const_cast<CoordArray&>(fullMesh->getCoords());          CoordArray& coordbase = const_cast<CoordArray&>(fullMesh->getCoords());
541          DBPutPointmesh(dbfile, "/pointmesh",          DBPutPointmesh(dbfile, "/pointmesh",
542                originalMesh->getNumDims(), &coordbase[0],                originalNodes->getNumDims(), &coordbase[0],
543                fullMesh->getNumNodes(), DB_FLOAT, NULL);                fullMesh->getNumNodes(), DB_FLOAT, NULL);
544      }      }
545    
# Line 579  bool ElementData::writeToSilo(DBfile* db Line 547  bool ElementData::writeToSilo(DBfile* db
547          return false;          return false;
548    
549      // decide whether to additionally write out the reduced mesh      // decide whether to additionally write out the reduced mesh
550      if (reducedCount > 0) {      if (reducedNumElements > 0) {
551          if (siloPath == "/")          reducedMesh->setSiloPath(siloPath);
552              siloMeshName = siloPath + reducedMesh->getName();          siloMeshNameStr = reducedMesh->getFullSiloName();
553          else          siloMeshName = siloMeshNameStr.c_str();
554              siloMeshName = siloPath + string("/") + reducedMesh->getName();          arraylen = reducedNumElements * reducedNodesPerElement;
         arraylen = reducedCount * reducedNodesPerElement;  
555          eltype = toSiloElementType(reducedType);          eltype = toSiloElementType(reducedType);
556          varName = string("Reduced") + name + string("_zones");          varName = string("Reduced") + name + string("_zones");
557          ret = DBPutZonelist2(dbfile, varName.c_str(), reducedCount,          ret = DBPutZonelist2(dbfile, varName.c_str(), reducedNumElements,
558                  originalMesh->getNumDims(), &reducedNodes[0], arraylen, 0, 0,                  originalNodes->getNumDims(), &reducedNodes[0], arraylen, 0, 0,
559                  numReducedGhostElements, &eltype, &reducedNodesPerElement,                  numReducedGhostElements, &eltype, &reducedNodesPerElement,
560                  &reducedCount, 1, NULL);                  &reducedNumElements, 1, NULL);
561          if (ret == 0) {          if (ret == 0) {
562              CoordArray& coordbase = const_cast<CoordArray&>(reducedMesh->getCoords());              CoordArray& coordbase = const_cast<CoordArray&>(reducedMesh->getCoords());
563              ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(),              ret = DBPutUcdmesh(dbfile, siloMeshName,
564                     originalMesh->getNumDims(), NULL, &coordbase[0],                     reducedMesh->getNumDims(), NULL, &coordbase[0],
565                     reducedMesh->getNumNodes(), reducedCount, varName.c_str(),                     reducedMesh->getNumNodes(), reducedNumElements, varName.c_str(),
566                     NULL, DB_FLOAT, NULL);                     NULL, DB_FLOAT, NULL);
567          }          }
568          if (ret != 0)          if (ret != 0)
569              return false;              return false;
570          numCells = reducedCount;          numCells = reducedNumElements;
571      } else {      } else {
572          numCells = count;          numCells = numElements;
573      }      }
     meshName = siloMeshName.c_str();  
574    
575      // finally, write out the element-centered variables on the correct mesh      // finally, write out the element-centered variables on the correct mesh
576      varName = name + string("_Color");      varName = name + string("_Color");
577      ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,      ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
578              (float*)&color[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);              (float*)&color[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
579      if (ret == 0) {      if (ret == 0) {
580          varName = name + string("_Id");          varName = name + string("_Id");
581          ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
582              (float*)&ID[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);              (float*)&ID[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
583      }      }
584      if (ret == 0) {      if (ret == 0) {
585          varName = name + string("_Owner");          varName = name + string("_Owner");
586          ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
587              (float*)&owner[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);              (float*)&owner[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
588      }      }
589      if (ret == 0) {      if (ret == 0) {
590          varName = name + string("_Tag");          varName = name + string("_Tag");
591          ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,          ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
592              (float*)&tag[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);              (float*)&tag[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
593      }      }
594        
595        if (name == "Elements") {
596            fullMesh->writeToSilo(dbfile);
597        }
598    
599      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
600      return (ret == 0);      return (ret == 0);
601    
602  #else // !HAVE_SILO  #else // !USE_SILO
603      return false;      return false;
604  #endif  #endif
605  }  }
# Line 771  FinleyElementInfo ElementData::getFinley Line 741  FinleyElementInfo ElementData::getFinley
741      return ret;      return ret;
742  }  }
743    
744  } // namespace EscriptReader  } // namespace escriptexport
745    
Legend:
Removed from v.2548  
changed lines
  Added in v.2810
  ViewVC Help
Powered by ViewVC 1.1.26