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

      1: #ifndef included_ALE_CartesianSieve_hh
      2: #define included_ALE_CartesianSieve_hh

      4: #ifndef  included_ALE_Mesh_hh
      5: #include <sieve/Mesh.hh>
      6: #endif

      8: namespace ALE {
      9:   namespace CartesianSieveDef {
     10:     static void nextCodeOdd(int a[], const int n, const int size) {
     11:       for(int d = size-2; d >= 0; --d) a[d+1] = a[d];
     12:       if (n == (n/2)*2) {
     13:         a[0] = 1;
     14:       } else {
     15:         a[0] = 0;
     16:       }
     17:     };
     18:     static void nextCodeEven(int a[], const int n, const int size) {
     19:       ALE::CartesianSieveDef::nextCodeOdd(a, n, size);
     20:       a[0] ^= 1;
     21:       for(int d = 1; d < size; ++d) a[d] ^= 0;
     22:     };
     23:     static void getGrayCode(const int n, const int size, int code[]) {
     24:       if (n == 1 << size) {
     25:         code[0] = -1;
     26:       } else if (n == 0) {
     27:         for(int d = 0; d < size; ++d) code[d] = 0;
     28:       } else if (n == 1) {
     29:         code[0] = 1;
     30:         for(int d = 1; d < size; ++d) code[d] = 0;
     31:       } else if (n == (n/2)*2) {
     32:         // Might have to shift pointer
     33:         ALE::CartesianSieveDef::getGrayCode(n/2, size, code);
     34:         ALE::CartesianSieveDef::nextCodeEven(code, n/2, size);
     35:       } else {
     36:         // Might have to shift pointer
     37:         ALE::CartesianSieveDef::getGrayCode((n-1)/2, size, code);
     38:         ALE::CartesianSieveDef::nextCodeOdd(code, (n-1)/2, size);
     39:       }
     40:     };
     41:     static void getIndices(const int dim, const int sizes[], const int& point, int indices[]) {
     42:       int p       = point;
     43:       int divisor = 1;

     45:       for(int d = 0; d < dim-1; ++d) {
     46:         divisor *= sizes[d];
     47:       }
     48:       //std::cout << "  Getting indices for point " << point << ":";
     49:       for(int d = dim-1; d >= 0; --d) {
     50:         indices[d] = p/divisor;
     51:         p -= divisor*indices[d];
     52:         if (d > 0) divisor /= sizes[d-1];
     53:       }
     54:       //for(int d = 0; d < dim; ++d) {
     55:       //  std::cout << " " << indices[d];
     56:       //}
     57:       //std::cout << std::endl;
     58:     };
     59:     static inline int getValue(const int dim, const int sizes[], const int indices[], const int code[], const bool forward = true) {
     60:       //std::cout << "Building value" << std::endl << "  code:";
     61:       //for(int d = 0; d < dim; ++d) std::cout << " " << code[d];
     62:       //std::cout << std::endl << "  indices:";
     63:       //for(int d = 0; d < dim; ++d) std::cout << " " << indices[d];
     64:       //std::cout << std::endl;
     65:       if (code[0] < 0) {
     66:         //std::cout << "  value " << -1 << std::endl;
     67:         return -1;
     68:       }
     69:       int value = indices[dim-1];
     70:       if (forward) {
     71:         if (code[dim-1]) value++;
     72:       } else {
     73:         if (!code[dim-1]) value--;
     74:       }
     75:       //std::cout << "  value " << value << std::endl;
     76:       for(int d = dim-2; d >= 0; --d) {
     77:         value = value*sizes[d] + indices[d];
     78:         if (forward) {
     79:           if (code[d]) value++;
     80:         } else {
     81:           if (!code[d]) value--;
     82:         }
     83:         //std::cout << "  value " << value << std::endl;
     84:       }
     85:       //std::cout << "  final value " << value << std::endl;
     86:       return value;
     87:     };

     90:     template <typename PointType_>
     91:     class PointSequence {
     92:     public:
     93:       typedef PointType_ point_type;
     94:       // We traverse the points in natural order
     95:       class iterator {
     96:       public:
     97:         // Standard iterator typedefs
     98:         typedef std::input_iterator_tag iterator_category;
     99:         typedef point_type              value_type;
    100:         typedef int                     difference_type;
    101:         typedef value_type*             pointer;
    102:         typedef value_type&             reference;
    103:       protected:
    104:         value_type _value;
    105:       public:
    106:         iterator(const value_type& start) : _value(start) {};
    107:         iterator(const iterator& iter) : _value(iter._value) {};
    108:         virtual ~iterator() {};
    109:       public:
    110:         virtual bool             operator==(const iterator& iter) const {return this->_value == iter._value;};
    111:         virtual bool             operator!=(const iterator& iter) const {return this->_value != iter._value;};
    112:         virtual const value_type operator*() const {return this->_value;};
    113:         virtual iterator&        operator++() {++this->_value; return *this;};
    114:         virtual iterator         operator++(int n) {iterator tmp(*this); ++this->_value; return tmp;};
    115:       };
    116:     protected:
    117:       point_type _start;
    118:       point_type _end;
    119:     public:
    120:       PointSequence(const point_type& start, const point_type& end) : _start(start), _end(end) {};
    121:       virtual ~PointSequence() {};
    122:     public:
    123:       virtual iterator begin() {return iterator(this->_start);};
    124:       virtual iterator end()   {return iterator(this->_end);};
    125:       virtual bool     empty() {return this->_start >= this->_end;};
    126:     };
    127:     // In one dimension, we have
    128:     //
    129:     //   v_k---c_k---v_{k+1}
    130:     //
    131:     // In two dimensions, we have
    132:     //
    133:     //            c_{k+m}             c_{k+m+1}
    134:     //
    135:     // v_{k+m+1}-----------v_{k+m+2}--
    136:     //    |                    |
    137:     //    |       c_k          |      c_{k+1}
    138:     //    |                    |
    139:     //   v_k---------------v_{k+1}----
    140:     //
    141:     // So for c_k the cone is (v_{k+(k/m)}, v_{k+(k/m)+1}, v_{k+m+(k/m)}, v_{k+m+(k/m)+1})
    142:     //
    143:     // In three dimensions, we have
    144:     //
    145:     //         v_{k+m+1+(m+1)*(n+1)} v_{k+m+1+(m+1)*(n+1)+1}
    146:     //
    147:     // v_{k+m+1}           v_{k+m+2}
    148:     //
    149:     //
    150:     //         v_{k+(m+1)*(n+1)}     v_{k+(m+1)*(n+1)+1}
    151:     //
    152:     //   v_k               v_{k+1}
    153:     //
    154:     // Suppose we break down the index k into a tuple (i,j), then 2d becomes
    155:     //
    156:     //   cone(c_{i,j}) = (v_{i,j}, v_{i+1,j}, v_{i,j+1}, v_{i+1,j+1})
    157:     //                 = (v_{k}, v_{k+1}, v_{k+m+1}, v_{k+m+2})
    158:     // Now for 2D
    159:     //   i = k%m
    160:     //   j = k/m
    161:     // and 3D
    162:     //   k = q/(m*n)
    163:     //   j = (q - m*n*k)/m
    164:     //   i = (q - m*n*k - m*j)
    165:     template <typename PointType_>
    166:     class ConeSequence {
    167:     public:
    168:       typedef PointType_ point_type;
    169:       // We traverse the points in natural order
    170:       class iterator {
    171:       public:
    172:         // Standard iterator typedefs
    173:         typedef std::input_iterator_tag iterator_category;
    174:         typedef PointType_              value_type;
    175:         typedef int                     difference_type;
    176:         typedef value_type*             pointer;
    177:         typedef value_type&             reference;
    178:       protected:
    179:         int               _dim;
    180:         const int        *_sizes;
    181:         const int        *_cellSizes;
    182:         int               _numCells;
    183:         int               _pos;
    184:         int              *_code;
    185:         value_type        _value;
    186:         value_type *_indices;
    187:       protected:
    188:         void init() {
    189:           this->_code = new int[this->_dim];
    190:           if (this->_pos == -2) {
    191:             this->_value = this->_indices[0];
    192:             ++this->_pos;
    193:           } else if (this->_pos == -1) {
    194:             this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_cellSizes, this->_indices, this->_code);
    195:           } else {
    196:             ALE::CartesianSieveDef::getGrayCode(this->_pos, this->_dim, this->_code);
    197:             this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code) + this->_numCells;
    198:           }
    199:         };
    200:       public:
    201:         iterator(const int dim, const int sizes[], const int cellSizes[], const int numCells, const value_type indices[], const int pos, const bool addSelf = false) : _dim(dim), _sizes(sizes), _cellSizes(cellSizes), _numCells(numCells), _pos(pos) {
    202:           this->_indices = new int[this->_dim];
    203:           for(int d = 0; d < this->_dim; ++d) this->_indices[d] = indices[d];
    204:           this->init();
    205:         };
    206:         iterator(const iterator& iter) : _dim(iter._dim), _sizes(iter._sizes), _cellSizes(iter._cellSizes), _numCells(iter._numCells), _pos(iter._pos) {
    207:           this->_indices = new int[this->_dim];
    208:           for(int d = 0; d < this->_dim; ++d) this->_indices[d] = iter._indices[d];
    209:           this->init();
    210:         };
    211:         virtual ~iterator() {
    212:           delete [] this->_code;
    213:           delete [] this->_indices;
    214:         };
    215:       public:
    216:         virtual bool              operator==(const iterator& iter) const {return this->_pos == iter._pos;};
    217:         virtual bool              operator!=(const iterator& iter) const {return this->_pos != iter._pos;};
    218:         virtual const value_type  operator*() const {return this->_value;};
    219:         virtual iterator&         operator++() {
    220:           ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
    221:           this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code) + this->_numCells;
    222:           return *this;
    223:         };
    224:         virtual iterator          operator++(int n) {
    225:           iterator tmp(*this);
    226:           ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
    227:           this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code) + this->_numCells;
    228:           return tmp;
    229:         };
    230:       };
    231:     protected:
    232:       int         _dim;
    233:       const int  *_sizes;
    234:       int        *_vertexSizes;
    235:       int         _numCells;
    236:       point_type  _cell;
    237:       int        *_indices;
    238:       bool        _addSelf;
    239:     public:
    240:       ConeSequence(const int dim, const int sizes[], const int numCells, const point_type& cell, const bool addSelf = false) : _dim(dim), _sizes(sizes), _numCells(numCells), _cell(cell), _addSelf(addSelf) {
    241:         this->_vertexSizes = new int[dim];
    242:         this->_indices     = new point_type[dim];
    243:         //this->getIndices(this->_cell, this->_indices);
    244:         for(int d = 0; d < this->_dim; ++d) this->_vertexSizes[d] = this->_sizes[d]+1;
    245:         ALE::CartesianSieveDef::getIndices(this->_dim, this->_sizes, this->_cell, this->_indices);
    246:       };
    247:       virtual ~ConeSequence() {
    248:         delete [] this->_vertexSizes;
    249:         delete [] this->_indices;
    250:       };
    251:     protected:
    252:       void getIndices(const point_type& cell, const int indices[]) {
    253:         point_type c       = cell;
    254:         int        divisor = 1;

    256:         for(int d = 0; d < this->_dim-1; ++d) {
    257:           divisor *= this->_sizes[d];
    258:         }
    259:         //std::cout << "  Got indices for cell " << cell << ":";
    260:         for(int d = this->_dim-1; d >= 0; --d) {
    261:           this->_indices[d] = c/divisor;
    262:           c -= divisor*this->_indices[d];
    263:           if (d > 0) divisor /= this->_sizes[d-1];
    264:         }
    265:         //for(int d = 0; d < this->_dim; ++d) {
    266:         //  std::cout << " " << this->_indices[d];
    267:         //}
    268:         //std::cout << std::endl;
    269:       };
    270:     public:
    271:       virtual iterator begin() {
    272:         int start = 0;

    274:         if (this->_addSelf) {
    275:           if (this->_cell >= this->_numCells) {
    276:             start = -2;
    277:             this->_indices[0] = this->_cell;
    278:           } else {
    279:             start = -1;
    280:           }
    281:         }
    282:         return iterator(this->_dim, this->_vertexSizes, this->_sizes, this->_numCells, this->_indices, start);
    283:       };
    284:       virtual iterator end() {
    285:         int end = 1 << this->_dim;

    287:         if (!this->_dim || (this->_cell >= this->_numCells)) {
    288:           end = 0;
    289:         }
    290:         return iterator(this->_dim, this->_vertexSizes, this->_sizes, this->_numCells, this->_indices, end);
    291:       };
    292:     };
    293:     template <typename PointType_>
    294:     class SupportSequence {
    295:     public:
    296:       typedef PointType_ point_type;
    297:       // We traverse the points in natural order
    298:       class iterator {
    299:       public:
    300:         // Standard iterator typedefs
    301:         typedef std::input_iterator_tag iterator_category;
    302:         typedef PointType_              value_type;
    303:         typedef int                     difference_type;
    304:         typedef value_type*             pointer;
    305:         typedef value_type&             reference;
    306:       protected:
    307:         int               _dim;
    308:         const int        *_sizes;
    309:         int               _pos;
    310:         int              *_code;
    311:         value_type        _value;
    312:         value_type *_indices;
    313:       protected:
    314:         bool validIndex(const int indices[], const int code[]) {
    315:           for(int d = 0; d < this->_dim; ++d) {
    316:             if ((code[d])  && (indices[d] >= this->_sizes[d])) return false;
    317:             if ((!code[d]) && (indices[d] < 1)) return false;
    318:           }
    319:           return true;
    320:         }
    321:         void init() {
    322:           this->_code = new int[this->_dim];
    323:           ALE::CartesianSieveDef::getGrayCode(this->_pos, this->_dim, this->_code);
    324:           while((this->_code[0] >= 0) && !this->validIndex(this->_indices, this->_code)) {
    325:             ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
    326:           }
    327:           this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code, false);
    328:         };
    329:       public:
    330:         iterator(const int dim, const int sizes[], const value_type indices[], const int pos) : _dim(dim), _sizes(sizes), _pos(pos) {
    331:           this->_indices = new int[this->_dim];
    332:           for(int d = 0; d < this->_dim; ++d) this->_indices[d] = indices[d];
    333:           this->init();
    334:         };
    335:         iterator(const iterator& iter) : _dim(iter._dim), _sizes(iter._sizes), _pos(iter._pos) {
    336:           this->_indices = new int[this->_dim];
    337:           for(int d = 0; d < this->_dim; ++d) this->_indices[d] = iter._indices[d];
    338:           this->init();
    339:         };
    340:         virtual ~iterator() {
    341:           delete [] this->_code;
    342:           delete [] this->_indices;
    343:         };
    344:       public:
    345:         virtual bool              operator==(const iterator& iter) const {return this->_pos == iter._pos;};
    346:         virtual bool              operator!=(const iterator& iter) const {return this->_pos != iter._pos;};
    347:         virtual const value_type  operator*() const {return this->_value;};
    348:         virtual iterator&         operator++() {
    349:           do {
    350:             ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
    351:           } while((this->_code[0] >= 0) && !this->validIndex(this->_indices, this->_code));
    352:           this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code, false);
    353:           return *this;
    354:         };
    355:         virtual iterator          operator++(int n) {
    356:           iterator tmp(*this);
    357:           do {
    358:             ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
    359:           } while((this->_code[0] >= 0) && !this->validIndex(this->_indices, this->_code));
    360:           this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code, false);
    361:           return tmp;
    362:         };
    363:       };
    364:     protected:
    365:       int         _dim;
    366:       const int  *_sizes;
    367:       int        *_vertexSizes;
    368:       int         _numCells;
    369:       point_type  _vertex;
    370:       int        *_indices;
    371:     public:
    372:       SupportSequence(const int dim, const int sizes[], const int numCells, const point_type& vertex) : _dim(dim), _sizes(sizes), _numCells(numCells), _vertex(vertex) {
    373:         this->_vertexSizes = new int[dim];
    374:         this->_indices     = new point_type[dim];
    375:         //this->getIndices(this->_vertex, this->_indices);
    376:         for(int d = 0; d < this->_dim; ++d) this->_vertexSizes[d] = this->_sizes[d]+1;
    377:         ALE::CartesianSieveDef::getIndices(this->_dim, this->_vertexSizes, this->_vertex - this->_numCells, this->_indices);
    378:       };
    379:       virtual ~SupportSequence() {
    380:         delete [] this->_vertexSizes;
    381:         delete [] this->_indices;
    382:       };
    383:     protected:
    384:       void getIndices(const point_type& vertex, const int indices[]) {
    385:         point_type v       = vertex - this->_numCells;
    386:         int        divisor = 1;

    388:         for(int d = 0; d < this->_dim-1; ++d) {
    389:           divisor *= this->_sizes[d]+1;
    390:         }
    391:         std::cout << "  Got indices for vertex " << vertex << ":";
    392:         for(int d = this->_dim-1; d >= 0; --d) {
    393:           this->_indices[d] = v/divisor;
    394:           v -= divisor*this->_indices[d];
    395:           if (d > 0) divisor /= this->_sizes[d-1]+1;
    396:         }
    397:         for(int d = 0; d < this->_dim; ++d) {
    398:           std::cout << " " << this->_indices[d];
    399:         }
    400:         std::cout << std::endl;
    401:       };
    402:     public:
    403:       virtual iterator begin() {return iterator(this->_dim, this->_sizes, this->_indices, 0);};
    404:       virtual iterator end()   {return iterator(this->_dim, this->_sizes, this->_indices, this->_dim ? 1 << this->_dim: 0);};
    405:     };
    406:   }
    407:   // We can do meets of two cells as empty if they do not intersect, and a k-dimensional face (with 2^k vertices)
    408:   //   if they have k indices in common
    409:   // We can do joins of two vertices. If they have have k indices in common, they span a d-k face, and thus have
    410:   //   2^{k} cells in the join.
    411:   // Note meet and join are both empty if any index differs by more than 1

    413:   // This is a purely Cartesian mesh
    414:   //   We will only represent cells and vertices, meaning the cap consists only of
    415:   //     vertices, and the base only of cells.
    416:   //   numCells:    m x n x ...
    417:   //   numVertices: m+1 x n+1 x ...
    418:   //   Points:
    419:   //     cells:       0     - numCells-1
    420:   //     vertices: numCells - numCells+numVertices-1
    421:   template <typename Point_>
    422:   class CartesianSieve : public ALE::ParallelObject {
    423:   public:
    424:     typedef Point_ point_type;
    425:     typedef CartesianSieveDef::PointSequence<point_type>   baseSequence;
    426:     typedef CartesianSieveDef::PointSequence<point_type>   capSequence;
    427:     typedef CartesianSieveDef::ConeSequence<point_type>    coneSequence;
    428:     typedef CartesianSieveDef::SupportSequence<point_type> supportSequence;
    429:     // Backward compatibility
    430:     typedef coneSequence coneArray;
    431:     typedef coneSequence coneSet;
    432:     typedef supportSequence supportArray;
    433:     typedef supportSequence supportSet;
    434:   protected:
    435:     int  _dim;
    436:     int *_sizes;
    437:     int  _numCells;
    438:     int  _numVertices;
    439:   public:
    440:     CartesianSieve(MPI_Comm comm, const int dim, const int numCells[], const int& debug = 0) : ParallelObject(comm, debug), _dim(dim) {
    441:       this->_sizes       = new int[dim];
    442:       this->_numCells    = 1;
    443:       this->_numVertices = 1;
    444:       for(int d = 0; d < dim; ++d) {
    445:         this->_sizes[d]     = numCells[d];
    446:         this->_numCells    *= numCells[d];
    447:         this->_numVertices *= numCells[d]+1;
    448:       }
    449:     };
    450:     virtual ~CartesianSieve() {
    451:       delete [] this->_sizes;
    452:     };
    453:   public:
    454:     int getDimension() const {return this->_dim;};
    455:     const int *getSizes() const {return this->_sizes;};
    456:     int getNumCells() const {return this->_numCells;};
    457:     int getNumVertices() const {return this->_numVertices;};
    458:   public:
    459:     // WARNING: Creating a lot of objects here
    460:     Obj<baseSequence> base() {
    461:       Obj<baseSequence> base = new baseSequence(0, this->_numCells);
    462:       return base;
    463:     };
    464:     Obj<baseSequence> leaves() {
    465:       return this->base();
    466:     };
    467:     // WARNING: Creating a lot of objects here
    468:     Obj<capSequence> cap() {
    469:        Obj<capSequence> cap = new capSequence(this->_numCells, this->_numCells+this->_numVertices);
    470:        return cap;
    471:     };
    472:     Obj<capSequence> roots() {
    473:        return this->cap();
    474:     };
    475:     // WARNING: Creating a lot of objects here
    476:     Obj<coneSequence> cone(const point_type& p) {
    477:       if ((p < 0) || (p >= this->_numCells)) {
    478:         Obj<coneSequence> cone = new coneSequence(0, this->_sizes, this->_numCells, 0);
    479:         return cone;
    480:       }
    481:       Obj<coneSequence> cone = new coneSequence(this->_dim, this->_sizes, this->_numCells, p);
    482:       return cone;
    483:     };
    484:     // WARNING: Creating a lot of objects here
    485:     Obj<supportSequence> support(const point_type& p) {
    486:       if ((p < this->_numCells) || (p >= this->_numCells+this->_numVertices)) {
    487:         Obj<supportSequence> support = new supportSequence(0, this->_sizes, this->_numCells, 0);
    488:         return support;
    489:       }
    490:       Obj<supportSequence> support = new supportSequence(this->_dim, this->_sizes, this->_numCells, p);
    491:       return support;
    492:     };
    493:     // WARNING: Creating a lot of objects here
    494:     Obj<coneSequence> closure(const point_type& p) {
    495:       Obj<coneSequence> cone = new coneSequence(this->_dim, this->_sizes, this->_numCells, p, true);
    496:       return cone;
    497:     };
    498:     // WARNING: Creating a lot of objects here
    499:     Obj<supportSequence> star(const point_type& p) {return this->support(p);};
    500:     // WARNING: Creating a lot of objects here
    501:     Obj<coneSequence> nCone(const point_type& p, const int n) {
    502:       if (n == 0) return this->cone(-1);
    503:       if (n != 1) throw ALE::Exception("Invalid height for nCone");
    504:       return this->cone(p);
    505:     };
    506:     // WARNING: Creating a lot of objects here
    507:     Obj<supportSequence> nSupport(const point_type& p, const int n) {
    508:       if (n == 0) return this->support(-1);
    509:       if (n != 1) throw ALE::Exception("Invalid depth for nSupport");
    510:       return this->support(p);
    511:     };
    512:   public:
    513:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
    514:         ostringstream txt;

    516:       if (comm == MPI_COMM_NULL) {
    517:         comm = this->comm();
    518:       }
    519:       if (name == "") {
    520:         PetscPrintf(comm, "viewing a CartesianSieve\n");
    521:       } else {
    522:         PetscPrintf(comm, "viewing CartesianSieve '%s'\n", name.c_str());
    523:       }
    524:       if(this->commRank() == 0) {
    525:         txt << "cap --> base:\n";
    526:       }
    527:       Obj<capSequence>  cap  = this->cap();
    528:       Obj<baseSequence> base = this->base();
    529:       if(cap->empty()) {
    530:         txt << "[" << this->commRank() << "]: empty" << std::endl;
    531:       }
    532:       for(typename capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
    533:         const Obj<supportSequence>&              supp    = this->support(*capi);
    534:         const typename supportSequence::iterator suppEnd = supp->end();

    536:         for(typename supportSequence::iterator suppi = supp->begin(); suppi != suppEnd; ++suppi) {
    537:           txt << "[" << this->commRank() << "]: " << *capi << "---->" << *suppi << std::endl;
    538:         }
    539:       }
    540:       PetscSynchronizedPrintf(this->comm(), txt.str().c_str());
    541:       PetscSynchronizedFlush(this->comm());
    542:       //
    543:       ostringstream txt1;
    544:       if(this->commRank() == 0) {
    545:         txt1 << "base --> cap:\n";
    546:       }
    547:       if(base->empty()) {
    548:         txt1 << "[" << this->commRank() << "]: empty" << std::endl;
    549:       }
    550:       for(typename baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
    551:         const Obj<coneSequence>&              cone    = this->cone(*basei);
    552:         const typename coneSequence::iterator coneEnd = cone->end();

    554:         for(typename coneSequence::iterator conei = cone->begin(); conei != coneEnd; ++conei) {
    555:           txt1 << "[" << this->commRank() << "]: " << *basei << "<----" << *conei << std::endl;
    556:         }
    557:       }
    558:       //
    559:       PetscSynchronizedPrintf(this->comm(), txt1.str().c_str());
    560:       PetscSynchronizedFlush(this->comm());
    561:       //
    562:       ostringstream txt2;
    563:       if(this->commRank() == 0) {
    564:         txt2 << "cap <point>:\n";
    565:       }
    566:       txt2 << "[" << this->commRank() << "]:  [";
    567:       for(typename capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
    568:         txt2 << " <" << *capi << ">";
    569:       }
    570:       txt2 << " ]" << std::endl;
    571:       //
    572:       PetscSynchronizedPrintf(this->comm(), txt2.str().c_str());
    573:       PetscSynchronizedFlush(this->comm());
    574:       //
    575:       ostringstream txt3;
    576:       if(this->commRank() == 0) {
    577:         txt3 << "base <point>:\n";
    578:       }
    579:       txt3 << "[" << this->commRank() << "]:  [";
    580:       for(typename baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
    581:         txt3 << " <" << *basei << ">";
    582:       }
    583:       txt3 << " ]" << std::endl;
    584:       //
    585:       PetscSynchronizedPrintf(this->comm(), txt3.str().c_str());
    586:       PetscSynchronizedFlush(this->comm());
    587:     };
    588:   };

    590:   // We do not just use Bundle, because we need to optimize labels
    591:   class CartesianBundle : public ALE::ParallelObject {
    592:   public:
    593:     typedef CartesianSieve<int>                                       sieve_type;
    594:     typedef sieve_type::point_type                                    point_type;
    595:     typedef Section<point_type, double>                               real_section_type;
    596:     typedef Section<point_type, int>                                  int_section_type;
    597:     typedef UniformSection<MinimalArrow<point_type, point_type>, int> arrow_section_type;
    598:     typedef std::map<std::string, Obj<arrow_section_type> >           arrow_sections_type;
    599:     typedef std::map<std::string, Obj<real_section_type> >            real_sections_type;
    600:     typedef std::map<std::string, Obj<int_section_type> >             int_sections_type;
    601:     typedef sieve_type::baseSequence                                  label_sequence;
    602:     typedef ALE::Sifter<int,point_type,point_type>                    send_overlap_type;
    603:     typedef ALE::Sifter<point_type,int,point_type>                    recv_overlap_type;
    604:   protected:
    605:     Obj<sieve_type>        _sieve;
    606:     int                    _maxHeight;
    607:     int                    _maxDepth;
    608:     arrow_sections_type    _arrowSections;
    609:     real_sections_type     _realSections;
    610:     int_sections_type      _intSections;
    611:     Obj<send_overlap_type> _sendOverlap;
    612:     Obj<recv_overlap_type> _recvOverlap;
    613:     Obj<send_overlap_type> _distSendOverlap;
    614:     Obj<recv_overlap_type> _distRecvOverlap;
    615:   public:
    616:     CartesianBundle(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1) {
    617:       this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
    618:       this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
    619:     };
    620:     virtual ~CartesianBundle() {};
    621:     public: // Accessors
    622:       void clear() {
    623:         this->_maxHeight = -1;
    624:         this->_maxDepth  = -1;
    625:       };
    626:     const Obj<sieve_type>& getSieve() const {return this->_sieve;};
    627:     void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
    628:     bool hasArrowSection(const std::string& name) const {
    629:       return this->_arrowSections.find(name) != this->_arrowSections.end();
    630:     };
    631:     const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
    632:       if (!this->hasArrowSection(name)) {
    633:         Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());

    635:         section->setName(name);
    636:         if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
    637:         this->_arrowSections[name] = section;
    638:       }
    639:       return this->_arrowSections[name];
    640:     };
    641:     void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
    642:       this->_arrowSections[name] = section;
    643:     };
    644:     Obj<std::set<std::string> > getArrowSections() const {
    645:       Obj<std::set<std::string> > names = std::set<std::string>();

    647:       for(arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
    648:         names->insert(s_iter->first);
    649:       }
    650:       return names;
    651:     };
    652:     bool hasRealSection(const std::string& name) const {
    653:       return this->_realSections.find(name) != this->_realSections.end();
    654:     };
    655:     const Obj<real_section_type>& getRealSection(const std::string& name) {
    656:       if (!this->hasRealSection(name)) {
    657:         Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());

    659:         section->setName(name);
    660:         if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
    661:         this->_realSections[name] = section;
    662:       }
    663:       return this->_realSections[name];
    664:     };
    665:     void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
    666:       this->_realSections[name] = section;
    667:     };
    668:     Obj<std::set<std::string> > getRealSections() const {
    669:       Obj<std::set<std::string> > names = std::set<std::string>();

    671:       for(real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
    672:         names->insert(s_iter->first);
    673:       }
    674:       return names;
    675:     };
    676:     bool hasIntSection(const std::string& name) const {
    677:       return this->_intSections.find(name) != this->_intSections.end();
    678:     };
    679:     const Obj<int_section_type>& getIntSection(const std::string& name) {
    680:       if (!this->hasIntSection(name)) {
    681:         Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());

    683:         section->setName(name);
    684:         if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
    685:         this->_intSections[name] = section;
    686:       }
    687:       return this->_intSections[name];
    688:     };
    689:     void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
    690:       this->_intSections[name] = section;
    691:     };
    692:     Obj<std::set<std::string> > getIntSections() const {
    693:       Obj<std::set<std::string> > names = std::set<std::string>();

    695:       for(int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
    696:         names->insert(s_iter->first);
    697:       }
    698:       return names;
    699:     };
    700:     const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
    701:     void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
    702:     const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
    703:     void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
    704:     const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
    705:     void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
    706:     const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
    707:     void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
    708:   public: // Stratification
    711:     virtual void stratify() {
    712:       ALE_LOG_EVENT_BEGIN;
    713:       this->_maxHeight = 1;
    714:       this->_maxDepth  = 1;
    715:       ALE_LOG_EVENT_END;
    716:     };
    717:     virtual int height() const {return this->_maxHeight;};
    718:     virtual int height(const point_type& point) {
    719:       const int numCells    = this->getSieve()->getNumCells();
    720:       const int numVertices = this->getSieve()->getNumVertices();
    721:       if ((point >= 0)        && (point < numCells))             return 0;
    722:       if ((point >= numCells) && (point < numCells+numVertices)) return 1;
    723:       return -1;
    724:     };
    725:     // WARNING: Creating a lot of objects here
    726:     virtual const Obj<label_sequence> heightStratum(int height) {
    727:       if (height == 0) return this->getSieve()->base();
    728:       if (height == 1) return this->getSieve()->cap();
    729:       throw ALE::Exception("Invalid height stratum");
    730:     };
    731:     virtual int depth() const {return this->_maxDepth;};
    732:     virtual int depth(const point_type& point) {
    733:       const int numCells    = this->getSieve()->getNumCells();
    734:       const int numVertices = this->getSieve()->getNumVertices();
    735:       if ((point >= 0)        && (point < numCells))             return 1;
    736:       if ((point >= numCells) && (point < numCells+numVertices)) return 0;
    737:       return -1;
    738:     };
    739:     // WARNING: Creating a lot of objects here
    740:     virtual const Obj<label_sequence> depthStratum(int depth) {
    741:       if (depth == 0) return this->getSieve()->cap();
    742:       if (depth == 1) return this->getSieve()->base();
    743:       throw ALE::Exception("Invalid depth stratum");
    744:     };
    745:     virtual const Obj<label_sequence> getLabelStratum(const std::string& name, int value) {
    746:       if (name == "height") {
    747:         return this->heightStratum(value);
    748:       } else if (name == "depth") {
    749:         return this->depthStratum(value);
    750:       }
    751:       throw ALE::Exception("Invalid label name");
    752:     }
    753:   public: // Viewers
    754:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
    755:       if (comm == MPI_COMM_NULL) {
    756:         comm = this->comm();
    757:       }
    758:       if (name == "") {
    759:         PetscPrintf(comm, "viewing a CartesianBundle\n");
    760:       } else {
    761:         PetscPrintf(comm, "viewing CartesianBundle '%s'\n", name.c_str());
    762:       }
    763:       PetscPrintf(comm, "  maximum height %d maximum depth %d\n", this->height(), this->depth());
    764:     };
    765:   public:
    766:     void constructOverlap() {};
    767:   };

    769:   class CartesianMesh : public CartesianBundle {
    770:   public:
    771:     typedef int                                point_type;
    772:     typedef CartesianBundle::sieve_type        sieve_type;
    773:     typedef CartesianBundle::label_sequence    label_sequence;
    774:     typedef CartesianBundle::send_overlap_type send_overlap_type;
    775:     typedef CartesianBundle::recv_overlap_type recv_overlap_type;
    776:     typedef sieve_type::coneSequence           coneSequence;
    777:     typedef sieve_type::supportSequence        supportSequence;
    778:   protected:
    779:     int                    _dim;
    780:     // Discretization
    781:     Obj<Discretization>    _discretization;
    782:     Obj<BoundaryCondition> _boundaryCondition;
    783:   public:
    784:     CartesianMesh(MPI_Comm comm, int dim, int debug = 0) : CartesianBundle(comm, debug), _dim(dim) {
    785:       this->_discretization    = new Discretization(comm, debug);
    786:       this->_boundaryCondition = new BoundaryCondition(comm, debug);
    787:     };
    788:     virtual ~CartesianMesh() {};
    789:   public: // Accessors
    790:     int getDimension() {return this->_dim;};
    791:     const Obj<Discretization>&    getDiscretization() {return this->_discretization;};
    792:     void setDiscretization(const Obj<Discretization>& discretization) {this->_discretization = discretization;};
    793:     const Obj<BoundaryCondition>& getBoundaryCondition() {return this->_boundaryCondition;};
    794:     void setBoundaryCondition(const Obj<BoundaryCondition>& boundaryCondition) {this->_boundaryCondition = boundaryCondition;};
    795:   public: // Discretization
    796: #if 0
    797:     void markBoundaryCells(const std::string& name) {
    798:       const Obj<label_type>&     label    = this->getLabel(name);
    799:       const Obj<label_sequence>& boundary = this->getLabelStratum(name, 1);
    800:       const Obj<sieve_type>&     sieve    = this->getSieve();

    802:       for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
    803:         if (this->height(*e_iter) == 1) {
    804:           const point_type cell = *sieve->support(*e_iter)->begin();

    806:           this->setValue(label, cell, 2);
    807:         }
    808:       }
    809:     };
    810: #endif
    811:     void setupField(const Obj<real_section_type>& s, const bool postponeGhosts = false) {
    812:       const std::string& name = this->_boundaryCondition->getLabelName();

    814:       for(int d = 0; d <= this->depth(); ++d) {
    815:         s->setFiberDimension(this->depthStratum(d), this->_discretization->getNumDof(d));
    816:       }
    817:       if (!name.empty()) {
    818:         const Obj<label_sequence>& boundary = this->getLabelStratum(name, 1);

    820:         for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
    821:           s->setFiberDimension(*e_iter, -this->_discretization->getNumDof(this->depth(*e_iter)));
    822:         }
    823:       }
    824:       s->allocatePoint();
    825:       if (!name.empty()) {
    826: #if 0
    827:         const Obj<label_sequence>&     boundaryCells = this->getLabelStratum(name, 2);
    828:         const Obj<real_section_type>&  coordinates   = this->getRealSection("coordinates");
    829:         real_section_type::value_type *values        = new real_section_type::value_type[this->sizeWithBC(s, *boundaryCells->begin())];
    830:         double                        *v0            = new double[this->getDimension()];
    831:         double                        *J             = new double[this->getDimension()*this->getDimension()];
    832:         double                         detJ;

    834:         for(label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
    835:           const Obj<coneArray>      closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
    836:           const coneArray::iterator end     = closure->end();
    837:           int                       v       = 0;

    839:           this->computeElementGeometry(coordinates, *c_iter, v0, J, NULL, detJ);
    840:           for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
    841:             const int cDim = s->getConstraintDimension(*cl_iter);

    843:             if (cDim) {
    844:               for(int d = 0; d < cDim; ++d, ++v) {
    845:                 values[v] = this->_boundaryCondition->integrateDual(v0, J, v);
    846:               }
    847:             } else {
    848:               const int dim = s->getFiberDimension(*cl_iter);

    850:               for(int d = 0; d < dim; ++d, ++v) {
    851:                 values[v] = 0.0;
    852:               }
    853:             }
    854:           }
    855:           this->updateAll(s, *c_iter, values);
    856:         }
    857:         delete [] values;
    858:         delete [] v0;
    859:         delete [] J;
    860: #endif
    861:       }
    862:     };
    863:   public: // Size traversal
    864:     template<typename Section_>
    865:     int size(const Obj<Section_>& section, const point_type& p) {
    866:       const typename Section_::chart_type& chart = section->getChart();
    867:       const Obj<coneSequence>              cone  = this->getSieve()->cone(p);
    868:       typename coneSequence::iterator      end   = cone->end();
    869:       int                                  size  = 0;

    871:       if (chart.count(p)) {
    872:         size += section->getConstrainedFiberDimension(p);
    873:       }
    874:       for(typename coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
    875:         if (chart.count(*c_iter)) {
    876:           size += section->getConstrainedFiberDimension(*c_iter);
    877:         }
    878:       }
    879:       return size;
    880:     }
    881:     template<typename Section_>
    882:     int sizeWithBC(const Obj<Section_>& section, const point_type& p) {
    883:       const typename Section_::chart_type& chart = section->getChart();
    884:       const Obj<coneSequence>              cone  = this->getSieve()->cone(p);
    885:       typename coneSequence::iterator      end   = cone->end();
    886:       int                                  size  = 0;

    888:       if (chart.count(p)) {
    889:         size += section->getFiberDimension(p);
    890:       }
    891:       for(typename coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
    892:         if (chart.count(*c_iter)) {
    893:           size += section->getFiberDimension(*c_iter);
    894:         }
    895:       }
    896:       return size;
    897:     }
    898:   public: // Allocation
    899:     template<typename Section_>
    900:     void allocate(const Obj<Section_>& section) {
    901:       section->allocatePoint();
    902:     }
    903:   public: // Retrieval traversal
    904:     // Return the values for the closure of this point
    905:     //   use a smart pointer?
    906:     template<typename Section_>
    907:     const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p) {
    908:       const int                       size   = this->sizeWithBC(section, p);
    909:       typename Section_::value_type  *values = section->getRawArray(size);
    910:       int                             j      = -1;

    912:       const int& dim = section->getFiberDimension(p);
    913:       const typename Section_::value_type *array = section->restrictPoint(p);

    915:       for(int i = 0; i < dim; ++i) {
    916:         values[++j] = array[i];
    917:       }
    918:       const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
    919:       typename sieve_type::coneSequence::iterator   end  = cone->end();

    921:       for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
    922:         const int& dim = section->getFiberDimension(*p_iter);

    924:         array = section->restrictPoint(*p_iter);
    925:         for(int i = 0; i < dim; ++i) {
    926:           values[++j] = array[i];
    927:         }
    928:       }
    929:       if (j != size-1) {
    930:         ostringstream txt;

    932:         txt << "Invalid restrictClosure to point " << p << std::endl;
    933:         txt << "  j " << j << " should be " << (size-1) << std::endl;
    934:         std::cout << txt.str();
    935:         throw ALE::Exception(txt.str().c_str());
    936:       }
    937:       return values;
    938:     }
    939:     template<typename Section_>
    940:     void update(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
    941:       int j = 0;

    943:       section->updatePoint(p, &v[j]);
    944:       j += section->getFiberDimension(p);
    945:       const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

    947:       for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
    948:         section->updatePoint(*p_iter, &v[j]);
    949:         j += section->getFiberDimension(*p_iter);
    950:       }
    951:     }
    952:     template<typename Section_>
    953:     void updateAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
    954:       int j = 0;

    956:       section->updateAddPoint(p, &v[j]);
    957:       j += section->getFiberDimension(p);
    958:       const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

    960:       for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
    961:         section->updateAddPoint(*p_iter, &v[j]);
    962:         j += section->getFiberDimension(*p_iter);
    963:       }
    964:     }
    965:     template<typename Section_>
    966:     void updateBC(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
    967:       int j = 0;

    969:       section->updatePointBC(p, &v[j]);
    970:       j += section->getFiberDimension(p);
    971:       const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

    973:       for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
    974:         section->updatePointBC(*p_iter, &v[j]);
    975:         j += section->getFiberDimension(*p_iter);
    976:       }
    977:     }
    978:     template<typename Section_>
    979:     void updateAll(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
    980:       int j = 0;

    982:       section->updatePointBC(p, &v[j]);
    983:       j += section->getFiberDimension(p);
    984:       const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

    986:       for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
    987:         section->updatePointAll(*p_iter, &v[j]);
    988:         j += section->getFiberDimension(*p_iter);
    989:       }
    990:     }
    991:   public: // Mesh geometry
    992:     void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
    993:       const double *coords = this->restrictClosure(coordinates, e);
    994:       const int     dim    = this->_dim;
    995:       const int     last   = 1 << dim;

    997:       if (v0) {
    998:         for(int d = 0; d < dim; d++) {
    999:           v0[d] = coords[d];
    1000:         }
    1001:       }
    1002:       for(int d = 0; d < dim; ++d) {
    1003:         if (J) {
    1004:           J[d] = 0.5*(coords[last*dim+d] - v0[d]);
    1005:         }
    1006:         if (invJ) {
    1007:           invJ[d] = 1.0/J[d];
    1008:         }
    1009:       }
    1010:     }
    1011:   public: // Viewers
    1012:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
    1013:       if (comm == MPI_COMM_NULL) {
    1014:         comm = this->comm();
    1015:       }
    1016:       if (name == "") {
    1017:         PetscPrintf(comm, "viewing a CartesianMesh\n");
    1018:       } else {
    1019:         PetscPrintf(comm, "viewing CartesianMesh '%s'\n", name.c_str());
    1020:       }
    1021:       this->getSieve()->view("mesh sieve");
    1022:     }
    1023:   };

    1025:   class CartesianMeshBuilder {
    1026:   public:
    1027:     static Obj<CartesianMesh> createCartesianMesh(const MPI_Comm comm, const int dim, const int numCells[], const int partitions[], const int debug = 0) {

    1029:       // Steal PETSc code that calculates partitions
    1030:       //   We could conceivably allow multiple patches per partition
    1031:       int size, rank;
    1032:       int totalPartitions = 1;

    1035:       MPI_Comm_size(comm, &size);CHKERRXX(ierr);
    1036:       MPI_Comm_rank(comm, &rank);CHKERRXX(ierr);
    1037:       for(int d = 0; d < dim; ++d) totalPartitions *= partitions[d];
    1038:       if (size != totalPartitions) throw ALE::Exception("Invalid partitioning");
    1039:       // Determine local sizes
    1040:       int *numLocalCells       = new int[dim];
    1041:       int *numNeighborVertices = new int[dim];
    1042:       int *numLocalVertices    = new int[dim];

    1044:       for(int d = 0; d < dim; ++d) {
    1045:         numLocalCells[d]    = numCells[d]/partitions[d] + (rank < numCells[d]%partitions[d]);
    1046:         numLocalVertices[d] = numLocalCells[d]+1;
    1047:       }
    1048:       // Create mesh
    1049:       const Obj<CartesianMesh>             mesh  = new CartesianMesh(comm, dim, debug);
    1050:       const Obj<CartesianMesh::sieve_type> sieve = new CartesianMesh::sieve_type(comm, dim, numLocalCells, debug);

    1052:       mesh->setSieve(sieve);
    1053:       mesh->stratify();
    1054:       delete [] numLocalCells;
    1055:       // Create overlaps
    1056:       //   We overlap anyone within 1 index of us for the entire boundary
    1057:       const Obj<CartesianMesh::send_overlap_type>& sendOverlap = mesh->getSendOverlap();
    1058:       const Obj<CartesianMesh::recv_overlap_type>& recvOverlap = mesh->getRecvOverlap();
    1059:       int *indices  = new int[dim];
    1060:       int *code     = new int[dim];
    1061:       int *zeroCode = new int[dim];
    1062:       int *vIndices = new int[dim];
    1063:       int numNeighborCells;

    1065:       ALE::CartesianSieveDef::getIndices(dim, partitions, rank, indices);
    1066:       for(int e = 0; e < dim; ++e) zeroCode[e] = 0;
    1067:       for(int d = 0; d < dim; ++d) {
    1068:         if (indices[d] > 0) {
    1069:           for(int e = 0; e < dim; ++e) code[e] = 1;
    1070:           code[d] = 0;
    1071:           int neighborRank = ALE::CartesianSieveDef::getValue(dim, partitions, indices, code, false);
    1072:           numNeighborCells = 1;
    1073:           for(int e = 0; e < dim; ++e) {
    1074:             numNeighborCells      *= (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]));
    1075:             numNeighborVertices[e] = (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]))+1;
    1076:           }

    1078:           // Add the whole d-1 face on the left edge of dimension d
    1079:           for(int v = sieve->getNumCells(); v < sieve->getNumCells()+sieve->getNumVertices(); ++v) {
    1080:             ALE::CartesianSieveDef::getIndices(dim, numLocalVertices, v - sieve->getNumCells(), vIndices);
    1081:             if (vIndices[d] == 0) {
    1082:               vIndices[d]   = numNeighborVertices[d]-1;
    1083:               int neighborV = ALE::CartesianSieveDef::getValue(dim, numNeighborVertices, vIndices, zeroCode) + numNeighborCells;

    1085:               sendOverlap->addCone(v, neighborRank, neighborV);
    1086:               recvOverlap->addCone(neighborRank, v, neighborV);
    1087:             }
    1088:           }
    1089:         }
    1090:         if (indices[d] < partitions[d]-1) {
    1091:           for(int e = 0; e < dim; ++e) code[e] = 0;
    1092:           code[d] = 1;
    1093:           int neighborRank = ALE::CartesianSieveDef::getValue(dim, partitions, indices, code);
    1094:           numNeighborCells = 1;
    1095:           for(int e = 0; e < dim; ++e) {
    1096:             numNeighborCells      *= (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]));
    1097:             numNeighborVertices[e] = (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]))+1;
    1098:           }

    1100:           // Add the whole d-1 face on the right edge of dimension d
    1101:           for(int v = sieve->getNumCells(); v < sieve->getNumCells()+sieve->getNumVertices(); ++v) {
    1102:             ALE::CartesianSieveDef::getIndices(dim, numLocalVertices, v - sieve->getNumCells(), vIndices);
    1103:             if (vIndices[d] == numLocalVertices[d]-1) {
    1104:               vIndices[d]   = 0;
    1105:               int neighborV = ALE::CartesianSieveDef::getValue(dim, numNeighborVertices, vIndices, zeroCode) + numNeighborCells;

    1107:               sendOverlap->addCone(v, neighborRank, neighborV);
    1108:               recvOverlap->addCone(neighborRank, v, neighborV);
    1109:             }
    1110:           }
    1111:         }
    1112:       }
    1113:       if (debug) {
    1114:         sendOverlap->view("Send Overlap");
    1115:         recvOverlap->view("Receive Overlap");
    1116:       }
    1117:       delete [] numLocalVertices;
    1118:       delete [] numNeighborVertices;
    1119:       delete [] indices;
    1120:       delete [] code;
    1121:       delete [] zeroCode;
    1122:       delete [] vIndices;
    1123:       // Create coordinates
    1124:       return mesh;
    1125:     };
    1126:   };
    1127: }

    1129: #endif