Actual source code: IField.hh
petsc-3.4.2 2013-07-02

      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(int comp = 1) {
    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->_comps.push_back(comp);
    602:       this->_spaces.push_back(space);
    603:       this->_bcs.push_back(bc);
    604:     };
    605:     Obj<IGeneralSection> getFibration(const int space) const {
    606:       Obj<IGeneralSection> field = new IGeneralSection(this->comm(), this->debug());
    607: //     Obj<atlas_type> _atlas;
    608: //     std::vector<Obj<atlas_type> > _spaces;
    609: //     Obj<bc_type>    _bc;
    610: //     std::vector<Obj<bc_type> >    _bcs;
    611:       field->setChart(this->getChart());
    612:       field->addSpace();
    613:       const chart_type& chart = this->getChart();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1158: #endif