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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3792 - (show annotations)
Wed Feb 1 06:16:25 2012 UTC (7 years, 1 month ago) by caltinay
File size: 19881 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/RipleyElements.h>
15 #include <weipa/NodeData.h>
16
17 #ifndef VISIT_PLUGIN
18 #include <ripley/RipleyDomain.h>
19 #endif
20
21 #include <iostream>
22
23 #if USE_SILO
24 #include <silo.h>
25 #endif
26
27
28 using namespace std;
29
30 namespace weipa {
31
32 //
33 // Constructor
34 //
35 RipleyElements::RipleyElements(const string& elementName, RipleyNodes_ptr nodeData)
36 : originalMesh(nodeData), name(elementName), numElements(0),
37 numGhostElements(0), nodesPerElement(0),
38 type(ZONETYPE_UNKNOWN)
39 {
40 nodeMesh.reset(new RipleyNodes(name));
41 }
42
43 //
44 // Copy constructor
45 //
46 RipleyElements::RipleyElements(const RipleyElements& e)
47 {
48 name = e.name;
49 numElements = e.numElements;
50 numGhostElements = e.numGhostElements;
51 type = e.type;
52 nodesPerElement = e.nodesPerElement;
53 originalMesh = e.originalMesh;
54 if (e.nodeMesh)
55 nodeMesh.reset(new RipleyNodes(*e.nodeMesh));
56 else
57 nodeMesh.reset(new RipleyNodes(name));
58
59 NperDim = e.NperDim;
60 nodes = e.nodes;
61 ID = e.ID;
62 tag = e.tag;
63 owner = e.owner;
64 }
65
66 //
67 //
68 //
69 bool RipleyElements::initFromRipley(const ripley::RipleyDomain* dom, int fsType)
70 {
71 #ifndef VISIT_PLUGIN
72 const pair<int,int> shape = dom->getDataShape(fsType);
73 const IntVec faces = dom->getNumFacesPerBoundary();
74 const IntVec NS = dom->getNumSubdivisionsPerDim();
75
76 numElements = shape.second;
77 if (fsType==ripley::Elements)
78 NperDim = dom->getNumElementsPerDim();
79 else
80 NperDim = faces;
81
82 if (numElements > 0) {
83 nodesPerElement = shape.first;
84 switch (nodesPerElement) {
85 case 2:
86 type = ZONETYPE_BEAM;
87 break;
88 case 4:
89 type = ZONETYPE_QUAD;
90 break;
91 case 8:
92 type = ZONETYPE_HEX;
93 break;
94 }
95 owner.assign(numElements, dom->getMPIRank());
96
97 const int* iPtr = dom->borrowSampleReferenceIDs(fsType);
98 ID.assign(iPtr, iPtr+numElements);
99
100 //iPtr = dom->borrowListOfTags(fsType);
101 tag.assign(iPtr, iPtr+numElements);
102
103 const IntVec NN = dom->getNumNodesPerDim();
104 nodes.clear();
105 if (dom->getDim() == 2) {
106 if (fsType==ripley::Elements) {
107 if (faces[0]==0) {
108 owner[0]=(faces[2]==0 ? dom->getMPIRank()-NS[0]-1
109 : dom->getMPIRank()-1);
110 for (int i=1; i<NperDim[1]; i++)
111 owner[i*NperDim[0]]=dom->getMPIRank()-1;
112 }
113 if (faces[2]==0) {
114 const int first=(faces[0]==0 ? 1 : 0);
115 for (int i=first; i<NperDim[0]; i++)
116 owner[i]=dom->getMPIRank()-NS[0];
117 }
118 int id=0;
119 for (int i=0; i<numElements; i++) {
120 nodes.push_back(id);
121 nodes.push_back(id+1);
122 nodes.push_back(id+1+NN[0]);
123 nodes.push_back(id+NN[0]);
124 id++;
125 if ((i+1)%NperDim[0]==0)
126 id++;
127 }
128 } else if (fsType==ripley::FaceElements) {
129 if (faces[0]==0) {
130 if (faces[2]>0)
131 owner[faces[1]]=dom->getMPIRank()-1;
132 if (faces[3]>0)
133 owner[faces[1]+faces[2]]=dom->getMPIRank()-1;
134 }
135 if (faces[2]==0) {
136 if (faces[0]>0)
137 owner[0]=dom->getMPIRank()-NS[0];
138 if (faces[1]>0)
139 owner[faces[0]]=dom->getMPIRank()-NS[0];
140 }
141
142 int id=0;
143 for (int i=0; i<NperDim[0]; i++) {
144 nodes.push_back(id);
145 nodes.push_back(id+NN[0]);
146 id+=NN[0];
147 }
148 id=NN[0]-1;
149 for (int i=0; i<NperDim[1]; i++) {
150 nodes.push_back(id);
151 nodes.push_back(id+NN[0]);
152 id+=NN[0];
153 }
154 id=0;
155 for (int i=0; i<NperDim[2]; i++) {
156 nodes.push_back(id);
157 nodes.push_back(id+1);
158 id++;
159 }
160 id=NN[0]*(NN[1]-1);
161 for (int i=0; i<NperDim[3]; i++) {
162 nodes.push_back(id);
163 nodes.push_back(id+1);
164 id++;
165 }
166 }
167 } else {
168 if (fsType==ripley::Elements) {
169 // ownership is not entirely correct but that is not critical.
170 // fix when there is time.
171 if (faces[1]==0) {
172 for (int k2=0; k2<NperDim[2]; k2++) {
173 for (int k1=0; k1<NperDim[1]; k1++) {
174 const int e=k2*NperDim[0]*NperDim[1]+(k1+1)*NperDim[0]-1;
175 owner[e]=dom->getMPIRank()+1;
176 }
177 }
178 }
179 if (faces[3]==0) {
180 for (int k2=0; k2<NperDim[2]; k2++) {
181 for (int k0=0; k0<NperDim[0]; k0++) {
182 const int e=(k2+1)*NperDim[0]*NperDim[1]-NperDim[0]+k0;
183 owner[e]=dom->getMPIRank()+NS[0];
184 }
185 }
186 }
187 if (faces[5]==0) {
188 for (int k1=0; k1<NperDim[1]; k1++) {
189 for (int k0=0; k0<NperDim[0]; k0++) {
190 const int e=k1*NperDim[0]+k0+NperDim[0]*NperDim[1]*(NperDim[2]-1);
191 owner[e]=dom->getMPIRank()+NS[0]*NS[1];
192 }
193 }
194 }
195 int id=0;
196 for (int i=0; i<numElements; i++) {
197 nodes.push_back(id+NN[0]*NN[1]);
198 nodes.push_back(id);
199 nodes.push_back(id+1);
200 nodes.push_back(id+NN[0]*NN[1]+1);
201
202 nodes.push_back(id+NN[0]*(NN[1]+1));
203 nodes.push_back(id+NN[0]);
204 nodes.push_back(id+1+NN[0]);
205 nodes.push_back(id+NN[0]*(NN[1]+1)+1);
206 id++;
207 if ((i+1)%NperDim[0]==0)
208 id++;
209 if ((i+1)%(NperDim[0]*NperDim[1])==0)
210 id+=NN[0];
211 }
212 } else if (fsType==ripley::FaceElements) {
213 // ownership is not entirely correct but that is not critical.
214 // fix when there is time.
215 const IntVec NE = dom->getNumElementsPerDim();
216 int offset=0;
217 if (faces[0]>0) {
218 if (faces[3]==0) {
219 for (int k2=0; k2<NE[2]; k2++)
220 owner[(k2+1)*NE[1]-1]=dom->getMPIRank()+NS[0];
221 }
222 if (faces[5]==0) {
223 for (int k1=0; k1<NE[1]; k1++)
224 owner[NE[1]*(NE[2]-1)+k1]=dom->getMPIRank()+NS[0]*NS[1];
225 }
226 for (int k2=0; k2<NE[2]; k2++) {
227 for (int k1=0; k1<NE[1]; k1++) {
228 const int first=k2*NN[0]*NN[1]+k1*NN[0];
229 nodes.push_back(first+NN[0]*NN[1]);
230 nodes.push_back(first);
231 nodes.push_back(first+NN[0]);
232 nodes.push_back(first+NN[0]*(NN[1]+1));
233 }
234 }
235 offset+=faces[0];
236 }
237 if (faces[1]>0) {
238 if (faces[3]==0) {
239 for (int k2=0; k2<NE[2]; k2++)
240 owner[offset+(k2+1)*NE[1]-1]=dom->getMPIRank()+NS[0]+1;
241 }
242 if (faces[5]==0) {
243 for (int k1=0; k1<NE[1]; k1++)
244 owner[offset+NE[1]*(NE[2]-1)+k1]=dom->getMPIRank()+NS[0]*NS[1]+1;
245 }
246 for (int k2=0; k2<NE[2]; k2++) {
247 for (int k1=0; k1<NE[1]; k1++) {
248 const int first=k2*NN[0]*NN[1]+(k1+1)*NN[0]-1;
249 nodes.push_back(first+NN[0]*NN[1]);
250 nodes.push_back(first);
251 nodes.push_back(first+NN[0]);
252 nodes.push_back(first+NN[0]*(NN[1]+1));
253 }
254 }
255 offset+=faces[1];
256 }
257 if (faces[2]>0) {
258 if (faces[1]==0) {
259 for (int k2=0; k2<NE[2]; k2++)
260 owner[offset+(k2+1)*NE[0]-1]=dom->getMPIRank()+1;
261 }
262 if (faces[5]==0) {
263 for (int k0=0; k0<NE[0]; k0++)
264 owner[offset+NE[0]*(NE[2]-1)+k0]=dom->getMPIRank()+NS[0]*NS[1];
265 }
266 for (int k2=0; k2<NE[2]; k2++) {
267 for (int k0=0; k0<NE[0]; k0++) {
268 const int first=k2*NN[0]*NN[1]+k0;
269 nodes.push_back(first+NN[0]*NN[1]);
270 nodes.push_back(first);
271 nodes.push_back(first+1);
272 nodes.push_back(first+1+NN[0]*NN[1]);
273 }
274 }
275 offset+=faces[2];
276 }
277 if (faces[3]>0) {
278 if (faces[1]==0) {
279 for (int k2=0; k2<NE[2]; k2++)
280 owner[offset+(k2+1)*NE[0]-1]=dom->getMPIRank()+NS[0]+1;
281 }
282 if (faces[5]==0) {
283 for (int k0=0; k0<NE[0]; k0++)
284 owner[offset+NE[0]*(NE[2]-1)+k0]=dom->getMPIRank()+NS[0]*(NS[1]+1);
285 }
286 for (int k2=0; k2<NE[2]; k2++) {
287 for (int k0=0; k0<NE[0]; k0++) {
288 const int first=(k2+1)*NN[0]*NN[1]-NN[0]+k0;
289 nodes.push_back(first+NN[0]*NN[1]);
290 nodes.push_back(first);
291 nodes.push_back(first+1);
292 nodes.push_back(first+1+NN[0]*NN[1]);
293 }
294 }
295 offset+=faces[3];
296 }
297 if (faces[4]>0) {
298 if (faces[1]==0) {
299 for (int k1=0; k1<NE[1]; k1++)
300 owner[offset+(k1+1)*NE[0]-1]=dom->getMPIRank()+1;
301 }
302 if (faces[3]==0) {
303 for (int k0=0; k0<NE[0]; k0++)
304 owner[offset+NE[0]*(NE[1]-1)+k0]=dom->getMPIRank()+NS[0];
305 }
306 for (int k1=0; k1<NE[1]; k1++) {
307 for (int k0=0; k0<NE[0]; k0++) {
308 const int first=k1*NN[0]+k0;
309 nodes.push_back(first);
310 nodes.push_back(first+1);
311 nodes.push_back(first+NN[0]+1);
312 nodes.push_back(first+NN[0]);
313 }
314 }
315 offset+=faces[4];
316 }
317 if (faces[5]>0) {
318 if (faces[1]==0) {
319 for (int k1=0; k1<NE[1]; k1++)
320 owner[offset+(k1+1)*NE[0]-1]=dom->getMPIRank()+NS[0]*NS[1]+1;
321 }
322 if (faces[3]==0) {
323 for (int k0=0; k0<NE[0]; k0++)
324 owner[offset+NE[0]*(NE[1]-1)+k0]=dom->getMPIRank()+NS[0]*(NS[1]+1);
325 }
326 for (int k1=0; k1<NE[1]; k1++) {
327 for (int k0=0; k0<NE[0]; k0++) {
328 const int first=NN[0]*NN[1]*(NN[2]-1)+k1*NN[0]+k0;
329 nodes.push_back(first);
330 nodes.push_back(first+1);
331 nodes.push_back(first+NN[0]+1);
332 nodes.push_back(first+NN[0]);
333 }
334 }
335 }
336 }
337 }
338
339 buildMeshes();
340 }
341 return true;
342
343 #else // VISIT_PLUGIN
344 return false;
345 #endif
346 }
347
348 StringVec RipleyElements::getMeshNames() const
349 {
350 StringVec res;
351 if (nodeMesh)
352 res.push_back(nodeMesh->getName());
353 return res;
354 }
355
356 StringVec RipleyElements::getVarNames() const
357 {
358 StringVec res;
359 res.push_back(name + string("_Id"));
360 res.push_back(name + string("_Owner"));
361 res.push_back(name + string("_Tag"));
362 return res;
363 }
364
365 const IntVec& RipleyElements::getVarDataByName(const string varName) const
366 {
367 if (varName == name+string("_Id"))
368 return ID;
369 if (varName == name+string("_Owner"))
370 return owner;
371 if (varName == name+string("_Tag"))
372 return tag;
373
374 throw "Invalid variable name";
375 }
376
377 void RipleyElements::reorderArray(IntVec& v, const IntVec& idx,
378 int elementsPerIndex)
379 {
380 IntVec newArray(v.size());
381 IntVec::iterator arrIt = newArray.begin();
382 IntVec::const_iterator idxIt;
383 if (elementsPerIndex == 1) {
384 for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
385 *arrIt++ = v[*idxIt];
386 }
387 } else {
388 for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
389 int i = *idxIt;
390 copy(&v[i*elementsPerIndex], &v[(i+1)*elementsPerIndex], arrIt);
391 arrIt += elementsPerIndex;
392 }
393 }
394 v.swap(newArray);
395 }
396
397 IntVec RipleyElements::prepareGhostIndices(int ownIndex)
398 {
399 IntVec indexArray;
400 numGhostElements = 0;
401
402 // move indices of "ghost zones" to the end to be able to reorder
403 // data accordingly
404 for (int i=0; i<numElements; i++) {
405 if (owner[i] == ownIndex)
406 indexArray.push_back(i);
407 }
408
409 for (int i=0; i<numElements; i++) {
410 if (owner[i] != ownIndex) {
411 numGhostElements++;
412 indexArray.push_back(i);
413 }
414 }
415 return indexArray;
416 }
417
418 void RipleyElements::reorderGhostZones(int ownIndex)
419 {
420 IntVec indexArray = prepareGhostIndices(ownIndex);
421
422 // move "ghost data" to the end of the arrays
423 if (numGhostElements > 0) {
424 reorderArray(nodes, indexArray, nodesPerElement);
425 reorderArray(owner, indexArray, 1);
426 reorderArray(ID, indexArray, 1);
427 reorderArray(tag, indexArray, 1);
428 }
429 }
430
431 void RipleyElements::removeGhostZones(int ownIndex)
432 {
433 reorderGhostZones(ownIndex);
434
435 if (numGhostElements > 0) {
436 numElements -= numGhostElements;
437 nodes.resize(numElements*nodesPerElement);
438 owner.resize(numElements);
439 ID.resize(numElements);
440 tag.resize(numElements);
441 numGhostElements = 0;
442 }
443 }
444
445 void RipleyElements::buildMeshes()
446 {
447 // build a new mesh containing only the required nodes
448 if (numElements > 0) {
449 if (nodeMesh && nodeMesh->getNumNodes() > 0) {
450 RipleyNodes_ptr newMesh(new RipleyNodes(nodeMesh, nodes, name));
451 nodeMesh.swap(newMesh);
452 } else {
453 nodeMesh.reset(new RipleyNodes(originalMesh, nodes, name));
454 }
455 #ifdef _DEBUG
456 cout << nodeMesh->getName() << " has " << nodeMesh->getNumNodes()
457 << " nodes and " << numElements << " elements" << endl;
458 #endif
459 }
460 }
461
462 void RipleyElements::writeConnectivityVTK(ostream& os)
463 {
464 if (numElements > 0) {
465 const IntVec& gNI = nodeMesh->getGlobalNodeIndices();
466 IntVec::const_iterator it;
467 int count = 1;
468 for (it=nodes.begin(); it!=nodes.end(); it++, count++) {
469 os << gNI[*it];
470 if (count % nodesPerElement == 0)
471 os << endl;
472 else
473 os << " ";
474 }
475 }
476 }
477
478 #if USE_SILO
479 inline int toSiloElementType(int type)
480 {
481 switch (type) {
482 case ZONETYPE_BEAM: return DB_ZONETYPE_BEAM;
483 case ZONETYPE_HEX: return DB_ZONETYPE_HEX;
484 case ZONETYPE_POLYGON: return DB_ZONETYPE_POLYGON;
485 case ZONETYPE_QUAD: return DB_ZONETYPE_QUAD;
486 }
487 return 0;
488 }
489 #endif
490
491 bool RipleyElements::writeToSilo(DBfile* dbfile, const string& siloPath,
492 const StringVec& labels,
493 const StringVec& units, bool writeMeshData)
494 {
495 #if USE_SILO
496 if (numElements == 0)
497 return true;
498
499 int ret;
500
501 if (siloPath != "") {
502 ret = DBSetDir(dbfile, siloPath.c_str());
503 if (ret != 0)
504 return false;
505 }
506
507 // write out the full mesh in any case
508 nodeMesh->setSiloPath(siloPath);
509 string siloMeshNameStr = nodeMesh->getFullSiloName();
510 const char* siloMeshName = siloMeshNameStr.c_str();
511 int arraylen = numElements * nodesPerElement;
512 int eltype = toSiloElementType(type);
513
514 string varName = name + string("_zones");
515 ret = DBPutZonelist2(dbfile, varName.c_str(), numElements,
516 nodeMesh->getNumDims(), &nodes[0], arraylen, 0, 0,
517 numGhostElements, &eltype, &nodesPerElement, &numElements, 1, NULL);
518
519 if (ret == 0) {
520 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
521 DBoptlist* optList = NULL;
522 int nOpts = labels.size()+units.size();
523 if (nOpts>0) {
524 optList = DBMakeOptlist(nOpts);
525 if (labels.size()>0)
526 DBAddOption(optList, DBOPT_XLABEL, (void*)labels[0].c_str());
527 if (labels.size()>1)
528 DBAddOption(optList, DBOPT_YLABEL, (void*)labels[1].c_str());
529 if (labels.size()>2)
530 DBAddOption(optList, DBOPT_ZLABEL, (void*)labels[2].c_str());
531 if (units.size()>0)
532 DBAddOption(optList, DBOPT_XUNITS, (void*)units[0].c_str());
533 if (units.size()>1)
534 DBAddOption(optList, DBOPT_YUNITS, (void*)units[1].c_str());
535 if (units.size()>2)
536 DBAddOption(optList, DBOPT_ZUNITS, (void*)units[2].c_str());
537 }
538 ret = DBPutUcdmesh(dbfile, siloMeshName,
539 nodeMesh->getNumDims(), NULL, &coordbase[0],
540 nodeMesh->getNumNodes(), numElements, varName.c_str(),
541 /*"facelist"*/NULL, DB_FLOAT, optList);
542
543 if (optList)
544 DBFreeOptlist(optList);
545 }
546
547 if (ret != 0)
548 return false;
549
550 // write out the element-centered variables if enabled
551 if (writeMeshData) {
552 varName = name + string("_Id");
553 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
554 (float*)&ID[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
555 NULL);
556 if (ret == 0) {
557 varName = name + string("_Owner");
558 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
559 (float*)&owner[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
560 NULL);
561 }
562 }
563
564 // "Elements" is a special case
565 if (writeMeshData && name == "Elements") {
566 nodeMesh->writeToSilo(dbfile);
567 }
568
569 return (ret == 0);
570
571 #else // !USE_SILO
572 return false;
573 #endif
574 }
575
576 } // namespace weipa
577

  ViewVC Help
Powered by ViewVC 1.1.26