Actual source code: IField.hh
petsc-3.4.5 2014-06-29
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(int comp = 1) {
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->_comps.push_back(comp);
602: this->_spaces.push_back(space);
603: this->_bcs.push_back(bc);
604: };
605: Obj<IGeneralSection> getFibration(const int space) const {
606: Obj<IGeneralSection> field = new IGeneralSection(this->comm(), this->debug());
607: // Obj<atlas_type> _atlas;
608: // std::vector<Obj<atlas_type> > _spaces;
609: // Obj<bc_type> _bc;
610: // std::vector<Obj<bc_type> > _bcs;
611: field->setChart(this->getChart());
612: field->addSpace();
613: const chart_type& chart = this->getChart();
615: // Copy sizes
616: const typename chart_type::const_iterator chartEnd = chart.end();
617: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chartEnd; ++c_iter) {
618: const int fDim = this->getFiberDimension(*c_iter, space);
619: const int cDim = this->getConstraintDimension(*c_iter, space);
621: if (fDim) {
622: field->setFiberDimension(*c_iter, fDim);
623: field->setFiberDimension(*c_iter, fDim, 0);
624: }
625: if (cDim) {
626: field->setConstraintDimension(*c_iter, cDim);
627: field->setConstraintDimension(*c_iter, cDim, 0);
628: }
629: }
630: field->allocateStorage();
631: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
632: const chart_type& newChart = field->getChart();
634: const typename chart_type::const_iterator newChartEnd = newChart.end();
635: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChartEnd; ++c_iter) {
636: const int cDim = field->getConstraintDimension(*c_iter);
637: const int dof[1] = {0};
639: if (cDim) {
640: field->setConstraintDof(*c_iter, dof);
641: }
642: }
643: // Copy offsets
644: newAtlas->setChart(newChart);
645: newAtlas->allocatePoint();
646: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChartEnd; ++c_iter) {
647: index_type idx;
649: idx.prefix = field->getFiberDimension(*c_iter);
650: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
651: for(int s = 0; s < space; ++s) {
652: idx.index += this->getFiberDimension(*c_iter, s);
653: }
654: newAtlas->addPoint(*c_iter);
655: newAtlas->updatePoint(*c_iter, &idx);
656: }
657: field->replaceStorage(this->_array, true, this->getStorageSize());
658: field->setAtlas(newAtlas);
659: return field;
660: };
661: };
663: class SectionSerializer {
664: public:
665: template<typename Point_, typename Value_>
666: static void writeSection(std::ofstream& fs, IConstantSection<Point_, Value_>& section) {
667: MPIMover<Value_> mover(section.comm());
669: if (section.commRank() == 0) {
670: // Write local
671: fs << section.getChart().min() << " " << section.getChart().max() << std::endl;
672: fs.precision(15);
673: fs << section.restrictPoint(section.getChart().min())[0] << " ";
674: fs << section.getDefaultValue() << std::endl;
675: // Receive and write remote
676: for(int p = 1; p < section.commSize(); ++p) {
677: PetscInt sizes[2];
678: Value_ values[2];
679: MPI_Status status;
682: MPI_Recv(sizes, 2, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
683: MPI_Recv(values, 2, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
684: fs << sizes[0] << " " << sizes[1] << std::endl;
685: fs.precision(15);
686: fs << values[0] << " " << values[1] << std::endl;
687: }
688: } else {
689: // Send remote
690: PetscInt sizes[2];
691: Value_ values[2];
694: sizes[0] = section.getChart().min();
695: sizes[1] = section.getChart().max();
696: values[0] = section.restrictPoint(section.getChart().min())[0];
697: values[1] = section.getDefaultValue();
698: MPI_Send(sizes, 2, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
699: MPI_Send(values, 2, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
700: }
701: }
702: template<typename Point_, typename Value_, int fiberDim>
703: static void writeSection(std::ofstream& fs, IUniformSection<Point_, Value_, fiberDim>& section) {
704: typedef typename IUniformSection<Point_, Value_, fiberDim>::index_type index_type;
705: typedef typename IUniformSection<Point_, Value_, fiberDim>::value_type value_type;
706: MPIMover<value_type> mover(section.comm());
707: index_type min = section.getChart().min();
708: index_type max = section.getChart().max();
710: // Write atlas
711: writeSection(fs, *section.getAtlas());
712: if (section.commRank() == 0) {
713: // Write local values
714: fs.precision(15);
715: for(index_type p = min; p < max; ++p) {
716: const value_type *values = section.restrictPoint(p);
718: for(int i = 0; i < fiberDim; ++i) {
719: fs << values[i] << std::endl;
720: }
721: }
722: // Write empty value
723: const value_type *defValue = section.getDefault();
725: for(int i = 0; i < fiberDim; ++i) {
726: if (i > 0) fs << " ";
727: fs << defValue[i];
728: }
729: fs << std::endl;
730: // Receive and write remote
731: for(int p = 1; p < section.commSize(); ++p) {
732: PetscInt size;
733: value_type *values;
734: value_type emptyValues[fiberDim];
735: MPI_Status status;
738: MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
739: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
740: MPI_Recv(values, size*fiberDim, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
741: for(PetscInt v = 0; v < size; ++v) {
742: fs << values[v] << std::endl;
743: }
744: PetscFree(values);CHKERRXX(ierr);
745: MPI_Recv(emptyValues, fiberDim, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
746: for(int i = 0; i < fiberDim; ++i) {
747: if (i > 0) fs << " ";
748: fs << emptyValues[i];
749: }
750: fs << std::endl;
751: }
752: } else {
753: // Send remote
754: PetscInt size = section.getChart().size();
755: PetscInt v = 0;
756: const value_type *defValue = section.getDefault();
757: value_type *values;
758: value_type emptyValues[fiberDim];
759: PetscErrorCode ierr;
761: MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
762: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
763: for(index_type p = min; p < max; ++p) {
764: const value_type *val = section.restrictPoint(p);
766: for(int i = 0; i < fiberDim; ++i, ++v) {
767: values[v] = val[i];
768: }
769: }
770: MPI_Send(values, size*fiberDim, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
771: for(int i = 0; i < fiberDim; ++i) {emptyValues[i] = ((value_type *) &defValue[i])[0];}
772: MPI_Send(emptyValues, fiberDim, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
773: }
774: }
775: template<typename Point_, typename Value_>
776: static void writeSection(std::ofstream& fs, ISection<Point_, Value_>& section) {
777: typedef typename ISection<Point_, Value_>::point_type point_type;
778: typedef typename ISection<Point_, Value_>::value_type value_type;
779: MPIMover<value_type> mover(section.comm());
780: point_type min = section.getChart().min();
781: point_type max = section.getChart().max();
783: // Write atlas
784: writeSection(fs, *section.getAtlas());
785: if (section.commRank() == 0) {
786: // Write local values
787: fs.precision(15);
788: for(point_type p = min; p < max; ++p) {
789: const int fiberDim = section.getFiberDimension(p);
790: const value_type *values = section.restrictPoint(p);
792: for(int i = 0; i < fiberDim; ++i) {
793: fs << values[i] << std::endl;
794: }
795: }
796: // Receive and write remote
797: for(int p = 1; p < section.commSize(); ++p) {
798: PetscInt size;
799: value_type *values;
800: MPI_Status status;
803: MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
804: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
805: MPI_Recv(values, size, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
806: for(PetscInt v = 0; v < size; ++v) {
807: fs << values[v] << std::endl;
808: }
809: PetscFree(values);CHKERRXX(ierr);
810: }
811: } else {
812: // Send remote
813: PetscInt size = section.size();
814: PetscInt v = 0;
815: value_type *values;
818: MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
819: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
820: for(point_type p = min; p < max; ++p) {
821: const int fiberDim = section.getFiberDimension(p);
822: const value_type *val = section.restrictPoint(p);
824: for(int i = 0; i < fiberDim; ++i, ++v) {
825: values[v] = val[i];
826: }
827: }
828: MPI_Send(values, size, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
829: }
830: }
831: template<typename Point_, typename Value_>
832: static void writeSection(std::ofstream& fs, IGeneralSection<Point_, Value_>& section) {
833: typedef typename IGeneralSection<Point_, Value_>::point_type point_type;
834: typedef typename IGeneralSection<Point_, Value_>::value_type value_type;
835: MPIMover<value_type> mover(section.comm());
836: point_type min = section.getChart().min();
837: point_type max = section.getChart().max();
839: // Write atlas
840: writeSection(fs, *section.getAtlas());
841: if (section.commRank() == 0) {
842: // Write local values
843: fs.precision(15);
844: for(point_type p = min; p < max; ++p) {
845: const int fiberDim = section.getFiberDimension(p);
846: const value_type *values = section.restrictPoint(p);
848: for(int i = 0; i < fiberDim; ++i) {
849: fs << values[i] << std::endl;
850: }
851: }
852: // Receive and write remote
853: for(int p = 1; p < section.commSize(); ++p) {
854: PetscInt size;
855: value_type *values;
856: MPI_Status status;
859: MPI_Recv(&size, 1, MPIU_INT, p, 1, section.comm(), &status);CHKERRXX(ierr);
860: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
861: MPI_Recv(values, size, mover.getType(), p, 1, section.comm(), &status);CHKERRXX(ierr);
862: for(PetscInt v = 0; v < size; ++v) {
863: fs << values[v] << std::endl;
864: }
865: PetscFree(values);CHKERRXX(ierr);
866: }
867: } else {
868: // Send remote
869: PetscInt size = section.sizeWithBC();
870: PetscInt v = 0;
871: value_type *values;
874: MPI_Send(&size, 1, MPIU_INT, 0, 1, section.comm());CHKERRXX(ierr);
875: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
876: for(point_type p = min; p < max; ++p) {
877: const int fiberDim = section.getFiberDimension(p);
878: const value_type *val = section.restrictPoint(p);
880: for(int i = 0; i < fiberDim; ++i, ++v) {
881: values[v] = val[i];
882: }
883: }
884: MPI_Send(values, size, mover.getType(), 0, 1, section.comm());CHKERRXX(ierr);
885: }
886: // Write BC
887: writeSection(fs, *section.getBC());
888: // Write spaces
889: // std::vector<Obj<atlas_type> > _spaces;
890: // std::vector<Obj<bc_type> > _bcs;
891: }
892: template<typename Point_, typename Value_>
893: static void loadSection(std::ifstream& fs, IConstantSection<Point_, Value_>& section) {
894: typedef typename IConstantSection<Point_, Value_>::index_type index_type;
895: typedef typename IConstantSection<Point_, Value_>::value_type value_type;
896: MPIMover<value_type> mover(section.comm());
897: index_type min, max;
898: value_type val;
900: if (section.commRank() == 0) {
901: // Load local
902: fs >> min;
903: fs >> max;
904: section.setChart(typename IConstantSection<Point_, Value_>::chart_type(min, max));
905: fs >> val;
906: section.updatePoint(min, &val);
907: fs >> val;
908: section.setDefaultValue(val);
909: // Load and send remote
910: for(int p = 1; p < section.commSize(); ++p) {
911: PetscInt sizes[2];
912: value_type values[2];
915: fs >> sizes[0];
916: fs >> sizes[1];
917: fs >> values[0];
918: fs >> values[1];
919: MPI_Send(sizes, 2, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
920: MPI_Send(values, 2, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
921: }
922: } else {
923: // Load remote
924: PetscInt sizes[2];
925: value_type values[2];
926: MPI_Status status;
929: MPI_Recv(sizes, 2, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
930: section.setChart(typename IConstantSection<Point_, Value_>::chart_type(sizes[0], sizes[1]));
931: MPI_Recv(values, 2, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
932: section.updatePoint(min, values);
933: section.setDefaultValue(values[1]);
934: }
935: }
936: template<typename Point_, typename Value_, int fiberDim>
937: static void loadSection(std::ifstream& fs, IUniformSection<Point_, Value_, fiberDim>& section) {
938: typedef typename IUniformSection<Point_, Value_, fiberDim>::index_type index_type;
939: typedef typename IUniformSection<Point_, Value_, fiberDim>::value_type value_type;
940: MPIMover<value_type> mover(section.comm());
941: // Load atlas
942: loadSection(fs, *section.getAtlas());
943: section.allocatePoint();
944: index_type min = section.getChart().min();
945: index_type max = section.getChart().max();
947: if (section.commRank() == 0) {
948: // Load local values
949: for(index_type p = min; p < max; ++p) {
950: value_type values[fiberDim];
952: for(int i = 0; i < fiberDim; ++i) {
953: typename IUniformSection<Point_, Value_, fiberDim>::value_type value;
955: fs >> value;
956: values[i] = value;
957: }
958: section.updatePoint(p, values);
959: }
960: // Load empty value
961: value_type defValue[fiberDim];
963: for(int i = 0; i < fiberDim; ++i) {
964: fs >> defValue[i];
965: }
966: section.setDefault(defValue);
967: // Load and send remote
968: for(int pr = 1; pr < section.commSize(); ++pr) {
969: PetscInt size = section.getChart().size();
970: PetscInt v = 0;
971: value_type *values;
972: value_type emptyValues[fiberDim];
975: MPI_Send(&size, 1, MPIU_INT, pr, 1, section.comm());CHKERRXX(ierr);
976: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
977: for(index_type p = min; p < max; ++p) {
978: for(int i = 0; i < fiberDim; ++i, ++v) {
979: fs >> values[v];
980: }
981: }
982: MPI_Send(values, size*fiberDim, mover.getType(), pr, 1, section.comm());CHKERRXX(ierr);
983: for(int i = 0; i < fiberDim; ++i) {
984: fs >> emptyValues[i];
985: }
986: MPI_Send(emptyValues, fiberDim, mover.getType(), pr, 1, section.comm());CHKERRXX(ierr);
987: }
988: } else {
989: // Load remote
990: PetscInt size;
991: value_type *values;
992: value_type emptyValues[fiberDim];
993: value_type pvalues[fiberDim];
994: MPI_Status status;
995: PetscInt v = 0;
998: assert(sizeof(value_type) <= sizeof(PetscScalar));
999: MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
1000: PetscMalloc(size*fiberDim * sizeof(value_type), &values);CHKERRXX(ierr);
1001: MPI_Recv(values, size*fiberDim, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1002: for(index_type p = min; p < max; ++p) {
1003: for(int i = 0; i < fiberDim; ++i, ++v) {
1004: pvalues[i] = values[v];
1005: }
1006: section.updatePoint(p, pvalues);
1007: }
1008: PetscFree(values);CHKERRXX(ierr);
1009: MPI_Recv(emptyValues, fiberDim, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1010: section.setDefault(emptyValues);
1011: }
1012: }
1013: template<typename Point_, typename Value_>
1014: static void loadSection(std::ifstream& fs, ISection<Point_, Value_>& section) {
1015: typedef typename ISection<Point_, Value_>::point_type point_type;
1016: typedef typename ISection<Point_, Value_>::value_type value_type;
1017: MPIMover<value_type> mover(section.comm());
1018: // Load atlas
1019: loadSection(fs, *section.getAtlas());
1020: section.allocatePoint();
1021: point_type min = section.getChart().min();
1022: point_type max = section.getChart().max();
1023: int maxDim = -1;
1025: if (section.commRank() == 0) {
1026: // Load local values
1027: for(point_type p = min; p < max; ++p) {
1028: maxDim = std::max(maxDim, section.getFiberDimension(p));
1029: }
1030: value_type *values = new value_type[maxDim];
1031: for(point_type p = min; p < max; ++p) {
1032: const int fiberDim = section.getFiberDimension(p);
1034: for(int i = 0; i < fiberDim; ++i) {
1035: fs >> values[i];
1036: }
1037: section.updatePoint(p, values);
1038: }
1039: delete [] values;
1040: // Load and send remote
1041: for(int p = 1; p < section.commSize(); ++p) {
1042: PetscInt size = section.size();
1043: value_type *values;
1046: MPI_Send(&size, 1, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
1047: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1048: for(PetscInt v = 0; v < size; ++v) {
1049: fs >> values[v];
1050: }
1051: MPI_Send(values, size, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
1052: }
1053: } else {
1054: // Load remote
1055: PetscInt size;
1056: value_type *values;
1057: MPI_Status status;
1058: PetscInt maxDim = 0;
1059: PetscInt off = 0;
1062: MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
1063: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1064: MPI_Recv(values, size, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1065: for(point_type p = min; p < max; ++p) {
1066: maxDim = std::max(maxDim, section.getFiberDimension(p));
1067: }
1068: value_type *pvalues = new value_type[maxDim];
1069: for(point_type p = min; p < max; ++p) {
1070: const int fiberDim = section.getFiberDimension(p);
1072: for(int i = 0; i < fiberDim; ++i, ++off) {
1073: pvalues[i] = values[off];
1074: }
1075: section.updatePoint(p, pvalues);
1076: }
1077: delete [] pvalues;
1078: PetscFree(values);CHKERRXX(ierr);
1079: }
1080: }
1081: template<typename Point_, typename Value_>
1082: static void loadSection(std::ifstream& fs, IGeneralSection<Point_, Value_>& section) {
1083: typedef typename IGeneralSection<Point_, Value_>::point_type point_type;
1084: typedef typename IGeneralSection<Point_, Value_>::value_type value_type;
1085: MPIMover<value_type> mover(section.comm());
1086: // Load atlas
1087: loadSection(fs, *section.getAtlas());
1088: section.allocatePoint();
1089: point_type min = section.getChart().min();
1090: point_type max = section.getChart().max();
1091: int maxDim = -1;
1093: if (section.commRank() == 0) {
1094: // Load local values
1095: for(point_type p = min; p < max; ++p) {
1096: maxDim = std::max(maxDim, section.getFiberDimension(p));
1097: }
1098: value_type *values = new value_type[maxDim];
1099: for(point_type p = min; p < max; ++p) {
1100: const int fiberDim = section.getFiberDimension(p);
1102: for(int i = 0; i < fiberDim; ++i) {
1103: fs >> values[i];
1104: }
1105: section.updatePoint(p, values);
1106: }
1107: delete [] values;
1108: // Load and send remote
1109: for(int p = 1; p < section.commSize(); ++p) {
1110: PetscInt size = section.sizeWithBC();
1111: value_type *values;
1114: MPI_Send(&size, 1, MPIU_INT, p, 1, section.comm());CHKERRXX(ierr);
1115: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1116: for(PetscInt v = 0; v < size; ++v) {
1117: fs >> values[v];
1118: }
1119: MPI_Send(values, size, mover.getType(), p, 1, section.comm());CHKERRXX(ierr);
1120: }
1121: } else {
1122: // Load remote
1123: PetscInt size;
1124: value_type *values;
1125: MPI_Status status;
1126: PetscInt maxDim = 0;
1127: PetscInt off = 0;
1130: assert(sizeof(value_type) <= sizeof(PetscScalar));
1131: MPI_Recv(&size, 1, MPIU_INT, 0, 1, section.comm(), &status);CHKERRXX(ierr);
1132: PetscMalloc(size * sizeof(value_type), &values);CHKERRXX(ierr);
1133: MPI_Recv(values, size, mover.getType(), 0, 1, section.comm(), &status);CHKERRXX(ierr);
1134: for(point_type p = min; p < max; ++p) {
1135: maxDim = std::max(maxDim, section.getFiberDimension(p));
1136: }
1137: value_type *pvalues = new value_type[maxDim];
1138: for(point_type p = min; p < max; ++p) {
1139: const int fiberDim = section.getFiberDimension(p);
1141: for(int i = 0; i < fiberDim; ++i, ++off) {
1142: pvalues[i] = values[off];
1143: }
1144: section.updatePoint(p, pvalues);
1145: }
1146: delete [] pvalues;
1147: PetscFree(values);CHKERRXX(ierr);
1148: }
1149: // Load BC
1150: loadSection(fs, *section.getBC());
1151: // Load spaces
1152: // std::vector<Obj<atlas_type> > _spaces;
1153: // std::vector<Obj<bc_type> > _bcs;
1154: }
1155: };
1156: }
1158: #endif