Actual source code: SieveBuilder.hh
petsc-3.4.5 2014-06-29
1: #ifndef included_ALE_SieveBuilder_hh
2: #define included_ALE_SieveBuilder_hh
4: #ifndef included_ALE_ALE_hh
5: #include <sieve/ALE.hh>
6: #endif
8: namespace ALE {
9: template<typename Bundle_>
10: class SieveBuilder {
11: public:
12: typedef Bundle_ bundle_type;
13: typedef typename bundle_type::sieve_type sieve_type;
14: typedef typename bundle_type::arrow_section_type arrow_section_type;
15: typedef std::vector<typename sieve_type::point_type> PointArray;
16: typedef std::pair<typename sieve_type::point_type, int> oPoint_type;
17: typedef std::vector<oPoint_type> oPointArray;
18: public:
19: static void buildHexFaces(Obj<sieve_type> sieve, Obj<arrow_section_type> orientation, int dim, std::map<int, int*>& curElement, std::map<int,PointArray>& bdVertices, std::map<int,oPointArray>& faces, typename sieve_type::point_type& cell, int& cellOrientation) {
20: int debug = sieve->debug();
22: if (debug > 1) {std::cout << " Building hex faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;}
23: faces[dim].clear();
24: if (dim > 3) {
25: throw ALE::Exception("Cannot do hexes of dimension greater than three");
26: } else if (dim > 2) {
27: int nodes[24] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 5, 4,
28: 1, 2, 6, 5, 2, 3, 7, 6, 3, 0, 4, 7};
30: for(int b = 0; b < 6; b++) {
31: typename sieve_type::point_type face;
32: int o = 1;
34: bdVertices[dim-1].clear();
35: for(int c = 0; c < 4; c++) {
36: bdVertices[dim-1].push_back(bdVertices[dim][nodes[b*4+c]]);
37: }
38: if (debug > 1) {std::cout << " boundary hex face " << b << std::endl;}
39: buildHexFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
40: if (debug > 1) {std::cout << " added face " << face << std::endl;}
41: faces[dim].push_back(oPoint_type(face, o));
42: }
43: } else if (dim > 1) {
44: int boundarySize = bdVertices[dim].size();
46: for(int b = 0; b < boundarySize; b++) {
47: typename sieve_type::point_type face;
48: int o = 1;
50: bdVertices[dim-1].clear();
51: for(int c = 0; c < 2; c++) {
52: bdVertices[dim-1].push_back(bdVertices[dim][(b+c)%boundarySize]);
53: }
54: if (debug > 1) {
55: std::cout << " boundary point " << bdVertices[dim][b] << std::endl;
56: std::cout << " boundary vertices";
57: for(int c = 0; c < (int) bdVertices[dim-1].size(); c++) {
58: std::cout << " " << bdVertices[dim-1][c];
59: }
60: std::cout << std::endl;
61: }
62: buildHexFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
63: if (debug > 1) {std::cout << " added face " << face << std::endl;}
64: faces[dim].push_back(oPoint_type(face, o));
65: }
66: } else {
67: if (debug > 1) {std::cout << " Just set faces to boundary in 1d" << std::endl;}
68: typename PointArray::iterator bd_iter = bdVertices[dim].begin();
69: faces[dim].push_back(oPoint_type(*bd_iter, 0));++bd_iter;
70: faces[dim].push_back(oPoint_type(*bd_iter, 0));
71: //faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
72: }
73: if (debug > 1) {
74: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
75: std::cout << " face point " << f_iter->first << " orientation " << f_iter->second << std::endl;
76: }
77: }
78: // We always create the toplevel, so we could short circuit somehow
79: // Should not have to loop here since the meet of just 2 boundary elements is an element
80: typename oPointArray::iterator f_itor = faces[dim].begin();
81: const typename sieve_type::point_type& start = f_itor->first;
82: const typename sieve_type::point_type& next = (++f_itor)->first;
83: Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);
85: if (preElement->size() > 0) {
86: cell = *preElement->begin();
88: const int size = faces[dim].size();
89: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
90: int wrap = size > 2 ? size-1 : 0;
91: int indA = 0, indB = 0;
93: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indA) {
94: if (start == *c_iter) break;
95: }
96: if (debug > 1) {std::cout << " pointA " << start << " indA " << indA << std::endl;}
97: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indB) {
98: if (next == *c_iter) break;
99: }
100: if (debug > 1) {std::cout << " pointB " << next << " indB " << indB << std::endl;}
101: if ((indB - indA == 1) || (indA - indB == wrap)) {
102: if (cellOrientation > 0) {
103: cellOrientation = indA+1;
104: } else {
105: if (dim == 1) {
106: cellOrientation = -2;
107: } else {
108: cellOrientation = -(indA+1);
109: }
110: }
111: } else if ((indA - indB == 1) || (indB - indA == wrap)) {
112: if (debug > 1) {std::cout << " reversing cell orientation" << std::endl;}
113: if (cellOrientation > 0) {
114: cellOrientation = -(indA+1);
115: } else {
116: if (dim == 1) {
117: cellOrientation = 1;
118: } else {
119: cellOrientation = indA+1;
120: }
121: }
122: } else {
123: throw ALE::Exception("Inconsistent orientation");
124: }
125: if (debug > 1) {std::cout << " Found old cell " << cell << " orientation " << cellOrientation << std::endl;}
126: } else {
127: int color = 0;
129: cell = typename sieve_type::point_type((*curElement[dim])++);
130: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
131: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(f_iter->first, cell);
133: sieve->addArrow(f_iter->first, cell, color++);
134: if (f_iter->second) {
135: orientation->addPoint(arrow);
136: orientation->updatePoint(arrow, &(f_iter->second));
137: if (debug > 1) {std::cout << " Orienting arrow (" << f_iter->first << ", " << cell << ") to " << f_iter->second << std::endl;}
138: }
139: }
140: if (cellOrientation > 0) {
141: cellOrientation = 1;
142: } else {
143: cellOrientation = -(dim+1);
144: }
145: if (debug > 1) {std::cout << " Added cell " << cell << " dim " << dim << std::endl;}
146: }
147: };
148: static void buildQuadraticHexFaces(Obj<sieve_type> sieve, Obj<arrow_section_type> orientation, int dim, std::map<int, int*>& curElement, std::map<int,PointArray>& bdVertices, std::map<int,oPointArray>& faces, typename sieve_type::point_type& cell, int& cellOrientation) {
149: int debug = sieve->debug();
151: if (debug > 1) {std::cout << " Building hex faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;}
152: faces[dim].clear();
153: if (dim > 3) {
154: throw ALE::Exception("Cannot do hexes of dimension greater than three");
155: } else if (dim > 2) {
156: const int faceSizeQuadHex = 9;
157: const int numFacesQuadHex = 6;
158: int nodes[54] = {
159: 3, 2, 1, 0, 10, 9, 8, 11, 24, // bottom
160: 4, 5, 6, 7, 12, 13, 14, 15, 25, // top
161: 0, 1, 5, 4, 8, 17, 12, 16, 22, // front
162: 1, 2, 6, 5, 9, 18, 13, 17, 21, // right
163: 2, 3, 7, 6, 10, 19, 14, 18, 23, // back
164: 3, 0, 4, 7, 11, 16, 15, 19, 20 // left
165: };
167: for(int b = 0; b < numFacesQuadHex; b++) {
168: typename sieve_type::point_type face;
169: int o = 1;
171: bdVertices[dim-1].clear();
172: for(int c = 0; c < faceSizeQuadHex; c++) {
173: bdVertices[dim-1].push_back(bdVertices[dim][nodes[b*faceSizeQuadHex+c]]);
174: }
175: if (debug > 1) {std::cout << " boundary hex face " << b << std::endl;}
176: buildQuadraticHexFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
177: if (debug > 1) {std::cout << " added face " << face << std::endl;}
178: faces[dim].push_back(oPoint_type(face, o));
179: }
180: } else if (dim > 1) {
181: const int edgeSizeQuadHex = 3;
182: const int numEdgesQuadHex = 4;
183: int nodes[12] = {
184: 0, 1, 4, // bottom
185: 1, 2, 5, // right
186: 2, 3, 6, // top
187: 3, 0, 7, // left
188: };
189: for(int b = 0; b < numEdgesQuadHex; b++) {
190: typename sieve_type::point_type face;
191: int o = 1;
193: bdVertices[dim-1].clear();
194: for(int c = 0; c < edgeSizeQuadHex; c++) {
195: bdVertices[dim-1].push_back(bdVertices[dim][nodes[b*edgeSizeQuadHex+c]]);
196: }
197: if (debug > 1) {
198: std::cout << " boundary point " << bdVertices[dim][b] << std::endl;
199: std::cout << " boundary vertices";
200: for(int c = 0; c < (int) bdVertices[dim-1].size(); c++) {
201: std::cout << " " << bdVertices[dim-1][c];
202: }
203: std::cout << std::endl;
204: }
205: buildQuadraticHexFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
206: if (debug > 1) {std::cout << " added face " << face << std::endl;}
207: faces[dim].push_back(oPoint_type(face, o));
208: }
209: } else {
210: if (debug > 1) {std::cout << " Just set faces to boundary in 1d" << std::endl;}
211: typename PointArray::iterator bdEnd = bdVertices[dim].end();
212: for (typename PointArray::iterator bd_iter = bdVertices[dim].begin(); bd_iter != bdEnd; ++bd_iter)
213: faces[dim].push_back(oPoint_type(*bd_iter, 0));
214: //faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
215: }
216: if (debug > 1) {
217: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
218: std::cout << " face point " << f_iter->first << " orientation " << f_iter->second << std::endl;
219: }
220: }
221: // We always create the toplevel, so we could short circuit somehow
222: // Should not have to loop here since the meet of just 2 boundary elements is an element
223: typename oPointArray::iterator f_itor = faces[dim].begin();
224: const typename sieve_type::point_type& start = f_itor->first;
225: const typename sieve_type::point_type& next = (++f_itor)->first;
226: Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);
228: if (preElement->size() > 0) {
229: cell = *preElement->begin();
231: const int size = faces[dim].size();
232: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
233: int wrap = size > 2 ? size-1 : 0;
234: int indA = 0, indB = 0;
236: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indA) {
237: if (start == *c_iter) break;
238: }
239: if (debug > 1) {std::cout << " pointA " << start << " indA " << indA << std::endl;}
240: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indB) {
241: if (next == *c_iter) break;
242: }
243: if (debug > 1) {std::cout << " pointB " << next << " indB " << indB << std::endl;}
244: if ((indB - indA == 1) || (indA - indB == wrap)) {
245: if (cellOrientation > 0) {
246: cellOrientation = indA+1;
247: } else {
248: if (dim == 1) {
249: cellOrientation = -2;
250: } else {
251: cellOrientation = -(indA+1);
252: }
253: }
254: } else if ((indA - indB == 1) || (indB - indA == wrap)) {
255: if (debug > 1) {std::cout << " reversing cell orientation" << std::endl;}
256: if (cellOrientation > 0) {
257: cellOrientation = -(indA+1);
258: } else {
259: if (dim == 1) {
260: cellOrientation = 1;
261: } else {
262: cellOrientation = indA+1;
263: }
264: }
265: } else {
266: throw ALE::Exception("Inconsistent orientation");
267: }
268: if (debug > 1) {std::cout << " Found old cell " << cell << " orientation " << cellOrientation << std::endl;}
269: } else {
270: int color = 0;
272: cell = typename sieve_type::point_type((*curElement[dim])++);
273: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
274: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(f_iter->first, cell);
276: sieve->addArrow(f_iter->first, cell, color++);
277: if (f_iter->second) {
278: orientation->addPoint(arrow);
279: orientation->updatePoint(arrow, &(f_iter->second));
280: if (debug > 1) {std::cout << " Orienting arrow (" << f_iter->first << ", " << cell << ") to " << f_iter->second << std::endl;}
281: }
282: }
283: if (cellOrientation > 0) {
284: cellOrientation = 1;
285: } else {
286: cellOrientation = -(dim+1);
287: }
288: if (debug > 1) {std::cout << " Added cell " << cell << " dim " << dim << std::endl;}
289: }
290: };
291: static void buildQuadraticTetFaces(Obj<sieve_type> sieve, Obj<arrow_section_type> orientation,
292: int dim,
293: std::map<int, int*>& curElement,
294: std::map<int,PointArray>& bdVertices,
295: std::map<int,oPointArray>& faces,
296: typename sieve_type::point_type& cell,
297: int& cellOrientation) {
298: int debug = sieve->debug();
300: if (debug > 1)
301: std::cout << " Building tet faces for cell " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
303: switch (dim) {
304: case 1: {
305: break;
306: } // case 1
307: case 2: {
308: assert(6 == bdVertices[dim].size());
310: // Edges on face
311: int color = 0;
312: int faceOrientation = 1;
314: // Bottom edge
315: typename sieve_type::point_type edge = bdVertices[dim][0];
316: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(edge, cell);
317: sieve->addArrow(edge, cell, color++);
318: orientation->addPoint(arrow);
319: if (bdVertices[dim][3] < bdVertices[dim][4])
320: faceOrientation = 1;
321: else
322: faceOrientation = -dim;
323: orientation->updatePoint(arrow, &faceOrientation);
324: if (debug > 1)
325: std::cout << " Adding bottom edge " << edge << " with orientation " << faceOrientation << std::endl;
327: // Right edge
328: faceOrientation = 1;
329: edge = bdVertices[dim][1];
330: arrow.source = edge;
331: sieve->addArrow(edge, cell, color++);
332: orientation->addPoint(arrow);
333: if (bdVertices[dim][4] < bdVertices[dim][5])
334: faceOrientation = 1;
335: else
336: faceOrientation = -dim;
337: orientation->updatePoint(arrow, &faceOrientation);
338: if (debug > 1)
339: std::cout << " Adding right edge " << edge << " with orientation " << faceOrientation << std::endl;
340:
341: // Left edge
342: faceOrientation = 1;
343: edge = bdVertices[dim][2];
344: arrow.source = edge;
345: sieve->addArrow(edge, cell, color++);
346: orientation->addPoint(arrow);
347: if (bdVertices[dim][5] < bdVertices[dim][3])
348: faceOrientation = 1;
349: else
350: faceOrientation = -dim;
351: orientation->updatePoint(arrow, &faceOrientation);
352: if (debug > 1)
353: std::cout << " Adding left edge " << edge << " with orientation " << faceOrientation << std::endl;
354:
355: // Vertices on edges
357: // Vertices on bottom edge
358: color = 0;
359: edge = bdVertices[dim][0];
360: typename sieve_type::point_type vertexA = bdVertices[dim][3];
361: typename sieve_type::point_type vertexB = bdVertices[dim][4];
362: sieve->addArrow(vertexA, edge, color++);
363: sieve->addArrow(vertexB, edge, color++);
364: if (debug > 1)
365: std::cout << " Adding vertices " << vertexA << " and " << vertexB << std::endl;
366:
367: // Vertices on right edge
368: color = 0;
369: edge = bdVertices[dim][1];
370: vertexA = bdVertices[dim][4];
371: vertexB = bdVertices[dim][5];
372: sieve->addArrow(vertexA, edge, color++);
373: sieve->addArrow(vertexB, edge, color++);
374: if (debug > 1)
375: std::cout << " Adding vertices " << vertexA << " and " << vertexB << std::endl;
376:
377: // Vertices on left edge
378: color = 0;
379: edge = bdVertices[dim][2];
380: vertexA = bdVertices[dim][5];
381: vertexB = bdVertices[dim][3];
382: sieve->addArrow(vertexA, edge, color++);
383: sieve->addArrow(vertexB, edge, color++);
384: if (debug > 1)
385: std::cout << " Adding vertices " << vertexA << " and " << vertexB << std::endl;
386:
387: break;
388: } // case 2
389: case 3: {
390: break;
391: } // case 3
392: default:
393: std::cerr << "Unknown dimension " << dim << std::endl;
394: assert(0);
395: throw ALE::Exception("Unknown dimension");
396: } // switch
398: };
399: static void buildFaces(Obj<sieve_type> sieve, Obj<arrow_section_type> orientation, int dim, std::map<int, int*>& curElement, std::map<int,PointArray>& bdVertices, std::map<int,oPointArray>& faces, typename sieve_type::point_type& cell, int& cellOrientation) {
400: int debug = sieve->debug();
402: if (debug > 1) {
403: if (cell >= 0) {
404: std::cout << " Building faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
405: } else {
406: std::cout << " Building faces for boundary of undetermined cell (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
407: }
408: }
409: if (dim == 0) return;
410: faces[dim].clear();
411: if (dim > 1) {
412: int b = 0;
413: // Use the cone construction
414: for(typename PointArray::iterator b_itor = bdVertices[dim].begin(); b_itor != bdVertices[dim].end(); ++b_itor, ++b) {
415: typename sieve_type::point_type face = -1;
416: int o = b%2 ? -bdVertices[dim].size() : 1;
418: bdVertices[dim-1].clear();
419: for(typename PointArray::iterator i_itor = bdVertices[dim].begin(); i_itor != bdVertices[dim].end(); ++i_itor) {
420: if (i_itor != b_itor) {
421: bdVertices[dim-1].push_back(*i_itor);
422: }
423: }
424: if (debug > 1) {std::cout << " boundary point " << *b_itor << std::endl;}
425: buildFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
426: if (debug > 1) {std::cout << " added face " << face << std::endl;}
427: faces[dim].push_back(oPoint_type(face, o));
428: }
429: } else {
430: if (debug > 1) {std::cout << " Just set faces to boundary in 1d" << std::endl;}
431: typename PointArray::iterator bd_iter = bdVertices[dim].begin();
432: faces[dim].push_back(oPoint_type(*bd_iter, 0));++bd_iter;
433: faces[dim].push_back(oPoint_type(*bd_iter, 0));
434: //faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
435: }
436: if (debug > 1) {
437: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
438: std::cout << " face point " << f_iter->first << " orientation " << f_iter->second << std::endl;
439: }
440: }
441: // We always create the toplevel, so we could short circuit somehow
442: // Should not have to loop here since the meet of just 2 boundary elements is an element
443: typename oPointArray::iterator f_itor = faces[dim].begin();
444: const typename sieve_type::point_type& start = f_itor->first;
445: const typename sieve_type::point_type& next = (++f_itor)->first;
446: Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);
448: if (preElement->size() > 0) {
449: cell = *preElement->begin();
451: const int size = faces[dim].size();
452: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
453: int wrap = size > 2 ? size-1 : 0;
454: int indA = 0, indB = 0;
456: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indA) {
457: if (start == *c_iter) break;
458: }
459: if (debug > 1) {std::cout << " pointA " << start << " indA " << indA << std::endl;}
460: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indB) {
461: if (next == *c_iter) break;
462: }
463: if (debug > 1) {std::cout << " pointB " << next << " indB " << indB << std::endl;}
464: if ((indB - indA == 1) || (indA - indB == wrap)) {
465: if (cellOrientation > 0) {
466: cellOrientation = indA+1;
467: } else {
468: if (dim == 1) {
469: cellOrientation = -2;
470: } else {
471: cellOrientation = -(indA+1);
472: }
473: }
474: } else if ((indA - indB == 1) || (indB - indA == wrap)) {
475: if (debug > 1) {std::cout << " reversing cell orientation" << std::endl;}
476: if (cellOrientation > 0) {
477: cellOrientation = -(indA+1);
478: } else {
479: if (dim == 1) {
480: cellOrientation = 1;
481: } else {
482: cellOrientation = indA+1;
483: }
484: }
485: } else {
486: throw ALE::Exception("Inconsistent orientation");
487: }
488: if (debug > 1) {std::cout << " Found old cell " << cell << " orientation " << cellOrientation << std::endl;}
489: } else {
490: int color = 0;
492: cell = typename sieve_type::point_type((*curElement[dim])++);
493: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
494: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(f_iter->first, cell);
496: sieve->addArrow(f_iter->first, cell, color++);
497: if (f_iter->second) {
498: orientation->addPoint(arrow);
499: orientation->updatePoint(arrow, &(f_iter->second));
500: if (debug > 1) {std::cout << " Orienting arrow (" << f_iter->first << ", " << cell << ") to " << f_iter->second << std::endl;}
501: }
502: }
503: if (cellOrientation > 0) {
504: cellOrientation = 1;
505: } else {
506: cellOrientation = -(dim+1);
507: }
508: if (debug > 1) {std::cout << " Added cell " << cell << " dim " << dim << " orientation " << cellOrientation << std::endl;}
509: }
510: };
511: #if 0
512: static void orientTriangle(const typename sieve_type::point_type cell, const int vertices[], const Obj<sieve_type>& sieve, const Obj<arrow_section_type>& orientation, int firstVertex[]) {
513: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
514: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
515: int debug = sieve->debug();
516: int corners = 3;
517: int e = 0;
519: if (debug > 1) {std::cout << "Orienting triangle " << cell << std::endl;}
520: for(typename sieve_type::traits::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter, ++e) {
521: if (debug > 1) {std::cout << " edge " << *p_iter << std::endl;}
522: const Obj<typename sieve_type::traits::coneSequence>& endpoints = sieve->cone(*p_iter);
523: typename sieve_type::traits::coneSequence::iterator vertex = endpoints->begin();
524: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(*p_iter, cell);
525: int indA, indB, value;
527: orientation->addPoint(arrow);
528: for(indA = 0; indA < corners; indA++) {if (*vertex == vertices[indA]) break;}
529: if (debug > 1) {std::cout << " vertexA " << *vertex << " indA " << indA <<std::endl;}
530: firstVertex[e] = *vertex;
531: ++vertex;
532: for(indB = 0; indB < corners; indB++) {if (*vertex == vertices[indB]) break;}
533: if (debug > 1) {std::cout << " vertexB " << *vertex << " indB " << indB <<std::endl;}
534: if ((indA == corners) || (indB == corners) || (indA == indB)) {throw ALE::Exception("Invalid edge endpoints");}
535: if ((indB - indA == 1) || (indA - indB == 2)) {
536: value = 1;
537: } else {
538: value = -1;
539: firstVertex[e] = *vertex;
540: }
541: if (debug > 1) {std::cout << " value " << value <<std::endl;}
542: orientation->updatePoint(arrow, &value);
543: }
544: };
545: #endif
548: // Build a topology from a connectivity description
549: // (0, 0) ... (0, numCells-1): dim-dimensional cells
550: // (0, numCells) ... (0, numVertices): vertices
551: // The other cells are numbered as they are requested
552: static void buildTopology(Obj<sieve_type> sieve, int dim, int numCells, int cells[], int numVertices, bool interpolate = true, int corners = -1, int firstVertex = -1, Obj<arrow_section_type> orientation = NULL, int firstCell = 0) {
553: if (interpolate && orientation.isNull()) {
554: throw ALE::Exception("Cannot interpolate mesh without providing an orientation Section");
555: }
556: ALE_LOG_EVENT_BEGIN;
557: if (sieve->commRank() != 0) {
558: ALE_LOG_EVENT_END;
559: return;
560: }
561: buildTopology_private(sieve, dim, numCells, cells, numVertices, interpolate, corners, firstVertex, orientation, firstCell);
562: ALE_LOG_EVENT_END;
563: };
564: static void buildTopology_private(Obj<sieve_type> sieve, int dim, int numCells, int cells[], int numVertices, bool interpolate = true, int corners = -1, int firstVertex = -1, Obj<arrow_section_type> orientation = NULL, int firstCell = 0) {
565: int debug = sieve->debug();
567: if (firstVertex < 0) firstVertex = numCells;
568: // Create a map from dimension to the current element number for that dimension
569: std::map<int,int*> curElement;
570: std::map<int,PointArray> bdVertices;
571: std::map<int,PointArray> faces;
572: std::map<int,oPointArray> oFaces;
573: int curCell = firstCell;
574: int curVertex = firstVertex;
575: int newElement = firstVertex > firstCell ? firstVertex + numVertices : firstCell + numCells;
576: int o = 1;
578: if (corners < 0) corners = dim+1;
579: curElement[0] = &curVertex;
580: curElement[dim] = &curCell;
581: for(int d = 1; d < dim; d++) {
582: curElement[d] = &newElement;
583: }
584: for(int c = 0; c < numCells; c++) {
585: typename sieve_type::point_type cell(c);
587: // Build the cell
588: if (interpolate) {
589: bdVertices[dim].clear();
590: for(int b = 0; b < corners; b++) {
591: // This ordering produces the same vertex order as the uninterpolated mesh
592: //typename sieve_type::point_type vertex(cells[c*corners+(b+corners-1)%corners]+firstVertex);
593: typename sieve_type::point_type vertex(cells[c*corners+b]+firstVertex);
595: if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
596: bdVertices[dim].push_back(vertex);
597: }
598: if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}
600: if ((2 == dim && 4 == corners) || (3 == dim && 8 == corners)) {
601: buildHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
602: } else if ((2 == dim && 9 == corners) || (3 == dim && 27 == corners)) {
603: buildQuadraticHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
604: } else if ((2 == dim && 6 == corners) || (3 == dim && 10 == corners)) {
605: buildQuadraticTetFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
606: } else {
607: buildFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
608: }
609: #if 0
610: if ((dim == 2) && (!orientation.isNull())) {
611: typename sieve_type::point_type eVertices[3];
612: typename sieve_type::point_type fVertices[3];
614: for(int v = 0; v < 3; ++v) {
615: fVertices[v] = cells[c*corners+v]+numCells;
616: }
617: orientTriangle(cell, fVertices, sieve, orientation, eVertices);
618: } else if ((dim == 3) && (!orientation.isNull())) {
619: // The order of vertices in cells[] orients the cell itself
620: if (debug > 1) {std::cout << "Orienting tetrahedron " << cell << std::endl;}
621: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
622: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
623: int f = 0;
625: for(typename sieve_type::traits::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter, ++f) {
626: if (debug > 1) {std::cout << " face " << *p_iter << std::endl;}
627: const Obj<typename sieve_type::traits::coneSequence>& fCone = sieve->cone(*p_iter);
628: const typename sieve_type::traits::coneSequence::iterator fEnd = fCone->end();
629: typename sieve_type::point_type fVertices[3];
630: typename sieve_type::point_type eVertices[3];
632: // We will choose the orientation such that the normals are outward
633: for(int v = 0, i = 0; v < corners; ++v) {
634: if (v == f) continue;
635: fVertices[i++] = cells[c*corners+v]+numCells;
636: }
637: if (f%2) {
638: int tmp = fVertices[0];
639: fVertices[0] = fVertices[1];
640: fVertices[1] = tmp;
641: }
642: orientTriangle(*p_iter, fVertices, sieve, orientation, eVertices);
643: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> fArrow(*p_iter, cell);
644: int indC, indD, indE, value;
646: orientation->addPoint(fArrow);
647: for(indC = 0; indC < corners; indC++) {if (eVertices[0] == fVertices[indC]) break;}
648: if (debug > 1) {std::cout << " vertexC " << eVertices[0] << " indC " << indC <<std::endl;}
649: for(indD = 0; indD < corners; indD++) {if (eVertices[1] == fVertices[indD]) break;}
650: if (debug > 1) {std::cout << " vertexD " << eVertices[1] << " indD " << indD <<std::endl;}
651: for(indE = 0; indE < corners; indE++) {if (eVertices[2] == fVertices[indE]) break;}
652: if (debug > 1) {std::cout << " vertexE " << eVertices[2] << " indE " << indE <<std::endl;}
653: if ((indC == corners) || (indD == corners) || (indE == corners) ||
654: (indC == indD) || (indD == indE) || (indE == indC)) {throw ALE::Exception("Invalid face corners");}
655: if ((indD - indC == 1) || (indC - indD == 2)) {
656: if (!((indE - indD == 1) || (indD - indE == 2)) ||
657: !((indC - indE == 1) || (indE - indC == 2))) {throw ALE::Exception("Invalid order");}
658: value = 1;
659: } else {
660: value = -1;
661: }
662: if (debug > 1) {std::cout << " value " << value <<std::endl;}
663: orientation->updatePoint(fArrow, &value);
664: orientation->view("Intermediate orientation");
665: }
666: }
667: #endif
668: } else {
669: for(int b = 0; b < corners; b++) {
670: sieve->addArrow(typename sieve_type::point_type(cells[c*corners+b]+firstVertex), cell, b);
671: }
672: if (debug) {
673: if (debug > 1) {
674: for(int b = 0; b < corners; b++) {
675: std::cout << " Adding vertex " << typename sieve_type::point_type(cells[c*corners+b]+firstVertex) << std::endl;
676: }
677: }
678: if ((numCells < 10000) || (c%1000 == 0)) {
679: std::cout << "Adding cell " << cell << " dim " << dim << std::endl;
680: }
681: }
682: }
683: }
684: };
685: static void buildCoordinates(const Obj<Bundle_>& bundle, const int embedDim, const PetscReal coords[], int numCells = -1) {
686: const Obj<typename Bundle_::real_section_type>& coordinates = bundle->getRealSection("coordinates");
687: const Obj<typename Bundle_::label_sequence>& vertices = bundle->depthStratum(0);
688: const int debug = bundle->debug();
690: if (numCells < 0) {
691: numCells = bundle->heightStratum(0)->size();
692: }
693: bundle->setupCoordinates(coordinates);
694: coordinates->setFiberDimension(vertices, embedDim);
695: bundle->allocate(coordinates);
696: for(typename Bundle_::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
697: coordinates->updatePoint(*v_iter, &(coords[(*v_iter - numCells)*embedDim]));
698: if (debug) {
699: if ((numCells < 10000) || ((*v_iter)%1000 == 0)) {
700: std::cout << "Adding coordinates for vertex " << *v_iter << std::endl;
701: }
702: }
703: }
704: };
707: // Build a topology from a connectivity description
708: // (0, 0) ... (0, numCells-1): dim-dimensional cells
709: // (0, numCells) ... (0, numVertices): vertices
710: // The other cells are numbered as they are requested
711: static void buildTopologyMultiple(Obj<sieve_type> sieve, int dim, int numCells, int cells[], int numVertices, bool interpolate = true, int corners = -1, int firstVertex = -1, Obj<arrow_section_type> orientation = NULL) {
712: int debug = sieve->debug();
714: ALE_LOG_EVENT_BEGIN;
715: int *cellOffset = new int[sieve->commSize()+1];
716: cellOffset[0] = 0;
717: MPI_Allgather(&numCells, 1, MPI_INT, &cellOffset[1], 1, MPI_INT, sieve->comm());
718: for(int p = 1; p <= sieve->commSize(); ++p) cellOffset[p] += cellOffset[p-1];
719: int *vertexOffset = new int[sieve->commSize()+1];
720: vertexOffset[0] = 0;
721: MPI_Allgather(&numVertices, 1, MPI_INT, &vertexOffset[1], 1, MPI_INT, sieve->comm());
722: for(int p = 1; p <= sieve->commSize(); ++p) vertexOffset[p] += vertexOffset[p-1];
723: if (firstVertex < 0) firstVertex = cellOffset[sieve->commSize()] + vertexOffset[sieve->commRank()];
724: // Estimate the number of intermediates as (V+C)*(dim-1)
725: // Should include a check for running over the namespace
726: // Create a map from dimension to the current element number for that dimension
727: std::map<int,int*> curElement;
728: std::map<int,PointArray> bdVertices;
729: std::map<int,PointArray> faces;
730: std::map<int,oPointArray> oFaces;
731: int curCell = cellOffset[sieve->commRank()];
732: int curVertex = firstVertex;
733: int newElement = firstVertex+vertexOffset[sieve->commSize()] + (cellOffset[sieve->commRank()] + vertexOffset[sieve->commRank()])*(dim-1);
734: int o = 1;
736: if (corners < 0) corners = dim+1;
737: curElement[0] = &curVertex;
738: curElement[dim] = &curCell;
739: for(int d = 1; d < dim; d++) {
740: curElement[d] = &newElement;
741: }
742: for(int c = 0; c < numCells; c++) {
743: typename sieve_type::point_type cell(c);
745: // Build the cell
746: if (interpolate) {
747: bdVertices[dim].clear();
748: for(int b = 0; b < corners; b++) {
749: // This ordering produces the same vertex order as the uninterpolated mesh
750: //typename sieve_type::point_type vertex(cells[c*corners+(b+corners-1)%corners]+firstVertex);
751: typename sieve_type::point_type vertex(cells[c*corners+b]+firstVertex);
753: if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
754: bdVertices[dim].push_back(vertex);
755: }
756: if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}
758: if ((2 == dim && 4 == corners) || (3 == dim && 8 == corners)) {
759: buildHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
760: } else if ((2 == dim && 9 == corners) || (3 == dim && 27 == corners)) {
761: buildQuadraticHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
762: } else if ((2 == dim && 6 == corners) || (3 == dim && 10 == corners)) {
763: buildQuadraticTetFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
764: } else {
765: buildFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
766: }
767: } else {
768: for(int b = 0; b < corners; b++) {
769: sieve->addArrow(typename sieve_type::point_type(cells[c*corners+b]+firstVertex), cell, b);
770: }
771: if (debug) {
772: if (debug > 1) {
773: for(int b = 0; b < corners; b++) {
774: std::cout << " Adding vertex " << typename sieve_type::point_type(cells[c*corners+b]+firstVertex) << std::endl;
775: }
776: }
777: if ((numCells < 10000) || (c%1000 == 0)) {
778: std::cout << "Adding cell " << cell << " dim " << dim << std::endl;
779: }
780: }
781: }
782: }
784: if (newElement >= firstVertex+vertexOffset[sieve->commSize()] + (cellOffset[sieve->commRank()+1] + vertexOffset[sieve->commRank()+1])*(dim-1)) {
785: throw ALE::Exception("Namespace violation during intermediate element construction");
786: }
787: delete [] cellOffset;
788: delete [] vertexOffset;
789: ALE_LOG_EVENT_END;
790: };
791: static void buildCoordinatesMultiple(const Obj<Bundle_>& bundle, const int embedDim, const double coords[], int numGlobalCells = -1) {
792: const Obj<typename Bundle_::real_section_type>& coordinates = bundle->getRealSection("coordinates");
793: const Obj<typename Bundle_::label_sequence>& vertices = bundle->depthStratum(0);
794: const int numCells = bundle->heightStratum(0)->size(), numVertices = vertices->size();
795: const int debug = bundle->debug();
796: int offset;
798: if (numGlobalCells < 0) {
799: MPI_Allreduce((void *) &numCells, &numGlobalCells, 1, MPI_INT, MPI_SUM, bundle->comm());
800: }
801: MPI_Scan((void *) &numVertices, &offset, 1, MPI_INT, MPI_SUM, bundle->comm());
802: offset += numGlobalCells - numVertices;
803: coordinates->setFiberDimension(vertices, embedDim);
804: bundle->allocate(coordinates);
805: for(typename Bundle_::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
806: coordinates->updatePoint(*v_iter, &(coords[(*v_iter - offset)*embedDim]));
807: if (debug) {
808: if ((numCells < 10000) || ((*v_iter)%1000 == 0)) {
809: std::cout << "Adding coordinates for vertex " << *v_iter << std::endl;
810: }
811: }
812: }
813: };
814: };
815: }
816: #endif