Actual source code: SieveBuilder.hh

  1: #ifndef included_ALE_SieveBuilder_hh
  2: #define included_ALE_SieveBuilder_hh

  4: #ifndef  included_ALE_ALE_hh
  5: #include <ALE.hh>
  6: #endif

  8: namespace ALE {
  9:   template<typename Bundle_>
 10:   class SieveBuilder {
 11:   public:
 12:     typedef Bundle_                                      bundle_type;
 13:     typedef typename bundle_type::sieve_type             sieve_type;
 14:     typedef typename bundle_type::arrow_section_type     arrow_section_type;
 15:     typedef std::vector<typename sieve_type::point_type> PointArray;
 16:     typedef std::pair<typename sieve_type::point_type, int> oPoint_type;
 17:     typedef std::vector<oPoint_type>                        oPointArray;
 18:   public:
 19:     static void buildHexFaces(Obj<sieve_type> sieve, Obj<arrow_section_type> orientation, int dim, std::map<int, int*>& curElement, std::map<int,PointArray>& bdVertices, std::map<int,oPointArray>& faces, typename sieve_type::point_type& cell, int& cellOrientation) {
 20:       int debug = sieve->debug();

 22:       if (debug > 1) {std::cout << "  Building hex faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;}
 23:       faces[dim].clear();
 24:       if (dim > 3) {
 25:         throw ALE::Exception("Cannot do hexes of dimension greater than three");
 26:       } else if (dim > 2) {
 27:         int nodes[24] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 5, 4,
 28:                          1, 2, 6, 5, 2, 3, 7, 6, 3, 0, 4, 7};

 30:         for(int b = 0; b < 6; b++) {
 31:           typename sieve_type::point_type face;
 32:           int o = 1;

 34:           bdVertices[dim-1].clear();
 35:           for(int c = 0; c < 4; c++) {
 36:             bdVertices[dim-1].push_back(bdVertices[dim][nodes[b*4+c]]);
 37:           }
 38:           if (debug > 1) {std::cout << "    boundary hex face " << b << std::endl;}
 39:           buildHexFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
 40:           if (debug > 1) {std::cout << "    added face " << face << std::endl;}
 41:           faces[dim].push_back(oPoint_type(face, o));
 42:         }
 43:       } else if (dim > 1) {
 44:         int boundarySize = bdVertices[dim].size();

 46:         for(int b = 0; b < boundarySize; b++) {
 47:           typename sieve_type::point_type face;
 48:           int o = 1;

 50:           bdVertices[dim-1].clear();
 51:           for(int c = 0; c < 2; c++) {
 52:             bdVertices[dim-1].push_back(bdVertices[dim][(b+c)%boundarySize]);
 53:           }
 54:           if (debug > 1) {
 55:             std::cout << "    boundary point " << bdVertices[dim][b] << std::endl;
 56:             std::cout << "      boundary vertices";
 57:             for(int c = 0; c < (int) bdVertices[dim-1].size(); c++) {
 58:               std::cout << " " << bdVertices[dim-1][c];
 59:             }
 60:             std::cout << std::endl;
 61:           }
 62:           buildHexFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
 63:           if (debug > 1) {std::cout << "    added face " << face << std::endl;}
 64:           faces[dim].push_back(oPoint_type(face, o));
 65:         }
 66:       } else {
 67:         if (debug > 1) {std::cout << "  Just set faces to boundary in 1d" << std::endl;}
 68:         typename PointArray::iterator bd_iter = bdVertices[dim].begin();
 69:         faces[dim].push_back(oPoint_type(*bd_iter, 0));++bd_iter;
 70:         faces[dim].push_back(oPoint_type(*bd_iter, 0));
 71:         //faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
 72:       }
 73:       if (debug > 1) {
 74:         for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
 75:           std::cout << "  face point " << f_iter->first << " orientation " << f_iter->second << std::endl;
 76:         }
 77:       }
 78:       // We always create the toplevel, so we could short circuit somehow
 79:       // Should not have to loop here since the meet of just 2 boundary elements is an element
 80:       typename oPointArray::iterator         f_itor = faces[dim].begin();
 81:       const typename sieve_type::point_type& start  = f_itor->first;
 82:       const typename sieve_type::point_type& next   = (++f_itor)->first;
 83:       Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);

 85:       if (preElement->size() > 0) {
 86:         cell = *preElement->begin();

 88:         const int size = faces[dim].size();
 89:         const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
 90:         int       wrap = size > 2 ? size-1 : 0;
 91:         int       indA = 0, indB = 0;

 93:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indA) {
 94:           if (start == *c_iter) break;
 95:         }
 96:         if (debug > 1) {std::cout << "    pointA " << start << " indA " << indA << std::endl;}
 97:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indB) {
 98:           if (next  == *c_iter) break;
 99:         }
100:         if (debug > 1) {std::cout << "    pointB " << next  << " indB " << indB << std::endl;}
101:         if ((indB - indA == 1) || (indA - indB == wrap)) {
102:           if (cellOrientation > 0) {
103:             cellOrientation = indA+1;
104:           } else {
105:             if (dim == 1) {
106:               cellOrientation = -2;
107:             } else {
108:               cellOrientation = -(indA+1);
109:             }
110:           }
111:         } else if ((indA - indB == 1) || (indB - indA == wrap)) {
112:           if (debug > 1) {std::cout << "      reversing cell orientation" << std::endl;}
113:           if (cellOrientation > 0) {
114:             cellOrientation = -(indA+1);
115:           } else {
116:             if (dim == 1) {
117:               cellOrientation = 1;
118:             } else {
119:               cellOrientation = indA+1;
120:             }
121:           }
122:         } else {
123:           throw ALE::Exception("Inconsistent orientation");
124:         }
125:         if (debug > 1) {std::cout << "  Found old cell " << cell << " orientation " << cellOrientation << std::endl;}
126:       } else {
127:         int color = 0;

129:         cell = typename sieve_type::point_type((*curElement[dim])++);
130:         for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
131:           MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(f_iter->first, cell);

133:           sieve->addArrow(f_iter->first, cell, color++);
134:           if (f_iter->second) {
135:             orientation->addPoint(arrow);
136:             orientation->updatePoint(arrow, &(f_iter->second));
137:             if (debug > 1) {std::cout << "    Orienting arrow (" << f_iter->first << ", " << cell << ") to " << f_iter->second << std::endl;}
138:           }
139:         }
140:         if (cellOrientation > 0) {
141:           cellOrientation = 1;
142:         } else {
143:           cellOrientation = -(dim+1);
144:         }
145:         if (debug > 1) {std::cout << "  Added cell " << cell << " dim " << dim << std::endl;}
146:       }
147:     };
148:     static void buildFaces(Obj<sieve_type> sieve, Obj<arrow_section_type> orientation, int dim, std::map<int, int*>& curElement, std::map<int,PointArray>& bdVertices, std::map<int,oPointArray>& faces, typename sieve_type::point_type& cell, int& cellOrientation) {
149:       int debug = sieve->debug();

151:       if (debug > 1) {
152:         if (cell >= 0) {
153:           std::cout << "  Building faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
154:         } else {
155:           std::cout << "  Building faces for boundary of undetermined cell (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
156:         }
157:       }
158:       if (dim == 0) return;
159:       faces[dim].clear();
160:       if (dim > 1) {
161:         int b = 0;
162:         // Use the cone construction
163:         for(typename PointArray::iterator b_itor = bdVertices[dim].begin(); b_itor != bdVertices[dim].end(); ++b_itor, ++b) {
164:           typename sieve_type::point_type face = -1;
165:           int                             o    = b%2 ? -bdVertices[dim].size() : 1;

167:           bdVertices[dim-1].clear();
168:           for(typename PointArray::iterator i_itor = bdVertices[dim].begin(); i_itor != bdVertices[dim].end(); ++i_itor) {
169:             if (i_itor != b_itor) {
170:               bdVertices[dim-1].push_back(*i_itor);
171:             }
172:           }
173:           if (debug > 1) {std::cout << "    boundary point " << *b_itor << std::endl;}
174:           buildFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
175:           if (debug > 1) {std::cout << "    added face " << face << std::endl;}
176:           faces[dim].push_back(oPoint_type(face, o));
177:         }
178:       } else {
179:         if (debug > 1) {std::cout << "  Just set faces to boundary in 1d" << std::endl;}
180:         typename PointArray::iterator bd_iter = bdVertices[dim].begin();
181:         faces[dim].push_back(oPoint_type(*bd_iter, 0));++bd_iter;
182:         faces[dim].push_back(oPoint_type(*bd_iter, 0));
183:         //faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
184:       }
185:       if (debug > 1) {
186:         for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
187:           std::cout << "  face point " << f_iter->first << " orientation " << f_iter->second << std::endl;
188:         }
189:       }
190:       // We always create the toplevel, so we could short circuit somehow
191:       // Should not have to loop here since the meet of just 2 boundary elements is an element
192:       typename oPointArray::iterator         f_itor = faces[dim].begin();
193:       const typename sieve_type::point_type& start  = f_itor->first;
194:       const typename sieve_type::point_type& next   = (++f_itor)->first;
195:       Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);

197:       if (preElement->size() > 0) {
198:         cell = *preElement->begin();

200:         const int size = faces[dim].size();
201:         const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
202:         int       wrap = size > 2 ? size-1 : 0;
203:         int       indA = 0, indB = 0;

205:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indA) {
206:           if (start == *c_iter) break;
207:         }
208:         if (debug > 1) {std::cout << "    pointA " << start << " indA " << indA << std::endl;}
209:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indB) {
210:           if (next  == *c_iter) break;
211:         }
212:         if (debug > 1) {std::cout << "    pointB " << next  << " indB " << indB << std::endl;}
213:         if ((indB - indA == 1) || (indA - indB == wrap)) {
214:           if (cellOrientation > 0) {
215:             cellOrientation = indA+1;
216:           } else {
217:             if (dim == 1) {
218:               cellOrientation = -2;
219:             } else {
220:               cellOrientation = -(indA+1);
221:             }
222:           }
223:         } else if ((indA - indB == 1) || (indB - indA == wrap)) {
224:           if (debug > 1) {std::cout << "      reversing cell orientation" << std::endl;}
225:           if (cellOrientation > 0) {
226:             cellOrientation = -(indA+1);
227:           } else {
228:             if (dim == 1) {
229:               cellOrientation = 1;
230:             } else {
231:               cellOrientation = indA+1;
232:             }
233:           }
234:         } else {
235:           throw ALE::Exception("Inconsistent orientation");
236:         }
237:         if (debug > 1) {std::cout << "  Found old cell " << cell << " orientation " << cellOrientation << std::endl;}
238:       } else {
239:         int color = 0;

241:         cell = typename sieve_type::point_type((*curElement[dim])++);
242:         for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
243:           MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(f_iter->first, cell);

245:           sieve->addArrow(f_iter->first, cell, color++);
246:           if (f_iter->second) {
247:             orientation->addPoint(arrow);
248:             orientation->updatePoint(arrow, &(f_iter->second));
249:             if (debug > 1) {std::cout << "    Orienting arrow (" << f_iter->first << ", " << cell << ") to " << f_iter->second << std::endl;}
250:           }
251:         }
252:         if (cellOrientation > 0) {
253:           cellOrientation = 1;
254:         } else {
255:           cellOrientation = -(dim+1);
256:         }
257:         if (debug > 1) {std::cout << "  Added cell " << cell << " dim " << dim << " orientation " << cellOrientation << std::endl;}
258:       }
259:     };
260: #if 0
261:     static void orientTriangle(const typename sieve_type::point_type cell, const int vertices[], const Obj<sieve_type>& sieve, const Obj<arrow_section_type>& orientation, int firstVertex[]) {
262:       const Obj<typename sieve_type::traits::coneSequence>&     cone = sieve->cone(cell);
263:       const typename sieve_type::traits::coneSequence::iterator end  = cone->end();
264:       int debug   = sieve->debug();
265:       int corners = 3;
266:       int e       = 0;

268:       if (debug > 1) {std::cout << "Orienting triangle " << cell << std::endl;}
269:       for(typename sieve_type::traits::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter, ++e) {
270:         if (debug > 1) {std::cout << "  edge " << *p_iter << std::endl;}
271:         const Obj<typename sieve_type::traits::coneSequence>& endpoints = sieve->cone(*p_iter);
272:         typename sieve_type::traits::coneSequence::iterator   vertex    = endpoints->begin();
273:         MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(*p_iter, cell);
274:         int                                                                           indA, indB, value;

276:         orientation->addPoint(arrow);
277:         for(indA = 0; indA < corners; indA++) {if (*vertex == vertices[indA]) break;}
278:         if (debug > 1) {std::cout << "    vertexA " << *vertex << " indA " << indA <<std::endl;}
279:         firstVertex[e] = *vertex;
280:         ++vertex;
281:         for(indB = 0; indB < corners; indB++) {if (*vertex == vertices[indB]) break;}
282:         if (debug > 1) {std::cout << "    vertexB " << *vertex << " indB " << indB <<std::endl;}
283:         if ((indA == corners) || (indB == corners) || (indA == indB)) {throw ALE::Exception("Invalid edge endpoints");}
284:         if ((indB - indA == 1) || (indA - indB == 2)) {
285:           value =  1;
286:         } else {
287:           value = -1;
288:           firstVertex[e] = *vertex;
289:         }
290:         if (debug > 1) {std::cout << "  value " << value <<std::endl;}
291:         orientation->updatePoint(arrow, &value);
292:       }
293:     };
294: #endif
297:     // Build a topology from a connectivity description
298:     //   (0, 0)        ... (0, numCells-1):  dim-dimensional cells
299:     //   (0, numCells) ... (0, numVertices): vertices
300:     // The other cells are numbered as they are requested
301:     static void buildTopology(Obj<sieve_type> sieve, int dim, int numCells, int cells[], int numVertices, bool interpolate = true, int corners = -1, int firstVertex = -1, Obj<arrow_section_type> orientation = NULL, int firstCell = 0) {
302:       int debug = sieve->debug();

304:       ALE_LOG_EVENT_BEGIN;
305:       if (sieve->commRank() != 0) {
306:         ALE_LOG_EVENT_END;
307:         return;
308:       }
309:       if (firstVertex < 0) firstVertex = numCells;
310:       // Create a map from dimension to the current element number for that dimension
311:       std::map<int,int*>       curElement;
312:       std::map<int,PointArray> bdVertices;
313:       std::map<int,PointArray> faces;
314:       std::map<int,oPointArray> oFaces;
315:       int                      curCell    = firstCell;
316:       int                      curVertex  = firstVertex;
317:       int                      newElement = firstVertex > firstCell ? firstVertex + numVertices : firstCell + numCells;
318:       int                      o          = 1;

320:       if (corners < 0) corners = dim+1;
321:       curElement[0]   = &curVertex;
322:       curElement[dim] = &curCell;
323:       for(int d = 1; d < dim; d++) {
324:         curElement[d] = &newElement;
325:       }
326:       for(int c = 0; c < numCells; c++) {
327:         typename sieve_type::point_type cell(c);

329:         // Build the cell
330:         if (interpolate) {
331:           bdVertices[dim].clear();
332:           for(int b = 0; b < corners; b++) {
333:             // This ordering produces the same vertex order as the uninterpolated mesh
334:             //typename sieve_type::point_type vertex(cells[c*corners+(b+corners-1)%corners]+firstVertex);
335:             typename sieve_type::point_type vertex(cells[c*corners+b]+firstVertex);

337:             if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
338:             bdVertices[dim].push_back(vertex);
339:           }
340:           if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}

342:           if (corners != dim+1) {
343:             buildHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
344:           } else {
345:             buildFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
346:           }
347: #if 0
348:           if ((dim == 2) && (!orientation.isNull())) {
349:             typename sieve_type::point_type eVertices[3];
350:             typename sieve_type::point_type fVertices[3];

352:             for(int v = 0; v < 3; ++v) {
353:               fVertices[v] = cells[c*corners+v]+numCells;
354:             }
355:             orientTriangle(cell, fVertices, sieve, orientation, eVertices);
356:           } else if ((dim == 3) && (!orientation.isNull())) {
357:             // The order of vertices in cells[] orients the cell itself
358:             if (debug > 1) {std::cout << "Orienting tetrahedron " << cell << std::endl;}
359:             const Obj<typename sieve_type::traits::coneSequence>&     cone = sieve->cone(cell);
360:             const typename sieve_type::traits::coneSequence::iterator end  = cone->end();
361:             int f = 0;

363:             for(typename sieve_type::traits::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter, ++f) {
364:               if (debug > 1) {std::cout << "  face " << *p_iter << std::endl;}
365:               const Obj<typename sieve_type::traits::coneSequence>&     fCone = sieve->cone(*p_iter);
366:               const typename sieve_type::traits::coneSequence::iterator fEnd  = fCone->end();
367:               typename sieve_type::point_type fVertices[3];
368:               typename sieve_type::point_type eVertices[3];

370:               // We will choose the orientation such that the normals are outward
371:               for(int v = 0, i = 0; v < corners; ++v) {
372:                 if (v == f) continue;
373:                 fVertices[i++] = cells[c*corners+v]+numCells;
374:               }
375:               if (f%2) {
376:                 int tmp      = fVertices[0];
377:                 fVertices[0] = fVertices[1];
378:                 fVertices[1] = tmp;
379:               }
380:               orientTriangle(*p_iter, fVertices, sieve, orientation, eVertices);
381:               MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> fArrow(*p_iter, cell);
382:               int                                                                           indC, indD, indE, value;

384:               orientation->addPoint(fArrow);
385:               for(indC = 0; indC < corners; indC++) {if (eVertices[0] == fVertices[indC]) break;}
386:               if (debug > 1) {std::cout << "    vertexC " << eVertices[0] << " indC " << indC <<std::endl;}
387:               for(indD = 0; indD < corners; indD++) {if (eVertices[1] == fVertices[indD]) break;}
388:               if (debug > 1) {std::cout << "    vertexD " << eVertices[1] << " indD " << indD <<std::endl;}
389:               for(indE = 0; indE < corners; indE++) {if (eVertices[2] == fVertices[indE]) break;}
390:               if (debug > 1) {std::cout << "    vertexE " << eVertices[2] << " indE " << indE <<std::endl;}
391:               if ((indC == corners) || (indD == corners) || (indE == corners) ||
392:                   (indC == indD) || (indD == indE) || (indE == indC)) {throw ALE::Exception("Invalid face corners");}
393:               if ((indD - indC == 1) || (indC - indD == 2)) {
394:                 if (!((indE - indD == 1) || (indD - indE == 2)) ||
395:                     !((indC - indE == 1) || (indE - indC == 2))) {throw ALE::Exception("Invalid order");}
396:                 value =  1;
397:               } else {
398:                 value = -1;
399:               }
400:               if (debug > 1) {std::cout << "  value " << value <<std::endl;}
401:               orientation->updatePoint(fArrow, &value);
402:               orientation->view("Intermediate orientation");
403:             }
404:           }
405: #endif
406:         } else {
407:           for(int b = 0; b < corners; b++) {
408:             sieve->addArrow(typename sieve_type::point_type(cells[c*corners+b]+firstVertex), cell, b);
409:           }
410:           if (debug) {
411:             if (debug > 1) {
412:               for(int b = 0; b < corners; b++) {
413:                 std::cout << "  Adding vertex " << typename sieve_type::point_type(cells[c*corners+b]+firstVertex) << std::endl;
414:               }
415:             }
416:             if ((numCells < 10000) || (c%1000 == 0)) {
417:               std::cout << "Adding cell " << cell << " dim " << dim << std::endl;
418:             }
419:           }
420:         }
421:       }
422:       ALE_LOG_EVENT_END;
423:     };
424:     static void buildCoordinates(const Obj<Bundle_>& bundle, const int embedDim, const double coords[]) {
425:       const Obj<typename Bundle_::real_section_type>& coordinates = bundle->getRealSection("coordinates");
426:       const Obj<typename Bundle_::label_sequence>&    vertices    = bundle->depthStratum(0);
427:       const int numCells = bundle->heightStratum(0)->size();
428:       const int debug    = bundle->debug();

430:       bundle->setupCoordinates(coordinates);
431:       coordinates->setFiberDimension(vertices, embedDim);
432:       bundle->allocate(coordinates);
433:       for(typename Bundle_::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
434:         coordinates->updatePoint(*v_iter, &(coords[(*v_iter - numCells)*embedDim]));
435:         if (debug) {
436:           if ((numCells < 10000) || ((*v_iter)%1000 == 0)) {
437:             std::cout << "Adding coordinates for vertex " << *v_iter << std::endl;
438:           }
439:         }
440:       }
441:     };
444:     // Build a topology from a connectivity description
445:     //   (0, 0)        ... (0, numCells-1):  dim-dimensional cells
446:     //   (0, numCells) ... (0, numVertices): vertices
447:     // The other cells are numbered as they are requested
448:     static void buildTopologyMultiple(Obj<sieve_type> sieve, int dim, int numCells, int cells[], int numVertices, bool interpolate = true, int corners = -1, int firstVertex = -1, Obj<arrow_section_type> orientation = NULL) {
449:       int debug = sieve->debug();

451:       ALE_LOG_EVENT_BEGIN;
452:       int *cellOffset = new int[sieve->commSize()+1];
453:       cellOffset[0] = 0;
454:       MPI_Allgather(&numCells, 1, MPI_INT, &cellOffset[1], 1, MPI_INT, sieve->comm());
455:       for(int p = 1; p <= sieve->commSize(); ++p) cellOffset[p] += cellOffset[p-1];
456:       int *vertexOffset = new int[sieve->commSize()+1];
457:       vertexOffset[0] = 0;
458:       MPI_Allgather(&numVertices, 1, MPI_INT, &vertexOffset[1], 1, MPI_INT, sieve->comm());
459:       for(int p = 1; p <= sieve->commSize(); ++p) vertexOffset[p] += vertexOffset[p-1];
460:       if (firstVertex < 0) firstVertex = cellOffset[sieve->commSize()] + vertexOffset[sieve->commRank()];
461:       // Estimate the number of intermediates as (V+C)*(dim-1)
462:       //   Should include a check for running over the namespace
463:       // Create a map from dimension to the current element number for that dimension
464:       std::map<int,int*>       curElement;
465:       std::map<int,PointArray> bdVertices;
466:       std::map<int,PointArray> faces;
467:       std::map<int,oPointArray> oFaces;
468:       int                      curCell    = cellOffset[sieve->commRank()];
469:       int                      curVertex  = firstVertex;
470:       int                      newElement = firstVertex+vertexOffset[sieve->commSize()] + (cellOffset[sieve->commRank()] + vertexOffset[sieve->commRank()])*(dim-1);
471:       int                      o          = 1;

473:       if (corners < 0) corners = dim+1;
474:       curElement[0]   = &curVertex;
475:       curElement[dim] = &curCell;
476:       for(int d = 1; d < dim; d++) {
477:         curElement[d] = &newElement;
478:       }
479:       for(int c = 0; c < numCells; c++) {
480:         typename sieve_type::point_type cell(c);

482:         // Build the cell
483:         if (interpolate) {
484:           bdVertices[dim].clear();
485:           for(int b = 0; b < corners; b++) {
486:             // This ordering produces the same vertex order as the uninterpolated mesh
487:             //typename sieve_type::point_type vertex(cells[c*corners+(b+corners-1)%corners]+firstVertex);
488:             typename sieve_type::point_type vertex(cells[c*corners+b]+firstVertex);

490:             if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
491:             bdVertices[dim].push_back(vertex);
492:           }
493:           if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}

495:           if (corners != dim+1) {
496:             buildHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
497:           } else {
498:             buildFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
499:           }
500:         } else {
501:           for(int b = 0; b < corners; b++) {
502:             sieve->addArrow(typename sieve_type::point_type(cells[c*corners+b]+firstVertex), cell, b);
503:           }
504:           if (debug) {
505:             if (debug > 1) {
506:               for(int b = 0; b < corners; b++) {
507:                 std::cout << "  Adding vertex " << typename sieve_type::point_type(cells[c*corners+b]+firstVertex) << std::endl;
508:               }
509:             }
510:             if ((numCells < 10000) || (c%1000 == 0)) {
511:               std::cout << "Adding cell " << cell << " dim " << dim << std::endl;
512:             }
513:           }
514:         }
515:       }

517:       if (newElement >= firstVertex+vertexOffset[sieve->commSize()] + (cellOffset[sieve->commRank()+1] + vertexOffset[sieve->commRank()+1])*(dim-1)) {
518:         throw ALE::Exception("Namespace violation during intermediate element construction");
519:       }
520:       delete [] cellOffset;
521:       delete [] vertexOffset;
522:       ALE_LOG_EVENT_END;
523:     };
524:     static void buildCoordinatesMultiple(const Obj<Bundle_>& bundle, const int embedDim, const double coords[]) {
525:       const Obj<typename Bundle_::real_section_type>& coordinates = bundle->getRealSection("coordinates");
526:       const Obj<typename Bundle_::label_sequence>&    vertices    = bundle->depthStratum(0);
527:       const int numCells = bundle->heightStratum(0)->size(), numVertices = vertices->size();
528:       const int debug    = bundle->debug();
529:       int       numGlobalCells, offset;

531:       MPI_Allreduce((void *) &numCells, &numGlobalCells, 1, MPI_INT, MPI_SUM, bundle->comm());
532:       MPI_Scan((void *) &numVertices, &offset, 1, MPI_INT, MPI_SUM, bundle->comm());
533:       offset += numGlobalCells - numVertices;
534:       coordinates->setFiberDimension(vertices, embedDim);
535:       bundle->allocate(coordinates);
536:       for(typename Bundle_::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
537:         coordinates->updatePoint(*v_iter, &(coords[(*v_iter - offset)*embedDim]));
538:         if (debug) {
539:           if ((numCells < 10000) || ((*v_iter)%1000 == 0)) {
540:             std::cout << "Adding coordinates for vertex " << *v_iter << std::endl;
541:           }
542:         }
543:       }
544:     };
545:   };
546: }
547: #endif