Actual source code: Field.hh

  1: #ifndef included_ALE_Field_hh
  2: #define included_ALE_Field_hh

  4: #ifndef  included_ALE_SieveAlgorithms_hh
  5: #include <sieve/SieveAlgorithms.hh>
  6: #endif


 10: // Sieve need point_type
 11: // Section need point_type and value_type
 12: //   size(), restrict(), update() need orientation which is a Bundle (circular)
 13: // Bundle is Sieve+Section

 15: // An AbstractSection is a mapping from Sieve points to sets of values
 16: //   This is our presentation of a section of a fibre bundle,
 17: //     in which the Topology is the base space, and
 18: //     the value sets are vectors in the fiber spaces
 19: //   The main interface to values is through restrict() and update()
 20: //     This retrieval uses Sieve closure()
 21: //     We should call these rawRestrict/rawUpdate
 22: //   The Section must also be able to set/report the dimension of each fiber
 23: //     for which we use another section we call an \emph{Atlas}
 24: //     For some storage schemes, we also want offsets to go with these dimensions
 25: //   We must have some way of limiting the points associated with values
 26: //     so each section must support a getPatch() call returning the points with associated fibers
 27: //     I was using the Chart for this
 28: //   The Section must be able to participate in \emph{completion}
 29: //     This means restrict to a provided overlap, and exchange in the restricted sections
 30: //     Completion does not use hierarchy, so we see the Topology as a DiscreteTopology
 31: namespace ALE {
 32:   template<typename Point_, typename Alloc_ = malloc_allocator<Point_> >
 33:   class DiscreteSieve {
 34:   public:
 35:     typedef Point_                              point_type;
 36:     typedef Alloc_                              alloc_type;
 37:     typedef std::vector<point_type, alloc_type> coneSequence;
 38:     typedef std::vector<point_type, alloc_type> coneSet;
 39:     typedef std::vector<point_type, alloc_type> coneArray;
 40:     typedef std::vector<point_type, alloc_type> supportSequence;
 41:     typedef std::vector<point_type, alloc_type> supportSet;
 42:     typedef std::vector<point_type, alloc_type> supportArray;
 43:     typedef std::set<point_type, std::less<point_type>, alloc_type>   points_type;
 44:     typedef points_type                                               baseSequence;
 45:     typedef points_type                                               capSequence;
 46:     typedef typename alloc_type::template rebind<points_type>::other  points_alloc_type;
 47:     typedef typename points_alloc_type::pointer                       points_ptr;
 48:     typedef typename alloc_type::template rebind<coneSequence>::other coneSequence_alloc_type;
 49:     typedef typename coneSequence_alloc_type::pointer                 coneSequence_ptr;
 50:   protected:
 51:     Obj<points_type>  _points;
 52:     Obj<coneSequence> _empty;
 53:     Obj<coneSequence> _return;
 54:     alloc_type        _allocator;
 55:     void _init() {
 56:       points_ptr pPoints = points_alloc_type(this->_allocator).allocate(1);
 57:       points_alloc_type(this->_allocator).construct(pPoints, points_type());
 58:       this->_points = Obj<points_type>(pPoints, sizeof(points_type));
 59:       ///this->_points = new points_type();
 60:       coneSequence_ptr pEmpty = coneSequence_alloc_type(this->_allocator).allocate(1);
 61:       coneSequence_alloc_type(this->_allocator).construct(pEmpty, coneSequence());
 62:       this->_empty = Obj<coneSequence>(pEmpty, sizeof(coneSequence));
 63:       ///this->_empty  = new coneSequence();
 64:       coneSequence_ptr pReturn = coneSequence_alloc_type(this->_allocator).allocate(1);
 65:       coneSequence_alloc_type(this->_allocator).construct(pReturn, coneSequence());
 66:       this->_return = Obj<coneSequence>(pReturn, sizeof(coneSequence));
 67:       ///this->_return = new coneSequence();
 68:     };
 69:   public:
 70:     DiscreteSieve() {
 71:       this->_init();
 72:     };
 73:     template<typename Input>
 74:     DiscreteSieve(const Obj<Input>& points) {
 75:       this->_init();
 76:       this->_points->insert(points->begin(), points->end());
 77:     }
 78:     virtual ~DiscreteSieve() {};
 79:   public:
 80:     void addPoint(const point_type& point) {
 81:       this->_points->insert(point);
 82:     }
 83:     template<typename Input>
 84:     void addPoints(const Obj<Input>& points) {
 85:       this->_points->insert(points->begin(), points->end());
 86:     }
 87:     const Obj<coneSequence>& cone(const point_type& p) {return this->_empty;}
 88:     template<typename Input>
 89:     const Obj<coneSequence>& cone(const Input& p) {return this->_empty;}
 90:     const Obj<coneSet>& nCone(const point_type& p, const int level) {
 91:       if (level == 0) {
 92:         return this->closure(p);
 93:       } else {
 94:         return this->_empty;
 95:       }
 96:     }
 97:     const Obj<coneArray>& closure(const point_type& p) {
 98:       this->_return->clear();
 99:       this->_return->push_back(p);
100:       return this->_return;
101:     }
102:     const Obj<supportSequence>& support(const point_type& p) {return this->_empty;}
103:     template<typename Input>
104:     const Obj<supportSequence>& support(const Input& p) {return this->_empty;}
105:     const Obj<supportSet>& nSupport(const point_type& p, const int level) {
106:       if (level == 0) {
107:         return this->star(p);
108:       } else {
109:         return this->_empty;
110:       }
111:     }
112:     const Obj<supportArray>& star(const point_type& p) {
113:       this->_return->clear();
114:       this->_return->push_back(p);
115:       return this->_return;
116:     }
117:     const Obj<capSequence>& roots() {return this->_points;}
118:     const Obj<capSequence>& cap() {return this->_points;}
119:     const Obj<baseSequence>& leaves() {return this->_points;}
120:     const Obj<baseSequence>& base() {return this->_points;}
121:     template<typename Color>
122:     void addArrow(const point_type& p, const point_type& q, const Color& color) {
123:       throw ALE::Exception("Cannot add an arrow to a DiscreteSieve");
124:     }
125:     void stratify() {};
126:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
127:       ostringstream txt;
128:       int rank;

130:       if (comm == MPI_COMM_NULL) {
131:         comm = MPI_COMM_SELF;
132:         rank = 0;
133:       } else {
134:         MPI_Comm_rank(comm, &rank);
135:       }
136:       if (name == "") {
137:         if(rank == 0) {
138:           txt << "viewing a DiscreteSieve" << std::endl;
139:         }
140:       } else {
141:         if(rank == 0) {
142:           txt << "viewing DiscreteSieve '" << name << "'" << std::endl;
143:         }
144:       }
145:       for(typename points_type::const_iterator p_iter = this->_points->begin(); p_iter != this->_points->end(); ++p_iter) {
146:         txt << "  Point " << *p_iter << std::endl;
147:       }
148:       PetscSynchronizedPrintf(comm, txt.str().c_str());
149:       PetscSynchronizedFlush(comm);
150:     };
151:   };
152:   // A ConstantSection is the simplest Section
153:   //   All fibers are dimension 1
154:   //   All values are equal to a constant
155:   //     We need no value storage and no communication for completion
156:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
157:   class ConstantSection : public ALE::ParallelObject {
158:   public:
159:     typedef Point_                                                  point_type;
160:     typedef Value_                                                  value_type;
161:     typedef Alloc_                                                  alloc_type;
162:     typedef std::set<point_type, std::less<point_type>, alloc_type> chart_type;
163:   protected:
164:     chart_type _chart;
165:     value_type _value[2];
166:   public:
167:     ConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
168:       _value[1] = 0;
169:     };
170:     ConstantSection(MPI_Comm comm, const value_type& value, const int debug) : ParallelObject(comm, debug) {
171:       _value[0] = value;
172:       _value[1] = value;
173:     };
174:     ConstantSection(MPI_Comm comm, const value_type& value, const value_type& defaultValue, const int debug) : ParallelObject(comm, debug) {
175:       _value[0] = value;
176:       _value[1] = defaultValue;
177:     };
178:   public: // Verifiers
179:     void checkPoint(const point_type& point) const {
180:       if (this->_chart.find(point) == this->_chart.end()) {
181:         ostringstream msg;
182:         msg << "Invalid section point " << point << std::endl;
183:         throw ALE::Exception(msg.str().c_str());
184:       }
185:     };
186:     void checkDimension(const int& dim) {
187:       if (dim != 1) {
188:         ostringstream msg;
189:         msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
190:         throw ALE::Exception(msg.str().c_str());
191:       }
192:     };
193:     bool hasPoint(const point_type& point) const {
194:       return this->_chart.count(point) > 0;
195:     };
196:   public: // Accessors
197:     const chart_type& getChart() {return this->_chart;};
198:     void addPoint(const point_type& point) {
199:       this->_chart.insert(point);
200:     };
201:     template<typename Points>
202:     void addPoint(const Obj<Points>& points) {
203:       this->_chart.insert(points->begin(), points->end());
204:     }
205:     template<typename Points>
206:     void addPoint(const Points& points) {
207:       this->_chart.insert(points.begin(), points.end());
208:     }
209: //     void addPoint(const std::set<point_type>& points) {
210: //       this->_chart.insert(points.begin(), points.end());
211: //     };
212:     value_type getDefaultValue() {return this->_value[1];};
213:     void setDefaultValue(const value_type value) {this->_value[1] = value;};
214:     void copy(const Obj<ConstantSection>& section) {
215:       const chart_type& chart = section->getChart();

217:       this->addPoint(chart);
218:       this->_value[0] = section->restrictSpace()[0];
219:       this->setDefaultValue(section->getDefaultValue());
220:     };
221:   public: // Sizes
222:     void clear() {
223:       this->_chart.clear();
224:     };
225:     int getFiberDimension(const point_type& p) const {
226:       if (this->hasPoint(p)) return 1;
227:       return 0;
228:     };
229:     void setFiberDimension(const point_type& p, int dim) {
230:       this->checkDimension(dim);
231:       this->addPoint(p);
232:     };
233:     template<typename Sequence>
234:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
235:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
236:         this->setFiberDimension(*p_iter, dim);
237:       }
238:     }
239:     void addFiberDimension(const point_type& p, int dim) {
240:       if (this->_chart.find(p) != this->_chart.end()) {
241:         ostringstream msg;
242:         msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
243:         throw ALE::Exception(msg.str().c_str());
244:       } else {
245:         this->setFiberDimension(p, dim);
246:       }
247:     }
248:     int size(const point_type& p) {return this->getFiberDimension(p);};
249:     void allocatePoint() {};
250:   public: // Restriction
251:     const value_type *restrictSpace() const {
252:       return this->_value;
253:     };
254:     const value_type *restrictPoint(const point_type& p) const {
255:       if (this->hasPoint(p)) {
256:         return this->_value;
257:       }
258:       return &this->_value[1];
259:     };
260:     void updatePoint(const point_type& p, const value_type v[]) {
261:       this->_value[0] = v[0];
262:     };
263:     void updateAddPoint(const point_type& p, const value_type v[]) {
264:       this->_value[0] += v[0];
265:     };
266:   public:
267:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
268:       ostringstream txt;
269:       int rank;

271:       if (comm == MPI_COMM_NULL) {
272:         comm = this->comm();
273:         rank = this->commRank();
274:       } else {
275:         MPI_Comm_rank(comm, &rank);
276:       }
277:       if (name == "") {
278:         if(rank == 0) {
279:           txt << "viewing a ConstantSection" << std::endl;
280:         }
281:       } else {
282:         if(rank == 0) {
283:           txt << "viewing ConstantSection '" << name << "'" << std::endl;
284:         }
285:       }
286:       txt <<"["<<this->commRank()<<"]: Value " << this->_value[0] << " Default Value " << this->_value[1] << std::endl;
287:       PetscSynchronizedPrintf(comm, txt.str().c_str());
288:       PetscSynchronizedFlush(comm);
289:     };
290:   };
291:   // A UniformSection often acts as an Atlas
292:   //   All fibers are the same dimension
293:   //     Note we can use a ConstantSection for this Atlas
294:   //   Each point may have a different vector
295:   //     Thus we need storage for values, and hence must implement completion
296:   template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
297:   class UniformSection : public ALE::ParallelObject {
298:   public:
299:     typedef Point_                                           point_type;
300:     typedef Value_                                           value_type;
301:     typedef Alloc_                                           alloc_type;
302:     typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
303:     typedef ConstantSection<point_type, int, point_alloc_type> atlas_type;
304:     typedef typename atlas_type::chart_type                  chart_type;
305:     typedef struct {value_type v[fiberDim];}                 fiber_type;
306:     typedef typename alloc_type::template rebind<std::pair<const point_type, fiber_type> >::other pair_alloc_type;
307:     typedef std::map<point_type, fiber_type, std::less<point_type>, pair_alloc_type>              values_type;
308:     typedef typename alloc_type::template rebind<atlas_type>::other                               atlas_alloc_type;
309:     typedef typename atlas_alloc_type::pointer                                                    atlas_ptr;
310:   protected:
311:     size_t          _contiguous_array_size;
312:     value_type     *_contiguous_array;
313:     Obj<atlas_type> _atlas;
314:     values_type     _array;
315:     fiber_type      _emptyValue;
316:     alloc_type      _allocator;
317:   public:
318:     UniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _contiguous_array_size(0), _contiguous_array(NULL) {
319:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
320:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
321:       this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
322:       int dim = fiberDim;
323:       this->_atlas->updatePoint(*this->_atlas->getChart().begin(), &dim);
324:       for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
325:     };
326:     UniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _contiguous_array_size(0), _contiguous_array(NULL), _atlas(atlas) {
327:       int dim = fiberDim;
328:       this->_atlas->update(*this->_atlas->getChart().begin(), &dim);
329:       for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
330:     };
331:     virtual ~UniformSection() {
332:       this->_atlas = NULL;
333:       if (this->_contiguous_array) {
334:         for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
335:         this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
336:       }
337:     };
338:   public:
339:     value_type *getRawArray(const int size) {
340:       static value_type *array   = NULL;
341:       static int         maxSize = 0;

343:       if (size > maxSize) {
344:         const value_type dummy(0);

346:         if (array) {
347:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
348:           this->_allocator.deallocate(array, maxSize);
349:           ///delete [] array;
350:         }
351:         maxSize = size;
352:         array   = this->_allocator.allocate(maxSize);
353:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
354:         ///array = new value_type[maxSize];
355:       };
356:       return array;
357:     };
358:   public: // Verifiers
359:     bool hasPoint(const point_type& point) {
360:       return this->_atlas->hasPoint(point);
361:     };
362:     void checkDimension(const int& dim) {
363:       if (dim != fiberDim) {
364:         ostringstream msg;
365:         msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
366:         throw ALE::Exception(msg.str().c_str());
367:       }
368:     };
369:   public: // Accessors
370:     const chart_type& getChart() {return this->_atlas->getChart();};
371:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
372:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
373:     void addPoint(const point_type& point) {
374:       this->setFiberDimension(point, fiberDim);
375:     };
376:     template<typename Points>
377:     void addPoint(const Obj<Points>& points) {
378:       for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
379:         this->setFiberDimension(*p_iter, fiberDim);
380:       }
381:     }
382:     void copy(const Obj<UniformSection>& section) {
383:       this->getAtlas()->copy(section->getAtlas());
384:       const chart_type& chart = section->getChart();

386:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
387:         this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
388:       }
389:     }
390:   public: // Sizes
391:     void clear() {
392:       this->_array.clear();
393:       this->_atlas->clear();
394:     }
395:     int getFiberDimension(const point_type& p) const {
396:       return this->_atlas->restrictPoint(p)[0];
397:     }
398:     void setFiberDimension(const point_type& p, int dim) {
399:       this->update();
400:       this->checkDimension(dim);
401:       this->_atlas->addPoint(p);
402:       this->_atlas->updatePoint(p, &dim);
403:     }
404:     template<typename Sequence>
405:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
406:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
407:         this->setFiberDimension(*p_iter, dim);
408:       }
409:     }
410:     void setFiberDimension(const std::set<point_type>& points, int dim) {
411:       for(typename std::set<point_type>::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
412:         this->setFiberDimension(*p_iter, dim);
413:       }
414:     };
415:     void addFiberDimension(const point_type& p, int dim) {
416:       if (this->hasPoint(p)) {
417:         ostringstream msg;
418:         msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
419:         throw ALE::Exception(msg.str().c_str());
420:       } else {
421:         this->setFiberDimension(p, dim);
422:       }
423:     };
424:     int size() const {
425:       const chart_type& points = this->_atlas->getChart();
426:       int               size   = 0;

428:       for(typename chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
429:         size += this->getFiberDimension(*p_iter);
430:       }
431:       return size;
432:     };
433:     int sizeWithBC() const {
434:       return this->size();
435:     };
436:     void allocatePoint() {};
437:   public: // Restriction
438:     const value_type *restrictSpace() {
439:       const chart_type& chart = this->getChart();
440:       const value_type  dummy = 0;
441:       int               k     = 0;

443:       if (chart.size() > this->_contiguous_array_size*fiberDim) {
444:         if (this->_contiguous_array) {
445:           for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
446:           this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
447:         }
448:         this->_contiguous_array_size = chart.size()*fiberDim;
449:         this->_contiguous_array = this->_allocator.allocate(this->_contiguous_array_size*fiberDim);
450:         for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.construct(this->_contiguous_array+i, dummy);}
451:       }
452:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
453:         const value_type *a = this->_array[*p_iter].v;

455:         for(int i = 0; i < fiberDim; ++i, ++k) {
456:           this->_contiguous_array[k] = a[i];
457:         }
458:       }
459:       return this->_contiguous_array;
460:     };
461:     void update() {
462:       if (this->_contiguous_array) {
463:         const chart_type& chart = this->getChart();
464:         int               k     = 0;

466:         for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
467:           value_type *a = this->_array[*p_iter].v;

469:           for(int i = 0; i < fiberDim; ++i, ++k) {
470:             a[i] = this->_contiguous_array[k];
471:           }
472:         }
473:         for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
474:         this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
475:         this->_contiguous_array_size = 0;
476:         this->_contiguous_array      = NULL;
477:       }
478:     };
479:     // Return only the values associated to this point, not its closure
480:     const value_type *restrictPoint(const point_type& p) {
481:       if (this->_array.find(p) == this->_array.end()) return this->_emptyValue.v;
482:       this->update();
483:       return this->_array[p].v;
484:     };
485:     // Update only the values associated to this point, not its closure
486:     void updatePoint(const point_type& p, const value_type v[]) {
487:       this->update();
488:       for(int i = 0; i < fiberDim; ++i) {
489:         this->_array[p].v[i] = v[i];
490:       }
491:     };
492:     // Update only the values associated to this point, not its closure
493:     void updateAddPoint(const point_type& p, const value_type v[]) {
494:       this->update();
495:       for(int i = 0; i < fiberDim; ++i) {
496:         this->_array[p].v[i] += v[i];
497:       }
498:     };
499:     void updatePointAll(const point_type& p, const value_type v[]) {
500:       this->updatePoint(p, v);
501:     };
502:   public:
503:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
504:       ostringstream txt;
505:       int rank;

507:       this->update();
508:       if (comm == MPI_COMM_NULL) {
509:         comm = this->comm();
510:         rank = this->commRank();
511:       } else {
512:         MPI_Comm_rank(comm, &rank);
513:       }
514:       if (name == "") {
515:         if(rank == 0) {
516:           txt << "viewing a UniformSection" << std::endl;
517:         }
518:       } else {
519:         if(rank == 0) {
520:           txt << "viewing UniformSection '" << name << "'" << std::endl;
521:         }
522:       }
523:       const typename atlas_type::chart_type& chart = this->_atlas->getChart();
524:       values_type&                           array = this->_array;

526:       for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
527:         const point_type&                     point = *p_iter;
528:         const typename atlas_type::value_type dim   = this->_atlas->restrictPoint(point)[0];

530:         if (dim != 0) {
531:           txt << "[" << this->commRank() << "]:   " << point << " dim " << dim << "  ";
532:           for(int i = 0; i < dim; i++) {
533:             txt << " " << array[point].v[i];
534:           }
535:           txt << std::endl;
536:         }
537:       }
538:       if (chart.size() == 0) {
539:         txt << "[" << this->commRank() << "]: empty" << std::endl;
540:       }
541:       PetscSynchronizedPrintf(comm, txt.str().c_str());
542:       PetscSynchronizedFlush(comm);
543:     };
544:   };

546:   // A FauxConstantSection is the simplest Section
547:   //   All fibers are dimension 1
548:   //   All values are equal to a constant
549:   //     We need no value storage and no communication for completion
550:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
551:   class FauxConstantSection : public ALE::ParallelObject {
552:   public:
553:     typedef Point_ point_type;
554:     typedef Value_ value_type;
555:     typedef Alloc_ alloc_type;
556:   protected:
557:     value_type _value;
558:   public:
559:     FauxConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
560:     FauxConstantSection(MPI_Comm comm, const value_type& value, const int debug) : ParallelObject(comm, debug), _value(value) {};
561:   public: // Verifiers
562:     void checkDimension(const int& dim) {
563:       if (dim != 1) {
564:         ostringstream msg;
565:         msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
566:         throw ALE::Exception(msg.str().c_str());
567:       }
568:     };
569:   public: // Accessors
570:     void addPoint(const point_type& point) {
571:     }
572:     template<typename Points>
573:     void addPoint(const Obj<Points>& points) {
574:     }
575:     template<typename Points>
576:     void addPoint(const Points& points) {
577:     }
578:     void copy(const Obj<FauxConstantSection>& section) {
579:       this->_value = section->restrictPoint(point_type())[0];
580:     }
581:   public: // Sizes
582:     void clear() {
583:     };
584:     int getFiberDimension(const point_type& p) const {
585:       return 1;
586:     };
587:     void setFiberDimension(const point_type& p, int dim) {
588:       this->checkDimension(dim);
589:       this->addPoint(p);
590:     };
591:     template<typename Sequence>
592:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
593:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
594:         this->setFiberDimension(*p_iter, dim);
595:       }
596:     }
597:     void addFiberDimension(const point_type& p, int dim) {
598:       this->setFiberDimension(p, dim);
599:     }
600:     int size(const point_type& p) {return this->getFiberDimension(p);}
601:   public: // Restriction
602:     const value_type *restrictPoint(const point_type& p) const {
603:       return &this->_value;
604:     };
605:     void updatePoint(const point_type& p, const value_type v[]) {
606:       this->_value = v[0];
607:     };
608:     void updateAddPoint(const point_type& p, const value_type v[]) {
609:       this->_value += v[0];
610:     };
611:   public:
612:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
613:       ostringstream txt;
614:       int rank;

616:       if (comm == MPI_COMM_NULL) {
617:         comm = this->comm();
618:         rank = this->commRank();
619:       } else {
620:         MPI_Comm_rank(comm, &rank);
621:       }
622:       if (name == "") {
623:         if(rank == 0) {
624:           txt << "viewing a FauxConstantSection" << std::endl;
625:         }
626:       } else {
627:         if(rank == 0) {
628:           txt << "viewing FauxConstantSection '" << name << "'" << std::endl;
629:         }
630:       }
631:       txt <<"["<<this->commRank()<<"]: Value " << this->_value << std::endl;
632:       PetscSynchronizedPrintf(comm, txt.str().c_str());
633:       PetscSynchronizedFlush(comm);
634:     };
635:   };
636:   // Make a simple set from the keys of a map
637:   template<typename Map>
638:   class SetFromMap {
639:   public:
640:     typedef typename Map::size_type size_type;
641:     class const_iterator {
642:     public:
643:       typedef typename Map::key_type  value_type;
644:       typedef typename Map::size_type size_type;
645:     protected:
646:       typename Map::const_iterator _iter;
647:     public:
648:       const_iterator(const typename Map::const_iterator& iter): _iter(iter) {};
649:       ~const_iterator() {};
650:     public:
651:       const_iterator& operator=(const const_iterator& iter) {this->_iter = iter._iter; return this->_iter;};
652:       bool operator==(const const_iterator& iter) const {return this->_iter == iter._iter;};
653:       bool operator!=(const const_iterator& iter) const {return this->_iter != iter._iter;};
654:       const_iterator& operator++() {++this->_iter; return *this;}
655:       const_iterator operator++(int) {
656:         const_iterator tmp(*this);
657:         ++(*this);
658:         return tmp;
659:       };
660:       const_iterator& operator--() {--this->_iter; return *this;}
661:       const_iterator operator--(int) {
662:         const_iterator tmp(*this);
663:         --(*this);
664:         return tmp;
665:       };
666:       value_type operator*() const {return this->_iter->first;};
667:     };
668:   protected:
669:     const Map& _map;
670:   public:
671:     SetFromMap(const Map& map): _map(map) {};
672:   public:
673:     const_iterator begin() const {return const_iterator(this->_map.begin());};
674:     const_iterator end()   const {return const_iterator(this->_map.end());};
675:     size_type      size()  const {return this->_map.size();};
676:   };
677:   // A NewUniformSection often acts as an Atlas
678:   //   All fibers are the same dimension
679:   //     Note we can use a ConstantSection for this Atlas
680:   //   Each point may have a different vector
681:   //     Thus we need storage for values, and hence must implement completion
682:   template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
683:   class NewUniformSection : public ALE::ParallelObject {
684:   public:
685:     typedef Point_                                           point_type;
686:     typedef Value_                                           value_type;
687:     typedef Alloc_                                           alloc_type;
688:     typedef typename alloc_type::template rebind<point_type>::other                               point_alloc_type;
689:     typedef FauxConstantSection<point_type, int, point_alloc_type>                                atlas_type;
690:     typedef struct {value_type v[fiberDim];}                                                      fiber_type;
691:     typedef typename alloc_type::template rebind<std::pair<const point_type, fiber_type> >::other pair_alloc_type;
692:     typedef std::map<point_type, fiber_type, std::less<point_type>, pair_alloc_type>              values_type;
693:     typedef SetFromMap<values_type>                                                               chart_type;
694:     typedef typename alloc_type::template rebind<atlas_type>::other                               atlas_alloc_type;
695:     typedef typename atlas_alloc_type::pointer                                                    atlas_ptr;
696:   protected:
697:     values_type     _array;
698:     chart_type      _chart;
699:     size_t          _contiguous_array_size;
700:     value_type     *_contiguous_array;
701:     Obj<atlas_type> _atlas;
702:     fiber_type      _emptyValue;
703:     alloc_type      _allocator;
704:   public:
705:     NewUniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _chart(_array), _contiguous_array_size(0), _contiguous_array(NULL) {
706:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
707:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
708:       this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
709:       int dim = fiberDim;
710:       this->_atlas->update(point_type(), &dim);
711:       for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
712:     };
713:     NewUniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _chart(_array), _contiguous_array_size(0), _contiguous_array(NULL), _atlas(atlas) {
714:       int dim = fiberDim;
715:       this->_atlas->update(point_type(), &dim);
716:       for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
717:     };
718:     ~NewUniformSection() {
719:       this->_atlas = NULL;
720:       if (this->_contiguous_array) {
721:         for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
722:         this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
723:       }
724:     };
725:   public:
726:     value_type *getRawArray(const int size) {
727:       static value_type *array   = NULL;
728:       static int         maxSize = 0;

730:       if (size > maxSize) {
731:         const value_type dummy(0);

733:         if (array) {
734:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
735:           this->_allocator.deallocate(array, maxSize);
736:           ///delete [] array;
737:         }
738:         maxSize = size;
739:         array   = this->_allocator.allocate(maxSize);
740:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
741:         ///array = new value_type[maxSize];
742:       };
743:       return array;
744:     };
745:   public: // Verifiers
746:     bool hasPoint(const point_type& point) {
747:       return (this->_array.find(point) != this->_array.end());
748:       ///return this->_atlas->hasPoint(point);
749:     };
750:     void checkDimension(const int& dim) {
751:       if (dim != fiberDim) {
752:         ostringstream msg;
753:         msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
754:         throw ALE::Exception(msg.str().c_str());
755:       }
756:     };
757:   public: // Accessors
758:     const chart_type& getChart() {return this->_chart;}
759:     const Obj<atlas_type>& getAtlas() {return this->_atlas;}
760:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;}
761:     void addPoint(const point_type& point) {
762:       this->setFiberDimension(point, fiberDim);
763:     }
764:     template<typename Points>
765:     void addPoint(const Obj<Points>& points) {
766:       for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
767:         this->setFiberDimension(*p_iter, fiberDim);
768:       }
769:     }
770:     void copy(const Obj<NewUniformSection>& section) {
771:       this->getAtlas()->copy(section->getAtlas());
772:       const chart_type& chart = section->getChart();

774:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
775:         this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
776:       }
777:     }
778:   public: // Sizes
779:     void clear() {
780:       this->_array.clear();
781:       this->_atlas->clear();
782:     };
783:     int getFiberDimension(const point_type& p) const {
784:       return fiberDim;
785:     };
786:     void setFiberDimension(const point_type& p, int dim) {
787:       this->update();
788:       this->checkDimension(dim);
789:       this->_atlas->addPoint(p);
790:       this->_atlas->updatePoint(p, &dim);
791:       this->_array[p] = fiber_type();
792:     };
793:     template<typename Sequence>
794:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
795:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
796:         this->setFiberDimension(*p_iter, dim);
797:       }
798:     }
799:     void setFiberDimension(const std::set<point_type>& points, int dim) {
800:       for(typename std::set<point_type>::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
801:         this->setFiberDimension(*p_iter, dim);
802:       }
803:     };
804:     void addFiberDimension(const point_type& p, int dim) {
805:       if (this->hasPoint(p)) {
806:         ostringstream msg;
807:         msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
808:         throw ALE::Exception(msg.str().c_str());
809:       } else {
810:         this->setFiberDimension(p, dim);
811:       }
812:     };
813:     int size() {
814:       const chart_type& points = this->getChart();
815:       int               size   = 0;

817:       for(typename chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
818:         size += this->getFiberDimension(*p_iter);
819:       }
820:       return size;
821:     };
822:     int sizeWithBC() {
823:       return this->size();
824:     };
825:     void allocatePoint() {};
826:   public: // Restriction
827:     // Return a pointer to the entire contiguous storage array
828:     const value_type *restrictSpace() {
829:       const chart_type& chart = this->getChart();
830:       const value_type  dummy = 0;
831:       int               k     = 0;

833:       if (chart.size() > this->_contiguous_array_size*fiberDim) {
834:         if (this->_contiguous_array) {
835:           for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
836:           this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
837:         }
838:         this->_contiguous_array_size = chart.size()*fiberDim;
839:         this->_contiguous_array = this->_allocator.allocate(this->_contiguous_array_size*fiberDim);
840:         for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.construct(this->_contiguous_array+i, dummy);}
841:       }
842:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
843:         const value_type *a = this->_array[*p_iter].v;

845:         for(int i = 0; i < fiberDim; ++i, ++k) {
846:           this->_contiguous_array[k] = a[i];
847:         }
848:       }
849:       return this->_contiguous_array;
850:     };
851:     void update() {
852:       if (this->_contiguous_array) {
853:         const chart_type& chart = this->getChart();
854:         int               k     = 0;

856:         for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
857:           value_type *a = this->_array[*p_iter].v;

859:           for(int i = 0; i < fiberDim; ++i, ++k) {
860:             a[i] = this->_contiguous_array[k];
861:           }
862:         }
863:         for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
864:         this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
865:         this->_contiguous_array_size = 0;
866:         this->_contiguous_array      = NULL;
867:       }
868:     };
869:     // Return only the values associated to this point, not its closure
870:     const value_type *restrictPoint(const point_type& p) {
871:       if (this->_array.find(p) == this->_array.end()) return this->_emptyValue.v;
872:       this->update();
873:       return this->_array[p].v;
874:     };
875:     // Update only the values associated to this point, not its closure
876:     void updatePoint(const point_type& p, const value_type v[]) {
877:       this->update();
878:       for(int i = 0; i < fiberDim; ++i) {
879:         this->_array[p].v[i] = v[i];
880:       }
881:     };
882:     // Update only the values associated to this point, not its closure
883:     void updateAddPoint(const point_type& p, const value_type v[]) {
884:       this->update();
885:       for(int i = 0; i < fiberDim; ++i) {
886:         this->_array[p].v[i] += v[i];
887:       }
888:     };
889:     void updatePointAll(const point_type& p, const value_type v[]) {
890:       this->updatePoint(p, v);
891:     };
892:   public:
893:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
894:       ostringstream txt;
895:       int rank;

897:       this->update();
898:       if (comm == MPI_COMM_NULL) {
899:         comm = this->comm();
900:         rank = this->commRank();
901:       } else {
902:         MPI_Comm_rank(comm, &rank);
903:       }
904:       if (name == "") {
905:         if(rank == 0) {
906:           txt << "viewing a NewUniformSection" << std::endl;
907:         }
908:       } else {
909:         if(rank == 0) {
910:           txt << "viewing NewUniformSection '" << name << "'" << std::endl;
911:         }
912:       }
913:       const chart_type& chart = this->getChart();
914:       values_type&      array = this->_array;

916:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
917:         const point_type& point = *p_iter;

919:         if (fiberDim != 0) {
920:           txt << "[" << this->commRank() << "]:   " << point << " dim " << fiberDim << "  ";
921:           for(int i = 0; i < fiberDim; i++) {
922:             txt << " " << array[point].v[i];
923:           }
924:           txt << std::endl;
925:         }
926:       }
927:       if (chart.size() == 0) {
928:         txt << "[" << this->commRank() << "]: empty" << std::endl;
929:       }
930:       PetscSynchronizedPrintf(comm, txt.str().c_str());
931:       PetscSynchronizedFlush(comm);
932:     };
933:   };
934:   // A Section is our most general construct (more general ones could be envisioned)
935:   //   The Atlas is a UniformSection of dimension 1 and value type Point
936:   //     to hold each fiber dimension and offsets into a contiguous patch array
937:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
938:            typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >
939:   class Section : public ALE::ParallelObject {
940:   public:
941:     typedef Point_                                                  point_type;
942:     typedef Value_                                                  value_type;
943:     typedef Alloc_                                                  alloc_type;
944:     typedef Atlas_                                                  atlas_type;
945:     typedef Point                                                   index_type;
946:     typedef typename atlas_type::chart_type                         chart_type;
947:     typedef value_type *                                            values_type;
948:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
949:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
950:   protected:
951:     Obj<atlas_type> _atlas;
952:     Obj<atlas_type> _atlasNew;
953:     values_type     _array;
954:     alloc_type      _allocator;
955:   public:
956:     Section(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
957:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
958:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
959:       this->_atlas    = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
960:       this->_atlasNew = NULL;
961:       this->_array    = NULL;
962:     };
963:     Section(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL), _array(NULL) {};
964:     virtual ~Section() {
965:       if (this->_array) {
966:         const int totalSize = this->sizeWithBC();

968:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
969:         this->_allocator.deallocate(this->_array, totalSize);
970:         ///delete [] this->_array;
971:         this->_array = NULL;
972:       }
973:     };
974:   public:
975:     value_type *getRawArray(const int size) {
976:       static value_type *array   = NULL;
977:       static int         maxSize = 0;

979:       if (size > maxSize) {
980:         const value_type dummy(0);

982:         if (array) {
983:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
984:           this->_allocator.deallocate(array, maxSize);
985:           ///delete [] array;
986:         }
987:         maxSize = size;
988:         array   = this->_allocator.allocate(maxSize);
989:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
990:         ///array = new value_type[maxSize];
991:       };
992:       return array;
993:     };
994:     int getStorageSize() const {
995:       return this->sizeWithBC();
996:     };
997:   public: // Verifiers
998:     bool hasPoint(const point_type& point) {
999:       return this->_atlas->hasPoint(point);
1000:     };
1001:   public: // Accessors
1002:     const chart_type& getChart() {return this->_atlas->getChart();};
1003:     void setChart(chart_type& chart) {};
1004:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
1005:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1006:     const Obj<atlas_type>& getNewAtlas() {return this->_atlasNew;};
1007:     void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1008:     const chart_type& getChart() const {return this->_atlas->getChart();};
1009:   public: // BC
1010:     // Returns the number of constraints on a point
1011:     int getConstraintDimension(const point_type& p) const {
1012:       return std::max(0, -this->_atlas->restrictPoint(p)->prefix);
1013:     }
1014:     // Set the number of constraints on a point
1015:     //   We only allow the entire point to be constrained, so these will be the
1016:     //   only dofs on the point
1017:     void setConstraintDimension(const point_type& p, const int numConstraints) {
1018:       this->setFiberDimension(p, -numConstraints);
1019:     };
1020:     void addConstraintDimension(const point_type& p, const int numConstraints) {
1021:       throw ALE::Exception("Variable constraint dimensions not available for this Section type");
1022:     };
1023:     void copyBC(const Obj<Section>& section) {
1024:       const chart_type& chart = this->getChart();

1026:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1027:         if (this->getConstraintDimension(*p_iter)) {
1028:           this->updatePointBC(*p_iter, section->restrictPoint(*p_iter));
1029:         }
1030:       }
1031:     };
1032:     void defaultConstraintDof() {};
1033:   public: // Sizes
1034:     void clear() {
1035:       const int totalSize = this->sizeWithBC();

1037:       for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1038:       this->_allocator.deallocate(this->_array, totalSize);
1039:       ///delete [] this->_array;
1040:       this->_array = NULL;
1041:       this->_atlas->clear();
1042:     };
1043:     // Return the total number of dofs on the point (free and constrained)
1044:     int getFiberDimension(const point_type& p) const {
1045:       return std::abs(this->_atlas->restrictPoint(p)->prefix);
1046:     };
1047:     int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
1048:       return std::abs(atlas->restrictPoint(p)->prefix);
1049:     };
1050:     // Set the total number of dofs on the point (free and constrained)
1051:     void setFiberDimension(const point_type& p, int dim) {
1052:       const index_type idx(dim, -1);
1053:       this->_atlas->addPoint(p);
1054:       this->_atlas->updatePoint(p, &idx);
1055:     };
1056:     template<typename Sequence>
1057:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
1058:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1059:         this->setFiberDimension(*p_iter, dim);
1060:       }
1061:     }
1062:     void addFiberDimension(const point_type& p, int dim) {
1063:       if (this->_atlas->hasPoint(p)) {
1064:         const index_type values(dim, 0);
1065:         this->_atlas->updateAddPoint(p, &values);
1066:       } else {
1067:         this->setFiberDimension(p, dim);
1068:       }
1069:     }
1070:     // Return the number of constrained dofs on this point
1071:     //   If constrained, this is equal to the fiber dimension
1072:     //   Otherwise, 0
1073:     int getConstrainedFiberDimension(const point_type& p) const {
1074:       return std::max((index_type::prefix_type) 0, this->_atlas->restrictPoint(p)->prefix);
1075:     };
1076:     // Return the total number of free dofs
1077:     int size() const {
1078:       const chart_type& points = this->getChart();
1079:       int size = 0;

1081:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1082:         size += this->getConstrainedFiberDimension(*p_iter);
1083:       }
1084:       return size;
1085:     };
1086:     // Return the total number of dofs (free and constrained)
1087:     int sizeWithBC() const {
1088:       const chart_type& points = this->getChart();
1089:       int size = 0;

1091:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1092:         size += this->getFiberDimension(*p_iter);
1093:       }
1094:       return size;
1095:     };
1096:   public: // Index retrieval
1097:     const typename index_type::index_type& getIndex(const point_type& p) {
1098:       return this->_atlas->restrictPoint(p)->index;
1099:     };
1100:     void setIndex(const point_type& p, const typename index_type::index_type& index) {
1101:       ((typename atlas_type::value_type *) this->_atlas->restrictPoint(p))->index = index;
1102:     };
1103:     void setIndexBC(const point_type& p, const typename index_type::index_type& index) {
1104:       this->setIndex(p, index);
1105:     };
1106:     void getIndices(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1107:       this->getIndices(p, this->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1108:     };
1109:     template<typename Order_>
1110:     void getIndices(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1111:       this->getIndices(p, order->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1112:     }
1113:     void getIndices(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1114:       const int& dim   = this->getFiberDimension(p);
1115:       const int& cDim  = this->getConstraintDimension(p);
1116:       const int  end   = start + dim;

1118:       if (!cDim) {
1119:         if (orientation >= 0) {
1120:           for(int i = start; i < end; ++i) {
1121:             indices[(*indx)++] = i;
1122:           }
1123:         } else {
1124:           for(int i = end-1; i >= start; --i) {
1125:             indices[(*indx)++] = i;
1126:           }
1127:         }
1128:       } else {
1129:         if (!freeOnly) {
1130:           for(int i = start; i < end; ++i) {
1131:             indices[(*indx)++] = -1;
1132:           }
1133:         }
1134:       }
1135:     };
1136:   public: // Allocation
1137:     void allocateStorage() {
1138:       const int totalSize = this->sizeWithBC();
1139:       const value_type dummy(0);

1141:       this->_array = this->_allocator.allocate(totalSize);
1142:       ///this->_array = new value_type[totalSize];
1143:       for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1144:       ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1145:     };
1146:     void replaceStorage(value_type *newArray) {
1147:       const int totalSize = this->sizeWithBC();

1149:       for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1150:       this->_allocator.deallocate(this->_array, totalSize);
1151:       ///delete [] this->_array;
1152:       this->_array    = newArray;
1153:       this->_atlasNew = NULL;
1154:     };
1155:     void addPoint(const point_type& point, const int dim) {
1156:       if (dim == 0) return;
1157:       if (this->_atlasNew.isNull()) {
1158:         this->_atlasNew = new atlas_type(this->comm(), this->debug());
1159:         this->_atlasNew->copy(this->_atlas);
1160:       }
1161:       const index_type idx(dim, -1);
1162:       this->_atlasNew->addPoint(point);
1163:       this->_atlasNew->updatePoint(point, &idx);
1164:     };
1165:     template<typename Sequence>
1166:     void addPoints(const Obj<Sequence>& points, const int dim) {
1167:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1168:         this->addPoint(*p_iter, dim);
1169:       }
1170:     }
1171:     void orderPoints(const Obj<atlas_type>& atlas){
1172:       const chart_type& chart    = this->getChart();
1173:       int               offset   = 0;
1174:       int               bcOffset = this->size();

1176:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1177:         typename atlas_type::value_type idx  = atlas->restrictPoint(*c_iter)[0];
1178:         const int&                      dim  = idx.prefix;

1180:         if (dim >= 0) {
1181:           idx.index = offset;
1182:           atlas->updatePoint(*c_iter, &idx);
1183:           offset += dim;
1184:         } else {
1185:           idx.index = bcOffset;
1186:           atlas->updatePoint(*c_iter, &idx);
1187:           bcOffset -= dim;
1188:         }
1189:       }
1190:     };
1191:     void allocatePoint() {
1192:       this->orderPoints(this->_atlas);
1193:       this->allocateStorage();
1194:     };
1195:   public: // Restriction and Update
1196:     // Zero entries
1197:     void zero() {
1198:       memset(this->_array, 0, this->size()* sizeof(value_type));
1199:     };
1200:     // Return a pointer to the entire contiguous storage array
1201:     const value_type *restrictSpace() {
1202:       return this->_array;
1203:     };
1204:     // Update the entire contiguous storage array
1205:     void update(const value_type v[]) {
1206:       const value_type *array = this->_array;
1207:       const int         size  = this->size();

1209:       for(int i = 0; i < size; i++) {
1210:         array[i] = v[i];
1211:       }
1212:     };
1213:     // Return the free values on a point
1214:     const value_type *restrictPoint(const point_type& p) {
1215:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1216:     };
1217:     // Update the free values on a point
1218:     void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1219:       const index_type& idx = this->_atlas->restrictPoint(p)[0];
1220:       value_type       *a   = &(this->_array[idx.index]);

1222:       if (orientation >= 0) {
1223:         for(int i = 0; i < idx.prefix; ++i) {
1224:           a[i] = v[i];
1225:         }
1226:       } else {
1227:         const int last = idx.prefix-1;

1229:         for(int i = 0; i < idx.prefix; ++i) {
1230:           a[i] = v[last-i];
1231:         }
1232:       }
1233:     };
1234:     // Update the free values on a point
1235:     void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1236:       const index_type& idx = this->_atlas->restrictPoint(p)[0];
1237:       value_type       *a   = &(this->_array[idx.index]);

1239:       if (orientation >= 0) {
1240:         for(int i = 0; i < idx.prefix; ++i) {
1241:           a[i] += v[i];
1242:         }
1243:       } else {
1244:         const int last = idx.prefix-1;

1246:         for(int i = 0; i < idx.prefix; ++i) {
1247:           a[i] += v[last-i];
1248:         }
1249:       }
1250:     };
1251:     // Update only the constrained dofs on a point
1252:     void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
1253:       const index_type& idx = this->_atlas->restrictPoint(p)[0];
1254:       const int         dim = -idx.prefix;
1255:       value_type       *a   = &(this->_array[idx.index]);

1257:       if (orientation >= 0) {
1258:         for(int i = 0; i < dim; ++i) {
1259:           a[i] = v[i];
1260:         }
1261:       } else {
1262:         const int last = dim-1;

1264:         for(int i = 0; i < dim; ++i) {
1265:           a[i] = v[last-i];
1266:         }
1267:       }
1268:     };
1269:     // Update all dofs on a point (free and constrained)
1270:     void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
1271:       const index_type& idx = this->_atlas->restrictPoint(p)[0];
1272:       const int         dim = std::abs(idx.prefix);
1273:       value_type       *a   = &(this->_array[idx.index]);

1275:       if (orientation >= 0) {
1276:         for(int i = 0; i < dim; ++i) {
1277:           a[i] = v[i];
1278:         }
1279:       } else {
1280:         const int last = dim-1;

1282:         for(int i = 0; i < dim; ++i) {
1283:           a[i] = v[last-i];
1284:         }
1285:       }
1286:     };
1287:   public:
1288:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
1289:       ostringstream txt;
1290:       int rank;

1292:       if (comm == MPI_COMM_NULL) {
1293:         comm = this->comm();
1294:         rank = this->commRank();
1295:       } else {
1296:         MPI_Comm_rank(comm, &rank);
1297:       }
1298:       if(rank == 0) {
1299:         txt << "viewing Section " << this->getName() << std::endl;
1300:         if (name != "") {
1301:           txt << ": " << name << "'";
1302:         }
1303:         txt << std::endl;
1304:       }
1305:       const typename atlas_type::chart_type& chart = this->_atlas->getChart();
1306:       const value_type                      *array = this->_array;

1308:       for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1309:         const point_type& point = *p_iter;
1310:         const index_type& idx   = this->_atlas->restrictPoint(point)[0];
1311:         const int         dim   = this->getFiberDimension(point);

1313:         if (idx.prefix != 0) {
1314:           txt << "[" << this->commRank() << "]:   " << point << " dim " << idx.prefix << " offset " << idx.index << "  ";
1315:           for(int i = 0; i < dim; i++) {
1316:             txt << " " << array[idx.index+i];
1317:           }
1318:           txt << std::endl;
1319:         }
1320:       }
1321:       if (chart.size() == 0) {
1322:         txt << "[" << this->commRank() << "]: empty" << std::endl;
1323:       }
1324:       PetscSynchronizedPrintf(comm, txt.str().c_str());
1325:       PetscSynchronizedFlush(comm);
1326:     };
1327:   public: // Fibrations
1328:     void setConstraintDof(const point_type& p, const int dofs[]) {};
1329:     int getNumSpaces() const {return 1;};
1330:     void addSpace() {};
1331:     int getFiberDimension(const point_type& p, const int space) {return this->getFiberDimension(p);};
1332:     void setFiberDimension(const point_type& p, int dim, const int space) {this->setFiberDimension(p, dim);};
1333:     template<typename Sequence>
1334:     void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
1335:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1336:         this->setFiberDimension(*p_iter, dim, space);
1337:       }
1338:     }
1339:     void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
1340:       this->setConstraintDimension(p, numConstraints);
1341:     }
1342:   };
1343:   // GeneralSection will support BC on a subset of unknowns on a point
1344:   //   We make a separate constraint Atlas to mark constrained dofs on a point
1345:   //   Storage will be contiguous by node, just as in Section
1346:   //     This allows fast restrict(p)
1347:   //     Then update() is accomplished by skipping constrained unknowns
1348:   //     We must eliminate restrictSpace() since it does not correspond to the constrained system
1349:   //   Numbering will have to be rewritten to calculate correct mappings
1350:   //     I think we can just generate multiple tuples per point
1351:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
1352:            typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
1353:            typename BCAtlas_ = Section<Point_, int, typename Alloc_::template rebind<int>::other> >
1354:   class GeneralSection : public ALE::ParallelObject {
1355:   public:
1356:     typedef Point_                                                  point_type;
1357:     typedef Value_                                                  value_type;
1358:     typedef Alloc_                                                  alloc_type;
1359:     typedef Atlas_                                                  atlas_type;
1360:     typedef BCAtlas_                                                bc_type;
1361:     typedef Point                                                   index_type;
1362:     typedef typename atlas_type::chart_type                         chart_type;
1363:     typedef value_type *                                            values_type;
1364:     typedef std::pair<const int *, const int *>                     customAtlasInd_type;
1365:     typedef std::pair<customAtlasInd_type, bool>                    customAtlas_type;
1366:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
1367:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
1368:     typedef typename alloc_type::template rebind<bc_type>::other    bc_alloc_type;
1369:     typedef typename bc_alloc_type::pointer                         bc_ptr;
1370:   protected:
1371:     // Describes layout of storage, point --> (# of values, offset into array)
1372:     Obj<atlas_type> _atlas;
1373:     // Spare atlas to allow dynamic updates
1374:     Obj<atlas_type> _atlasNew;
1375:     // Storage
1376:     values_type     _array;
1377:     bool            _sharedStorage;
1378:     int             _sharedStorageSize;
1379:     // A section that maps points to sets of constrained local dofs
1380:     Obj<bc_type>    _bc;
1381:     alloc_type      _allocator;
1382:     // Fibration structures
1383:     //   _spaces is a set of atlases which describe the layout of each in the storage of this section
1384:     //   _bcs is the same as _bc, but for each field
1385:     std::vector<Obj<atlas_type> > _spaces;
1386:     std::vector<Obj<bc_type> >    _bcs;
1387:     // Optimization
1388:     std::vector<customAtlas_type> _customRestrictAtlas;
1389:     std::vector<customAtlas_type> _customUpdateAtlas;
1390:   public:
1391:     GeneralSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
1392:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
1393:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
1394:       this->_atlas         = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
1395:       bc_ptr pBC           = bc_alloc_type(this->_allocator).allocate(1);
1396:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
1397:       this->_bc            = Obj<bc_type>(pBC, sizeof(bc_type));
1398:       this->_atlasNew      = NULL;
1399:       this->_array         = NULL;
1400:       this->_sharedStorage = false;
1401:     };
1402:     GeneralSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
1403:       bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
1404:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(this->comm(), this->debug()));
1405:       this->_bc  = Obj<bc_type>(pBC, sizeof(bc_type));
1406:     };
1407:     ~GeneralSection() {
1408:       if (this->_array && !this->_sharedStorage) {
1409:         const int totalSize = this->sizeWithBC();

1411:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1412:         this->_allocator.deallocate(this->_array, totalSize);
1413:         ///delete [] this->_array;
1414:         this->_array = NULL;
1415:       }
1416:       for(std::vector<customAtlas_type>::iterator a_iter = this->_customRestrictAtlas.begin(); a_iter != this->_customRestrictAtlas.end(); ++a_iter) {
1417:         if (a_iter->second) {
1418:           delete [] a_iter->first.first;
1419:           delete [] a_iter->first.second;
1420:         }
1421:       }
1422:       for(std::vector<customAtlas_type>::iterator a_iter = this->_customUpdateAtlas.begin(); a_iter != this->_customUpdateAtlas.end(); ++a_iter) {
1423:         if (a_iter->second) {
1424:           delete [] a_iter->first.first;
1425:           delete [] a_iter->first.second;
1426:         }
1427:       }
1428:     };
1429:   public:
1430:     value_type *getRawArray(const int size) {
1431:       // Put in a sentinel value that deallocates the array
1432:       static value_type *array   = NULL;
1433:       static int         maxSize = 0;

1435:       if (size > maxSize) {
1436:         const value_type dummy(0);

1438:         if (array) {
1439:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
1440:           this->_allocator.deallocate(array, maxSize);
1441:           ///delete [] array;
1442:         }
1443:         maxSize = size;
1444:         array   = this->_allocator.allocate(maxSize);
1445:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
1446:         ///array = new value_type[maxSize];
1447:       };
1448:       return array;
1449:     };
1450:     int getStorageSize() const {
1451:       if (this->_sharedStorage) {
1452:         return this->_sharedStorageSize;
1453:       }
1454:       return this->sizeWithBC();
1455:     };
1456:     bool sharedStorage() const {return this->_sharedStorage;};
1457:   public: // Verifiers
1458:     bool hasPoint(const point_type& point) {
1459:       return this->_atlas->hasPoint(point);
1460:     };
1461:   public: // Accessors
1462:     const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
1463:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1464:     const Obj<atlas_type>& getNewAtlas() const {return this->_atlasNew;};
1465:     void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1466:     const Obj<bc_type>& getBC() const {return this->_bc;};
1467:     void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
1468:     const chart_type& getChart() const {return this->_atlas->getChart();};
1469:     void setChart(const chart_type& chart) {throw ALE::Exception("setChart() for GeneralSection is invalid");};
1470:   public: // BC
1471:     // Returns the number of constraints on a point
1472:     int getConstraintDimension(const point_type& p) const {
1473:       if (!this->_bc->hasPoint(p)) return 0;
1474:       return this->_bc->getFiberDimension(p);
1475:     }
1476:     // Set the number of constraints on a point
1477:     void setConstraintDimension(const point_type& p, const int numConstraints) {
1478:       this->_bc->setFiberDimension(p, numConstraints);
1479:     }
1480:     // Increment the number of constraints on a point
1481:     void addConstraintDimension(const point_type& p, const int numConstraints) {
1482:       this->_bc->addFiberDimension(p, numConstraints);
1483:     }
1484:     // Return the local dofs which are constrained on a point
1485:     const int *getConstraintDof(const point_type& p) const {
1486:       return this->_bc->restrictPoint(p);
1487:     }
1488:     // Set the local dofs which are constrained on a point
1489:     void setConstraintDof(const point_type& p, const int dofs[]) {
1490:       this->_bc->updatePoint(p, dofs);
1491:     }
1492:     template<typename OtherSection>
1493:     void copyBC(const Obj<OtherSection>& section) {
1494:       this->setBC(section->getBC());
1495:       const chart_type& chart = this->getChart();

1497:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1498:         if (this->getConstraintDimension(*p_iter)) {
1499:           this->updatePointBCFull(*p_iter, section->restrictPoint(*p_iter));
1500:         }
1501:       }
1502:       this->copySpaces(section);
1503:     }
1504:     void defaultConstraintDof() {
1505:       const chart_type& chart = this->getChart();
1506:       int size = 0;

1508:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1509:         size = std::max(size, this->getConstraintDimension(*p_iter));
1510:       }
1511:       int *dofs = new int[size];
1512:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1513:         const int cDim = this->getConstraintDimension(*p_iter);

1515:         if (cDim) {
1516:           for(int d = 0; d < cDim; ++d) {
1517:             dofs[d] = d;
1518:           }
1519:           this->_bc->updatePoint(*p_iter, dofs);
1520:         }
1521:       }
1522:       delete [] dofs;
1523:     };
1524:   public: // Sizes
1525:     void clear() {
1526:       if (!this->_sharedStorage) {
1527:         const int totalSize = this->sizeWithBC();

1529:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1530:         this->_allocator.deallocate(this->_array, totalSize);
1531:         ///delete [] this->_array;
1532:       }
1533:       this->_array = NULL;
1534:       this->_atlas->clear();
1535:       this->_bc->clear();
1536:     };
1537:     // Return the total number of dofs on the point (free and constrained)
1538:     int getFiberDimension(const point_type& p) const {
1539:       return this->_atlas->restrictPoint(p)->prefix;
1540:     };
1541:     int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
1542:       return atlas->restrictPoint(p)->prefix;
1543:     };
1544:     // Set the total number of dofs on the point (free and constrained)
1545:     void setFiberDimension(const point_type& p, int dim) {
1546:       const index_type idx(dim, -1);
1547:       this->_atlas->addPoint(p);
1548:       this->_atlas->updatePoint(p, &idx);
1549:     };
1550:     template<typename Sequence>
1551:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
1552:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1553:         this->setFiberDimension(*p_iter, dim);
1554:       }
1555:     }
1556:     void addFiberDimension(const point_type& p, int dim) {
1557:       if (this->_atlas->hasPoint(p)) {
1558:         const index_type values(dim, 0);
1559:         this->_atlas->updateAddPoint(p, &values);
1560:       } else {
1561:         this->setFiberDimension(p, dim);
1562:       }
1563:     };
1564:     // Return the number of constrained dofs on this point
1565:     int getConstrainedFiberDimension(const point_type& p) const {
1566:       return this->getFiberDimension(p) - this->getConstraintDimension(p);
1567:     };
1568:     // Return the total number of free dofs
1569:     int size() const {
1570:       const chart_type& points = this->getChart();
1571:       int               size   = 0;

1573:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1574:         size += this->getConstrainedFiberDimension(*p_iter);
1575:       }
1576:       return size;
1577:     };
1578:     // Return the total number of dofs (free and constrained)
1579:     int sizeWithBC() const {
1580:       const chart_type& points = this->getChart();
1581:       int               size   = 0;

1583:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1584:         size += this->getFiberDimension(*p_iter);
1585:       }
1586:       return size;
1587:     };
1588:   public: // Index retrieval
1589:     const typename index_type::index_type& getIndex(const point_type& p) {
1590:       return this->_atlas->restrictPoint(p)->index;
1591:     };
1592:     void setIndex(const point_type& p, const typename index_type::index_type& index) {
1593:       ((typename atlas_type::value_type *) this->_atlas->restrictPoint(p))->index = index;
1594:     };
1595:     void setIndexBC(const point_type& p, const typename index_type::index_type& index) {};
1596:     void getIndicesRaw(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1) {
1597:       this->getIndicesRaw(p, this->getIndex(p), indices, indx, orientation);
1598:     };
1599:     template<typename Order_>
1600:     void getIndicesRaw(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1) {
1601:       this->getIndicesRaw(p, order->getIndex(p), indices, indx, orientation);
1602:     }
1603:     void getIndices(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = true) {
1604:       this->getIndices(p, this->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1605:     };
1606:     template<typename Order_>
1607:     void getIndices(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1608:       this->getIndices(p, order->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1609:     }
1610:     void getIndicesRaw(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation) {
1611:       if (orientation >= 0) {
1612:         const int& dim = this->getFiberDimension(p);
1613:         const int  end = start + dim;

1615:         for(int i = start; i < end; ++i) {
1616:           indices[(*indx)++] = i;
1617:         }
1618:       } else {
1619:         const int numSpaces = this->getNumSpaces();
1620:         int offset = start;

1622:         for(int space = 0; space < numSpaces; ++space) {
1623:           const int& dim = this->getFiberDimension(p, space);

1625:           for(int i = offset+dim-1; i >= offset; --i) {
1626:             indices[(*indx)++] = i;
1627:           }
1628:           offset += dim;
1629:         }
1630:         if (!numSpaces) {
1631:           const int& dim = this->getFiberDimension(p);

1633:           for(int i = offset+dim-1; i >= offset; --i) {
1634:             indices[(*indx)++] = i;
1635:           }
1636:           offset += dim;
1637:         }
1638:       }
1639:     }
1640:     /*
1641:       p           - The Sieve point
1642:       start       - Offset for the set of indices on this point
1643:       indices     - Storage for the indices
1644:       indx        - Pointer to an offset into the indices argument (to allow composition of index sets)
1645:       orientation - A negative argument reverses the indices
1646:       freeOnly        - If false, include constrained dofs with negative indices, otherwise leave them out
1647:       skipConstraints - If true, when a constrained dof is encountered, increment the index (even though it is not placed in indices[])
1648:      */
1649:     void getIndices(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1650:       const int& cDim = this->getConstraintDimension(p);

1652:       if (!cDim) {
1653:         this->getIndicesRaw(p, start, indices, indx, orientation);
1654:       } else {
1655:         if (orientation >= 0) {
1656:           const int&                          dim  = this->getFiberDimension(p);
1657:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1658:           int                                 cInd = 0;

1660:           /* i is the returned index, k is the local dof number */
1661:           for(int i = start, k = 0; k < dim; ++k) {
1662:             if ((cInd < cDim) && (k == cDof[cInd])) {
1663:               /* Constrained dof */
1664:               if (!freeOnly) indices[(*indx)++] = -(k+1);
1665:               if (skipConstraints) ++i;
1666:               ++cInd;
1667:             } else {
1668:               /* Unconstrained dof */
1669:               indices[(*indx)++] = i++;
1670:             }
1671:           }
1672:         } else {
1673:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1674:           int                                 offset  = 0;
1675:           int                                 cOffset = 0;
1676:           int                                 j       = -1;

1678:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1679:             const int  dim = this->getFiberDimension(p, space);
1680:             const int tDim = this->getConstrainedFiberDimension(p, space);
1681:             int       cInd = (dim - tDim)-1;

1683:             j += dim;
1684:             for(int i = 0, k = start+tDim+offset; i < dim; ++i, --j) {
1685:               if ((cInd >= 0) && (j == cDof[cInd+cOffset])) {
1686:                 if (!freeOnly) indices[(*indx)++] = -(offset+i+1);
1687:                 if (skipConstraints) --k;
1688:                 --cInd;
1689:               } else {
1690:                 indices[(*indx)++] = --k;
1691:               }
1692:             }
1693:             j       += dim;
1694:             offset  += dim;
1695:             cOffset += dim - tDim;
1696:           }
1697:         }
1698:       }
1699:     };
1700:   public: // Allocation
1701:     void allocateStorage() {
1702:       const int totalSize = this->sizeWithBC();
1703:       const value_type dummy(0) ;

1705:       this->_array             = this->_allocator.allocate(totalSize);
1706:       ///this->_array             = new value_type[totalSize];
1707:       this->_sharedStorage     = false;
1708:       this->_sharedStorageSize = 0;
1709:       for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1710:       ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1711:       this->_bc->allocatePoint();
1712:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = this->_bcs.begin(); b_iter != this->_bcs.end(); ++b_iter) {
1713:         (*b_iter)->allocatePoint();;
1714:       }
1715:     };
1716:     void replaceStorage(value_type *newArray, const bool sharedStorage = false, const int sharedStorageSize = 0) {
1717:       if (this->_array && !this->_sharedStorage) {
1718:         const int totalSize = this->sizeWithBC();

1720:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1721:         this->_allocator.deallocate(this->_array, totalSize);
1722:         ///delete [] this->_array;
1723:       }
1724:       this->_array             = newArray;
1725:       this->_sharedStorage     = sharedStorage;
1726:       this->_sharedStorageSize = sharedStorageSize;
1727:       this->_atlas             = this->_atlasNew;
1728:       this->_atlasNew          = NULL;
1729:     };
1730:     // DANGEROUS
1731:     void setStorage(value_type *newArray) {
1732:       this->_array = newArray;
1733:     };
1734:     void addPoint(const point_type& point, const int dim) {
1735:       if (dim == 0) return;
1736:       if (this->_atlasNew.isNull()) {
1737:         this->_atlasNew = new atlas_type(this->comm(), this->debug());
1738:         this->_atlasNew->copy(this->_atlas);
1739:       }
1740:       const index_type idx(dim, -1);
1741:       this->_atlasNew->addPoint(point);
1742:       this->_atlasNew->updatePoint(point, &idx);
1743:     };
1744:     void orderPoints(const Obj<atlas_type>& atlas){
1745:       const chart_type& chart  = this->getChart();
1746:       int               offset = 0;

1748:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1749:         typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
1750:         const int&                      dim = idx.prefix;

1752:         idx.index = offset;
1753:         atlas->updatePoint(*c_iter, &idx);
1754:         offset += dim;
1755:       }
1756:     };
1757:     void allocatePoint() {
1758:       this->orderPoints(this->_atlas);
1759:       this->allocateStorage();
1760:     };
1761:   public: // Restriction and Update
1762:     // Zero entries
1763:     void zero() {
1764:       this->set((value_type) 0.0);
1765:     };
1766:     void zeroWithBC() {
1767:       memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1768:     };
1769:     void set(const value_type value) {
1770:       //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1771:       const chart_type& chart = this->getChart();

1773:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1774:         value_type *array = (value_type *) this->restrictPoint(*c_iter);
1775:         const int&  dim   = this->getFiberDimension(*c_iter);
1776:         const int&  cDim  = this->getConstraintDimension(*c_iter);

1778:         if (!cDim) {
1779:           for(int i = 0; i < dim; ++i) {
1780:             array[i] = value;
1781:           }
1782:         } else {
1783:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1784:           int                                 cInd = 0;

1786:           for(int i = 0; i < dim; ++i) {
1787:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1788:             array[i] = value;
1789:           }
1790:         }
1791:       }
1792:     };
1793:     // Add two sections and put the result in a third
1794:     void add(const Obj<GeneralSection>& x, const Obj<GeneralSection>& y) {
1795:       // Check atlases
1796:       const chart_type& chart = this->getChart();

1798:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1799:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
1800:         const value_type *xArray = x->restrictPoint(*c_iter);
1801:         const value_type *yArray = y->restrictPoint(*c_iter);
1802:         const int&        dim    = this->getFiberDimension(*c_iter);
1803:         const int&        cDim   = this->getConstraintDimension(*c_iter);

1805:         if (!cDim) {
1806:           for(int i = 0; i < dim; ++i) {
1807:             array[i] = xArray[i] + yArray[i];
1808:           }
1809:         } else {
1810:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1811:           int                                 cInd = 0;

1813:           for(int i = 0; i < dim; ++i) {
1814:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1815:             array[i] = xArray[i] + yArray[i];
1816:           }
1817:         }
1818:       }
1819:     };
1820:     // this = this + alpha * x
1821:     template<typename OtherSection>
1822:     void axpy(const value_type alpha, const Obj<OtherSection>& x) {
1823:       // Check atlases
1824:       const chart_type& chart = this->getChart();

1826:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1827:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
1828:         const value_type *xArray = x->restrictPoint(*c_iter);
1829:         const int&        dim    = this->getFiberDimension(*c_iter);
1830:         const int&        cDim   = this->getConstraintDimension(*c_iter);

1832:         if (!cDim) {
1833:           for(int i = 0; i < dim; ++i) {
1834:             array[i] += alpha*xArray[i];
1835:           }
1836:         } else {
1837:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1838:           int                                 cInd = 0;

1840:           for(int i = 0; i < dim; ++i) {
1841:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1842:             array[i] += alpha*xArray[i];
1843:           }
1844:         }
1845:       }
1846:     }
1847:     // Return the free values on a point
1848:     const value_type *restrictSpace() const {
1849:       return this->_array;
1850:     }
1851:     // Return the free values on a point
1852:     const value_type *restrictPoint(const point_type& p) const {
1853:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1854:     }
1855:     void restrictPoint(const point_type& p, value_type values[], const int size) const {
1856:       assert(this->_atlas->restrictPoint(p)[0].prefix == size);
1857:       const value_type *v = &(this->_array[this->_atlas->restrictPoint(p)[0].index]);

1859:       for(int i = 0; i < size; ++i) {
1860:         values[i] = v[i];
1861:       }
1862:     };
1863:     // Update the free values on a point
1864:     //   Takes a full array and ignores constrained values
1865:     void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1866:       value_type *array = (value_type *) this->restrictPoint(p);
1867:       const int&  cDim  = this->getConstraintDimension(p);

1869:       if (!cDim) {
1870:         if (orientation >= 0) {
1871:           const int& dim = this->getFiberDimension(p);

1873:           for(int i = 0; i < dim; ++i) {
1874:             array[i] = v[i];
1875:           }
1876:         } else {
1877:           int offset = 0;
1878:           int j      = -1;

1880:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1881:             const int& dim = this->getFiberDimension(p, space);

1883:             for(int i = dim-1; i >= 0; --i) {
1884:               array[++j] = v[i+offset];
1885:             }
1886:             offset += dim;
1887:           }
1888:         }
1889:       } else {
1890:         if (orientation >= 0) {
1891:           const int&                          dim  = this->getFiberDimension(p);
1892:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1893:           int                                 cInd = 0;

1895:           for(int i = 0; i < dim; ++i) {
1896:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1897:             array[i] = v[i];
1898:           }
1899:         } else {
1900:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1901:           int                                 offset  = 0;
1902:           int                                 cOffset = 0;
1903:           int                                 j       = 0;

1905:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1906:             const int  dim = this->getFiberDimension(p, space);
1907:             const int tDim = this->getConstrainedFiberDimension(p, space);
1908:             const int sDim = dim - tDim;
1909:             int       cInd = 0;

1911:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1912:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1913:               array[j] = v[k];
1914:             }
1915:             offset  += dim;
1916:             cOffset += dim - tDim;
1917:           }
1918:         }
1919:       }
1920:     };
1921:     // Update the free values on a point
1922:     //   Takes a full array and ignores constrained values
1923:     void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1924:       value_type *array = (value_type *) this->restrictPoint(p);
1925:       const int&  cDim  = this->getConstraintDimension(p);

1927:       if (!cDim) {
1928:         if (orientation >= 0) {
1929:           const int& dim = this->getFiberDimension(p);

1931:           for(int i = 0; i < dim; ++i) {
1932:             array[i] += v[i];
1933:           }
1934:         } else {
1935:           int offset = 0;
1936:           int j      = -1;

1938:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1939:             const int& dim = this->getFiberDimension(p, space);

1941:             for(int i = dim-1; i >= 0; --i) {
1942:               array[++j] += v[i+offset];
1943:             }
1944:             offset += dim;
1945:           }
1946:         }
1947:       } else {
1948:         if (orientation >= 0) {
1949:           const int&                          dim  = this->getFiberDimension(p);
1950:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1951:           int                                 cInd = 0;

1953:           for(int i = 0; i < dim; ++i) {
1954:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1955:             array[i] += v[i];
1956:           }
1957:         } else {
1958:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1959:           int                                 offset  = 0;
1960:           int                                 cOffset = 0;
1961:           int                                 j       = 0;

1963:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1964:             const int  dim = this->getFiberDimension(p, space);
1965:             const int tDim = this->getConstrainedFiberDimension(p, space);
1966:             const int sDim = dim - tDim;
1967:             int       cInd = 0;

1969:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1970:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1971:               array[j] += v[k];
1972:             }
1973:             offset  += dim;
1974:             cOffset += dim - tDim;
1975:           }
1976:         }
1977:       }
1978:     };
1979:     // Update the free values on a point
1980:     //   Takes ONLY unconstrained values
1981:     void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1982:       value_type *array = (value_type *) this->restrictPoint(p);
1983:       const int&  cDim  = this->getConstraintDimension(p);

1985:       if (!cDim) {
1986:         if (orientation >= 0) {
1987:           const int& dim = this->getFiberDimension(p);

1989:           for(int i = 0; i < dim; ++i) {
1990:             array[i] = v[i];
1991:           }
1992:         } else {
1993:           int offset = 0;
1994:           int j      = -1;

1996:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1997:             const int& dim = this->getFiberDimension(p, space);

1999:             for(int i = dim-1; i >= 0; --i) {
2000:               array[++j] = v[i+offset];
2001:             }
2002:             offset += dim;
2003:           }
2004:         }
2005:       } else {
2006:         if (orientation >= 0) {
2007:           const int&                          dim  = this->getFiberDimension(p);
2008:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2009:           int                                 cInd = 0;

2011:           for(int i = 0, k = -1; i < dim; ++i) {
2012:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2013:             array[i] = v[++k];
2014:           }
2015:         } else {
2016:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2017:           int                                 offset  = 0;
2018:           int                                 cOffset = 0;
2019:           int                                 j       = 0;

2021:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2022:             const int  dim = this->getFiberDimension(p, space);
2023:             const int tDim = this->getConstrainedFiberDimension(p, space);
2024:             const int sDim = dim - tDim;
2025:             int       cInd = 0;

2027:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2028:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2029:               array[j] = v[--k];
2030:             }
2031:             offset  += dim;
2032:             cOffset += dim - tDim;
2033:           }
2034:         }
2035:       }
2036:     };
2037:     // Update the free values on a point
2038:     //   Takes ONLY unconstrained values
2039:     void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2040:       value_type *array = (value_type *) this->restrictPoint(p);
2041:       const int&  cDim  = this->getConstraintDimension(p);

2043:       if (!cDim) {
2044:         if (orientation >= 0) {
2045:           const int& dim = this->getFiberDimension(p);

2047:           for(int i = 0; i < dim; ++i) {
2048:             array[i] += v[i];
2049:           }
2050:         } else {
2051:           int offset = 0;
2052:           int j      = -1;

2054:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2055:             const int& dim = this->getFiberDimension(p, space);

2057:             for(int i = dim-1; i >= 0; --i) {
2058:               array[++j] += v[i+offset];
2059:             }
2060:             offset += dim;
2061:           }
2062:         }
2063:       } else {
2064:         if (orientation >= 0) {
2065:           const int&                          dim  = this->getFiberDimension(p);
2066:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2067:           int                                 cInd = 0;

2069:           for(int i = 0, k = -1; i < dim; ++i) {
2070:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2071:             array[i] += v[++k];
2072:           }
2073:         } else {
2074:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2075:           int                                 offset  = 0;
2076:           int                                 cOffset = 0;
2077:           int                                 j       = 0;

2079:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2080:             const int  dim = this->getFiberDimension(p, space);
2081:             const int tDim = this->getConstrainedFiberDimension(p, space);
2082:             const int sDim = dim - tDim;
2083:             int       cInd = 0;

2085:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2086:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2087:               array[j] += v[--k];
2088:             }
2089:             offset  += dim;
2090:             cOffset += dim - tDim;
2091:           }
2092:         }
2093:       }
2094:     };
2095:     // Update only the constrained dofs on a point
2096:     //   This takes an array with ONLY bc values
2097:     void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2098:       value_type *array = (value_type *) this->restrictPoint(p);
2099:       const int&  cDim  = this->getConstraintDimension(p);

2101:       if (cDim) {
2102:         if (orientation >= 0) {
2103:           const int&                          dim  = this->getFiberDimension(p);
2104:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2105:           int                                 cInd = 0;

2107:           for(int i = 0; i < dim; ++i) {
2108:             if (cInd == cDim) break;
2109:             if (i == cDof[cInd]) {
2110:               array[i] = v[cInd];
2111:               ++cInd;
2112:             }
2113:           }
2114:         } else {
2115:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2116:           int                                 cOffset = 0;
2117:           int                                 j       = 0;

2119:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2120:             const int  dim = this->getFiberDimension(p, space);
2121:             const int tDim = this->getConstrainedFiberDimension(p, space);
2122:             int       cInd = 0;

2124:             for(int i = 0; i < dim; ++i, ++j) {
2125:               if (cInd < 0) break;
2126:               if (j == cDof[cInd+cOffset]) {
2127:                 array[j] = v[cInd+cOffset];
2128:                 ++cInd;
2129:               }
2130:             }
2131:             cOffset += dim - tDim;
2132:           }
2133:         }
2134:       }
2135:     };
2136:     // Update only the constrained dofs on a point
2137:     //   This takes an array with ALL values, not just BC
2138:     void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2139:       value_type *array = (value_type *) this->restrictPoint(p);
2140:       const int&  cDim  = this->getConstraintDimension(p);

2142:       if (cDim) {
2143:         if (orientation >= 0) {
2144:           const int&                          dim  = this->getFiberDimension(p);
2145:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2146:           int                                 cInd = 0;

2148:           for(int i = 0; i < dim; ++i) {
2149:             if (cInd == cDim) break;
2150:             if (i == cDof[cInd]) {
2151:               array[i] = v[i];
2152:               ++cInd;
2153:             }
2154:           }
2155:         } else {
2156:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2157:           int                                 offset  = 0;
2158:           int                                 cOffset = 0;
2159:           int                                 j       = 0;

2161:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2162:             const int  dim = this->getFiberDimension(p, space);
2163:             const int tDim = this->getConstrainedFiberDimension(p, space);
2164:             int       cInd = 0;

2166:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2167:               if (cInd < 0) break;
2168:               if (j == cDof[cInd+cOffset]) {
2169:                 array[j] = v[k];
2170:                 ++cInd;
2171:               }
2172:             }
2173:             offset  += dim;
2174:             cOffset += dim - tDim;
2175:           }
2176:         }
2177:       }
2178:     };
2179:     // Update all dofs on a point (free and constrained)
2180:     void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
2181:       value_type *array = (value_type *) this->restrictPoint(p);

2183:       if (orientation >= 0) {
2184:         const int& dim = this->getFiberDimension(p);

2186:         for(int i = 0; i < dim; ++i) {
2187:           array[i] = v[i];
2188:         }
2189:       } else {
2190:         int offset = 0;
2191:         int j      = -1;

2193:         for(int space = 0; space < this->getNumSpaces(); ++space) {
2194:           const int& dim = this->getFiberDimension(p, space);

2196:           for(int i = dim-1; i >= 0; --i) {
2197:             array[++j] = v[i+offset];
2198:           }
2199:           offset += dim;
2200:         }
2201:       }
2202:     };
2203:     // Update all dofs on a point (free and constrained)
2204:     void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
2205:       value_type *array = (value_type *) this->restrictPoint(p);

2207:       if (orientation >= 0) {
2208:         const int& dim = this->getFiberDimension(p);

2210:         for(int i = 0; i < dim; ++i) {
2211:           array[i] += v[i];
2212:         }
2213:       } else {
2214:         int offset = 0;
2215:         int j      = -1;

2217:         for(int space = 0; space < this->getNumSpaces(); ++space) {
2218:           const int& dim = this->getFiberDimension(p, space);

2220:           for(int i = dim-1; i >= 0; --i) {
2221:             array[++j] += v[i+offset];
2222:           }
2223:           offset += dim;
2224:         }
2225:       }
2226:     };
2227:   public: // Fibrations
2228:     int getNumSpaces() const {return this->_spaces.size();};
2229:     const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
2230:     const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
2231:     void addSpace() {
2232:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
2233:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
2234:       this->_spaces.push_back(space);
2235:       this->_bcs.push_back(bc);
2236:     };
2237:     int getFiberDimension(const point_type& p, const int space) const {
2238:       return this->_spaces[space]->restrictPoint(p)->prefix;
2239:     };
2240:     void setFiberDimension(const point_type& p, int dim, const int space) {
2241:       const index_type idx(dim, -1);
2242:       this->_spaces[space]->addPoint(p);
2243:       this->_spaces[space]->updatePoint(p, &idx);
2244:     };
2245:     template<typename Sequence>
2246:     void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
2247:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2248:         this->setFiberDimension(*p_iter, dim, space);
2249:       }
2250:     }
2251:     int getConstraintDimension(const point_type& p, const int space) const {
2252:       if (!this->_bcs[space]->hasPoint(p)) return 0;
2253:       return this->_bcs[space]->getFiberDimension(p);
2254:     }
2255:     void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2256:       this->_bcs[space]->setFiberDimension(p, numConstraints);
2257:     }
2258:     void addConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2259:       this->_bcs[space]->addFiberDimension(p, numConstraints);
2260:     }
2261:     int getConstrainedFiberDimension(const point_type& p, const int space) const {
2262:       return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
2263:     }
2264:     // Return the local dofs which are constrained on a point
2265:     const int *getConstraintDof(const point_type& p, const int space) const {
2266:       return this->_bcs[space]->restrictPoint(p);
2267:     }
2268:     // Set the local dofs which are constrained on a point
2269:     void setConstraintDof(const point_type& p, const int dofs[], const int space) {
2270:       this->_bcs[space]->updatePoint(p, dofs);
2271:     }
2272:     // Return the total number of free dofs
2273:     int size(const int space) const {
2274:       const chart_type& points = this->getChart();
2275:       int               size   = 0;

2277:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2278:         size += this->getConstrainedFiberDimension(*p_iter, space);
2279:       }
2280:       return size;
2281:     };
2282:     template<typename OtherSection>
2283:     void copySpaces(const Obj<OtherSection>& section) {
2284:       const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
2285:       const std::vector<Obj<bc_type> >&    bcs    = section->getBCs();

2287:       this->_spaces.clear();
2288:       for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
2289:         this->_spaces.push_back(*s_iter);
2290:       }
2291:       this->_bcs.clear();
2292:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
2293:         this->_bcs.push_back(*b_iter);
2294:       }
2295:     }
2296:     Obj<GeneralSection> getFibration(const int space) const {
2297:       Obj<GeneralSection> field = new GeneralSection(this->comm(), this->debug());
2298: //     Obj<atlas_type> _atlas;
2299: //     std::vector<Obj<atlas_type> > _spaces;
2300: //     Obj<bc_type>    _bc;
2301: //     std::vector<Obj<bc_type> >    _bcs;
2302:       field->addSpace();
2303:       const chart_type& chart = this->getChart();

2305:       // Copy sizes
2306:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2307:         const int fDim = this->getFiberDimension(*c_iter, space);
2308:         const int cDim = this->getConstraintDimension(*c_iter, space);

2310:         if (fDim) {
2311:           field->setFiberDimension(*c_iter, fDim);
2312:           field->setFiberDimension(*c_iter, fDim, 0);
2313:         }
2314:         if (cDim) {
2315:           field->setConstraintDimension(*c_iter, cDim);
2316:           field->setConstraintDimension(*c_iter, cDim, 0);
2317:         }
2318:       }
2319:       field->allocateStorage();
2320:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
2321:       const chart_type& newChart = field->getChart();

2323:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2324:         const int cDim   = field->getConstraintDimension(*c_iter);
2325:         //const int dof[1] = {0};

2327:         if (cDim) {
2328:           field->setConstraintDof(*c_iter, this->getConstraintDof(*c_iter, space));
2329:         }
2330:       }
2331:       // Copy offsets
2332:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2333:         index_type idx;

2335:         idx.prefix = field->getFiberDimension(*c_iter);
2336:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
2337:         for(int s = 0; s < space; ++s) {
2338:           idx.index += this->getFiberDimension(*c_iter, s);
2339:         }
2340:         newAtlas->addPoint(*c_iter);
2341:         newAtlas->updatePoint(*c_iter, &idx);
2342:       }
2343:       field->replaceStorage(this->_array, true, this->getStorageSize());
2344:       field->setAtlas(newAtlas);
2345:       return field;
2346:     };
2347:   public: // Optimization
2348:     void getCustomRestrictAtlas(const int tag, const int *offsets[], const int *indices[]) {
2349:       *offsets = this->_customRestrictAtlas[tag].first.first;
2350:       *indices = this->_customRestrictAtlas[tag].first.second;
2351:     };
2352:     void getCustomUpdateAtlas(const int tag, const int *offsets[], const int *indices[]) {
2353:       *offsets = this->_customUpdateAtlas[tag].first.first;
2354:       *indices = this->_customUpdateAtlas[tag].first.second;
2355:     };
2356:     // This returns the tag assigned to the atlas
2357:     int setCustomAtlas(const int restrictOffsets[], const int restrictIndices[], const int updateOffsets[], const int updateIndices[], bool autoFree = true) {
2358:       this->_customRestrictAtlas.push_back(customAtlas_type(customAtlasInd_type(restrictOffsets, restrictIndices), autoFree));
2359:       this->_customUpdateAtlas.push_back(customAtlas_type(customAtlasInd_type(updateOffsets, updateIndices), autoFree));
2360:       return this->_customUpdateAtlas.size()-1;
2361:     };
2362:     int copyCustomAtlas(const Obj<GeneralSection>& section, const int tag) {
2363:       const int *rOffsets, *rIndices, *uOffsets, *uIndices;

2365:       section->getCustomRestrictAtlas(tag, &rOffsets, &rIndices);
2366:       section->getCustomUpdateAtlas(tag, &uOffsets, &uIndices);
2367:       return this->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices, false);
2368:     };
2369:   public:
2370:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
2371:       ostringstream txt;
2372:       int rank;

2374:       if (comm == MPI_COMM_NULL) {
2375:         comm = this->comm();
2376:         rank = this->commRank();
2377:       } else {
2378:         MPI_Comm_rank(comm, &rank);
2379:       }
2380:       if (name == "") {
2381:         if(rank == 0) {
2382:           txt << "viewing a GeneralSection" << std::endl;
2383:         }
2384:       } else {
2385:         if (rank == 0) {
2386:           txt << "viewing GeneralSection '" << name << "'" << std::endl;
2387:         }
2388:       }
2389:       if (rank == 0) {
2390:         txt << "  Fields: " << this->getNumSpaces() << std::endl;
2391:       }
2392:       const chart_type& chart = this->getChart();

2394:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2395:         const value_type *array = this->restrictPoint(*p_iter);
2396:         const int&        dim   = this->getFiberDimension(*p_iter);

2398:         if (dim != 0) {
2399:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << "  ";
2400:           for(int i = 0; i < dim; i++) {
2401:             txt << " " << array[i];
2402:           }
2403:           const int& dim = this->getConstraintDimension(*p_iter);

2405:           if (dim) {
2406:             const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);

2408:             txt << " constrained";
2409:             for(int i = 0; i < dim; ++i) {
2410:               txt << " " << bcArray[i];
2411:             }
2412:           }
2413:           txt << std::endl;
2414:         }
2415:       }
2416:       if (chart.size() == 0) {
2417:         txt << "[" << this->commRank() << "]: empty" << std::endl;
2418:       }
2419:       PetscSynchronizedPrintf(comm, txt.str().c_str());
2420:       PetscSynchronizedFlush(comm);
2421:     };
2422:   };
2423:   // FEMSection will support vector BC on a subset of unknowns on a point
2424:   //   We make a separate constraint Section to hold the transformation and projection operators
2425:   //   Storage will be contiguous by node, just as in Section
2426:   //     This allows fast restrict(p)
2427:   //     Then update() is accomplished by projecting constrained unknowns
2428:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
2429:            typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
2430:            typename BCAtlas_ = UniformSection<Point_, int, 1, typename Alloc_::template rebind<int>::other>,
2431:            typename BC_ = Section<Point_, double, typename Alloc_::template rebind<double>::other> >
2432:   class FEMSection : public ALE::ParallelObject {
2433:   public:
2434:     typedef Point_                                                  point_type;
2435:     typedef Value_                                                  value_type;
2436:     typedef Alloc_                                                  alloc_type;
2437:     typedef Atlas_                                                  atlas_type;
2438:     typedef BCAtlas_                                                bc_atlas_type;
2439:     typedef BC_                                                     bc_type;
2440:     typedef Point                                                   index_type;
2441:     typedef typename atlas_type::chart_type                         chart_type;
2442:     typedef value_type *                                            values_type;
2443:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
2444:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
2445:     typedef typename alloc_type::template rebind<bc_type>::other    bc_atlas_alloc_type;
2446:     typedef typename bc_atlas_alloc_type::pointer                   bc_atlas_ptr;
2447:     typedef typename alloc_type::template rebind<bc_type>::other    bc_alloc_type;
2448:     typedef typename bc_alloc_type::pointer                         bc_ptr;
2449:   protected:
2450:     Obj<atlas_type>    _atlas;
2451:     Obj<bc_atlas_type> _bc_atlas;
2452:     Obj<bc_type>       _bc;
2453:     values_type        _array;
2454:     bool               _sharedStorage;
2455:     int                _sharedStorageSize;
2456:     alloc_type         _allocator;
2457:     std::vector<Obj<atlas_type> > _spaces;
2458:     std::vector<Obj<bc_type> >    _bcs;
2459:   public:
2460:     FEMSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
2461:       atlas_ptr pAtlas      = atlas_alloc_type(this->_allocator).allocate(1);
2462:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
2463:       this->_atlas          = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
2464:       bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2465:       bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(comm, debug));
2466:       this->_bc_atlas       = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2467:       bc_ptr pBC            = bc_alloc_type(this->_allocator).allocate(1);
2468:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
2469:       this->_bc             = Obj<bc_type>(pBC, sizeof(bc_type));
2470:       this->_array          = NULL;
2471:       this->_sharedStorage  = false;
2472:     };
2473:     FEMSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
2474:       bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2475:       bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(this->comm(), this->debug()));
2476:       this->_bc_atlas       = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2477:       bc_ptr pBC            = bc_alloc_type(this->_allocator).allocate(1);
2478:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(this->comm(), this->debug()));
2479:       this->_bc             = Obj<bc_type>(pBC, sizeof(bc_type));
2480:     };
2481:     ~FEMSection() {
2482:       if (this->_array && !this->_sharedStorage) {
2483:         const int totalSize = this->sizeWithBC();

2485:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2486:         this->_allocator.deallocate(this->_array, totalSize);
2487:         this->_array = NULL;
2488:       }
2489:     };
2490:   public:
2491:     value_type *getRawArray(const int size) {
2492:       // Put in a sentinel value that deallocates the array
2493:       static value_type *array   = NULL;
2494:       static int         maxSize = 0;

2496:       if (size > maxSize) {
2497:         const value_type dummy(0);

2499:         if (array) {
2500:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
2501:           this->_allocator.deallocate(array, maxSize);
2502:         }
2503:         maxSize = size;
2504:         array   = this->_allocator.allocate(maxSize);
2505:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
2506:       };
2507:       return array;
2508:     };
2509:     int getStorageSize() const {
2510:       if (this->_sharedStorage) {
2511:         return this->_sharedStorageSize;
2512:       }
2513:       return this->sizeWithBC();
2514:     };
2515:   public: // Verifiers
2516:     bool hasPoint(const point_type& point) {
2517:       return this->_atlas->hasPoint(point);
2518:     };
2519:   public: // Accessors
2520:     const chart_type& getChart() const {return this->_atlas->getChart();};
2521:     const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
2522:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
2523:     const Obj<bc_type>& getBC() const {return this->_bc;};
2524:     void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
2525:   public: // BC
2526:     // Returns the number of constraints on a point
2527:     int getConstraintDimension(const point_type& p) const {
2528:       if (!this->_bc_atlas->hasPoint(p)) return 0;
2529:       return this->_bc_atlas->restrictPoint(p)[0];
2530:     };
2531:     // Set the number of constraints on a point
2532:     void setConstraintDimension(const point_type& p, const int numConstraints) {
2533:       this->_bc_atlas->updatePoint(p, &numConstraints);
2534:     };
2535:     // Increment the number of constraints on a point
2536:     void addConstraintDimension(const point_type& p, const int numConstraints) {
2537:       this->_bc_atlas->updatePointAdd(p, &numConstraints);
2538:     };
2539:     const int *getConstraintDof(const point_type& p) const {
2540:       return this->_bc->restrictPoint(p);
2541:     }
2542:   public: // Sizes
2543:     void clear() {
2544:       if (!this->_sharedStorage) {
2545:         const int totalSize = this->sizeWithBC();

2547:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2548:         this->_allocator.deallocate(this->_array, totalSize);
2549:       }
2550:       this->_array = NULL;
2551:       this->_atlas->clear();
2552:       this->_bc_atlas->clear();
2553:       this->_bc->clear();
2554:     };
2555:     // Return the total number of dofs on the point (free and constrained)
2556:     int getFiberDimension(const point_type& p) const {
2557:       return this->_atlas->restrictPoint(p)->prefix;
2558:     };
2559:     int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
2560:       return atlas->restrictPoint(p)->prefix;
2561:     };
2562:     // Set the total number of dofs on the point (free and constrained)
2563:     void setFiberDimension(const point_type& p, int dim) {
2564:       const index_type idx(dim, -1);
2565:       this->_atlas->addPoint(p);
2566:       this->_atlas->updatePoint(p, &idx);
2567:     };
2568:     template<typename Sequence>
2569:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
2570:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2571:         this->setFiberDimension(*p_iter, dim);
2572:       }
2573:     }
2574:     void addFiberDimension(const point_type& p, int dim) {
2575:       if (this->_atlas->hasPoint(p)) {
2576:         const index_type values(dim, 0);
2577:         this->_atlas->updateAddPoint(p, &values);
2578:       } else {
2579:         this->setFiberDimension(p, dim);
2580:       }
2581:     };
2582:     // Return the number of constrained dofs on this point
2583:     int getConstrainedFiberDimension(const point_type& p) const {
2584:       return this->getFiberDimension(p) - this->getConstraintDimension(p);
2585:     };
2586:     // Return the total number of free dofs
2587:     int size() const {
2588:       const chart_type& points = this->getChart();
2589:       int               size   = 0;

2591:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2592:         size += this->getConstrainedFiberDimension(*p_iter);
2593:       }
2594:       return size;
2595:     };
2596:     // Return the total number of dofs (free and constrained)
2597:     int sizeWithBC() const {
2598:       const chart_type& points = this->getChart();
2599:       int               size   = 0;

2601:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2602:         size += this->getFiberDimension(*p_iter);
2603:       }
2604:       return size;
2605:     };
2606:   public: // Allocation
2607:     void allocateStorage() {
2608:       const int totalSize = this->sizeWithBC();
2609:       const value_type dummy(0) ;

2611:       this->_array             = this->_allocator.allocate(totalSize);
2612:       this->_sharedStorage     = false;
2613:       this->_sharedStorageSize = 0;
2614:       for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
2615:       this->_bc_atlas->allocatePoint();
2616:       this->_bc->allocatePoint();
2617:     };
2618:     void orderPoints(const Obj<atlas_type>& atlas){
2619:       const chart_type& chart  = this->getChart();
2620:       int               offset = 0;

2622:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2623:         typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
2624:         const int&                      dim = idx.prefix;

2626:         idx.index = offset;
2627:         atlas->updatePoint(*c_iter, &idx);
2628:         offset += dim;
2629:       }
2630:     };
2631:     void allocatePoint() {
2632:       this->orderPoints(this->_atlas);
2633:       this->allocateStorage();
2634:     };
2635:   public: // Restriction and Update
2636:     // Zero entries
2637:     void zero() {
2638:       this->set(0.0);
2639:     };
2640:     void set(const value_type value) {
2641:       //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
2642:       const chart_type& chart = this->getChart();

2644:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2645:         value_type *array = (value_type *) this->restrictPoint(*c_iter);
2646:         const int&  dim   = this->getFiberDimension(*c_iter);
2647:         const int&  cDim  = this->getConstraintDimension(*c_iter);

2649:         if (!cDim) {
2650:           for(int i = 0; i < dim; ++i) {
2651:             array[i] = value;
2652:           }
2653:         } else {
2654:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2655:           int                                 cInd = 0;

2657:           for(int i = 0; i < dim; ++i) {
2658:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2659:             array[i] = value;
2660:           }
2661:         }
2662:       }
2663:     };
2664:     // Add two sections and put the result in a third
2665:     void add(const Obj<FEMSection>& x, const Obj<FEMSection>& y) {
2666:       // Check atlases
2667:       const chart_type& chart = this->getChart();

2669:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2670:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
2671:         const value_type *xArray = x->restrictPoint(*c_iter);
2672:         const value_type *yArray = y->restrictPoint(*c_iter);
2673:         const int&        dim    = this->getFiberDimension(*c_iter);
2674:         const int&        cDim   = this->getConstraintDimension(*c_iter);

2676:         if (!cDim) {
2677:           for(int i = 0; i < dim; ++i) {
2678:             array[i] = xArray[i] + yArray[i];
2679:           }
2680:         } else {
2681:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2682:           int                                 cInd = 0;

2684:           for(int i = 0; i < dim; ++i) {
2685:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2686:             array[i] = xArray[i] + yArray[i];
2687:           }
2688:         }
2689:       }
2690:     };
2691:     // this = this + alpha * x
2692:     void axpy(const value_type alpha, const Obj<FEMSection>& x) {
2693:       // Check atlases
2694:       const chart_type& chart = this->getChart();

2696:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2697:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
2698:         const value_type *xArray = x->restrictPoint(*c_iter);
2699:         const int&        dim    = this->getFiberDimension(*c_iter);
2700:         const int&        cDim   = this->getConstraintDimension(*c_iter);

2702:         if (!cDim) {
2703:           for(int i = 0; i < dim; ++i) {
2704:             array[i] += alpha*xArray[i];
2705:           }
2706:         } else {
2707:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2708:           int                                 cInd = 0;

2710:           for(int i = 0; i < dim; ++i) {
2711:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2712:             array[i] += alpha*xArray[i];
2713:           }
2714:         }
2715:       }
2716:     };
2717:     // Return the free values on a point
2718:     const value_type *restrictSpace() const {
2719:       return this->_array;
2720:     };
2721:     // Return the free values on a point
2722:     const value_type *restrictPoint(const point_type& p) const {
2723:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
2724:     };
2725:     // Update the free values on a point
2726:     //   Takes a full array and ignores constrained values
2727:     void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2728:       value_type *array = (value_type *) this->restrictPoint(p);
2729:       const int&  cDim  = this->getConstraintDimension(p);

2731:       if (!cDim) {
2732:         if (orientation >= 0) {
2733:           const int& dim = this->getFiberDimension(p);

2735:           for(int i = 0; i < dim; ++i) {
2736:             array[i] = v[i];
2737:           }
2738:         } else {
2739:           int offset = 0;
2740:           int j      = -1;

2742:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2743:             const int& dim = this->getFiberDimension(p, space);

2745:             for(int i = dim-1; i >= 0; --i) {
2746:               array[++j] = v[i+offset];
2747:             }
2748:             offset += dim;
2749:           }
2750:         }
2751:       } else {
2752:         if (orientation >= 0) {
2753:           const int&                          dim  = this->getFiberDimension(p);
2754:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2755:           int                                 cInd = 0;

2757:           for(int i = 0; i < dim; ++i) {
2758:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2759:             array[i] = v[i];
2760:           }
2761:         } else {
2762:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2763:           int                                 offset  = 0;
2764:           int                                 cOffset = 0;
2765:           int                                 j       = 0;

2767:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2768:             const int  dim = this->getFiberDimension(p, space);
2769:             const int tDim = this->getConstrainedFiberDimension(p, space);
2770:             const int sDim = dim - tDim;
2771:             int       cInd = 0;

2773:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2774:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2775:               array[j] = v[k];
2776:             }
2777:             offset  += dim;
2778:             cOffset += dim - tDim;
2779:           }
2780:         }
2781:       }
2782:     };
2783:     // Update the free values on a point
2784:     //   Takes a full array and ignores constrained values
2785:     void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2786:       value_type *array = (value_type *) this->restrictPoint(p);
2787:       const int&  cDim  = this->getConstraintDimension(p);

2789:       if (!cDim) {
2790:         if (orientation >= 0) {
2791:           const int& dim = this->getFiberDimension(p);

2793:           for(int i = 0; i < dim; ++i) {
2794:             array[i] += v[i];
2795:           }
2796:         } else {
2797:           int offset = 0;
2798:           int j      = -1;

2800:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2801:             const int& dim = this->getFiberDimension(p, space);

2803:             for(int i = dim-1; i >= 0; --i) {
2804:               array[++j] += v[i+offset];
2805:             }
2806:             offset += dim;
2807:           }
2808:         }
2809:       } else {
2810:         if (orientation >= 0) {
2811:           const int&                          dim  = this->getFiberDimension(p);
2812:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2813:           int                                 cInd = 0;

2815:           for(int i = 0; i < dim; ++i) {
2816:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2817:             array[i] += v[i];
2818:           }
2819:         } else {
2820:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2821:           int                                 offset  = 0;
2822:           int                                 cOffset = 0;
2823:           int                                 j       = 0;

2825:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2826:             const int  dim = this->getFiberDimension(p, space);
2827:             const int tDim = this->getConstrainedFiberDimension(p, space);
2828:             const int sDim = dim - tDim;
2829:             int       cInd = 0;

2831:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2832:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2833:               array[j] += v[k];
2834:             }
2835:             offset  += dim;
2836:             cOffset += dim - tDim;
2837:           }
2838:         }
2839:       }
2840:     };
2841:     // Update the free values on a point
2842:     //   Takes ONLY unconstrained values
2843:     void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2844:       value_type *array = (value_type *) this->restrictPoint(p);
2845:       const int&  cDim  = this->getConstraintDimension(p);

2847:       if (!cDim) {
2848:         if (orientation >= 0) {
2849:           const int& dim = this->getFiberDimension(p);

2851:           for(int i = 0; i < dim; ++i) {
2852:             array[i] = v[i];
2853:           }
2854:         } else {
2855:           int offset = 0;
2856:           int j      = -1;

2858:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2859:             const int& dim = this->getFiberDimension(p, space);

2861:             for(int i = dim-1; i >= 0; --i) {
2862:               array[++j] = v[i+offset];
2863:             }
2864:             offset += dim;
2865:           }
2866:         }
2867:       } else {
2868:         if (orientation >= 0) {
2869:           const int&                          dim  = this->getFiberDimension(p);
2870:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2871:           int                                 cInd = 0;

2873:           for(int i = 0, k = -1; i < dim; ++i) {
2874:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2875:             array[i] = v[++k];
2876:           }
2877:         } else {
2878:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2879:           int                                 offset  = 0;
2880:           int                                 cOffset = 0;
2881:           int                                 j       = 0;

2883:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2884:             const int  dim = this->getFiberDimension(p, space);
2885:             const int tDim = this->getConstrainedFiberDimension(p, space);
2886:             const int sDim = dim - tDim;
2887:             int       cInd = 0;

2889:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2890:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2891:               array[j] = v[--k];
2892:             }
2893:             offset  += dim;
2894:             cOffset += dim - tDim;
2895:           }
2896:         }
2897:       }
2898:     };
2899:     // Update the free values on a point
2900:     //   Takes ONLY unconstrained values
2901:     void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2902:       value_type *array = (value_type *) this->restrictPoint(p);
2903:       const int&  cDim  = this->getConstraintDimension(p);

2905:       if (!cDim) {
2906:         if (orientation >= 0) {
2907:           const int& dim = this->getFiberDimension(p);

2909:           for(int i = 0; i < dim; ++i) {
2910:             array[i] += v[i];
2911:           }
2912:         } else {
2913:           int offset = 0;
2914:           int j      = -1;

2916:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2917:             const int& dim = this->getFiberDimension(p, space);

2919:             for(int i = dim-1; i >= 0; --i) {
2920:               array[++j] += v[i+offset];
2921:             }
2922:             offset += dim;
2923:           }
2924:         }
2925:       } else {
2926:         if (orientation >= 0) {
2927:           const int&                          dim  = this->getFiberDimension(p);
2928:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2929:           int                                 cInd = 0;

2931:           for(int i = 0, k = -1; i < dim; ++i) {
2932:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2933:             array[i] += v[++k];
2934:           }
2935:         } else {
2936:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2937:           int                                 offset  = 0;
2938:           int                                 cOffset = 0;
2939:           int                                 j       = 0;

2941:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2942:             const int  dim = this->getFiberDimension(p, space);
2943:             const int tDim = this->getConstrainedFiberDimension(p, space);
2944:             const int sDim = dim - tDim;
2945:             int       cInd = 0;

2947:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2948:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2949:               array[j] += v[--k];
2950:             }
2951:             offset  += dim;
2952:             cOffset += dim - tDim;
2953:           }
2954:         }
2955:       }
2956:     };
2957:     // Update only the constrained dofs on a point
2958:     //   This takes an array with ONLY bc values
2959:     void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2960:       value_type *array = (value_type *) this->restrictPoint(p);
2961:       const int&  cDim  = this->getConstraintDimension(p);

2963:       if (cDim) {
2964:         if (orientation >= 0) {
2965:           const int&                          dim  = this->getFiberDimension(p);
2966:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2967:           int                                 cInd = 0;

2969:           for(int i = 0; i < dim; ++i) {
2970:             if (cInd == cDim) break;
2971:             if (i == cDof[cInd]) {
2972:               array[i] = v[cInd];
2973:               ++cInd;
2974:             }
2975:           }
2976:         } else {
2977:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2978:           int                                 cOffset = 0;
2979:           int                                 j       = 0;

2981:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2982:             const int  dim = this->getFiberDimension(p, space);
2983:             const int tDim = this->getConstrainedFiberDimension(p, space);
2984:             int       cInd = 0;

2986:             for(int i = 0; i < dim; ++i, ++j) {
2987:               if (cInd < 0) break;
2988:               if (j == cDof[cInd+cOffset]) {
2989:                 array[j] = v[cInd+cOffset];
2990:                 ++cInd;
2991:               }
2992:             }
2993:             cOffset += dim - tDim;
2994:           }
2995:         }
2996:       }
2997:     };
2998:     // Update only the constrained dofs on a point
2999:     //   This takes an array with ALL values, not just BC
3000:     void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
3001:       value_type *array = (value_type *) this->restrictPoint(p);
3002:       const int&  cDim  = this->getConstraintDimension(p);

3004:       if (cDim) {
3005:         if (orientation >= 0) {
3006:           const int&                          dim  = this->getFiberDimension(p);
3007:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
3008:           int                                 cInd = 0;

3010:           for(int i = 0; i < dim; ++i) {
3011:             if (cInd == cDim) break;
3012:             if (i == cDof[cInd]) {
3013:               array[i] = v[i];
3014:               ++cInd;
3015:             }
3016:           }
3017:         } else {
3018:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
3019:           int                                 offset  = 0;
3020:           int                                 cOffset = 0;
3021:           int                                 j       = 0;

3023:           for(int space = 0; space < this->getNumSpaces(); ++space) {
3024:             const int  dim = this->getFiberDimension(p, space);
3025:             const int tDim = this->getConstrainedFiberDimension(p, space);
3026:             int       cInd = 0;

3028:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
3029:               if (cInd < 0) break;
3030:               if (j == cDof[cInd+cOffset]) {
3031:                 array[j] = v[k];
3032:                 ++cInd;
3033:               }
3034:             }
3035:             offset  += dim;
3036:             cOffset += dim - tDim;
3037:           }
3038:         }
3039:       }
3040:     };
3041:     // Update all dofs on a point (free and constrained)
3042:     void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
3043:       value_type *array = (value_type *) this->restrictPoint(p);

3045:       if (orientation >= 0) {
3046:         const int& dim = this->getFiberDimension(p);

3048:         for(int i = 0; i < dim; ++i) {
3049:           array[i] = v[i];
3050:         }
3051:       } else {
3052:         int offset = 0;
3053:         int j      = -1;

3055:         for(int space = 0; space < this->getNumSpaces(); ++space) {
3056:           const int& dim = this->getFiberDimension(p, space);

3058:           for(int i = dim-1; i >= 0; --i) {
3059:             array[++j] = v[i+offset];
3060:           }
3061:           offset += dim;
3062:         }
3063:       }
3064:     };
3065:     // Update all dofs on a point (free and constrained)
3066:     void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
3067:       value_type *array = (value_type *) this->restrictPoint(p);

3069:       if (orientation >= 0) {
3070:         const int& dim = this->getFiberDimension(p);

3072:         for(int i = 0; i < dim; ++i) {
3073:           array[i] += v[i];
3074:         }
3075:       } else {
3076:         int offset = 0;
3077:         int j      = -1;

3079:         for(int space = 0; space < this->getNumSpaces(); ++space) {
3080:           const int& dim = this->getFiberDimension(p, space);

3082:           for(int i = dim-1; i >= 0; --i) {
3083:             array[++j] += v[i+offset];
3084:           }
3085:           offset += dim;
3086:         }
3087:       }
3088:     };
3089:   public: // Fibrations
3090:     int getNumSpaces() const {return this->_spaces.size();};
3091:     const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
3092:     const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
3093:     void addSpace() {
3094:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
3095:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
3096:       this->_spaces.push_back(space);
3097:       this->_bcs.push_back(bc);
3098:     };
3099:     int getFiberDimension(const point_type& p, const int space) const {
3100:       return this->_spaces[space]->restrictPoint(p)->prefix;
3101:     };
3102:     void setFiberDimension(const point_type& p, int dim, const int space) {
3103:       const index_type idx(dim, -1);
3104:       this->_spaces[space]->addPoint(p);
3105:       this->_spaces[space]->updatePoint(p, &idx);
3106:     };
3107:     template<typename Sequence>
3108:     void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
3109:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
3110:         this->setFiberDimension(*p_iter, dim, space);
3111:       }
3112:     }
3113:     int getConstraintDimension(const point_type& p, const int space) const {
3114:       if (!this->_bcs[space]->hasPoint(p)) return 0;
3115:       return this->_bcs[space]->getFiberDimension(p);
3116:     };
3117:     void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
3118:       this->_bcs[space]->setFiberDimension(p, numConstraints);
3119:     };
3120:     int getConstrainedFiberDimension(const point_type& p, const int space) const {
3121:       return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
3122:     };
3123:     void copyFibration(const Obj<FEMSection>& section) {
3124:       const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
3125:       const std::vector<Obj<bc_type> >&    bcs    = section->getBCs();

3127:       this->_spaces.clear();
3128:       for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
3129:         this->_spaces.push_back(*s_iter);
3130:       }
3131:       this->_bcs.clear();
3132:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
3133:         this->_bcs.push_back(*b_iter);
3134:       }
3135:     };
3136:     Obj<FEMSection> getFibration(const int space) const {
3137:       Obj<FEMSection> field = new FEMSection(this->comm(), this->debug());
3138: //     Obj<atlas_type> _atlas;
3139: //     std::vector<Obj<atlas_type> > _spaces;
3140: //     Obj<bc_type>    _bc;
3141: //     std::vector<Obj<bc_type> >    _bcs;
3142:       field->addSpace();
3143:       const chart_type& chart = this->getChart();

3145:       // Copy sizes
3146:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3147:         const int fDim = this->getFiberDimension(*c_iter, space);
3148:         const int cDim = this->getConstraintDimension(*c_iter, space);

3150:         if (fDim) {
3151:           field->setFiberDimension(*c_iter, fDim);
3152:           field->setFiberDimension(*c_iter, fDim, 0);
3153:         }
3154:         if (cDim) {
3155:           field->setConstraintDimension(*c_iter, cDim);
3156:           field->setConstraintDimension(*c_iter, cDim, 0);
3157:         }
3158:       }
3159:       field->allocateStorage();
3160:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
3161:       const chart_type& newChart = field->getChart();

3163:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3164:         const int cDim   = field->getConstraintDimension(*c_iter);
3165:         const int dof[1] = {0};

3167:         if (cDim) {
3168:           field->setConstraintDof(*c_iter, dof);
3169:         }
3170:       }
3171:       // Copy offsets
3172:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3173:         index_type idx;

3175:         idx.prefix = field->getFiberDimension(*c_iter);
3176:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
3177:         for(int s = 0; s < space; ++s) {
3178:           idx.index += this->getFiberDimension(*c_iter, s);
3179:         }
3180:         newAtlas->addPoint(*c_iter);
3181:         newAtlas->updatePoint(*c_iter, &idx);
3182:       }
3183:       field->replaceStorage(this->_array, true, this->getStorageSize());
3184:       field->setAtlas(newAtlas);
3185:       return field;
3186:     };
3187:   public:
3188:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3189:       ostringstream txt;
3190:       int rank;

3192:       if (comm == MPI_COMM_NULL) {
3193:         comm = this->comm();
3194:         rank = this->commRank();
3195:       } else {
3196:         MPI_Comm_rank(comm, &rank);
3197:       }
3198:       if (name == "") {
3199:         if(rank == 0) {
3200:           txt << "viewing a FEMSection" << std::endl;
3201:         }
3202:       } else {
3203:         if (rank == 0) {
3204:           txt << "viewing FEMSection '" << name << "'" << std::endl;
3205:         }
3206:       }
3207:       if (rank == 0) {
3208:         txt << "  Fields: " << this->getNumSpaces() << std::endl;
3209:       }
3210:       const chart_type& chart = this->getChart();

3212:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
3213:         const value_type *array = this->restrictPoint(*p_iter);
3214:         const int&        dim   = this->getFiberDimension(*p_iter);

3216:         if (dim != 0) {
3217:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << "  ";
3218:           for(int i = 0; i < dim; i++) {
3219:             txt << " " << array[i];
3220:           }
3221:           const int& dim = this->getConstraintDimension(*p_iter);

3223:           if (dim) {
3224:             const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);

3226:             txt << " constrained";
3227:             for(int i = 0; i < dim; ++i) {
3228:               txt << " " << bcArray[i];
3229:             }
3230:           }
3231:           txt << std::endl;
3232:         }
3233:       }
3234:       if (chart.size() == 0) {
3235:         txt << "[" << this->commRank() << "]: empty" << std::endl;
3236:       }
3237:       PetscSynchronizedPrintf(comm, txt.str().c_str());
3238:       PetscSynchronizedFlush(comm);
3239:     };
3240:   };
3241:   // A Field combines several sections
3242:   template<typename Overlap_, typename Patch_, typename Section_>
3243:   class Field : public ALE::ParallelObject {
3244:   public:
3245:     typedef Overlap_                                 overlap_type;
3246:     typedef Patch_                                   patch_type;
3247:     typedef Section_                                 section_type;
3248:     typedef typename section_type::point_type        point_type;
3249:     typedef typename section_type::chart_type        chart_type;
3250:     typedef typename section_type::value_type        value_type;
3251:     typedef std::map<patch_type, Obj<section_type> > sheaf_type;
3252:     typedef enum {SEND, RECEIVE}                     request_type;
3253:     typedef std::map<patch_type, MPI_Request>        requests_type;
3254:   protected:
3255:     sheaf_type    _sheaf;
3256:     int           _tag;
3257:     MPI_Datatype  _datatype;
3258:     requests_type _requests;
3259:   public:
3260:     Field(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
3261:       this->_tag      = this->getNewTag();
3262:       this->_datatype = this->getMPIDatatype();
3263:     };
3264:     Field(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _tag(tag) {
3265:       this->_datatype = this->getMPIDatatype();
3266:     };
3267:     virtual ~Field() {};
3268:   protected:
3269:     MPI_Datatype getMPIDatatype() {
3270:       if (sizeof(value_type) == 4) {
3271:         return MPI_INT;
3272:       } else if (sizeof(value_type) == 8) {
3273:         return MPI_DOUBLE;
3274:       } else if (sizeof(value_type) == 28) {
3275:         int          blen[2];
3276:         MPI_Aint     indices[2];
3277:         MPI_Datatype oldtypes[2], newtype;
3278:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
3279:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3280:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3281:         MPI_Type_commit(&newtype);
3282:         return newtype;
3283:       } else if (sizeof(value_type) == 32) {
3284:         int          blen[2];
3285:         MPI_Aint     indices[2];
3286:         MPI_Datatype oldtypes[2], newtype;
3287:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_DOUBLE;
3288:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3289:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3290:         MPI_Type_commit(&newtype);
3291:         return newtype;
3292:       }
3293:       throw ALE::Exception("Cannot determine MPI type for value type");
3294:     };
3295:     int getNewTag() {
3296:       static int tagKeyval = MPI_KEYVAL_INVALID;
3297:       int *tagvalp = NULL, *maxval, flg;

3299:       if (tagKeyval == MPI_KEYVAL_INVALID) {
3300:         tagvalp = (int *) malloc(sizeof(int));
3301:         MPI_Keyval_create(MPI_NULL_COPY_FN, DMMesh_DelTag, &tagKeyval, (void *) NULL);
3302:         MPI_Attr_put(this->_comm, tagKeyval, tagvalp);
3303:         tagvalp[0] = 0;
3304:       }
3305:       MPI_Attr_get(this->_comm, tagKeyval, (void **) &tagvalp, &flg);
3306:       if (tagvalp[0] < 1) {
3307:         MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
3308:         tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
3309:       }
3310:       if (this->debug()) {
3311:         std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
3312:       }
3313:       return tagvalp[0]--;
3314:     };
3315:   public: // Verifiers
3316:     void checkPatch(const patch_type& patch) const {
3317:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3318:         ostringstream msg;
3319:         msg << "Invalid field patch " << patch << std::endl;
3320:         throw ALE::Exception(msg.str().c_str());
3321:       }
3322:     };
3323:     bool hasPatch(const patch_type& patch) {
3324:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3325:         return false;
3326:       }
3327:       return true;
3328:     };
3329:   public: // Accessors
3330:     int getTag() const {return this->_tag;};
3331:     void setTag(const int tag) {this->_tag = tag;};
3332:     Obj<section_type>& getSection(const patch_type& patch) {
3333:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3334:         this->_sheaf[patch] = new section_type(this->comm(), this->debug());
3335:       }
3336:       return this->_sheaf[patch];
3337:     };
3338:     void setSection(const patch_type& patch, const Obj<section_type>& section) {this->_sheaf[patch] = section;};
3339:     const sheaf_type& getPatches() {
3340:       return this->_sheaf;
3341:     };
3342:     void clear() {
3343:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3344:         p_iter->second->clear();
3345:       }
3346:     };
3347:   public: //  Adapter
3348:     template<typename Topology_>
3349:     void setTopology(const Obj<Topology_>& topology) {
3350:       const typename Topology_::sheaf_type& patches = topology->getPatches();

3352:       for(typename Topology_::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3353:         int                      rank    = p_iter->first;
3354:         const Obj<section_type>& section = this->getSection(rank);
3355:         const Obj<typename Topology_::sieve_type::baseSequence>& base = p_iter->second->base();

3357:         for(typename Topology_::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
3358:           section->setFiberDimension(*b_iter, 1);
3359:         }
3360:       }
3361:     }
3362:     void allocate() {
3363:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3364:         p_iter->second->allocatePoint();
3365:       }
3366:     }
3367:   public: // Communication
3368:     void construct(const int size) {
3369:       const sheaf_type& patches = this->getPatches();

3371:       for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3372:         const patch_type         rank    = p_iter->first;
3373:         const Obj<section_type>& section = this->getSection(rank);
3374:         const chart_type&        chart   = section->getChart();

3376:         for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3377:           section->setFiberDimension(*c_iter, size);
3378:         }
3379:       }
3380:     };
3381:     template<typename Sizer>
3382:     void construct(const Obj<Sizer>& sizer) {
3383:       const sheaf_type& patches = this->getPatches();

3385:       for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3386:         const patch_type         rank    = p_iter->first;
3387:         const Obj<section_type>& section = this->getSection(rank);
3388:         const chart_type&        chart   = section->getChart();
3389: 
3390:         for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3391:           section->setFiberDimension(*c_iter, *(sizer->getSection(rank)->restrictPoint(*c_iter)));
3392:         }
3393:       }
3394:     }
3395:     void constructCommunication(const request_type& requestType) {
3396:       const sheaf_type& patches = this->getPatches();

3398:       for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3399:         const patch_type         patch   = p_iter->first;
3400:         const Obj<section_type>& section = this->getSection(patch);
3401:         MPI_Request              request;

3403:         if (requestType == RECEIVE) {
3404:           if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data(" << section->size() << ") from " << patch << " tag " << this->_tag << std::endl;}
3405:           MPI_Recv_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3406:         } else {
3407:           if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << section->size() << ") to " << patch << " tag " << this->_tag << std::endl;}
3408:           MPI_Send_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3409:         }
3410:         this->_requests[patch] = request;
3411:       }
3412:     };
3413:     void startCommunication() {
3414:       const sheaf_type& patches = this->getPatches();

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

3419:         MPI_Start(&request);
3420:       }
3421:     };
3422:     void endCommunication() {
3423:       const sheaf_type& patches = this->getPatches();
3424:       MPI_Status status;

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

3429:         MPI_Wait(&request, &status);
3430:       }
3431:     };
3432:   public:
3433:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3434:       ostringstream txt;
3435:       int rank;

3437:       if (comm == MPI_COMM_NULL) {
3438:         comm = this->comm();
3439:         rank = this->commRank();
3440:       } else {
3441:         MPI_Comm_rank(comm, &rank);
3442:       }
3443:       if (name == "") {
3444:         if(rank == 0) {
3445:           txt << "viewing a Field" << std::endl;
3446:         }
3447:       } else {
3448:         if(rank == 0) {
3449:           txt << "viewing Field '" << name << "'" << std::endl;
3450:         }
3451:       }
3452:       PetscSynchronizedPrintf(comm, txt.str().c_str());
3453:       PetscSynchronizedFlush(comm);
3454:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3455:         ostringstream txt1;

3457:         txt1 << "[" << this->commRank() << "]: Patch " << p_iter->first << std::endl;
3458:         PetscSynchronizedPrintf(comm, txt1.str().c_str());
3459:         PetscSynchronizedFlush(comm);
3460:         p_iter->second->view("field section", comm);
3461:       }
3462:     };
3463:   };
3464: }
3465: #endif