/[escript]/trunk/weipa/src/FinleyNodes.cpp
ViewVC logotype

Contents of /trunk/weipa/src/FinleyNodes.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 12078 byte(s)
Merging dudley and scons updates from branches

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 <weipa/FinleyNodes.h>
15
16 #ifndef VISIT_PLUGIN
17 extern "C" {
18 #include <dudley/Mesh.h>
19 #include <dudley/NodeFile.h>
20 #include <finley/Mesh.h>
21 #include <finley/NodeFile.h>
22 }
23 #endif
24
25 #if USE_NETCDF
26 #include <netcdfcpp.h>
27 #endif
28
29 #if USE_SILO
30 #include <silo.h>
31 #endif
32
33 using namespace std;
34
35 namespace weipa {
36
37 //
38 // Constructor with name
39 //
40 FinleyNodes::FinleyNodes(const string& meshName) :
41 numDims(0), numNodes(0), name(meshName)
42 {
43 }
44
45 //
46 //
47 //
48 FinleyNodes::FinleyNodes(FinleyNodes_ptr fullNodes, IntVec& requiredNodes,
49 const string& meshName) :
50 name(meshName)
51 {
52 numDims = fullNodes->numDims;
53 nodeDist = fullNodes->nodeDist;
54
55 // first: find the unique set of required nodes and their IDs while
56 // updating the contents of requiredNodes at the same time
57 // requiredNodes contains node indices (not IDs!)
58 IntVec::iterator it;
59 IndexMap indexMap; // maps old index to new index
60 size_t newIndex = 0;
61
62 for (it = requiredNodes.begin(); it != requiredNodes.end(); it++) {
63 IndexMap::iterator res = indexMap.find(*it);
64 if (res == indexMap.end()) {
65 nodeID.push_back(fullNodes->nodeID[*it]);
66 nodeTag.push_back(fullNodes->nodeTag[*it]);
67 nodeGDOF.push_back(fullNodes->nodeGDOF[*it]);
68 nodeGNI.push_back(fullNodes->nodeGNI[*it]);
69 nodeGRDFI.push_back(fullNodes->nodeGRDFI[*it]);
70 nodeGRNI.push_back(fullNodes->nodeGRNI[*it]);
71 indexMap[*it] = newIndex;
72 *it = newIndex++;
73 } else {
74 *it = res->second;
75 }
76 }
77
78 // second: now that we know how many nodes we need use the map to fill
79 // the coordinates
80 numNodes = newIndex;
81 for (int dim=0; dim<numDims; dim++) {
82 const float* origC = fullNodes->coords[dim];
83 float* c = new float[numNodes];
84 coords.push_back(c);
85 IndexMap::const_iterator mIt;
86 for (mIt = indexMap.begin(); mIt != indexMap.end(); mIt++) {
87 c[mIt->second] = origC[mIt->first];
88 }
89 }
90 }
91
92 //
93 // Copy constructor
94 //
95 FinleyNodes::FinleyNodes(const FinleyNodes& m)
96 {
97 numDims = m.numDims;
98 numNodes = m.numNodes;
99 nodeID = m.nodeID;
100 nodeTag = m.nodeTag;
101 nodeGDOF = m.nodeGDOF;
102 nodeGNI = m.nodeGNI;
103 nodeGRDFI = m.nodeGRDFI;
104 nodeGRNI = m.nodeGRNI;
105 nodeDist = m.nodeDist;
106 name = m.name;
107 for (int i=0; i<numDims; i++) {
108 float* c = new float[numNodes];
109 copy(m.coords[i], m.coords[i]+numNodes, c);
110 coords.push_back(c);
111 }
112 }
113
114 //
115 //
116 //
117 FinleyNodes::~FinleyNodes()
118 {
119 CoordArray::iterator it;
120 for (it = coords.begin(); it != coords.end(); it++)
121 delete[] *it;
122 }
123
124 //
125 //
126 //
127 bool FinleyNodes::initFromDudley(const Dudley_NodeFile* dudleyFile)
128 {
129 #ifndef VISIT_PLUGIN
130 numDims = dudleyFile->numDim;
131 numNodes = dudleyFile->numNodes;
132
133 int mpisize = dudleyFile->MPIInfo->size;
134 int* iPtr = dudleyFile->nodesDistribution->first_component;
135 nodeDist.clear();
136 nodeDist.insert(nodeDist.end(), mpisize+1, 0);
137 copy(iPtr, iPtr+mpisize+1, nodeDist.begin());
138
139 CoordArray::iterator it;
140 for (it = coords.begin(); it != coords.end(); it++)
141 delete[] *it;
142 coords.clear();
143 nodeID.clear();
144 nodeTag.clear();
145 nodeGDOF.clear();
146 nodeGNI.clear();
147 nodeGRDFI.clear();
148 nodeGRNI.clear();
149
150 if (numNodes > 0) {
151 for (int i=0; i<numDims; i++) {
152 double* srcPtr = dudleyFile->Coordinates + i;
153 float* c = new float[numNodes];
154 coords.push_back(c);
155 for (int j=0; j<numNodes; j++, srcPtr+=numDims) {
156 *c++ = (float) *srcPtr;
157 }
158 }
159
160 iPtr = dudleyFile->Id;
161 nodeID.insert(nodeID.end(), numNodes, 0);
162 copy(iPtr, iPtr+numNodes, nodeID.begin());
163
164 iPtr = dudleyFile->Tag;
165 nodeTag.insert(nodeTag.end(), numNodes, 0);
166 copy(iPtr, iPtr+numNodes, nodeTag.begin());
167
168 iPtr = dudleyFile->globalDegreesOfFreedom;
169 nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
170 copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
171
172 iPtr = dudleyFile->globalNodesIndex;
173 nodeGNI.insert(nodeGNI.end(), numNodes, 0);
174 copy(iPtr, iPtr+numNodes, nodeGNI.begin());
175
176 iPtr = dudleyFile->globalReducedDOFIndex;
177 nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
178 copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
179
180 iPtr = dudleyFile->globalReducedNodesIndex;
181 nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
182 copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
183
184 }
185 return true;
186 #else // VISIT_PLUGIN
187 return false;
188 #endif
189 }
190
191 //
192 //
193 //
194 bool FinleyNodes::initFromFinley(const Finley_NodeFile* finleyFile)
195 {
196 #ifndef VISIT_PLUGIN
197 numDims = finleyFile->numDim;
198 numNodes = finleyFile->numNodes;
199
200 int mpisize = finleyFile->MPIInfo->size;
201 int* iPtr = finleyFile->nodesDistribution->first_component;
202 nodeDist.clear();
203 nodeDist.insert(nodeDist.end(), mpisize+1, 0);
204 copy(iPtr, iPtr+mpisize+1, nodeDist.begin());
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 if (numNodes > 0) {
218 for (int i=0; i<numDims; i++) {
219 double* srcPtr = finleyFile->Coordinates + i;
220 float* c = new float[numNodes];
221 coords.push_back(c);
222 for (int j=0; j<numNodes; j++, srcPtr+=numDims) {
223 *c++ = (float) *srcPtr;
224 }
225 }
226
227 iPtr = finleyFile->Id;
228 nodeID.insert(nodeID.end(), numNodes, 0);
229 copy(iPtr, iPtr+numNodes, nodeID.begin());
230
231 iPtr = finleyFile->Tag;
232 nodeTag.insert(nodeTag.end(), numNodes, 0);
233 copy(iPtr, iPtr+numNodes, nodeTag.begin());
234
235 iPtr = finleyFile->globalDegreesOfFreedom;
236 nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
237 copy(iPtr, iPtr+numNodes, nodeGDOF.begin());
238
239 iPtr = finleyFile->globalNodesIndex;
240 nodeGNI.insert(nodeGNI.end(), numNodes, 0);
241 copy(iPtr, iPtr+numNodes, nodeGNI.begin());
242
243 iPtr = finleyFile->globalReducedDOFIndex;
244 nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
245 copy(iPtr, iPtr+numNodes, nodeGRDFI.begin());
246
247 iPtr = finleyFile->globalReducedNodesIndex;
248 nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
249 copy(iPtr, iPtr+numNodes, nodeGRNI.begin());
250
251 }
252 return true;
253 #else // VISIT_PLUGIN
254 return false;
255 #endif
256 }
257
258 //
259 //
260 //
261 bool FinleyNodes::readFromNc(NcFile* ncFile)
262 {
263 #if USE_NETCDF
264 NcAtt* att;
265 NcVar* var;
266
267 att = ncFile->get_att("numDim");
268 numDims = att->as_int(0);
269
270 att = ncFile->get_att("numNodes");
271 numNodes = att->as_int(0);
272
273 att = ncFile->get_att("mpi_size");
274 int mpisize = att->as_int(0);
275
276 nodeDist.clear();
277 nodeDist.insert(nodeDist.end(), mpisize+1, 0);
278 var = ncFile->get_var("Nodes_NodeDistribution");
279 var->get(&nodeDist[0], mpisize+1);
280
281 CoordArray::iterator it;
282 for (it = coords.begin(); it != coords.end(); it++)
283 delete[] *it;
284 coords.clear();
285 nodeID.clear();
286 nodeTag.clear();
287 nodeGDOF.clear();
288 nodeGNI.clear();
289 nodeGRDFI.clear();
290 nodeGRNI.clear();
291
292 // Only attempt to read further if there are any nodes.
293 // Having no nodes is not an error.
294 if (numNodes > 0) {
295 var = ncFile->get_var("Nodes_Coordinates");
296 for (int i=0; i<numDims; i++) {
297 float* c = new float[numNodes];
298 var->set_cur(0, i);
299 var->get(c, numNodes, 1);
300 coords.push_back(c);
301 }
302
303 nodeID.insert(nodeID.end(), numNodes, 0);
304 var = ncFile->get_var("Nodes_Id");
305 var->get(&nodeID[0], numNodes);
306
307 nodeTag.insert(nodeTag.end(), numNodes, 0);
308 var = ncFile->get_var("Nodes_Tag");
309 var->get(&nodeTag[0], numNodes);
310
311 nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
312 var = ncFile->get_var("Nodes_gDOF");
313 var->get(&nodeGDOF[0], numNodes);
314
315 nodeGNI.insert(nodeGNI.end(), numNodes, 0);
316 var = ncFile->get_var("Nodes_gNI");
317 var->get(&nodeGNI[0], numNodes);
318
319 nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
320 var = ncFile->get_var("Nodes_grDfI");
321 var->get(&nodeGRDFI[0], numNodes);
322
323 nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
324 var = ncFile->get_var("Nodes_grNI");
325 var->get(&nodeGRNI[0], numNodes);
326 }
327
328 return true;
329 #else // !USE_NETCDF
330 return false;
331 #endif
332 }
333
334 //
335 //
336 //
337 const IntVec& FinleyNodes::getVarDataByName(const string& name) const
338 {
339 if (name == "Nodes_Id")
340 return nodeID;
341 else if (name == "Nodes_Tag")
342 return nodeTag;
343 else if (name == "Nodes_gDOF")
344 return nodeGDOF;
345 else if (name == "Nodes_gNI")
346 return nodeGNI;
347 else if (name == "Nodes_grDfI")
348 return nodeGRDFI;
349 else if (name == "Nodes_grNI")
350 return nodeGRNI;
351 else
352 throw "Invalid variable name";
353 }
354
355 //
356 //
357 //
358 StringVec FinleyNodes::getVarNames() const
359 {
360 StringVec res;
361 res.push_back("Nodes_Id");
362 res.push_back("Nodes_Tag");
363 res.push_back("Nodes_gDOF");
364 res.push_back("Nodes_gNI");
365 res.push_back("Nodes_grDfI");
366 res.push_back("Nodes_grNI");
367 return res;
368 }
369
370 //
371 //
372 //
373 int FinleyNodes::getGlobalNumNodes() const
374 {
375 int ret=0;
376 if (!nodeDist.empty())
377 ret = nodeDist[nodeDist.size()-1];
378 return ret;
379 }
380
381 //
382 //
383 //
384 void FinleyNodes::writeCoordinatesVTK(ostream& os, int ownIndex)
385 {
386 if (numNodes > 0) {
387 int firstId = nodeDist[ownIndex];
388 int lastId = nodeDist[ownIndex+1];
389 for (size_t i=0; i<numNodes; i++) {
390 if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
391 os << coords[0][i] << " " << coords[1][i] << " ";
392 if (numDims == 3)
393 os << coords[2][i];
394 else
395 os << 0.;
396 os << endl;
397 }
398 }
399 }
400 }
401
402 //
403 //
404 //
405 bool FinleyNodes::writeToSilo(DBfile* dbfile)
406 {
407 #if USE_SILO
408 if (numNodes == 0)
409 return true;
410
411 int ret;
412
413 if (siloPath != "") {
414 ret = DBSetDir(dbfile, siloPath.c_str());
415 if (ret != 0)
416 return false;
417 }
418 string siloMeshName = getFullSiloName();
419
420 // Write node-centered variables
421 ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
422 (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
423
424 if (ret == 0)
425 ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
426 (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
427 DB_NODECENT, NULL);
428 if (ret == 0)
429 ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
430 (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
431 DB_NODECENT, NULL);
432 if (ret == 0)
433 ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
434 (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
435 DB_NODECENT, NULL);
436 if (ret == 0)
437 ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
438 (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
439 DB_NODECENT, NULL);
440 if (ret == 0)
441 ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
442 (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
443 DB_NODECENT, NULL);
444
445 DBSetDir(dbfile, "/");
446 return (ret == 0);
447
448 #else // !USE_SILO
449 return false;
450 #endif
451 }
452
453 } // namespace weipa
454

  ViewVC Help
Powered by ViewVC 1.1.26