1 |
|
2 |
/******************************************************* |
3 |
* |
4 |
* Copyright (c) 2003-2009 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 |
throw "Invalid variable name"; |
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 |
|