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