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

Contents of /trunk/dataexporter/src/NodeData.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2810 - (show annotations)
Mon Dec 7 04:13:49 2009 UTC (9 years, 9 months ago) by caltinay
File size: 10382 byte(s)
Reincarnation of the escriptreader as a more flexible escriptexport library:
 - can be initialised with instances of escript::Data and Finley_Mesh
 - is now MPI aware at the EscriptDataset level including Silo writer
 - now uses boost shared pointers
The lib is currently only used by the escriptconvert tool but this is going to
change...

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 // Copy constructor
89 //
90 NodeData::NodeData(const NodeData& m)
91 {
92 numDims = m.numDims;
93 numNodes = m.numNodes;
94 nodeID = m.nodeID;
95 nodeTag = m.nodeTag;
96 nodeGDOF = m.nodeGDOF;
97 nodeGNI = m.nodeGNI;
98 nodeGRDFI = m.nodeGRDFI;
99 nodeGRNI = m.nodeGRNI;
100 nodeDist = m.nodeDist;
101 name = m.name;
102 for (int i=0; i<numDims; i++) {
103 float* c = new float[numNodes];
104 copy(m.coords[i], m.coords[i]+numNodes, c);
105 coords.push_back(c);
106 }
107 }
108
109 //
110 //
111 //
112 NodeData::~NodeData()
113 {
114 CoordArray::iterator it;
115 for (it = coords.begin(); it != coords.end(); it++)
116 delete[] *it;
117 }
118
119 //
120 //
121 //
122 bool NodeData::initFromFinley(const Finley_NodeFile* finleyFile)
123 {
124 numDims = finleyFile->numDim;
125 numNodes = finleyFile->numNodes;
126
127 CoordArray::iterator it;
128 for (it = coords.begin(); it != coords.end(); it++)
129 delete[] *it;
130 coords.clear();
131
132 if (numNodes > 0) {
133 for (int i=0; i<numDims; i++) {
134 double* srcPtr = finleyFile->Coordinates + i;
135 float* c = new float[numNodes];
136 coords.push_back(c);
137 for (int j=0; j<numNodes; j++, srcPtr+=numDims) {
138 *c++ = (float) *srcPtr;
139 }
140 }
141
142 int* iPtr;
143
144 iPtr = finleyFile->Id;
145 nodeID.clear();
146 nodeID.insert(nodeID.end(), numNodes, 0);
147 copy(iPtr, iPtr+numNodes, nodeID.begin());
148
149 iPtr = finleyFile->Tag;
150 nodeTag.clear();
151 nodeTag.insert(nodeTag.end(), numNodes, 0);
152 copy(iPtr, iPtr+numNodes, nodeTag.begin());
153
154 iPtr = finleyFile->globalDegreesOfFreedom;
155 nodeGDOF.clear();
156 nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
157 copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
158
159 iPtr = finleyFile->globalNodesIndex;
160 nodeGNI.clear();
161 nodeGNI.insert(nodeGNI.end(), numNodes, 0);
162 copy(iPtr, iPtr+numNodes, nodeGNI.begin());
163
164 iPtr = finleyFile->globalReducedDOFIndex;
165 nodeGRDFI.clear();
166 nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
167 copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
168
169 iPtr = finleyFile->globalReducedNodesIndex;
170 nodeGRNI.clear();
171 nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
172 copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
173
174 int mpisize = finleyFile->MPIInfo->size;
175 iPtr = finleyFile->nodesDistribution->first_component;
176 nodeDist.clear();
177 nodeDist.insert(nodeDist.end(), mpisize+1, 0);
178 copy(iPtr, iPtr+mpisize+1, nodeDist.begin());
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 CoordArray::iterator it;
202 for (it = coords.begin(); it != coords.end(); it++)
203 delete[] *it;
204 coords.clear();
205 var = ncFile->get_var("Nodes_Coordinates");
206 for (int i=0; i<numDims; i++) {
207 float* c = new float[numNodes];
208 var->set_cur(0, i);
209 var->get(c, numNodes, 1);
210 coords.push_back(c);
211 }
212
213 nodeID.clear();
214 nodeID.insert(nodeID.end(), numNodes, 0);
215 var = ncFile->get_var("Nodes_Id");
216 var->get(&nodeID[0], numNodes);
217
218 nodeTag.clear();
219 nodeTag.insert(nodeTag.end(), numNodes, 0);
220 var = ncFile->get_var("Nodes_Tag");
221 var->get(&nodeTag[0], numNodes);
222
223 nodeGDOF.clear();
224 nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
225 var = ncFile->get_var("Nodes_gDOF");
226 var->get(&nodeGDOF[0], numNodes);
227
228 nodeGNI.clear();
229 nodeGNI.insert(nodeGNI.end(), numNodes, 0);
230 var = ncFile->get_var("Nodes_gNI");
231 var->get(&nodeGNI[0], numNodes);
232
233 nodeGRDFI.clear();
234 nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
235 var = ncFile->get_var("Nodes_grDfI");
236 var->get(&nodeGRDFI[0], numNodes);
237
238 nodeGRNI.clear();
239 nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
240 var = ncFile->get_var("Nodes_grNI");
241 var->get(&nodeGRNI[0], numNodes);
242
243 nodeDist.clear();
244 nodeDist.insert(nodeDist.end(), mpisize+1, 0);
245 var = ncFile->get_var("Nodes_NodeDistribution");
246 var->get(&nodeDist[0], mpisize+1);
247
248 return true;
249 #else // !USE_NETCDF
250 return false;
251 #endif
252 }
253
254 //
255 //
256 //
257 const IntVec& NodeData::getVarDataByName(const std::string& name) const
258 {
259 if (name == "Nodes_Id")
260 return nodeID;
261 else if (name == "Nodes_Tag")
262 return nodeTag;
263 else if (name == "Nodes_gDOF")
264 return nodeGDOF;
265 else if (name == "Nodes_gNI")
266 return nodeGNI;
267 else if (name == "Nodes_grDfI")
268 return nodeGRDFI;
269 else if (name == "Nodes_grNI")
270 return nodeGRNI;
271 else
272 throw "Invalid variable name";
273 }
274
275 //
276 //
277 //
278 StringVec NodeData::getVarNames() const
279 {
280 StringVec res;
281 if (numNodes > 0) {
282 res.push_back("Nodes_Id");
283 res.push_back("Nodes_Tag");
284 res.push_back("Nodes_gDOF");
285 res.push_back("Nodes_gNI");
286 res.push_back("Nodes_grDfI");
287 res.push_back("Nodes_grNI");
288 }
289 return res;
290 }
291
292 //
293 //
294 //
295 void NodeData::removeGhostNodes(int ownIndex)
296 {
297 if (nodeDist.empty() || ownIndex > nodeDist.size()-1)
298 return;
299
300 int firstId = nodeDist[ownIndex];
301 int lastId = nodeDist[ownIndex+1];
302
303 numNodes = lastId-firstId;
304
305 CoordArray newCoords;
306 CoordArray::iterator it;
307 for (int i=0; i<numDims; i++) {
308 float* c = new float[numNodes];
309 newCoords.push_back(c);
310 }
311
312 IntVec newNodeID, newNodeTag;
313 IntVec newNodeGDOF, newNodeGNI, newNodeGRDFI, newNodeGRNI;
314
315 IndexMap nodeID2idx = getIndexMap();
316
317 for (int i=firstId; i<lastId; i++) {
318 int idx = nodeID2idx[i];
319 for (int j=0; j<numDims; j++) {
320 newCoords[j][i-firstId] = coords[j][idx];
321 }
322 newNodeID.push_back(i);
323 newNodeTag.push_back(nodeTag[idx]);
324 newNodeGDOF.push_back(nodeGDOF[idx]);
325 newNodeGNI.push_back(nodeGNI[idx]);
326 newNodeGRDFI.push_back(nodeGRDFI[idx]);
327 newNodeGRNI.push_back(nodeGRNI[idx]);
328 }
329
330 for (it = coords.begin(); it != coords.end(); it++)
331 delete[] *it;
332
333 coords = newCoords;
334 nodeID = newNodeID;
335 nodeTag = newNodeTag;
336 nodeGDOF = newNodeGDOF;
337 nodeGNI = newNodeGNI;
338 nodeGRDFI = newNodeGRDFI;
339 nodeGRNI = newNodeGRNI;
340 }
341
342 //
343 //
344 //
345 bool NodeData::writeToSilo(DBfile* dbfile)
346 {
347 #if USE_SILO
348 if (numNodes == 0)
349 return true;
350
351 int ret;
352
353 if (siloPath != "") {
354 ret = DBSetDir(dbfile, siloPath.c_str());
355 if (ret != 0)
356 return false;
357 }
358 string siloMeshName = getFullSiloName();
359
360 // Write node-centered variables
361 ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
362 (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
363
364 if (ret == 0)
365 ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
366 (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
367 DB_NODECENT, NULL);
368 if (ret == 0)
369 ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
370 (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
371 DB_NODECENT, NULL);
372 if (ret == 0)
373 ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
374 (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
375 DB_NODECENT, NULL);
376 if (ret == 0)
377 ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
378 (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
379 DB_NODECENT, NULL);
380 if (ret == 0)
381 ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
382 (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
383 DB_NODECENT, NULL);
384
385 DBSetDir(dbfile, "/");
386 return (ret == 0);
387
388 #else // !USE_SILO
389 return false;
390 #endif
391 }
392
393 } // namespace escriptexport
394

  ViewVC Help
Powered by ViewVC 1.1.26