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 <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: public: // Verifiers
1457: bool hasPoint(const point_type& point) {
1458: return this->_atlas->hasPoint(point);
1459: };
1460: public: // Accessors
1461: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
1462: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1463: const Obj<atlas_type>& getNewAtlas() const {return this->_atlasNew;};
1464: void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1465: const Obj<bc_type>& getBC() const {return this->_bc;};
1466: void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
1467: const chart_type& getChart() const {return this->_atlas->getChart();};
1468: void setChart(const chart_type& chart) {throw ALE::Exception("setChart() for GeneralSection is invalid");};
1469: public: // BC
1470: // Returns the number of constraints on a point
1471: int getConstraintDimension(const point_type& p) const {
1472: if (!this->_bc->hasPoint(p)) return 0;
1473: return this->_bc->getFiberDimension(p);
1474: }
1475: // Set the number of constraints on a point
1476: void setConstraintDimension(const point_type& p, const int numConstraints) {
1477: this->_bc->setFiberDimension(p, numConstraints);
1478: }
1479: // Increment the number of constraints on a point
1480: void addConstraintDimension(const point_type& p, const int numConstraints) {
1481: this->_bc->addFiberDimension(p, numConstraints);
1482: }
1483: // Return the local dofs which are constrained on a point
1484: const int *getConstraintDof(const point_type& p) const {
1485: return this->_bc->restrictPoint(p);
1486: }
1487: // Set the local dofs which are constrained on a point
1488: void setConstraintDof(const point_type& p, const int dofs[]) {
1489: this->_bc->updatePoint(p, dofs);
1490: }
1491: template<typename OtherSection>
1492: void copyBC(const Obj<OtherSection>& section) {
1493: this->setBC(section->getBC());
1494: const chart_type& chart = this->getChart();
1496: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1497: if (this->getConstraintDimension(*p_iter)) {
1498: this->updatePointBCFull(*p_iter, section->restrictPoint(*p_iter));
1499: }
1500: }
1501: this->copyFibration(section);
1502: }
1503: void defaultConstraintDof() {
1504: const chart_type& chart = this->getChart();
1505: int size = 0;
1507: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1508: size = std::max(size, this->getConstraintDimension(*p_iter));
1509: }
1510: int *dofs = new int[size];
1511: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1512: const int cDim = this->getConstraintDimension(*p_iter);
1514: if (cDim) {
1515: for(int d = 0; d < cDim; ++d) {
1516: dofs[d] = d;
1517: }
1518: this->_bc->updatePoint(*p_iter, dofs);
1519: }
1520: }
1521: delete [] dofs;
1522: };
1523: public: // Sizes
1524: void clear() {
1525: if (!this->_sharedStorage) {
1526: const int totalSize = this->sizeWithBC();
1528: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1529: this->_allocator.deallocate(this->_array, totalSize);
1530: ///delete [] this->_array;
1531: }
1532: this->_array = NULL;
1533: this->_atlas->clear();
1534: this->_bc->clear();
1535: };
1536: // Return the total number of dofs on the point (free and constrained)
1537: int getFiberDimension(const point_type& p) const {
1538: return this->_atlas->restrictPoint(p)->prefix;
1539: };
1540: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
1541: return atlas->restrictPoint(p)->prefix;
1542: };
1543: // Set the total number of dofs on the point (free and constrained)
1544: void setFiberDimension(const point_type& p, int dim) {
1545: const index_type idx(dim, -1);
1546: this->_atlas->addPoint(p);
1547: this->_atlas->updatePoint(p, &idx);
1548: };
1549: template<typename Sequence>
1550: void setFiberDimension(const Obj<Sequence>& points, int dim) {
1551: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1552: this->setFiberDimension(*p_iter, dim);
1553: }
1554: }
1555: void addFiberDimension(const point_type& p, int dim) {
1556: if (this->_atlas->hasPoint(p)) {
1557: const index_type values(dim, 0);
1558: this->_atlas->updateAddPoint(p, &values);
1559: } else {
1560: this->setFiberDimension(p, dim);
1561: }
1562: };
1563: // Return the number of constrained dofs on this point
1564: int getConstrainedFiberDimension(const point_type& p) const {
1565: return this->getFiberDimension(p) - this->getConstraintDimension(p);
1566: };
1567: // Return the total number of free dofs
1568: int size() const {
1569: const chart_type& points = this->getChart();
1570: int size = 0;
1572: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1573: size += this->getConstrainedFiberDimension(*p_iter);
1574: }
1575: return size;
1576: };
1577: // Return the total number of dofs (free and constrained)
1578: int sizeWithBC() const {
1579: const chart_type& points = this->getChart();
1580: int size = 0;
1582: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1583: size += this->getFiberDimension(*p_iter);
1584: }
1585: return size;
1586: };
1587: public: // Index retrieval
1588: const typename index_type::index_type& getIndex(const point_type& p) {
1589: return this->_atlas->restrictPoint(p)->index;
1590: };
1591: void setIndex(const point_type& p, const typename index_type::index_type& index) {
1592: ((typename atlas_type::value_type *) this->_atlas->restrictPoint(p))->index = index;
1593: };
1594: void setIndexBC(const point_type& p, const typename index_type::index_type& index) {};
1595: void getIndices(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = true) {
1596: this->getIndices(p, this->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1597: };
1598: template<typename Order_>
1599: 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) {
1600: this->getIndices(p, order->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1601: }
1602: void getIndicesRaw(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation) {
1603: if (orientation >= 0) {
1604: const int& dim = this->getFiberDimension(p);
1605: const int end = start + dim;
1607: for(int i = start; i < end; ++i) {
1608: indices[(*indx)++] = i;
1609: }
1610: } else {
1611: const int numSpaces = this->getNumSpaces();
1612: int offset = start;
1614: for(int space = 0; space < numSpaces; ++space) {
1615: const int& dim = this->getFiberDimension(p, space);
1617: for(int i = offset+dim-1; i >= offset; --i) {
1618: indices[(*indx)++] = i;
1619: }
1620: offset += dim;
1621: }
1622: if (!numSpaces) {
1623: const int& dim = this->getFiberDimension(p);
1625: for(int i = offset+dim-1; i >= offset; --i) {
1626: indices[(*indx)++] = i;
1627: }
1628: offset += dim;
1629: }
1630: }
1631: }
1632: 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) {
1633: const int& cDim = this->getConstraintDimension(p);
1635: if (!cDim) {
1636: this->getIndicesRaw(p, start, indices, indx, orientation);
1637: } else {
1638: if (orientation >= 0) {
1639: const int& dim = this->getFiberDimension(p);
1640: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1641: int cInd = 0;
1643: for(int i = start, k = 0; k < dim; ++k) {
1644: if ((cInd < cDim) && (k == cDof[cInd])) {
1645: if (!freeOnly) indices[(*indx)++] = -(k+1);
1646: if (skipConstraints) ++i;
1647: ++cInd;
1648: } else {
1649: indices[(*indx)++] = i++;
1650: }
1651: }
1652: } else {
1653: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1654: int offset = 0;
1655: int cOffset = 0;
1656: int j = -1;
1658: for(int space = 0; space < this->getNumSpaces(); ++space) {
1659: const int dim = this->getFiberDimension(p, space);
1660: const int tDim = this->getConstrainedFiberDimension(p, space);
1661: int cInd = (dim - tDim)-1;
1663: j += dim;
1664: for(int i = 0, k = start+tDim+offset; i < dim; ++i, --j) {
1665: if ((cInd >= 0) && (j == cDof[cInd+cOffset])) {
1666: if (!freeOnly) indices[(*indx)++] = -(offset+i+1);
1667: if (skipConstraints) --k;
1668: --cInd;
1669: } else {
1670: indices[(*indx)++] = --k;
1671: }
1672: }
1673: j += dim;
1674: offset += dim;
1675: cOffset += dim - tDim;
1676: }
1677: }
1678: }
1679: };
1680: public: // Allocation
1681: void allocateStorage() {
1682: const int totalSize = this->sizeWithBC();
1683: const value_type dummy(0) ;
1685: this->_array = this->_allocator.allocate(totalSize);
1686: ///this->_array = new value_type[totalSize];
1687: this->_sharedStorage = false;
1688: this->_sharedStorageSize = 0;
1689: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1690: ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1691: this->_bc->allocatePoint();
1692: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = this->_bcs.begin(); b_iter != this->_bcs.end(); ++b_iter) {
1693: (*b_iter)->allocatePoint();;
1694: }
1695: };
1696: void replaceStorage(value_type *newArray, const bool sharedStorage = false, const int sharedStorageSize = 0) {
1697: if (this->_array && !this->_sharedStorage) {
1698: const int totalSize = this->sizeWithBC();
1700: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1701: this->_allocator.deallocate(this->_array, totalSize);
1702: ///delete [] this->_array;
1703: }
1704: this->_array = newArray;
1705: this->_sharedStorage = sharedStorage;
1706: this->_sharedStorageSize = sharedStorageSize;
1707: this->_atlas = this->_atlasNew;
1708: this->_atlasNew = NULL;
1709: };
1710: void addPoint(const point_type& point, const int dim) {
1711: if (dim == 0) return;
1712: if (this->_atlasNew.isNull()) {
1713: this->_atlasNew = new atlas_type(this->comm(), this->debug());
1714: this->_atlasNew->copy(this->_atlas);
1715: }
1716: const index_type idx(dim, -1);
1717: this->_atlasNew->addPoint(point);
1718: this->_atlasNew->updatePoint(point, &idx);
1719: };
1720: void orderPoints(const Obj<atlas_type>& atlas){
1721: const chart_type& chart = this->getChart();
1722: int offset = 0;
1724: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1725: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
1726: const int& dim = idx.prefix;
1728: idx.index = offset;
1729: atlas->updatePoint(*c_iter, &idx);
1730: offset += dim;
1731: }
1732: };
1733: void allocatePoint() {
1734: this->orderPoints(this->_atlas);
1735: this->allocateStorage();
1736: };
1737: public: // Restriction and Update
1738: // Zero entries
1739: void zero() {
1740: this->set(0.0);
1741: };
1742: void zeroWithBC() {
1743: memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1744: };
1745: void set(const value_type value) {
1746: //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1747: const chart_type& chart = this->getChart();
1749: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1750: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1751: const int& dim = this->getFiberDimension(*c_iter);
1752: const int& cDim = this->getConstraintDimension(*c_iter);
1754: if (!cDim) {
1755: for(int i = 0; i < dim; ++i) {
1756: array[i] = value;
1757: }
1758: } else {
1759: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1760: int cInd = 0;
1762: for(int i = 0; i < dim; ++i) {
1763: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1764: array[i] = value;
1765: }
1766: }
1767: }
1768: };
1769: // Add two sections and put the result in a third
1770: void add(const Obj<GeneralSection>& x, const Obj<GeneralSection>& y) {
1771: // Check atlases
1772: const chart_type& chart = this->getChart();
1774: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1775: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1776: const value_type *xArray = x->restrictPoint(*c_iter);
1777: const value_type *yArray = y->restrictPoint(*c_iter);
1778: const int& dim = this->getFiberDimension(*c_iter);
1779: const int& cDim = this->getConstraintDimension(*c_iter);
1781: if (!cDim) {
1782: for(int i = 0; i < dim; ++i) {
1783: array[i] = xArray[i] + yArray[i];
1784: }
1785: } else {
1786: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1787: int cInd = 0;
1789: for(int i = 0; i < dim; ++i) {
1790: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1791: array[i] = xArray[i] + yArray[i];
1792: }
1793: }
1794: }
1795: };
1796: // this = this + alpha * x
1797: template<typename OtherSection>
1798: void axpy(const value_type alpha, const Obj<OtherSection>& x) {
1799: // Check atlases
1800: const chart_type& chart = this->getChart();
1802: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1803: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1804: const value_type *xArray = x->restrictPoint(*c_iter);
1805: const int& dim = this->getFiberDimension(*c_iter);
1806: const int& cDim = this->getConstraintDimension(*c_iter);
1808: if (!cDim) {
1809: for(int i = 0; i < dim; ++i) {
1810: array[i] += alpha*xArray[i];
1811: }
1812: } else {
1813: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1814: int cInd = 0;
1816: for(int i = 0; i < dim; ++i) {
1817: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1818: array[i] += alpha*xArray[i];
1819: }
1820: }
1821: }
1822: }
1823: // Return the free values on a point
1824: const value_type *restrictSpace() const {
1825: return this->_array;
1826: }
1827: // Return the free values on a point
1828: const value_type *restrictPoint(const point_type& p) const {
1829: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1830: }
1831: void restrictPoint(const point_type& p, value_type values[], const int size) const {
1832: assert(this->_atlas->restrictPoint(p)[0].prefix == size);
1833: const value_type *v = &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1835: for(int i = 0; i < size; ++i) {
1836: values[i] = v[i];
1837: }
1838: };
1839: // Update the free values on a point
1840: // Takes a full array and ignores constrained values
1841: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1842: value_type *array = (value_type *) this->restrictPoint(p);
1843: const int& cDim = this->getConstraintDimension(p);
1845: if (!cDim) {
1846: if (orientation >= 0) {
1847: const int& dim = this->getFiberDimension(p);
1849: for(int i = 0; i < dim; ++i) {
1850: array[i] = v[i];
1851: }
1852: } else {
1853: int offset = 0;
1854: int j = -1;
1856: for(int space = 0; space < this->getNumSpaces(); ++space) {
1857: const int& dim = this->getFiberDimension(p, space);
1859: for(int i = dim-1; i >= 0; --i) {
1860: array[++j] = v[i+offset];
1861: }
1862: offset += dim;
1863: }
1864: }
1865: } else {
1866: if (orientation >= 0) {
1867: const int& dim = this->getFiberDimension(p);
1868: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1869: int cInd = 0;
1871: for(int i = 0; i < dim; ++i) {
1872: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1873: array[i] = v[i];
1874: }
1875: } else {
1876: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1877: int offset = 0;
1878: int cOffset = 0;
1879: int j = 0;
1881: for(int space = 0; space < this->getNumSpaces(); ++space) {
1882: const int dim = this->getFiberDimension(p, space);
1883: const int tDim = this->getConstrainedFiberDimension(p, space);
1884: const int sDim = dim - tDim;
1885: int cInd = 0;
1887: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1888: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1889: array[j] = v[k];
1890: }
1891: offset += dim;
1892: cOffset += dim - tDim;
1893: }
1894: }
1895: }
1896: };
1897: // Update the free values on a point
1898: // Takes a full array and ignores constrained values
1899: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1900: value_type *array = (value_type *) this->restrictPoint(p);
1901: const int& cDim = this->getConstraintDimension(p);
1903: if (!cDim) {
1904: if (orientation >= 0) {
1905: const int& dim = this->getFiberDimension(p);
1907: for(int i = 0; i < dim; ++i) {
1908: array[i] += v[i];
1909: }
1910: } else {
1911: int offset = 0;
1912: int j = -1;
1914: for(int space = 0; space < this->getNumSpaces(); ++space) {
1915: const int& dim = this->getFiberDimension(p, space);
1917: for(int i = dim-1; i >= 0; --i) {
1918: array[++j] += v[i+offset];
1919: }
1920: offset += dim;
1921: }
1922: }
1923: } else {
1924: if (orientation >= 0) {
1925: const int& dim = this->getFiberDimension(p);
1926: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1927: int cInd = 0;
1929: for(int i = 0; i < dim; ++i) {
1930: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1931: array[i] += v[i];
1932: }
1933: } else {
1934: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1935: int offset = 0;
1936: int cOffset = 0;
1937: int j = 0;
1939: for(int space = 0; space < this->getNumSpaces(); ++space) {
1940: const int dim = this->getFiberDimension(p, space);
1941: const int tDim = this->getConstrainedFiberDimension(p, space);
1942: const int sDim = dim - tDim;
1943: int cInd = 0;
1945: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1946: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1947: array[j] += v[k];
1948: }
1949: offset += dim;
1950: cOffset += dim - tDim;
1951: }
1952: }
1953: }
1954: };
1955: // Update the free values on a point
1956: // Takes ONLY unconstrained values
1957: void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1958: value_type *array = (value_type *) this->restrictPoint(p);
1959: const int& cDim = this->getConstraintDimension(p);
1961: if (!cDim) {
1962: if (orientation >= 0) {
1963: const int& dim = this->getFiberDimension(p);
1965: for(int i = 0; i < dim; ++i) {
1966: array[i] = v[i];
1967: }
1968: } else {
1969: int offset = 0;
1970: int j = -1;
1972: for(int space = 0; space < this->getNumSpaces(); ++space) {
1973: const int& dim = this->getFiberDimension(p, space);
1975: for(int i = dim-1; i >= 0; --i) {
1976: array[++j] = v[i+offset];
1977: }
1978: offset += dim;
1979: }
1980: }
1981: } else {
1982: if (orientation >= 0) {
1983: const int& dim = this->getFiberDimension(p);
1984: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1985: int cInd = 0;
1987: for(int i = 0, k = -1; i < dim; ++i) {
1988: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1989: array[i] = v[++k];
1990: }
1991: } else {
1992: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1993: int offset = 0;
1994: int cOffset = 0;
1995: int j = 0;
1997: for(int space = 0; space < this->getNumSpaces(); ++space) {
1998: const int dim = this->getFiberDimension(p, space);
1999: const int tDim = this->getConstrainedFiberDimension(p, space);
2000: const int sDim = dim - tDim;
2001: int cInd = 0;
2003: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2004: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2005: array[j] = v[--k];
2006: }
2007: offset += dim;
2008: cOffset += dim - tDim;
2009: }
2010: }
2011: }
2012: };
2013: // Update the free values on a point
2014: // Takes ONLY unconstrained values
2015: void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2016: value_type *array = (value_type *) this->restrictPoint(p);
2017: const int& cDim = this->getConstraintDimension(p);
2019: if (!cDim) {
2020: if (orientation >= 0) {
2021: const int& dim = this->getFiberDimension(p);
2023: for(int i = 0; i < dim; ++i) {
2024: array[i] += v[i];
2025: }
2026: } else {
2027: int offset = 0;
2028: int j = -1;
2030: for(int space = 0; space < this->getNumSpaces(); ++space) {
2031: const int& dim = this->getFiberDimension(p, space);
2033: for(int i = dim-1; i >= 0; --i) {
2034: array[++j] += v[i+offset];
2035: }
2036: offset += dim;
2037: }
2038: }
2039: } else {
2040: if (orientation >= 0) {
2041: const int& dim = this->getFiberDimension(p);
2042: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2043: int cInd = 0;
2045: for(int i = 0, k = -1; i < dim; ++i) {
2046: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2047: array[i] += v[++k];
2048: }
2049: } else {
2050: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2051: int offset = 0;
2052: int cOffset = 0;
2053: int j = 0;
2055: for(int space = 0; space < this->getNumSpaces(); ++space) {
2056: const int dim = this->getFiberDimension(p, space);
2057: const int tDim = this->getConstrainedFiberDimension(p, space);
2058: const int sDim = dim - tDim;
2059: int cInd = 0;
2061: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2062: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2063: array[j] += v[--k];
2064: }
2065: offset += dim;
2066: cOffset += dim - tDim;
2067: }
2068: }
2069: }
2070: };
2071: // Update only the constrained dofs on a point
2072: // This takes an array with ONLY bc values
2073: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2074: value_type *array = (value_type *) this->restrictPoint(p);
2075: const int& cDim = this->getConstraintDimension(p);
2077: if (cDim) {
2078: if (orientation >= 0) {
2079: const int& dim = this->getFiberDimension(p);
2080: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2081: int cInd = 0;
2083: for(int i = 0; i < dim; ++i) {
2084: if (cInd == cDim) break;
2085: if (i == cDof[cInd]) {
2086: array[i] = v[cInd];
2087: ++cInd;
2088: }
2089: }
2090: } else {
2091: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2092: int cOffset = 0;
2093: int j = 0;
2095: for(int space = 0; space < this->getNumSpaces(); ++space) {
2096: const int dim = this->getFiberDimension(p, space);
2097: const int tDim = this->getConstrainedFiberDimension(p, space);
2098: int cInd = 0;
2100: for(int i = 0; i < dim; ++i, ++j) {
2101: if (cInd < 0) break;
2102: if (j == cDof[cInd+cOffset]) {
2103: array[j] = v[cInd+cOffset];
2104: ++cInd;
2105: }
2106: }
2107: cOffset += dim - tDim;
2108: }
2109: }
2110: }
2111: };
2112: // Update only the constrained dofs on a point
2113: // This takes an array with ALL values, not just BC
2114: void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2115: value_type *array = (value_type *) this->restrictPoint(p);
2116: const int& cDim = this->getConstraintDimension(p);
2118: if (cDim) {
2119: if (orientation >= 0) {
2120: const int& dim = this->getFiberDimension(p);
2121: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2122: int cInd = 0;
2124: for(int i = 0; i < dim; ++i) {
2125: if (cInd == cDim) break;
2126: if (i == cDof[cInd]) {
2127: array[i] = v[i];
2128: ++cInd;
2129: }
2130: }
2131: } else {
2132: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2133: int offset = 0;
2134: int cOffset = 0;
2135: int j = 0;
2137: for(int space = 0; space < this->getNumSpaces(); ++space) {
2138: const int dim = this->getFiberDimension(p, space);
2139: const int tDim = this->getConstrainedFiberDimension(p, space);
2140: int cInd = 0;
2142: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2143: if (cInd < 0) break;
2144: if (j == cDof[cInd+cOffset]) {
2145: array[j] = v[k];
2146: ++cInd;
2147: }
2148: }
2149: offset += dim;
2150: cOffset += dim - tDim;
2151: }
2152: }
2153: }
2154: };
2155: // Update all dofs on a point (free and constrained)
2156: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
2157: value_type *array = (value_type *) this->restrictPoint(p);
2159: if (orientation >= 0) {
2160: const int& dim = this->getFiberDimension(p);
2162: for(int i = 0; i < dim; ++i) {
2163: array[i] = v[i];
2164: }
2165: } else {
2166: int offset = 0;
2167: int j = -1;
2169: for(int space = 0; space < this->getNumSpaces(); ++space) {
2170: const int& dim = this->getFiberDimension(p, space);
2172: for(int i = dim-1; i >= 0; --i) {
2173: array[++j] = v[i+offset];
2174: }
2175: offset += dim;
2176: }
2177: }
2178: };
2179: // Update all dofs on a point (free and constrained)
2180: void updatePointAllAdd(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: public: // Fibrations
2204: int getNumSpaces() const {return this->_spaces.size();};
2205: const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
2206: const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
2207: void addSpace() {
2208: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
2209: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
2210: this->_spaces.push_back(space);
2211: this->_bcs.push_back(bc);
2212: };
2213: int getFiberDimension(const point_type& p, const int space) const {
2214: return this->_spaces[space]->restrictPoint(p)->prefix;
2215: };
2216: void setFiberDimension(const point_type& p, int dim, const int space) {
2217: const index_type idx(dim, -1);
2218: this->_spaces[space]->addPoint(p);
2219: this->_spaces[space]->updatePoint(p, &idx);
2220: };
2221: template<typename Sequence>
2222: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
2223: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2224: this->setFiberDimension(*p_iter, dim, space);
2225: }
2226: }
2227: int getConstraintDimension(const point_type& p, const int space) const {
2228: if (!this->_bcs[space]->hasPoint(p)) return 0;
2229: return this->_bcs[space]->getFiberDimension(p);
2230: }
2231: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2232: this->_bcs[space]->setFiberDimension(p, numConstraints);
2233: }
2234: void addConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2235: this->_bcs[space]->addFiberDimension(p, numConstraints);
2236: }
2237: int getConstrainedFiberDimension(const point_type& p, const int space) const {
2238: return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
2239: }
2240: // Return the local dofs which are constrained on a point
2241: const int *getConstraintDof(const point_type& p, const int space) const {
2242: return this->_bcs[space]->restrictPoint(p);
2243: }
2244: // Set the local dofs which are constrained on a point
2245: void setConstraintDof(const point_type& p, const int dofs[], const int space) {
2246: this->_bcs[space]->updatePoint(p, dofs);
2247: }
2248: // Return the total number of free dofs
2249: int size(const int space) const {
2250: const chart_type& points = this->getChart();
2251: int size = 0;
2253: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2254: size += this->getConstrainedFiberDimension(*p_iter, space);
2255: }
2256: return size;
2257: };
2258: template<typename OtherSection>
2259: void copyFibration(const Obj<OtherSection>& section) {
2260: const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
2261: const std::vector<Obj<bc_type> >& bcs = section->getBCs();
2263: this->_spaces.clear();
2264: for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
2265: this->_spaces.push_back(*s_iter);
2266: }
2267: this->_bcs.clear();
2268: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
2269: this->_bcs.push_back(*b_iter);
2270: }
2271: }
2272: Obj<GeneralSection> getFibration(const int space) const {
2273: Obj<GeneralSection> field = new GeneralSection(this->comm(), this->debug());
2274: // Obj<atlas_type> _atlas;
2275: // std::vector<Obj<atlas_type> > _spaces;
2276: // Obj<bc_type> _bc;
2277: // std::vector<Obj<bc_type> > _bcs;
2278: field->addSpace();
2279: const chart_type& chart = this->getChart();
2281: // Copy sizes
2282: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2283: const int fDim = this->getFiberDimension(*c_iter, space);
2284: const int cDim = this->getConstraintDimension(*c_iter, space);
2286: if (fDim) {
2287: field->setFiberDimension(*c_iter, fDim);
2288: field->setFiberDimension(*c_iter, fDim, 0);
2289: }
2290: if (cDim) {
2291: field->setConstraintDimension(*c_iter, cDim);
2292: field->setConstraintDimension(*c_iter, cDim, 0);
2293: }
2294: }
2295: field->allocateStorage();
2296: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
2297: const chart_type& newChart = field->getChart();
2299: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2300: const int cDim = field->getConstraintDimension(*c_iter);
2301: const int dof[1] = {0};
2303: if (cDim) {
2304: field->setConstraintDof(*c_iter, this->getConstraintDof(*c_iter, space));
2305: }
2306: }
2307: // Copy offsets
2308: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2309: index_type idx;
2311: idx.prefix = field->getFiberDimension(*c_iter);
2312: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
2313: for(int s = 0; s < space; ++s) {
2314: idx.index += this->getFiberDimension(*c_iter, s);
2315: }
2316: newAtlas->addPoint(*c_iter);
2317: newAtlas->updatePoint(*c_iter, &idx);
2318: }
2319: field->replaceStorage(this->_array, true, this->getStorageSize());
2320: field->setAtlas(newAtlas);
2321: return field;
2322: };
2323: public: // Optimization
2324: void getCustomRestrictAtlas(const int tag, const int *offsets[], const int *indices[]) {
2325: *offsets = this->_customRestrictAtlas[tag].first.first;
2326: *indices = this->_customRestrictAtlas[tag].first.second;
2327: };
2328: void getCustomUpdateAtlas(const int tag, const int *offsets[], const int *indices[]) {
2329: *offsets = this->_customUpdateAtlas[tag].first.first;
2330: *indices = this->_customUpdateAtlas[tag].first.second;
2331: };
2332: // This returns the tag assigned to the atlas
2333: int setCustomAtlas(const int restrictOffsets[], const int restrictIndices[], const int updateOffsets[], const int updateIndices[], bool autoFree = true) {
2334: this->_customRestrictAtlas.push_back(customAtlas_type(customAtlasInd_type(restrictOffsets, restrictIndices), autoFree));
2335: this->_customUpdateAtlas.push_back(customAtlas_type(customAtlasInd_type(updateOffsets, updateIndices), autoFree));
2336: return this->_customUpdateAtlas.size()-1;
2337: };
2338: int copyCustomAtlas(const Obj<GeneralSection>& section, const int tag) {
2339: const int *rOffsets, *rIndices, *uOffsets, *uIndices;
2341: section->getCustomRestrictAtlas(tag, &rOffsets, &rIndices);
2342: section->getCustomUpdateAtlas(tag, &uOffsets, &uIndices);
2343: return this->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices, false);
2344: };
2345: public:
2346: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
2347: ostringstream txt;
2348: int rank;
2350: if (comm == MPI_COMM_NULL) {
2351: comm = this->comm();
2352: rank = this->commRank();
2353: } else {
2354: MPI_Comm_rank(comm, &rank);
2355: }
2356: if (name == "") {
2357: if(rank == 0) {
2358: txt << "viewing a GeneralSection" << std::endl;
2359: }
2360: } else {
2361: if (rank == 0) {
2362: txt << "viewing GeneralSection '" << name << "'" << std::endl;
2363: }
2364: }
2365: if (rank == 0) {
2366: txt << " Fields: " << this->getNumSpaces() << std::endl;
2367: }
2368: const chart_type& chart = this->getChart();
2370: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2371: const value_type *array = this->restrictPoint(*p_iter);
2372: const int& dim = this->getFiberDimension(*p_iter);
2374: if (dim != 0) {
2375: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << " ";
2376: for(int i = 0; i < dim; i++) {
2377: txt << " " << array[i];
2378: }
2379: const int& dim = this->getConstraintDimension(*p_iter);
2381: if (dim) {
2382: const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);
2384: txt << " constrained";
2385: for(int i = 0; i < dim; ++i) {
2386: txt << " " << bcArray[i];
2387: }
2388: }
2389: txt << std::endl;
2390: }
2391: }
2392: if (chart.size() == 0) {
2393: txt << "[" << this->commRank() << "]: empty" << std::endl;
2394: }
2395: PetscSynchronizedPrintf(comm, txt.str().c_str());
2396: PetscSynchronizedFlush(comm);
2397: };
2398: };
2399: // FEMSection will support vector BC on a subset of unknowns on a point
2400: // We make a separate constraint Section to hold the transformation and projection operators
2401: // Storage will be contiguous by node, just as in Section
2402: // This allows fast restrict(p)
2403: // Then update() is accomplished by projecting constrained unknowns
2404: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
2405: typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
2406: typename BCAtlas_ = UniformSection<Point_, int, 1, typename Alloc_::template rebind<int>::other>,
2407: typename BC_ = Section<Point_, double, typename Alloc_::template rebind<double>::other> >
2408: class FEMSection : public ALE::ParallelObject {
2409: public:
2410: typedef Point_ point_type;
2411: typedef Value_ value_type;
2412: typedef Alloc_ alloc_type;
2413: typedef Atlas_ atlas_type;
2414: typedef BCAtlas_ bc_atlas_type;
2415: typedef BC_ bc_type;
2416: typedef Point index_type;
2417: typedef typename atlas_type::chart_type chart_type;
2418: typedef value_type * values_type;
2419: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
2420: typedef typename atlas_alloc_type::pointer atlas_ptr;
2421: typedef typename alloc_type::template rebind<bc_type>::other bc_atlas_alloc_type;
2422: typedef typename bc_atlas_alloc_type::pointer bc_atlas_ptr;
2423: typedef typename alloc_type::template rebind<bc_type>::other bc_alloc_type;
2424: typedef typename bc_alloc_type::pointer bc_ptr;
2425: protected:
2426: Obj<atlas_type> _atlas;
2427: Obj<bc_atlas_type> _bc_atlas;
2428: Obj<bc_type> _bc;
2429: values_type _array;
2430: bool _sharedStorage;
2431: int _sharedStorageSize;
2432: alloc_type _allocator;
2433: std::vector<Obj<atlas_type> > _spaces;
2434: std::vector<Obj<bc_type> > _bcs;
2435: public:
2436: FEMSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
2437: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
2438: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
2439: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
2440: bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2441: bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(comm, debug));
2442: this->_bc_atlas = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2443: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
2444: bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
2445: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
2446: this->_array = NULL;
2447: this->_sharedStorage = false;
2448: };
2449: FEMSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
2450: bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2451: bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(this->comm(), this->debug()));
2452: this->_bc_atlas = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2453: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
2454: bc_alloc_type(this->_allocator).construct(pBC, bc_type(this->comm(), this->debug()));
2455: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
2456: };
2457: ~FEMSection() {
2458: if (this->_array && !this->_sharedStorage) {
2459: const int totalSize = this->sizeWithBC();
2461: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2462: this->_allocator.deallocate(this->_array, totalSize);
2463: this->_array = NULL;
2464: }
2465: };
2466: public:
2467: value_type *getRawArray(const int size) {
2468: // Put in a sentinel value that deallocates the array
2469: static value_type *array = NULL;
2470: static int maxSize = 0;
2472: if (size > maxSize) {
2473: const value_type dummy(0);
2475: if (array) {
2476: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
2477: this->_allocator.deallocate(array, maxSize);
2478: }
2479: maxSize = size;
2480: array = this->_allocator.allocate(maxSize);
2481: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
2482: };
2483: return array;
2484: };
2485: int getStorageSize() const {
2486: if (this->_sharedStorage) {
2487: return this->_sharedStorageSize;
2488: }
2489: return this->sizeWithBC();
2490: };
2491: public: // Verifiers
2492: bool hasPoint(const point_type& point) {
2493: return this->_atlas->hasPoint(point);
2494: };
2495: public: // Accessors
2496: const chart_type& getChart() const {return this->_atlas->getChart();};
2497: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
2498: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
2499: const Obj<bc_type>& getBC() const {return this->_bc;};
2500: void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
2501: public: // BC
2502: // Returns the number of constraints on a point
2503: int getConstraintDimension(const point_type& p) const {
2504: if (!this->_bc_atlas->hasPoint(p)) return 0;
2505: return this->_bc_atlas->restrictPoint(p)[0];
2506: };
2507: // Set the number of constraints on a point
2508: void setConstraintDimension(const point_type& p, const int numConstraints) {
2509: this->_bc_atlas->updatePoint(p, &numConstraints);
2510: };
2511: // Increment the number of constraints on a point
2512: void addConstraintDimension(const point_type& p, const int numConstraints) {
2513: this->_bc_atlas->updatePointAdd(p, &numConstraints);
2514: };
2515: const int *getConstraintDof(const point_type& p) const {
2516: return this->_bc->restrictPoint(p);
2517: }
2518: public: // Sizes
2519: void clear() {
2520: if (!this->_sharedStorage) {
2521: const int totalSize = this->sizeWithBC();
2523: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2524: this->_allocator.deallocate(this->_array, totalSize);
2525: }
2526: this->_array = NULL;
2527: this->_atlas->clear();
2528: this->_bc_atlas->clear();
2529: this->_bc->clear();
2530: };
2531: // Return the total number of dofs on the point (free and constrained)
2532: int getFiberDimension(const point_type& p) const {
2533: return this->_atlas->restrictPoint(p)->prefix;
2534: };
2535: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
2536: return atlas->restrictPoint(p)->prefix;
2537: };
2538: // Set the total number of dofs on the point (free and constrained)
2539: void setFiberDimension(const point_type& p, int dim) {
2540: const index_type idx(dim, -1);
2541: this->_atlas->addPoint(p);
2542: this->_atlas->updatePoint(p, &idx);
2543: };
2544: template<typename Sequence>
2545: void setFiberDimension(const Obj<Sequence>& points, int dim) {
2546: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2547: this->setFiberDimension(*p_iter, dim);
2548: }
2549: }
2550: void addFiberDimension(const point_type& p, int dim) {
2551: if (this->_atlas->hasPoint(p)) {
2552: const index_type values(dim, 0);
2553: this->_atlas->updateAddPoint(p, &values);
2554: } else {
2555: this->setFiberDimension(p, dim);
2556: }
2557: };
2558: // Return the number of constrained dofs on this point
2559: int getConstrainedFiberDimension(const point_type& p) const {
2560: return this->getFiberDimension(p) - this->getConstraintDimension(p);
2561: };
2562: // Return the total number of free dofs
2563: int size() const {
2564: const chart_type& points = this->getChart();
2565: int size = 0;
2567: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2568: size += this->getConstrainedFiberDimension(*p_iter);
2569: }
2570: return size;
2571: };
2572: // Return the total number of dofs (free and constrained)
2573: int sizeWithBC() const {
2574: const chart_type& points = this->getChart();
2575: int size = 0;
2577: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2578: size += this->getFiberDimension(*p_iter);
2579: }
2580: return size;
2581: };
2582: public: // Allocation
2583: void allocateStorage() {
2584: const int totalSize = this->sizeWithBC();
2585: const value_type dummy(0) ;
2587: this->_array = this->_allocator.allocate(totalSize);
2588: this->_sharedStorage = false;
2589: this->_sharedStorageSize = 0;
2590: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
2591: this->_bc_atlas->allocatePoint();
2592: this->_bc->allocatePoint();
2593: };
2594: void orderPoints(const Obj<atlas_type>& atlas){
2595: const chart_type& chart = this->getChart();
2596: int offset = 0;
2598: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2599: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
2600: const int& dim = idx.prefix;
2602: idx.index = offset;
2603: atlas->updatePoint(*c_iter, &idx);
2604: offset += dim;
2605: }
2606: };
2607: void allocatePoint() {
2608: this->orderPoints(this->_atlas);
2609: this->allocateStorage();
2610: };
2611: public: // Restriction and Update
2612: // Zero entries
2613: void zero() {
2614: this->set(0.0);
2615: };
2616: void set(const value_type value) {
2617: //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
2618: const chart_type& chart = this->getChart();
2620: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2621: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2622: const int& dim = this->getFiberDimension(*c_iter);
2623: const int& cDim = this->getConstraintDimension(*c_iter);
2625: if (!cDim) {
2626: for(int i = 0; i < dim; ++i) {
2627: array[i] = value;
2628: }
2629: } else {
2630: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2631: int cInd = 0;
2633: for(int i = 0; i < dim; ++i) {
2634: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2635: array[i] = value;
2636: }
2637: }
2638: }
2639: };
2640: // Add two sections and put the result in a third
2641: void add(const Obj<FEMSection>& x, const Obj<FEMSection>& y) {
2642: // Check atlases
2643: const chart_type& chart = this->getChart();
2645: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2646: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2647: const value_type *xArray = x->restrictPoint(*c_iter);
2648: const value_type *yArray = y->restrictPoint(*c_iter);
2649: const int& dim = this->getFiberDimension(*c_iter);
2650: const int& cDim = this->getConstraintDimension(*c_iter);
2652: if (!cDim) {
2653: for(int i = 0; i < dim; ++i) {
2654: array[i] = xArray[i] + yArray[i];
2655: }
2656: } else {
2657: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2658: int cInd = 0;
2660: for(int i = 0; i < dim; ++i) {
2661: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2662: array[i] = xArray[i] + yArray[i];
2663: }
2664: }
2665: }
2666: };
2667: // this = this + alpha * x
2668: void axpy(const value_type alpha, const Obj<FEMSection>& x) {
2669: // Check atlases
2670: const chart_type& chart = this->getChart();
2672: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2673: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2674: const value_type *xArray = x->restrictPoint(*c_iter);
2675: const int& dim = this->getFiberDimension(*c_iter);
2676: const int& cDim = this->getConstraintDimension(*c_iter);
2678: if (!cDim) {
2679: for(int i = 0; i < dim; ++i) {
2680: array[i] += alpha*xArray[i];
2681: }
2682: } else {
2683: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2684: int cInd = 0;
2686: for(int i = 0; i < dim; ++i) {
2687: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2688: array[i] += alpha*xArray[i];
2689: }
2690: }
2691: }
2692: };
2693: // Return the free values on a point
2694: const value_type *restrictSpace() const {
2695: return this->_array;
2696: };
2697: // Return the free values on a point
2698: const value_type *restrictPoint(const point_type& p) const {
2699: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
2700: };
2701: // Update the free values on a point
2702: // Takes a full array and ignores constrained values
2703: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2704: value_type *array = (value_type *) this->restrictPoint(p);
2705: const int& cDim = this->getConstraintDimension(p);
2707: if (!cDim) {
2708: if (orientation >= 0) {
2709: const int& dim = this->getFiberDimension(p);
2711: for(int i = 0; i < dim; ++i) {
2712: array[i] = v[i];
2713: }
2714: } else {
2715: int offset = 0;
2716: int j = -1;
2718: for(int space = 0; space < this->getNumSpaces(); ++space) {
2719: const int& dim = this->getFiberDimension(p, space);
2721: for(int i = dim-1; i >= 0; --i) {
2722: array[++j] = v[i+offset];
2723: }
2724: offset += dim;
2725: }
2726: }
2727: } else {
2728: if (orientation >= 0) {
2729: const int& dim = this->getFiberDimension(p);
2730: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2731: int cInd = 0;
2733: for(int i = 0; i < dim; ++i) {
2734: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2735: array[i] = v[i];
2736: }
2737: } else {
2738: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2739: int offset = 0;
2740: int cOffset = 0;
2741: int j = 0;
2743: for(int space = 0; space < this->getNumSpaces(); ++space) {
2744: const int dim = this->getFiberDimension(p, space);
2745: const int tDim = this->getConstrainedFiberDimension(p, space);
2746: const int sDim = dim - tDim;
2747: int cInd = 0;
2749: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2750: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2751: array[j] = v[k];
2752: }
2753: offset += dim;
2754: cOffset += dim - tDim;
2755: }
2756: }
2757: }
2758: };
2759: // Update the free values on a point
2760: // Takes a full array and ignores constrained values
2761: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2762: value_type *array = (value_type *) this->restrictPoint(p);
2763: const int& cDim = this->getConstraintDimension(p);
2765: if (!cDim) {
2766: if (orientation >= 0) {
2767: const int& dim = this->getFiberDimension(p);
2769: for(int i = 0; i < dim; ++i) {
2770: array[i] += v[i];
2771: }
2772: } else {
2773: int offset = 0;
2774: int j = -1;
2776: for(int space = 0; space < this->getNumSpaces(); ++space) {
2777: const int& dim = this->getFiberDimension(p, space);
2779: for(int i = dim-1; i >= 0; --i) {
2780: array[++j] += v[i+offset];
2781: }
2782: offset += dim;
2783: }
2784: }
2785: } else {
2786: if (orientation >= 0) {
2787: const int& dim = this->getFiberDimension(p);
2788: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2789: int cInd = 0;
2791: for(int i = 0; i < dim; ++i) {
2792: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2793: array[i] += v[i];
2794: }
2795: } else {
2796: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2797: int offset = 0;
2798: int cOffset = 0;
2799: int j = 0;
2801: for(int space = 0; space < this->getNumSpaces(); ++space) {
2802: const int dim = this->getFiberDimension(p, space);
2803: const int tDim = this->getConstrainedFiberDimension(p, space);
2804: const int sDim = dim - tDim;
2805: int cInd = 0;
2807: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2808: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2809: array[j] += v[k];
2810: }
2811: offset += dim;
2812: cOffset += dim - tDim;
2813: }
2814: }
2815: }
2816: };
2817: // Update the free values on a point
2818: // Takes ONLY unconstrained values
2819: void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2820: value_type *array = (value_type *) this->restrictPoint(p);
2821: const int& cDim = this->getConstraintDimension(p);
2823: if (!cDim) {
2824: if (orientation >= 0) {
2825: const int& dim = this->getFiberDimension(p);
2827: for(int i = 0; i < dim; ++i) {
2828: array[i] = v[i];
2829: }
2830: } else {
2831: int offset = 0;
2832: int j = -1;
2834: for(int space = 0; space < this->getNumSpaces(); ++space) {
2835: const int& dim = this->getFiberDimension(p, space);
2837: for(int i = dim-1; i >= 0; --i) {
2838: array[++j] = v[i+offset];
2839: }
2840: offset += dim;
2841: }
2842: }
2843: } else {
2844: if (orientation >= 0) {
2845: const int& dim = this->getFiberDimension(p);
2846: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2847: int cInd = 0;
2849: for(int i = 0, k = -1; i < dim; ++i) {
2850: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2851: array[i] = v[++k];
2852: }
2853: } else {
2854: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2855: int offset = 0;
2856: int cOffset = 0;
2857: int j = 0;
2859: for(int space = 0; space < this->getNumSpaces(); ++space) {
2860: const int dim = this->getFiberDimension(p, space);
2861: const int tDim = this->getConstrainedFiberDimension(p, space);
2862: const int sDim = dim - tDim;
2863: int cInd = 0;
2865: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2866: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2867: array[j] = v[--k];
2868: }
2869: offset += dim;
2870: cOffset += dim - tDim;
2871: }
2872: }
2873: }
2874: };
2875: // Update the free values on a point
2876: // Takes ONLY unconstrained values
2877: void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2878: value_type *array = (value_type *) this->restrictPoint(p);
2879: const int& cDim = this->getConstraintDimension(p);
2881: if (!cDim) {
2882: if (orientation >= 0) {
2883: const int& dim = this->getFiberDimension(p);
2885: for(int i = 0; i < dim; ++i) {
2886: array[i] += v[i];
2887: }
2888: } else {
2889: int offset = 0;
2890: int j = -1;
2892: for(int space = 0; space < this->getNumSpaces(); ++space) {
2893: const int& dim = this->getFiberDimension(p, space);
2895: for(int i = dim-1; i >= 0; --i) {
2896: array[++j] += v[i+offset];
2897: }
2898: offset += dim;
2899: }
2900: }
2901: } else {
2902: if (orientation >= 0) {
2903: const int& dim = this->getFiberDimension(p);
2904: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2905: int cInd = 0;
2907: for(int i = 0, k = -1; i < dim; ++i) {
2908: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2909: array[i] += v[++k];
2910: }
2911: } else {
2912: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2913: int offset = 0;
2914: int cOffset = 0;
2915: int j = 0;
2917: for(int space = 0; space < this->getNumSpaces(); ++space) {
2918: const int dim = this->getFiberDimension(p, space);
2919: const int tDim = this->getConstrainedFiberDimension(p, space);
2920: const int sDim = dim - tDim;
2921: int cInd = 0;
2923: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2924: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2925: array[j] += v[--k];
2926: }
2927: offset += dim;
2928: cOffset += dim - tDim;
2929: }
2930: }
2931: }
2932: };
2933: // Update only the constrained dofs on a point
2934: // This takes an array with ONLY bc values
2935: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2936: value_type *array = (value_type *) this->restrictPoint(p);
2937: const int& cDim = this->getConstraintDimension(p);
2939: if (cDim) {
2940: if (orientation >= 0) {
2941: const int& dim = this->getFiberDimension(p);
2942: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2943: int cInd = 0;
2945: for(int i = 0; i < dim; ++i) {
2946: if (cInd == cDim) break;
2947: if (i == cDof[cInd]) {
2948: array[i] = v[cInd];
2949: ++cInd;
2950: }
2951: }
2952: } else {
2953: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2954: int cOffset = 0;
2955: int j = 0;
2957: for(int space = 0; space < this->getNumSpaces(); ++space) {
2958: const int dim = this->getFiberDimension(p, space);
2959: const int tDim = this->getConstrainedFiberDimension(p, space);
2960: int cInd = 0;
2962: for(int i = 0; i < dim; ++i, ++j) {
2963: if (cInd < 0) break;
2964: if (j == cDof[cInd+cOffset]) {
2965: array[j] = v[cInd+cOffset];
2966: ++cInd;
2967: }
2968: }
2969: cOffset += dim - tDim;
2970: }
2971: }
2972: }
2973: };
2974: // Update only the constrained dofs on a point
2975: // This takes an array with ALL values, not just BC
2976: void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2977: value_type *array = (value_type *) this->restrictPoint(p);
2978: const int& cDim = this->getConstraintDimension(p);
2980: if (cDim) {
2981: if (orientation >= 0) {
2982: const int& dim = this->getFiberDimension(p);
2983: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2984: int cInd = 0;
2986: for(int i = 0; i < dim; ++i) {
2987: if (cInd == cDim) break;
2988: if (i == cDof[cInd]) {
2989: array[i] = v[i];
2990: ++cInd;
2991: }
2992: }
2993: } else {
2994: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2995: int offset = 0;
2996: int cOffset = 0;
2997: int j = 0;
2999: for(int space = 0; space < this->getNumSpaces(); ++space) {
3000: const int dim = this->getFiberDimension(p, space);
3001: const int tDim = this->getConstrainedFiberDimension(p, space);
3002: int cInd = 0;
3004: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
3005: if (cInd < 0) break;
3006: if (j == cDof[cInd+cOffset]) {
3007: array[j] = v[k];
3008: ++cInd;
3009: }
3010: }
3011: offset += dim;
3012: cOffset += dim - tDim;
3013: }
3014: }
3015: }
3016: };
3017: // Update all dofs on a point (free and constrained)
3018: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
3019: value_type *array = (value_type *) this->restrictPoint(p);
3021: if (orientation >= 0) {
3022: const int& dim = this->getFiberDimension(p);
3024: for(int i = 0; i < dim; ++i) {
3025: array[i] = v[i];
3026: }
3027: } else {
3028: int offset = 0;
3029: int j = -1;
3031: for(int space = 0; space < this->getNumSpaces(); ++space) {
3032: const int& dim = this->getFiberDimension(p, space);
3034: for(int i = dim-1; i >= 0; --i) {
3035: array[++j] = v[i+offset];
3036: }
3037: offset += dim;
3038: }
3039: }
3040: };
3041: // Update all dofs on a point (free and constrained)
3042: void updatePointAllAdd(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: public: // Fibrations
3066: int getNumSpaces() const {return this->_spaces.size();};
3067: const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
3068: const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
3069: void addSpace() {
3070: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
3071: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
3072: this->_spaces.push_back(space);
3073: this->_bcs.push_back(bc);
3074: };
3075: int getFiberDimension(const point_type& p, const int space) const {
3076: return this->_spaces[space]->restrictPoint(p)->prefix;
3077: };
3078: void setFiberDimension(const point_type& p, int dim, const int space) {
3079: const index_type idx(dim, -1);
3080: this->_spaces[space]->addPoint(p);
3081: this->_spaces[space]->updatePoint(p, &idx);
3082: };
3083: template<typename Sequence>
3084: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
3085: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
3086: this->setFiberDimension(*p_iter, dim, space);
3087: }
3088: }
3089: int getConstraintDimension(const point_type& p, const int space) const {
3090: if (!this->_bcs[space]->hasPoint(p)) return 0;
3091: return this->_bcs[space]->getFiberDimension(p);
3092: };
3093: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
3094: this->_bcs[space]->setFiberDimension(p, numConstraints);
3095: };
3096: int getConstrainedFiberDimension(const point_type& p, const int space) const {
3097: return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
3098: };
3099: void copyFibration(const Obj<FEMSection>& section) {
3100: const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
3101: const std::vector<Obj<bc_type> >& bcs = section->getBCs();
3103: this->_spaces.clear();
3104: for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
3105: this->_spaces.push_back(*s_iter);
3106: }
3107: this->_bcs.clear();
3108: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
3109: this->_bcs.push_back(*b_iter);
3110: }
3111: };
3112: Obj<FEMSection> getFibration(const int space) const {
3113: Obj<FEMSection> field = new FEMSection(this->comm(), this->debug());
3114: // Obj<atlas_type> _atlas;
3115: // std::vector<Obj<atlas_type> > _spaces;
3116: // Obj<bc_type> _bc;
3117: // std::vector<Obj<bc_type> > _bcs;
3118: field->addSpace();
3119: const chart_type& chart = this->getChart();
3121: // Copy sizes
3122: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3123: const int fDim = this->getFiberDimension(*c_iter, space);
3124: const int cDim = this->getConstraintDimension(*c_iter, space);
3126: if (fDim) {
3127: field->setFiberDimension(*c_iter, fDim);
3128: field->setFiberDimension(*c_iter, fDim, 0);
3129: }
3130: if (cDim) {
3131: field->setConstraintDimension(*c_iter, cDim);
3132: field->setConstraintDimension(*c_iter, cDim, 0);
3133: }
3134: }
3135: field->allocateStorage();
3136: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
3137: const chart_type& newChart = field->getChart();
3139: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3140: const int cDim = field->getConstraintDimension(*c_iter);
3141: const int dof[1] = {0};
3143: if (cDim) {
3144: field->setConstraintDof(*c_iter, dof);
3145: }
3146: }
3147: // Copy offsets
3148: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3149: index_type idx;
3151: idx.prefix = field->getFiberDimension(*c_iter);
3152: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
3153: for(int s = 0; s < space; ++s) {
3154: idx.index += this->getFiberDimension(*c_iter, s);
3155: }
3156: newAtlas->addPoint(*c_iter);
3157: newAtlas->updatePoint(*c_iter, &idx);
3158: }
3159: field->replaceStorage(this->_array, true, this->getStorageSize());
3160: field->setAtlas(newAtlas);
3161: return field;
3162: };
3163: public:
3164: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3165: ostringstream txt;
3166: int rank;
3168: if (comm == MPI_COMM_NULL) {
3169: comm = this->comm();
3170: rank = this->commRank();
3171: } else {
3172: MPI_Comm_rank(comm, &rank);
3173: }
3174: if (name == "") {
3175: if(rank == 0) {
3176: txt << "viewing a FEMSection" << std::endl;
3177: }
3178: } else {
3179: if (rank == 0) {
3180: txt << "viewing FEMSection '" << name << "'" << std::endl;
3181: }
3182: }
3183: if (rank == 0) {
3184: txt << " Fields: " << this->getNumSpaces() << std::endl;
3185: }
3186: const chart_type& chart = this->getChart();
3188: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
3189: const value_type *array = this->restrictPoint(*p_iter);
3190: const int& dim = this->getFiberDimension(*p_iter);
3192: if (dim != 0) {
3193: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << " ";
3194: for(int i = 0; i < dim; i++) {
3195: txt << " " << array[i];
3196: }
3197: const int& dim = this->getConstraintDimension(*p_iter);
3199: if (dim) {
3200: const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);
3202: txt << " constrained";
3203: for(int i = 0; i < dim; ++i) {
3204: txt << " " << bcArray[i];
3205: }
3206: }
3207: txt << std::endl;
3208: }
3209: }
3210: if (chart.size() == 0) {
3211: txt << "[" << this->commRank() << "]: empty" << std::endl;
3212: }
3213: PetscSynchronizedPrintf(comm, txt.str().c_str());
3214: PetscSynchronizedFlush(comm);
3215: };
3216: };
3217: // A Field combines several sections
3218: template<typename Overlap_, typename Patch_, typename Section_>
3219: class Field : public ALE::ParallelObject {
3220: public:
3221: typedef Overlap_ overlap_type;
3222: typedef Patch_ patch_type;
3223: typedef Section_ section_type;
3224: typedef typename section_type::point_type point_type;
3225: typedef typename section_type::chart_type chart_type;
3226: typedef typename section_type::value_type value_type;
3227: typedef std::map<patch_type, Obj<section_type> > sheaf_type;
3228: typedef enum {SEND, RECEIVE} request_type;
3229: typedef std::map<patch_type, MPI_Request> requests_type;
3230: protected:
3231: sheaf_type _sheaf;
3232: int _tag;
3233: MPI_Datatype _datatype;
3234: requests_type _requests;
3235: public:
3236: Field(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
3237: this->_tag = this->getNewTag();
3238: this->_datatype = this->getMPIDatatype();
3239: };
3240: Field(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _tag(tag) {
3241: this->_datatype = this->getMPIDatatype();
3242: };
3243: virtual ~Field() {};
3244: protected:
3245: MPI_Datatype getMPIDatatype() {
3246: if (sizeof(value_type) == 4) {
3247: return MPI_INT;
3248: } else if (sizeof(value_type) == 8) {
3249: return MPI_DOUBLE;
3250: } else if (sizeof(value_type) == 28) {
3251: int blen[2];
3252: MPI_Aint indices[2];
3253: MPI_Datatype oldtypes[2], newtype;
3254: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
3255: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3256: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3257: MPI_Type_commit(&newtype);
3258: return newtype;
3259: } else if (sizeof(value_type) == 32) {
3260: int blen[2];
3261: MPI_Aint indices[2];
3262: MPI_Datatype oldtypes[2], newtype;
3263: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_DOUBLE;
3264: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3265: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3266: MPI_Type_commit(&newtype);
3267: return newtype;
3268: }
3269: throw ALE::Exception("Cannot determine MPI type for value type");
3270: };
3271: int getNewTag() {
3272: static int tagKeyval = MPI_KEYVAL_INVALID;
3273: int *tagvalp = NULL, *maxval, flg;
3275: if (tagKeyval == MPI_KEYVAL_INVALID) {
3276: tagvalp = (int *) malloc(sizeof(int));
3277: MPI_Keyval_create(MPI_NULL_COPY_FN, Mesh_DelTag, &tagKeyval, (void *) NULL);
3278: MPI_Attr_put(this->_comm, tagKeyval, tagvalp);
3279: tagvalp[0] = 0;
3280: }
3281: MPI_Attr_get(this->_comm, tagKeyval, (void **) &tagvalp, &flg);
3282: if (tagvalp[0] < 1) {
3283: MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
3284: tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
3285: }
3286: if (this->debug()) {
3287: std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
3288: }
3289: return tagvalp[0]--;
3290: };
3291: public: // Verifiers
3292: void checkPatch(const patch_type& patch) const {
3293: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3294: ostringstream msg;
3295: msg << "Invalid field patch " << patch << std::endl;
3296: throw ALE::Exception(msg.str().c_str());
3297: }
3298: };
3299: bool hasPatch(const patch_type& patch) {
3300: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3301: return false;
3302: }
3303: return true;
3304: };
3305: public: // Accessors
3306: int getTag() const {return this->_tag;};
3307: void setTag(const int tag) {this->_tag = tag;};
3308: Obj<section_type>& getSection(const patch_type& patch) {
3309: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3310: this->_sheaf[patch] = new section_type(this->comm(), this->debug());
3311: }
3312: return this->_sheaf[patch];
3313: };
3314: void setSection(const patch_type& patch, const Obj<section_type>& section) {this->_sheaf[patch] = section;};
3315: const sheaf_type& getPatches() {
3316: return this->_sheaf;
3317: };
3318: void clear() {
3319: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3320: p_iter->second->clear();
3321: }
3322: };
3323: public: // Adapter
3324: template<typename Topology_>
3325: void setTopology(const Obj<Topology_>& topology) {
3326: const typename Topology_::sheaf_type& patches = topology->getPatches();
3328: for(typename Topology_::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3329: int rank = p_iter->first;
3330: const Obj<section_type>& section = this->getSection(rank);
3331: const Obj<typename Topology_::sieve_type::baseSequence>& base = p_iter->second->base();
3333: for(typename Topology_::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
3334: section->setFiberDimension(*b_iter, 1);
3335: }
3336: }
3337: }
3338: void allocate() {
3339: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3340: p_iter->second->allocatePoint();
3341: }
3342: }
3343: public: // Communication
3344: void construct(const int size) {
3345: const sheaf_type& patches = this->getPatches();
3347: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3348: const patch_type rank = p_iter->first;
3349: const Obj<section_type>& section = this->getSection(rank);
3350: const chart_type& chart = section->getChart();
3352: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3353: section->setFiberDimension(*c_iter, size);
3354: }
3355: }
3356: };
3357: template<typename Sizer>
3358: void construct(const Obj<Sizer>& sizer) {
3359: const sheaf_type& patches = this->getPatches();
3361: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3362: const patch_type rank = p_iter->first;
3363: const Obj<section_type>& section = this->getSection(rank);
3364: const chart_type& chart = section->getChart();
3365:
3366: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3367: section->setFiberDimension(*c_iter, *(sizer->getSection(rank)->restrictPoint(*c_iter)));
3368: }
3369: }
3370: }
3371: void constructCommunication(const request_type& requestType) {
3372: const sheaf_type& patches = this->getPatches();
3374: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3375: const patch_type patch = p_iter->first;
3376: const Obj<section_type>& section = this->getSection(patch);
3377: MPI_Request request;
3379: if (requestType == RECEIVE) {
3380: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data(" << section->size() << ") from " << patch << " tag " << this->_tag << std::endl;}
3381: MPI_Recv_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3382: } else {
3383: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << section->size() << ") to " << patch << " tag " << this->_tag << std::endl;}
3384: MPI_Send_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3385: }
3386: this->_requests[patch] = request;
3387: }
3388: };
3389: void startCommunication() {
3390: const sheaf_type& patches = this->getPatches();
3392: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3393: MPI_Request request = this->_requests[p_iter->first];
3395: MPI_Start(&request);
3396: }
3397: };
3398: void endCommunication() {
3399: const sheaf_type& patches = this->getPatches();
3400: MPI_Status status;
3402: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3403: MPI_Request request = this->_requests[p_iter->first];
3405: MPI_Wait(&request, &status);
3406: }
3407: };
3408: public:
3409: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3410: ostringstream txt;
3411: int rank;
3413: if (comm == MPI_COMM_NULL) {
3414: comm = this->comm();
3415: rank = this->commRank();
3416: } else {
3417: MPI_Comm_rank(comm, &rank);
3418: }
3419: if (name == "") {
3420: if(rank == 0) {
3421: txt << "viewing a Field" << std::endl;
3422: }
3423: } else {
3424: if(rank == 0) {
3425: txt << "viewing Field '" << name << "'" << std::endl;
3426: }
3427: }
3428: PetscSynchronizedPrintf(comm, txt.str().c_str());
3429: PetscSynchronizedFlush(comm);
3430: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3431: ostringstream txt1;
3433: txt1 << "[" << this->commRank() << "]: Patch " << p_iter->first << std::endl;
3434: PetscSynchronizedPrintf(comm, txt1.str().c_str());
3435: PetscSynchronizedFlush(comm);
3436: p_iter->second->view("field section", comm);
3437: }
3438: };
3439: };
3440: }
3441: #endif