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

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

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

trunk/tools/libescriptreader/src/escriptreader/Mesh.cpp revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC trunk/dataexporter/src/NodeData.cpp revision 2810 by caltinay, Mon Dec 7 04:13:49 2009 UTC
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
14  //  #include <escriptexport/NodeData.h>
15  // Mesh.cpp  
16  //  extern "C" {
17  #include <escriptreader/Mesh.h>  #include <finley/Mesh.h>
18  #include <escriptreader/ElementData.h>  #include <finley/NodeFile.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    
29  using namespace std;  using namespace std;
30    
31  namespace EscriptReader {  namespace escriptexport {
32    
33  //  //
34    // Constructor with name
35  //  //
36  //  NodeData::NodeData(const string& meshName) :
37  Mesh::Mesh(CoordArray c, int nDims, int nNodes) :      numDims(0), numNodes(0), name(meshName)
     coords(c), numDims(nDims), numNodes(nNodes), name("Mesh")  
38  {  {
39  }  }
40    
41  //  //
42  //  //
43  //  //
44  Mesh::Mesh(const Mesh& m)  NodeData::NodeData(NodeData_ptr fullNodes, IntVec& requiredNodes,
45                       const string& meshName) :
46        name(meshName)
47    {
48        numDims = fullNodes->numDims;
49        nodeDist = fullNodes->nodeDist;
50    
51        // first: find the unique set of required nodes and their IDs while
52        // updating the contents of requiredNodes at the same time
53        // requiredNodes contains node indices (not IDs!)
54        IntVec::iterator it;
55        IndexMap indexMap; // maps old index to new index
56        size_t newIndex = 0;
57    
58        for (it = requiredNodes.begin(); it != requiredNodes.end(); it++) {
59            IndexMap::iterator res = indexMap.find(*it);
60            if (res == indexMap.end()) {
61                nodeID.push_back(fullNodes->nodeID[*it]);
62                nodeTag.push_back(fullNodes->nodeTag[*it]);
63                nodeGDOF.push_back(fullNodes->nodeGDOF[*it]);
64                nodeGNI.push_back(fullNodes->nodeGNI[*it]);
65                nodeGRDFI.push_back(fullNodes->nodeGRDFI[*it]);
66                nodeGRNI.push_back(fullNodes->nodeGRNI[*it]);
67                indexMap[*it] = newIndex;
68                *it = newIndex++;
69            } else {
70                *it = res->second;
71            }
72        }
73    
74        // second: now that we know how many nodes we need use the map to fill
75        // the coordinates
76        numNodes = newIndex;
77        for (int dim=0; dim<numDims; dim++) {
78            const float* origC = fullNodes->coords[dim];
79            float* c = new float[numNodes];
80            coords.push_back(c);
81            IndexMap::const_iterator mIt;
82            for (mIt = indexMap.begin(); mIt != indexMap.end(); mIt++)
83                c[mIt->second] = origC[mIt->first];
84        }
85    }
86    
87    //
88    // Copy constructor
89    //
90    NodeData::NodeData(const NodeData& m)
91  {  {
92      numDims = m.numDims;      numDims = m.numDims;
93      numNodes = m.numNodes;      numNodes = m.numNodes;
94      nodeID = m.nodeID;      nodeID = m.nodeID;
95      nodeID2idx = m.nodeID2idx;      nodeTag = m.nodeTag;
96        nodeGDOF = m.nodeGDOF;
97        nodeGNI = m.nodeGNI;
98        nodeGRDFI = m.nodeGRDFI;
99        nodeGRNI = m.nodeGRNI;
100        nodeDist = m.nodeDist;
101      name = m.name;      name = m.name;
     siloPath = m.siloPath;  
102      for (int i=0; i<numDims; i++) {      for (int i=0; i<numDims; i++) {
103          float* c = new float[numNodes];          float* c = new float[numNodes];
104          copy(m.coords[i], m.coords[i]+numNodes, c);          copy(m.coords[i], m.coords[i]+numNodes, c);
# Line 54  Mesh::Mesh(const Mesh& m) Line 109  Mesh::Mesh(const Mesh& m)
109  //  //
110  //  //
111  //  //
112  Mesh::~Mesh()  NodeData::~NodeData()
113  {  {
114      CoordArray::iterator it;      CoordArray::iterator it;
115      for (it = coords.begin(); it != coords.end(); it++)      for (it = coords.begin(); it != coords.end(); it++)
# Line 64  Mesh::~Mesh() Line 119  Mesh::~Mesh()
119  //  //
120  //  //
121  //  //
122  bool Mesh::readFromNc(const string& ncFile)  bool NodeData::initFromFinley(const Finley_NodeFile* finleyFile)
123  {  {
124      NcError ncerr(NcError::silent_nonfatal);      numDims = finleyFile->numDim;
125      NcFile* input;      numNodes = finleyFile->numNodes;
126      NcAtt* att;  
127      NcVar* var;      CoordArray::iterator it;
128        for (it = coords.begin(); it != coords.end(); it++)
129            delete[] *it;
130        coords.clear();
131    
132        if (numNodes > 0) {
133            for (int i=0; i<numDims; i++) {
134                double* srcPtr = finleyFile->Coordinates + i;
135                float* c = new float[numNodes];
136                coords.push_back(c);
137                for (int j=0; j<numNodes; j++, srcPtr+=numDims) {
138                    *c++ = (float) *srcPtr;
139                }
140            }
141    
142            int* iPtr;
143    
144      input = new NcFile(ncFile.c_str());          iPtr = finleyFile->Id;
145      if (!input->is_valid()) {          nodeID.clear();
146          cerr << "Could not open input file " << ncFile.c_str() << ".\n";          nodeID.insert(nodeID.end(), numNodes, 0);
147          delete input;          copy(iPtr, iPtr+numNodes, nodeID.begin());
148          return false;  
149            iPtr = finleyFile->Tag;
150            nodeTag.clear();
151            nodeTag.insert(nodeTag.end(), numNodes, 0);
152            copy(iPtr, iPtr+numNodes, nodeTag.begin());
153    
154            iPtr = finleyFile->globalDegreesOfFreedom;
155            nodeGDOF.clear();
156            nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
157            copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
158    
159            iPtr = finleyFile->globalNodesIndex;
160            nodeGNI.clear();
161            nodeGNI.insert(nodeGNI.end(), numNodes, 0);
162            copy(iPtr, iPtr+numNodes, nodeGNI.begin());
163    
164            iPtr = finleyFile->globalReducedDOFIndex;
165            nodeGRDFI.clear();
166            nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
167            copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
168    
169            iPtr = finleyFile->globalReducedNodesIndex;
170            nodeGRNI.clear();
171            nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
172            copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
173    
174            int mpisize = finleyFile->MPIInfo->size;
175            iPtr = finleyFile->nodesDistribution->first_component;
176            nodeDist.clear();
177            nodeDist.insert(nodeDist.end(), mpisize+1, 0);
178            copy(iPtr, iPtr+mpisize+1, nodeDist.begin());
179      }      }
180        return true;
181    }
182    
183      att = input->get_att("numDim");  //
184    //
185    //
186    bool NodeData::readFromNc(NcFile* ncFile)
187    {
188    #if USE_NETCDF
189        NcAtt* att;
190        NcVar* var;
191    
192        att = ncFile->get_att("numDim");
193      numDims = att->as_int(0);      numDims = att->as_int(0);
194    
195      att = input->get_att("numNodes");      att = ncFile->get_att("numNodes");
196      numNodes = att->as_int(0);      numNodes = att->as_int(0);
197    
198        att = ncFile->get_att("mpi_size");
199        int mpisize = att->as_int(0);
200    
201        CoordArray::iterator it;
202        for (it = coords.begin(); it != coords.end(); it++)
203            delete[] *it;
204      coords.clear();      coords.clear();
205      var = input->get_var("Nodes_Coordinates");      var = ncFile->get_var("Nodes_Coordinates");
206      for (int i=0; i<numDims; i++) {      for (int i=0; i<numDims; i++) {
207          float* c = new float[numNodes];          float* c = new float[numNodes];
208          var->set_cur(0, i);          var->set_cur(0, i);
# Line 95  bool Mesh::readFromNc(const string& ncFi Line 212  bool Mesh::readFromNc(const string& ncFi
212    
213      nodeID.clear();      nodeID.clear();
214      nodeID.insert(nodeID.end(), numNodes, 0);      nodeID.insert(nodeID.end(), numNodes, 0);
215      var = input->get_var("Nodes_Id");      var = ncFile->get_var("Nodes_Id");
216      var->get(&nodeID[0], numNodes);      var->get(&nodeID[0], numNodes);
217    
218      buildIndexMap();      nodeTag.clear();
219        nodeTag.insert(nodeTag.end(), numNodes, 0);
220        var = ncFile->get_var("Nodes_Tag");
221        var->get(&nodeTag[0], numNodes);
222    
223        nodeGDOF.clear();
224        nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
225        var = ncFile->get_var("Nodes_gDOF");
226        var->get(&nodeGDOF[0], numNodes);
227    
228        nodeGNI.clear();
229        nodeGNI.insert(nodeGNI.end(), numNodes, 0);
230        var = ncFile->get_var("Nodes_gNI");
231        var->get(&nodeGNI[0], numNodes);
232    
233        nodeGRDFI.clear();
234        nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
235        var = ncFile->get_var("Nodes_grDfI");
236        var->get(&nodeGRDFI[0], numNodes);
237    
238        nodeGRNI.clear();
239        nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
240        var = ncFile->get_var("Nodes_grNI");
241        var->get(&nodeGRNI[0], numNodes);
242    
243        nodeDist.clear();
244        nodeDist.insert(nodeDist.end(), mpisize+1, 0);
245        var = ncFile->get_var("Nodes_NodeDistribution");
246        var->get(&nodeDist[0], mpisize+1);
247    
     delete input;  
248      return true;      return true;
249    #else // !USE_NETCDF
250        return false;
251    #endif
252    }
253    
254    //
255    //
256    //
257    const IntVec& NodeData::getVarDataByName(const std::string& name) const
258    {
259        if (name == "Nodes_Id")
260            return nodeID;
261        else if (name == "Nodes_Tag")
262            return nodeTag;
263        else if (name == "Nodes_gDOF")
264            return nodeGDOF;
265        else if (name == "Nodes_gNI")
266            return nodeGNI;
267        else if (name == "Nodes_grDfI")
268            return nodeGRDFI;
269        else if (name == "Nodes_grNI")
270            return nodeGRNI;
271        else
272            throw "Invalid variable name";
273    }
274    
275    //
276    //
277    //
278    StringVec NodeData::getVarNames() const
279    {
280        StringVec res;
281        if (numNodes > 0) {
282            res.push_back("Nodes_Id");
283            res.push_back("Nodes_Tag");
284            res.push_back("Nodes_gDOF");
285            res.push_back("Nodes_gNI");
286            res.push_back("Nodes_grDfI");
287            res.push_back("Nodes_grNI");
288        }
289        return res;
290  }  }
291    
292  //  //
293  //  //
294  //  //
295  bool Mesh::writeToSilo(DBfile* dbfile, const string& pathInSilo)  void NodeData::removeGhostNodes(int ownIndex)
296  {  {
297  #if HAVE_SILO      if (nodeDist.empty() || ownIndex > nodeDist.size()-1)
298            return;
299    
300        int firstId = nodeDist[ownIndex];
301        int lastId = nodeDist[ownIndex+1];
302    
303        numNodes = lastId-firstId;
304    
305        CoordArray newCoords;
306        CoordArray::iterator it;
307        for (int i=0; i<numDims; i++) {
308            float* c = new float[numNodes];
309            newCoords.push_back(c);
310        }
311    
312        IntVec newNodeID, newNodeTag;
313        IntVec newNodeGDOF, newNodeGNI, newNodeGRDFI, newNodeGRNI;
314    
315        IndexMap nodeID2idx = getIndexMap();
316        
317        for (int i=firstId; i<lastId; i++) {
318            int idx = nodeID2idx[i];
319            for (int j=0; j<numDims; j++) {
320                newCoords[j][i-firstId] = coords[j][idx];
321            }
322            newNodeID.push_back(i);
323            newNodeTag.push_back(nodeTag[idx]);
324            newNodeGDOF.push_back(nodeGDOF[idx]);
325            newNodeGNI.push_back(nodeGNI[idx]);
326            newNodeGRDFI.push_back(nodeGRDFI[idx]);
327            newNodeGRNI.push_back(nodeGRNI[idx]);
328        }
329    
330        for (it = coords.begin(); it != coords.end(); it++)
331            delete[] *it;
332    
333        coords = newCoords;
334        nodeID = newNodeID;
335        nodeTag = newNodeTag;
336        nodeGDOF = newNodeGDOF;
337        nodeGNI = newNodeGNI;
338        nodeGRDFI = newNodeGRDFI;
339        nodeGRNI = newNodeGRNI;
340    }
341    
342    //
343    //
344    //
345    bool NodeData::writeToSilo(DBfile* dbfile)
346    {
347    #if USE_SILO
348        if (numNodes == 0)
349            return true;
350    
351      int ret;      int ret;
352    
353      if (pathInSilo != "") {      if (siloPath != "") {
         siloPath = pathInSilo;  
354          ret = DBSetDir(dbfile, siloPath.c_str());          ret = DBSetDir(dbfile, siloPath.c_str());
355          if (ret != 0)          if (ret != 0)
356              return false;              return false;
     } else {  
         siloPath = "/";  
357      }      }
358          string siloMeshName = getFullSiloName();
359      // Write node-centered variable  
360      ret = DBPutUcdvar1(dbfile, "Nodes_Id", getFullSiloName().c_str(),      // Write node-centered variables
361        ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
362              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
363    
364        if (ret == 0)
365            ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
366                    (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
367                    DB_NODECENT, NULL);
368        if (ret == 0)
369            ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
370                    (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
371                    DB_NODECENT, NULL);
372        if (ret == 0)
373            ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
374                    (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
375                    DB_NODECENT, NULL);
376        if (ret == 0)
377            ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
378                    (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
379                    DB_NODECENT, NULL);
380        if (ret == 0)
381            ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
382                    (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
383                    DB_NODECENT, NULL);
384    
385      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
386      return (ret == 0);      return (ret == 0);
387    
388  #else // !HAVE_SILO  #else // !USE_SILO
389      return false;      return false;
390  #endif  #endif
391  }  }
392    
393  } // namespace EscriptReader  } // namespace escriptexport
394    

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

  ViewVC Help
Powered by ViewVC 1.1.26