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->_atlas->clear();
393:       this->_array.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;};
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->_values.find(point) != this->_values.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->_atlas->clear();
781:       this->_array.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() const {
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() const {
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:     const 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:       this->_atlas->clear();
1038:       for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1039:       this->_allocator.deallocate(this->_array, totalSize);
1040:       ///delete [] this->_array;
1041:       this->_array = NULL;
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:     Obj<atlas_type> _atlas;
1372:     Obj<atlas_type> _atlasNew;
1373:     values_type     _array;
1374:     bool            _sharedStorage;
1375:     int             _sharedStorageSize;
1376:     Obj<bc_type>    _bc;
1377:     alloc_type      _allocator;
1378:     std::vector<Obj<atlas_type> > _spaces;
1379:     std::vector<Obj<bc_type> >    _bcs;
1380:     // Optimization
1381:     std::vector<customAtlas_type> _customRestrictAtlas;
1382:     std::vector<customAtlas_type> _customUpdateAtlas;
1383:   public:
1384:     GeneralSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
1385:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
1386:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
1387:       this->_atlas         = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
1388:       bc_ptr pBC           = bc_alloc_type(this->_allocator).allocate(1);
1389:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
1390:       this->_bc            = Obj<bc_type>(pBC, sizeof(bc_type));
1391:       this->_atlasNew      = NULL;
1392:       this->_array         = NULL;
1393:       this->_sharedStorage = false;
1394:     };
1395:     GeneralSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
1396:       bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
1397:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
1398:       this->_bc  = Obj<bc_type>(pBC, sizeof(bc_type));
1399:     };
1400:     ~GeneralSection() {
1401:       if (this->_array && !this->_sharedStorage) {
1402:         const int totalSize = this->sizeWithBC();

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

1428:       if (size > maxSize) {
1429:         const value_type dummy(0);

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

1489:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1490:         if (this->getConstraintDimension(*p_iter)) {
1491:           this->updatePointBCFull(*p_iter, section->restrictPoint(*p_iter));
1492:         }
1493:       }
1494:       this->copyFibration(section);
1495:     };
1496:     void defaultConstraintDof() {
1497:       const chart_type& chart = this->getChart();
1498:       int size = 0;

1500:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1501:         size = std::max(size, this->getConstraintDimension(*p_iter));
1502:       }
1503:       int *dofs = new int[size];
1504:       for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1505:         const int cDim = this->getConstraintDimension(*p_iter);

1507:         if (cDim) {
1508:           for(int d = 0; d < cDim; ++d) {
1509:             dofs[d] = d;
1510:           }
1511:           this->_bc->updatePoint(*p_iter, dofs);
1512:         }
1513:       }
1514:       delete [] dofs;
1515:     };
1516:   public: // Sizes
1517:     void clear() {
1518:       this->_atlas->clear();
1519:       if (!this->_sharedStorage) {
1520:         const int totalSize = this->sizeWithBC();

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

1565:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1566:         size += this->getConstrainedFiberDimension(*p_iter);
1567:       }
1568:       return size;
1569:     };
1570:     // Return the total number of dofs (free and constrained)
1571:     int sizeWithBC() const {
1572:       const chart_type& points = this->getChart();
1573:       int               size   = 0;

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

1600:         for(int i = start; i < end; ++i) {
1601:           indices[(*indx)++] = i;
1602:         }
1603:       } else {
1604:         const int numSpaces = this->getNumSpaces();
1605:         int offset = start;

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

1610:           for(int i = offset+dim-1; i >= offset; --i) {
1611:             indices[(*indx)++] = i;
1612:           }
1613:           offset += dim;
1614:         }
1615:         if (!numSpaces) {
1616:           const int& dim = this->getFiberDimension(p);

1618:           for(int i = offset+dim-1; i >= offset; --i) {
1619:             indices[(*indx)++] = i;
1620:           }
1621:           offset += dim;
1622:         }
1623:       }
1624:     }
1625:     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) {
1626:       const int& cDim = this->getConstraintDimension(p);

1628:       if (!cDim) {
1629:         this->getIndicesRaw(p, start, indices, indx, orientation);
1630:       } else {
1631:         if (orientation >= 0) {
1632:           const int&                          dim  = this->getFiberDimension(p);
1633:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1634:           int                                 cInd = 0;

1636:           for(int i = start, k = 0; k < dim; ++k) {
1637:             if ((cInd < cDim) && (k == cDof[cInd])) {
1638:               if (!freeOnly) indices[(*indx)++] = -(k+1);
1639:               if (skipConstraints) ++i;
1640:               ++cInd;
1641:             } else {
1642:               indices[(*indx)++] = i++;
1643:             }
1644:           }
1645:         } else {
1646:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1647:           int                                 offset  = 0;
1648:           int                                 cOffset = 0;
1649:           int                                 j       = -1;

1651:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1652:             const int  dim = this->getFiberDimension(p, space);
1653:             const int tDim = this->getConstrainedFiberDimension(p, space);
1654:             int       cInd = (dim - tDim)-1;

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

1678:       this->_array             = this->_allocator.allocate(totalSize);
1679:       ///this->_array             = new value_type[totalSize];
1680:       this->_sharedStorage     = false;
1681:       this->_sharedStorageSize = 0;
1682:       for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1683:       ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1684:       this->_bc->allocatePoint();
1685:     };
1686:     void replaceStorage(value_type *newArray, const bool sharedStorage = false, const int sharedStorageSize = 0) {
1687:       if (this->_array && !this->_sharedStorage) {
1688:         const int totalSize = this->sizeWithBC();

1690:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1691:         this->_allocator.deallocate(this->_array, totalSize);
1692:         ///delete [] this->_array;
1693:       }
1694:       this->_array             = newArray;
1695:       this->_sharedStorage     = sharedStorage;
1696:       this->_sharedStorageSize = sharedStorageSize;
1697:       this->_atlas             = this->_atlasNew;
1698:       this->_atlasNew          = NULL;
1699:     };
1700:     void addPoint(const point_type& point, const int dim) {
1701:       if (dim == 0) return;
1702:       if (this->_atlasNew.isNull()) {
1703:         this->_atlasNew = new atlas_type(this->comm(), this->debug());
1704:         this->_atlasNew->copy(this->_atlas);
1705:       }
1706:       const index_type idx(dim, -1);
1707:       this->_atlasNew->addPoint(point);
1708:       this->_atlasNew->updatePoint(point, &idx);
1709:     };
1710:     void orderPoints(const Obj<atlas_type>& atlas){
1711:       const chart_type& chart  = this->getChart();
1712:       int               offset = 0;

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

1718:         idx.index = offset;
1719:         atlas->updatePoint(*c_iter, &idx);
1720:         offset += dim;
1721:       }
1722:     };
1723:     void allocatePoint() {
1724:       this->orderPoints(this->_atlas);
1725:       this->allocateStorage();
1726:     };
1727:   public: // Restriction and Update
1728:     // Zero entries
1729:     void zero() {
1730:       this->set(0.0);
1731:     };
1732:     void set(const value_type value) {
1733:       //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1734:       const chart_type& chart = this->getChart();

1736:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1737:         value_type *array = (value_type *) this->restrictPoint(*c_iter);
1738:         const int&  dim   = this->getFiberDimension(*c_iter);
1739:         const int&  cDim  = this->getConstraintDimension(*c_iter);

1741:         if (!cDim) {
1742:           for(int i = 0; i < dim; ++i) {
1743:             array[i] = value;
1744:           }
1745:         } else {
1746:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1747:           int                                 cInd = 0;

1749:           for(int i = 0; i < dim; ++i) {
1750:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1751:             array[i] = value;
1752:           }
1753:         }
1754:       }
1755:     };
1756:     // Add two sections and put the result in a third
1757:     void add(const Obj<GeneralSection>& x, const Obj<GeneralSection>& y) {
1758:       // Check atlases
1759:       const chart_type& chart = this->getChart();

1761:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1762:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
1763:         const value_type *xArray = x->restrictPoint(*c_iter);
1764:         const value_type *yArray = y->restrictPoint(*c_iter);
1765:         const int&        dim    = this->getFiberDimension(*c_iter);
1766:         const int&        cDim   = this->getConstraintDimension(*c_iter);

1768:         if (!cDim) {
1769:           for(int i = 0; i < dim; ++i) {
1770:             array[i] = xArray[i] + yArray[i];
1771:           }
1772:         } else {
1773:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1774:           int                                 cInd = 0;

1776:           for(int i = 0; i < dim; ++i) {
1777:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1778:             array[i] = xArray[i] + yArray[i];
1779:           }
1780:         }
1781:       }
1782:     };
1783:     // this = this + alpha * x
1784:     template<typename OtherSection>
1785:     void axpy(const value_type alpha, const Obj<OtherSection>& x) {
1786:       // Check atlases
1787:       const chart_type& chart = this->getChart();

1789:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1790:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
1791:         const value_type *xArray = x->restrictPoint(*c_iter);
1792:         const int&        dim    = this->getFiberDimension(*c_iter);
1793:         const int&        cDim   = this->getConstraintDimension(*c_iter);

1795:         if (!cDim) {
1796:           for(int i = 0; i < dim; ++i) {
1797:             array[i] += alpha*xArray[i];
1798:           }
1799:         } else {
1800:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1801:           int                                 cInd = 0;

1803:           for(int i = 0; i < dim; ++i) {
1804:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1805:             array[i] += alpha*xArray[i];
1806:           }
1807:         }
1808:       }
1809:     };
1810:     // Return the free values on a point
1811:     const value_type *restrictSpace() const {
1812:       return this->_array;
1813:     };
1814:     // Return the free values on a point
1815:     const value_type *restrictPoint(const point_type& p) const {
1816:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1817:     };
1818:     void restrictPoint(const point_type& p, value_type values[], const int size) const {
1819:       assert(this->_atlas->restrictPoint(p)[0].prefix == size);
1820:       const value_type *v = &(this->_array[this->_atlas->restrictPoint(p)[0].index]);

1822:       for(int i = 0; i < size; ++i) {
1823:         values[i] = v[i];
1824:       }
1825:     };
1826:     // Update the free values on a point
1827:     //   Takes a full array and ignores constrained values
1828:     void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1829:       value_type *array = (value_type *) this->restrictPoint(p);
1830:       const int&  cDim  = this->getConstraintDimension(p);

1832:       if (!cDim) {
1833:         if (orientation >= 0) {
1834:           const int& dim = this->getFiberDimension(p);

1836:           for(int i = 0; i < dim; ++i) {
1837:             array[i] = v[i];
1838:           }
1839:         } else {
1840:           int offset = 0;
1841:           int j      = -1;

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

1846:             for(int i = dim-1; i >= 0; --i) {
1847:               array[++j] = v[i+offset];
1848:             }
1849:             offset += dim;
1850:           }
1851:         }
1852:       } else {
1853:         if (orientation >= 0) {
1854:           const int&                          dim  = this->getFiberDimension(p);
1855:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1856:           int                                 cInd = 0;

1858:           for(int i = 0; i < dim; ++i) {
1859:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1860:             array[i] = v[i];
1861:           }
1862:         } else {
1863:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1864:           int                                 offset  = 0;
1865:           int                                 cOffset = 0;
1866:           int                                 j       = 0;

1868:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1869:             const int  dim = this->getFiberDimension(p, space);
1870:             const int tDim = this->getConstrainedFiberDimension(p, space);
1871:             const int sDim = dim - tDim;
1872:             int       cInd = 0;

1874:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1875:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1876:               array[j] = v[k];
1877:             }
1878:             offset  += dim;
1879:             cOffset += dim - tDim;
1880:           }
1881:         }
1882:       }
1883:     };
1884:     // Update the free values on a point
1885:     //   Takes a full array and ignores constrained values
1886:     void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1887:       value_type *array = (value_type *) this->restrictPoint(p);
1888:       const int&  cDim  = this->getConstraintDimension(p);

1890:       if (!cDim) {
1891:         if (orientation >= 0) {
1892:           const int& dim = this->getFiberDimension(p);

1894:           for(int i = 0; i < dim; ++i) {
1895:             array[i] += v[i];
1896:           }
1897:         } else {
1898:           int offset = 0;
1899:           int j      = -1;

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

1904:             for(int i = dim-1; i >= 0; --i) {
1905:               array[++j] += v[i+offset];
1906:             }
1907:             offset += dim;
1908:           }
1909:         }
1910:       } else {
1911:         if (orientation >= 0) {
1912:           const int&                          dim  = this->getFiberDimension(p);
1913:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1914:           int                                 cInd = 0;

1916:           for(int i = 0; i < dim; ++i) {
1917:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1918:             array[i] += v[i];
1919:           }
1920:         } else {
1921:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1922:           int                                 offset  = 0;
1923:           int                                 cOffset = 0;
1924:           int                                 j       = 0;

1926:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1927:             const int  dim = this->getFiberDimension(p, space);
1928:             const int tDim = this->getConstrainedFiberDimension(p, space);
1929:             const int sDim = dim - tDim;
1930:             int       cInd = 0;

1932:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1933:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1934:               array[j] += v[k];
1935:             }
1936:             offset  += dim;
1937:             cOffset += dim - tDim;
1938:           }
1939:         }
1940:       }
1941:     };
1942:     // Update the free values on a point
1943:     //   Takes ONLY unconstrained values
1944:     void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1945:       value_type *array = (value_type *) this->restrictPoint(p);
1946:       const int&  cDim  = this->getConstraintDimension(p);

1948:       if (!cDim) {
1949:         if (orientation >= 0) {
1950:           const int& dim = this->getFiberDimension(p);

1952:           for(int i = 0; i < dim; ++i) {
1953:             array[i] = v[i];
1954:           }
1955:         } else {
1956:           int offset = 0;
1957:           int j      = -1;

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

1962:             for(int i = dim-1; i >= 0; --i) {
1963:               array[++j] = v[i+offset];
1964:             }
1965:             offset += dim;
1966:           }
1967:         }
1968:       } else {
1969:         if (orientation >= 0) {
1970:           const int&                          dim  = this->getFiberDimension(p);
1971:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1972:           int                                 cInd = 0;

1974:           for(int i = 0, k = -1; i < dim; ++i) {
1975:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1976:             array[i] = v[++k];
1977:           }
1978:         } else {
1979:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
1980:           int                                 offset  = 0;
1981:           int                                 cOffset = 0;
1982:           int                                 j       = 0;

1984:           for(int space = 0; space < this->getNumSpaces(); ++space) {
1985:             const int  dim = this->getFiberDimension(p, space);
1986:             const int tDim = this->getConstrainedFiberDimension(p, space);
1987:             const int sDim = dim - tDim;
1988:             int       cInd = 0;

1990:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
1991:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1992:               array[j] = v[--k];
1993:             }
1994:             offset  += dim;
1995:             cOffset += dim - tDim;
1996:           }
1997:         }
1998:       }
1999:     };
2000:     // Update the free values on a point
2001:     //   Takes ONLY unconstrained values
2002:     void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2003:       value_type *array = (value_type *) this->restrictPoint(p);
2004:       const int&  cDim  = this->getConstraintDimension(p);

2006:       if (!cDim) {
2007:         if (orientation >= 0) {
2008:           const int& dim = this->getFiberDimension(p);

2010:           for(int i = 0; i < dim; ++i) {
2011:             array[i] += v[i];
2012:           }
2013:         } else {
2014:           int offset = 0;
2015:           int j      = -1;

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

2020:             for(int i = dim-1; i >= 0; --i) {
2021:               array[++j] += v[i+offset];
2022:             }
2023:             offset += dim;
2024:           }
2025:         }
2026:       } else {
2027:         if (orientation >= 0) {
2028:           const int&                          dim  = this->getFiberDimension(p);
2029:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2030:           int                                 cInd = 0;

2032:           for(int i = 0, k = -1; i < dim; ++i) {
2033:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2034:             array[i] += v[++k];
2035:           }
2036:         } else {
2037:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2038:           int                                 offset  = 0;
2039:           int                                 cOffset = 0;
2040:           int                                 j       = 0;

2042:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2043:             const int  dim = this->getFiberDimension(p, space);
2044:             const int tDim = this->getConstrainedFiberDimension(p, space);
2045:             const int sDim = dim - tDim;
2046:             int       cInd = 0;

2048:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2049:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2050:               array[j] += v[--k];
2051:             }
2052:             offset  += dim;
2053:             cOffset += dim - tDim;
2054:           }
2055:         }
2056:       }
2057:     };
2058:     // Update only the constrained dofs on a point
2059:     //   This takes an array with ONLY bc values
2060:     void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2061:       value_type *array = (value_type *) this->restrictPoint(p);
2062:       const int&  cDim  = this->getConstraintDimension(p);

2064:       if (cDim) {
2065:         if (orientation >= 0) {
2066:           const int&                          dim  = this->getFiberDimension(p);
2067:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2068:           int                                 cInd = 0;

2070:           for(int i = 0; i < dim; ++i) {
2071:             if (cInd == cDim) break;
2072:             if (i == cDof[cInd]) {
2073:               array[i] = v[cInd];
2074:               ++cInd;
2075:             }
2076:           }
2077:         } else {
2078:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2079:           int                                 cOffset = 0;
2080:           int                                 j       = 0;

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

2087:             for(int i = 0; i < dim; ++i, ++j) {
2088:               if (cInd < 0) break;
2089:               if (j == cDof[cInd+cOffset]) {
2090:                 array[j] = v[cInd+cOffset];
2091:                 ++cInd;
2092:               }
2093:             }
2094:             cOffset += dim - tDim;
2095:           }
2096:         }
2097:       }
2098:     };
2099:     // Update only the constrained dofs on a point
2100:     //   This takes an array with ALL values, not just BC
2101:     void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2102:       value_type *array = (value_type *) this->restrictPoint(p);
2103:       const int&  cDim  = this->getConstraintDimension(p);

2105:       if (cDim) {
2106:         if (orientation >= 0) {
2107:           const int&                          dim  = this->getFiberDimension(p);
2108:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2109:           int                                 cInd = 0;

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

2124:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2125:             const int  dim = this->getFiberDimension(p, space);
2126:             const int tDim = this->getConstrainedFiberDimension(p, space);
2127:             int       cInd = 0;

2129:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2130:               if (cInd < 0) break;
2131:               if (j == cDof[cInd+cOffset]) {
2132:                 array[j] = v[k];
2133:                 ++cInd;
2134:               }
2135:             }
2136:             offset  += dim;
2137:             cOffset += dim - tDim;
2138:           }
2139:         }
2140:       }
2141:     };
2142:     // Update all dofs on a point (free and constrained)
2143:     void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
2144:       value_type *array = (value_type *) this->restrictPoint(p);

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

2149:         for(int i = 0; i < dim; ++i) {
2150:           array[i] = v[i];
2151:         }
2152:       } else {
2153:         int offset = 0;
2154:         int j      = -1;

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

2159:           for(int i = dim-1; i >= 0; --i) {
2160:             array[++j] = v[i+offset];
2161:           }
2162:           offset += dim;
2163:         }
2164:       }
2165:     };
2166:     // Update all dofs on a point (free and constrained)
2167:     void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
2168:       value_type *array = (value_type *) this->restrictPoint(p);

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

2173:         for(int i = 0; i < dim; ++i) {
2174:           array[i] += v[i];
2175:         }
2176:       } else {
2177:         int offset = 0;
2178:         int j      = -1;

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

2183:           for(int i = dim-1; i >= 0; --i) {
2184:             array[++j] += v[i+offset];
2185:           }
2186:           offset += dim;
2187:         }
2188:       }
2189:     };
2190:   public: // Fibrations
2191:     int getNumSpaces() const {return this->_spaces.size();};
2192:     const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
2193:     const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
2194:     void addSpace() {
2195:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
2196:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
2197:       this->_spaces.push_back(space);
2198:       this->_bcs.push_back(bc);
2199:     };
2200:     int getFiberDimension(const point_type& p, const int space) const {
2201:       return this->_spaces[space]->restrictPoint(p)->prefix;
2202:     };
2203:     void setFiberDimension(const point_type& p, int dim, const int space) {
2204:       const index_type idx(dim, -1);
2205:       this->_spaces[space]->addPoint(p);
2206:       this->_spaces[space]->updatePoint(p, &idx);
2207:     };
2208:     template<typename Sequence>
2209:     void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
2210:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2211:         this->setFiberDimension(*p_iter, dim, space);
2212:       }
2213:     };
2214:     const int getConstraintDimension(const point_type& p, const int space) const {
2215:       if (!this->_bcs[space]->hasPoint(p)) return 0;
2216:       return this->_bcs[space]->getFiberDimension(p);
2217:     };
2218:     void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2219:       this->_bcs[space]->setFiberDimension(p, numConstraints);
2220:     };
2221:     int getConstrainedFiberDimension(const point_type& p, const int space) const {
2222:       return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
2223:     };
2224:     template<typename OtherSection>
2225:     void copyFibration(const Obj<OtherSection>& section) {
2226:       const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
2227:       const std::vector<Obj<bc_type> >&    bcs    = section->getBCs();

2229:       this->_spaces.clear();
2230:       for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
2231:         this->_spaces.push_back(*s_iter);
2232:       }
2233:       this->_bcs.clear();
2234:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
2235:         this->_bcs.push_back(*b_iter);
2236:       }
2237:     };
2238:     Obj<GeneralSection> getFibration(const int space) const {
2239:       Obj<GeneralSection> field = new GeneralSection(this->comm(), this->debug());
2240: //     Obj<atlas_type> _atlas;
2241: //     std::vector<Obj<atlas_type> > _spaces;
2242: //     Obj<bc_type>    _bc;
2243: //     std::vector<Obj<bc_type> >    _bcs;
2244:       field->addSpace();
2245:       const chart_type& chart = this->getChart();

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

2252:         if (fDim) {
2253:           field->setFiberDimension(*c_iter, fDim);
2254:           field->setFiberDimension(*c_iter, fDim, 0);
2255:         }
2256:         if (cDim) {
2257:           field->setConstraintDimension(*c_iter, cDim);
2258:           field->setConstraintDimension(*c_iter, cDim, 0);
2259:         }
2260:       }
2261:       field->allocateStorage();
2262:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
2263:       const chart_type& newChart = field->getChart();

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

2269:         if (cDim) {
2270:           field->setConstraintDof(*c_iter, dof);
2271:         }
2272:       }
2273:       // Copy offsets
2274:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2275:         index_type idx;

2277:         idx.prefix = field->getFiberDimension(*c_iter);
2278:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
2279:         for(int s = 0; s < space; ++s) {
2280:           idx.index += this->getFiberDimension(*c_iter, s);
2281:         }
2282:         newAtlas->addPoint(*c_iter);
2283:         newAtlas->updatePoint(*c_iter, &idx);
2284:       }
2285:       field->replaceStorage(this->_array, true, this->getStorageSize());
2286:       field->setAtlas(newAtlas);
2287:       return field;
2288:     };
2289:   public: // Optimization
2290:     void getCustomRestrictAtlas(const int tag, const int *offsets[], const int *indices[]) {
2291:       *offsets = this->_customRestrictAtlas[tag].first.first;
2292:       *indices = this->_customRestrictAtlas[tag].first.second;
2293:     };
2294:     void getCustomUpdateAtlas(const int tag, const int *offsets[], const int *indices[]) {
2295:       *offsets = this->_customUpdateAtlas[tag].first.first;
2296:       *indices = this->_customUpdateAtlas[tag].first.second;
2297:     };
2298:     // This returns the tag assigned to the atlas
2299:     int setCustomAtlas(const int restrictOffsets[], const int restrictIndices[], const int updateOffsets[], const int updateIndices[], bool autoFree = true) {
2300:       this->_customRestrictAtlas.push_back(customAtlas_type(customAtlasInd_type(restrictOffsets, restrictIndices), autoFree));
2301:       this->_customUpdateAtlas.push_back(customAtlas_type(customAtlasInd_type(updateOffsets, updateIndices), autoFree));
2302:       return this->_customUpdateAtlas.size()-1;
2303:     };
2304:     int copyCustomAtlas(const Obj<GeneralSection>& section, const int tag) {
2305:       const int *rOffsets, *rIndices, *uOffsets, *uIndices;

2307:       section->getCustomRestrictAtlas(tag, &rOffsets, &rIndices);
2308:       section->getCustomUpdateAtlas(tag, &uOffsets, &uIndices);
2309:       return this->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices, false);
2310:     };
2311:   public:
2312:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
2313:       ostringstream txt;
2314:       int rank;

2316:       if (comm == MPI_COMM_NULL) {
2317:         comm = this->comm();
2318:         rank = this->commRank();
2319:       } else {
2320:         MPI_Comm_rank(comm, &rank);
2321:       }
2322:       if (name == "") {
2323:         if(rank == 0) {
2324:           txt << "viewing a GeneralSection" << std::endl;
2325:         }
2326:       } else {
2327:         if (rank == 0) {
2328:           txt << "viewing GeneralSection '" << name << "'" << std::endl;
2329:         }
2330:       }
2331:       if (rank == 0) {
2332:         txt << "  Fields: " << this->getNumSpaces() << std::endl;
2333:       }
2334:       const chart_type& chart = this->getChart();

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

2340:         if (dim != 0) {
2341:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << "  ";
2342:           for(int i = 0; i < dim; i++) {
2343:             txt << " " << array[i];
2344:           }
2345:           const int& dim = this->getConstraintDimension(*p_iter);

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

2350:             txt << " constrained";
2351:             for(int i = 0; i < dim; ++i) {
2352:               txt << " " << bcArray[i];
2353:             }
2354:           }
2355:           txt << std::endl;
2356:         }
2357:       }
2358:       if (chart.size() == 0) {
2359:         txt << "[" << this->commRank() << "]: empty" << std::endl;
2360:       }
2361:       PetscSynchronizedPrintf(comm, txt.str().c_str());
2362:       PetscSynchronizedFlush(comm);
2363:     };
2364:   };
2365:   // FEMSection will support vector BC on a subset of unknowns on a point
2366:   //   We make a separate constraint Section to hold the transformation and projection operators
2367:   //   Storage will be contiguous by node, just as in Section
2368:   //     This allows fast restrict(p)
2369:   //     Then update() is accomplished by projecting constrained unknowns
2370:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
2371:            typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
2372:            typename BCAtlas_ = UniformSection<Point_, int, 1, typename Alloc_::template rebind<int>::other>,
2373:            typename BC_ = Section<Point_, double, typename Alloc_::template rebind<double>::other> >
2374:   class FEMSection : public ALE::ParallelObject {
2375:   public:
2376:     typedef Point_                                                  point_type;
2377:     typedef Value_                                                  value_type;
2378:     typedef Alloc_                                                  alloc_type;
2379:     typedef Atlas_                                                  atlas_type;
2380:     typedef BCAtlas_                                                bc_atlas_type;
2381:     typedef BC_                                                     bc_type;
2382:     typedef Point                                                   index_type;
2383:     typedef typename atlas_type::chart_type                         chart_type;
2384:     typedef value_type *                                            values_type;
2385:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
2386:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
2387:     typedef typename alloc_type::template rebind<bc_type>::other    bc_atlas_alloc_type;
2388:     typedef typename bc_atlas_alloc_type::pointer                   bc_atlas_ptr;
2389:     typedef typename alloc_type::template rebind<bc_type>::other    bc_alloc_type;
2390:     typedef typename bc_alloc_type::pointer                         bc_ptr;
2391:   protected:
2392:     Obj<atlas_type>    _atlas;
2393:     Obj<bc_atlas_type> _bc_atlas;
2394:     Obj<bc_type>       _bc;
2395:     values_type        _array;
2396:     bool               _sharedStorage;
2397:     int                _sharedStorageSize;
2398:     alloc_type         _allocator;
2399:     std::vector<Obj<atlas_type> > _spaces;
2400:     std::vector<Obj<bc_type> >    _bcs;
2401:   public:
2402:     FEMSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
2403:       atlas_ptr pAtlas      = atlas_alloc_type(this->_allocator).allocate(1);
2404:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
2405:       this->_atlas          = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
2406:       bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2407:       bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(comm, debug));
2408:       this->_bc_atlas       = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2409:       bc_ptr pBC            = bc_alloc_type(this->_allocator).allocate(1);
2410:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
2411:       this->_bc             = Obj<bc_type>(pBC, sizeof(bc_type));
2412:       this->_array          = NULL;
2413:       this->_sharedStorage  = false;
2414:     };
2415:     FEMSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
2416:       bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2417:       bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(comm, debug));
2418:       this->_bc_atlas       = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2419:       bc_ptr pBC            = bc_alloc_type(this->_allocator).allocate(1);
2420:       bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
2421:       this->_bc             = Obj<bc_type>(pBC, sizeof(bc_type));
2422:     };
2423:     ~FEMSection() {
2424:       if (this->_array && !this->_sharedStorage) {
2425:         const int totalSize = this->sizeWithBC();

2427:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2428:         this->_allocator.deallocate(this->_array, totalSize);
2429:         this->_array = NULL;
2430:       }
2431:     };
2432:   public:
2433:     value_type *getRawArray(const int size) {
2434:       // Put in a sentinel value that deallocates the array
2435:       static value_type *array   = NULL;
2436:       static int         maxSize = 0;

2438:       if (size > maxSize) {
2439:         const value_type dummy(0);

2441:         if (array) {
2442:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
2443:           this->_allocator.deallocate(array, maxSize);
2444:         }
2445:         maxSize = size;
2446:         array   = this->_allocator.allocate(maxSize);
2447:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
2448:       };
2449:       return array;
2450:     };
2451:     int getStorageSize() const {
2452:       if (this->_sharedStorage) {
2453:         return this->_sharedStorageSize;
2454:       }
2455:       return this->sizeWithBC();
2456:     };
2457:   public: // Verifiers
2458:     bool hasPoint(const point_type& point) {
2459:       return this->_atlas->hasPoint(point);
2460:     };
2461:   public: // Accessors
2462:     const chart_type& getChart() const {return this->_atlas->getChart();};
2463:     const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
2464:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
2465:     const Obj<bc_type>& getBC() const {return this->_bc;};
2466:     void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
2467:   public: // BC
2468:     // Returns the number of constraints on a point
2469:     const int getConstraintDimension(const point_type& p) const {
2470:       if (!this->_bc_atlas->hasPoint(p)) return 0;
2471:       return this->_bc_atlas->restrictPoint(p)[0];
2472:     };
2473:     // Set the number of constraints on a point
2474:     void setConstraintDimension(const point_type& p, const int numConstraints) {
2475:       this->_bc_atlas->updatePoint(p, &numConstraints);
2476:     };
2477:     // Increment the number of constraints on a point
2478:     void addConstraintDimension(const point_type& p, const int numConstraints) {
2479:       this->_bc_atlas->updatePointAdd(p, &numConstraints);
2480:     };
2481:   public: // Sizes
2482:     void clear() {
2483:       this->_atlas->clear();
2484:       if (!this->_sharedStorage) {
2485:         const int totalSize = this->sizeWithBC();

2487:         for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2488:         this->_allocator.deallocate(this->_array, totalSize);
2489:       }
2490:       this->_array = NULL;
2491:       this->_bc_atlas->clear();
2492:       this->_bc->clear();
2493:     };
2494:     // Return the total number of dofs on the point (free and constrained)
2495:     int getFiberDimension(const point_type& p) const {
2496:       return this->_atlas->restrictPoint(p)->prefix;
2497:     };
2498:     int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
2499:       return atlas->restrictPoint(p)->prefix;
2500:     };
2501:     // Set the total number of dofs on the point (free and constrained)
2502:     void setFiberDimension(const point_type& p, int dim) {
2503:       const index_type idx(dim, -1);
2504:       this->_atlas->addPoint(p);
2505:       this->_atlas->updatePoint(p, &idx);
2506:     };
2507:     template<typename Sequence>
2508:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
2509:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2510:         this->setFiberDimension(*p_iter, dim);
2511:       }
2512:     };
2513:     void addFiberDimension(const point_type& p, int dim) {
2514:       if (this->_atlas->hasPoint(p)) {
2515:         const index_type values(dim, 0);
2516:         this->_atlas->updateAddPoint(p, &values);
2517:       } else {
2518:         this->setFiberDimension(p, dim);
2519:       }
2520:     };
2521:     // Return the number of constrained dofs on this point
2522:     int getConstrainedFiberDimension(const point_type& p) const {
2523:       return this->getFiberDimension(p) - this->getConstraintDimension(p);
2524:     };
2525:     // Return the total number of free dofs
2526:     int size() const {
2527:       const chart_type& points = this->getChart();
2528:       int               size   = 0;

2530:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2531:         size += this->getConstrainedFiberDimension(*p_iter);
2532:       }
2533:       return size;
2534:     };
2535:     // Return the total number of dofs (free and constrained)
2536:     int sizeWithBC() const {
2537:       const chart_type& points = this->getChart();
2538:       int               size   = 0;

2540:       for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2541:         size += this->getFiberDimension(*p_iter);
2542:       }
2543:       return size;
2544:     };
2545:   public: // Allocation
2546:     void allocateStorage() {
2547:       const int totalSize = this->sizeWithBC();
2548:       const value_type dummy(0) ;

2550:       this->_array             = this->_allocator.allocate(totalSize);
2551:       this->_sharedStorage     = false;
2552:       this->_sharedStorageSize = 0;
2553:       for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
2554:       this->_bc_atlas->allocatePoint();
2555:       this->_bc->allocatePoint();
2556:     };
2557:     void orderPoints(const Obj<atlas_type>& atlas){
2558:       const chart_type& chart  = this->getChart();
2559:       int               offset = 0;

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

2565:         idx.index = offset;
2566:         atlas->updatePoint(*c_iter, &idx);
2567:         offset += dim;
2568:       }
2569:     };
2570:     void allocatePoint() {
2571:       this->orderPoints(this->_atlas);
2572:       this->allocateStorage();
2573:     };
2574:   public: // Restriction and Update
2575:     // Zero entries
2576:     void zero() {
2577:       this->set(0.0);
2578:     };
2579:     void set(const value_type value) {
2580:       //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
2581:       const chart_type& chart = this->getChart();

2583:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2584:         value_type *array = (value_type *) this->restrictPoint(*c_iter);
2585:         const int&  dim   = this->getFiberDimension(*c_iter);
2586:         const int&  cDim  = this->getConstraintDimension(*c_iter);

2588:         if (!cDim) {
2589:           for(int i = 0; i < dim; ++i) {
2590:             array[i] = value;
2591:           }
2592:         } else {
2593:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2594:           int                                 cInd = 0;

2596:           for(int i = 0; i < dim; ++i) {
2597:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2598:             array[i] = value;
2599:           }
2600:         }
2601:       }
2602:     };
2603:     // Add two sections and put the result in a third
2604:     void add(const Obj<FEMSection>& x, const Obj<FEMSection>& y) {
2605:       // Check atlases
2606:       const chart_type& chart = this->getChart();

2608:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2609:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
2610:         const value_type *xArray = x->restrictPoint(*c_iter);
2611:         const value_type *yArray = y->restrictPoint(*c_iter);
2612:         const int&        dim    = this->getFiberDimension(*c_iter);
2613:         const int&        cDim   = this->getConstraintDimension(*c_iter);

2615:         if (!cDim) {
2616:           for(int i = 0; i < dim; ++i) {
2617:             array[i] = xArray[i] + yArray[i];
2618:           }
2619:         } else {
2620:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2621:           int                                 cInd = 0;

2623:           for(int i = 0; i < dim; ++i) {
2624:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2625:             array[i] = xArray[i] + yArray[i];
2626:           }
2627:         }
2628:       }
2629:     };
2630:     // this = this + alpha * x
2631:     void axpy(const value_type alpha, const Obj<FEMSection>& x) {
2632:       // Check atlases
2633:       const chart_type& chart = this->getChart();

2635:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2636:         value_type       *array  = (value_type *) this->restrictPoint(*c_iter);
2637:         const value_type *xArray = x->restrictPoint(*c_iter);
2638:         const int&        dim    = this->getFiberDimension(*c_iter);
2639:         const int&        cDim   = this->getConstraintDimension(*c_iter);

2641:         if (!cDim) {
2642:           for(int i = 0; i < dim; ++i) {
2643:             array[i] += alpha*xArray[i];
2644:           }
2645:         } else {
2646:           const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2647:           int                                 cInd = 0;

2649:           for(int i = 0; i < dim; ++i) {
2650:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2651:             array[i] += alpha*xArray[i];
2652:           }
2653:         }
2654:       }
2655:     };
2656:     // Return the free values on a point
2657:     const value_type *restrictSpace() const {
2658:       return this->_array;
2659:     };
2660:     // Return the free values on a point
2661:     const value_type *restrictPoint(const point_type& p) const {
2662:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
2663:     };
2664:     // Update the free values on a point
2665:     //   Takes a full array and ignores constrained values
2666:     void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2667:       value_type *array = (value_type *) this->restrictPoint(p);
2668:       const int&  cDim  = this->getConstraintDimension(p);

2670:       if (!cDim) {
2671:         if (orientation >= 0) {
2672:           const int& dim = this->getFiberDimension(p);

2674:           for(int i = 0; i < dim; ++i) {
2675:             array[i] = v[i];
2676:           }
2677:         } else {
2678:           int offset = 0;
2679:           int j      = -1;

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

2684:             for(int i = dim-1; i >= 0; --i) {
2685:               array[++j] = v[i+offset];
2686:             }
2687:             offset += dim;
2688:           }
2689:         }
2690:       } else {
2691:         if (orientation >= 0) {
2692:           const int&                          dim  = this->getFiberDimension(p);
2693:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2694:           int                                 cInd = 0;

2696:           for(int i = 0; i < dim; ++i) {
2697:             if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2698:             array[i] = v[i];
2699:           }
2700:         } else {
2701:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2702:           int                                 offset  = 0;
2703:           int                                 cOffset = 0;
2704:           int                                 j       = 0;

2706:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2707:             const int  dim = this->getFiberDimension(p, space);
2708:             const int tDim = this->getConstrainedFiberDimension(p, space);
2709:             const int sDim = dim - tDim;
2710:             int       cInd = 0;

2712:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2713:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2714:               array[j] = v[k];
2715:             }
2716:             offset  += dim;
2717:             cOffset += dim - tDim;
2718:           }
2719:         }
2720:       }
2721:     };
2722:     // Update the free values on a point
2723:     //   Takes a full array and ignores constrained values
2724:     void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2725:       value_type *array = (value_type *) this->restrictPoint(p);
2726:       const int&  cDim  = this->getConstraintDimension(p);

2728:       if (!cDim) {
2729:         if (orientation >= 0) {
2730:           const int& dim = this->getFiberDimension(p);

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

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

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

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

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

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

2786:       if (!cDim) {
2787:         if (orientation >= 0) {
2788:           const int& dim = this->getFiberDimension(p);

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

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

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

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

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

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

2844:       if (!cDim) {
2845:         if (orientation >= 0) {
2846:           const int& dim = this->getFiberDimension(p);

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

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

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

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

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

2886:             for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2887:               if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2888:               array[j] += v[--k];
2889:             }
2890:             offset  += dim;
2891:             cOffset += dim - tDim;
2892:           }
2893:         }
2894:       }
2895:     };
2896:     // Update only the constrained dofs on a point
2897:     //   This takes an array with ONLY bc values
2898:     void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2899:       value_type *array = (value_type *) this->restrictPoint(p);
2900:       const int&  cDim  = this->getConstraintDimension(p);

2902:       if (cDim) {
2903:         if (orientation >= 0) {
2904:           const int&                          dim  = this->getFiberDimension(p);
2905:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2906:           int                                 cInd = 0;

2908:           for(int i = 0; i < dim; ++i) {
2909:             if (cInd == cDim) break;
2910:             if (i == cDof[cInd]) {
2911:               array[i] = v[cInd];
2912:               ++cInd;
2913:             }
2914:           }
2915:         } else {
2916:           const typename bc_type::value_type *cDof    = this->getConstraintDof(p);
2917:           int                                 cOffset = 0;
2918:           int                                 j       = 0;

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

2925:             for(int i = 0; i < dim; ++i, ++j) {
2926:               if (cInd < 0) break;
2927:               if (j == cDof[cInd+cOffset]) {
2928:                 array[j] = v[cInd+cOffset];
2929:                 ++cInd;
2930:               }
2931:             }
2932:             cOffset += dim - tDim;
2933:           }
2934:         }
2935:       }
2936:     };
2937:     // Update only the constrained dofs on a point
2938:     //   This takes an array with ALL values, not just BC
2939:     void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2940:       value_type *array = (value_type *) this->restrictPoint(p);
2941:       const int&  cDim  = this->getConstraintDimension(p);

2943:       if (cDim) {
2944:         if (orientation >= 0) {
2945:           const int&                          dim  = this->getFiberDimension(p);
2946:           const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2947:           int                                 cInd = 0;

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

2962:           for(int space = 0; space < this->getNumSpaces(); ++space) {
2963:             const int  dim = this->getFiberDimension(p, space);
2964:             const int tDim = this->getConstrainedFiberDimension(p, space);
2965:             int       cInd = 0;

2967:             for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2968:               if (cInd < 0) break;
2969:               if (j == cDof[cInd+cOffset]) {
2970:                 array[j] = v[k];
2971:                 ++cInd;
2972:               }
2973:             }
2974:             offset  += dim;
2975:             cOffset += dim - tDim;
2976:           }
2977:         }
2978:       }
2979:     };
2980:     // Update all dofs on a point (free and constrained)
2981:     void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
2982:       value_type *array = (value_type *) this->restrictPoint(p);

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

2987:         for(int i = 0; i < dim; ++i) {
2988:           array[i] = v[i];
2989:         }
2990:       } else {
2991:         int offset = 0;
2992:         int j      = -1;

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

2997:           for(int i = dim-1; i >= 0; --i) {
2998:             array[++j] = v[i+offset];
2999:           }
3000:           offset += dim;
3001:         }
3002:       }
3003:     };
3004:     // Update all dofs on a point (free and constrained)
3005:     void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
3006:       value_type *array = (value_type *) this->restrictPoint(p);

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

3011:         for(int i = 0; i < dim; ++i) {
3012:           array[i] += v[i];
3013:         }
3014:       } else {
3015:         int offset = 0;
3016:         int j      = -1;

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

3021:           for(int i = dim-1; i >= 0; --i) {
3022:             array[++j] += v[i+offset];
3023:           }
3024:           offset += dim;
3025:         }
3026:       }
3027:     };
3028:   public: // Fibrations
3029:     int getNumSpaces() const {return this->_spaces.size();};
3030:     const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
3031:     const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
3032:     void addSpace() {
3033:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
3034:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
3035:       this->_spaces.push_back(space);
3036:       this->_bcs.push_back(bc);
3037:     };
3038:     int getFiberDimension(const point_type& p, const int space) const {
3039:       return this->_spaces[space]->restrictPoint(p)->prefix;
3040:     };
3041:     void setFiberDimension(const point_type& p, int dim, const int space) {
3042:       const index_type idx(dim, -1);
3043:       this->_spaces[space]->addPoint(p);
3044:       this->_spaces[space]->updatePoint(p, &idx);
3045:     };
3046:     template<typename Sequence>
3047:     void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
3048:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
3049:         this->setFiberDimension(*p_iter, dim, space);
3050:       }
3051:     };
3052:     const int getConstraintDimension(const point_type& p, const int space) const {
3053:       if (!this->_bcs[space]->hasPoint(p)) return 0;
3054:       return this->_bcs[space]->getFiberDimension(p);
3055:     };
3056:     void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
3057:       this->_bcs[space]->setFiberDimension(p, numConstraints);
3058:     };
3059:     int getConstrainedFiberDimension(const point_type& p, const int space) const {
3060:       return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
3061:     };
3062:     void copyFibration(const Obj<FEMSection>& section) {
3063:       const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
3064:       const std::vector<Obj<bc_type> >&    bcs    = section->getBCs();

3066:       this->_spaces.clear();
3067:       for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
3068:         this->_spaces.push_back(*s_iter);
3069:       }
3070:       this->_bcs.clear();
3071:       for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
3072:         this->_bcs.push_back(*b_iter);
3073:       }
3074:     };
3075:     Obj<FEMSection> getFibration(const int space) const {
3076:       Obj<FEMSection> field = new FEMSection(this->comm(), this->debug());
3077: //     Obj<atlas_type> _atlas;
3078: //     std::vector<Obj<atlas_type> > _spaces;
3079: //     Obj<bc_type>    _bc;
3080: //     std::vector<Obj<bc_type> >    _bcs;
3081:       field->addSpace();
3082:       const chart_type& chart = this->getChart();

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

3089:         if (fDim) {
3090:           field->setFiberDimension(*c_iter, fDim);
3091:           field->setFiberDimension(*c_iter, fDim, 0);
3092:         }
3093:         if (cDim) {
3094:           field->setConstraintDimension(*c_iter, cDim);
3095:           field->setConstraintDimension(*c_iter, cDim, 0);
3096:         }
3097:       }
3098:       field->allocateStorage();
3099:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
3100:       const chart_type& newChart = field->getChart();

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

3106:         if (cDim) {
3107:           field->setConstraintDof(*c_iter, dof);
3108:         }
3109:       }
3110:       // Copy offsets
3111:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3112:         index_type idx;

3114:         idx.prefix = field->getFiberDimension(*c_iter);
3115:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
3116:         for(int s = 0; s < space; ++s) {
3117:           idx.index += this->getFiberDimension(*c_iter, s);
3118:         }
3119:         newAtlas->addPoint(*c_iter);
3120:         newAtlas->updatePoint(*c_iter, &idx);
3121:       }
3122:       field->replaceStorage(this->_array, true, this->getStorageSize());
3123:       field->setAtlas(newAtlas);
3124:       return field;
3125:     };
3126:   public:
3127:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3128:       ostringstream txt;
3129:       int rank;

3131:       if (comm == MPI_COMM_NULL) {
3132:         comm = this->comm();
3133:         rank = this->commRank();
3134:       } else {
3135:         MPI_Comm_rank(comm, &rank);
3136:       }
3137:       if (name == "") {
3138:         if(rank == 0) {
3139:           txt << "viewing a FEMSection" << std::endl;
3140:         }
3141:       } else {
3142:         if (rank == 0) {
3143:           txt << "viewing FEMSection '" << name << "'" << std::endl;
3144:         }
3145:       }
3146:       if (rank == 0) {
3147:         txt << "  Fields: " << this->getNumSpaces() << std::endl;
3148:       }
3149:       const chart_type& chart = this->getChart();

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

3155:         if (dim != 0) {
3156:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << "  ";
3157:           for(int i = 0; i < dim; i++) {
3158:             txt << " " << array[i];
3159:           }
3160:           const int& dim = this->getConstraintDimension(*p_iter);

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

3165:             txt << " constrained";
3166:             for(int i = 0; i < dim; ++i) {
3167:               txt << " " << bcArray[i];
3168:             }
3169:           }
3170:           txt << std::endl;
3171:         }
3172:       }
3173:       if (chart.size() == 0) {
3174:         txt << "[" << this->commRank() << "]: empty" << std::endl;
3175:       }
3176:       PetscSynchronizedPrintf(comm, txt.str().c_str());
3177:       PetscSynchronizedFlush(comm);
3178:     };
3179:   };
3180:   // A Field combines several sections
3181:   template<typename Overlap_, typename Patch_, typename Section_>
3182:   class Field : public ALE::ParallelObject {
3183:   public:
3184:     typedef Overlap_                                 overlap_type;
3185:     typedef Patch_                                   patch_type;
3186:     typedef Section_                                 section_type;
3187:     typedef typename section_type::point_type        point_type;
3188:     typedef typename section_type::chart_type        chart_type;
3189:     typedef typename section_type::value_type        value_type;
3190:     typedef std::map<patch_type, Obj<section_type> > sheaf_type;
3191:     typedef enum {SEND, RECEIVE}                     request_type;
3192:     typedef std::map<patch_type, MPI_Request>        requests_type;
3193:   protected:
3194:     sheaf_type    _sheaf;
3195:     int           _tag;
3196:     MPI_Datatype  _datatype;
3197:     requests_type _requests;
3198:   public:
3199:     Field(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
3200:       this->_tag      = this->getNewTag();
3201:       this->_datatype = this->getMPIDatatype();
3202:     };
3203:     Field(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _tag(tag) {
3204:       this->_datatype = this->getMPIDatatype();
3205:     };
3206:     virtual ~Field() {};
3207:   protected:
3208:     MPI_Datatype getMPIDatatype() {
3209:       if (sizeof(value_type) == 4) {
3210:         return MPI_INT;
3211:       } else if (sizeof(value_type) == 8) {
3212:         return MPI_DOUBLE;
3213:       } else if (sizeof(value_type) == 28) {
3214:         int          blen[2];
3215:         MPI_Aint     indices[2];
3216:         MPI_Datatype oldtypes[2], newtype;
3217:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
3218:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3219:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3220:         MPI_Type_commit(&newtype);
3221:         return newtype;
3222:       } else if (sizeof(value_type) == 32) {
3223:         int          blen[2];
3224:         MPI_Aint     indices[2];
3225:         MPI_Datatype oldtypes[2], newtype;
3226:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_DOUBLE;
3227:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3228:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3229:         MPI_Type_commit(&newtype);
3230:         return newtype;
3231:       }
3232:       throw ALE::Exception("Cannot determine MPI type for value type");
3233:     };
3234:     int getNewTag() {
3235:       static int tagKeyval = MPI_KEYVAL_INVALID;
3236:       int *tagvalp = NULL, *maxval, flg;

3238:       if (tagKeyval == MPI_KEYVAL_INVALID) {
3239:         tagvalp = (int *) malloc(sizeof(int));
3240:         MPI_Keyval_create(MPI_NULL_COPY_FN, Mesh_DelTag, &tagKeyval, (void *) NULL);
3241:         MPI_Attr_put(this->_comm, tagKeyval, tagvalp);
3242:         tagvalp[0] = 0;
3243:       }
3244:       MPI_Attr_get(this->_comm, tagKeyval, (void **) &tagvalp, &flg);
3245:       if (tagvalp[0] < 1) {
3246:         MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
3247:         tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
3248:       }
3249:       if (this->debug()) {
3250:         std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
3251:       }
3252:       return tagvalp[0]--;
3253:     };
3254:   public: // Verifiers
3255:     void checkPatch(const patch_type& patch) const {
3256:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3257:         ostringstream msg;
3258:         msg << "Invalid field patch " << patch << std::endl;
3259:         throw ALE::Exception(msg.str().c_str());
3260:       }
3261:     };
3262:     bool hasPatch(const patch_type& patch) {
3263:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3264:         return false;
3265:       }
3266:       return true;
3267:     };
3268:   public: // Accessors
3269:     int getTag() const {return this->_tag;};
3270:     void setTag(const int tag) {this->_tag = tag;};
3271:     Obj<section_type>& getSection(const patch_type& patch) {
3272:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3273:         this->_sheaf[patch] = new section_type(this->comm(), this->debug());
3274:       }
3275:       return this->_sheaf[patch];
3276:     };
3277:     void setSection(const patch_type& patch, const Obj<section_type>& section) {this->_sheaf[patch] = section;};
3278:     const sheaf_type& getPatches() {
3279:       return this->_sheaf;
3280:     };
3281:     void clear() {
3282:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3283:         p_iter->second->clear();
3284:       }
3285:     };
3286:   public: //  Adapter
3287:     template<typename Topology_>
3288:     void setTopology(const Obj<Topology_>& topology) {
3289:       const typename Topology_::sheaf_type& patches = topology->getPatches();

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

3296:         for(typename Topology_::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
3297:           section->setFiberDimension(*b_iter, 1);
3298:         }
3299:       }
3300:     };
3301:     void allocate() {
3302:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3303:         p_iter->second->allocatePoint();
3304:       }
3305:     };
3306:   public: // Communication
3307:     void construct(const int size) {
3308:       const sheaf_type& patches = this->getPatches();

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

3315:         for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3316:           section->setFiberDimension(*c_iter, size);
3317:         }
3318:       }
3319:     };
3320:     template<typename Sizer>
3321:     void construct(const Obj<Sizer>& sizer) {
3322:       const sheaf_type& patches = this->getPatches();

3324:       for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3325:         const patch_type         rank    = p_iter->first;
3326:         const Obj<section_type>& section = this->getSection(rank);
3327:         const chart_type&        chart   = section->getChart();
3328: 
3329:         for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3330:           section->setFiberDimension(*c_iter, *(sizer->getSection(rank)->restrictPoint(*c_iter)));
3331:         }
3332:       }
3333:     };
3334:     void constructCommunication(const request_type& requestType) {
3335:       const sheaf_type& patches = this->getPatches();

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

3342:         if (requestType == RECEIVE) {
3343:           if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data(" << section->size() << ") from " << patch << " tag " << this->_tag << std::endl;}
3344:           MPI_Recv_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3345:         } else {
3346:           if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << section->size() << ") to " << patch << " tag " << this->_tag << std::endl;}
3347:           MPI_Send_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3348:         }
3349:         this->_requests[patch] = request;
3350:       }
3351:     };
3352:     void startCommunication() {
3353:       const sheaf_type& patches = this->getPatches();

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

3358:         MPI_Start(&request);
3359:       }
3360:     };
3361:     void endCommunication() {
3362:       const sheaf_type& patches = this->getPatches();
3363:       MPI_Status status;

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

3368:         MPI_Wait(&request, &status);
3369:       }
3370:     };
3371:   public:
3372:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3373:       ostringstream txt;
3374:       int rank;

3376:       if (comm == MPI_COMM_NULL) {
3377:         comm = this->comm();
3378:         rank = this->commRank();
3379:       } else {
3380:         MPI_Comm_rank(comm, &rank);
3381:       }
3382:       if (name == "") {
3383:         if(rank == 0) {
3384:           txt << "viewing a Field" << std::endl;
3385:         }
3386:       } else {
3387:         if(rank == 0) {
3388:           txt << "viewing Field '" << name << "'" << std::endl;
3389:         }
3390:       }
3391:       PetscSynchronizedPrintf(comm, txt.str().c_str());
3392:       PetscSynchronizedFlush(comm);
3393:       for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3394:         ostringstream txt1;

3396:         txt1 << "[" << this->commRank() << "]: Patch " << p_iter->first << std::endl;
3397:         PetscSynchronizedPrintf(comm, txt1.str().c_str());
3398:         PetscSynchronizedFlush(comm);
3399:         p_iter->second->view("field section", comm);
3400:       }
3401:     };
3402:   };
3403: }
3404: #endif