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:       ALE_LOG_EVENT_BEGIN;
303:       if (sieve->commRank() != 0) {
304:         ALE_LOG_EVENT_END;
305:         return;
306:       }
307:       buildTopology_private(sieve, dim, numCells, cells, numVertices, interpolate, corners, firstVertex, orientation, firstCell);
308:       ALE_LOG_EVENT_END;
309:     };
310:     static void buildTopology_private(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) {
311:       int debug = sieve->debug();

313:       if (firstVertex < 0) firstVertex = numCells;
314:       // Create a map from dimension to the current element number for that dimension
315:       std::map<int,int*>       curElement;
316:       std::map<int,PointArray> bdVertices;
317:       std::map<int,PointArray> faces;
318:       std::map<int,oPointArray> oFaces;
319:       int                      curCell    = firstCell;
320:       int                      curVertex  = firstVertex;
321:       int                      newElement = firstVertex > firstCell ? firstVertex + numVertices : firstCell + numCells;
322:       int                      o          = 1;

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

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

341:             if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
342:             bdVertices[dim].push_back(vertex);
343:           }
344:           if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}

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

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

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

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

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

433:       bundle->setupCoordinates(coordinates);
434:       coordinates->setFiberDimension(vertices, embedDim);
435:       bundle->allocate(coordinates);
436:       for(typename Bundle_::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
437:         coordinates->updatePoint(*v_iter, &(coords[(*v_iter - numCells)*embedDim]));
438:         if (debug) {
439:           if ((numCells < 10000) || ((*v_iter)%1000 == 0)) {
440:             std::cout << "Adding coordinates for vertex " << *v_iter << std::endl;
441:           }
442:         }
443:       }
444:     };
447:     // Build a topology from a connectivity description
448:     //   (0, 0)        ... (0, numCells-1):  dim-dimensional cells
449:     //   (0, numCells) ... (0, numVertices): vertices
450:     // The other cells are numbered as they are requested
451:     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) {
452:       int debug = sieve->debug();

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

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

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

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

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

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

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