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 |
#include <escriptexport/NodeData.h> |
15 |
|
16 |
extern "C" { |
17 |
#include <finley/Mesh.h> |
18 |
#include <finley/NodeFile.h> |
19 |
} |
20 |
|
21 |
#if USE_NETCDF |
22 |
#include <netcdf.hh> |
23 |
#endif |
24 |
|
25 |
#if USE_SILO |
26 |
#include <silo.h> |
27 |
#endif |
28 |
|
29 |
using namespace std; |
30 |
|
31 |
namespace escriptexport { |
32 |
|
33 |
// |
34 |
// Constructor with name |
35 |
// |
36 |
NodeData::NodeData(const string& meshName) : |
37 |
numDims(0), numNodes(0), name(meshName) |
38 |
{ |
39 |
} |
40 |
|
41 |
// |
42 |
// |
43 |
// |
44 |
NodeData::NodeData(NodeData_ptr fullNodes, IntVec& requiredNodes, |
45 |
const string& meshName) : |
46 |
name(meshName) |
47 |
{ |
48 |
numDims = fullNodes->numDims; |
49 |
nodeDist = fullNodes->nodeDist; |
50 |
|
51 |
// first: find the unique set of required nodes and their IDs while |
52 |
// updating the contents of requiredNodes at the same time |
53 |
// requiredNodes contains node indices (not IDs!) |
54 |
IntVec::iterator it; |
55 |
IndexMap indexMap; // maps old index to new index |
56 |
size_t newIndex = 0; |
57 |
|
58 |
for (it = requiredNodes.begin(); it != requiredNodes.end(); it++) { |
59 |
IndexMap::iterator res = indexMap.find(*it); |
60 |
if (res == indexMap.end()) { |
61 |
nodeID.push_back(fullNodes->nodeID[*it]); |
62 |
nodeTag.push_back(fullNodes->nodeTag[*it]); |
63 |
nodeGDOF.push_back(fullNodes->nodeGDOF[*it]); |
64 |
nodeGNI.push_back(fullNodes->nodeGNI[*it]); |
65 |
nodeGRDFI.push_back(fullNodes->nodeGRDFI[*it]); |
66 |
nodeGRNI.push_back(fullNodes->nodeGRNI[*it]); |
67 |
indexMap[*it] = newIndex; |
68 |
*it = newIndex++; |
69 |
} else { |
70 |
*it = res->second; |
71 |
} |
72 |
} |
73 |
|
74 |
// second: now that we know how many nodes we need use the map to fill |
75 |
// the coordinates |
76 |
numNodes = newIndex; |
77 |
for (int dim=0; dim<numDims; dim++) { |
78 |
const float* origC = fullNodes->coords[dim]; |
79 |
float* c = new float[numNodes]; |
80 |
coords.push_back(c); |
81 |
IndexMap::const_iterator mIt; |
82 |
for (mIt = indexMap.begin(); mIt != indexMap.end(); mIt++) { |
83 |
c[mIt->second] = origC[mIt->first]; |
84 |
} |
85 |
} |
86 |
} |
87 |
|
88 |
// |
89 |
// Copy constructor |
90 |
// |
91 |
NodeData::NodeData(const NodeData& m) |
92 |
{ |
93 |
numDims = m.numDims; |
94 |
numNodes = m.numNodes; |
95 |
nodeID = m.nodeID; |
96 |
nodeTag = m.nodeTag; |
97 |
nodeGDOF = m.nodeGDOF; |
98 |
nodeGNI = m.nodeGNI; |
99 |
nodeGRDFI = m.nodeGRDFI; |
100 |
nodeGRNI = m.nodeGRNI; |
101 |
nodeDist = m.nodeDist; |
102 |
name = m.name; |
103 |
for (int i=0; i<numDims; i++) { |
104 |
float* c = new float[numNodes]; |
105 |
copy(m.coords[i], m.coords[i]+numNodes, c); |
106 |
coords.push_back(c); |
107 |
} |
108 |
} |
109 |
|
110 |
// |
111 |
// |
112 |
// |
113 |
NodeData::~NodeData() |
114 |
{ |
115 |
CoordArray::iterator it; |
116 |
for (it = coords.begin(); it != coords.end(); it++) |
117 |
delete[] *it; |
118 |
} |
119 |
|
120 |
// |
121 |
// |
122 |
// |
123 |
bool NodeData::initFromFinley(const Finley_NodeFile* finleyFile) |
124 |
{ |
125 |
numDims = finleyFile->numDim; |
126 |
numNodes = finleyFile->numNodes; |
127 |
|
128 |
int mpisize = finleyFile->MPIInfo->size; |
129 |
int* iPtr = finleyFile->nodesDistribution->first_component; |
130 |
nodeDist.clear(); |
131 |
nodeDist.insert(nodeDist.end(), mpisize+1, 0); |
132 |
copy(iPtr, iPtr+mpisize+1, nodeDist.begin()); |
133 |
|
134 |
CoordArray::iterator it; |
135 |
for (it = coords.begin(); it != coords.end(); it++) |
136 |
delete[] *it; |
137 |
coords.clear(); |
138 |
nodeID.clear(); |
139 |
nodeTag.clear(); |
140 |
nodeGDOF.clear(); |
141 |
nodeGNI.clear(); |
142 |
nodeGRDFI.clear(); |
143 |
nodeGRNI.clear(); |
144 |
|
145 |
if (numNodes > 0) { |
146 |
for (int i=0; i<numDims; i++) { |
147 |
double* srcPtr = finleyFile->Coordinates + i; |
148 |
float* c = new float[numNodes]; |
149 |
coords.push_back(c); |
150 |
for (int j=0; j<numNodes; j++, srcPtr+=numDims) { |
151 |
*c++ = (float) *srcPtr; |
152 |
} |
153 |
} |
154 |
|
155 |
iPtr = finleyFile->Id; |
156 |
nodeID.insert(nodeID.end(), numNodes, 0); |
157 |
copy(iPtr, iPtr+numNodes, nodeID.begin()); |
158 |
|
159 |
iPtr = finleyFile->Tag; |
160 |
nodeTag.insert(nodeTag.end(), numNodes, 0); |
161 |
copy(iPtr, iPtr+numNodes, nodeTag.begin()); |
162 |
|
163 |
iPtr = finleyFile->globalDegreesOfFreedom; |
164 |
nodeGDOF.insert(nodeGDOF.end(), numNodes, 0); |
165 |
copy(iPtr, iPtr+numNodes, nodeGDOF.begin()); |
166 |
|
167 |
iPtr = finleyFile->globalNodesIndex; |
168 |
nodeGNI.insert(nodeGNI.end(), numNodes, 0); |
169 |
copy(iPtr, iPtr+numNodes, nodeGNI.begin()); |
170 |
|
171 |
iPtr = finleyFile->globalReducedDOFIndex; |
172 |
nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0); |
173 |
copy(iPtr, iPtr+numNodes, nodeGRDFI.begin()); |
174 |
|
175 |
iPtr = finleyFile->globalReducedNodesIndex; |
176 |
nodeGRNI.insert(nodeGRNI.end(), numNodes, 0); |
177 |
copy(iPtr, iPtr+numNodes, nodeGRNI.begin()); |
178 |
|
179 |
} |
180 |
return true; |
181 |
} |
182 |
|
183 |
// |
184 |
// |
185 |
// |
186 |
bool NodeData::readFromNc(NcFile* ncFile) |
187 |
{ |
188 |
#if USE_NETCDF |
189 |
NcAtt* att; |
190 |
NcVar* var; |
191 |
|
192 |
att = ncFile->get_att("numDim"); |
193 |
numDims = att->as_int(0); |
194 |
|
195 |
att = ncFile->get_att("numNodes"); |
196 |
numNodes = att->as_int(0); |
197 |
|
198 |
att = ncFile->get_att("mpi_size"); |
199 |
int mpisize = att->as_int(0); |
200 |
|
201 |
nodeDist.clear(); |
202 |
nodeDist.insert(nodeDist.end(), mpisize+1, 0); |
203 |
var = ncFile->get_var("Nodes_NodeDistribution"); |
204 |
var->get(&nodeDist[0], mpisize+1); |
205 |
|
206 |
CoordArray::iterator it; |
207 |
for (it = coords.begin(); it != coords.end(); it++) |
208 |
delete[] *it; |
209 |
coords.clear(); |
210 |
nodeID.clear(); |
211 |
nodeTag.clear(); |
212 |
nodeGDOF.clear(); |
213 |
nodeGNI.clear(); |
214 |
nodeGRDFI.clear(); |
215 |
nodeGRNI.clear(); |
216 |
|
217 |
// Only attempt to read further if there are any nodes. |
218 |
// Having no nodes is not an error. |
219 |
if (numNodes > 0) { |
220 |
var = ncFile->get_var("Nodes_Coordinates"); |
221 |
for (int i=0; i<numDims; i++) { |
222 |
float* c = new float[numNodes]; |
223 |
var->set_cur(0, i); |
224 |
var->get(c, numNodes, 1); |
225 |
coords.push_back(c); |
226 |
} |
227 |
|
228 |
nodeID.insert(nodeID.end(), numNodes, 0); |
229 |
var = ncFile->get_var("Nodes_Id"); |
230 |
var->get(&nodeID[0], numNodes); |
231 |
|
232 |
nodeTag.insert(nodeTag.end(), numNodes, 0); |
233 |
var = ncFile->get_var("Nodes_Tag"); |
234 |
var->get(&nodeTag[0], numNodes); |
235 |
|
236 |
nodeGDOF.insert(nodeGDOF.end(), numNodes, 0); |
237 |
var = ncFile->get_var("Nodes_gDOF"); |
238 |
var->get(&nodeGDOF[0], numNodes); |
239 |
|
240 |
nodeGNI.insert(nodeGNI.end(), numNodes, 0); |
241 |
var = ncFile->get_var("Nodes_gNI"); |
242 |
var->get(&nodeGNI[0], numNodes); |
243 |
|
244 |
nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0); |
245 |
var = ncFile->get_var("Nodes_grDfI"); |
246 |
var->get(&nodeGRDFI[0], numNodes); |
247 |
|
248 |
nodeGRNI.insert(nodeGRNI.end(), numNodes, 0); |
249 |
var = ncFile->get_var("Nodes_grNI"); |
250 |
var->get(&nodeGRNI[0], numNodes); |
251 |
} |
252 |
|
253 |
return true; |
254 |
#else // !USE_NETCDF |
255 |
return false; |
256 |
#endif |
257 |
} |
258 |
|
259 |
// |
260 |
// |
261 |
// |
262 |
const IntVec& NodeData::getVarDataByName(const string& name) const |
263 |
{ |
264 |
if (name == "Nodes_Id") |
265 |
return nodeID; |
266 |
else if (name == "Nodes_Tag") |
267 |
return nodeTag; |
268 |
else if (name == "Nodes_gDOF") |
269 |
return nodeGDOF; |
270 |
else if (name == "Nodes_gNI") |
271 |
return nodeGNI; |
272 |
else if (name == "Nodes_grDfI") |
273 |
return nodeGRDFI; |
274 |
else if (name == "Nodes_grNI") |
275 |
return nodeGRNI; |
276 |
else |
277 |
throw "Invalid variable name"; |
278 |
} |
279 |
|
280 |
// |
281 |
// |
282 |
// |
283 |
StringVec NodeData::getVarNames() const |
284 |
{ |
285 |
StringVec res; |
286 |
if (numNodes > 0) { |
287 |
res.push_back("Nodes_Id"); |
288 |
res.push_back("Nodes_Tag"); |
289 |
res.push_back("Nodes_gDOF"); |
290 |
res.push_back("Nodes_gNI"); |
291 |
res.push_back("Nodes_grDfI"); |
292 |
res.push_back("Nodes_grNI"); |
293 |
} |
294 |
return res; |
295 |
} |
296 |
|
297 |
// |
298 |
// |
299 |
// |
300 |
bool NodeData::writeToSilo(DBfile* dbfile) |
301 |
{ |
302 |
#if USE_SILO |
303 |
if (numNodes == 0) |
304 |
return true; |
305 |
|
306 |
int ret; |
307 |
|
308 |
if (siloPath != "") { |
309 |
ret = DBSetDir(dbfile, siloPath.c_str()); |
310 |
if (ret != 0) |
311 |
return false; |
312 |
} |
313 |
string siloMeshName = getFullSiloName(); |
314 |
|
315 |
// Write node-centered variables |
316 |
ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(), |
317 |
(float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL); |
318 |
|
319 |
if (ret == 0) |
320 |
ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(), |
321 |
(float*)&nodeTag[0], numNodes, NULL, 0, DB_INT, |
322 |
DB_NODECENT, NULL); |
323 |
if (ret == 0) |
324 |
ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(), |
325 |
(float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT, |
326 |
DB_NODECENT, NULL); |
327 |
if (ret == 0) |
328 |
ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(), |
329 |
(float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT, |
330 |
DB_NODECENT, NULL); |
331 |
if (ret == 0) |
332 |
ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(), |
333 |
(float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT, |
334 |
DB_NODECENT, NULL); |
335 |
if (ret == 0) |
336 |
ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(), |
337 |
(float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT, |
338 |
DB_NODECENT, NULL); |
339 |
|
340 |
DBSetDir(dbfile, "/"); |
341 |
return (ret == 0); |
342 |
|
343 |
#else // !USE_SILO |
344 |
return false; |
345 |
#endif |
346 |
} |
347 |
|
348 |
} // namespace escriptexport |
349 |
|