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

      1: #ifndef included_ALE_Selection_hh
      2: #define included_ALE_Selection_hh

      4: #ifndef  included_ALE_SieveAlgorithms_hh
      5: #include <sieve/SieveAlgorithms.hh>
      6: #endif

      8: #ifndef  included_ALE_SieveBuilder_hh
      9: #include <sieve/SieveBuilder.hh>
     10: #endif

     12: #ifndef  included_ALE_Mesh_hh
     13: #include <sieve/Mesh.hh>
     14: #endif

     16: namespace ALE {
     17:   template<typename Mesh_>
     18:   class Selection {
     19:   public:
     20:     typedef ALE::Mesh<PetscInt,PetscScalar>      FlexMesh;
     21:     typedef Mesh_                                mesh_type;
     22:     typedef typename mesh_type::sieve_type       sieve_type;
     23:     typedef typename mesh_type::point_type       point_type;
     24:     typedef typename mesh_type::int_section_type int_section_type;
     25:     typedef std::set<point_type>                 PointSet;
     26:     typedef std::vector<point_type>              PointArray;
     27:     typedef std::pair<typename sieve_type::point_type, int> oPoint_type;
     28:     typedef std::vector<oPoint_type>                        oPointArray;
     29:     typedef typename ALE::SieveAlg<mesh_type>    sieveAlg;
     30:   protected:
     31:     template<typename Sieve, typename FaceType>
     32:     class FaceVisitor {
     33:     public:
     34:       typedef typename Sieve::point_type point_type;
     35:       typedef typename Sieve::arrow_type arrow_type;
     36:     protected:
     37:       FaceType& face;
     38:       PointArray&     origVertices;
     39:       PointArray&     faceVertices;
     40:       int            *indices;
     41:       const int       debug;
     42:       int             oppositeVertex;
     43:       int             v;
     44:     public:
     45:       FaceVisitor(FaceType& f, PointArray& oV, PointArray& fV, int *i, const int debug) : face(f), origVertices(oV), faceVertices(fV), indices(i), debug(debug), oppositeVertex(-1), v(0) {};
     46:       void visitPoint(const point_type& point) {
     47:         if (face->find(point) != face->end()) {
     48:           if (debug) std::cout << "    vertex " << point << std::endl;
     49:           indices[origVertices.size()] = v;
     50:           origVertices.insert(origVertices.end(), point);
     51:         } else {
     52:           if (debug) std::cout << "    vertex " << point << std::endl;
     53:           oppositeVertex = v;
     54:         }
     55:         ++v;
     56:       };
     57:       void visitArrow(const arrow_type&) {};
     58:     public:
     59:       int getOppositeVertex() {return this->oppositeVertex;};
     60:     };
     61:   public:
     62:     template<typename Processor>
     63:     static void subsets(const PointArray& v, const int size, Processor& processor, Obj<PointArray> *out = NULL, const int min = 0) {
     64:       if (size == 0) {
     65:         processor(*out);
     66:         return;
     67:       }
     68:       if (min == 0) {
     69:         out  = new Obj<PointArray>();
     70:         *out = new PointArray();
     71:       }
     72:       for(int i = min; i < (int) v.size(); ++i) {
     73:         (*out)->push_back(v[i]);
     74:         subsets(v, size-1, processor, out, i+1);
     75:         (*out)->pop_back();
     76:       }
     77:       if (min == 0) {delete out;}
     78:     };
     79:     template<typename Mesh>
     80:     static int numFaceVertices(const Obj<Mesh>& mesh) {
     81:       return numFaceVertices(mesh, mesh->getNumCellCorners());
     82:     };
     83:     template<typename Mesh>
     84:     static int numFaceVertices(const Obj<Mesh>& mesh, const unsigned int numCorners) {
     85:       //unsigned int numCorners = mesh->getNumCellCorners(cell, depth);
     86:       const    int cellDim          = mesh->getDimension();
     87:       unsigned int _numFaceVertices = 0;

     89:       switch (cellDim) {
     90:       case 0 :
     91:         _numFaceVertices = 0;
     92:         break;
     93:       case 1 :
     94:         _numFaceVertices = 1;
     95:         break;
     96:       case 2:
     97:         switch (numCorners) {
     98:         case 3 : // triangle
     99:           _numFaceVertices = 2; // Edge has 2 vertices
    100:           break;
    101:         case 4 : // quadrilateral
    102:           _numFaceVertices = 2; // Edge has 2 vertices
    103:           break;
    104:         case 6 : // quadratic triangle
    105:           _numFaceVertices = 3; // Edge has 3 vertices
    106:           break;
    107:         case 9 : // quadratic quadrilateral
    108:           _numFaceVertices = 3; // Edge has 3 vertices
    109:           break;
    110:         default :
    111:           throw ALE::Exception("Invalid number of face corners");
    112:         }
    113:         break;
    114:       case 3:
    115:         switch (numCorners) {
    116:         case 4 : // tetradehdron
    117:           _numFaceVertices = 3; // Face has 3 vertices
    118:           break;
    119:         case 8 : // hexahedron
    120:           _numFaceVertices = 4; // Face has 4 vertices
    121:           break;
    122:         case 10 : // quadratic tetrahedron
    123:           _numFaceVertices = 6; // Face has 6 vertices
    124:           break;
    125:         case 27 : // quadratic hexahedron
    126:           _numFaceVertices = 9; // Face has 9 vertices
    127:           break;
    128:         default :
    129:           throw ALE::Exception("Invalid number of face corners");
    130:         }
    131:         break;
    132:       default:
    133:         throw ALE::Exception("Invalid cell dimension");
    134:       }
    135:       return _numFaceVertices;
    136:     };
    137:     // We need this method because we do not use interpolated sieves
    138:     //   - Without interpolation, we cannot say what vertex collections are
    139:     //     faces, and how they are oriented
    140:     //   - Now we read off the list of face vertices IN THE ORDER IN WHICH
    141:     //     THEY APPEAR IN THE CELL
    142:     //   - This leads to simple algorithms for simplices and quads to check
    143:     //     orientation since these sets are always valid faces
    144:     //   - This is not true with hexes, so we just sort and check explicit cases
    145:     //   - This means for hexes that we have to alter the vertex container as well
    146:     template<typename Mesh>
    147:     static bool faceOrientation(const point_type& cell, const Obj<Mesh>& mesh, const int numCorners,
    148:                                 const int indices[], const int oppositeVertex, PointArray *origVertices, PointArray *faceVertices) {
    149:       const int cellDim   = mesh->getDimension();
    150:       const int debug     = mesh->debug();
    151:       bool      posOrient = false;

    153:       if (debug) std::cout << "cellDim: " << cellDim << ", numCorners: " << numCorners << std::endl;

    155:       if (cellDim == numCorners-1) {
    156:         // Simplices
    157:         posOrient = !(oppositeVertex%2);
    158:       } else if (cellDim == 1 && numCorners == 3) {
    159:         posOrient = true;
    160:       } else if (cellDim == 2 && numCorners == 4) {
    161:         // Quads
    162:         if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
    163:           posOrient = true;
    164:         } else if ((indices[0] == 3) && (indices[1] == 0)) {
    165:           posOrient = true;
    166:         } else {
    167:           if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
    168:             posOrient = false;
    169:           } else {
    170:             throw ALE::Exception("Invalid quad crossedge");
    171:           }
    172:         }
    173:       } else if (cellDim == 2 && numCorners == 6) {
    174:         // Quadratic triangle (I hate this)
    175:         // Edges are determined by the first 2 vertices (corners of edges)
    176:         const int faceSizeTri = 3;
    177:         int  sortedIndices[3];
    178:         bool found = false;
    179:         int faceVerticesTriSorted[9] = {
    180:           0, 3,  4, // bottom
    181:           1, 4,  5, // right
    182:           2, 3,  5, // left
    183:         };
    184:         int faceVerticesTri[9] = {
    185:           0, 3,  4, // bottom
    186:           1, 4,  5, // right
    187:           2, 5,  3, // left
    188:         };

    190:         for(int i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
    191:         std::sort(sortedIndices, sortedIndices+faceSizeTri);
    192:         for (int iFace=0; iFace < 4; ++iFace) {
    193:           const int ii = iFace*faceSizeTri;
    194:           if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
    195:               (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
    196:             if (debug) {
    197:               if (iFace == 0) std::cout << "Bottom edge" << std::endl;
    198:               else if (iFace == 1) std::cout << "Right edge" << std::endl;
    199:               else if (iFace == 2) std::cout << "Left edge" << std::endl;
    200:             }  // if
    201:             for (int fVertex=0; fVertex < faceSizeTri; ++fVertex)
    202:               for (int cVertex=0; cVertex < faceSizeTri; ++cVertex)
    203:                 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
    204:                   faceVertices->push_back((*origVertices)[cVertex]);
    205:                   break;
    206:                 } // if
    207:             found = true;
    208:             break;
    209:           } // if
    210:         } // for
    211:         if (!found) {throw ALE::Exception("Invalid tri crossface");}
    212:         return true;
    213:       } else if (cellDim == 2 && numCorners == 9) {
    214:         // Quadratic quad (I hate this)
    215:         // Edges are determined by the first 2 vertices (corners of edges)
    216:         const int faceSizeQuad = 3;
    217:         int  sortedIndices[3];
    218:         bool found = false;
    219:         int faceVerticesQuadSorted[12] = {
    220:           0, 1,  4, // bottom
    221:           1, 2,  5, // right
    222:           2, 3,  6, // top
    223:           0, 3,  7, // left
    224:         };
    225:         int faceVerticesQuad[12] = {
    226:           0, 1,  4, // bottom
    227:           1, 2,  5, // right
    228:           2, 3,  6, // top
    229:           3, 0,  7, // left
    230:         };

    232:         for(int i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
    233:         std::sort(sortedIndices, sortedIndices+faceSizeQuad);
    234:         for (int iFace=0; iFace < 4; ++iFace) {
    235:           const int ii = iFace*faceSizeQuad;
    236:           if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
    237:               (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
    238:             if (debug) {
    239:               if (iFace == 0) std::cout << "Bottom edge" << std::endl;
    240:               else if (iFace == 1) std::cout << "Right edge" << std::endl;
    241:               else if (iFace == 2) std::cout << "Top edge" << std::endl;
    242:               else if (iFace == 3) std::cout << "Left edge" << std::endl;
    243:             }  // if
    244:             for (int fVertex=0; fVertex < faceSizeQuad; ++fVertex)
    245:               for (int cVertex=0; cVertex < faceSizeQuad; ++cVertex)
    246:                 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
    247:                   faceVertices->push_back((*origVertices)[cVertex]);
    248:                   break;
    249:                 } // if
    250:             found = true;
    251:             break;
    252:           } // if
    253:         } // for
    254:         if (!found) {throw ALE::Exception("Invalid quad crossface");}
    255:         return true;
    256:       } else if (cellDim == 3 && numCorners == 8) {
    257:         // Hexes
    258:         //   A hex is two oriented quads with the normal of the first
    259:         //   pointing up at the second.
    260:         //
    261:         //     7---6
    262:         //    /|  /|
    263:         //   4---5 |
    264:         //   | 3-|-2
    265:         //   |/  |/
    266:         //   0---1
    267:         //
    268:         // Faces are determined by the first 4 vertices (corners of faces)
    269:         const int faceSizeHex = 4;
    270:         int  sortedIndices[4];
    271:         bool found = false;
    272:         int faceVerticesHexSorted[24] = {
    273:           0, 1, 2, 3,  // bottom
    274:           4, 5, 6, 7,  // top
    275:           0, 1, 4, 5,  // front
    276:           1, 2, 5, 6,  // right
    277:           2, 3, 6, 7,  // back
    278:           0, 3, 4, 7,  // left
    279:         };
    280:         int faceVerticesHex[24] = {
    281:           3, 2, 1, 0,  // bottom
    282:           4, 5, 6, 7,  // top
    283:           0, 1, 5, 4,  // front
    284:           1, 2, 6, 5,  // right
    285:           2, 3, 7, 6,  // back
    286:           3, 0, 4, 7,  // left
    287:         };

    289:         for(int i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
    290:         std::sort(sortedIndices, sortedIndices+faceSizeHex);
    291:         for (int iFace=0; iFace < 6; ++iFace) {
    292:           const int ii = iFace*faceSizeHex;
    293:           if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
    294:               (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
    295:               (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
    296:               (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
    297:             if (debug) {
    298:               if (iFace == 0) std::cout << "Bottom quad" << std::endl;
    299:               else if (iFace == 1) std::cout << "Top quad" << std::endl;
    300:               else if (iFace == 2) std::cout << "Front quad" << std::endl;
    301:               else if (iFace == 3) std::cout << "Right quad" << std::endl;
    302:               else if (iFace == 4) std::cout << "Back quad" << std::endl;
    303:               else if (iFace == 5) std::cout << "Left quad" << std::endl;
    304:             }  // if
    305:             for (int fVertex=0; fVertex < faceSizeHex; ++fVertex)
    306:               for (int cVertex=0; cVertex < faceSizeHex; ++cVertex)
    307:                 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
    308:                   faceVertices->push_back((*origVertices)[cVertex]);
    309:                   break;
    310:                 } // if
    311:             found = true;
    312:             break;
    313:           } // if
    314:         } // for
    315:         if (!found) {throw ALE::Exception("Invalid hex crossface");}
    316:         return true;
    317:       } else if (cellDim == 3 && numCorners == 10) {
    318:         // Quadratic tet
    319:         // Faces are determined by the first 3 vertices (corners of faces)
    320:         const int faceSizeTet = 6;
    321:         int  sortedIndices[6];
    322:         bool found = false;
    323:         int faceVerticesTetSorted[24] = {
    324:           0, 1, 2,  6, 7, 8, // bottom
    325:           0, 3, 4,  6, 7, 9,  // front
    326:           1, 4, 5,  7, 8, 9,  // right
    327:           2, 3, 5,  6, 8, 9,  // left
    328:         };
    329:         int faceVerticesTet[24] = {
    330:           0, 1, 2,  6, 7, 8, // bottom
    331:           0, 4, 3,  6, 7, 9,  // front
    332:           1, 5, 4,  7, 8, 9,  // right
    333:           2, 3, 5,  8, 6, 9,  // left
    334:         };

    336:         for(int i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
    337:         std::sort(sortedIndices, sortedIndices+faceSizeTet);
    338:         for (int iFace=0; iFace < 6; ++iFace) {
    339:           const int ii = iFace*faceSizeTet;
    340:           if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
    341:               (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
    342:               (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
    343:               (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
    344:             if (debug) {
    345:               if (iFace == 0) std::cout << "Bottom tri" << std::endl;
    346:               else if (iFace == 1) std::cout << "Front tri" << std::endl;
    347:               else if (iFace == 2) std::cout << "Right tri" << std::endl;
    348:               else if (iFace == 3) std::cout << "Left tri" << std::endl;
    349:             }  // if
    350:             for (int fVertex=0; fVertex < faceSizeTet; ++fVertex)
    351:               for (int cVertex=0; cVertex < faceSizeTet; ++cVertex)
    352:                 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
    353:                   faceVertices->push_back((*origVertices)[cVertex]);
    354:                   break;
    355:                 } // if
    356:             found = true;
    357:             break;
    358:           } // if
    359:         } // for
    360:         if (!found) {throw ALE::Exception("Invalid tet crossface");}
    361:         return true;
    362:       } else if (cellDim == 3 && numCorners == 27) {
    363:         // Quadratic hexes (I hate this)
    364:         //   A hex is two oriented quads with the normal of the first
    365:         //   pointing up at the second.
    366:         //
    367:         //     7---6
    368:         //    /|  /|
    369:         //   4---5 |
    370:         //   | 3-|-2
    371:         //   |/  |/
    372:         //   0---1
    373:         //
    374:         // Faces are determined by the first 4 vertices (corners of faces)
    375:         const int faceSizeQuadHex = 9;
    376:         int  sortedIndices[9];
    377:         bool found = false;
    378:         int faceVerticesQuadHexSorted[54] = {
    379:           0, 1, 2, 3,  8, 9, 10, 11,  24, // bottom
    380:           4, 5, 6, 7,  12, 13, 14, 15,  25, // top
    381:           0, 1, 4, 5,  8, 12, 16, 17,  22, // front
    382:           1, 2, 5, 6,  9, 13, 17, 18,  21, // right
    383:           2, 3, 6, 7,  10, 14, 18, 19,  23, // back
    384:           0, 3, 4, 7,  11, 15, 16, 19,  20, // left
    385:         };
    386:         int faceVerticesQuadHex[54] = {
    387:           3, 2, 1, 0,  10, 9, 8, 11,  24, // bottom
    388:           4, 5, 6, 7,  12, 13, 14, 15,  25, // top
    389:           0, 1, 5, 4,  8, 17, 12, 16,  22, // front
    390:           1, 2, 6, 5,  9, 18, 13, 17,  21, // right
    391:           2, 3, 7, 6,  10, 19, 14, 18,  23, // back
    392:           3, 0, 4, 7,  11, 16, 15, 19,  20 // left
    393:         };

    395:         for (int i=0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
    396:         std::sort(sortedIndices, sortedIndices+faceSizeQuadHex);
    397:         for (int iFace=0; iFace < 6; ++iFace) {
    398:           const int ii = iFace*faceSizeQuadHex;
    399:           if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
    400:               (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
    401:               (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
    402:               (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
    403:             if (debug) {
    404:               if (iFace == 0) std::cout << "Bottom quad" << std::endl;
    405:               else if (iFace == 1) std::cout << "Top quad" << std::endl;
    406:               else if (iFace == 2) std::cout << "Front quad" << std::endl;
    407:               else if (iFace == 3) std::cout << "Right quad" << std::endl;
    408:               else if (iFace == 4) std::cout << "Back quad" << std::endl;
    409:               else if (iFace == 5) std::cout << "Left quad" << std::endl;
    410:             }  // if
    411:             for (int fVertex=0; fVertex < faceSizeQuadHex; ++fVertex)
    412:               for (int cVertex=0; cVertex < faceSizeQuadHex; ++cVertex)
    413:                 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
    414:                   faceVertices->push_back((*origVertices)[cVertex]);
    415:                   break;
    416:                 } // if
    417:             found = true;
    418:             break;
    419:           } // if
    420:         } // for
    421:         if (!found) {throw ALE::Exception("Invalid hex crossface");}
    422:         return true;
    423:       } else {
    424:         throw ALE::Exception("Unknown cell type for faceOrientation().");
    425:       }
    426:       if (!posOrient) {
    427:         if (debug) std::cout << "  Reversing initial face orientation" << std::endl;
    428:         faceVertices->insert(faceVertices->end(), (*origVertices).rbegin(), (*origVertices).rend());
    429:       } else {
    430:         if (debug) std::cout << "  Keeping initial face orientation" << std::endl;
    431:         faceVertices->insert(faceVertices->end(), (*origVertices).begin(), (*origVertices).end());
    432:       }
    433:       return posOrient;
    434:     };
    435:     // Given a cell and a face, as a set of vertices,
    436:     //   return the oriented face, as a set of vertices, in faceVertices
    437:     // The orientation is such that the face normal points out of the cell
    438:     template<typename Mesh, typename FaceType>
    439:     static bool getOrientedFace(const Obj<Mesh>& mesh, const point_type& cell, FaceType face,
    440:                                 const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
    441:     {
    442:       FaceVisitor<typename Mesh::sieve_type,FaceType> v(face, *origVertices, *faceVertices, indices, mesh->debug());

    444:       origVertices->clear();
    445:       faceVertices->clear();
    446:       mesh->getSieve()->cone(cell, v);
    447:       return faceOrientation(cell, mesh, numCorners, indices, v.getOppositeVertex(), origVertices, faceVertices);
    448:     };
    449:     template<typename FaceType>
    450:     static bool getOrientedFace(const Obj<FlexMesh>& mesh, const point_type& cell, FaceType face,
    451:                                 const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
    452:     {
    453:       const Obj<typename FlexMesh::sieve_type::traits::coneSequence>& cone = mesh->getSieve()->cone(cell);
    454:       const int debug = mesh->debug();
    455:       int       v     = 0;
    456:       int       oppositeVertex;

    458:       origVertices->clear();
    459:       faceVertices->clear();
    460:       for(typename FlexMesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter, ++v) {
    461:         if (face->find(*v_iter) != face->end()) {
    462:           if (debug) std::cout << "    vertex " << *v_iter << std::endl;
    463:           indices[origVertices->size()] = v;
    464:           origVertices->insert(origVertices->end(), *v_iter);
    465:         } else {
    466:           if (debug) std::cout << "    vertex " << *v_iter << std::endl;
    467:           oppositeVertex = v;
    468:         }
    469:       }
    470:       return faceOrientation(cell, mesh, numCorners, indices, oppositeVertex, origVertices, faceVertices);
    471:     };
    472:     template<typename Sieve>
    473:     static void insertFace(const Obj<mesh_type>& mesh, const Obj<Sieve>& subSieve, const Obj<PointSet>& face, point_type& f,
    474:                            const point_type& cell, const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
    475:     {
    476:       const Obj<typename Sieve::supportSet> preFace = subSieve->nJoin1(face);
    477:       const int                             debug   = subSieve->debug();

    479:       if (preFace->size() > 1) {
    480:         throw ALE::Exception("Invalid fault sieve: Multiple faces from vertex set");
    481:       } else if (preFace->size() == 1) {
    482:         // Add the other cell neighbor for this face
    483:         subSieve->addArrow(*preFace->begin(), cell);
    484:       } else if (preFace->size() == 0) {
    485:         if (debug) std::cout << "  Orienting face " << f << std::endl;
    486:         try {
    487:           getOrientedFace(mesh, cell, face, numCorners, indices, origVertices, faceVertices);
    488:           if (debug) std::cout << "  Adding face " << f << std::endl;
    489:           int color = 0;
    490:           for(typename PointArray::const_iterator f_iter = faceVertices->begin(); f_iter != faceVertices->end(); ++f_iter) {
    491:             if (debug) std::cout << "    vertex " << *f_iter << std::endl;
    492:             subSieve->addArrow(*f_iter, f, color++);
    493:           }
    494:           subSieve->addArrow(f, cell);
    495:           f++;
    496:         } catch (ALE::Exception e) {
    497:           if (debug) std::cout << "  Did not add invalid face " << f << std::endl;
    498:         }
    499:       }
    500:     };
    501:   public:
    502:     class FaceInserter {
    503: #if 0
    504:     public:
    505:       typedef Mesh_                                mesh_type;
    506:       typedef typename mesh_type::sieve_type       sieve_type;
    507:       typedef typename mesh_type::point_type       point_type;
    508:       typedef std::set<point_type>                 PointSet;
    509:       typedef std::vector<point_type>              PointArray;
    510: #endif
    511:     protected:
    512:       const Obj<mesh_type>  mesh;
    513:       const Obj<sieve_type> sieve;
    514:       const Obj<sieve_type> subSieve;
    515:       point_type&           f;
    516:       const point_type      cell;
    517:       const int             numCorners;
    518:       int                  *indices;
    519:       PointArray           *origVertices;
    520:       PointArray           *faceVertices;
    521:       PointSet             *subCells;
    522:       const int             debug;
    523:     public:
    524:       FaceInserter(const Obj<mesh_type>& mesh, const Obj<sieve_type>& sieve, const Obj<sieve_type>& subSieve, point_type& f, const point_type& cell, const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices, PointSet *subCells) : mesh(mesh), sieve(sieve), subSieve(subSieve), f(f), cell(cell), numCorners(numCorners), indices(indices), origVertices(origVertices), faceVertices(faceVertices), subCells(subCells), debug(mesh->debug()) {};
    525:       virtual ~FaceInserter() {};
    526:     public:
    527:       void operator()(const Obj<PointArray>& face) {
    528:         const Obj<typename sieve_type::supportSet> sievePreFace = sieve->nJoin1(face);

    530:         if (sievePreFace->size() == 1) {
    531:           if (debug) std::cout << "  Contains a boundary face on the submesh" << std::endl;
    532:           PointSet faceSet(face->begin(), face->end());
    533:           ALE::Selection<mesh_type>::insertFace(mesh, subSieve, faceSet, f, cell, numCorners, indices, origVertices, faceVertices);
    534:           subCells->insert(cell);
    535:         }
    536:       };
    537:     };
    538:     template<typename Sieve>
    539:     class FaceInserterV {
    540:     protected:
    541:       const Obj<mesh_type>&  mesh;
    542:       const Obj<sieve_type>& sieve;
    543:       const Obj<Sieve>&      subSieve;
    544:       point_type&            f;
    545:       const point_type       cell;
    546:       const int              numCorners;
    547:       int                   *indices;
    548:       PointArray            *origVertices;
    549:       PointArray            *faceVertices;
    550:       PointSet              *subCells;
    551:       const int              debug;
    552:     public:
    553:       FaceInserterV(const Obj<mesh_type>& mesh, const Obj<sieve_type>& sieve, const Obj<Sieve>& subSieve, point_type& f, const point_type& cell, const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices, PointSet *subCells) : mesh(mesh), sieve(sieve), subSieve(subSieve), f(f), cell(cell), numCorners(numCorners), indices(indices), origVertices(origVertices), faceVertices(faceVertices), subCells(subCells), debug(mesh->debug()) {};
    554:       virtual ~FaceInserterV() {};
    555:     public:
    556:       void operator()(const Obj<PointArray>& face) {
    557:         ISieveVisitor::PointRetriever<sieve_type> jV(sieve->getMaxSupportSize());

    559:         sieve->join(*face, jV);
    560:         if (jV.getSize() == 1) {
    561:           if (debug) std::cout << "  Contains a boundary face on the submesh" << std::endl;
    562:           PointSet faceSet(face->begin(), face->end());
    563:           ALE::Selection<mesh_type>::insertFace(mesh, subSieve, faceSet, f, cell, numCorners, indices, origVertices, faceVertices);
    564:           subCells->insert(cell);
    565:         }
    566:       };
    567:     };
    568:   protected:
    569:     static int binomial(const int i, const int j) {
    570:       assert(j <= i);
    571:       assert(i < 34);
    572:       if (j == 0) {
    573:         return 1;
    574:       } else if (j == i) {
    575:         return 1;
    576:       } else {
    577:         return binomial(i-1, j) + binomial(i-1, j-1);
    578:       }
    579:     };
    580:   public:
    581:     // This takes in a section and creates a submesh from the vertices in the section chart
    582:     //   This is a hyperplane of one dimension lower than the mesh
    583:     static Obj<mesh_type> submesh_uninterpolated(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1, const bool boundaryFaces = true) {
    584:       // A challenge here is to coordinate the extra numbering of new faces
    585:       //   In serial, it is enough to number after the last point:
    586:       //     Use sieve->base()->size() + sieve->cap()->size(), or determine the greatest point
    587:       //   In parallel, there are two steps:
    588:       //     1) Use the serial result, and reduce either with add (for size) or max (for greatest point)
    589:       //     2) Determine how many faces will be created on each process
    590:       //        This will be bounded by C(numCorners, faceSize)*submeshCells
    591:       //        Thus it looks like we should first accumulate submeshCells, and then create faces
    592:       typedef typename mesh_type::label_type        label_type;
    593:       typedef typename int_section_type::chart_type chart_type;
    594:       const int                  dim        = (dimension > 0) ? dimension : mesh->getDimension()-1;
    595:       const Obj<sieve_type>&     sieve      = mesh->getSieve();
    596:       Obj<mesh_type>             submesh    = new mesh_type(mesh->comm(), dim, mesh->debug());
    597:       Obj<sieve_type>            subSieve   = new sieve_type(mesh->comm(), mesh->debug());
    598:       const bool                 censor     = mesh->hasLabel("censored depth");
    599:       const Obj<label_type>&     depthLabel = censor ? mesh->getLabel("censored depth") : mesh->getLabel("depth");
    600:       const int                  depth      = mesh->depth();
    601:       const int                  height     = mesh->height();
    602:       const chart_type&          chart      = label->getChart();
    603:       const int                  numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), depth)->size();
    604:       const int                  faceSize   = numFaceVertices(mesh);
    605:       Obj<PointSet>              face       = new PointSet();
    606:       int                        f          = sieve->base()->size() + sieve->cap()->size();
    607:       const int                  debug      = mesh->debug();
    608:       int                       *indices    = new int[faceSize];
    609:       PointArray                 origVertices, faceVertices;
    610:       PointSet                   submeshVertices, submeshCells;

    613:       const typename chart_type::iterator chartEnd = chart.end();
    614:       for(typename chart_type::iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
    615:         if (label->getFiberDimension(*c_iter)) submeshVertices.insert(*c_iter);
    616:       }
    617:       const typename PointSet::const_iterator svBegin = submeshVertices.begin();
    618:       const typename PointSet::const_iterator svEnd   = submeshVertices.end();

    620:       for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
    621:         const Obj<typename sieveAlg::supportArray>&     cells  = sieveAlg::nSupport(mesh, *sv_iter, depth);
    622:         const typename sieveAlg::supportArray::iterator cBegin = cells->begin();
    623:         const typename sieveAlg::supportArray::iterator cEnd   = cells->end();

    625:         if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
    626:         for(typename sieveAlg::supportArray::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
    627:           if (debug) std::cout << "  Checking cell " << *c_iter << std::endl;
    628:           if (submeshCells.find(*c_iter) != submeshCells.end()) continue;
    629:           if (censor && (!mesh->getValue(depthLabel, *c_iter))) continue;
    630:           const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *c_iter, height);
    631:           const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
    632:           const typename sieveAlg::coneArray::iterator vEnd   = cone->end();

    634:           face->clear();
    635:           for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
    636:             if (submeshVertices.find(*v_iter) != svEnd) {
    637:               if (debug) std::cout << "    contains submesh vertex " << *v_iter << std::endl;
    638:               face->insert(face->end(), *v_iter);
    639:             }
    640:           }
    641:           if ((int) face->size() > faceSize) {
    642:             if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
    643:             if (debug) std::cout << "  Has all boundary faces on the submesh" << std::endl;
    644:             submeshCells.insert(*c_iter);
    645:           }
    646:           if ((int) face->size() == faceSize) {
    647:             if (debug) std::cout << "  Contains a face on the submesh" << std::endl;
    648:             submeshCells.insert(*c_iter);
    649:           }
    650:         }
    651:       }
    652:       if (mesh->commSize() > 1) {
    653:         int localF     = f;
    654:         int localFaces = binomial(numCorners, faceSize)*submeshCells.size();
    655:         int maxFaces;

    657:         MPI_Allreduce(&localF, &f, 1, MPI_INT, MPI_SUM, mesh->comm());
    658:         //     2) Determine how many faces will be created on each process
    659:         //        This will be bounded by faceSize*submeshCells
    660:         //        Thus it looks like we should first accumulate submeshCells, and then create faces
    661:         MPI_Scan(&localFaces, &maxFaces, 1, MPI_INT, MPI_SUM, mesh->comm());
    662:         f += maxFaces - localFaces;
    663:       }
    664:       for(typename PointSet::const_iterator c_iter = submeshCells.begin(); c_iter != submeshCells.end(); ++c_iter) {
    665:         if (debug) std::cout << "  Processing submesh cell " << *c_iter << std::endl;
    666:         const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *c_iter, height);
    667:         const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
    668:         const typename sieveAlg::coneArray::iterator vEnd   = cone->end();

    670:         face->clear();
    671:         for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
    672:           if (submeshVertices.find(*v_iter) != svEnd) {
    673:             if (debug) std::cout << "    contains submesh vertex " << *v_iter << std::endl;
    674:             face->insert(face->end(), *v_iter);
    675:           }
    676:         }
    677:         if ((int) face->size() > faceSize) {
    678:           // Here we allow a set of vertices to lie completely on a boundary cell (like a corner tetrahedron)
    679:           //   We have to take all the faces, and discard those in the interior
    680:           FaceInserter inserter(mesh, sieve, subSieve, f, *c_iter, numCorners, indices, &origVertices, &faceVertices, &submeshCells);
    681:           PointArray   faceVec(face->begin(), face->end());

    683:           subsets(faceVec, faceSize, inserter);
    684:         }
    685:         if ((int) face->size() == faceSize) {
    686:           insertFace(mesh, subSieve, face, f, *c_iter, numCorners, indices, &origVertices, &faceVertices);
    687:         }
    688:       }
    689:       delete [] indices;
    690:       submesh->setSieve(subSieve);
    691:       submesh->stratify();
    692:       if (debug) submesh->view("Submesh");
    693:       return submesh;
    694:     };
    695:     // This takes in a section and creates a submesh from the vertices in the section chart
    696:     //   This is a hyperplane of one dimension lower than the mesh
    697:     static Obj<mesh_type> submesh_interpolated(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1, const bool boundaryFaces = true) {
    698:       const int debug  = mesh->debug();
    699:       const int depth  = mesh->depth();
    700:       const int height = mesh->height();
    701:       const typename int_section_type::chart_type&                chart        = label->getChart();
    702:       const typename int_section_type::chart_type::const_iterator chartEnd     = chart.end();
    703:       const Obj<PointSet>                                         submeshFaces = new PointSet();
    704:       PointSet submeshVertices;

    706:       for(typename int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
    707:         //assert(!mesh->depth(*c_iter));
    708:         submeshVertices.insert(*c_iter);
    709:       }
    710:       const typename PointSet::const_iterator svBegin = submeshVertices.begin();
    711:       const typename PointSet::const_iterator svEnd   = submeshVertices.end();

    713:       for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
    714:         const Obj<typename sieveAlg::supportArray>& faces = sieveAlg::nSupport(mesh, *sv_iter, depth-1);
    715:         const typename sieveAlg::supportArray::iterator fBegin = faces->begin();
    716:         const typename sieveAlg::supportArray::iterator fEnd   = faces->end();

    718:         if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
    719:         for(typename sieveAlg::supportArray::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
    720:           if (debug) std::cout << "  Checking face " << *f_iter << std::endl;
    721:           if (submeshFaces->find(*f_iter) != submeshFaces->end()) continue;
    722:           const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *f_iter, height-1);
    723:           const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
    724:           const typename sieveAlg::coneArray::iterator vEnd   = cone->end();
    725:           bool                                         found  = true;

    727:           for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
    728:             if (submeshVertices.find(*v_iter) != svEnd) {
    729:               if (debug) std::cout << "    contains submesh vertex " << *v_iter << std::endl;
    730:             } else {
    731:               found = false;
    732:             }
    733:           }
    734:           if (found) {
    735:             if (boundaryFaces) {throw ALE::Exception("Not finished: should check that it is a boundary face");}
    736:             if (debug) std::cout << "  Is a face on the submesh" << std::endl;
    737:             submeshFaces->insert(*f_iter);
    738:           }
    739:         }
    740:       }
    741:       return submesh(mesh, submeshFaces, mesh->getDimension()-1);
    742:     };
    743:     template<typename output_mesh_type>
    744:     static Obj<output_mesh_type> submeshV_uninterpolated(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1, const bool boundaryFaces = true) {
    745:       typedef typename mesh_type::label_type        label_type;
    746:       typedef typename int_section_type::chart_type chart_type;
    747:       const int                           dim        = (dimension > 0) ? dimension : mesh->getDimension()-1;
    748:       const Obj<sieve_type>&              sieve      = mesh->getSieve();
    749:       Obj<FlexMesh>                       submesh    = new FlexMesh(mesh->comm(), dim, mesh->debug());
    750:       Obj<typename FlexMesh::sieve_type>  subSieve   = new typename FlexMesh::sieve_type(mesh->comm(), mesh->debug());
    751:       const bool                          censor     = mesh->hasLabel("censored depth");
    752:       const Obj<label_type>&              depthLabel = censor ? mesh->getLabel("censored depth") : mesh->getLabel("depth");
    753:       const chart_type&                   chart      = label->getChart();
    754:       const int                           numCorners = mesh->getNumCellCorners();
    755:       const int                           faceSize   = numFaceVertices(mesh);
    756:       Obj<PointSet>                       face       = new PointSet();
    757:       int                                 f          = sieve->getBaseSize() + sieve->getCapSize();
    758:       const int                           debug      = mesh->debug();
    759:       int                                *indices    = new int[faceSize];
    760:       PointArray                          origVertices, faceVertices;
    761:       PointSet                            submeshVertices, submeshCells;

    763:       const typename chart_type::const_iterator chartEnd = chart.end();
    764:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
    765:         if (label->getFiberDimension(*c_iter)) submeshVertices.insert(*c_iter);
    766:       }
    767:       const typename PointSet::const_iterator svBegin = submeshVertices.begin();
    768:       const typename PointSet::const_iterator svEnd   = submeshVertices.end();
    769:       typename ISieveVisitor::PointRetriever<sieve_type> sV(sieve->getMaxSupportSize());
    770:       typename ISieveVisitor::PointRetriever<sieve_type> cV(sieve->getMaxConeSize());

    772:       for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
    773:         sieve->support(*sv_iter, sV);
    774:         const int         numCells = sV.getSize();
    775:         const point_type *cells    = sV.getPoints();

    777:         if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
    778:         for(int c = 0; c < numCells; ++c) {
    779:           if (debug) std::cout << "  Checking cell " << cells[c] << std::endl;
    780:           if (submeshCells.find(cells[c]) != submeshCells.end()) continue;
    781:           if (censor && (!mesh->getValue(depthLabel, cells[c]))) continue;
    782:           sieve->cone(cells[c], cV);
    783:           const int         numVertices = cV.getSize();
    784:           const point_type *vertices    = cV.getPoints();

    786:           face->clear();
    787:           for(int v = 0; v < numVertices; ++v) {
    788:             if (submeshVertices.find(vertices[v]) != svEnd) {
    789:               if (debug) std::cout << "    contains submesh vertex " << vertices[v] << std::endl;
    790:               face->insert(face->end(), vertices[v]);
    791:             }
    792:           }
    793:           if ((int) face->size() > faceSize) {
    794:             if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
    795:             if (debug) std::cout << "  Has all boundary faces on the submesh" << std::endl;
    796:             submeshCells.insert(cells[c]);
    797:           }
    798:           if ((int) face->size() == faceSize) {
    799:             if (debug) std::cout << "  Contains a face on the submesh" << std::endl;
    800:             submeshCells.insert(cells[c]);
    801:           }
    802:           cV.clear();
    803:         }
    804:         sV.clear();
    805:       }
    806:       if (mesh->commSize() > 1) {
    807:         int localF     = f;
    808:         int localFaces = binomial(numCorners, faceSize)*submeshCells.size();
    809:         int maxFaces;

    811:         MPI_Allreduce(&localF, &f, 1, MPI_INT, MPI_SUM, mesh->comm());
    812:         //     2) Determine how many faces will be created on each process
    813:         //        This will be bounded by faceSize*submeshCells
    814:         //        Thus it looks like we should first accumulate submeshCells, and then create faces
    815:         MPI_Scan(&localFaces, &maxFaces, 1, MPI_INT, MPI_SUM, mesh->comm());
    816:         f += maxFaces - localFaces;
    817:       }
    818:       for(typename PointSet::const_iterator c_iter = submeshCells.begin(); c_iter != submeshCells.end(); ++c_iter) {
    819:         if (debug) std::cout << "  Processing submesh cell " << *c_iter << std::endl;
    820:         sieve->cone(*c_iter, cV);
    821:         const int         numVertices = cV.getSize();
    822:         const point_type *vertices    = cV.getPoints();

    824:         face->clear();
    825:         for(int v = 0; v < numVertices; ++v) {
    826:           if (submeshVertices.find(vertices[v]) != svEnd) {
    827:             if (debug) std::cout << "    contains submesh vertex " << vertices[v] << std::endl;
    828:             face->insert(face->end(), vertices[v]);
    829:           }
    830:         }
    831:         if ((int) face->size() > faceSize) {
    832:           if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
    833:           // Here we allow a set of vertices to lie completely on a boundary cell (like a corner tetrahedron)
    834:           //   We have to take all the faces, and discard those in the interior
    835:           FaceInserterV<FlexMesh::sieve_type> inserter(mesh, sieve, subSieve, f, *c_iter, numCorners, indices, &origVertices, &faceVertices, &submeshCells);
    836:           PointArray                          faceVec(face->begin(), face->end());

    838:           subsets(faceVec, faceSize, inserter);
    839:         }
    840:         if ((int) face->size() == faceSize) {
    841:           insertFace(mesh, subSieve, face, f, *c_iter, numCorners, indices, &origVertices, &faceVertices);
    842:         }
    843:         cV.clear();
    844:       }
    845:       delete [] indices;
    846:       submesh->setSieve(subSieve);
    847:       submesh->stratify();
    848:       if (debug) submesh->view("Submesh");

    850:       Obj<output_mesh_type> isubmesh = new output_mesh_type(submesh->comm(), submesh->getDimension(), submesh->debug());
    851:       Obj<typename output_mesh_type::sieve_type> isieve = new typename output_mesh_type::sieve_type(submesh->comm(), submesh->debug());
    852:       std::map<typename output_mesh_type::point_type,typename output_mesh_type::point_type> renumbering;
    853:       isubmesh->setSieve(isieve);
    854:       ALE::ISieveConverter::convertMesh(*submesh, *isubmesh, renumbering, false);
    855:       return isubmesh;
    856:     };
    857:   public:
    858:     // This takes in a section and creates a submesh from the vertices in the section chart
    859:     //   This is a hyperplane of one dimension lower than the mesh
    860:     static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
    861:       const int dim   = mesh->getDimension();
    862:       const int depth = mesh->depth();

    864:       if (dim == depth) {
    865:         return submesh_interpolated(mesh, label, dimension, false);
    866:       } else if (depth == 1) {
    867:         return submesh_uninterpolated(mesh, label, dimension);
    868:       }
    869:       throw ALE::Exception("Cannot handle partially interpolated meshes");
    870:     };
    871:     template<typename output_mesh_type>
    872:     static Obj<output_mesh_type> submeshV(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
    873:       const int dim   = mesh->getDimension();
    874:       const int depth = mesh->depth();

    876: #if 0
    877:       if (dim == depth) {
    878:         //return submesh_interpolated(mesh, label, dimension, false);
    879:         throw ALE::Exception("Cannot handle interpolated meshes");
    880:       } else if (depth == 1) {
    881:         return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
    882:       }
    883: #else
    884:       if (depth == 1) {
    885:         return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
    886:       } else if (dim == depth) {
    887:         //return submesh_interpolated(mesh, label, dimension, false);
    888:         throw ALE::Exception("Cannot handle interpolated meshes");
    889:       }
    890: #endif
    891:       throw ALE::Exception("Cannot handle partially interpolated meshes");
    892:     };
    893:     // This creates a submesh consisting of the union of the closures of the given points
    894:     //   This mesh is the same dimension as in the input mesh
    895:     template<typename Points>
    896:     static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<Points>& points, const int dim = -1) {
    897:       const Obj<sieve_type>& sieve     = mesh->getSieve();
    898:       Obj<mesh_type>         newMesh   = new mesh_type(mesh->comm(), dim >= 0 ? dim : mesh->getDimension(), mesh->debug());
    899:       Obj<sieve_type>        newSieve  = new sieve_type(mesh->comm(), mesh->debug());
    900:       Obj<PointSet>          newPoints = new PointSet();
    901:       Obj<PointSet>          modPoints = new PointSet();
    902:       Obj<PointSet>          tmpPoints;

    904:       newMesh->setSieve(newSieve);
    905:       for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
    906:         newPoints->insert(*p_iter);
    907:         do {
    908:           modPoints->clear();
    909:           for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
    910:             const Obj<typename sieve_type::traits::coneSequence>&     cone = sieve->cone(*np_iter);
    911:             const typename sieve_type::traits::coneSequence::iterator end  = cone->end();
    912:             int c = 0;

    914:             for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter, ++c) {
    915:               newSieve->addArrow(*c_iter, *np_iter, c);
    916:             }
    917:             modPoints->insert(cone->begin(), cone->end());
    918:           }
    919:           tmpPoints = newPoints;
    920:           newPoints = modPoints;
    921:           modPoints = tmpPoints;
    922:         } while(newPoints->size());
    923:         newPoints->insert(*p_iter);
    924:         do {
    925:           modPoints->clear();
    926:           for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
    927:             const Obj<typename sieve_type::traits::supportSequence>&     support = sieve->support(*np_iter);
    928:             const typename sieve_type::traits::supportSequence::iterator end     = support->end();
    929:             int s = 0;

    931:             for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != end; ++s_iter, ++s) {
    932:               newSieve->addArrow(*np_iter, *s_iter, s);
    933:             }
    934:             modPoints->insert(support->begin(), support->end());
    935:           }
    936:           tmpPoints = newPoints;
    937:           newPoints = modPoints;
    938:           modPoints = tmpPoints;
    939:         } while(newPoints->size());
    940:       }
    941:       newMesh->stratify();
    942:       newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
    943:       newMesh->setArrowSection("orientation", mesh->getArrowSection("orientation"));
    944:       return newMesh;
    945:     };
    946:   protected:
    947:     static Obj<mesh_type> boundary_uninterpolated(const Obj<mesh_type>& mesh) {
    948:       throw ALE::Exception("Not yet implemented");
    949:       const Obj<typename mesh_type::label_sequence>&     cells    = mesh->heightStratum(0);
    950:       const Obj<sieve_type>&                             sieve    = mesh->getSieve();
    951:       const typename mesh_type::label_sequence::iterator cBegin   = cells->begin();
    952:       const typename mesh_type::label_sequence::iterator cEnd     = cells->end();
    953:       const int                                          faceSize = numFaceVertices(mesh);

    955:       for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
    956:         const Obj<typename sieve_type::traits::coneSequence>&     vertices = sieve->cone(*c_iter);
    957:         const typename sieve_type::traits::coneSequence::iterator vBegin   = vertices->begin();
    958:         const typename sieve_type::traits::coneSequence::iterator vEnd     = vertices->end();
    959:         //PointArray cell(vertices->begin(), vertices->end());

    961:         // For each vertex, gather
    962:         for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
    963:           const Obj<typename sieve_type::traits::supportSequence>&     neighbors = sieve->support(*v_iter);
    964:           const typename sieve_type::traits::supportSequence::iterator nBegin    = neighbors->begin();
    965:           const typename sieve_type::traits::supportSequence::iterator nEnd      = neighbors->end();

    967:           for(typename sieve_type::traits::supportSequence::iterator n_iter = nBegin; n_iter != nEnd; ++n_iter) {
    968:             const Obj<typename sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, 1);

    970:             if (preFace->size() == faceSize) {
    971:             }
    972:           }
    973:         }
    974:         // For each face
    975:         // - determine if its legal

    977:         // - determine if its part of a neighboring cell
    978:         // - if not, its a boundary face
    979:         //subsets(cell, faceSize, inserter);
    980:       }
    981:     };
    982:     static void addClosure(const Obj<sieve_type>& sieveA, const Obj<sieve_type>& sieveB, const point_type& p, const int depth = 1) {
    983:       Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
    984:       Obj<typename sieve_type::coneSet> next    = new typename sieve_type::coneSet();
    985:       Obj<typename sieve_type::coneSet> tmp;

    987:       current->insert(p);
    988:       while(current->size()) {
    989:         for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
    990:           const Obj<typename sieve_type::traits::coneSequence>&     cone  = sieveA->cone(*p_iter);
    991:           const typename sieve_type::traits::coneSequence::iterator begin = cone->begin();
    992:           const typename sieve_type::traits::coneSequence::iterator end   = cone->end();

    994:           for(typename sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
    995:             sieveB->addArrow(*c_iter, *p_iter, c_iter.color());
    996:             next->insert(*c_iter);
    997:           }
    998:         }
    999:         tmp = current; current = next; next = tmp;
    1000:         next->clear();
    1001:       }
    1002:       if (!depth) {
    1003:         const Obj<typename sieve_type::traits::supportSequence>&     support = sieveA->support(p);
    1004:         const typename sieve_type::traits::supportSequence::iterator begin   = support->begin();
    1005:         const typename sieve_type::traits::supportSequence::iterator end     = support->end();

    1007:         for(typename sieve_type::traits::supportSequence::iterator s_iter = begin; s_iter != end; ++s_iter) {
    1008:           sieveB->addArrow(p, *s_iter, s_iter.color());
    1009:           next->insert(*s_iter);
    1010:         }
    1011:       }
    1012:     };
    1013:     static Obj<mesh_type> boundary_interpolated(const Obj<mesh_type>& mesh, const int faceHeight = 1) {
    1014:       Obj<mesh_type>                                     newMesh  = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
    1015:       Obj<sieve_type>                                    newSieve = new sieve_type(mesh->comm(), mesh->debug());
    1016:       const Obj<sieve_type>&                             sieve    = mesh->getSieve();
    1017:       const Obj<typename mesh_type::label_sequence>&     faces    = mesh->heightStratum(faceHeight);
    1018:       const typename mesh_type::label_sequence::iterator fBegin   = faces->begin();
    1019:       const typename mesh_type::label_sequence::iterator fEnd     = faces->end();
    1020:       const int                                          depth    = faceHeight - mesh->depth();

    1022:       for(typename mesh_type::label_sequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
    1023:         const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*f_iter);

    1025:         if (support->size() == 1) {
    1026:           addClosure(sieve, newSieve, *f_iter, depth);
    1027:         }
    1028:       }
    1029:       newMesh->setSieve(newSieve);
    1030:       newMesh->stratify();
    1031:       return newMesh;
    1032:     };
    1033:     template<typename SieveTypeA, typename SieveTypeB>
    1034:     static void addClosureV(const Obj<SieveTypeA>& sieveA, const Obj<SieveTypeB>& sieveB, const point_type& p, const int depth = 1) {
    1035:       typedef std::set<typename SieveTypeA::point_type> coneSet;
    1036:       ALE::ISieveVisitor::PointRetriever<SieveTypeA> cV(std::max(1, sieveA->getMaxConeSize()));
    1037:       Obj<coneSet> current = new coneSet();
    1038:       Obj<coneSet> next    = new coneSet();
    1039:       Obj<coneSet> tmp;

    1041:       current->insert(p);
    1042:       while(current->size()) {
    1043:         for(typename coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
    1044:           sieveA->cone(*p_iter, cV);
    1045:           const typename SieveTypeA::point_type *cone = cV.getPoints();

    1047:           for(int c = 0; c < (int) cV.getSize(); ++c) {
    1048:             sieveB->addArrow(cone[c], *p_iter);
    1049:             next->insert(cone[c]);
    1050:           }
    1051:           cV.clear();
    1052:         }
    1053:         tmp = current; current = next; next = tmp;
    1054:         next->clear();
    1055:       }
    1056:       if (!depth) {
    1057:         ALE::ISieveVisitor::PointRetriever<SieveTypeA> sV(std::max(1, sieveA->getMaxSupportSize()));

    1059:         sieveA->support(p, sV);
    1060:         const typename SieveTypeA::point_type *support = sV.getPoints();

    1062:         for(int s = 0; s < (int) sV.getSize(); ++s) {
    1063:           sieveB->addArrow(p, support[s]);
    1064:         }
    1065:         sV.clear();
    1066:       }
    1067:     };
    1068:   public:
    1069:     static Obj<mesh_type> boundary(const Obj<mesh_type>& mesh) {
    1070:       const int dim   = mesh->getDimension();
    1071:       const int depth = mesh->depth();

    1073:       if (dim == depth) {
    1074:         return boundary_interpolated(mesh);
    1075:       } else if (depth == dim+1) {
    1076:         return boundary_interpolated(mesh, 2);
    1077:       } else if (depth == 1) {
    1078:         return boundary_uninterpolated(mesh);
    1079:       } else if (depth == -1) {
    1080:         Obj<mesh_type>  newMesh  = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
    1081:         Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());

    1083:         newMesh->setSieve(newSieve);
    1084:         newMesh->stratify();
    1085:         return newMesh;
    1086:       }
    1087:       throw ALE::Exception("Cannot handle partially interpolated meshes");
    1088:     };
    1089:     template<typename MeshTypeQ>
    1090:     static Obj<FlexMesh> boundaryV_uninterpolated(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
    1091:         throw ALE::Exception("Cannot handle uninterpolated meshes");
    1092:     };
    1093:     // This method takes in an interpolated mesh, and returns the boundary
    1094:     template<typename MeshTypeQ>
    1095:     static Obj<FlexMesh> boundaryV_interpolated(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
    1096:       Obj<FlexMesh>                                       newMesh  = new FlexMesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
    1097:       Obj<typename FlexMesh::sieve_type>                  newSieve = new typename FlexMesh::sieve_type(mesh->comm(), mesh->debug());
    1098:       const Obj<typename MeshTypeQ::sieve_type>&          sieve    = mesh->getSieve();
    1099:       const Obj<typename MeshTypeQ::label_sequence>&      faces    = mesh->heightStratum(faceHeight);
    1100:       const typename MeshTypeQ::label_sequence::iterator  fBegin   = faces->begin();
    1101:       const typename MeshTypeQ::label_sequence::iterator  fEnd     = faces->end();
    1102:       const int                                           depth    = faceHeight - mesh->depth();
    1103:       ALE::ISieveVisitor::PointRetriever<sieve_type>      sV(std::max(1, sieve->getMaxSupportSize()));

    1105:       for(typename MeshTypeQ::label_sequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
    1106:         sieve->support(*f_iter, sV);

    1108:         if (sV.getSize() == 1) {
    1109:           addClosureV(sieve, newSieve, *f_iter, depth);
    1110:         }
    1111:         sV.clear();
    1112:       }
    1113:       newMesh->setSieve(newSieve);
    1114:       newMesh->stratify();
    1115:       return newMesh;
    1116:     };
    1117:     template<typename MeshTypeQ>
    1118:     static Obj<FlexMesh> boundaryV(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
    1119:       const int dim   = mesh->getDimension();
    1120:       const int depth = mesh->depth();

    1122:       if (dim == depth) {
    1123:         return boundaryV_interpolated(mesh);
    1124:       } else if (depth == dim+1) {
    1125:         return boundaryV_interpolated(mesh, 2);
    1126:       } else if (depth == 1) {
    1127:         throw ALE::Exception("Cannot handle uninterpolated meshes");
    1128:       } else if (depth == -1) {
    1129:         Obj<mesh_type>  newMesh  = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
    1130:         Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());

    1132:         newMesh->setSieve(newSieve);
    1133:         newMesh->stratify();
    1134:         return newMesh;
    1135:       }
    1136:       throw ALE::Exception("Cannot handle partially interpolated meshes");
    1137:     };
    1138:   public:
    1139:     static Obj<mesh_type> interpolateMesh(const Obj<mesh_type>& mesh) {
    1140:       const Obj<sieve_type>                              sieve       = mesh->getSieve();
    1141:       const int  dim         = mesh->getDimension();
    1142:       const int  numVertices = mesh->depthStratum(0)->size();
    1143:       const Obj<typename mesh_type::label_sequence>&     cells       = mesh->heightStratum(0);
    1144:       const int  numCells    = cells->size();
    1145:       const int  firstVertex = numCells;
    1146:       const int  debug       = sieve->debug();
    1147:       Obj<mesh_type>                                     newMesh     = new mesh_type(mesh->comm(), dim, mesh->debug());
    1148:       Obj<sieve_type>                                    newSieve    = new sieve_type(mesh->comm(), mesh->debug());
    1149:       const Obj<typename mesh_type::arrow_section_type>& orientation = newMesh->getArrowSection("orientation");

    1151:       newMesh->setSieve(newSieve);
    1152:       // Create a map from dimension to the current element number for that dimension
    1153:       std::map<int,point_type*> curElement;
    1154:       std::map<int,PointArray>  bdVertices;
    1155:       std::map<int,PointArray>  faces;
    1156:       std::map<int,oPointArray> oFaces;
    1157:       int                       curCell    = 0;
    1158:       int                       curVertex  = firstVertex;
    1159:       int                       newElement = firstVertex+numVertices;
    1160:       int                       o;

    1162:       curElement[0]   = &curVertex;
    1163:       curElement[dim] = &curCell;
    1164:       for(int d = 1; d < dim; d++) {
    1165:         curElement[d] = &newElement;
    1166:       }
    1167:       typename mesh_type::label_sequence::iterator cBegin = cells->begin();
    1168:       typename mesh_type::label_sequence::iterator cEnd   = cells->end();

    1170:       for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
    1171:         typename sieve_type::point_type                           cell    = *c_iter;
    1172:         const Obj<typename sieve_type::traits::coneSequence>&     cone    = sieve->cone(cell);
    1173:         const int                                                 corners = cone->size();
    1174:         const typename sieve_type::traits::coneSequence::iterator vBegin  = cone->begin();
    1175:         const typename sieve_type::traits::coneSequence::iterator vEnd    = cone->end();

    1177:         // Build the cell
    1178:         bdVertices[dim].clear();
    1179:         for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
    1180:           typename sieve_type::point_type vertex(*v_iter);

    1182:           if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
    1183:           bdVertices[dim].push_back(vertex);
    1184:         }
    1185:         if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}

    1187:         if (corners != dim+1) {
    1188:           ALE::SieveBuilder<mesh_type>::buildHexFaces(newSieve, dim, curElement, bdVertices, faces, cell);
    1189:           // Will need to handle cohesive cells here
    1190:         } else {
    1191:           ALE::SieveBuilder<mesh_type>::buildFaces(newSieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
    1192:         }
    1193:       }
    1194:       newMesh->stratify();
    1195:       newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
    1196:       return newMesh;
    1197:     };
    1198:   };
    1199: }

    1201: #if 0
    1202: namespace ALE {
    1203:   class MySelection {
    1204:   public:
    1205:     typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
    1206:     typedef ALE::Selection<ALE::Mesh> selection;
    1207:     typedef ALE::Mesh::sieve_type sieve_type;
    1208:     typedef ALE::Mesh::int_section_type int_section_type;
    1209:     typedef ALE::Mesh::real_section_type real_section_type;
    1210:     typedef std::set<ALE::Mesh::point_type> PointSet;
    1211:     typedef std::vector<ALE::Mesh::point_type> PointArray;
    1212:   public:
    1213:     template<class InputPoints>
    1214:     static bool _compatibleOrientation(const Obj<Mesh>& mesh,
    1215:                                        const ALE::Mesh::point_type& p,
    1216:                                        const ALE::Mesh::point_type& q,
    1217:                                        const int numFaultCorners,
    1218:                                        const int faultFaceSize,
    1219:                                        const int faultDepth,
    1220:                                        const Obj<InputPoints>& points,
    1221:                                        int indices[],
    1222:                                        PointArray *origVertices,
    1223:                                        PointArray *faceVertices,
    1224:                                        PointArray *neighborVertices)
    1225:     {
    1226:       typedef ALE::Selection<ALE::Mesh> selection;
    1227:       const int debug = mesh->debug();
    1228:       bool compatible;

    1230:       bool eOrient = selection::getOrientedFace(mesh, p, points, numFaultCorners, indices, origVertices, faceVertices);
    1231:       bool nOrient = selection::getOrientedFace(mesh, q, points, numFaultCorners, indices, origVertices, neighborVertices);

    1233:       if (faultFaceSize > 1) {
    1234:         if (debug) {
    1235:           for(PointArray::iterator v_iter = faceVertices->begin(); v_iter != faceVertices->end(); ++v_iter) {
    1236:             std::cout << "  face vertex " << *v_iter << std::endl;
    1237:           }
    1238:           for(PointArray::iterator v_iter = neighborVertices->begin(); v_iter != neighborVertices->end(); ++v_iter) {
    1239:             std::cout << "  neighbor vertex " << *v_iter << std::endl;
    1240:           }
    1241:         }
    1242:         compatible = !(*faceVertices->begin() == *neighborVertices->begin());
    1243:       } else {
    1244:         compatible = !(nOrient == eOrient);
    1245:       }
    1246:       return compatible;
    1247:     };
    1248:     static void _replaceCell(const Obj<sieve_type>& sieve,
    1249:                              const ALE::Mesh::point_type cell,
    1250:                              std::map<int,int> *vertexRenumber,
    1251:                              const int debug)
    1252:     {
    1253:       bool       replace = false;
    1254:       PointArray newVertices;

    1256:       const Obj<sieve_type::traits::coneSequence>& cCone = sieve->cone(cell);

    1258:       for(sieve_type::traits::coneSequence::iterator v_iter = cCone->begin(); v_iter != cCone->end(); ++v_iter) {
    1259:         if (vertexRenumber->find(*v_iter) != vertexRenumber->end()) {
    1260:           if (debug) std::cout << "    vertex " << (*vertexRenumber)[*v_iter] << std::endl;
    1261:           newVertices.insert(newVertices.end(), (*vertexRenumber)[*v_iter]);
    1262:           replace = true;
    1263:         } else {
    1264:           if (debug) std::cout << "    vertex " << *v_iter << std::endl;
    1265:           newVertices.insert(newVertices.end(), *v_iter);
    1266:         } // if/else
    1267:       } // for
    1268:       if (replace) {
    1269:         if (debug) std::cout << "  Replacing cell " << cell << std::endl;
    1270:         sieve->clearCone(cell);
    1271:         int color = 0;
    1272:         for(PointArray::const_iterator v_iter = newVertices.begin(); v_iter != newVertices.end(); ++v_iter) {
    1273:           sieve->addArrow(*v_iter, cell, color++);
    1274:         } // for
    1275:       }
    1276:     };
    1277:     template<class InputPoints>
    1278:     static void _computeCensoredDepth(const Obj<ALE::Mesh>& mesh,
    1279:                                       const Obj<ALE::Mesh::label_type>& depth,
    1280:                                       const Obj<ALE::Mesh::sieve_type>& sieve,
    1281:                                       const Obj<InputPoints>& points,
    1282:                                       const ALE::Mesh::point_type& firstCohesiveCell,
    1283:                                       const Obj<std::set<ALE::Mesh::point_type> >& modifiedPoints)
    1284:     {
    1285:       modifiedPoints->clear();

    1287:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
    1288:         if (*p_iter >= firstCohesiveCell) continue;
    1289:         // Compute the max depth of the points in the cone of p, and add 1
    1290:         int d0 = mesh->getValue(depth, *p_iter, -1);
    1291:         int d1 = mesh->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;

    1293:         if(d1 != d0) {
    1294:           mesh->setValue(depth, *p_iter, d1);
    1295:           modifiedPoints->insert(*p_iter);
    1296:         }
    1297:       }
    1298:       // FIX: We would like to avoid the copy here with support()
    1299:       if(modifiedPoints->size() > 0) {
    1300:         _computeCensoredDepth(mesh, depth, sieve, sieve->support(modifiedPoints), firstCohesiveCell, modifiedPoints);
    1301:       }
    1302:     };
    1303:     static void create(const Obj<Mesh>& mesh, Obj<Mesh> fault, const Obj<Mesh::int_section_type>& groupField) {
    1304:       static PetscLogEvent CreateFaultMesh_Event = 0, OrientFaultMesh_Event = 0, AddCohesivePoints_Event = 0, SplitMesh_Event = 0;

    1306:       if (!CreateFaultMesh_Event) {
    1307:         PetscLogEventRegister("CreateFaultMesh", 0,&CreateFaultMesh_Event);
    1308:       }
    1309:       if (!OrientFaultMesh_Event) {
    1310:         PetscLogEventRegister("OrientFaultMesh", 0,&OrientFaultMesh_Event);
    1311:       }
    1312:       if (!AddCohesivePoints_Event) {
    1313:         PetscLogEventRegister("AddCohesivePoints", 0,&AddCohesivePoints_Event);
    1314:       }
    1315:       if (!SplitMesh_Event) {
    1316:         PetscLogEventRegister("SplitMesh", 0,&SplitMesh_Event);
    1317:       }

    1319:       const Obj<sieve_type>& sieve = mesh->getSieve();
    1320:       const int  debug      = mesh->debug();
    1321:       int        numCorners = 0;    // The number of vertices in a mesh cell
    1322:       int        faceSize   = 0;    // The number of vertices in a mesh face
    1323:       int       *indices    = NULL; // The indices of a face vertex set in a cell
    1324:       PointArray origVertices;
    1325:       PointArray faceVertices;
    1326:       PointArray neighborVertices;
    1327:       const bool constraintCell = false;

    1329:       if (!mesh->commRank()) {
    1330:         numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), mesh->depth())->size();
    1331:         faceSize   = selection::numFaceVertices(mesh);
    1332:         indices    = new int[faceSize];
    1333:       }

    1335:       //int f = sieve->base()->size() + sieve->cap()->size();
    1336:       //ALE::Obj<PointSet> face = new PointSet();

    1338:       // Create a sieve which captures the fault
    1339:       PetscLogEventBegin(CreateFaultMesh_Event,0,0,0,0);
    1340:       fault = ALE::Selection<ALE::Mesh>::submesh(mesh, groupField);
    1341:       if (debug) {fault->view("Fault mesh");}
    1342:       PetscLogEventEnd(CreateFaultMesh_Event,0,0,0,0);
    1343:       // Orient the fault sieve
    1344:       PetscLogEventBegin(OrientFaultMesh_Event,0,0,0,0);
    1345:       const Obj<sieve_type>&                faultSieve = fault->getSieve();
    1346:       const ALE::Obj<Mesh::label_sequence>& fFaces     = fault->heightStratum(1);
    1347:       int faultDepth      = fault->depth()-1; // Depth of fault cells
    1348:       int numFaultCorners = 0; // The number of vertices in a fault cell

    1350:       if (!fault->commRank()) {
    1351:         numFaultCorners = faultSieve->nCone(*fFaces->begin(), faultDepth)->size();
    1352:         if (debug) std::cout << "  Fault corners " << numFaultCorners << std::endl;
    1353:         assert(numFaultCorners == faceSize);
    1354:       }
    1355:       PetscLogEventEnd(OrientFaultMesh_Event,0,0,0,0);

    1357:       // Add new shadow vertices and possibly Lagrange multipler vertices
    1358:       PetscLogEventBegin(AddCohesivePoints_Event,0,0,0,0);
    1359:       const ALE::Obj<Mesh::label_sequence>&   fVertices  = fault->depthStratum(0);
    1360:       const ALE::Obj<std::set<std::string> >& groupNames = mesh->getIntSections();
    1361:       Mesh::point_type newPoint = sieve->base()->size() + sieve->cap()->size();
    1362:       std::map<int,int> vertexRenumber;

    1364:       for(Mesh::label_sequence::iterator v_iter = fVertices->begin(); v_iter != fVertices->end(); ++v_iter, ++newPoint) {
    1365:         vertexRenumber[*v_iter] = newPoint;
    1366:         if (debug) {std::cout << "Duplicating " << *v_iter << " to " << vertexRenumber[*v_iter] << std::endl;}

    1368:         // Add shadow and constraint vertices (if they exist) to group
    1369:         // associated with fault
    1370:         groupField->addPoint(newPoint, 1);
    1371:         if (constraintCell) groupField->addPoint(newPoint+1, 1);

    1373:         // Add shadow vertices to other groups, don't add constraint
    1374:         // vertices (if they exist) because we don't want BC, etc to act
    1375:         // on constraint vertices
    1376:         for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
    1377:           const ALE::Obj<int_section_type>& group = mesh->getIntSection(*name);
    1378:           if (group->hasPoint(*v_iter)) group->addPoint(newPoint, 1);
    1379:         }
    1380:         if (constraintCell) newPoint++;
    1381:       }
    1382:       for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
    1383:         mesh->reallocate(mesh->getIntSection(*name));
    1384:       }

    1386:       // Split the mesh along the fault sieve and create cohesive elements
    1387:       const Obj<ALE::Mesh::label_sequence>&     faces       = fault->depthStratum(1);
    1388:       const Obj<ALE::Mesh::arrow_section_type>& orientation = mesh->getArrowSection("orientation");
    1389:       int firstCohesiveCell = newPoint;
    1390:       PointSet replaceCells;
    1391:       PointSet noReplaceCells;

    1393:       for(ALE::Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter, ++newPoint) {
    1394:         if (debug) std::cout << "Considering fault face " << *f_iter << std::endl;
    1395:         const ALE::Obj<sieve_type::traits::supportSequence>& cells = faultSieve->support(*f_iter);
    1396:         const ALE::Mesh::arrow_section_type::point_type arrow(*cells->begin(), *f_iter);
    1397:         bool reversed = orientation->restrictPoint(arrow)[0] < 0;
    1398:         const ALE::Mesh::point_type cell = *cells->begin();

    1400:         if (debug) std::cout << "  Checking orientation against cell " << cell << std::endl;
    1401:         if (numFaultCorners == 2) reversed = orientation->restrictPoint(arrow)[0] == -2;
    1402:         if (reversed) {
    1403:           replaceCells.insert(cell);
    1404:           noReplaceCells.insert(*(++cells->begin()));
    1405:         } else {
    1406:           replaceCells.insert(*(++cells->begin()));
    1407:           noReplaceCells.insert(cell);
    1408:         }
    1409:         //selection::getOrientedFace(mesh, cell, &vertexRenumber, numCorners, indices, &origVertices, &faceVertices);
    1410:         //const Obj<sieve_type::coneArray> faceCone = faultSieve->nCone(*f_iter, faultDepth);

    1412:         // Adding cohesive cell (not interpolated)
    1413:         const Obj<sieve_type::coneArray>&     fCone  = faultSieve->nCone(*f_iter, faultDepth);
    1414:         const sieve_type::coneArray::iterator fBegin = fCone->begin();
    1415:         const sieve_type::coneArray::iterator fEnd   = fCone->end();
    1416:         int color = 0;

    1418:         if (debug) {std::cout << "  Creating cohesive cell " << newPoint << std::endl;}
    1419:         for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
    1420:           if (debug) {std::cout << "    vertex " << *v_iter << std::endl;}
    1421:           sieve->addArrow(*v_iter, newPoint, color++);
    1422:         }
    1423:         for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
    1424:           if (debug) {std::cout << "    shadow vertex " << vertexRenumber[*v_iter] << std::endl;}
    1425:           sieve->addArrow(vertexRenumber[*v_iter], newPoint, color++);
    1426:         }
    1427:       }
    1428:       PetscLogEventEnd(AddCohesivePoints_Event,0,0,0,0);
    1429:       // Replace all cells with a vertex on the fault that share a face with this one, or with one that does
    1430:       PetscLogEventBegin(SplitMesh_Event,0,0,0,0);
    1431:       const int_section_type::chart_type&          chart    = groupField->getChart();
    1432:       const int_section_type::chart_type::iterator chartEnd = chart.end();

    1434:       for(PointSet::const_iterator v_iter = chart.begin(); v_iter != chartEnd; ++v_iter) {
    1435:         bool modified = true;

    1437:         while(modified) {
    1438:           modified = false;
    1439:           const Obj<sieve_type::traits::supportSequence>&     neighbors = sieve->support(*v_iter);
    1440:           const sieve_type::traits::supportSequence::iterator end       = neighbors->end();

    1442:           for(sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != end; ++n_iter) {
    1443:             if (replaceCells.find(*n_iter)   != replaceCells.end())   continue;
    1444:             if (noReplaceCells.find(*n_iter) != noReplaceCells.end()) continue;
    1445:             if (*n_iter >= firstCohesiveCell) continue;
    1446:             if (debug) std::cout << "  Checking fault neighbor " << *n_iter << std::endl;
    1447:             // If neighbors shares a faces with anyone in replaceCells, then add
    1448:             for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
    1449:               const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, mesh->depth());

    1451:               if ((int) preFace->size() == faceSize) {
    1452:                 if (debug) std::cout << "    Scheduling " << *n_iter << " for replacement" << std::endl;
    1453:                 replaceCells.insert(*n_iter);
    1454:                 modified = true;
    1455:                 break;
    1456:               }
    1457:             }
    1458:           }
    1459:         }
    1460:       }
    1461:       for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
    1462:         _replaceCell(sieve, *c_iter, &vertexRenumber, debug);
    1463:       }
    1464:       if (!fault->commRank()) delete [] indices;
    1465:       mesh->stratify();
    1466:       const ALE::Obj<Mesh::label_type>& label          = mesh->createLabel(std::string("censored depth"));
    1467:       const ALE::Obj<PointSet>          modifiedPoints = new PointSet();
    1468:       _computeCensoredDepth(mesh, label, mesh->getSieve(), mesh->getSieve()->roots(), firstCohesiveCell, modifiedPoints);
    1469:       if (debug) mesh->view("Mesh with Cohesive Elements");

    1471:       // Fix coordinates
    1472:       const Obj<real_section_type>&         coordinates = mesh->getRealSection("coordinates");
    1473:       const Obj<ALE::Mesh::label_sequence>& fVertices2  = fault->depthStratum(0);

    1475:       for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
    1476:         coordinates->addPoint(vertexRenumber[*v_iter], coordinates->getFiberDimension(*v_iter));
    1477:         if (constraintCell) {
    1478:           coordinates->addPoint(vertexRenumber[*v_iter]+1, coordinates->getFiberDimension(*v_iter));
    1479:         }
    1480:       }
    1481:       mesh->reallocate(coordinates);
    1482:       for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
    1483:         coordinates->updatePoint(vertexRenumber[*v_iter], coordinates->restrictPoint(*v_iter));
    1484:         if (constraintCell) {
    1485:         coordinates->updatePoint(vertexRenumber[*v_iter]+1, coordinates->restrictPoint(*v_iter));
    1486:         }
    1487:       }
    1488:       if (debug) coordinates->view("Coordinates with shadow vertices");
    1489:       PetscLogEventEnd(SplitMesh_Event,0,0,0,0);
    1490:     };
    1491:   };
    1492: };
    1493: #endif

    1495: #endif