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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3792 - (show annotations)
Wed Feb 1 06:16:25 2012 UTC (7 years, 6 months ago) by caltinay
File size: 15139 byte(s)
Merged ripley rectangular domain into trunk.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2012 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/RipleyNodes.h>
15
16 #ifndef VISIT_PLUGIN
17 #include <ripley/RipleyDomain.h>
18 #endif
19
20 #if USE_SILO
21 #include <silo.h>
22 #endif
23
24 using namespace std;
25
26 namespace weipa {
27
28 //
29 // Constructor with name
30 //
31 RipleyNodes::RipleyNodes(const string& meshName) :
32 numDims(0), numNodes(0), globalNumNodes(0), name(meshName)
33 {
34 }
35
36 //
37 //
38 //
39 RipleyNodes::RipleyNodes(RipleyNodes_ptr fullNodes, IntVec& requiredNodes,
40 const string& meshName) :
41 name(meshName)
42 {
43 numDims = fullNodes->numDims;
44 nodeDist = fullNodes->nodeDist;
45 globalNumNodes = fullNodes->globalNumNodes;
46
47 // first: find the unique set of required nodes and their IDs while
48 // updating the contents of requiredNodes at the same time
49 // requiredNodes contains node indices (not IDs!)
50 IntVec::iterator it;
51 IndexMap indexMap; // maps old index to new index
52 size_t newIndex = 0;
53
54 for (it = requiredNodes.begin(); it != requiredNodes.end(); it++) {
55 IndexMap::iterator res = indexMap.find(*it);
56 if (res == indexMap.end()) {
57 nodeID.push_back(fullNodes->nodeID[*it]);
58 nodeGNI.push_back(fullNodes->nodeGNI[*it]);
59 nodeTag.push_back(fullNodes->nodeTag[*it]);
60 indexMap[*it] = newIndex;
61 *it = newIndex++;
62 } else {
63 *it = res->second;
64 }
65 }
66
67 // second: now that we know how many nodes we need use the map to fill
68 // the coordinates
69 numNodes = newIndex;
70 for (int dim=0; dim<numDims; dim++) {
71 const float* origC = fullNodes->coords[dim];
72 float* c = new float[numNodes];
73 coords.push_back(c);
74 IndexMap::const_iterator mIt;
75 for (mIt = indexMap.begin(); mIt != indexMap.end(); mIt++) {
76 c[mIt->second] = origC[mIt->first];
77 }
78 }
79 }
80
81 //
82 // Copy constructor
83 //
84 RipleyNodes::RipleyNodes(const RipleyNodes& m)
85 {
86 numDims = m.numDims;
87 numNodes = m.numNodes;
88 globalNumNodes = m.globalNumNodes;
89 nodeID = m.nodeID;
90 nodeTag = m.nodeTag;
91 nodeDist = m.nodeDist;
92 name = m.name;
93 for (int i=0; i<numDims; i++) {
94 float* c = new float[numNodes];
95 copy(m.coords[i], m.coords[i]+numNodes, c);
96 coords.push_back(c);
97 }
98 }
99
100 //
101 //
102 //
103 RipleyNodes::~RipleyNodes()
104 {
105 CoordArray::iterator it;
106 for (it = coords.begin(); it != coords.end(); it++)
107 delete[] *it;
108 }
109
110 //
111 //
112 //
113 bool RipleyNodes::initFromRipley(const ripley::RipleyDomain* dom)
114 {
115 #ifndef VISIT_PLUGIN
116 CoordArray::iterator it;
117 for (it = coords.begin(); it != coords.end(); it++)
118 delete[] *it;
119 coords.clear();
120 nodeID.clear();
121 nodeTag.clear();
122
123 numDims = dom->getDim();
124 globalNumNodes = dom->getNumDataPointsGlobal();
125 pair<int,int> shape = dom->getDataShape(ripley::Nodes);
126 numNodes = shape.second;
127 nodeDist = dom->getNodeDistribution();
128 nodeGNI.assign(numNodes, 0);
129
130 if (numNodes > 0) {
131 for (int d=0; d<numDims; d++) {
132 float* c = new float[numNodes];
133 coords.push_back(c);
134 }
135 const IntVec NN = dom->getNumNodesPerDim();
136 const IntVec faces = dom->getNumFacesPerBoundary();
137 const IntVec NS = dom->getNumSubdivisionsPerDim();
138
139 if (numDims==2) {
140 pair<double,double> xx=dom->getFirstCoordAndSpacing(0);
141 pair<double,double> yy=dom->getFirstCoordAndSpacing(1);
142 const int left=(faces[0]==0 ? 1 : 0);
143 const int right=(faces[1]==0 ? NN[0]-2 : NN[0]-1);
144 const int bottom=(faces[2]==0 ? 1 : 0);
145 const int top=(faces[3]==0 ? NN[1]-2 : NN[1]-1);
146 const int ownN0=right-left+1;
147 const int ownN1=top-bottom+1;
148 const int rx=dom->getMPIRank()%NS[0];
149 const int ry=dom->getMPIRank()/NS[0];
150 for (int i1=0; i1<NN[1]; i1++) {
151 for (int i0=0; i0<NN[0]; i0++) {
152 coords[0][i0+NN[0]*i1] = (float)(xx.first+i0*xx.second);
153 coords[1][i0+NN[0]*i1] = (float)(yy.first+i1*yy.second);
154 }
155 }
156 for (int i1=-1; i1<2; i1++) {
157 for (int i0=-1; i0<2; i0++) {
158 // location of neighbour rank
159 const int nrx=rx+i0;
160 const int nry=ry+i1;
161 if (nrx>=0 && nry>=0 && nrx<NS[0] && nry<NS[1]) {
162 // index of first node on neighbour rank
163 const int first=nodeDist[nry*NS[0]+nrx];
164 if (i0==0 && i1==0) {
165 // own nodes
166 for (int y=bottom; y<=top; y++)
167 for (int x=left; x<=right; x++)
168 nodeGNI[x+y*NN[0]]=first+x-left+(y-bottom)*ownN0;
169 } else if (i0==0) {
170 // top or bottom
171 const int myFirst=(i1==-1 ? 0 : NN[0]*(NN[1]-1));
172 const int nFirst=(i1==-1 ? first+ownN0*(ownN1-1) : first);
173 for (int i=left; i<=right; i++)
174 nodeGNI[myFirst+i]=nFirst+i-left;
175 } else if (i1==0) {
176 // left or right
177 const int myFirst=(i0==-1 ? 0 : NN[0]-1);
178 const int nFirst=(i0==-1 ? first+ownN0-1 : first);
179 for (int i=bottom; i<=top; i++)
180 nodeGNI[myFirst+i*NN[0]]= nFirst+(i-bottom)*ownN0;
181 } else {
182 // corner
183 const int mIdx=(i0+1)/2*(NN[0]-1)+(i1+1)/2*(NN[0]*(NN[1]-1));
184 const int nIdx=(1-i0)/2*(ownN0-1)+(1-i1)/2*(ownN0*(ownN1-1));
185 nodeGNI[mIdx]=first+nIdx;
186 }
187 }
188 }
189 }
190
191 } else {
192 pair<double,double> xx=dom->getFirstCoordAndSpacing(0);
193 pair<double,double> yy=dom->getFirstCoordAndSpacing(1);
194 pair<double,double> zz=dom->getFirstCoordAndSpacing(2);
195 const int left=(faces[0]==0 ? 1 : 0);
196 const int right=(faces[1]==0 ? NN[0]-2 : NN[0]-1);
197 const int bottom=(faces[2]==0 ? 1 : 0);
198 const int top=(faces[3]==0 ? NN[1]-2 : NN[1]-1);
199 const int front=(faces[4]==0 ? 1 : 0);
200 const int back=(faces[5]==0 ? NN[2]-2 : NN[2]-1);
201 const int ownN0=right-left+1;
202 const int ownN1=top-bottom+1;
203 const int ownN2=back-front+1;
204 const int rx=dom->getMPIRank()%NS[0];
205 const int ry=dom->getMPIRank()%(NS[0]*NS[1])/NS[0];
206 const int rz=dom->getMPIRank()/(NS[0]*NS[1]);
207 for (int i2=0; i2<NN[2]; i2++) {
208 for (int i1=0; i1<NN[1]; i1++) {
209 for (int i0=0; i0<NN[0]; i0++) {
210 const int index = i0+NN[0]*i1+NN[0]*NN[1]*i2;
211 coords[0][index] = (float)(xx.first+i0*xx.second);
212 coords[1][index] = (float)(yy.first+i1*yy.second);
213 coords[2][index] = (float)(zz.first+i2*zz.second);
214 }
215 }
216 }
217 for (int i2=-1; i2<2; i2++) {
218 for (int i1=-1; i1<2; i1++) {
219 for (int i0=-1; i0<2; i0++) {
220 // location of neighbour rank
221 const int nrx=rx+i0;
222 const int nry=ry+i1;
223 const int nrz=rz+i2;
224 if (nrx>=0 && nry>=0 && nrz>=0
225 && nrx<NS[0] && nry<NS[1] && nrz<NS[2]) {
226 // index of first node on neighbour rank
227 const int first=nodeDist[nrz*NS[0]*NS[1]+nry*NS[0]+nrx];
228 if (i0==0 && i1==0 && i2==0) {
229 // own nodes
230 for (int z=front; z<=back; z++)
231 for (int y=bottom; y<=top; y++)
232 for (int x=left; x<=right; x++)
233 nodeGNI[x+y*NN[0]+z*NN[0]*NN[1]]=
234 first+x-left+(y-bottom)*ownN0
235 +(z-front)*ownN0*ownN1;
236 } else if (i0==0 && i1==0) {
237 // front or back plane
238 for (int y=bottom; y<=top; y++) {
239 const int myFirst=(i2==-1 ? y*NN[0] : NN[0]*NN[1]*(NN[2]-1)+y*NN[0]);
240 const int nFirst=(i2==-1 ?
241 first+ownN0*ownN1*(ownN2-1)+(y-bottom)*ownN0 : first+(y-bottom)*ownN0);
242 for (int x=left; x<=right; x++)
243 nodeGNI[myFirst+x]=nFirst+x-left;
244 }
245 } else if (i0==0 && i2==0) {
246 // top or bottom plane
247 for (int z=front; z<=back; z++) {
248 const int myFirst=(i1==-1 ? z*NN[0]*NN[1] : NN[0]*((z+1)*NN[1]-1));
249 const int nFirst=(i1==-1 ?
250 first+ownN0*((z-front+1)*ownN1-1) : first+(z-front)*ownN0*ownN1);
251 for (int x=left; x<=right; x++)
252 nodeGNI[myFirst+x]=nFirst+x-left;
253 }
254 } else if (i1==0 && i2==0) {
255 // left or right plane
256 for (int z=front; z<=back; z++) {
257 const int myFirst=(i0==-1 ? z*NN[0]*NN[1] : NN[0]*(1+z*NN[1])-1);
258 const int nFirst=(i0==-1 ?
259 first+ownN0*(1+(z-front)*ownN1)-1 : first+(z-front)*ownN0*ownN1);
260 for (int y=bottom; y<=top; y++)
261 nodeGNI[myFirst+y*NN[0]]=nFirst+(y-bottom)*ownN0;
262 }
263 } else if (i0==0) {
264 // edge in x direction
265 const int myFirst=(i1+1)/2*NN[0]*(NN[1]-1)
266 +(i2+1)/2*NN[0]*NN[1]*(NN[2]-1);
267 const int nFirst=first+(1-i1)/2*ownN0*(ownN1-1)
268 +(1-i2)/2*ownN0*ownN1*(ownN2-1);
269 for (int i=left; i<=right; i++)
270 nodeGNI[myFirst+i]=nFirst+i-left;
271 } else if (i1==0) {
272 // edge in y direction
273 const int myFirst=(i0+1)/2*(NN[0]-1)
274 +(i2+1)/2*NN[0]*NN[1]*(NN[2]-1);
275 const int nFirst=first+(1-i0)/2*(ownN0-1)
276 +(1-i2)/2*ownN0*ownN1*(ownN2-1);
277 for (int i=bottom; i<=top; i++)
278 nodeGNI[myFirst+i*NN[0]]=nFirst+(i-bottom)*ownN0;
279 } else if (i2==0) {
280 // edge in z direction
281 const int myFirst=(i0+1)/2*(NN[0]-1)
282 +(i1+1)/2*NN[0]*(NN[1]-1);
283 const int nFirst=first+(1-i0)/2*(ownN0-1)
284 +(1-i1)/2*ownN0*(ownN1-1);
285 for (int i=front; i<=back; i++)
286 nodeGNI[myFirst+i*NN[0]*NN[1]]=nFirst+(i-front)*ownN0*ownN1;
287 } else {
288 // corner
289 const int mIdx=(i0+1)/2*(NN[0]-1)
290 +(i1+1)/2*NN[0]*(NN[1]-1)
291 +(i2+1)/2*NN[0]*NN[1]*(NN[2]-1);
292 const int nIdx=(1-i0)/2*(ownN0-1)
293 +(1-i1)/2*ownN0*(ownN1-1)
294 +(1-i2)/2*ownN0*ownN1*(ownN2-1);
295 nodeGNI[mIdx]=first+nIdx;
296 }
297 }
298 }
299 }
300 }
301 }
302
303 const int* iPtr = dom->borrowSampleReferenceIDs(ripley::Nodes);
304 nodeID.assign(iPtr, iPtr+numNodes);
305
306 //iPtr = dom->borrowListOfTags(ripley::Nodes);
307 nodeTag.assign(iPtr, iPtr+numNodes);
308 }
309
310 return true;
311 #else // VISIT_PLUGIN
312 return false;
313 #endif
314 }
315
316 //
317 //
318 //
319 const IntVec& RipleyNodes::getVarDataByName(const string& name) const
320 {
321 if (name == "Nodes_Id")
322 return nodeID;
323 if (name == "Nodes_Tag")
324 return nodeTag;
325
326 throw "Invalid variable name";
327 }
328
329 //
330 //
331 //
332 StringVec RipleyNodes::getVarNames() const
333 {
334 StringVec res;
335 res.push_back("Nodes_Id");
336 res.push_back("Nodes_Tag");
337 return res;
338 }
339
340 //
341 //
342 //
343 void RipleyNodes::writeCoordinatesVTK(ostream& os, int ownIndex)
344 {
345 if (numNodes > 0) {
346 int firstId = nodeDist[ownIndex];
347 int lastId = nodeDist[ownIndex+1];
348 for (size_t i=0; i<numNodes; i++) {
349 if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
350 os << coords[0][i] << " " << coords[1][i] << " ";
351 if (numDims == 3)
352 os << coords[2][i];
353 else
354 os << 0.;
355 os << endl;
356 }
357 }
358 }
359 }
360
361 //
362 //
363 //
364 bool RipleyNodes::writeToSilo(DBfile* dbfile)
365 {
366 #if USE_SILO
367 if (numNodes == 0)
368 return true;
369
370 int ret;
371
372 if (siloPath != "") {
373 ret = DBSetDir(dbfile, siloPath.c_str());
374 if (ret != 0)
375 return false;
376 }
377 string siloMeshName = getFullSiloName();
378
379 // Write node-centered variables
380 ret = DBPutUcdvar1(dbfile, "Nodes_Id", siloMeshName.c_str(),
381 (float*)&nodeID[0], numNodes, NULL, 0, DB_INT, DB_NODECENT, NULL);
382
383 if (ret == 0)
384 ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
385 (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
386 DB_NODECENT, NULL);
387
388 DBSetDir(dbfile, "/");
389 return (ret == 0);
390
391 #else // !USE_SILO
392 return false;
393 #endif
394 }
395
396 } // namespace weipa
397

  ViewVC Help
Powered by ViewVC 1.1.26