Actual source code: IField.hh
1: #ifndef included_ALE_IField_hh
2: #define included_ALE_IField_hh
4: #ifndef included_ALE_BasicCommunication_hh
5: #include <sieve/BasicCommunication.hh>
6: #endif
8: #ifndef included_ALE_Field_hh
9: #include <sieve/Field.hh>
10: #endif
12: #ifndef included_ALE_ISieve_hh
13: #include <sieve/ISieve.hh>
14: #endif
16: // An ISection (or IntervalSection) is a section over a simple interval domain
17: namespace ALE {
18: // An IConstantSection is the simplest ISection
19: // All fibers are dimension 1
20: // All values are equal to a constant
21: // We need no value storage and no communication for completion
22: // The default value is returned for any point not in the domain
23: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
24: class IConstantSection : public ALE::ParallelObject {
25: public:
26: typedef Point_ point_type;
27: typedef Value_ value_type;
28: typedef Alloc_ alloc_type;
29: typedef Interval<point_type, alloc_type> chart_type;
30: typedef point_type index_type;
31: protected:
32: chart_type _chart;
33: value_type _value[2]; // Value and default value
34: public:
35: IConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
36: _value[1] = 0;
37: };
38: 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) {
39: _value[0] = value;
40: _value[1] = value;
41: };
42: 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) {
43: _value[0] = value;
44: _value[1] = defaultValue;
45: };
46: public: // Verifiers
47: void checkPoint(const point_type& point) const {
48: this->_chart.checkPoint(point);
49: };
50: void checkDimension(const int& dim) {
51: if (dim != 1) {
52: ostringstream msg;
53: msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
54: throw ALE::Exception(msg.str().c_str());
55: }
56: };
57: bool hasPoint(const point_type& point) const {
58: return this->_chart.hasPoint(point);
59: };
60: public: // Accessors
61: const chart_type& getChart() const {return this->_chart;};
62: void setChart(const chart_type& chart) {this->_chart = chart;};
63: void addPoint(const point_type& point) {
64: this->checkPoint(point);
65: };
66: template<typename Points>
67: void addPoint(const Obj<Points>& points) {
68: const typename Points::const_iterator pointsEnd = points->end();
69: for(typename Points::const_iterator p_iter = points->begin(); p_iter != pointsEnd; ++p_iter) {
70: this->checkPoint(*p_iter);
71: }
72: }
73: template<typename Points>
74: void addPoint(const Points& points) {
75: const typename Points::const_iterator pointsEnd = points.end();
76: for(typename Points::const_iterator p_iter = points.begin(); p_iter != pointsEnd; ++p_iter) {
77: this->checkPoint(*p_iter);
78: }
79: }
80: value_type getDefaultValue() {return this->_value[1];};
81: void setDefaultValue(const value_type value) {this->_value[1] = value;};
82: void copy(const Obj<IConstantSection>& section) {
83: const chart_type& chart = section->getChart();
85: this->_chart = chart;
86: this->_value[0] = section->restrictPoint(*chart.begin())[0];
87: this->_value[1] = section->restrictPoint(*chart.begin())[1];
88: }
89: public: // Sizes
90: ///void clear() {};
91: int getFiberDimension(const point_type& p) const {
92: if (this->hasPoint(p)) return 1;
93: return 0;
94: }
95: void setFiberDimension(const point_type& p, int dim) {
96: this->checkDimension(dim);
97: this->addPoint(p);
98: }
99: template<typename Sequence>
100: void setFiberDimension(const Obj<Sequence>& points, int dim) {
101: const typename Sequence::iterator pointsEnd = points->end();
102: for(typename Sequence::iterator p_iter = points->begin(); p_iter != pointsEnd; ++p_iter) {
103: this->setFiberDimension(*p_iter, dim);
104: }
105: }
106: void addFiberDimension(const point_type& p, int dim) {
107: if (this->hasPoint(p)) {
108: ostringstream msg;
109: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
110: throw ALE::Exception(msg.str().c_str());
111: } else {
112: this->setFiberDimension(p, dim);
113: }
114: }
115: int size(const point_type& p) {return this->getFiberDimension(p);}
116: public: // Restriction
117: void clear() {};
118: const value_type *restrictSpace() const {
119: return this->_value;
120: };
121: const value_type *restrictPoint(const point_type& p) const {
122: if (this->hasPoint(p)) {
123: return this->_value;
124: }
125: return &this->_value[1];
126: };
127: void updatePoint(const point_type&, const value_type v[]) {
128: this->_value[0] = v[0];
129: };
130: void updateAddPoint(const point_type&, const value_type v[]) {
131: this->_value[0] += v[0];
132: };
133: public:
134: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
135: ostringstream txt;
136: int rank;
138: if (comm == MPI_COMM_NULL) {
139: comm = this->comm();
140: rank = this->commRank();
141: } else {
142: MPI_Comm_rank(comm, &rank);
143: }
144: if (name == "") {
145: if(rank == 0) {
146: txt << "viewing an IConstantSection" << std::endl;
147: }
148: } else {
149: if(rank == 0) {
150: txt << "viewing IConstantSection '" << name << "'" << std::endl;
151: }
152: }
153: txt <<"["<<this->commRank()<<"]: chart " << this->_chart << std::endl;
154: txt <<"["<<this->commRank()<<"]: Value " << this->_value[0] << " Default Value " << this->_value[1] << std::endl;
155: PetscSynchronizedPrintf(comm, txt.str().c_str());
156: PetscSynchronizedFlush(comm);
157: };
158: };
160: // An IUniformSection often acts as an Atlas
161: // All fibers are the same dimension
162: // Note we can use a IConstantSection for this Atlas
163: // Each point may have a different vector
164: // Thus we need storage for values, and hence must implement completion
165: template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
166: class IUniformSection : public ALE::ParallelObject {
167: public:
168: typedef Point_ point_type;
169: typedef Value_ value_type;
170: typedef Alloc_ alloc_type;
171: typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
172: typedef IConstantSection<point_type, int, point_alloc_type> atlas_type;
173: typedef typename atlas_type::chart_type chart_type;
174: typedef point_type index_type;
175: typedef struct {value_type v[fiberDim];} fiber_type;
176: typedef value_type* values_type;
177: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
178: typedef typename atlas_alloc_type::pointer atlas_ptr;
179: protected:
180: Obj<atlas_type> _atlas;
181: values_type _array;
182: fiber_type _emptyValue;
183: alloc_type _allocator;
184: public:
185: IUniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
186: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
187: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
188: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
189: this->_array = NULL;
190: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
191: };
192: IUniformSection(MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : ParallelObject(comm, debug) {
193: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
194: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, min, max, fiberDim, debug));
195: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
196: this->_array = NULL;
197: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
198: };
199: IUniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas) {
200: int dim = fiberDim;
201: this->_atlas->update(*this->_atlas->getChart().begin(), &dim);
202: this->_array = NULL;
203: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
204: };
205: virtual ~IUniformSection() {
206: if (this->_array) {
207: for(int i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.destroy(this->_array+i);}
208: this->_array += this->getChart().min()*fiberDim;
209: this->_allocator.deallocate(this->_array, this->sizeWithBC());
210: this->_array = NULL;
211: this->_atlas = NULL;
212: }
213: };
214: public:
215: value_type *getRawArray(const int size) {
216: static value_type *array = NULL;
217: static int maxSize = 0;
219: if (size > maxSize) {
220: const value_type dummy(0);
222: if (array) {
223: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
224: this->_allocator.deallocate(array, maxSize);
225: }
226: maxSize = size;
227: array = this->_allocator.allocate(maxSize);
228: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
229: };
230: return array;
231: };
232: public: // Verifiers
233: bool hasPoint(const point_type& point) const {
234: return this->_atlas->hasPoint(point);
235: };
236: void checkDimension(const int& dim) {
237: if (dim != fiberDim) {
238: ostringstream msg;
239: msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
240: throw ALE::Exception(msg.str().c_str());
241: }
242: };
243: public: // Accessors
244: const chart_type& getChart() const {return this->_atlas->getChart();};
245: void setChart(const chart_type& chart) {
246: this->_atlas->setChart(chart);
247: int dim = fiberDim;
248: this->_atlas->updatePoint(*this->getChart().begin(), &dim);
249: };
250: bool resizeChart(const chart_type& chart) {
251: if ((chart.min() >= this->getChart().min()) && (chart.max() <= this->getChart().max())) return false;
252: this->setChart(chart);
253: return true;
254: };
255: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
256: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
257: void addPoint(const point_type& point) {
258: this->setFiberDimension(point, fiberDim);
259: };
260: template<typename Points>
261: void addPoint(const Obj<Points>& points) {
262: const typename Points::const_iterator pointsEnd = points.end();
263: for(typename Points::iterator p_iter = points->begin(); p_iter != pointsEnd; ++p_iter) {
264: this->setFiberDimension(*p_iter, fiberDim);
265: }
266: }
267: void copy(const Obj<IUniformSection>& section) {
268: this->getAtlas()->copy(section->getAtlas());
269: const chart_type& chart = section->getChart();
271: const typename chart_type::const_iterator chartEnd = chart.end();
272: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
273: this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
274: }
275: }
276: const value_type *getDefault() const {return this->_emptyValue.v;}
277: void setDefault(const value_type v[]) {for(int i = 0; i < fiberDim; ++i) {this->_emptyValue.v[i] = v[i];}}
278: public: // Sizes
279: void clear() {
280: this->zero();
281: this->_atlas->clear();
282: }
283: int getFiberDimension(const point_type& p) const {
284: return this->_atlas->restrictPoint(p)[0];
285: }
286: void setFiberDimension(const point_type& p, int dim) {
287: this->checkDimension(dim);
288: this->_atlas->addPoint(p);
289: this->_atlas->updatePoint(p, &dim);
290: }
291: template<typename Sequence>
292: void setFiberDimension(const Obj<Sequence>& points, int dim) {
293: const typename Sequence::iterator pointsEnd = points->end();
294: for(typename Sequence::iterator p_iter = points->begin(); p_iter != pointsEnd; ++p_iter) {
295: this->setFiberDimension(*p_iter, dim);
296: }
297: }
298: void setFiberDimension(const std::set<point_type>& points, int dim) {
299: const typename std::set<point_type>::const_iterator pointsEnd = points.end();
300: for(typename std::set<point_type>::iterator p_iter = points.begin(); p_iter != pointsEnd; ++p_iter) {
301: this->setFiberDimension(*p_iter, dim);
302: }
303: };
304: void addFiberDimension(const point_type& p, int dim) {
305: if (this->hasPoint(p)) {
306: ostringstream msg;
307: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
308: throw ALE::Exception(msg.str().c_str());
309: } else {
310: this->setFiberDimension(p, dim);
311: }
312: };
313: int size() const {return this->_atlas->getChart().size()*fiberDim;};
314: int sizeWithBC() const {return this->size();};
315: void allocatePoint() {
316: this->_array = this->_allocator.allocate(this->sizeWithBC());
317: this->_array -= this->getChart().min()*fiberDim;
318: for(index_type i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);}
319: };
320: bool reallocatePoint(const chart_type& chart, values_type *oldData = NULL) {
321: const chart_type oldChart = this->getChart();
322: const int oldSize = this->sizeWithBC();
323: values_type oldArray = this->_array;
324: if (!this->resizeChart(chart)) return false;
325: const int size = this->sizeWithBC();
327: this->_array = this->_allocator.allocate(size);
328: this->_array -= this->getChart().min()*fiberDim;
329: for(index_type i = this->getChart().min()*fiberDim; i < this->getChart().max()*fiberDim; ++i) {this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);}
330: for(int i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {
331: this->_array[i] = oldArray[i];
332: }
333: if (!oldData) {
334: for(index_type i = oldChart.min()*fiberDim; i < oldChart.max()*fiberDim; ++i) {this->_allocator.destroy(oldArray+i);}
335: oldArray += this->getChart().min()*fiberDim;
336: this->_allocator.deallocate(oldArray, oldSize);
337: ///std::cout << "Freed IUniformSection data" << std::endl;
338: } else {
339: ///std::cout << "Did not free IUniformSection data" << std::endl;
340: *oldData = oldArray;
341: }
342: return true;
343: };
344: template<typename Iterator, typename Extractor>
345: bool reallocatePoint(const Iterator& begin, const Iterator& end, const Extractor& extractor) {
346: point_type min = this->getChart().min();
347: point_type max = this->getChart().max()-1;
349: for(Iterator p_iter = begin; p_iter != end; ++p_iter) {
350: min = std::min(extractor(*p_iter), min);
351: max = std::max(extractor(*p_iter), max);
352: }
353: return reallocatePoint(chart_type(min, max+1));
354: }
355: public: // Restriction
356: // Zero entries
357: void zero() {
358: memset(this->_array+(this->getChart().min()*fiberDim), 0, this->sizeWithBC()* sizeof(value_type));
359: };
360: // Return a pointer to the entire contiguous storage array
361: const values_type& restrictSpace() const {
362: return this->_array;
363: };
364: // Return only the values associated to this point, not its closure
365: const value_type *restrictPoint(const point_type& p) const {
366: if (!this->hasPoint(p)) return this->_emptyValue.v;
367: return &this->_array[p*fiberDim];
368: };
369: // Update only the values associated to this point, not its closure
370: void updatePoint(const point_type& p, const value_type v[]) {
371: for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
372: this->_array[idx] = v[i];
373: }
374: };
375: // Update only the values associated to this point, not its closure
376: void updateAddPoint(const point_type& p, const value_type v[]) {
377: for(int i = 0, idx = p*fiberDim; i < fiberDim; ++i, ++idx) {
378: this->_array[idx] += v[i];
379: }
380: };
381: void updatePointAll(const point_type& p, const value_type v[]) {
382: this->updatePoint(p, v);
383: };
384: public:
385: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
386: ostringstream txt;
387: int rank;
389: if (comm == MPI_COMM_NULL) {
390: comm = this->comm();
391: rank = this->commRank();
392: } else {
393: MPI_Comm_rank(comm, &rank);
394: }
395: if (name == "") {
396: if(rank == 0) {
397: txt << "viewing an IUniformSection" << std::endl;
398: }
399: } else {
400: if(rank == 0) {
401: txt << "viewing IUniformSection '" << name << "'" << std::endl;
402: }
403: }
404: const typename atlas_type::chart_type& chart = this->_atlas->getChart();
405: values_type array = this->_array;
407: const typename atlas_type::chart_type::const_iterator chartEnd = chart.end();
408: for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chartEnd; ++p_iter) {
409: const int idx = (*p_iter)*fiberDim;
411: if (fiberDim != 0) {
412: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << fiberDim << " ";
413: for(int i = 0; i < fiberDim; i++) {
414: txt << " " << array[idx+i];
415: }
416: txt << std::endl;
417: }
418: }
419: if (chart.size() == 0) {
420: txt << "[" << this->commRank() << "]: empty" << std::endl;
421: }
422: PetscSynchronizedPrintf(comm, txt.str().c_str());
423: PetscSynchronizedFlush(comm);
424: };
425: };
426: // An ISection allows variable fiber sizes per point
427: // The Atlas is a UniformSection of dimension 1 and value type Point
428: // to hold each fiber dimension and offsets into a contiguous array
429: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_> >
430: class ISection : public Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > {
431: public:
432: typedef Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> > base;
433: typedef typename base::point_type point_type;
434: typedef typename base::value_type value_type;
435: typedef typename base::alloc_type alloc_type;
436: typedef typename base::index_type index_type;
437: typedef typename base::atlas_type atlas_type;
438: typedef typename base::chart_type chart_type;
439: typedef typename base::values_type values_type;
440: typedef typename base::atlas_alloc_type atlas_alloc_type;
441: typedef typename base::atlas_ptr atlas_ptr;
442: public:
443: ISection(MPI_Comm comm, const int debug = 0) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(comm, debug) {
444: };
445: 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) {
446: this->_atlas->setChart(chart_type(min, max));
447: this->_atlas->allocatePoint();
448: };
449: ISection(const Obj<atlas_type>& atlas) : Section<Point_, Value_, Alloc_, IUniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >(atlas) {};
450: virtual ~ISection() {};
451: public:
452: void setChart(const chart_type& chart) {
453: this->_atlas->setChart(chart);
454: this->_atlas->allocatePoint();
455: };
456: bool resizeChart(const chart_type& chart) {
457: if (!this->_atlas->reallocatePoint(chart)) return false;
458: return true;
459: };
460: bool reallocatePoint(const chart_type& chart) {
461: typedef typename atlas_type::alloc_type atlas_alloc_type;
462: const chart_type oldChart = this->getChart();
463: const int oldSize = this->sizeWithBC();
464: const values_type oldArray = this->_array;
465: const int oldAtlasSize = this->_atlas->sizeWithBC();
466: typename atlas_type::values_type oldAtlasArray;
467: if (!this->_atlas->reallocatePoint(chart, &oldAtlasArray)) return false;
469: this->orderPoints(this->_atlas);
470: this->allocateStorage();
471: for(int i = oldChart.min(); i < oldChart.max(); ++i) {
472: const typename atlas_type::value_type& idx = this->_atlas->restrictPoint(i)[0];
473: const int dim = idx.prefix;
474: const int off = idx.index;
476: for(int d = 0; d < dim; ++d) {
477: this->_array[off+d] = oldArray[oldAtlasArray[i].index+d];
478: }
479: }
480: for(int i = 0; i < oldSize; ++i) {this->_allocator.destroy(oldArray+i);}
481: this->_allocator.deallocate(oldArray, oldSize);
482: for(int i = oldChart.min(); i < oldChart.max(); ++i) {atlas_alloc_type(this->_allocator).destroy(oldAtlasArray+i);}
483: oldAtlasArray += oldChart.min();
484: atlas_alloc_type(this->_allocator).deallocate(oldAtlasArray, oldAtlasSize);
485: ///std::cout << "In ISection, Freed IUniformSection data" << std::endl;
486: return true;
487: };
488: public:
489: // Return the free values on a point
490: // This is overridden, because the one in Section cannot be const due to problem in the interface with UniformSection
491: const value_type *restrictPoint(const point_type& p) const {
492: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
493: };
494: };
495: // IGeneralSection will support BC on a subset of unknowns on a point
496: // We use a separate constraint Atlas to mark constrained dofs on a point
497: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_> >
498: 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> > {
499: public:
500: 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;
501: typedef typename base::point_type point_type;
502: typedef typename base::value_type value_type;
503: typedef typename base::alloc_type alloc_type;
504: typedef typename base::index_type index_type;
505: typedef typename base::atlas_type atlas_type;
506: typedef typename base::bc_type bc_type;
507: typedef typename base::chart_type chart_type;
508: typedef typename base::values_type values_type;
509: typedef typename base::atlas_alloc_type atlas_alloc_type;
510: typedef typename base::atlas_ptr atlas_ptr;
511: typedef typename base::bc_alloc_type bc_alloc_type;
512: typedef typename base::bc_ptr bc_ptr;
513: typedef std::pair<point_type,int> newpoint_type;
514: protected:
515: std::set<newpoint_type> newPoints;
516: public:
517: 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) {};
518: 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) {
519: this->_atlas->setChart(chart_type(min, max));
520: this->_atlas->allocatePoint();
521: this->_bc->setChart(chart_type(min, max));
522: };
523: 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) {
524: this->_bc->setChart(atlas->getChart());
525: };
526: ~IGeneralSection() {};
527: public:
528: void setChart(const chart_type& chart) {
529: this->_atlas->setChart(chart);
530: this->_atlas->allocatePoint();
531: this->_bc->setChart(chart);
532: ///this->_bc->getAtlas()->allocatePoint();
533: for(int s = 0; s < (int) this->_spaces.size(); ++s) {
534: this->_spaces[s]->setChart(chart);
535: this->_spaces[s]->allocatePoint();
536: this->_bcs[s]->setChart(chart);
537: ///this->_bcs[s]->getAtlas()->allocatePoint();
538: }
539: };
540: public:
541: bool hasNewPoints() {return this->newPoints.size() > 0;};
542: const std::set<newpoint_type>& getNewPoints() {return this->newPoints;};
543: void addPoint(const point_type& point, const int dim) {
544: if (dim == 0) return;
545: this->newPoints.insert(newpoint_type(point, dim));
546: };
547: // Returns true if the chart was changed
548: bool resizeChart(const chart_type& chart) {
549: if (!this->_atlas->reallocatePoint(chart)) return false;
550: this->_bc->reallocatePoint(chart);
551: for(int s = 0; s < (int) this->_spaces.size(); ++s) {
552: this->_spaces[s]->reallocatePoint(chart);
553: this->_bcs[s]->reallocatePoint(chart);
554: }
555: return true;
556: };
557: // Returns true if the chart was changed
558: bool reallocatePoint(const chart_type& chart) {
559: typedef typename alloc_type::template rebind<typename atlas_type::value_type>::other atlas_alloc_type;
560: const chart_type oldChart = this->getChart();
561: const int oldSize = this->sizeWithBC();
562: const values_type oldArray = this->_array;
563: const int oldAtlasSize = this->_atlas->sizeWithBC();
564: typename atlas_type::values_type oldAtlasArray;
565: if (!this->_atlas->reallocatePoint(chart, &oldAtlasArray)) return false;
566: this->_bc->reallocatePoint(chart);
567: for(int s = 0; s < (int) this->_spaces.size(); ++s) {
568: this->_spaces[s]->reallocatePoint(chart);
569: this->_bcs[s]->reallocatePoint(chart);
570: }
571: const typename std::set<newpoint_type>::const_iterator newPointsEnd = this->newPoints.end();
572: for(typename std::set<newpoint_type>::const_iterator p_iter = this->newPoints.begin(); p_iter != newPointsEnd; ++p_iter) {
573: this->setFiberDimension(p_iter->first, p_iter->second);
574: }
575: this->orderPoints(this->_atlas);
576: this->allocateStorage();
577: for(int i = oldChart.min(); i < oldChart.max(); ++i) {
578: const typename atlas_type::value_type& idx = this->_atlas->restrictPoint(i)[0];
579: const int dim = idx.prefix;
580: const int off = idx.index;
582: for(int d = 0; d < dim; ++d) {
583: this->_array[off+d] = oldArray[oldAtlasArray[i].index+d];
584: }
585: }
586: for(int i = 0; i < oldSize; ++i) {this->_allocator.destroy(oldArray+i);}
587: this->_allocator.deallocate(oldArray, oldSize);
588: for(int i = oldChart.min(); i < oldChart.max(); ++i) {atlas_alloc_type(this->_allocator).destroy(oldAtlasArray+i);}
589: oldAtlasArray += oldChart.min();
590: atlas_alloc_type(this->_allocator).deallocate(oldAtlasArray, oldAtlasSize);
591: this->newPoints.clear();
592: return true;
593: };
594: public:
595: void addSpace() {
596: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
597: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
598: space->setChart(this->_atlas->getChart());
599: space->allocatePoint();
600: bc->setChart(this->_bc->getChart());
601: this->_spaces.push_back(space);
602: this->_bcs.push_back(bc);
603: };
604: Obj<IGeneralSection> getFibration(const int space) const {
605: Obj<IGeneralSection> field = new IGeneralSection(this->comm(), this->debug());
606: // Obj<atlas_type> _atlas;
607: // std::vector<Obj<atlas_type> > _spaces;
608: // Obj<bc_type> _bc;
609: // std::vector<Obj<bc_type> > _bcs;
610: field->setChart(this->getChart());
611: field->addSpace();
612: const chart_type& chart = this->getChart();
614: // Copy sizes
615: const typename chart_type::const_iterator chartEnd = chart.end();
616: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
617: const int fDim = this->getFiberDimension(*c_iter, space);
618: const int cDim = this->getConstraintDimension(*c_iter, space);
620: if (fDim) {
621: field->setFiberDimension(*c_iter, fDim);
622: field->setFiberDimension(*c_iter, fDim, 0);
623: }
624: if (cDim) {
625: field->setConstraintDimension(*c_iter, cDim);
626: field->setConstraintDimension(*c_iter, cDim, 0);
627: }
628: }
629: field->allocateStorage();
630: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
631: const chart_type& newChart = field->getChart();
633: const typename chart_type::const_iterator newChartEnd = newChart.end();
634: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChartEnd; ++c_iter) {
635: const int cDim = field->getConstraintDimension(*c_iter);
636: const int dof[1] = {0};
638: if (cDim) {
639: field->setConstraintDof(*c_iter, dof);
640: }
641: }
642: // Copy offsets
643: newAtlas->setChart(newChart);
644: newAtlas->allocatePoint();
645: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChartEnd; ++c_iter) {
646: index_type idx;
648: idx.prefix = field->getFiberDimension(*c_iter);
649: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
650: for(int s = 0; s < space; ++s) {
651: idx.index += this->getFiberDimension(*c_iter, s);
652: }
653: newAtlas->addPoint(*c_iter);
654: newAtlas->updatePoint(*c_iter, &idx);
655: }
656: field->replaceStorage(this->_array, true, this->getStorageSize());
657: field->setAtlas(newAtlas);
658: return field;
659: };
660: };
662: class SectionSerializer {
663: public:
664: template<typename Point_, typename Value_>
665: static void writeSection(std::ofstream& fs, IConstantSection<Point_, Value_>& section) {
666: MPIMover<Value_> mover(section.comm());
668: if (section.commRank() == 0) {
669: // Write local
670: fs << section.getChart().min() << " " << section.getChart().max() << std::endl;
671: fs.precision(15);
672: fs << section.restrictPoint(section.getChart().min())[0] << " ";
673: fs << section.getDefaultValue() << std::endl;
674: // Receive and write remote
675: for(int p = 1; p < section.commSize(); ++p) {
676: PetscInt sizes[2];
677: Value_ values[2];
678: MPI_Status status;
681: MPI_Recv(sizes, 2, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
682: MPI_Recv(values, 2, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
683: fs << sizes[0] << " " << sizes[1] << std::endl;
684: fs.precision(15);
685: fs << values[0] << " " << values[1] << std::endl;
686: }
687: } else {
688: // Send remote
689: PetscInt sizes[2];
690: Value_ values[2];
693: sizes[0] = section.getChart().min();
694: sizes[1] = section.getChart().max();
695: values[0] = section.restrictPoint(section.getChart().min())[0];
696: values[1] = section.getDefaultValue();
697: MPI_Send(sizes, 2, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
698: MPI_Send(values, 2, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
699: }
700: }
701: template<typename Point_, typename Value_, int fiberDim>
702: static void writeSection(std::ofstream& fs, IUniformSection<Point_, Value_, fiberDim>& section) {
703: typedef typename IUniformSection<Point_, Value_, fiberDim>::index_type index_type;
704: typedef typename IUniformSection<Point_, Value_, fiberDim>::value_type value_type;
705: MPIMover<value_type> mover(section.comm());
706: index_type min = section.getChart().min();
707: index_type max = section.getChart().max();
709: // Write atlas
710: writeSection(fs, *section.getAtlas());
711: if (section.commRank() == 0) {
712: // Write local values
713: fs.precision(15);
714: for(index_type p = min; p < max; ++p) {
715: const value_type *values = section.restrictPoint(p);
717: for(int i = 0; i < fiberDim; ++i) {
718: fs << values[i] << std::endl;
719: }
720: }
721: // Write empty value
722: const value_type *defValue = section.getDefault();
724: for(int i = 0; i < fiberDim; ++i) {
725: if (i > 0) fs << " ";
726: fs << defValue[i];
727: }
728: fs << std::endl;
729: // Receive and write remote
730: for(int p = 1; p < section.commSize(); ++p) {
731: PetscInt size;
732: value_type *values;
733: value_type emptyValues[fiberDim];
734: MPI_Status status;
737: MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
738: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
739: MPI_Recv(values, size*fiberDim, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
740: for(PetscInt v = 0; v < size; ++v) {
741: fs << values[v] << std::endl;
742: }
743: PetscFree(values);CHKERRXX(ierr);
744: MPI_Recv(emptyValues, fiberDim, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
745: for(int i = 0; i < fiberDim; ++i) {
746: if (i > 0) fs << " ";
747: fs << emptyValues[i];
748: }
749: fs << std::endl;
750: }
751: } else {
752: // Send remote
753: PetscInt size = section.getChart().size();
754: PetscInt v = 0;
755: const value_type *defValue = section.getDefault();
756: value_type *values;
757: value_type emptyValues[fiberDim];
758: PetscErrorCode ierr;
760: MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
761: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
762: for(index_type p = min; p < max; ++p) {
763: const value_type *val = section.restrictPoint(p);
765: for(int i = 0; i < fiberDim; ++i, ++v) {
766: values[v] = val[i];
767: }
768: }
769: MPI_Send(values, size*fiberDim, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
770: for(int i = 0; i < fiberDim; ++i) {emptyValues[i] = ((value_type *) &defValue[i])[0];}
771: MPI_Send(emptyValues, fiberDim, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
772: }
773: }
774: template<typename Point_, typename Value_>
775: static void writeSection(std::ofstream& fs, ISection<Point_, Value_>& section) {
776: typedef typename ISection<Point_, Value_>::point_type point_type;
777: typedef typename ISection<Point_, Value_>::value_type value_type;
778: MPIMover<value_type> mover(section.comm());
779: point_type min = section.getChart().min();
780: point_type max = section.getChart().max();
782: // Write atlas
783: writeSection(fs, *section.getAtlas());
784: if (section.commRank() == 0) {
785: // Write local values
786: fs.precision(15);
787: for(point_type p = min; p < max; ++p) {
788: const int fiberDim = section.getFiberDimension(p);
789: const value_type *values = section.restrictPoint(p);
791: for(int i = 0; i < fiberDim; ++i) {
792: fs << values[i] << std::endl;
793: }
794: }
795: // Receive and write remote
796: for(int p = 1; p < section.commSize(); ++p) {
797: PetscInt size;
798: value_type *values;
799: MPI_Status status;
802: MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
803: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
804: MPI_Recv(values, size, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
805: for(PetscInt v = 0; v < size; ++v) {
806: fs << values[v] << std::endl;
807: }
808: PetscFree(values);CHKERRXX(ierr);
809: }
810: } else {
811: // Send remote
812: PetscInt size = section.size();
813: PetscInt v = 0;
814: value_type *values;
817: MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
818: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
819: for(point_type p = min; p < max; ++p) {
820: const int fiberDim = section.getFiberDimension(p);
821: const value_type *val = section.restrictPoint(p);
823: for(int i = 0; i < fiberDim; ++i, ++v) {
824: values[v] = val[i];
825: }
826: }
827: MPI_Send(values, size, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
828: }
829: }
830: template<typename Point_, typename Value_>
831: static void writeSection(std::ofstream& fs, IGeneralSection<Point_, Value_>& section) {
832: typedef typename IGeneralSection<Point_, Value_>::point_type point_type;
833: typedef typename IGeneralSection<Point_, Value_>::value_type value_type;
834: MPIMover<value_type> mover(section.comm());
835: point_type min = section.getChart().min();
836: point_type max = section.getChart().max();
838: // Write atlas
839: writeSection(fs, *section.getAtlas());
840: if (section.commRank() == 0) {
841: // Write local values
842: fs.precision(15);
843: for(point_type p = min; p < max; ++p) {
844: const int fiberDim = section.getFiberDimension(p);
845: const value_type *values = section.restrictPoint(p);
847: for(int i = 0; i < fiberDim; ++i) {
848: fs << values[i] << std::endl;
849: }
850: }
851: // Receive and write remote
852: for(int p = 1; p < section.commSize(); ++p) {
853: PetscInt size;
854: value_type *values;
855: MPI_Status status;
858: MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
859: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
860: MPI_Recv(values, size, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
861: for(PetscInt v = 0; v < size; ++v) {
862: fs << values[v] << std::endl;
863: }
864: PetscFree(values);CHKERRXX(ierr);
865: }
866: } else {
867: // Send remote
868: PetscInt size = section.sizeWithBC();
869: PetscInt v = 0;
870: value_type *values;
873: MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
874: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
875: for(point_type p = min; p < max; ++p) {
876: const int fiberDim = section.getFiberDimension(p);
877: const value_type *val = section.restrictPoint(p);
879: for(int i = 0; i < fiberDim; ++i, ++v) {
880: values[v] = val[i];
881: }
882: }
883: MPI_Send(values, size, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
884: }
885: // Write BC
886: writeSection(fs, *section.getBC());
887: // Write spaces
888: // std::vector<Obj<atlas_type> > _spaces;
889: // std::vector<Obj<bc_type> > _bcs;
890: }
891: template<typename Point_, typename Value_>
892: static void loadSection(std::ifstream& fs, IConstantSection<Point_, Value_>& section) {
893: typedef typename IConstantSection<Point_, Value_>::index_type index_type;
894: typedef typename IConstantSection<Point_, Value_>::value_type value_type;
895: MPIMover<value_type> mover(section.comm());
896: index_type min, max;
897: value_type val;
899: if (section.commRank() == 0) {
900: // Load local
901: fs >> min;
902: fs >> max;
903: section.setChart(typename IConstantSection<Point_, Value_>::chart_type(min, max));
904: fs >> val;
905: section.updatePoint(min, &val);
906: fs >> val;
907: section.setDefaultValue(val);
908: // Load and send remote
909: for(int p = 1; p < section.commSize(); ++p) {
910: PetscInt sizes[2];
911: value_type values[2];
914: fs >> sizes[0];
915: fs >> sizes[1];
916: fs >> values[0];
917: fs >> values[1];
918: MPI_Send(sizes, 2, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
919: MPI_Send(values, 2, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
920: }
921: } else {
922: // Load remote
923: PetscInt sizes[2];
924: value_type values[2];
925: MPI_Status status;
928: MPI_Recv(sizes, 2, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
929: section.setChart(typename IConstantSection<Point_, Value_>::chart_type(sizes[0], sizes[1]));
930: MPI_Recv(values, 2, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
931: section.updatePoint(min, values);
932: section.setDefaultValue(values[1]);
933: }
934: }
935: template<typename Point_, typename Value_, int fiberDim>
936: static void loadSection(std::ifstream& fs, IUniformSection<Point_, Value_, fiberDim>& section) {
937: typedef typename IUniformSection<Point_, Value_, fiberDim>::index_type index_type;
938: typedef typename IUniformSection<Point_, Value_, fiberDim>::value_type value_type;
939: MPIMover<value_type> mover(section.comm());
940: // Load atlas
941: loadSection(fs, *section.getAtlas());
942: section.allocatePoint();
943: index_type min = section.getChart().min();
944: index_type max = section.getChart().max();
946: if (section.commRank() == 0) {
947: // Load local values
948: for(index_type p = min; p < max; ++p) {
949: value_type values[fiberDim];
951: for(int i = 0; i < fiberDim; ++i) {
952: typename IUniformSection<Point_, Value_, fiberDim>::value_type value;
954: fs >> value;
955: values[i] = value;
956: }
957: section.updatePoint(p, values);
958: }
959: // Load empty value
960: value_type defValue[fiberDim];
962: for(int i = 0; i < fiberDim; ++i) {
963: fs >> defValue[i];
964: }
965: section.setDefault(defValue);
966: // Load and send remote
967: for(int pr = 1; pr < section.commSize(); ++pr) {
968: PetscInt size = section.getChart().size();
969: PetscInt v = 0;
970: value_type *values;
971: value_type emptyValues[fiberDim];
974: MPI_Send(&size, 1, MPIU_INT, pr, 1, section.comm());CHKERRXX(ierr);
975: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
976: for(index_type p = min; p < max; ++p) {
977: for(int i = 0; i < fiberDim; ++i, ++v) {
978: fs >> values[v];
979: }
980: }
981: MPI_Send(values, size*fiberDim, mover.getType(), pr, 1, section.comm());CHKERRXX(ierr);
982: for(int i = 0; i < fiberDim; ++i) {
983: fs >> emptyValues[i];
984: }
985: MPI_Send(emptyValues, fiberDim, mover.getType(), pr, 1, section.comm());CHKERRXX(ierr);
986: }
987: } else {
988: // Load remote
989: PetscInt size;
990: value_type *values;
991: value_type emptyValues[fiberDim];
992: value_type pvalues[fiberDim];
993: MPI_Status status;
994: PetscInt v = 0;
997: assert(sizeof(value_type) <= sizeof(PetscScalar));
998: MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
999: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
1000: MPI_Recv(values, size*fiberDim, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1001: for(index_type p = min; p < max; ++p) {
1002: for(int i = 0; i < fiberDim; ++i, ++v) {
1003: pvalues[i] = values[v];
1004: }
1005: section.updatePoint(p, pvalues);
1006: }
1007: PetscFree(values);CHKERRXX(ierr);
1008: MPI_Recv(emptyValues, fiberDim, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1009: section.setDefault(emptyValues);
1010: }
1011: }
1012: template<typename Point_, typename Value_>
1013: static void loadSection(std::ifstream& fs, ISection<Point_, Value_>& section) {
1014: typedef typename ISection<Point_, Value_>::point_type point_type;
1015: typedef typename ISection<Point_, Value_>::value_type value_type;
1016: MPIMover<value_type> mover(section.comm());
1017: // Load atlas
1018: loadSection(fs, *section.getAtlas());
1019: section.allocatePoint();
1020: point_type min = section.getChart().min();
1021: point_type max = section.getChart().max();
1022: int maxDim = -1;
1024: if (section.commRank() == 0) {
1025: // Load local values
1026: for(point_type p = min; p < max; ++p) {
1027: maxDim = std::max(maxDim, section.getFiberDimension(p));
1028: }
1029: value_type *values = new value_type[maxDim];
1030: for(point_type p = min; p < max; ++p) {
1031: const int fiberDim = section.getFiberDimension(p);
1033: for(int i = 0; i < fiberDim; ++i) {
1034: fs >> values[i];
1035: }
1036: section.updatePoint(p, values);
1037: }
1038: delete [] values;
1039: // Load and send remote
1040: for(int p = 1; p < section.commSize(); ++p) {
1041: PetscInt size = section.size();
1042: value_type *values;
1045: MPI_Send(&size, 1, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
1046: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1047: for(PetscInt v = 0; v < size; ++v) {
1048: fs >> values[v];
1049: }
1050: MPI_Send(values, size, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
1051: }
1052: } else {
1053: // Load remote
1054: PetscInt size;
1055: value_type *values;
1056: MPI_Status status;
1057: PetscInt maxDim = 0;
1058: PetscInt off = 0;
1061: MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
1062: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1063: MPI_Recv(values, size, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1064: for(point_type p = min; p < max; ++p) {
1065: maxDim = std::max(maxDim, section.getFiberDimension(p));
1066: }
1067: value_type *pvalues = new value_type[maxDim];
1068: for(point_type p = min; p < max; ++p) {
1069: const int fiberDim = section.getFiberDimension(p);
1071: for(int i = 0; i < fiberDim; ++i, ++off) {
1072: pvalues[i] = values[off];
1073: }
1074: section.updatePoint(p, pvalues);
1075: }
1076: delete [] pvalues;
1077: PetscFree(values);CHKERRXX(ierr);
1078: }
1079: }
1080: template<typename Point_, typename Value_>
1081: static void loadSection(std::ifstream& fs, IGeneralSection<Point_, Value_>& section) {
1082: typedef typename IGeneralSection<Point_, Value_>::point_type point_type;
1083: typedef typename IGeneralSection<Point_, Value_>::value_type value_type;
1084: MPIMover<value_type> mover(section.comm());
1085: // Load atlas
1086: loadSection(fs, *section.getAtlas());
1087: section.allocatePoint();
1088: point_type min = section.getChart().min();
1089: point_type max = section.getChart().max();
1090: int maxDim = -1;
1092: if (section.commRank() == 0) {
1093: // Load local values
1094: for(point_type p = min; p < max; ++p) {
1095: maxDim = std::max(maxDim, section.getFiberDimension(p));
1096: }
1097: value_type *values = new value_type[maxDim];
1098: for(point_type p = min; p < max; ++p) {
1099: const int fiberDim = section.getFiberDimension(p);
1101: for(int i = 0; i < fiberDim; ++i) {
1102: fs >> values[i];
1103: }
1104: section.updatePoint(p, values);
1105: }
1106: delete [] values;
1107: // Load and send remote
1108: for(int p = 1; p < section.commSize(); ++p) {
1109: PetscInt size = section.sizeWithBC();
1110: value_type *values;
1113: MPI_Send(&size, 1, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
1114: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1115: for(PetscInt v = 0; v < size; ++v) {
1116: fs >> values[v];
1117: }
1118: MPI_Send(values, size, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
1119: }
1120: } else {
1121: // Load remote
1122: PetscInt size;
1123: value_type *values;
1124: MPI_Status status;
1125: PetscInt maxDim = 0;
1126: PetscInt off = 0;
1129: assert(sizeof(value_type) <= sizeof(PetscScalar));
1130: MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
1131: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1132: MPI_Recv(values, size, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1133: for(point_type p = min; p < max; ++p) {
1134: maxDim = std::max(maxDim, section.getFiberDimension(p));
1135: }
1136: value_type *pvalues = new value_type[maxDim];
1137: for(point_type p = min; p < max; ++p) {
1138: const int fiberDim = section.getFiberDimension(p);
1140: for(int i = 0; i < fiberDim; ++i, ++off) {
1141: pvalues[i] = values[off];
1142: }
1143: section.updatePoint(p, pvalues);
1144: }
1145: delete [] pvalues;
1146: PetscFree(values);CHKERRXX(ierr);
1147: }
1148: // Load BC
1149: loadSection(fs, *section.getBC());
1150: // Load spaces
1151: // std::vector<Obj<atlas_type> > _spaces;
1152: // std::vector<Obj<bc_type> > _bcs;
1153: }
1154: };
1155: }
1157: #endif