Actual source code: Field.hh
1: #ifndef included_ALE_Field_hh
2: #define included_ALE_Field_hh
4: #ifndef included_ALE_SieveAlgorithms_hh
5: #include <sieve/SieveAlgorithms.hh>
6: #endif
10: // Sieve need point_type
11: // Section need point_type and value_type
12: // size(), restrict(), update() need orientation which is a Bundle (circular)
13: // Bundle is Sieve+Section
15: // An AbstractSection is a mapping from Sieve points to sets of values
16: // This is our presentation of a section of a fibre bundle,
17: // in which the Topology is the base space, and
18: // the value sets are vectors in the fiber spaces
19: // The main interface to values is through restrict() and update()
20: // This retrieval uses Sieve closure()
21: // We should call these rawRestrict/rawUpdate
22: // The Section must also be able to set/report the dimension of each fiber
23: // for which we use another section we call an \emph{Atlas}
24: // For some storage schemes, we also want offsets to go with these dimensions
25: // We must have some way of limiting the points associated with values
26: // so each section must support a getPatch() call returning the points with associated fibers
27: // I was using the Chart for this
28: // The Section must be able to participate in \emph{completion}
29: // This means restrict to a provided overlap, and exchange in the restricted sections
30: // Completion does not use hierarchy, so we see the Topology as a DiscreteTopology
31: namespace ALE {
32: template<typename Point_, typename Alloc_ = malloc_allocator<Point_> >
33: class DiscreteSieve {
34: public:
35: typedef Point_ point_type;
36: typedef Alloc_ alloc_type;
37: typedef std::vector<point_type, alloc_type> coneSequence;
38: typedef std::vector<point_type, alloc_type> coneSet;
39: typedef std::vector<point_type, alloc_type> coneArray;
40: typedef std::vector<point_type, alloc_type> supportSequence;
41: typedef std::vector<point_type, alloc_type> supportSet;
42: typedef std::vector<point_type, alloc_type> supportArray;
43: typedef std::set<point_type, std::less<point_type>, alloc_type> points_type;
44: typedef points_type baseSequence;
45: typedef points_type capSequence;
46: typedef typename alloc_type::template rebind<points_type>::other points_alloc_type;
47: typedef typename points_alloc_type::pointer points_ptr;
48: typedef typename alloc_type::template rebind<coneSequence>::other coneSequence_alloc_type;
49: typedef typename coneSequence_alloc_type::pointer coneSequence_ptr;
50: protected:
51: Obj<points_type> _points;
52: Obj<coneSequence> _empty;
53: Obj<coneSequence> _return;
54: alloc_type _allocator;
55: void _init() {
56: points_ptr pPoints = points_alloc_type(this->_allocator).allocate(1);
57: points_alloc_type(this->_allocator).construct(pPoints, points_type());
58: this->_points = Obj<points_type>(pPoints, sizeof(points_type));
59: ///this->_points = new points_type();
60: coneSequence_ptr pEmpty = coneSequence_alloc_type(this->_allocator).allocate(1);
61: coneSequence_alloc_type(this->_allocator).construct(pEmpty, coneSequence());
62: this->_empty = Obj<coneSequence>(pEmpty, sizeof(coneSequence));
63: ///this->_empty = new coneSequence();
64: coneSequence_ptr pReturn = coneSequence_alloc_type(this->_allocator).allocate(1);
65: coneSequence_alloc_type(this->_allocator).construct(pReturn, coneSequence());
66: this->_return = Obj<coneSequence>(pReturn, sizeof(coneSequence));
67: ///this->_return = new coneSequence();
68: };
69: public:
70: DiscreteSieve() {
71: this->_init();
72: };
73: template<typename Input>
74: DiscreteSieve(const Obj<Input>& points) {
75: this->_init();
76: this->_points->insert(points->begin(), points->end());
77: }
78: virtual ~DiscreteSieve() {};
79: public:
80: void addPoint(const point_type& point) {
81: this->_points->insert(point);
82: }
83: template<typename Input>
84: void addPoints(const Obj<Input>& points) {
85: this->_points->insert(points->begin(), points->end());
86: }
87: const Obj<coneSequence>& cone(const point_type& p) {return this->_empty;}
88: template<typename Input>
89: const Obj<coneSequence>& cone(const Input& p) {return this->_empty;}
90: const Obj<coneSet>& nCone(const point_type& p, const int level) {
91: if (level == 0) {
92: return this->closure(p);
93: } else {
94: return this->_empty;
95: }
96: }
97: const Obj<coneArray>& closure(const point_type& p) {
98: this->_return->clear();
99: this->_return->push_back(p);
100: return this->_return;
101: }
102: const Obj<supportSequence>& support(const point_type& p) {return this->_empty;}
103: template<typename Input>
104: const Obj<supportSequence>& support(const Input& p) {return this->_empty;}
105: const Obj<supportSet>& nSupport(const point_type& p, const int level) {
106: if (level == 0) {
107: return this->star(p);
108: } else {
109: return this->_empty;
110: }
111: }
112: const Obj<supportArray>& star(const point_type& p) {
113: this->_return->clear();
114: this->_return->push_back(p);
115: return this->_return;
116: }
117: const Obj<capSequence>& roots() {return this->_points;}
118: const Obj<capSequence>& cap() {return this->_points;}
119: const Obj<baseSequence>& leaves() {return this->_points;}
120: const Obj<baseSequence>& base() {return this->_points;}
121: template<typename Color>
122: void addArrow(const point_type& p, const point_type& q, const Color& color) {
123: throw ALE::Exception("Cannot add an arrow to a DiscreteSieve");
124: }
125: void stratify() {};
126: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
127: ostringstream txt;
128: int rank;
130: if (comm == MPI_COMM_NULL) {
131: comm = MPI_COMM_SELF;
132: rank = 0;
133: } else {
134: MPI_Comm_rank(comm, &rank);
135: }
136: if (name == "") {
137: if(rank == 0) {
138: txt << "viewing a DiscreteSieve" << std::endl;
139: }
140: } else {
141: if(rank == 0) {
142: txt << "viewing DiscreteSieve '" << name << "'" << std::endl;
143: }
144: }
145: for(typename points_type::const_iterator p_iter = this->_points->begin(); p_iter != this->_points->end(); ++p_iter) {
146: txt << " Point " << *p_iter << std::endl;
147: }
148: PetscSynchronizedPrintf(comm, txt.str().c_str());
149: PetscSynchronizedFlush(comm);
150: };
151: };
152: // A ConstantSection is the simplest Section
153: // All fibers are dimension 1
154: // All values are equal to a constant
155: // We need no value storage and no communication for completion
156: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
157: class ConstantSection : public ALE::ParallelObject {
158: public:
159: typedef Point_ point_type;
160: typedef Value_ value_type;
161: typedef Alloc_ alloc_type;
162: typedef std::set<point_type, std::less<point_type>, alloc_type> chart_type;
163: protected:
164: chart_type _chart;
165: value_type _value[2];
166: public:
167: ConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
168: _value[1] = 0;
169: };
170: ConstantSection(MPI_Comm comm, const value_type& value, const int debug) : ParallelObject(comm, debug) {
171: _value[0] = value;
172: _value[1] = value;
173: };
174: ConstantSection(MPI_Comm comm, const value_type& value, const value_type& defaultValue, const int debug) : ParallelObject(comm, debug) {
175: _value[0] = value;
176: _value[1] = defaultValue;
177: };
178: public: // Verifiers
179: void checkPoint(const point_type& point) const {
180: if (this->_chart.find(point) == this->_chart.end()) {
181: ostringstream msg;
182: msg << "Invalid section point " << point << std::endl;
183: throw ALE::Exception(msg.str().c_str());
184: }
185: };
186: void checkDimension(const int& dim) {
187: if (dim != 1) {
188: ostringstream msg;
189: msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
190: throw ALE::Exception(msg.str().c_str());
191: }
192: };
193: bool hasPoint(const point_type& point) const {
194: return this->_chart.count(point) > 0;
195: };
196: public: // Accessors
197: const chart_type& getChart() {return this->_chart;};
198: void addPoint(const point_type& point) {
199: this->_chart.insert(point);
200: };
201: template<typename Points>
202: void addPoint(const Obj<Points>& points) {
203: this->_chart.insert(points->begin(), points->end());
204: }
205: template<typename Points>
206: void addPoint(const Points& points) {
207: this->_chart.insert(points.begin(), points.end());
208: }
209: // void addPoint(const std::set<point_type>& points) {
210: // this->_chart.insert(points.begin(), points.end());
211: // };
212: value_type getDefaultValue() {return this->_value[1];};
213: void setDefaultValue(const value_type value) {this->_value[1] = value;};
214: void copy(const Obj<ConstantSection>& section) {
215: const chart_type& chart = section->getChart();
217: this->addPoint(chart);
218: this->_value[0] = section->restrictSpace()[0];
219: this->setDefaultValue(section->getDefaultValue());
220: };
221: public: // Sizes
222: void clear() {
223: this->_chart.clear();
224: };
225: int getFiberDimension(const point_type& p) const {
226: if (this->hasPoint(p)) return 1;
227: return 0;
228: };
229: void setFiberDimension(const point_type& p, int dim) {
230: this->checkDimension(dim);
231: this->addPoint(p);
232: };
233: template<typename Sequence>
234: void setFiberDimension(const Obj<Sequence>& points, int dim) {
235: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
236: this->setFiberDimension(*p_iter, dim);
237: }
238: }
239: void addFiberDimension(const point_type& p, int dim) {
240: if (this->_chart.find(p) != this->_chart.end()) {
241: ostringstream msg;
242: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
243: throw ALE::Exception(msg.str().c_str());
244: } else {
245: this->setFiberDimension(p, dim);
246: }
247: }
248: int size(const point_type& p) {return this->getFiberDimension(p);};
249: void allocatePoint() {};
250: public: // Restriction
251: const value_type *restrictSpace() const {
252: return this->_value;
253: };
254: const value_type *restrictPoint(const point_type& p) const {
255: if (this->hasPoint(p)) {
256: return this->_value;
257: }
258: return &this->_value[1];
259: };
260: void updatePoint(const point_type& p, const value_type v[]) {
261: this->_value[0] = v[0];
262: };
263: void updateAddPoint(const point_type& p, const value_type v[]) {
264: this->_value[0] += v[0];
265: };
266: public:
267: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
268: ostringstream txt;
269: int rank;
271: if (comm == MPI_COMM_NULL) {
272: comm = this->comm();
273: rank = this->commRank();
274: } else {
275: MPI_Comm_rank(comm, &rank);
276: }
277: if (name == "") {
278: if(rank == 0) {
279: txt << "viewing a ConstantSection" << std::endl;
280: }
281: } else {
282: if(rank == 0) {
283: txt << "viewing ConstantSection '" << name << "'" << std::endl;
284: }
285: }
286: txt <<"["<<this->commRank()<<"]: Value " << this->_value[0] << " Default Value " << this->_value[1] << std::endl;
287: PetscSynchronizedPrintf(comm, txt.str().c_str());
288: PetscSynchronizedFlush(comm);
289: };
290: };
291: // A UniformSection often acts as an Atlas
292: // All fibers are the same dimension
293: // Note we can use a ConstantSection for this Atlas
294: // Each point may have a different vector
295: // Thus we need storage for values, and hence must implement completion
296: template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
297: class UniformSection : public ALE::ParallelObject {
298: public:
299: typedef Point_ point_type;
300: typedef Value_ value_type;
301: typedef Alloc_ alloc_type;
302: typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
303: typedef ConstantSection<point_type, int, point_alloc_type> atlas_type;
304: typedef typename atlas_type::chart_type chart_type;
305: typedef struct {value_type v[fiberDim];} fiber_type;
306: typedef typename alloc_type::template rebind<std::pair<const point_type, fiber_type> >::other pair_alloc_type;
307: typedef std::map<point_type, fiber_type, std::less<point_type>, pair_alloc_type> values_type;
308: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
309: typedef typename atlas_alloc_type::pointer atlas_ptr;
310: protected:
311: size_t _contiguous_array_size;
312: value_type *_contiguous_array;
313: Obj<atlas_type> _atlas;
314: values_type _array;
315: fiber_type _emptyValue;
316: alloc_type _allocator;
317: public:
318: UniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _contiguous_array_size(0), _contiguous_array(NULL) {
319: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
320: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
321: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
322: int dim = fiberDim;
323: this->_atlas->updatePoint(*this->_atlas->getChart().begin(), &dim);
324: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
325: };
326: UniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _contiguous_array_size(0), _contiguous_array(NULL), _atlas(atlas) {
327: int dim = fiberDim;
328: this->_atlas->update(*this->_atlas->getChart().begin(), &dim);
329: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
330: };
331: virtual ~UniformSection() {
332: this->_atlas = NULL;
333: if (this->_contiguous_array) {
334: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
335: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
336: }
337: };
338: public:
339: value_type *getRawArray(const int size) {
340: static value_type *array = NULL;
341: static int maxSize = 0;
343: if (size > maxSize) {
344: const value_type dummy(0);
346: if (array) {
347: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
348: this->_allocator.deallocate(array, maxSize);
349: ///delete [] array;
350: }
351: maxSize = size;
352: array = this->_allocator.allocate(maxSize);
353: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
354: ///array = new value_type[maxSize];
355: };
356: return array;
357: };
358: public: // Verifiers
359: bool hasPoint(const point_type& point) {
360: return this->_atlas->hasPoint(point);
361: };
362: void checkDimension(const int& dim) {
363: if (dim != fiberDim) {
364: ostringstream msg;
365: msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
366: throw ALE::Exception(msg.str().c_str());
367: }
368: };
369: public: // Accessors
370: const chart_type& getChart() {return this->_atlas->getChart();};
371: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
372: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
373: void addPoint(const point_type& point) {
374: this->setFiberDimension(point, fiberDim);
375: };
376: template<typename Points>
377: void addPoint(const Obj<Points>& points) {
378: for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
379: this->setFiberDimension(*p_iter, fiberDim);
380: }
381: }
382: void copy(const Obj<UniformSection>& section) {
383: this->getAtlas()->copy(section->getAtlas());
384: const chart_type& chart = section->getChart();
386: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
387: this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
388: }
389: }
390: public: // Sizes
391: void clear() {
392: this->_array.clear();
393: this->_atlas->clear();
394: }
395: int getFiberDimension(const point_type& p) const {
396: return this->_atlas->restrictPoint(p)[0];
397: }
398: void setFiberDimension(const point_type& p, int dim) {
399: this->update();
400: this->checkDimension(dim);
401: this->_atlas->addPoint(p);
402: this->_atlas->updatePoint(p, &dim);
403: }
404: template<typename Sequence>
405: void setFiberDimension(const Obj<Sequence>& points, int dim) {
406: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
407: this->setFiberDimension(*p_iter, dim);
408: }
409: }
410: void setFiberDimension(const std::set<point_type>& points, int dim) {
411: for(typename std::set<point_type>::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
412: this->setFiberDimension(*p_iter, dim);
413: }
414: };
415: void addFiberDimension(const point_type& p, int dim) {
416: if (this->hasPoint(p)) {
417: ostringstream msg;
418: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
419: throw ALE::Exception(msg.str().c_str());
420: } else {
421: this->setFiberDimension(p, dim);
422: }
423: };
424: int size() const {
425: const chart_type& points = this->_atlas->getChart();
426: int size = 0;
428: for(typename chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
429: size += this->getFiberDimension(*p_iter);
430: }
431: return size;
432: };
433: int sizeWithBC() const {
434: return this->size();
435: };
436: void allocatePoint() {};
437: public: // Restriction
438: const value_type *restrictSpace() {
439: const chart_type& chart = this->getChart();
440: const value_type dummy = 0;
441: int k = 0;
443: if (chart.size() > this->_contiguous_array_size*fiberDim) {
444: if (this->_contiguous_array) {
445: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
446: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
447: }
448: this->_contiguous_array_size = chart.size()*fiberDim;
449: this->_contiguous_array = this->_allocator.allocate(this->_contiguous_array_size*fiberDim);
450: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.construct(this->_contiguous_array+i, dummy);}
451: }
452: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
453: const value_type *a = this->_array[*p_iter].v;
455: for(int i = 0; i < fiberDim; ++i, ++k) {
456: this->_contiguous_array[k] = a[i];
457: }
458: }
459: return this->_contiguous_array;
460: };
461: void update() {
462: if (this->_contiguous_array) {
463: const chart_type& chart = this->getChart();
464: int k = 0;
466: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
467: value_type *a = this->_array[*p_iter].v;
469: for(int i = 0; i < fiberDim; ++i, ++k) {
470: a[i] = this->_contiguous_array[k];
471: }
472: }
473: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
474: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
475: this->_contiguous_array_size = 0;
476: this->_contiguous_array = NULL;
477: }
478: };
479: // Return only the values associated to this point, not its closure
480: const value_type *restrictPoint(const point_type& p) {
481: if (this->_array.find(p) == this->_array.end()) return this->_emptyValue.v;
482: this->update();
483: return this->_array[p].v;
484: };
485: // Update only the values associated to this point, not its closure
486: void updatePoint(const point_type& p, const value_type v[]) {
487: this->update();
488: for(int i = 0; i < fiberDim; ++i) {
489: this->_array[p].v[i] = v[i];
490: }
491: };
492: // Update only the values associated to this point, not its closure
493: void updateAddPoint(const point_type& p, const value_type v[]) {
494: this->update();
495: for(int i = 0; i < fiberDim; ++i) {
496: this->_array[p].v[i] += v[i];
497: }
498: };
499: void updatePointAll(const point_type& p, const value_type v[]) {
500: this->updatePoint(p, v);
501: };
502: public:
503: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
504: ostringstream txt;
505: int rank;
507: this->update();
508: if (comm == MPI_COMM_NULL) {
509: comm = this->comm();
510: rank = this->commRank();
511: } else {
512: MPI_Comm_rank(comm, &rank);
513: }
514: if (name == "") {
515: if(rank == 0) {
516: txt << "viewing a UniformSection" << std::endl;
517: }
518: } else {
519: if(rank == 0) {
520: txt << "viewing UniformSection '" << name << "'" << std::endl;
521: }
522: }
523: const typename atlas_type::chart_type& chart = this->_atlas->getChart();
524: values_type& array = this->_array;
526: for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
527: const point_type& point = *p_iter;
528: const typename atlas_type::value_type dim = this->_atlas->restrictPoint(point)[0];
530: if (dim != 0) {
531: txt << "[" << this->commRank() << "]: " << point << " dim " << dim << " ";
532: for(int i = 0; i < dim; i++) {
533: txt << " " << array[point].v[i];
534: }
535: txt << std::endl;
536: }
537: }
538: if (chart.size() == 0) {
539: txt << "[" << this->commRank() << "]: empty" << std::endl;
540: }
541: PetscSynchronizedPrintf(comm, txt.str().c_str());
542: PetscSynchronizedFlush(comm);
543: };
544: };
546: // A FauxConstantSection is the simplest Section
547: // All fibers are dimension 1
548: // All values are equal to a constant
549: // We need no value storage and no communication for completion
550: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
551: class FauxConstantSection : public ALE::ParallelObject {
552: public:
553: typedef Point_ point_type;
554: typedef Value_ value_type;
555: typedef Alloc_ alloc_type;
556: protected:
557: value_type _value;
558: public:
559: FauxConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
560: FauxConstantSection(MPI_Comm comm, const value_type& value, const int debug) : ParallelObject(comm, debug), _value(value) {};
561: public: // Verifiers
562: void checkDimension(const int& dim) {
563: if (dim != 1) {
564: ostringstream msg;
565: msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
566: throw ALE::Exception(msg.str().c_str());
567: }
568: };
569: public: // Accessors
570: void addPoint(const point_type& point) {
571: }
572: template<typename Points>
573: void addPoint(const Obj<Points>& points) {
574: }
575: template<typename Points>
576: void addPoint(const Points& points) {
577: }
578: void copy(const Obj<FauxConstantSection>& section) {
579: this->_value = section->restrictPoint(point_type())[0];
580: }
581: public: // Sizes
582: void clear() {
583: };
584: int getFiberDimension(const point_type& p) const {
585: return 1;
586: };
587: void setFiberDimension(const point_type& p, int dim) {
588: this->checkDimension(dim);
589: this->addPoint(p);
590: };
591: template<typename Sequence>
592: void setFiberDimension(const Obj<Sequence>& points, int dim) {
593: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
594: this->setFiberDimension(*p_iter, dim);
595: }
596: }
597: void addFiberDimension(const point_type& p, int dim) {
598: this->setFiberDimension(p, dim);
599: }
600: int size(const point_type& p) {return this->getFiberDimension(p);}
601: public: // Restriction
602: const value_type *restrictPoint(const point_type& p) const {
603: return &this->_value;
604: };
605: void updatePoint(const point_type& p, const value_type v[]) {
606: this->_value = v[0];
607: };
608: void updateAddPoint(const point_type& p, const value_type v[]) {
609: this->_value += v[0];
610: };
611: public:
612: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
613: ostringstream txt;
614: int rank;
616: if (comm == MPI_COMM_NULL) {
617: comm = this->comm();
618: rank = this->commRank();
619: } else {
620: MPI_Comm_rank(comm, &rank);
621: }
622: if (name == "") {
623: if(rank == 0) {
624: txt << "viewing a FauxConstantSection" << std::endl;
625: }
626: } else {
627: if(rank == 0) {
628: txt << "viewing FauxConstantSection '" << name << "'" << std::endl;
629: }
630: }
631: txt <<"["<<this->commRank()<<"]: Value " << this->_value << std::endl;
632: PetscSynchronizedPrintf(comm, txt.str().c_str());
633: PetscSynchronizedFlush(comm);
634: };
635: };
636: // Make a simple set from the keys of a map
637: template<typename Map>
638: class SetFromMap {
639: public:
640: typedef typename Map::size_type size_type;
641: class const_iterator {
642: public:
643: typedef typename Map::key_type value_type;
644: typedef typename Map::size_type size_type;
645: protected:
646: typename Map::const_iterator _iter;
647: public:
648: const_iterator(const typename Map::const_iterator& iter): _iter(iter) {};
649: ~const_iterator() {};
650: public:
651: const_iterator& operator=(const const_iterator& iter) {this->_iter = iter._iter; return this->_iter;};
652: bool operator==(const const_iterator& iter) const {return this->_iter == iter._iter;};
653: bool operator!=(const const_iterator& iter) const {return this->_iter != iter._iter;};
654: const_iterator& operator++() {++this->_iter; return *this;}
655: const_iterator operator++(int) {
656: const_iterator tmp(*this);
657: ++(*this);
658: return tmp;
659: };
660: const_iterator& operator--() {--this->_iter; return *this;}
661: const_iterator operator--(int) {
662: const_iterator tmp(*this);
663: --(*this);
664: return tmp;
665: };
666: value_type operator*() const {return this->_iter->first;};
667: };
668: protected:
669: const Map& _map;
670: public:
671: SetFromMap(const Map& map): _map(map) {};
672: public:
673: const_iterator begin() const {return const_iterator(this->_map.begin());};
674: const_iterator end() const {return const_iterator(this->_map.end());};
675: size_type size() const {return this->_map.size();};
676: };
677: // A NewUniformSection often acts as an Atlas
678: // All fibers are the same dimension
679: // Note we can use a ConstantSection for this Atlas
680: // Each point may have a different vector
681: // Thus we need storage for values, and hence must implement completion
682: template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
683: class NewUniformSection : public ALE::ParallelObject {
684: public:
685: typedef Point_ point_type;
686: typedef Value_ value_type;
687: typedef Alloc_ alloc_type;
688: typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
689: typedef FauxConstantSection<point_type, int, point_alloc_type> atlas_type;
690: typedef struct {value_type v[fiberDim];} fiber_type;
691: typedef typename alloc_type::template rebind<std::pair<const point_type, fiber_type> >::other pair_alloc_type;
692: typedef std::map<point_type, fiber_type, std::less<point_type>, pair_alloc_type> values_type;
693: typedef SetFromMap<values_type> chart_type;
694: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
695: typedef typename atlas_alloc_type::pointer atlas_ptr;
696: protected:
697: values_type _array;
698: chart_type _chart;
699: size_t _contiguous_array_size;
700: value_type *_contiguous_array;
701: Obj<atlas_type> _atlas;
702: fiber_type _emptyValue;
703: alloc_type _allocator;
704: public:
705: NewUniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _chart(_array), _contiguous_array_size(0), _contiguous_array(NULL) {
706: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
707: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
708: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
709: int dim = fiberDim;
710: this->_atlas->update(point_type(), &dim);
711: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
712: };
713: NewUniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _chart(_array), _contiguous_array_size(0), _contiguous_array(NULL), _atlas(atlas) {
714: int dim = fiberDim;
715: this->_atlas->update(point_type(), &dim);
716: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
717: };
718: ~NewUniformSection() {
719: this->_atlas = NULL;
720: if (this->_contiguous_array) {
721: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
722: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
723: }
724: };
725: public:
726: value_type *getRawArray(const int size) {
727: static value_type *array = NULL;
728: static int maxSize = 0;
730: if (size > maxSize) {
731: const value_type dummy(0);
733: if (array) {
734: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
735: this->_allocator.deallocate(array, maxSize);
736: ///delete [] array;
737: }
738: maxSize = size;
739: array = this->_allocator.allocate(maxSize);
740: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
741: ///array = new value_type[maxSize];
742: };
743: return array;
744: };
745: public: // Verifiers
746: bool hasPoint(const point_type& point) {
747: return (this->_array.find(point) != this->_array.end());
748: ///return this->_atlas->hasPoint(point);
749: };
750: void checkDimension(const int& dim) {
751: if (dim != fiberDim) {
752: ostringstream msg;
753: msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
754: throw ALE::Exception(msg.str().c_str());
755: }
756: };
757: public: // Accessors
758: const chart_type& getChart() {return this->_chart;}
759: const Obj<atlas_type>& getAtlas() {return this->_atlas;}
760: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;}
761: void addPoint(const point_type& point) {
762: this->setFiberDimension(point, fiberDim);
763: }
764: template<typename Points>
765: void addPoint(const Obj<Points>& points) {
766: for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
767: this->setFiberDimension(*p_iter, fiberDim);
768: }
769: }
770: void copy(const Obj<NewUniformSection>& section) {
771: this->getAtlas()->copy(section->getAtlas());
772: const chart_type& chart = section->getChart();
774: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
775: this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
776: }
777: }
778: public: // Sizes
779: void clear() {
780: this->_array.clear();
781: this->_atlas->clear();
782: };
783: int getFiberDimension(const point_type& p) const {
784: return fiberDim;
785: };
786: void setFiberDimension(const point_type& p, int dim) {
787: this->update();
788: this->checkDimension(dim);
789: this->_atlas->addPoint(p);
790: this->_atlas->updatePoint(p, &dim);
791: this->_array[p] = fiber_type();
792: };
793: template<typename Sequence>
794: void setFiberDimension(const Obj<Sequence>& points, int dim) {
795: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
796: this->setFiberDimension(*p_iter, dim);
797: }
798: }
799: void setFiberDimension(const std::set<point_type>& points, int dim) {
800: for(typename std::set<point_type>::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
801: this->setFiberDimension(*p_iter, dim);
802: }
803: };
804: void addFiberDimension(const point_type& p, int dim) {
805: if (this->hasPoint(p)) {
806: ostringstream msg;
807: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
808: throw ALE::Exception(msg.str().c_str());
809: } else {
810: this->setFiberDimension(p, dim);
811: }
812: };
813: int size() {
814: const chart_type& points = this->getChart();
815: int size = 0;
817: for(typename chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
818: size += this->getFiberDimension(*p_iter);
819: }
820: return size;
821: };
822: int sizeWithBC() {
823: return this->size();
824: };
825: void allocatePoint() {};
826: public: // Restriction
827: // Return a pointer to the entire contiguous storage array
828: const value_type *restrictSpace() {
829: const chart_type& chart = this->getChart();
830: const value_type dummy = 0;
831: int k = 0;
833: if (chart.size() > this->_contiguous_array_size*fiberDim) {
834: if (this->_contiguous_array) {
835: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
836: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
837: }
838: this->_contiguous_array_size = chart.size()*fiberDim;
839: this->_contiguous_array = this->_allocator.allocate(this->_contiguous_array_size*fiberDim);
840: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.construct(this->_contiguous_array+i, dummy);}
841: }
842: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
843: const value_type *a = this->_array[*p_iter].v;
845: for(int i = 0; i < fiberDim; ++i, ++k) {
846: this->_contiguous_array[k] = a[i];
847: }
848: }
849: return this->_contiguous_array;
850: };
851: void update() {
852: if (this->_contiguous_array) {
853: const chart_type& chart = this->getChart();
854: int k = 0;
856: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
857: value_type *a = this->_array[*p_iter].v;
859: for(int i = 0; i < fiberDim; ++i, ++k) {
860: a[i] = this->_contiguous_array[k];
861: }
862: }
863: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
864: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
865: this->_contiguous_array_size = 0;
866: this->_contiguous_array = NULL;
867: }
868: };
869: // Return only the values associated to this point, not its closure
870: const value_type *restrictPoint(const point_type& p) {
871: if (this->_array.find(p) == this->_array.end()) return this->_emptyValue.v;
872: this->update();
873: return this->_array[p].v;
874: };
875: // Update only the values associated to this point, not its closure
876: void updatePoint(const point_type& p, const value_type v[]) {
877: this->update();
878: for(int i = 0; i < fiberDim; ++i) {
879: this->_array[p].v[i] = v[i];
880: }
881: };
882: // Update only the values associated to this point, not its closure
883: void updateAddPoint(const point_type& p, const value_type v[]) {
884: this->update();
885: for(int i = 0; i < fiberDim; ++i) {
886: this->_array[p].v[i] += v[i];
887: }
888: };
889: void updatePointAll(const point_type& p, const value_type v[]) {
890: this->updatePoint(p, v);
891: };
892: public:
893: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
894: ostringstream txt;
895: int rank;
897: this->update();
898: if (comm == MPI_COMM_NULL) {
899: comm = this->comm();
900: rank = this->commRank();
901: } else {
902: MPI_Comm_rank(comm, &rank);
903: }
904: if (name == "") {
905: if(rank == 0) {
906: txt << "viewing a NewUniformSection" << std::endl;
907: }
908: } else {
909: if(rank == 0) {
910: txt << "viewing NewUniformSection '" << name << "'" << std::endl;
911: }
912: }
913: const chart_type& chart = this->getChart();
914: values_type& array = this->_array;
916: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
917: const point_type& point = *p_iter;
919: if (fiberDim != 0) {
920: txt << "[" << this->commRank() << "]: " << point << " dim " << fiberDim << " ";
921: for(int i = 0; i < fiberDim; i++) {
922: txt << " " << array[point].v[i];
923: }
924: txt << std::endl;
925: }
926: }
927: if (chart.size() == 0) {
928: txt << "[" << this->commRank() << "]: empty" << std::endl;
929: }
930: PetscSynchronizedPrintf(comm, txt.str().c_str());
931: PetscSynchronizedFlush(comm);
932: };
933: };
934: // A Section is our most general construct (more general ones could be envisioned)
935: // The Atlas is a UniformSection of dimension 1 and value type Point
936: // to hold each fiber dimension and offsets into a contiguous patch array
937: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
938: typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >
939: class Section : public ALE::ParallelObject {
940: public:
941: typedef Point_ point_type;
942: typedef Value_ value_type;
943: typedef Alloc_ alloc_type;
944: typedef Atlas_ atlas_type;
945: typedef Point index_type;
946: typedef typename atlas_type::chart_type chart_type;
947: typedef value_type * values_type;
948: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
949: typedef typename atlas_alloc_type::pointer atlas_ptr;
950: protected:
951: Obj<atlas_type> _atlas;
952: Obj<atlas_type> _atlasNew;
953: values_type _array;
954: alloc_type _allocator;
955: public:
956: Section(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
957: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
958: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
959: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
960: this->_atlasNew = NULL;
961: this->_array = NULL;
962: };
963: Section(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL), _array(NULL) {};
964: virtual ~Section() {
965: if (this->_array) {
966: const int totalSize = this->sizeWithBC();
968: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
969: this->_allocator.deallocate(this->_array, totalSize);
970: ///delete [] this->_array;
971: this->_array = NULL;
972: }
973: };
974: public:
975: value_type *getRawArray(const int size) {
976: static value_type *array = NULL;
977: static int maxSize = 0;
979: if (size > maxSize) {
980: const value_type dummy(0);
982: if (array) {
983: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
984: this->_allocator.deallocate(array, maxSize);
985: ///delete [] array;
986: }
987: maxSize = size;
988: array = this->_allocator.allocate(maxSize);
989: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
990: ///array = new value_type[maxSize];
991: };
992: return array;
993: };
994: int getStorageSize() const {
995: return this->sizeWithBC();
996: };
997: public: // Verifiers
998: bool hasPoint(const point_type& point) {
999: return this->_atlas->hasPoint(point);
1000: };
1001: public: // Accessors
1002: const chart_type& getChart() {return this->_atlas->getChart();};
1003: void setChart(chart_type& chart) {};
1004: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
1005: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1006: const Obj<atlas_type>& getNewAtlas() {return this->_atlasNew;};
1007: void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1008: const chart_type& getChart() const {return this->_atlas->getChart();};
1009: public: // BC
1010: // Returns the number of constraints on a point
1011: int getConstraintDimension(const point_type& p) const {
1012: return std::max(0, -this->_atlas->restrictPoint(p)->prefix);
1013: }
1014: // Set the number of constraints on a point
1015: // We only allow the entire point to be constrained, so these will be the
1016: // only dofs on the point
1017: void setConstraintDimension(const point_type& p, const int numConstraints) {
1018: this->setFiberDimension(p, -numConstraints);
1019: };
1020: void addConstraintDimension(const point_type& p, const int numConstraints) {
1021: throw ALE::Exception("Variable constraint dimensions not available for this Section type");
1022: };
1023: void copyBC(const Obj<Section>& section) {
1024: const chart_type& chart = this->getChart();
1026: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1027: if (this->getConstraintDimension(*p_iter)) {
1028: this->updatePointBC(*p_iter, section->restrictPoint(*p_iter));
1029: }
1030: }
1031: };
1032: void defaultConstraintDof() {};
1033: public: // Sizes
1034: void clear() {
1035: const int totalSize = this->sizeWithBC();
1037: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1038: this->_allocator.deallocate(this->_array, totalSize);
1039: ///delete [] this->_array;
1040: this->_array = NULL;
1041: this->_atlas->clear();
1042: };
1043: // Return the total number of dofs on the point (free and constrained)
1044: int getFiberDimension(const point_type& p) const {
1045: return std::abs(this->_atlas->restrictPoint(p)->prefix);
1046: };
1047: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
1048: return std::abs(atlas->restrictPoint(p)->prefix);
1049: };
1050: // Set the total number of dofs on the point (free and constrained)
1051: void setFiberDimension(const point_type& p, int dim) {
1052: const index_type idx(dim, -1);
1053: this->_atlas->addPoint(p);
1054: this->_atlas->updatePoint(p, &idx);
1055: };
1056: template<typename Sequence>
1057: void setFiberDimension(const Obj<Sequence>& points, int dim) {
1058: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1059: this->setFiberDimension(*p_iter, dim);
1060: }
1061: }
1062: void addFiberDimension(const point_type& p, int dim) {
1063: if (this->_atlas->hasPoint(p)) {
1064: const index_type values(dim, 0);
1065: this->_atlas->updateAddPoint(p, &values);
1066: } else {
1067: this->setFiberDimension(p, dim);
1068: }
1069: }
1070: // Return the number of constrained dofs on this point
1071: // If constrained, this is equal to the fiber dimension
1072: // Otherwise, 0
1073: int getConstrainedFiberDimension(const point_type& p) const {
1074: return std::max((index_type::prefix_type) 0, this->_atlas->restrictPoint(p)->prefix);
1075: };
1076: // Return the total number of free dofs
1077: int size() const {
1078: const chart_type& points = this->getChart();
1079: int size = 0;
1081: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1082: size += this->getConstrainedFiberDimension(*p_iter);
1083: }
1084: return size;
1085: };
1086: // Return the total number of dofs (free and constrained)
1087: int sizeWithBC() const {
1088: const chart_type& points = this->getChart();
1089: int size = 0;
1091: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1092: size += this->getFiberDimension(*p_iter);
1093: }
1094: return size;
1095: };
1096: public: // Index retrieval
1097: const typename index_type::index_type& getIndex(const point_type& p) {
1098: return this->_atlas->restrictPoint(p)->index;
1099: };
1100: void setIndex(const point_type& p, const typename index_type::index_type& index) {
1101: ((typename atlas_type::value_type *) this->_atlas->restrictPoint(p))->index = index;
1102: };
1103: void setIndexBC(const point_type& p, const typename index_type::index_type& index) {
1104: this->setIndex(p, index);
1105: };
1106: void getIndices(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1107: this->getIndices(p, this->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1108: };
1109: template<typename Order_>
1110: void getIndices(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1111: this->getIndices(p, order->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1112: }
1113: void getIndices(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1114: const int& dim = this->getFiberDimension(p);
1115: const int& cDim = this->getConstraintDimension(p);
1116: const int end = start + dim;
1118: if (!cDim) {
1119: if (orientation >= 0) {
1120: for(int i = start; i < end; ++i) {
1121: indices[(*indx)++] = i;
1122: }
1123: } else {
1124: for(int i = end-1; i >= start; --i) {
1125: indices[(*indx)++] = i;
1126: }
1127: }
1128: } else {
1129: if (!freeOnly) {
1130: for(int i = start; i < end; ++i) {
1131: indices[(*indx)++] = -1;
1132: }
1133: }
1134: }
1135: };
1136: public: // Allocation
1137: void allocateStorage() {
1138: const int totalSize = this->sizeWithBC();
1139: const value_type dummy(0);
1141: this->_array = this->_allocator.allocate(totalSize);
1142: ///this->_array = new value_type[totalSize];
1143: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1144: ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1145: };
1146: void replaceStorage(value_type *newArray) {
1147: const int totalSize = this->sizeWithBC();
1149: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1150: this->_allocator.deallocate(this->_array, totalSize);
1151: ///delete [] this->_array;
1152: this->_array = newArray;
1153: this->_atlasNew = NULL;
1154: };
1155: void addPoint(const point_type& point, const int dim) {
1156: if (dim == 0) return;
1157: if (this->_atlasNew.isNull()) {
1158: this->_atlasNew = new atlas_type(this->comm(), this->debug());
1159: this->_atlasNew->copy(this->_atlas);
1160: }
1161: const index_type idx(dim, -1);
1162: this->_atlasNew->addPoint(point);
1163: this->_atlasNew->updatePoint(point, &idx);
1164: };
1165: template<typename Sequence>
1166: void addPoints(const Obj<Sequence>& points, const int dim) {
1167: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1168: this->addPoint(*p_iter, dim);
1169: }
1170: }
1171: void orderPoints(const Obj<atlas_type>& atlas){
1172: const chart_type& chart = this->getChart();
1173: int offset = 0;
1174: int bcOffset = this->size();
1176: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1177: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
1178: const int& dim = idx.prefix;
1180: if (dim >= 0) {
1181: idx.index = offset;
1182: atlas->updatePoint(*c_iter, &idx);
1183: offset += dim;
1184: } else {
1185: idx.index = bcOffset;
1186: atlas->updatePoint(*c_iter, &idx);
1187: bcOffset -= dim;
1188: }
1189: }
1190: };
1191: void allocatePoint() {
1192: this->orderPoints(this->_atlas);
1193: this->allocateStorage();
1194: };
1195: public: // Restriction and Update
1196: // Zero entries
1197: void zero() {
1198: memset(this->_array, 0, this->size()* sizeof(value_type));
1199: };
1200: // Return a pointer to the entire contiguous storage array
1201: const value_type *restrictSpace() {
1202: return this->_array;
1203: };
1204: // Update the entire contiguous storage array
1205: void update(const value_type v[]) {
1206: const value_type *array = this->_array;
1207: const int size = this->size();
1209: for(int i = 0; i < size; i++) {
1210: array[i] = v[i];
1211: }
1212: };
1213: // Return the free values on a point
1214: const value_type *restrictPoint(const point_type& p) {
1215: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1216: };
1217: // Update the free values on a point
1218: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1219: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1220: value_type *a = &(this->_array[idx.index]);
1222: if (orientation >= 0) {
1223: for(int i = 0; i < idx.prefix; ++i) {
1224: a[i] = v[i];
1225: }
1226: } else {
1227: const int last = idx.prefix-1;
1229: for(int i = 0; i < idx.prefix; ++i) {
1230: a[i] = v[last-i];
1231: }
1232: }
1233: };
1234: // Update the free values on a point
1235: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1236: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1237: value_type *a = &(this->_array[idx.index]);
1239: if (orientation >= 0) {
1240: for(int i = 0; i < idx.prefix; ++i) {
1241: a[i] += v[i];
1242: }
1243: } else {
1244: const int last = idx.prefix-1;
1246: for(int i = 0; i < idx.prefix; ++i) {
1247: a[i] += v[last-i];
1248: }
1249: }
1250: };
1251: // Update only the constrained dofs on a point
1252: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
1253: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1254: const int dim = -idx.prefix;
1255: value_type *a = &(this->_array[idx.index]);
1257: if (orientation >= 0) {
1258: for(int i = 0; i < dim; ++i) {
1259: a[i] = v[i];
1260: }
1261: } else {
1262: const int last = dim-1;
1264: for(int i = 0; i < dim; ++i) {
1265: a[i] = v[last-i];
1266: }
1267: }
1268: };
1269: // Update all dofs on a point (free and constrained)
1270: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
1271: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1272: const int dim = std::abs(idx.prefix);
1273: value_type *a = &(this->_array[idx.index]);
1275: if (orientation >= 0) {
1276: for(int i = 0; i < dim; ++i) {
1277: a[i] = v[i];
1278: }
1279: } else {
1280: const int last = dim-1;
1282: for(int i = 0; i < dim; ++i) {
1283: a[i] = v[last-i];
1284: }
1285: }
1286: };
1287: public:
1288: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
1289: ostringstream txt;
1290: int rank;
1292: if (comm == MPI_COMM_NULL) {
1293: comm = this->comm();
1294: rank = this->commRank();
1295: } else {
1296: MPI_Comm_rank(comm, &rank);
1297: }
1298: if(rank == 0) {
1299: txt << "viewing Section " << this->getName() << std::endl;
1300: if (name != "") {
1301: txt << ": " << name << "'";
1302: }
1303: txt << std::endl;
1304: }
1305: const typename atlas_type::chart_type& chart = this->_atlas->getChart();
1306: const value_type *array = this->_array;
1308: for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1309: const point_type& point = *p_iter;
1310: const index_type& idx = this->_atlas->restrictPoint(point)[0];
1311: const int dim = this->getFiberDimension(point);
1313: if (idx.prefix != 0) {
1314: txt << "[" << this->commRank() << "]: " << point << " dim " << idx.prefix << " offset " << idx.index << " ";
1315: for(int i = 0; i < dim; i++) {
1316: txt << " " << array[idx.index+i];
1317: }
1318: txt << std::endl;
1319: }
1320: }
1321: if (chart.size() == 0) {
1322: txt << "[" << this->commRank() << "]: empty" << std::endl;
1323: }
1324: PetscSynchronizedPrintf(comm, txt.str().c_str());
1325: PetscSynchronizedFlush(comm);
1326: };
1327: public: // Fibrations
1328: void setConstraintDof(const point_type& p, const int dofs[]) {};
1329: int getNumSpaces() const {return 1;};
1330: void addSpace() {};
1331: int getFiberDimension(const point_type& p, const int space) {return this->getFiberDimension(p);};
1332: void setFiberDimension(const point_type& p, int dim, const int space) {this->setFiberDimension(p, dim);};
1333: template<typename Sequence>
1334: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
1335: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1336: this->setFiberDimension(*p_iter, dim, space);
1337: }
1338: }
1339: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
1340: this->setConstraintDimension(p, numConstraints);
1341: }
1342: };
1343: // GeneralSection will support BC on a subset of unknowns on a point
1344: // We make a separate constraint Atlas to mark constrained dofs on a point
1345: // Storage will be contiguous by node, just as in Section
1346: // This allows fast restrict(p)
1347: // Then update() is accomplished by skipping constrained unknowns
1348: // We must eliminate restrictSpace() since it does not correspond to the constrained system
1349: // Numbering will have to be rewritten to calculate correct mappings
1350: // I think we can just generate multiple tuples per point
1351: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
1352: typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
1353: typename BCAtlas_ = Section<Point_, int, typename Alloc_::template rebind<int>::other> >
1354: class GeneralSection : public ALE::ParallelObject {
1355: public:
1356: typedef Point_ point_type;
1357: typedef Value_ value_type;
1358: typedef Alloc_ alloc_type;
1359: typedef Atlas_ atlas_type;
1360: typedef BCAtlas_ bc_type;
1361: typedef Point index_type;
1362: typedef typename atlas_type::chart_type chart_type;
1363: typedef value_type * values_type;
1364: typedef std::pair<const int *, const int *> customAtlasInd_type;
1365: typedef std::pair<customAtlasInd_type, bool> customAtlas_type;
1366: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
1367: typedef typename atlas_alloc_type::pointer atlas_ptr;
1368: typedef typename alloc_type::template rebind<bc_type>::other bc_alloc_type;
1369: typedef typename bc_alloc_type::pointer bc_ptr;
1370: protected:
1371: // Describes layout of storage, point --> (# of values, offset into array)
1372: Obj<atlas_type> _atlas;
1373: // Spare atlas to allow dynamic updates
1374: Obj<atlas_type> _atlasNew;
1375: // Storage
1376: values_type _array;
1377: bool _sharedStorage;
1378: int _sharedStorageSize;
1379: // A section that maps points to sets of constrained local dofs
1380: Obj<bc_type> _bc;
1381: alloc_type _allocator;
1382: // Fibration structures
1383: // _spaces is a set of atlases which describe the layout of each in the storage of this section
1384: // _bcs is the same as _bc, but for each field
1385: std::vector<Obj<atlas_type> > _spaces;
1386: std::vector<Obj<bc_type> > _bcs;
1387: // Optimization
1388: std::vector<customAtlas_type> _customRestrictAtlas;
1389: std::vector<customAtlas_type> _customUpdateAtlas;
1390: public:
1391: GeneralSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
1392: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
1393: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
1394: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
1395: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
1396: bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
1397: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
1398: this->_atlasNew = NULL;
1399: this->_array = NULL;
1400: this->_sharedStorage = false;
1401: };
1402: GeneralSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
1403: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
1404: bc_alloc_type(this->_allocator).construct(pBC, bc_type(this->comm(), this->debug()));
1405: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
1406: };
1407: ~GeneralSection() {
1408: if (this->_array && !this->_sharedStorage) {
1409: const int totalSize = this->sizeWithBC();
1411: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1412: this->_allocator.deallocate(this->_array, totalSize);
1413: ///delete [] this->_array;
1414: this->_array = NULL;
1415: }
1416: for(std::vector<customAtlas_type>::iterator a_iter = this->_customRestrictAtlas.begin(); a_iter != this->_customRestrictAtlas.end(); ++a_iter) {
1417: if (a_iter->second) {
1418: delete [] a_iter->first.first;
1419: delete [] a_iter->first.second;
1420: }
1421: }
1422: for(std::vector<customAtlas_type>::iterator a_iter = this->_customUpdateAtlas.begin(); a_iter != this->_customUpdateAtlas.end(); ++a_iter) {
1423: if (a_iter->second) {
1424: delete [] a_iter->first.first;
1425: delete [] a_iter->first.second;
1426: }
1427: }
1428: };
1429: public:
1430: value_type *getRawArray(const int size) {
1431: // Put in a sentinel value that deallocates the array
1432: static value_type *array = NULL;
1433: static int maxSize = 0;
1435: if (size > maxSize) {
1436: const value_type dummy(0);
1438: if (array) {
1439: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
1440: this->_allocator.deallocate(array, maxSize);
1441: ///delete [] array;
1442: }
1443: maxSize = size;
1444: array = this->_allocator.allocate(maxSize);
1445: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
1446: ///array = new value_type[maxSize];
1447: };
1448: return array;
1449: };
1450: int getStorageSize() const {
1451: if (this->_sharedStorage) {
1452: return this->_sharedStorageSize;
1453: }
1454: return this->sizeWithBC();
1455: };
1456: bool sharedStorage() const {return this->_sharedStorage;};
1457: public: // Verifiers
1458: bool hasPoint(const point_type& point) {
1459: return this->_atlas->hasPoint(point);
1460: };
1461: public: // Accessors
1462: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
1463: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1464: const Obj<atlas_type>& getNewAtlas() const {return this->_atlasNew;};
1465: void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1466: const Obj<bc_type>& getBC() const {return this->_bc;};
1467: void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
1468: const chart_type& getChart() const {return this->_atlas->getChart();};
1469: void setChart(const chart_type& chart) {throw ALE::Exception("setChart() for GeneralSection is invalid");};
1470: public: // BC
1471: // Returns the number of constraints on a point
1472: int getConstraintDimension(const point_type& p) const {
1473: if (!this->_bc->hasPoint(p)) return 0;
1474: return this->_bc->getFiberDimension(p);
1475: }
1476: // Set the number of constraints on a point
1477: void setConstraintDimension(const point_type& p, const int numConstraints) {
1478: this->_bc->setFiberDimension(p, numConstraints);
1479: }
1480: // Increment the number of constraints on a point
1481: void addConstraintDimension(const point_type& p, const int numConstraints) {
1482: this->_bc->addFiberDimension(p, numConstraints);
1483: }
1484: // Return the local dofs which are constrained on a point
1485: const int *getConstraintDof(const point_type& p) const {
1486: return this->_bc->restrictPoint(p);
1487: }
1488: // Set the local dofs which are constrained on a point
1489: void setConstraintDof(const point_type& p, const int dofs[]) {
1490: this->_bc->updatePoint(p, dofs);
1491: }
1492: template<typename OtherSection>
1493: void copyBC(const Obj<OtherSection>& section) {
1494: this->setBC(section->getBC());
1495: const chart_type& chart = this->getChart();
1497: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1498: if (this->getConstraintDimension(*p_iter)) {
1499: this->updatePointBCFull(*p_iter, section->restrictPoint(*p_iter));
1500: }
1501: }
1502: this->copySpaces(section);
1503: }
1504: void defaultConstraintDof() {
1505: const chart_type& chart = this->getChart();
1506: int size = 0;
1508: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1509: size = std::max(size, this->getConstraintDimension(*p_iter));
1510: }
1511: int *dofs = new int[size];
1512: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1513: const int cDim = this->getConstraintDimension(*p_iter);
1515: if (cDim) {
1516: for(int d = 0; d < cDim; ++d) {
1517: dofs[d] = d;
1518: }
1519: this->_bc->updatePoint(*p_iter, dofs);
1520: }
1521: }
1522: delete [] dofs;
1523: };
1524: public: // Sizes
1525: void clear() {
1526: if (!this->_sharedStorage) {
1527: const int totalSize = this->sizeWithBC();
1529: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1530: this->_allocator.deallocate(this->_array, totalSize);
1531: ///delete [] this->_array;
1532: }
1533: this->_array = NULL;
1534: this->_atlas->clear();
1535: this->_bc->clear();
1536: };
1537: // Return the total number of dofs on the point (free and constrained)
1538: int getFiberDimension(const point_type& p) const {
1539: return this->_atlas->restrictPoint(p)->prefix;
1540: };
1541: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
1542: return atlas->restrictPoint(p)->prefix;
1543: };
1544: // Set the total number of dofs on the point (free and constrained)
1545: void setFiberDimension(const point_type& p, int dim) {
1546: const index_type idx(dim, -1);
1547: this->_atlas->addPoint(p);
1548: this->_atlas->updatePoint(p, &idx);
1549: };
1550: template<typename Sequence>
1551: void setFiberDimension(const Obj<Sequence>& points, int dim) {
1552: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1553: this->setFiberDimension(*p_iter, dim);
1554: }
1555: }
1556: void addFiberDimension(const point_type& p, int dim) {
1557: if (this->_atlas->hasPoint(p)) {
1558: const index_type values(dim, 0);
1559: this->_atlas->updateAddPoint(p, &values);
1560: } else {
1561: this->setFiberDimension(p, dim);
1562: }
1563: };
1564: // Return the number of constrained dofs on this point
1565: int getConstrainedFiberDimension(const point_type& p) const {
1566: return this->getFiberDimension(p) - this->getConstraintDimension(p);
1567: };
1568: // Return the total number of free dofs
1569: int size() const {
1570: const chart_type& points = this->getChart();
1571: int size = 0;
1573: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1574: size += this->getConstrainedFiberDimension(*p_iter);
1575: }
1576: return size;
1577: };
1578: // Return the total number of dofs (free and constrained)
1579: int sizeWithBC() const {
1580: const chart_type& points = this->getChart();
1581: int size = 0;
1583: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1584: size += this->getFiberDimension(*p_iter);
1585: }
1586: return size;
1587: };
1588: public: // Index retrieval
1589: const typename index_type::index_type& getIndex(const point_type& p) {
1590: return this->_atlas->restrictPoint(p)->index;
1591: };
1592: void setIndex(const point_type& p, const typename index_type::index_type& index) {
1593: ((typename atlas_type::value_type *) this->_atlas->restrictPoint(p))->index = index;
1594: };
1595: void setIndexBC(const point_type& p, const typename index_type::index_type& index) {};
1596: void getIndicesRaw(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1) {
1597: this->getIndicesRaw(p, this->getIndex(p), indices, indx, orientation);
1598: };
1599: template<typename Order_>
1600: void getIndicesRaw(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1) {
1601: this->getIndicesRaw(p, order->getIndex(p), indices, indx, orientation);
1602: }
1603: void getIndices(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = true) {
1604: this->getIndices(p, this->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1605: };
1606: template<typename Order_>
1607: void getIndices(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1608: this->getIndices(p, order->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1609: }
1610: void getIndicesRaw(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation) {
1611: if (orientation >= 0) {
1612: const int& dim = this->getFiberDimension(p);
1613: const int end = start + dim;
1615: for(int i = start; i < end; ++i) {
1616: indices[(*indx)++] = i;
1617: }
1618: } else {
1619: const int numSpaces = this->getNumSpaces();
1620: int offset = start;
1622: for(int space = 0; space < numSpaces; ++space) {
1623: const int& dim = this->getFiberDimension(p, space);
1625: for(int i = offset+dim-1; i >= offset; --i) {
1626: indices[(*indx)++] = i;
1627: }
1628: offset += dim;
1629: }
1630: if (!numSpaces) {
1631: const int& dim = this->getFiberDimension(p);
1633: for(int i = offset+dim-1; i >= offset; --i) {
1634: indices[(*indx)++] = i;
1635: }
1636: offset += dim;
1637: }
1638: }
1639: }
1640: /*
1641: p - The Sieve point
1642: start - Offset for the set of indices on this point
1643: indices - Storage for the indices
1644: indx - Pointer to an offset into the indices argument (to allow composition of index sets)
1645: orientation - A negative argument reverses the indices
1646: freeOnly - If false, include constrained dofs with negative indices, otherwise leave them out
1647: skipConstraints - If true, when a constrained dof is encountered, increment the index (even though it is not placed in indices[])
1648: */
1649: void getIndices(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1650: const int& cDim = this->getConstraintDimension(p);
1652: if (!cDim) {
1653: this->getIndicesRaw(p, start, indices, indx, orientation);
1654: } else {
1655: if (orientation >= 0) {
1656: const int& dim = this->getFiberDimension(p);
1657: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1658: int cInd = 0;
1660: /* i is the returned index, k is the local dof number */
1661: for(int i = start, k = 0; k < dim; ++k) {
1662: if ((cInd < cDim) && (k == cDof[cInd])) {
1663: /* Constrained dof */
1664: if (!freeOnly) indices[(*indx)++] = -(k+1);
1665: if (skipConstraints) ++i;
1666: ++cInd;
1667: } else {
1668: /* Unconstrained dof */
1669: indices[(*indx)++] = i++;
1670: }
1671: }
1672: } else {
1673: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1674: int offset = 0;
1675: int cOffset = 0;
1676: int j = -1;
1678: for(int space = 0; space < this->getNumSpaces(); ++space) {
1679: const int dim = this->getFiberDimension(p, space);
1680: const int tDim = this->getConstrainedFiberDimension(p, space);
1681: int cInd = (dim - tDim)-1;
1683: j += dim;
1684: for(int i = 0, k = start+tDim+offset; i < dim; ++i, --j) {
1685: if ((cInd >= 0) && (j == cDof[cInd+cOffset])) {
1686: if (!freeOnly) indices[(*indx)++] = -(offset+i+1);
1687: if (skipConstraints) --k;
1688: --cInd;
1689: } else {
1690: indices[(*indx)++] = --k;
1691: }
1692: }
1693: j += dim;
1694: offset += dim;
1695: cOffset += dim - tDim;
1696: }
1697: }
1698: }
1699: };
1700: public: // Allocation
1701: void allocateStorage() {
1702: const int totalSize = this->sizeWithBC();
1703: const value_type dummy(0) ;
1705: this->_array = this->_allocator.allocate(totalSize);
1706: ///this->_array = new value_type[totalSize];
1707: this->_sharedStorage = false;
1708: this->_sharedStorageSize = 0;
1709: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1710: ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1711: this->_bc->allocatePoint();
1712: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = this->_bcs.begin(); b_iter != this->_bcs.end(); ++b_iter) {
1713: (*b_iter)->allocatePoint();;
1714: }
1715: };
1716: void replaceStorage(value_type *newArray, const bool sharedStorage = false, const int sharedStorageSize = 0) {
1717: if (this->_array && !this->_sharedStorage) {
1718: const int totalSize = this->sizeWithBC();
1720: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1721: this->_allocator.deallocate(this->_array, totalSize);
1722: ///delete [] this->_array;
1723: }
1724: this->_array = newArray;
1725: this->_sharedStorage = sharedStorage;
1726: this->_sharedStorageSize = sharedStorageSize;
1727: this->_atlas = this->_atlasNew;
1728: this->_atlasNew = NULL;
1729: };
1730: // DANGEROUS
1731: void setStorage(value_type *newArray) {
1732: this->_array = newArray;
1733: };
1734: void addPoint(const point_type& point, const int dim) {
1735: if (dim == 0) return;
1736: if (this->_atlasNew.isNull()) {
1737: this->_atlasNew = new atlas_type(this->comm(), this->debug());
1738: this->_atlasNew->copy(this->_atlas);
1739: }
1740: const index_type idx(dim, -1);
1741: this->_atlasNew->addPoint(point);
1742: this->_atlasNew->updatePoint(point, &idx);
1743: };
1744: void orderPoints(const Obj<atlas_type>& atlas){
1745: const chart_type& chart = this->getChart();
1746: int offset = 0;
1748: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1749: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
1750: const int& dim = idx.prefix;
1752: idx.index = offset;
1753: atlas->updatePoint(*c_iter, &idx);
1754: offset += dim;
1755: }
1756: };
1757: void allocatePoint() {
1758: this->orderPoints(this->_atlas);
1759: this->allocateStorage();
1760: };
1761: public: // Restriction and Update
1762: // Zero entries
1763: void zero() {
1764: this->set((value_type) 0.0);
1765: };
1766: void zeroWithBC() {
1767: memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1768: };
1769: void set(const value_type value) {
1770: //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1771: const chart_type& chart = this->getChart();
1773: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1774: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1775: const int& dim = this->getFiberDimension(*c_iter);
1776: const int& cDim = this->getConstraintDimension(*c_iter);
1778: if (!cDim) {
1779: for(int i = 0; i < dim; ++i) {
1780: array[i] = value;
1781: }
1782: } else {
1783: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1784: int cInd = 0;
1786: for(int i = 0; i < dim; ++i) {
1787: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1788: array[i] = value;
1789: }
1790: }
1791: }
1792: };
1793: // Add two sections and put the result in a third
1794: void add(const Obj<GeneralSection>& x, const Obj<GeneralSection>& y) {
1795: // Check atlases
1796: const chart_type& chart = this->getChart();
1798: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1799: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1800: const value_type *xArray = x->restrictPoint(*c_iter);
1801: const value_type *yArray = y->restrictPoint(*c_iter);
1802: const int& dim = this->getFiberDimension(*c_iter);
1803: const int& cDim = this->getConstraintDimension(*c_iter);
1805: if (!cDim) {
1806: for(int i = 0; i < dim; ++i) {
1807: array[i] = xArray[i] + yArray[i];
1808: }
1809: } else {
1810: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1811: int cInd = 0;
1813: for(int i = 0; i < dim; ++i) {
1814: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1815: array[i] = xArray[i] + yArray[i];
1816: }
1817: }
1818: }
1819: };
1820: // this = this + alpha * x
1821: template<typename OtherSection>
1822: void axpy(const value_type alpha, const Obj<OtherSection>& x) {
1823: // Check atlases
1824: const chart_type& chart = this->getChart();
1826: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1827: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1828: const value_type *xArray = x->restrictPoint(*c_iter);
1829: const int& dim = this->getFiberDimension(*c_iter);
1830: const int& cDim = this->getConstraintDimension(*c_iter);
1832: if (!cDim) {
1833: for(int i = 0; i < dim; ++i) {
1834: array[i] += alpha*xArray[i];
1835: }
1836: } else {
1837: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1838: int cInd = 0;
1840: for(int i = 0; i < dim; ++i) {
1841: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1842: array[i] += alpha*xArray[i];
1843: }
1844: }
1845: }
1846: }
1847: // Return the free values on a point
1848: const value_type *restrictSpace() const {
1849: return this->_array;
1850: }
1851: // Return the free values on a point
1852: const value_type *restrictPoint(const point_type& p) const {
1853: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1854: }
1855: void restrictPoint(const point_type& p, value_type values[], const int size) const {
1856: assert(this->_atlas->restrictPoint(p)[0].prefix == size);
1857: const value_type *v = &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1859: for(int i = 0; i < size; ++i) {
1860: values[i] = v[i];
1861: }
1862: };
1863: // Update the free values on a point
1864: // Takes a full array and ignores constrained values
1865: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1866: value_type *array = (value_type *) this->restrictPoint(p);
1867: const int& cDim = this->getConstraintDimension(p);
1869: if (!cDim) {
1870: if (orientation >= 0) {
1871: const int& dim = this->getFiberDimension(p);
1873: for(int i = 0; i < dim; ++i) {
1874: array[i] = v[i];
1875: }
1876: } else {
1877: int offset = 0;
1878: int j = -1;
1880: for(int space = 0; space < this->getNumSpaces(); ++space) {
1881: const int& dim = this->getFiberDimension(p, space);
1883: for(int i = dim-1; i >= 0; --i) {
1884: array[++j] = v[i+offset];
1885: }
1886: offset += dim;
1887: }
1888: }
1889: } else {
1890: if (orientation >= 0) {
1891: const int& dim = this->getFiberDimension(p);
1892: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1893: int cInd = 0;
1895: for(int i = 0; i < dim; ++i) {
1896: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1897: array[i] = v[i];
1898: }
1899: } else {
1900: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1901: int offset = 0;
1902: int cOffset = 0;
1903: int j = 0;
1905: for(int space = 0; space < this->getNumSpaces(); ++space) {
1906: const int dim = this->getFiberDimension(p, space);
1907: const int tDim = this->getConstrainedFiberDimension(p, space);
1908: const int sDim = dim - tDim;
1909: int cInd = 0;
1911: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1912: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1913: array[j] = v[k];
1914: }
1915: offset += dim;
1916: cOffset += dim - tDim;
1917: }
1918: }
1919: }
1920: };
1921: // Update the free values on a point
1922: // Takes a full array and ignores constrained values
1923: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1924: value_type *array = (value_type *) this->restrictPoint(p);
1925: const int& cDim = this->getConstraintDimension(p);
1927: if (!cDim) {
1928: if (orientation >= 0) {
1929: const int& dim = this->getFiberDimension(p);
1931: for(int i = 0; i < dim; ++i) {
1932: array[i] += v[i];
1933: }
1934: } else {
1935: int offset = 0;
1936: int j = -1;
1938: for(int space = 0; space < this->getNumSpaces(); ++space) {
1939: const int& dim = this->getFiberDimension(p, space);
1941: for(int i = dim-1; i >= 0; --i) {
1942: array[++j] += v[i+offset];
1943: }
1944: offset += dim;
1945: }
1946: }
1947: } else {
1948: if (orientation >= 0) {
1949: const int& dim = this->getFiberDimension(p);
1950: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1951: int cInd = 0;
1953: for(int i = 0; i < dim; ++i) {
1954: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1955: array[i] += v[i];
1956: }
1957: } else {
1958: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1959: int offset = 0;
1960: int cOffset = 0;
1961: int j = 0;
1963: for(int space = 0; space < this->getNumSpaces(); ++space) {
1964: const int dim = this->getFiberDimension(p, space);
1965: const int tDim = this->getConstrainedFiberDimension(p, space);
1966: const int sDim = dim - tDim;
1967: int cInd = 0;
1969: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1970: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1971: array[j] += v[k];
1972: }
1973: offset += dim;
1974: cOffset += dim - tDim;
1975: }
1976: }
1977: }
1978: };
1979: // Update the free values on a point
1980: // Takes ONLY unconstrained values
1981: void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1982: value_type *array = (value_type *) this->restrictPoint(p);
1983: const int& cDim = this->getConstraintDimension(p);
1985: if (!cDim) {
1986: if (orientation >= 0) {
1987: const int& dim = this->getFiberDimension(p);
1989: for(int i = 0; i < dim; ++i) {
1990: array[i] = v[i];
1991: }
1992: } else {
1993: int offset = 0;
1994: int j = -1;
1996: for(int space = 0; space < this->getNumSpaces(); ++space) {
1997: const int& dim = this->getFiberDimension(p, space);
1999: for(int i = dim-1; i >= 0; --i) {
2000: array[++j] = v[i+offset];
2001: }
2002: offset += dim;
2003: }
2004: }
2005: } else {
2006: if (orientation >= 0) {
2007: const int& dim = this->getFiberDimension(p);
2008: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2009: int cInd = 0;
2011: for(int i = 0, k = -1; i < dim; ++i) {
2012: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2013: array[i] = v[++k];
2014: }
2015: } else {
2016: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2017: int offset = 0;
2018: int cOffset = 0;
2019: int j = 0;
2021: for(int space = 0; space < this->getNumSpaces(); ++space) {
2022: const int dim = this->getFiberDimension(p, space);
2023: const int tDim = this->getConstrainedFiberDimension(p, space);
2024: const int sDim = dim - tDim;
2025: int cInd = 0;
2027: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2028: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2029: array[j] = v[--k];
2030: }
2031: offset += dim;
2032: cOffset += dim - tDim;
2033: }
2034: }
2035: }
2036: };
2037: // Update the free values on a point
2038: // Takes ONLY unconstrained values
2039: void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2040: value_type *array = (value_type *) this->restrictPoint(p);
2041: const int& cDim = this->getConstraintDimension(p);
2043: if (!cDim) {
2044: if (orientation >= 0) {
2045: const int& dim = this->getFiberDimension(p);
2047: for(int i = 0; i < dim; ++i) {
2048: array[i] += v[i];
2049: }
2050: } else {
2051: int offset = 0;
2052: int j = -1;
2054: for(int space = 0; space < this->getNumSpaces(); ++space) {
2055: const int& dim = this->getFiberDimension(p, space);
2057: for(int i = dim-1; i >= 0; --i) {
2058: array[++j] += v[i+offset];
2059: }
2060: offset += dim;
2061: }
2062: }
2063: } else {
2064: if (orientation >= 0) {
2065: const int& dim = this->getFiberDimension(p);
2066: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2067: int cInd = 0;
2069: for(int i = 0, k = -1; i < dim; ++i) {
2070: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2071: array[i] += v[++k];
2072: }
2073: } else {
2074: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2075: int offset = 0;
2076: int cOffset = 0;
2077: int j = 0;
2079: for(int space = 0; space < this->getNumSpaces(); ++space) {
2080: const int dim = this->getFiberDimension(p, space);
2081: const int tDim = this->getConstrainedFiberDimension(p, space);
2082: const int sDim = dim - tDim;
2083: int cInd = 0;
2085: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2086: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2087: array[j] += v[--k];
2088: }
2089: offset += dim;
2090: cOffset += dim - tDim;
2091: }
2092: }
2093: }
2094: };
2095: // Update only the constrained dofs on a point
2096: // This takes an array with ONLY bc values
2097: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2098: value_type *array = (value_type *) this->restrictPoint(p);
2099: const int& cDim = this->getConstraintDimension(p);
2101: if (cDim) {
2102: if (orientation >= 0) {
2103: const int& dim = this->getFiberDimension(p);
2104: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2105: int cInd = 0;
2107: for(int i = 0; i < dim; ++i) {
2108: if (cInd == cDim) break;
2109: if (i == cDof[cInd]) {
2110: array[i] = v[cInd];
2111: ++cInd;
2112: }
2113: }
2114: } else {
2115: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2116: int cOffset = 0;
2117: int j = 0;
2119: for(int space = 0; space < this->getNumSpaces(); ++space) {
2120: const int dim = this->getFiberDimension(p, space);
2121: const int tDim = this->getConstrainedFiberDimension(p, space);
2122: int cInd = 0;
2124: for(int i = 0; i < dim; ++i, ++j) {
2125: if (cInd < 0) break;
2126: if (j == cDof[cInd+cOffset]) {
2127: array[j] = v[cInd+cOffset];
2128: ++cInd;
2129: }
2130: }
2131: cOffset += dim - tDim;
2132: }
2133: }
2134: }
2135: };
2136: // Update only the constrained dofs on a point
2137: // This takes an array with ALL values, not just BC
2138: void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2139: value_type *array = (value_type *) this->restrictPoint(p);
2140: const int& cDim = this->getConstraintDimension(p);
2142: if (cDim) {
2143: if (orientation >= 0) {
2144: const int& dim = this->getFiberDimension(p);
2145: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2146: int cInd = 0;
2148: for(int i = 0; i < dim; ++i) {
2149: if (cInd == cDim) break;
2150: if (i == cDof[cInd]) {
2151: array[i] = v[i];
2152: ++cInd;
2153: }
2154: }
2155: } else {
2156: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2157: int offset = 0;
2158: int cOffset = 0;
2159: int j = 0;
2161: for(int space = 0; space < this->getNumSpaces(); ++space) {
2162: const int dim = this->getFiberDimension(p, space);
2163: const int tDim = this->getConstrainedFiberDimension(p, space);
2164: int cInd = 0;
2166: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2167: if (cInd < 0) break;
2168: if (j == cDof[cInd+cOffset]) {
2169: array[j] = v[k];
2170: ++cInd;
2171: }
2172: }
2173: offset += dim;
2174: cOffset += dim - tDim;
2175: }
2176: }
2177: }
2178: };
2179: // Update all dofs on a point (free and constrained)
2180: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
2181: value_type *array = (value_type *) this->restrictPoint(p);
2183: if (orientation >= 0) {
2184: const int& dim = this->getFiberDimension(p);
2186: for(int i = 0; i < dim; ++i) {
2187: array[i] = v[i];
2188: }
2189: } else {
2190: int offset = 0;
2191: int j = -1;
2193: for(int space = 0; space < this->getNumSpaces(); ++space) {
2194: const int& dim = this->getFiberDimension(p, space);
2196: for(int i = dim-1; i >= 0; --i) {
2197: array[++j] = v[i+offset];
2198: }
2199: offset += dim;
2200: }
2201: }
2202: };
2203: // Update all dofs on a point (free and constrained)
2204: void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
2205: value_type *array = (value_type *) this->restrictPoint(p);
2207: if (orientation >= 0) {
2208: const int& dim = this->getFiberDimension(p);
2210: for(int i = 0; i < dim; ++i) {
2211: array[i] += v[i];
2212: }
2213: } else {
2214: int offset = 0;
2215: int j = -1;
2217: for(int space = 0; space < this->getNumSpaces(); ++space) {
2218: const int& dim = this->getFiberDimension(p, space);
2220: for(int i = dim-1; i >= 0; --i) {
2221: array[++j] += v[i+offset];
2222: }
2223: offset += dim;
2224: }
2225: }
2226: };
2227: public: // Fibrations
2228: int getNumSpaces() const {return this->_spaces.size();};
2229: const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
2230: const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
2231: void addSpace() {
2232: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
2233: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
2234: this->_spaces.push_back(space);
2235: this->_bcs.push_back(bc);
2236: };
2237: int getFiberDimension(const point_type& p, const int space) const {
2238: return this->_spaces[space]->restrictPoint(p)->prefix;
2239: };
2240: void setFiberDimension(const point_type& p, int dim, const int space) {
2241: const index_type idx(dim, -1);
2242: this->_spaces[space]->addPoint(p);
2243: this->_spaces[space]->updatePoint(p, &idx);
2244: };
2245: template<typename Sequence>
2246: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
2247: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2248: this->setFiberDimension(*p_iter, dim, space);
2249: }
2250: }
2251: int getConstraintDimension(const point_type& p, const int space) const {
2252: if (!this->_bcs[space]->hasPoint(p)) return 0;
2253: return this->_bcs[space]->getFiberDimension(p);
2254: }
2255: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2256: this->_bcs[space]->setFiberDimension(p, numConstraints);
2257: }
2258: void addConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2259: this->_bcs[space]->addFiberDimension(p, numConstraints);
2260: }
2261: int getConstrainedFiberDimension(const point_type& p, const int space) const {
2262: return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
2263: }
2264: // Return the local dofs which are constrained on a point
2265: const int *getConstraintDof(const point_type& p, const int space) const {
2266: return this->_bcs[space]->restrictPoint(p);
2267: }
2268: // Set the local dofs which are constrained on a point
2269: void setConstraintDof(const point_type& p, const int dofs[], const int space) {
2270: this->_bcs[space]->updatePoint(p, dofs);
2271: }
2272: // Return the total number of free dofs
2273: int size(const int space) const {
2274: const chart_type& points = this->getChart();
2275: int size = 0;
2277: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2278: size += this->getConstrainedFiberDimension(*p_iter, space);
2279: }
2280: return size;
2281: };
2282: template<typename OtherSection>
2283: void copySpaces(const Obj<OtherSection>& section) {
2284: const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
2285: const std::vector<Obj<bc_type> >& bcs = section->getBCs();
2287: this->_spaces.clear();
2288: for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
2289: this->_spaces.push_back(*s_iter);
2290: }
2291: this->_bcs.clear();
2292: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
2293: this->_bcs.push_back(*b_iter);
2294: }
2295: }
2296: Obj<GeneralSection> getFibration(const int space) const {
2297: Obj<GeneralSection> field = new GeneralSection(this->comm(), this->debug());
2298: // Obj<atlas_type> _atlas;
2299: // std::vector<Obj<atlas_type> > _spaces;
2300: // Obj<bc_type> _bc;
2301: // std::vector<Obj<bc_type> > _bcs;
2302: field->addSpace();
2303: const chart_type& chart = this->getChart();
2305: // Copy sizes
2306: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2307: const int fDim = this->getFiberDimension(*c_iter, space);
2308: const int cDim = this->getConstraintDimension(*c_iter, space);
2310: if (fDim) {
2311: field->setFiberDimension(*c_iter, fDim);
2312: field->setFiberDimension(*c_iter, fDim, 0);
2313: }
2314: if (cDim) {
2315: field->setConstraintDimension(*c_iter, cDim);
2316: field->setConstraintDimension(*c_iter, cDim, 0);
2317: }
2318: }
2319: field->allocateStorage();
2320: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
2321: const chart_type& newChart = field->getChart();
2323: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2324: const int cDim = field->getConstraintDimension(*c_iter);
2325: //const int dof[1] = {0};
2327: if (cDim) {
2328: field->setConstraintDof(*c_iter, this->getConstraintDof(*c_iter, space));
2329: }
2330: }
2331: // Copy offsets
2332: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2333: index_type idx;
2335: idx.prefix = field->getFiberDimension(*c_iter);
2336: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
2337: for(int s = 0; s < space; ++s) {
2338: idx.index += this->getFiberDimension(*c_iter, s);
2339: }
2340: newAtlas->addPoint(*c_iter);
2341: newAtlas->updatePoint(*c_iter, &idx);
2342: }
2343: field->replaceStorage(this->_array, true, this->getStorageSize());
2344: field->setAtlas(newAtlas);
2345: return field;
2346: };
2347: public: // Optimization
2348: void getCustomRestrictAtlas(const int tag, const int *offsets[], const int *indices[]) {
2349: *offsets = this->_customRestrictAtlas[tag].first.first;
2350: *indices = this->_customRestrictAtlas[tag].first.second;
2351: };
2352: void getCustomUpdateAtlas(const int tag, const int *offsets[], const int *indices[]) {
2353: *offsets = this->_customUpdateAtlas[tag].first.first;
2354: *indices = this->_customUpdateAtlas[tag].first.second;
2355: };
2356: // This returns the tag assigned to the atlas
2357: int setCustomAtlas(const int restrictOffsets[], const int restrictIndices[], const int updateOffsets[], const int updateIndices[], bool autoFree = true) {
2358: this->_customRestrictAtlas.push_back(customAtlas_type(customAtlasInd_type(restrictOffsets, restrictIndices), autoFree));
2359: this->_customUpdateAtlas.push_back(customAtlas_type(customAtlasInd_type(updateOffsets, updateIndices), autoFree));
2360: return this->_customUpdateAtlas.size()-1;
2361: };
2362: int copyCustomAtlas(const Obj<GeneralSection>& section, const int tag) {
2363: const int *rOffsets, *rIndices, *uOffsets, *uIndices;
2365: section->getCustomRestrictAtlas(tag, &rOffsets, &rIndices);
2366: section->getCustomUpdateAtlas(tag, &uOffsets, &uIndices);
2367: return this->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices, false);
2368: };
2369: public:
2370: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
2371: ostringstream txt;
2372: int rank;
2374: if (comm == MPI_COMM_NULL) {
2375: comm = this->comm();
2376: rank = this->commRank();
2377: } else {
2378: MPI_Comm_rank(comm, &rank);
2379: }
2380: if (name == "") {
2381: if(rank == 0) {
2382: txt << "viewing a GeneralSection" << std::endl;
2383: }
2384: } else {
2385: if (rank == 0) {
2386: txt << "viewing GeneralSection '" << name << "'" << std::endl;
2387: }
2388: }
2389: if (rank == 0) {
2390: txt << " Fields: " << this->getNumSpaces() << std::endl;
2391: }
2392: const chart_type& chart = this->getChart();
2394: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2395: const value_type *array = this->restrictPoint(*p_iter);
2396: const int& dim = this->getFiberDimension(*p_iter);
2398: if (dim != 0) {
2399: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << " ";
2400: for(int i = 0; i < dim; i++) {
2401: txt << " " << array[i];
2402: }
2403: const int& dim = this->getConstraintDimension(*p_iter);
2405: if (dim) {
2406: const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);
2408: txt << " constrained";
2409: for(int i = 0; i < dim; ++i) {
2410: txt << " " << bcArray[i];
2411: }
2412: }
2413: txt << std::endl;
2414: }
2415: }
2416: if (chart.size() == 0) {
2417: txt << "[" << this->commRank() << "]: empty" << std::endl;
2418: }
2419: PetscSynchronizedPrintf(comm, txt.str().c_str());
2420: PetscSynchronizedFlush(comm);
2421: };
2422: };
2423: // FEMSection will support vector BC on a subset of unknowns on a point
2424: // We make a separate constraint Section to hold the transformation and projection operators
2425: // Storage will be contiguous by node, just as in Section
2426: // This allows fast restrict(p)
2427: // Then update() is accomplished by projecting constrained unknowns
2428: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
2429: typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
2430: typename BCAtlas_ = UniformSection<Point_, int, 1, typename Alloc_::template rebind<int>::other>,
2431: typename BC_ = Section<Point_, double, typename Alloc_::template rebind<double>::other> >
2432: class FEMSection : public ALE::ParallelObject {
2433: public:
2434: typedef Point_ point_type;
2435: typedef Value_ value_type;
2436: typedef Alloc_ alloc_type;
2437: typedef Atlas_ atlas_type;
2438: typedef BCAtlas_ bc_atlas_type;
2439: typedef BC_ bc_type;
2440: typedef Point index_type;
2441: typedef typename atlas_type::chart_type chart_type;
2442: typedef value_type * values_type;
2443: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
2444: typedef typename atlas_alloc_type::pointer atlas_ptr;
2445: typedef typename alloc_type::template rebind<bc_type>::other bc_atlas_alloc_type;
2446: typedef typename bc_atlas_alloc_type::pointer bc_atlas_ptr;
2447: typedef typename alloc_type::template rebind<bc_type>::other bc_alloc_type;
2448: typedef typename bc_alloc_type::pointer bc_ptr;
2449: protected:
2450: Obj<atlas_type> _atlas;
2451: Obj<bc_atlas_type> _bc_atlas;
2452: Obj<bc_type> _bc;
2453: values_type _array;
2454: bool _sharedStorage;
2455: int _sharedStorageSize;
2456: alloc_type _allocator;
2457: std::vector<Obj<atlas_type> > _spaces;
2458: std::vector<Obj<bc_type> > _bcs;
2459: public:
2460: FEMSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
2461: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
2462: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
2463: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
2464: bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2465: bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(comm, debug));
2466: this->_bc_atlas = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2467: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
2468: bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
2469: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
2470: this->_array = NULL;
2471: this->_sharedStorage = false;
2472: };
2473: FEMSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
2474: bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2475: bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(this->comm(), this->debug()));
2476: this->_bc_atlas = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2477: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
2478: bc_alloc_type(this->_allocator).construct(pBC, bc_type(this->comm(), this->debug()));
2479: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
2480: };
2481: ~FEMSection() {
2482: if (this->_array && !this->_sharedStorage) {
2483: const int totalSize = this->sizeWithBC();
2485: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2486: this->_allocator.deallocate(this->_array, totalSize);
2487: this->_array = NULL;
2488: }
2489: };
2490: public:
2491: value_type *getRawArray(const int size) {
2492: // Put in a sentinel value that deallocates the array
2493: static value_type *array = NULL;
2494: static int maxSize = 0;
2496: if (size > maxSize) {
2497: const value_type dummy(0);
2499: if (array) {
2500: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
2501: this->_allocator.deallocate(array, maxSize);
2502: }
2503: maxSize = size;
2504: array = this->_allocator.allocate(maxSize);
2505: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
2506: };
2507: return array;
2508: };
2509: int getStorageSize() const {
2510: if (this->_sharedStorage) {
2511: return this->_sharedStorageSize;
2512: }
2513: return this->sizeWithBC();
2514: };
2515: public: // Verifiers
2516: bool hasPoint(const point_type& point) {
2517: return this->_atlas->hasPoint(point);
2518: };
2519: public: // Accessors
2520: const chart_type& getChart() const {return this->_atlas->getChart();};
2521: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
2522: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
2523: const Obj<bc_type>& getBC() const {return this->_bc;};
2524: void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
2525: public: // BC
2526: // Returns the number of constraints on a point
2527: int getConstraintDimension(const point_type& p) const {
2528: if (!this->_bc_atlas->hasPoint(p)) return 0;
2529: return this->_bc_atlas->restrictPoint(p)[0];
2530: };
2531: // Set the number of constraints on a point
2532: void setConstraintDimension(const point_type& p, const int numConstraints) {
2533: this->_bc_atlas->updatePoint(p, &numConstraints);
2534: };
2535: // Increment the number of constraints on a point
2536: void addConstraintDimension(const point_type& p, const int numConstraints) {
2537: this->_bc_atlas->updatePointAdd(p, &numConstraints);
2538: };
2539: const int *getConstraintDof(const point_type& p) const {
2540: return this->_bc->restrictPoint(p);
2541: }
2542: public: // Sizes
2543: void clear() {
2544: if (!this->_sharedStorage) {
2545: const int totalSize = this->sizeWithBC();
2547: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2548: this->_allocator.deallocate(this->_array, totalSize);
2549: }
2550: this->_array = NULL;
2551: this->_atlas->clear();
2552: this->_bc_atlas->clear();
2553: this->_bc->clear();
2554: };
2555: // Return the total number of dofs on the point (free and constrained)
2556: int getFiberDimension(const point_type& p) const {
2557: return this->_atlas->restrictPoint(p)->prefix;
2558: };
2559: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
2560: return atlas->restrictPoint(p)->prefix;
2561: };
2562: // Set the total number of dofs on the point (free and constrained)
2563: void setFiberDimension(const point_type& p, int dim) {
2564: const index_type idx(dim, -1);
2565: this->_atlas->addPoint(p);
2566: this->_atlas->updatePoint(p, &idx);
2567: };
2568: template<typename Sequence>
2569: void setFiberDimension(const Obj<Sequence>& points, int dim) {
2570: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2571: this->setFiberDimension(*p_iter, dim);
2572: }
2573: }
2574: void addFiberDimension(const point_type& p, int dim) {
2575: if (this->_atlas->hasPoint(p)) {
2576: const index_type values(dim, 0);
2577: this->_atlas->updateAddPoint(p, &values);
2578: } else {
2579: this->setFiberDimension(p, dim);
2580: }
2581: };
2582: // Return the number of constrained dofs on this point
2583: int getConstrainedFiberDimension(const point_type& p) const {
2584: return this->getFiberDimension(p) - this->getConstraintDimension(p);
2585: };
2586: // Return the total number of free dofs
2587: int size() const {
2588: const chart_type& points = this->getChart();
2589: int size = 0;
2591: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2592: size += this->getConstrainedFiberDimension(*p_iter);
2593: }
2594: return size;
2595: };
2596: // Return the total number of dofs (free and constrained)
2597: int sizeWithBC() const {
2598: const chart_type& points = this->getChart();
2599: int size = 0;
2601: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2602: size += this->getFiberDimension(*p_iter);
2603: }
2604: return size;
2605: };
2606: public: // Allocation
2607: void allocateStorage() {
2608: const int totalSize = this->sizeWithBC();
2609: const value_type dummy(0) ;
2611: this->_array = this->_allocator.allocate(totalSize);
2612: this->_sharedStorage = false;
2613: this->_sharedStorageSize = 0;
2614: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
2615: this->_bc_atlas->allocatePoint();
2616: this->_bc->allocatePoint();
2617: };
2618: void orderPoints(const Obj<atlas_type>& atlas){
2619: const chart_type& chart = this->getChart();
2620: int offset = 0;
2622: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2623: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
2624: const int& dim = idx.prefix;
2626: idx.index = offset;
2627: atlas->updatePoint(*c_iter, &idx);
2628: offset += dim;
2629: }
2630: };
2631: void allocatePoint() {
2632: this->orderPoints(this->_atlas);
2633: this->allocateStorage();
2634: };
2635: public: // Restriction and Update
2636: // Zero entries
2637: void zero() {
2638: this->set(0.0);
2639: };
2640: void set(const value_type value) {
2641: //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
2642: const chart_type& chart = this->getChart();
2644: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2645: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2646: const int& dim = this->getFiberDimension(*c_iter);
2647: const int& cDim = this->getConstraintDimension(*c_iter);
2649: if (!cDim) {
2650: for(int i = 0; i < dim; ++i) {
2651: array[i] = value;
2652: }
2653: } else {
2654: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2655: int cInd = 0;
2657: for(int i = 0; i < dim; ++i) {
2658: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2659: array[i] = value;
2660: }
2661: }
2662: }
2663: };
2664: // Add two sections and put the result in a third
2665: void add(const Obj<FEMSection>& x, const Obj<FEMSection>& y) {
2666: // Check atlases
2667: const chart_type& chart = this->getChart();
2669: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2670: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2671: const value_type *xArray = x->restrictPoint(*c_iter);
2672: const value_type *yArray = y->restrictPoint(*c_iter);
2673: const int& dim = this->getFiberDimension(*c_iter);
2674: const int& cDim = this->getConstraintDimension(*c_iter);
2676: if (!cDim) {
2677: for(int i = 0; i < dim; ++i) {
2678: array[i] = xArray[i] + yArray[i];
2679: }
2680: } else {
2681: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2682: int cInd = 0;
2684: for(int i = 0; i < dim; ++i) {
2685: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2686: array[i] = xArray[i] + yArray[i];
2687: }
2688: }
2689: }
2690: };
2691: // this = this + alpha * x
2692: void axpy(const value_type alpha, const Obj<FEMSection>& x) {
2693: // Check atlases
2694: const chart_type& chart = this->getChart();
2696: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2697: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2698: const value_type *xArray = x->restrictPoint(*c_iter);
2699: const int& dim = this->getFiberDimension(*c_iter);
2700: const int& cDim = this->getConstraintDimension(*c_iter);
2702: if (!cDim) {
2703: for(int i = 0; i < dim; ++i) {
2704: array[i] += alpha*xArray[i];
2705: }
2706: } else {
2707: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2708: int cInd = 0;
2710: for(int i = 0; i < dim; ++i) {
2711: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2712: array[i] += alpha*xArray[i];
2713: }
2714: }
2715: }
2716: };
2717: // Return the free values on a point
2718: const value_type *restrictSpace() const {
2719: return this->_array;
2720: };
2721: // Return the free values on a point
2722: const value_type *restrictPoint(const point_type& p) const {
2723: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
2724: };
2725: // Update the free values on a point
2726: // Takes a full array and ignores constrained values
2727: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2728: value_type *array = (value_type *) this->restrictPoint(p);
2729: const int& cDim = this->getConstraintDimension(p);
2731: if (!cDim) {
2732: if (orientation >= 0) {
2733: const int& dim = this->getFiberDimension(p);
2735: for(int i = 0; i < dim; ++i) {
2736: array[i] = v[i];
2737: }
2738: } else {
2739: int offset = 0;
2740: int j = -1;
2742: for(int space = 0; space < this->getNumSpaces(); ++space) {
2743: const int& dim = this->getFiberDimension(p, space);
2745: for(int i = dim-1; i >= 0; --i) {
2746: array[++j] = v[i+offset];
2747: }
2748: offset += dim;
2749: }
2750: }
2751: } else {
2752: if (orientation >= 0) {
2753: const int& dim = this->getFiberDimension(p);
2754: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2755: int cInd = 0;
2757: for(int i = 0; i < dim; ++i) {
2758: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2759: array[i] = v[i];
2760: }
2761: } else {
2762: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2763: int offset = 0;
2764: int cOffset = 0;
2765: int j = 0;
2767: for(int space = 0; space < this->getNumSpaces(); ++space) {
2768: const int dim = this->getFiberDimension(p, space);
2769: const int tDim = this->getConstrainedFiberDimension(p, space);
2770: const int sDim = dim - tDim;
2771: int cInd = 0;
2773: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2774: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2775: array[j] = v[k];
2776: }
2777: offset += dim;
2778: cOffset += dim - tDim;
2779: }
2780: }
2781: }
2782: };
2783: // Update the free values on a point
2784: // Takes a full array and ignores constrained values
2785: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2786: value_type *array = (value_type *) this->restrictPoint(p);
2787: const int& cDim = this->getConstraintDimension(p);
2789: if (!cDim) {
2790: if (orientation >= 0) {
2791: const int& dim = this->getFiberDimension(p);
2793: for(int i = 0; i < dim; ++i) {
2794: array[i] += v[i];
2795: }
2796: } else {
2797: int offset = 0;
2798: int j = -1;
2800: for(int space = 0; space < this->getNumSpaces(); ++space) {
2801: const int& dim = this->getFiberDimension(p, space);
2803: for(int i = dim-1; i >= 0; --i) {
2804: array[++j] += v[i+offset];
2805: }
2806: offset += dim;
2807: }
2808: }
2809: } else {
2810: if (orientation >= 0) {
2811: const int& dim = this->getFiberDimension(p);
2812: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2813: int cInd = 0;
2815: for(int i = 0; i < dim; ++i) {
2816: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2817: array[i] += v[i];
2818: }
2819: } else {
2820: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2821: int offset = 0;
2822: int cOffset = 0;
2823: int j = 0;
2825: for(int space = 0; space < this->getNumSpaces(); ++space) {
2826: const int dim = this->getFiberDimension(p, space);
2827: const int tDim = this->getConstrainedFiberDimension(p, space);
2828: const int sDim = dim - tDim;
2829: int cInd = 0;
2831: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2832: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2833: array[j] += v[k];
2834: }
2835: offset += dim;
2836: cOffset += dim - tDim;
2837: }
2838: }
2839: }
2840: };
2841: // Update the free values on a point
2842: // Takes ONLY unconstrained values
2843: void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2844: value_type *array = (value_type *) this->restrictPoint(p);
2845: const int& cDim = this->getConstraintDimension(p);
2847: if (!cDim) {
2848: if (orientation >= 0) {
2849: const int& dim = this->getFiberDimension(p);
2851: for(int i = 0; i < dim; ++i) {
2852: array[i] = v[i];
2853: }
2854: } else {
2855: int offset = 0;
2856: int j = -1;
2858: for(int space = 0; space < this->getNumSpaces(); ++space) {
2859: const int& dim = this->getFiberDimension(p, space);
2861: for(int i = dim-1; i >= 0; --i) {
2862: array[++j] = v[i+offset];
2863: }
2864: offset += dim;
2865: }
2866: }
2867: } else {
2868: if (orientation >= 0) {
2869: const int& dim = this->getFiberDimension(p);
2870: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2871: int cInd = 0;
2873: for(int i = 0, k = -1; i < dim; ++i) {
2874: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2875: array[i] = v[++k];
2876: }
2877: } else {
2878: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2879: int offset = 0;
2880: int cOffset = 0;
2881: int j = 0;
2883: for(int space = 0; space < this->getNumSpaces(); ++space) {
2884: const int dim = this->getFiberDimension(p, space);
2885: const int tDim = this->getConstrainedFiberDimension(p, space);
2886: const int sDim = dim - tDim;
2887: int cInd = 0;
2889: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2890: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2891: array[j] = v[--k];
2892: }
2893: offset += dim;
2894: cOffset += dim - tDim;
2895: }
2896: }
2897: }
2898: };
2899: // Update the free values on a point
2900: // Takes ONLY unconstrained values
2901: void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2902: value_type *array = (value_type *) this->restrictPoint(p);
2903: const int& cDim = this->getConstraintDimension(p);
2905: if (!cDim) {
2906: if (orientation >= 0) {
2907: const int& dim = this->getFiberDimension(p);
2909: for(int i = 0; i < dim; ++i) {
2910: array[i] += v[i];
2911: }
2912: } else {
2913: int offset = 0;
2914: int j = -1;
2916: for(int space = 0; space < this->getNumSpaces(); ++space) {
2917: const int& dim = this->getFiberDimension(p, space);
2919: for(int i = dim-1; i >= 0; --i) {
2920: array[++j] += v[i+offset];
2921: }
2922: offset += dim;
2923: }
2924: }
2925: } else {
2926: if (orientation >= 0) {
2927: const int& dim = this->getFiberDimension(p);
2928: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2929: int cInd = 0;
2931: for(int i = 0, k = -1; i < dim; ++i) {
2932: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2933: array[i] += v[++k];
2934: }
2935: } else {
2936: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2937: int offset = 0;
2938: int cOffset = 0;
2939: int j = 0;
2941: for(int space = 0; space < this->getNumSpaces(); ++space) {
2942: const int dim = this->getFiberDimension(p, space);
2943: const int tDim = this->getConstrainedFiberDimension(p, space);
2944: const int sDim = dim - tDim;
2945: int cInd = 0;
2947: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2948: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2949: array[j] += v[--k];
2950: }
2951: offset += dim;
2952: cOffset += dim - tDim;
2953: }
2954: }
2955: }
2956: };
2957: // Update only the constrained dofs on a point
2958: // This takes an array with ONLY bc values
2959: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2960: value_type *array = (value_type *) this->restrictPoint(p);
2961: const int& cDim = this->getConstraintDimension(p);
2963: if (cDim) {
2964: if (orientation >= 0) {
2965: const int& dim = this->getFiberDimension(p);
2966: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2967: int cInd = 0;
2969: for(int i = 0; i < dim; ++i) {
2970: if (cInd == cDim) break;
2971: if (i == cDof[cInd]) {
2972: array[i] = v[cInd];
2973: ++cInd;
2974: }
2975: }
2976: } else {
2977: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2978: int cOffset = 0;
2979: int j = 0;
2981: for(int space = 0; space < this->getNumSpaces(); ++space) {
2982: const int dim = this->getFiberDimension(p, space);
2983: const int tDim = this->getConstrainedFiberDimension(p, space);
2984: int cInd = 0;
2986: for(int i = 0; i < dim; ++i, ++j) {
2987: if (cInd < 0) break;
2988: if (j == cDof[cInd+cOffset]) {
2989: array[j] = v[cInd+cOffset];
2990: ++cInd;
2991: }
2992: }
2993: cOffset += dim - tDim;
2994: }
2995: }
2996: }
2997: };
2998: // Update only the constrained dofs on a point
2999: // This takes an array with ALL values, not just BC
3000: void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
3001: value_type *array = (value_type *) this->restrictPoint(p);
3002: const int& cDim = this->getConstraintDimension(p);
3004: if (cDim) {
3005: if (orientation >= 0) {
3006: const int& dim = this->getFiberDimension(p);
3007: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
3008: int cInd = 0;
3010: for(int i = 0; i < dim; ++i) {
3011: if (cInd == cDim) break;
3012: if (i == cDof[cInd]) {
3013: array[i] = v[i];
3014: ++cInd;
3015: }
3016: }
3017: } else {
3018: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
3019: int offset = 0;
3020: int cOffset = 0;
3021: int j = 0;
3023: for(int space = 0; space < this->getNumSpaces(); ++space) {
3024: const int dim = this->getFiberDimension(p, space);
3025: const int tDim = this->getConstrainedFiberDimension(p, space);
3026: int cInd = 0;
3028: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
3029: if (cInd < 0) break;
3030: if (j == cDof[cInd+cOffset]) {
3031: array[j] = v[k];
3032: ++cInd;
3033: }
3034: }
3035: offset += dim;
3036: cOffset += dim - tDim;
3037: }
3038: }
3039: }
3040: };
3041: // Update all dofs on a point (free and constrained)
3042: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
3043: value_type *array = (value_type *) this->restrictPoint(p);
3045: if (orientation >= 0) {
3046: const int& dim = this->getFiberDimension(p);
3048: for(int i = 0; i < dim; ++i) {
3049: array[i] = v[i];
3050: }
3051: } else {
3052: int offset = 0;
3053: int j = -1;
3055: for(int space = 0; space < this->getNumSpaces(); ++space) {
3056: const int& dim = this->getFiberDimension(p, space);
3058: for(int i = dim-1; i >= 0; --i) {
3059: array[++j] = v[i+offset];
3060: }
3061: offset += dim;
3062: }
3063: }
3064: };
3065: // Update all dofs on a point (free and constrained)
3066: void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
3067: value_type *array = (value_type *) this->restrictPoint(p);
3069: if (orientation >= 0) {
3070: const int& dim = this->getFiberDimension(p);
3072: for(int i = 0; i < dim; ++i) {
3073: array[i] += v[i];
3074: }
3075: } else {
3076: int offset = 0;
3077: int j = -1;
3079: for(int space = 0; space < this->getNumSpaces(); ++space) {
3080: const int& dim = this->getFiberDimension(p, space);
3082: for(int i = dim-1; i >= 0; --i) {
3083: array[++j] += v[i+offset];
3084: }
3085: offset += dim;
3086: }
3087: }
3088: };
3089: public: // Fibrations
3090: int getNumSpaces() const {return this->_spaces.size();};
3091: const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
3092: const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
3093: void addSpace() {
3094: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
3095: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
3096: this->_spaces.push_back(space);
3097: this->_bcs.push_back(bc);
3098: };
3099: int getFiberDimension(const point_type& p, const int space) const {
3100: return this->_spaces[space]->restrictPoint(p)->prefix;
3101: };
3102: void setFiberDimension(const point_type& p, int dim, const int space) {
3103: const index_type idx(dim, -1);
3104: this->_spaces[space]->addPoint(p);
3105: this->_spaces[space]->updatePoint(p, &idx);
3106: };
3107: template<typename Sequence>
3108: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
3109: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
3110: this->setFiberDimension(*p_iter, dim, space);
3111: }
3112: }
3113: int getConstraintDimension(const point_type& p, const int space) const {
3114: if (!this->_bcs[space]->hasPoint(p)) return 0;
3115: return this->_bcs[space]->getFiberDimension(p);
3116: };
3117: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
3118: this->_bcs[space]->setFiberDimension(p, numConstraints);
3119: };
3120: int getConstrainedFiberDimension(const point_type& p, const int space) const {
3121: return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
3122: };
3123: void copyFibration(const Obj<FEMSection>& section) {
3124: const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
3125: const std::vector<Obj<bc_type> >& bcs = section->getBCs();
3127: this->_spaces.clear();
3128: for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
3129: this->_spaces.push_back(*s_iter);
3130: }
3131: this->_bcs.clear();
3132: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
3133: this->_bcs.push_back(*b_iter);
3134: }
3135: };
3136: Obj<FEMSection> getFibration(const int space) const {
3137: Obj<FEMSection> field = new FEMSection(this->comm(), this->debug());
3138: // Obj<atlas_type> _atlas;
3139: // std::vector<Obj<atlas_type> > _spaces;
3140: // Obj<bc_type> _bc;
3141: // std::vector<Obj<bc_type> > _bcs;
3142: field->addSpace();
3143: const chart_type& chart = this->getChart();
3145: // Copy sizes
3146: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3147: const int fDim = this->getFiberDimension(*c_iter, space);
3148: const int cDim = this->getConstraintDimension(*c_iter, space);
3150: if (fDim) {
3151: field->setFiberDimension(*c_iter, fDim);
3152: field->setFiberDimension(*c_iter, fDim, 0);
3153: }
3154: if (cDim) {
3155: field->setConstraintDimension(*c_iter, cDim);
3156: field->setConstraintDimension(*c_iter, cDim, 0);
3157: }
3158: }
3159: field->allocateStorage();
3160: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
3161: const chart_type& newChart = field->getChart();
3163: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3164: const int cDim = field->getConstraintDimension(*c_iter);
3165: const int dof[1] = {0};
3167: if (cDim) {
3168: field->setConstraintDof(*c_iter, dof);
3169: }
3170: }
3171: // Copy offsets
3172: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3173: index_type idx;
3175: idx.prefix = field->getFiberDimension(*c_iter);
3176: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
3177: for(int s = 0; s < space; ++s) {
3178: idx.index += this->getFiberDimension(*c_iter, s);
3179: }
3180: newAtlas->addPoint(*c_iter);
3181: newAtlas->updatePoint(*c_iter, &idx);
3182: }
3183: field->replaceStorage(this->_array, true, this->getStorageSize());
3184: field->setAtlas(newAtlas);
3185: return field;
3186: };
3187: public:
3188: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3189: ostringstream txt;
3190: int rank;
3192: if (comm == MPI_COMM_NULL) {
3193: comm = this->comm();
3194: rank = this->commRank();
3195: } else {
3196: MPI_Comm_rank(comm, &rank);
3197: }
3198: if (name == "") {
3199: if(rank == 0) {
3200: txt << "viewing a FEMSection" << std::endl;
3201: }
3202: } else {
3203: if (rank == 0) {
3204: txt << "viewing FEMSection '" << name << "'" << std::endl;
3205: }
3206: }
3207: if (rank == 0) {
3208: txt << " Fields: " << this->getNumSpaces() << std::endl;
3209: }
3210: const chart_type& chart = this->getChart();
3212: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
3213: const value_type *array = this->restrictPoint(*p_iter);
3214: const int& dim = this->getFiberDimension(*p_iter);
3216: if (dim != 0) {
3217: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << " ";
3218: for(int i = 0; i < dim; i++) {
3219: txt << " " << array[i];
3220: }
3221: const int& dim = this->getConstraintDimension(*p_iter);
3223: if (dim) {
3224: const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);
3226: txt << " constrained";
3227: for(int i = 0; i < dim; ++i) {
3228: txt << " " << bcArray[i];
3229: }
3230: }
3231: txt << std::endl;
3232: }
3233: }
3234: if (chart.size() == 0) {
3235: txt << "[" << this->commRank() << "]: empty" << std::endl;
3236: }
3237: PetscSynchronizedPrintf(comm, txt.str().c_str());
3238: PetscSynchronizedFlush(comm);
3239: };
3240: };
3241: // A Field combines several sections
3242: template<typename Overlap_, typename Patch_, typename Section_>
3243: class Field : public ALE::ParallelObject {
3244: public:
3245: typedef Overlap_ overlap_type;
3246: typedef Patch_ patch_type;
3247: typedef Section_ section_type;
3248: typedef typename section_type::point_type point_type;
3249: typedef typename section_type::chart_type chart_type;
3250: typedef typename section_type::value_type value_type;
3251: typedef std::map<patch_type, Obj<section_type> > sheaf_type;
3252: typedef enum {SEND, RECEIVE} request_type;
3253: typedef std::map<patch_type, MPI_Request> requests_type;
3254: protected:
3255: sheaf_type _sheaf;
3256: int _tag;
3257: MPI_Datatype _datatype;
3258: requests_type _requests;
3259: public:
3260: Field(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
3261: this->_tag = this->getNewTag();
3262: this->_datatype = this->getMPIDatatype();
3263: };
3264: Field(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _tag(tag) {
3265: this->_datatype = this->getMPIDatatype();
3266: };
3267: virtual ~Field() {};
3268: protected:
3269: MPI_Datatype getMPIDatatype() {
3270: if (sizeof(value_type) == 4) {
3271: return MPI_INT;
3272: } else if (sizeof(value_type) == 8) {
3273: return MPI_DOUBLE;
3274: } else if (sizeof(value_type) == 28) {
3275: int blen[2];
3276: MPI_Aint indices[2];
3277: MPI_Datatype oldtypes[2], newtype;
3278: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
3279: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3280: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3281: MPI_Type_commit(&newtype);
3282: return newtype;
3283: } else if (sizeof(value_type) == 32) {
3284: int blen[2];
3285: MPI_Aint indices[2];
3286: MPI_Datatype oldtypes[2], newtype;
3287: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_DOUBLE;
3288: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3289: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3290: MPI_Type_commit(&newtype);
3291: return newtype;
3292: }
3293: throw ALE::Exception("Cannot determine MPI type for value type");
3294: };
3295: int getNewTag() {
3296: static int tagKeyval = MPI_KEYVAL_INVALID;
3297: int *tagvalp = NULL, *maxval, flg;
3299: if (tagKeyval == MPI_KEYVAL_INVALID) {
3300: tagvalp = (int *) malloc(sizeof(int));
3301: MPI_Keyval_create(MPI_NULL_COPY_FN, DMMesh_DelTag, &tagKeyval, (void *) NULL);
3302: MPI_Attr_put(this->_comm, tagKeyval, tagvalp);
3303: tagvalp[0] = 0;
3304: }
3305: MPI_Attr_get(this->_comm, tagKeyval, (void **) &tagvalp, &flg);
3306: if (tagvalp[0] < 1) {
3307: MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
3308: tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
3309: }
3310: if (this->debug()) {
3311: std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
3312: }
3313: return tagvalp[0]--;
3314: };
3315: public: // Verifiers
3316: void checkPatch(const patch_type& patch) const {
3317: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3318: ostringstream msg;
3319: msg << "Invalid field patch " << patch << std::endl;
3320: throw ALE::Exception(msg.str().c_str());
3321: }
3322: };
3323: bool hasPatch(const patch_type& patch) {
3324: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3325: return false;
3326: }
3327: return true;
3328: };
3329: public: // Accessors
3330: int getTag() const {return this->_tag;};
3331: void setTag(const int tag) {this->_tag = tag;};
3332: Obj<section_type>& getSection(const patch_type& patch) {
3333: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3334: this->_sheaf[patch] = new section_type(this->comm(), this->debug());
3335: }
3336: return this->_sheaf[patch];
3337: };
3338: void setSection(const patch_type& patch, const Obj<section_type>& section) {this->_sheaf[patch] = section;};
3339: const sheaf_type& getPatches() {
3340: return this->_sheaf;
3341: };
3342: void clear() {
3343: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3344: p_iter->second->clear();
3345: }
3346: };
3347: public: // Adapter
3348: template<typename Topology_>
3349: void setTopology(const Obj<Topology_>& topology) {
3350: const typename Topology_::sheaf_type& patches = topology->getPatches();
3352: for(typename Topology_::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3353: int rank = p_iter->first;
3354: const Obj<section_type>& section = this->getSection(rank);
3355: const Obj<typename Topology_::sieve_type::baseSequence>& base = p_iter->second->base();
3357: for(typename Topology_::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
3358: section->setFiberDimension(*b_iter, 1);
3359: }
3360: }
3361: }
3362: void allocate() {
3363: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3364: p_iter->second->allocatePoint();
3365: }
3366: }
3367: public: // Communication
3368: void construct(const int size) {
3369: const sheaf_type& patches = this->getPatches();
3371: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3372: const patch_type rank = p_iter->first;
3373: const Obj<section_type>& section = this->getSection(rank);
3374: const chart_type& chart = section->getChart();
3376: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3377: section->setFiberDimension(*c_iter, size);
3378: }
3379: }
3380: };
3381: template<typename Sizer>
3382: void construct(const Obj<Sizer>& sizer) {
3383: const sheaf_type& patches = this->getPatches();
3385: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3386: const patch_type rank = p_iter->first;
3387: const Obj<section_type>& section = this->getSection(rank);
3388: const chart_type& chart = section->getChart();
3389:
3390: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3391: section->setFiberDimension(*c_iter, *(sizer->getSection(rank)->restrictPoint(*c_iter)));
3392: }
3393: }
3394: }
3395: void constructCommunication(const request_type& requestType) {
3396: const sheaf_type& patches = this->getPatches();
3398: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3399: const patch_type patch = p_iter->first;
3400: const Obj<section_type>& section = this->getSection(patch);
3401: MPI_Request request;
3403: if (requestType == RECEIVE) {
3404: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data(" << section->size() << ") from " << patch << " tag " << this->_tag << std::endl;}
3405: MPI_Recv_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3406: } else {
3407: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << section->size() << ") to " << patch << " tag " << this->_tag << std::endl;}
3408: MPI_Send_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3409: }
3410: this->_requests[patch] = request;
3411: }
3412: };
3413: void startCommunication() {
3414: const sheaf_type& patches = this->getPatches();
3416: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3417: MPI_Request request = this->_requests[p_iter->first];
3419: MPI_Start(&request);
3420: }
3421: };
3422: void endCommunication() {
3423: const sheaf_type& patches = this->getPatches();
3424: MPI_Status status;
3426: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3427: MPI_Request request = this->_requests[p_iter->first];
3429: MPI_Wait(&request, &status);
3430: }
3431: };
3432: public:
3433: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3434: ostringstream txt;
3435: int rank;
3437: if (comm == MPI_COMM_NULL) {
3438: comm = this->comm();
3439: rank = this->commRank();
3440: } else {
3441: MPI_Comm_rank(comm, &rank);
3442: }
3443: if (name == "") {
3444: if(rank == 0) {
3445: txt << "viewing a Field" << std::endl;
3446: }
3447: } else {
3448: if(rank == 0) {
3449: txt << "viewing Field '" << name << "'" << std::endl;
3450: }
3451: }
3452: PetscSynchronizedPrintf(comm, txt.str().c_str());
3453: PetscSynchronizedFlush(comm);
3454: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3455: ostringstream txt1;
3457: txt1 << "[" << this->commRank() << "]: Patch " << p_iter->first << std::endl;
3458: PetscSynchronizedPrintf(comm, txt1.str().c_str());
3459: PetscSynchronizedFlush(comm);
3460: p_iter->second->view("field section", comm);
3461: }
3462: };
3463: };
3464: }
3465: #endif