Actual source code: CoSieve.hh

  1: #ifndef included_ALE_CoSieve_hh
  2: #define included_ALE_CoSieve_hh

  4: #ifndef  included_ALE_Sieve_hh
  5: #include <Sieve.hh>
  6: #endif


 10: namespace ALE {
 11:   class ParallelObject {
 12:   protected:
 13:     int         _debug;
 14:     MPI_Comm    _comm;
 15:     int         _commRank;
 16:     int         _commSize;
 17:     std::string _name;
 18:     PetscObject _petscObj;
 19:     void __init(MPI_Comm comm) {
 20:       static PetscCookie objType = -1;
 21:       //const char        *id_name = ALE::getClassName<T>();
 22:       const char        *id_name = "ParallelObject";
 23:       PetscErrorCode     ierr;

 25:       if (objType < 0) {
 26:         PetscLogClassRegister(&objType, id_name);CHKERROR(ierr, "Error in PetscLogClassRegister");
 27:       }
 28:       this->_comm = comm;
 29:       MPI_Comm_rank(this->_comm, &this->_commRank); CHKERROR(ierr, "Error in MPI_Comm_rank");
 30:       MPI_Comm_size(this->_comm, &this->_commSize); CHKERROR(ierr, "Error in MPI_Comm_size");
 31:       PetscObjectCreateGeneric(this->_comm, objType, id_name, &this->_petscObj);CHKERROR(ierr, "Error in PetscObjectCreate");
 32:       //ALE::restoreClassName<T>(id_name);
 33:     };
 34:   public:
 35:     ParallelObject(MPI_Comm comm = PETSC_COMM_SELF, const int debug = 0) : _debug(debug), _petscObj(NULL) {__init(comm);}
 36:     virtual ~ParallelObject() {
 37:       if (this->_petscObj) {
 38:         PetscErrorCode PetscObjectDestroy(this->_petscObj); CHKERROR(ierr, "Failed in PetscObjectDestroy");
 39:         this->_petscObj = NULL;
 40:       }
 41:     };
 42:   public:
 43:     int                debug()    const {return this->_debug;};
 44:     void               setDebug(const int debug) {this->_debug = debug;};
 45:     MPI_Comm           comm()     const {return this->_comm;};
 46:     int                commSize() const {return this->_commSize;};
 47:     int                commRank() const {return this->_commRank;}
 48:     PetscObject        petscObj() const {return this->_petscObj;};
 49:     const std::string& getName() const {return this->_name;};
 50:     void               setName(const std::string& name) {this->_name = name;};
 51:   };

 53:   namespace New {
 54:     template<typename Sieve_>
 55:     class SieveBuilder {
 56:     public:
 57:       typedef Sieve_                                       sieve_type;
 58:       typedef std::vector<typename sieve_type::point_type> PointArray;
 59:     public:
 60:       static void buildHexFaces(Obj<sieve_type> sieve, int dim, std::map<int, int*>& curElement, std::map<int,PointArray>& bdVertices, std::map<int,PointArray>& faces, typename sieve_type::point_type& cell) {
 61:         int debug = sieve->debug();

 63:         if (debug > 1) {std::cout << "  Building hex faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;}
 64:         if (dim > 3) {
 65:           throw ALE::Exception("Cannot do hexes of dimension greater than three");
 66:         } else if (dim > 2) {
 67:           int nodes[24] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 5, 4,
 68:                            1, 2, 6, 5, 2, 3, 7, 6, 3, 0, 4, 7};

 70:           for(int b = 0; b < 6; b++) {
 71:             typename sieve_type::point_type face;

 73:             for(int c = 0; c < 4; c++) {
 74:               bdVertices[dim-1].push_back(bdVertices[dim][nodes[b*4+c]]);
 75:             }
 76:             if (debug > 1) {std::cout << "    boundary hex face " << b << std::endl;}
 77:             buildHexFaces(sieve, dim-1, curElement, bdVertices, faces, face);
 78:             if (debug > 1) {std::cout << "    added face " << face << std::endl;}
 79:             faces[dim].push_back(face);
 80:           }
 81:         } else if (dim > 1) {
 82:           int boundarySize = bdVertices[dim].size();

 84:           for(int b = 0; b < boundarySize; b++) {
 85:             typename sieve_type::point_type face;

 87:             for(int c = 0; c < 2; c++) {
 88:               bdVertices[dim-1].push_back(bdVertices[dim][(b+c)%boundarySize]);
 89:             }
 90:             if (debug > 1) {std::cout << "    boundary point " << bdVertices[dim][b] << std::endl;}
 91:             buildHexFaces(sieve, dim-1, curElement, bdVertices, faces, face);
 92:             if (debug > 1) {std::cout << "    added face " << face << std::endl;}
 93:             faces[dim].push_back(face);
 94:           }
 95:         } else {
 96:           if (debug > 1) {std::cout << "  Just set faces to boundary in 1d" << std::endl;}
 97:           faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
 98:         }
 99:         if (debug > 1) {
100:           for(typename PointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
101:             std::cout << "  face point " << *f_iter << std::endl;
102:           }
103:         }
104:         // We always create the toplevel, so we could short circuit somehow
105:         // Should not have to loop here since the meet of just 2 boundary elements is an element
106:         typename PointArray::iterator          f_itor = faces[dim].begin();
107:         const typename sieve_type::point_type& start  = *f_itor;
108:         const typename sieve_type::point_type& next   = *(++f_itor);
109:         Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);

111:         if (preElement->size() > 0) {
112:           cell = *preElement->begin();
113:           if (debug > 1) {std::cout << "  Found old cell " << cell << std::endl;}
114:         } else {
115:           int color = 0;

117:           cell = typename sieve_type::point_type((*curElement[dim])++);
118:           for(typename PointArray::iterator f_itor = faces[dim].begin(); f_itor != faces[dim].end(); ++f_itor) {
119:             sieve->addArrow(*f_itor, cell, color++);
120:           }
121:           if (debug > 1) {std::cout << "  Added cell " << cell << " dim " << dim << std::endl;}
122:         }
123:       };
124:       static void buildFaces(Obj<sieve_type> sieve, int dim, std::map<int, int*>& curElement, std::map<int,PointArray>& bdVertices, std::map<int,PointArray>& faces, typename sieve_type::point_type& cell) {
125:         int debug = sieve->debug();

127:         if (debug > 1) {
128:           if (cell >= 0) {
129:             std::cout << "  Building faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
130:           } else {
131:             std::cout << "  Building faces for boundary of undetermined cell (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
132:           }
133:         }
134:         faces[dim].clear();
135:         if (dim > 1) {
136:           // Use the cone construction
137:           for(typename PointArray::iterator b_itor = bdVertices[dim].begin(); b_itor != bdVertices[dim].end(); ++b_itor) {
138:             typename sieve_type::point_type face   = -1;

140:             bdVertices[dim-1].clear();
141:             for(typename PointArray::iterator i_itor = bdVertices[dim].begin(); i_itor != bdVertices[dim].end(); ++i_itor) {
142:               if (i_itor != b_itor) {
143:                 bdVertices[dim-1].push_back(*i_itor);
144:               }
145:             }
146:             if (debug > 1) {std::cout << "    boundary point " << *b_itor << std::endl;}
147:             buildFaces(sieve, dim-1, curElement, bdVertices, faces, face);
148:             if (debug > 1) {std::cout << "    added face " << face << std::endl;}
149:             faces[dim].push_back(face);
150:           }
151:         } else {
152:           if (debug > 1) {std::cout << "  Just set faces to boundary in 1d" << std::endl;}
153:           faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
154:         }
155:         if (debug > 1) {
156:           for(typename PointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
157:             std::cout << "  face point " << *f_iter << std::endl;
158:           }
159:         }
160:         // We always create the toplevel, so we could short circuit somehow
161:         // Should not have to loop here since the meet of just 2 boundary elements is an element
162:         typename PointArray::iterator          f_itor = faces[dim].begin();
163:         const typename sieve_type::point_type& start  = *f_itor;
164:         const typename sieve_type::point_type& next   = *(++f_itor);
165:         Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);

167:         if (preElement->size() > 0) {
168:           cell = *preElement->begin();
169:           if (debug > 1) {std::cout << "  Found old cell " << cell << std::endl;}
170:         } else {
171:           int color = 0;

173:           cell = typename sieve_type::point_type((*curElement[dim])++);
174:           for(typename PointArray::iterator f_itor = faces[dim].begin(); f_itor != faces[dim].end(); ++f_itor) {
175:             sieve->addArrow(*f_itor, cell, color++);
176:           }
177:           if (debug > 1) {std::cout << "  Added cell " << cell << " dim " << dim << std::endl;}
178:         }
179:       };

183:       // Build a topology from a connectivity description
184:       //   (0, 0)        ... (0, numCells-1):  dim-dimensional cells
185:       //   (0, numCells) ... (0, numVertices): vertices
186:       // The other cells are numbered as they are requested
187:       static void buildTopology(Obj<sieve_type> sieve, int dim, int numCells, int cells[], int numVertices, bool interpolate = true, int corners = -1, int firstVertex = -1) {
188:         int debug = sieve->debug();

190:         ALE_LOG_EVENT_BEGIN;
191:         if (sieve->commRank() != 0) {
192:           ALE_LOG_EVENT_END;
193:           return;
194:         }
195:         if (firstVertex < 0) firstVertex = numCells;
196:         // Create a map from dimension to the current element number for that dimension
197:         std::map<int,int*>       curElement;
198:         std::map<int,PointArray> bdVertices;
199:         std::map<int,PointArray> faces;
200:         int                      curCell    = 0;
201:         int                      curVertex  = firstVertex;
202:         int                      newElement = firstVertex+numVertices;

204:         if (corners < 0) corners = dim+1;
205:         curElement[0]   = &curVertex;
206:         curElement[dim] = &curCell;
207:         for(int d = 1; d < dim; d++) {
208:           curElement[d] = &newElement;
209:         }
210:         for(int c = 0; c < numCells; c++) {
211:           typename sieve_type::point_type cell(c);

213:           // Build the cell
214:           if (interpolate) {
215:             bdVertices[dim].clear();
216:             for(int b = 0; b < corners; b++) {
217:               typename sieve_type::point_type vertex(cells[c*corners+b]+firstVertex);

219:               if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
220:               bdVertices[dim].push_back(vertex);
221:             }
222:             if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}

224:             if (corners != dim+1) {
225:               buildHexFaces(sieve, dim, curElement, bdVertices, faces, cell);
226:             } else {
227:               buildFaces(sieve, dim, curElement, bdVertices, faces, cell);
228:             }
229:           } else {
230:             for(int b = 0; b < corners; b++) {
231:               sieve->addArrow(typename sieve_type::point_type(cells[c*corners+b]+firstVertex), cell, b);
232:             }
233:             if (debug) {
234:               if (debug > 1) {
235:                 for(int b = 0; b < corners; b++) {
236:                   std::cout << "  Adding vertex " << typename sieve_type::point_type(cells[c*corners+b]+firstVertex) << std::endl;
237:                 }
238:               }
239:               std::cout << "Adding cell " << cell << " dim " << dim << std::endl;
240:             }
241:           }
242:         }
243:         ALE_LOG_EVENT_END;
244:       };
245:       template<typename Section>
246:       static void buildCoordinates(const Obj<Section>& coords, const int embedDim, const double coordinates[]) {
247:         const typename Section::patch_type                          patch    = 0;
248:         const Obj<typename Section::topology_type::label_sequence>& vertices = coords->getTopology()->depthStratum(patch, 0);
249:         const int numCells = coords->getTopology()->heightStratum(patch, 0)->size();

251:         coords->setName("coordinates");
252:         coords->setFiberDimensionByDepth(patch, 0, embedDim);
253:         coords->allocate();
254:         for(typename Section::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
255:           coords->update(patch, *v_iter, &(coordinates[(*v_iter - numCells)*embedDim]));
256:         }
257:       };
258:     };

260:     // A Topology is a collection of Sieves
261:     //   Each Sieve has a label, which we call a \emph{patch}
262:     //   The collection itself we call a \emph{sheaf}
263:     //   The main operation we provide in Topology is the creation of a \emph{label}
264:     //     A label is a bidirectional mapping of Sieve points to integers, implemented with a Sifter
265:     template<typename Patch_, typename Sieve_>
266:     class Topology : public ALE::ParallelObject {
267:     public:
268:       typedef Patch_                                                patch_type;
269:       typedef Sieve_                                                sieve_type;
270:       typedef typename sieve_type::point_type                       point_type;
271:       typedef typename std::map<patch_type, Obj<sieve_type> >       sheaf_type;
272:       typedef typename ALE::Sifter<int, point_type, int>            patch_label_type;
273:       typedef typename std::map<patch_type, Obj<patch_label_type> > label_type;
274:       typedef typename std::map<patch_type, int>                    max_label_type;
275:       typedef typename std::map<const std::string, label_type>      labels_type;
276:       typedef typename patch_label_type::supportSequence            label_sequence;
277:       typedef typename std::set<point_type>                         point_set_type;
278:       typedef typename ALE::Sifter<int,point_type,point_type>       send_overlap_type;
279:       typedef typename ALE::Sifter<point_type,int,point_type>       recv_overlap_type;
280:     protected:
281:       sheaf_type     _sheaf;
282:       labels_type    _labels;
283:       int            _maxHeight;
284:       max_label_type _maxHeights;
285:       int            _maxDepth;
286:       max_label_type _maxDepths;
287:       bool           _calculatedOverlap;
288:       Obj<send_overlap_type> _sendOverlap;
289:       Obj<recv_overlap_type> _recvOverlap;
290:       Obj<send_overlap_type> _distSendOverlap;
291:       Obj<recv_overlap_type> _distRecvOverlap;
292:       // Work space
293:       Obj<point_set_type>    _modifiedPoints;
294:     public:
295:       Topology(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1), _calculatedOverlap(false) {
296:         this->_sendOverlap    = new send_overlap_type(this->comm(), this->debug());
297:         this->_recvOverlap    = new recv_overlap_type(this->comm(), this->debug());
298:         this->_modifiedPoints = new point_set_type();
299:       };
300:     public: // Verifiers
301:       void checkPatch(const patch_type& patch) {
302:         if (this->_sheaf.find(patch) == this->_sheaf.end()) {
303:           ostringstream msg;
304:           msg << "Invalid topology patch: " << patch << std::endl;
305:           throw ALE::Exception(msg.str().c_str());
306:         }
307:       };
308:       void checkLabel(const std::string& name, const patch_type& patch) {
309:         this->checkPatch(patch);
310:         if ((this->_labels.find(name) == this->_labels.end()) || (this->_labels[name].find(patch) == this->_labels[name].end())) {
311:           ostringstream msg;
312:           msg << "Invalid label name: " << name << " for patch " << patch << std::endl;
313:           throw ALE::Exception(msg.str().c_str());
314:         }
315:       };
316:       bool hasPatch(const patch_type& patch) {
317:         if (this->_sheaf.find(patch) != this->_sheaf.end()) {
318:           return true;
319:         }
320:         return false;
321:       };
322:       bool hasLabel(const std::string& name, const patch_type& patch) {
323:         if ((this->_labels.find(name) != this->_labels.end()) && (this->_labels[name].find(patch) != this->_labels[name].end())) {
324:           return true;
325:         }
326:         return false;
327:       };
328:     public: // Accessors
329:       const Obj<sieve_type>& getPatch(const patch_type& patch) {
330:         this->checkPatch(patch);
331:         return this->_sheaf[patch];
332:       };
333:       void setPatch(const patch_type& patch, const Obj<sieve_type>& sieve) {
334:         this->_sheaf[patch] = sieve;
335:       };
336:       int getValue (const Obj<patch_label_type>& label, const point_type& point, const int defValue = 0) {
337:         const Obj<typename patch_label_type::coneSequence>& cone = label->cone(point);

339:         if (cone->size() == 0) return defValue;
340:         return *cone->begin();
341:       };
342:       template<typename InputPoints>
343:       int getMaxValue (const Obj<patch_label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
344:         int maxValue = defValue;

346:         for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
347:           maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
348:         }
349:         return maxValue;
350:       };
351:       void setValue(const Obj<patch_label_type>& label, const point_type& point, const int value) {
352:         label->setCone(value, point);
353:       };
354:       const Obj<patch_label_type>& createLabel(const patch_type& patch, const std::string& name) {
355:         this->checkPatch(patch);
356:         this->_labels[name][patch] = new patch_label_type(this->comm(), this->debug());
357:         return this->_labels[name][patch];
358:       };
359:       const Obj<patch_label_type>& getLabel(const patch_type& patch, const std::string& name) {
360:         this->checkLabel(name, patch);
361:         return this->_labels[name][patch];
362:       };
363:       const Obj<label_sequence>& getLabelStratum(const patch_type& patch, const std::string& name, int label) {
364:         this->checkLabel(name, patch);
365:         return this->_labels[name][patch]->support(label);
366:       };
367:       const sheaf_type& getPatches() {
368:         return this->_sheaf;
369:       };
370:       const labels_type& getLabels() {
371:         return this->_labels;
372:       };
373:       void clear() {
374:         this->_sheaf.clear();
375:         this->_labels.clear();
376:         this->_maxHeight = -1;
377:         this->_maxHeights.clear();
378:         this->_maxDepth = -1;
379:         this->_maxDepths.clear();
380:       };
381:       const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
382:       void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
383:       const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
384:       void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
385:       const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
386:       void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
387:       const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
388:       void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
389:     public:
390:       template<class InputPoints>
391:       void computeHeight(const Obj<patch_label_type>& height, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxHeight) {
392:         this->_modifiedPoints->clear();

394:         for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
395:           // Compute the max height of the points in the support of p, and add 1
396:           int h0 = this->getValue(height, *p_iter, -1);
397:           int h1 = this->getMaxValue(height, sieve->support(*p_iter), -1) + 1;

399:           if(h1 != h0) {
400:             this->setValue(height, *p_iter, h1);
401:             if (h1 > maxHeight) maxHeight = h1;
402:             this->_modifiedPoints->insert(*p_iter);
403:           }
404:         }
405:         // FIX: We would like to avoid the copy here with cone()
406:         if(this->_modifiedPoints->size() > 0) {
407:           this->computeHeight(height, sieve, sieve->cone(this->_modifiedPoints), maxHeight);
408:         }
409:       };
410:       void computeHeights() {
411:         const std::string name("height");

413:         this->_maxHeight = -1;
414:         for(typename sheaf_type::iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
415:           const Obj<patch_label_type>& label = this->createLabel(s_iter->first, name);

417:           this->_maxHeights[s_iter->first] = -1;
418:           this->computeHeight(label, s_iter->second, s_iter->second->leaves(), this->_maxHeights[s_iter->first]);
419:           if (this->_maxHeights[s_iter->first] > this->_maxHeight) this->_maxHeight = this->_maxHeights[s_iter->first];
420:         }
421:       };
422:       int height() const {return this->_maxHeight;};
423:       int height(const patch_type& patch) {
424:         this->checkPatch(patch);
425:         return this->_maxHeights[patch];
426:       };
427:       int height(const patch_type& patch, const point_type& point) {
428:         return this->getValue(this->_labels["height"][patch], point, -1);
429:       };
430:       const Obj<label_sequence>& heightStratum(const patch_type& patch, int height) {
431:         return this->getLabelStratum(patch, "height", height);
432:       };
433:       template<class InputPoints>
434:       void computeDepth(const Obj<patch_label_type>& depth, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxDepth) {
435:         this->_modifiedPoints->clear();

437:         for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
438:           // Compute the max depth of the points in the cone of p, and add 1
439:           int d0 = this->getValue(depth, *p_iter, -1);
440:           int d1 = this->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;

442:           if(d1 != d0) {
443:             this->setValue(depth, *p_iter, d1);
444:             if (d1 > maxDepth) maxDepth = d1;
445:             this->_modifiedPoints->insert(*p_iter);
446:           }
447:         }
448:         // FIX: We would like to avoid the copy here with support()
449:         if(this->_modifiedPoints->size() > 0) {
450:           this->computeDepth(depth, sieve, sieve->support(this->_modifiedPoints), maxDepth);
451:         }
452:       };
453:       void computeDepths() {
454:         const std::string name("depth");

456:         this->_maxDepth = -1;
457:         for(typename sheaf_type::iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
458:           const Obj<patch_label_type>& label = this->createLabel(s_iter->first, name);

460:           this->_maxDepths[s_iter->first] = -1;
461:           this->computeDepth(label, s_iter->second, s_iter->second->roots(), this->_maxDepths[s_iter->first]);
462:           if (this->_maxDepths[s_iter->first] > this->_maxDepth) this->_maxDepth = this->_maxDepths[s_iter->first];
463:         }
464:       };
465:       int depth() const {return this->_maxDepth;};
466:       int depth(const patch_type& patch) {
467:         this->checkPatch(patch);
468:         return this->_maxDepths[patch];
469:       };
470:       int depth(const patch_type& patch, const point_type& point) {
471:         return this->getValue(this->_labels["depth"][patch], point, -1);
472:       };
473:       const Obj<label_sequence>& depthStratum(const patch_type& patch, int depth) {
474:         return this->getLabelStratum(patch, "depth", depth);
475:       };
478:       void stratify() {
479:         ALE_LOG_EVENT_BEGIN;
480:         this->computeHeights();
481:         this->computeDepths();
482:         ALE_LOG_EVENT_END;
483:       };
484:     public: // Viewers
485:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
486:         if (comm == MPI_COMM_NULL) {
487:           comm = this->comm();
488:         }
489:         if (name == "") {
490:           PetscPrintf(comm, "viewing a Topology\n");
491:         } else {
492:           PetscPrintf(comm, "viewing Topology '%s'\n", name.c_str());
493:         }
494:         PetscPrintf(comm, "  maximum height %d maximum depth %d\n", this->height(), this->depth());
495:         for(typename sheaf_type::const_iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
496:           ostringstream txt;

498:           txt << "Patch " << s_iter->first;
499:           s_iter->second->view(txt.str().c_str(), comm);
500:           PetscPrintf(comm, "  maximum height %d maximum depth %d\n", this->height(s_iter->first), this->depth(s_iter->first));
501:         }
502:         for(typename labels_type::const_iterator l_iter = this->_labels.begin(); l_iter != this->_labels.end(); ++l_iter) {
503:           PetscPrintf(comm, "  label %s constructed\n", l_iter->first.c_str());
504:         }
505:       };
506:     public:
507:       void constructOverlap(const patch_type& patch) {
508:         if (this->_calculatedOverlap) return;
509:         this->constructOverlap(this->getPatch(patch)->base(), this->_sendOverlap, this->_recvOverlap);
510:         this->constructOverlap(this->getPatch(patch)->cap(), this->_sendOverlap, this->_recvOverlap);
511:         if (this->debug()) {
512:           this->_sendOverlap->view("Send overlap");
513:           this->_recvOverlap->view("Receive overlap");
514:         }
515:         this->_calculatedOverlap = true;
516:       };
517:       template<typename Sequence>
518:       void constructOverlap(const Obj<Sequence>& points, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
519:         point_type *sendBuf = new point_type[points->size()];
520:         int         size    = 0;
521:         for(typename Sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
522:           sendBuf[size++] = *l_iter;
523:         }
524:         int *sizes   = new int[this->commSize()];   // The number of points coming from each process
525:         int *offsets = new int[this->commSize()+1]; // Prefix sums for sizes
526:         int *oldOffs = new int[this->commSize()+1]; // Temporary storage
527:         point_type *remotePoints = NULL;            // The points from each process
528:         int        *remoteRanks  = NULL;            // The rank and number of overlap points of each process that overlaps another

530:         // Change to Allgather() for the correct binning algorithm
531:         MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, this->comm());
532:         if (this->commRank() == 0) {
533:           offsets[0] = 0;
534:           for(int p = 1; p <= this->commSize(); p++) {
535:             offsets[p] = offsets[p-1] + sizes[p-1];
536:           }
537:           remotePoints = new point_type[offsets[this->commSize()]];
538:         }
539:         MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, this->comm());
540:         std::map<int, std::map<int, std::set<point_type> > > overlapInfo; // Maps (p,q) to their set of overlap points

542:         if (this->commRank() == 0) {
543:           for(int p = 0; p < this->commSize(); p++) {
544:             std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]]);
545:           }
546:           for(int p = 0; p <= this->commSize(); p++) {
547:             oldOffs[p] = offsets[p];
548:           }
549:           for(int p = 0; p < this->commSize(); p++) {
550:             for(int q = p+1; q < this->commSize(); q++) {
551:               std::set_intersection(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
552:                                     &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
553:                                     std::insert_iterator<std::set<point_type> >(overlapInfo[p][q], overlapInfo[p][q].begin()));
554:               overlapInfo[q][p] = overlapInfo[p][q];
555:             }
556:             sizes[p]     = overlapInfo[p].size()*2;
557:             offsets[p+1] = offsets[p] + sizes[p];
558:           }
559:           remoteRanks = new int[offsets[this->commSize()]];
560:           int       k = 0;
561:           for(int p = 0; p < this->commSize(); p++) {
562:             for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
563:               remoteRanks[k*2]   = r_iter->first;
564:               remoteRanks[k*2+1] = r_iter->second.size();
565:               k++;
566:             }
567:           }
568:         }
569:         int numOverlaps;                          // The number of processes overlapping this process
570:         MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, this->comm());
571:         int *overlapRanks = new int[numOverlaps]; // The rank and overlap size for each overlapping process
572:         MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, this->comm());
573:         point_type *sendPoints = NULL;            // The points to send to each process
574:         if (this->commRank() == 0) {
575:           for(int p = 0, k = 0; p < this->commSize(); p++) {
576:             sizes[p] = 0;
577:             for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
578:               sizes[p] += remoteRanks[k*2+1];
579:               k++;
580:             }
581:             offsets[p+1] = offsets[p] + sizes[p];
582:           }
583:           sendPoints = new point_type[offsets[this->commSize()]];
584:           for(int p = 0, k = 0; p < this->commSize(); p++) {
585:             for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
586:               int rank = r_iter->first;
587:               for(typename std::set<point_type>::iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
588:                 sendPoints[k++] = *p_iter;
589:               }
590:             }
591:           }
592:         }
593:         int numOverlapPoints = 0;
594:         for(int r = 0; r < numOverlaps/2; r++) {
595:           numOverlapPoints += overlapRanks[r*2+1];
596:         }
597:         point_type *overlapPoints = new point_type[numOverlapPoints];
598:         MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints, MPI_INT, 0, this->comm());

600:         for(int r = 0, k = 0; r < numOverlaps/2; r++) {
601:           int rank = overlapRanks[r*2];

603:           for(int p = 0; p < overlapRanks[r*2+1]; p++) {
604:             point_type point = overlapPoints[k++];

606:             sendOverlap->addArrow(point, rank, point);
607:             recvOverlap->addArrow(rank, point, point);
608:           }
609:         }

611:         delete [] overlapPoints;
612:         delete [] overlapRanks;
613:         delete [] sizes;
614:         delete [] offsets;
615:         delete [] oldOffs;
616:         if (this->commRank() == 0) {
617:           delete [] remoteRanks;
618:           delete [] remotePoints;
619:           delete [] sendPoints;
620:         }
621:       };
622:     };

624:     // An AbstractSection is a mapping from Sieve points to sets of values
625:     //   This is our presentation of a section of a fibre bundle,
626:     //     in which the Topology is the base space, and
627:     //     the value sets are vectors in the fiber spaces
628:     //   The main interface to values is through restrict() and update()
629:     //     This retrieval uses Sieve closure()
630:     //     We should call these rawRestrict/rawUpdate
631:     //   The Section must also be able to set/report the dimension of each fiber
632:     //     for which we use another section we call an \emph{Atlas}
633:     //     For some storage schemes, we also want offsets to go with these dimensions
634:     //   We must have some way of limiting the points associated with values
635:     //     so each section must support a getPatch() call returning the points with associated fibers
636:     //     I was using the Chart for this
637:     //   The Section must be able to participate in \emph{completion}
638:     //     This means restrict to a provided overlap, and exchange in the restricted sections
639:     //     Completion does not use hierarchy, so we see the Topology as a DiscreteTopology

641:     // A ConstantSection is the simplest Section
642:     //   All fibers are dimension 1
643:     //   All values are equal to a constant
644:     //     We need no value storage and no communication for completion
645:     template<typename Topology_, typename Value_>
646:     class NewConstantSection : public ALE::ParallelObject {
647:     public:
648:       typedef Topology_                          topology_type;
649:       typedef typename topology_type::patch_type patch_type;
650:       typedef typename topology_type::sieve_type sieve_type;
651:       typedef typename topology_type::point_type point_type;
652:       typedef std::set<point_type>               chart_type;
653:       typedef std::map<patch_type, chart_type>   atlas_type;
654:       typedef Value_                             value_type;
655:     protected:
656:       Obj<topology_type> _topology;
657:       atlas_type         _atlas;
658:       chart_type         _emptyChart;
659:       value_type         _value;
660:       value_type         _defaultValue;
661:     public:
662:       NewConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _defaultValue(0) {
663:         this->_topology = new topology_type(comm, debug);
664:       };
665:       NewConstantSection(const Obj<topology_type>& topology) : ParallelObject(topology->comm(), topology->debug()), _topology(topology) {};
666:       NewConstantSection(const Obj<topology_type>& topology, const value_type& value) : ParallelObject(topology->comm(), topology->debug()), _topology(topology), _value(value), _defaultValue(value) {};
667:       NewConstantSection(const Obj<topology_type>& topology, const value_type& value, const value_type& defaultValue) : ParallelObject(topology->comm(), topology->debug()), _topology(topology), _value(value), _defaultValue(defaultValue) {};
668:     public: // Verifiers
669:       void checkPatch(const patch_type& patch) const {
670:         this->_topology->checkPatch(patch);
671:         if (this->_atlas.find(patch) == this->_atlas.end()) {
672:           ostringstream msg;
673:           msg << "Invalid atlas patch " << patch << std::endl;
674:           throw ALE::Exception(msg.str().c_str());
675:         }
676:       };
677:       void checkPoint(const patch_type& patch, const point_type& point) const {
678:         this->checkPatch(patch);
679:         if (this->_atlas.find(patch)->second.find(point) == this->_atlas.find(patch)->second.end()) {
680:           ostringstream msg;
681:           msg << "Invalid section point " << point << std::endl;
682:           throw ALE::Exception(msg.str().c_str());
683:         }
684:       };
685:       void checkDimension(const int& dim) {
686:         if (dim != 1) {
687:           ostringstream msg;
688:           msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
689:           throw ALE::Exception(msg.str().c_str());
690:         }
691:       };
692:       bool hasPatch(const patch_type& patch) {
693:         if (this->_atlas.find(patch) == this->_atlas.end()) {
694:           return false;
695:         }
696:         return true;
697:       };
698:       bool hasPoint(const patch_type& patch, const point_type& point) const {
699:         this->checkPatch(patch);
700:         return this->_atlas.find(patch)->second.count(point) > 0;
701:       };
702:     public: // Accessors
703:       const Obj<topology_type>& getTopology() const {return this->_topology;};
704:       void setTopology(const Obj<topology_type>& topology) {this->_topology = topology;};
705:       const chart_type& getPatch(const patch_type& patch) {
706:         if (this->hasPatch(patch)) {
707:           return this->_atlas[patch];
708:         }
709:         return this->_emptyChart;
710:       };
711:       void updatePatch(const patch_type& patch, const point_type& point) {
712:         this->_atlas[patch].insert(point);
713:       };
714:       template<typename Points>
715:       void updatePatch(const patch_type& patch, const Obj<Points>& points) {
716:         this->_atlas[patch].insert(points->begin(), points->end());
717:       };
718:       value_type getDefaultValue() {return this->_defaultValue;};
719:       void setDefaultValue(const value_type value) {this->_defaultValue = value;};
720:     public: // Sizes
721:       void clear() {
722:         this->_atlas.clear();
723:       };
724:       int getFiberDimension(const patch_type& patch, const point_type& p) const {
725:         if (this->hasPoint(patch, p)) return 1;
726:         return 0;
727:       };
728:       void setFiberDimension(const patch_type& patch, const point_type& p, int dim) {
729:         this->checkDimension(dim);
730:         this->updatePatch(patch, p);
731:       };
732:       template<typename Sequence>
733:       void setFiberDimension(const patch_type& patch, const Obj<Sequence>& points, int dim) {
734:         for(typename topology_type::label_sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
735:           this->setFiberDimension(patch, *p_iter, dim);
736:         }
737:       };
738:       void addFiberDimension(const patch_type& patch, const point_type& p, int dim) {
739:         if (this->hasPatch(patch) && (this->_atlas[patch].find(p) != this->_atlas[patch].end())) {
740:           ostringstream msg;
741:           msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
742:           throw ALE::Exception(msg.str().c_str());
743:         } else {
744:           this->setFiberDimension(patch, p, dim);
745:         }
746:       };
747:       void setFiberDimensionByDepth(const patch_type& patch, int depth, int dim) {
748:         this->setFiberDimension(patch, this->_topology->getLabelStratum(patch, "depth", depth), dim);
749:       };
750:       void setFiberDimensionByHeight(const patch_type& patch, int height, int dim) {
751:         this->setFiberDimension(patch, this->_topology->getLabelStratum(patch, "height", height), dim);
752:       };
753:       int size(const patch_type& patch) {return this->_atlas[patch].size();};
754:       int size(const patch_type& patch, const point_type& p) {return this->getFiberDimension(patch, p);};
755:     public: // Restriction
756:       const value_type *restrict(const patch_type& patch, const point_type& p) const {
757:         //std::cout <<"["<<this->commRank()<<"]: Constant restrict ("<<patch<<","<<p<<") from " << std::endl;
758:         //for(typename chart_type::iterator c_iter = this->_atlas.find(patch)->second.begin(); c_iter != this->_atlas.find(patch)->second.end(); ++c_iter) {
759:         //  std::cout <<"["<<this->commRank()<<"]:   point " << *c_iter << std::endl;
760:         //}
761:         if (this->hasPoint(patch, p)) {
762:           return &this->_value;
763:         }
764:         return &this->_defaultValue;
765:       };
766:       const value_type *restrictPoint(const patch_type& patch, const point_type& p) const {return this->restrict(patch, p);};
767:       void update(const patch_type& patch, const point_type& p, const value_type v[]) {
768:         this->checkPatch(patch);
769:         this->_value = v[0];
770:       };
771:       void updatePoint(const patch_type& patch, const point_type& p, const value_type v[]) {return this->update(patch, p, v);};
772:       void updateAdd(const patch_type& patch, const point_type& p, const value_type v[]) {
773:         this->checkPatch(patch);
774:         this->_value += v[0];
775:       };
776:       void updateAddPoint(const patch_type& patch, const point_type& p, const value_type v[]) {return this->updateAdd(patch, p, v);};
777:     public:
778:       void copy(const Obj<NewConstantSection>& section) {
779:         const typename topology_type::sheaf_type& patches = this->_topology->getPatches();

781:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
782:           const patch_type& patch = p_iter->first;
783:           if (!section->hasPatch(patch)) continue;
784:           const chart_type& chart = section->getPatch(patch);

786:           for(typename chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
787:             this->updatePatch(patch, *c_iter);
788:           }
789:           this->_value = section->restrict(patch, *chart.begin())[0];
790:         }
791:       };
792:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
793:         ostringstream txt;
794:         int rank;

796:         if (comm == MPI_COMM_NULL) {
797:           comm = this->comm();
798:           rank = this->commRank();
799:         } else {
800:           MPI_Comm_rank(comm, &rank);
801:         }
802:         if (name == "") {
803:           if(rank == 0) {
804:             txt << "viewing a ConstantSection" << std::endl;
805:           }
806:         } else {
807:           if(rank == 0) {
808:             txt << "viewing ConstantSection '" << name << "'" << std::endl;
809:           }
810:         }
811:         const typename topology_type::sheaf_type& sheaf = this->_topology->getPatches();

813:         for(typename topology_type::sheaf_type::const_iterator p_iter = sheaf.begin(); p_iter != sheaf.end(); ++p_iter) {
814:           txt <<"["<<this->commRank()<<"]: Patch " << p_iter->first << std::endl;
815:           txt <<"["<<this->commRank()<<"]:   Value " << this->_value << std::endl;
816:         }
817:         PetscSynchronizedPrintf(comm, txt.str().c_str());
818:         PetscSynchronizedFlush(comm);
819:       };
820:     };

822:     // A UniformSection often acts as an Atlas
823:     //   All fibers are the same dimension
824:     //     Note we can use a ConstantSection for this Atlas
825:     //   Each point may have a different vector
826:     //     Thus we need storage for values, and hence must implement completion
827:     template<typename Topology_, typename Value_, int fiberDim = 1>
828:     class UniformSection : public ALE::ParallelObject {
829:     public:
830:       typedef Topology_                              topology_type;
831:       typedef typename topology_type::patch_type     patch_type;
832:       typedef typename topology_type::sieve_type     sieve_type;
833:       typedef typename topology_type::point_type     point_type;
834:       typedef NewConstantSection<topology_type, int> atlas_type;
835:       typedef typename atlas_type::chart_type        chart_type;
836:       typedef Value_                                 value_type;
837:       typedef struct {value_type v[fiberDim];}       fiber_type;
838:       typedef std::map<point_type, fiber_type>       array_type;
839:       typedef std::map<patch_type, array_type>       values_type;
840:     protected:
841:       Obj<atlas_type> _atlas;
842:       values_type     _arrays;
843:     public:
844:       UniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
845:         this->_atlas = new atlas_type(comm, debug);
846:       };
847:       UniformSection(const Obj<topology_type>& topology) : ParallelObject(topology->comm(), topology->debug()) {
848:         this->_atlas = new atlas_type(topology, fiberDim, 0);
849:       };
850:       UniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas) {};
851:     protected:
852:       value_type *getRawArray(const int size) {
853:         static value_type *array   = NULL;
854:         static int         maxSize = 0;

856:         if (size > maxSize) {
857:           maxSize = size;
858:           if (array) delete [] array;
859:           array = new value_type[maxSize];
860:         };
861:         return array;
862:       };
863:     public: // Verifiers
864:       void checkPatch(const patch_type& patch) {
865:         this->_atlas->checkPatch(patch);
866:         if (this->_arrays.find(patch) == this->_arrays.end()) {
867:           ostringstream msg;
868:           msg << "Invalid section patch: " << patch << std::endl;
869:           throw ALE::Exception(msg.str().c_str());
870:         }
871:       };
872:       bool hasPatch(const patch_type& patch) {
873:         return this->_atlas->hasPatch(patch);
874:       };
875:       bool hasPoint(const patch_type& patch, const point_type& point) {
876:         return this->_atlas->hasPoint(patch, point);
877:       };
878:       void checkDimension(const int& dim) {
879:         if (dim != fiberDim) {
880:           ostringstream msg;
881:           msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
882:           throw ALE::Exception(msg.str().c_str());
883:         }
884:       };
885:     public: // Accessors
886:       const Obj<atlas_type>& getAtlas() {return this->_atlas;};
887:       void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
888:       const Obj<topology_type>& getTopology() {return this->_atlas->getTopology();};
889:       void setTopology(const Obj<topology_type>& topology) {this->_atlas->setTopology(topology);};
890:       const chart_type& getPatch(const patch_type& patch) {
891:         return this->_atlas->getPatch(patch);
892:       };
893:       void updatePatch(const patch_type& patch, const point_type& point) {
894:         this->setFiberDimension(patch, point, 1);
895:       };
896:       template<typename Points>
897:       void updatePatch(const patch_type& patch, const Obj<Points>& points) {
898:         for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
899:           this->setFiberDimension(patch, *p_iter, 1);
900:         }
901:       };
902:       void copy(const Obj<UniformSection<Topology_, Value_, fiberDim> >& section) {
903:         this->getAtlas()->copy(section->getAtlas());
904:         const typename topology_type::sheaf_type& sheaf = section->getTopology()->getPatches();

906:         for(typename topology_type::sheaf_type::const_iterator s_iter = sheaf.begin(); s_iter != sheaf.end(); ++s_iter) {
907:           const patch_type& patch = s_iter->first;
908:           if (!section->hasPatch(patch)) continue;
909:           const chart_type& chart = section->getPatch(patch);

911:           for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
912:             this->update(s_iter->first, *c_iter, section->restrict(s_iter->first, *c_iter));
913:           }
914:         }
915:       };
916:     public: // Sizes
917:       void clear() {
918:         this->_atlas->clear();
919:         this->_arrays.clear();
920:       };
921:       int getFiberDimension(const patch_type& patch, const point_type& p) const {
922:         // Could check for non-existence here
923:         return this->_atlas->restrictPoint(patch, p)[0];
924:       };
925:       void setFiberDimension(const patch_type& patch, const point_type& p, int dim) {
926:         this->checkDimension(dim);
927:         this->_atlas->updatePatch(patch, p);
928:         this->_atlas->updatePoint(patch, p, &dim);
929:       };
930:       template<typename Sequence>
931:       void setFiberDimension(const patch_type& patch, const Obj<Sequence>& points, int dim) {
932:         for(typename topology_type::label_sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
933:           this->setFiberDimension(patch, *p_iter, dim);
934:         }
935:       };
936:       void addFiberDimension(const patch_type& patch, const point_type& p, int dim) {
937:         if (this->hasPatch(patch) && (this->_atlas[patch].find(p) != this->_atlas[patch].end())) {
938:           ostringstream msg;
939:           msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
940:           throw ALE::Exception(msg.str().c_str());
941:         } else {
942:           this->setFiberDimension(patch, p, dim);
943:         }
944:       };
945:       void setFiberDimensionByDepth(const patch_type& patch, int depth, int dim) {
946:         this->setFiberDimension(patch, this->getTopology()->getLabelStratum(patch, "depth", depth), dim);
947:       };
948:       void setFiberDimensionByHeight(const patch_type& patch, int height, int dim) {
949:         this->setFiberDimension(patch, this->getTopology()->getLabelStratum(patch, "height", height), dim);
950:       };
951:       int size(const patch_type& patch) {
952:         const typename atlas_type::chart_type& points = this->_atlas->getPatch(patch);
953:         int size = 0;

955:         for(typename atlas_type::chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
956:           size += this->getFiberDimension(patch, *p_iter);
957:         }
958:         return size;
959:       };
960:       int size(const patch_type& patch, const point_type& p) {
961:         const typename atlas_type::chart_type&  points  = this->_atlas->getPatch(patch);
962:         const Obj<typename sieve_type::coneSet> closure = this->getTopology()->getPatch(patch)->closure(p);
963:         typename sieve_type::coneSet::iterator  end     = closure->end();
964:         int size = 0;

966:         for(typename sieve_type::coneSet::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
967:           if (points.count(*c_iter)) {
968:             size += this->getFiberDimension(patch, *c_iter);
969:           }
970:         }
971:         return size;
972:       };
973:       void orderPatches() {};
974:     public: // Restriction
975:       const array_type& restrict(const patch_type& patch) {
976:         this->checkPatch(patch);
977:         return this->_arrays[patch];
978:       };
979:       // Return the values for the closure of this point
980:       //   use a smart pointer?
981:       const value_type *restrict(const patch_type& patch, const point_type& p) {
982:         this->checkPatch(patch);
983:         const chart_type& chart = this->getPatch(patch);
984:         array_type& array  = this->_arrays[patch];
985:         const int   size   = this->size(patch, p);
986:         value_type *values = this->getRawArray(size);
987:         int         j      = -1;

989:         // We could actually ask for the height of the individual point
990:         if (this->getTopology()->height(patch) < 2) {
991:           // Only avoids the copy of closure()
992:           const int& dim = this->_atlas->restrict(patch, p)[0];

994:           if (chart.count(p)) {
995:             for(int i = 0; i < dim; ++i) {
996:               values[++j] = array[p].v[i];
997:             }
998:           }
999:           // Should be closure()
1000:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1001:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1003:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1004:             if (chart.count(*p_iter)) {
1005:               const int& dim = this->_atlas->restrict(patch, *p_iter)[0];

1007:               for(int i = 0; i < dim; ++i) {
1008:                 values[++j] = array[*p_iter].v[i];
1009:               }
1010:             }
1011:           }
1012:         } else {
1013:           throw ALE::Exception("Not yet implemented for interpolated sieves");
1014:         }
1015:         if (j != size-1) {
1016:           ostringstream txt;

1018:           txt << "Invalid restrict to point " << p << std::endl;
1019:           txt << "  j " << j << " should be " << (size-1) << std::endl;
1020:           std::cout << txt.str();
1021:           throw ALE::Exception(txt.str().c_str());
1022:         }
1023:         return values;
1024:       };
1025:       void update(const patch_type& patch, const point_type& p, const value_type v[]) {
1026:         this->_atlas->checkPatch(patch);
1027:         const chart_type& chart = this->getPatch(patch);
1028:         array_type& array = this->_arrays[patch];
1029:         int         j     = -1;

1031:         if (this->getTopology()->height(patch) < 2) {
1032:           // Only avoids the copy of closure()
1033:           const int& dim = this->_atlas->restrict(patch, p)[0];

1035:           if (chart.count(p)) {
1036:             for(int i = 0; i < dim; ++i) {
1037:               array[p].v[i] = v[++j];
1038:             }
1039:           }
1040:           // Should be closure()
1041:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);

1043:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
1044:             if (chart.count(*p_iter)) {
1045:               const int& dim = this->_atlas->restrict(patch, *p_iter)[0];

1047:               for(int i = 0; i < dim; ++i) {
1048:                 array[*p_iter].v[i] = v[++j];
1049:               }
1050:             }
1051:           }
1052:         } else {
1053:           throw ALE::Exception("Not yet implemented for interpolated sieves");
1054:         }
1055:       };
1056:       void updateAdd(const patch_type& patch, const point_type& p, const value_type v[]) {
1057:         this->_atlas->checkPatch(patch);
1058:         const chart_type& chart = this->getPatch(patch);
1059:         array_type& array = this->_arrays[patch];
1060:         int         j     = -1;

1062:         if (this->getTopology()->height(patch) < 2) {
1063:           // Only avoids the copy of closure()
1064:           const int& dim = this->_atlas->restrict(patch, p)[0];

1066:           if (chart.count(p)) {
1067:             for(int i = 0; i < dim; ++i) {
1068:               array[p].v[i] += v[++j];
1069:             }
1070:           }
1071:           // Should be closure()
1072:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);

1074:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
1075:             if (chart.count(*p_iter)) {
1076:               const int& dim = this->_atlas->restrict(patch, *p_iter)[0];

1078:               for(int i = 0; i < dim; ++i) {
1079:                 array[*p_iter].v[i] += v[++j];
1080:               }
1081:             }
1082:           }
1083:         } else {
1084:           throw ALE::Exception("Not yet implemented for interpolated sieves");
1085:         }
1086:       };
1087:       // Return only the values associated to this point, not its closure
1088:       const value_type *restrictPoint(const patch_type& patch, const point_type& p) {
1089:         this->checkPatch(patch);
1090:         return this->_arrays[patch][p].v;
1091:       };
1092:       // Update only the values associated to this point, not its closure
1093:       void updatePoint(const patch_type& patch, const point_type& p, const value_type v[]) {
1094:         this->_atlas->checkPatch(patch);
1095:         for(int i = 0; i < fiberDim; ++i) {
1096:           this->_arrays[patch][p].v[i] = v[i];
1097:         }
1098:       };
1099:       // Update only the values associated to this point, not its closure
1100:       void updateAddPoint(const patch_type& patch, const point_type& p, const value_type v[]) {
1101:         this->_atlas->checkPatch(patch);
1102:         for(int i = 0; i < fiberDim; ++i) {
1103:           this->_arrays[patch][p].v[i] += v[i];
1104:         }
1105:       };
1106:     public:
1107:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1108:         ostringstream txt;
1109:         int rank;

1111:         if (comm == MPI_COMM_NULL) {
1112:           comm = this->comm();
1113:           rank = this->commRank();
1114:         } else {
1115:           MPI_Comm_rank(comm, &rank);
1116:         }
1117:         if (name == "") {
1118:           if(rank == 0) {
1119:             txt << "viewing a UniformSection" << std::endl;
1120:           }
1121:         } else {
1122:           if(rank == 0) {
1123:             txt << "viewing UniformSection '" << name << "'" << std::endl;
1124:           }
1125:         }
1126:         for(typename values_type::const_iterator a_iter = this->_arrays.begin(); a_iter != this->_arrays.end(); ++a_iter) {
1127:           const patch_type& patch = a_iter->first;
1128:           array_type&       array = this->_arrays[patch];

1130:           txt << "[" << this->commRank() << "]: Patch " << patch << std::endl;
1131:           const typename atlas_type::chart_type& chart = this->_atlas->getPatch(patch);

1133:           for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1134:             const point_type&                     point = *p_iter;
1135:             const typename atlas_type::value_type dim   = this->_atlas->restrict(patch, point)[0];

1137:             if (dim != 0) {
1138:               txt << "[" << this->commRank() << "]:   " << point << " dim " << dim << "  ";
1139:               for(int i = 0; i < dim; i++) {
1140:                 txt << " " << array[point].v[i];
1141:               }
1142:               txt << std::endl;
1143:             }
1144:           }
1145:         }
1146:         PetscSynchronizedPrintf(comm, txt.str().c_str());
1147:         PetscSynchronizedFlush(comm);
1148:       };
1149:     };

1151:     // A Section is our most general construct (more general ones could be envisioned)
1152:     //   The Atlas is a UniformSection of dimension 1 and value type Point
1153:     //     to hold each fiber dimension and offsets into a contiguous patch array
1154:     template<typename Topology_, typename Value_>
1155:     class Section : public ALE::ParallelObject {
1156:     public:
1157:       typedef Topology_                                 topology_type;
1158:       typedef typename topology_type::patch_type        patch_type;
1159:       typedef typename topology_type::sieve_type        sieve_type;
1160:       typedef typename topology_type::point_type        point_type;
1161:       typedef ALE::Point                                index_type;
1162:       typedef UniformSection<topology_type, index_type> atlas_type;
1163:       typedef typename atlas_type::chart_type           chart_type;
1164:       typedef Value_                                    value_type;
1165:       typedef value_type *                              array_type;
1166:       typedef std::map<patch_type, array_type>          values_type;
1167:       typedef std::vector<index_type>                   IndexArray;
1168:     protected:
1169:       Obj<atlas_type> _atlas;
1170:       Obj<atlas_type> _atlasNew;
1171:       values_type     _arrays;
1172:       Obj<IndexArray> _indexArray;
1173:     public:
1174:       Section(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
1175:         this->_atlas      = new atlas_type(comm, debug);
1176:         this->_atlasNew   = NULL;
1177:         this->_indexArray = new IndexArray();
1178:       };
1179:       Section(const Obj<topology_type>& topology) : ParallelObject(topology->comm(), topology->debug()), _atlasNew(NULL) {
1180:         this->_atlas      = new atlas_type(topology);
1181:         this->_indexArray = new IndexArray();
1182:       };
1183:       Section(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL) {
1184:         this->_indexArray = new IndexArray();
1185:       };
1186:       virtual ~Section() {
1187:         for(typename values_type::iterator a_iter = this->_arrays.begin(); a_iter != this->_arrays.end(); ++a_iter) {
1188:           delete [] a_iter->second;
1189:           a_iter->second = NULL;
1190:         }
1191:       };
1192:     protected:
1193:       value_type *getRawArray(const int size) {
1194:         static value_type *array   = NULL;
1195:         static int         maxSize = 0;

1197:         if (size > maxSize) {
1198:           maxSize = size;
1199:           if (array) delete [] array;
1200:           array = new value_type[maxSize];
1201:         };
1202:         return array;
1203:       };
1204:     public: // Verifiers
1205:       void checkPatch(const patch_type& patch) {
1206:         this->_atlas->checkPatch(patch);
1207:         if (this->_arrays.find(patch) == this->_arrays.end()) {
1208:           ostringstream msg;
1209:           msg << "Invalid section patch: " << patch << std::endl;
1210:           throw ALE::Exception(msg.str().c_str());
1211:         }
1212:       };
1213:       bool hasPatch(const patch_type& patch) {
1214:         return this->_atlas->hasPatch(patch);
1215:       };
1216:     public: // Accessors
1217:       const Obj<atlas_type>& getAtlas() {return this->_atlas;};
1218:       void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1219:       const Obj<topology_type>& getTopology() {return this->_atlas->getTopology();};
1220:       void setTopology(const Obj<topology_type>& topology) {this->_atlas->setTopology(topology);};
1221:       const chart_type& getPatch(const patch_type& patch) {
1222:         return this->_atlas->getPatch(patch);
1223:       };
1224:       bool hasPoint(const patch_type& patch, const point_type& point) {
1225:         return this->_atlas->hasPoint(patch, point);
1226:       };
1227:     public: // Sizes
1228:       void clear() {
1229:         this->_atlas->clear();
1230:         this->_arrays.clear();
1231:       };
1232:       int getFiberDimension(const patch_type& patch, const point_type& p) const {
1233:         // Could check for non-existence here
1234:         return this->_atlas->restrictPoint(patch, p)->prefix;
1235:       };
1236:       int getFiberDimension(const Obj<atlas_type>& atlas, const patch_type& patch, const point_type& p) const {
1237:         // Could check for non-existence here
1238:         return atlas->restrictPoint(patch, p)->prefix;
1239:       };
1240:       void setFiberDimension(const patch_type& patch, const point_type& p, int dim) {
1241:         const index_type idx(dim, -1);
1242:         this->_atlas->updatePatch(patch, p);
1243:         this->_atlas->updatePoint(patch, p, &idx);
1244:       };
1245:       template<typename Sequence>
1246:       void setFiberDimension(const patch_type& patch, const Obj<Sequence>& points, int dim) {
1247:         for(typename topology_type::label_sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1248:           this->setFiberDimension(patch, *p_iter, dim);
1249:         }
1250:       };
1251:       void addFiberDimension(const patch_type& patch, const point_type& p, int dim) {
1252:         if (this->_atlas->hasPatch(patch) && this->_atlas->hasPoint(patch, p)) {
1253:           const index_type values(dim, 0);
1254:           this->_atlas->updateAddPoint(patch, p, &values);
1255:         } else {
1256:           this->setFiberDimension(patch, p, dim);
1257:         }
1258:       };
1259:       void setFiberDimensionByDepth(const patch_type& patch, int depth, int dim) {
1260:         this->setFiberDimension(patch, this->getTopology()->getLabelStratum(patch, "depth", depth), dim);
1261:       };
1262:       void setFiberDimensionByHeight(const patch_type& patch, int height, int dim) {
1263:         this->setFiberDimension(patch, this->getTopology()->getLabelStratum(patch, "height", height), dim);
1264:       };
1265:       int size(const patch_type& patch) {
1266:         const typename atlas_type::chart_type& points = this->_atlas->getPatch(patch);
1267:         int size = 0;

1269:         for(typename atlas_type::chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1270:           size += std::max(0, this->getFiberDimension(patch, *p_iter));
1271:         }
1272:         return size;
1273:       };
1274:       int sizeWithBC(const patch_type& patch) {
1275:         const typename atlas_type::chart_type& points = this->_atlas->getPatch(patch);
1276:         int size = 0;

1278:         for(typename atlas_type::chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1279:           size += std::abs(this->getFiberDimension(patch, *p_iter));
1280:         }
1281:         return size;
1282:       };
1283:       int size(const patch_type& patch, const point_type& p) {
1284:         const typename atlas_type::chart_type&  points  = this->_atlas->getPatch(patch);
1285:         const Obj<typename sieve_type::coneSet> closure = this->getTopology()->getPatch(patch)->closure(p);
1286:         typename sieve_type::coneSet::iterator  end     = closure->end();
1287:         int size = 0;

1289:         for(typename sieve_type::coneSet::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
1290:           if (points.count(*c_iter)) {
1291:             size += std::max(0, this->getFiberDimension(patch, *c_iter));
1292:           }
1293:         }
1294:         return size;
1295:       };
1296:       int sizeWithBC(const patch_type& patch, const point_type& p) {
1297:         const typename atlas_type::chart_type&  points  = this->_atlas->getPatch(patch);
1298:         const Obj<typename sieve_type::coneSet> closure = this->getTopology()->getPatch(patch)->closure(p);
1299:         typename sieve_type::coneSet::iterator  end     = closure->end();
1300:         int size = 0;

1302:         for(typename sieve_type::coneSet::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
1303:           if (points.count(*c_iter)) {
1304:             size += std::abs(this->getFiberDimension(patch, *c_iter));
1305:           }
1306:         }
1307:         return size;
1308:       };
1309:       int size(const Obj<atlas_type>& atlas, const patch_type& patch) {
1310:         const typename atlas_type::chart_type& points = atlas->getPatch(patch);
1311:         int size = 0;

1313:         for(typename atlas_type::chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1314:           size += std::max(0, this->getFiberDimension(atlas, patch, *p_iter));
1315:         }
1316:         return size;
1317:       };
1318:     public: // Index retrieval
1319:       const index_type& getIndex(const patch_type& patch, const point_type& p) {
1320:         this->checkPatch(patch);
1321:         return this->_atlas->restrictPoint(patch, p)[0];
1322:       };
1323:       template<typename Numbering>
1324:       const index_type getIndex(const patch_type& patch, const point_type& p, const Obj<Numbering>& numbering) {
1325:         this->checkPatch(patch);
1326:         return index_type(this->getFiberDimension(patch, p), numbering->getIndex(p));
1327:       };
1328:       const Obj<IndexArray>& getIndices(const patch_type& patch, const point_type& p, const int level = -1) {
1329:         this->_indexArray->clear();

1331:         if (level == 0) {
1332:           this->_indexArray->push_back(this->getIndex(patch, p));
1333:         } else if ((level == 1) || (this->getTopology()->height(patch) == 1)) {
1334:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1335:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1337:           this->_indexArray->push_back(this->getIndex(patch, p));
1338:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1339:             this->_indexArray->push_back(this->getIndex(patch, *p_iter));
1340:           }
1341:         } else if (level == -1) {
1342:           const Obj<typename sieve_type::coneSet> closure = this->getTopology()->getPatch(patch)->closure(p);
1343:           typename sieve_type::coneSet::iterator  end     = closure->end();

1345:           for(typename sieve_type::coneSet::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
1346:             this->_indexArray->push_back(this->getIndex(patch, *p_iter));
1347:           }
1348:         } else {
1349:           const Obj<typename sieve_type::coneArray> cone = this->getTopology()->getPatch(patch)->nCone(p, level);
1350:           typename sieve_type::coneArray::iterator  end  = cone->end();

1352:           for(typename sieve_type::coneArray::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1353:             this->_indexArray->push_back(this->getIndex(patch, *p_iter));
1354:           }
1355:         }
1356:         return this->_indexArray;
1357:       };
1358:       template<typename Numbering>
1359:       const Obj<IndexArray>& getIndices(const patch_type& patch, const point_type& p, const Obj<Numbering>& numbering, const int level = -1) {
1360:         this->_indexArray->clear();

1362:         if (level == 0) {
1363:           this->_indexArray->push_back(this->getIndex(patch, p, numbering));
1364:         } else if ((level == 1) || (this->getTopology()->height(patch) == 1)) {
1365:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1366:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1368:           this->_indexArray->push_back(this->getIndex(patch, p, numbering));
1369:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1370:             this->_indexArray->push_back(this->getIndex(patch, *p_iter, numbering));
1371:           }
1372:         } else if (level == -1) {
1373:           const Obj<typename sieve_type::coneSet> closure = this->getTopology()->getPatch(patch)->closure(p);
1374:           typename sieve_type::coneSet::iterator  end     = closure->end();

1376:           for(typename sieve_type::coneSet::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
1377:             this->_indexArray->push_back(this->getIndex(patch, *p_iter, numbering));
1378:           }
1379:         } else {
1380:           const Obj<typename sieve_type::coneArray> cone = this->getTopology()->getPatch(patch)->nCone(p, level);
1381:           typename sieve_type::coneArray::iterator  end  = cone->end();

1383:           for(typename sieve_type::coneArray::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1384:             this->_indexArray->push_back(this->getIndex(patch, *p_iter, numbering));
1385:           }
1386:         }
1387:         return this->_indexArray;
1388:       };
1389:     public: // Allocation
1390:       void orderPoint(const Obj<atlas_type>& atlas, const Obj<sieve_type>& sieve, const patch_type& patch, const point_type& point, int& offset, int& bcOffset) {
1391:         const Obj<typename sieve_type::coneSequence>& cone = sieve->cone(point);
1392:         typename sieve_type::coneSequence::iterator   end  = cone->end();
1393:         index_type                                    idx  = atlas->restrictPoint(patch, point)[0];
1394:         const int&                                    dim  = idx.prefix;
1395:         const index_type                              defaultIdx(0, -1);

1397:         if (atlas->getPatch(patch).count(point) == 0) {
1398:           idx = defaultIdx;
1399:         }
1400:         if (idx.index == -1) {
1401:           for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
1402:             if (this->_debug > 1) {std::cout << "    Recursing to " << *c_iter << std::endl;}
1403:             this->orderPoint(atlas, sieve, patch, *c_iter, offset, bcOffset);
1404:           }
1405:           if (dim > 0) {
1406:             if (this->_debug > 1) {std::cout << "  Ordering point " << point << " at " << offset << std::endl;}
1407:             idx.index = offset;
1408:             atlas->updatePoint(patch, point, &idx);
1409:             offset += dim;
1410:           } else if (dim < 0) {
1411:             if (this->_debug > 1) {std::cout << "  Ordering boundary point " << point << " at " << bcOffset << std::endl;}
1412:             idx.index = bcOffset;
1413:             atlas->updatePoint(patch, point, &idx);
1414:             bcOffset += dim;
1415:           }
1416:         }
1417:       }
1418:       void orderPatch(const Obj<atlas_type>& atlas, const patch_type& patch, int& offset, int& bcOffset) {
1419:         const typename atlas_type::chart_type& chart = atlas->getPatch(patch);

1421:         for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1422:           if (this->_debug > 1) {std::cout << "Ordering closure of point " << *p_iter << std::endl;}
1423:           this->orderPoint(atlas, this->getTopology()->getPatch(patch), patch, *p_iter, offset, bcOffset);
1424:         }
1425:         for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1426:           index_type idx = atlas->restrictPoint(patch, *p_iter)[0];
1427:           const int& dim = idx.prefix;

1429:           if (dim < 0) {
1430:             if (this->_debug > 1) {std::cout << "Correcting boundary offset of point " << *p_iter << std::endl;}
1431:             idx.index = offset - (idx.index+2);
1432:             atlas->updatePoint(patch, *p_iter, &idx);
1433:           }
1434:         }
1435:       };
1436:       void orderPatches(const Obj<atlas_type>& atlas) {
1437:         const typename topology_type::sheaf_type& patches = this->getTopology()->getPatches();

1439:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
1440:           if (this->_debug > 1) {std::cout << "Ordering patch " << p_iter->first << std::endl;}
1441:           int offset = 0, bcOffset = -2;

1443:           if (!atlas->hasPatch(p_iter->first)) continue;
1444:           this->orderPatch(atlas, p_iter->first, offset, bcOffset);
1445:         }
1446:       };
1447:       void orderPatches() {
1448:         this->orderPatches(this->_atlas);
1449:       };
1450:       void allocateStorage() {
1451:         const typename topology_type::sheaf_type& patches = this->getTopology()->getPatches();

1453:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
1454:           if (!this->_atlas->hasPatch(p_iter->first)) continue;
1455:           this->_arrays[p_iter->first] = new value_type[this->sizeWithBC(p_iter->first)];
1456:           PetscMemzero(this->_arrays[p_iter->first], this->sizeWithBC(p_iter->first) * sizeof(value_type));
1457:         }
1458:       };
1459:       void allocate() {
1460:         this->orderPatches();
1461:         this->allocateStorage();
1462:       };
1463:       void addPoint(const patch_type& patch, const point_type& point, const int dim) {
1464:         if (dim == 0) return;
1465:         //const typename atlas_type::chart_type& chart = this->_atlas->getPatch(patch);

1467:         //if (chart.find(point) == chart.end()) {
1468:         if (this->_atlasNew.isNull()) {
1469:           this->_atlasNew = new atlas_type(this->getTopology());
1470:           this->_atlasNew->copy(this->_atlas);
1471:         }
1472:         const index_type idx(dim, -1);
1473:         this->_atlasNew->updatePatch(patch, point);
1474:         this->_atlasNew->update(patch, point, &idx);
1475:       };
1476:       void reallocate() {
1477:         if (this->_atlasNew.isNull()) return;
1478:         const typename topology_type::sheaf_type& patches = this->getTopology()->getPatches();

1480:         // Since copy() preserves offsets, we must reinitialize them before ordering
1481:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
1482:           const patch_type&                      patch = p_iter->first;
1483:           const typename atlas_type::chart_type& chart = this->_atlasNew->getPatch(patch);
1484:           index_type                             defaultIdx(0, -1);

1486:           for(typename atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1487:             defaultIdx.prefix = this->_atlasNew->restrict(patch, *c_iter)[0].prefix;
1488:             this->_atlasNew->update(patch, *c_iter, &defaultIdx);
1489:           }
1490:         }
1491:         this->orderPatches(this->_atlasNew);
1492:         // Copy over existing values
1493:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
1494:           const patch_type&                      patch    = p_iter->first;
1495:           value_type                            *newArray = new value_type[this->size(this->_atlasNew, patch)];

1497:           if (!this->_atlas->hasPatch(patch)) {
1498:             this->_arrays[patch] = newArray;
1499:             continue;
1500:           }
1501:           const typename atlas_type::chart_type& chart    = this->_atlas->getPatch(patch);
1502:           const value_type                      *array    = this->_arrays[patch];

1504:           for(typename atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1505:             const index_type& idx       = this->_atlas->restrict(patch, *c_iter)[0];
1506:             const int         size      = idx.prefix;
1507:             const int         offset    = idx.index;
1508:             const int&        newOffset = this->_atlasNew->restrict(patch, *c_iter)[0].index;

1510:             for(int i = 0; i < size; ++i) {
1511:               newArray[newOffset+i] = array[offset+i];
1512:             }
1513:           }
1514:           delete [] this->_arrays[patch];
1515:           this->_arrays[patch] = newArray;
1516:         }
1517:         this->_atlas    = this->_atlasNew;
1518:         this->_atlasNew = NULL;
1519:       };
1520:     public: // Restriction and Update
1521:       // Zero entries
1522:       void zero(const patch_type& patch) {
1523:         this->checkPatch(patch);
1524:         memset(this->_arrays[patch], 0, this->size(patch)* sizeof(value_type));
1525:       };
1526:       // Return a pointer to the entire contiguous storage array
1527:       const value_type *restrict(const patch_type& patch) {
1528:         this->checkPatch(patch);
1529:         return this->_arrays[patch];
1530:       };
1531:       // Update the entire contiguous storage array
1532:       void update(const patch_type& patch, const value_type v[]) {
1533:         const value_type *array = this->_arrays[patch];
1534:         const int         size  = this->size(patch);

1536:         for(int i = 0; i < size; i++) {
1537:           array[i] = v[i];
1538:         }
1539:       };
1540:       // Return the values for the closure of this point
1541:       //   use a smart pointer?
1542:       const value_type *restrict(const patch_type& patch, const point_type& p) {
1543:         this->checkPatch(patch);
1544:         const value_type *a      = this->_arrays[patch];
1545:         const int         size   = this->sizeWithBC(patch, p);
1546:         value_type       *values = this->getRawArray(size);
1547:         int               j      = -1;

1549:         if (this->getTopology()->height(patch) < 2) {
1550:           // Avoids the copy of both
1551:           //   points  in topology->closure()
1552:           //   indices in _atlas->restrict()
1553:           const index_type& pInd = this->_atlas->restrictPoint(patch, p)[0];

1555:           for(int i = pInd.index; i < std::abs(pInd.prefix) + pInd.index; ++i) {
1556:             values[++j] = a[i];
1557:           }
1558:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1559:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1561:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1562:             const index_type& ind    = this->_atlas->restrictPoint(patch, *p_iter)[0];
1563:             const int&        start  = ind.index;
1564:             const int&        length = std::abs(ind.prefix);

1566:             for(int i = start; i < start + length; ++i) {
1567:               values[++j] = a[i];
1568:             }
1569:           }
1570:         } else {
1571:           const Obj<IndexArray>& ind = this->getIndices(patch, p);

1573:           for(typename IndexArray::iterator i_iter = ind->begin(); i_iter != ind->end(); ++i_iter) {
1574:             const int& start  = i_iter->index;
1575:             const int& length = std::abs(i_iter->prefix);

1577:             for(int i = start; i < start + length; ++i) {
1578:               values[++j] = a[i];
1579:             }
1580:           }
1581:         }
1582:         if (j != size-1) {
1583:           ostringstream txt;

1585:           txt << "Invalid restrict to point " << p << std::endl;
1586:           txt << "  j " << j << " should be " << (size-1) << std::endl;
1587:           std::cout << txt.str();
1588:           throw ALE::Exception(txt.str().c_str());
1589:         }
1590:         return values;
1591:       };
1592:       // Update the values for the closure of this point
1593:       void update(const patch_type& patch, const point_type& p, const value_type v[]) {
1594:         this->checkPatch(patch);
1595:         value_type *a = this->_arrays[patch];
1596:         int         j = -1;

1598:         if (this->getTopology()->height(patch) < 2) {
1599:           // Avoids the copy of both
1600:           //   points  in topology->closure()
1601:           //   indices in _atlas->restrict()
1602:           const index_type& pInd = this->_atlas->restrictPoint(patch, p)[0];

1604:           for(int i = pInd.index; i < pInd.prefix + pInd.index; ++i) {
1605:             a[i] = v[++j];
1606:           }
1607:           j += std::max(0, -pInd.prefix);
1608:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1609:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1611:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1612:             const index_type& ind    = this->_atlas->restrictPoint(patch, *p_iter)[0];
1613:             const int&        start  = ind.index;
1614:             const int&        length = ind.prefix;

1616:             for(int i = start; i < start + length; ++i) {
1617:               a[i] = v[++j];
1618:             }
1619:             j += std::max(0, -length);
1620:           }
1621:         } else {
1622:           const Obj<IndexArray>& ind = this->getIndices(patch, p);

1624:           for(typename IndexArray::iterator i_iter = ind->begin(); i_iter != ind->end(); ++i_iter) {
1625:             const int& start  = i_iter->index;
1626:             const int& length = i_iter->prefix;

1628:             for(int i = start; i < start + length; ++i) {
1629:               a[i] = v[++j];
1630:             }
1631:             j += std::max(0, -length);
1632:           }
1633:         }
1634:       };
1635:       // Update the values for the closure of this point
1636:       void updateAdd(const patch_type& patch, const point_type& p, const value_type v[]) {
1637:         this->checkPatch(patch);
1638:         value_type *a = this->_arrays[patch];
1639:         int         j = -1;

1641:         if (this->getTopology()->height(patch) < 2) {
1642:           // Avoids the copy of both
1643:           //   points  in topology->closure()
1644:           //   indices in _atlas->restrict()
1645:           const index_type& pInd = this->_atlas->restrictPoint(patch, p)[0];

1647:           for(int i = pInd.index; i < pInd.prefix + pInd.index; ++i) {
1648:             a[i] += v[++j];
1649:           }
1650:           j += std::max(0, -pInd.prefix);
1651:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1652:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1654:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1655:             const index_type& ind    = this->_atlas->restrictPoint(patch, *p_iter)[0];
1656:             const int&        start  = ind.index;
1657:             const int&        length = ind.prefix;

1659:             for(int i = start; i < start + length; ++i) {
1660:               a[i] += v[++j];
1661:             }
1662:             j += std::max(0, -length);
1663:           }
1664:         } else {
1665:           const Obj<IndexArray>& ind = this->getIndices(patch, p);

1667:           for(typename IndexArray::iterator i_iter = ind->begin(); i_iter != ind->end(); ++i_iter) {
1668:             const int& start  = i_iter->index;
1669:             const int& length = i_iter->prefix;

1671:             for(int i = start; i < start + length; ++i) {
1672:               a[i] += v[++j];
1673:             }
1674:             j += std::max(0, -length);
1675:           }
1676:         }
1677:       };
1678:       // Update the values for the closure of this point
1679:       template<typename Input>
1680:       void update(const patch_type& patch, const point_type& p, const Obj<Input>& v) {
1681:         this->checkPatch(patch);
1682:         value_type *a = this->_arrays[patch];

1684:         if (this->getTopology()->height(patch) == 1) {
1685:           // Only avoids the copy
1686:           const index_type& pInd = this->_atlas->restrictPoint(patch, p)[0];
1687:           typename Input::iterator v_iter = v->begin();
1688:           typename Input::iterator v_end  = v->end();

1690:           for(int i = pInd.index; i < pInd.prefix + pInd.index; ++i) {
1691:             a[i] = *v_iter;
1692:             ++v_iter;
1693:           }
1694:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1695:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1697:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1698:             const index_type& ind    = this->_atlas->restrictPoint(patch, *p_iter)[0];
1699:             const int&        start  = ind.index;
1700:             const int&        length = ind.prefix;

1702:             for(int i = start; i < start + length; ++i) {
1703:               a[i] = *v_iter;
1704:               ++v_iter;
1705:             }
1706:           }
1707:         } else {
1708:           const Obj<IndexArray>& ind = this->getIndices(patch, p);
1709:           typename Input::iterator v_iter = v->begin();
1710:           typename Input::iterator v_end  = v->end();

1712:           for(typename IndexArray::iterator i_iter = ind->begin(); i_iter != ind->end(); ++i_iter) {
1713:             const int& start  = i_iter->index;
1714:             const int& length = i_iter->prefix;

1716:             for(int i = start; i < start + length; ++i) {
1717:               a[i] = *v_iter;
1718:               ++v_iter;
1719:             }
1720:           }
1721:         }
1722:       };
1723:       void updateBC(const patch_type& patch, const point_type& p, const value_type v[]) {
1724:         this->checkPatch(patch);
1725:         value_type *a = this->_arrays[patch];
1726:         int         j = -1;

1728:         if (this->getTopology()->height(patch) < 2) {
1729:           // Avoids the copy of both
1730:           //   points  in topology->closure()
1731:           //   indices in _atlas->restrict()
1732:           const index_type& pInd = this->_atlas->restrictPoint(patch, p)[0];

1734:           for(int i = pInd.index; i < std::abs(pInd.prefix) + pInd.index; ++i) {
1735:             a[i] = v[++j];
1736:           }
1737:           const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
1738:           typename sieve_type::coneSequence::iterator   end  = cone->end();

1740:           for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
1741:             const index_type& ind    = this->_atlas->restrictPoint(patch, *p_iter)[0];
1742:             const int&        start  = ind.index;
1743:             const int&        length = std::abs(ind.prefix);

1745:             for(int i = start; i < start + length; ++i) {
1746:               a[i] = v[++j];
1747:             }
1748:           }
1749:         } else {
1750:           const Obj<IndexArray>& ind = this->getIndices(patch, p);

1752:           for(typename IndexArray::iterator i_iter = ind->begin(); i_iter != ind->end(); ++i_iter) {
1753:             const int& start  = i_iter->index;
1754:             const int& length = std::abs(i_iter->prefix);

1756:             for(int i = start; i < start + length; ++i) {
1757:               a[i] = v[++j];
1758:             }
1759:           }
1760:         }
1761:       };
1762:       // Return only the values associated to this point, not its closure
1763:       const value_type *restrictPoint(const patch_type& patch, const point_type& p) {
1764:         this->checkPatch(patch);
1765:         return &(this->_arrays[patch][this->_atlas->restrictPoint(patch, p)[0].index]);
1766:       };
1767:       // Update only the values associated to this point, not its closure
1768:       void updatePoint(const patch_type& patch, const point_type& p, const value_type v[]) {
1769:         this->checkPatch(patch);
1770:         const index_type& idx = this->_atlas->restrictPoint(patch, p)[0];
1771:         value_type       *a   = &(this->_arrays[patch][idx.index]);

1773:         for(int i = 0; i < idx.prefix; ++i) {
1774:           a[i] = v[i];
1775:         }
1776:       };
1777:       // Update only the values associated to this point, not its closure
1778:       void updateAddPoint(const patch_type& patch, const point_type& p, const value_type v[]) {
1779:         this->checkPatch(patch);
1780:         const index_type& idx = this->_atlas->restrictPoint(patch, p)[0];
1781:         value_type       *a   = &(this->_arrays[patch][idx.index]);

1783:         for(int i = 0; i < idx.prefix; ++i) {
1784:           a[i] += v[i];
1785:         }
1786:       };
1787:       void updatePointBC(const patch_type& patch, const point_type& p, const value_type v[]) {
1788:         this->checkPatch(patch);
1789:         const index_type& idx = this->_atlas->restrictPoint(patch, p)[0];
1790:         value_type       *a   = &(this->_arrays[patch][idx.index]);

1792:         for(int i = 0; i < std::abs(idx.prefix); ++i) {
1793:           a[i] = v[i];
1794:         }
1795:       };
1796:     public: // BC
1797:       void copyBC(const Obj<Section>& section) {
1798:         const typename topology_type::sheaf_type& patches = this->getTopology()->getPatches();

1800:         for(typename topology_type::sheaf_type::const_iterator patch_iter = patches.begin(); patch_iter != patches.end(); ++patch_iter) {
1801:           const typename atlas_type::chart_type& chart = this->_atlas->getPatch(patch_iter->first);

1803:           for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1804:             const index_type& idx = this->_atlas->restrictPoint(patch_iter->first, *p_iter)[0];

1806:             if (idx.prefix < 0) {
1807:               this->updatePointBC(patch_iter->first, *p_iter, section->restrictPoint(patch_iter->first, *p_iter));
1808:             }
1809:           }
1810:         }
1811:       };
1812:     public:
1813:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
1814:         ostringstream txt;
1815:         int rank;

1817:         if (comm == MPI_COMM_NULL) {
1818:           comm = this->comm();
1819:           rank = this->commRank();
1820:         } else {
1821:           MPI_Comm_rank(comm, &rank);
1822:         }
1823:         if (name == "") {
1824:           if(rank == 0) {
1825:             txt << "viewing a Section" << std::endl;
1826:           }
1827:         } else {
1828:           if(rank == 0) {
1829:             txt << "viewing Section '" << name << "'" << std::endl;
1830:           }
1831:         }
1832:         for(typename values_type::const_iterator a_iter = this->_arrays.begin(); a_iter != this->_arrays.end(); ++a_iter) {
1833:           const patch_type&  patch = a_iter->first;
1834:           const value_type  *array = a_iter->second;

1836:           txt << "[" << this->commRank() << "]: Patch " << patch << std::endl;
1837:           const typename atlas_type::chart_type& chart = this->_atlas->getPatch(patch);

1839:           for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1840:             const point_type& point = *p_iter;
1841:             const index_type& idx   = this->_atlas->restrict(patch, point)[0];

1843:             if (idx.prefix != 0) {
1844:               txt << "[" << this->commRank() << "]:   " << point << " dim " << idx.prefix << " offset " << idx.index << "  ";
1845:               for(int i = 0; i < std::abs(idx.prefix); i++) {
1846:                 txt << " " << array[idx.index+i];
1847:               }
1848:               txt << std::endl;
1849:             }
1850:           }
1851:         }
1852:         if (this->_arrays.empty()) {
1853:           txt << "[" << this->commRank() << "]: empty" << std::endl;
1854:         }
1855:         PetscSynchronizedPrintf(comm, txt.str().c_str());
1856:         PetscSynchronizedFlush(comm);
1857:       };
1858:     };

1860:     // An Overlap is a Sifter describing the overlap of two Sieves
1861:     //   Each arrow is local point ---(remote point)---> remote rank right now
1862:     //     For XSifter, this should change to (local patch, local point) ---> (remote rank, remote patch, remote point)

1864:     template<typename Topology_, typename Value_>
1865:     class ConstantSection : public ALE::ParallelObject {
1866:     public:
1867:       typedef Topology_                          topology_type;
1868:       typedef typename topology_type::patch_type patch_type;
1869:       typedef typename topology_type::sieve_type sieve_type;
1870:       typedef typename topology_type::point_type point_type;
1871:       typedef Value_                             value_type;
1872:     protected:
1873:       Obj<topology_type> _topology;
1874:       const value_type   _value;
1875:     public:
1876:       ConstantSection(MPI_Comm comm, const value_type value, const int debug = 0) : ParallelObject(comm, debug), _value(value) {
1877:         this->_topology = new topology_type(comm, debug);
1878:       };
1879:       ConstantSection(const Obj<topology_type>& topology, const value_type value) : ParallelObject(topology->comm(), topology->debug()), _topology(topology), _value(value) {};
1880:       virtual ~ConstantSection() {};
1881:     public:
1882:       bool hasPoint(const patch_type& patch, const point_type& point) const {return true;};
1883:       const value_type *restrict(const patch_type& patch) {return &this->_value;};
1884:       // This should return something the size of the closure
1885:       const value_type *restrict(const patch_type& patch, const point_type& p) {return &this->_value;};
1886:       const value_type *restrictPoint(const patch_type& patch, const point_type& p) {return &this->_value;};
1887:       void update(const patch_type& patch, const point_type& p, const value_type v[]) {
1888:         throw ALE::Exception("Cannot update a ConstantSection");
1889:       };
1890:       void updateAdd(const patch_type& patch, const point_type& p, const value_type v[]) {
1891:         throw ALE::Exception("Cannot update a ConstantSection");
1892:       };
1893:       void updatePoint(const patch_type& patch, const point_type& p, const value_type v[]) {
1894:         throw ALE::Exception("Cannot update a ConstantSection");
1895:       };
1896:       template<typename Input>
1897:       void update(const patch_type& patch, const point_type& p, const Obj<Input>& v) {
1898:         throw ALE::Exception("Cannot update a ConstantSection");
1899:       };
1900:     public:
1901:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
1902:         ostringstream txt;
1903:         int rank;

1905:         if (comm == MPI_COMM_NULL) {
1906:           comm = this->comm();
1907:           rank = this->commRank();
1908:         } else {
1909:           MPI_Comm_rank(comm, &rank);
1910:         }
1911:         if (name == "") {
1912:           if(rank == 0) {
1913:             txt << "viewing a ConstantSection with value " << this->_value << std::endl;
1914:           }
1915:         } else {
1916:           if(rank == 0) {
1917:             txt << "viewing ConstantSection '" << name << "' with value " << this->_value << std::endl;
1918:           }
1919:         }
1920:         PetscSynchronizedPrintf(comm, txt.str().c_str());
1921:         PetscSynchronizedFlush(comm);
1922:       };
1923:     };

1925:     template<typename Point_>
1926:     class DiscreteSieve {
1927:     public:
1928:       typedef Point_                  point_type;
1929:       typedef std::vector<point_type> coneSequence;
1930:       typedef std::vector<point_type> coneSet;
1931:       typedef std::vector<point_type> coneArray;
1932:       typedef std::vector<point_type> supportSequence;
1933:       typedef std::vector<point_type> supportSet;
1934:       typedef std::vector<point_type> supportArray;
1935:       typedef std::set<point_type>    points_type;
1936:       typedef points_type             baseSequence;
1937:       typedef points_type             capSequence;
1938:     protected:
1939:       Obj<points_type>  _points;
1940:       Obj<coneSequence> _empty;
1941:       Obj<coneSequence> _return;
1942:       void _init() {
1943:         this->_points = new points_type();
1944:         this->_empty  = new coneSequence();
1945:         this->_return = new coneSequence();
1946:       };
1947:     public:
1948:       DiscreteSieve() {
1949:         this->_init();
1950:       };
1951:       template<typename Input>
1952:       DiscreteSieve(const Obj<Input>& points) {
1953:         this->_init();
1954:         this->_points->insert(points->begin(), points->end());
1955:       };
1956:       virtual ~DiscreteSieve() {};
1957:     public:
1958:       void addPoint(const point_type& point) {
1959:         this->_points->insert(point);
1960:       };
1961:       template<typename Input>
1962:       void addPoints(const Obj<Input>& points) {
1963:         this->_points->insert(points->begin(), points->end());
1964:       };
1965:       const Obj<coneSequence>& cone(const point_type& p) {return this->_empty;};
1966:       template<typename Input>
1967:       const Obj<coneSequence>& cone(const Input& p) {return this->_empty;};
1968:       const Obj<coneSet>& nCone(const point_type& p, const int level) {
1969:         if (level == 0) {
1970:           return this->closure(p);
1971:         } else {
1972:           return this->_empty;
1973:         }
1974:       };
1975:       const Obj<coneArray>& closure(const point_type& p) {
1976:         this->_return->clear();
1977:         this->_return->push_back(p);
1978:         return this->_return;
1979:       };
1980:       const Obj<supportSequence>& support(const point_type& p) {return this->_empty;};
1981:       template<typename Input>
1982:       const Obj<supportSequence>& support(const Input& p) {return this->_empty;};
1983:       const Obj<supportSet>& nSupport(const point_type& p, const int level) {
1984:         if (level == 0) {
1985:           return this->star(p);
1986:         } else {
1987:           return this->_empty;
1988:         }
1989:       };
1990:       const Obj<supportArray>& star(const point_type& p) {
1991:         this->_return->clear();
1992:         this->_return->push_back(p);
1993:         return this->_return;
1994:       };
1995:       const Obj<capSequence>& roots() {return this->_points;};
1996:       const Obj<capSequence>& cap() {return this->_points;};
1997:       const Obj<baseSequence>& leaves() {return this->_points;};
1998:       const Obj<baseSequence>& base() {return this->_points;};
1999:       template<typename Color>
2000:       void addArrow(const point_type& p, const point_type& q, const Color& color) {
2001:         throw ALE::Exception("Cannot add an arrow to a DiscreteSieve");
2002:       };
2003:       void stratify() {};
2004:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
2005:         ostringstream txt;
2006:         int rank;

2008:         if (comm == MPI_COMM_NULL) {
2009:           comm = MPI_COMM_SELF;
2010:           rank = 0;
2011:         } else {
2012:           MPI_Comm_rank(comm, &rank);
2013:         }
2014:         if (name == "") {
2015:           if(rank == 0) {
2016:             txt << "viewing a DiscreteSieve" << std::endl;
2017:           }
2018:         } else {
2019:           if(rank == 0) {
2020:             txt << "viewing DiscreteSieve '" << name << "'" << std::endl;
2021:           }
2022:         }
2023:         for(typename points_type::const_iterator p_iter = this->_points->begin(); p_iter != this->_points->end(); ++p_iter) {
2024:           txt << "  Point " << *p_iter << std::endl;
2025:         }
2026:         PetscSynchronizedPrintf(comm, txt.str().c_str());
2027:         PetscSynchronizedFlush(comm);
2028:       };
2029:     };


2032:     template<typename Overlap_, typename Topology_, typename Value_>
2033:     class OverlapValues : public Section<Topology_, Value_> {
2034:     public:
2035:       typedef Overlap_                          overlap_type;
2036:       typedef Section<Topology_, Value_>        base_type;
2037:       typedef typename base_type::patch_type    patch_type;
2038:       typedef typename base_type::topology_type topology_type;
2039:       typedef typename base_type::chart_type    chart_type;
2040:       typedef typename base_type::atlas_type    atlas_type;
2041:       typedef typename base_type::value_type    value_type;
2042:       typedef enum {SEND, RECEIVE}              request_type;
2043:       typedef std::map<patch_type, MPI_Request> requests_type;
2044:     protected:
2045:       int           _tag;
2046:       MPI_Datatype  _datatype;
2047:       requests_type _requests;
2048:     public:
2049:       OverlapValues(MPI_Comm comm, const int debug = 0) : Section<Topology_, Value_>(comm, debug) {
2050:         this->_tag      = this->getNewTag();
2051:         this->_datatype = this->getMPIDatatype();
2052:       };
2053:       OverlapValues(MPI_Comm comm, const int tag, const int debug) : Section<Topology_, Value_>(comm, debug), _tag(tag) {
2054:         this->_datatype = this->getMPIDatatype();
2055:       };
2056:       virtual ~OverlapValues() {};
2057:     protected:
2058:       MPI_Datatype getMPIDatatype() {
2059:         if (sizeof(value_type) == 4) {
2060:           return MPI_INT;
2061:         } else if (sizeof(value_type) == 8) {
2062:           return MPI_DOUBLE;
2063:         } else if (sizeof(value_type) == 28) {
2064:           int          blen[2];
2065:           MPI_Aint     indices[2];
2066:           MPI_Datatype oldtypes[2], newtype;
2067:           blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
2068:           blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
2069:           MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
2070:           MPI_Type_commit(&newtype);
2071:           return newtype;
2072:         } else if (sizeof(value_type) == 32) {
2073:           int          blen[2];
2074:           MPI_Aint     indices[2];
2075:           MPI_Datatype oldtypes[2], newtype;
2076:           blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_DOUBLE;
2077:           blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
2078:           MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
2079:           MPI_Type_commit(&newtype);
2080:           return newtype;
2081:         }
2082:         throw ALE::Exception("Cannot determine MPI type for value type");
2083:       };
2084:       int getNewTag() {
2085:         static int tagKeyval = MPI_KEYVAL_INVALID;
2086:         int *tagvalp = NULL, *maxval, flg;

2088:         if (tagKeyval == MPI_KEYVAL_INVALID) {
2089:           PetscMalloc(sizeof(int), &tagvalp);
2090:           MPI_Keyval_create(MPI_NULL_COPY_FN, Petsc_DelTag, &tagKeyval, (void *) NULL);
2091:           MPI_Attr_put(this->_comm, tagKeyval, tagvalp);
2092:           tagvalp[0] = 0;
2093:         }
2094:         MPI_Attr_get(this->_comm, tagKeyval, (void **) &tagvalp, &flg);
2095:         if (tagvalp[0] < 1) {
2096:           MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
2097:           tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
2098:         }
2099:         //std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
2100:         return tagvalp[0]--;
2101:       };
2102:     public: // Accessors
2103:       int getTag() const {return this->_tag;};
2104:       void setTag(const int tag) {this->_tag = tag;};
2105:     public:
2106:       void construct(const int size) {
2107:         const typename topology_type::sheaf_type& patches = this->getTopology()->getPatches();

2109:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
2110:           const Obj<typename topology_type::sieve_type::baseSequence>& base = p_iter->second->base();
2111:           int                                  rank = p_iter->first;

2113:           for(typename topology_type::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2114:             this->setFiberDimension(rank, *b_iter, size);
2115:           }
2116:         }
2117:       };
2118:       template<typename Sizer>
2119:       void construct(const Obj<Sizer>& sizer) {
2120:         const typename topology_type::sheaf_type& patches = this->getTopology()->getPatches();

2122:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
2123:           const Obj<typename topology_type::sieve_type::baseSequence>& base = p_iter->second->base();
2124:           int                                  rank = p_iter->first;

2126:           for(typename topology_type::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2127:             this->setFiberDimension(rank, *b_iter, *(sizer->restrict(rank, *b_iter)));
2128:           }
2129:         }
2130:       };
2131:       void constructCommunication(const request_type& requestType) {
2132:         const typename topology_type::sheaf_type& patches = this->getAtlas()->getTopology()->getPatches();

2134:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
2135:           const patch_type patch = p_iter->first;
2136:           MPI_Request request;

2138:           if (requestType == RECEIVE) {
2139:             if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data(" << this->size(patch) << ") from " << patch << " tag " << this->_tag << std::endl;}
2140:             MPI_Recv_init(this->_arrays[patch], this->size(patch), this->_datatype, patch, this->_tag, this->_comm, &request);
2141:           } else {
2142:             if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << this->size(patch) << ") to " << patch << " tag " << this->_tag << std::endl;}
2143:             MPI_Send_init(this->_arrays[patch], this->size(patch), this->_datatype, patch, this->_tag, this->_comm, &request);
2144:           }
2145:           this->_requests[patch] = request;
2146:         }
2147:       };
2148:       void startCommunication() {
2149:         const typename topology_type::sheaf_type& patches = this->getAtlas()->getTopology()->getPatches();

2151:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
2152:           MPI_Request request = this->_requests[p_iter->first];

2154:           MPI_Start(&request);
2155:         }
2156:       };
2157:       void endCommunication() {
2158:         const typename topology_type::sheaf_type& patches = this->getAtlas()->getTopology()->getPatches();
2159:         MPI_Status status;

2161:         for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
2162:           MPI_Request request = this->_requests[p_iter->first];

2164:           MPI_Wait(&request, &status);
2165:         }
2166:       };
2167:     };
2168:   }
2169: }

2171: #endif