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