/[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 2187 by caltinay, Tue Dec 23 04:13:15 2008 UTC trunk/dataexporter/src/NodeData.cpp revision 2880 by caltinay, Thu Jan 28 01:21:30 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 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    //
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 54  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;
116        for (it = coords.begin(); it != coords.end(); it++)
117            delete[] *it;
118    }
119    
120    //
121    //
122    //
123    bool NodeData::initFromFinley(const Finley_NodeFile* finleyFile)
124  {  {
125        numDims = finleyFile->numDim;
126        numNodes = finleyFile->numNodes;
127    
128        int mpisize = finleyFile->MPIInfo->size;
129        int* iPtr = finleyFile->nodesDistribution->first_component;
130        nodeDist.clear();
131        nodeDist.insert(nodeDist.end(), mpisize+1, 0);
132        copy(iPtr, iPtr+mpisize+1, nodeDist.begin());
133    
134      CoordArray::iterator it;      CoordArray::iterator it;
135      for (it = coords.begin(); it != coords.end(); it++)      for (it = coords.begin(); it != coords.end(); it++)
136          delete[] *it;          delete[] *it;
137        coords.clear();
138        nodeID.clear();
139        nodeTag.clear();
140        nodeGDOF.clear();
141        nodeGNI.clear();
142        nodeGRDFI.clear();
143        nodeGRNI.clear();
144    
145        if (numNodes > 0) {
146            for (int i=0; i<numDims; i++) {
147                double* srcPtr = finleyFile->Coordinates + i;
148                float* c = new float[numNodes];
149                coords.push_back(c);
150                for (int j=0; j<numNodes; j++, srcPtr+=numDims) {
151                    *c++ = (float) *srcPtr;
152                }
153            }
154    
155            iPtr = finleyFile->Id;
156            nodeID.insert(nodeID.end(), numNodes, 0);
157            copy(iPtr, iPtr+numNodes, nodeID.begin());
158    
159            iPtr = finleyFile->Tag;
160            nodeTag.insert(nodeTag.end(), numNodes, 0);
161            copy(iPtr, iPtr+numNodes, nodeTag.begin());
162    
163            iPtr = finleyFile->globalDegreesOfFreedom;
164            nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
165            copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
166    
167            iPtr = finleyFile->globalNodesIndex;
168            nodeGNI.insert(nodeGNI.end(), numNodes, 0);
169            copy(iPtr, iPtr+numNodes, nodeGNI.begin());
170    
171            iPtr = finleyFile->globalReducedDOFIndex;
172            nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
173            copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
174    
175            iPtr = finleyFile->globalReducedNodesIndex;
176            nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
177            copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
178    
179        }
180        return true;
181  }  }
182    
183  //  //
184  //  //
185  //  //
186  bool Mesh::readFromNc(const string& ncFile)  bool NodeData::readFromNc(NcFile* ncFile)
187  {  {
188      NcError ncerr(NcError::silent_nonfatal);  #if USE_NETCDF
     NcFile* input;  
189      NcAtt* att;      NcAtt* att;
190      NcVar* var;      NcVar* var;
191    
192      input = new NcFile(ncFile.c_str());      att = ncFile->get_att("numDim");
     if (!input->is_valid()) {  
         cerr << "Could not open input file " << ncFile.c_str() << ".\n";  
         delete input;  
         return false;  
     }  
   
     att = input->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        nodeDist.clear();
202        nodeDist.insert(nodeDist.end(), mpisize+1, 0);
203        var = ncFile->get_var("Nodes_NodeDistribution");
204        var->get(&nodeDist[0], mpisize+1);
205    
206        CoordArray::iterator it;
207        for (it = coords.begin(); it != coords.end(); it++)
208            delete[] *it;
209      coords.clear();      coords.clear();
210      var = input->get_var("Nodes_Coordinates");      nodeID.clear();
211      for (int i=0; i<numDims; i++) {      nodeTag.clear();
212          float* c = new float[numNodes];      nodeGDOF.clear();
213          var->set_cur(0, i);      nodeGNI.clear();
214          var->get(c, numNodes, 1);      nodeGRDFI.clear();
215          coords.push_back(c);      nodeGRNI.clear();
216    
217        // Only attempt to read further if there are any nodes.
218        // Having no nodes is not an error.
219        if (numNodes > 0) {
220            var = ncFile->get_var("Nodes_Coordinates");
221            for (int i=0; i<numDims; i++) {
222                float* c = new float[numNodes];
223                var->set_cur(0, i);
224                var->get(c, numNodes, 1);
225                coords.push_back(c);
226            }
227    
228            nodeID.insert(nodeID.end(), numNodes, 0);
229            var = ncFile->get_var("Nodes_Id");
230            var->get(&nodeID[0], numNodes);
231    
232            nodeTag.insert(nodeTag.end(), numNodes, 0);
233            var = ncFile->get_var("Nodes_Tag");
234            var->get(&nodeTag[0], numNodes);
235    
236            nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
237            var = ncFile->get_var("Nodes_gDOF");
238            var->get(&nodeGDOF[0], numNodes);
239    
240            nodeGNI.insert(nodeGNI.end(), numNodes, 0);
241            var = ncFile->get_var("Nodes_gNI");
242            var->get(&nodeGNI[0], numNodes);
243    
244            nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
245            var = ncFile->get_var("Nodes_grDfI");
246            var->get(&nodeGRDFI[0], numNodes);
247    
248            nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
249            var = ncFile->get_var("Nodes_grNI");
250            var->get(&nodeGRNI[0], numNodes);
251      }      }
252    
253      nodeID.clear();      return true;
254      nodeID.insert(nodeID.end(), numNodes, 0);  #else // !USE_NETCDF
255      var = input->get_var("Nodes_Id");      return false;
256      var->get(&nodeID[0], numNodes);  #endif
257    }
258    
259      buildIndexMap();  //
260    //
261    //
262    const IntVec& NodeData::getVarDataByName(const string& name) const
263    {
264        if (name == "Nodes_Id")
265            return nodeID;
266        else if (name == "Nodes_Tag")
267            return nodeTag;
268        else if (name == "Nodes_gDOF")
269            return nodeGDOF;
270        else if (name == "Nodes_gNI")
271            return nodeGNI;
272        else if (name == "Nodes_grDfI")
273            return nodeGRDFI;
274        else if (name == "Nodes_grNI")
275            return nodeGRNI;
276        else
277            throw "Invalid variable name";
278    }
279    
280      delete input;  //
281      return true;  //
282    //
283    StringVec NodeData::getVarNames() const
284    {
285        StringVec res;
286        res.push_back("Nodes_Id");
287        res.push_back("Nodes_Tag");
288        res.push_back("Nodes_gDOF");
289        res.push_back("Nodes_gNI");
290        res.push_back("Nodes_grDfI");
291        res.push_back("Nodes_grNI");
292        return res;
293  }  }
294    
295  //  //
296  //  //
297  //  //
298  bool Mesh::writeToSilo(DBfile* dbfile, const string& pathInSilo)  bool NodeData::writeToSilo(DBfile* dbfile)
299  {  {
300  #if HAVE_SILO  #if USE_SILO
301        if (numNodes == 0)
302            return true;
303    
304      int ret;      int ret;
305    
306      if (pathInSilo != "") {      if (siloPath != "") {
         siloPath = pathInSilo;  
307          ret = DBSetDir(dbfile, siloPath.c_str());          ret = DBSetDir(dbfile, siloPath.c_str());
308          if (ret != 0)          if (ret != 0)
309              return false;              return false;
     } else {  
         siloPath = "/";  
310      }      }
311          string siloMeshName = getFullSiloName();
312      // Write node-centered variable  
313      ret = DBPutUcdvar1(dbfile, "Nodes_Id", getFullSiloName().c_str(),      // Write node-centered variables
314        ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
315              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
316    
317        if (ret == 0)
318            ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
319                    (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
320                    DB_NODECENT, NULL);
321        if (ret == 0)
322            ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
323                    (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
324                    DB_NODECENT, NULL);
325        if (ret == 0)
326            ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
327                    (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
328                    DB_NODECENT, NULL);
329        if (ret == 0)
330            ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
331                    (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
332                    DB_NODECENT, NULL);
333        if (ret == 0)
334            ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
335                    (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
336                    DB_NODECENT, NULL);
337    
338      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
339      return (ret == 0);      return (ret == 0);
340    
341  #else // !HAVE_SILO  #else // !USE_SILO
342      return false;      return false;
343  #endif  #endif
344  }  }
345    
346  } // namespace EscriptReader  } // namespace escriptexport
347    

Legend:
Removed from v.2187  
changed lines
  Added in v.2880

  ViewVC Help
Powered by ViewVC 1.1.26