Actual source code: IField.hh
1: #ifndef included_ALE_IField_hh
2: #define included_ALE_IField_hh
4: #ifndef included_ALE_Field_hh
5: #include <Field.hh>
6: #endif
8: #ifndef included_ALE_ISieve_hh
9: #include <ISieve.hh>
10: #endif
12: // An ISection (or IntervalSection) is a section over a simple interval domain
13: namespace ALE {
14: // An IConstantSection is the simplest ISection
15: // All fibers are dimension 1
16: // All values are equal to a constant
17: // We need no value storage and no communication for completion
18: // The default value is returned for any point not in the domain
19: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
20: class IConstantSection : public ALE::ParallelObject {
21: public:
22: typedef Point_ point_type;
23: typedef Value_ value_type;
24: typedef Alloc_ alloc_type;
25: typedef Interval<point_type, alloc_type> chart_type;
26: protected:
27: chart_type _chart;
28: value_type _value[2]; // Value and default value
29: public:
30: IConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
31: _value[1] = 0;
32: };
33: IConstantSection(MPI_Comm comm, const point_type& min, const point_type& max, const value_type& value, const int debug) : ParallelObject(comm, debug), _chart(min, max) {
34: _value[0] = value;
35: _value[1] = value;
36: };
37: IConstantSection(MPI_Comm comm, const point_type& min, const point_type& max, const value_type& value, const value_type& defaultValue, const int debug) : ParallelObject(comm, debug), _chart(min, max) {
38: _value[0] = value;
39: _value[1] = defaultValue;
40: };
41: public: // Verifiers
42: void checkPoint(const point_type& point) const {
43: this->_chart.checkPoint(point);
44: };
45: void checkDimension(const int& dim) {
46: if (dim != 1) {
47: ostringstream msg;
48: msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
49: throw ALE::Exception(msg.str().c_str());
50: }
51: };
52: bool hasPoint(const point_type& point) const {
53: return this->_chart.hasPoint(point);
54: };
55: public: // Accessors
56: const chart_type& getChart() const {return this->_chart;};
57: void setChart(const chart_type& chart) {this->_chart = chart;};
58: void addPoint(const point_type& point) {
59: this->checkPoint(point);
60: };
61: template<typename Points>
62: void addPoint(const Obj<Points>& points) {
63: for(typename Points::const_iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
64: this->checkPoint(*p_iter);
65: }
66: };
67: template<typename Points>
68: void addPoint(const Points& points) {
69: for(typename Points::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
70: this->checkPoint(*p_iter);
71: }
72: };
73: value_type getDefaultValue() {return this->_value[1];};
74: void setDefaultValue(const value_type value) {this->_value[1] = value;};
75: void copy(const Obj<IConstantSection>& section) {
76: const chart_type& chart = section->getChart();
78: this->_chart = chart;
79: this->_value[0] = section->restrictPoint(*chart.begin())[0];
80: this->_value[1] = section->restrictPoint(*chart.begin())[1];
81: };
82: public: // Sizes
83: ///void clear() {};
84: int getFiberDimension(const point_type& p) const {
85: if (this->hasPoint(p)) return 1;
86: return 0;
87: };
88: void setFiberDimension(const point_type& p, int dim) {
89: this->checkDimension(dim);
90: this->addPoint(p);
91: };
92: template<typename Sequence>
93: void setFiberDimension(const Obj<Sequence>& points, int dim) {
94: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
95: this->setFiberDimension(*p_iter, dim);
96: }
97: };
98: void addFiberDimension(const point_type& p, int dim) {
99: if (this->hasPoint(p)) {
100: ostringstream msg;
101: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
102: throw ALE::Exception(msg.str().c_str());
103: } else {
104: this->setFiberDimension(p, dim);
105: }
106: };
107: int size(const point_type& p) {return this->getFiberDimension(p);};
108: public: // Restriction
109: void clear() {};
110: const value_type *restrictSpace() const {
111: return this->_value;
112: };
113: const value_type *restrictPoint(const point_type& p) const {
114: if (this->hasPoint(p)) {
115: return this->_value;
116: }
117: return &this->_value[1];
118: };
119: void updatePoint(const point_type& p, const value_type v[]) {
120: this->_value[0] = v[0];
121: };
122: void updateAddPoint(const point_type& p, const value_type v[]) {
123: this->_value[0] += v[0];
124: };
125: public:
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 = this->comm();
132: rank = this->commRank();
133: } else {
134: MPI_Comm_rank(comm, &rank);
135: }
136: if (name == "") {
137: if(rank == 0) {
138: txt << "viewing an IConstantSection" << std::endl;
139: }
140: } else {
141: if(rank == 0) {
142: txt << "viewing IConstantSection '" << name << "'" << std::endl;
143: }
144: }
145: txt <<"["<<this->commRank()<<"]: chart " << this->_chart << std::endl;
146: txt <<"["<<this->commRank()<<"]: Value " << this->_value[0] << " Default Value " << this->_value[1] << std::endl;
147: PetscSynchronizedPrintf(comm, txt.str().c_str());
148: PetscSynchronizedFlush(comm);
149: };
150: };
152: // An IUniformSection often acts as an Atlas
153: // All fibers are the same dimension
154: // Note we can use a IConstantSection for this Atlas
155: // Each point may have a different vector
156: // Thus we need storage for values, and hence must implement completion
157: template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
158: class IUniformSection : public ALE::ParallelObject {
159: public:
160: typedef Point_ point_type;
161: typedef Value_ value_type;
162: typedef Alloc_ alloc_type;
163: typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
164: typedef IConstantSection<point_type, int, point_alloc_type> atlas_type;
165: typedef typename atlas_type::chart_type chart_type;
166: typedef point_type index_type;
167: typedef struct {value_type v[fiberDim];} fiber_type;
168: typedef value_type* values_type;
169: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
170: typedef typename atlas_alloc_type::pointer atlas_ptr;
171: protected:
172: Obj<atlas_type> _atlas;
173: values_type _array;
174: fiber_type _emptyValue;
175: alloc_type _allocator;
176: public:
177: IUniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
178: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
179: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
180: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
181: this->_array = NULL;
182: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
183: };
184: IUniformSection(MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : ParallelObject(comm, debug) {
185: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
186: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, min, max, fiberDim, debug));
187: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
188: this->_array = NULL;
189: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
190: };
191: IUniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas) {
192: int dim = fiberDim;
193: this->_atlas->update(*this->_atlas->getChart().begin(), &dim);
194: this->_array = NULL;
195: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
196: };
197: virtual ~IUniformSection() {
198: if (this->_array) {
199: for(int i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.destroy(this->_array+i);}
200: this->_array += this->getChart().min()*fiberDim;
201: this->_allocator.deallocate(this->_array, this->sizeWithBC());
202: this->_array = NULL;
203: this->_atlas = NULL;
204: }
205: };
206: public:
207: value_type *getRawArray(const int size) {
208: static value_type *array = NULL;
209: static int maxSize = 0;
211: if (size > maxSize) {
212: const value_type dummy(0);
214: if (array) {
215: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
216: this->_allocator.deallocate(array, maxSize);
217: }
218: maxSize = size;
219: array = this->_allocator.allocate(maxSize);
220: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
221: };
222: return array;
223: };
224: public: // Verifiers
225: bool hasPoint(const point_type& point) const {
226: return this->_atlas->hasPoint(point);
227: };
228: void checkDimension(const int& dim) {
229: if (dim != fiberDim) {
230: ostringstream msg;
231: msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
232: throw ALE::Exception(msg.str().c_str());
233: }
234: };
235: public: // Accessors
236: const chart_type& getChart() const {return this->_atlas->getChart();};
237: void setChart(const chart_type& chart) {
238: this->_atlas->setChart(chart);
239: int dim = fiberDim;
240: this->_atlas->updatePoint(*this->getChart().begin(), &dim);
241: };
242: bool resizeChart(const chart_type& chart) {
243: if ((chart.min() >= this->getChart().min()) && (chart.max() <= this->getChart().max())) return false;
244: this->setChart(chart);
245: return true;
246: };
247: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
248: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
249: void addPoint(const point_type& point) {
250: this->setFiberDimension(point, fiberDim);
251: };
252: template<typename Points>
253: void addPoint(const Obj<Points>& points) {
254: for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
255: this->setFiberDimension(*p_iter, fiberDim);
256: }
257: };
258: void copy(const Obj<IUniformSection>& section) {
259: this->getAtlas()->copy(section->getAtlas());
260: const chart_type& chart = section->getChart();
262: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
263: this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
264: }
265: };
266: const value_type *getDefault() const {return this->_emptyValue;};
267: void setDefault(const value_type v[]) {for(int i = 0; i < fiberDim; ++i) {this->_emptyValue.v[i] = v[i];}};
268: public: // Sizes
269: void clear() {
270: this->_atlas->clear();
271: };
272: int getFiberDimension(const point_type& p) const {
273: return this->_atlas->restrictPoint(p)[0];
274: };
275: void setFiberDimension(const point_type& p, int dim) {
276: this->checkDimension(dim);
277: this->_atlas->addPoint(p);
278: this->_atlas->updatePoint(p, &dim);
279: };
280: template<typename Sequence>
281: void setFiberDimension(const Obj<Sequence>& points, int dim) {
282: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
283: this->setFiberDimension(*p_iter, dim);
284: }
285: };
286: void setFiberDimension(const std::set<point_type>& points, int dim) {
287: for(typename std::set<point_type>::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
288: this->setFiberDimension(*p_iter, dim);
289: }
290: };
291: void addFiberDimension(const point_type& p, int dim) {
292: if (this->hasPoint(p)) {
293: ostringstream msg;
294: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
295: throw ALE::Exception(msg.str().c_str());
296: } else {
297: this->setFiberDimension(p, dim);
298: }
299: };
300: int size() const {return this->_atlas->getChart().size()*fiberDim;};
301: int sizeWithBC() const {return this->size();};
302: void allocatePoint() {
303: this->_array = this->_allocator.allocate(this->sizeWithBC());
304: this->_array -= this->getChart().min()*fiberDim;
305: for(index_type i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);}
306: };
307: bool reallocatePoint(const chart_type& chart, values_type *oldData = NULL) {
308: const chart_type oldChart = this->getChart();
309: const int oldSize = this->sizeWithBC();
310: values_type oldArray = this->_array;
311: if (!this->resizeChart(chart)) return false;
312: const int size = this->sizeWithBC();
314: this->_array = this->_allocator.allocate(size);
315: this->_array -= this->getChart().min()*fiberDim;
316: for(index_type i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);}
317: for(int i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {
318: this->_array[i] = oldArray[i];
319: }
320: if (!oldData) {
321: for(index_type i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {this->_allocator.destroy(oldArray+i);}
322: oldArray += this->getChart().min()*fiberDim;
323: this->_allocator.deallocate(oldArray, oldSize);
324: ///std::cout << "Freed IUniformSection data" << std::endl;
325: } else {
326: ///std::cout << "Did not free IUniformSection data" << std::endl;
327: *oldData = oldArray;
328: }
329: return true;
330: };
331: template<typename Iterator, typename Extractor>
332: bool reallocatePoint(const Iterator& begin, const Iterator& end, const Extractor& extractor) {
333: point_type min = this->getChart().min();
334: point_type max = this->getChart().max()-1;
336: for(Iterator p_iter = begin; p_iter != end; ++p_iter) {
337: min = std::min(extractor(*p_iter), min);
338: max = std::max(extractor(*p_iter), max);
339: }
340: return reallocatePoint(chart_type(min, max+1));
341: };
342: public: // Restriction
343: // Return a pointer to the entire contiguous storage array
344: const values_type& restrictSpace() const {
345: return this->_array;
346: };
347: // Return only the values associated to this point, not its closure
348: const value_type *restrictPoint(const point_type& p) const {
349: if (!this->hasPoint(p)) return this->_emptyValue.v;
350: return &this->_array[p*fiberDim];
351: };
352: // Update only the values associated to this point, not its closure
353: void updatePoint(const point_type& p, const value_type v[]) {
354: for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
355: this->_array[idx] = v[i];
356: }
357: };
358: // Update only the values associated to this point, not its closure
359: void updateAddPoint(const point_type& p, const value_type v[]) {
360: for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
361: this->_array[idx] += v[i];
362: }
363: };
364: void updatePointAll(const point_type& p, const value_type v[]) {
365: this->updatePoint(p, v);
366: };
367: public:
368: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
369: ostringstream txt;
370: int rank;
372: if (comm == MPI_COMM_NULL) {
373: comm = this->comm();
374: rank = this->commRank();
375: } else {
376: MPI_Comm_rank(comm, &rank);
377: }
378: if (name == "") {
379: if(rank == 0) {
380: txt << "viewing an IUniformSection" << std::endl;
381: }
382: } else {
383: if(rank == 0) {
384: txt << "viewing IUniformSection '" << name << "'" << std::endl;
385: }
386: }
387: const typename atlas_type::chart_type& chart = this->_atlas->getChart();
388: values_type array = this->_array;
390: for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
391: const int idx = (*p_iter)*fiberDim;
393: if (fiberDim != 0) {
394: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << fiberDim << " ";
395: for(int i = 0; i < fiberDim; i++) {
396: txt << " " << array[idx+i];
397: }
398: txt << std::endl;
399: }
400: }
401: if (chart.size() == 0) {
402: txt << "[" << this->commRank() << "]: empty" << std::endl;
403: }
404: PetscSynchronizedPrintf(comm, txt.str().c_str());
405: PetscSynchronizedFlush(comm);
406: };
407: };
408: // An ISection allows variable fiber sizes per point
409: // The Atlas is a UniformSection of dimension 1 and value type Point
410: // to hold each fiber dimension and offsets into a contiguous array
411: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_> >
412: class ISection : public Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > {
413: public:
414: typedef Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > base;
415: typedef typename base::point_type point_type;
416: typedef typename base::value_type value_type;
417: typedef typename base::alloc_type alloc_type;
418: typedef typename base::index_type index_type;
419: typedef typename base::atlas_type atlas_type;
420: typedef typename base::chart_type chart_type;
421: typedef typename base::values_type values_type;
422: typedef typename base::atlas_alloc_type atlas_alloc_type;
423: typedef typename base::atlas_ptr atlas_ptr;
424: public:
425: ISection(MPI_Comm comm, const int debug = 0) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(comm, debug) {
426: };
427: ISection(MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(comm, debug) {
428: this->_atlas->setChart(chart_type(min, max));
429: this->_atlas->allocatePoint();
430: };
431: ISection(const Obj<atlas_type>& atlas) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(atlas) {};
432: virtual ~ISection() {};
433: public:
434: void setChart(const chart_type& chart) {
435: this->_atlas->setChart(chart);
436: this->_atlas->allocatePoint();
437: };
438: bool resizeChart(const chart_type& chart) {
439: if (!this->_atlas->reallocatePoint(chart)) return false;
440: return true;
441: };
442: bool reallocatePoint(const chart_type& chart) {
443: typedef typename atlas_type::alloc_type atlas_alloc_type;
444: const chart_type oldChart = this->getChart();
445: const int oldSize = this->sizeWithBC();
446: const values_type oldArray = this->_array;
447: const int oldAtlasSize = this->_atlas->sizeWithBC();
448: typename atlas_type::values_type oldAtlasArray;
449: if (!this->_atlas->reallocatePoint(chart, &oldAtlasArray)) return false;
451: this->orderPoints(this->_atlas);
452: this->allocateStorage();
453: for(int i = oldChart.min(); i < oldChart.max(); ++i) {
454: const typename atlas_type::value_type& idx = this->_atlas->restrictPoint(i)[0];
455: const int dim = idx.prefix;
456: const int off = idx.index;
458: for(int d = 0; d < dim; ++d) {
459: this->_array[off+d] = oldArray[oldAtlasArray[i].index+d];
460: }
461: }
462: for(int i = 0; i < oldSize; ++i) {this->_allocator.destroy(oldArray+i);}
463: this->_allocator.deallocate(oldArray, oldSize);
464: for(int i = oldChart.min(); i < oldChart.max(); ++i) {atlas_alloc_type(this->_allocator).destroy(oldAtlasArray+i);}
465: oldAtlasArray += oldChart.min();
466: atlas_alloc_type(this->_allocator).deallocate(oldAtlasArray, oldAtlasSize);
467: ///std::cout << "In ISection, Freed IUniformSection data" << std::endl;
468: };
469: public:
470: // Return the free values on a point
471: // This is overridden, because the one in Section cannot be const due to problem in the interface with UniformSection
472: const value_type *restrictPoint(const point_type& p) const {
473: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
474: };
475: };
476: // IGeneralSection will support BC on a subset of unknowns on a point
477: // We use a separate constraint Atlas to mark constrained dofs on a point
478: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_> >
479: class IGeneralSection : public GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> > {
480: public:
481: typedef GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> > base;
482: typedef typename base::point_type point_type;
483: typedef typename base::value_type value_type;
484: typedef typename base::alloc_type alloc_type;
485: typedef typename base::index_type index_type;
486: typedef typename base::atlas_type atlas_type;
487: typedef typename base::bc_type bc_type;
488: typedef typename base::chart_type chart_type;
489: typedef typename base::values_type values_type;
490: typedef typename base::atlas_alloc_type atlas_alloc_type;
491: typedef typename base::atlas_ptr atlas_ptr;
492: typedef typename base::bc_alloc_type bc_alloc_type;
493: typedef typename base::bc_ptr bc_ptr;
494: typedef std::pair<point_type,int> newpoint_type;
495: protected:
496: std::set<newpoint_type> newPoints;
497: public:
498: IGeneralSection(MPI_Comm comm, const int debug = 0) : GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> >(comm, debug) {};
499: IGeneralSection(MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> >(comm, debug) {
500: this->_atlas->setChart(chart_type(min, max));
501: this->_atlas->allocatePoint();
502: this->_bc->setChart(chart_type(min, max));
503: };
504: IGeneralSection(const Obj<atlas_type>& atlas) : GeneralSection<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>, ISection<Point_, int, typename Alloc_::template rebind<int>::other> >(atlas) {
505: this->_bc->setChart(atlas->getChart());
506: };
507: ~IGeneralSection() {};
508: public:
509: void setChart(const chart_type& chart) {
510: this->_atlas->setChart(chart);
511: this->_atlas->allocatePoint();
512: this->_bc->setChart(chart);
513: ///this->_bc->getAtlas()->allocatePoint();
514: for(int s = 0; s < (int) this->_spaces.size(); ++s) {
515: this->_spaces[s]->setChart(chart);
516: this->_spaces[s]->allocatePoint();
517: this->_bcs[s]->setChart(chart);
518: ///this->_bcs[s]->getAtlas()->allocatePoint();
519: }
520: };
521: public:
522: bool hasNewPoints() {return this->newPoints.size() > 0;};
523: const std::set<newpoint_type>& getNewPoints() {return this->newPoints;};
524: void addPoint(const point_type& point, const int dim) {
525: if (dim == 0) return;
526: this->newPoints.insert(newpoint_type(point, dim));
527: };
528: // Returns true if the chart was changed
529: bool resizeChart(const chart_type& chart) {
530: if (!this->_atlas->reallocatePoint(chart)) return false;
531: this->_bc->reallocatePoint(chart);
532: for(int s = 0; s < (int) this->_spaces.size(); ++s) {
533: this->_spaces[s]->reallocatePoint(chart);
534: this->_bcs[s]->reallocatePoint(chart);
535: }
536: return true;
537: };
538: // Returns true if the chart was changed
539: bool reallocatePoint(const chart_type& chart) {
540: typedef typename alloc_type::template rebind<typename atlas_type::value_type>::other atlas_alloc_type;
541: const chart_type oldChart = this->getChart();
542: const int oldSize = this->sizeWithBC();
543: const values_type oldArray = this->_array;
544: const int oldAtlasSize = this->_atlas->sizeWithBC();
545: typename atlas_type::values_type oldAtlasArray;
546: if (!this->_atlas->reallocatePoint(chart, &oldAtlasArray)) return false;
547: this->_bc->reallocatePoint(chart);
548: for(int s = 0; s < (int) this->_spaces.size(); ++s) {
549: this->_spaces[s]->reallocatePoint(chart);
550: this->_bcs[s]->reallocatePoint(chart);
551: }
552: for(typename std::set<newpoint_type>::const_iterator p_iter = this->newPoints.begin(); p_iter != this->newPoints.end(); ++p_iter) {
553: this->setFiberDimension(p_iter->first, p_iter->second);
554: }
555: this->orderPoints(this->_atlas);
556: this->allocateStorage();
557: for(int i = oldChart.min(); i < oldChart.max(); ++i) {
558: const typename atlas_type::value_type& idx = this->_atlas->restrictPoint(i)[0];
559: const int dim = idx.prefix;
560: const int off = idx.index;
562: for(int d = 0; d < dim; ++d) {
563: this->_array[off+d] = oldArray[oldAtlasArray[i].index+d];
564: }
565: }
566: for(int i = 0; i < oldSize; ++i) {this->_allocator.destroy(oldArray+i);}
567: this->_allocator.deallocate(oldArray, oldSize);
568: for(int i = oldChart.min(); i < oldChart.max(); ++i) {atlas_alloc_type(this->_allocator).destroy(oldAtlasArray+i);}
569: oldAtlasArray += oldChart.min();
570: atlas_alloc_type(this->_allocator).deallocate(oldAtlasArray, oldAtlasSize);
571: this->newPoints.clear();
572: return true;
573: };
574: public:
575: void addSpace() {
576: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
577: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
578: space->setChart(this->_atlas->getChart());
579: space->allocatePoint();
580: bc->setChart(this->_bc->getChart());
581: this->_spaces.push_back(space);
582: this->_bcs.push_back(bc);
583: };
584: Obj<IGeneralSection> getFibration(const int space) const {
585: Obj<IGeneralSection> field = new IGeneralSection(this->comm(), this->debug());
586: // Obj<atlas_type> _atlas;
587: // std::vector<Obj<atlas_type> > _spaces;
588: // Obj<bc_type> _bc;
589: // std::vector<Obj<bc_type> > _bcs;
590: field->setChart(this->getChart());
591: field->addSpace();
592: const chart_type& chart = this->getChart();
594: // Copy sizes
595: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
596: const int fDim = this->getFiberDimension(*c_iter, space);
597: const int cDim = this->getConstraintDimension(*c_iter, space);
599: if (fDim) {
600: field->setFiberDimension(*c_iter, fDim);
601: field->setFiberDimension(*c_iter, fDim, 0);
602: }
603: if (cDim) {
604: field->setConstraintDimension(*c_iter, cDim);
605: field->setConstraintDimension(*c_iter, cDim, 0);
606: }
607: }
608: field->allocateStorage();
609: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
610: const chart_type& newChart = field->getChart();
612: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
613: const int cDim = field->getConstraintDimension(*c_iter);
614: const int dof[1] = {0};
616: if (cDim) {
617: field->setConstraintDof(*c_iter, dof);
618: }
619: }
620: // Copy offsets
621: newAtlas->setChart(newChart);
622: newAtlas->allocatePoint();
623: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
624: index_type idx;
626: idx.prefix = field->getFiberDimension(*c_iter);
627: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
628: for(int s = 0; s < space; ++s) {
629: idx.index += this->getFiberDimension(*c_iter, s);
630: }
631: newAtlas->addPoint(*c_iter);
632: newAtlas->updatePoint(*c_iter, &idx);
633: }
634: field->replaceStorage(this->_array, true, this->getStorageSize());
635: field->setAtlas(newAtlas);
636: return field;
637: };
638: };
639: }
641: #endif