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