Actual source code: Selection.hh
petsc-3.3-p7 2013-05-11
1: #ifndef included_ALE_Selection_hh
2: #define included_ALE_Selection_hh
4: #ifndef included_ALE_SieveAlgorithms_hh
5: #include <sieve/SieveAlgorithms.hh>
6: #endif
8: #ifndef included_ALE_SieveBuilder_hh
9: #include <sieve/SieveBuilder.hh>
10: #endif
12: #ifndef included_ALE_Mesh_hh
13: #include <sieve/Mesh.hh>
14: #endif
16: namespace ALE {
17: template<typename Mesh_>
18: class Selection {
19: public:
20: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
21: typedef Mesh_ mesh_type;
22: typedef typename mesh_type::sieve_type sieve_type;
23: typedef typename mesh_type::point_type point_type;
24: typedef typename mesh_type::int_section_type int_section_type;
25: typedef std::set<point_type> PointSet;
26: typedef std::vector<point_type> PointArray;
27: typedef std::pair<typename sieve_type::point_type, int> oPoint_type;
28: typedef std::vector<oPoint_type> oPointArray;
29: typedef typename ALE::SieveAlg<mesh_type> sieveAlg;
30: protected:
31: template<typename Sieve, typename FaceType>
32: class FaceVisitor {
33: public:
34: typedef typename Sieve::point_type point_type;
35: typedef typename Sieve::arrow_type arrow_type;
36: protected:
37: const FaceType& face;
38: PointArray& origVertices;
39: PointArray& faceVertices;
40: int *indices;
41: const int debug;
42: int oppositeVertex;
43: int v;
44: public:
45: FaceVisitor(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) {};
46: void visitPoint(const point_type& point) {
47: if (face->find(point) != face->end()) {
48: if (debug) std::cout << " vertex " << point << std::endl;
49: indices[origVertices.size()] = v;
50: origVertices.insert(origVertices.end(), point);
51: } else {
52: if (debug) std::cout << " vertex " << point << std::endl;
53: oppositeVertex = v;
54: }
55: ++v;
56: };
57: void visitArrow(const arrow_type&) {};
58: public:
59: int getOppositeVertex() {return this->oppositeVertex;};
60: };
61: public:
62: template<typename Processor>
63: static void subsets(const PointArray& v, const int size, Processor& processor, Obj<PointArray> *out = NULL, const int min = 0) {
64: if (size == 0) {
65: processor(*out);
66: return;
67: }
68: if (min == 0) {
69: out = new Obj<PointArray>();
70: *out = new PointArray();
71: }
72: for(int i = min; i < (int) v.size(); ++i) {
73: (*out)->push_back(v[i]);
74: subsets(v, size-1, processor, out, i+1);
75: (*out)->pop_back();
76: }
77: if (min == 0) {delete out;}
78: };
79: template<typename Mesh>
80: static int numFaceVertices(const Obj<Mesh>& mesh) {
81: return numFaceVertices(mesh, mesh->getNumCellCorners());
82: };
83: template<typename Mesh>
84: static int numFaceVertices(const Obj<Mesh>& mesh, const unsigned int numCorners) {
85: //unsigned int numCorners = mesh->getNumCellCorners(cell, depth);
86: const int cellDim = mesh->getDimension();
87: unsigned int _numFaceVertices = 0;
89: switch (cellDim) {
90: case 0 :
91: _numFaceVertices = 0;
92: break;
93: case 1 :
94: _numFaceVertices = 1;
95: break;
96: case 2:
97: switch (numCorners) {
98: case 3 : // triangle
99: _numFaceVertices = 2; // Edge has 2 vertices
100: break;
101: case 4 : // quadrilateral
102: _numFaceVertices = 2; // Edge has 2 vertices
103: break;
104: case 6 : // quadratic triangle
105: _numFaceVertices = 3; // Edge has 3 vertices
106: break;
107: case 9 : // quadratic quadrilateral
108: _numFaceVertices = 3; // Edge has 3 vertices
109: break;
110: default :
111: throw ALE::Exception("Invalid number of face corners");
112: }
113: break;
114: case 3:
115: switch (numCorners) {
116: case 4 : // tetradehdron
117: _numFaceVertices = 3; // Face has 3 vertices
118: break;
119: case 8 : // hexahedron
120: _numFaceVertices = 4; // Face has 4 vertices
121: break;
122: case 10 : // quadratic tetrahedron
123: _numFaceVertices = 6; // Face has 6 vertices
124: break;
125: case 27 : // quadratic hexahedron
126: _numFaceVertices = 9; // Face has 9 vertices
127: break;
128: default :
129: throw ALE::Exception("Invalid number of face corners");
130: }
131: break;
132: default:
133: throw ALE::Exception("Invalid cell dimension");
134: }
135: return _numFaceVertices;
136: };
137: // We need this method because we do not use interpolated sieves
138: // - Without interpolation, we cannot say what vertex collections are
139: // faces, and how they are oriented
140: // - Now we read off the list of face vertices IN THE ORDER IN WHICH
141: // THEY APPEAR IN THE CELL
142: // - This leads to simple algorithms for simplices and quads to check
143: // orientation since these sets are always valid faces
144: // - This is not true with hexes, so we just sort and check explicit cases
145: // - This means for hexes that we have to alter the vertex container as well
146: template<typename Mesh>
147: static bool faceOrientation(const point_type& cell, const Obj<Mesh>& mesh, const int numCorners,
148: const int indices[], const int oppositeVertex, PointArray *origVertices, PointArray *faceVertices) {
149: const int cellDim = mesh->getDimension();
150: const int debug = mesh->debug();
151: bool posOrient = false;
153: if (debug) std::cout << "cellDim: " << cellDim << ", numCorners: " << numCorners << std::endl;
155: if (cellDim == numCorners-1) {
156: // Simplices
157: posOrient = !(oppositeVertex%2);
158: } else if (cellDim == 1 && numCorners == 3) {
159: posOrient = true;
160: } else if (cellDim == 2 && numCorners == 4) {
161: // Quads
162: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
163: posOrient = true;
164: } else if ((indices[0] == 3) && (indices[1] == 0)) {
165: posOrient = true;
166: } else {
167: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
168: posOrient = false;
169: } else {
170: throw ALE::Exception("Invalid quad crossedge");
171: }
172: }
173: } else if (cellDim == 2 && numCorners == 6) {
174: // Quadratic triangle (I hate this)
175: // Edges are determined by the first 2 vertices (corners of edges)
176: const int faceSizeTri = 3;
177: int sortedIndices[3];
178: bool found = false;
179: int faceVerticesTriSorted[9] = {
180: 0, 3, 4, // bottom
181: 1, 4, 5, // right
182: 2, 3, 5, // left
183: };
184: int faceVerticesTri[9] = {
185: 0, 3, 4, // bottom
186: 1, 4, 5, // right
187: 2, 5, 3, // left
188: };
190: for(int i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
191: std::sort(sortedIndices, sortedIndices+faceSizeTri);
192: for (int iFace=0; iFace < 4; ++iFace) {
193: const int ii = iFace*faceSizeTri;
194: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
195: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
196: if (debug) {
197: if (iFace == 0) std::cout << "Bottom edge" << std::endl;
198: else if (iFace == 1) std::cout << "Right edge" << std::endl;
199: else if (iFace == 2) std::cout << "Left edge" << std::endl;
200: } // if
201: for (int fVertex=0; fVertex < faceSizeTri; ++fVertex)
202: for (int cVertex=0; cVertex < faceSizeTri; ++cVertex)
203: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
204: faceVertices->push_back((*origVertices)[cVertex]);
205: break;
206: } // if
207: found = true;
208: break;
209: } // if
210: } // for
211: if (!found) {throw ALE::Exception("Invalid tri crossface");}
212: return true;
213: } else if (cellDim == 2 && numCorners == 9) {
214: // Quadratic quad (I hate this)
215: // Edges are determined by the first 2 vertices (corners of edges)
216: const int faceSizeQuad = 3;
217: int sortedIndices[3];
218: bool found = false;
219: int faceVerticesQuadSorted[12] = {
220: 0, 1, 4, // bottom
221: 1, 2, 5, // right
222: 2, 3, 6, // top
223: 0, 3, 7, // left
224: };
225: int faceVerticesQuad[12] = {
226: 0, 1, 4, // bottom
227: 1, 2, 5, // right
228: 2, 3, 6, // top
229: 3, 0, 7, // left
230: };
232: for(int i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
233: std::sort(sortedIndices, sortedIndices+faceSizeQuad);
234: for (int iFace=0; iFace < 4; ++iFace) {
235: const int ii = iFace*faceSizeQuad;
236: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
237: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
238: if (debug) {
239: if (iFace == 0) std::cout << "Bottom edge" << std::endl;
240: else if (iFace == 1) std::cout << "Right edge" << std::endl;
241: else if (iFace == 2) std::cout << "Top edge" << std::endl;
242: else if (iFace == 3) std::cout << "Left edge" << std::endl;
243: } // if
244: for (int fVertex=0; fVertex < faceSizeQuad; ++fVertex)
245: for (int cVertex=0; cVertex < faceSizeQuad; ++cVertex)
246: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
247: faceVertices->push_back((*origVertices)[cVertex]);
248: break;
249: } // if
250: found = true;
251: break;
252: } // if
253: } // for
254: if (!found) {throw ALE::Exception("Invalid quad crossface");}
255: return true;
256: } else if (cellDim == 3 && numCorners == 8) {
257: // Hexes
258: // A hex is two oriented quads with the normal of the first
259: // pointing up at the second.
260: //
261: // 7---6
262: // /| /|
263: // 4---5 |
264: // | 3-|-2
265: // |/ |/
266: // 0---1
267: //
268: // Faces are determined by the first 4 vertices (corners of faces)
269: const int faceSizeHex = 4;
270: int sortedIndices[4];
271: bool found = false;
272: int faceVerticesHexSorted[24] = {
273: 0, 1, 2, 3, // bottom
274: 4, 5, 6, 7, // top
275: 0, 1, 4, 5, // front
276: 1, 2, 5, 6, // right
277: 2, 3, 6, 7, // back
278: 0, 3, 4, 7, // left
279: };
280: int faceVerticesHex[24] = {
281: 3, 2, 1, 0, // bottom
282: 4, 5, 6, 7, // top
283: 0, 1, 5, 4, // front
284: 1, 2, 6, 5, // right
285: 2, 3, 7, 6, // back
286: 3, 0, 4, 7, // left
287: };
289: for(int i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
290: std::sort(sortedIndices, sortedIndices+faceSizeHex);
291: for (int iFace=0; iFace < 6; ++iFace) {
292: const int ii = iFace*faceSizeHex;
293: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
294: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
295: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
296: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
297: if (debug) {
298: if (iFace == 0) std::cout << "Bottom quad" << std::endl;
299: else if (iFace == 1) std::cout << "Top quad" << std::endl;
300: else if (iFace == 2) std::cout << "Front quad" << std::endl;
301: else if (iFace == 3) std::cout << "Right quad" << std::endl;
302: else if (iFace == 4) std::cout << "Back quad" << std::endl;
303: else if (iFace == 5) std::cout << "Left quad" << std::endl;
304: } // if
305: for (int fVertex=0; fVertex < faceSizeHex; ++fVertex)
306: for (int cVertex=0; cVertex < faceSizeHex; ++cVertex)
307: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
308: faceVertices->push_back((*origVertices)[cVertex]);
309: break;
310: } // if
311: found = true;
312: break;
313: } // if
314: } // for
315: if (!found) {throw ALE::Exception("Invalid hex crossface");}
316: return true;
317: } else if (cellDim == 3 && numCorners == 10) {
318: // Quadratic tet
319: // Faces are determined by the first 3 vertices (corners of faces)
320: const int faceSizeTet = 6;
321: int sortedIndices[6];
322: bool found = false;
323: int faceVerticesTetSorted[24] = {
324: 0, 1, 2, 6, 7, 8, // bottom
325: 0, 3, 4, 6, 7, 9, // front
326: 1, 4, 5, 7, 8, 9, // right
327: 2, 3, 5, 6, 8, 9, // left
328: };
329: int faceVerticesTet[24] = {
330: 0, 1, 2, 6, 7, 8, // bottom
331: 0, 4, 3, 6, 7, 9, // front
332: 1, 5, 4, 7, 8, 9, // right
333: 2, 3, 5, 8, 6, 9, // left
334: };
336: for(int i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
337: std::sort(sortedIndices, sortedIndices+faceSizeTet);
338: for (int iFace=0; iFace < 6; ++iFace) {
339: const int ii = iFace*faceSizeTet;
340: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
341: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
342: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
343: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
344: if (debug) {
345: if (iFace == 0) std::cout << "Bottom tri" << std::endl;
346: else if (iFace == 1) std::cout << "Front tri" << std::endl;
347: else if (iFace == 2) std::cout << "Right tri" << std::endl;
348: else if (iFace == 3) std::cout << "Left tri" << std::endl;
349: } // if
350: for (int fVertex=0; fVertex < faceSizeTet; ++fVertex)
351: for (int cVertex=0; cVertex < faceSizeTet; ++cVertex)
352: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
353: faceVertices->push_back((*origVertices)[cVertex]);
354: break;
355: } // if
356: found = true;
357: break;
358: } // if
359: } // for
360: if (!found) {throw ALE::Exception("Invalid tet crossface");}
361: return true;
362: } else if (cellDim == 3 && numCorners == 27) {
363: // Quadratic hexes (I hate this)
364: // A hex is two oriented quads with the normal of the first
365: // pointing up at the second.
366: //
367: // 7---6
368: // /| /|
369: // 4---5 |
370: // | 3-|-2
371: // |/ |/
372: // 0---1
373: //
374: // Faces are determined by the first 4 vertices (corners of faces)
375: const int faceSizeQuadHex = 9;
376: int sortedIndices[9];
377: bool found = false;
378: int faceVerticesQuadHexSorted[54] = {
379: 0, 1, 2, 3, 8, 9, 10, 11, 24, // bottom
380: 4, 5, 6, 7, 12, 13, 14, 15, 25, // top
381: 0, 1, 4, 5, 8, 12, 16, 17, 22, // front
382: 1, 2, 5, 6, 9, 13, 17, 18, 21, // right
383: 2, 3, 6, 7, 10, 14, 18, 19, 23, // back
384: 0, 3, 4, 7, 11, 15, 16, 19, 20, // left
385: };
386: int faceVerticesQuadHex[54] = {
387: 3, 2, 1, 0, 10, 9, 8, 11, 24, // bottom
388: 4, 5, 6, 7, 12, 13, 14, 15, 25, // top
389: 0, 1, 5, 4, 8, 17, 12, 16, 22, // front
390: 1, 2, 6, 5, 9, 18, 13, 17, 21, // right
391: 2, 3, 7, 6, 10, 19, 14, 18, 23, // back
392: 3, 0, 4, 7, 11, 16, 15, 19, 20 // left
393: };
395: for (int i=0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
396: std::sort(sortedIndices, sortedIndices+faceSizeQuadHex);
397: for (int iFace=0; iFace < 6; ++iFace) {
398: const int ii = iFace*faceSizeQuadHex;
399: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
400: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
401: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
402: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
403: if (debug) {
404: if (iFace == 0) std::cout << "Bottom quad" << std::endl;
405: else if (iFace == 1) std::cout << "Top quad" << std::endl;
406: else if (iFace == 2) std::cout << "Front quad" << std::endl;
407: else if (iFace == 3) std::cout << "Right quad" << std::endl;
408: else if (iFace == 4) std::cout << "Back quad" << std::endl;
409: else if (iFace == 5) std::cout << "Left quad" << std::endl;
410: } // if
411: for (int fVertex=0; fVertex < faceSizeQuadHex; ++fVertex)
412: for (int cVertex=0; cVertex < faceSizeQuadHex; ++cVertex)
413: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
414: faceVertices->push_back((*origVertices)[cVertex]);
415: break;
416: } // if
417: found = true;
418: break;
419: } // if
420: } // for
421: if (!found) {throw ALE::Exception("Invalid hex crossface");}
422: return true;
423: } else {
424: throw ALE::Exception("Unknown cell type for faceOrientation().");
425: }
426: if (!posOrient) {
427: if (debug) std::cout << " Reversing initial face orientation" << std::endl;
428: faceVertices->insert(faceVertices->end(), (*origVertices).rbegin(), (*origVertices).rend());
429: } else {
430: if (debug) std::cout << " Keeping initial face orientation" << std::endl;
431: faceVertices->insert(faceVertices->end(), (*origVertices).begin(), (*origVertices).end());
432: }
433: return posOrient;
434: };
435: // Given a cell and a face, as a set of vertices,
436: // return the oriented face, as a set of vertices, in faceVertices
437: // The orientation is such that the face normal points out of the cell
438: template<typename Mesh, typename FaceType>
439: static bool getOrientedFace(const Obj<Mesh>& mesh, const point_type& cell, FaceType face,
440: const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
441: {
442: FaceVisitor<typename Mesh::sieve_type,FaceType> v(face, *origVertices, *faceVertices, indices, mesh->debug());
444: origVertices->clear();
445: faceVertices->clear();
446: mesh->getSieve()->cone(cell, v);
447: return faceOrientation(cell, mesh, numCorners, indices, v.getOppositeVertex(), origVertices, faceVertices);
448: };
449: template<typename FaceType>
450: static bool getOrientedFace(const Obj<FlexMesh>& mesh, const point_type& cell, FaceType face,
451: const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
452: {
453: const Obj<typename FlexMesh::sieve_type::traits::coneSequence>& cone = mesh->getSieve()->cone(cell);
454: const int debug = mesh->debug();
455: int v = 0;
456: int oppositeVertex;
458: origVertices->clear();
459: faceVertices->clear();
460: for(typename FlexMesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter, ++v) {
461: if (face->find(*v_iter) != face->end()) {
462: if (debug) std::cout << " vertex " << *v_iter << std::endl;
463: indices[origVertices->size()] = v;
464: origVertices->insert(origVertices->end(), *v_iter);
465: } else {
466: if (debug) std::cout << " vertex " << *v_iter << std::endl;
467: oppositeVertex = v;
468: }
469: }
470: return faceOrientation(cell, mesh, numCorners, indices, oppositeVertex, origVertices, faceVertices);
471: };
472: template<typename Sieve>
473: static void insertFace(const Obj<mesh_type>& mesh, const Obj<Sieve>& subSieve, const Obj<PointSet>& face, point_type& f,
474: const point_type& cell, const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices)
475: {
476: const Obj<typename Sieve::supportSet> preFace = subSieve->nJoin1(face);
477: const int debug = subSieve->debug();
479: if (preFace->size() > 1) {
480: throw ALE::Exception("Invalid fault sieve: Multiple faces from vertex set");
481: } else if (preFace->size() == 1) {
482: // Add the other cell neighbor for this face
483: subSieve->addArrow(*preFace->begin(), cell);
484: } else if (preFace->size() == 0) {
485: if (debug) std::cout << " Orienting face " << f << std::endl;
486: try {
487: getOrientedFace(mesh, cell, face, numCorners, indices, origVertices, faceVertices);
488: if (debug) std::cout << " Adding face " << f << std::endl;
489: int color = 0;
490: for(typename PointArray::const_iterator f_iter = faceVertices->begin(); f_iter != faceVertices->end(); ++f_iter) {
491: if (debug) std::cout << " vertex " << *f_iter << std::endl;
492: subSieve->addArrow(*f_iter, f, color++);
493: }
494: subSieve->addArrow(f, cell);
495: f++;
496: } catch (ALE::Exception e) {
497: if (debug) std::cout << " Did not add invalid face " << f << std::endl;
498: }
499: }
500: };
501: public:
502: class FaceInserter {
503: #if 0
504: public:
505: typedef Mesh_ mesh_type;
506: typedef typename mesh_type::sieve_type sieve_type;
507: typedef typename mesh_type::point_type point_type;
508: typedef std::set<point_type> PointSet;
509: typedef std::vector<point_type> PointArray;
510: #endif
511: protected:
512: const Obj<mesh_type> mesh;
513: const Obj<sieve_type> sieve;
514: const Obj<sieve_type> subSieve;
515: point_type& f;
516: const point_type cell;
517: const int numCorners;
518: int *indices;
519: PointArray *origVertices;
520: PointArray *faceVertices;
521: PointSet *subCells;
522: const int debug;
523: public:
524: FaceInserter(const Obj<mesh_type>& mesh, const Obj<sieve_type>& sieve, const Obj<sieve_type>& subSieve, point_type& f, const point_type& cell, const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices, PointSet *subCells) : mesh(mesh), sieve(sieve), subSieve(subSieve), f(f), cell(cell), numCorners(numCorners), indices(indices), origVertices(origVertices), faceVertices(faceVertices), subCells(subCells), debug(mesh->debug()) {};
525: virtual ~FaceInserter() {};
526: public:
527: void operator()(const Obj<PointArray>& face) {
528: const Obj<typename sieve_type::supportSet> sievePreFace = sieve->nJoin1(face);
530: if (sievePreFace->size() == 1) {
531: if (debug) std::cout << " Contains a boundary face on the submesh" << std::endl;
532: PointSet faceSet(face->begin(), face->end());
533: ALE::Selection<mesh_type>::insertFace(mesh, subSieve, faceSet, f, cell, numCorners, indices, origVertices, faceVertices);
534: subCells->insert(cell);
535: }
536: };
537: };
538: template<typename Sieve>
539: class FaceInserterV {
540: protected:
541: const Obj<mesh_type>& mesh;
542: const Obj<sieve_type>& sieve;
543: const Obj<Sieve>& subSieve;
544: point_type& f;
545: const point_type cell;
546: const int numCorners;
547: int *indices;
548: PointArray *origVertices;
549: PointArray *faceVertices;
550: PointSet *subCells;
551: const int debug;
552: public:
553: FaceInserterV(const Obj<mesh_type>& mesh, const Obj<sieve_type>& sieve, const Obj<Sieve>& subSieve, point_type& f, const point_type& cell, const int numCorners, int indices[], PointArray *origVertices, PointArray *faceVertices, PointSet *subCells) : mesh(mesh), sieve(sieve), subSieve(subSieve), f(f), cell(cell), numCorners(numCorners), indices(indices), origVertices(origVertices), faceVertices(faceVertices), subCells(subCells), debug(mesh->debug()) {};
554: virtual ~FaceInserterV() {};
555: public:
556: void operator()(const Obj<PointArray>& face) {
557: ISieveVisitor::PointRetriever<sieve_type> jV(sieve->getMaxSupportSize());
559: sieve->join(*face, jV);
560: if (jV.getSize() == 1) {
561: if (debug) std::cout << " Contains a boundary face on the submesh" << std::endl;
562: PointSet faceSet(face->begin(), face->end());
563: ALE::Selection<mesh_type>::insertFace(mesh, subSieve, faceSet, f, cell, numCorners, indices, origVertices, faceVertices);
564: subCells->insert(cell);
565: }
566: };
567: };
568: protected:
569: static int binomial(const int i, const int j) {
570: assert(j <= i);
571: assert(i < 34);
572: if (j == 0) {
573: return 1;
574: } else if (j == i) {
575: return 1;
576: } else {
577: return binomial(i-1, j) + binomial(i-1, j-1);
578: }
579: };
580: public:
581: // This takes in a section and creates a submesh from the vertices in the section chart
582: // This is a hyperplane of one dimension lower than the mesh
583: static Obj<mesh_type> submesh_uninterpolated(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1, const bool boundaryFaces = true) {
584: // A challenge here is to coordinate the extra numbering of new faces
585: // In serial, it is enough to number after the last point:
586: // Use sieve->base()->size() + sieve->cap()->size(), or determine the greatest point
587: // In parallel, there are two steps:
588: // 1) Use the serial result, and reduce either with add (for size) or max (for greatest point)
589: // 2) Determine how many faces will be created on each process
590: // This will be bounded by C(numCorners, faceSize)*submeshCells
591: // Thus it looks like we should first accumulate submeshCells, and then create faces
592: typedef typename mesh_type::label_type label_type;
593: typedef typename int_section_type::chart_type chart_type;
594: const int dim = (dimension > 0) ? dimension : mesh->getDimension()-1;
595: const Obj<sieve_type>& sieve = mesh->getSieve();
596: Obj<mesh_type> submesh = new mesh_type(mesh->comm(), dim, mesh->debug());
597: Obj<sieve_type> subSieve = new sieve_type(mesh->comm(), mesh->debug());
598: const bool censor = mesh->hasLabel("censored depth");
599: const Obj<label_type>& depthLabel = censor ? mesh->getLabel("censored depth") : mesh->getLabel("depth");
600: const int depth = mesh->depth();
601: const int height = mesh->height();
602: const chart_type& chart = label->getChart();
603: const int numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), depth)->size();
604: const int faceSize = numFaceVertices(mesh);
605: Obj<PointSet> face = new PointSet();
606: int f = sieve->base()->size() + sieve->cap()->size();
607: const int debug = mesh->debug();
608: int *indices = new int[faceSize];
609: PointArray origVertices, faceVertices;
610: PointSet submeshVertices, submeshCells;
613: const typename chart_type::iterator chartEnd = chart.end();
614: for(typename chart_type::iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
615: if (label->getFiberDimension(*c_iter)) submeshVertices.insert(*c_iter);
616: }
617: const typename PointSet::const_iterator svBegin = submeshVertices.begin();
618: const typename PointSet::const_iterator svEnd = submeshVertices.end();
620: for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
621: const Obj<typename sieveAlg::supportArray>& cells = sieveAlg::nSupport(mesh, *sv_iter, depth);
622: const typename sieveAlg::supportArray::iterator cBegin = cells->begin();
623: const typename sieveAlg::supportArray::iterator cEnd = cells->end();
625: if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
626: for(typename sieveAlg::supportArray::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
627: if (debug) std::cout << " Checking cell " << *c_iter << std::endl;
628: if (submeshCells.find(*c_iter) != submeshCells.end()) continue;
629: if (censor && (!mesh->getValue(depthLabel, *c_iter))) continue;
630: const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *c_iter, height);
631: const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
632: const typename sieveAlg::coneArray::iterator vEnd = cone->end();
634: face->clear();
635: for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
636: if (submeshVertices.find(*v_iter) != svEnd) {
637: if (debug) std::cout << " contains submesh vertex " << *v_iter << std::endl;
638: face->insert(face->end(), *v_iter);
639: }
640: }
641: if ((int) face->size() > faceSize) {
642: if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
643: if (debug) std::cout << " Has all boundary faces on the submesh" << std::endl;
644: submeshCells.insert(*c_iter);
645: }
646: if ((int) face->size() == faceSize) {
647: if (debug) std::cout << " Contains a face on the submesh" << std::endl;
648: submeshCells.insert(*c_iter);
649: }
650: }
651: }
652: if (mesh->commSize() > 1) {
653: int localF = f;
654: int localFaces = binomial(numCorners, faceSize)*submeshCells.size();
655: int maxFaces;
657: MPI_Allreduce(&localF, &f, 1, MPI_INT, MPI_SUM, mesh->comm());
658: // 2) Determine how many faces will be created on each process
659: // This will be bounded by faceSize*submeshCells
660: // Thus it looks like we should first accumulate submeshCells, and then create faces
661: MPI_Scan(&localFaces, &maxFaces, 1, MPI_INT, MPI_SUM, mesh->comm());
662: f += maxFaces - localFaces;
663: }
664: for(typename PointSet::const_iterator c_iter = submeshCells.begin(); c_iter != submeshCells.end(); ++c_iter) {
665: if (debug) std::cout << " Processing submesh cell " << *c_iter << std::endl;
666: const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *c_iter, height);
667: const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
668: const typename sieveAlg::coneArray::iterator vEnd = cone->end();
670: face->clear();
671: for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
672: if (submeshVertices.find(*v_iter) != svEnd) {
673: if (debug) std::cout << " contains submesh vertex " << *v_iter << std::endl;
674: face->insert(face->end(), *v_iter);
675: }
676: }
677: if ((int) face->size() > faceSize) {
678: // Here we allow a set of vertices to lie completely on a boundary cell (like a corner tetrahedron)
679: // We have to take all the faces, and discard those in the interior
680: FaceInserter inserter(mesh, sieve, subSieve, f, *c_iter, numCorners, indices, &origVertices, &faceVertices, &submeshCells);
681: PointArray faceVec(face->begin(), face->end());
683: subsets(faceVec, faceSize, inserter);
684: }
685: if ((int) face->size() == faceSize) {
686: insertFace(mesh, subSieve, face, f, *c_iter, numCorners, indices, &origVertices, &faceVertices);
687: }
688: }
689: delete [] indices;
690: submesh->setSieve(subSieve);
691: submesh->stratify();
692: if (debug) submesh->view("Submesh");
693: return submesh;
694: };
695: // This takes in a section and creates a submesh from the vertices in the section chart
696: // This is a hyperplane of one dimension lower than the mesh
697: static Obj<mesh_type> submesh_interpolated(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1, const bool boundaryFaces = true) {
698: const int debug = mesh->debug();
699: const int depth = mesh->depth();
700: const int height = mesh->height();
701: const typename int_section_type::chart_type& chart = label->getChart();
702: const typename int_section_type::chart_type::const_iterator chartEnd = chart.end();
703: const Obj<PointSet> submeshFaces = new PointSet();
704: PointSet submeshVertices;
706: for(typename int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
707: //assert(!mesh->depth(*c_iter));
708: submeshVertices.insert(*c_iter);
709: }
710: const typename PointSet::const_iterator svBegin = submeshVertices.begin();
711: const typename PointSet::const_iterator svEnd = submeshVertices.end();
713: for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
714: const Obj<typename sieveAlg::supportArray>& faces = sieveAlg::nSupport(mesh, *sv_iter, depth-1);
715: const typename sieveAlg::supportArray::iterator fBegin = faces->begin();
716: const typename sieveAlg::supportArray::iterator fEnd = faces->end();
717:
718: if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
719: for(typename sieveAlg::supportArray::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
720: if (debug) std::cout << " Checking face " << *f_iter << std::endl;
721: if (submeshFaces->find(*f_iter) != submeshFaces->end()) continue;
722: const Obj<typename sieveAlg::coneArray>& cone = sieveAlg::nCone(mesh, *f_iter, height-1);
723: const typename sieveAlg::coneArray::iterator vBegin = cone->begin();
724: const typename sieveAlg::coneArray::iterator vEnd = cone->end();
725: bool found = true;
727: for(typename sieveAlg::coneArray::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
728: if (submeshVertices.find(*v_iter) != svEnd) {
729: if (debug) std::cout << " contains submesh vertex " << *v_iter << std::endl;
730: } else {
731: found = false;
732: }
733: }
734: if (found) {
735: if (boundaryFaces) {throw ALE::Exception("Not finished: should check that it is a boundary face");}
736: if (debug) std::cout << " Is a face on the submesh" << std::endl;
737: submeshFaces->insert(*f_iter);
738: }
739: }
740: }
741: return submesh(mesh, submeshFaces, mesh->getDimension()-1);
742: };
743: template<typename output_mesh_type>
744: static Obj<output_mesh_type> submeshV_uninterpolated(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1, const bool boundaryFaces = true) {
745: typedef typename mesh_type::label_type label_type;
746: typedef typename int_section_type::chart_type chart_type;
747: const int dim = (dimension > 0) ? dimension : mesh->getDimension()-1;
748: const Obj<sieve_type>& sieve = mesh->getSieve();
749: Obj<FlexMesh> submesh = new FlexMesh(mesh->comm(), dim, mesh->debug());
750: Obj<typename FlexMesh::sieve_type> subSieve = new typename FlexMesh::sieve_type(mesh->comm(), mesh->debug());
751: const bool censor = mesh->hasLabel("censored depth");
752: const Obj<label_type>& depthLabel = censor ? mesh->getLabel("censored depth") : mesh->getLabel("depth");
753: const chart_type& chart = label->getChart();
754: const int numCorners = mesh->getNumCellCorners();
755: const int faceSize = numFaceVertices(mesh);
756: Obj<PointSet> face = new PointSet();
757: int f = sieve->getBaseSize() + sieve->getCapSize();
758: const int debug = mesh->debug();
759: int *indices = new int[faceSize];
760: PointArray origVertices, faceVertices;
761: PointSet submeshVertices, submeshCells;
763: const typename chart_type::const_iterator chartEnd = chart.end();
764: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
765: if (label->getFiberDimension(*c_iter)) submeshVertices.insert(*c_iter);
766: }
767: const typename PointSet::const_iterator svBegin = submeshVertices.begin();
768: const typename PointSet::const_iterator svEnd = submeshVertices.end();
769: typename ISieveVisitor::PointRetriever<sieve_type> sV(sieve->getMaxSupportSize());
770: typename ISieveVisitor::PointRetriever<sieve_type> cV(sieve->getMaxConeSize());
772: for(typename PointSet::const_iterator sv_iter = svBegin; sv_iter != svEnd; ++sv_iter) {
773: sieve->support(*sv_iter, sV);
774: const int numCells = sV.getSize();
775: const point_type *cells = sV.getPoints();
776:
777: if (debug) std::cout << "Checking submesh vertex " << *sv_iter << std::endl;
778: for(int c = 0; c < numCells; ++c) {
779: if (debug) std::cout << " Checking cell " << cells[c] << std::endl;
780: if (submeshCells.find(cells[c]) != submeshCells.end()) continue;
781: if (censor && (!mesh->getValue(depthLabel, cells[c]))) continue;
782: sieve->cone(cells[c], cV);
783: const int numVertices = cV.getSize();
784: const point_type *vertices = cV.getPoints();
786: face->clear();
787: for(int v = 0; v < numVertices; ++v) {
788: if (submeshVertices.find(vertices[v]) != svEnd) {
789: if (debug) std::cout << " contains submesh vertex " << vertices[v] << std::endl;
790: face->insert(face->end(), vertices[v]);
791: }
792: }
793: if ((int) face->size() > faceSize) {
794: if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
795: if (debug) std::cout << " Has all boundary faces on the submesh" << std::endl;
796: submeshCells.insert(cells[c]);
797: }
798: if ((int) face->size() == faceSize) {
799: if (debug) std::cout << " Contains a face on the submesh" << std::endl;
800: submeshCells.insert(cells[c]);
801: }
802: cV.clear();
803: }
804: sV.clear();
805: }
806: if (mesh->commSize() > 1) {
807: int localF = f;
808: int localFaces = binomial(numCorners, faceSize)*submeshCells.size();
809: int maxFaces;
811: MPI_Allreduce(&localF, &f, 1, MPI_INT, MPI_SUM, mesh->comm());
812: // 2) Determine how many faces will be created on each process
813: // This will be bounded by faceSize*submeshCells
814: // Thus it looks like we should first accumulate submeshCells, and then create faces
815: MPI_Scan(&localFaces, &maxFaces, 1, MPI_INT, MPI_SUM, mesh->comm());
816: f += maxFaces - localFaces;
817: }
818: for(typename PointSet::const_iterator c_iter = submeshCells.begin(); c_iter != submeshCells.end(); ++c_iter) {
819: if (debug) std::cout << " Processing submesh cell " << *c_iter << std::endl;
820: sieve->cone(*c_iter, cV);
821: const int numVertices = cV.getSize();
822: const point_type *vertices = cV.getPoints();
824: face->clear();
825: for(int v = 0; v < numVertices; ++v) {
826: if (submeshVertices.find(vertices[v]) != svEnd) {
827: if (debug) std::cout << " contains submesh vertex " << vertices[v] << std::endl;
828: face->insert(face->end(), vertices[v]);
829: }
830: }
831: if ((int) face->size() > faceSize) {
832: if (!boundaryFaces) throw ALE::Exception("Invalid fault mesh: Too many vertices of an element on the fault");
833: // Here we allow a set of vertices to lie completely on a boundary cell (like a corner tetrahedron)
834: // We have to take all the faces, and discard those in the interior
835: FaceInserterV<FlexMesh::sieve_type> inserter(mesh, sieve, subSieve, f, *c_iter, numCorners, indices, &origVertices, &faceVertices, &submeshCells);
836: PointArray faceVec(face->begin(), face->end());
838: subsets(faceVec, faceSize, inserter);
839: }
840: if ((int) face->size() == faceSize) {
841: insertFace(mesh, subSieve, face, f, *c_iter, numCorners, indices, &origVertices, &faceVertices);
842: }
843: cV.clear();
844: }
845: delete [] indices;
846: submesh->setSieve(subSieve);
847: submesh->stratify();
848: if (debug) submesh->view("Submesh");
850: Obj<output_mesh_type> isubmesh = new output_mesh_type(submesh->comm(), submesh->getDimension(), submesh->debug());
851: Obj<typename output_mesh_type::sieve_type> isieve = new typename output_mesh_type::sieve_type(submesh->comm(), submesh->debug());
852: std::map<typename output_mesh_type::point_type,typename output_mesh_type::point_type> renumbering;
853: isubmesh->setSieve(isieve);
854: ALE::ISieveConverter::convertMesh(*submesh, *isubmesh, renumbering, false);
855: return isubmesh;
856: };
857: public:
858: // This takes in a section and creates a submesh from the vertices in the section chart
859: // This is a hyperplane of one dimension lower than the mesh
860: static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
861: const int dim = mesh->getDimension();
862: const int depth = mesh->depth();
864: if (dim == depth) {
865: return submesh_interpolated(mesh, label, dimension, false);
866: } else if (depth == 1) {
867: return submesh_uninterpolated(mesh, label, dimension);
868: }
869: throw ALE::Exception("Cannot handle partially interpolated meshes");
870: };
871: template<typename output_mesh_type>
872: static Obj<output_mesh_type> submeshV(const Obj<mesh_type>& mesh, const Obj<int_section_type>& label, const int dimension = -1) {
873: const int dim = mesh->getDimension();
874: const int depth = mesh->depth();
876: #if 0
877: if (dim == depth) {
878: //return submesh_interpolated(mesh, label, dimension, false);
879: throw ALE::Exception("Cannot handle interpolated meshes");
880: } else if (depth == 1) {
881: return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
882: }
883: #else
884: if (depth == 1) {
885: return submeshV_uninterpolated<output_mesh_type>(mesh, label, dimension);
886: } else if (dim == depth) {
887: //return submesh_interpolated(mesh, label, dimension, false);
888: throw ALE::Exception("Cannot handle interpolated meshes");
889: }
890: #endif
891: throw ALE::Exception("Cannot handle partially interpolated meshes");
892: };
893: // This creates a submesh consisting of the union of the closures of the given points
894: // This mesh is the same dimension as in the input mesh
895: template<typename Points>
896: static Obj<mesh_type> submesh(const Obj<mesh_type>& mesh, const Obj<Points>& points, const int dim = -1) {
897: const Obj<sieve_type>& sieve = mesh->getSieve();
898: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), dim >= 0 ? dim : mesh->getDimension(), mesh->debug());
899: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
900: Obj<PointSet> newPoints = new PointSet();
901: Obj<PointSet> modPoints = new PointSet();
902: Obj<PointSet> tmpPoints;
904: newMesh->setSieve(newSieve);
905: for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
906: newPoints->insert(*p_iter);
907: do {
908: modPoints->clear();
909: for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
910: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*np_iter);
911: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
912: int c = 0;
914: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter, ++c) {
915: newSieve->addArrow(*c_iter, *np_iter, c);
916: }
917: modPoints->insert(cone->begin(), cone->end());
918: }
919: tmpPoints = newPoints;
920: newPoints = modPoints;
921: modPoints = tmpPoints;
922: } while(newPoints->size());
923: newPoints->insert(*p_iter);
924: do {
925: modPoints->clear();
926: for(typename PointSet::iterator np_iter = newPoints->begin(); np_iter != newPoints->end(); ++np_iter) {
927: const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*np_iter);
928: const typename sieve_type::traits::supportSequence::iterator end = support->end();
929: int s = 0;
931: for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != end; ++s_iter, ++s) {
932: newSieve->addArrow(*np_iter, *s_iter, s);
933: }
934: modPoints->insert(support->begin(), support->end());
935: }
936: tmpPoints = newPoints;
937: newPoints = modPoints;
938: modPoints = tmpPoints;
939: } while(newPoints->size());
940: }
941: newMesh->stratify();
942: newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
943: newMesh->setArrowSection("orientation", mesh->getArrowSection("orientation"));
944: return newMesh;
945: };
946: protected:
947: static Obj<mesh_type> boundary_uninterpolated(const Obj<mesh_type>& mesh) {
948: throw ALE::Exception("Not yet implemented");
949: const Obj<typename mesh_type::label_sequence>& cells = mesh->heightStratum(0);
950: const Obj<sieve_type>& sieve = mesh->getSieve();
951: const typename mesh_type::label_sequence::iterator cBegin = cells->begin();
952: const typename mesh_type::label_sequence::iterator cEnd = cells->end();
953: const int faceSize = numFaceVertices(mesh);
955: for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
956: const Obj<typename sieve_type::traits::coneSequence>& vertices = sieve->cone(*c_iter);
957: const typename sieve_type::traits::coneSequence::iterator vBegin = vertices->begin();
958: const typename sieve_type::traits::coneSequence::iterator vEnd = vertices->end();
959: //PointArray cell(vertices->begin(), vertices->end());
961: // For each vertex, gather
962: for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
963: const Obj<typename sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
964: const typename sieve_type::traits::supportSequence::iterator nBegin = neighbors->begin();
965: const typename sieve_type::traits::supportSequence::iterator nEnd = neighbors->end();
967: for(typename sieve_type::traits::supportSequence::iterator n_iter = nBegin; n_iter != nEnd; ++n_iter) {
968: const Obj<typename sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, 1);
970: if (preFace->size() == faceSize) {
971: }
972: }
973: }
974: // For each face
975: // - determine if its legal
976:
977: // - determine if its part of a neighboring cell
978: // - if not, its a boundary face
979: //subsets(cell, faceSize, inserter);
980: }
981: };
982: static void addClosure(const Obj<sieve_type>& sieveA, const Obj<sieve_type>& sieveB, const point_type& p, const int depth = 1) {
983: Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
984: Obj<typename sieve_type::coneSet> next = new typename sieve_type::coneSet();
985: Obj<typename sieve_type::coneSet> tmp;
987: current->insert(p);
988: while(current->size()) {
989: for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
990: const Obj<typename sieve_type::traits::coneSequence>& cone = sieveA->cone(*p_iter);
991: const typename sieve_type::traits::coneSequence::iterator begin = cone->begin();
992: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
994: for(typename sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
995: sieveB->addArrow(*c_iter, *p_iter, c_iter.color());
996: next->insert(*c_iter);
997: }
998: }
999: tmp = current; current = next; next = tmp;
1000: next->clear();
1001: }
1002: if (!depth) {
1003: const Obj<typename sieve_type::traits::supportSequence>& support = sieveA->support(p);
1004: const typename sieve_type::traits::supportSequence::iterator begin = support->begin();
1005: const typename sieve_type::traits::supportSequence::iterator end = support->end();
1006:
1007: for(typename sieve_type::traits::supportSequence::iterator s_iter = begin; s_iter != end; ++s_iter) {
1008: sieveB->addArrow(p, *s_iter, s_iter.color());
1009: next->insert(*s_iter);
1010: }
1011: }
1012: };
1013: static Obj<mesh_type> boundary_interpolated(const Obj<mesh_type>& mesh, const int faceHeight = 1) {
1014: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
1015: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
1016: const Obj<sieve_type>& sieve = mesh->getSieve();
1017: const Obj<typename mesh_type::label_sequence>& faces = mesh->heightStratum(faceHeight);
1018: const typename mesh_type::label_sequence::iterator fBegin = faces->begin();
1019: const typename mesh_type::label_sequence::iterator fEnd = faces->end();
1020: const int depth = faceHeight - mesh->depth();
1022: for(typename mesh_type::label_sequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
1023: const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*f_iter);
1025: if (support->size() == 1) {
1026: addClosure(sieve, newSieve, *f_iter, depth);
1027: }
1028: }
1029: newMesh->setSieve(newSieve);
1030: newMesh->stratify();
1031: return newMesh;
1032: };
1033: template<typename SieveTypeA, typename SieveTypeB>
1034: static void addClosureV(const Obj<SieveTypeA>& sieveA, const Obj<SieveTypeB>& sieveB, const point_type& p, const int depth = 1) {
1035: typedef std::set<typename SieveTypeA::point_type> coneSet;
1036: ALE::ISieveVisitor::PointRetriever<SieveTypeA> cV(std::max(1, sieveA->getMaxConeSize()));
1037: Obj<coneSet> current = new coneSet();
1038: Obj<coneSet> next = new coneSet();
1039: Obj<coneSet> tmp;
1041: current->insert(p);
1042: while(current->size()) {
1043: for(typename coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1044: sieveA->cone(*p_iter, cV);
1045: const typename SieveTypeA::point_type *cone = cV.getPoints();
1047: for(int c = 0; c < (int) cV.getSize(); ++c) {
1048: sieveB->addArrow(cone[c], *p_iter);
1049: next->insert(cone[c]);
1050: }
1051: cV.clear();
1052: }
1053: tmp = current; current = next; next = tmp;
1054: next->clear();
1055: }
1056: if (!depth) {
1057: ALE::ISieveVisitor::PointRetriever<SieveTypeA> sV(std::max(1, sieveA->getMaxSupportSize()));
1059: sieveA->support(p, sV);
1060: const typename SieveTypeA::point_type *support = sV.getPoints();
1061:
1062: for(int s = 0; s < (int) sV.getSize(); ++s) {
1063: sieveB->addArrow(p, support[s]);
1064: }
1065: sV.clear();
1066: }
1067: };
1068: public:
1069: static Obj<mesh_type> boundary(const Obj<mesh_type>& mesh) {
1070: const int dim = mesh->getDimension();
1071: const int depth = mesh->depth();
1073: if (dim == depth) {
1074: return boundary_interpolated(mesh);
1075: } else if (depth == dim+1) {
1076: return boundary_interpolated(mesh, 2);
1077: } else if (depth == 1) {
1078: return boundary_uninterpolated(mesh);
1079: } else if (depth == -1) {
1080: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
1081: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
1083: newMesh->setSieve(newSieve);
1084: newMesh->stratify();
1085: return newMesh;
1086: }
1087: throw ALE::Exception("Cannot handle partially interpolated meshes");
1088: };
1089: template<typename MeshTypeQ>
1090: static Obj<FlexMesh> boundaryV_uninterpolated(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
1091: throw ALE::Exception("Cannot handle uninterpolated meshes");
1092: };
1093: // This method takes in an interpolated mesh, and returns the boundary
1094: template<typename MeshTypeQ>
1095: static Obj<FlexMesh> boundaryV_interpolated(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
1096: Obj<FlexMesh> newMesh = new FlexMesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
1097: Obj<typename FlexMesh::sieve_type> newSieve = new typename FlexMesh::sieve_type(mesh->comm(), mesh->debug());
1098: const Obj<typename MeshTypeQ::sieve_type>& sieve = mesh->getSieve();
1099: const Obj<typename MeshTypeQ::label_sequence>& faces = mesh->heightStratum(faceHeight);
1100: const typename MeshTypeQ::label_sequence::iterator fBegin = faces->begin();
1101: const typename MeshTypeQ::label_sequence::iterator fEnd = faces->end();
1102: const int depth = faceHeight - mesh->depth();
1103: ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
1105: for(typename MeshTypeQ::label_sequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
1106: sieve->support(*f_iter, sV);
1108: if (sV.getSize() == 1) {
1109: addClosureV(sieve, newSieve, *f_iter, depth);
1110: }
1111: sV.clear();
1112: }
1113: newMesh->setSieve(newSieve);
1114: newMesh->stratify();
1115: return newMesh;
1116: };
1117: template<typename MeshTypeQ>
1118: static Obj<FlexMesh> boundaryV(const Obj<MeshTypeQ>& mesh, const int faceHeight = 1) {
1119: const int dim = mesh->getDimension();
1120: const int depth = mesh->depth();
1122: if (dim == depth) {
1123: return boundaryV_interpolated(mesh);
1124: } else if (depth == dim+1) {
1125: return boundaryV_interpolated(mesh, 2);
1126: } else if (depth == 1) {
1127: throw ALE::Exception("Cannot handle uninterpolated meshes");
1128: } else if (depth == -1) {
1129: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), mesh->getDimension()-1, mesh->debug());
1130: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
1132: newMesh->setSieve(newSieve);
1133: newMesh->stratify();
1134: return newMesh;
1135: }
1136: throw ALE::Exception("Cannot handle partially interpolated meshes");
1137: };
1138: public:
1139: static Obj<mesh_type> interpolateMesh(const Obj<mesh_type>& mesh) {
1140: const Obj<sieve_type> sieve = mesh->getSieve();
1141: const int dim = mesh->getDimension();
1142: const int numVertices = mesh->depthStratum(0)->size();
1143: const Obj<typename mesh_type::label_sequence>& cells = mesh->heightStratum(0);
1144: const int numCells = cells->size();
1145: const int firstVertex = numCells;
1146: const int debug = sieve->debug();
1147: Obj<mesh_type> newMesh = new mesh_type(mesh->comm(), dim, mesh->debug());
1148: Obj<sieve_type> newSieve = new sieve_type(mesh->comm(), mesh->debug());
1149: const Obj<typename mesh_type::arrow_section_type>& orientation = newMesh->getArrowSection("orientation");
1151: newMesh->setSieve(newSieve);
1152: // Create a map from dimension to the current element number for that dimension
1153: std::map<int,point_type*> curElement;
1154: std::map<int,PointArray> bdVertices;
1155: std::map<int,PointArray> faces;
1156: std::map<int,oPointArray> oFaces;
1157: int curCell = 0;
1158: int curVertex = firstVertex;
1159: int newElement = firstVertex+numVertices;
1160: int o;
1162: curElement[0] = &curVertex;
1163: curElement[dim] = &curCell;
1164: for(int d = 1; d < dim; d++) {
1165: curElement[d] = &newElement;
1166: }
1167: typename mesh_type::label_sequence::iterator cBegin = cells->begin();
1168: typename mesh_type::label_sequence::iterator cEnd = cells->end();
1170: for(typename mesh_type::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
1171: typename sieve_type::point_type cell = *c_iter;
1172: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
1173: const int corners = cone->size();
1174: const typename sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
1175: const typename sieve_type::traits::coneSequence::iterator vEnd = cone->end();
1177: // Build the cell
1178: bdVertices[dim].clear();
1179: for(typename sieve_type::traits::coneSequence::iterator v_iter = vBegin; v_iter != vEnd; ++v_iter) {
1180: typename sieve_type::point_type vertex(*v_iter);
1182: if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
1183: bdVertices[dim].push_back(vertex);
1184: }
1185: if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}
1187: if (corners != dim+1) {
1188: ALE::SieveBuilder<mesh_type>::buildHexFaces(newSieve, dim, curElement, bdVertices, faces, cell);
1189: // Will need to handle cohesive cells here
1190: } else {
1191: ALE::SieveBuilder<mesh_type>::buildFaces(newSieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
1192: }
1193: }
1194: newMesh->stratify();
1195: newMesh->setRealSection("coordinates", mesh->getRealSection("coordinates"));
1196: return newMesh;
1197: };
1198: };
1199: }
1201: #if 0
1202: namespace ALE {
1203: class MySelection {
1204: public:
1205: typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
1206: typedef ALE::Selection<ALE::Mesh> selection;
1207: typedef ALE::Mesh::sieve_type sieve_type;
1208: typedef ALE::Mesh::int_section_type int_section_type;
1209: typedef ALE::Mesh::real_section_type real_section_type;
1210: typedef std::set<ALE::Mesh::point_type> PointSet;
1211: typedef std::vector<ALE::Mesh::point_type> PointArray;
1212: public:
1213: template<class InputPoints>
1214: static bool _compatibleOrientation(const Obj<Mesh>& mesh,
1215: const ALE::Mesh::point_type& p,
1216: const ALE::Mesh::point_type& q,
1217: const int numFaultCorners,
1218: const int faultFaceSize,
1219: const int faultDepth,
1220: const Obj<InputPoints>& points,
1221: int indices[],
1222: PointArray *origVertices,
1223: PointArray *faceVertices,
1224: PointArray *neighborVertices)
1225: {
1226: typedef ALE::Selection<ALE::Mesh> selection;
1227: const int debug = mesh->debug();
1228: bool compatible;
1230: bool eOrient = selection::getOrientedFace(mesh, p, points, numFaultCorners, indices, origVertices, faceVertices);
1231: bool nOrient = selection::getOrientedFace(mesh, q, points, numFaultCorners, indices, origVertices, neighborVertices);
1233: if (faultFaceSize > 1) {
1234: if (debug) {
1235: for(PointArray::iterator v_iter = faceVertices->begin(); v_iter != faceVertices->end(); ++v_iter) {
1236: std::cout << " face vertex " << *v_iter << std::endl;
1237: }
1238: for(PointArray::iterator v_iter = neighborVertices->begin(); v_iter != neighborVertices->end(); ++v_iter) {
1239: std::cout << " neighbor vertex " << *v_iter << std::endl;
1240: }
1241: }
1242: compatible = !(*faceVertices->begin() == *neighborVertices->begin());
1243: } else {
1244: compatible = !(nOrient == eOrient);
1245: }
1246: return compatible;
1247: };
1248: static void _replaceCell(const Obj<sieve_type>& sieve,
1249: const ALE::Mesh::point_type cell,
1250: std::map<int,int> *vertexRenumber,
1251: const int debug)
1252: {
1253: bool replace = false;
1254: PointArray newVertices;
1256: const Obj<sieve_type::traits::coneSequence>& cCone = sieve->cone(cell);
1258: for(sieve_type::traits::coneSequence::iterator v_iter = cCone->begin(); v_iter != cCone->end(); ++v_iter) {
1259: if (vertexRenumber->find(*v_iter) != vertexRenumber->end()) {
1260: if (debug) std::cout << " vertex " << (*vertexRenumber)[*v_iter] << std::endl;
1261: newVertices.insert(newVertices.end(), (*vertexRenumber)[*v_iter]);
1262: replace = true;
1263: } else {
1264: if (debug) std::cout << " vertex " << *v_iter << std::endl;
1265: newVertices.insert(newVertices.end(), *v_iter);
1266: } // if/else
1267: } // for
1268: if (replace) {
1269: if (debug) std::cout << " Replacing cell " << cell << std::endl;
1270: sieve->clearCone(cell);
1271: int color = 0;
1272: for(PointArray::const_iterator v_iter = newVertices.begin(); v_iter != newVertices.end(); ++v_iter) {
1273: sieve->addArrow(*v_iter, cell, color++);
1274: } // for
1275: }
1276: };
1277: template<class InputPoints>
1278: static void _computeCensoredDepth(const Obj<ALE::Mesh>& mesh,
1279: const Obj<ALE::Mesh::label_type>& depth,
1280: const Obj<ALE::Mesh::sieve_type>& sieve,
1281: const Obj<InputPoints>& points,
1282: const ALE::Mesh::point_type& firstCohesiveCell,
1283: const Obj<std::set<ALE::Mesh::point_type> >& modifiedPoints)
1284: {
1285: modifiedPoints->clear();
1287: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1288: if (*p_iter >= firstCohesiveCell) continue;
1289: // Compute the max depth of the points in the cone of p, and add 1
1290: int d0 = mesh->getValue(depth, *p_iter, -1);
1291: int d1 = mesh->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;
1293: if(d1 != d0) {
1294: mesh->setValue(depth, *p_iter, d1);
1295: modifiedPoints->insert(*p_iter);
1296: }
1297: }
1298: // FIX: We would like to avoid the copy here with support()
1299: if(modifiedPoints->size() > 0) {
1300: _computeCensoredDepth(mesh, depth, sieve, sieve->support(modifiedPoints), firstCohesiveCell, modifiedPoints);
1301: }
1302: };
1303: static void create(const Obj<Mesh>& mesh, Obj<Mesh> fault, const Obj<Mesh::int_section_type>& groupField) {
1304: static PetscLogEvent CreateFaultMesh_Event = 0, OrientFaultMesh_Event = 0, AddCohesivePoints_Event = 0, SplitMesh_Event = 0;
1306: if (!CreateFaultMesh_Event) {
1307: PetscLogEventRegister("CreateFaultMesh", 0,&CreateFaultMesh_Event);
1308: }
1309: if (!OrientFaultMesh_Event) {
1310: PetscLogEventRegister("OrientFaultMesh", 0,&OrientFaultMesh_Event);
1311: }
1312: if (!AddCohesivePoints_Event) {
1313: PetscLogEventRegister("AddCohesivePoints", 0,&AddCohesivePoints_Event);
1314: }
1315: if (!SplitMesh_Event) {
1316: PetscLogEventRegister("SplitMesh", 0,&SplitMesh_Event);
1317: }
1319: const Obj<sieve_type>& sieve = mesh->getSieve();
1320: const int debug = mesh->debug();
1321: int numCorners = 0; // The number of vertices in a mesh cell
1322: int faceSize = 0; // The number of vertices in a mesh face
1323: int *indices = NULL; // The indices of a face vertex set in a cell
1324: PointArray origVertices;
1325: PointArray faceVertices;
1326: PointArray neighborVertices;
1327: const bool constraintCell = false;
1329: if (!mesh->commRank()) {
1330: numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), mesh->depth())->size();
1331: faceSize = selection::numFaceVertices(mesh);
1332: indices = new int[faceSize];
1333: }
1335: //int f = sieve->base()->size() + sieve->cap()->size();
1336: //ALE::Obj<PointSet> face = new PointSet();
1337:
1338: // Create a sieve which captures the fault
1339: PetscLogEventBegin(CreateFaultMesh_Event,0,0,0,0);
1340: fault = ALE::Selection<ALE::Mesh>::submesh(mesh, groupField);
1341: if (debug) {fault->view("Fault mesh");}
1342: PetscLogEventEnd(CreateFaultMesh_Event,0,0,0,0);
1343: // Orient the fault sieve
1344: PetscLogEventBegin(OrientFaultMesh_Event,0,0,0,0);
1345: const Obj<sieve_type>& faultSieve = fault->getSieve();
1346: const ALE::Obj<Mesh::label_sequence>& fFaces = fault->heightStratum(1);
1347: int faultDepth = fault->depth()-1; // Depth of fault cells
1348: int numFaultCorners = 0; // The number of vertices in a fault cell
1350: if (!fault->commRank()) {
1351: numFaultCorners = faultSieve->nCone(*fFaces->begin(), faultDepth)->size();
1352: if (debug) std::cout << " Fault corners " << numFaultCorners << std::endl;
1353: assert(numFaultCorners == faceSize);
1354: }
1355: PetscLogEventEnd(OrientFaultMesh_Event,0,0,0,0);
1357: // Add new shadow vertices and possibly Lagrange multipler vertices
1358: PetscLogEventBegin(AddCohesivePoints_Event,0,0,0,0);
1359: const ALE::Obj<Mesh::label_sequence>& fVertices = fault->depthStratum(0);
1360: const ALE::Obj<std::set<std::string> >& groupNames = mesh->getIntSections();
1361: Mesh::point_type newPoint = sieve->base()->size() + sieve->cap()->size();
1362: std::map<int,int> vertexRenumber;
1363:
1364: for(Mesh::label_sequence::iterator v_iter = fVertices->begin(); v_iter != fVertices->end(); ++v_iter, ++newPoint) {
1365: vertexRenumber[*v_iter] = newPoint;
1366: if (debug) {std::cout << "Duplicating " << *v_iter << " to " << vertexRenumber[*v_iter] << std::endl;}
1368: // Add shadow and constraint vertices (if they exist) to group
1369: // associated with fault
1370: groupField->addPoint(newPoint, 1);
1371: if (constraintCell) groupField->addPoint(newPoint+1, 1);
1373: // Add shadow vertices to other groups, don't add constraint
1374: // vertices (if they exist) because we don't want BC, etc to act
1375: // on constraint vertices
1376: for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
1377: const ALE::Obj<int_section_type>& group = mesh->getIntSection(*name);
1378: if (group->hasPoint(*v_iter)) group->addPoint(newPoint, 1);
1379: }
1380: if (constraintCell) newPoint++;
1381: }
1382: for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
1383: mesh->reallocate(mesh->getIntSection(*name));
1384: }
1386: // Split the mesh along the fault sieve and create cohesive elements
1387: const Obj<ALE::Mesh::label_sequence>& faces = fault->depthStratum(1);
1388: const Obj<ALE::Mesh::arrow_section_type>& orientation = mesh->getArrowSection("orientation");
1389: int firstCohesiveCell = newPoint;
1390: PointSet replaceCells;
1391: PointSet noReplaceCells;
1393: for(ALE::Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter, ++newPoint) {
1394: if (debug) std::cout << "Considering fault face " << *f_iter << std::endl;
1395: const ALE::Obj<sieve_type::traits::supportSequence>& cells = faultSieve->support(*f_iter);
1396: const ALE::Mesh::arrow_section_type::point_type arrow(*cells->begin(), *f_iter);
1397: bool reversed = orientation->restrictPoint(arrow)[0] < 0;
1398: const ALE::Mesh::point_type cell = *cells->begin();
1400: if (debug) std::cout << " Checking orientation against cell " << cell << std::endl;
1401: if (numFaultCorners == 2) reversed = orientation->restrictPoint(arrow)[0] == -2;
1402: if (reversed) {
1403: replaceCells.insert(cell);
1404: noReplaceCells.insert(*(++cells->begin()));
1405: } else {
1406: replaceCells.insert(*(++cells->begin()));
1407: noReplaceCells.insert(cell);
1408: }
1409: //selection::getOrientedFace(mesh, cell, &vertexRenumber, numCorners, indices, &origVertices, &faceVertices);
1410: //const Obj<sieve_type::coneArray> faceCone = faultSieve->nCone(*f_iter, faultDepth);
1412: // Adding cohesive cell (not interpolated)
1413: const Obj<sieve_type::coneArray>& fCone = faultSieve->nCone(*f_iter, faultDepth);
1414: const sieve_type::coneArray::iterator fBegin = fCone->begin();
1415: const sieve_type::coneArray::iterator fEnd = fCone->end();
1416: int color = 0;
1418: if (debug) {std::cout << " Creating cohesive cell " << newPoint << std::endl;}
1419: for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
1420: if (debug) {std::cout << " vertex " << *v_iter << std::endl;}
1421: sieve->addArrow(*v_iter, newPoint, color++);
1422: }
1423: for(sieve_type::coneArray::iterator v_iter = fBegin; v_iter != fEnd; ++v_iter) {
1424: if (debug) {std::cout << " shadow vertex " << vertexRenumber[*v_iter] << std::endl;}
1425: sieve->addArrow(vertexRenumber[*v_iter], newPoint, color++);
1426: }
1427: }
1428: PetscLogEventEnd(AddCohesivePoints_Event,0,0,0,0);
1429: // Replace all cells with a vertex on the fault that share a face with this one, or with one that does
1430: PetscLogEventBegin(SplitMesh_Event,0,0,0,0);
1431: const int_section_type::chart_type& chart = groupField->getChart();
1432: const int_section_type::chart_type::iterator chartEnd = chart.end();
1434: for(PointSet::const_iterator v_iter = chart.begin(); v_iter != chartEnd; ++v_iter) {
1435: bool modified = true;
1437: while(modified) {
1438: modified = false;
1439: const Obj<sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
1440: const sieve_type::traits::supportSequence::iterator end = neighbors->end();
1442: for(sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != end; ++n_iter) {
1443: if (replaceCells.find(*n_iter) != replaceCells.end()) continue;
1444: if (noReplaceCells.find(*n_iter) != noReplaceCells.end()) continue;
1445: if (*n_iter >= firstCohesiveCell) continue;
1446: if (debug) std::cout << " Checking fault neighbor " << *n_iter << std::endl;
1447: // If neighbors shares a faces with anyone in replaceCells, then add
1448: for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
1449: const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, mesh->depth());
1451: if ((int) preFace->size() == faceSize) {
1452: if (debug) std::cout << " Scheduling " << *n_iter << " for replacement" << std::endl;
1453: replaceCells.insert(*n_iter);
1454: modified = true;
1455: break;
1456: }
1457: }
1458: }
1459: }
1460: }
1461: for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
1462: _replaceCell(sieve, *c_iter, &vertexRenumber, debug);
1463: }
1464: if (!fault->commRank()) delete [] indices;
1465: mesh->stratify();
1466: const ALE::Obj<Mesh::label_type>& label = mesh->createLabel(std::string("censored depth"));
1467: const ALE::Obj<PointSet> modifiedPoints = new PointSet();
1468: _computeCensoredDepth(mesh, label, mesh->getSieve(), mesh->getSieve()->roots(), firstCohesiveCell, modifiedPoints);
1469: if (debug) mesh->view("Mesh with Cohesive Elements");
1471: // Fix coordinates
1472: const Obj<real_section_type>& coordinates = mesh->getRealSection("coordinates");
1473: const Obj<ALE::Mesh::label_sequence>& fVertices2 = fault->depthStratum(0);
1475: for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
1476: coordinates->addPoint(vertexRenumber[*v_iter], coordinates->getFiberDimension(*v_iter));
1477: if (constraintCell) {
1478: coordinates->addPoint(vertexRenumber[*v_iter]+1, coordinates->getFiberDimension(*v_iter));
1479: }
1480: }
1481: mesh->reallocate(coordinates);
1482: for(ALE::Mesh::label_sequence::iterator v_iter = fVertices2->begin(); v_iter != fVertices2->end(); ++v_iter) {
1483: coordinates->updatePoint(vertexRenumber[*v_iter], coordinates->restrictPoint(*v_iter));
1484: if (constraintCell) {
1485: coordinates->updatePoint(vertexRenumber[*v_iter]+1, coordinates->restrictPoint(*v_iter));
1486: }
1487: }
1488: if (debug) coordinates->view("Coordinates with shadow vertices");
1489: PetscLogEventEnd(SplitMesh_Event,0,0,0,0);
1490: };
1491: };
1492: };
1493: #endif
1495: #endif