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 <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:   public: // Verifiers
1457:     bool hasPoint(const point_type& point) {
1458:       return this->_atlas->hasPoint(point);
1459:     };
1460:   public: // Accessors
1461:     const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
1462:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1463:     const Obj<atlas_type>& getNewAtlas() const {return this->_atlasNew;};
1464:     void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1465:     const Obj<bc_type>& getBC() const {return this->_bc;};
1466:     void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
1467:     const chart_type& getChart() const {return this->_atlas->getChart();};
1468:     void setChart(const chart_type& chart) {throw ALE::Exception("setChart() for GeneralSection is invalid");};
1469:   public: // BC
1470:     // Returns the number of constraints on a point
1471:     int getConstraintDimension(const point_type& p) const {
1472:       if (!this->_bc->hasPoint(p)) return 0;
1473:       return this->_bc->getFiberDimension(p);
1474:     }
1475:     // Set the number of constraints on a point
1476:     void setConstraintDimension(const point_type& p, const int numConstraints) {
1477:       this->_bc->setFiberDimension(p, numConstraints);
1478:     }
1479:     // Increment the number of constraints on a point
1480:     void addConstraintDimension(const point_type& p, const int numConstraints) {
1481:       this->_bc->addFiberDimension(p, numConstraints);
1482:     }
1483:     // Return the local dofs which are constrained on a point
1484:     const int *getConstraintDof(const point_type& p) const {
1485:       return this->_bc->restrictPoint(p);
1486:     }
1487:     // Set the local dofs which are constrained on a point
1488:     void setConstraintDof(const point_type& p, const int dofs[]) {
1489:       this->_bc->updatePoint(p, dofs);
1490:     }
1491:     template<typename OtherSection>
1492:     void copyBC(const Obj<OtherSection>& section) {
1493:       this->setBC(section->getBC());
1494:       const chart_type& chart = this->getChart();

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

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

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

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

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

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

1607:         for(int i = start; i < end; ++i) {
1608:           indices[(*indx)++] = i;
1609:         }
1610:       } else {
1611:         const int numSpaces = this->getNumSpaces();
1612:         int offset = start;

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

1617:           for(int i = offset+dim-1; i >= offset; --i) {
1618:             indices[(*indx)++] = i;
1619:           }
1620:           offset += dim;
1621:         }
1622:         if (!numSpaces) {
1623:           const int& dim = this->getFiberDimension(p);

1625:           for(int i = offset+dim-1; i >= offset; --i) {
1626:             indices[(*indx)++] = i;
1627:           }
1628:           offset += dim;
1629:         }
1630:       }
1631:     }
1632:     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) {
1633:       const int& cDim = this->getConstraintDimension(p);

1635:       if (!cDim) {
1636:         this->getIndicesRaw(p, start, indices, indx, orientation);
1637:       } else {
1638:         if (orientation >= 0) {
1639:           const int&                          dim  = this->getFiberDimension(p);
1640:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1641:           int                                 cInd = 0;

1643:           for(int i = start, k = 0; k < dim; ++k) {
1644:             if ((cInd < cDim) && (k == cDof[cInd])) {
1645:               if (!freeOnly) indices[(*indx)++] = -(k+1);
1646:               if (skipConstraints) ++i;
1647:               ++cInd;
1648:             } else {
1649:               indices[(*indx)++] = i++;
1650:             }
1651:           }
1652:         } else {
1653:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1654:           int                                 offset  = 0;
1655:           int                                 cOffset = 0;
1656:           int                                 j       = -1;

1658:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1659:             const int  dim = this->getFiberDimension(p, space);
1660:             const int tDim = this->getConstrainedFiberDimension(p, space);
1661:             int       cInd = (dim - tDim)-1;

1663:             j += dim;
1664:             for(int i = 0, k = start+tDim+offset; i < dim; ++i, --j) {
1665:               if ((cInd >= 0) && (j == cDof[cInd+cOffset])) {
1666:                 if (!freeOnly) indices[(*indx)++] = -(offset+i+1);
1667:                 if (skipConstraints) --k;
1668:                 --cInd;
1669:               } else {
1670:                 indices[(*indx)++] = --k;
1671:               }
1672:             }
1673:             j       += dim;
1674:             offset  += dim;
1675:             cOffset += dim - tDim;
1676:           }
1677:         }
1678:       }
1679:     };
1680:   public: // Allocation
1681:     void allocateStorage() {
1682:       const int totalSize = this->sizeWithBC();
1683:       const value_type dummy(0) ;

1685:       this->_array             = this->_allocator.allocate(totalSize);
1686:       ///this->_array             = new value_type[totalSize];
1687:       this->_sharedStorage     = false;
1688:       this->_sharedStorageSize = 0;
1689:       for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1690:       ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1691:       this->_bc->allocatePoint();
1692:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = this->_bcs.begin(); b_iter != this->_bcs.end(); ++b_iter) {
1693:         (*b_iter)->allocatePoint();;
1694:       }
1695:     };
1696:     void replaceStorage(value_type *newArray, const bool sharedStorage = false, const int sharedStorageSize = 0) {
1697:       if (this->_array && !this->_sharedStorage) {
1698:         const int totalSize = this->sizeWithBC();

1700:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1701:         this->_allocator.deallocate(this->_array, totalSize);
1702:         ///delete [] this->_array;
1703:       }
1704:       this->_array             = newArray;
1705:       this->_sharedStorage     = sharedStorage;
1706:       this->_sharedStorageSize = sharedStorageSize;
1707:       this->_atlas             = this->_atlasNew;
1708:       this->_atlasNew          = NULL;
1709:     };
1710:     void addPoint(const point_type& point, const int dim) {
1711:       if (dim == 0) return;
1712:       if (this->_atlasNew.isNull()) {
1713:         this->_atlasNew = new atlas_type(this->comm(), this->debug());
1714:         this->_atlasNew->copy(this->_atlas);
1715:       }
1716:       const index_type idx(dim, -1);
1717:       this->_atlasNew->addPoint(point);
1718:       this->_atlasNew->updatePoint(point, &idx);
1719:     };
1720:     void orderPoints(const Obj<atlas_type>& atlas){
1721:       const chart_type& chart  = this->getChart();
1722:       int               offset = 0;

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

1728:         idx.index = offset;
1729:         atlas->updatePoint(*c_iter, &idx);
1730:         offset += dim;
1731:       }
1732:     };
1733:     void allocatePoint() {
1734:       this->orderPoints(this->_atlas);
1735:       this->allocateStorage();
1736:     };
1737:   public: // Restriction and Update
1738:     // Zero entries
1739:     void zero() {
1740:       this->set(0.0);
1741:     };
1742:     void zeroWithBC() {
1743:       memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1744:     };
1745:     void set(const value_type value) {
1746:       //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1747:       const chart_type& chart = this->getChart();

1749:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1750:         value_type *array = (value_type *) this->restrictPoint(*c_iter);
1751:         const int&  dim   = this->getFiberDimension(*c_iter);
1752:         const int&  cDim  = this->getConstraintDimension(*c_iter);

1754:         if (!cDim) {
1755:           for(int i = 0; i < dim; ++i) {
1756:             array[i] = value;
1757:           }
1758:         } else {
1759:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1760:           int                                 cInd = 0;

1762:           for(int i = 0; i < dim; ++i) {
1763:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1764:             array[i] = value;
1765:           }
1766:         }
1767:       }
1768:     };
1769:     // Add two sections and put the result in a third
1770:     void add(const Obj<GeneralSection>& x, const Obj<GeneralSection>& y) {
1771:       // Check atlases
1772:       const chart_type& chart = this->getChart();

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

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

1789:           for(int i = 0; i < dim; ++i) {
1790:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1791:             array[i] = xArray[i] + yArray[i];
1792:           }
1793:         }
1794:       }
1795:     };
1796:     // this = this + alpha * x
1797:     template<typename OtherSection>
1798:     void axpy(const value_type alpha, const Obj<OtherSection>& x) {
1799:       // Check atlases
1800:       const chart_type& chart = this->getChart();

1802:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1803:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
1804:         const value_type *xArray = x->restrictPoint(*c_iter);
1805:         const int&        dim    = this->getFiberDimension(*c_iter);
1806:         const int&        cDim   = this->getConstraintDimension(*c_iter);

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

1816:           for(int i = 0; i < dim; ++i) {
1817:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1818:             array[i] += alpha*xArray[i];
1819:           }
1820:         }
1821:       }
1822:     }
1823:     // Return the free values on a point
1824:     const value_type *restrictSpace() const {
1825:       return this->_array;
1826:     }
1827:     // Return the free values on a point
1828:     const value_type *restrictPoint(const point_type& p) const {
1829:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1830:     }
1831:     void restrictPoint(const point_type& p, value_type values[], const int size) const {
1832:       assert(this->_atlas->restrictPoint(p)[0].prefix == size);
1833:       const value_type *v = &(this->_array[this->_atlas->restrictPoint(p)[0].index]);

1835:       for(int i = 0; i < size; ++i) {
1836:         values[i] = v[i];
1837:       }
1838:     };
1839:     // Update the free values on a point
1840:     //   Takes a full array and ignores constrained values
1841:     void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1842:       value_type *array = (value_type *) this->restrictPoint(p);
1843:       const int&  cDim  = this->getConstraintDimension(p);

1845:       if (!cDim) {
1846:         if (orientation >= 0) {
1847:           const int& dim = this->getFiberDimension(p);

1849:           for(int i = 0; i < dim; ++i) {
1850:             array[i] = v[i];
1851:           }
1852:         } else {
1853:           int offset = 0;
1854:           int j      = -1;

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

1859:             for(int i = dim-1; i >= 0; --i) {
1860:               array[++j] = v[i+offset];
1861:             }
1862:             offset += dim;
1863:           }
1864:         }
1865:       } else {
1866:         if (orientation >= 0) {
1867:           const int&                          dim  = this->getFiberDimension(p);
1868:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1869:           int                                 cInd = 0;

1871:           for(int i = 0; i < dim; ++i) {
1872:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1873:             array[i] = v[i];
1874:           }
1875:         } else {
1876:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1877:           int                                 offset  = 0;
1878:           int                                 cOffset = 0;
1879:           int                                 j       = 0;

1881:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1882:             const int  dim = this->getFiberDimension(p, space);
1883:             const int tDim = this->getConstrainedFiberDimension(p, space);
1884:             const int sDim = dim - tDim;
1885:             int       cInd = 0;

1887:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1888:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1889:               array[j] = v[k];
1890:             }
1891:             offset  += dim;
1892:             cOffset += dim - tDim;
1893:           }
1894:         }
1895:       }
1896:     };
1897:     // Update the free values on a point
1898:     //   Takes a full array and ignores constrained values
1899:     void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1900:       value_type *array = (value_type *) this->restrictPoint(p);
1901:       const int&  cDim  = this->getConstraintDimension(p);

1903:       if (!cDim) {
1904:         if (orientation >= 0) {
1905:           const int& dim = this->getFiberDimension(p);

1907:           for(int i = 0; i < dim; ++i) {
1908:             array[i] += v[i];
1909:           }
1910:         } else {
1911:           int offset = 0;
1912:           int j      = -1;

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

1917:             for(int i = dim-1; i >= 0; --i) {
1918:               array[++j] += v[i+offset];
1919:             }
1920:             offset += dim;
1921:           }
1922:         }
1923:       } else {
1924:         if (orientation >= 0) {
1925:           const int&                          dim  = this->getFiberDimension(p);
1926:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1927:           int                                 cInd = 0;

1929:           for(int i = 0; i < dim; ++i) {
1930:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1931:             array[i] += v[i];
1932:           }
1933:         } else {
1934:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1935:           int                                 offset  = 0;
1936:           int                                 cOffset = 0;
1937:           int                                 j       = 0;

1939:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1940:             const int  dim = this->getFiberDimension(p, space);
1941:             const int tDim = this->getConstrainedFiberDimension(p, space);
1942:             const int sDim = dim - tDim;
1943:             int       cInd = 0;

1945:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1946:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1947:               array[j] += v[k];
1948:             }
1949:             offset  += dim;
1950:             cOffset += dim - tDim;
1951:           }
1952:         }
1953:       }
1954:     };
1955:     // Update the free values on a point
1956:     //   Takes ONLY unconstrained values
1957:     void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1958:       value_type *array = (value_type *) this->restrictPoint(p);
1959:       const int&  cDim  = this->getConstraintDimension(p);

1961:       if (!cDim) {
1962:         if (orientation >= 0) {
1963:           const int& dim = this->getFiberDimension(p);

1965:           for(int i = 0; i < dim; ++i) {
1966:             array[i] = v[i];
1967:           }
1968:         } else {
1969:           int offset = 0;
1970:           int j      = -1;

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

1975:             for(int i = dim-1; i >= 0; --i) {
1976:               array[++j] = v[i+offset];
1977:             }
1978:             offset += dim;
1979:           }
1980:         }
1981:       } else {
1982:         if (orientation >= 0) {
1983:           const int&                          dim  = this->getFiberDimension(p);
1984:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1985:           int                                 cInd = 0;

1987:           for(int i = 0, k = -1; i < dim; ++i) {
1988:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1989:             array[i] = v[++k];
1990:           }
1991:         } else {
1992:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1993:           int                                 offset  = 0;
1994:           int                                 cOffset = 0;
1995:           int                                 j       = 0;

1997:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1998:             const int  dim = this->getFiberDimension(p, space);
1999:             const int tDim = this->getConstrainedFiberDimension(p, space);
2000:             const int sDim = dim - tDim;
2001:             int       cInd = 0;

2003:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2004:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2005:               array[j] = v[--k];
2006:             }
2007:             offset  += dim;
2008:             cOffset += dim - tDim;
2009:           }
2010:         }
2011:       }
2012:     };
2013:     // Update the free values on a point
2014:     //   Takes ONLY unconstrained values
2015:     void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2016:       value_type *array = (value_type *) this->restrictPoint(p);
2017:       const int&  cDim  = this->getConstraintDimension(p);

2019:       if (!cDim) {
2020:         if (orientation >= 0) {
2021:           const int& dim = this->getFiberDimension(p);

2023:           for(int i = 0; i < dim; ++i) {
2024:             array[i] += v[i];
2025:           }
2026:         } else {
2027:           int offset = 0;
2028:           int j      = -1;

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

2033:             for(int i = dim-1; i >= 0; --i) {
2034:               array[++j] += v[i+offset];
2035:             }
2036:             offset += dim;
2037:           }
2038:         }
2039:       } else {
2040:         if (orientation >= 0) {
2041:           const int&                          dim  = this->getFiberDimension(p);
2042:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2043:           int                                 cInd = 0;

2045:           for(int i = 0, k = -1; i < dim; ++i) {
2046:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2047:             array[i] += v[++k];
2048:           }
2049:         } else {
2050:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2051:           int                                 offset  = 0;
2052:           int                                 cOffset = 0;
2053:           int                                 j       = 0;

2055:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2056:             const int  dim = this->getFiberDimension(p, space);
2057:             const int tDim = this->getConstrainedFiberDimension(p, space);
2058:             const int sDim = dim - tDim;
2059:             int       cInd = 0;

2061:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2062:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2063:               array[j] += v[--k];
2064:             }
2065:             offset  += dim;
2066:             cOffset += dim - tDim;
2067:           }
2068:         }
2069:       }
2070:     };
2071:     // Update only the constrained dofs on a point
2072:     //   This takes an array with ONLY bc values
2073:     void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2074:       value_type *array = (value_type *) this->restrictPoint(p);
2075:       const int&  cDim  = this->getConstraintDimension(p);

2077:       if (cDim) {
2078:         if (orientation >= 0) {
2079:           const int&                          dim  = this->getFiberDimension(p);
2080:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2081:           int                                 cInd = 0;

2083:           for(int i = 0; i < dim; ++i) {
2084:             if (cInd == cDim) break;
2085:             if (i == cDof[cInd]) {
2086:               array[i] = v[cInd];
2087:               ++cInd;
2088:             }
2089:           }
2090:         } else {
2091:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2092:           int                                 cOffset = 0;
2093:           int                                 j       = 0;

2095:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2096:             const int  dim = this->getFiberDimension(p, space);
2097:             const int tDim = this->getConstrainedFiberDimension(p, space);
2098:             int       cInd = 0;

2100:             for(int i = 0; i < dim; ++i, ++j) {
2101:               if (cInd < 0) break;
2102:               if (j == cDof[cInd+cOffset]) {
2103:                 array[j] = v[cInd+cOffset];
2104:                 ++cInd;
2105:               }
2106:             }
2107:             cOffset += dim - tDim;
2108:           }
2109:         }
2110:       }
2111:     };
2112:     // Update only the constrained dofs on a point
2113:     //   This takes an array with ALL values, not just BC
2114:     void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2115:       value_type *array = (value_type *) this->restrictPoint(p);
2116:       const int&  cDim  = this->getConstraintDimension(p);

2118:       if (cDim) {
2119:         if (orientation >= 0) {
2120:           const int&                          dim  = this->getFiberDimension(p);
2121:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2122:           int                                 cInd = 0;

2124:           for(int i = 0; i < dim; ++i) {
2125:             if (cInd == cDim) break;
2126:             if (i == cDof[cInd]) {
2127:               array[i] = v[i];
2128:               ++cInd;
2129:             }
2130:           }
2131:         } else {
2132:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2133:           int                                 offset  = 0;
2134:           int                                 cOffset = 0;
2135:           int                                 j       = 0;

2137:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2138:             const int  dim = this->getFiberDimension(p, space);
2139:             const int tDim = this->getConstrainedFiberDimension(p, space);
2140:             int       cInd = 0;

2142:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2143:               if (cInd < 0) break;
2144:               if (j == cDof[cInd+cOffset]) {
2145:                 array[j] = v[k];
2146:                 ++cInd;
2147:               }
2148:             }
2149:             offset  += dim;
2150:             cOffset += dim - tDim;
2151:           }
2152:         }
2153:       }
2154:     };
2155:     // Update all dofs on a point (free and constrained)
2156:     void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
2157:       value_type *array = (value_type *) this->restrictPoint(p);

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

2162:         for(int i = 0; i < dim; ++i) {
2163:           array[i] = v[i];
2164:         }
2165:       } else {
2166:         int offset = 0;
2167:         int j      = -1;

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

2172:           for(int i = dim-1; i >= 0; --i) {
2173:             array[++j] = v[i+offset];
2174:           }
2175:           offset += dim;
2176:         }
2177:       }
2178:     };
2179:     // Update all dofs on a point (free and constrained)
2180:     void updatePointAllAdd(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:   public: // Fibrations
2204:     int getNumSpaces() const {return this->_spaces.size();};
2205:     const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
2206:     const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
2207:     void addSpace() {
2208:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
2209:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
2210:       this->_spaces.push_back(space);
2211:       this->_bcs.push_back(bc);
2212:     };
2213:     int getFiberDimension(const point_type& p, const int space) const {
2214:       return this->_spaces[space]->restrictPoint(p)->prefix;
2215:     };
2216:     void setFiberDimension(const point_type& p, int dim, const int space) {
2217:       const index_type idx(dim, -1);
2218:       this->_spaces[space]->addPoint(p);
2219:       this->_spaces[space]->updatePoint(p, &idx);
2220:     };
2221:     template<typename Sequence>
2222:     void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
2223:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2224:         this->setFiberDimension(*p_iter, dim, space);
2225:       }
2226:     }
2227:     int getConstraintDimension(const point_type& p, const int space) const {
2228:       if (!this->_bcs[space]->hasPoint(p)) return 0;
2229:       return this->_bcs[space]->getFiberDimension(p);
2230:     }
2231:     void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2232:       this->_bcs[space]->setFiberDimension(p, numConstraints);
2233:     }
2234:     void addConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2235:       this->_bcs[space]->addFiberDimension(p, numConstraints);
2236:     }
2237:     int getConstrainedFiberDimension(const point_type& p, const int space) const {
2238:       return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
2239:     }
2240:     // Return the local dofs which are constrained on a point
2241:     const int *getConstraintDof(const point_type& p, const int space) const {
2242:       return this->_bcs[space]->restrictPoint(p);
2243:     }
2244:     // Set the local dofs which are constrained on a point
2245:     void setConstraintDof(const point_type& p, const int dofs[], const int space) {
2246:       this->_bcs[space]->updatePoint(p, dofs);
2247:     }
2248:     // Return the total number of free dofs
2249:     int size(const int space) const {
2250:       const chart_type& points = this->getChart();
2251:       int               size   = 0;

2253:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2254:         size += this->getConstrainedFiberDimension(*p_iter, space);
2255:       }
2256:       return size;
2257:     };
2258:     template<typename OtherSection>
2259:     void copyFibration(const Obj<OtherSection>& section) {
2260:       const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
2261:       const std::vector<Obj<bc_type> >&    bcs    = section->getBCs();

2263:       this->_spaces.clear();
2264:       for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
2265:         this->_spaces.push_back(*s_iter);
2266:       }
2267:       this->_bcs.clear();
2268:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
2269:         this->_bcs.push_back(*b_iter);
2270:       }
2271:     }
2272:     Obj<GeneralSection> getFibration(const int space) const {
2273:       Obj<GeneralSection> field = new GeneralSection(this->comm(), this->debug());
2274: //     Obj<atlas_type> _atlas;
2275: //     std::vector<Obj<atlas_type> > _spaces;
2276: //     Obj<bc_type>    _bc;
2277: //     std::vector<Obj<bc_type> >    _bcs;
2278:       field->addSpace();
2279:       const chart_type& chart = this->getChart();

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

2286:         if (fDim) {
2287:           field->setFiberDimension(*c_iter, fDim);
2288:           field->setFiberDimension(*c_iter, fDim, 0);
2289:         }
2290:         if (cDim) {
2291:           field->setConstraintDimension(*c_iter, cDim);
2292:           field->setConstraintDimension(*c_iter, cDim, 0);
2293:         }
2294:       }
2295:       field->allocateStorage();
2296:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
2297:       const chart_type& newChart = field->getChart();

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

2303:         if (cDim) {
2304:           field->setConstraintDof(*c_iter, this->getConstraintDof(*c_iter, space));
2305:         }
2306:       }
2307:       // Copy offsets
2308:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2309:         index_type idx;

2311:         idx.prefix = field->getFiberDimension(*c_iter);
2312:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
2313:         for(int s = 0; s < space; ++s) {
2314:           idx.index += this->getFiberDimension(*c_iter, s);
2315:         }
2316:         newAtlas->addPoint(*c_iter);
2317:         newAtlas->updatePoint(*c_iter, &idx);
2318:       }
2319:       field->replaceStorage(this->_array, true, this->getStorageSize());
2320:       field->setAtlas(newAtlas);
2321:       return field;
2322:     };
2323:   public: // Optimization
2324:     void getCustomRestrictAtlas(const int tag, const int *offsets[], const int *indices[]) {
2325:       *offsets = this->_customRestrictAtlas[tag].first.first;
2326:       *indices = this->_customRestrictAtlas[tag].first.second;
2327:     };
2328:     void getCustomUpdateAtlas(const int tag, const int *offsets[], const int *indices[]) {
2329:       *offsets = this->_customUpdateAtlas[tag].first.first;
2330:       *indices = this->_customUpdateAtlas[tag].first.second;
2331:     };
2332:     // This returns the tag assigned to the atlas
2333:     int setCustomAtlas(const int restrictOffsets[], const int restrictIndices[], const int updateOffsets[], const int updateIndices[], bool autoFree = true) {
2334:       this->_customRestrictAtlas.push_back(customAtlas_type(customAtlasInd_type(restrictOffsets, restrictIndices), autoFree));
2335:       this->_customUpdateAtlas.push_back(customAtlas_type(customAtlasInd_type(updateOffsets, updateIndices), autoFree));
2336:       return this->_customUpdateAtlas.size()-1;
2337:     };
2338:     int copyCustomAtlas(const Obj<GeneralSection>& section, const int tag) {
2339:       const int *rOffsets, *rIndices, *uOffsets, *uIndices;

2341:       section->getCustomRestrictAtlas(tag, &rOffsets, &rIndices);
2342:       section->getCustomUpdateAtlas(tag, &uOffsets, &uIndices);
2343:       return this->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices, false);
2344:     };
2345:   public:
2346:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
2347:       ostringstream txt;
2348:       int rank;

2350:       if (comm == MPI_COMM_NULL) {
2351:         comm = this->comm();
2352:         rank = this->commRank();
2353:       } else {
2354:         MPI_Comm_rank(comm, &rank);
2355:       }
2356:       if (name == "") {
2357:         if(rank == 0) {
2358:           txt << "viewing a GeneralSection" << std::endl;
2359:         }
2360:       } else {
2361:         if (rank == 0) {
2362:           txt << "viewing GeneralSection '" << name << "'" << std::endl;
2363:         }
2364:       }
2365:       if (rank == 0) {
2366:         txt << "  Fields: " << this->getNumSpaces() << std::endl;
2367:       }
2368:       const chart_type& chart = this->getChart();

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

2374:         if (dim != 0) {
2375:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << "  ";
2376:           for(int i = 0; i < dim; i++) {
2377:             txt << " " << array[i];
2378:           }
2379:           const int& dim = this->getConstraintDimension(*p_iter);

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

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

2461:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2462:         this->_allocator.deallocate(this->_array, totalSize);
2463:         this->_array = NULL;
2464:       }
2465:     };
2466:   public:
2467:     value_type *getRawArray(const int size) {
2468:       // Put in a sentinel value that deallocates the array
2469:       static value_type *array   = NULL;
2470:       static int         maxSize = 0;

2472:       if (size > maxSize) {
2473:         const value_type dummy(0);

2475:         if (array) {
2476:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
2477:           this->_allocator.deallocate(array, maxSize);
2478:         }
2479:         maxSize = size;
2480:         array   = this->_allocator.allocate(maxSize);
2481:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
2482:       };
2483:       return array;
2484:     };
2485:     int getStorageSize() const {
2486:       if (this->_sharedStorage) {
2487:         return this->_sharedStorageSize;
2488:       }
2489:       return this->sizeWithBC();
2490:     };
2491:   public: // Verifiers
2492:     bool hasPoint(const point_type& point) {
2493:       return this->_atlas->hasPoint(point);
2494:     };
2495:   public: // Accessors
2496:     const chart_type& getChart() const {return this->_atlas->getChart();};
2497:     const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
2498:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
2499:     const Obj<bc_type>& getBC() const {return this->_bc;};
2500:     void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
2501:   public: // BC
2502:     // Returns the number of constraints on a point
2503:     int getConstraintDimension(const point_type& p) const {
2504:       if (!this->_bc_atlas->hasPoint(p)) return 0;
2505:       return this->_bc_atlas->restrictPoint(p)[0];
2506:     };
2507:     // Set the number of constraints on a point
2508:     void setConstraintDimension(const point_type& p, const int numConstraints) {
2509:       this->_bc_atlas->updatePoint(p, &numConstraints);
2510:     };
2511:     // Increment the number of constraints on a point
2512:     void addConstraintDimension(const point_type& p, const int numConstraints) {
2513:       this->_bc_atlas->updatePointAdd(p, &numConstraints);
2514:     };
2515:     const int *getConstraintDof(const point_type& p) const {
2516:       return this->_bc->restrictPoint(p);
2517:     }
2518:   public: // Sizes
2519:     void clear() {
2520:       if (!this->_sharedStorage) {
2521:         const int totalSize = this->sizeWithBC();

2523:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2524:         this->_allocator.deallocate(this->_array, totalSize);
2525:       }
2526:       this->_array = NULL;
2527:       this->_atlas->clear();
2528:       this->_bc_atlas->clear();
2529:       this->_bc->clear();
2530:     };
2531:     // Return the total number of dofs on the point (free and constrained)
2532:     int getFiberDimension(const point_type& p) const {
2533:       return this->_atlas->restrictPoint(p)->prefix;
2534:     };
2535:     int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
2536:       return atlas->restrictPoint(p)->prefix;
2537:     };
2538:     // Set the total number of dofs on the point (free and constrained)
2539:     void setFiberDimension(const point_type& p, int dim) {
2540:       const index_type idx(dim, -1);
2541:       this->_atlas->addPoint(p);
2542:       this->_atlas->updatePoint(p, &idx);
2543:     };
2544:     template<typename Sequence>
2545:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
2546:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2547:         this->setFiberDimension(*p_iter, dim);
2548:       }
2549:     }
2550:     void addFiberDimension(const point_type& p, int dim) {
2551:       if (this->_atlas->hasPoint(p)) {
2552:         const index_type values(dim, 0);
2553:         this->_atlas->updateAddPoint(p, &values);
2554:       } else {
2555:         this->setFiberDimension(p, dim);
2556:       }
2557:     };
2558:     // Return the number of constrained dofs on this point
2559:     int getConstrainedFiberDimension(const point_type& p) const {
2560:       return this->getFiberDimension(p) - this->getConstraintDimension(p);
2561:     };
2562:     // Return the total number of free dofs
2563:     int size() const {
2564:       const chart_type& points = this->getChart();
2565:       int               size   = 0;

2567:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2568:         size += this->getConstrainedFiberDimension(*p_iter);
2569:       }
2570:       return size;
2571:     };
2572:     // Return the total number of dofs (free and constrained)
2573:     int sizeWithBC() const {
2574:       const chart_type& points = this->getChart();
2575:       int               size   = 0;

2577:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2578:         size += this->getFiberDimension(*p_iter);
2579:       }
2580:       return size;
2581:     };
2582:   public: // Allocation
2583:     void allocateStorage() {
2584:       const int totalSize = this->sizeWithBC();
2585:       const value_type dummy(0) ;

2587:       this->_array             = this->_allocator.allocate(totalSize);
2588:       this->_sharedStorage     = false;
2589:       this->_sharedStorageSize = 0;
2590:       for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
2591:       this->_bc_atlas->allocatePoint();
2592:       this->_bc->allocatePoint();
2593:     };
2594:     void orderPoints(const Obj<atlas_type>& atlas){
2595:       const chart_type& chart  = this->getChart();
2596:       int               offset = 0;

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

2602:         idx.index = offset;
2603:         atlas->updatePoint(*c_iter, &idx);
2604:         offset += dim;
2605:       }
2606:     };
2607:     void allocatePoint() {
2608:       this->orderPoints(this->_atlas);
2609:       this->allocateStorage();
2610:     };
2611:   public: // Restriction and Update
2612:     // Zero entries
2613:     void zero() {
2614:       this->set(0.0);
2615:     };
2616:     void set(const value_type value) {
2617:       //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
2618:       const chart_type& chart = this->getChart();

2620:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2621:         value_type *array = (value_type *) this->restrictPoint(*c_iter);
2622:         const int&  dim   = this->getFiberDimension(*c_iter);
2623:         const int&  cDim  = this->getConstraintDimension(*c_iter);

2625:         if (!cDim) {
2626:           for(int i = 0; i < dim; ++i) {
2627:             array[i] = value;
2628:           }
2629:         } else {
2630:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2631:           int                                 cInd = 0;

2633:           for(int i = 0; i < dim; ++i) {
2634:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2635:             array[i] = value;
2636:           }
2637:         }
2638:       }
2639:     };
2640:     // Add two sections and put the result in a third
2641:     void add(const Obj<FEMSection>& x, const Obj<FEMSection>& y) {
2642:       // Check atlases
2643:       const chart_type& chart = this->getChart();

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

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

2660:           for(int i = 0; i < dim; ++i) {
2661:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2662:             array[i] = xArray[i] + yArray[i];
2663:           }
2664:         }
2665:       }
2666:     };
2667:     // this = this + alpha * x
2668:     void axpy(const value_type alpha, const Obj<FEMSection>& x) {
2669:       // Check atlases
2670:       const chart_type& chart = this->getChart();

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

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

2686:           for(int i = 0; i < dim; ++i) {
2687:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2688:             array[i] += alpha*xArray[i];
2689:           }
2690:         }
2691:       }
2692:     };
2693:     // Return the free values on a point
2694:     const value_type *restrictSpace() const {
2695:       return this->_array;
2696:     };
2697:     // Return the free values on a point
2698:     const value_type *restrictPoint(const point_type& p) const {
2699:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
2700:     };
2701:     // Update the free values on a point
2702:     //   Takes a full array and ignores constrained values
2703:     void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2704:       value_type *array = (value_type *) this->restrictPoint(p);
2705:       const int&  cDim  = this->getConstraintDimension(p);

2707:       if (!cDim) {
2708:         if (orientation >= 0) {
2709:           const int& dim = this->getFiberDimension(p);

2711:           for(int i = 0; i < dim; ++i) {
2712:             array[i] = v[i];
2713:           }
2714:         } else {
2715:           int offset = 0;
2716:           int j      = -1;

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

2721:             for(int i = dim-1; i >= 0; --i) {
2722:               array[++j] = v[i+offset];
2723:             }
2724:             offset += dim;
2725:           }
2726:         }
2727:       } else {
2728:         if (orientation >= 0) {
2729:           const int&                          dim  = this->getFiberDimension(p);
2730:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2731:           int                                 cInd = 0;

2733:           for(int i = 0; i < dim; ++i) {
2734:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2735:             array[i] = v[i];
2736:           }
2737:         } else {
2738:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2739:           int                                 offset  = 0;
2740:           int                                 cOffset = 0;
2741:           int                                 j       = 0;

2743:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2744:             const int  dim = this->getFiberDimension(p, space);
2745:             const int tDim = this->getConstrainedFiberDimension(p, space);
2746:             const int sDim = dim - tDim;
2747:             int       cInd = 0;

2749:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2750:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2751:               array[j] = v[k];
2752:             }
2753:             offset  += dim;
2754:             cOffset += dim - tDim;
2755:           }
2756:         }
2757:       }
2758:     };
2759:     // Update the free values on a point
2760:     //   Takes a full array and ignores constrained values
2761:     void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2762:       value_type *array = (value_type *) this->restrictPoint(p);
2763:       const int&  cDim  = this->getConstraintDimension(p);

2765:       if (!cDim) {
2766:         if (orientation >= 0) {
2767:           const int& dim = this->getFiberDimension(p);

2769:           for(int i = 0; i < dim; ++i) {
2770:             array[i] += v[i];
2771:           }
2772:         } else {
2773:           int offset = 0;
2774:           int j      = -1;

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

2779:             for(int i = dim-1; i >= 0; --i) {
2780:               array[++j] += v[i+offset];
2781:             }
2782:             offset += dim;
2783:           }
2784:         }
2785:       } else {
2786:         if (orientation >= 0) {
2787:           const int&                          dim  = this->getFiberDimension(p);
2788:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2789:           int                                 cInd = 0;

2791:           for(int i = 0; i < dim; ++i) {
2792:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2793:             array[i] += v[i];
2794:           }
2795:         } else {
2796:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2797:           int                                 offset  = 0;
2798:           int                                 cOffset = 0;
2799:           int                                 j       = 0;

2801:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2802:             const int  dim = this->getFiberDimension(p, space);
2803:             const int tDim = this->getConstrainedFiberDimension(p, space);
2804:             const int sDim = dim - tDim;
2805:             int       cInd = 0;

2807:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2808:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2809:               array[j] += v[k];
2810:             }
2811:             offset  += dim;
2812:             cOffset += dim - tDim;
2813:           }
2814:         }
2815:       }
2816:     };
2817:     // Update the free values on a point
2818:     //   Takes ONLY unconstrained values
2819:     void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2820:       value_type *array = (value_type *) this->restrictPoint(p);
2821:       const int&  cDim  = this->getConstraintDimension(p);

2823:       if (!cDim) {
2824:         if (orientation >= 0) {
2825:           const int& dim = this->getFiberDimension(p);

2827:           for(int i = 0; i < dim; ++i) {
2828:             array[i] = v[i];
2829:           }
2830:         } else {
2831:           int offset = 0;
2832:           int j      = -1;

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

2837:             for(int i = dim-1; i >= 0; --i) {
2838:               array[++j] = v[i+offset];
2839:             }
2840:             offset += dim;
2841:           }
2842:         }
2843:       } else {
2844:         if (orientation >= 0) {
2845:           const int&                          dim  = this->getFiberDimension(p);
2846:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2847:           int                                 cInd = 0;

2849:           for(int i = 0, k = -1; i < dim; ++i) {
2850:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2851:             array[i] = v[++k];
2852:           }
2853:         } else {
2854:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2855:           int                                 offset  = 0;
2856:           int                                 cOffset = 0;
2857:           int                                 j       = 0;

2859:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2860:             const int  dim = this->getFiberDimension(p, space);
2861:             const int tDim = this->getConstrainedFiberDimension(p, space);
2862:             const int sDim = dim - tDim;
2863:             int       cInd = 0;

2865:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2866:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2867:               array[j] = v[--k];
2868:             }
2869:             offset  += dim;
2870:             cOffset += dim - tDim;
2871:           }
2872:         }
2873:       }
2874:     };
2875:     // Update the free values on a point
2876:     //   Takes ONLY unconstrained values
2877:     void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2878:       value_type *array = (value_type *) this->restrictPoint(p);
2879:       const int&  cDim  = this->getConstraintDimension(p);

2881:       if (!cDim) {
2882:         if (orientation >= 0) {
2883:           const int& dim = this->getFiberDimension(p);

2885:           for(int i = 0; i < dim; ++i) {
2886:             array[i] += v[i];
2887:           }
2888:         } else {
2889:           int offset = 0;
2890:           int j      = -1;

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

2895:             for(int i = dim-1; i >= 0; --i) {
2896:               array[++j] += v[i+offset];
2897:             }
2898:             offset += dim;
2899:           }
2900:         }
2901:       } else {
2902:         if (orientation >= 0) {
2903:           const int&                          dim  = this->getFiberDimension(p);
2904:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2905:           int                                 cInd = 0;

2907:           for(int i = 0, k = -1; i < dim; ++i) {
2908:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2909:             array[i] += v[++k];
2910:           }
2911:         } else {
2912:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2913:           int                                 offset  = 0;
2914:           int                                 cOffset = 0;
2915:           int                                 j       = 0;

2917:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2918:             const int  dim = this->getFiberDimension(p, space);
2919:             const int tDim = this->getConstrainedFiberDimension(p, space);
2920:             const int sDim = dim - tDim;
2921:             int       cInd = 0;

2923:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2924:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2925:               array[j] += v[--k];
2926:             }
2927:             offset  += dim;
2928:             cOffset += dim - tDim;
2929:           }
2930:         }
2931:       }
2932:     };
2933:     // Update only the constrained dofs on a point
2934:     //   This takes an array with ONLY bc values
2935:     void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2936:       value_type *array = (value_type *) this->restrictPoint(p);
2937:       const int&  cDim  = this->getConstraintDimension(p);

2939:       if (cDim) {
2940:         if (orientation >= 0) {
2941:           const int&                          dim  = this->getFiberDimension(p);
2942:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2943:           int                                 cInd = 0;

2945:           for(int i = 0; i < dim; ++i) {
2946:             if (cInd == cDim) break;
2947:             if (i == cDof[cInd]) {
2948:               array[i] = v[cInd];
2949:               ++cInd;
2950:             }
2951:           }
2952:         } else {
2953:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2954:           int                                 cOffset = 0;
2955:           int                                 j       = 0;

2957:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2958:             const int  dim = this->getFiberDimension(p, space);
2959:             const int tDim = this->getConstrainedFiberDimension(p, space);
2960:             int       cInd = 0;

2962:             for(int i = 0; i < dim; ++i, ++j) {
2963:               if (cInd < 0) break;
2964:               if (j == cDof[cInd+cOffset]) {
2965:                 array[j] = v[cInd+cOffset];
2966:                 ++cInd;
2967:               }
2968:             }
2969:             cOffset += dim - tDim;
2970:           }
2971:         }
2972:       }
2973:     };
2974:     // Update only the constrained dofs on a point
2975:     //   This takes an array with ALL values, not just BC
2976:     void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2977:       value_type *array = (value_type *) this->restrictPoint(p);
2978:       const int&  cDim  = this->getConstraintDimension(p);

2980:       if (cDim) {
2981:         if (orientation >= 0) {
2982:           const int&                          dim  = this->getFiberDimension(p);
2983:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2984:           int                                 cInd = 0;

2986:           for(int i = 0; i < dim; ++i) {
2987:             if (cInd == cDim) break;
2988:             if (i == cDof[cInd]) {
2989:               array[i] = v[i];
2990:               ++cInd;
2991:             }
2992:           }
2993:         } else {
2994:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2995:           int                                 offset  = 0;
2996:           int                                 cOffset = 0;
2997:           int                                 j       = 0;

2999:           for(int space = 0; space < this->getNumSpaces(); ++space) {
3000:             const int  dim = this->getFiberDimension(p, space);
3001:             const int tDim = this->getConstrainedFiberDimension(p, space);
3002:             int       cInd = 0;

3004:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
3005:               if (cInd < 0) break;
3006:               if (j == cDof[cInd+cOffset]) {
3007:                 array[j] = v[k];
3008:                 ++cInd;
3009:               }
3010:             }
3011:             offset  += dim;
3012:             cOffset += dim - tDim;
3013:           }
3014:         }
3015:       }
3016:     };
3017:     // Update all dofs on a point (free and constrained)
3018:     void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
3019:       value_type *array = (value_type *) this->restrictPoint(p);

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

3024:         for(int i = 0; i < dim; ++i) {
3025:           array[i] = v[i];
3026:         }
3027:       } else {
3028:         int offset = 0;
3029:         int j      = -1;

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

3034:           for(int i = dim-1; i >= 0; --i) {
3035:             array[++j] = v[i+offset];
3036:           }
3037:           offset += dim;
3038:         }
3039:       }
3040:     };
3041:     // Update all dofs on a point (free and constrained)
3042:     void updatePointAllAdd(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:   public: // Fibrations
3066:     int getNumSpaces() const {return this->_spaces.size();};
3067:     const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
3068:     const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
3069:     void addSpace() {
3070:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
3071:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
3072:       this->_spaces.push_back(space);
3073:       this->_bcs.push_back(bc);
3074:     };
3075:     int getFiberDimension(const point_type& p, const int space) const {
3076:       return this->_spaces[space]->restrictPoint(p)->prefix;
3077:     };
3078:     void setFiberDimension(const point_type& p, int dim, const int space) {
3079:       const index_type idx(dim, -1);
3080:       this->_spaces[space]->addPoint(p);
3081:       this->_spaces[space]->updatePoint(p, &idx);
3082:     };
3083:     template<typename Sequence>
3084:     void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
3085:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
3086:         this->setFiberDimension(*p_iter, dim, space);
3087:       }
3088:     }
3089:     int getConstraintDimension(const point_type& p, const int space) const {
3090:       if (!this->_bcs[space]->hasPoint(p)) return 0;
3091:       return this->_bcs[space]->getFiberDimension(p);
3092:     };
3093:     void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
3094:       this->_bcs[space]->setFiberDimension(p, numConstraints);
3095:     };
3096:     int getConstrainedFiberDimension(const point_type& p, const int space) const {
3097:       return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
3098:     };
3099:     void copyFibration(const Obj<FEMSection>& section) {
3100:       const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
3101:       const std::vector<Obj<bc_type> >&    bcs    = section->getBCs();

3103:       this->_spaces.clear();
3104:       for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
3105:         this->_spaces.push_back(*s_iter);
3106:       }
3107:       this->_bcs.clear();
3108:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
3109:         this->_bcs.push_back(*b_iter);
3110:       }
3111:     };
3112:     Obj<FEMSection> getFibration(const int space) const {
3113:       Obj<FEMSection> field = new FEMSection(this->comm(), this->debug());
3114: //     Obj<atlas_type> _atlas;
3115: //     std::vector<Obj<atlas_type> > _spaces;
3116: //     Obj<bc_type>    _bc;
3117: //     std::vector<Obj<bc_type> >    _bcs;
3118:       field->addSpace();
3119:       const chart_type& chart = this->getChart();

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

3126:         if (fDim) {
3127:           field->setFiberDimension(*c_iter, fDim);
3128:           field->setFiberDimension(*c_iter, fDim, 0);
3129:         }
3130:         if (cDim) {
3131:           field->setConstraintDimension(*c_iter, cDim);
3132:           field->setConstraintDimension(*c_iter, cDim, 0);
3133:         }
3134:       }
3135:       field->allocateStorage();
3136:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
3137:       const chart_type& newChart = field->getChart();

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

3143:         if (cDim) {
3144:           field->setConstraintDof(*c_iter, dof);
3145:         }
3146:       }
3147:       // Copy offsets
3148:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3149:         index_type idx;

3151:         idx.prefix = field->getFiberDimension(*c_iter);
3152:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
3153:         for(int s = 0; s < space; ++s) {
3154:           idx.index += this->getFiberDimension(*c_iter, s);
3155:         }
3156:         newAtlas->addPoint(*c_iter);
3157:         newAtlas->updatePoint(*c_iter, &idx);
3158:       }
3159:       field->replaceStorage(this->_array, true, this->getStorageSize());
3160:       field->setAtlas(newAtlas);
3161:       return field;
3162:     };
3163:   public:
3164:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3165:       ostringstream txt;
3166:       int rank;

3168:       if (comm == MPI_COMM_NULL) {
3169:         comm = this->comm();
3170:         rank = this->commRank();
3171:       } else {
3172:         MPI_Comm_rank(comm, &rank);
3173:       }
3174:       if (name == "") {
3175:         if(rank == 0) {
3176:           txt << "viewing a FEMSection" << std::endl;
3177:         }
3178:       } else {
3179:         if (rank == 0) {
3180:           txt << "viewing FEMSection '" << name << "'" << std::endl;
3181:         }
3182:       }
3183:       if (rank == 0) {
3184:         txt << "  Fields: " << this->getNumSpaces() << std::endl;
3185:       }
3186:       const chart_type& chart = this->getChart();

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

3192:         if (dim != 0) {
3193:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << "  ";
3194:           for(int i = 0; i < dim; i++) {
3195:             txt << " " << array[i];
3196:           }
3197:           const int& dim = this->getConstraintDimension(*p_iter);

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

3202:             txt << " constrained";
3203:             for(int i = 0; i < dim; ++i) {
3204:               txt << " " << bcArray[i];
3205:             }
3206:           }
3207:           txt << std::endl;
3208:         }
3209:       }
3210:       if (chart.size() == 0) {
3211:         txt << "[" << this->commRank() << "]: empty" << std::endl;
3212:       }
3213:       PetscSynchronizedPrintf(comm, txt.str().c_str());
3214:       PetscSynchronizedFlush(comm);
3215:     };
3216:   };
3217:   // A Field combines several sections
3218:   template<typename Overlap_, typename Patch_, typename Section_>
3219:   class Field : public ALE::ParallelObject {
3220:   public:
3221:     typedef Overlap_                                 overlap_type;
3222:     typedef Patch_                                   patch_type;
3223:     typedef Section_                                 section_type;
3224:     typedef typename section_type::point_type        point_type;
3225:     typedef typename section_type::chart_type        chart_type;
3226:     typedef typename section_type::value_type        value_type;
3227:     typedef std::map<patch_type, Obj<section_type> > sheaf_type;
3228:     typedef enum {SEND, RECEIVE}                     request_type;
3229:     typedef std::map<patch_type, MPI_Request>        requests_type;
3230:   protected:
3231:     sheaf_type    _sheaf;
3232:     int           _tag;
3233:     MPI_Datatype  _datatype;
3234:     requests_type _requests;
3235:   public:
3236:     Field(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
3237:       this->_tag      = this->getNewTag();
3238:       this->_datatype = this->getMPIDatatype();
3239:     };
3240:     Field(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _tag(tag) {
3241:       this->_datatype = this->getMPIDatatype();
3242:     };
3243:     virtual ~Field() {};
3244:   protected:
3245:     MPI_Datatype getMPIDatatype() {
3246:       if (sizeof(value_type) == 4) {
3247:         return MPI_INT;
3248:       } else if (sizeof(value_type) == 8) {
3249:         return MPI_DOUBLE;
3250:       } else if (sizeof(value_type) == 28) {
3251:         int          blen[2];
3252:         MPI_Aint     indices[2];
3253:         MPI_Datatype oldtypes[2], newtype;
3254:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
3255:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3256:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3257:         MPI_Type_commit(&newtype);
3258:         return newtype;
3259:       } else if (sizeof(value_type) == 32) {
3260:         int          blen[2];
3261:         MPI_Aint     indices[2];
3262:         MPI_Datatype oldtypes[2], newtype;
3263:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_DOUBLE;
3264:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3265:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3266:         MPI_Type_commit(&newtype);
3267:         return newtype;
3268:       }
3269:       throw ALE::Exception("Cannot determine MPI type for value type");
3270:     };
3271:     int getNewTag() {
3272:       static int tagKeyval = MPI_KEYVAL_INVALID;
3273:       int *tagvalp = NULL, *maxval, flg;

3275:       if (tagKeyval == MPI_KEYVAL_INVALID) {
3276:         tagvalp = (int *) malloc(sizeof(int));
3277:         MPI_Keyval_create(MPI_NULL_COPY_FN, Mesh_DelTag, &tagKeyval, (void *) NULL);
3278:         MPI_Attr_put(this->_comm, tagKeyval, tagvalp);
3279:         tagvalp[0] = 0;
3280:       }
3281:       MPI_Attr_get(this->_comm, tagKeyval, (void **) &tagvalp, &flg);
3282:       if (tagvalp[0] < 1) {
3283:         MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
3284:         tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
3285:       }
3286:       if (this->debug()) {
3287:         std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
3288:       }
3289:       return tagvalp[0]--;
3290:     };
3291:   public: // Verifiers
3292:     void checkPatch(const patch_type& patch) const {
3293:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3294:         ostringstream msg;
3295:         msg << "Invalid field patch " << patch << std::endl;
3296:         throw ALE::Exception(msg.str().c_str());
3297:       }
3298:     };
3299:     bool hasPatch(const patch_type& patch) {
3300:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3301:         return false;
3302:       }
3303:       return true;
3304:     };
3305:   public: // Accessors
3306:     int getTag() const {return this->_tag;};
3307:     void setTag(const int tag) {this->_tag = tag;};
3308:     Obj<section_type>& getSection(const patch_type& patch) {
3309:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3310:         this->_sheaf[patch] = new section_type(this->comm(), this->debug());
3311:       }
3312:       return this->_sheaf[patch];
3313:     };
3314:     void setSection(const patch_type& patch, const Obj<section_type>& section) {this->_sheaf[patch] = section;};
3315:     const sheaf_type& getPatches() {
3316:       return this->_sheaf;
3317:     };
3318:     void clear() {
3319:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3320:         p_iter->second->clear();
3321:       }
3322:     };
3323:   public: //  Adapter
3324:     template<typename Topology_>
3325:     void setTopology(const Obj<Topology_>& topology) {
3326:       const typename Topology_::sheaf_type& patches = topology->getPatches();

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

3333:         for(typename Topology_::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
3334:           section->setFiberDimension(*b_iter, 1);
3335:         }
3336:       }
3337:     }
3338:     void allocate() {
3339:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3340:         p_iter->second->allocatePoint();
3341:       }
3342:     }
3343:   public: // Communication
3344:     void construct(const int size) {
3345:       const sheaf_type& patches = this->getPatches();

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

3352:         for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3353:           section->setFiberDimension(*c_iter, size);
3354:         }
3355:       }
3356:     };
3357:     template<typename Sizer>
3358:     void construct(const Obj<Sizer>& sizer) {
3359:       const sheaf_type& patches = this->getPatches();

3361:       for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3362:         const patch_type         rank    = p_iter->first;
3363:         const Obj<section_type>& section = this->getSection(rank);
3364:         const chart_type&        chart   = section->getChart();
3365: 
3366:         for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3367:           section->setFiberDimension(*c_iter, *(sizer->getSection(rank)->restrictPoint(*c_iter)));
3368:         }
3369:       }
3370:     }
3371:     void constructCommunication(const request_type& requestType) {
3372:       const sheaf_type& patches = this->getPatches();

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

3379:         if (requestType == RECEIVE) {
3380:           if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data(" << section->size() << ") from " << patch << " tag " << this->_tag << std::endl;}
3381:           MPI_Recv_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3382:         } else {
3383:           if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << section->size() << ") to " << patch << " tag " << this->_tag << std::endl;}
3384:           MPI_Send_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3385:         }
3386:         this->_requests[patch] = request;
3387:       }
3388:     };
3389:     void startCommunication() {
3390:       const sheaf_type& patches = this->getPatches();

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

3395:         MPI_Start(&request);
3396:       }
3397:     };
3398:     void endCommunication() {
3399:       const sheaf_type& patches = this->getPatches();
3400:       MPI_Status status;

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

3405:         MPI_Wait(&request, &status);
3406:       }
3407:     };
3408:   public:
3409:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3410:       ostringstream txt;
3411:       int rank;

3413:       if (comm == MPI_COMM_NULL) {
3414:         comm = this->comm();
3415:         rank = this->commRank();
3416:       } else {
3417:         MPI_Comm_rank(comm, &rank);
3418:       }
3419:       if (name == "") {
3420:         if(rank == 0) {
3421:           txt << "viewing a Field" << std::endl;
3422:         }
3423:       } else {
3424:         if(rank == 0) {
3425:           txt << "viewing Field '" << name << "'" << std::endl;
3426:         }
3427:       }
3428:       PetscSynchronizedPrintf(comm, txt.str().c_str());
3429:       PetscSynchronizedFlush(comm);
3430:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3431:         ostringstream txt1;

3433:         txt1 << "[" << this->commRank() << "]: Patch " << p_iter->first << std::endl;
3434:         PetscSynchronizedPrintf(comm, txt1.str().c_str());
3435:         PetscSynchronizedFlush(comm);
3436:         p_iter->second->view("field section", comm);
3437:       }
3438:     };
3439:   };
3440: }
3441: #endif