/[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 2183 by caltinay, Fri Dec 19 03:52:50 2008 UTC trunk/dataexporter/src/NodeData.cpp revision 2834 by caltinay, Thu Jan 7 06:06:56 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 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 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 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    //
89    // Copy constructor
90    //
91    NodeData::NodeData(const NodeData& m)
92  {  {
93      numDims = m.numDims;      numDims = m.numDims;
94      numNodes = m.numNodes;      numNodes = m.numNodes;
95      nodeID = m.nodeID;      nodeID = m.nodeID;
96      nodeID2idx = m.nodeID2idx;      nodeTag = m.nodeTag;
97        nodeGDOF = m.nodeGDOF;
98        nodeGNI = m.nodeGNI;
99        nodeGRDFI = m.nodeGRDFI;
100        nodeGRNI = m.nodeGRNI;
101        nodeDist = m.nodeDist;
102      name = m.name;      name = m.name;
     siloPath = m.siloPath;  
103      for (int i=0; i<numDims; i++) {      for (int i=0; i<numDims; i++) {
104          float* c = new float[numNodes];          float* c = new float[numNodes];
105          copy(m.coords[i], m.coords[i]+numNodes, c);          copy(m.coords[i], m.coords[i]+numNodes, c);
# Line 52  Mesh::Mesh(const Mesh& m) Line 110  Mesh::Mesh(const Mesh& m)
110  //  //
111  //  //
112  //  //
113  Mesh::~Mesh()  NodeData::~NodeData()
114  {  {
115      CoordArray::iterator it;      CoordArray::iterator it;
116      for (it = coords.begin(); it != coords.end(); it++)      for (it = coords.begin(); it != coords.end(); it++)
# Line 62  Mesh::~Mesh() Line 120  Mesh::~Mesh()
120  //  //
121  //  //
122  //  //
123  bool Mesh::readFromNc(const string& ncFile)  bool NodeData::initFromFinley(const Finley_NodeFile* finleyFile)
124  {  {
125      NcError ncerr(NcError::silent_nonfatal);      numDims = finleyFile->numDim;
126      NcFile* input;      numNodes = finleyFile->numNodes;
127      NcAtt* att;  
128      NcVar* var;      CoordArray::iterator it;
129        for (it = coords.begin(); it != coords.end(); it++)
130            delete[] *it;
131        coords.clear();
132    
133        if (numNodes > 0) {
134            for (int i=0; i<numDims; i++) {
135                double* srcPtr = finleyFile->Coordinates + i;
136                float* c = new float[numNodes];
137                coords.push_back(c);
138                for (int j=0; j<numNodes; j++, srcPtr+=numDims) {
139                    *c++ = (float) *srcPtr;
140                }
141            }
142    
143            int* iPtr;
144    
145      input = new NcFile(ncFile.c_str());          iPtr = finleyFile->Id;
146      if (!input->is_valid()) {          nodeID.clear();
147          cerr << "Could not open input file " << ncFile.c_str() << ".\n";          nodeID.insert(nodeID.end(), numNodes, 0);
148          delete input;          copy(iPtr, iPtr+numNodes, nodeID.begin());
149          return false;  
150            iPtr = finleyFile->Tag;
151            nodeTag.clear();
152            nodeTag.insert(nodeTag.end(), numNodes, 0);
153            copy(iPtr, iPtr+numNodes, nodeTag.begin());
154    
155            iPtr = finleyFile->globalDegreesOfFreedom;
156            nodeGDOF.clear();
157            nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
158            copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
159    
160            iPtr = finleyFile->globalNodesIndex;
161            nodeGNI.clear();
162            nodeGNI.insert(nodeGNI.end(), numNodes, 0);
163            copy(iPtr, iPtr+numNodes, nodeGNI.begin());
164    
165            iPtr = finleyFile->globalReducedDOFIndex;
166            nodeGRDFI.clear();
167            nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
168            copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
169    
170            iPtr = finleyFile->globalReducedNodesIndex;
171            nodeGRNI.clear();
172            nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
173            copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
174    
175            int mpisize = finleyFile->MPIInfo->size;
176            iPtr = finleyFile->nodesDistribution->first_component;
177            nodeDist.clear();
178            nodeDist.insert(nodeDist.end(), mpisize+1, 0);
179            copy(iPtr, iPtr+mpisize+1, nodeDist.begin());
180      }      }
181        return true;
182    }
183    
184      att = input->get_att("numDim");  //
185    //
186    //
187    bool NodeData::readFromNc(NcFile* ncFile)
188    {
189    #if USE_NETCDF
190        NcAtt* att;
191        NcVar* var;
192    
193        att = ncFile->get_att("numDim");
194      numDims = att->as_int(0);      numDims = att->as_int(0);
195    
196      att = input->get_att("numNodes");      att = ncFile->get_att("numNodes");
197      numNodes = att->as_int(0);      numNodes = att->as_int(0);
198    
199        att = ncFile->get_att("mpi_size");
200        int mpisize = att->as_int(0);
201    
202        CoordArray::iterator it;
203        for (it = coords.begin(); it != coords.end(); it++)
204            delete[] *it;
205      coords.clear();      coords.clear();
206      var = input->get_var("Nodes_Coordinates");      var = ncFile->get_var("Nodes_Coordinates");
207      for (int i=0; i<numDims; i++) {      for (int i=0; i<numDims; i++) {
208          float* c = new float[numNodes];          float* c = new float[numNodes];
209          var->set_cur(0, i);          var->set_cur(0, i);
# Line 93  bool Mesh::readFromNc(const string& ncFi Line 213  bool Mesh::readFromNc(const string& ncFi
213    
214      nodeID.clear();      nodeID.clear();
215      nodeID.insert(nodeID.end(), numNodes, 0);      nodeID.insert(nodeID.end(), numNodes, 0);
216      var = input->get_var("Nodes_Id");      var = ncFile->get_var("Nodes_Id");
217      var->get(&nodeID[0], numNodes);      var->get(&nodeID[0], numNodes);
218    
219      buildIndexMap();      nodeTag.clear();
220        nodeTag.insert(nodeTag.end(), numNodes, 0);
221        var = ncFile->get_var("Nodes_Tag");
222        var->get(&nodeTag[0], numNodes);
223    
224        nodeGDOF.clear();
225        nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
226        var = ncFile->get_var("Nodes_gDOF");
227        var->get(&nodeGDOF[0], numNodes);
228    
229        nodeGNI.clear();
230        nodeGNI.insert(nodeGNI.end(), numNodes, 0);
231        var = ncFile->get_var("Nodes_gNI");
232        var->get(&nodeGNI[0], numNodes);
233    
234        nodeGRDFI.clear();
235        nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
236        var = ncFile->get_var("Nodes_grDfI");
237        var->get(&nodeGRDFI[0], numNodes);
238    
239        nodeGRNI.clear();
240        nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
241        var = ncFile->get_var("Nodes_grNI");
242        var->get(&nodeGRNI[0], numNodes);
243    
244        nodeDist.clear();
245        nodeDist.insert(nodeDist.end(), mpisize+1, 0);
246        var = ncFile->get_var("Nodes_NodeDistribution");
247        var->get(&nodeDist[0], mpisize+1);
248    
     delete input;  
249      return true;      return true;
250    #else // !USE_NETCDF
251        return false;
252    #endif
253  }  }
254    
255  //  //
256  //  //
257  //  //
258  bool Mesh::writeToSilo(DBfile* dbfile, const string& pathInSilo)  const IntVec& NodeData::getVarDataByName(const string& name) const
259  {  {
260  #if HAVE_SILO      if (name == "Nodes_Id")
261            return nodeID;
262        else if (name == "Nodes_Tag")
263            return nodeTag;
264        else if (name == "Nodes_gDOF")
265            return nodeGDOF;
266        else if (name == "Nodes_gNI")
267            return nodeGNI;
268        else if (name == "Nodes_grDfI")
269            return nodeGRDFI;
270        else if (name == "Nodes_grNI")
271            return nodeGRNI;
272        else
273            throw "Invalid variable name";
274    }
275    
276    //
277    //
278    //
279    StringVec NodeData::getVarNames() const
280    {
281        StringVec res;
282        if (numNodes > 0) {
283            res.push_back("Nodes_Id");
284            res.push_back("Nodes_Tag");
285            res.push_back("Nodes_gDOF");
286            res.push_back("Nodes_gNI");
287            res.push_back("Nodes_grDfI");
288            res.push_back("Nodes_grNI");
289        }
290        return res;
291    }
292    
293    //
294    //
295    //
296    void NodeData::removeGhostNodes(int ownIndex)
297    {
298        if (nodeDist.empty() || ownIndex > nodeDist.size()-1)
299            return;
300    
301        int firstId = nodeDist[ownIndex];
302        int lastId = nodeDist[ownIndex+1];
303    
304        // no ghost nodes
305        if (lastId-firstId == numNodes)
306            return;
307    
308        // we have at most lastId-firstId nodes, it could be less however if
309        // nodes were culled already
310        numNodes = lastId-firstId;
311    
312        CoordArray newCoords;
313        CoordArray::iterator it;
314        for (int i=0; i<numDims; i++) {
315            float* c = new float[numNodes];
316            newCoords.push_back(c);
317        }
318    
319        IntVec newNodeID, newNodeTag;
320        IntVec newNodeGDOF, newNodeGNI, newNodeGRDFI, newNodeGRNI;
321    
322        IndexMap nodeID2idx = getIndexMap();
323        int destIdx = 0;
324        for (int i=firstId; i<lastId; i++) {
325            IndexMap::iterator it = nodeID2idx.find(i);
326            if (it == nodeID2idx.end()) {
327                continue;
328            }
329            int idx = it->second;
330            for (int dim=0; dim<numDims; dim++) {
331                newCoords[dim][destIdx] = coords[dim][idx];
332            }
333            destIdx++;
334            newNodeID.push_back(i);
335            newNodeTag.push_back(nodeTag[idx]);
336            newNodeGDOF.push_back(nodeGDOF[idx]);
337            newNodeGNI.push_back(nodeGNI[idx]);
338            newNodeGRDFI.push_back(nodeGRDFI[idx]);
339            newNodeGRNI.push_back(nodeGRNI[idx]);
340        }
341    
342        numNodes = destIdx;
343    
344        for (it = coords.begin(); it != coords.end(); it++)
345            delete[] *it;
346    
347        coords = newCoords;
348        nodeID = newNodeID;
349        nodeTag = newNodeTag;
350        nodeGDOF = newNodeGDOF;
351        nodeGNI = newNodeGNI;
352        nodeGRDFI = newNodeGRDFI;
353        nodeGRNI = newNodeGRNI;
354    }
355    
356    //
357    //
358    //
359    bool NodeData::writeToSilo(DBfile* dbfile)
360    {
361    #if USE_SILO
362        if (numNodes == 0)
363            return true;
364    
365      int ret;      int ret;
366    
367      if (pathInSilo != "") {      if (siloPath != "") {
         siloPath = pathInSilo;  
368          ret = DBSetDir(dbfile, siloPath.c_str());          ret = DBSetDir(dbfile, siloPath.c_str());
369          if (ret != 0)          if (ret != 0)
370              return false;              return false;
     } else {  
         siloPath = "/";  
371      }      }
372          string siloMeshName = getFullSiloName();
373      // Write node-centered variable  
374      ret = DBPutUcdvar1(dbfile, "Nodes_Id", getFullSiloName().c_str(),      // Write node-centered variables
375        ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
376              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
377    
378        if (ret == 0)
379            ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
380                    (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
381                    DB_NODECENT, NULL);
382        if (ret == 0)
383            ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
384                    (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
385                    DB_NODECENT, NULL);
386        if (ret == 0)
387            ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
388                    (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
389                    DB_NODECENT, NULL);
390        if (ret == 0)
391            ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
392                    (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
393                    DB_NODECENT, NULL);
394        if (ret == 0)
395            ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
396                    (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
397                    DB_NODECENT, NULL);
398    
399      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
400      return (ret == 0);      return (ret == 0);
401    
402  #else // !HAVE_SILO  #else // !USE_SILO
403      return false;      return false;
404  #endif  #endif
405  }  }
406    
407    } // namespace escriptexport
408    

Legend:
Removed from v.2183  
changed lines
  Added in v.2834

  ViewVC Help
Powered by ViewVC 1.1.26