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