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 interpolates 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->view("Submesh before stratify");
743: submesh->stratify();
744: if (debug) submesh->view("Submesh");
746: Obj<output_mesh_type> isubmesh = new output_mesh_type(submesh->comm(), submesh->getDimension(), submesh->debug());
747: Obj<typename output_mesh_type::sieve_type> isieve = new typename output_mesh_type::sieve_type(submesh->comm(), submesh->debug());
748: std::map<typename output_mesh_type::point_type,typename output_mesh_type::point_type> renumbering;
749: isubmesh->setSieve(isieve);
750: ALE::ISieveConverter::convertMesh(*submesh, *isubmesh, renumbering, false);
751: return isubmesh;
752: };
753: public:
754: // This takes in a section and creates a submesh from the vertices in the section chart
755: // This is a hyperplane of one dimension lower than the mesh
756: static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
757: const int dim = mesh->getDimension();
758: const int depth = mesh->depth();
760: if (dim == depth) {
761: return submesh_interpolated(mesh, label, dimension, false);
762: } else if (depth == 1) {
763: return submesh_uninterpolated(mesh, label, dimension);
764: }
765: throw ALE::Exception("Cannot handle partially interpolated meshes");
766: };
767: template<typename output_mesh_type>
768: static Obj<output_mesh_type> submeshV(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
769: const int dim = mesh->getDimension();
770: const int depth = mesh->depth();
772: #if 0
773: if (dim == depth) {
774: //return submesh_interpolated(mesh, label, dimension, false);
775: throw ALE::Exception("Cannot handle interpolated meshes");
776: } else if (depth == 1) {
777: return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
778: }
779: #else
780: if (depth == 1) {
781: return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
782: } else if (dim == depth) {
783: //return submesh_interpolated(mesh, label, dimension, false);
784: throw ALE::Exception("Cannot handle interpolated meshes");
785: }
786: #endif
787: throw ALE::Exception("Cannot handle partially interpolated meshes");
788: };
789: // This creates a submesh consisting of the union of the closures of the given points
790: // This mesh is the same dimension as in the input mesh
791: template<typename Points>
792: static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<Points>& points, const int dim = -1) {
793: const Obj<sieve_type>& sieve = mesh->getSieve();
794: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), dim >= 0 ? dim : mesh->getDimension(), mesh->debug());
795: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
796: Obj<PointSet> newPoints = new PointSet();
797: Obj<PointSet> modPoints = new PointSet();
798: Obj<PointSet> tmpPoints;
800: newMesh->setSieve(newSieve);
801: for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
802: newPoints->insert(*p_iter);
803: do {
804: modPoints->clear();
805: for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
806: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*np_iter);
807: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
808: int c = 0;
810: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter, ++c) {
811: newSieve->addArrow(*c_iter, *np_iter, c);
812: }
813: modPoints->insert(cone->begin(), cone->end());
814: }
815: tmpPoints = newPoints;
816: newPoints = modPoints;
817: modPoints = tmpPoints;
818: } while(newPoints->size());
819: newPoints->insert(*p_iter);
820: do {
821: modPoints->clear();
822: for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
823: const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*np_iter);
824: const typename sieve_type::traits::supportSequence::iterator end = support->end();
825: int s = 0;
827: for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != end; ++s_iter, ++s) {
828: newSieve->addArrow(*np_iter, *s_iter, s);
829: }
830: modPoints->insert(support->begin(), support->end());
831: }
832: tmpPoints = newPoints;
833: newPoints = modPoints;
834: modPoints = tmpPoints;
835: } while(newPoints->size());
836: }
837: newMesh->stratify();
838: newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
839: newMesh->setArrowSection("orientation", mesh->getArrowSection("orientation"));
840: return newMesh;
841: };
842: protected:
843: static Obj<mesh_type> boundary_uninterpolated(const Obj<mesh_type>& mesh) {
844: throw ALE::Exception("Not yet implemented");
845: const Obj<typename mesh_type::label_sequence>& cells = mesh->heightStratum(0);
846: const Obj<sieve_type>& sieve = mesh->getSieve();
847: const typename mesh_type::label_sequence::iterator cBegin = cells->begin();
848: const typename mesh_type::label_sequence::iterator cEnd = cells->end();
849: const int faceSize = numFaceVertices(mesh);
851: for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
852: const Obj<typename sieve_type::traits::coneSequence>& vertices = sieve->cone(*c_iter);
853: const typename sieve_type::traits::coneSequence::iterator vBegin = vertices->begin();
854: const typename sieve_type::traits::coneSequence::iterator vEnd = vertices->end();
855: //PointArray cell(vertices->begin(), vertices->end());
857: // For each vertex, gather
858: for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
859: const Obj<typename sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
860: const typename sieve_type::traits::supportSequence::iterator nBegin = neighbors->begin();
861: const typename sieve_type::traits::supportSequence::iterator nEnd = neighbors->end();
863: for(typename sieve_type::traits::supportSequence::iterator n_iter = nBegin; n_iter != nEnd; ++n_iter) {
864: const Obj<typename sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, 1);
866: if (preFace->size() == faceSize) {
867: }
868: }
869: }
870: // For each face
871: // - determine if its legal
872:
873: // - determine if its part of a neighboring cell
874: // - if not, its a boundary face
875: //subsets(cell, faceSize, inserter);
876: }
877: };
878: static void addClosure(const Obj<sieve_type>& sieveA, const Obj<sieve_type>& sieveB, const point_type& p, const int depth = 1) {
879: Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
880: Obj<typename sieve_type::coneSet> next = new typename sieve_type::coneSet();
881: Obj<typename sieve_type::coneSet> tmp;
883: current->insert(p);
884: while(current->size()) {
885: for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
886: const Obj<typename sieve_type::traits::coneSequence>& cone = sieveA->cone(*p_iter);
887: const typename sieve_type::traits::coneSequence::iterator begin = cone->begin();
888: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
890: for(typename sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
891: sieveB->addArrow(*c_iter, *p_iter, c_iter.color());
892: next->insert(*c_iter);
893: }
894: }
895: tmp = current; current = next; next = tmp;
896: next->clear();
897: }
898: if (!depth) {
899: const Obj<typename sieve_type::traits::supportSequence>& support = sieveA->support(p);
900: const typename sieve_type::traits::supportSequence::iterator begin = support->begin();
901: const typename sieve_type::traits::supportSequence::iterator end = support->end();
902:
903: for(typename sieve_type::traits::supportSequence::iterator s_iter = begin; s_iter != end; ++s_iter) {
904: sieveB->addArrow(p, *s_iter, s_iter.color());
905: next->insert(*s_iter);
906: }
907: }
908: };
909: static Obj<mesh_type> boundary_interpolated(const Obj<mesh_type>& mesh, const int faceHeight = 1) {
910: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
911: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
912: const Obj<sieve_type>& sieve = mesh->getSieve();
913: const Obj<typename mesh_type::label_sequence>& faces = mesh->heightStratum(faceHeight);
914: const typename mesh_type::label_sequence::iterator fBegin = faces->begin();
915: const typename mesh_type::label_sequence::iterator fEnd = faces->end();
916: const int depth = faceHeight - mesh->depth();
918: for(typename mesh_type::label_sequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
919: const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*f_iter);
921: if (support->size() == 1) {
922: addClosure(sieve, newSieve, *f_iter, depth);
923: }
924: }
925: newMesh->setSieve(newSieve);
926: newMesh->stratify();
927: return newMesh;
928: };
929: template<typename SieveType>
930: static void addClosureV(const Obj<SieveType>& sieveA, const Obj<SieveType>& sieveB, const point_type& p, const int depth = 1) {
931: typedef std::set<typename SieveType::point_type> coneSet;
932: ALE::ISieveVisitor::PointRetriever<SieveType> cV(std::max(1, sieveA->getMaxConeSize()));
933: Obj<coneSet> current = new coneSet();
934: Obj<coneSet> next = new coneSet();
935: Obj<coneSet> tmp;
937: current->insert(p);
938: while(current->size()) {
939: for(typename coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
940: sieveA->cone(*p_iter, cV);
941: const Mesh::point_type *cone = cV.getPoints();
943: for(int c = 0; c < cV.getSize(); ++c) {
944: sieveB->addArrow(cone[c], *p_iter);
945: next->insert(cone[c]);
946: }
947: }
948: tmp = current; current = next; next = tmp;
949: next->clear();
950: }
951: if (!depth) {
952: ALE::ISieveVisitor::PointRetriever<SieveType> sV(std::max(1, sieveA->getMaxSupportSize()));
954: sieveA->support(p, sV);
955: const typename SieveType::point_type *support = sV.getPoints();
956:
957: for(int s = 0; s < sV.getSize(); ++s) {
958: sieveB->addArrow(p, support[s]);
959: }
960: }
961: };
962: public:
963: static Obj<mesh_type> boundary(const Obj<mesh_type>& mesh) {
964: const int dim = mesh->getDimension();
965: const int depth = mesh->depth();
967: if (dim == depth) {
968: return boundary_interpolated(mesh);
969: } else if (depth == dim+1) {
970: return boundary_interpolated(mesh, 2);
971: } else if (depth == 1) {
972: return boundary_uninterpolated(mesh);
973: } else if (depth == -1) {
974: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
975: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
977: newMesh->setSieve(newSieve);
978: newMesh->stratify();
979: return newMesh;
980: }
981: throw ALE::Exception("Cannot handle partially interpolated meshes");
982: };
983: template<typename MeshTypeQ>
984: static Obj<MeshTypeQ> boundaryV(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
985: Obj<MeshTypeQ> newMesh = new MeshTypeQ(mesh->comm(), mesh->getDimension()-1, mesh->debug());
986: Obj<typename MeshTypeQ::sieve_type> newSieve = new typename MeshTypeQ::sieve_type(mesh->comm(), mesh->debug());
987: const Obj<typename MeshTypeQ::sieve_type>& sieve = mesh->getSieve();
988: const Obj<typename MeshTypeQ::label_sequence>& faces = mesh->heightStratum(faceHeight);
989: const typename MeshTypeQ::label_sequence::iterator fBegin = faces->begin();
990: const typename MeshTypeQ::label_sequence::iterator fEnd = faces->end();
991: const int depth = faceHeight - mesh->depth();
992: ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
994: for(typename MeshTypeQ::label_sequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
995: sieve->support(*f_iter, sV);
997: if (sV.getSize() == 1) {
998: addClosureV(sieve, newSieve, *f_iter, depth);
999: }
1000: sV.clear();
1001: }
1002: newMesh->setSieve(newSieve);
1003: newMesh->stratify();
1004: return newMesh;
1005: };
1006: public:
1007: static Obj<mesh_type> interpolateMesh(const Obj<mesh_type>& mesh) {
1008: const Obj<sieve_type> sieve = mesh->getSieve();
1009: const int dim = mesh->getDimension();
1010: const int numVertices = mesh->depthStratum(0)->size();
1011: const Obj<typename mesh_type::label_sequence>& cells = mesh->heightStratum(0);
1012: const int numCells = cells->size();
1013: const int corners = sieve->cone(*cells->begin())->size();
1014: const int firstVertex = numCells;
1015: const int debug = sieve->debug();
1016: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), dim, mesh->debug());
1017: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
1018: const Obj<typename mesh_type::arrow_section_type>& orientation = newMesh->getArrowSection("orientation");
1020: newMesh->setSieve(newSieve);
1021: // Create a map from dimension to the current element number for that dimension
1022: std::map<int,point_type*> curElement;
1023: std::map<int,PointArray> bdVertices;
1024: std::map<int,PointArray> faces;
1025: std::map<int,oPointArray> oFaces;
1026: int curCell = 0;
1027: int curVertex = firstVertex;
1028: int newElement = firstVertex+numVertices;
1029: int o;
1031: curElement[0] = &curVertex;
1032: curElement[dim] = &curCell;
1033: for(int d = 1; d < dim; d++) {
1034: curElement[d] = &newElement;
1035: }
1036: typename mesh_type::label_sequence::iterator cBegin = cells->begin();
1037: typename mesh_type::label_sequence::iterator cEnd = cells->end();
1039: for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
1040: typename sieve_type::point_type cell = *c_iter;
1041: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
1042: const typename sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
1043: const typename sieve_type::traits::coneSequence::iterator vEnd = cone->end();
1045: // Build the cell
1046: bdVertices[dim].clear();
1047: for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
1048: typename sieve_type::point_type vertex(*v_iter);
1050: if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
1051: bdVertices[dim].push_back(vertex);
1052: }
1053: if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}
1055: if (corners != dim+1) {
1056: ALE::SieveBuilder<mesh_type>::buildHexFaces(newSieve, dim, curElement, bdVertices, faces, cell);
1057: } else {
1058: ALE::SieveBuilder<mesh_type>::buildFaces(newSieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
1059: }
1060: }
1061: newMesh->stratify();
1062: newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
1063: return newMesh;
1064: };
1065: };
1068: }
1070: #if 0
1071: namespace ALE {
1072: class MySelection {
1073: public:
1074: typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
1075: typedef ALE::Selection<ALE::Mesh> selection;
1076: typedef ALE::Mesh::sieve_type sieve_type;
1077: typedef ALE::Mesh::int_section_type int_section_type;
1078: typedef ALE::Mesh::real_section_type real_section_type;
1079: typedef std::set<ALE::Mesh::point_type> PointSet;
1080: typedef std::vector<ALE::Mesh::point_type> PointArray;
1081: public:
1082: template<class InputPoints>
1083: static bool _compatibleOrientation(const Obj<Mesh>& mesh,
1084: const ALE::Mesh::point_type& p,
1085: const ALE::Mesh::point_type& q,
1086: const int numFaultCorners,
1087: const int faultFaceSize,
1088: const int faultDepth,
1089: const Obj<InputPoints>& points,
1090: int indices[],
1091: PointArray *origVertices,
1092: PointArray *faceVertices,
1093: PointArray *neighborVertices)
1094: {
1095: typedef ALE::Selection<ALE::Mesh> selection;
1096: const int debug = mesh->debug();
1097: bool compatible;
1099: bool eOrient = selection::getOrientedFace(mesh, p, points, numFaultCorners, indices, origVertices, faceVertices);
1100: bool nOrient = selection::getOrientedFace(mesh, q, points, numFaultCorners, indices, origVertices, neighborVertices);
1102: if (faultFaceSize > 1) {
1103: if (debug) {
1104: for(PointArray::iterator v_iter = faceVertices->begin(); v_iter != faceVertices->end(); ++v_iter) {
1105: std::cout << " face vertex " << *v_iter << std::endl;
1106: }
1107: for(PointArray::iterator v_iter = neighborVertices->begin(); v_iter != neighborVertices->end(); ++v_iter) {
1108: std::cout << " neighbor vertex " << *v_iter << std::endl;
1109: }
1110: }
1111: compatible = !(*faceVertices->begin() == *neighborVertices->begin());
1112: } else {
1113: compatible = !(nOrient == eOrient);
1114: }
1115: return compatible;
1116: };
1117: static void _replaceCell(const Obj<sieve_type>& sieve,
1118: const ALE::Mesh::point_type cell,
1119: std::map<int,int> *vertexRenumber,
1120: const int debug)
1121: {
1122: bool replace = false;
1123: PointArray newVertices;
1125: const Obj<sieve_type::traits::coneSequence>& cCone = sieve->cone(cell);
1127: for(sieve_type::traits::coneSequence::iterator v_iter = cCone->begin(); v_iter != cCone->end(); ++v_iter) {
1128: if (vertexRenumber->find(*v_iter) != vertexRenumber->end()) {
1129: if (debug) std::cout << " vertex " << (*vertexRenumber)[*v_iter] << std::endl;
1130: newVertices.insert(newVertices.end(), (*vertexRenumber)[*v_iter]);
1131: replace = true;
1132: } else {
1133: if (debug) std::cout << " vertex " << *v_iter << std::endl;
1134: newVertices.insert(newVertices.end(), *v_iter);
1135: } // if/else
1136: } // for
1137: if (replace) {
1138: if (debug) std::cout << " Replacing cell " << cell << std::endl;
1139: sieve->clearCone(cell);
1140: int color = 0;
1141: for(PointArray::const_iterator v_iter = newVertices.begin(); v_iter != newVertices.end(); ++v_iter) {
1142: sieve->addArrow(*v_iter, cell, color++);
1143: } // for
1144: }
1145: };
1146: template<class InputPoints>
1147: static void _computeCensoredDepth(const Obj<ALE::Mesh>& mesh,
1148: const Obj<ALE::Mesh::label_type>& depth,
1149: const Obj<ALE::Mesh::sieve_type>& sieve,
1150: const Obj<InputPoints>& points,
1151: const ALE::Mesh::point_type& firstCohesiveCell,
1152: const Obj<std::set<ALE::Mesh::point_type> >& modifiedPoints)
1153: {
1154: modifiedPoints->clear();
1156: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1157: if (*p_iter >= firstCohesiveCell) continue;
1158: // Compute the max depth of the points in the cone of p, and add 1
1159: int d0 = mesh->getValue(depth, *p_iter, -1);
1160: int d1 = mesh->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;
1162: if(d1 != d0) {
1163: mesh->setValue(depth, *p_iter, d1);
1164: modifiedPoints->insert(*p_iter);
1165: }
1166: }
1167: // FIX: We would like to avoid the copy here with support()
1168: if(modifiedPoints->size() > 0) {
1169: _computeCensoredDepth(mesh, depth, sieve, sieve->support(modifiedPoints), firstCohesiveCell, modifiedPoints);
1170: }
1171: };
1172: static void create(const Obj<Mesh>& mesh, Obj<Mesh> fault, const Obj<Mesh::int_section_type>& groupField) {
1173: static PetscLogEvent CreateFaultMesh_Event = 0, OrientFaultMesh_Event = 0, AddCohesivePoints_Event = 0, SplitMesh_Event = 0;
1175: if (!CreateFaultMesh_Event) {
1176: PetscLogEventRegister("CreateFaultMesh", 0,&CreateFaultMesh_Event);
1177: }
1178: if (!OrientFaultMesh_Event) {
1179: PetscLogEventRegister("OrientFaultMesh", 0,&OrientFaultMesh_Event);
1180: }
1181: if (!AddCohesivePoints_Event) {
1182: PetscLogEventRegister("AddCohesivePoints", 0,&AddCohesivePoints_Event);
1183: }
1184: if (!SplitMesh_Event) {
1185: PetscLogEventRegister("SplitMesh", 0,&SplitMesh_Event);
1186: }
1188: const Obj<sieve_type>& sieve = mesh->getSieve();
1189: const int debug = mesh->debug();
1190: int numCorners = 0; // The number of vertices in a mesh cell
1191: int faceSize = 0; // The number of vertices in a mesh face
1192: int *indices = NULL; // The indices of a face vertex set in a cell
1193: PointArray origVertices;
1194: PointArray faceVertices;
1195: PointArray neighborVertices;
1196: const bool constraintCell = false;
1198: if (!mesh->commRank()) {
1199: numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), mesh->depth())->size();
1200: faceSize = selection::numFaceVertices(mesh);
1201: indices = new int[faceSize];
1202: }
1204: //int f = sieve->base()->size() + sieve->cap()->size();
1205: //ALE::Obj<PointSet> face = new PointSet();
1206:
1207: // Create a sieve which captures the fault
1208: PetscLogEventBegin(CreateFaultMesh_Event,0,0,0,0);
1209: fault = ALE::Selection<ALE::Mesh>::submesh(mesh, groupField);
1210: if (debug) {fault->view("Fault mesh");}
1211: PetscLogEventEnd(CreateFaultMesh_Event,0,0,0,0);
1212: // Orient the fault sieve
1213: PetscLogEventBegin(OrientFaultMesh_Event,0,0,0,0);
1214: const Obj<sieve_type>& faultSieve = fault->getSieve();
1215: const ALE::Obj<Mesh::label_sequence>& fFaces = fault->heightStratum(1);
1216: int faultDepth = fault->depth()-1; // Depth of fault cells
1217: int numFaultCorners = 0; // The number of vertices in a fault cell
1219: if (!fault->commRank()) {
1220: numFaultCorners = faultSieve->nCone(*fFaces->begin(), faultDepth)->size();
1221: if (debug) std::cout << " Fault corners " << numFaultCorners << std::endl;
1222: assert(numFaultCorners == faceSize);
1223: }
1224: PetscLogEventEnd(OrientFaultMesh_Event,0,0,0,0);
1226: // Add new shadow vertices and possibly Lagrange multipler vertices
1227: PetscLogEventBegin(AddCohesivePoints_Event,0,0,0,0);
1228: const ALE::Obj<Mesh::label_sequence>& fVertices = fault->depthStratum(0);
1229: const ALE::Obj<std::set<std::string> >& groupNames = mesh->getIntSections();
1230: Mesh::point_type newPoint = sieve->base()->size() + sieve->cap()->size();
1231: std::map<int,int> vertexRenumber;
1232:
1233: for(Mesh::label_sequence::iterator v_iter = fVertices->begin(); v_iter != fVertices->end(); ++v_iter, ++newPoint) {
1234: vertexRenumber[*v_iter] = newPoint;
1235: if (debug) {std::cout << "Duplicating " << *v_iter << " to " << vertexRenumber[*v_iter] << std::endl;}
1237: // Add shadow and constraint vertices (if they exist) to group
1238: // associated with fault
1239: groupField->addPoint(newPoint, 1);
1240: if (constraintCell) groupField->addPoint(newPoint+1, 1);
1242: // Add shadow vertices to other groups, don't add constraint
1243: // vertices (if they exist) because we don't want BC, etc to act
1244: // on constraint vertices
1245: for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
1246: const ALE::Obj<int_section_type>& group = mesh->getIntSection(*name);
1247: if (group->hasPoint(*v_iter)) group->addPoint(newPoint, 1);
1248: }
1249: if (constraintCell) newPoint++;
1250: }
1251: for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
1252: mesh->reallocate(mesh->getIntSection(*name));
1253: }
1255: // Split the mesh along the fault sieve and create cohesive elements
1256: const Obj<ALE::Mesh::label_sequence>& faces = fault->depthStratum(1);
1257: const Obj<ALE::Mesh::arrow_section_type>& orientation = mesh->getArrowSection("orientation");
1258: int firstCohesiveCell = newPoint;
1259: PointSet replaceCells;
1260: PointSet noReplaceCells;
1262: for(ALE::Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter, ++newPoint) {
1263: if (debug) std::cout << "Considering fault face " << *f_iter << std::endl;
1264: const ALE::Obj<sieve_type::traits::supportSequence>& cells = faultSieve->support(*f_iter);
1265: const ALE::Mesh::arrow_section_type::point_type arrow(*cells->begin(), *f_iter);
1266: bool reversed = orientation->restrictPoint(arrow)[0] < 0;
1267: const ALE::Mesh::point_type cell = *cells->begin();
1269: if (debug) std::cout << " Checking orientation against cell " << cell << std::endl;
1270: if (numFaultCorners == 2) reversed = orientation->restrictPoint(arrow)[0] == -2;
1271: if (reversed) {
1272: replaceCells.insert(cell);
1273: noReplaceCells.insert(*(++cells->begin()));
1274: } else {
1275: replaceCells.insert(*(++cells->begin()));
1276: noReplaceCells.insert(cell);
1277: }
1278: //selection::getOrientedFace(mesh, cell, &vertexRenumber, numCorners, indices, &origVertices, &faceVertices);
1279: //const Obj<sieve_type::coneArray> faceCone = faultSieve->nCone(*f_iter, faultDepth);
1281: // Adding cohesive cell (not interpolated)
1282: const Obj<sieve_type::coneArray>& fCone = faultSieve->nCone(*f_iter, faultDepth);
1283: const sieve_type::coneArray::iterator fBegin = fCone->begin();
1284: const sieve_type::coneArray::iterator fEnd = fCone->end();
1285: int color = 0;
1287: if (debug) {std::cout << " Creating cohesive cell " << newPoint << std::endl;}
1288: for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
1289: if (debug) {std::cout << " vertex " << *v_iter << std::endl;}
1290: sieve->addArrow(*v_iter, newPoint, color++);
1291: }
1292: for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
1293: if (debug) {std::cout << " shadow vertex " << vertexRenumber[*v_iter] << std::endl;}
1294: sieve->addArrow(vertexRenumber[*v_iter], newPoint, color++);
1295: }
1296: }
1297: PetscLogEventEnd(AddCohesivePoints_Event,0,0,0,0);
1298: // Replace all cells with a vertex on the fault that share a face with this one, or with one that does
1299: PetscLogEventBegin(SplitMesh_Event,0,0,0,0);
1300: const int_section_type::chart_type& chart = groupField->getChart();
1301: const int_section_type::chart_type::iterator chartEnd = chart.end();
1303: for(PointSet::const_iterator v_iter = chart.begin(); v_iter != chartEnd; ++v_iter) {
1304: bool modified = true;
1306: while(modified) {
1307: modified = false;
1308: const Obj<sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
1309: const sieve_type::traits::supportSequence::iterator end = neighbors->end();
1311: for(sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != end; ++n_iter) {
1312: if (replaceCells.find(*n_iter) != replaceCells.end()) continue;
1313: if (noReplaceCells.find(*n_iter) != noReplaceCells.end()) continue;
1314: if (*n_iter >= firstCohesiveCell) continue;
1315: if (debug) std::cout << " Checking fault neighbor " << *n_iter << std::endl;
1316: // If neighbors shares a faces with anyone in replaceCells, then add
1317: for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
1318: const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, mesh->depth());
1320: if ((int) preFace->size() == faceSize) {
1321: if (debug) std::cout << " Scheduling " << *n_iter << " for replacement" << std::endl;
1322: replaceCells.insert(*n_iter);
1323: modified = true;
1324: break;
1325: }
1326: }
1327: }
1328: }
1329: }
1330: for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
1331: _replaceCell(sieve, *c_iter, &vertexRenumber, debug);
1332: }
1333: if (!fault->commRank()) delete [] indices;
1334: mesh->stratify();
1335: const ALE::Obj<Mesh::label_type>& label = mesh->createLabel(std::string("censored depth"));
1336: const ALE::Obj<PointSet> modifiedPoints = new PointSet();
1337: _computeCensoredDepth(mesh, label, mesh->getSieve(), mesh->getSieve()->roots(), firstCohesiveCell, modifiedPoints);
1338: if (debug) mesh->view("Mesh with Cohesive Elements");
1340: // Fix coordinates
1341: const Obj<real_section_type>& coordinates = mesh->getRealSection("coordinates");
1342: const Obj<ALE::Mesh::label_sequence>& fVertices2 = fault->depthStratum(0);
1344: for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
1345: coordinates->addPoint(vertexRenumber[*v_iter], coordinates->getFiberDimension(*v_iter));
1346: if (constraintCell) {
1347: coordinates->addPoint(vertexRenumber[*v_iter]+1, coordinates->getFiberDimension(*v_iter));
1348: }
1349: }
1350: mesh->reallocate(coordinates);
1351: for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
1352: coordinates->updatePoint(vertexRenumber[*v_iter], coordinates->restrictPoint(*v_iter));
1353: if (constraintCell) {
1354: coordinates->updatePoint(vertexRenumber[*v_iter]+1, coordinates->restrictPoint(*v_iter));
1355: }
1356: }
1357: if (debug) coordinates->view("Coordinates with shadow vertices");
1358: PetscLogEventEnd(SplitMesh_Event,0,0,0,0);
1359: };
1360: };
1361: };
1362: #endif
1364: #endif