Actual source code: Selection.hh

  1: #ifndef included_ALE_Selection_hh
  2: #define included_ALE_Selection_hh

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

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

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

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

136:       // Simplices
137:       if (cellDim == numCorners-1) {
138:         posOrient = !(oppositeVertex%2);
139:       } else if (cellDim == 2) {
140:         // Quads
141:         if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
142:           posOrient = true;
143:         } else if ((indices[0] == 3) && (indices[1] == 0)) {
144:           posOrient = true;
145:         } else {
146:           if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
147:             posOrient = false;
148:           } else {
149:             throw ALE::Exception("Invalid quad crossedge");
150:           }
151:         }
152:       } else if (cellDim == 3) {
153:         // Hexes
154:         //   A hex is two oriented quads with the normal of the first
155:         //   pointing up at the second.
156:         //
157:         //     7---6
158:         //    /|  /|
159:         //   4---5 |
160:         //   | 3-|-2
161:         //   |/  |/
162:         //   0---1
163:         int  sortedIndices[4];
164:         bool found = false;

166:         for(int i = 0; i < 4; ++i) sortedIndices[i] = indices[i];
167:         std::sort(sortedIndices, sortedIndices+4);
168:         // Case 1: Bottom quad
169:         if ((sortedIndices[0] == 0) && (sortedIndices[1] == 1) && (sortedIndices[2] == 2) && (sortedIndices[3] == 3)) {
170:           if (debug) std::cout << "Bottom quad" << std::endl;
171:           for(int i = 0; i < 4; ++i) {
172:             if (indices[i] == 3) {
173:               faceVertices->push_back((*origVertices)[i]); break;
174:             }
175:           }
176:           for(int i = 0; i < 4; ++i) {
177:             if (indices[i] == 2) {
178:               faceVertices->push_back((*origVertices)[i]); break;
179:             }
180:           }
181:           for(int i = 0; i < 4; ++i) {
182:             if (indices[i] == 1) {
183:               faceVertices->push_back((*origVertices)[i]); break;
184:             }
185:           }
186:           for(int i = 0; i < 4; ++i) {
187:             if (indices[i] == 0) {
188:               faceVertices->push_back((*origVertices)[i]); break;
189:             }
190:           }
191:           found = true;
192:         }
193:         // Case 2: Top quad
194:         if ((sortedIndices[0] == 4) && (sortedIndices[1] == 5) && (sortedIndices[2] == 6) && (sortedIndices[3] == 7)) {
195:           if (debug) std::cout << "Top quad" << std::endl;
196:           for(int i = 0; i < 4; ++i) {
197:             if (indices[i] == 5) {
198:               faceVertices->push_back((*origVertices)[i]); break;
199:             }
200:           }
201:           for(int i = 0; i < 4; ++i) {
202:             if (indices[i] == 6) {
203:               faceVertices->push_back((*origVertices)[i]); break;
204:             }
205:           }
206:           for(int i = 0; i < 4; ++i) {
207:             if (indices[i] == 7) {
208:               faceVertices->push_back((*origVertices)[i]); break;
209:             }
210:           }
211:           for(int i = 0; i < 4; ++i) {
212:             if (indices[i] == 4) {
213:               faceVertices->push_back((*origVertices)[i]); break;
214:             }
215:           }
216:           found = true;
217:         }
218:         // Case 3: Front quad
219:         if ((sortedIndices[0] == 0) && (sortedIndices[1] == 1) && (sortedIndices[2] == 4) && (sortedIndices[3] == 5)) {
220:           if (debug) std::cout << "Front quad" << std::endl;
221:           for(int i = 0; i < 4; ++i) {
222:             if (indices[i] == 1) {
223:               faceVertices->push_back((*origVertices)[i]); break;
224:             }
225:           }
226:           for(int i = 0; i < 4; ++i) {
227:             if (indices[i] == 5) {
228:               faceVertices->push_back((*origVertices)[i]); break;
229:             }
230:           }
231:           for(int i = 0; i < 4; ++i) {
232:             if (indices[i] == 4) {
233:               faceVertices->push_back((*origVertices)[i]); break;
234:             }
235:           }
236:           for(int i = 0; i < 4; ++i) {
237:             if (indices[i] == 0) {
238:               faceVertices->push_back((*origVertices)[i]); break;
239:             }
240:           }
241:           found = true;
242:         }
243:         // Case 4: Back quad
244:         if ((sortedIndices[0] == 2) && (sortedIndices[1] == 3) && (sortedIndices[2] == 6) && (sortedIndices[3] == 7)) {
245:           if (debug) std::cout << "Back quad" << std::endl;
246:           for(int i = 0; i < 4; ++i) {
247:             if (indices[i] == 7) {
248:               faceVertices->push_back((*origVertices)[i]); break;
249:             }
250:           }
251:           for(int i = 0; i < 4; ++i) {
252:             if (indices[i] == 6) {
253:               faceVertices->push_back((*origVertices)[i]); break;
254:             }
255:           }
256:           for(int i = 0; i < 4; ++i) {
257:             if (indices[i] == 2) {
258:               faceVertices->push_back((*origVertices)[i]); break;
259:             }
260:           }
261:           for(int i = 0; i < 4; ++i) {
262:             if (indices[i] == 3) {
263:               faceVertices->push_back((*origVertices)[i]); break;
264:             }
265:           }
266:           found = true;
267:         }
268:         // Case 5: Right quad
269:         if ((sortedIndices[0] == 1) && (sortedIndices[1] == 2) && (sortedIndices[2] == 5) && (sortedIndices[3] == 6)) {
270:           if (debug) std::cout << "Right quad" << std::endl;
271:           for(int i = 0; i < 4; ++i) {
272:             if (indices[i] == 2) {
273:               faceVertices->push_back((*origVertices)[i]); break;
274:             }
275:           }
276:           for(int i = 0; i < 4; ++i) {
277:             if (indices[i] == 6) {
278:               faceVertices->push_back((*origVertices)[i]); break;
279:             }
280:           }
281:           for(int i = 0; i < 4; ++i) {
282:             if (indices[i] == 5) {
283:               faceVertices->push_back((*origVertices)[i]); break;
284:             }
285:           }
286:           for(int i = 0; i < 4; ++i) {
287:             if (indices[i] == 1) {
288:               faceVertices->push_back((*origVertices)[i]); break;
289:             }
290:           }
291:           found = true;
292:         }
293:         // Case 6: Left quad
294:         if ((sortedIndices[0] == 0) && (sortedIndices[1] == 3) && (sortedIndices[2] == 4) && (sortedIndices[3] == 7)) {
295:           if (debug) std::cout << "Left quad" << std::endl;
296:           for(int i = 0; i < 4; ++i) {
297:             if (indices[i] == 4) {
298:               faceVertices->push_back((*origVertices)[i]); break;
299:             }
300:           }
301:           for(int i = 0; i < 4; ++i) {
302:             if (indices[i] == 7) {
303:               faceVertices->push_back((*origVertices)[i]); break;
304:             }
305:           }
306:           for(int i = 0; i < 4; ++i) {
307:             if (indices[i] == 3) {
308:               faceVertices->push_back((*origVertices)[i]); break;
309:             }
310:           }
311:           for(int i = 0; i < 4; ++i) {
312:             if (indices[i] == 0) {
313:               faceVertices->push_back((*origVertices)[i]); break;
314:             }
315:           }
316:           found = true;
317:         }
318:         if (!found) {throw ALE::Exception("Invalid hex crossface");}
319:         return true;
320:       }
321:       if (!posOrient) {
322:         if (debug) std::cout << "  Reversing initial face orientation" << std::endl;
323:         faceVertices->insert(faceVertices->end(), (*origVertices).rbegin(), (*origVertices).rend());
324:       } else {
325:         if (debug) std::cout << "  Keeping initial face orientation" << std::endl;
326:         faceVertices->insert(faceVertices->end(), (*origVertices).begin(), (*origVertices).end());
327:       }
328:       return posOrient;
329:     };
330:     // Given a cell and a face, as a set of vertices,
331:     //   return the oriented face, as a set of vertices, in faceVertices
332:     // The orientation is such that the face normal points out of the cell
333:     template<typename Mesh, typename FaceType>
334:     static bool getOrientedFace(const Obj<Mesh>& mesh, const point_type& cell, FaceType face,
335:                                 const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
336:     {
337:       FaceVisitor<typename Mesh::sieve_type,FaceType> v(face, *origVertices, *faceVertices, indices, mesh->debug());

339:       origVertices->clear();
340:       faceVertices->clear();
341:       mesh->getSieve()->cone(cell, v);
342:       return faceOrientation(cell, mesh, numCorners, indices, v.getOppositeVertex(), origVertices, faceVertices);
343:     };
344:     template<typename FaceType>
345:     static bool getOrientedFace(const Obj<ALE::Mesh>& mesh, const point_type& cell, FaceType face,
346:                                 const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
347:     {
348:       const Obj<typename ALE::Mesh::sieve_type::traits::coneSequence>& cone = mesh->getSieve()->cone(cell);
349:       const int debug = mesh->debug();
350:       int       v     = 0;
351:       int       oppositeVertex;

353:       origVertices->clear();
354:       faceVertices->clear();
355:       for(typename ALE::Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter, ++v) {
356:         if (face->find(*v_iter) != face->end()) {
357:           if (debug) std::cout << "    vertex " << *v_iter << std::endl;
358:           indices[origVertices->size()] = v;
359:           origVertices->insert(origVertices->end(), *v_iter);
360:         } else {
361:           if (debug) std::cout << "    vertex " << *v_iter << std::endl;
362:           oppositeVertex = v;
363:         }
364:       }
365:       return faceOrientation(cell, mesh, numCorners, indices, oppositeVertex, origVertices, faceVertices);
366:     };
367:     template<typename Sieve>
368:     static void insertFace(const Obj<mesh_type>& mesh, const Obj<Sieve>& subSieve, const Obj<PointSet>& face, point_type& f,
369:                            const point_type& cell, const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
370:     {
371:       const Obj<typename Sieve::supportSet> preFace = subSieve->nJoin1(face);
372:       const int                             debug   = subSieve->debug();

374:       if (preFace->size() > 1) {
375:         throw ALE::Exception("Invalid fault sieve: Multiple faces from vertex set");
376:       } else if (preFace->size() == 1) {
377:         // Add the other cell neighbor for this face
378:         subSieve->addArrow(*preFace->begin(), cell);
379:       } else if (preFace->size() == 0) {
380:         if (debug) std::cout << "  Orienting face " << f << std::endl;
381:         try {
382:           getOrientedFace(mesh, cell, face, numCorners, indices, origVertices, faceVertices);
383:           if (debug) std::cout << "  Adding face " << f << std::endl;
384:           int color = 0;
385:           for(typename PointArray::const_iterator f_iter = faceVertices->begin(); f_iter != faceVertices->end(); ++f_iter) {
386:             if (debug) std::cout << "    vertex " << *f_iter << std::endl;
387:             subSieve->addArrow(*f_iter, f, color++);
388:           }
389:           subSieve->addArrow(f, cell);
390:           f++;
391:         } catch (ALE::Exception e) {
392:           if (debug) std::cout << "  Did not add invalid face " << f << std::endl;
393:         }
394:       }
395:     };
396:   public:
397:     class FaceInserter {
398: #if 0
399:     public:
400:       typedef Mesh_                                mesh_type;
401:       typedef typename mesh_type::sieve_type       sieve_type;
402:       typedef typename mesh_type::point_type       point_type;
403:       typedef std::set<point_type>                 PointSet;
404:       typedef std::vector<point_type>              PointArray;
405: #endif
406:     protected:
407:       const Obj<mesh_type>  mesh;
408:       const Obj<sieve_type> sieve;
409:       const Obj<sieve_type> subSieve;
410:       point_type&           f;
411:       const point_type      cell;
412:       const int             numCorners;
413:       int                  *indices;
414:       PointArray           *origVertices;
415:       PointArray           *faceVertices;
416:       PointSet             *subCells;
417:       const int             debug;
418:     public:
419:       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()) {};
420:       virtual ~FaceInserter() {};
421:     public:
422:       void operator()(const Obj<PointArray>& face) {
423:         const Obj<typename sieve_type::supportSet> sievePreFace = sieve->nJoin1(face);

425:         if (sievePreFace->size() == 1) {
426:           if (debug) std::cout << "  Contains a boundary face on the submesh" << std::endl;
427:           PointSet faceSet(face->begin(), face->end());
428:           ALE::Selection<mesh_type>::insertFace(mesh, subSieve, faceSet, f, cell, numCorners, indices, origVertices, faceVertices);
429:           subCells->insert(cell);
430:         }
431:       };
432:     };
433:     template<typename Sieve>
434:     class FaceInserterV {
435:     protected:
436:       const Obj<mesh_type>&  mesh;
437:       const Obj<sieve_type>& sieve;
438:       const Obj<Sieve>&      subSieve;
439:       point_type&            f;
440:       const point_type       cell;
441:       const int              numCorners;
442:       int                   *indices;
443:       PointArray            *origVertices;
444:       PointArray            *faceVertices;
445:       PointSet              *subCells;
446:       const int              debug;
447:     public:
448:       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()) {};
449:       virtual ~FaceInserterV() {};
450:     public:
451:       void operator()(const Obj<PointArray>& face) {
452:         ISieveVisitor::PointRetriever<sieve_type> jV(sieve->getMaxSupportSize());

454:         sieve->join(*face, jV);
455:         if (jV.getSize() == 1) {
456:           if (debug) std::cout << "  Contains a boundary face on the submesh" << std::endl;
457:           PointSet faceSet(face->begin(), face->end());
458:           ALE::Selection<mesh_type>::insertFace(mesh, subSieve, faceSet, f, cell, numCorners, indices, origVertices, faceVertices);
459:           subCells->insert(cell);
460:         }
461:       };
462:     };
463:   protected:
464:     static int binomial(const int i, const int j) {
465:       assert(j <= i);
466:       assert(i < 34);
467:       if (j == 0) {
468:         return 1;
469:       } else if (j == i) {
470:         return 1;
471:       } else {
472:         return binomial(i-1, j) + binomial(i-1, j-1);
473:       }
474:     };
475:   public:
476:     // This takes in a section and creates a submesh from the vertices in the section chart
477:     //   This is a hyperplane of one dimension lower than the mesh
478:     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) {
479:       // A challenge here is to coordinate the extra numbering of new faces
480:       //   In serial, it is enough to number after the last point:
481:       //     Use sieve->base()->size() + sieve->cap()->size(), or determine the greatest point
482:       //   In parallel, there are two steps:
483:       //     1) Use the serial result, and reduce either with add (for size) or max (for greatest point)
484:       //     2) Determine how many faces will be created on each process
485:       //        This will be bounded by C(numCorners, faceSize)*submeshCells
486:       //        Thus it looks like we should first accumulate submeshCells, and then create faces
487:       typedef typename mesh_type::label_type        label_type;
488:       typedef typename int_section_type::chart_type chart_type;
489:       const int                  dim        = (dimension > 0) ? dimension : mesh->getDimension()-1;
490:       const Obj<sieve_type>&     sieve      = mesh->getSieve();
491:       Obj<mesh_type>             submesh    = new mesh_type(mesh->comm(), dim, mesh->debug());
492:       Obj<sieve_type>            subSieve   = new sieve_type(mesh->comm(), mesh->debug());
493:       const bool                 censor     = mesh->hasLabel("censored depth");
494:       const Obj<label_type>&     depthLabel = censor ? mesh->getLabel("censored depth") : mesh->getLabel("depth");
495:       const int                  depth      = mesh->depth();
496:       const int                  height     = mesh->height();
497:       const chart_type&          chart      = label->getChart();
498:       const int                  numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), depth)->size();
499:       const int                  faceSize   = numFaceVertices(mesh);
500:       Obj<PointSet>              face       = new PointSet();
501:       int                        f          = sieve->base()->size() + sieve->cap()->size();
502:       const int                  debug      = mesh->debug();
503:       int                       *indices    = new int[faceSize];
504:       PointArray                 origVertices, faceVertices;
505:       PointSet                   submeshVertices, submeshCells;


508:       const typename chart_type::iterator chartEnd = chart.end();
509:       for(typename chart_type::iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
510:         if (label->getFiberDimension(*c_iter)) submeshVertices.insert(*c_iter);
511:       }
512:       const typename PointSet::const_iterator svBegin = submeshVertices.begin();
513:       const typename PointSet::const_iterator svEnd   = submeshVertices.end();

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

520:         if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
521:         for(typename sieveAlg::supportArray::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
522:           if (debug) std::cout << "  Checking cell " << *c_iter << std::endl;
523:           if (submeshCells.find(*c_iter) != submeshCells.end())        continue;
524:           if (censor && (!mesh->getValue(depthLabel, *c_iter))) continue;
525:           const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *c_iter, height);
526:           const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
527:           const typename sieveAlg::coneArray::iterator vEnd   = cone->end();

529:           face->clear();
530:           for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
531:             if (submeshVertices.find(*v_iter) != svEnd) {
532:               if (debug) std::cout << "    contains submesh vertex " << *v_iter << std::endl;
533:               face->insert(face->end(), *v_iter);
534:             }
535:           }
536:           if ((int) face->size() > faceSize) {
537:             if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
538:             if (debug) std::cout << "  Has all boundary faces on the submesh" << std::endl;
539:             submeshCells.insert(*c_iter);
540:           }
541:           if ((int) face->size() == faceSize) {
542:             if (debug) std::cout << "  Contains a face on the submesh" << std::endl;
543:             submeshCells.insert(*c_iter);
544:           }
545:         }
546:       }
547:       if (mesh->commSize() > 1) {
548:         int localF     = f;
549:         int localFaces = binomial(numCorners, faceSize)*submeshCells.size();
550:         int maxFaces;

552:         MPI_Allreduce(&localF, &f, 1, MPI_INT, MPI_SUM, mesh->comm());
553:         //     2) Determine how many faces will be created on each process
554:         //        This will be bounded by faceSize*submeshCells
555:         //        Thus it looks like we should first accumulate submeshCells, and then create faces
556:         MPI_Scan(&localFaces, &maxFaces, 1, MPI_INT, MPI_SUM, mesh->comm());
557:         f += maxFaces - localFaces;
558:       }
559:       for(typename PointSet::const_iterator c_iter = submeshCells.begin(); c_iter != submeshCells.end(); ++c_iter) {
560:         if (debug) std::cout << "  Processing submesh cell " << *c_iter << std::endl;
561:         const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *c_iter, height);
562:         const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
563:         const typename sieveAlg::coneArray::iterator vEnd   = cone->end();

565:         face->clear();
566:         for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
567:           if (submeshVertices.find(*v_iter) != svEnd) {
568:             if (debug) std::cout << "    contains submesh vertex " << *v_iter << std::endl;
569:             face->insert(face->end(), *v_iter);
570:           }
571:         }
572:         if ((int) face->size() > faceSize) {
573:           // Here we allow a set of vertices to lie completely on a boundary cell (like a corner tetrahedron)
574:           //   We have to take all the faces, and discard those in the interior
575:           FaceInserter inserter(mesh, sieve, subSieve, f, *c_iter, numCorners, indices, &origVertices, &faceVertices, &submeshCells);
576:           PointArray   faceVec(face->begin(), face->end());

578:           subsets(faceVec, faceSize, inserter);
579:         }
580:         if ((int) face->size() == faceSize) {
581:           insertFace(mesh, subSieve, face, f, *c_iter, numCorners, indices, &origVertices, &faceVertices);
582:         }
583:       }
584:       delete [] indices;
585:       submesh->setSieve(subSieve);
586:       submesh->stratify();
587:       if (debug) submesh->view("Submesh");
588:       return submesh;
589:     };
590:     // This takes in a section and creates a submesh from the vertices in the section chart
591:     //   This is a hyperplane of one dimension lower than the mesh
592:     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) {
593:       const int debug  = mesh->debug();
594:       const int depth  = mesh->depth();
595:       const int height = mesh->height();
596:       const typename int_section_type::chart_type&                chart        = label->getChart();
597:       const typename int_section_type::chart_type::const_iterator chartEnd     = chart.end();
598:       const Obj<PointSet>                                         submeshFaces = new PointSet();
599:       PointSet submeshVertices;

601:       for(typename int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
602:         //assert(!mesh->depth(*c_iter));
603:         submeshVertices.insert(*c_iter);
604:       }
605:       const typename PointSet::const_iterator svBegin = submeshVertices.begin();
606:       const typename PointSet::const_iterator svEnd   = submeshVertices.end();

608:       for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
609:         const Obj<typename sieveAlg::supportArray>& faces = sieveAlg::nSupport(mesh, *sv_iter, depth-1);
610:         const typename sieveAlg::supportArray::iterator fBegin = faces->begin();
611:         const typename sieveAlg::supportArray::iterator fEnd   = faces->end();
612: 
613:         if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
614:         for(typename sieveAlg::supportArray::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
615:           if (debug) std::cout << "  Checking face " << *f_iter << std::endl;
616:           if (submeshFaces->find(*f_iter) != submeshFaces->end())        continue;
617:           const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *f_iter, height-1);
618:           const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
619:           const typename sieveAlg::coneArray::iterator vEnd   = cone->end();
620:           bool                                         found  = true;

622:           for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
623:             if (submeshVertices.find(*v_iter) != svEnd) {
624:               if (debug) std::cout << "    contains submesh vertex " << *v_iter << std::endl;
625:             } else {
626:               found = false;
627:             }
628:           }
629:           if (found) {
630:             if (boundaryFaces) {throw ALE::Exception("Not finished: should check that it is a boundary face");}
631:             if (debug) std::cout << "  Is a face on the submesh" << std::endl;
632:             submeshFaces->insert(*f_iter);
633:           }
634:         }
635:       }
636:       return submesh(mesh, submeshFaces, mesh->getDimension()-1);
637:     };
638:     template<typename output_mesh_type>
639:     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) {
640:       typedef typename mesh_type::label_type        label_type;
641:       typedef typename int_section_type::chart_type chart_type;
642:       const int                           dim        = (dimension > 0) ? dimension : mesh->getDimension()-1;
643:       const Obj<sieve_type>&              sieve      = mesh->getSieve();
644:       Obj<ALE::Mesh>                      submesh    = new ALE::Mesh(mesh->comm(), dim, mesh->debug());
645:       Obj<typename ALE::Mesh::sieve_type> subSieve   = new typename ALE::Mesh::sieve_type(mesh->comm(), mesh->debug());
646:       const bool                          censor     = mesh->hasLabel("censored depth");
647:       const Obj<label_type>&              depthLabel = censor ? mesh->getLabel("censored depth") : mesh->getLabel("depth");
648:       const chart_type&                   chart      = label->getChart();
649:       const int                           numCorners = mesh->getNumCellCorners();
650:       const int                           faceSize   = numFaceVertices(mesh);
651:       Obj<PointSet>                       face       = new PointSet();
652:       int                                 f          = sieve->getBaseSize() + sieve->getCapSize();
653:       const int                           debug      = mesh->debug();
654:       int                                *indices    = new int[faceSize];
655:       PointArray                          origVertices, faceVertices;
656:       PointSet                            submeshVertices, submeshCells;

658:       const typename chart_type::const_iterator chartEnd = chart.end();
659:       for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
660:         if (label->getFiberDimension(*c_iter)) submeshVertices.insert(*c_iter);
661:       }
662:       const typename PointSet::const_iterator svBegin = submeshVertices.begin();
663:       const typename PointSet::const_iterator svEnd   = submeshVertices.end();
664:       typename ISieveVisitor::PointRetriever<sieve_type> sV(sieve->getMaxSupportSize());
665:       typename ISieveVisitor::PointRetriever<sieve_type> cV(sieve->getMaxConeSize());

667:       for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
668:         sieve->support(*sv_iter, sV);
669:         const int         numCells = sV.getSize();
670:         const point_type *cells    = sV.getPoints();
671: 
672:         if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
673:         for(int c = 0; c < numCells; ++c) {
674:           if (debug) std::cout << "  Checking cell " << cells[c] << std::endl;
675:           if (submeshCells.find(cells[c]) != submeshCells.end()) continue;
676:           if (censor && (!mesh->getValue(depthLabel, cells[c]))) continue;
677:           sieve->cone(cells[c], cV);
678:           const int         numVertices = cV.getSize();
679:           const point_type *vertices    = cV.getPoints();

681:           face->clear();
682:           for(int v = 0; v < numVertices; ++v) {
683:             if (submeshVertices.find(vertices[v]) != svEnd) {
684:               if (debug) std::cout << "    contains submesh vertex " << vertices[v] << std::endl;
685:               face->insert(face->end(), vertices[v]);
686:             }
687:           }
688:           if ((int) face->size() > faceSize) {
689:             if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
690:             if (debug) std::cout << "  Has all boundary faces on the submesh" << std::endl;
691:             submeshCells.insert(cells[c]);
692:           }
693:           if ((int) face->size() == faceSize) {
694:             if (debug) std::cout << "  Contains a face on the submesh" << std::endl;
695:             submeshCells.insert(cells[c]);
696:           }
697:           cV.clear();
698:         }
699:         sV.clear();
700:       }
701:       if (mesh->commSize() > 1) {
702:         int localF     = f;
703:         int localFaces = binomial(numCorners, faceSize)*submeshCells.size();
704:         int maxFaces;

706:         MPI_Allreduce(&localF, &f, 1, MPI_INT, MPI_SUM, mesh->comm());
707:         //     2) Determine how many faces will be created on each process
708:         //        This will be bounded by faceSize*submeshCells
709:         //        Thus it looks like we should first accumulate submeshCells, and then create faces
710:         MPI_Scan(&localFaces, &maxFaces, 1, MPI_INT, MPI_SUM, mesh->comm());
711:         f += maxFaces - localFaces;
712:       }
713:       for(typename PointSet::const_iterator c_iter = submeshCells.begin(); c_iter != submeshCells.end(); ++c_iter) {
714:         if (debug) std::cout << "  Processing submesh cell " << *c_iter << std::endl;
715:         sieve->cone(*c_iter, cV);
716:         const int         numVertices = cV.getSize();
717:         const point_type *vertices    = cV.getPoints();

719:         face->clear();
720:         for(int v = 0; v < numVertices; ++v) {
721:           if (submeshVertices.find(vertices[v]) != svEnd) {
722:             if (debug) std::cout << "    contains submesh vertex " << vertices[v] << std::endl;
723:             face->insert(face->end(), vertices[v]);
724:           }
725:         }
726:         if ((int) face->size() > faceSize) {
727:           if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
728:           // Here we allow a set of vertices to lie completely on a boundary cell (like a corner tetrahedron)
729:           //   We have to take all the faces, and discard those in the interior
730:           FaceInserterV<ALE::Mesh::sieve_type> inserter(mesh, sieve, subSieve, f, *c_iter, numCorners, indices, &origVertices, &faceVertices, &submeshCells);
731:           PointArray                           faceVec(face->begin(), face->end());

733:           subsets(faceVec, faceSize, inserter);
734:         }
735:         if ((int) face->size() == faceSize) {
736:           insertFace(mesh, subSieve, face, f, *c_iter, numCorners, indices, &origVertices, &faceVertices);
737:         }
738:         cV.clear();
739:       }
740:       delete [] indices;
741:       submesh->setSieve(subSieve);
742:       submesh->stratify();
743:       if (debug) submesh->view("Submesh");

745:       Obj<output_mesh_type> isubmesh = new output_mesh_type(submesh->comm(), submesh->getDimension(), submesh->debug());
746:       Obj<typename output_mesh_type::sieve_type> isieve = new typename output_mesh_type::sieve_type(submesh->comm(), submesh->debug());
747:       std::map<typename output_mesh_type::point_type,typename output_mesh_type::point_type> renumbering;
748:       isubmesh->setSieve(isieve);
749:       ALE::ISieveConverter::convertMesh(*submesh, *isubmesh, renumbering, false);
750:       return isubmesh;
751:     };
752:   public:
753:     // This takes in a section and creates a submesh from the vertices in the section chart
754:     //   This is a hyperplane of one dimension lower than the mesh
755:     static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
756:       const int dim   = mesh->getDimension();
757:       const int depth = mesh->depth();

759:       if (dim == depth) {
760:         return submesh_interpolated(mesh, label, dimension, false);
761:       } else if (depth == 1) {
762:         return submesh_uninterpolated(mesh, label, dimension);
763:       }
764:       throw ALE::Exception("Cannot handle partially interpolated meshes");
765:     };
766:     template<typename output_mesh_type>
767:     static Obj<output_mesh_type> submeshV(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
768:       const int dim   = mesh->getDimension();
769:       const int depth = mesh->depth();

771: #if 0
772:       if (dim == depth) {
773:         //return submesh_interpolated(mesh, label, dimension, false);
774:         throw ALE::Exception("Cannot handle interpolated meshes");
775:       } else if (depth == 1) {
776:         return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
777:       }
778: #else
779:       if (depth == 1) {
780:         return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
781:       } else if (dim == depth) {
782:         //return submesh_interpolated(mesh, label, dimension, false);
783:         throw ALE::Exception("Cannot handle interpolated meshes");
784:       }
785: #endif
786:       throw ALE::Exception("Cannot handle partially interpolated meshes");
787:     };
788:     // This creates a submesh consisting of the union of the closures of the given points
789:     //   This mesh is the same dimension as in the input mesh
790:     template<typename Points>
791:     static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<Points>& points, const int dim = -1) {
792:       const Obj<sieve_type>& sieve     = mesh->getSieve();
793:       Obj<mesh_type>         newMesh   = new mesh_type(mesh->comm(), dim >= 0 ? dim : mesh->getDimension(), mesh->debug());
794:       Obj<sieve_type>        newSieve  = new sieve_type(mesh->comm(), mesh->debug());
795:       Obj<PointSet>          newPoints = new PointSet();
796:       Obj<PointSet>          modPoints = new PointSet();
797:       Obj<PointSet>          tmpPoints;

799:       newMesh->setSieve(newSieve);
800:       for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
801:         newPoints->insert(*p_iter);
802:         do {
803:           modPoints->clear();
804:           for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
805:             const Obj<typename sieve_type::traits::coneSequence>&     cone = sieve->cone(*np_iter);
806:             const typename sieve_type::traits::coneSequence::iterator end  = cone->end();
807:             int c = 0;

809:             for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter, ++c) {
810:               newSieve->addArrow(*c_iter, *np_iter, c);
811:             }
812:             modPoints->insert(cone->begin(), cone->end());
813:           }
814:           tmpPoints = newPoints;
815:           newPoints = modPoints;
816:           modPoints = tmpPoints;
817:         } while(newPoints->size());
818:         newPoints->insert(*p_iter);
819:         do {
820:           modPoints->clear();
821:           for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
822:             const Obj<typename sieve_type::traits::supportSequence>&     support = sieve->support(*np_iter);
823:             const typename sieve_type::traits::supportSequence::iterator end     = support->end();
824:             int s = 0;

826:             for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != end; ++s_iter, ++s) {
827:               newSieve->addArrow(*np_iter, *s_iter, s);
828:             }
829:             modPoints->insert(support->begin(), support->end());
830:           }
831:           tmpPoints = newPoints;
832:           newPoints = modPoints;
833:           modPoints = tmpPoints;
834:         } while(newPoints->size());
835:       }
836:       newMesh->stratify();
837:       newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
838:       newMesh->setArrowSection("orientation", mesh->getArrowSection("orientation"));
839:       return newMesh;
840:     };
841:   protected:
842:     static Obj<mesh_type> boundary_uninterpolated(const Obj<mesh_type>& mesh) {
843:       throw ALE::Exception("Not yet implemented");
844:       const Obj<typename mesh_type::label_sequence>&     cells    = mesh->heightStratum(0);
845:       const Obj<sieve_type>&                             sieve    = mesh->getSieve();
846:       const typename mesh_type::label_sequence::iterator cBegin   = cells->begin();
847:       const typename mesh_type::label_sequence::iterator cEnd     = cells->end();
848:       const int                                          faceSize = numFaceVertices(mesh);

850:       for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
851:         const Obj<typename sieve_type::traits::coneSequence>&     vertices = sieve->cone(*c_iter);
852:         const typename sieve_type::traits::coneSequence::iterator vBegin   = vertices->begin();
853:         const typename sieve_type::traits::coneSequence::iterator vEnd     = vertices->end();
854:         //PointArray cell(vertices->begin(), vertices->end());

856:         // For each vertex, gather
857:         for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
858:           const Obj<typename sieve_type::traits::supportSequence>&     neighbors = sieve->support(*v_iter);
859:           const typename sieve_type::traits::supportSequence::iterator nBegin    = neighbors->begin();
860:           const typename sieve_type::traits::supportSequence::iterator nEnd      = neighbors->end();

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

865:             if (preFace->size() == faceSize) {
866:             }
867:           }
868:         }
869:         // For each face
870:         // - determine if its legal
871: 
872:         // - determine if its part of a neighboring cell
873:         // - if not, its a boundary face
874:         //subsets(cell, faceSize, inserter);
875:       }
876:     };
877:     static void addClosure(const Obj<sieve_type>& sieveA, const Obj<sieve_type>& sieveB, const point_type& p, const int depth = 1) {
878:       Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
879:       Obj<typename sieve_type::coneSet> next    = new typename sieve_type::coneSet();
880:       Obj<typename sieve_type::coneSet> tmp;

882:       current->insert(p);
883:       while(current->size()) {
884:         for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
885:           const Obj<typename sieve_type::traits::coneSequence>&     cone  = sieveA->cone(*p_iter);
886:           const typename sieve_type::traits::coneSequence::iterator begin = cone->begin();
887:           const typename sieve_type::traits::coneSequence::iterator end   = cone->end();

889:           for(typename sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
890:             sieveB->addArrow(*c_iter, *p_iter, c_iter.color());
891:             next->insert(*c_iter);
892:           }
893:         }
894:         tmp = current; current = next; next = tmp;
895:         next->clear();
896:       }
897:       if (!depth) {
898:         const Obj<typename sieve_type::traits::supportSequence>&     support = sieveA->support(p);
899:         const typename sieve_type::traits::supportSequence::iterator begin   = support->begin();
900:         const typename sieve_type::traits::supportSequence::iterator end     = support->end();
901: 
902:         for(typename sieve_type::traits::supportSequence::iterator s_iter = begin; s_iter != end; ++s_iter) {
903:           sieveB->addArrow(p, *s_iter, s_iter.color());
904:           next->insert(*s_iter);
905:         }
906:       }
907:     };
908:     static Obj<mesh_type> boundary_interpolated(const Obj<mesh_type>& mesh, const int faceHeight = 1) {
909:       Obj<mesh_type>                                     newMesh  = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
910:       Obj<sieve_type>                                    newSieve = new sieve_type(mesh->comm(), mesh->debug());
911:       const Obj<sieve_type>&                             sieve    = mesh->getSieve();
912:       const Obj<typename mesh_type::label_sequence>&     faces    = mesh->heightStratum(faceHeight);
913:       const typename mesh_type::label_sequence::iterator fBegin   = faces->begin();
914:       const typename mesh_type::label_sequence::iterator fEnd     = faces->end();
915:       const int                                          depth    = faceHeight - mesh->depth();

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

920:         if (support->size() == 1) {
921:           addClosure(sieve, newSieve, *f_iter, depth);
922:         }
923:       }
924:       newMesh->setSieve(newSieve);
925:       newMesh->stratify();
926:       return newMesh;
927:     };
928:     template<typename SieveTypeA, typename SieveTypeB>
929:     static void addClosureV(const Obj<SieveTypeA>& sieveA, const Obj<SieveTypeB>& sieveB, const point_type& p, const int depth = 1) {
930:       typedef std::set<typename SieveTypeA::point_type> coneSet;
931:       ALE::ISieveVisitor::PointRetriever<SieveTypeA> cV(std::max(1, sieveA->getMaxConeSize()));
932:       Obj<coneSet> current = new coneSet();
933:       Obj<coneSet> next    = new coneSet();
934:       Obj<coneSet> tmp;

936:       current->insert(p);
937:       while(current->size()) {
938:         for(typename coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
939:           sieveA->cone(*p_iter, cV);
940:           const typename SieveTypeA::point_type *cone = cV.getPoints();

942:           for(int c = 0; c < (int) cV.getSize(); ++c) {
943:             sieveB->addArrow(cone[c], *p_iter);
944:             next->insert(cone[c]);
945:           }
946:           cV.clear();
947:         }
948:         tmp = current; current = next; next = tmp;
949:         next->clear();
950:       }
951:       if (!depth) {
952:         ALE::ISieveVisitor::PointRetriever<SieveTypeA> sV(std::max(1, sieveA->getMaxSupportSize()));

954:         sieveA->support(p, sV);
955:         const typename SieveTypeA::point_type *support = sV.getPoints();
956: 
957:         for(int s = 0; s < (int) sV.getSize(); ++s) {
958:           sieveB->addArrow(p, support[s]);
959:         }
960:         sV.clear();
961:       }
962:     };
963:   public:
964:     static Obj<mesh_type> boundary(const Obj<mesh_type>& mesh) {
965:       const int dim   = mesh->getDimension();
966:       const int depth = mesh->depth();

968:       if (dim == depth) {
969:         return boundary_interpolated(mesh);
970:       } else if (depth == dim+1) {
971:         return boundary_interpolated(mesh, 2);
972:       } else if (depth == 1) {
973:         return boundary_uninterpolated(mesh);
974:       } else if (depth == -1) {
975:         Obj<mesh_type>  newMesh  = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
976:         Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());

978:         newMesh->setSieve(newSieve);
979:         newMesh->stratify();
980:         return newMesh;
981:       }
982:       throw ALE::Exception("Cannot handle partially interpolated meshes");
983:     };
984:     template<typename MeshTypeQ>
985:     static Obj<ALE::Mesh> boundaryV_uninterpolated(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
986:         throw ALE::Exception("Cannot handle uninterpolated meshes");
987:     };
988:     // This method takes in an interpolated mesh, and returns the boundary
989:     template<typename MeshTypeQ>
990:     static Obj<ALE::Mesh> boundaryV_interpolated(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
991:       Obj<ALE::Mesh>                                      newMesh  = new ALE::Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
992:       Obj<typename ALE::Mesh::sieve_type>                 newSieve = new typename ALE::Mesh::sieve_type(mesh->comm(), mesh->debug());
993:       const Obj<typename MeshTypeQ::sieve_type>&          sieve    = mesh->getSieve();
994:       const Obj<typename MeshTypeQ::label_sequence>&      faces    = mesh->heightStratum(faceHeight);
995:       const typename MeshTypeQ::label_sequence::iterator  fBegin   = faces->begin();
996:       const typename MeshTypeQ::label_sequence::iterator  fEnd     = faces->end();
997:       const int                                           depth    = faceHeight - mesh->depth();
998:       ALE::ISieveVisitor::PointRetriever<sieve_type>      sV(std::max(1, sieve->getMaxSupportSize()));

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

1003:         if (sV.getSize() == 1) {
1004:           addClosureV(sieve, newSieve, *f_iter, depth);
1005:         }
1006:         sV.clear();
1007:       }
1008:       newMesh->setSieve(newSieve);
1009:       newMesh->stratify();
1010:       return newMesh;
1011:     };
1012:     template<typename MeshTypeQ>
1013:     static Obj<ALE::Mesh> boundaryV(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
1014:       const int dim   = mesh->getDimension();
1015:       const int depth = mesh->depth();

1017:       if (dim == depth) {
1018:         return boundaryV_interpolated(mesh);
1019:       } else if (depth == dim+1) {
1020:         return boundaryV_interpolated(mesh, 2);
1021:       } else if (depth == 1) {
1022:         throw ALE::Exception("Cannot handle uninterpolated meshes");
1023:       } else if (depth == -1) {
1024:         Obj<mesh_type>  newMesh  = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
1025:         Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());

1027:         newMesh->setSieve(newSieve);
1028:         newMesh->stratify();
1029:         return newMesh;
1030:       }
1031:       throw ALE::Exception("Cannot handle partially interpolated meshes");
1032:     };
1033:   public:
1034:     static Obj<mesh_type> interpolateMesh(const Obj<mesh_type>& mesh) {
1035:       const Obj<sieve_type>                              sieve       = mesh->getSieve();
1036:       const int  dim         = mesh->getDimension();
1037:       const int  numVertices = mesh->depthStratum(0)->size();
1038:       const Obj<typename mesh_type::label_sequence>&     cells       = mesh->heightStratum(0);
1039:       const int  numCells    = cells->size();
1040:       const int  corners     = sieve->cone(*cells->begin())->size();
1041:       const int  firstVertex = numCells;
1042:       const int  debug       = sieve->debug();
1043:       Obj<mesh_type>                                     newMesh     = new mesh_type(mesh->comm(), dim, mesh->debug());
1044:       Obj<sieve_type>                                    newSieve    = new sieve_type(mesh->comm(), mesh->debug());
1045:       const Obj<typename mesh_type::arrow_section_type>& orientation = newMesh->getArrowSection("orientation");

1047:       newMesh->setSieve(newSieve);
1048:       // Create a map from dimension to the current element number for that dimension
1049:       std::map<int,point_type*> curElement;
1050:       std::map<int,PointArray>  bdVertices;
1051:       std::map<int,PointArray>  faces;
1052:       std::map<int,oPointArray> oFaces;
1053:       int                       curCell    = 0;
1054:       int                       curVertex  = firstVertex;
1055:       int                       newElement = firstVertex+numVertices;
1056:       int                       o;

1058:       curElement[0]   = &curVertex;
1059:       curElement[dim] = &curCell;
1060:       for(int d = 1; d < dim; d++) {
1061:         curElement[d] = &newElement;
1062:       }
1063:       typename mesh_type::label_sequence::iterator cBegin = cells->begin();
1064:       typename mesh_type::label_sequence::iterator cEnd   = cells->end();

1066:       for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
1067:         typename sieve_type::point_type                           cell   = *c_iter;
1068:         const Obj<typename sieve_type::traits::coneSequence>&     cone   = sieve->cone(cell);
1069:         const typename sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
1070:         const typename sieve_type::traits::coneSequence::iterator vEnd   = cone->end();

1072:         // Build the cell
1073:         bdVertices[dim].clear();
1074:         for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
1075:           typename sieve_type::point_type vertex(*v_iter);

1077:           if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
1078:           bdVertices[dim].push_back(vertex);
1079:         }
1080:         if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}

1082:         if (corners != dim+1) {
1083:           ALE::SieveBuilder<mesh_type>::buildHexFaces(newSieve, dim, curElement, bdVertices, faces, cell);
1084:         } else {
1085:           ALE::SieveBuilder<mesh_type>::buildFaces(newSieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
1086:         }
1087:       }
1088:       newMesh->stratify();
1089:       newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
1090:       return newMesh;
1091:     };
1092:   };
1093: }

1095: #if 0
1096: namespace ALE {
1097:   class MySelection {
1098:   public:
1099:     typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
1100:     typedef ALE::Selection<ALE::Mesh> selection;
1101:     typedef ALE::Mesh::sieve_type sieve_type;
1102:     typedef ALE::Mesh::int_section_type int_section_type;
1103:     typedef ALE::Mesh::real_section_type real_section_type;
1104:     typedef std::set<ALE::Mesh::point_type> PointSet;
1105:     typedef std::vector<ALE::Mesh::point_type> PointArray;
1106:   public:
1107:     template<class InputPoints>
1108:     static bool _compatibleOrientation(const Obj<Mesh>& mesh,
1109:                                        const ALE::Mesh::point_type& p,
1110:                                        const ALE::Mesh::point_type& q,
1111:                                        const int numFaultCorners,
1112:                                        const int faultFaceSize,
1113:                                        const int faultDepth,
1114:                                        const Obj<InputPoints>& points,
1115:                                        int indices[],
1116:                                        PointArray *origVertices,
1117:                                        PointArray *faceVertices,
1118:                                        PointArray *neighborVertices)
1119:     {
1120:       typedef ALE::Selection<ALE::Mesh> selection;
1121:       const int debug = mesh->debug();
1122:       bool compatible;

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

1127:       if (faultFaceSize > 1) {
1128:         if (debug) {
1129:           for(PointArray::iterator v_iter = faceVertices->begin(); v_iter != faceVertices->end(); ++v_iter) {
1130:             std::cout << "  face vertex " << *v_iter << std::endl;
1131:           }
1132:           for(PointArray::iterator v_iter = neighborVertices->begin(); v_iter != neighborVertices->end(); ++v_iter) {
1133:             std::cout << "  neighbor vertex " << *v_iter << std::endl;
1134:           }
1135:         }
1136:         compatible = !(*faceVertices->begin() == *neighborVertices->begin());
1137:       } else {
1138:         compatible = !(nOrient == eOrient);
1139:       }
1140:       return compatible;
1141:     };
1142:     static void _replaceCell(const Obj<sieve_type>& sieve,
1143:                              const ALE::Mesh::point_type cell,
1144:                              std::map<int,int> *vertexRenumber,
1145:                              const int debug)
1146:     {
1147:       bool       replace = false;
1148:       PointArray newVertices;

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

1152:       for(sieve_type::traits::coneSequence::iterator v_iter = cCone->begin(); v_iter != cCone->end(); ++v_iter) {
1153:         if (vertexRenumber->find(*v_iter) != vertexRenumber->end()) {
1154:           if (debug) std::cout << "    vertex " << (*vertexRenumber)[*v_iter] << std::endl;
1155:           newVertices.insert(newVertices.end(), (*vertexRenumber)[*v_iter]);
1156:           replace = true;
1157:         } else {
1158:           if (debug) std::cout << "    vertex " << *v_iter << std::endl;
1159:           newVertices.insert(newVertices.end(), *v_iter);
1160:         } // if/else
1161:       } // for
1162:       if (replace) {
1163:         if (debug) std::cout << "  Replacing cell " << cell << std::endl;
1164:         sieve->clearCone(cell);
1165:         int color = 0;
1166:         for(PointArray::const_iterator v_iter = newVertices.begin(); v_iter != newVertices.end(); ++v_iter) {
1167:           sieve->addArrow(*v_iter, cell, color++);
1168:         } // for
1169:       }
1170:     };
1171:     template<class InputPoints>
1172:     static void _computeCensoredDepth(const Obj<ALE::Mesh>& mesh,
1173:                                       const Obj<ALE::Mesh::label_type>& depth,
1174:                                       const Obj<ALE::Mesh::sieve_type>& sieve,
1175:                                       const Obj<InputPoints>& points,
1176:                                       const ALE::Mesh::point_type& firstCohesiveCell,
1177:                                       const Obj<std::set<ALE::Mesh::point_type> >& modifiedPoints)
1178:     {
1179:       modifiedPoints->clear();

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

1187:         if(d1 != d0) {
1188:           mesh->setValue(depth, *p_iter, d1);
1189:           modifiedPoints->insert(*p_iter);
1190:         }
1191:       }
1192:       // FIX: We would like to avoid the copy here with support()
1193:       if(modifiedPoints->size() > 0) {
1194:         _computeCensoredDepth(mesh, depth, sieve, sieve->support(modifiedPoints), firstCohesiveCell, modifiedPoints);
1195:       }
1196:     };
1197:     static void create(const Obj<Mesh>& mesh, Obj<Mesh> fault, const Obj<Mesh::int_section_type>& groupField) {
1198:       static PetscLogEvent CreateFaultMesh_Event = 0, OrientFaultMesh_Event = 0, AddCohesivePoints_Event = 0, SplitMesh_Event = 0;

1200:       if (!CreateFaultMesh_Event) {
1201:         PetscLogEventRegister("CreateFaultMesh", 0,&CreateFaultMesh_Event);
1202:       }
1203:       if (!OrientFaultMesh_Event) {
1204:         PetscLogEventRegister("OrientFaultMesh", 0,&OrientFaultMesh_Event);
1205:       }
1206:       if (!AddCohesivePoints_Event) {
1207:         PetscLogEventRegister("AddCohesivePoints", 0,&AddCohesivePoints_Event);
1208:       }
1209:       if (!SplitMesh_Event) {
1210:         PetscLogEventRegister("SplitMesh", 0,&SplitMesh_Event);
1211:       }

1213:       const Obj<sieve_type>& sieve = mesh->getSieve();
1214:       const int  debug      = mesh->debug();
1215:       int        numCorners = 0;    // The number of vertices in a mesh cell
1216:       int        faceSize   = 0;    // The number of vertices in a mesh face
1217:       int       *indices    = NULL; // The indices of a face vertex set in a cell
1218:       PointArray origVertices;
1219:       PointArray faceVertices;
1220:       PointArray neighborVertices;
1221:       const bool constraintCell = false;

1223:       if (!mesh->commRank()) {
1224:         numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), mesh->depth())->size();
1225:         faceSize   = selection::numFaceVertices(mesh);
1226:         indices    = new int[faceSize];
1227:       }

1229:       //int f = sieve->base()->size() + sieve->cap()->size();
1230:       //ALE::Obj<PointSet> face = new PointSet();
1231: 
1232:       // Create a sieve which captures the fault
1233:       PetscLogEventBegin(CreateFaultMesh_Event,0,0,0,0);
1234:       fault = ALE::Selection<ALE::Mesh>::submesh(mesh, groupField);
1235:       if (debug) {fault->view("Fault mesh");}
1236:       PetscLogEventEnd(CreateFaultMesh_Event,0,0,0,0);
1237:       // Orient the fault sieve
1238:       PetscLogEventBegin(OrientFaultMesh_Event,0,0,0,0);
1239:       const Obj<sieve_type>&                faultSieve = fault->getSieve();
1240:       const ALE::Obj<Mesh::label_sequence>& fFaces     = fault->heightStratum(1);
1241:       int faultDepth      = fault->depth()-1; // Depth of fault cells
1242:       int numFaultCorners = 0; // The number of vertices in a fault cell

1244:       if (!fault->commRank()) {
1245:         numFaultCorners = faultSieve->nCone(*fFaces->begin(), faultDepth)->size();
1246:         if (debug) std::cout << "  Fault corners " << numFaultCorners << std::endl;
1247:         assert(numFaultCorners == faceSize);
1248:       }
1249:       PetscLogEventEnd(OrientFaultMesh_Event,0,0,0,0);

1251:       // Add new shadow vertices and possibly Lagrange multipler vertices
1252:       PetscLogEventBegin(AddCohesivePoints_Event,0,0,0,0);
1253:       const ALE::Obj<Mesh::label_sequence>&   fVertices  = fault->depthStratum(0);
1254:       const ALE::Obj<std::set<std::string> >& groupNames = mesh->getIntSections();
1255:       Mesh::point_type newPoint = sieve->base()->size() + sieve->cap()->size();
1256:       std::map<int,int> vertexRenumber;
1257: 
1258:       for(Mesh::label_sequence::iterator v_iter = fVertices->begin(); v_iter != fVertices->end(); ++v_iter, ++newPoint) {
1259:         vertexRenumber[*v_iter] = newPoint;
1260:         if (debug) {std::cout << "Duplicating " << *v_iter << " to " << vertexRenumber[*v_iter] << std::endl;}

1262:         // Add shadow and constraint vertices (if they exist) to group
1263:         // associated with fault
1264:         groupField->addPoint(newPoint, 1);
1265:         if (constraintCell) groupField->addPoint(newPoint+1, 1);

1267:         // Add shadow vertices to other groups, don't add constraint
1268:         // vertices (if they exist) because we don't want BC, etc to act
1269:         // on constraint vertices
1270:         for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
1271:           const ALE::Obj<int_section_type>& group = mesh->getIntSection(*name);
1272:           if (group->hasPoint(*v_iter)) group->addPoint(newPoint, 1);
1273:         }
1274:         if (constraintCell) newPoint++;
1275:       }
1276:       for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
1277:         mesh->reallocate(mesh->getIntSection(*name));
1278:       }

1280:       // Split the mesh along the fault sieve and create cohesive elements
1281:       const Obj<ALE::Mesh::label_sequence>&     faces       = fault->depthStratum(1);
1282:       const Obj<ALE::Mesh::arrow_section_type>& orientation = mesh->getArrowSection("orientation");
1283:       int firstCohesiveCell = newPoint;
1284:       PointSet replaceCells;
1285:       PointSet noReplaceCells;

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

1294:         if (debug) std::cout << "  Checking orientation against cell " << cell << std::endl;
1295:         if (numFaultCorners == 2) reversed = orientation->restrictPoint(arrow)[0] == -2;
1296:         if (reversed) {
1297:           replaceCells.insert(cell);
1298:           noReplaceCells.insert(*(++cells->begin()));
1299:         } else {
1300:           replaceCells.insert(*(++cells->begin()));
1301:           noReplaceCells.insert(cell);
1302:         }
1303:         //selection::getOrientedFace(mesh, cell, &vertexRenumber, numCorners, indices, &origVertices, &faceVertices);
1304:         //const Obj<sieve_type::coneArray> faceCone = faultSieve->nCone(*f_iter, faultDepth);

1306:         // Adding cohesive cell (not interpolated)
1307:         const Obj<sieve_type::coneArray>&     fCone  = faultSieve->nCone(*f_iter, faultDepth);
1308:         const sieve_type::coneArray::iterator fBegin = fCone->begin();
1309:         const sieve_type::coneArray::iterator fEnd   = fCone->end();
1310:         int color = 0;

1312:         if (debug) {std::cout << "  Creating cohesive cell " << newPoint << std::endl;}
1313:         for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
1314:           if (debug) {std::cout << "    vertex " << *v_iter << std::endl;}
1315:           sieve->addArrow(*v_iter, newPoint, color++);
1316:         }
1317:         for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
1318:           if (debug) {std::cout << "    shadow vertex " << vertexRenumber[*v_iter] << std::endl;}
1319:           sieve->addArrow(vertexRenumber[*v_iter], newPoint, color++);
1320:         }
1321:       }
1322:       PetscLogEventEnd(AddCohesivePoints_Event,0,0,0,0);
1323:       // Replace all cells with a vertex on the fault that share a face with this one, or with one that does
1324:       PetscLogEventBegin(SplitMesh_Event,0,0,0,0);
1325:       const int_section_type::chart_type&          chart    = groupField->getChart();
1326:       const int_section_type::chart_type::iterator chartEnd = chart.end();

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

1331:         while(modified) {
1332:           modified = false;
1333:           const Obj<sieve_type::traits::supportSequence>&     neighbors = sieve->support(*v_iter);
1334:           const sieve_type::traits::supportSequence::iterator end       = neighbors->end();

1336:           for(sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != end; ++n_iter) {
1337:             if (replaceCells.find(*n_iter)   != replaceCells.end())   continue;
1338:             if (noReplaceCells.find(*n_iter) != noReplaceCells.end()) continue;
1339:             if (*n_iter >= firstCohesiveCell) continue;
1340:             if (debug) std::cout << "  Checking fault neighbor " << *n_iter << std::endl;
1341:             // If neighbors shares a faces with anyone in replaceCells, then add
1342:             for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
1343:               const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, mesh->depth());

1345:               if ((int) preFace->size() == faceSize) {
1346:                 if (debug) std::cout << "    Scheduling " << *n_iter << " for replacement" << std::endl;
1347:                 replaceCells.insert(*n_iter);
1348:                 modified = true;
1349:                 break;
1350:               }
1351:             }
1352:           }
1353:         }
1354:       }
1355:       for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
1356:         _replaceCell(sieve, *c_iter, &vertexRenumber, debug);
1357:       }
1358:       if (!fault->commRank()) delete [] indices;
1359:       mesh->stratify();
1360:       const ALE::Obj<Mesh::label_type>& label          = mesh->createLabel(std::string("censored depth"));
1361:       const ALE::Obj<PointSet>          modifiedPoints = new PointSet();
1362:       _computeCensoredDepth(mesh, label, mesh->getSieve(), mesh->getSieve()->roots(), firstCohesiveCell, modifiedPoints);
1363:       if (debug) mesh->view("Mesh with Cohesive Elements");

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

1369:       for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
1370:         coordinates->addPoint(vertexRenumber[*v_iter], coordinates->getFiberDimension(*v_iter));
1371:         if (constraintCell) {
1372:           coordinates->addPoint(vertexRenumber[*v_iter]+1, coordinates->getFiberDimension(*v_iter));
1373:         }
1374:       }
1375:       mesh->reallocate(coordinates);
1376:       for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
1377:         coordinates->updatePoint(vertexRenumber[*v_iter], coordinates->restrictPoint(*v_iter));
1378:         if (constraintCell) {
1379:         coordinates->updatePoint(vertexRenumber[*v_iter]+1, coordinates->restrictPoint(*v_iter));
1380:         }
1381:       }
1382:       if (debug) coordinates->view("Coordinates with shadow vertices");
1383:       PetscLogEventEnd(SplitMesh_Event,0,0,0,0);
1384:     };
1385:   };
1386: };
1387: #endif

1389: #endif