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

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

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

trunk/tools/libescriptreader/src/escriptreader/MeshWithElements.cpp revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC trunk/dataexporter/src/FinleyMesh.cpp revision 2810 by caltinay, Mon Dec 7 04:13:49 2009 UTC
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
14  //  #include <escriptexport/FinleyMesh.h>
15  // MeshWithElements.cpp  #include <escriptexport/ElementData.h>
16  //  #include <escriptexport/NodeData.h>
17  #include <escriptreader/MeshWithElements.h>  
18  #include <escriptreader/ElementData.h>  #include <finley/CppAdapter/MeshAdapter.h>
19    extern "C" {
20    #include <finley/Mesh.h>
21    }
22    
23    #include <iostream>
24    
25    #if USE_NETCDF
26  #include <netcdf.hh>  #include <netcdf.hh>
27  #if HAVE_SILO  #endif
28    
29    #if USE_SILO
30  #include <silo.h>  #include <silo.h>
31  #endif  #endif
32    
33  using namespace std;  using namespace std;
34    
35  namespace EscriptReader {  namespace escriptexport {
36    
37  //  //
38  //  //
39  //  //
40  MeshWithElements::MeshWithElements() : Mesh(),  FinleyMesh::FinleyMesh() :
41      cells(NULL),      initialized(false)
     faces(NULL),  
     contacts(NULL),  
     points(NULL)  
42  {  {
     name = "Elements";  
43  }  }
44    
45  //  //
46  //  //
47  //  //
48  MeshWithElements::MeshWithElements(const MeshWithElements& m) :  FinleyMesh::FinleyMesh(const FinleyMesh& m)
     Mesh(m)  
49  {  {
50      nodeTag = m.nodeTag;      nodes = NodeData_ptr(new NodeData(*m.nodes));
51      nodeGDOF = m.nodeGDOF;      cells = ElementData_ptr(new ElementData(*m.cells));
52      nodeGNI = m.nodeGNI;      faces = ElementData_ptr(new ElementData(*m.faces));
53      nodeGRDFI = m.nodeGRDFI;      contacts = ElementData_ptr(new ElementData(*m.contacts));
54      nodeGRNI = m.nodeGRNI;      initialized = m.initialized;
     cells = new ElementData(*m.cells);  
     faces = new ElementData(*m.faces);  
     contacts = new ElementData(*m.contacts);  
     points = new ElementData(*m.points);  
55  }  }
56    
57  //  //
58  //  //
59  //  //
60  MeshWithElements::~MeshWithElements()  FinleyMesh::~FinleyMesh()
61  {  {
62      delete cells;      cleanup();
63      delete faces;  }
64      delete contacts;  
65      delete points;  //
66    //
67    //
68    void FinleyMesh::cleanup()
69    {
70        nodes.reset();
71        cells.reset();
72        faces.reset();
73        contacts.reset();
74        initialized = false;
75    }
76    
77    //
78    //
79    //
80    bool FinleyMesh::initFromEscript(escript::const_Domain_ptr escriptDomain)
81    {
82        cleanup();
83    
84        finleyMesh = dynamic_cast<const finley::MeshAdapter*>(escriptDomain.get())->getFinley_Mesh();
85        if (!finleyMesh) {
86            return false;
87        }
88    
89        nodes = NodeData_ptr(new NodeData("Elements"));
90        cells = ElementData_ptr(new ElementData("Elements", nodes));
91        faces = ElementData_ptr(new ElementData("FaceElements", nodes));
92        contacts = ElementData_ptr(new ElementData("ContactElements", nodes));
93    
94        if (nodes->initFromFinley(finleyMesh->Nodes) &&
95                cells->initFromFinley(finleyMesh->Elements) &&
96                faces->initFromFinley(finleyMesh->FaceElements) &&
97                contacts->initFromFinley(finleyMesh->ContactElements)) {
98            initialized = true;
99        }
100    
101        return initialized;
102    }
103    
104    //
105    // Reads mesh and element data from NetCDF file with given name
106    //
107    bool FinleyMesh::initFromNetCDF(const string& filename)
108    {
109        cleanup();
110        
111    #if USE_NETCDF
112        NcError ncerr(NcError::silent_nonfatal);
113        NcFile* input;
114    
115        input = new NcFile(filename.c_str());
116        if (!input->is_valid()) {
117            cerr << "Could not open input file " << filename << "." << endl;
118            delete input;
119            return false;
120        }
121    
122        nodes = NodeData_ptr(new NodeData("Elements"));
123        if (!nodes->readFromNc(input))
124            return false;
125    
126        // Read all element types
127        cells = ElementData_ptr(new ElementData("Elements", nodes));
128        cells->readFromNc(input);
129        faces = ElementData_ptr(new ElementData("FaceElements", nodes));
130        faces->readFromNc(input);
131        contacts = ElementData_ptr(new ElementData("ContactElements", nodes));
132        contacts->readFromNc(input);
133    
134        delete input;
135        initialized = true;
136    #endif
137    
138        return initialized;
139    }
140    
141    //
142    //
143    //
144    NodeData_ptr FinleyMesh::getMeshForFinleyFS(int functionSpace) const
145    {
146        NodeData_ptr result;
147    
148        if (!initialized)
149            return result;
150    
151        ElementData_ptr elements = getElementsForFinleyFS(functionSpace);
152        if (elements != NULL) {
153            if (functionSpace == FINLEY_NODES)
154                result = elements->getNodeMesh();
155            else {
156                if (elements->getReducedNumElements() > 0)
157                    result = elements->getReducedNodeMesh();
158                else
159                    result = elements->getNodeMesh();
160            }
161        }
162    
163        return result;
164    }
165    
166    //
167    //
168    //
169    ElementData_ptr FinleyMesh::getElementsForFinleyFS(int functionSpace) const
170    {
171        ElementData_ptr result;
172    
173        if (!initialized) {
174            return result;
175        }
176    
177        switch (functionSpace) {
178            case FINLEY_NODES:
179            case FINLEY_REDUCED_NODES:
180            case FINLEY_REDUCED_ELEMENTS:
181            case FINLEY_ELEMENTS:
182                result = cells;
183                break;
184    
185            case FINLEY_REDUCED_FACE_ELEMENTS:
186            case FINLEY_FACE_ELEMENTS:
187                result = faces;
188                break;
189    
190            case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
191            case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
192            case FINLEY_CONTACT_ELEMENTS_1:
193            case FINLEY_CONTACT_ELEMENTS_2:
194                result = contacts;
195                break;
196    
197            default:
198                cerr << "Unsupported function space type " << functionSpace
199                    << "!" << endl;
200        }
201        return result;
202  }  }
203    
204  //  //
205  // Returns a vector of strings containing Silo mesh names that have been  // Returns a vector of strings containing Silo mesh names that have been
206  // written  // written
207  //  //
208  StringVec MeshWithElements::getMeshNames() const  StringVec FinleyMesh::getMeshNames() const
209  {  {
210      StringVec res;      StringVec res;
211      res.push_back(name);      if (initialized) {
212      StringVec tmpVec;          StringVec tmpVec;
213      tmpVec = cells->getMeshNames();          tmpVec = cells->getMeshNames();
214      res.insert(res.end(), tmpVec.begin(), tmpVec.end());          res.insert(res.end(), tmpVec.begin(), tmpVec.end());
215      tmpVec = faces->getMeshNames();          tmpVec = faces->getMeshNames();
216      res.insert(res.end(), tmpVec.begin(), tmpVec.end());          res.insert(res.end(), tmpVec.begin(), tmpVec.end());
217      tmpVec = contacts->getMeshNames();          tmpVec = contacts->getMeshNames();
218      res.insert(res.end(), tmpVec.begin(), tmpVec.end());          res.insert(res.end(), tmpVec.begin(), tmpVec.end());
219      tmpVec = points->getMeshNames();      }
     res.insert(res.end(), tmpVec.begin(), tmpVec.end());  
220      return res;      return res;
221  }  }
222    
223  Mesh* MeshWithElements::getMeshByName(const string name) const  //
224    //
225    //
226    NodeData_ptr FinleyMesh::getMeshByName(const string& name) const
227  {  {
228      Mesh* ret = NULL;      NodeData_ptr ret;
229        if (!initialized) {
230            return ret;
231        }
232      if (name == "Elements")      if (name == "Elements")
233          ret = cells->fullMesh;          ret = cells->getNodeMesh();
234      else if (name == "ReducedElements")      else if (name == "ReducedElements")
235          ret = cells->reducedMesh;          ret = cells->getReducedNodeMesh();
236      else if (name == "FaceElements")      else if (name == "FaceElements")
237          ret = faces->fullMesh;          ret = faces->getNodeMesh();
238      else if (name == "ReducedFaceElements")      else if (name == "ReducedFaceElements")
239          ret = faces->reducedMesh;          ret = faces->getReducedNodeMesh();
240      else if (name == "ContactElements")      else if (name == "ContactElements")
241          ret = contacts->fullMesh;          ret = contacts->getNodeMesh();
242      else if (name == "ReducedContactElements")      else if (name == "ReducedContactElements")
243          ret = contacts->reducedMesh;          ret = contacts->getReducedNodeMesh();
     else if (name == "Points")  
         ret = points->fullMesh;  
244    
245      return ret;      return ret;
246  }  }
# Line 110  Mesh* MeshWithElements::getMeshByName(co Line 249  Mesh* MeshWithElements::getMeshByName(co
249  // Returns a vector of strings containing Silo variable names that have  // Returns a vector of strings containing Silo variable names that have
250  // been written  // been written
251  //  //
252  StringVec MeshWithElements::getVarNames() const  StringVec FinleyMesh::getVarNames() const
253  {  {
254      StringVec res;      StringVec res;
255      res.push_back("Nodes_Id");  
256      res.push_back("Nodes_Tag");      if (initialized) {
257      res.push_back("Nodes_gDOF");          res = nodes->getVarNames();
258      res.push_back("Nodes_gNI");          StringVec tmpVec = cells->getVarNames();
259      res.push_back("Nodes_grDfI");          res.insert(res.end(), tmpVec.begin(), tmpVec.end());
260      res.push_back("Nodes_grNI");          tmpVec = faces->getVarNames();
261            res.insert(res.end(), tmpVec.begin(), tmpVec.end());
262      StringVec tmpVec;          tmpVec = contacts->getVarNames();
263      tmpVec = cells->getVarNames();          res.insert(res.end(), tmpVec.begin(), tmpVec.end());
264      res.insert(res.end(), tmpVec.begin(), tmpVec.end());      }
     tmpVec = faces->getVarNames();  
     res.insert(res.end(), tmpVec.begin(), tmpVec.end());  
     tmpVec = contacts->getVarNames();  
     res.insert(res.end(), tmpVec.begin(), tmpVec.end());  
     tmpVec = points->getVarNames();  
     res.insert(res.end(), tmpVec.begin(), tmpVec.end());  
265    
266      return res;      return res;
267  }  }
268    
269  const IntVec& MeshWithElements::getVarDataByName(const string name) const  //
270    //
271    //
272    const IntVec& FinleyMesh::getVarDataByName(const string& name) const
273  {  {
274      if (name == "Nodes_Id")      if (!initialized) {
275          return nodeID;          throw "Mesh not initialized";
276      else if (name == "Nodes_Tag")      }
277          return nodeTag;      
278      else if (name == "Nodes_gDOF")      if (name.compare(0, 6, "Nodes_") == 0)
279          return nodeGDOF;          return nodes->getVarDataByName(name);
     else if (name == "Nodes_gNI")  
         return nodeGNI;  
     else if (name == "Nodes_grDfI")  
         return nodeGRDFI;  
     else if (name == "Nodes_grNI")  
         return nodeGRNI;  
280      else if (name.compare(0, 9, "Elements_") == 0)      else if (name.compare(0, 9, "Elements_") == 0)
281          return cells->getVarDataByName(name);          return cells->getVarDataByName(name);
282      else if (name.compare(0, 13, "FaceElements_") == 0)      else if (name.compare(0, 13, "FaceElements_") == 0)
283          return faces->getVarDataByName(name);          return faces->getVarDataByName(name);
284      else if (name.compare(0, 16, "ContactElements_") == 0)      else if (name.compare(0, 16, "ContactElements_") == 0)
285          return contacts->getVarDataByName(name);          return contacts->getVarDataByName(name);
     else if (name.compare(0, 7, "Points_") == 0)  
         return points->getVarDataByName(name);  
286      else      else
287          throw "Invalid variable name";          throw "Invalid variable name";
288  }  }
289    
290  ElementData* MeshWithElements::getElementsByName(const string name) const  //
291    //
292    //
293    ElementData_ptr FinleyMesh::getElementsByName(const string& name) const
294  {  {
295      ElementData* ret = NULL;      ElementData_ptr ret;
296      if (name == "Elements" || name == "ReducedElements")      if (name == "Elements" || name == "ReducedElements")
297          ret = cells;          ret = cells;
298      else if (name == "FaceElements" || name == "ReducedFaceElements")      else if (name == "FaceElements" || name == "ReducedFaceElements")
299          ret = faces;          ret = faces;
300      else if (name == "ContactElements" || name == "ReducedContactElements")      else if (name == "ContactElements" || name == "ReducedContactElements")
301          ret = contacts;          ret = contacts;
     else if (name == "Points")  
         ret = points;  
302    
303      return ret;      return ret;
304  }  }
305    
306  //  //
 // Reads mesh and element data from NetCDF file with given name  
 //  
 bool MeshWithElements::readFromNc(const string& filename)  
 {  
     if (!Mesh::readFromNc(filename))  
         return false;  
   
     NcError ncerr(NcError::silent_nonfatal);  
     NcFile* input;  
     NcVar* var;  
   
     input = new NcFile(filename.c_str());  
     if (!input->is_valid()) {  
         cerr << "Could not open input file " << filename << "." << endl;  
         delete input;  
         return false;  
     }  
   
     nodeTag.clear();  
     nodeTag.insert(nodeTag.end(), numNodes, 0);  
     var = input->get_var("Nodes_Tag");  
     var->get(&nodeTag[0], numNodes);  
   
     nodeGDOF.clear();  
     nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);  
     var = input->get_var("Nodes_gDOF");  
     var->get(&nodeGDOF[0], numNodes);  
   
     nodeGNI.clear();  
     nodeGNI.insert(nodeGNI.end(), numNodes, 0);  
     var = input->get_var("Nodes_gNI");  
     var->get(&nodeGNI[0], numNodes);  
   
     nodeGRDFI.clear();  
     nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);  
     var = input->get_var("Nodes_grDfI");  
     var->get(&nodeGRDFI[0], numNodes);  
   
     nodeGRNI.clear();  
     nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);  
     var = input->get_var("Nodes_grNI");  
     var->get(&nodeGRNI[0], numNodes);  
   
     // Read all element types  
     cells = new ElementData("Elements", this);  
     cells->readFromNc(input);  
     faces = new ElementData("FaceElements", this);  
     faces->readFromNc(input);  
     contacts = new ElementData("ContactElements", this);  
     contacts->readFromNc(input);  
     points = new ElementData("Points", this);  
     points->readFromNc(input);  
   
     delete input;  
     return true;  
 }  
   
 //  
307  //  //
308  //  //
309  void MeshWithElements::handleGhostZones(int ownIndex)  void FinleyMesh::reorderGhostZones(int ownIndex)
310  {  {
311      cells->handleGhostZones(ownIndex);      if (initialized) {
312      faces->handleGhostZones(ownIndex);          cells->reorderGhostZones(ownIndex);
313      contacts->handleGhostZones(ownIndex);          faces->reorderGhostZones(ownIndex);
314      points->handleGhostZones(ownIndex);          contacts->reorderGhostZones(ownIndex);
315  #ifdef _DEBUG  #ifdef _DEBUG
316      cout << "block " << ownIndex << " has " << cells->getGhostCount()          cout << "block " << ownIndex << " has " << cells->getGhostCount()
317           << " (" << cells->getReducedGhostCount() << ") ghost zones," << endl;               << " (reduced=" << cells->getReducedGhostCount()
318      cout << faces->getGhostCount() << " (" << faces->getReducedGhostCount()               << ") ghost zones," << endl;
319           << ") ghost faces," << endl;          cout << "\t" << faces->getGhostCount() << " ("
320      cout << contacts->getGhostCount() << " ("               << faces->getReducedGhostCount() << ") ghost faces," << endl;
321           << contacts->getReducedGhostCount() << ") ghost contacts," << endl;          cout << "\t" << contacts->getGhostCount() << " ("
322      cout << points->getGhostCount() << " (" << points->getReducedGhostCount()               << contacts->getReducedGhostCount() << ") ghost contacts." << endl;
          << ") ghost points." << endl;  
323  #endif  #endif
324        }
325  }  }
326    
327  //  //
328  //  //
329  //  //
330  void MeshWithElements::removeGhostZones()  void FinleyMesh::removeGhostZones(int ownIndex)
331  {  {
332      cells->removeGhostZones();      if (initialized) {
333      faces->removeGhostZones();          cells->removeGhostZones(ownIndex);
334      contacts->removeGhostZones();          faces->removeGhostZones(ownIndex);
335      points->removeGhostZones();          contacts->removeGhostZones(ownIndex);
336  #ifdef _DEBUG  #ifdef _DEBUG
337      cout << "After removing ghost zones there are" << endl;          cout << "After removing ghost zones there are" << endl;
338      cout << "    " << numNodes << " Nodes, ";          cout << "    " << nodes->getNumNodes() << " Nodes, ";
339      cout << cells->getCount() << " Elements, ";          cout << cells->getCount() << " Elements, ";
340      cout << faces->getCount() << " Face elements, ";          cout << faces->getCount() << " Face elements, ";
341      cout << contacts->getCount() << " Contact elements, ";          cout << contacts->getCount() << " Contact elements left." << endl;
     cout << points->getCount() << " Points left." << endl;  
342  #endif  #endif
343        }
344  }  }
345    
346  //  //
347  //  //
348  //  //
349  bool MeshWithElements::writeToSilo(DBfile* dbfile, const string& pathInSilo)  bool FinleyMesh::writeToSilo(DBfile* dbfile, const string& pathInSilo)
350  {  {
351  #if HAVE_SILO  #if USE_SILO
352      int ret;      // Write nodes, elements and mesh variables
353        if (!initialized ||
354      // Write meshes and zone-centered variables              !cells->writeToSilo(dbfile, pathInSilo) ||
     if (!cells->writeToSilo(dbfile, pathInSilo) ||  
355              !faces->writeToSilo(dbfile, pathInSilo) ||              !faces->writeToSilo(dbfile, pathInSilo) ||
356              !contacts->writeToSilo(dbfile, pathInSilo) ||              !contacts->writeToSilo(dbfile, pathInSilo))
             !points->writeToSilo(dbfile, pathInSilo))  
357          return false;          return false;
358    
359      if (!Mesh::writeToSilo(dbfile, pathInSilo))      siloPath = pathInSilo;
360          return false;      return true;
       
     if (pathInSilo != "") {  
         ret = DBSetDir(dbfile, pathInSilo.c_str());  
         if (ret != 0)  
             return false;  
     }  
   
     string siloMeshName = getFullSiloName();  
   
     // Write node-centered variables  
     ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),  
                 (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,  
                 DB_NODECENT, NULL);  
     if (ret == 0)  
         ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),  
                 (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,  
                 DB_NODECENT, NULL);  
     if (ret == 0)  
         ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),  
                 (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,  
                 DB_NODECENT, NULL);  
     if (ret == 0)  
         ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),  
                 (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,  
                 DB_NODECENT, NULL);  
     if (ret == 0)  
         ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),  
                 (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,  
                 DB_NODECENT, NULL);  
   
     DBSetDir(dbfile, "/");  
     return (ret == 0);  
361    
362  #else // !HAVE_SILO  #else // !USE_SILO
363      return false;      return false;
364  #endif  #endif
365  }  }
366    
367  } // namespace EscriptReader  } // namespace escriptexport
368    

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

  ViewVC Help
Powered by ViewVC 1.1.26