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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2548 - (hide annotations)
Mon Jul 20 06:20:06 2009 UTC (11 years, 4 months ago) by jfenwick
Original Path: trunk/tools/libescriptreader/src/escriptreader/MeshWithElements.cpp
File size: 9253 byte(s)
Updating copyright notices
1 caltinay 2183
2     /*******************************************************
3     *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 caltinay 2183 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     //
15     // MeshWithElements.cpp
16     //
17     #include <escriptreader/MeshWithElements.h>
18     #include <escriptreader/ElementData.h>
19     #include <netcdf.hh>
20     #if HAVE_SILO
21     #include <silo.h>
22     #endif
23    
24     using namespace std;
25    
26 caltinay 2187 namespace EscriptReader {
27    
28 caltinay 2183 //
29     //
30     //
31 caltinay 2194 MeshWithElements::MeshWithElements() : Mesh(),
32     cells(NULL),
33     faces(NULL),
34     contacts(NULL),
35     points(NULL)
36 caltinay 2183 {
37     name = "Elements";
38     }
39    
40     //
41     //
42     //
43     MeshWithElements::MeshWithElements(const MeshWithElements& m) :
44     Mesh(m)
45     {
46     nodeTag = m.nodeTag;
47     nodeGDOF = m.nodeGDOF;
48     nodeGNI = m.nodeGNI;
49     nodeGRDFI = m.nodeGRDFI;
50     nodeGRNI = m.nodeGRNI;
51     cells = new ElementData(*m.cells);
52     faces = new ElementData(*m.faces);
53     contacts = new ElementData(*m.contacts);
54     points = new ElementData(*m.points);
55     }
56    
57     //
58     //
59     //
60     MeshWithElements::~MeshWithElements()
61     {
62     delete cells;
63     delete faces;
64     delete contacts;
65     delete points;
66     }
67    
68     //
69     // Returns a vector of strings containing Silo mesh names that have been
70     // written
71     //
72     StringVec MeshWithElements::getMeshNames() const
73     {
74     StringVec res;
75     res.push_back(name);
76     StringVec tmpVec;
77     tmpVec = cells->getMeshNames();
78     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
79     tmpVec = faces->getMeshNames();
80     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
81     tmpVec = contacts->getMeshNames();
82     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
83     tmpVec = points->getMeshNames();
84     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
85     return res;
86     }
87    
88     Mesh* MeshWithElements::getMeshByName(const string name) const
89     {
90     Mesh* ret = NULL;
91     if (name == "Elements")
92     ret = cells->fullMesh;
93     else if (name == "ReducedElements")
94     ret = cells->reducedMesh;
95     else if (name == "FaceElements")
96     ret = faces->fullMesh;
97     else if (name == "ReducedFaceElements")
98     ret = faces->reducedMesh;
99     else if (name == "ContactElements")
100     ret = contacts->fullMesh;
101     else if (name == "ReducedContactElements")
102     ret = contacts->reducedMesh;
103     else if (name == "Points")
104     ret = points->fullMesh;
105    
106     return ret;
107     }
108    
109     //
110     // Returns a vector of strings containing Silo variable names that have
111     // been written
112     //
113     StringVec MeshWithElements::getVarNames() const
114     {
115     StringVec res;
116     res.push_back("Nodes_Id");
117     res.push_back("Nodes_Tag");
118     res.push_back("Nodes_gDOF");
119     res.push_back("Nodes_gNI");
120     res.push_back("Nodes_grDfI");
121     res.push_back("Nodes_grNI");
122    
123     StringVec tmpVec;
124     tmpVec = cells->getVarNames();
125     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
126     tmpVec = faces->getVarNames();
127     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
128     tmpVec = contacts->getVarNames();
129     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
130     tmpVec = points->getVarNames();
131     res.insert(res.end(), tmpVec.begin(), tmpVec.end());
132    
133     return res;
134     }
135    
136     const IntVec& MeshWithElements::getVarDataByName(const string name) const
137     {
138     if (name == "Nodes_Id")
139     return nodeID;
140     else if (name == "Nodes_Tag")
141     return nodeTag;
142     else if (name == "Nodes_gDOF")
143     return nodeGDOF;
144     else if (name == "Nodes_gNI")
145     return nodeGNI;
146     else if (name == "Nodes_grDfI")
147     return nodeGRDFI;
148     else if (name == "Nodes_grNI")
149     return nodeGRNI;
150     else if (name.compare(0, 9, "Elements_") == 0)
151     return cells->getVarDataByName(name);
152     else if (name.compare(0, 13, "FaceElements_") == 0)
153     return faces->getVarDataByName(name);
154     else if (name.compare(0, 16, "ContactElements_") == 0)
155     return contacts->getVarDataByName(name);
156     else if (name.compare(0, 7, "Points_") == 0)
157     return points->getVarDataByName(name);
158     else
159 caltinay 2201 throw "Invalid variable name";
160 caltinay 2183 }
161    
162     ElementData* MeshWithElements::getElementsByName(const string name) const
163     {
164     ElementData* ret = NULL;
165     if (name == "Elements" || name == "ReducedElements")
166     ret = cells;
167     else if (name == "FaceElements" || name == "ReducedFaceElements")
168     ret = faces;
169     else if (name == "ContactElements" || name == "ReducedContactElements")
170     ret = contacts;
171     else if (name == "Points")
172     ret = points;
173    
174     return ret;
175     }
176    
177     //
178     // Reads mesh and element data from NetCDF file with given name
179     //
180     bool MeshWithElements::readFromNc(const string& filename)
181     {
182     if (!Mesh::readFromNc(filename))
183     return false;
184    
185     NcError ncerr(NcError::silent_nonfatal);
186     NcFile* input;
187     NcVar* var;
188    
189     input = new NcFile(filename.c_str());
190     if (!input->is_valid()) {
191     cerr << "Could not open input file " << filename << "." << endl;
192     delete input;
193     return false;
194     }
195    
196     nodeTag.clear();
197     nodeTag.insert(nodeTag.end(), numNodes, 0);
198     var = input->get_var("Nodes_Tag");
199     var->get(&nodeTag[0], numNodes);
200    
201     nodeGDOF.clear();
202     nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
203     var = input->get_var("Nodes_gDOF");
204     var->get(&nodeGDOF[0], numNodes);
205    
206     nodeGNI.clear();
207     nodeGNI.insert(nodeGNI.end(), numNodes, 0);
208     var = input->get_var("Nodes_gNI");
209     var->get(&nodeGNI[0], numNodes);
210    
211     nodeGRDFI.clear();
212     nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
213     var = input->get_var("Nodes_grDfI");
214     var->get(&nodeGRDFI[0], numNodes);
215    
216     nodeGRNI.clear();
217     nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
218     var = input->get_var("Nodes_grNI");
219     var->get(&nodeGRNI[0], numNodes);
220    
221     // Read all element types
222     cells = new ElementData("Elements", this);
223     cells->readFromNc(input);
224     faces = new ElementData("FaceElements", this);
225     faces->readFromNc(input);
226     contacts = new ElementData("ContactElements", this);
227     contacts->readFromNc(input);
228     points = new ElementData("Points", this);
229     points->readFromNc(input);
230    
231     delete input;
232     return true;
233     }
234    
235     //
236     //
237     //
238     void MeshWithElements::handleGhostZones(int ownIndex)
239     {
240     cells->handleGhostZones(ownIndex);
241     faces->handleGhostZones(ownIndex);
242     contacts->handleGhostZones(ownIndex);
243     points->handleGhostZones(ownIndex);
244     #ifdef _DEBUG
245     cout << "block " << ownIndex << " has " << cells->getGhostCount()
246     << " (" << cells->getReducedGhostCount() << ") ghost zones," << endl;
247     cout << faces->getGhostCount() << " (" << faces->getReducedGhostCount()
248     << ") ghost faces," << endl;
249     cout << contacts->getGhostCount() << " ("
250     << contacts->getReducedGhostCount() << ") ghost contacts," << endl;
251     cout << points->getGhostCount() << " (" << points->getReducedGhostCount()
252     << ") ghost points." << endl;
253     #endif
254     }
255    
256     //
257     //
258     //
259     void MeshWithElements::removeGhostZones()
260     {
261     cells->removeGhostZones();
262     faces->removeGhostZones();
263     contacts->removeGhostZones();
264     points->removeGhostZones();
265     #ifdef _DEBUG
266     cout << "After removing ghost zones there are" << endl;
267     cout << " " << numNodes << " Nodes, ";
268     cout << cells->getCount() << " Elements, ";
269     cout << faces->getCount() << " Face elements, ";
270     cout << contacts->getCount() << " Contact elements, ";
271     cout << points->getCount() << " Points left." << endl;
272     #endif
273     }
274    
275     //
276     //
277     //
278     bool MeshWithElements::writeToSilo(DBfile* dbfile, const string& pathInSilo)
279     {
280     #if HAVE_SILO
281     int ret;
282    
283     // Write meshes and zone-centered variables
284     if (!cells->writeToSilo(dbfile, pathInSilo) ||
285     !faces->writeToSilo(dbfile, pathInSilo) ||
286     !contacts->writeToSilo(dbfile, pathInSilo) ||
287     !points->writeToSilo(dbfile, pathInSilo))
288     return false;
289    
290     if (!Mesh::writeToSilo(dbfile, pathInSilo))
291     return false;
292    
293     if (pathInSilo != "") {
294     ret = DBSetDir(dbfile, pathInSilo.c_str());
295     if (ret != 0)
296     return false;
297     }
298    
299     string siloMeshName = getFullSiloName();
300    
301     // Write node-centered variables
302     ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
303     (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
304     DB_NODECENT, NULL);
305     if (ret == 0)
306     ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
307     (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
308     DB_NODECENT, NULL);
309     if (ret == 0)
310     ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
311     (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
312     DB_NODECENT, NULL);
313     if (ret == 0)
314     ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
315     (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
316     DB_NODECENT, NULL);
317     if (ret == 0)
318     ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
319     (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
320     DB_NODECENT, NULL);
321    
322     DBSetDir(dbfile, "/");
323     return (ret == 0);
324    
325     #else // !HAVE_SILO
326     return false;
327     #endif
328     }
329    
330 caltinay 2187 } // namespace EscriptReader
331    

  ViewVC Help
Powered by ViewVC 1.1.26