Actual source code: CartesianSieve.hh
petsc-3.3-p7 2013-05-11
1: #ifndef included_ALE_CartesianSieve_hh
2: #define included_ALE_CartesianSieve_hh
4: #ifndef included_ALE_Mesh_hh
5: #include <sieve/Mesh.hh>
6: #endif
8: namespace ALE {
9: namespace CartesianSieveDef {
10: static void nextCodeOdd(int a[], const int n, const int size) {
11: for(int d = size-2; d >= 0; --d) a[d+1] = a[d];
12: if (n == (n/2)*2) {
13: a[0] = 1;
14: } else {
15: a[0] = 0;
16: }
17: };
18: static void nextCodeEven(int a[], const int n, const int size) {
19: ALE::CartesianSieveDef::nextCodeOdd(a, n, size);
20: a[0] ^= 1;
21: for(int d = 1; d < size; ++d) a[d] ^= 0;
22: };
23: static void getGrayCode(const int n, const int size, int code[]) {
24: if (n == 1 << size) {
25: code[0] = -1;
26: } else if (n == 0) {
27: for(int d = 0; d < size; ++d) code[d] = 0;
28: } else if (n == 1) {
29: code[0] = 1;
30: for(int d = 1; d < size; ++d) code[d] = 0;
31: } else if (n == (n/2)*2) {
32: // Might have to shift pointer
33: ALE::CartesianSieveDef::getGrayCode(n/2, size, code);
34: ALE::CartesianSieveDef::nextCodeEven(code, n/2, size);
35: } else {
36: // Might have to shift pointer
37: ALE::CartesianSieveDef::getGrayCode((n-1)/2, size, code);
38: ALE::CartesianSieveDef::nextCodeOdd(code, (n-1)/2, size);
39: }
40: };
41: static void getIndices(const int dim, const int sizes[], const int& point, int indices[]) {
42: int p = point;
43: int divisor = 1;
45: for(int d = 0; d < dim-1; ++d) {
46: divisor *= sizes[d];
47: }
48: //std::cout << " Getting indices for point " << point << ":";
49: for(int d = dim-1; d >= 0; --d) {
50: indices[d] = p/divisor;
51: p -= divisor*indices[d];
52: if (d > 0) divisor /= sizes[d-1];
53: }
54: //for(int d = 0; d < dim; ++d) {
55: // std::cout << " " << indices[d];
56: //}
57: //std::cout << std::endl;
58: };
59: static inline int getValue(const int dim, const int sizes[], const int indices[], const int code[], const bool forward = true) {
60: //std::cout << "Building value" << std::endl << " code:";
61: //for(int d = 0; d < dim; ++d) std::cout << " " << code[d];
62: //std::cout << std::endl << " indices:";
63: //for(int d = 0; d < dim; ++d) std::cout << " " << indices[d];
64: //std::cout << std::endl;
65: if (code[0] < 0) {
66: //std::cout << " value " << -1 << std::endl;
67: return -1;
68: }
69: int value = indices[dim-1];
70: if (forward) {
71: if (code[dim-1]) value++;
72: } else {
73: if (!code[dim-1]) value--;
74: }
75: //std::cout << " value " << value << std::endl;
76: for(int d = dim-2; d >= 0; --d) {
77: value = value*sizes[d] + indices[d];
78: if (forward) {
79: if (code[d]) value++;
80: } else {
81: if (!code[d]) value--;
82: }
83: //std::cout << " value " << value << std::endl;
84: }
85: //std::cout << " final value " << value << std::endl;
86: return value;
87: };
90: template <typename PointType_>
91: class PointSequence {
92: public:
93: typedef PointType_ point_type;
94: // We traverse the points in natural order
95: class iterator {
96: public:
97: // Standard iterator typedefs
98: typedef std::input_iterator_tag iterator_category;
99: typedef point_type value_type;
100: typedef int difference_type;
101: typedef value_type* pointer;
102: typedef value_type& reference;
103: protected:
104: value_type _value;
105: public:
106: iterator(const value_type& start) : _value(start) {};
107: iterator(const iterator& iter) : _value(iter._value) {};
108: virtual ~iterator() {};
109: public:
110: virtual bool operator==(const iterator& iter) const {return this->_value == iter._value;};
111: virtual bool operator!=(const iterator& iter) const {return this->_value != iter._value;};
112: virtual const value_type operator*() const {return this->_value;};
113: virtual iterator& operator++() {++this->_value; return *this;};
114: virtual iterator operator++(int n) {iterator tmp(*this); ++this->_value; return tmp;};
115: };
116: protected:
117: point_type _start;
118: point_type _end;
119: public:
120: PointSequence(const point_type& start, const point_type& end) : _start(start), _end(end) {};
121: virtual ~PointSequence() {};
122: public:
123: virtual iterator begin() {return iterator(this->_start);};
124: virtual iterator end() {return iterator(this->_end);};
125: virtual bool empty() {return this->_start >= this->_end;};
126: };
127: // In one dimension, we have
128: //
129: // v_k---c_k---v_{k+1}
130: //
131: // In two dimensions, we have
132: //
133: // c_{k+m} c_{k+m+1}
134: //
135: // v_{k+m+1}-----------v_{k+m+2}--
136: // | |
137: // | c_k | c_{k+1}
138: // | |
139: // v_k---------------v_{k+1}----
140: //
141: // So for c_k the cone is (v_{k+(k/m)}, v_{k+(k/m)+1}, v_{k+m+(k/m)}, v_{k+m+(k/m)+1})
142: //
143: // In three dimensions, we have
144: //
145: // v_{k+m+1+(m+1)*(n+1)} v_{k+m+1+(m+1)*(n+1)+1}
146: //
147: // v_{k+m+1} v_{k+m+2}
148: //
149: //
150: // v_{k+(m+1)*(n+1)} v_{k+(m+1)*(n+1)+1}
151: //
152: // v_k v_{k+1}
153: //
154: // Suppose we break down the index k into a tuple (i,j), then 2d becomes
155: //
156: // cone(c_{i,j}) = (v_{i,j}, v_{i+1,j}, v_{i,j+1}, v_{i+1,j+1})
157: // = (v_{k}, v_{k+1}, v_{k+m+1}, v_{k+m+2})
158: // Now for 2D
159: // i = k%m
160: // j = k/m
161: // and 3D
162: // k = q/(m*n)
163: // j = (q - m*n*k)/m
164: // i = (q - m*n*k - m*j)
165: template <typename PointType_>
166: class ConeSequence {
167: public:
168: typedef PointType_ point_type;
169: // We traverse the points in natural order
170: class iterator {
171: public:
172: // Standard iterator typedefs
173: typedef std::input_iterator_tag iterator_category;
174: typedef PointType_ value_type;
175: typedef int difference_type;
176: typedef value_type* pointer;
177: typedef value_type& reference;
178: protected:
179: int _dim;
180: const int *_sizes;
181: const int *_cellSizes;
182: int _numCells;
183: int _pos;
184: int *_code;
185: value_type _value;
186: value_type *_indices;
187: protected:
188: void init() {
189: this->_code = new int[this->_dim];
190: if (this->_pos == -2) {
191: this->_value = this->_indices[0];
192: ++this->_pos;
193: } else if (this->_pos == -1) {
194: this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_cellSizes, this->_indices, this->_code);
195: } else {
196: ALE::CartesianSieveDef::getGrayCode(this->_pos, this->_dim, this->_code);
197: this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code) + this->_numCells;
198: }
199: };
200: public:
201: iterator(const int dim, const int sizes[], const int cellSizes[], const int numCells, const value_type indices[], const int pos, const bool addSelf = false) : _dim(dim), _sizes(sizes), _cellSizes(cellSizes), _numCells(numCells), _pos(pos) {
202: this->_indices = new int[this->_dim];
203: for(int d = 0; d < this->_dim; ++d) this->_indices[d] = indices[d];
204: this->init();
205: };
206: iterator(const iterator& iter) : _dim(iter._dim), _sizes(iter._sizes), _cellSizes(iter._cellSizes), _numCells(iter._numCells), _pos(iter._pos) {
207: this->_indices = new int[this->_dim];
208: for(int d = 0; d < this->_dim; ++d) this->_indices[d] = iter._indices[d];
209: this->init();
210: };
211: virtual ~iterator() {
212: delete [] this->_code;
213: delete [] this->_indices;
214: };
215: public:
216: virtual bool operator==(const iterator& iter) const {return this->_pos == iter._pos;};
217: virtual bool operator!=(const iterator& iter) const {return this->_pos != iter._pos;};
218: virtual const value_type operator*() const {return this->_value;};
219: virtual iterator& operator++() {
220: ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
221: this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code) + this->_numCells;
222: return *this;
223: };
224: virtual iterator operator++(int n) {
225: iterator tmp(*this);
226: ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
227: this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code) + this->_numCells;
228: return tmp;
229: };
230: };
231: protected:
232: int _dim;
233: const int *_sizes;
234: int *_vertexSizes;
235: int _numCells;
236: point_type _cell;
237: int *_indices;
238: bool _addSelf;
239: public:
240: ConeSequence(const int dim, const int sizes[], const int numCells, const point_type& cell, const bool addSelf = false) : _dim(dim), _sizes(sizes), _numCells(numCells), _cell(cell), _addSelf(addSelf) {
241: this->_vertexSizes = new int[dim];
242: this->_indices = new point_type[dim];
243: //this->getIndices(this->_cell, this->_indices);
244: for(int d = 0; d < this->_dim; ++d) this->_vertexSizes[d] = this->_sizes[d]+1;
245: ALE::CartesianSieveDef::getIndices(this->_dim, this->_sizes, this->_cell, this->_indices);
246: };
247: virtual ~ConeSequence() {
248: delete [] this->_vertexSizes;
249: delete [] this->_indices;
250: };
251: protected:
252: void getIndices(const point_type& cell, const int indices[]) {
253: point_type c = cell;
254: int divisor = 1;
256: for(int d = 0; d < this->_dim-1; ++d) {
257: divisor *= this->_sizes[d];
258: }
259: //std::cout << " Got indices for cell " << cell << ":";
260: for(int d = this->_dim-1; d >= 0; --d) {
261: this->_indices[d] = c/divisor;
262: c -= divisor*this->_indices[d];
263: if (d > 0) divisor /= this->_sizes[d-1];
264: }
265: //for(int d = 0; d < this->_dim; ++d) {
266: // std::cout << " " << this->_indices[d];
267: //}
268: //std::cout << std::endl;
269: };
270: public:
271: virtual iterator begin() {
272: int start = 0;
274: if (this->_addSelf) {
275: if (this->_cell >= this->_numCells) {
276: start = -2;
277: this->_indices[0] = this->_cell;
278: } else {
279: start = -1;
280: }
281: }
282: return iterator(this->_dim, this->_vertexSizes, this->_sizes, this->_numCells, this->_indices, start);
283: };
284: virtual iterator end() {
285: int end = 1 << this->_dim;
287: if (!this->_dim || (this->_cell >= this->_numCells)) {
288: end = 0;
289: }
290: return iterator(this->_dim, this->_vertexSizes, this->_sizes, this->_numCells, this->_indices, end);
291: };
292: };
293: template <typename PointType_>
294: class SupportSequence {
295: public:
296: typedef PointType_ point_type;
297: // We traverse the points in natural order
298: class iterator {
299: public:
300: // Standard iterator typedefs
301: typedef std::input_iterator_tag iterator_category;
302: typedef PointType_ value_type;
303: typedef int difference_type;
304: typedef value_type* pointer;
305: typedef value_type& reference;
306: protected:
307: int _dim;
308: const int *_sizes;
309: int _pos;
310: int *_code;
311: value_type _value;
312: value_type *_indices;
313: protected:
314: bool validIndex(const int indices[], const int code[]) {
315: for(int d = 0; d < this->_dim; ++d) {
316: if ((code[d]) && (indices[d] >= this->_sizes[d])) return false;
317: if ((!code[d]) && (indices[d] < 1)) return false;
318: }
319: return true;
320: }
321: void init() {
322: this->_code = new int[this->_dim];
323: ALE::CartesianSieveDef::getGrayCode(this->_pos, this->_dim, this->_code);
324: while((this->_code[0] >= 0) && !this->validIndex(this->_indices, this->_code)) {
325: ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
326: }
327: this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code, false);
328: };
329: public:
330: iterator(const int dim, const int sizes[], const value_type indices[], const int pos) : _dim(dim), _sizes(sizes), _pos(pos) {
331: this->_indices = new int[this->_dim];
332: for(int d = 0; d < this->_dim; ++d) this->_indices[d] = indices[d];
333: this->init();
334: };
335: iterator(const iterator& iter) : _dim(iter._dim), _sizes(iter._sizes), _pos(iter._pos) {
336: this->_indices = new int[this->_dim];
337: for(int d = 0; d < this->_dim; ++d) this->_indices[d] = iter._indices[d];
338: this->init();
339: };
340: virtual ~iterator() {
341: delete [] this->_code;
342: delete [] this->_indices;
343: };
344: public:
345: virtual bool operator==(const iterator& iter) const {return this->_pos == iter._pos;};
346: virtual bool operator!=(const iterator& iter) const {return this->_pos != iter._pos;};
347: virtual const value_type operator*() const {return this->_value;};
348: virtual iterator& operator++() {
349: do {
350: ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
351: } while((this->_code[0] >= 0) && !this->validIndex(this->_indices, this->_code));
352: this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code, false);
353: return *this;
354: };
355: virtual iterator operator++(int n) {
356: iterator tmp(*this);
357: do {
358: ALE::CartesianSieveDef::getGrayCode(++this->_pos, this->_dim, this->_code);
359: } while((this->_code[0] >= 0) && !this->validIndex(this->_indices, this->_code));
360: this->_value = ALE::CartesianSieveDef::getValue(this->_dim, this->_sizes, this->_indices, this->_code, false);
361: return tmp;
362: };
363: };
364: protected:
365: int _dim;
366: const int *_sizes;
367: int *_vertexSizes;
368: int _numCells;
369: point_type _vertex;
370: int *_indices;
371: public:
372: SupportSequence(const int dim, const int sizes[], const int numCells, const point_type& vertex) : _dim(dim), _sizes(sizes), _numCells(numCells), _vertex(vertex) {
373: this->_vertexSizes = new int[dim];
374: this->_indices = new point_type[dim];
375: //this->getIndices(this->_vertex, this->_indices);
376: for(int d = 0; d < this->_dim; ++d) this->_vertexSizes[d] = this->_sizes[d]+1;
377: ALE::CartesianSieveDef::getIndices(this->_dim, this->_vertexSizes, this->_vertex - this->_numCells, this->_indices);
378: };
379: virtual ~SupportSequence() {
380: delete [] this->_vertexSizes;
381: delete [] this->_indices;
382: };
383: protected:
384: void getIndices(const point_type& vertex, const int indices[]) {
385: point_type v = vertex - this->_numCells;
386: int divisor = 1;
388: for(int d = 0; d < this->_dim-1; ++d) {
389: divisor *= this->_sizes[d]+1;
390: }
391: std::cout << " Got indices for vertex " << vertex << ":";
392: for(int d = this->_dim-1; d >= 0; --d) {
393: this->_indices[d] = v/divisor;
394: v -= divisor*this->_indices[d];
395: if (d > 0) divisor /= this->_sizes[d-1]+1;
396: }
397: for(int d = 0; d < this->_dim; ++d) {
398: std::cout << " " << this->_indices[d];
399: }
400: std::cout << std::endl;
401: };
402: public:
403: virtual iterator begin() {return iterator(this->_dim, this->_sizes, this->_indices, 0);};
404: virtual iterator end() {return iterator(this->_dim, this->_sizes, this->_indices, this->_dim ? 1 << this->_dim: 0);};
405: };
406: }
407: // We can do meets of two cells as empty if they do not intersect, and a k-dimensional face (with 2^k vertices)
408: // if they have k indices in common
409: // We can do joins of two vertices. If they have have k indices in common, they span a d-k face, and thus have
410: // 2^{k} cells in the join.
411: // Note meet and join are both empty if any index differs by more than 1
413: // This is a purely Cartesian mesh
414: // We will only represent cells and vertices, meaning the cap consists only of
415: // vertices, and the base only of cells.
416: // numCells: m x n x ...
417: // numVertices: m+1 x n+1 x ...
418: // Points:
419: // cells: 0 - numCells-1
420: // vertices: numCells - numCells+numVertices-1
421: template <typename Point_>
422: class CartesianSieve : public ALE::ParallelObject {
423: public:
424: typedef Point_ point_type;
425: typedef CartesianSieveDef::PointSequence<point_type> baseSequence;
426: typedef CartesianSieveDef::PointSequence<point_type> capSequence;
427: typedef CartesianSieveDef::ConeSequence<point_type> coneSequence;
428: typedef CartesianSieveDef::SupportSequence<point_type> supportSequence;
429: // Backward compatibility
430: typedef coneSequence coneArray;
431: typedef coneSequence coneSet;
432: typedef supportSequence supportArray;
433: typedef supportSequence supportSet;
434: protected:
435: int _dim;
436: int *_sizes;
437: int _numCells;
438: int _numVertices;
439: public:
440: CartesianSieve(MPI_Comm comm, const int dim, const int numCells[], const int& debug = 0) : ParallelObject(comm, debug), _dim(dim) {
441: this->_sizes = new int[dim];
442: this->_numCells = 1;
443: this->_numVertices = 1;
444: for(int d = 0; d < dim; ++d) {
445: this->_sizes[d] = numCells[d];
446: this->_numCells *= numCells[d];
447: this->_numVertices *= numCells[d]+1;
448: }
449: };
450: virtual ~CartesianSieve() {
451: delete [] this->_sizes;
452: };
453: public:
454: int getDimension() const {return this->_dim;};
455: const int *getSizes() const {return this->_sizes;};
456: int getNumCells() const {return this->_numCells;};
457: int getNumVertices() const {return this->_numVertices;};
458: public:
459: // WARNING: Creating a lot of objects here
460: Obj<baseSequence> base() {
461: Obj<baseSequence> base = new baseSequence(0, this->_numCells);
462: return base;
463: };
464: Obj<baseSequence> leaves() {
465: return this->base();
466: };
467: // WARNING: Creating a lot of objects here
468: Obj<capSequence> cap() {
469: Obj<capSequence> cap = new capSequence(this->_numCells, this->_numCells+this->_numVertices);
470: return cap;
471: };
472: Obj<capSequence> roots() {
473: return this->cap();
474: };
475: // WARNING: Creating a lot of objects here
476: Obj<coneSequence> cone(const point_type& p) {
477: if ((p < 0) || (p >= this->_numCells)) {
478: Obj<coneSequence> cone = new coneSequence(0, this->_sizes, this->_numCells, 0);
479: return cone;
480: }
481: Obj<coneSequence> cone = new coneSequence(this->_dim, this->_sizes, this->_numCells, p);
482: return cone;
483: };
484: // WARNING: Creating a lot of objects here
485: Obj<supportSequence> support(const point_type& p) {
486: if ((p < this->_numCells) || (p >= this->_numCells+this->_numVertices)) {
487: Obj<supportSequence> support = new supportSequence(0, this->_sizes, this->_numCells, 0);
488: return support;
489: }
490: Obj<supportSequence> support = new supportSequence(this->_dim, this->_sizes, this->_numCells, p);
491: return support;
492: };
493: // WARNING: Creating a lot of objects here
494: Obj<coneSequence> closure(const point_type& p) {
495: Obj<coneSequence> cone = new coneSequence(this->_dim, this->_sizes, this->_numCells, p, true);
496: return cone;
497: };
498: // WARNING: Creating a lot of objects here
499: Obj<supportSequence> star(const point_type& p) {return this->support(p);};
500: // WARNING: Creating a lot of objects here
501: Obj<coneSequence> nCone(const point_type& p, const int n) {
502: if (n == 0) return this->cone(-1);
503: if (n != 1) throw ALE::Exception("Invalid height for nCone");
504: return this->cone(p);
505: };
506: // WARNING: Creating a lot of objects here
507: Obj<supportSequence> nSupport(const point_type& p, const int n) {
508: if (n == 0) return this->support(-1);
509: if (n != 1) throw ALE::Exception("Invalid depth for nSupport");
510: return this->support(p);
511: };
512: public:
513: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
514: ostringstream txt;
516: if (comm == MPI_COMM_NULL) {
517: comm = this->comm();
518: }
519: if (name == "") {
520: PetscPrintf(comm, "viewing a CartesianSieve\n");
521: } else {
522: PetscPrintf(comm, "viewing CartesianSieve '%s'\n", name.c_str());
523: }
524: if(this->commRank() == 0) {
525: txt << "cap --> base:\n";
526: }
527: Obj<capSequence> cap = this->cap();
528: Obj<baseSequence> base = this->base();
529: if(cap->empty()) {
530: txt << "[" << this->commRank() << "]: empty" << std::endl;
531: }
532: for(typename capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
533: const Obj<supportSequence>& supp = this->support(*capi);
534: const typename supportSequence::iterator suppEnd = supp->end();
536: for(typename supportSequence::iterator suppi = supp->begin(); suppi != suppEnd; ++suppi) {
537: txt << "[" << this->commRank() << "]: " << *capi << "---->" << *suppi << std::endl;
538: }
539: }
540: PetscSynchronizedPrintf(this->comm(), txt.str().c_str());
541: PetscSynchronizedFlush(this->comm());
542: //
543: ostringstream txt1;
544: if(this->commRank() == 0) {
545: txt1 << "base --> cap:\n";
546: }
547: if(base->empty()) {
548: txt1 << "[" << this->commRank() << "]: empty" << std::endl;
549: }
550: for(typename baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
551: const Obj<coneSequence>& cone = this->cone(*basei);
552: const typename coneSequence::iterator coneEnd = cone->end();
554: for(typename coneSequence::iterator conei = cone->begin(); conei != coneEnd; ++conei) {
555: txt1 << "[" << this->commRank() << "]: " << *basei << "<----" << *conei << std::endl;
556: }
557: }
558: //
559: PetscSynchronizedPrintf(this->comm(), txt1.str().c_str());
560: PetscSynchronizedFlush(this->comm());
561: //
562: ostringstream txt2;
563: if(this->commRank() == 0) {
564: txt2 << "cap <point>:\n";
565: }
566: txt2 << "[" << this->commRank() << "]: [";
567: for(typename capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
568: txt2 << " <" << *capi << ">";
569: }
570: txt2 << " ]" << std::endl;
571: //
572: PetscSynchronizedPrintf(this->comm(), txt2.str().c_str());
573: PetscSynchronizedFlush(this->comm());
574: //
575: ostringstream txt3;
576: if(this->commRank() == 0) {
577: txt3 << "base <point>:\n";
578: }
579: txt3 << "[" << this->commRank() << "]: [";
580: for(typename baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
581: txt3 << " <" << *basei << ">";
582: }
583: txt3 << " ]" << std::endl;
584: //
585: PetscSynchronizedPrintf(this->comm(), txt3.str().c_str());
586: PetscSynchronizedFlush(this->comm());
587: };
588: };
590: // We do not just use Bundle, because we need to optimize labels
591: class CartesianBundle : public ALE::ParallelObject {
592: public:
593: typedef CartesianSieve<int> sieve_type;
594: typedef sieve_type::point_type point_type;
595: typedef Section<point_type, double> real_section_type;
596: typedef Section<point_type, int> int_section_type;
597: typedef UniformSection<MinimalArrow<point_type, point_type>, int> arrow_section_type;
598: typedef std::map<std::string, Obj<arrow_section_type> > arrow_sections_type;
599: typedef std::map<std::string, Obj<real_section_type> > real_sections_type;
600: typedef std::map<std::string, Obj<int_section_type> > int_sections_type;
601: typedef sieve_type::baseSequence label_sequence;
602: typedef ALE::Sifter<int,point_type,point_type> send_overlap_type;
603: typedef ALE::Sifter<point_type,int,point_type> recv_overlap_type;
604: protected:
605: Obj<sieve_type> _sieve;
606: int _maxHeight;
607: int _maxDepth;
608: arrow_sections_type _arrowSections;
609: real_sections_type _realSections;
610: int_sections_type _intSections;
611: Obj<send_overlap_type> _sendOverlap;
612: Obj<recv_overlap_type> _recvOverlap;
613: Obj<send_overlap_type> _distSendOverlap;
614: Obj<recv_overlap_type> _distRecvOverlap;
615: public:
616: CartesianBundle(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1) {
617: this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
618: this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
619: };
620: virtual ~CartesianBundle() {};
621: public: // Accessors
622: void clear() {
623: this->_maxHeight = -1;
624: this->_maxDepth = -1;
625: };
626: const Obj<sieve_type>& getSieve() const {return this->_sieve;};
627: void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
628: bool hasArrowSection(const std::string& name) const {
629: return this->_arrowSections.find(name) != this->_arrowSections.end();
630: };
631: const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
632: if (!this->hasArrowSection(name)) {
633: Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());
635: section->setName(name);
636: if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
637: this->_arrowSections[name] = section;
638: }
639: return this->_arrowSections[name];
640: };
641: void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
642: this->_arrowSections[name] = section;
643: };
644: Obj<std::set<std::string> > getArrowSections() const {
645: Obj<std::set<std::string> > names = std::set<std::string>();
647: for(arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
648: names->insert(s_iter->first);
649: }
650: return names;
651: };
652: bool hasRealSection(const std::string& name) const {
653: return this->_realSections.find(name) != this->_realSections.end();
654: };
655: const Obj<real_section_type>& getRealSection(const std::string& name) {
656: if (!this->hasRealSection(name)) {
657: Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());
659: section->setName(name);
660: if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
661: this->_realSections[name] = section;
662: }
663: return this->_realSections[name];
664: };
665: void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
666: this->_realSections[name] = section;
667: };
668: Obj<std::set<std::string> > getRealSections() const {
669: Obj<std::set<std::string> > names = std::set<std::string>();
671: for(real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
672: names->insert(s_iter->first);
673: }
674: return names;
675: };
676: bool hasIntSection(const std::string& name) const {
677: return this->_intSections.find(name) != this->_intSections.end();
678: };
679: const Obj<int_section_type>& getIntSection(const std::string& name) {
680: if (!this->hasIntSection(name)) {
681: Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());
683: section->setName(name);
684: if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
685: this->_intSections[name] = section;
686: }
687: return this->_intSections[name];
688: };
689: void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
690: this->_intSections[name] = section;
691: };
692: Obj<std::set<std::string> > getIntSections() const {
693: Obj<std::set<std::string> > names = std::set<std::string>();
695: for(int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
696: names->insert(s_iter->first);
697: }
698: return names;
699: };
700: const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
701: void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
702: const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
703: void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
704: const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
705: void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
706: const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
707: void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
708: public: // Stratification
711: virtual void stratify() {
712: ALE_LOG_EVENT_BEGIN;
713: this->_maxHeight = 1;
714: this->_maxDepth = 1;
715: ALE_LOG_EVENT_END;
716: };
717: virtual int height() const {return this->_maxHeight;};
718: virtual int height(const point_type& point) {
719: const int numCells = this->getSieve()->getNumCells();
720: const int numVertices = this->getSieve()->getNumVertices();
721: if ((point >= 0) && (point < numCells)) return 0;
722: if ((point >= numCells) && (point < numCells+numVertices)) return 1;
723: return -1;
724: };
725: // WARNING: Creating a lot of objects here
726: virtual const Obj<label_sequence> heightStratum(int height) {
727: if (height == 0) return this->getSieve()->base();
728: if (height == 1) return this->getSieve()->cap();
729: throw ALE::Exception("Invalid height stratum");
730: };
731: virtual int depth() const {return this->_maxDepth;};
732: virtual int depth(const point_type& point) {
733: const int numCells = this->getSieve()->getNumCells();
734: const int numVertices = this->getSieve()->getNumVertices();
735: if ((point >= 0) && (point < numCells)) return 1;
736: if ((point >= numCells) && (point < numCells+numVertices)) return 0;
737: return -1;
738: };
739: // WARNING: Creating a lot of objects here
740: virtual const Obj<label_sequence> depthStratum(int depth) {
741: if (depth == 0) return this->getSieve()->cap();
742: if (depth == 1) return this->getSieve()->base();
743: throw ALE::Exception("Invalid depth stratum");
744: };
745: virtual const Obj<label_sequence> getLabelStratum(const std::string& name, int value) {
746: if (name == "height") {
747: return this->heightStratum(value);
748: } else if (name == "depth") {
749: return this->depthStratum(value);
750: }
751: throw ALE::Exception("Invalid label name");
752: }
753: public: // Viewers
754: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
755: if (comm == MPI_COMM_NULL) {
756: comm = this->comm();
757: }
758: if (name == "") {
759: PetscPrintf(comm, "viewing a CartesianBundle\n");
760: } else {
761: PetscPrintf(comm, "viewing CartesianBundle '%s'\n", name.c_str());
762: }
763: PetscPrintf(comm, " maximum height %d maximum depth %d\n", this->height(), this->depth());
764: };
765: public:
766: void constructOverlap() {};
767: };
769: class CartesianMesh : public CartesianBundle {
770: public:
771: typedef int point_type;
772: typedef CartesianBundle::sieve_type sieve_type;
773: typedef CartesianBundle::label_sequence label_sequence;
774: typedef CartesianBundle::send_overlap_type send_overlap_type;
775: typedef CartesianBundle::recv_overlap_type recv_overlap_type;
776: typedef sieve_type::coneSequence coneSequence;
777: typedef sieve_type::supportSequence supportSequence;
778: protected:
779: int _dim;
780: // Discretization
781: Obj<Discretization> _discretization;
782: Obj<BoundaryCondition> _boundaryCondition;
783: public:
784: CartesianMesh(MPI_Comm comm, int dim, int debug = 0) : CartesianBundle(comm, debug), _dim(dim) {
785: this->_discretization = new Discretization(comm, debug);
786: this->_boundaryCondition = new BoundaryCondition(comm, debug);
787: };
788: virtual ~CartesianMesh() {};
789: public: // Accessors
790: int getDimension() {return this->_dim;};
791: const Obj<Discretization>& getDiscretization() {return this->_discretization;};
792: void setDiscretization(const Obj<Discretization>& discretization) {this->_discretization = discretization;};
793: const Obj<BoundaryCondition>& getBoundaryCondition() {return this->_boundaryCondition;};
794: void setBoundaryCondition(const Obj<BoundaryCondition>& boundaryCondition) {this->_boundaryCondition = boundaryCondition;};
795: public: // Discretization
796: #if 0
797: void markBoundaryCells(const std::string& name) {
798: const Obj<label_type>& label = this->getLabel(name);
799: const Obj<label_sequence>& boundary = this->getLabelStratum(name, 1);
800: const Obj<sieve_type>& sieve = this->getSieve();
802: for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
803: if (this->height(*e_iter) == 1) {
804: const point_type cell = *sieve->support(*e_iter)->begin();
806: this->setValue(label, cell, 2);
807: }
808: }
809: };
810: #endif
811: void setupField(const Obj<real_section_type>& s, const bool postponeGhosts = false) {
812: const std::string& name = this->_boundaryCondition->getLabelName();
814: for(int d = 0; d <= this->depth(); ++d) {
815: s->setFiberDimension(this->depthStratum(d), this->_discretization->getNumDof(d));
816: }
817: if (!name.empty()) {
818: const Obj<label_sequence>& boundary = this->getLabelStratum(name, 1);
820: for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
821: s->setFiberDimension(*e_iter, -this->_discretization->getNumDof(this->depth(*e_iter)));
822: }
823: }
824: s->allocatePoint();
825: if (!name.empty()) {
826: #if 0
827: const Obj<label_sequence>& boundaryCells = this->getLabelStratum(name, 2);
828: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
829: real_section_type::value_type *values = new real_section_type::value_type[this->sizeWithBC(s, *boundaryCells->begin())];
830: double *v0 = new double[this->getDimension()];
831: double *J = new double[this->getDimension()*this->getDimension()];
832: double detJ;
834: for(label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
835: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
836: const coneArray::iterator end = closure->end();
837: int v = 0;
839: this->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
840: for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
841: const int cDim = s->getConstraintDimension(*cl_iter);
843: if (cDim) {
844: for(int d = 0; d < cDim; ++d, ++v) {
845: values[v] = this->_boundaryCondition->integrateDual(v0, J, v);
846: }
847: } else {
848: const int dim = s->getFiberDimension(*cl_iter);
850: for(int d = 0; d < dim; ++d, ++v) {
851: values[v] = 0.0;
852: }
853: }
854: }
855: this->updateAll(s, *c_iter, values);
856: }
857: delete [] values;
858: delete [] v0;
859: delete [] J;
860: #endif
861: }
862: };
863: public: // Size traversal
864: template<typename Section_>
865: int size(const Obj<Section_>& section, const point_type& p) {
866: const typename Section_::chart_type& chart = section->getChart();
867: const Obj<coneSequence> cone = this->getSieve()->cone(p);
868: typename coneSequence::iterator end = cone->end();
869: int size = 0;
871: if (chart.count(p)) {
872: size += section->getConstrainedFiberDimension(p);
873: }
874: for(typename coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
875: if (chart.count(*c_iter)) {
876: size += section->getConstrainedFiberDimension(*c_iter);
877: }
878: }
879: return size;
880: }
881: template<typename Section_>
882: int sizeWithBC(const Obj<Section_>& section, const point_type& p) {
883: const typename Section_::chart_type& chart = section->getChart();
884: const Obj<coneSequence> cone = this->getSieve()->cone(p);
885: typename coneSequence::iterator end = cone->end();
886: int size = 0;
888: if (chart.count(p)) {
889: size += section->getFiberDimension(p);
890: }
891: for(typename coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
892: if (chart.count(*c_iter)) {
893: size += section->getFiberDimension(*c_iter);
894: }
895: }
896: return size;
897: }
898: public: // Allocation
899: template<typename Section_>
900: void allocate(const Obj<Section_>& section) {
901: section->allocatePoint();
902: }
903: public: // Retrieval traversal
904: // Return the values for the closure of this point
905: // use a smart pointer?
906: template<typename Section_>
907: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p) {
908: const int size = this->sizeWithBC(section, p);
909: typename Section_::value_type *values = section->getRawArray(size);
910: int j = -1;
912: const int& dim = section->getFiberDimension(p);
913: const typename Section_::value_type *array = section->restrictPoint(p);
915: for(int i = 0; i < dim; ++i) {
916: values[++j] = array[i];
917: }
918: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
919: typename sieve_type::coneSequence::iterator end = cone->end();
921: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
922: const int& dim = section->getFiberDimension(*p_iter);
924: array = section->restrictPoint(*p_iter);
925: for(int i = 0; i < dim; ++i) {
926: values[++j] = array[i];
927: }
928: }
929: if (j != size-1) {
930: ostringstream txt;
932: txt << "Invalid restrictClosure to point " << p << std::endl;
933: txt << " j " << j << " should be " << (size-1) << std::endl;
934: std::cout << txt.str();
935: throw ALE::Exception(txt.str().c_str());
936: }
937: return values;
938: }
939: template<typename Section_>
940: void update(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
941: int j = 0;
943: section->updatePoint(p, &v[j]);
944: j += section->getFiberDimension(p);
945: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
947: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
948: section->updatePoint(*p_iter, &v[j]);
949: j += section->getFiberDimension(*p_iter);
950: }
951: }
952: template<typename Section_>
953: void updateAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
954: int j = 0;
956: section->updateAddPoint(p, &v[j]);
957: j += section->getFiberDimension(p);
958: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
960: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
961: section->updateAddPoint(*p_iter, &v[j]);
962: j += section->getFiberDimension(*p_iter);
963: }
964: }
965: template<typename Section_>
966: void updateBC(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
967: int j = 0;
969: section->updatePointBC(p, &v[j]);
970: j += section->getFiberDimension(p);
971: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
973: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
974: section->updatePointBC(*p_iter, &v[j]);
975: j += section->getFiberDimension(*p_iter);
976: }
977: }
978: template<typename Section_>
979: void updateAll(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
980: int j = 0;
982: section->updatePointBC(p, &v[j]);
983: j += section->getFiberDimension(p);
984: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
986: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
987: section->updatePointAll(*p_iter, &v[j]);
988: j += section->getFiberDimension(*p_iter);
989: }
990: }
991: public: // Mesh geometry
992: void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
993: const double *coords = this->restrictClosure(coordinates, e);
994: const int dim = this->_dim;
995: const int last = 1 << dim;
997: if (v0) {
998: for(int d = 0; d < dim; d++) {
999: v0[d] = coords[d];
1000: }
1001: }
1002: for(int d = 0; d < dim; ++d) {
1003: if (J) {
1004: J[d] = 0.5*(coords[last*dim+d] - v0[d]);
1005: }
1006: if (invJ) {
1007: invJ[d] = 1.0/J[d];
1008: }
1009: }
1010: }
1011: public: // Viewers
1012: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1013: if (comm == MPI_COMM_NULL) {
1014: comm = this->comm();
1015: }
1016: if (name == "") {
1017: PetscPrintf(comm, "viewing a CartesianMesh\n");
1018: } else {
1019: PetscPrintf(comm, "viewing CartesianMesh '%s'\n", name.c_str());
1020: }
1021: this->getSieve()->view("mesh sieve");
1022: }
1023: };
1025: class CartesianMeshBuilder {
1026: public:
1027: static Obj<CartesianMesh> createCartesianMesh(const MPI_Comm comm, const int dim, const int numCells[], const int partitions[], const int debug = 0) {
1029: // Steal PETSc code that calculates partitions
1030: // We could conceivably allow multiple patches per partition
1031: int size, rank;
1032: int totalPartitions = 1;
1035: MPI_Comm_size(comm, &size);CHKERRXX(ierr);
1036: MPI_Comm_rank(comm, &rank);CHKERRXX(ierr);
1037: for(int d = 0; d < dim; ++d) totalPartitions *= partitions[d];
1038: if (size != totalPartitions) throw ALE::Exception("Invalid partitioning");
1039: // Determine local sizes
1040: int *numLocalCells = new int[dim];
1041: int *numNeighborVertices = new int[dim];
1042: int *numLocalVertices = new int[dim];
1044: for(int d = 0; d < dim; ++d) {
1045: numLocalCells[d] = numCells[d]/partitions[d] + (rank < numCells[d]%partitions[d]);
1046: numLocalVertices[d] = numLocalCells[d]+1;
1047: }
1048: // Create mesh
1049: const Obj<CartesianMesh> mesh = new CartesianMesh(comm, dim, debug);
1050: const Obj<CartesianMesh::sieve_type> sieve = new CartesianMesh::sieve_type(comm, dim, numLocalCells, debug);
1052: mesh->setSieve(sieve);
1053: mesh->stratify();
1054: delete [] numLocalCells;
1055: // Create overlaps
1056: // We overlap anyone within 1 index of us for the entire boundary
1057: const Obj<CartesianMesh::send_overlap_type>& sendOverlap = mesh->getSendOverlap();
1058: const Obj<CartesianMesh::recv_overlap_type>& recvOverlap = mesh->getRecvOverlap();
1059: int *indices = new int[dim];
1060: int *code = new int[dim];
1061: int *zeroCode = new int[dim];
1062: int *vIndices = new int[dim];
1063: int numNeighborCells;
1065: ALE::CartesianSieveDef::getIndices(dim, partitions, rank, indices);
1066: for(int e = 0; e < dim; ++e) zeroCode[e] = 0;
1067: for(int d = 0; d < dim; ++d) {
1068: if (indices[d] > 0) {
1069: for(int e = 0; e < dim; ++e) code[e] = 1;
1070: code[d] = 0;
1071: int neighborRank = ALE::CartesianSieveDef::getValue(dim, partitions, indices, code, false);
1072: numNeighborCells = 1;
1073: for(int e = 0; e < dim; ++e) {
1074: numNeighborCells *= (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]));
1075: numNeighborVertices[e] = (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]))+1;
1076: }
1078: // Add the whole d-1 face on the left edge of dimension d
1079: for(int v = sieve->getNumCells(); v < sieve->getNumCells()+sieve->getNumVertices(); ++v) {
1080: ALE::CartesianSieveDef::getIndices(dim, numLocalVertices, v - sieve->getNumCells(), vIndices);
1081: if (vIndices[d] == 0) {
1082: vIndices[d] = numNeighborVertices[d]-1;
1083: int neighborV = ALE::CartesianSieveDef::getValue(dim, numNeighborVertices, vIndices, zeroCode) + numNeighborCells;
1085: sendOverlap->addCone(v, neighborRank, neighborV);
1086: recvOverlap->addCone(neighborRank, v, neighborV);
1087: }
1088: }
1089: }
1090: if (indices[d] < partitions[d]-1) {
1091: for(int e = 0; e < dim; ++e) code[e] = 0;
1092: code[d] = 1;
1093: int neighborRank = ALE::CartesianSieveDef::getValue(dim, partitions, indices, code);
1094: numNeighborCells = 1;
1095: for(int e = 0; e < dim; ++e) {
1096: numNeighborCells *= (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]));
1097: numNeighborVertices[e] = (numCells[e]/partitions[e] + (neighborRank < numCells[e]%partitions[e]))+1;
1098: }
1100: // Add the whole d-1 face on the right edge of dimension d
1101: for(int v = sieve->getNumCells(); v < sieve->getNumCells()+sieve->getNumVertices(); ++v) {
1102: ALE::CartesianSieveDef::getIndices(dim, numLocalVertices, v - sieve->getNumCells(), vIndices);
1103: if (vIndices[d] == numLocalVertices[d]-1) {
1104: vIndices[d] = 0;
1105: int neighborV = ALE::CartesianSieveDef::getValue(dim, numNeighborVertices, vIndices, zeroCode) + numNeighborCells;
1107: sendOverlap->addCone(v, neighborRank, neighborV);
1108: recvOverlap->addCone(neighborRank, v, neighborV);
1109: }
1110: }
1111: }
1112: }
1113: if (debug) {
1114: sendOverlap->view("Send Overlap");
1115: recvOverlap->view("Receive Overlap");
1116: }
1117: delete [] numLocalVertices;
1118: delete [] numNeighborVertices;
1119: delete [] indices;
1120: delete [] code;
1121: delete [] zeroCode;
1122: delete [] vIndices;
1123: // Create coordinates
1124: return mesh;
1125: };
1126: };
1127: }
1129: #endif