Actual source code: IField.hh

  1: #ifndef included_ALE_IField_hh
  2: #define included_ALE_IField_hh

  4: #ifndef  included_ALE_Field_hh
  5: #include <Field.hh>
  6: #endif

  8: #ifndef  included_ALE_ISieve_hh
  9: #include <ISieve.hh>
 10: #endif

 12: // An ISection (or IntervalSection) is a section over a simple interval domain
 13: namespace ALE {
 14:   // An IConstantSection is the simplest ISection
 15:   //   All fibers are dimension 1
 16:   //   All values are equal to a constant
 17:   //     We need no value storage and no communication for completion
 18:   //     The default value is returned for any point not in the domain
 19:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
 20:   class IConstantSection : public ALE::ParallelObject {
 21:   public:
 22:     typedef Point_ point_type;
 23:     typedef Value_ value_type;
 24:     typedef Alloc_ alloc_type;
 25:     typedef Interval<point_type, alloc_type> chart_type;
 26:   protected:
 27:     chart_type _chart;
 28:     value_type _value[2]; // Value and default value
 29:   public:
 30:     IConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
 31:       _value[1] = 0;
 32:     };
 33:     IConstantSection(MPI_Comm comm, const point_type& min, const point_type& max, const value_type& value, const int debug) : ParallelObject(comm, debug), _chart(min, max) {
 34:       _value[0] = value;
 35:       _value[1] = value;
 36:     };
 37:     IConstantSection(MPI_Comm comm, const point_type& min, const point_type& max, const value_type& value, const value_type& defaultValue, const int debug) : ParallelObject(comm, debug), _chart(min, max) {
 38:       _value[0] = value;
 39:       _value[1] = defaultValue;
 40:     };
 41:   public: // Verifiers
 42:     void checkPoint(const point_type& point) const {
 43:       this->_chart.checkPoint(point);
 44:     };
 45:     void checkDimension(const int& dim) {
 46:       if (dim != 1) {
 47:         ostringstream msg;
 48:         msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
 49:         throw ALE::Exception(msg.str().c_str());
 50:       }
 51:     };
 52:     bool hasPoint(const point_type& point) const {
 53:       return this->_chart.hasPoint(point);
 54:     };
 55:   public: // Accessors
 56:     const chart_type& getChart() const {return this->_chart;};
 57:     void setChart(const chart_type& chart) {this->_chart = chart;};
 58:     void addPoint(const point_type& point) {
 59:       this->checkPoint(point);
 60:     };
 61:     template<typename Points>
 62:     void addPoint(const Obj<Points>& points) {
 63:       for(typename Points::const_iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
 64:         this->checkPoint(*p_iter);
 65:       }
 66:     };
 67:     template<typename Points>
 68:     void addPoint(const Points& points) {
 69:       for(typename Points::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
 70:         this->checkPoint(*p_iter);
 71:       }
 72:     };
 73:     value_type getDefaultValue() {return this->_value[1];};
 74:     void setDefaultValue(const value_type value) {this->_value[1] = value;};
 75:     void copy(const Obj<IConstantSection>& section) {
 76:       const chart_type& chart = section->getChart();

 78:       this->_chart = chart;
 79:       this->_value[0] = section->restrictPoint(*chart.begin())[0];
 80:       this->_value[1] = section->restrictPoint(*chart.begin())[1];
 81:     };
 82:   public: // Sizes
 83:     ///void clear() {};
 84:     int getFiberDimension(const point_type& p) const {
 85:       if (this->hasPoint(p)) return 1;
 86:       return 0;
 87:     };
 88:     void setFiberDimension(const point_type& p, int dim) {
 89:       this->checkDimension(dim);
 90:       this->addPoint(p);
 91:     };
 92:     template<typename Sequence>
 93:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
 94:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
 95:         this->setFiberDimension(*p_iter, dim);
 96:       }
 97:     };
 98:     void addFiberDimension(const point_type& p, int dim) {
 99:       if (this->hasPoint(p)) {
100:         ostringstream msg;
101:         msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
102:         throw ALE::Exception(msg.str().c_str());
103:       } else {
104:         this->setFiberDimension(p, dim);
105:       }
106:     };
107:     int size(const point_type& p) {return this->getFiberDimension(p);};
108:   public: // Restriction
109:     void clear() {};
110:     const value_type *restrictSpace() const {
111:       return this->_value;
112:     };
113:     const value_type *restrictPoint(const point_type& p) const {
114:       if (this->hasPoint(p)) {
115:         return this->_value;
116:       }
117:       return &this->_value[1];
118:     };
119:     void updatePoint(const point_type& p, const value_type v[]) {
120:       this->_value[0] = v[0];
121:     };
122:     void updateAddPoint(const point_type& p, const value_type v[]) {
123:       this->_value[0] += v[0];
124:     };
125:   public:
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 = this->comm();
132:         rank = this->commRank();
133:       } else {
134:         MPI_Comm_rank(comm, &rank);
135:       }
136:       if (name == "") {
137:         if(rank == 0) {
138:           txt << "viewing an IConstantSection" << std::endl;
139:         }
140:       } else {
141:         if(rank == 0) {
142:           txt << "viewing IConstantSection '" << name << "'" << std::endl;
143:         }
144:       }
145:       txt <<"["<<this->commRank()<<"]: chart " << this->_chart << std::endl;
146:       txt <<"["<<this->commRank()<<"]: Value " << this->_value[0] << " Default Value " << this->_value[1] << std::endl;
147:       PetscSynchronizedPrintf(comm, txt.str().c_str());
148:       PetscSynchronizedFlush(comm);
149:     };
150:   };

152:   // An IUniformSection often acts as an Atlas
153:   //   All fibers are the same dimension
154:   //     Note we can use a IConstantSection for this Atlas
155:   //   Each point may have a different vector
156:   //     Thus we need storage for values, and hence must implement completion
157:   template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
158:   class IUniformSection : public ALE::ParallelObject {
159:   public:
160:     typedef Point_                                                  point_type;
161:     typedef Value_                                                  value_type;
162:     typedef Alloc_                                                  alloc_type;
163:     typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
164:     typedef IConstantSection<point_type, int, point_alloc_type>     atlas_type;
165:     typedef typename atlas_type::chart_type                         chart_type;
166:     typedef point_type                                              index_type;
167:     typedef struct {value_type v[fiberDim];}                        fiber_type;
168:     typedef value_type*                                             values_type;
169:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
170:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
171:   protected:
172:     Obj<atlas_type> _atlas;
173:     values_type     _array;
174:     fiber_type      _emptyValue;
175:     alloc_type      _allocator;
176:   public:
177:     IUniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
178:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
179:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
180:       this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
181:       this->_array = NULL;
182:       for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
183:     };
184:     IUniformSection(MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : ParallelObject(comm, debug) {
185:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
186:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, min, max, fiberDim, debug));
187:       this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
188:       this->_array = NULL;
189:       for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
190:     };
191:     IUniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas) {
192:       int dim = fiberDim;
193:       this->_atlas->update(*this->_atlas->getChart().begin(), &dim);
194:       this->_array = NULL;
195:       for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
196:     };
197:     virtual ~IUniformSection() {
198:       if (this->_array) {
199:         for(int i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.destroy(this->_array+i);}
200:         this->_array += this->getChart().min()*fiberDim;
201:         this->_allocator.deallocate(this->_array, this->sizeWithBC());
202:         this->_array = NULL;
203:         this->_atlas = NULL;
204:       }
205:     };
206:   public:
207:     value_type *getRawArray(const int size) {
208:       static value_type *array   = NULL;
209:       static int         maxSize = 0;

211:       if (size > maxSize) {
212:         const value_type dummy(0);

214:         if (array) {
215:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
216:           this->_allocator.deallocate(array, maxSize);
217:         }
218:         maxSize = size;
219:         array   = this->_allocator.allocate(maxSize);
220:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
221:       };
222:       return array;
223:     };
224:   public: // Verifiers
225:     bool hasPoint(const point_type& point) const {
226:       return this->_atlas->hasPoint(point);
227:     };
228:     void checkDimension(const int& dim) {
229:       if (dim != fiberDim) {
230:         ostringstream msg;
231:         msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
232:         throw ALE::Exception(msg.str().c_str());
233:       }
234:     };
235:   public: // Accessors
236:     const chart_type& getChart() const {return this->_atlas->getChart();};
237:     void setChart(const chart_type& chart) {
238:       this->_atlas->setChart(chart);
239:       int dim = fiberDim;
240:       this->_atlas->updatePoint(*this->getChart().begin(), &dim);
241:     };
242:     bool resizeChart(const chart_type& chart) {
243:       if ((chart.min() >= this->getChart().min()) && (chart.max() <= this->getChart().max())) return false;
244:       this->setChart(chart);
245:       return true;
246:     };
247:     const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
248:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
249:     void addPoint(const point_type& point) {
250:       this->setFiberDimension(point, fiberDim);
251:     };
252:     template<typename Points>
253:     void addPoint(const Obj<Points>& points) {
254:       for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
255:         this->setFiberDimension(*p_iter, fiberDim);
256:       }
257:     };
258:     void copy(const Obj<IUniformSection>& section) {
259:       this->getAtlas()->copy(section->getAtlas());
260:       const chart_type& chart = section->getChart();

262:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
263:         this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
264:       }
265:     };
266:     const value_type *getDefault() const {return this->_emptyValue;};
267:     void setDefault(const value_type v[]) {for(int i = 0; i < fiberDim; ++i) {this->_emptyValue.v[i] = v[i];}};
268:   public: // Sizes
269:     void clear() {
270:       this->_atlas->clear();
271:     };
272:     int getFiberDimension(const point_type& p) const {
273:       return this->_atlas->restrictPoint(p)[0];
274:     };
275:     void setFiberDimension(const point_type& p, int dim) {
276:       this->checkDimension(dim);
277:       this->_atlas->addPoint(p);
278:       this->_atlas->updatePoint(p, &dim);
279:     };
280:     template<typename Sequence>
281:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
282:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
283:         this->setFiberDimension(*p_iter, dim);
284:       }
285:     };
286:     void setFiberDimension(const std::set<point_type>& points, int dim) {
287:       for(typename std::set<point_type>::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
288:         this->setFiberDimension(*p_iter, dim);
289:       }
290:     };
291:     void addFiberDimension(const point_type& p, int dim) {
292:       if (this->hasPoint(p)) {
293:         ostringstream msg;
294:         msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
295:         throw ALE::Exception(msg.str().c_str());
296:       } else {
297:         this->setFiberDimension(p, dim);
298:       }
299:     };
300:     int size() const {return this->_atlas->getChart().size()*fiberDim;};
301:     int sizeWithBC() const {return this->size();};
302:     void allocatePoint() {
303:       this->_array = this->_allocator.allocate(this->sizeWithBC());
304:       this->_array -= this->getChart().min()*fiberDim;
305:       for(index_type i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);}
306:     };
307:     bool reallocatePoint(const chart_type& chart, values_type *oldData = NULL) {
308:       const chart_type  oldChart = this->getChart();
309:       const int         oldSize  = this->sizeWithBC();
310:       values_type       oldArray = this->_array;
311:       if (!this->resizeChart(chart)) return false;
312:       const int         size     = this->sizeWithBC();

314:       this->_array = this->_allocator.allocate(size);
315:       this->_array -= this->getChart().min()*fiberDim;
316:       for(index_type i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);}
317:       for(int i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {
318:         this->_array[i] = oldArray[i];
319:       }
320:       if (!oldData) {
321:         for(index_type i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {this->_allocator.destroy(oldArray+i);}
322:         oldArray += this->getChart().min()*fiberDim;
323:         this->_allocator.deallocate(oldArray, oldSize);
324:         ///std::cout << "Freed IUniformSection data" << std::endl;
325:       } else {
326:         ///std::cout << "Did not free IUniformSection data" << std::endl;
327:         *oldData = oldArray;
328:       }
329:       return true;
330:     };
331:     template<typename Iterator, typename Extractor>
332:     bool reallocatePoint(const Iterator& begin, const Iterator& end, const Extractor& extractor) {
333:       point_type min = this->getChart().min();
334:       point_type max = this->getChart().max()-1;

336:       for(Iterator p_iter = begin; p_iter != end; ++p_iter) {
337:         min = std::min(extractor(*p_iter), min);
338:         max = std::max(extractor(*p_iter), max);
339:       }
340:       return reallocatePoint(chart_type(min, max+1));
341:     };
342:   public: // Restriction
343:     // Return a pointer to the entire contiguous storage array
344:     const values_type& restrictSpace() const {
345:       return this->_array;
346:     };
347:     // Return only the values associated to this point, not its closure
348:     const value_type *restrictPoint(const point_type& p) const {
349:       if (!this->hasPoint(p)) return this->_emptyValue.v;
350:       return &this->_array[p*fiberDim];
351:     };
352:     // Update only the values associated to this point, not its closure
353:     void updatePoint(const point_type& p, const value_type v[]) {
354:       for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
355:         this->_array[idx] = v[i];
356:       }
357:     };
358:     // Update only the values associated to this point, not its closure
359:     void updateAddPoint(const point_type& p, const value_type v[]) {
360:       for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
361:         this->_array[idx] += v[i];
362:       }
363:     };
364:     void updatePointAll(const point_type& p, const value_type v[]) {
365:       this->updatePoint(p, v);
366:     };
367:   public:
368:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
369:       ostringstream txt;
370:       int rank;

372:       if (comm == MPI_COMM_NULL) {
373:         comm = this->comm();
374:         rank = this->commRank();
375:       } else {
376:         MPI_Comm_rank(comm, &rank);
377:       }
378:       if (name == "") {
379:         if(rank == 0) {
380:           txt << "viewing an IUniformSection" << std::endl;
381:         }
382:       } else {
383:         if(rank == 0) {
384:           txt << "viewing IUniformSection '" << name << "'" << std::endl;
385:         }
386:       }
387:       const typename atlas_type::chart_type& chart = this->_atlas->getChart();
388:       values_type                            array = this->_array;

390:       for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
391:         const int idx = (*p_iter)*fiberDim;

393:         if (fiberDim != 0) {
394:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << fiberDim << "  ";
395:           for(int i = 0; i < fiberDim; i++) {
396:             txt << " " << array[idx+i];
397:           }
398:           txt << std::endl;
399:         }
400:       }
401:       if (chart.size() == 0) {
402:         txt << "[" << this->commRank() << "]: empty" << std::endl;
403:       }
404:       PetscSynchronizedPrintf(comm, txt.str().c_str());
405:       PetscSynchronizedFlush(comm);
406:     };
407:   };
408:   // An ISection allows variable fiber sizes per point
409:   //   The Atlas is a UniformSection of dimension 1 and value type Point
410:   //     to hold each fiber dimension and offsets into a contiguous array
411:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_> >
412:   class ISection : public Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > {
413:   public:
414:     typedef Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > base;
415:     typedef typename base::point_type       point_type;
416:     typedef typename base::value_type       value_type;
417:     typedef typename base::alloc_type       alloc_type;
418:     typedef typename base::index_type       index_type;
419:     typedef typename base::atlas_type       atlas_type;
420:     typedef typename base::chart_type       chart_type;
421:     typedef typename base::values_type      values_type;
422:     typedef typename base::atlas_alloc_type atlas_alloc_type;
423:     typedef typename base::atlas_ptr        atlas_ptr;
424:   public:
425:     ISection(MPI_Comm comm, const int debug = 0) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(comm, debug) {
426:     };
427:     ISection(MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(comm, debug) {
428:       this->_atlas->setChart(chart_type(min, max));
429:       this->_atlas->allocatePoint();
430:     };
431:     ISection(const Obj<atlas_type>& atlas) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(atlas) {};
432:     virtual ~ISection() {};
433:   public:
434:     void setChart(const chart_type& chart) {
435:       this->_atlas->setChart(chart);
436:       this->_atlas->allocatePoint();
437:     };
438:     bool resizeChart(const chart_type& chart) {
439:       if (!this->_atlas->reallocatePoint(chart)) return false;
440:       return true;
441:     };
442:     bool reallocatePoint(const chart_type& chart) {
443:       typedef typename atlas_type::alloc_type atlas_alloc_type;
444:       const chart_type        oldChart = this->getChart();
445:       const int               oldSize  = this->sizeWithBC();
446:       const values_type       oldArray = this->_array;
447:       const int               oldAtlasSize = this->_atlas->sizeWithBC();
448:       typename atlas_type::values_type oldAtlasArray;
449:       if (!this->_atlas->reallocatePoint(chart, &oldAtlasArray)) return false;

451:       this->orderPoints(this->_atlas);
452:       this->allocateStorage();
453:       for(int i = oldChart.min(); i < oldChart.max(); ++i) {
454:         const typename atlas_type::value_type& idx = this->_atlas->restrictPoint(i)[0];
455:         const int                              dim = idx.prefix;
456:         const int                              off = idx.index;

458:         for(int d = 0; d < dim; ++d) {
459:           this->_array[off+d] = oldArray[oldAtlasArray[i].index+d];
460:         }
461:       }
462:       for(int i = 0; i < oldSize; ++i) {this->_allocator.destroy(oldArray+i);}
463:       this->_allocator.deallocate(oldArray, oldSize);
464:       for(int i = oldChart.min(); i < oldChart.max(); ++i) {atlas_alloc_type(this->_allocator).destroy(oldAtlasArray+i);}
465:       oldAtlasArray += oldChart.min();
466:       atlas_alloc_type(this->_allocator).deallocate(oldAtlasArray, oldAtlasSize);
467:       ///std::cout << "In ISection, Freed IUniformSection data" << std::endl;
468:     };
469:   public:
470:     // Return the free values on a point
471:     //   This is overridden, because the one in Section cannot be const due to problem in the interface with UniformSection
472:     const value_type *restrictPoint(const point_type& p) const {
473:       return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
474:     };
475:   };
476:   // IGeneralSection will support BC on a subset of unknowns on a point
477:   //   We use a separate constraint Atlas to mark constrained dofs on a point
478:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_> >
479:   class IGeneralSection : public GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> > {
480:   public:
481:     typedef GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> > base;
482:     typedef typename base::point_type       point_type;
483:     typedef typename base::value_type       value_type;
484:     typedef typename base::alloc_type       alloc_type;
485:     typedef typename base::index_type       index_type;
486:     typedef typename base::atlas_type       atlas_type;
487:     typedef typename base::bc_type          bc_type;
488:     typedef typename base::chart_type       chart_type;
489:     typedef typename base::values_type      values_type;
490:     typedef typename base::atlas_alloc_type atlas_alloc_type;
491:     typedef typename base::atlas_ptr        atlas_ptr;
492:     typedef typename base::bc_alloc_type    bc_alloc_type;
493:     typedef typename base::bc_ptr           bc_ptr;
494:     typedef std::pair<point_type,int>       newpoint_type;
495:   protected:
496:     std::set<newpoint_type> newPoints;
497:   public:
498:     IGeneralSection(MPI_Comm comm, const int debug = 0) : GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> >(comm, debug) {};
499:     IGeneralSection(MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> >(comm, debug) {
500:       this->_atlas->setChart(chart_type(min, max));
501:       this->_atlas->allocatePoint();
502:       this->_bc->setChart(chart_type(min, max));
503:     };
504:     IGeneralSection(const Obj<atlas_type>& atlas) : GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> >(atlas) {
505:       this->_bc->setChart(atlas->getChart());
506:     };
507:     ~IGeneralSection() {};
508:   public:
509:     void setChart(const chart_type& chart) {
510:       this->_atlas->setChart(chart);
511:       this->_atlas->allocatePoint();
512:       this->_bc->setChart(chart);
513:       ///this->_bc->getAtlas()->allocatePoint();
514:       for(int s = 0; s < (int) this->_spaces.size(); ++s) {
515:         this->_spaces[s]->setChart(chart);
516:         this->_spaces[s]->allocatePoint();
517:         this->_bcs[s]->setChart(chart);
518:         ///this->_bcs[s]->getAtlas()->allocatePoint();
519:       }
520:     };
521:   public:
522:     bool hasNewPoints() {return this->newPoints.size() > 0;};
523:     const std::set<newpoint_type>& getNewPoints() {return this->newPoints;};
524:     void addPoint(const point_type& point, const int dim) {
525:       if (dim == 0) return;
526:       this->newPoints.insert(newpoint_type(point, dim));
527:     };
528:     // Returns true if the chart was changed
529:     bool resizeChart(const chart_type& chart) {
530:       if (!this->_atlas->reallocatePoint(chart)) return false;
531:       this->_bc->reallocatePoint(chart);
532:       for(int s = 0; s < (int) this->_spaces.size(); ++s) {
533:         this->_spaces[s]->reallocatePoint(chart);
534:         this->_bcs[s]->reallocatePoint(chart);
535:       }
536:       return true;
537:     };
538:     // Returns true if the chart was changed
539:     bool reallocatePoint(const chart_type& chart) {
540:       typedef typename alloc_type::template rebind<typename atlas_type::value_type>::other atlas_alloc_type;
541:       const chart_type        oldChart = this->getChart();
542:       const int               oldSize  = this->sizeWithBC();
543:       const values_type       oldArray = this->_array;
544:       const int               oldAtlasSize = this->_atlas->sizeWithBC();
545:       typename atlas_type::values_type oldAtlasArray;
546:       if (!this->_atlas->reallocatePoint(chart, &oldAtlasArray)) return false;
547:       this->_bc->reallocatePoint(chart);
548:       for(int s = 0; s < (int) this->_spaces.size(); ++s) {
549:         this->_spaces[s]->reallocatePoint(chart);
550:         this->_bcs[s]->reallocatePoint(chart);
551:       }
552:       for(typename std::set<newpoint_type>::const_iterator p_iter = this->newPoints.begin(); p_iter != this->newPoints.end(); ++p_iter) {
553:         this->setFiberDimension(p_iter->first, p_iter->second);
554:       }
555:       this->orderPoints(this->_atlas);
556:       this->allocateStorage();
557:       for(int i = oldChart.min(); i < oldChart.max(); ++i) {
558:         const typename atlas_type::value_type& idx = this->_atlas->restrictPoint(i)[0];
559:         const int                              dim = idx.prefix;
560:         const int                              off = idx.index;

562:         for(int d = 0; d < dim; ++d) {
563:           this->_array[off+d] = oldArray[oldAtlasArray[i].index+d];
564:         }
565:       }
566:       for(int i = 0; i < oldSize; ++i) {this->_allocator.destroy(oldArray+i);}
567:       this->_allocator.deallocate(oldArray, oldSize);
568:       for(int i = oldChart.min(); i < oldChart.max(); ++i) {atlas_alloc_type(this->_allocator).destroy(oldAtlasArray+i);}
569:       oldAtlasArray += oldChart.min();
570:       atlas_alloc_type(this->_allocator).deallocate(oldAtlasArray, oldAtlasSize);
571:       this->newPoints.clear();
572:       return true;
573:     };
574:   public:
575:     void addSpace() {
576:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
577:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
578:       space->setChart(this->_atlas->getChart());
579:       space->allocatePoint();
580:       bc->setChart(this->_bc->getChart());
581:       this->_spaces.push_back(space);
582:       this->_bcs.push_back(bc);
583:     };
584:     Obj<IGeneralSection> getFibration(const int space) const {
585:       Obj<IGeneralSection> field = new IGeneralSection(this->comm(), this->debug());
586: //     Obj<atlas_type> _atlas;
587: //     std::vector<Obj<atlas_type> > _spaces;
588: //     Obj<bc_type>    _bc;
589: //     std::vector<Obj<bc_type> >    _bcs;
590:       field->setChart(this->getChart());
591:       field->addSpace();
592:       const chart_type& chart = this->getChart();

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

599:         if (fDim) {
600:           field->setFiberDimension(*c_iter, fDim);
601:           field->setFiberDimension(*c_iter, fDim, 0);
602:         }
603:         if (cDim) {
604:           field->setConstraintDimension(*c_iter, cDim);
605:           field->setConstraintDimension(*c_iter, cDim, 0);
606:         }
607:       }
608:       field->allocateStorage();
609:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
610:       const chart_type& newChart = field->getChart();

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

616:         if (cDim) {
617:           field->setConstraintDof(*c_iter, dof);
618:         }
619:       }
620:       // Copy offsets
621:       newAtlas->setChart(newChart);
622:       newAtlas->allocatePoint();
623:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
624:         index_type idx;

626:         idx.prefix = field->getFiberDimension(*c_iter);
627:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
628:         for(int s = 0; s < space; ++s) {
629:           idx.index += this->getFiberDimension(*c_iter, s);
630:         }
631:         newAtlas->addPoint(*c_iter);
632:         newAtlas->updatePoint(*c_iter, &idx);
633:       }
634:       field->replaceStorage(this->_array, true, this->getStorageSize());
635:       field->setAtlas(newAtlas);
636:       return field;
637:     };
638:   };
639: }

641: #endif