/[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 2853 by caltinay, Mon Jan 18 01:46:05 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;      CoordArray::iterator it;
116      for (it = coords.begin(); it != coords.end(); it++)      for (it = coords.begin(); it != coords.end(); it++)
# Line 64  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;      int mpisize = finleyFile->MPIInfo->size;
129        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;
135        for (it = coords.begin(); it != coords.end(); it++)
136            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            int* iPtr;
156    
157      input = new NcFile(ncFile.c_str());          iPtr = finleyFile->Id;
158      if (!input->is_valid()) {          nodeID.insert(nodeID.end(), numNodes, 0);
159          cerr << "Could not open input file " << ncFile.c_str() << ".\n";          copy(iPtr, iPtr+numNodes, nodeID.begin());
160          delete input;  
161          return false;          iPtr = finleyFile->Tag;
162            nodeTag.insert(nodeTag.end(), numNodes, 0);
163            copy(iPtr, iPtr+numNodes, nodeTag.begin());
164    
165            iPtr = finleyFile->globalDegreesOfFreedom;
166            nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
167            copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
168    
169            iPtr = finleyFile->globalNodesIndex;
170            nodeGNI.insert(nodeGNI.end(), numNodes, 0);
171            copy(iPtr, iPtr+numNodes, nodeGNI.begin());
172    
173            iPtr = finleyFile->globalReducedDOFIndex;
174            nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
175            copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
176    
177            iPtr = finleyFile->globalReducedNodesIndex;
178            nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
179            copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
180    
181      }      }
182        return true;
183    }
184    
185      att = input->get_att("numDim");  //
186    //
187    //
188    bool NodeData::readFromNc(NcFile* ncFile)
189    {
190    #if USE_NETCDF
191        NcAtt* att;
192        NcVar* var;
193    
194        att = ncFile->get_att("numDim");
195      numDims = att->as_int(0);      numDims = att->as_int(0);
196    
197      att = input->get_att("numNodes");      att = ncFile->get_att("numNodes");
198      numNodes = att->as_int(0);      numNodes = att->as_int(0);
199    
200        att = ncFile->get_att("mpi_size");
201        int mpisize = att->as_int(0);
202    
203        nodeDist.clear();
204        nodeDist.insert(nodeDist.end(), mpisize+1, 0);
205        var = ncFile->get_var("Nodes_NodeDistribution");
206        var->get(&nodeDist[0], mpisize+1);
207    
208        CoordArray::iterator it;
209        for (it = coords.begin(); it != coords.end(); it++)
210            delete[] *it;
211      coords.clear();      coords.clear();
212      var = input->get_var("Nodes_Coordinates");      nodeID.clear();
213      for (int i=0; i<numDims; i++) {      nodeTag.clear();
214          float* c = new float[numNodes];      nodeGDOF.clear();
215          var->set_cur(0, i);      nodeGNI.clear();
216          var->get(c, numNodes, 1);      nodeGRDFI.clear();
217          coords.push_back(c);      nodeGRNI.clear();
218    
219        // Only attempt to read further if there are any nodes.
220        // Having no nodes is not an error.
221        if (numNodes > 0) {
222            var = ncFile->get_var("Nodes_Coordinates");
223            for (int i=0; i<numDims; i++) {
224                float* c = new float[numNodes];
225                var->set_cur(0, i);
226                var->get(c, numNodes, 1);
227                coords.push_back(c);
228            }
229    
230            nodeID.insert(nodeID.end(), numNodes, 0);
231            var = ncFile->get_var("Nodes_Id");
232            var->get(&nodeID[0], numNodes);
233    
234            nodeTag.insert(nodeTag.end(), numNodes, 0);
235            var = ncFile->get_var("Nodes_Tag");
236            var->get(&nodeTag[0], numNodes);
237    
238            nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
239            var = ncFile->get_var("Nodes_gDOF");
240            var->get(&nodeGDOF[0], numNodes);
241    
242            nodeGNI.insert(nodeGNI.end(), numNodes, 0);
243            var = ncFile->get_var("Nodes_gNI");
244            var->get(&nodeGNI[0], numNodes);
245    
246            nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
247            var = ncFile->get_var("Nodes_grDfI");
248            var->get(&nodeGRDFI[0], numNodes);
249    
250            nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
251            var = ncFile->get_var("Nodes_grNI");
252            var->get(&nodeGRNI[0], numNodes);
253      }      }
254    
255      nodeID.clear();      return true;
256      nodeID.insert(nodeID.end(), numNodes, 0);  #else // !USE_NETCDF
257      var = input->get_var("Nodes_Id");      return false;
258      var->get(&nodeID[0], numNodes);  #endif
259    }
260    
261      buildIndexMap();  //
262    //
263    //
264    const IntVec& NodeData::getVarDataByName(const string& name) const
265    {
266        if (name == "Nodes_Id")
267            return nodeID;
268        else if (name == "Nodes_Tag")
269            return nodeTag;
270        else if (name == "Nodes_gDOF")
271            return nodeGDOF;
272        else if (name == "Nodes_gNI")
273            return nodeGNI;
274        else if (name == "Nodes_grDfI")
275            return nodeGRDFI;
276        else if (name == "Nodes_grNI")
277            return nodeGRNI;
278        else
279            throw "Invalid variable name";
280    }
281    
282      delete input;  //
283      return true;  //
284    //
285    StringVec NodeData::getVarNames() const
286    {
287        StringVec res;
288        if (numNodes > 0) {
289            res.push_back("Nodes_Id");
290            res.push_back("Nodes_Tag");
291            res.push_back("Nodes_gDOF");
292            res.push_back("Nodes_gNI");
293            res.push_back("Nodes_grDfI");
294            res.push_back("Nodes_grNI");
295        }
296        return res;
297  }  }
298    
299  //  //
300  //  //
301  //  //
302  bool Mesh::writeToSilo(DBfile* dbfile, const string& pathInSilo)  bool NodeData::writeToSilo(DBfile* dbfile)
303  {  {
304  #if HAVE_SILO  #if USE_SILO
305        if (numNodes == 0)
306            return true;
307    
308      int ret;      int ret;
309    
310      if (pathInSilo != "") {      if (siloPath != "") {
         siloPath = pathInSilo;  
311          ret = DBSetDir(dbfile, siloPath.c_str());          ret = DBSetDir(dbfile, siloPath.c_str());
312          if (ret != 0)          if (ret != 0)
313              return false;              return false;
     } else {  
         siloPath = "/";  
314      }      }
315          string siloMeshName = getFullSiloName();
316      // Write node-centered variable  
317      ret = DBPutUcdvar1(dbfile, "Nodes_Id", getFullSiloName().c_str(),      // Write node-centered variables
318        ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
319              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);              (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
320    
321        if (ret == 0)
322            ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
323                    (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
324                    DB_NODECENT, NULL);
325        if (ret == 0)
326            ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
327                    (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
328                    DB_NODECENT, NULL);
329        if (ret == 0)
330            ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
331                    (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
332                    DB_NODECENT, NULL);
333        if (ret == 0)
334            ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
335                    (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
336                    DB_NODECENT, NULL);
337        if (ret == 0)
338            ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
339                    (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
340                    DB_NODECENT, NULL);
341    
342      DBSetDir(dbfile, "/");      DBSetDir(dbfile, "/");
343      return (ret == 0);      return (ret == 0);
344    
345  #else // !HAVE_SILO  #else // !USE_SILO
346      return false;      return false;
347  #endif  #endif
348  }  }
349    
350  } // namespace EscriptReader  } // namespace escriptexport
351    
Legend:
Removed from v.2187  
changed lines
  Added in v.2853
  ViewVC Help
Powered by ViewVC 1.1.26