Actual source code: Field.hh
petsc-3.4.5 2014-06-29
1: #ifndef included_ALE_Field_hh
2: #define included_ALE_Field_hh
4: #ifndef included_ALE_SieveAlgorithms_hh
5: #include <sieve/SieveAlgorithms.hh>
6: #endif
8: extern "C" PetscMPIInt DMMesh_DelTag(MPI_Comm comm,PetscMPIInt keyval,void* attr_val,void* extra_state);
10: // Sieve need point_type
11: // Section need point_type and value_type
12: // size(), restrict(), update() need orientation which is a Bundle (circular)
13: // Bundle is Sieve+Section
15: // An AbstractSection is a mapping from Sieve points to sets of values
16: // This is our presentation of a section of a fibre bundle,
17: // in which the Topology is the base space, and
18: // the value sets are vectors in the fiber spaces
19: // The main interface to values is through restrict() and update()
20: // This retrieval uses Sieve closure()
21: // We should call these rawRestrict/rawUpdate
22: // The Section must also be able to set/report the dimension of each fiber
23: // for which we use another section we call an \emph{Atlas}
24: // For some storage schemes, we also want offsets to go with these dimensions
25: // We must have some way of limiting the points associated with values
26: // so each section must support a getPatch() call returning the points with associated fibers
27: // I was using the Chart for this
28: // The Section must be able to participate in \emph{completion}
29: // This means restrict to a provided overlap, and exchange in the restricted sections
30: // Completion does not use hierarchy, so we see the Topology as a DiscreteTopology
31: namespace ALE {
32: template<typename Point_, typename Alloc_ = malloc_allocator<Point_> >
33: class DiscreteSieve {
34: public:
35: typedef Point_ point_type;
36: typedef Alloc_ alloc_type;
37: typedef std::vector<point_type, alloc_type> coneSequence;
38: typedef std::vector<point_type, alloc_type> coneSet;
39: typedef std::vector<point_type, alloc_type> coneArray;
40: typedef std::vector<point_type, alloc_type> supportSequence;
41: typedef std::vector<point_type, alloc_type> supportSet;
42: typedef std::vector<point_type, alloc_type> supportArray;
43: typedef std::set<point_type, std::less<point_type>, alloc_type> points_type;
44: typedef points_type baseSequence;
45: typedef points_type capSequence;
46: typedef typename alloc_type::template rebind<points_type>::other points_alloc_type;
47: typedef typename points_alloc_type::pointer points_ptr;
48: typedef typename alloc_type::template rebind<coneSequence>::other coneSequence_alloc_type;
49: typedef typename coneSequence_alloc_type::pointer coneSequence_ptr;
50: protected:
51: Obj<points_type> _points;
52: Obj<coneSequence> _empty;
53: Obj<coneSequence> _return;
54: alloc_type _allocator;
55: void _init() {
56: points_ptr pPoints = points_alloc_type(this->_allocator).allocate(1);
57: points_alloc_type(this->_allocator).construct(pPoints, points_type());
58: this->_points = Obj<points_type>(pPoints, sizeof(points_type));
59: ///this->_points = new points_type();
60: coneSequence_ptr pEmpty = coneSequence_alloc_type(this->_allocator).allocate(1);
61: coneSequence_alloc_type(this->_allocator).construct(pEmpty, coneSequence());
62: this->_empty = Obj<coneSequence>(pEmpty, sizeof(coneSequence));
63: ///this->_empty = new coneSequence();
64: coneSequence_ptr pReturn = coneSequence_alloc_type(this->_allocator).allocate(1);
65: coneSequence_alloc_type(this->_allocator).construct(pReturn, coneSequence());
66: this->_return = Obj<coneSequence>(pReturn, sizeof(coneSequence));
67: ///this->_return = new coneSequence();
68: };
69: public:
70: DiscreteSieve() {
71: this->_init();
72: };
73: template<typename Input>
74: DiscreteSieve(const Obj<Input>& points) {
75: this->_init();
76: this->_points->insert(points->begin(), points->end());
77: }
78: virtual ~DiscreteSieve() {};
79: public:
80: void addPoint(const point_type& point) {
81: this->_points->insert(point);
82: }
83: template<typename Input>
84: void addPoints(const Obj<Input>& points) {
85: this->_points->insert(points->begin(), points->end());
86: }
87: const Obj<coneSequence>& cone(const point_type& p) {return this->_empty;}
88: template<typename Input>
89: const Obj<coneSequence>& cone(const Input& p) {return this->_empty;}
90: const Obj<coneSet>& nCone(const point_type& p, const int level) {
91: if (level == 0) {
92: return this->closure(p);
93: } else {
94: return this->_empty;
95: }
96: }
97: const Obj<coneArray>& closure(const point_type& p) {
98: this->_return->clear();
99: this->_return->push_back(p);
100: return this->_return;
101: }
102: const Obj<supportSequence>& support(const point_type& p) {return this->_empty;}
103: template<typename Input>
104: const Obj<supportSequence>& support(const Input& p) {return this->_empty;}
105: const Obj<supportSet>& nSupport(const point_type& p, const int level) {
106: if (level == 0) {
107: return this->star(p);
108: } else {
109: return this->_empty;
110: }
111: }
112: const Obj<supportArray>& star(const point_type& p) {
113: this->_return->clear();
114: this->_return->push_back(p);
115: return this->_return;
116: }
117: const Obj<capSequence>& roots() {return this->_points;}
118: const Obj<capSequence>& cap() {return this->_points;}
119: const Obj<baseSequence>& leaves() {return this->_points;}
120: const Obj<baseSequence>& base() {return this->_points;}
121: template<typename Color>
122: void addArrow(const point_type& p, const point_type& q, const Color& color) {
123: throw ALE::Exception("Cannot add an arrow to a DiscreteSieve");
124: }
125: void stratify() {};
126: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
127: ostringstream txt;
128: int rank;
130: if (comm == MPI_COMM_NULL) {
131: comm = MPI_COMM_SELF;
132: rank = 0;
133: } else {
134: MPI_Comm_rank(comm, &rank);
135: }
136: if (name == "") {
137: if(rank == 0) {
138: txt << "viewing a DiscreteSieve" << std::endl;
139: }
140: } else {
141: if(rank == 0) {
142: txt << "viewing DiscreteSieve '" << name << "'" << std::endl;
143: }
144: }
145: for(typename points_type::const_iterator p_iter = this->_points->begin(); p_iter != this->_points->end(); ++p_iter) {
146: txt << " Point " << *p_iter << std::endl;
147: }
148: PetscSynchronizedPrintf(comm, txt.str().c_str());
149: PetscSynchronizedFlush(comm);
150: };
151: };
152: // A ConstantSection is the simplest Section
153: // All fibers are dimension 1
154: // All values are equal to a constant
155: // We need no value storage and no communication for completion
156: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
157: class ConstantSection : public ALE::ParallelObject {
158: public:
159: typedef Point_ point_type;
160: typedef Value_ value_type;
161: typedef Alloc_ alloc_type;
162: typedef std::set<point_type, std::less<point_type>, alloc_type> chart_type;
163: protected:
164: chart_type _chart;
165: value_type _value[2];
166: public:
167: ConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
168: _value[1] = 0;
169: };
170: ConstantSection(MPI_Comm comm, const value_type& value, const int debug) : ParallelObject(comm, debug) {
171: _value[0] = value;
172: _value[1] = value;
173: };
174: ConstantSection(MPI_Comm comm, const value_type& value, const value_type& defaultValue, const int debug) : ParallelObject(comm, debug) {
175: _value[0] = value;
176: _value[1] = defaultValue;
177: };
178: public: // Verifiers
179: void checkPoint(const point_type& point) const {
180: if (this->_chart.find(point) == this->_chart.end()) {
181: ostringstream msg;
182: msg << "Invalid section point " << point << std::endl;
183: throw ALE::Exception(msg.str().c_str());
184: }
185: };
186: void checkDimension(const int& dim) {
187: if (dim != 1) {
188: ostringstream msg;
189: msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
190: throw ALE::Exception(msg.str().c_str());
191: }
192: };
193: bool hasPoint(const point_type& point) const {
194: return this->_chart.count(point) > 0;
195: };
196: public: // Accessors
197: const chart_type& getChart() {return this->_chart;};
198: void addPoint(const point_type& point) {
199: this->_chart.insert(point);
200: };
201: template<typename Points>
202: void addPoint(const Obj<Points>& points) {
203: this->_chart.insert(points->begin(), points->end());
204: }
205: template<typename Points>
206: void addPoint(const Points& points) {
207: this->_chart.insert(points.begin(), points.end());
208: }
209: // void addPoint(const std::set<point_type>& points) {
210: // this->_chart.insert(points.begin(), points.end());
211: // };
212: value_type getDefaultValue() {return this->_value[1];};
213: void setDefaultValue(const value_type value) {this->_value[1] = value;};
214: void copy(const Obj<ConstantSection>& section) {
215: const chart_type& chart = section->getChart();
217: this->addPoint(chart);
218: this->_value[0] = section->restrictSpace()[0];
219: this->setDefaultValue(section->getDefaultValue());
220: };
221: public: // Sizes
222: void clear() {
223: this->_chart.clear();
224: };
225: int getFiberDimension(const point_type& p) const {
226: if (this->hasPoint(p)) return 1;
227: return 0;
228: };
229: void setFiberDimension(const point_type& p, int dim) {
230: this->checkDimension(dim);
231: this->addPoint(p);
232: };
233: template<typename Sequence>
234: void setFiberDimension(const Obj<Sequence>& points, int dim) {
235: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
236: this->setFiberDimension(*p_iter, dim);
237: }
238: }
239: void addFiberDimension(const point_type& p, int dim) {
240: if (this->_chart.find(p) != this->_chart.end()) {
241: ostringstream msg;
242: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed 1" << std::endl;
243: throw ALE::Exception(msg.str().c_str());
244: } else {
245: this->setFiberDimension(p, dim);
246: }
247: }
248: int size(const point_type& p) {return this->getFiberDimension(p);};
249: void allocatePoint() {};
250: public: // Restriction
251: const value_type *restrictSpace() const {
252: return this->_value;
253: };
254: const value_type *restrictPoint(const point_type& p) const {
255: if (this->hasPoint(p)) {
256: return this->_value;
257: }
258: return &this->_value[1];
259: };
260: void updatePoint(const point_type& p, const value_type v[]) {
261: this->_value[0] = v[0];
262: };
263: void updateAddPoint(const point_type& p, const value_type v[]) {
264: this->_value[0] += v[0];
265: };
266: public:
267: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
268: ostringstream txt;
269: int rank;
271: if (comm == MPI_COMM_NULL) {
272: comm = this->comm();
273: rank = this->commRank();
274: } else {
275: MPI_Comm_rank(comm, &rank);
276: }
277: if (name == "") {
278: if(rank == 0) {
279: txt << "viewing a ConstantSection" << std::endl;
280: }
281: } else {
282: if(rank == 0) {
283: txt << "viewing ConstantSection '" << name << "'" << std::endl;
284: }
285: }
286: txt <<"["<<this->commRank()<<"]: Value " << this->_value[0] << " Default Value " << this->_value[1] << std::endl;
287: PetscSynchronizedPrintf(comm, txt.str().c_str());
288: PetscSynchronizedFlush(comm);
289: };
290: };
291: // A UniformSection often acts as an Atlas
292: // All fibers are the same dimension
293: // Note we can use a ConstantSection for this Atlas
294: // Each point may have a different vector
295: // Thus we need storage for values, and hence must implement completion
296: template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
297: class UniformSection : public ALE::ParallelObject {
298: public:
299: typedef Point_ point_type;
300: typedef Value_ value_type;
301: typedef Alloc_ alloc_type;
302: typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
303: typedef ConstantSection<point_type, int, point_alloc_type> atlas_type;
304: typedef typename atlas_type::chart_type chart_type;
305: typedef struct {value_type v[fiberDim];} fiber_type;
306: typedef typename alloc_type::template rebind<std::pair<const point_type, fiber_type> >::other pair_alloc_type;
307: typedef std::map<point_type, fiber_type, std::less<point_type>, pair_alloc_type> values_type;
308: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
309: typedef typename atlas_alloc_type::pointer atlas_ptr;
310: protected:
311: size_t _contiguous_array_size;
312: value_type *_contiguous_array;
313: Obj<atlas_type> _atlas;
314: values_type _array;
315: fiber_type _emptyValue;
316: alloc_type _allocator;
317: public:
318: UniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _contiguous_array_size(0), _contiguous_array(NULL) {
319: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
320: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
321: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
322: int dim = fiberDim;
323: this->_atlas->updatePoint(*this->_atlas->getChart().begin(), &dim);
324: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
325: };
326: UniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _contiguous_array_size(0), _contiguous_array(NULL), _atlas(atlas) {
327: int dim = fiberDim;
328: this->_atlas->update(*this->_atlas->getChart().begin(), &dim);
329: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
330: };
331: virtual ~UniformSection() {
332: this->_atlas = NULL;
333: if (this->_contiguous_array) {
334: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
335: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
336: }
337: };
338: public:
339: value_type *getRawArray(const int size) {
340: static value_type *array = NULL;
341: static int maxSize = 0;
343: if (size > maxSize) {
344: const value_type dummy(0);
346: if (array) {
347: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
348: this->_allocator.deallocate(array, maxSize);
349: ///delete [] array;
350: }
351: maxSize = size;
352: array = this->_allocator.allocate(maxSize);
353: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
354: ///array = new value_type[maxSize];
355: };
356: return array;
357: };
358: public: // Verifiers
359: bool hasPoint(const point_type& point) {
360: return this->_atlas->hasPoint(point);
361: };
362: void checkDimension(const int& dim) {
363: if (dim != fiberDim) {
364: ostringstream msg;
365: msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
366: throw ALE::Exception(msg.str().c_str());
367: }
368: };
369: public: // Accessors
370: const chart_type& getChart() {return this->_atlas->getChart();};
371: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
372: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
373: void addPoint(const point_type& point) {
374: this->setFiberDimension(point, fiberDim);
375: };
376: template<typename Points>
377: void addPoint(const Obj<Points>& points) {
378: for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
379: this->setFiberDimension(*p_iter, fiberDim);
380: }
381: }
382: void copy(const Obj<UniformSection>& section) {
383: this->getAtlas()->copy(section->getAtlas());
384: const chart_type& chart = section->getChart();
386: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
387: this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
388: }
389: }
390: public: // Sizes
391: void clear() {
392: this->_array.clear();
393: this->_atlas->clear();
394: }
395: int getFiberDimension(const point_type& p) const {
396: return this->_atlas->restrictPoint(p)[0];
397: }
398: void setFiberDimension(const point_type& p, int dim) {
399: this->update();
400: this->checkDimension(dim);
401: this->_atlas->addPoint(p);
402: this->_atlas->updatePoint(p, &dim);
403: }
404: template<typename Sequence>
405: void setFiberDimension(const Obj<Sequence>& points, int dim) {
406: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
407: this->setFiberDimension(*p_iter, dim);
408: }
409: }
410: void setFiberDimension(const std::set<point_type>& points, int dim) {
411: for(typename std::set<point_type>::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
412: this->setFiberDimension(*p_iter, dim);
413: }
414: };
415: void addFiberDimension(const point_type& p, int dim) {
416: if (this->hasPoint(p)) {
417: ostringstream msg;
418: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
419: throw ALE::Exception(msg.str().c_str());
420: } else {
421: this->setFiberDimension(p, dim);
422: }
423: };
424: int size() const {
425: const chart_type& points = this->_atlas->getChart();
426: int size = 0;
428: for(typename chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
429: size += this->getFiberDimension(*p_iter);
430: }
431: return size;
432: };
433: int sizeWithBC() const {
434: return this->size();
435: };
436: void allocatePoint() {};
437: public: // Restriction
438: const value_type *restrictSpace() {
439: const chart_type& chart = this->getChart();
440: const value_type dummy = 0;
441: int k = 0;
443: if (chart.size() > this->_contiguous_array_size*fiberDim) {
444: if (this->_contiguous_array) {
445: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
446: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
447: }
448: this->_contiguous_array_size = chart.size()*fiberDim;
449: this->_contiguous_array = this->_allocator.allocate(this->_contiguous_array_size*fiberDim);
450: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.construct(this->_contiguous_array+i, dummy);}
451: }
452: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
453: const value_type *a = this->_array[*p_iter].v;
455: for(int i = 0; i < fiberDim; ++i, ++k) {
456: this->_contiguous_array[k] = a[i];
457: }
458: }
459: return this->_contiguous_array;
460: };
461: void update() {
462: if (this->_contiguous_array) {
463: const chart_type& chart = this->getChart();
464: int k = 0;
466: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
467: value_type *a = this->_array[*p_iter].v;
469: for(int i = 0; i < fiberDim; ++i, ++k) {
470: a[i] = this->_contiguous_array[k];
471: }
472: }
473: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
474: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
475: this->_contiguous_array_size = 0;
476: this->_contiguous_array = NULL;
477: }
478: };
479: // Return only the values associated to this point, not its closure
480: const value_type *restrictPoint(const point_type& p) {
481: if (this->_array.find(p) == this->_array.end()) return this->_emptyValue.v;
482: this->update();
483: return this->_array[p].v;
484: };
485: // Update only the values associated to this point, not its closure
486: void updatePoint(const point_type& p, const value_type v[]) {
487: this->update();
488: for(int i = 0; i < fiberDim; ++i) {
489: this->_array[p].v[i] = v[i];
490: }
491: };
492: // Update only the values associated to this point, not its closure
493: void updateAddPoint(const point_type& p, const value_type v[]) {
494: this->update();
495: for(int i = 0; i < fiberDim; ++i) {
496: this->_array[p].v[i] += v[i];
497: }
498: };
499: void updatePointAll(const point_type& p, const value_type v[]) {
500: this->updatePoint(p, v);
501: };
502: public:
503: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
504: ostringstream txt;
505: int rank;
507: this->update();
508: if (comm == MPI_COMM_NULL) {
509: comm = this->comm();
510: rank = this->commRank();
511: } else {
512: MPI_Comm_rank(comm, &rank);
513: }
514: if (name == "") {
515: if(rank == 0) {
516: txt << "viewing a UniformSection" << std::endl;
517: }
518: } else {
519: if(rank == 0) {
520: txt << "viewing UniformSection '" << name << "'" << std::endl;
521: }
522: }
523: const typename atlas_type::chart_type& chart = this->_atlas->getChart();
524: values_type& array = this->_array;
526: for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
527: const point_type& point = *p_iter;
528: const typename atlas_type::value_type dim = this->_atlas->restrictPoint(point)[0];
530: if (dim != 0) {
531: txt << "[" << this->commRank() << "]: " << point << " dim " << dim << " ";
532: for(int i = 0; i < dim; i++) {
533: txt << " " << array[point].v[i];
534: }
535: txt << std::endl;
536: }
537: }
538: if (chart.size() == 0) {
539: txt << "[" << this->commRank() << "]: empty" << std::endl;
540: }
541: PetscSynchronizedPrintf(comm, txt.str().c_str());
542: PetscSynchronizedFlush(comm);
543: };
544: };
546: // A FauxConstantSection is the simplest Section
547: // All fibers are dimension 1
548: // All values are equal to a constant
549: // We need no value storage and no communication for completion
550: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Point_> >
551: class FauxConstantSection : public ALE::ParallelObject {
552: public:
553: typedef Point_ point_type;
554: typedef Value_ value_type;
555: typedef Alloc_ alloc_type;
556: protected:
557: value_type _value;
558: public:
559: FauxConstantSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
560: FauxConstantSection(MPI_Comm comm, const value_type& value, const int debug) : ParallelObject(comm, debug), _value(value) {};
561: public: // Verifiers
562: void checkDimension(const int& dim) {
563: if (dim != 1) {
564: ostringstream msg;
565: msg << "Invalid fiber dimension " << dim << " must be 1" << std::endl;
566: throw ALE::Exception(msg.str().c_str());
567: }
568: };
569: public: // Accessors
570: void addPoint(const point_type& point) {
571: }
572: template<typename Points>
573: void addPoint(const Obj<Points>& points) {
574: }
575: template<typename Points>
576: void addPoint(const Points& points) {
577: }
578: void copy(const Obj<FauxConstantSection>& section) {
579: this->_value = section->restrictPoint(point_type())[0];
580: }
581: public: // Sizes
582: void clear() {
583: };
584: int getFiberDimension(const point_type& p) const {
585: return 1;
586: };
587: void setFiberDimension(const point_type& p, int dim) {
588: this->checkDimension(dim);
589: this->addPoint(p);
590: };
591: template<typename Sequence>
592: void setFiberDimension(const Obj<Sequence>& points, int dim) {
593: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
594: this->setFiberDimension(*p_iter, dim);
595: }
596: }
597: void addFiberDimension(const point_type& p, int dim) {
598: this->setFiberDimension(p, dim);
599: }
600: int size(const point_type& p) {return this->getFiberDimension(p);}
601: public: // Restriction
602: const value_type *restrictPoint(const point_type& p) const {
603: return &this->_value;
604: };
605: void updatePoint(const point_type& p, const value_type v[]) {
606: this->_value = v[0];
607: };
608: void updateAddPoint(const point_type& p, const value_type v[]) {
609: this->_value += v[0];
610: };
611: public:
612: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
613: ostringstream txt;
614: int rank;
616: if (comm == MPI_COMM_NULL) {
617: comm = this->comm();
618: rank = this->commRank();
619: } else {
620: MPI_Comm_rank(comm, &rank);
621: }
622: if (name == "") {
623: if(rank == 0) {
624: txt << "viewing a FauxConstantSection" << std::endl;
625: }
626: } else {
627: if(rank == 0) {
628: txt << "viewing FauxConstantSection '" << name << "'" << std::endl;
629: }
630: }
631: txt <<"["<<this->commRank()<<"]: Value " << this->_value << std::endl;
632: PetscSynchronizedPrintf(comm, txt.str().c_str());
633: PetscSynchronizedFlush(comm);
634: };
635: };
636: // Make a simple set from the keys of a map
637: template<typename Map>
638: class SetFromMap {
639: public:
640: typedef typename Map::size_type size_type;
641: class const_iterator {
642: public:
643: typedef typename Map::key_type value_type;
644: typedef typename Map::size_type size_type;
645: protected:
646: typename Map::const_iterator _iter;
647: public:
648: const_iterator(const typename Map::const_iterator& iter): _iter(iter) {};
649: ~const_iterator() {};
650: public:
651: const_iterator& operator=(const const_iterator& iter) {this->_iter = iter._iter; return this->_iter;};
652: bool operator==(const const_iterator& iter) const {return this->_iter == iter._iter;};
653: bool operator!=(const const_iterator& iter) const {return this->_iter != iter._iter;};
654: const_iterator& operator++() {++this->_iter; return *this;}
655: const_iterator operator++(int) {
656: const_iterator tmp(*this);
657: ++(*this);
658: return tmp;
659: };
660: const_iterator& operator--() {--this->_iter; return *this;}
661: const_iterator operator--(int) {
662: const_iterator tmp(*this);
663: --(*this);
664: return tmp;
665: };
666: value_type operator*() const {return this->_iter->first;};
667: };
668: protected:
669: const Map& _map;
670: public:
671: SetFromMap(const Map& map): _map(map) {};
672: public:
673: const_iterator begin() const {return const_iterator(this->_map.begin());};
674: const_iterator end() const {return const_iterator(this->_map.end());};
675: size_type size() const {return this->_map.size();};
676: };
677: // A NewUniformSection often acts as an Atlas
678: // All fibers are the same dimension
679: // Note we can use a ConstantSection for this Atlas
680: // Each point may have a different vector
681: // Thus we need storage for values, and hence must implement completion
682: template<typename Point_, typename Value_, int fiberDim = 1, typename Alloc_ = malloc_allocator<Value_> >
683: class NewUniformSection : public ALE::ParallelObject {
684: public:
685: typedef Point_ point_type;
686: typedef Value_ value_type;
687: typedef Alloc_ alloc_type;
688: typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
689: typedef FauxConstantSection<point_type, int, point_alloc_type> atlas_type;
690: typedef struct {value_type v[fiberDim];} fiber_type;
691: typedef typename alloc_type::template rebind<std::pair<const point_type, fiber_type> >::other pair_alloc_type;
692: typedef std::map<point_type, fiber_type, std::less<point_type>, pair_alloc_type> values_type;
693: typedef SetFromMap<values_type> chart_type;
694: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
695: typedef typename atlas_alloc_type::pointer atlas_ptr;
696: protected:
697: values_type _array;
698: chart_type _chart;
699: size_t _contiguous_array_size;
700: value_type *_contiguous_array;
701: Obj<atlas_type> _atlas;
702: fiber_type _emptyValue;
703: alloc_type _allocator;
704: public:
705: NewUniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _chart(_array), _contiguous_array_size(0), _contiguous_array(NULL) {
706: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
707: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
708: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
709: int dim = fiberDim;
710: this->_atlas->update(point_type(), &dim);
711: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
712: };
713: NewUniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _chart(_array), _contiguous_array_size(0), _contiguous_array(NULL), _atlas(atlas) {
714: int dim = fiberDim;
715: this->_atlas->update(point_type(), &dim);
716: for(int i = 0; i < fiberDim; ++i) this->_emptyValue.v[i] = value_type();
717: };
718: ~NewUniformSection() {
719: this->_atlas = NULL;
720: if (this->_contiguous_array) {
721: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
722: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
723: }
724: };
725: public:
726: value_type *getRawArray(const int size) {
727: static value_type *array = NULL;
728: static int maxSize = 0;
730: if (size > maxSize) {
731: const value_type dummy(0);
733: if (array) {
734: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
735: this->_allocator.deallocate(array, maxSize);
736: ///delete [] array;
737: }
738: maxSize = size;
739: array = this->_allocator.allocate(maxSize);
740: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
741: ///array = new value_type[maxSize];
742: };
743: return array;
744: };
745: public: // Verifiers
746: bool hasPoint(const point_type& point) {
747: return (this->_array.find(point) != this->_array.end());
748: ///return this->_atlas->hasPoint(point);
749: };
750: void checkDimension(const int& dim) {
751: if (dim != fiberDim) {
752: ostringstream msg;
753: msg << "Invalid fiber dimension " << dim << " must be " << fiberDim << std::endl;
754: throw ALE::Exception(msg.str().c_str());
755: }
756: };
757: public: // Accessors
758: const chart_type& getChart() {return this->_chart;}
759: const Obj<atlas_type>& getAtlas() {return this->_atlas;}
760: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;}
761: void addPoint(const point_type& point) {
762: this->setFiberDimension(point, fiberDim);
763: }
764: template<typename Points>
765: void addPoint(const Obj<Points>& points) {
766: for(typename Points::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
767: this->setFiberDimension(*p_iter, fiberDim);
768: }
769: }
770: void copy(const Obj<NewUniformSection>& section) {
771: this->getAtlas()->copy(section->getAtlas());
772: const chart_type& chart = section->getChart();
774: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
775: this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
776: }
777: }
778: public: // Sizes
779: void clear() {
780: this->_array.clear();
781: this->_atlas->clear();
782: };
783: int getFiberDimension(const point_type& p) const {
784: return fiberDim;
785: };
786: void setFiberDimension(const point_type& p, int dim) {
787: this->update();
788: this->checkDimension(dim);
789: this->_atlas->addPoint(p);
790: this->_atlas->updatePoint(p, &dim);
791: this->_array[p] = fiber_type();
792: };
793: template<typename Sequence>
794: void setFiberDimension(const Obj<Sequence>& points, int dim) {
795: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
796: this->setFiberDimension(*p_iter, dim);
797: }
798: }
799: void setFiberDimension(const std::set<point_type>& points, int dim) {
800: for(typename std::set<point_type>::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
801: this->setFiberDimension(*p_iter, dim);
802: }
803: };
804: void addFiberDimension(const point_type& p, int dim) {
805: if (this->hasPoint(p)) {
806: ostringstream msg;
807: msg << "Invalid addition to fiber dimension " << dim << " cannot exceed " << fiberDim << std::endl;
808: throw ALE::Exception(msg.str().c_str());
809: } else {
810: this->setFiberDimension(p, dim);
811: }
812: };
813: int size() {
814: const chart_type& points = this->getChart();
815: int size = 0;
817: for(typename chart_type::iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
818: size += this->getFiberDimension(*p_iter);
819: }
820: return size;
821: };
822: int sizeWithBC() {
823: return this->size();
824: };
825: void allocatePoint() {};
826: public: // Restriction
827: // Return a pointer to the entire contiguous storage array
828: const value_type *restrictSpace() {
829: const chart_type& chart = this->getChart();
830: const value_type dummy = 0;
831: int k = 0;
833: if (chart.size() > this->_contiguous_array_size*fiberDim) {
834: if (this->_contiguous_array) {
835: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
836: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
837: }
838: this->_contiguous_array_size = chart.size()*fiberDim;
839: this->_contiguous_array = this->_allocator.allocate(this->_contiguous_array_size*fiberDim);
840: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.construct(this->_contiguous_array+i, dummy);}
841: }
842: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
843: const value_type *a = this->_array[*p_iter].v;
845: for(int i = 0; i < fiberDim; ++i, ++k) {
846: this->_contiguous_array[k] = a[i];
847: }
848: }
849: return this->_contiguous_array;
850: };
851: void update() {
852: if (this->_contiguous_array) {
853: const chart_type& chart = this->getChart();
854: int k = 0;
856: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
857: value_type *a = this->_array[*p_iter].v;
859: for(int i = 0; i < fiberDim; ++i, ++k) {
860: a[i] = this->_contiguous_array[k];
861: }
862: }
863: for(size_t i = 0; i < this->_contiguous_array_size; ++i) {this->_allocator.destroy(this->_contiguous_array+i);}
864: this->_allocator.deallocate(this->_contiguous_array, this->_contiguous_array_size);
865: this->_contiguous_array_size = 0;
866: this->_contiguous_array = NULL;
867: }
868: };
869: // Return only the values associated to this point, not its closure
870: const value_type *restrictPoint(const point_type& p) {
871: if (this->_array.find(p) == this->_array.end()) return this->_emptyValue.v;
872: this->update();
873: return this->_array[p].v;
874: };
875: // Update only the values associated to this point, not its closure
876: void updatePoint(const point_type& p, const value_type v[]) {
877: this->update();
878: for(int i = 0; i < fiberDim; ++i) {
879: this->_array[p].v[i] = v[i];
880: }
881: };
882: // Update only the values associated to this point, not its closure
883: void updateAddPoint(const point_type& p, const value_type v[]) {
884: this->update();
885: for(int i = 0; i < fiberDim; ++i) {
886: this->_array[p].v[i] += v[i];
887: }
888: };
889: void updatePointAll(const point_type& p, const value_type v[]) {
890: this->updatePoint(p, v);
891: };
892: public:
893: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
894: ostringstream txt;
895: int rank;
897: this->update();
898: if (comm == MPI_COMM_NULL) {
899: comm = this->comm();
900: rank = this->commRank();
901: } else {
902: MPI_Comm_rank(comm, &rank);
903: }
904: if (name == "") {
905: if(rank == 0) {
906: txt << "viewing a NewUniformSection" << std::endl;
907: }
908: } else {
909: if(rank == 0) {
910: txt << "viewing NewUniformSection '" << name << "'" << std::endl;
911: }
912: }
913: const chart_type& chart = this->getChart();
914: values_type& array = this->_array;
916: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
917: const point_type& point = *p_iter;
919: if (fiberDim != 0) {
920: txt << "[" << this->commRank() << "]: " << point << " dim " << fiberDim << " ";
921: for(int i = 0; i < fiberDim; i++) {
922: txt << " " << array[point].v[i];
923: }
924: txt << std::endl;
925: }
926: }
927: if (chart.size() == 0) {
928: txt << "[" << this->commRank() << "]: empty" << std::endl;
929: }
930: PetscSynchronizedPrintf(comm, txt.str().c_str());
931: PetscSynchronizedFlush(comm);
932: };
933: };
934: // A Section is our most general construct (more general ones could be envisioned)
935: // The Atlas is a UniformSection of dimension 1 and value type Point
936: // to hold each fiber dimension and offsets into a contiguous patch array
937: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
938: typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other> >
939: class Section : public ALE::ParallelObject {
940: public:
941: typedef Point_ point_type;
942: typedef Value_ value_type;
943: typedef Alloc_ alloc_type;
944: typedef Atlas_ atlas_type;
945: typedef Point index_type;
946: typedef typename atlas_type::chart_type chart_type;
947: typedef value_type * values_type;
948: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
949: typedef typename atlas_alloc_type::pointer atlas_ptr;
950: protected:
951: Obj<atlas_type> _atlas;
952: Obj<atlas_type> _atlasNew;
953: values_type _array;
954: alloc_type _allocator;
955: public:
956: Section(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
957: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
958: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
959: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
960: this->_atlasNew = NULL;
961: this->_array = NULL;
962: };
963: Section(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL), _array(NULL) {};
964: virtual ~Section() {
965: if (this->_array) {
966: const int totalSize = this->sizeWithBC();
968: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
969: this->_allocator.deallocate(this->_array, totalSize);
970: ///delete [] this->_array;
971: this->_array = NULL;
972: }
973: };
974: public:
975: value_type *getRawArray(const int size) {
976: static value_type *array = NULL;
977: static int maxSize = 0;
979: if (size > maxSize) {
980: const value_type dummy(0);
982: if (array) {
983: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
984: this->_allocator.deallocate(array, maxSize);
985: ///delete [] array;
986: }
987: maxSize = size;
988: array = this->_allocator.allocate(maxSize);
989: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
990: ///array = new value_type[maxSize];
991: };
992: return array;
993: };
994: int getStorageSize() const {
995: return this->sizeWithBC();
996: };
997: public: // Verifiers
998: bool hasPoint(const point_type& point) {
999: return this->_atlas->hasPoint(point);
1000: };
1001: public: // Accessors
1002: const chart_type& getChart() {return this->_atlas->getChart();};
1003: void setChart(chart_type& chart) {};
1004: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
1005: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1006: const Obj<atlas_type>& getNewAtlas() {return this->_atlasNew;};
1007: void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1008: const chart_type& getChart() const {return this->_atlas->getChart();};
1009: public: // BC
1010: // Returns the number of constraints on a point
1011: int getConstraintDimension(const point_type& p) const {
1012: return std::max(0, -this->_atlas->restrictPoint(p)->prefix);
1013: }
1014: // Set the number of constraints on a point
1015: // We only allow the entire point to be constrained, so these will be the
1016: // only dofs on the point
1017: void setConstraintDimension(const point_type& p, const int numConstraints) {
1018: this->setFiberDimension(p, -numConstraints);
1019: };
1020: void addConstraintDimension(const point_type& p, const int numConstraints) {
1021: throw ALE::Exception("Variable constraint dimensions not available for this Section type");
1022: };
1023: void copyBC(const Obj<Section>& section) {
1024: const chart_type& chart = this->getChart();
1026: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1027: if (this->getConstraintDimension(*p_iter)) {
1028: this->updatePointBC(*p_iter, section->restrictPoint(*p_iter));
1029: }
1030: }
1031: };
1032: void defaultConstraintDof() {};
1033: public: // Sizes
1034: void clear() {
1035: const int totalSize = this->sizeWithBC();
1037: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1038: this->_allocator.deallocate(this->_array, totalSize);
1039: ///delete [] this->_array;
1040: this->_array = NULL;
1041: this->_atlas->clear();
1042: };
1043: // Return the total number of dofs on the point (free and constrained)
1044: int getFiberDimension(const point_type& p) const {
1045: return std::abs(this->_atlas->restrictPoint(p)->prefix);
1046: };
1047: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
1048: return std::abs(atlas->restrictPoint(p)->prefix);
1049: };
1050: // Set the total number of dofs on the point (free and constrained)
1051: void setFiberDimension(const point_type& p, int dim) {
1052: const index_type idx(dim, -1);
1053: this->_atlas->addPoint(p);
1054: this->_atlas->updatePoint(p, &idx);
1055: };
1056: template<typename Sequence>
1057: void setFiberDimension(const Obj<Sequence>& points, int dim) {
1058: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1059: this->setFiberDimension(*p_iter, dim);
1060: }
1061: }
1062: void addFiberDimension(const point_type& p, int dim) {
1063: if (this->_atlas->hasPoint(p)) {
1064: const index_type values(dim, 0);
1065: this->_atlas->updateAddPoint(p, &values);
1066: } else {
1067: this->setFiberDimension(p, dim);
1068: }
1069: }
1070: // Return the number of constrained dofs on this point
1071: // If constrained, this is equal to the fiber dimension
1072: // Otherwise, 0
1073: int getConstrainedFiberDimension(const point_type& p) const {
1074: return std::max((index_type::prefix_type) 0, this->_atlas->restrictPoint(p)->prefix);
1075: };
1076: // Return the total number of free dofs
1077: int size() const {
1078: const chart_type& points = this->getChart();
1079: int size = 0;
1081: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1082: size += this->getConstrainedFiberDimension(*p_iter);
1083: }
1084: return size;
1085: };
1086: // Return the total number of dofs (free and constrained)
1087: int sizeWithBC() const {
1088: const chart_type& points = this->getChart();
1089: int size = 0;
1091: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1092: size += this->getFiberDimension(*p_iter);
1093: }
1094: return size;
1095: };
1096: public: // Index retrieval
1097: const typename index_type::index_type& getIndex(const point_type& p) {
1098: return this->_atlas->restrictPoint(p)->index;
1099: };
1100: void setIndex(const point_type& p, const typename index_type::index_type& index) {
1101: ((typename atlas_type::value_type *) this->_atlas->restrictPoint(p))->index = index;
1102: };
1103: void setIndexBC(const point_type& p, const typename index_type::index_type& index) {
1104: this->setIndex(p, index);
1105: };
1106: void getIndices(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1107: this->getIndices(p, this->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1108: };
1109: template<typename Order_>
1110: void getIndices(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1111: this->getIndices(p, order->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1112: }
1113: void getIndices(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1114: const int& dim = this->getFiberDimension(p);
1115: const int& cDim = this->getConstraintDimension(p);
1116: const int end = start + dim;
1118: if (!cDim) {
1119: if (orientation >= 0) {
1120: for(int i = start; i < end; ++i) {
1121: indices[(*indx)++] = i;
1122: }
1123: } else {
1124: for(int i = end-1; i >= start; --i) {
1125: indices[(*indx)++] = i;
1126: }
1127: }
1128: } else {
1129: if (!freeOnly) {
1130: for(int i = start; i < end; ++i) {
1131: indices[(*indx)++] = -1;
1132: }
1133: }
1134: }
1135: };
1136: public: // Allocation
1137: void allocateStorage() {
1138: const int totalSize = this->sizeWithBC();
1139: const value_type dummy(0);
1141: this->_array = this->_allocator.allocate(totalSize);
1142: ///this->_array = new value_type[totalSize];
1143: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1144: ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1145: };
1146: void replaceStorage(value_type *newArray) {
1147: const int totalSize = this->sizeWithBC();
1149: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1150: this->_allocator.deallocate(this->_array, totalSize);
1151: ///delete [] this->_array;
1152: this->_array = newArray;
1153: this->_atlasNew = NULL;
1154: };
1155: void addPoint(const point_type& point, const int dim) {
1156: if (dim == 0) return;
1157: if (this->_atlasNew.isNull()) {
1158: this->_atlasNew = new atlas_type(this->comm(), this->debug());
1159: this->_atlasNew->copy(this->_atlas);
1160: }
1161: const index_type idx(dim, -1);
1162: this->_atlasNew->addPoint(point);
1163: this->_atlasNew->updatePoint(point, &idx);
1164: };
1165: template<typename Sequence>
1166: void addPoints(const Obj<Sequence>& points, const int dim) {
1167: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1168: this->addPoint(*p_iter, dim);
1169: }
1170: }
1171: void orderPoints(const Obj<atlas_type>& atlas){
1172: const chart_type& chart = this->getChart();
1173: int offset = 0;
1174: int bcOffset = this->size();
1176: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1177: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
1178: const int& dim = idx.prefix;
1180: if (dim >= 0) {
1181: idx.index = offset;
1182: atlas->updatePoint(*c_iter, &idx);
1183: offset += dim;
1184: } else {
1185: idx.index = bcOffset;
1186: atlas->updatePoint(*c_iter, &idx);
1187: bcOffset -= dim;
1188: }
1189: }
1190: };
1191: void allocatePoint() {
1192: this->orderPoints(this->_atlas);
1193: this->allocateStorage();
1194: };
1195: public: // Restriction and Update
1196: // Zero entries
1197: void zero() {
1198: memset(this->_array, 0, this->size()* sizeof(value_type));
1199: };
1200: // Return a pointer to the entire contiguous storage array
1201: const value_type *restrictSpace() {
1202: return this->_array;
1203: };
1204: // Update the entire contiguous storage array
1205: void update(const value_type v[]) {
1206: const value_type *array = this->_array;
1207: const int size = this->size();
1209: for(int i = 0; i < size; i++) {
1210: array[i] = v[i];
1211: }
1212: };
1213: // Return the free values on a point
1214: const value_type *restrictPoint(const point_type& p) {
1215: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1216: };
1217: // Update the free values on a point
1218: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1219: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1220: value_type *a = &(this->_array[idx.index]);
1222: if (orientation >= 0) {
1223: for(int i = 0; i < idx.prefix; ++i) {
1224: a[i] = v[i];
1225: }
1226: } else {
1227: const int last = idx.prefix-1;
1229: for(int i = 0; i < idx.prefix; ++i) {
1230: a[i] = v[last-i];
1231: }
1232: }
1233: };
1234: // Update the free values on a point
1235: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1236: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1237: value_type *a = &(this->_array[idx.index]);
1239: if (orientation >= 0) {
1240: for(int i = 0; i < idx.prefix; ++i) {
1241: a[i] += v[i];
1242: }
1243: } else {
1244: const int last = idx.prefix-1;
1246: for(int i = 0; i < idx.prefix; ++i) {
1247: a[i] += v[last-i];
1248: }
1249: }
1250: };
1251: // Update only the constrained dofs on a point
1252: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
1253: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1254: const int dim = -idx.prefix;
1255: value_type *a = &(this->_array[idx.index]);
1257: if (orientation >= 0) {
1258: for(int i = 0; i < dim; ++i) {
1259: a[i] = v[i];
1260: }
1261: } else {
1262: const int last = dim-1;
1264: for(int i = 0; i < dim; ++i) {
1265: a[i] = v[last-i];
1266: }
1267: }
1268: };
1269: // Update all dofs on a point (free and constrained)
1270: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
1271: const index_type& idx = this->_atlas->restrictPoint(p)[0];
1272: const int dim = std::abs(idx.prefix);
1273: value_type *a = &(this->_array[idx.index]);
1275: if (orientation >= 0) {
1276: for(int i = 0; i < dim; ++i) {
1277: a[i] = v[i];
1278: }
1279: } else {
1280: const int last = dim-1;
1282: for(int i = 0; i < dim; ++i) {
1283: a[i] = v[last-i];
1284: }
1285: }
1286: };
1287: public:
1288: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
1289: ostringstream txt;
1290: int rank;
1292: if (comm == MPI_COMM_NULL) {
1293: comm = this->comm();
1294: rank = this->commRank();
1295: } else {
1296: MPI_Comm_rank(comm, &rank);
1297: }
1298: if(rank == 0) {
1299: txt << "viewing Section " << this->getName() << std::endl;
1300: if (name != "") {
1301: txt << ": " << name << "'";
1302: }
1303: txt << std::endl;
1304: }
1305: const typename atlas_type::chart_type& chart = this->_atlas->getChart();
1306: const value_type *array = this->_array;
1308: for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1309: const point_type& point = *p_iter;
1310: const index_type& idx = this->_atlas->restrictPoint(point)[0];
1311: const int dim = this->getFiberDimension(point);
1313: if (idx.prefix != 0) {
1314: txt << "[" << this->commRank() << "]: " << point << " dim " << idx.prefix << " offset " << idx.index << " ";
1315: for(int i = 0; i < dim; i++) {
1316: txt << " " << array[idx.index+i];
1317: }
1318: txt << std::endl;
1319: }
1320: }
1321: if (chart.size() == 0) {
1322: txt << "[" << this->commRank() << "]: empty" << std::endl;
1323: }
1324: PetscSynchronizedPrintf(comm, txt.str().c_str());
1325: PetscSynchronizedFlush(comm);
1326: };
1327: public: // Fibrations
1328: void setConstraintDof(const point_type& p, const int dofs[]) {};
1329: int getNumSpaces() const {return 1;};
1330: void addSpace() {};
1331: int getFiberDimension(const point_type& p, const int space) {return this->getFiberDimension(p);};
1332: void setFiberDimension(const point_type& p, int dim, const int space) {this->setFiberDimension(p, dim);};
1333: template<typename Sequence>
1334: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
1335: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1336: this->setFiberDimension(*p_iter, dim, space);
1337: }
1338: }
1339: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
1340: this->setConstraintDimension(p, numConstraints);
1341: }
1342: };
1343: // GeneralSection will support BC on a subset of unknowns on a point
1344: // We make a separate constraint Atlas to mark constrained dofs on a point
1345: // Storage will be contiguous by node, just as in Section
1346: // This allows fast restrict(p)
1347: // Then update() is accomplished by skipping constrained unknowns
1348: // We must eliminate restrictSpace() since it does not correspond to the constrained system
1349: // Numbering will have to be rewritten to calculate correct mappings
1350: // I think we can just generate multiple tuples per point
1351: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
1352: typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
1353: typename BCAtlas_ = Section<Point_, int, typename Alloc_::template rebind<int>::other> >
1354: class GeneralSection : public ALE::ParallelObject {
1355: public:
1356: typedef Point_ point_type;
1357: typedef Value_ value_type;
1358: typedef Alloc_ alloc_type;
1359: typedef Atlas_ atlas_type;
1360: typedef BCAtlas_ bc_type;
1361: typedef Point index_type;
1362: typedef typename atlas_type::chart_type chart_type;
1363: typedef value_type * values_type;
1364: typedef std::pair<const int *, const int *> customAtlasInd_type;
1365: typedef std::pair<customAtlasInd_type, bool> customAtlas_type;
1366: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
1367: typedef typename atlas_alloc_type::pointer atlas_ptr;
1368: typedef typename alloc_type::template rebind<bc_type>::other bc_alloc_type;
1369: typedef typename bc_alloc_type::pointer bc_ptr;
1370: protected:
1371: // Describes layout of storage, point --> (# of values, offset into array)
1372: Obj<atlas_type> _atlas;
1373: // Spare atlas to allow dynamic updates
1374: Obj<atlas_type> _atlasNew;
1375: // Storage
1376: values_type _array;
1377: bool _sharedStorage;
1378: int _sharedStorageSize;
1379: // A section that maps points to sets of constrained local dofs
1380: Obj<bc_type> _bc;
1381: alloc_type _allocator;
1382: // Fibration structures
1383: // _spaces is a set of atlases which describe the layout of each in the storage of this section
1384: // _bcs is the same as _bc, but for each field
1385: std::vector<int > _comps;
1386: std::vector<Obj<atlas_type> > _spaces;
1387: std::vector<Obj<bc_type> > _bcs;
1388: // Optimization
1389: std::vector<customAtlas_type> _customRestrictAtlas;
1390: std::vector<customAtlas_type> _customUpdateAtlas;
1391: public:
1392: GeneralSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
1393: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
1394: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
1395: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
1396: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
1397: bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
1398: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
1399: this->_atlasNew = NULL;
1400: this->_array = NULL;
1401: this->_sharedStorage = false;
1402: };
1403: GeneralSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
1404: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
1405: bc_alloc_type(this->_allocator).construct(pBC, bc_type(this->comm(), this->debug()));
1406: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
1407: };
1408: ~GeneralSection() {
1409: if (this->_array && !this->_sharedStorage) {
1410: const int totalSize = this->sizeWithBC();
1412: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1413: this->_allocator.deallocate(this->_array, totalSize);
1414: ///delete [] this->_array;
1415: this->_array = NULL;
1416: }
1417: for(std::vector<customAtlas_type>::iterator a_iter = this->_customRestrictAtlas.begin(); a_iter != this->_customRestrictAtlas.end(); ++a_iter) {
1418: if (a_iter->second) {
1419: delete [] a_iter->first.first;
1420: delete [] a_iter->first.second;
1421: }
1422: }
1423: for(std::vector<customAtlas_type>::iterator a_iter = this->_customUpdateAtlas.begin(); a_iter != this->_customUpdateAtlas.end(); ++a_iter) {
1424: if (a_iter->second) {
1425: delete [] a_iter->first.first;
1426: delete [] a_iter->first.second;
1427: }
1428: }
1429: };
1430: public:
1431: value_type *getRawArray(const int size) {
1432: // Put in a sentinel value that deallocates the array
1433: static value_type *array = NULL;
1434: static int maxSize = 0;
1436: if (size > maxSize) {
1437: const value_type dummy(0);
1439: if (array) {
1440: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
1441: this->_allocator.deallocate(array, maxSize);
1442: ///delete [] array;
1443: }
1444: maxSize = size;
1445: array = this->_allocator.allocate(maxSize);
1446: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
1447: ///array = new value_type[maxSize];
1448: };
1449: return array;
1450: };
1451: int getStorageSize() const {
1452: if (this->_sharedStorage) {
1453: return this->_sharedStorageSize;
1454: }
1455: return this->sizeWithBC();
1456: };
1457: bool sharedStorage() const {return this->_sharedStorage;};
1458: public: // Verifiers
1459: bool hasPoint(const point_type& point) {
1460: return this->_atlas->hasPoint(point);
1461: };
1462: public: // Accessors
1463: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
1464: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
1465: const Obj<atlas_type>& getNewAtlas() const {return this->_atlasNew;};
1466: void setNewAtlas(const Obj<atlas_type>& atlas) {this->_atlasNew = atlas;};
1467: const Obj<bc_type>& getBC() const {return this->_bc;};
1468: void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
1469: const chart_type& getChart() const {return this->_atlas->getChart();};
1470: void setChart(const chart_type& chart) {throw ALE::Exception("setChart() for GeneralSection is invalid");};
1471: public: // BC
1472: // Returns the number of constraints on a point
1473: int getConstraintDimension(const point_type& p) const {
1474: if (!this->_bc->hasPoint(p)) return 0;
1475: return this->_bc->getFiberDimension(p);
1476: }
1477: // Set the number of constraints on a point
1478: void setConstraintDimension(const point_type& p, const int numConstraints) {
1479: this->_bc->setFiberDimension(p, numConstraints);
1480: }
1481: // Increment the number of constraints on a point
1482: void addConstraintDimension(const point_type& p, const int numConstraints) {
1483: this->_bc->addFiberDimension(p, numConstraints);
1484: }
1485: // Return the local dofs which are constrained on a point
1486: const int *getConstraintDof(const point_type& p) const {
1487: return this->_bc->restrictPoint(p);
1488: }
1489: // Set the local dofs which are constrained on a point
1490: void setConstraintDof(const point_type& p, const int dofs[]) {
1491: this->_bc->updatePoint(p, dofs);
1492: }
1493: template<typename OtherSection>
1494: void copyBC(const Obj<OtherSection>& section) {
1495: this->setBC(section->getBC());
1496: const chart_type& chart = this->getChart();
1498: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1499: if (this->getConstraintDimension(*p_iter)) {
1500: this->updatePointBCFull(*p_iter, section->restrictPoint(*p_iter));
1501: }
1502: }
1503: this->copySpaces(section);
1504: }
1505: void defaultConstraintDof() {
1506: const chart_type& chart = this->getChart();
1507: int size = 0;
1509: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1510: size = std::max(size, this->getConstraintDimension(*p_iter));
1511: }
1512: int *dofs = new int[size];
1513: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1514: const int cDim = this->getConstraintDimension(*p_iter);
1516: if (cDim) {
1517: for(int d = 0; d < cDim; ++d) {
1518: dofs[d] = d;
1519: }
1520: this->_bc->updatePoint(*p_iter, dofs);
1521: }
1522: }
1523: delete [] dofs;
1524: };
1525: public: // Sizes
1526: void clear() {
1527: if (!this->_sharedStorage) {
1528: const int totalSize = this->sizeWithBC();
1530: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1531: this->_allocator.deallocate(this->_array, totalSize);
1532: ///delete [] this->_array;
1533: }
1534: this->_array = NULL;
1535: this->_atlas->clear();
1536: this->_bc->clear();
1537: };
1538: // Return the total number of dofs on the point (free and constrained)
1539: int getFiberDimension(const point_type& p) const {
1540: return this->_atlas->restrictPoint(p)->prefix;
1541: };
1542: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
1543: return atlas->restrictPoint(p)->prefix;
1544: };
1545: // Set the total number of dofs on the point (free and constrained)
1546: void setFiberDimension(const point_type& p, int dim) {
1547: const index_type idx(dim, -1);
1548: this->_atlas->addPoint(p);
1549: this->_atlas->updatePoint(p, &idx);
1550: };
1551: template<typename Sequence>
1552: void setFiberDimension(const Obj<Sequence>& points, int dim) {
1553: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1554: this->setFiberDimension(*p_iter, dim);
1555: }
1556: }
1557: void addFiberDimension(const point_type& p, int dim) {
1558: if (this->_atlas->hasPoint(p)) {
1559: const index_type values(dim, 0);
1560: this->_atlas->updateAddPoint(p, &values);
1561: } else {
1562: this->setFiberDimension(p, dim);
1563: }
1564: };
1565: // Return the number of constrained dofs on this point
1566: int getConstrainedFiberDimension(const point_type& p) const {
1567: return this->getFiberDimension(p) - this->getConstraintDimension(p);
1568: };
1569: // Return the total number of free dofs
1570: int size() const {
1571: const chart_type& points = this->getChart();
1572: int size = 0;
1574: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1575: size += this->getConstrainedFiberDimension(*p_iter);
1576: }
1577: return size;
1578: };
1579: // Return the total number of dofs (free and constrained)
1580: int sizeWithBC() const {
1581: const chart_type& points = this->getChart();
1582: int size = 0;
1584: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1585: size += this->getFiberDimension(*p_iter);
1586: }
1587: return size;
1588: };
1589: public: // Index retrieval
1590: const typename index_type::index_type& getIndex(const point_type& p) {
1591: return this->_atlas->restrictPoint(p)->index;
1592: };
1593: void setIndex(const point_type& p, const typename index_type::index_type& index) {
1594: ((typename atlas_type::value_type *) this->_atlas->restrictPoint(p))->index = index;
1595: };
1596: void setIndexBC(const point_type& p, const typename index_type::index_type& index) {};
1597: void getIndicesRaw(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1) {
1598: this->getIndicesRaw(p, this->getIndex(p), indices, indx, orientation);
1599: };
1600: template<typename Order_>
1601: void getIndicesRaw(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1) {
1602: this->getIndicesRaw(p, order->getIndex(p), indices, indx, orientation);
1603: }
1604: void getIndices(const point_type& p, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = true) {
1605: this->getIndices(p, this->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1606: };
1607: template<typename Order_>
1608: void getIndices(const point_type& p, const Obj<Order_>& order, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1609: this->getIndices(p, order->getIndex(p), indices, indx, orientation, freeOnly, skipConstraints);
1610: }
1611: void getIndicesRaw(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation) {
1612: if (orientation >= 0) {
1613: const int& dim = this->getFiberDimension(p);
1614: const int end = start + dim;
1616: for(int i = start; i < end; ++i) {
1617: indices[(*indx)++] = i;
1618: }
1619: } else {
1620: const int numSpaces = this->getNumSpaces();
1621: int offset = start;
1623: for(int space = 0; space < numSpaces; ++space) {
1624: const int& dim = this->getFiberDimension(p, space);
1626: for(int i = offset+dim-1; i >= offset; --i) {
1627: indices[(*indx)++] = i;
1628: }
1629: offset += dim;
1630: }
1631: if (!numSpaces) {
1632: const int& dim = this->getFiberDimension(p);
1634: for(int i = offset+dim-1; i >= offset; --i) {
1635: indices[(*indx)++] = i;
1636: }
1637: offset += dim;
1638: }
1639: }
1640: }
1641: /*
1642: p - The Sieve point
1643: start - Offset for the set of indices on this point
1644: indices - Storage for the indices
1645: indx - Pointer to an offset into the indices argument (to allow composition of index sets)
1646: orientation - A negative argument reverses the indices
1647: freeOnly - If false, include constrained dofs with negative indices, otherwise leave them out
1648: skipConstraints - If true, when a constrained dof is encountered, increment the index (even though it is not placed in indices[])
1649: */
1650: void getIndices(const point_type& p, const int start, PetscInt indices[], PetscInt *indx, const int orientation = 1, const bool freeOnly = false, const bool skipConstraints = false) {
1651: const int& cDim = this->getConstraintDimension(p);
1653: if (!cDim) {
1654: this->getIndicesRaw(p, start, indices, indx, orientation);
1655: } else {
1656: if (orientation >= 0) {
1657: const int& dim = this->getFiberDimension(p);
1658: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1659: int cInd = 0;
1661: /* i is the returned index, k is the local dof number */
1662: for(int i = start, k = 0; k < dim; ++k) {
1663: if ((cInd < cDim) && (k == cDof[cInd])) {
1664: /* Constrained dof */
1665: if (!freeOnly) indices[(*indx)++] = -(k+1);
1666: if (skipConstraints) ++i;
1667: ++cInd;
1668: } else {
1669: /* Unconstrained dof */
1670: indices[(*indx)++] = i++;
1671: }
1672: }
1673: } else {
1674: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1675: int offset = 0;
1676: int cOffset = 0;
1677: int j = -1;
1679: for(int space = 0; space < this->getNumSpaces(); ++space) {
1680: const int dim = this->getFiberDimension(p, space);
1681: const int tDim = this->getConstrainedFiberDimension(p, space);
1682: int cInd = (dim - tDim)-1;
1684: j += dim;
1685: for(int i = 0, k = start+tDim+offset; i < dim; ++i, --j) {
1686: if ((cInd >= 0) && (j == cDof[cInd+cOffset])) {
1687: if (!freeOnly) indices[(*indx)++] = -(offset+i+1);
1688: if (skipConstraints) --k;
1689: --cInd;
1690: } else {
1691: indices[(*indx)++] = --k;
1692: }
1693: }
1694: j += dim;
1695: offset += dim;
1696: cOffset += dim - tDim;
1697: }
1698: }
1699: }
1700: };
1701: public: // Allocation
1702: void allocateStorage() {
1703: const int totalSize = this->sizeWithBC();
1704: const value_type dummy(0) ;
1706: this->_array = this->_allocator.allocate(totalSize);
1707: ///this->_array = new value_type[totalSize];
1708: this->_sharedStorage = false;
1709: this->_sharedStorageSize = 0;
1710: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
1711: ///PetscMemzero(this->_array, totalSize * sizeof(value_type));
1712: this->_bc->allocatePoint();
1713: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = this->_bcs.begin(); b_iter != this->_bcs.end(); ++b_iter) {
1714: (*b_iter)->allocatePoint();;
1715: }
1716: };
1717: void replaceStorage(value_type *newArray, const bool sharedStorage = false, const int sharedStorageSize = 0) {
1718: if (this->_array && !this->_sharedStorage) {
1719: const int totalSize = this->sizeWithBC();
1721: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
1722: this->_allocator.deallocate(this->_array, totalSize);
1723: ///delete [] this->_array;
1724: }
1725: this->_array = newArray;
1726: this->_sharedStorage = sharedStorage;
1727: this->_sharedStorageSize = sharedStorageSize;
1728: this->_atlas = this->_atlasNew;
1729: this->_atlasNew = NULL;
1730: };
1731: // DANGEROUS
1732: void setStorage(value_type *newArray) {
1733: this->_array = newArray;
1734: };
1735: void addPoint(const point_type& point, const int dim) {
1736: if (dim == 0) return;
1737: if (this->_atlasNew.isNull()) {
1738: this->_atlasNew = new atlas_type(this->comm(), this->debug());
1739: this->_atlasNew->copy(this->_atlas);
1740: }
1741: const index_type idx(dim, -1);
1742: this->_atlasNew->addPoint(point);
1743: this->_atlasNew->updatePoint(point, &idx);
1744: };
1745: void orderPoints(const Obj<atlas_type>& atlas){
1746: const chart_type& chart = this->getChart();
1747: int offset = 0;
1749: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1750: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
1751: const int& dim = idx.prefix;
1753: idx.index = offset;
1754: atlas->updatePoint(*c_iter, &idx);
1755: offset += dim;
1756: }
1757: };
1758: void allocatePoint() {
1759: this->orderPoints(this->_atlas);
1760: this->allocateStorage();
1761: };
1762: public: // Restriction and Update
1763: // Zero entries
1764: void zero() {
1765: this->set((value_type) 0.0);
1766: };
1767: void zeroWithBC() {
1768: memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1769: };
1770: void set(const value_type value) {
1771: //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
1772: const chart_type& chart = this->getChart();
1774: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1775: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1776: const int& dim = this->getFiberDimension(*c_iter);
1777: const int& cDim = this->getConstraintDimension(*c_iter);
1779: if (!cDim) {
1780: for(int i = 0; i < dim; ++i) {
1781: array[i] = value;
1782: }
1783: } else {
1784: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1785: int cInd = 0;
1787: for(int i = 0; i < dim; ++i) {
1788: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1789: array[i] = value;
1790: }
1791: }
1792: }
1793: };
1794: // Add two sections and put the result in a third
1795: void add(const Obj<GeneralSection>& x, const Obj<GeneralSection>& y) {
1796: // Check atlases
1797: const chart_type& chart = this->getChart();
1799: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1800: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1801: const value_type *xArray = x->restrictPoint(*c_iter);
1802: const value_type *yArray = y->restrictPoint(*c_iter);
1803: const int& dim = this->getFiberDimension(*c_iter);
1804: const int& cDim = this->getConstraintDimension(*c_iter);
1806: if (!cDim) {
1807: for(int i = 0; i < dim; ++i) {
1808: array[i] = xArray[i] + yArray[i];
1809: }
1810: } else {
1811: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1812: int cInd = 0;
1814: for(int i = 0; i < dim; ++i) {
1815: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1816: array[i] = xArray[i] + yArray[i];
1817: }
1818: }
1819: }
1820: };
1821: // this = this + alpha * x
1822: template<typename OtherSection>
1823: void axpy(const value_type alpha, const Obj<OtherSection>& x) {
1824: // Check atlases
1825: const chart_type& chart = this->getChart();
1827: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1828: value_type *array = (value_type *) this->restrictPoint(*c_iter);
1829: const value_type *xArray = x->restrictPoint(*c_iter);
1830: const int& dim = this->getFiberDimension(*c_iter);
1831: const int& cDim = this->getConstraintDimension(*c_iter);
1833: if (!cDim) {
1834: for(int i = 0; i < dim; ++i) {
1835: array[i] += alpha*xArray[i];
1836: }
1837: } else {
1838: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
1839: int cInd = 0;
1841: for(int i = 0; i < dim; ++i) {
1842: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1843: array[i] += alpha*xArray[i];
1844: }
1845: }
1846: }
1847: }
1848: // Return the free values on a point
1849: const value_type *restrictSpace() const {
1850: return this->_array;
1851: }
1852: // Return the free values on a point
1853: const value_type *restrictPoint(const point_type& p) const {
1854: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1855: }
1856: void restrictPoint(const point_type& p, value_type values[], const int size) const {
1857: assert(this->_atlas->restrictPoint(p)[0].prefix == size);
1858: const value_type *v = &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
1860: for(int i = 0; i < size; ++i) {
1861: values[i] = v[i];
1862: }
1863: };
1864: // Update the free values on a point
1865: // Takes a full array and ignores constrained values
1866: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1867: value_type *array = (value_type *) this->restrictPoint(p);
1868: const int& cDim = this->getConstraintDimension(p);
1870: if (!cDim) {
1871: if (orientation >= 0) {
1872: const int& dim = this->getFiberDimension(p);
1874: for(int i = 0; i < dim; ++i) {
1875: array[i] = v[i];
1876: }
1877: } else {
1878: int offset = 0;
1879: int j = -1;
1881: for(int space = 0; space < this->getNumSpaces(); ++space) {
1882: const int& dim = this->getFiberDimension(p, space);
1884: for(int i = dim-1; i >= 0; --i) {
1885: array[++j] = v[i+offset];
1886: }
1887: offset += dim;
1888: }
1889: }
1890: } else {
1891: if (orientation >= 0) {
1892: const int& dim = this->getFiberDimension(p);
1893: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1894: int cInd = 0;
1896: for(int i = 0; i < dim; ++i) {
1897: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1898: array[i] = v[i];
1899: }
1900: } else {
1901: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1902: int offset = 0;
1903: int cOffset = 0;
1904: int j = 0;
1906: for(int space = 0; space < this->getNumSpaces(); ++space) {
1907: const int dim = this->getFiberDimension(p, space);
1908: const int tDim = this->getConstrainedFiberDimension(p, space);
1909: const int sDim = dim - tDim;
1910: int cInd = 0;
1912: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1913: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1914: array[j] = v[k];
1915: }
1916: offset += dim;
1917: cOffset += dim - tDim;
1918: }
1919: }
1920: }
1921: };
1922: // Update the free values on a point
1923: // Takes a full array and ignores constrained values
1924: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
1925: value_type *array = (value_type *) this->restrictPoint(p);
1926: const int& cDim = this->getConstraintDimension(p);
1928: if (!cDim) {
1929: if (orientation >= 0) {
1930: const int& dim = this->getFiberDimension(p);
1932: for(int i = 0; i < dim; ++i) {
1933: array[i] += v[i];
1934: }
1935: } else {
1936: int offset = 0;
1937: int j = -1;
1939: for(int space = 0; space < this->getNumSpaces(); ++space) {
1940: const int& dim = this->getFiberDimension(p, space);
1942: for(int i = dim-1; i >= 0; --i) {
1943: array[++j] += v[i+offset];
1944: }
1945: offset += dim;
1946: }
1947: }
1948: } else {
1949: if (orientation >= 0) {
1950: const int& dim = this->getFiberDimension(p);
1951: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1952: int cInd = 0;
1954: for(int i = 0; i < dim; ++i) {
1955: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1956: array[i] += v[i];
1957: }
1958: } else {
1959: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
1960: int offset = 0;
1961: int cOffset = 0;
1962: int j = 0;
1964: for(int space = 0; space < this->getNumSpaces(); ++space) {
1965: const int dim = this->getFiberDimension(p, space);
1966: const int tDim = this->getConstrainedFiberDimension(p, space);
1967: const int sDim = dim - tDim;
1968: int cInd = 0;
1970: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1971: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1972: array[j] += v[k];
1973: }
1974: offset += dim;
1975: cOffset += dim - tDim;
1976: }
1977: }
1978: }
1979: };
1980: // Update the free values on a point
1981: // Takes ONLY unconstrained values
1982: void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
1983: value_type *array = (value_type *) this->restrictPoint(p);
1984: const int& cDim = this->getConstraintDimension(p);
1986: if (!cDim) {
1987: if (orientation >= 0) {
1988: const int& dim = this->getFiberDimension(p);
1990: for(int i = 0; i < dim; ++i) {
1991: array[i] = v[i];
1992: }
1993: } else {
1994: int offset = 0;
1995: int j = -1;
1997: for(int space = 0; space < this->getNumSpaces(); ++space) {
1998: const int& dim = this->getFiberDimension(p, space);
2000: for(int i = dim-1; i >= 0; --i) {
2001: array[++j] = v[i+offset];
2002: }
2003: offset += dim;
2004: }
2005: }
2006: } else {
2007: if (orientation >= 0) {
2008: const int& dim = this->getFiberDimension(p);
2009: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2010: int cInd = 0;
2012: for(int i = 0, k = -1; i < dim; ++i) {
2013: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2014: array[i] = v[++k];
2015: }
2016: } else {
2017: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2018: int offset = 0;
2019: int cOffset = 0;
2020: int j = 0;
2022: for(int space = 0; space < this->getNumSpaces(); ++space) {
2023: const int dim = this->getFiberDimension(p, space);
2024: const int tDim = this->getConstrainedFiberDimension(p, space);
2025: const int sDim = dim - tDim;
2026: int cInd = 0;
2028: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2029: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2030: array[j] = v[--k];
2031: }
2032: offset += dim;
2033: cOffset += dim - tDim;
2034: }
2035: }
2036: }
2037: };
2038: // Update the free values on a point
2039: // Takes ONLY unconstrained values
2040: void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2041: value_type *array = (value_type *) this->restrictPoint(p);
2042: const int& cDim = this->getConstraintDimension(p);
2044: if (!cDim) {
2045: if (orientation >= 0) {
2046: const int& dim = this->getFiberDimension(p);
2048: for(int i = 0; i < dim; ++i) {
2049: array[i] += v[i];
2050: }
2051: } else {
2052: int offset = 0;
2053: int j = -1;
2055: for(int space = 0; space < this->getNumSpaces(); ++space) {
2056: const int& dim = this->getFiberDimension(p, space);
2058: for(int i = dim-1; i >= 0; --i) {
2059: array[++j] += v[i+offset];
2060: }
2061: offset += dim;
2062: }
2063: }
2064: } else {
2065: if (orientation >= 0) {
2066: const int& dim = this->getFiberDimension(p);
2067: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2068: int cInd = 0;
2070: for(int i = 0, k = -1; i < dim; ++i) {
2071: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2072: array[i] += v[++k];
2073: }
2074: } else {
2075: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2076: int offset = 0;
2077: int cOffset = 0;
2078: int j = 0;
2080: for(int space = 0; space < this->getNumSpaces(); ++space) {
2081: const int dim = this->getFiberDimension(p, space);
2082: const int tDim = this->getConstrainedFiberDimension(p, space);
2083: const int sDim = dim - tDim;
2084: int cInd = 0;
2086: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2087: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2088: array[j] += v[--k];
2089: }
2090: offset += dim;
2091: cOffset += dim - tDim;
2092: }
2093: }
2094: }
2095: };
2096: // Update only the constrained dofs on a point
2097: // This takes an array with ONLY bc values
2098: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2099: value_type *array = (value_type *) this->restrictPoint(p);
2100: const int& cDim = this->getConstraintDimension(p);
2102: if (cDim) {
2103: if (orientation >= 0) {
2104: const int& dim = this->getFiberDimension(p);
2105: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2106: int cInd = 0;
2108: for(int i = 0; i < dim; ++i) {
2109: if (cInd == cDim) break;
2110: if (i == cDof[cInd]) {
2111: array[i] = v[cInd];
2112: ++cInd;
2113: }
2114: }
2115: } else {
2116: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2117: int cOffset = 0;
2118: int j = 0;
2120: for(int space = 0; space < this->getNumSpaces(); ++space) {
2121: const int dim = this->getFiberDimension(p, space);
2122: const int tDim = this->getConstrainedFiberDimension(p, space);
2123: int cInd = 0;
2125: for(int i = 0; i < dim; ++i, ++j) {
2126: if (cInd < 0) break;
2127: if (j == cDof[cInd+cOffset]) {
2128: array[j] = v[cInd+cOffset];
2129: ++cInd;
2130: }
2131: }
2132: cOffset += dim - tDim;
2133: }
2134: }
2135: }
2136: };
2137: // Update only the constrained dofs on a point
2138: // This takes an array with ALL values, not just BC
2139: void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
2140: value_type *array = (value_type *) this->restrictPoint(p);
2141: const int& cDim = this->getConstraintDimension(p);
2143: if (cDim) {
2144: if (orientation >= 0) {
2145: const int& dim = this->getFiberDimension(p);
2146: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2147: int cInd = 0;
2149: for(int i = 0; i < dim; ++i) {
2150: if (cInd == cDim) break;
2151: if (i == cDof[cInd]) {
2152: array[i] = v[i];
2153: ++cInd;
2154: }
2155: }
2156: } else {
2157: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2158: int offset = 0;
2159: int cOffset = 0;
2160: int j = 0;
2162: for(int space = 0; space < this->getNumSpaces(); ++space) {
2163: const int dim = this->getFiberDimension(p, space);
2164: const int tDim = this->getConstrainedFiberDimension(p, space);
2165: int cInd = 0;
2167: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2168: if (cInd < 0) break;
2169: if (j == cDof[cInd+cOffset]) {
2170: array[j] = v[k];
2171: ++cInd;
2172: }
2173: }
2174: offset += dim;
2175: cOffset += dim - tDim;
2176: }
2177: }
2178: }
2179: };
2180: // Update all dofs on a point (free and constrained)
2181: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
2182: value_type *array = (value_type *) this->restrictPoint(p);
2184: if (orientation >= 0) {
2185: const int& dim = this->getFiberDimension(p);
2187: for(int i = 0; i < dim; ++i) {
2188: array[i] = v[i];
2189: }
2190: } else {
2191: int offset = 0;
2192: int j = -1;
2194: for(int space = 0; space < this->getNumSpaces(); ++space) {
2195: const int& dim = this->getFiberDimension(p, space);
2197: for(int i = dim-1; i >= 0; --i) {
2198: array[++j] = v[i+offset];
2199: }
2200: offset += dim;
2201: }
2202: }
2203: };
2204: // Update all dofs on a point (free and constrained)
2205: void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
2206: value_type *array = (value_type *) this->restrictPoint(p);
2208: if (orientation >= 0) {
2209: const int& dim = this->getFiberDimension(p);
2211: for(int i = 0; i < dim; ++i) {
2212: array[i] += v[i];
2213: }
2214: } else {
2215: int offset = 0;
2216: int j = -1;
2218: for(int space = 0; space < this->getNumSpaces(); ++space) {
2219: const int& dim = this->getFiberDimension(p, space);
2221: for(int i = dim-1; i >= 0; --i) {
2222: array[++j] += v[i+offset];
2223: }
2224: offset += dim;
2225: }
2226: }
2227: };
2228: public: // Fibrations
2229: int getNumSpaces() const {return this->_spaces.size();};
2230: const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
2231: const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
2232: void addSpace(int comp = 1) {
2233: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
2234: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
2235: this->_comps.push_back(comp);
2236: this->_spaces.push_back(space);
2237: this->_bcs.push_back(bc);
2238: };
2239: int getSpaceComponents(const int space) {return this->_comps[space];};
2240: int getFiberDimension(const point_type& p, const int space) const {
2241: return this->_spaces[space]->restrictPoint(p)->prefix;
2242: };
2243: void setFiberDimension(const point_type& p, int dim, const int space) {
2244: const index_type idx(dim, -1);
2245: this->_spaces[space]->addPoint(p);
2246: this->_spaces[space]->updatePoint(p, &idx);
2247: };
2248: template<typename Sequence>
2249: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
2250: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2251: this->setFiberDimension(*p_iter, dim, space);
2252: }
2253: }
2254: int getConstraintDimension(const point_type& p, const int space) const {
2255: if (!this->_bcs[space]->hasPoint(p)) return 0;
2256: return this->_bcs[space]->getFiberDimension(p);
2257: }
2258: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2259: this->_bcs[space]->setFiberDimension(p, numConstraints);
2260: }
2261: void addConstraintDimension(const point_type& p, const int numConstraints, const int space) {
2262: this->_bcs[space]->addFiberDimension(p, numConstraints);
2263: }
2264: int getConstrainedFiberDimension(const point_type& p, const int space) const {
2265: return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
2266: }
2267: // Return the local dofs which are constrained on a point
2268: const int *getConstraintDof(const point_type& p, const int space) const {
2269: return this->_bcs[space]->restrictPoint(p);
2270: }
2271: // Set the local dofs which are constrained on a point
2272: void setConstraintDof(const point_type& p, const int dofs[], const int space) {
2273: this->_bcs[space]->updatePoint(p, dofs);
2274: }
2275: // Return the total number of free dofs
2276: int size(const int space) const {
2277: const chart_type& points = this->getChart();
2278: int size = 0;
2280: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2281: size += this->getConstrainedFiberDimension(*p_iter, space);
2282: }
2283: return size;
2284: };
2285: template<typename OtherSection>
2286: void copySpaces(const Obj<OtherSection>& section) {
2287: const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
2288: const std::vector<Obj<bc_type> >& bcs = section->getBCs();
2290: this->_spaces.clear();
2291: for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
2292: this->_spaces.push_back(*s_iter);
2293: }
2294: this->_bcs.clear();
2295: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
2296: this->_bcs.push_back(*b_iter);
2297: }
2298: }
2299: Obj<GeneralSection> getFibration(const int space) const {
2300: Obj<GeneralSection> field = new GeneralSection(this->comm(), this->debug());
2301: // Obj<atlas_type> _atlas;
2302: // std::vector<Obj<atlas_type> > _spaces;
2303: // Obj<bc_type> _bc;
2304: // std::vector<Obj<bc_type> > _bcs;
2305: field->addSpace();
2306: const chart_type& chart = this->getChart();
2308: // Copy sizes
2309: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2310: const int fDim = this->getFiberDimension(*c_iter, space);
2311: const int cDim = this->getConstraintDimension(*c_iter, space);
2313: if (fDim) {
2314: field->setFiberDimension(*c_iter, fDim);
2315: field->setFiberDimension(*c_iter, fDim, 0);
2316: }
2317: if (cDim) {
2318: field->setConstraintDimension(*c_iter, cDim);
2319: field->setConstraintDimension(*c_iter, cDim, 0);
2320: }
2321: }
2322: field->allocateStorage();
2323: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
2324: const chart_type& newChart = field->getChart();
2326: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2327: const int cDim = field->getConstraintDimension(*c_iter);
2328: //const int dof[1] = {0};
2330: if (cDim) {
2331: field->setConstraintDof(*c_iter, this->getConstraintDof(*c_iter, space));
2332: }
2333: }
2334: // Copy offsets
2335: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
2336: index_type idx;
2338: idx.prefix = field->getFiberDimension(*c_iter);
2339: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
2340: for(int s = 0; s < space; ++s) {
2341: idx.index += this->getFiberDimension(*c_iter, s);
2342: }
2343: newAtlas->addPoint(*c_iter);
2344: newAtlas->updatePoint(*c_iter, &idx);
2345: }
2346: field->replaceStorage(this->_array, true, this->getStorageSize());
2347: field->setAtlas(newAtlas);
2348: return field;
2349: };
2350: public: // Optimization
2351: void getCustomRestrictAtlas(const int tag, const int *offsets[], const int *indices[]) {
2352: *offsets = this->_customRestrictAtlas[tag].first.first;
2353: *indices = this->_customRestrictAtlas[tag].first.second;
2354: };
2355: void getCustomUpdateAtlas(const int tag, const int *offsets[], const int *indices[]) {
2356: *offsets = this->_customUpdateAtlas[tag].first.first;
2357: *indices = this->_customUpdateAtlas[tag].first.second;
2358: };
2359: // This returns the tag assigned to the atlas
2360: int setCustomAtlas(const int restrictOffsets[], const int restrictIndices[], const int updateOffsets[], const int updateIndices[], bool autoFree = true) {
2361: this->_customRestrictAtlas.push_back(customAtlas_type(customAtlasInd_type(restrictOffsets, restrictIndices), autoFree));
2362: this->_customUpdateAtlas.push_back(customAtlas_type(customAtlasInd_type(updateOffsets, updateIndices), autoFree));
2363: return this->_customUpdateAtlas.size()-1;
2364: };
2365: int copyCustomAtlas(const Obj<GeneralSection>& section, const int tag) {
2366: const int *rOffsets, *rIndices, *uOffsets, *uIndices;
2368: section->getCustomRestrictAtlas(tag, &rOffsets, &rIndices);
2369: section->getCustomUpdateAtlas(tag, &uOffsets, &uIndices);
2370: return this->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices, false);
2371: };
2372: public:
2373: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
2374: ostringstream txt;
2375: int rank;
2377: if (comm == MPI_COMM_NULL) {
2378: comm = this->comm();
2379: rank = this->commRank();
2380: } else {
2381: MPI_Comm_rank(comm, &rank);
2382: }
2383: if (name == "") {
2384: if(rank == 0) {
2385: txt << "viewing a GeneralSection" << std::endl;
2386: }
2387: } else {
2388: if (rank == 0) {
2389: txt << "viewing GeneralSection '" << name << "'" << std::endl;
2390: }
2391: }
2392: if (rank == 0) {
2393: txt << " Fields: " << this->getNumSpaces() << std::endl;
2394: }
2395: const chart_type& chart = this->getChart();
2397: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2398: const value_type *array = this->restrictPoint(*p_iter);
2399: const int& dim = this->getFiberDimension(*p_iter);
2401: if (dim != 0) {
2402: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << " ";
2403: for(int i = 0; i < dim; i++) {
2404: txt << " " << array[i];
2405: }
2406: const int& dim = this->getConstraintDimension(*p_iter);
2408: if (dim) {
2409: const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);
2411: txt << " constrained";
2412: for(int i = 0; i < dim; ++i) {
2413: txt << " " << bcArray[i];
2414: }
2415: }
2416: txt << std::endl;
2417: }
2418: }
2419: if (chart.size() == 0) {
2420: txt << "[" << this->commRank() << "]: empty" << std::endl;
2421: }
2422: PetscSynchronizedPrintf(comm, txt.str().c_str());
2423: PetscSynchronizedFlush(comm);
2424: };
2425: };
2426: // FEMSection will support vector BC on a subset of unknowns on a point
2427: // We make a separate constraint Section to hold the transformation and projection operators
2428: // Storage will be contiguous by node, just as in Section
2429: // This allows fast restrict(p)
2430: // Then update() is accomplished by projecting constrained unknowns
2431: template<typename Point_, typename Value_, typename Alloc_ = malloc_allocator<Value_>,
2432: typename Atlas_ = UniformSection<Point_, Point, 1, typename Alloc_::template rebind<Point>::other>,
2433: typename BCAtlas_ = UniformSection<Point_, int, 1, typename Alloc_::template rebind<int>::other>,
2434: typename BC_ = Section<Point_, double, typename Alloc_::template rebind<double>::other> >
2435: class FEMSection : public ALE::ParallelObject {
2436: public:
2437: typedef Point_ point_type;
2438: typedef Value_ value_type;
2439: typedef Alloc_ alloc_type;
2440: typedef Atlas_ atlas_type;
2441: typedef BCAtlas_ bc_atlas_type;
2442: typedef BC_ bc_type;
2443: typedef Point index_type;
2444: typedef typename atlas_type::chart_type chart_type;
2445: typedef value_type * values_type;
2446: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
2447: typedef typename atlas_alloc_type::pointer atlas_ptr;
2448: typedef typename alloc_type::template rebind<bc_type>::other bc_atlas_alloc_type;
2449: typedef typename bc_atlas_alloc_type::pointer bc_atlas_ptr;
2450: typedef typename alloc_type::template rebind<bc_type>::other bc_alloc_type;
2451: typedef typename bc_alloc_type::pointer bc_ptr;
2452: protected:
2453: Obj<atlas_type> _atlas;
2454: Obj<bc_atlas_type> _bc_atlas;
2455: Obj<bc_type> _bc;
2456: values_type _array;
2457: bool _sharedStorage;
2458: int _sharedStorageSize;
2459: alloc_type _allocator;
2460: std::vector<Obj<atlas_type> > _spaces;
2461: std::vector<Obj<bc_type> > _bcs;
2462: public:
2463: FEMSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
2464: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
2465: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
2466: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
2467: bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2468: bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(comm, debug));
2469: this->_bc_atlas = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2470: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
2471: bc_alloc_type(this->_allocator).construct(pBC, bc_type(comm, debug));
2472: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
2473: this->_array = NULL;
2474: this->_sharedStorage = false;
2475: };
2476: FEMSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _array(NULL), _sharedStorage(false), _sharedStorageSize(0) {
2477: bc_atlas_ptr pBCAtlas = bc_atlas_alloc_type(this->_allocator).allocate(1);
2478: bc_atlas_alloc_type(this->_allocator).construct(pBCAtlas, bc_atlas_type(this->comm(), this->debug()));
2479: this->_bc_atlas = Obj<bc_atlas_type>(pBCAtlas, sizeof(bc_atlas_type));
2480: bc_ptr pBC = bc_alloc_type(this->_allocator).allocate(1);
2481: bc_alloc_type(this->_allocator).construct(pBC, bc_type(this->comm(), this->debug()));
2482: this->_bc = Obj<bc_type>(pBC, sizeof(bc_type));
2483: };
2484: ~FEMSection() {
2485: if (this->_array && !this->_sharedStorage) {
2486: const int totalSize = this->sizeWithBC();
2488: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2489: this->_allocator.deallocate(this->_array, totalSize);
2490: this->_array = NULL;
2491: }
2492: };
2493: public:
2494: value_type *getRawArray(const int size) {
2495: // Put in a sentinel value that deallocates the array
2496: static value_type *array = NULL;
2497: static int maxSize = 0;
2499: if (size > maxSize) {
2500: const value_type dummy(0);
2502: if (array) {
2503: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
2504: this->_allocator.deallocate(array, maxSize);
2505: }
2506: maxSize = size;
2507: array = this->_allocator.allocate(maxSize);
2508: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
2509: };
2510: return array;
2511: };
2512: int getStorageSize() const {
2513: if (this->_sharedStorage) {
2514: return this->_sharedStorageSize;
2515: }
2516: return this->sizeWithBC();
2517: };
2518: public: // Verifiers
2519: bool hasPoint(const point_type& point) {
2520: return this->_atlas->hasPoint(point);
2521: };
2522: public: // Accessors
2523: const chart_type& getChart() const {return this->_atlas->getChart();};
2524: const Obj<atlas_type>& getAtlas() const {return this->_atlas;};
2525: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
2526: const Obj<bc_type>& getBC() const {return this->_bc;};
2527: void setBC(const Obj<bc_type>& bc) {this->_bc = bc;};
2528: public: // BC
2529: // Returns the number of constraints on a point
2530: int getConstraintDimension(const point_type& p) const {
2531: if (!this->_bc_atlas->hasPoint(p)) return 0;
2532: return this->_bc_atlas->restrictPoint(p)[0];
2533: };
2534: // Set the number of constraints on a point
2535: void setConstraintDimension(const point_type& p, const int numConstraints) {
2536: this->_bc_atlas->updatePoint(p, &numConstraints);
2537: };
2538: // Increment the number of constraints on a point
2539: void addConstraintDimension(const point_type& p, const int numConstraints) {
2540: this->_bc_atlas->updatePointAdd(p, &numConstraints);
2541: };
2542: const int *getConstraintDof(const point_type& p) const {
2543: return this->_bc->restrictPoint(p);
2544: }
2545: public: // Sizes
2546: void clear() {
2547: if (!this->_sharedStorage) {
2548: const int totalSize = this->sizeWithBC();
2550: for(int i = 0; i < totalSize; ++i) {this->_allocator.destroy(this->_array+i);}
2551: this->_allocator.deallocate(this->_array, totalSize);
2552: }
2553: this->_array = NULL;
2554: this->_atlas->clear();
2555: this->_bc_atlas->clear();
2556: this->_bc->clear();
2557: };
2558: // Return the total number of dofs on the point (free and constrained)
2559: int getFiberDimension(const point_type& p) const {
2560: return this->_atlas->restrictPoint(p)->prefix;
2561: };
2562: int getFiberDimension(const Obj<atlas_type>& atlas, const point_type& p) const {
2563: return atlas->restrictPoint(p)->prefix;
2564: };
2565: // Set the total number of dofs on the point (free and constrained)
2566: void setFiberDimension(const point_type& p, int dim) {
2567: const index_type idx(dim, -1);
2568: this->_atlas->addPoint(p);
2569: this->_atlas->updatePoint(p, &idx);
2570: };
2571: template<typename Sequence>
2572: void setFiberDimension(const Obj<Sequence>& points, int dim) {
2573: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
2574: this->setFiberDimension(*p_iter, dim);
2575: }
2576: }
2577: void addFiberDimension(const point_type& p, int dim) {
2578: if (this->_atlas->hasPoint(p)) {
2579: const index_type values(dim, 0);
2580: this->_atlas->updateAddPoint(p, &values);
2581: } else {
2582: this->setFiberDimension(p, dim);
2583: }
2584: };
2585: // Return the number of constrained dofs on this point
2586: int getConstrainedFiberDimension(const point_type& p) const {
2587: return this->getFiberDimension(p) - this->getConstraintDimension(p);
2588: };
2589: // Return the total number of free dofs
2590: int size() const {
2591: const chart_type& points = this->getChart();
2592: int size = 0;
2594: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2595: size += this->getConstrainedFiberDimension(*p_iter);
2596: }
2597: return size;
2598: };
2599: // Return the total number of dofs (free and constrained)
2600: int sizeWithBC() const {
2601: const chart_type& points = this->getChart();
2602: int size = 0;
2604: for(typename chart_type::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2605: size += this->getFiberDimension(*p_iter);
2606: }
2607: return size;
2608: };
2609: public: // Allocation
2610: void allocateStorage() {
2611: const int totalSize = this->sizeWithBC();
2612: const value_type dummy(0) ;
2614: this->_array = this->_allocator.allocate(totalSize);
2615: this->_sharedStorage = false;
2616: this->_sharedStorageSize = 0;
2617: for(int i = 0; i < totalSize; ++i) {this->_allocator.construct(this->_array+i, dummy);}
2618: this->_bc_atlas->allocatePoint();
2619: this->_bc->allocatePoint();
2620: };
2621: void orderPoints(const Obj<atlas_type>& atlas){
2622: const chart_type& chart = this->getChart();
2623: int offset = 0;
2625: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2626: typename atlas_type::value_type idx = atlas->restrictPoint(*c_iter)[0];
2627: const int& dim = idx.prefix;
2629: idx.index = offset;
2630: atlas->updatePoint(*c_iter, &idx);
2631: offset += dim;
2632: }
2633: };
2634: void allocatePoint() {
2635: this->orderPoints(this->_atlas);
2636: this->allocateStorage();
2637: };
2638: public: // Restriction and Update
2639: // Zero entries
2640: void zero() {
2641: this->set(0.0);
2642: };
2643: void set(const value_type value) {
2644: //memset(this->_array, 0, this->sizeWithBC()* sizeof(value_type));
2645: const chart_type& chart = this->getChart();
2647: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2648: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2649: const int& dim = this->getFiberDimension(*c_iter);
2650: const int& cDim = this->getConstraintDimension(*c_iter);
2652: if (!cDim) {
2653: for(int i = 0; i < dim; ++i) {
2654: array[i] = value;
2655: }
2656: } else {
2657: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2658: int cInd = 0;
2660: for(int i = 0; i < dim; ++i) {
2661: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2662: array[i] = value;
2663: }
2664: }
2665: }
2666: };
2667: // Add two sections and put the result in a third
2668: void add(const Obj<FEMSection>& x, const Obj<FEMSection>& y) {
2669: // Check atlases
2670: const chart_type& chart = this->getChart();
2672: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2673: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2674: const value_type *xArray = x->restrictPoint(*c_iter);
2675: const value_type *yArray = y->restrictPoint(*c_iter);
2676: const int& dim = this->getFiberDimension(*c_iter);
2677: const int& cDim = this->getConstraintDimension(*c_iter);
2679: if (!cDim) {
2680: for(int i = 0; i < dim; ++i) {
2681: array[i] = xArray[i] + yArray[i];
2682: }
2683: } else {
2684: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2685: int cInd = 0;
2687: for(int i = 0; i < dim; ++i) {
2688: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2689: array[i] = xArray[i] + yArray[i];
2690: }
2691: }
2692: }
2693: };
2694: // this = this + alpha * x
2695: void axpy(const value_type alpha, const Obj<FEMSection>& x) {
2696: // Check atlases
2697: const chart_type& chart = this->getChart();
2699: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2700: value_type *array = (value_type *) this->restrictPoint(*c_iter);
2701: const value_type *xArray = x->restrictPoint(*c_iter);
2702: const int& dim = this->getFiberDimension(*c_iter);
2703: const int& cDim = this->getConstraintDimension(*c_iter);
2705: if (!cDim) {
2706: for(int i = 0; i < dim; ++i) {
2707: array[i] += alpha*xArray[i];
2708: }
2709: } else {
2710: const typename bc_type::value_type *cDof = this->getConstraintDof(*c_iter);
2711: int cInd = 0;
2713: for(int i = 0; i < dim; ++i) {
2714: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2715: array[i] += alpha*xArray[i];
2716: }
2717: }
2718: }
2719: };
2720: // Return the free values on a point
2721: const value_type *restrictSpace() const {
2722: return this->_array;
2723: };
2724: // Return the free values on a point
2725: const value_type *restrictPoint(const point_type& p) const {
2726: return &(this->_array[this->_atlas->restrictPoint(p)[0].index]);
2727: };
2728: // Update the free values on a point
2729: // Takes a full array and ignores constrained values
2730: void updatePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2731: value_type *array = (value_type *) this->restrictPoint(p);
2732: const int& cDim = this->getConstraintDimension(p);
2734: if (!cDim) {
2735: if (orientation >= 0) {
2736: const int& dim = this->getFiberDimension(p);
2738: for(int i = 0; i < dim; ++i) {
2739: array[i] = v[i];
2740: }
2741: } else {
2742: int offset = 0;
2743: int j = -1;
2745: for(int space = 0; space < this->getNumSpaces(); ++space) {
2746: const int& dim = this->getFiberDimension(p, space);
2748: for(int i = dim-1; i >= 0; --i) {
2749: array[++j] = v[i+offset];
2750: }
2751: offset += dim;
2752: }
2753: }
2754: } else {
2755: if (orientation >= 0) {
2756: const int& dim = this->getFiberDimension(p);
2757: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2758: int cInd = 0;
2760: for(int i = 0; i < dim; ++i) {
2761: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2762: array[i] = v[i];
2763: }
2764: } else {
2765: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2766: int offset = 0;
2767: int cOffset = 0;
2768: int j = 0;
2770: for(int space = 0; space < this->getNumSpaces(); ++space) {
2771: const int dim = this->getFiberDimension(p, space);
2772: const int tDim = this->getConstrainedFiberDimension(p, space);
2773: const int sDim = dim - tDim;
2774: int cInd = 0;
2776: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2777: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2778: array[j] = v[k];
2779: }
2780: offset += dim;
2781: cOffset += dim - tDim;
2782: }
2783: }
2784: }
2785: };
2786: // Update the free values on a point
2787: // Takes a full array and ignores constrained values
2788: void updateAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2789: value_type *array = (value_type *) this->restrictPoint(p);
2790: const int& cDim = this->getConstraintDimension(p);
2792: if (!cDim) {
2793: if (orientation >= 0) {
2794: const int& dim = this->getFiberDimension(p);
2796: for(int i = 0; i < dim; ++i) {
2797: array[i] += v[i];
2798: }
2799: } else {
2800: int offset = 0;
2801: int j = -1;
2803: for(int space = 0; space < this->getNumSpaces(); ++space) {
2804: const int& dim = this->getFiberDimension(p, space);
2806: for(int i = dim-1; i >= 0; --i) {
2807: array[++j] += v[i+offset];
2808: }
2809: offset += dim;
2810: }
2811: }
2812: } else {
2813: if (orientation >= 0) {
2814: const int& dim = this->getFiberDimension(p);
2815: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2816: int cInd = 0;
2818: for(int i = 0; i < dim; ++i) {
2819: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2820: array[i] += v[i];
2821: }
2822: } else {
2823: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2824: int offset = 0;
2825: int cOffset = 0;
2826: int j = 0;
2828: for(int space = 0; space < this->getNumSpaces(); ++space) {
2829: const int dim = this->getFiberDimension(p, space);
2830: const int tDim = this->getConstrainedFiberDimension(p, space);
2831: const int sDim = dim - tDim;
2832: int cInd = 0;
2834: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2835: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2836: array[j] += v[k];
2837: }
2838: offset += dim;
2839: cOffset += dim - tDim;
2840: }
2841: }
2842: }
2843: };
2844: // Update the free values on a point
2845: // Takes ONLY unconstrained values
2846: void updateFreePoint(const point_type& p, const value_type v[], const int orientation = 1) {
2847: value_type *array = (value_type *) this->restrictPoint(p);
2848: const int& cDim = this->getConstraintDimension(p);
2850: if (!cDim) {
2851: if (orientation >= 0) {
2852: const int& dim = this->getFiberDimension(p);
2854: for(int i = 0; i < dim; ++i) {
2855: array[i] = v[i];
2856: }
2857: } else {
2858: int offset = 0;
2859: int j = -1;
2861: for(int space = 0; space < this->getNumSpaces(); ++space) {
2862: const int& dim = this->getFiberDimension(p, space);
2864: for(int i = dim-1; i >= 0; --i) {
2865: array[++j] = v[i+offset];
2866: }
2867: offset += dim;
2868: }
2869: }
2870: } else {
2871: if (orientation >= 0) {
2872: const int& dim = this->getFiberDimension(p);
2873: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2874: int cInd = 0;
2876: for(int i = 0, k = -1; i < dim; ++i) {
2877: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2878: array[i] = v[++k];
2879: }
2880: } else {
2881: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2882: int offset = 0;
2883: int cOffset = 0;
2884: int j = 0;
2886: for(int space = 0; space < this->getNumSpaces(); ++space) {
2887: const int dim = this->getFiberDimension(p, space);
2888: const int tDim = this->getConstrainedFiberDimension(p, space);
2889: const int sDim = dim - tDim;
2890: int cInd = 0;
2892: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2893: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2894: array[j] = v[--k];
2895: }
2896: offset += dim;
2897: cOffset += dim - tDim;
2898: }
2899: }
2900: }
2901: };
2902: // Update the free values on a point
2903: // Takes ONLY unconstrained values
2904: void updateFreeAddPoint(const point_type& p, const value_type v[], const int orientation = 1) {
2905: value_type *array = (value_type *) this->restrictPoint(p);
2906: const int& cDim = this->getConstraintDimension(p);
2908: if (!cDim) {
2909: if (orientation >= 0) {
2910: const int& dim = this->getFiberDimension(p);
2912: for(int i = 0; i < dim; ++i) {
2913: array[i] += v[i];
2914: }
2915: } else {
2916: int offset = 0;
2917: int j = -1;
2919: for(int space = 0; space < this->getNumSpaces(); ++space) {
2920: const int& dim = this->getFiberDimension(p, space);
2922: for(int i = dim-1; i >= 0; --i) {
2923: array[++j] += v[i+offset];
2924: }
2925: offset += dim;
2926: }
2927: }
2928: } else {
2929: if (orientation >= 0) {
2930: const int& dim = this->getFiberDimension(p);
2931: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2932: int cInd = 0;
2934: for(int i = 0, k = -1; i < dim; ++i) {
2935: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2936: array[i] += v[++k];
2937: }
2938: } else {
2939: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2940: int offset = 0;
2941: int cOffset = 0;
2942: int j = 0;
2944: for(int space = 0; space < this->getNumSpaces(); ++space) {
2945: const int dim = this->getFiberDimension(p, space);
2946: const int tDim = this->getConstrainedFiberDimension(p, space);
2947: const int sDim = dim - tDim;
2948: int cInd = 0;
2950: for(int i = 0, k = tDim+offset-1; i < dim; ++i, ++j) {
2951: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2952: array[j] += v[--k];
2953: }
2954: offset += dim;
2955: cOffset += dim - tDim;
2956: }
2957: }
2958: }
2959: };
2960: // Update only the constrained dofs on a point
2961: // This takes an array with ONLY bc values
2962: void updatePointBC(const point_type& p, const value_type v[], const int orientation = 1) {
2963: value_type *array = (value_type *) this->restrictPoint(p);
2964: const int& cDim = this->getConstraintDimension(p);
2966: if (cDim) {
2967: if (orientation >= 0) {
2968: const int& dim = this->getFiberDimension(p);
2969: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2970: int cInd = 0;
2972: for(int i = 0; i < dim; ++i) {
2973: if (cInd == cDim) break;
2974: if (i == cDof[cInd]) {
2975: array[i] = v[cInd];
2976: ++cInd;
2977: }
2978: }
2979: } else {
2980: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
2981: int cOffset = 0;
2982: int j = 0;
2984: for(int space = 0; space < this->getNumSpaces(); ++space) {
2985: const int dim = this->getFiberDimension(p, space);
2986: const int tDim = this->getConstrainedFiberDimension(p, space);
2987: int cInd = 0;
2989: for(int i = 0; i < dim; ++i, ++j) {
2990: if (cInd < 0) break;
2991: if (j == cDof[cInd+cOffset]) {
2992: array[j] = v[cInd+cOffset];
2993: ++cInd;
2994: }
2995: }
2996: cOffset += dim - tDim;
2997: }
2998: }
2999: }
3000: };
3001: // Update only the constrained dofs on a point
3002: // This takes an array with ALL values, not just BC
3003: void updatePointBCFull(const point_type& p, const value_type v[], const int orientation = 1) {
3004: value_type *array = (value_type *) this->restrictPoint(p);
3005: const int& cDim = this->getConstraintDimension(p);
3007: if (cDim) {
3008: if (orientation >= 0) {
3009: const int& dim = this->getFiberDimension(p);
3010: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
3011: int cInd = 0;
3013: for(int i = 0; i < dim; ++i) {
3014: if (cInd == cDim) break;
3015: if (i == cDof[cInd]) {
3016: array[i] = v[i];
3017: ++cInd;
3018: }
3019: }
3020: } else {
3021: const typename bc_type::value_type *cDof = this->getConstraintDof(p);
3022: int offset = 0;
3023: int cOffset = 0;
3024: int j = 0;
3026: for(int space = 0; space < this->getNumSpaces(); ++space) {
3027: const int dim = this->getFiberDimension(p, space);
3028: const int tDim = this->getConstrainedFiberDimension(p, space);
3029: int cInd = 0;
3031: for(int i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
3032: if (cInd < 0) break;
3033: if (j == cDof[cInd+cOffset]) {
3034: array[j] = v[k];
3035: ++cInd;
3036: }
3037: }
3038: offset += dim;
3039: cOffset += dim - tDim;
3040: }
3041: }
3042: }
3043: };
3044: // Update all dofs on a point (free and constrained)
3045: void updatePointAll(const point_type& p, const value_type v[], const int orientation = 1) {
3046: value_type *array = (value_type *) this->restrictPoint(p);
3048: if (orientation >= 0) {
3049: const int& dim = this->getFiberDimension(p);
3051: for(int i = 0; i < dim; ++i) {
3052: array[i] = v[i];
3053: }
3054: } else {
3055: int offset = 0;
3056: int j = -1;
3058: for(int space = 0; space < this->getNumSpaces(); ++space) {
3059: const int& dim = this->getFiberDimension(p, space);
3061: for(int i = dim-1; i >= 0; --i) {
3062: array[++j] = v[i+offset];
3063: }
3064: offset += dim;
3065: }
3066: }
3067: };
3068: // Update all dofs on a point (free and constrained)
3069: void updatePointAllAdd(const point_type& p, const value_type v[], const int orientation = 1) {
3070: value_type *array = (value_type *) this->restrictPoint(p);
3072: if (orientation >= 0) {
3073: const int& dim = this->getFiberDimension(p);
3075: for(int i = 0; i < dim; ++i) {
3076: array[i] += v[i];
3077: }
3078: } else {
3079: int offset = 0;
3080: int j = -1;
3082: for(int space = 0; space < this->getNumSpaces(); ++space) {
3083: const int& dim = this->getFiberDimension(p, space);
3085: for(int i = dim-1; i >= 0; --i) {
3086: array[++j] += v[i+offset];
3087: }
3088: offset += dim;
3089: }
3090: }
3091: };
3092: public: // Fibrations
3093: int getNumSpaces() const {return this->_spaces.size();};
3094: const std::vector<Obj<atlas_type> >& getSpaces() {return this->_spaces;};
3095: const std::vector<Obj<bc_type> >& getBCs() {return this->_bcs;};
3096: void addSpace() {
3097: Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
3098: Obj<bc_type> bc = new bc_type(this->comm(), this->debug());
3099: this->_spaces.push_back(space);
3100: this->_bcs.push_back(bc);
3101: };
3102: int getFiberDimension(const point_type& p, const int space) const {
3103: return this->_spaces[space]->restrictPoint(p)->prefix;
3104: };
3105: void setFiberDimension(const point_type& p, int dim, const int space) {
3106: const index_type idx(dim, -1);
3107: this->_spaces[space]->addPoint(p);
3108: this->_spaces[space]->updatePoint(p, &idx);
3109: };
3110: template<typename Sequence>
3111: void setFiberDimension(const Obj<Sequence>& points, int dim, const int space) {
3112: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
3113: this->setFiberDimension(*p_iter, dim, space);
3114: }
3115: }
3116: int getConstraintDimension(const point_type& p, const int space) const {
3117: if (!this->_bcs[space]->hasPoint(p)) return 0;
3118: return this->_bcs[space]->getFiberDimension(p);
3119: };
3120: void setConstraintDimension(const point_type& p, const int numConstraints, const int space) {
3121: this->_bcs[space]->setFiberDimension(p, numConstraints);
3122: };
3123: int getConstrainedFiberDimension(const point_type& p, const int space) const {
3124: return this->getFiberDimension(p, space) - this->getConstraintDimension(p, space);
3125: };
3126: void copyFibration(const Obj<FEMSection>& section) {
3127: const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
3128: const std::vector<Obj<bc_type> >& bcs = section->getBCs();
3130: this->_spaces.clear();
3131: for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter = spaces.begin(); s_iter != spaces.end(); ++s_iter) {
3132: this->_spaces.push_back(*s_iter);
3133: }
3134: this->_bcs.clear();
3135: for(typename std::vector<Obj<bc_type> >::const_iterator b_iter = bcs.begin(); b_iter != bcs.end(); ++b_iter) {
3136: this->_bcs.push_back(*b_iter);
3137: }
3138: };
3139: Obj<FEMSection> getFibration(const int space) const {
3140: Obj<FEMSection> field = new FEMSection(this->comm(), this->debug());
3141: // Obj<atlas_type> _atlas;
3142: // std::vector<Obj<atlas_type> > _spaces;
3143: // Obj<bc_type> _bc;
3144: // std::vector<Obj<bc_type> > _bcs;
3145: field->addSpace();
3146: const chart_type& chart = this->getChart();
3148: // Copy sizes
3149: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3150: const int fDim = this->getFiberDimension(*c_iter, space);
3151: const int cDim = this->getConstraintDimension(*c_iter, space);
3153: if (fDim) {
3154: field->setFiberDimension(*c_iter, fDim);
3155: field->setFiberDimension(*c_iter, fDim, 0);
3156: }
3157: if (cDim) {
3158: field->setConstraintDimension(*c_iter, cDim);
3159: field->setConstraintDimension(*c_iter, cDim, 0);
3160: }
3161: }
3162: field->allocateStorage();
3163: Obj<atlas_type> newAtlas = new atlas_type(this->comm(), this->debug());
3164: const chart_type& newChart = field->getChart();
3166: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3167: const int cDim = field->getConstraintDimension(*c_iter);
3168: const int dof[1] = {0};
3170: if (cDim) {
3171: field->setConstraintDof(*c_iter, dof);
3172: }
3173: }
3174: // Copy offsets
3175: for(typename chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
3176: index_type idx;
3178: idx.prefix = field->getFiberDimension(*c_iter);
3179: idx.index = this->_atlas->restrictPoint(*c_iter)[0].index;
3180: for(int s = 0; s < space; ++s) {
3181: idx.index += this->getFiberDimension(*c_iter, s);
3182: }
3183: newAtlas->addPoint(*c_iter);
3184: newAtlas->updatePoint(*c_iter, &idx);
3185: }
3186: field->replaceStorage(this->_array, true, this->getStorageSize());
3187: field->setAtlas(newAtlas);
3188: return field;
3189: };
3190: public:
3191: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3192: ostringstream txt;
3193: int rank;
3195: if (comm == MPI_COMM_NULL) {
3196: comm = this->comm();
3197: rank = this->commRank();
3198: } else {
3199: MPI_Comm_rank(comm, &rank);
3200: }
3201: if (name == "") {
3202: if(rank == 0) {
3203: txt << "viewing a FEMSection" << std::endl;
3204: }
3205: } else {
3206: if (rank == 0) {
3207: txt << "viewing FEMSection '" << name << "'" << std::endl;
3208: }
3209: }
3210: if (rank == 0) {
3211: txt << " Fields: " << this->getNumSpaces() << std::endl;
3212: }
3213: const chart_type& chart = this->getChart();
3215: for(typename chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
3216: const value_type *array = this->restrictPoint(*p_iter);
3217: const int& dim = this->getFiberDimension(*p_iter);
3219: if (dim != 0) {
3220: txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << dim << " offset " << this->_atlas->restrictPoint(*p_iter)->index << " ";
3221: for(int i = 0; i < dim; i++) {
3222: txt << " " << array[i];
3223: }
3224: const int& dim = this->getConstraintDimension(*p_iter);
3226: if (dim) {
3227: const typename bc_type::value_type *bcArray = this->_bc->restrictPoint(*p_iter);
3229: txt << " constrained";
3230: for(int i = 0; i < dim; ++i) {
3231: txt << " " << bcArray[i];
3232: }
3233: }
3234: txt << std::endl;
3235: }
3236: }
3237: if (chart.size() == 0) {
3238: txt << "[" << this->commRank() << "]: empty" << std::endl;
3239: }
3240: PetscSynchronizedPrintf(comm, txt.str().c_str());
3241: PetscSynchronizedFlush(comm);
3242: };
3243: };
3244: // A Field combines several sections
3245: template<typename Overlap_, typename Patch_, typename Section_>
3246: class Field : public ALE::ParallelObject {
3247: public:
3248: typedef Overlap_ overlap_type;
3249: typedef Patch_ patch_type;
3250: typedef Section_ section_type;
3251: typedef typename section_type::point_type point_type;
3252: typedef typename section_type::chart_type chart_type;
3253: typedef typename section_type::value_type value_type;
3254: typedef std::map<patch_type, Obj<section_type> > sheaf_type;
3255: typedef enum {SEND, RECEIVE} request_type;
3256: typedef std::map<patch_type, MPI_Request> requests_type;
3257: protected:
3258: sheaf_type _sheaf;
3259: int _tag;
3260: MPI_Datatype _datatype;
3261: requests_type _requests;
3262: public:
3263: Field(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
3264: this->_tag = this->getNewTag();
3265: this->_datatype = this->getMPIDatatype();
3266: };
3267: Field(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _tag(tag) {
3268: this->_datatype = this->getMPIDatatype();
3269: };
3270: virtual ~Field() {};
3271: protected:
3272: MPI_Datatype getMPIDatatype() {
3273: if (sizeof(value_type) == 4) {
3274: return MPI_INT;
3275: } else if (sizeof(value_type) == 8) {
3276: return MPI_DOUBLE;
3277: } else if (sizeof(value_type) == 28) {
3278: int blen[2];
3279: MPI_Aint indices[2];
3280: MPI_Datatype oldtypes[2], newtype;
3281: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
3282: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3283: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3284: MPI_Type_commit(&newtype);
3285: return newtype;
3286: } else if (sizeof(value_type) == 32) {
3287: int blen[2];
3288: MPI_Aint indices[2];
3289: MPI_Datatype oldtypes[2], newtype;
3290: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_DOUBLE;
3291: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
3292: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
3293: MPI_Type_commit(&newtype);
3294: return newtype;
3295: }
3296: throw ALE::Exception("Cannot determine MPI type for value type");
3297: };
3298: int getNewTag() {
3299: static int tagKeyval = MPI_KEYVAL_INVALID;
3300: int *tagvalp = NULL, *maxval, flg;
3302: if (tagKeyval == MPI_KEYVAL_INVALID) {
3303: tagvalp = (int *) malloc(sizeof(int));
3304: MPI_Keyval_create(MPI_NULL_COPY_FN, DMMesh_DelTag, &tagKeyval, (void *) NULL);
3305: MPI_Attr_put(this->_comm, tagKeyval, tagvalp);
3306: tagvalp[0] = 0;
3307: }
3308: MPI_Attr_get(this->_comm, tagKeyval, (void **) &tagvalp, &flg);
3309: if (tagvalp[0] < 1) {
3310: MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
3311: tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
3312: }
3313: if (this->debug()) {
3314: std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
3315: }
3316: return tagvalp[0]--;
3317: };
3318: public: // Verifiers
3319: void checkPatch(const patch_type& patch) const {
3320: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3321: ostringstream msg;
3322: msg << "Invalid field patch " << patch << std::endl;
3323: throw ALE::Exception(msg.str().c_str());
3324: }
3325: };
3326: bool hasPatch(const patch_type& patch) {
3327: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3328: return false;
3329: }
3330: return true;
3331: };
3332: public: // Accessors
3333: int getTag() const {return this->_tag;};
3334: void setTag(const int tag) {this->_tag = tag;};
3335: Obj<section_type>& getSection(const patch_type& patch) {
3336: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
3337: this->_sheaf[patch] = new section_type(this->comm(), this->debug());
3338: }
3339: return this->_sheaf[patch];
3340: };
3341: void setSection(const patch_type& patch, const Obj<section_type>& section) {this->_sheaf[patch] = section;};
3342: const sheaf_type& getPatches() {
3343: return this->_sheaf;
3344: };
3345: void clear() {
3346: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3347: p_iter->second->clear();
3348: }
3349: };
3350: public: // Adapter
3351: template<typename Topology_>
3352: void setTopology(const Obj<Topology_>& topology) {
3353: const typename Topology_::sheaf_type& patches = topology->getPatches();
3355: for(typename Topology_::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3356: int rank = p_iter->first;
3357: const Obj<section_type>& section = this->getSection(rank);
3358: const Obj<typename Topology_::sieve_type::baseSequence>& base = p_iter->second->base();
3360: for(typename Topology_::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
3361: section->setFiberDimension(*b_iter, 1);
3362: }
3363: }
3364: }
3365: void allocate() {
3366: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3367: p_iter->second->allocatePoint();
3368: }
3369: }
3370: public: // Communication
3371: void construct(const int size) {
3372: const sheaf_type& patches = this->getPatches();
3374: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3375: const patch_type rank = p_iter->first;
3376: const Obj<section_type>& section = this->getSection(rank);
3377: const chart_type& chart = section->getChart();
3379: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3380: section->setFiberDimension(*c_iter, size);
3381: }
3382: }
3383: };
3384: template<typename Sizer>
3385: void construct(const Obj<Sizer>& sizer) {
3386: const sheaf_type& patches = this->getPatches();
3388: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3389: const patch_type rank = p_iter->first;
3390: const Obj<section_type>& section = this->getSection(rank);
3391: const chart_type& chart = section->getChart();
3393: for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
3394: section->setFiberDimension(*c_iter, *(sizer->getSection(rank)->restrictPoint(*c_iter)));
3395: }
3396: }
3397: }
3398: void constructCommunication(const request_type& requestType) {
3399: const sheaf_type& patches = this->getPatches();
3401: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3402: const patch_type patch = p_iter->first;
3403: const Obj<section_type>& section = this->getSection(patch);
3404: MPI_Request request;
3406: if (requestType == RECEIVE) {
3407: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data(" << section->size() << ") from " << patch << " tag " << this->_tag << std::endl;}
3408: MPI_Recv_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3409: } else {
3410: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << section->size() << ") to " << patch << " tag " << this->_tag << std::endl;}
3411: MPI_Send_init((void *) section->restrictSpace(), section->size(), this->_datatype, patch, this->_tag, this->comm(), &request);
3412: }
3413: this->_requests[patch] = request;
3414: }
3415: };
3416: void startCommunication() {
3417: const sheaf_type& patches = this->getPatches();
3419: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3420: MPI_Request request = this->_requests[p_iter->first];
3422: MPI_Start(&request);
3423: }
3424: };
3425: void endCommunication() {
3426: const sheaf_type& patches = this->getPatches();
3427: MPI_Status status;
3429: for(typename sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
3430: MPI_Request request = this->_requests[p_iter->first];
3432: MPI_Wait(&request, &status);
3433: }
3434: };
3435: public:
3436: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
3437: ostringstream txt;
3438: int rank;
3440: if (comm == MPI_COMM_NULL) {
3441: comm = this->comm();
3442: rank = this->commRank();
3443: } else {
3444: MPI_Comm_rank(comm, &rank);
3445: }
3446: if (name == "") {
3447: if(rank == 0) {
3448: txt << "viewing a Field" << std::endl;
3449: }
3450: } else {
3451: if(rank == 0) {
3452: txt << "viewing Field '" << name << "'" << std::endl;
3453: }
3454: }
3455: PetscSynchronizedPrintf(comm, txt.str().c_str());
3456: PetscSynchronizedFlush(comm);
3457: for(typename sheaf_type::const_iterator p_iter = this->_sheaf.begin(); p_iter != this->_sheaf.end(); ++p_iter) {
3458: ostringstream txt1;
3460: txt1 << "[" << this->commRank() << "]: Patch " << p_iter->first << std::endl;
3461: PetscSynchronizedPrintf(comm, txt1.str().c_str());
3462: PetscSynchronizedFlush(comm);
3463: p_iter->second->view("field section", comm);
3464: }
3465: };
3466: };
3467: }
3468: #endif