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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2886 - (show annotations)
Thu Jan 28 05:39:23 2010 UTC (9 years, 7 months ago) by caltinay
File size: 9935 byte(s)
New version of saveVTK within dataexporter. Not pythonised yet since it needs
more testing and release 3.1 will not ship dataexporter.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 res.push_back("Nodes_Id");
287 res.push_back("Nodes_Tag");
288 res.push_back("Nodes_gDOF");
289 res.push_back("Nodes_gNI");
290 res.push_back("Nodes_grDfI");
291 res.push_back("Nodes_grNI");
292 return res;
293 }
294
295 //
296 //
297 //
298 int NodeData::getGlobalNumNodes() const
299 {
300 int ret=0;
301 if (!nodeDist.empty())
302 ret = nodeDist[nodeDist.size()-1];
303 return ret;
304 }
305
306 //
307 //
308 //
309 void NodeData::writeCoordinatesVTK(ostream& os, int ownIndex)
310 {
311 if (numNodes > 0) {
312 int firstId = nodeDist[ownIndex];
313 int lastId = nodeDist[ownIndex+1];
314 for (size_t i=0; i<numNodes; i++) {
315 if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
316 os << coords[0][i] << " " << coords[1][i] << " ";
317 if (numDims == 3)
318 os << coords[2][i];
319 else
320 os << 0.;
321 os << endl;
322 }
323 }
324 }
325 }
326
327 //
328 //
329 //
330 bool NodeData::writeToSilo(DBfile* dbfile)
331 {
332 #if USE_SILO
333 if (numNodes == 0)
334 return true;
335
336 int ret;
337
338 if (siloPath != "") {
339 ret = DBSetDir(dbfile, siloPath.c_str());
340 if (ret != 0)
341 return false;
342 }
343 string siloMeshName = getFullSiloName();
344
345 // Write node-centered variables
346 ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
347 (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
348
349 if (ret == 0)
350 ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
351 (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
352 DB_NODECENT, NULL);
353 if (ret == 0)
354 ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
355 (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
356 DB_NODECENT, NULL);
357 if (ret == 0)
358 ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
359 (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
360 DB_NODECENT, NULL);
361 if (ret == 0)
362 ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
363 (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
364 DB_NODECENT, NULL);
365 if (ret == 0)
366 ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
367 (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
368 DB_NODECENT, NULL);
369
370 DBSetDir(dbfile, "/");
371 return (ret == 0);
372
373 #else // !USE_SILO
374 return false;
375 #endif
376 }
377
378 } // namespace escriptexport
379

  ViewVC Help
Powered by ViewVC 1.1.26