/[escript]/trunk/tools/libescriptreader/src/escriptreader/MeshWithElements.cpp
ViewVC logotype

Contents of /trunk/tools/libescriptreader/src/escriptreader/MeshWithElements.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2194 - (show annotations)
Tue Jan 6 01:23:37 2009 UTC (10 years, 7 months ago) by caltinay
File size: 9247 byte(s)
escript2silo: Allow different meshes for different timesteps.
escriptreader: Fixed segfaults for some error conditions.


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * 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 namespace EscriptReader {
27
28 //
29 //
30 //
31 MeshWithElements::MeshWithElements() : Mesh(),
32 cells(NULL),
33 faces(NULL),
34 contacts(NULL),
35 points(NULL)
36 {
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 return *(IntVec*)(NULL);
160 }
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 } // namespace EscriptReader
331

  ViewVC Help
Powered by ViewVC 1.1.26