Actual source code: IField.hh

  1: #ifndef included_ALE_IField_hh
  2: #define included_ALE_IField_hh

  4: #ifndef  included_ALE_BasicCommunication_hh
  5: #include <sieve/BasicCommunication.hh>
  6: #endif

  8: #ifndef  included_ALE_Field_hh
  9: #include <sieve/Field.hh>
 10: #endif

 12: #ifndef  included_ALE_ISieve_hh
 13: #include <sieve/ISieve.hh>
 14: #endif

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

 85:       this->_chart = chart;
 86:       this->_value[0] = section->restrictPoint(*chart.begin())[0];
 87:       this->_value[1] = section->restrictPoint(*chart.begin())[1];
 88:     }
 89:   public: // Sizes
 90:     ///void clear() {};
 91:     int getFiberDimension(const point_type& p) const {
 92:       if (this->hasPoint(p)) return 1;
 93:       return 0;
 94:     }
 95:     void setFiberDimension(const point_type& p, int dim) {
 96:       this->checkDimension(dim);
 97:       this->addPoint(p);
 98:     }
 99:     template<typename Sequence>
100:     void setFiberDimension(const Obj<Sequence>& points, int dim) {
101:       const typename Sequence::iterator pointsEnd = points->end();
102:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != pointsEnd; ++p_iter) {
103:         this->setFiberDimension(*p_iter, dim);
104:       }
105:     }
106:     void addFiberDimension(const point_type& p, int dim) {
107:       if (this->hasPoint(p)) {
108:         ostringstream msg;
109:         msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
110:         throw ALE::Exception(msg.str().c_str());
111:       } else {
112:         this->setFiberDimension(p, dim);
113:       }
114:     }
115:     int size(const point_type& p) {return this->getFiberDimension(p);}
116:   public: // Restriction
117:     void clear() {};
118:     const value_type *restrictSpace() const {
119:       return this->_value;
120:     };
121:     const value_type *restrictPoint(const point_type& p) const {
122:       if (this->hasPoint(p)) {
123:         return this->_value;
124:       }
125:       return &this->_value[1];
126:     };
127:     void updatePoint(const point_type&, const value_type v[]) {
128:       this->_value[0] = v[0];
129:     };
130:     void updateAddPoint(const point_type&, const value_type v[]) {
131:       this->_value[0] += v[0];
132:     };
133:   public:
134:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
135:       ostringstream txt;
136:       int rank;

138:       if (comm == MPI_COMM_NULL) {
139:         comm = this->comm();
140:         rank = this->commRank();
141:       } else {
142:         MPI_Comm_rank(comm, &rank);
143:       }
144:       if (name == "") {
145:         if(rank == 0) {
146:           txt << "viewing an IConstantSection" << std::endl;
147:         }
148:       } else {
149:         if(rank == 0) {
150:           txt << "viewing IConstantSection '" << name << "'" << std::endl;
151:         }
152:       }
153:       txt <<"["<<this->commRank()<<"]: chart " << this->_chart << std::endl;
154:       txt <<"["<<this->commRank()<<"]: Value " << this->_value[0] << " Default Value " << this->_value[1] << std::endl;
155:       PetscSynchronizedPrintf(comm, txt.str().c_str());
156:       PetscSynchronizedFlush(comm);
157:     };
158:   };

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

219:       if (size > maxSize) {
220:         const value_type dummy(0);

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

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

327:       this->_array = this->_allocator.allocate(size);
328:       this->_array -= this->getChart().min()*fiberDim;
329:       for(index_type i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);}
330:       for(int i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {
331:         this->_array[i] = oldArray[i];
332:       }
333:       if (!oldData) {
334:         for(index_type i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {this->_allocator.destroy(oldArray+i);}
335:         oldArray += this->getChart().min()*fiberDim;
336:         this->_allocator.deallocate(oldArray, oldSize);
337:         ///std::cout << "Freed IUniformSection data" << std::endl;
338:       } else {
339:         ///std::cout << "Did not free IUniformSection data" << std::endl;
340:         *oldData = oldArray;
341:       }
342:       return true;
343:     };
344:     template<typename Iterator, typename Extractor>
345:     bool reallocatePoint(const Iterator& begin, const Iterator& end, const Extractor& extractor) {
346:       point_type min = this->getChart().min();
347:       point_type max = this->getChart().max()-1;

349:       for(Iterator p_iter = begin; p_iter != end; ++p_iter) {
350:         min = std::min(extractor(*p_iter), min);
351:         max = std::max(extractor(*p_iter), max);
352:       }
353:       return reallocatePoint(chart_type(min, max+1));
354:     }
355:   public: // Restriction
356:     // Zero entries
357:     void zero() {
358:       memset(this->_array+(this->getChart().min()*fiberDim), 0, this->sizeWithBC()* sizeof(value_type));
359:     };
360:     // Return a pointer to the entire contiguous storage array
361:     const values_type& restrictSpace() const {
362:       return this->_array;
363:     };
364:     // Return only the values associated to this point, not its closure
365:     const value_type *restrictPoint(const point_type& p) const {
366:       if (!this->hasPoint(p)) return this->_emptyValue.v;
367:       return &this->_array[p*fiberDim];
368:     };
369:     // Update only the values associated to this point, not its closure
370:     void updatePoint(const point_type& p, const value_type v[]) {
371:       for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
372:         this->_array[idx] = v[i];
373:       }
374:     };
375:     // Update only the values associated to this point, not its closure
376:     void updateAddPoint(const point_type& p, const value_type v[]) {
377:       for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
378:         this->_array[idx] += v[i];
379:       }
380:     };
381:     void updatePointAll(const point_type& p, const value_type v[]) {
382:       this->updatePoint(p, v);
383:     };
384:   public:
385:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
386:       ostringstream txt;
387:       int rank;

389:       if (comm == MPI_COMM_NULL) {
390:         comm = this->comm();
391:         rank = this->commRank();
392:       } else {
393:         MPI_Comm_rank(comm, &rank);
394:       }
395:       if (name == "") {
396:         if(rank == 0) {
397:           txt << "viewing an IUniformSection" << std::endl;
398:         }
399:       } else {
400:         if(rank == 0) {
401:           txt << "viewing IUniformSection '" << name << "'" << std::endl;
402:         }
403:       }
404:       const typename atlas_type::chart_type& chart = this->_atlas->getChart();
405:       values_type                            array = this->_array;

407:       const typename atlas_type::chart_type::const_iterator chartEnd = chart.end();
408:       for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chartEnd; ++p_iter) {
409:         const int idx = (*p_iter)*fiberDim;

411:         if (fiberDim != 0) {
412:           txt << "[" << this->commRank() << "]:   " << *p_iter << " dim " << fiberDim << "  ";
413:           for(int i = 0; i < fiberDim; i++) {
414:             txt << " " << array[idx+i];
415:           }
416:           txt << std::endl;
417:         }
418:       }
419:       if (chart.size() == 0) {
420:         txt << "[" << this->commRank() << "]: empty" << std::endl;
421:       }
422:       PetscSynchronizedPrintf(comm, txt.str().c_str());
423:       PetscSynchronizedFlush(comm);
424:     };
425:   };
426:   // An ISection allows variable fiber sizes per point
427:   //   The Atlas is a UniformSection of dimension 1 and value type Point
428:   //     to hold each fiber dimension and offsets into a contiguous array
429:   template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_> >
430:   class ISection : public Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > {
431:   public:
432:     typedef Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > base;
433:     typedef typename base::point_type       point_type;
434:     typedef typename base::value_type       value_type;
435:     typedef typename base::alloc_type       alloc_type;
436:     typedef typename base::index_type       index_type;
437:     typedef typename base::atlas_type       atlas_type;
438:     typedef typename base::chart_type       chart_type;
439:     typedef typename base::values_type      values_type;
440:     typedef typename base::atlas_alloc_type atlas_alloc_type;
441:     typedef typename base::atlas_ptr        atlas_ptr;
442:   public:
443:     ISection(MPI_Comm comm, const int debug = 0) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(comm, debug) {
444:     };
445:     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) {
446:       this->_atlas->setChart(chart_type(min, max));
447:       this->_atlas->allocatePoint();
448:     };
449:     ISection(const Obj<atlas_type>& atlas) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(atlas) {};
450:     virtual ~ISection() {};
451:   public:
452:     void setChart(const chart_type& chart) {
453:       this->_atlas->setChart(chart);
454:       this->_atlas->allocatePoint();
455:     };
456:     bool resizeChart(const chart_type& chart) {
457:       if (!this->_atlas->reallocatePoint(chart)) return false;
458:       return true;
459:     };
460:     bool reallocatePoint(const chart_type& chart) {
461:       typedef typename atlas_type::alloc_type atlas_alloc_type;
462:       const chart_type        oldChart = this->getChart();
463:       const int               oldSize  = this->sizeWithBC();
464:       const values_type       oldArray = this->_array;
465:       const int               oldAtlasSize = this->_atlas->sizeWithBC();
466:       typename atlas_type::values_type oldAtlasArray;
467:       if (!this->_atlas->reallocatePoint(chart, &oldAtlasArray)) return false;

469:       this->orderPoints(this->_atlas);
470:       this->allocateStorage();
471:       for(int i = oldChart.min(); i < oldChart.max(); ++i) {
472:         const typename atlas_type::value_type& idx = this->_atlas->restrictPoint(i)[0];
473:         const int                              dim = idx.prefix;
474:         const int                              off = idx.index;

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

582:         for(int d = 0; d < dim; ++d) {
583:           this->_array[off+d] = oldArray[oldAtlasArray[i].index+d];
584:         }
585:       }
586:       for(int i = 0; i < oldSize; ++i) {this->_allocator.destroy(oldArray+i);}
587:       this->_allocator.deallocate(oldArray, oldSize);
588:       for(int i = oldChart.min(); i < oldChart.max(); ++i) {atlas_alloc_type(this->_allocator).destroy(oldAtlasArray+i);}
589:       oldAtlasArray += oldChart.min();
590:       atlas_alloc_type(this->_allocator).deallocate(oldAtlasArray, oldAtlasSize);
591:       this->newPoints.clear();
592:       return true;
593:     };
594:   public:
595:     void addSpace() {
596:       Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
597:       Obj<bc_type>    bc    = new bc_type(this->comm(), this->debug());
598:       space->setChart(this->_atlas->getChart());
599:       space->allocatePoint();
600:       bc->setChart(this->_bc->getChart());
601:       this->_spaces.push_back(space);
602:       this->_bcs.push_back(bc);
603:     };
604:     Obj<IGeneralSection> getFibration(const int space) const {
605:       Obj<IGeneralSection> field = new IGeneralSection(this->comm(), this->debug());
606: //     Obj<atlas_type> _atlas;
607: //     std::vector<Obj<atlas_type> > _spaces;
608: //     Obj<bc_type>    _bc;
609: //     std::vector<Obj<bc_type> >    _bcs;
610:       field->setChart(this->getChart());
611:       field->addSpace();
612:       const chart_type& chart = this->getChart();

614:       // Copy sizes
615:       const typename chart_type::const_iterator chartEnd = chart.end();
616:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
617:         const int fDim = this->getFiberDimension(*c_iter, space);
618:         const int cDim = this->getConstraintDimension(*c_iter, space);

620:         if (fDim) {
621:           field->setFiberDimension(*c_iter, fDim);
622:           field->setFiberDimension(*c_iter, fDim, 0);
623:         }
624:         if (cDim) {
625:           field->setConstraintDimension(*c_iter, cDim);
626:           field->setConstraintDimension(*c_iter, cDim, 0);
627:         }
628:       }
629:       field->allocateStorage();
630:       Obj<atlas_type>   newAtlas = new atlas_type(this->comm(), this->debug());
631:       const chart_type& newChart = field->getChart();

633:       const typename chart_type::const_iterator newChartEnd = newChart.end();
634:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChartEnd; ++c_iter) {
635:         const int cDim   = field->getConstraintDimension(*c_iter);
636:         const int dof[1] = {0};

638:         if (cDim) {
639:           field->setConstraintDof(*c_iter, dof);
640:         }
641:       }
642:       // Copy offsets
643:       newAtlas->setChart(newChart);
644:       newAtlas->allocatePoint();
645:       for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChartEnd; ++c_iter) {
646:         index_type idx;

648:         idx.prefix = field->getFiberDimension(*c_iter);
649:         idx.index  = this->_atlas->restrictPoint(*c_iter)[0].index;
650:         for(int s = 0; s < space; ++s) {
651:           idx.index += this->getFiberDimension(*c_iter, s);
652:         }
653:         newAtlas->addPoint(*c_iter);
654:         newAtlas->updatePoint(*c_iter, &idx);
655:       }
656:       field->replaceStorage(this->_array, true, this->getStorageSize());
657:       field->setAtlas(newAtlas);
658:       return field;
659:     };
660:   };

662:   class SectionSerializer {
663:   public:
664:     template<typename Point_, typename Value_>
665:     static void writeSection(std::ofstream& fs, IConstantSection<Point_, Value_>& section) {
666:       MPIMover<Value_> mover(section.comm());

668:       if (section.commRank() == 0) {
669:         // Write local
670:         fs << section.getChart().min() << " " << section.getChart().max() << std::endl;
671:         fs.precision(15);
672:         fs << section.restrictPoint(section.getChart().min())[0] << " ";
673:         fs << section.getDefaultValue() << std::endl;
674:         // Receive and write remote
675:         for(int p = 1; p < section.commSize(); ++p) {
676:           PetscInt       sizes[2];
677:           Value_         values[2];
678:           MPI_Status     status;

681:           MPI_Recv(sizes,  2, MPIU_INT,        p, 1, section.comm(), &status);CHKERRXX(ierr);
682:           MPI_Recv(values, 2, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
683:           fs << sizes[0] << " " << sizes[1] << std::endl;
684:           fs.precision(15);
685:           fs << values[0] << " " << values[1] << std::endl;
686:         }
687:       } else {
688:         // Send remote
689:         PetscInt       sizes[2];
690:         Value_         values[2];

693:         sizes[0]  = section.getChart().min();
694:         sizes[1]  = section.getChart().max();
695:         values[0] = section.restrictPoint(section.getChart().min())[0];
696:         values[1] = section.getDefaultValue();
697:         MPI_Send(sizes,  2, MPIU_INT,        0, 1, section.comm());CHKERRXX(ierr);
698:         MPI_Send(values, 2, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
699:       }
700:     }
701:     template<typename Point_, typename Value_, int fiberDim>
702:     static void writeSection(std::ofstream& fs, IUniformSection<Point_, Value_, fiberDim>& section) {
703:       typedef typename IUniformSection<Point_, Value_, fiberDim>::index_type index_type;
704:       typedef typename IUniformSection<Point_, Value_, fiberDim>::value_type value_type;
705:       MPIMover<value_type> mover(section.comm());
706:       index_type min = section.getChart().min();
707:       index_type max = section.getChart().max();

709:       // Write atlas
710:       writeSection(fs, *section.getAtlas());
711:       if (section.commRank() == 0) {
712:         // Write local values
713:         fs.precision(15);
714:         for(index_type p = min; p < max; ++p) {
715:           const value_type *values = section.restrictPoint(p);

717:           for(int i = 0; i < fiberDim; ++i) {
718:             fs << values[i] << std::endl;
719:           }
720:         }
721:         // Write empty value
722:         const value_type *defValue = section.getDefault();

724:         for(int i = 0; i < fiberDim; ++i) {
725:           if (i > 0) fs << " ";
726:           fs << defValue[i];
727:         }
728:         fs << std::endl;
729:         // Receive and write remote
730:         for(int p = 1; p < section.commSize(); ++p) {
731:           PetscInt       size;
732:           value_type    *values;
733:           value_type     emptyValues[fiberDim];
734:           MPI_Status     status;

737:           MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
738:           PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
739:           MPI_Recv(values, size*fiberDim, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
740:           for(PetscInt v = 0; v < size; ++v) {
741:             fs << values[v] << std::endl;
742:           }
743:           PetscFree(values);CHKERRXX(ierr);
744:           MPI_Recv(emptyValues, fiberDim, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
745:           for(int i = 0; i < fiberDim; ++i) {
746:             if (i > 0) fs << " ";
747:             fs << emptyValues[i];
748:           }
749:           fs << std::endl;
750:         }
751:       } else {
752:         // Send remote
753:         PetscInt          size = section.getChart().size();
754:         PetscInt          v    = 0;
755:         const value_type *defValue = section.getDefault();
756:         value_type       *values;
757:         value_type        emptyValues[fiberDim];
758:         PetscErrorCode    ierr;

760:         MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
761:         PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
762:         for(index_type p = min; p < max; ++p) {
763:           const value_type *val = section.restrictPoint(p);

765:           for(int i = 0; i < fiberDim; ++i, ++v) {
766:             values[v] = val[i];
767:           }
768:         }
769:         MPI_Send(values, size*fiberDim, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
770:         for(int i = 0; i < fiberDim; ++i) {emptyValues[i] = ((value_type *) &defValue[i])[0];}
771:         MPI_Send(emptyValues, fiberDim, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
772:       }
773:     }
774:     template<typename Point_, typename Value_>
775:     static void writeSection(std::ofstream& fs, ISection<Point_, Value_>& section) {
776:       typedef typename ISection<Point_, Value_>::point_type point_type;
777:       typedef typename ISection<Point_, Value_>::value_type value_type;
778:       MPIMover<value_type> mover(section.comm());
779:       point_type min = section.getChart().min();
780:       point_type max = section.getChart().max();

782:       // Write atlas
783:       writeSection(fs, *section.getAtlas());
784:       if (section.commRank() == 0) {
785:       // Write local values
786:         fs.precision(15);
787:         for(point_type p = min; p < max; ++p) {
788:           const int         fiberDim = section.getFiberDimension(p);
789:           const value_type *values   = section.restrictPoint(p);

791:           for(int i = 0; i < fiberDim; ++i) {
792:             fs << values[i] << std::endl;
793:           }
794:         }
795:         // Receive and write remote
796:         for(int p = 1; p < section.commSize(); ++p) {
797:           PetscInt       size;
798:           value_type    *values;
799:           MPI_Status     status;

802:           MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
803:           PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
804:           MPI_Recv(values, size, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
805:           for(PetscInt v = 0; v < size; ++v) {
806:             fs << values[v] << std::endl;
807:           }
808:           PetscFree(values);CHKERRXX(ierr);
809:         }
810:       } else {
811:         // Send remote
812:         PetscInt       size = section.size();
813:         PetscInt       v    = 0;
814:         value_type    *values;

817:         MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
818:         PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
819:         for(point_type p = min; p < max; ++p) {
820:           const int         fiberDim = section.getFiberDimension(p);
821:           const value_type *val      = section.restrictPoint(p);

823:           for(int i = 0; i < fiberDim; ++i, ++v) {
824:             values[v] = val[i];
825:           }
826:         }
827:         MPI_Send(values, size, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
828:       }
829:     }
830:     template<typename Point_, typename Value_>
831:     static void writeSection(std::ofstream& fs, IGeneralSection<Point_, Value_>& section) {
832:       typedef typename IGeneralSection<Point_, Value_>::point_type point_type;
833:       typedef typename IGeneralSection<Point_, Value_>::value_type value_type;
834:       MPIMover<value_type> mover(section.comm());
835:       point_type min = section.getChart().min();
836:       point_type max = section.getChart().max();

838:       // Write atlas
839:       writeSection(fs, *section.getAtlas());
840:       if (section.commRank() == 0) {
841:         // Write local values
842:         fs.precision(15);
843:         for(point_type p = min; p < max; ++p) {
844:           const int         fiberDim = section.getFiberDimension(p);
845:           const value_type *values   = section.restrictPoint(p);

847:           for(int i = 0; i < fiberDim; ++i) {
848:             fs << values[i] << std::endl;
849:           }
850:         }
851:         // Receive and write remote
852:         for(int p = 1; p < section.commSize(); ++p) {
853:           PetscInt       size;
854:           value_type    *values;
855:           MPI_Status     status;

858:           MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
859:           PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
860:           MPI_Recv(values, size, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
861:           for(PetscInt v = 0; v < size; ++v) {
862:             fs << values[v] << std::endl;
863:           }
864:           PetscFree(values);CHKERRXX(ierr);
865:         }
866:       } else {
867:         // Send remote
868:         PetscInt       size = section.sizeWithBC();
869:         PetscInt       v    = 0;
870:         value_type    *values;

873:         MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
874:         PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
875:         for(point_type p = min; p < max; ++p) {
876:           const int         fiberDim = section.getFiberDimension(p);
877:           const value_type *val      = section.restrictPoint(p);

879:           for(int i = 0; i < fiberDim; ++i, ++v) {
880:             values[v] = val[i];
881:           }
882:         }
883:         MPI_Send(values, size, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
884:       }
885:       // Write BC
886:       writeSection(fs, *section.getBC());
887:       // Write spaces
888:       //   std::vector<Obj<atlas_type> > _spaces;
889:       //   std::vector<Obj<bc_type> >    _bcs;
890:     }
891:     template<typename Point_, typename Value_>
892:     static void loadSection(std::ifstream& fs, IConstantSection<Point_, Value_>& section) {
893:       typedef typename IConstantSection<Point_, Value_>::index_type index_type;
894:       typedef typename IConstantSection<Point_, Value_>::value_type value_type;
895:       MPIMover<value_type> mover(section.comm());
896:       index_type min, max;
897:       value_type val;

899:       if (section.commRank() == 0) {
900:         // Load local
901:         fs >> min;
902:         fs >> max;
903:         section.setChart(typename IConstantSection<Point_, Value_>::chart_type(min, max));
904:         fs >> val;
905:         section.updatePoint(min, &val);
906:         fs >> val;
907:         section.setDefaultValue(val);
908:         // Load and send remote
909:         for(int p = 1; p < section.commSize(); ++p) {
910:           PetscInt       sizes[2];
911:           value_type     values[2];

914:           fs >> sizes[0];
915:           fs >> sizes[1];
916:           fs >> values[0];
917:           fs >> values[1];
918:           MPI_Send(sizes,  2, MPIU_INT,        p, 1, section.comm());CHKERRXX(ierr);
919:           MPI_Send(values, 2, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
920:         }
921:       } else {
922:         // Load remote
923:         PetscInt       sizes[2];
924:         value_type     values[2];
925:         MPI_Status     status;

928:         MPI_Recv(sizes,  2, MPIU_INT,    0, 1, section.comm(), &status);CHKERRXX(ierr);
929:         section.setChart(typename IConstantSection<Point_, Value_>::chart_type(sizes[0], sizes[1]));
930:         MPI_Recv(values, 2, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
931:         section.updatePoint(min, values);
932:         section.setDefaultValue(values[1]);
933:       }
934:     }
935:     template<typename Point_, typename Value_, int fiberDim>
936:     static void loadSection(std::ifstream& fs, IUniformSection<Point_, Value_, fiberDim>& section) {
937:       typedef typename IUniformSection<Point_, Value_, fiberDim>::index_type index_type;
938:       typedef typename IUniformSection<Point_, Value_, fiberDim>::value_type value_type;
939:       MPIMover<value_type> mover(section.comm());
940:       // Load atlas
941:       loadSection(fs, *section.getAtlas());
942:       section.allocatePoint();
943:       index_type min = section.getChart().min();
944:       index_type max = section.getChart().max();

946:       if (section.commRank() == 0) {
947:         // Load local values
948:         for(index_type p = min; p < max; ++p) {
949:           value_type values[fiberDim];

951:           for(int i = 0; i < fiberDim; ++i) {
952:             typename IUniformSection<Point_, Value_, fiberDim>::value_type value;

954:             fs >> value;
955:             values[i] = value;
956:           }
957:           section.updatePoint(p, values);
958:         }
959:         // Load empty value
960:         value_type defValue[fiberDim];

962:         for(int i = 0; i < fiberDim; ++i) {
963:           fs >> defValue[i];
964:         }
965:         section.setDefault(defValue);
966:         // Load and send remote
967:         for(int pr = 1; pr < section.commSize(); ++pr) {
968:           PetscInt       size = section.getChart().size();
969:           PetscInt       v    = 0;
970:           value_type    *values;
971:           value_type     emptyValues[fiberDim];

974:           MPI_Send(&size, 1, MPIU_INT, pr, 1, section.comm());CHKERRXX(ierr);
975:           PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
976:           for(index_type p = min; p < max; ++p) {
977:             for(int i = 0; i < fiberDim; ++i, ++v) {
978:               fs >> values[v];
979:             }
980:           }
981:           MPI_Send(values, size*fiberDim, mover.getType(), pr, 1, section.comm());CHKERRXX(ierr);
982:           for(int i = 0; i < fiberDim; ++i) {
983:             fs >> emptyValues[i];
984:           }
985:           MPI_Send(emptyValues, fiberDim, mover.getType(), pr, 1, section.comm());CHKERRXX(ierr);
986:         }
987:       } else {
988:         // Load remote
989:         PetscInt       size;
990:         value_type    *values;
991:         value_type     emptyValues[fiberDim];
992:         value_type     pvalues[fiberDim];
993:         MPI_Status     status;
994:         PetscInt       v = 0;

997:         assert(sizeof(value_type) <= sizeof(PetscScalar));
998:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
999:         PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
1000:         MPI_Recv(values, size*fiberDim, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1001:         for(index_type p = min; p < max; ++p) {
1002:           for(int i = 0; i < fiberDim; ++i, ++v) {
1003:             pvalues[i] = values[v];
1004:           }
1005:           section.updatePoint(p, pvalues);
1006:         }
1007:         PetscFree(values);CHKERRXX(ierr);
1008:         MPI_Recv(emptyValues, fiberDim, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1009:         section.setDefault(emptyValues);
1010:       }
1011:     }
1012:     template<typename Point_, typename Value_>
1013:     static void loadSection(std::ifstream& fs, ISection<Point_, Value_>& section) {
1014:       typedef typename ISection<Point_, Value_>::point_type point_type;
1015:       typedef typename ISection<Point_, Value_>::value_type value_type;
1016:       MPIMover<value_type> mover(section.comm());
1017:       // Load atlas
1018:       loadSection(fs, *section.getAtlas());
1019:       section.allocatePoint();
1020:       point_type min    = section.getChart().min();
1021:       point_type max    = section.getChart().max();
1022:       int        maxDim = -1;

1024:       if (section.commRank() == 0) {
1025:         // Load local values
1026:         for(point_type p = min; p < max; ++p) {
1027:           maxDim = std::max(maxDim, section.getFiberDimension(p));
1028:         }
1029:         value_type *values = new value_type[maxDim];
1030:         for(point_type p = min; p < max; ++p) {
1031:           const int fiberDim = section.getFiberDimension(p);

1033:           for(int i = 0; i < fiberDim; ++i) {
1034:             fs >> values[i];
1035:           }
1036:           section.updatePoint(p, values);
1037:         }
1038:         delete [] values;
1039:         // Load and send remote
1040:         for(int p = 1; p < section.commSize(); ++p) {
1041:           PetscInt       size = section.size();
1042:           value_type    *values;

1045:           MPI_Send(&size, 1, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
1046:           PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1047:           for(PetscInt v = 0; v < size; ++v) {
1048:             fs >> values[v];
1049:           }
1050:           MPI_Send(values, size, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
1051:         }
1052:       } else {
1053:         // Load remote
1054:         PetscInt       size;
1055:         value_type    *values;
1056:         MPI_Status     status;
1057:         PetscInt       maxDim = 0;
1058:         PetscInt       off    = 0;

1061:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
1062:         PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1063:         MPI_Recv(values, size, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1064:         for(point_type p = min; p < max; ++p) {
1065:           maxDim = std::max(maxDim, section.getFiberDimension(p));
1066:         }
1067:         value_type *pvalues = new value_type[maxDim];
1068:         for(point_type p = min; p < max; ++p) {
1069:           const int fiberDim = section.getFiberDimension(p);

1071:           for(int i = 0; i < fiberDim; ++i, ++off) {
1072:             pvalues[i] = values[off];
1073:           }
1074:           section.updatePoint(p, pvalues);
1075:         }
1076:         delete [] pvalues;
1077:         PetscFree(values);CHKERRXX(ierr);
1078:       }
1079:     }
1080:     template<typename Point_, typename Value_>
1081:     static void loadSection(std::ifstream& fs, IGeneralSection<Point_, Value_>& section) {
1082:       typedef typename IGeneralSection<Point_, Value_>::point_type point_type;
1083:       typedef typename IGeneralSection<Point_, Value_>::value_type value_type;
1084:       MPIMover<value_type> mover(section.comm());
1085:       // Load atlas
1086:       loadSection(fs, *section.getAtlas());
1087:       section.allocatePoint();
1088:       point_type min    = section.getChart().min();
1089:       point_type max    = section.getChart().max();
1090:       int        maxDim = -1;

1092:       if (section.commRank() == 0) {
1093:         // Load local values
1094:         for(point_type p = min; p < max; ++p) {
1095:           maxDim = std::max(maxDim, section.getFiberDimension(p));
1096:         }
1097:         value_type *values = new value_type[maxDim];
1098:         for(point_type p = min; p < max; ++p) {
1099:           const int fiberDim = section.getFiberDimension(p);

1101:           for(int i = 0; i < fiberDim; ++i) {
1102:             fs >> values[i];
1103:           }
1104:           section.updatePoint(p, values);
1105:         }
1106:         delete [] values;
1107:         // Load and send remote
1108:         for(int p = 1; p < section.commSize(); ++p) {
1109:           PetscInt       size = section.sizeWithBC();
1110:           value_type    *values;

1113:           MPI_Send(&size, 1, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
1114:           PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1115:           for(PetscInt v = 0; v < size; ++v) {
1116:             fs >> values[v];
1117:           }
1118:           MPI_Send(values, size, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
1119:         }
1120:       } else {
1121:         // Load remote
1122:         PetscInt       size;
1123:         value_type    *values;
1124:         MPI_Status     status;
1125:         PetscInt       maxDim = 0;
1126:         PetscInt       off    = 0;

1129:         assert(sizeof(value_type) <= sizeof(PetscScalar));
1130:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
1131:         PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1132:         MPI_Recv(values, size, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1133:         for(point_type p = min; p < max; ++p) {
1134:           maxDim = std::max(maxDim, section.getFiberDimension(p));
1135:         }
1136:         value_type *pvalues = new value_type[maxDim];
1137:         for(point_type p = min; p < max; ++p) {
1138:           const int fiberDim = section.getFiberDimension(p);

1140:           for(int i = 0; i < fiberDim; ++i, ++off) {
1141:             pvalues[i] = values[off];
1142:           }
1143:           section.updatePoint(p, pvalues);
1144:         }
1145:         delete [] pvalues;
1146:         PetscFree(values);CHKERRXX(ierr);
1147:       }
1148:       // Load BC
1149:       loadSection(fs, *section.getBC());
1150:       // Load spaces
1151:       //   std::vector<Obj<atlas_type> > _spaces;
1152:       //   std::vector<Obj<bc_type> >    _bcs;
1153:     }
1154:   };
1155: }

1157: #endif