Actual source code: SieveBuilder.hh
1: #ifndef included_ALE_SieveBuilder_hh
2: #define included_ALE_SieveBuilder_hh
4: #ifndef included_ALE_ALE_hh
5: #include <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 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) {
149: int debug = sieve->debug();
151: if (debug > 1) {
152: if (cell >= 0) {
153: std::cout << " Building faces for boundary of " << cell << " (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
154: } else {
155: std::cout << " Building faces for boundary of undetermined cell (size " << bdVertices[dim].size() << "), dim " << dim << std::endl;
156: }
157: }
158: if (dim == 0) return;
159: faces[dim].clear();
160: if (dim > 1) {
161: int b = 0;
162: // Use the cone construction
163: for(typename PointArray::iterator b_itor = bdVertices[dim].begin(); b_itor != bdVertices[dim].end(); ++b_itor, ++b) {
164: typename sieve_type::point_type face = -1;
165: int o = b%2 ? -bdVertices[dim].size() : 1;
167: bdVertices[dim-1].clear();
168: for(typename PointArray::iterator i_itor = bdVertices[dim].begin(); i_itor != bdVertices[dim].end(); ++i_itor) {
169: if (i_itor != b_itor) {
170: bdVertices[dim-1].push_back(*i_itor);
171: }
172: }
173: if (debug > 1) {std::cout << " boundary point " << *b_itor << std::endl;}
174: buildFaces(sieve, orientation, dim-1, curElement, bdVertices, faces, face, o);
175: if (debug > 1) {std::cout << " added face " << face << std::endl;}
176: faces[dim].push_back(oPoint_type(face, o));
177: }
178: } else {
179: if (debug > 1) {std::cout << " Just set faces to boundary in 1d" << std::endl;}
180: typename PointArray::iterator bd_iter = bdVertices[dim].begin();
181: faces[dim].push_back(oPoint_type(*bd_iter, 0));++bd_iter;
182: faces[dim].push_back(oPoint_type(*bd_iter, 0));
183: //faces[dim].insert(faces[dim].end(), bdVertices[dim].begin(), bdVertices[dim].end());
184: }
185: if (debug > 1) {
186: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
187: std::cout << " face point " << f_iter->first << " orientation " << f_iter->second << std::endl;
188: }
189: }
190: // We always create the toplevel, so we could short circuit somehow
191: // Should not have to loop here since the meet of just 2 boundary elements is an element
192: typename oPointArray::iterator f_itor = faces[dim].begin();
193: const typename sieve_type::point_type& start = f_itor->first;
194: const typename sieve_type::point_type& next = (++f_itor)->first;
195: Obj<typename sieve_type::supportSet> preElement = sieve->nJoin(start, next, 1);
197: if (preElement->size() > 0) {
198: cell = *preElement->begin();
200: const int size = faces[dim].size();
201: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
202: int wrap = size > 2 ? size-1 : 0;
203: int indA = 0, indB = 0;
205: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indA) {
206: if (start == *c_iter) break;
207: }
208: if (debug > 1) {std::cout << " pointA " << start << " indA " << indA << std::endl;}
209: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++indB) {
210: if (next == *c_iter) break;
211: }
212: if (debug > 1) {std::cout << " pointB " << next << " indB " << indB << std::endl;}
213: if ((indB - indA == 1) || (indA - indB == wrap)) {
214: if (cellOrientation > 0) {
215: cellOrientation = indA+1;
216: } else {
217: if (dim == 1) {
218: cellOrientation = -2;
219: } else {
220: cellOrientation = -(indA+1);
221: }
222: }
223: } else if ((indA - indB == 1) || (indB - indA == wrap)) {
224: if (debug > 1) {std::cout << " reversing cell orientation" << std::endl;}
225: if (cellOrientation > 0) {
226: cellOrientation = -(indA+1);
227: } else {
228: if (dim == 1) {
229: cellOrientation = 1;
230: } else {
231: cellOrientation = indA+1;
232: }
233: }
234: } else {
235: throw ALE::Exception("Inconsistent orientation");
236: }
237: if (debug > 1) {std::cout << " Found old cell " << cell << " orientation " << cellOrientation << std::endl;}
238: } else {
239: int color = 0;
241: cell = typename sieve_type::point_type((*curElement[dim])++);
242: for(typename oPointArray::iterator f_iter = faces[dim].begin(); f_iter != faces[dim].end(); ++f_iter) {
243: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(f_iter->first, cell);
245: sieve->addArrow(f_iter->first, cell, color++);
246: if (f_iter->second) {
247: orientation->addPoint(arrow);
248: orientation->updatePoint(arrow, &(f_iter->second));
249: if (debug > 1) {std::cout << " Orienting arrow (" << f_iter->first << ", " << cell << ") to " << f_iter->second << std::endl;}
250: }
251: }
252: if (cellOrientation > 0) {
253: cellOrientation = 1;
254: } else {
255: cellOrientation = -(dim+1);
256: }
257: if (debug > 1) {std::cout << " Added cell " << cell << " dim " << dim << " orientation " << cellOrientation << std::endl;}
258: }
259: };
260: #if 0
261: 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[]) {
262: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
263: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
264: int debug = sieve->debug();
265: int corners = 3;
266: int e = 0;
268: if (debug > 1) {std::cout << "Orienting triangle " << cell << std::endl;}
269: for(typename sieve_type::traits::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter, ++e) {
270: if (debug > 1) {std::cout << " edge " << *p_iter << std::endl;}
271: const Obj<typename sieve_type::traits::coneSequence>& endpoints = sieve->cone(*p_iter);
272: typename sieve_type::traits::coneSequence::iterator vertex = endpoints->begin();
273: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> arrow(*p_iter, cell);
274: int indA, indB, value;
276: orientation->addPoint(arrow);
277: for(indA = 0; indA < corners; indA++) {if (*vertex == vertices[indA]) break;}
278: if (debug > 1) {std::cout << " vertexA " << *vertex << " indA " << indA <<std::endl;}
279: firstVertex[e] = *vertex;
280: ++vertex;
281: for(indB = 0; indB < corners; indB++) {if (*vertex == vertices[indB]) break;}
282: if (debug > 1) {std::cout << " vertexB " << *vertex << " indB " << indB <<std::endl;}
283: if ((indA == corners) || (indB == corners) || (indA == indB)) {throw ALE::Exception("Invalid edge endpoints");}
284: if ((indB - indA == 1) || (indA - indB == 2)) {
285: value = 1;
286: } else {
287: value = -1;
288: firstVertex[e] = *vertex;
289: }
290: if (debug > 1) {std::cout << " value " << value <<std::endl;}
291: orientation->updatePoint(arrow, &value);
292: }
293: };
294: #endif
297: // Build a topology from a connectivity description
298: // (0, 0) ... (0, numCells-1): dim-dimensional cells
299: // (0, numCells) ... (0, numVertices): vertices
300: // The other cells are numbered as they are requested
301: 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) {
302: ALE_LOG_EVENT_BEGIN;
303: if (sieve->commRank() != 0) {
304: ALE_LOG_EVENT_END;
305: return;
306: }
307: buildTopology_private(sieve, dim, numCells, cells, numVertices, interpolate, corners, firstVertex, orientation, firstCell);
308: ALE_LOG_EVENT_END;
309: };
310: 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) {
311: int debug = sieve->debug();
313: if (firstVertex < 0) firstVertex = numCells;
314: // Create a map from dimension to the current element number for that dimension
315: std::map<int,int*> curElement;
316: std::map<int,PointArray> bdVertices;
317: std::map<int,PointArray> faces;
318: std::map<int,oPointArray> oFaces;
319: int curCell = firstCell;
320: int curVertex = firstVertex;
321: int newElement = firstVertex > firstCell ? firstVertex + numVertices : firstCell + numCells;
322: int o = 1;
324: if (corners < 0) corners = dim+1;
325: curElement[0] = &curVertex;
326: curElement[dim] = &curCell;
327: for(int d = 1; d < dim; d++) {
328: curElement[d] = &newElement;
329: }
330: for(int c = 0; c < numCells; c++) {
331: typename sieve_type::point_type cell(c);
333: // Build the cell
334: if (interpolate) {
335: bdVertices[dim].clear();
336: for(int b = 0; b < corners; b++) {
337: // This ordering produces the same vertex order as the uninterpolated mesh
338: //typename sieve_type::point_type vertex(cells[c*corners+(b+corners-1)%corners]+firstVertex);
339: typename sieve_type::point_type vertex(cells[c*corners+b]+firstVertex);
341: if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
342: bdVertices[dim].push_back(vertex);
343: }
344: if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}
346: if (corners != dim+1) {
347: buildHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
348: } else {
349: buildFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
350: }
351: #if 0
352: if ((dim == 2) && (!orientation.isNull())) {
353: typename sieve_type::point_type eVertices[3];
354: typename sieve_type::point_type fVertices[3];
356: for(int v = 0; v < 3; ++v) {
357: fVertices[v] = cells[c*corners+v]+numCells;
358: }
359: orientTriangle(cell, fVertices, sieve, orientation, eVertices);
360: } else if ((dim == 3) && (!orientation.isNull())) {
361: // The order of vertices in cells[] orients the cell itself
362: if (debug > 1) {std::cout << "Orienting tetrahedron " << cell << std::endl;}
363: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(cell);
364: const typename sieve_type::traits::coneSequence::iterator end = cone->end();
365: int f = 0;
367: for(typename sieve_type::traits::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter, ++f) {
368: if (debug > 1) {std::cout << " face " << *p_iter << std::endl;}
369: const Obj<typename sieve_type::traits::coneSequence>& fCone = sieve->cone(*p_iter);
370: const typename sieve_type::traits::coneSequence::iterator fEnd = fCone->end();
371: typename sieve_type::point_type fVertices[3];
372: typename sieve_type::point_type eVertices[3];
374: // We will choose the orientation such that the normals are outward
375: for(int v = 0, i = 0; v < corners; ++v) {
376: if (v == f) continue;
377: fVertices[i++] = cells[c*corners+v]+numCells;
378: }
379: if (f%2) {
380: int tmp = fVertices[0];
381: fVertices[0] = fVertices[1];
382: fVertices[1] = tmp;
383: }
384: orientTriangle(*p_iter, fVertices, sieve, orientation, eVertices);
385: MinimalArrow<typename sieve_type::point_type,typename sieve_type::point_type> fArrow(*p_iter, cell);
386: int indC, indD, indE, value;
388: orientation->addPoint(fArrow);
389: for(indC = 0; indC < corners; indC++) {if (eVertices[0] == fVertices[indC]) break;}
390: if (debug > 1) {std::cout << " vertexC " << eVertices[0] << " indC " << indC <<std::endl;}
391: for(indD = 0; indD < corners; indD++) {if (eVertices[1] == fVertices[indD]) break;}
392: if (debug > 1) {std::cout << " vertexD " << eVertices[1] << " indD " << indD <<std::endl;}
393: for(indE = 0; indE < corners; indE++) {if (eVertices[2] == fVertices[indE]) break;}
394: if (debug > 1) {std::cout << " vertexE " << eVertices[2] << " indE " << indE <<std::endl;}
395: if ((indC == corners) || (indD == corners) || (indE == corners) ||
396: (indC == indD) || (indD == indE) || (indE == indC)) {throw ALE::Exception("Invalid face corners");}
397: if ((indD - indC == 1) || (indC - indD == 2)) {
398: if (!((indE - indD == 1) || (indD - indE == 2)) ||
399: !((indC - indE == 1) || (indE - indC == 2))) {throw ALE::Exception("Invalid order");}
400: value = 1;
401: } else {
402: value = -1;
403: }
404: if (debug > 1) {std::cout << " value " << value <<std::endl;}
405: orientation->updatePoint(fArrow, &value);
406: orientation->view("Intermediate orientation");
407: }
408: }
409: #endif
410: } else {
411: for(int b = 0; b < corners; b++) {
412: sieve->addArrow(typename sieve_type::point_type(cells[c*corners+b]+firstVertex), cell, b);
413: }
414: if (debug) {
415: if (debug > 1) {
416: for(int b = 0; b < corners; b++) {
417: std::cout << " Adding vertex " << typename sieve_type::point_type(cells[c*corners+b]+firstVertex) << std::endl;
418: }
419: }
420: if ((numCells < 10000) || (c%1000 == 0)) {
421: std::cout << "Adding cell " << cell << " dim " << dim << std::endl;
422: }
423: }
424: }
425: }
426: };
427: static void buildCoordinates(const Obj<Bundle_>& bundle, const int embedDim, const double coords[]) {
428: const Obj<typename Bundle_::real_section_type>& coordinates = bundle->getRealSection("coordinates");
429: const Obj<typename Bundle_::label_sequence>& vertices = bundle->depthStratum(0);
430: const int numCells = bundle->heightStratum(0)->size();
431: const int debug = bundle->debug();
433: bundle->setupCoordinates(coordinates);
434: coordinates->setFiberDimension(vertices, embedDim);
435: bundle->allocate(coordinates);
436: for(typename Bundle_::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
437: coordinates->updatePoint(*v_iter, &(coords[(*v_iter - numCells)*embedDim]));
438: if (debug) {
439: if ((numCells < 10000) || ((*v_iter)%1000 == 0)) {
440: std::cout << "Adding coordinates for vertex " << *v_iter << std::endl;
441: }
442: }
443: }
444: };
447: // Build a topology from a connectivity description
448: // (0, 0) ... (0, numCells-1): dim-dimensional cells
449: // (0, numCells) ... (0, numVertices): vertices
450: // The other cells are numbered as they are requested
451: 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) {
452: int debug = sieve->debug();
454: ALE_LOG_EVENT_BEGIN;
455: int *cellOffset = new int[sieve->commSize()+1];
456: cellOffset[0] = 0;
457: MPI_Allgather(&numCells, 1, MPI_INT, &cellOffset[1], 1, MPI_INT, sieve->comm());
458: for(int p = 1; p <= sieve->commSize(); ++p) cellOffset[p] += cellOffset[p-1];
459: int *vertexOffset = new int[sieve->commSize()+1];
460: vertexOffset[0] = 0;
461: MPI_Allgather(&numVertices, 1, MPI_INT, &vertexOffset[1], 1, MPI_INT, sieve->comm());
462: for(int p = 1; p <= sieve->commSize(); ++p) vertexOffset[p] += vertexOffset[p-1];
463: if (firstVertex < 0) firstVertex = cellOffset[sieve->commSize()] + vertexOffset[sieve->commRank()];
464: // Estimate the number of intermediates as (V+C)*(dim-1)
465: // Should include a check for running over the namespace
466: // Create a map from dimension to the current element number for that dimension
467: std::map<int,int*> curElement;
468: std::map<int,PointArray> bdVertices;
469: std::map<int,PointArray> faces;
470: std::map<int,oPointArray> oFaces;
471: int curCell = cellOffset[sieve->commRank()];
472: int curVertex = firstVertex;
473: int newElement = firstVertex+vertexOffset[sieve->commSize()] + (cellOffset[sieve->commRank()] + vertexOffset[sieve->commRank()])*(dim-1);
474: int o = 1;
476: if (corners < 0) corners = dim+1;
477: curElement[0] = &curVertex;
478: curElement[dim] = &curCell;
479: for(int d = 1; d < dim; d++) {
480: curElement[d] = &newElement;
481: }
482: for(int c = 0; c < numCells; c++) {
483: typename sieve_type::point_type cell(c);
485: // Build the cell
486: if (interpolate) {
487: bdVertices[dim].clear();
488: for(int b = 0; b < corners; b++) {
489: // This ordering produces the same vertex order as the uninterpolated mesh
490: //typename sieve_type::point_type vertex(cells[c*corners+(b+corners-1)%corners]+firstVertex);
491: typename sieve_type::point_type vertex(cells[c*corners+b]+firstVertex);
493: if (debug > 1) {std::cout << "Adding boundary vertex " << vertex << std::endl;}
494: bdVertices[dim].push_back(vertex);
495: }
496: if (debug) {std::cout << "cell " << cell << " num boundary vertices " << bdVertices[dim].size() << std::endl;}
498: if (corners != dim+1) {
499: buildHexFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
500: } else {
501: buildFaces(sieve, orientation, dim, curElement, bdVertices, oFaces, cell, o);
502: }
503: } else {
504: for(int b = 0; b < corners; b++) {
505: sieve->addArrow(typename sieve_type::point_type(cells[c*corners+b]+firstVertex), cell, b);
506: }
507: if (debug) {
508: if (debug > 1) {
509: for(int b = 0; b < corners; b++) {
510: std::cout << " Adding vertex " << typename sieve_type::point_type(cells[c*corners+b]+firstVertex) << std::endl;
511: }
512: }
513: if ((numCells < 10000) || (c%1000 == 0)) {
514: std::cout << "Adding cell " << cell << " dim " << dim << std::endl;
515: }
516: }
517: }
518: }
520: if (newElement >= firstVertex+vertexOffset[sieve->commSize()] + (cellOffset[sieve->commRank()+1] + vertexOffset[sieve->commRank()+1])*(dim-1)) {
521: throw ALE::Exception("Namespace violation during intermediate element construction");
522: }
523: delete [] cellOffset;
524: delete [] vertexOffset;
525: ALE_LOG_EVENT_END;
526: };
527: static void buildCoordinatesMultiple(const Obj<Bundle_>& bundle, const int embedDim, const double coords[]) {
528: const Obj<typename Bundle_::real_section_type>& coordinates = bundle->getRealSection("coordinates");
529: const Obj<typename Bundle_::label_sequence>& vertices = bundle->depthStratum(0);
530: const int numCells = bundle->heightStratum(0)->size(), numVertices = vertices->size();
531: const int debug = bundle->debug();
532: int numGlobalCells, offset;
534: MPI_Allreduce((void *) &numCells, &numGlobalCells, 1, MPI_INT, MPI_SUM, bundle->comm());
535: MPI_Scan((void *) &numVertices, &offset, 1, MPI_INT, MPI_SUM, bundle->comm());
536: offset += numGlobalCells - numVertices;
537: coordinates->setFiberDimension(vertices, embedDim);
538: bundle->allocate(coordinates);
539: for(typename Bundle_::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
540: coordinates->updatePoint(*v_iter, &(coords[(*v_iter - offset)*embedDim]));
541: if (debug) {
542: if ((numCells < 10000) || ((*v_iter)%1000 == 0)) {
543: std::cout << "Adding coordinates for vertex " << *v_iter << std::endl;
544: }
545: }
546: }
547: };
548: };
549: }
550: #endif