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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2834 - (show annotations)
Thu Jan 7 06:06:56 2010 UTC (9 years, 8 months ago) by caltinay
File size: 10757 byte(s)
Reduced elements are now handled separate from the main elements within the
data exporter. This simplifies usage considerably.

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 CoordArray::iterator it;
129 for (it = coords.begin(); it != coords.end(); it++)
130 delete[] *it;
131 coords.clear();
132
133 if (numNodes > 0) {
134 for (int i=0; i<numDims; i++) {
135 double* srcPtr = finleyFile->Coordinates + i;
136 float* c = new float[numNodes];
137 coords.push_back(c);
138 for (int j=0; j<numNodes; j++, srcPtr+=numDims) {
139 *c++ = (float) *srcPtr;
140 }
141 }
142
143 int* iPtr;
144
145 iPtr = finleyFile->Id;
146 nodeID.clear();
147 nodeID.insert(nodeID.end(), numNodes, 0);
148 copy(iPtr, iPtr+numNodes, nodeID.begin());
149
150 iPtr = finleyFile->Tag;
151 nodeTag.clear();
152 nodeTag.insert(nodeTag.end(), numNodes, 0);
153 copy(iPtr, iPtr+numNodes, nodeTag.begin());
154
155 iPtr = finleyFile->globalDegreesOfFreedom;
156 nodeGDOF.clear();
157 nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
158 copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
159
160 iPtr = finleyFile->globalNodesIndex;
161 nodeGNI.clear();
162 nodeGNI.insert(nodeGNI.end(), numNodes, 0);
163 copy(iPtr, iPtr+numNodes, nodeGNI.begin());
164
165 iPtr = finleyFile->globalReducedDOFIndex;
166 nodeGRDFI.clear();
167 nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
168 copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
169
170 iPtr = finleyFile->globalReducedNodesIndex;
171 nodeGRNI.clear();
172 nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
173 copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
174
175 int mpisize = finleyFile->MPIInfo->size;
176 iPtr = finleyFile->nodesDistribution->first_component;
177 nodeDist.clear();
178 nodeDist.insert(nodeDist.end(), mpisize+1, 0);
179 copy(iPtr, iPtr+mpisize+1, nodeDist.begin());
180 }
181 return true;
182 }
183
184 //
185 //
186 //
187 bool NodeData::readFromNc(NcFile* ncFile)
188 {
189 #if USE_NETCDF
190 NcAtt* att;
191 NcVar* var;
192
193 att = ncFile->get_att("numDim");
194 numDims = att->as_int(0);
195
196 att = ncFile->get_att("numNodes");
197 numNodes = att->as_int(0);
198
199 att = ncFile->get_att("mpi_size");
200 int mpisize = att->as_int(0);
201
202 CoordArray::iterator it;
203 for (it = coords.begin(); it != coords.end(); it++)
204 delete[] *it;
205 coords.clear();
206 var = ncFile->get_var("Nodes_Coordinates");
207 for (int i=0; i<numDims; i++) {
208 float* c = new float[numNodes];
209 var->set_cur(0, i);
210 var->get(c, numNodes, 1);
211 coords.push_back(c);
212 }
213
214 nodeID.clear();
215 nodeID.insert(nodeID.end(), numNodes, 0);
216 var = ncFile->get_var("Nodes_Id");
217 var->get(&nodeID[0], numNodes);
218
219 nodeTag.clear();
220 nodeTag.insert(nodeTag.end(), numNodes, 0);
221 var = ncFile->get_var("Nodes_Tag");
222 var->get(&nodeTag[0], numNodes);
223
224 nodeGDOF.clear();
225 nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
226 var = ncFile->get_var("Nodes_gDOF");
227 var->get(&nodeGDOF[0], numNodes);
228
229 nodeGNI.clear();
230 nodeGNI.insert(nodeGNI.end(), numNodes, 0);
231 var = ncFile->get_var("Nodes_gNI");
232 var->get(&nodeGNI[0], numNodes);
233
234 nodeGRDFI.clear();
235 nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
236 var = ncFile->get_var("Nodes_grDfI");
237 var->get(&nodeGRDFI[0], numNodes);
238
239 nodeGRNI.clear();
240 nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
241 var = ncFile->get_var("Nodes_grNI");
242 var->get(&nodeGRNI[0], numNodes);
243
244 nodeDist.clear();
245 nodeDist.insert(nodeDist.end(), mpisize+1, 0);
246 var = ncFile->get_var("Nodes_NodeDistribution");
247 var->get(&nodeDist[0], mpisize+1);
248
249 return true;
250 #else // !USE_NETCDF
251 return false;
252 #endif
253 }
254
255 //
256 //
257 //
258 const IntVec& NodeData::getVarDataByName(const string& name) const
259 {
260 if (name == "Nodes_Id")
261 return nodeID;
262 else if (name == "Nodes_Tag")
263 return nodeTag;
264 else if (name == "Nodes_gDOF")
265 return nodeGDOF;
266 else if (name == "Nodes_gNI")
267 return nodeGNI;
268 else if (name == "Nodes_grDfI")
269 return nodeGRDFI;
270 else if (name == "Nodes_grNI")
271 return nodeGRNI;
272 else
273 throw "Invalid variable name";
274 }
275
276 //
277 //
278 //
279 StringVec NodeData::getVarNames() const
280 {
281 StringVec res;
282 if (numNodes > 0) {
283 res.push_back("Nodes_Id");
284 res.push_back("Nodes_Tag");
285 res.push_back("Nodes_gDOF");
286 res.push_back("Nodes_gNI");
287 res.push_back("Nodes_grDfI");
288 res.push_back("Nodes_grNI");
289 }
290 return res;
291 }
292
293 //
294 //
295 //
296 void NodeData::removeGhostNodes(int ownIndex)
297 {
298 if (nodeDist.empty() || ownIndex > nodeDist.size()-1)
299 return;
300
301 int firstId = nodeDist[ownIndex];
302 int lastId = nodeDist[ownIndex+1];
303
304 // no ghost nodes
305 if (lastId-firstId == numNodes)
306 return;
307
308 // we have at most lastId-firstId nodes, it could be less however if
309 // nodes were culled already
310 numNodes = lastId-firstId;
311
312 CoordArray newCoords;
313 CoordArray::iterator it;
314 for (int i=0; i<numDims; i++) {
315 float* c = new float[numNodes];
316 newCoords.push_back(c);
317 }
318
319 IntVec newNodeID, newNodeTag;
320 IntVec newNodeGDOF, newNodeGNI, newNodeGRDFI, newNodeGRNI;
321
322 IndexMap nodeID2idx = getIndexMap();
323 int destIdx = 0;
324 for (int i=firstId; i<lastId; i++) {
325 IndexMap::iterator it = nodeID2idx.find(i);
326 if (it == nodeID2idx.end()) {
327 continue;
328 }
329 int idx = it->second;
330 for (int dim=0; dim<numDims; dim++) {
331 newCoords[dim][destIdx] = coords[dim][idx];
332 }
333 destIdx++;
334 newNodeID.push_back(i);
335 newNodeTag.push_back(nodeTag[idx]);
336 newNodeGDOF.push_back(nodeGDOF[idx]);
337 newNodeGNI.push_back(nodeGNI[idx]);
338 newNodeGRDFI.push_back(nodeGRDFI[idx]);
339 newNodeGRNI.push_back(nodeGRNI[idx]);
340 }
341
342 numNodes = destIdx;
343
344 for (it = coords.begin(); it != coords.end(); it++)
345 delete[] *it;
346
347 coords = newCoords;
348 nodeID = newNodeID;
349 nodeTag = newNodeTag;
350 nodeGDOF = newNodeGDOF;
351 nodeGNI = newNodeGNI;
352 nodeGRDFI = newNodeGRDFI;
353 nodeGRNI = newNodeGRNI;
354 }
355
356 //
357 //
358 //
359 bool NodeData::writeToSilo(DBfile* dbfile)
360 {
361 #if USE_SILO
362 if (numNodes == 0)
363 return true;
364
365 int ret;
366
367 if (siloPath != "") {
368 ret = DBSetDir(dbfile, siloPath.c_str());
369 if (ret != 0)
370 return false;
371 }
372 string siloMeshName = getFullSiloName();
373
374 // Write node-centered variables
375 ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
376 (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
377
378 if (ret == 0)
379 ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
380 (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
381 DB_NODECENT, NULL);
382 if (ret == 0)
383 ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
384 (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
385 DB_NODECENT, NULL);
386 if (ret == 0)
387 ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
388 (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
389 DB_NODECENT, NULL);
390 if (ret == 0)
391 ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
392 (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
393 DB_NODECENT, NULL);
394 if (ret == 0)
395 ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
396 (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
397 DB_NODECENT, NULL);
398
399 DBSetDir(dbfile, "/");
400 return (ret == 0);
401
402 #else // !USE_SILO
403 return false;
404 #endif
405 }
406
407 } // namespace escriptexport
408

  ViewVC Help
Powered by ViewVC 1.1.26