Actual source code: Mesh.hh
1: #ifndef included_ALE_Mesh_hh
2: #define included_ALE_Mesh_hh
4: #include <list>
5: #include <valarray>
7: #ifndef included_ALE_Numbering_hh
8: #include <Numbering.hh>
9: #endif
11: #ifndef included_ALE_INumbering_hh
12: #include <INumbering.hh>
13: #endif
15: #ifndef included_ALE_Field_hh
16: #include <Field.hh>
17: #endif
19: #ifndef included_ALE_IField_hh
20: #include <IField.hh>
21: #endif
23: #ifndef included_ALE_SieveBuilder_hh
24: #include <SieveBuilder.hh>
25: #endif
27: #ifndef included_ALE_LabelSifter_hh
28: #include <LabelSifter.hh>
29: #endif
31: #ifndef included_ALE_Partitioner_hh
32: #include <Partitioner.hh>
33: #endif
35: namespace ALE {
36: class indexSet : public std::valarray<int> {
37: public:
38: inline bool
39: operator<(const indexSet& __x) {
40: if (__x.size() != this->size()) return __x.size() < this->size();
41: for(unsigned int i = 0; i < __x.size(); ++i) {
42: if (__x[i] == (*this)[i]) continue;
43: return __x[i] < (*this)[i];
44: }
45: return false;
46: }
47: };
48: inline bool
49: operator<(const indexSet& __x, const indexSet& __y) {
50: if (__x.size() != __y.size()) return __x.size() < __y.size();
51: for(unsigned int i = 0; i < __x.size(); ++i) {
52: if (__x[i] == __y[i]) continue;
53: return __x[i] < __y[i];
54: }
55: return false;
56: };
57: inline bool
58: operator<=(const indexSet& __x, const indexSet& __y) {
59: if (__x.size() != __y.size()) return __x.size() < __y.size();
60: for(unsigned int i = 0; i < __x.size(); ++i) {
61: if (__x[i] == __y[i]) continue;
62: return __x[i] < __y[i];
63: }
64: return true;
65: };
66: inline bool
67: operator==(const indexSet& __x, const indexSet& __y) {
68: if (__x.size() != __y.size()) return false;
69: for(unsigned int i = 0; i < __x.size(); ++i) {
70: if (__x[i] != __y[i]) return false;
71: }
72: return true;
73: };
74: inline bool
75: operator!=(const indexSet& __x, const indexSet& __y) {
76: if (__x.size() != __y.size()) return true;
77: for(unsigned int i = 0; i < __x.size(); ++i) {
78: if (__x[i] != __y[i]) return true;
79: }
80: return false;
81: };
83: template<typename Sieve_,
84: typename RealSection_ = Section<typename Sieve_::point_type, double>,
85: typename IntSection_ = Section<typename Sieve_::point_type, int>,
86: typename ArrowSection_ = UniformSection<MinimalArrow<typename Sieve_::point_type, typename Sieve_::point_type>, int> >
87: class Bundle : public ALE::ParallelObject {
88: public:
89: typedef Sieve_ sieve_type;
90: typedef RealSection_ real_section_type;
91: typedef IntSection_ int_section_type;
92: typedef ArrowSection_ arrow_section_type;
93: typedef Bundle<Sieve_,RealSection_,IntSection_,ArrowSection_> this_type;
94: typedef typename sieve_type::point_type point_type;
95: typedef malloc_allocator<point_type> alloc_type;
96: typedef typename ALE::LabelSifter<int, point_type> label_type;
97: typedef typename std::map<const std::string, Obj<label_type> > labels_type;
98: typedef typename label_type::supportSequence label_sequence;
99: typedef std::map<std::string, Obj<arrow_section_type> > arrow_sections_type;
100: typedef std::map<std::string, Obj<real_section_type> > real_sections_type;
101: typedef std::map<std::string, Obj<int_section_type> > int_sections_type;
102: typedef ALE::Point index_type;
103: typedef std::pair<index_type, int> oIndex_type;
104: typedef std::vector<oIndex_type> oIndexArray;
105: typedef std::pair<int *, int> indices_type;
106: typedef NumberingFactory<this_type> MeshNumberingFactory;
107: typedef typename ALE::Partitioner<>::part_type rank_type;
108: typedef typename ALE::Sifter<point_type,rank_type,point_type> send_overlap_type;
109: typedef typename ALE::Sifter<rank_type,point_type,point_type> recv_overlap_type;
110: typedef typename MeshNumberingFactory::numbering_type numbering_type;
111: typedef typename MeshNumberingFactory::order_type order_type;
112: typedef std::map<point_type, point_type> renumbering_type;
113: typedef typename ALE::SieveAlg<this_type> sieve_alg_type;
114: typedef typename sieve_alg_type::coneArray coneArray;
115: typedef typename sieve_alg_type::orientedConeArray oConeArray;
116: typedef typename sieve_alg_type::supportArray supportArray;
117: protected:
118: Obj<sieve_type> _sieve;
119: labels_type _labels;
120: int _maxHeight;
121: int _maxDepth;
122: arrow_sections_type _arrowSections;
123: real_sections_type _realSections;
124: int_sections_type _intSections;
125: Obj<oIndexArray> _indexArray;
126: Obj<MeshNumberingFactory> _factory;
127: bool _calculatedOverlap;
128: Obj<send_overlap_type> _sendOverlap;
129: Obj<recv_overlap_type> _recvOverlap;
130: Obj<send_overlap_type> _distSendOverlap;
131: Obj<recv_overlap_type> _distRecvOverlap;
132: renumbering_type _renumbering; // Maps global points to local points
133: // Work space
134: Obj<std::set<point_type> > _modifiedPoints;
135: public:
136: Bundle(MPI_Comm comm, int debug = 0) : ALE::ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1) {
137: this->_indexArray = new oIndexArray();
138: this->_modifiedPoints = new std::set<point_type>();
139: this->_factory = MeshNumberingFactory::singleton(this->comm(), this->debug());
140: this->_calculatedOverlap = false;
141: this->_sendOverlap = new send_overlap_type(comm, debug);
142: this->_recvOverlap = new recv_overlap_type(comm, debug);
143: };
144: Bundle(const Obj<sieve_type>& sieve) : ALE::ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _maxHeight(-1), _maxDepth(-1) {
145: this->_indexArray = new oIndexArray();
146: this->_modifiedPoints = new std::set<point_type>();
147: this->_factory = MeshNumberingFactory::singleton(this->comm(), this->debug());
148: this->_calculatedOverlap = false;
149: this->_sendOverlap = new send_overlap_type(comm, debug);
150: this->_recvOverlap = new recv_overlap_type(comm, debug);
151: };
152: virtual ~Bundle() {};
153: public: // Verifiers
154: bool hasLabel(const std::string& name) {
155: if (this->_labels.find(name) != this->_labels.end()) {
156: return true;
157: }
158: return false;
159: };
160: void checkLabel(const std::string& name) {
161: if (!this->hasLabel(name)) {
162: ostringstream msg;
163: msg << "Invalid label name: " << name << std::endl;
164: throw ALE::Exception(msg.str().c_str());
165: }
166: };
167: public: // Accessors
168: const Obj<sieve_type>& getSieve() const {return this->_sieve;};
169: void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
170: bool hasArrowSection(const std::string& name) const {
171: return this->_arrowSections.find(name) != this->_arrowSections.end();
172: };
173: const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
174: if (!this->hasArrowSection(name)) {
175: Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());
177: section->setName(name);
178: if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
179: this->_arrowSections[name] = section;
180: }
181: return this->_arrowSections[name];
182: };
183: void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
184: this->_arrowSections[name] = section;
185: };
186: Obj<std::set<std::string> > getArrowSections() const {
187: Obj<std::set<std::string> > names = std::set<std::string>();
189: for(typename arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
190: names->insert(s_iter->first);
191: }
192: return names;
193: };
194: bool hasRealSection(const std::string& name) const {
195: return this->_realSections.find(name) != this->_realSections.end();
196: };
197: const Obj<real_section_type>& getRealSection(const std::string& name) {
198: if (!this->hasRealSection(name)) {
199: Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());
201: section->setName(name);
202: if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
203: this->_realSections[name] = section;
204: }
205: return this->_realSections[name];
206: };
207: void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
208: this->_realSections[name] = section;
209: };
210: Obj<std::set<std::string> > getRealSections() const {
211: Obj<std::set<std::string> > names = std::set<std::string>();
213: for(typename real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
214: names->insert(s_iter->first);
215: }
216: return names;
217: };
218: bool hasIntSection(const std::string& name) const {
219: return this->_intSections.find(name) != this->_intSections.end();
220: };
221: const Obj<int_section_type>& getIntSection(const std::string& name) {
222: if (!this->hasIntSection(name)) {
223: Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());
225: section->setName(name);
226: if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
227: this->_intSections[name] = section;
228: }
229: return this->_intSections[name];
230: };
231: void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
232: this->_intSections[name] = section;
233: };
234: Obj<std::set<std::string> > getIntSections() const {
235: Obj<std::set<std::string> > names = std::set<std::string>();
237: for(typename int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
238: names->insert(s_iter->first);
239: }
240: return names;
241: };
242: const Obj<MeshNumberingFactory>& getFactory() const {return this->_factory;};
243: bool getCalculatedOverlap() const {return this->_calculatedOverlap;};
244: void setCalculatedOverlap(const bool calc) {this->_calculatedOverlap = calc;};
245: const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
246: void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
247: const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
248: void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
249: const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
250: void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
251: const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
252: void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
253: renumbering_type& getRenumbering() {return this->_renumbering;};
254: public: // Labels
255: int getValue (const Obj<label_type>& label, const point_type& point, const int defValue = 0) {
256: const Obj<typename label_type::coneSequence>& cone = label->cone(point);
258: if (cone->size() == 0) return defValue;
259: return *cone->begin();
260: };
261: void setValue(const Obj<label_type>& label, const point_type& point, const int value) {
262: label->setCone(value, point);
263: };
264: template<typename InputPoints>
265: int getMaxValue (const Obj<label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
266: int maxValue = defValue;
268: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
269: maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
270: }
271: return maxValue;
272: };
273: const Obj<label_type>& createLabel(const std::string& name) {
274: this->_labels[name] = new label_type(this->comm(), this->debug());
275: return this->_labels[name];
276: };
277: const Obj<label_type>& getLabel(const std::string& name) {
278: this->checkLabel(name);
279: return this->_labels[name];
280: };
281: void setLabel(const std::string& name, const Obj<label_type>& label) {
282: this->_labels[name] = label;
283: };
284: const labels_type& getLabels() {
285: return this->_labels;
286: };
287: virtual const Obj<label_sequence>& getLabelStratum(const std::string& name, int value) {
288: this->checkLabel(name);
289: return this->_labels[name]->support(value);
290: };
291: public: // Stratification
292: template<class InputPoints>
293: void computeHeight(const Obj<label_type>& height, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxHeight) {
294: this->_modifiedPoints->clear();
296: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
297: // Compute the max height of the points in the support of p, and add 1
298: int h0 = this->getValue(height, *p_iter, -1);
299: int h1 = this->getMaxValue(height, sieve->support(*p_iter), -1) + 1;
301: if(h1 != h0) {
302: this->setValue(height, *p_iter, h1);
303: if (h1 > maxHeight) maxHeight = h1;
304: this->_modifiedPoints->insert(*p_iter);
305: }
306: }
307: // FIX: We would like to avoid the copy here with cone()
308: if(this->_modifiedPoints->size() > 0) {
309: this->computeHeight(height, sieve, sieve->cone(this->_modifiedPoints), maxHeight);
310: }
311: };
312: void computeHeights() {
313: const Obj<label_type>& label = this->createLabel(std::string("height"));
315: this->_maxHeight = -1;
316: this->computeHeight(label, this->_sieve, this->_sieve->leaves(), this->_maxHeight);
317: };
318: virtual int height() const {return this->_maxHeight;};
319: virtual int height(const point_type& point) {
320: return this->getValue(this->_labels["height"], point, -1);
321: };
322: virtual const Obj<label_sequence>& heightStratum(int height) {
323: return this->getLabelStratum("height", height);
324: };
325: void setHeight(const Obj<label_type>& label) {
326: this->_labels["height"] = label;
327: const Obj<typename label_type::traits::capSequence> cap = label->cap();
329: for(typename label_type::traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
330: this->_maxHeight = std::max(this->_maxHeight, *c_iter);
331: }
332: };
333: template<class InputPoints>
334: void computeDepth(const Obj<label_type>& depth, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxDepth) {
335: this->_modifiedPoints->clear();
337: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
338: // Compute the max depth of the points in the cone of p, and add 1
339: int d0 = this->getValue(depth, *p_iter, -1);
340: int d1 = this->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;
342: if(d1 != d0) {
343: this->setValue(depth, *p_iter, d1);
344: if (d1 > maxDepth) maxDepth = d1;
345: this->_modifiedPoints->insert(*p_iter);
346: }
347: }
348: // FIX: We would like to avoid the copy here with support()
349: if(this->_modifiedPoints->size() > 0) {
350: this->computeDepth(depth, sieve, sieve->support(this->_modifiedPoints), maxDepth);
351: }
352: };
353: void computeDepths() {
354: const Obj<label_type>& label = this->createLabel(std::string("depth"));
356: this->_maxDepth = -1;
357: this->computeDepth(label, this->_sieve, this->_sieve->roots(), this->_maxDepth);
358: };
359: virtual int depth() const {return this->_maxDepth;};
360: virtual int depth(const point_type& point) {
361: return this->getValue(this->_labels["depth"], point, -1);
362: };
363: virtual const Obj<label_sequence>& depthStratum(int depth) {
364: return this->getLabelStratum("depth", depth);
365: };
368: virtual void stratify() {
369: ALE_LOG_EVENT_BEGIN;
370: this->computeHeights();
371: this->computeDepths();
372: ALE_LOG_EVENT_END;
373: };
374: public: // Size traversal
375: template<typename Section_>
376: int size(const Obj<Section_>& section, const point_type& p) {
377: const typename Section_::chart_type& chart = section->getChart();
378: int size = 0;
380: if (this->height() < 2) {
381: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
382: typename sieve_type::coneSequence::iterator end = cone->end();
384: if (chart.count(p)) {
385: size += section->getConstrainedFiberDimension(p);
386: }
387: for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
388: if (chart.count(*c_iter)) {
389: size += section->getConstrainedFiberDimension(*c_iter);
390: }
391: }
392: } else {
393: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), p);
394: typename coneArray::iterator end = closure->end();
396: for(typename coneArray::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
397: if (chart.count(*c_iter)) {
398: size += section->getConstrainedFiberDimension(*c_iter);
399: }
400: }
401: }
402: return size;
403: };
404: template<typename Section_>
405: int sizeWithBC(const Obj<Section_>& section, const point_type& p) {
406: const typename Section_::chart_type& chart = section->getChart();
407: int size = 0;
409: if (this->height() < 2) {
410: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
411: typename sieve_type::coneSequence::iterator end = cone->end();
413: if (chart.count(p)) {
414: size += section->getFiberDimension(p);
415: }
416: for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
417: if (chart.count(*c_iter)) {
418: size += section->getFiberDimension(*c_iter);
419: }
420: }
421: } else {
422: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), p);
423: typename coneArray::iterator end = closure->end();
425: for(typename coneArray::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
426: if (chart.count(*c_iter)) {
427: size += section->getFiberDimension(*c_iter);
428: }
429: }
430: }
431: return size;
432: };
433: protected:
434: int *getIndexArray(const int size) {
435: static int *array = NULL;
436: static int maxSize = 0;
438: if (size > maxSize) {
439: maxSize = size;
440: if (array) delete [] array;
441: array = new int[maxSize];
442: };
443: return array;
444: };
445: public: // Index traversal
446: void expandInterval(const index_type& interval, PetscInt indices[], PetscInt *indx) {
447: const int end = interval.prefix + interval.index;
449: for(int i = interval.index; i < end; ++i) {
450: indices[(*indx)++] = i;
451: }
452: };
453: void expandInterval(const index_type& interval, const int orientation, PetscInt indices[], PetscInt *indx) {
454: if (orientation >= 0) {
455: for(int i = 0; i < interval.prefix; ++i) {
456: indices[(*indx)++] = interval.index + i;
457: }
458: } else {
459: for(int i = interval.prefix-1; i >= 0; --i) {
460: indices[(*indx)++] = interval.index + i;
461: }
462: }
463: for(int i = 0; i < -interval.prefix; ++i) {
464: indices[(*indx)++] = -1;
465: }
466: };
467: void expandIntervals(Obj<oIndexArray> intervals, PetscInt *indices) {
468: int k = 0;
470: for(typename oIndexArray::iterator i_iter = intervals->begin(); i_iter != intervals->end(); ++i_iter) {
471: this->expandInterval(i_iter->first, i_iter->second, indices, &k);
472: }
473: }
474: template<typename Section_>
475: const indices_type getIndicesRaw(const Obj<Section_>& section, const point_type& p) {
476: int *indexArray = NULL;
477: int size = 0;
479: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
480: typename oConeArray::iterator begin = closure->begin();
481: typename oConeArray::iterator end = closure->end();
483: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
484: size += section->getFiberDimension(p_iter->first);
485: }
486: indexArray = this->getIndexArray(size);
487: int k = 0;
488: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
489: section->getIndicesRaw(p_iter->first, section->getIndex(p_iter->first), indexArray, &k, p_iter->second);
490: }
491: return indices_type(indexArray, size);
492: };
493: template<typename Section_>
494: const indices_type getIndices(const Obj<Section_>& section, const point_type& p, const int level = -1) {
495: int *indexArray = NULL;
496: int size = 0;
498: if (level == 0) {
499: size += section->getFiberDimension(p);
500: indexArray = this->getIndexArray(size);
501: int k = 0;
503: section->getIndices(p, indexArray, &k);
504: } else if ((level == 1) || (this->height() == 1)) {
505: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
506: typename sieve_type::coneSequence::iterator end = cone->end();
508: size += section->getFiberDimension(p);
509: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
510: size += section->getFiberDimension(*p_iter);
511: }
512: indexArray = this->getIndexArray(size);
513: int k = 0;
515: section->getIndices(p, indexArray, &k);
516: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
517: section->getIndices(*p_iter, indexArray, &k);
518: }
519: } else if (level == -1) {
520: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
521: typename oConeArray::iterator begin = closure->begin();
522: typename oConeArray::iterator end = closure->end();
524: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
525: size += section->getFiberDimension(p_iter->first);
526: }
527: indexArray = this->getIndexArray(size);
528: int k = 0;
529: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
530: section->getIndices(p_iter->first, indexArray, &k, p_iter->second);
531: }
532: } else {
533: throw ALE::Exception("Bundle has not yet implemented getIndices() for an arbitrary level");
534: }
535: if (this->debug()) {
536: for(int i = 0; i < size; ++i) {
537: printf("[%d]index %d: %d\n", this->commRank(), i, indexArray[i]);
538: }
539: }
540: return indices_type(indexArray, size);
541: };
542: template<typename Section_, typename Numbering>
543: const indices_type getIndices(const Obj<Section_>& section, const point_type& p, const Obj<Numbering>& numbering, const int level = -1) {
544: int *indexArray = NULL;
545: int size = 0;
547: if (level == 0) {
548: size += section->getFiberDimension(p);
549: indexArray = this->getIndexArray(size);
550: int k = 0;
552: section->getIndices(p, numbering, indexArray, &k);
553: } else if ((level == 1) || (this->height() == 1)) {
554: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
555: typename sieve_type::coneSequence::iterator end = cone->end();
557: size += section->getFiberDimension(p);
558: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
559: size += section->getFiberDimension(*p_iter);
560: }
561: indexArray = this->getIndexArray(size);
562: int k = 0;
564: section->getIndices(p, numbering, indexArray, &k);
565: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
566: section->getIndices(*p_iter, numbering, indexArray, &k);
567: }
568: } else if (level == -1) {
569: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
570: typename oConeArray::iterator end = closure->end();
572: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
573: size += section->getFiberDimension(p_iter->first);
574: }
575: indexArray = this->getIndexArray(size);
576: int k = 0;
577: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
578: section->getIndices(p_iter->first, numbering, indexArray, &k, p_iter->second);
579: }
580: } else {
581: throw ALE::Exception("Bundle has not yet implemented getIndices() for an arbitrary level");
582: }
583: return indices_type(indexArray, size);
584: };
585: public: // Retrieval traversal
586: // Return the values for the closure of this point
587: // use a smart pointer?
588: template<typename Section_>
589: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p) {
590: const int size = this->sizeWithBC(section, p);
591: return this->restrictClosure(section, p, section->getRawArray(size), size);
592: };
593: template<typename Section_>
594: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p, typename Section_::value_type *values, const int valuesSize) {
595: const int size = this->sizeWithBC(section, p);
596: int j = -1;
597: if (valuesSize < size) throw ALE::Exception("Input array too small");
599: // We could actually ask for the height of the individual point
600: if (this->height() < 2) {
601: const int& dim = section->getFiberDimension(p);
602: const typename Section_::value_type *array = section->restrictPoint(p);
604: for(int i = 0; i < dim; ++i) {
605: values[++j] = array[i];
606: }
607: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
608: typename sieve_type::coneSequence::iterator end = cone->end();
610: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
611: const int& dim = section->getFiberDimension(*p_iter);
613: array = section->restrictPoint(*p_iter);
614: for(int i = 0; i < dim; ++i) {
615: values[++j] = array[i];
616: }
617: }
618: } else {
619: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
620: typename oConeArray::iterator end = closure->end();
622: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
623: const typename Section_::value_type *array = section->restrictPoint(p_iter->first);
624: const int& dim = section->getFiberDimension(p_iter->first);
626: if (p_iter->second >= 0) {
627: for(int i = 0; i < dim; ++i) {
628: values[++j] = array[i];
629: }
630: } else {
631: for(int i = dim-1; i >= 0; --i) {
632: values[++j] = array[i];
633: }
634: }
635: }
636: }
637: if (j != size-1) {
638: ostringstream txt;
640: txt << "Invalid restrict to point " << p << std::endl;
641: txt << " j " << j << " should be " << (size-1) << std::endl;
642: std::cout << txt.str();
643: throw ALE::Exception(txt.str().c_str());
644: }
645: return values;
646: };
647: template<typename Section_>
648: const typename Section_::value_type *restrictNew(const Obj<Section_>& section, const point_type& p) {
649: const int size = this->sizeWithBC(section, p);
650: return this->restrictNew(section, p, section->getRawArray(size), size);
651: };
652: template<typename Section_>
653: const typename Section_::value_type *restrictNew(const Obj<Section_>& section, const point_type& p, typename Section_::value_type *values, const int valuesSize) {
654: const int size = this->sizeWithBC(section, p);
655: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
656: typename oConeArray::iterator end = closure->end();
657: int j = -1;
658: if (valuesSize < size) throw ALE::Exception("Input array too small");
660: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
661: const typename Section_::value_type *array = section->restrictPoint(p_iter->first);
663: if (p_iter->second >= 0) {
664: const int& dim = section->getFiberDimension(p_iter->first);
666: for(int i = 0; i < dim; ++i) {
667: values[++j] = array[i];
668: }
669: } else {
670: int offset = 0;
672: for(int space = 0; space < section->getNumSpaces(); ++space) {
673: const int& dim = section->getFiberDimension(p_iter->first, space);
675: for(int i = dim-1; i >= 0; --i) {
676: values[++j] = array[i+offset];
677: }
678: offset += dim;
679: }
680: }
681: }
682: if (j != size-1) {
683: ostringstream txt;
685: txt << "Invalid restrict to point " << p << std::endl;
686: txt << " j " << j << " should be " << (size-1) << std::endl;
687: std::cout << txt.str();
688: throw ALE::Exception(txt.str().c_str());
689: }
690: return values;
691: };
692: template<typename Section_>
693: void update(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
694: int j = 0;
696: if (this->height() < 2) {
697: section->updatePoint(p, &v[j]);
698: j += section->getFiberDimension(p);
699: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
701: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
702: section->updatePoint(*p_iter, &v[j]);
703: j += section->getFiberDimension(*p_iter);
704: }
705: } else {
706: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
707: typename oConeArray::iterator end = closure->end();
709: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
710: section->updatePoint(p_iter->first, &v[j], p_iter->second);
711: j += section->getFiberDimension(p_iter->first);
712: }
713: }
714: };
715: template<typename Section_>
716: void updateAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
717: int j = 0;
719: if (this->height() < 2) {
720: section->updateAddPoint(p, &v[j]);
721: j += section->getFiberDimension(p);
722: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
724: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
725: section->updateAddPoint(*p_iter, &v[j]);
726: j += section->getFiberDimension(*p_iter);
727: }
728: } else {
729: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
730: typename oConeArray::iterator end = closure->end();
732: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
733: section->updateAddPoint(p_iter->first, &v[j], p_iter->second);
734: j += section->getFiberDimension(p_iter->first);
735: }
736: }
737: };
738: template<typename Section_>
739: void updateBC(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
740: int j = 0;
742: if (this->height() < 2) {
743: section->updatePointBC(p, &v[j]);
744: j += section->getFiberDimension(p);
745: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
747: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
748: section->updatePointBC(*p_iter, &v[j]);
749: j += section->getFiberDimension(*p_iter);
750: }
751: } else {
752: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
753: typename oConeArray::iterator end = closure->end();
755: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
756: section->updatePointBC(p_iter->first, &v[j], p_iter->second);
757: j += section->getFiberDimension(p_iter->first);
758: }
759: }
760: };
761: template<typename Section_>
762: void updateAll(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
763: int j = 0;
765: if (this->height() < 2) {
766: section->updatePointAll(p, &v[j]);
767: j += section->getFiberDimension(p);
768: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
770: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
771: section->updatePointAll(*p_iter, &v[j]);
772: j += section->getFiberDimension(*p_iter);
773: }
774: } else {
775: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
776: typename oConeArray::iterator end = closure->end();
778: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
779: section->updatePointAll(p_iter->first, &v[j], p_iter->second);
780: j += section->getFiberDimension(p_iter->first);
781: }
782: }
783: };
784: template<typename Section_>
785: void updateAllAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
786: int j = 0;
788: if (this->height() < 2) {
789: section->updatePointAllAdd(p, &v[j]);
790: j += section->getFiberDimension(p);
791: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
793: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
794: section->updatePointAllAdd(*p_iter, &v[j]);
795: j += section->getFiberDimension(*p_iter);
796: }
797: } else {
798: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
799: typename oConeArray::iterator end = closure->end();
801: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
802: section->updatePointAllAdd(p_iter->first, &v[j], p_iter->second);
803: j += section->getFiberDimension(p_iter->first);
804: }
805: }
806: };
807: public: // Optimization
808: // Calculate a custom atlas for the given traversal
809: // This returns the tag value assigned to the traversal
810: template<typename Section_, typename Sequence_>
811: int calculateCustomAtlas(const Obj<Section_>& section, const Obj<Sequence_>& points) {
812: const typename Sequence_::iterator begin = points->begin();
813: const typename Sequence_::iterator end = points->end();
814: const int num = points->size();
815: int *rOffsets = new int[num+1];
816: int *rIndices;
817: int *uOffsets = new int[num+1];
818: int *uIndices;
819: int p;
821: p = 0;
822: rOffsets[p] = 0;
823: uOffsets[p] = 0;
824: for(typename Sequence_::iterator p_iter = begin; p_iter != end; ++p_iter, ++p) {
825: rOffsets[p+1] = rOffsets[p] + this->sizeWithBC(section, *p_iter);
826: uOffsets[p+1] = rOffsets[p+1];
827: //uOffsets[p+1] = uOffsets[p] + this->size(section, *p_iter);
828: }
829: rIndices = new int[rOffsets[p]];
830: uIndices = new int[uOffsets[p]];
831: p = 0;
832: for(typename Sequence_::iterator p_iter = begin; p_iter != end; ++p_iter, ++p) {
833: const indices_type rIdx = this->getIndicesRaw(section, *p_iter);
834: for(int i = 0, k = rOffsets[p]; k < rOffsets[p+1]; ++i, ++k) rIndices[k] = rIdx.first[i];
836: const indices_type uIdx = this->getIndices(section, *p_iter);
837: for(int i = 0, k = uOffsets[p]; k < uOffsets[p+1]; ++i, ++k) uIndices[k] = uIdx.first[i];
838: }
839: return section->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices);
840: };
841: template<typename Section_>
842: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p) {
843: const int *offsets, *indices;
845: section->getCustomRestrictAtlas(tag, &offsets, &indices);
846: const int size = offsets[p+1] - offsets[p];
847: return this->restrictClosure(section, tag, p, section->getRawArray(size), offsets, indices);
848: };
849: template<typename Section_>
850: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p, typename Section_::value_type *values, const int valuesSize) {
851: const int *offsets, *indices;
853: section->getCustomRestrictAtlas(tag, &offsets, &indices);
854: const int size = offsets[p+1] - offsets[p];
855: if (valuesSize < size) {throw ALE::Exception("Input array too small");}
856: return this->restrictClosure(section, tag, p, values, offsets, indices);
857: };
858: template<typename Section_>
859: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p, typename Section_::value_type *values, const int offsets[], const int indices[]) {
860: const typename Section_::value_type *array = section->restrictSpace();
862: const int size = offsets[p+1] - offsets[p];
863: for(int j = 0, k = offsets[p]; j < size; ++j, ++k) {
864: values[j] = array[indices[k]];
865: }
866: return values;
867: };
868: template<typename Section_>
869: void updateAdd(const Obj<Section_>& section, const int tag, const int p, const typename Section_::value_type values[]) {
870: typename Section_::value_type *array = (typename Section_::value_type *) section->restrictSpace();
871: const int *offsets, *indices;
873: section->getCustomUpdateAtlas(tag, &offsets, &indices);
874: const int size = offsets[p+1] - offsets[p];
875: for(int j = 0, k = offsets[p]; j < size; ++j, ++k) {
876: if (indices[k] < 0) continue;
877: array[indices[k]] += values[j];
878: }
879: };
880: public: // Allocation
881: template<typename Section_>
882: void allocate(const Obj<Section_>& section, const Obj<send_overlap_type>& sendOverlap = NULL) {
883: bool doGhosts = !sendOverlap.isNull();
885: this->_factory->orderPatch(section, this->getSieve(), sendOverlap);
886: if (doGhosts) {
887: if (this->_debug > 1) {std::cout << "Ordering patch for ghosts" << std::endl;}
888: const typename Section_::chart_type& points = section->getChart();
889: typename Section_::index_type::index_type offset = 0;
891: for(typename Section_::chart_type::const_iterator point = points.begin(); point != points.end(); ++point) {
892: const typename Section_::index_type& idx = section->getIndex(*point);
894: offset = std::max(offset, idx.index + std::abs(idx.prefix));
895: }
896: this->_factory->orderPatch(section, this->getSieve(), NULL, offset);
897: if (offset != section->sizeWithBC()) throw ALE::Exception("Inconsistent array sizes in section");
898: }
899: section->allocateStorage();
900: };
901: template<typename Section_>
902: void reallocate(const Obj<Section_>& section) {
903: if (section->getNewAtlas().isNull()) return;
904: // Since copy() preserves offsets, we must reinitialize them before ordering
905: const Obj<typename Section_::atlas_type> atlas = section->getAtlas();
906: const Obj<typename Section_::atlas_type>& newAtlas = section->getNewAtlas();
907: const typename Section_::atlas_type::chart_type& chart = newAtlas->getChart();
908: const typename Section_::atlas_type::chart_type& oldChart = atlas->getChart();
909: int newSize = 0;
910: typename Section_::index_type defaultIdx(0, -1);
912: for(typename Section_::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
913: defaultIdx.prefix = newAtlas->restrictPoint(*c_iter)[0].prefix;
914: newAtlas->updatePoint(*c_iter, &defaultIdx);
915: newSize += defaultIdx.prefix;
916: }
917: section->setAtlas(newAtlas);
918: this->_factory->orderPatch(section, this->getSieve());
919: // Copy over existing values
920: typedef typename alloc_type::template rebind<typename Section_::value_type>::other value_alloc_type;
921: value_alloc_type value_allocator;
922: typename Section_::value_type *newArray = value_allocator.allocate(newSize);
923: for(int i = 0; i < newSize; ++i) {value_allocator.construct(newArray+i, typename Section_::value_type());}
924: ///typename Section_::value_type *newArray = new typename Section_::value_type[newSize];
925: const typename Section_::value_type *array = section->restrictSpace();
927: for(typename Section_::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
928: const int& dim = section->getFiberDimension(*c_iter);
929: const int& offset = atlas->restrictPoint(*c_iter)->index;
930: const int& newOffset = newAtlas->restrictPoint(*c_iter)->index;
932: for(int i = 0; i < dim; ++i) {
933: newArray[newOffset+i] = array[offset+i];
934: }
935: }
936: section->replaceStorage(newArray);
937: };
938: public: // Overlap
939: template<typename Sequence>
940: void constructOverlap(const Obj<Sequence>& points, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
941: point_type *sendBuf = new point_type[points->size()];
942: int size = 0;
943: for(typename Sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
944: sendBuf[size++] = *l_iter;
945: }
946: int *sizes = new int[this->commSize()]; // The number of points coming from each process
947: int *offsets = new int[this->commSize()+1]; // Prefix sums for sizes
948: int *oldOffs = new int[this->commSize()+1]; // Temporary storage
949: point_type *remotePoints = NULL; // The points from each process
950: int *remoteRanks = NULL; // The rank and number of overlap points of each process that overlaps another
952: // Change to Allgather() for the correct binning algorithm
953: MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, this->comm());
954: if (this->commRank() == 0) {
955: offsets[0] = 0;
956: for(int p = 1; p <= this->commSize(); p++) {
957: offsets[p] = offsets[p-1] + sizes[p-1];
958: }
959: remotePoints = new point_type[offsets[this->commSize()]];
960: }
961: MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, this->comm());
962: delete [] sendBuf;
963: std::map<int, std::map<int, std::set<point_type> > > overlapInfo; // Maps (p,q) to their set of overlap points
965: if (this->commRank() == 0) {
966: for(int p = 0; p < this->commSize(); p++) {
967: std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]]);
968: }
969: for(int p = 0; p <= this->commSize(); p++) {
970: oldOffs[p] = offsets[p];
971: }
972: for(int p = 0; p < this->commSize(); p++) {
973: for(int q = p+1; q < this->commSize(); q++) {
974: std::set_intersection(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
975: &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
976: std::insert_iterator<std::set<point_type> >(overlapInfo[p][q], overlapInfo[p][q].begin()));
977: overlapInfo[q][p] = overlapInfo[p][q];
978: }
979: sizes[p] = overlapInfo[p].size()*2;
980: offsets[p+1] = offsets[p] + sizes[p];
981: }
982: remoteRanks = new int[offsets[this->commSize()]];
983: int k = 0;
984: for(int p = 0; p < this->commSize(); p++) {
985: for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
986: remoteRanks[k*2] = r_iter->first;
987: remoteRanks[k*2+1] = r_iter->second.size();
988: k++;
989: }
990: }
991: }
992: int numOverlaps; // The number of processes overlapping this process
993: MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, this->comm());
994: int *overlapRanks = new int[numOverlaps]; // The rank and overlap size for each overlapping process
995: MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, this->comm());
996: point_type *sendPoints = NULL; // The points to send to each process
997: if (this->commRank() == 0) {
998: for(int p = 0, k = 0; p < this->commSize(); p++) {
999: sizes[p] = 0;
1000: for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
1001: sizes[p] += remoteRanks[k*2+1];
1002: k++;
1003: }
1004: offsets[p+1] = offsets[p] + sizes[p];
1005: }
1006: sendPoints = new point_type[offsets[this->commSize()]];
1007: for(int p = 0, k = 0; p < this->commSize(); p++) {
1008: for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
1009: int rank = r_iter->first;
1010: for(typename std::set<point_type>::iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
1011: sendPoints[k++] = *p_iter;
1012: }
1013: }
1014: }
1015: }
1016: int numOverlapPoints = 0;
1017: for(int r = 0; r < numOverlaps/2; r++) {
1018: numOverlapPoints += overlapRanks[r*2+1];
1019: }
1020: point_type *overlapPoints = new point_type[numOverlapPoints];
1021: MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints, MPI_INT, 0, this->comm());
1023: for(int r = 0, k = 0; r < numOverlaps/2; r++) {
1024: int rank = overlapRanks[r*2];
1026: for(int p = 0; p < overlapRanks[r*2+1]; p++) {
1027: point_type point = overlapPoints[k++];
1029: sendOverlap->addArrow(point, rank, point);
1030: recvOverlap->addArrow(rank, point, point);
1031: }
1032: }
1034: delete [] overlapPoints;
1035: delete [] overlapRanks;
1036: delete [] sizes;
1037: delete [] offsets;
1038: delete [] oldOffs;
1039: if (this->commRank() == 0) {
1040: delete [] remoteRanks;
1041: delete [] remotePoints;
1042: delete [] sendPoints;
1043: }
1044: };
1045: void constructOverlap() {
1046: if (this->_calculatedOverlap) return;
1047: this->constructOverlap(this->getSieve()->base(), this->getSendOverlap(), this->getRecvOverlap());
1048: this->constructOverlap(this->getSieve()->cap(), this->getSendOverlap(), this->getRecvOverlap());
1049: if (this->debug()) {
1050: this->_sendOverlap->view("Send overlap");
1051: this->_recvOverlap->view("Receive overlap");
1052: }
1053: this->_calculatedOverlap = true;
1054: };
1055: };
1056: class BoundaryCondition : public ALE::ParallelObject {
1057: public:
1058: typedef double (*function_type)(const double []);
1059: typedef double (*integrator_type)(const double [], const double [], const int, function_type);
1060: protected:
1061: std::string _labelName;
1062: int _marker;
1063: function_type _func;
1064: integrator_type _integrator;
1065: public:
1066: BoundaryCondition(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
1067: ~BoundaryCondition() {};
1068: public:
1069: const std::string& getLabelName() const {return this->_labelName;};
1070: void setLabelName(const std::string& name) {this->_labelName = name;};
1071: int getMarker() const {return this->_marker;};
1072: void setMarker(const int marker) {this->_marker = marker;};
1073: function_type getFunction() const {return this->_func;};
1074: void setFunction(function_type func) {this->_func = func;};
1075: integrator_type getDualIntegrator() const {return this->_integrator;};
1076: void setDualIntegrator(integrator_type integrator) {this->_integrator = integrator;};
1077: public:
1078: double evaluate(const double coords[]) const {return this->_func(coords);};
1079: double integrateDual(const double v0[], const double J[], const int dualIndex) const {return this->_integrator(v0, J, dualIndex, this->_func);};
1080: };
1081: class Discretization : public ALE::ParallelObject {
1082: typedef std::map<std::string, Obj<BoundaryCondition> > boundaryConditions_type;
1083: protected:
1084: boundaryConditions_type _boundaryConditions;
1085: Obj<BoundaryCondition> _exactSolution;
1086: std::map<int,int> _dim2dof;
1087: std::map<int,int> _dim2class;
1088: int _quadSize;
1089: const double *_points;
1090: const double *_weights;
1091: int _basisSize;
1092: const double *_basis;
1093: const double *_basisDer;
1094: const int *_indices;
1095: std::map<int, const int *> _exclusionIndices;
1096: public:
1097: Discretization(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _quadSize(0), _points(NULL), _weights(NULL), _basisSize(0), _basis(NULL), _basisDer(NULL), _indices(NULL) {};
1098: virtual ~Discretization() {
1099: if (this->_indices) {delete [] this->_indices;}
1100: for(std::map<int, const int *>::iterator i_iter = _exclusionIndices.begin(); i_iter != _exclusionIndices.end(); ++i_iter) {
1101: delete [] i_iter->second;
1102: }
1103: };
1104: public:
1105: const bool hasBoundaryCondition() {return (this->_boundaryConditions.find("default") != this->_boundaryConditions.end());};
1106: const Obj<BoundaryCondition>& getBoundaryCondition() {return this->getBoundaryCondition("default");};
1107: void setBoundaryCondition(const Obj<BoundaryCondition>& boundaryCondition) {this->setBoundaryCondition("default", boundaryCondition);};
1108: const Obj<BoundaryCondition>& getBoundaryCondition(const std::string& name) {return this->_boundaryConditions[name];};
1109: void setBoundaryCondition(const std::string& name, const Obj<BoundaryCondition>& boundaryCondition) {this->_boundaryConditions[name] = boundaryCondition;};
1110: Obj<std::set<std::string> > getBoundaryConditions() const {
1111: Obj<std::set<std::string> > names = std::set<std::string>();
1113: for(boundaryConditions_type::const_iterator d_iter = this->_boundaryConditions.begin(); d_iter != this->_boundaryConditions.end(); ++d_iter) {
1114: names->insert(d_iter->first);
1115: }
1116: return names;
1117: };
1118: const Obj<BoundaryCondition>& getExactSolution() {return this->_exactSolution;};
1119: void setExactSolution(const Obj<BoundaryCondition>& exactSolution) {this->_exactSolution = exactSolution;};
1120: const int getQuadratureSize() {return this->_quadSize;};
1121: void setQuadratureSize(const int size) {this->_quadSize = size;};
1122: const double *getQuadraturePoints() {return this->_points;};
1123: void setQuadraturePoints(const double *points) {this->_points = points;};
1124: const double *getQuadratureWeights() {return this->_weights;};
1125: void setQuadratureWeights(const double *weights) {this->_weights = weights;};
1126: const int getBasisSize() {return this->_basisSize;};
1127: void setBasisSize(const int size) {this->_basisSize = size;};
1128: const double *getBasis() {return this->_basis;};
1129: void setBasis(const double *basis) {this->_basis = basis;};
1130: const double *getBasisDerivatives() {return this->_basisDer;};
1131: void setBasisDerivatives(const double *basisDer) {this->_basisDer = basisDer;};
1132: int getNumDof(const int dim) {return this->_dim2dof[dim];};
1133: void setNumDof(const int dim, const int numDof) {this->_dim2dof[dim] = numDof;};
1134: int getDofClass(const int dim) {return this->_dim2class[dim];};
1135: void setDofClass(const int dim, const int dofClass) {this->_dim2class[dim] = dofClass;};
1136: public:
1137: const int *getIndices() {return this->_indices;};
1138: const int *getIndices(const int marker) {
1139: if (!marker) return this->getIndices();
1140: return this->_exclusionIndices[marker];
1141: };
1142: void setIndices(const int *indices) {this->_indices = indices;};
1143: void setIndices(const int *indices, const int marker) {
1144: if (!marker) this->_indices = indices;
1145: this->_exclusionIndices[marker] = indices;
1146: };
1147: template<typename Bundle>
1148: int sizeV(Bundle& mesh) {
1149: typedef typename ISieveVisitor::PointRetriever<typename Bundle::sieve_type> Visitor;
1150: Visitor pV((int) pow((double) mesh.getSieve()->getMaxConeSize(), mesh.depth())+1, true);
1151: ISieveTraversal<typename Bundle::sieve_type>::orientedClosure(*mesh.getSieve(), *mesh.heightStratum(0)->begin(), pV);
1152: const typename Visitor::point_type *oPoints = pV.getPoints();
1153: const int oSize = pV.getSize();
1154: int size = 0;
1156: for(int cl = 0; cl < oSize; ++cl) {
1157: size += this->_dim2dof[mesh.depth(oPoints[cl])];
1158: }
1159: return size;
1160: };
1161: template<typename Bundle>
1162: int size(const Obj<Bundle>& mesh) {
1163: const Obj<typename Bundle::label_sequence>& cells = mesh->heightStratum(0);
1164: const Obj<typename Bundle::coneArray> closure = ALE::SieveAlg<Bundle>::closure(mesh, *cells->begin());
1165: const typename Bundle::coneArray::iterator end = closure->end();
1166: int size = 0;
1168: for(typename Bundle::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1169: size += this->_dim2dof[mesh->depth(*cl_iter)];
1170: }
1171: return size;
1172: };
1173: };
1174: }
1176: namespace ALE {
1177: template<typename Sieve_,
1178: typename RealSection_ = Section<typename Sieve_::point_type, double>,
1179: typename IntSection_ = Section<typename Sieve_::point_type, int>,
1180: typename Label_ = LabelSifter<int, typename Sieve_::point_type>,
1181: typename ArrowSection_ = UniformSection<MinimalArrow<typename Sieve_::point_type, typename Sieve_::point_type>, int> >
1182: class IBundle : public ALE::ParallelObject {
1183: public:
1184: typedef Sieve_ sieve_type;
1185: typedef RealSection_ real_section_type;
1186: typedef IntSection_ int_section_type;
1187: typedef ArrowSection_ arrow_section_type;
1188: typedef IBundle<Sieve_,RealSection_,IntSection_,Label_,ArrowSection_> this_type;
1189: typedef typename sieve_type::point_type point_type;
1190: typedef malloc_allocator<point_type> alloc_type;
1191: typedef Label_ label_type;
1192: typedef typename std::map<const std::string, Obj<label_type> > labels_type;
1193: typedef typename label_type::supportSequence label_sequence;
1194: typedef std::map<std::string, Obj<arrow_section_type> > arrow_sections_type;
1195: typedef std::map<std::string, Obj<real_section_type> > real_sections_type;
1196: typedef std::map<std::string, Obj<int_section_type> > int_sections_type;
1197: typedef ALE::Point index_type;
1198: typedef std::pair<index_type, int> oIndex_type;
1199: typedef std::vector<oIndex_type> oIndexArray;
1200: typedef std::pair<int *, int> indices_type;
1201: typedef NumberingFactory<this_type> MeshNumberingFactory;
1202: typedef typename ALE::Partitioner<>::part_type rank_type;
1203: typedef typename ALE::Sifter<point_type,rank_type,point_type> send_overlap_type;
1204: typedef typename ALE::Sifter<rank_type,point_type,point_type> recv_overlap_type;
1205: typedef typename MeshNumberingFactory::numbering_type numbering_type;
1206: typedef typename MeshNumberingFactory::order_type order_type;
1207: typedef std::map<point_type, point_type> renumbering_type;
1208: // These should go away
1209: typedef typename ALE::SieveAlg<this_type> sieve_alg_type;
1210: typedef typename sieve_alg_type::coneArray coneArray;
1211: typedef typename sieve_alg_type::orientedConeArray oConeArray;
1212: typedef typename sieve_alg_type::supportArray supportArray;
1213: public:
1214: class LabelVisitor {
1215: protected:
1216: label_type& label;
1217: int defaultValue;
1218: int value;
1219: public:
1220: LabelVisitor(label_type& l, const int defValue) : label(l), defaultValue(defValue), value(defValue) {};
1221: int getLabelValue(const point_type& point) const {
1222: const Obj<typename label_type::coneSequence>& cone = this->label.cone(point);
1224: if (cone->size() == 0) return this->defaultValue;
1225: return *cone->begin();
1226: };
1227: void setLabelValue(const point_type& point, const int value) {
1228: this->label.setCone(value, point);
1229: };
1230: int getValue() const {return this->value;};
1231: };
1232: class MaxConeVisitor : public LabelVisitor {
1233: public:
1234: MaxConeVisitor(label_type& l, const int defValue) : LabelVisitor(l, defValue) {};
1235: void visitPoint(const typename sieve_type::point_type& point) {};
1236: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1237: this->value = std::max(this->value, this->getLabelValue(arrow.source));
1238: };
1239: };
1240: class MaxSupportVisitor : public LabelVisitor {
1241: public:
1242: MaxSupportVisitor(label_type& l, const int defValue) : LabelVisitor(l, defValue) {};
1243: void visitPoint(const typename sieve_type::point_type& point) {};
1244: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1245: this->value = std::max(this->value, this->getLabelValue(arrow.target));
1246: };
1247: };
1248: class HeightVisitor {
1249: protected:
1250: const sieve_type& sieve;
1251: label_type& height;
1252: int maxHeight;
1253: std::set<typename sieve_type::point_type> modifiedPoints;
1254: public:
1255: HeightVisitor(const sieve_type& s, label_type& h) : sieve(s), height(h), maxHeight(-1) {};
1256: void visitPoint(const typename sieve_type::point_type& point) {
1257: MaxSupportVisitor v(height, -1);
1259: // Compute the max height of the points in the support of p, and add 1
1260: this->sieve.support(point, v);
1261: const int h0 = v.getLabelValue(point);
1262: const int h1 = v.getValue() + 1;
1264: if(h1 != h0) {
1265: v.setLabelValue(point, h1);
1266: if (h1 > this->maxHeight) this->maxHeight = h1;
1267: this->modifiedPoints.insert(point);
1268: }
1269: };
1270: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1271: this->visitPoint(arrow.source);
1272: };
1273: int getMaxHeight() const {return this->maxHeight;};
1274: bool isModified() const {return this->modifiedPoints.size() > 0;};
1275: const std::set<typename sieve_type::point_type>& getModifiedPoints() const {return this->modifiedPoints;};
1276: void clear() {this->modifiedPoints.clear();};
1277: };
1278: class DepthVisitor {
1279: public:
1280: typedef typename sieve_type::point_type point_type;
1281: protected:
1282: const sieve_type& sieve;
1283: label_type& depth;
1284: int maxDepth;
1285: const point_type limitPoint;
1286: std::set<point_type> modifiedPoints;
1287: public:
1288: DepthVisitor(const sieve_type& s, label_type& d) : sieve(s), depth(d), maxDepth(-1), limitPoint(sieve.getChart().max()+1) {};
1289: DepthVisitor(const sieve_type& s, const point_type& limit, label_type& d) : sieve(s), depth(d), maxDepth(-1), limitPoint(limit) {};
1290: void visitPoint(const point_type& point) {
1291: if (point >= this->limitPoint) return;
1292: MaxConeVisitor v(depth, -1);
1294: // Compute the max height of the points in the support of p, and add 1
1295: this->sieve.cone(point, v);
1296: const int d0 = v.getLabelValue(point);
1297: const int d1 = v.getValue() + 1;
1299: if(d1 != d0) {
1300: v.setLabelValue(point, d1);
1301: if (d1 > this->maxDepth) this->maxDepth = d1;
1302: this->modifiedPoints.insert(point);
1303: }
1304: };
1305: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1306: this->visitPoint(arrow.target);
1307: };
1308: int getMaxDepth() const {return this->maxDepth;};
1309: bool isModified() const {return this->modifiedPoints.size() > 0;};
1310: const std::set<typename sieve_type::point_type>& getModifiedPoints() const {return this->modifiedPoints;};
1311: void clear() {this->modifiedPoints.clear();};
1312: };
1313: protected:
1314: Obj<sieve_type> _sieve;
1315: labels_type _labels;
1316: int _maxHeight;
1317: int _maxDepth;
1318: arrow_sections_type _arrowSections;
1319: real_sections_type _realSections;
1320: int_sections_type _intSections;
1321: Obj<oIndexArray> _indexArray;
1322: Obj<MeshNumberingFactory> _factory;
1323: bool _calculatedOverlap;
1324: Obj<send_overlap_type> _sendOverlap;
1325: Obj<recv_overlap_type> _recvOverlap;
1326: Obj<send_overlap_type> _distSendOverlap;
1327: Obj<recv_overlap_type> _distRecvOverlap;
1328: renumbering_type _renumbering;
1329: // Work space
1330: Obj<std::set<point_type> > _modifiedPoints;
1331: public:
1332: IBundle(MPI_Comm comm, int debug = 0) : ALE::ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1) {
1333: this->_indexArray = new oIndexArray();
1334: this->_modifiedPoints = new std::set<point_type>();
1335: this->_factory = MeshNumberingFactory::singleton(this->comm(), this->debug());
1336: this->_calculatedOverlap = false;
1337: this->_sendOverlap = new send_overlap_type(comm, debug);
1338: this->_recvOverlap = new recv_overlap_type(comm, debug);
1339: };
1340: IBundle(const Obj<sieve_type>& sieve) : ALE::ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _maxHeight(-1), _maxDepth(-1) {
1341: this->_indexArray = new oIndexArray();
1342: this->_modifiedPoints = new std::set<point_type>();
1343: this->_factory = MeshNumberingFactory::singleton(this->comm(), this->debug());
1344: this->_calculatedOverlap = false;
1345: this->_sendOverlap = new send_overlap_type(comm, debug);
1346: this->_recvOverlap = new recv_overlap_type(comm, debug);
1347: };
1348: virtual ~IBundle() {};
1349: public: // Verifiers
1350: bool hasLabel(const std::string& name) {
1351: if (this->_labels.find(name) != this->_labels.end()) {
1352: return true;
1353: }
1354: return false;
1355: };
1356: void checkLabel(const std::string& name) {
1357: if (!this->hasLabel(name)) {
1358: ostringstream msg;
1359: msg << "Invalid label name: " << name << std::endl;
1360: throw ALE::Exception(msg.str().c_str());
1361: }
1362: };
1363: public: // Accessors
1364: const Obj<sieve_type>& getSieve() const {return this->_sieve;};
1365: void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
1366: bool hasArrowSection(const std::string& name) const {
1367: return this->_arrowSections.find(name) != this->_arrowSections.end();
1368: };
1369: const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
1370: if (!this->hasArrowSection(name)) {
1371: Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());
1373: section->setName(name);
1374: if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
1375: this->_arrowSections[name] = section;
1376: }
1377: return this->_arrowSections[name];
1378: };
1379: void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
1380: this->_arrowSections[name] = section;
1381: };
1382: Obj<std::set<std::string> > getArrowSections() const {
1383: Obj<std::set<std::string> > names = std::set<std::string>();
1385: for(typename arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
1386: names->insert(s_iter->first);
1387: }
1388: return names;
1389: };
1390: bool hasRealSection(const std::string& name) const {
1391: return this->_realSections.find(name) != this->_realSections.end();
1392: };
1393: const Obj<real_section_type>& getRealSection(const std::string& name) {
1394: if (!this->hasRealSection(name)) {
1395: Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());
1397: section->setName(name);
1398: if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
1399: this->_realSections[name] = section;
1400: }
1401: return this->_realSections[name];
1402: };
1403: void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
1404: this->_realSections[name] = section;
1405: };
1406: Obj<std::set<std::string> > getRealSections() const {
1407: Obj<std::set<std::string> > names = std::set<std::string>();
1409: for(typename real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
1410: names->insert(s_iter->first);
1411: }
1412: return names;
1413: };
1414: bool hasIntSection(const std::string& name) const {
1415: return this->_intSections.find(name) != this->_intSections.end();
1416: };
1417: const Obj<int_section_type>& getIntSection(const std::string& name) {
1418: if (!this->hasIntSection(name)) {
1419: Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());
1421: section->setName(name);
1422: if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
1423: this->_intSections[name] = section;
1424: }
1425: return this->_intSections[name];
1426: };
1427: void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
1428: this->_intSections[name] = section;
1429: };
1430: Obj<std::set<std::string> > getIntSections() const {
1431: Obj<std::set<std::string> > names = std::set<std::string>();
1433: for(typename int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
1434: names->insert(s_iter->first);
1435: }
1436: return names;
1437: };
1438: const Obj<MeshNumberingFactory>& getFactory() const {return this->_factory;};
1439: renumbering_type& getRenumbering() {return this->_renumbering;};
1440: public: // Labels
1441: int getValue (const Obj<label_type>& label, const point_type& point, const int defValue = 0) {
1442: const Obj<typename label_type::coneSequence>& cone = label->cone(point);
1444: if (cone->size() == 0) return defValue;
1445: return *cone->begin();
1446: };
1447: void setValue(const Obj<label_type>& label, const point_type& point, const int value) {
1448: label->setCone(value, point);
1449: };
1450: void addValue(const Obj<label_type>& label, const point_type& point, const int value) {
1451: label->addCone(value, point);
1452: };
1453: template<typename InputPoints>
1454: int getMaxValue (const Obj<label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
1455: int maxValue = defValue;
1457: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1458: maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
1459: }
1460: return maxValue;
1461: };
1462: const Obj<label_type>& createLabel(const std::string& name) {
1463: this->_labels[name] = new label_type(this->comm(), this->debug());
1464: return this->_labels[name];
1465: };
1466: const Obj<label_type>& getLabel(const std::string& name) {
1467: this->checkLabel(name);
1468: return this->_labels[name];
1469: };
1470: void setLabel(const std::string& name, const Obj<label_type>& label) {
1471: this->_labels[name] = label;
1472: };
1473: const labels_type& getLabels() {
1474: return this->_labels;
1475: };
1476: virtual const Obj<label_sequence>& getLabelStratum(const std::string& name, int value) {
1477: this->checkLabel(name);
1478: return this->_labels[name]->support(value);
1479: };
1480: public: // Stratification
1481: void computeHeights() {
1482: const Obj<label_type>& label = this->createLabel(std::string("height"));
1483: HeightVisitor h(*this->_sieve, *label);
1485: #ifdef IMESH_NEW_LABELS
1486: label->setChart(this->getSieve()->getChart());
1487: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setConeSize(p, 1);}
1488: if (label->getChart().size()) {label->setSupportSize(0, label->getChart().size());}
1489: label->allocate();
1490: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setCone(-1, p);}
1491: #endif
1492: this->_sieve->leaves(h);
1493: while(h.isModified()) {
1494: // FIX: Avoid the copy here somehow by fixing the traversal
1495: std::vector<point_type> modifiedPoints(h.getModifiedPoints().begin(), h.getModifiedPoints().end());
1497: h.clear();
1498: this->_sieve->cone(modifiedPoints, h);
1499: }
1500: #ifdef IMESH_NEW_LABELS
1501: // Recalculate supportOffsets and populate support
1502: label->recalculateLabel();
1503: #endif
1504: this->_maxHeight = h.getMaxHeight();
1505: };
1506: virtual int height() const {return this->_maxHeight;};
1507: virtual void setHeight(const int height) {this->_maxHeight = height;};
1508: virtual int height(const point_type& point) {
1509: return this->getValue(this->_labels["height"], point, -1);
1510: };
1511: virtual void setHeight(const point_type& point, const int height) {
1512: return this->setValue(this->_labels["height"], point, height);
1513: };
1514: virtual const Obj<label_sequence>& heightStratum(int height) {
1515: return this->getLabelStratum("height", height);
1516: };
1517: void computeDepths() {
1518: const Obj<label_type>& label = this->createLabel(std::string("depth"));
1519: DepthVisitor d(*this->_sieve, *label);
1521: #ifdef IMESH_NEW_LABELS
1522: label->setChart(this->getSieve()->getChart());
1523: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setConeSize(p, 1);}
1524: if (label->getChart().size()) {label->setSupportSize(0, label->getChart().size());}
1525: label->allocate();
1526: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setCone(-1, p);}
1527: #endif
1528: this->_sieve->roots(d);
1529: while(d.isModified()) {
1530: // FIX: Avoid the copy here somehow by fixing the traversal
1531: std::vector<point_type> modifiedPoints(d.getModifiedPoints().begin(), d.getModifiedPoints().end());
1533: d.clear();
1534: this->_sieve->support(modifiedPoints, d);
1535: }
1536: #ifdef IMESH_NEW_LABELS
1537: // Recalculate supportOffsets and populate support
1538: label->recalculateLabel();
1539: #endif
1540: this->_maxDepth = d.getMaxDepth();
1541: };
1542: virtual int depth() const {return this->_maxDepth;};
1543: virtual void setDepth(const int depth) {this->_maxDepth = depth;};
1544: virtual int depth(const point_type& point) {
1545: return this->getValue(this->_labels["depth"], point, -1);
1546: };
1547: virtual void setDepth(const point_type& point, const int depth) {
1548: return this->setValue(this->_labels["depth"], point, depth);
1549: };
1550: virtual const Obj<label_sequence>& depthStratum(int depth) {
1551: return this->getLabelStratum("depth", depth);
1552: };
1553: void stratify() {
1554: this->computeHeights();
1555: this->computeDepths();
1556: };
1557: protected:
1558: template<typename Value>
1559: static bool lt1(const Value& a, const Value& b) {
1560: return a.first < b.first;
1561: };
1562: public: // Allocation
1563: template<typename Section_>
1564: void reallocate(const Obj<Section_>& section) {
1565: if (!section->hasNewPoints()) return;
1566: typename Section_::chart_type newChart(std::min(std::min_element(section->getNewPoints().begin(), section->getNewPoints().end(), lt1<typename Section_::newpoint_type>)->first, section->getChart().min()),
1567: std::max(std::max_element(section->getNewPoints().begin(), section->getNewPoints().end(), lt1<typename Section_::newpoint_type>)->first, section->getChart().max()-1)+1);
1568: section->reallocatePoint(newChart);
1569: };
1570: };
1571: #ifdef IMESH_NEW_LABELS
1572: template<typename Label_ = IFSieve<int> >
1573: #else
1574: template<typename Label_ = LabelSifter<int, int> >
1575: #endif
1576: class IMesh : public IBundle<IFSieve<int>, IGeneralSection<int, double>, IGeneralSection<int, int>, Label_> {
1577: public:
1578: typedef IBundle<IFSieve<int>, IGeneralSection<int, double>, IGeneralSection<int, int>, Label_> base_type;
1579: typedef typename base_type::sieve_type sieve_type;
1580: typedef typename sieve_type::point_type point_type;
1581: typedef typename base_type::alloc_type alloc_type;
1582: typedef typename base_type::label_type label_type;
1583: typedef typename base_type::label_sequence label_sequence;
1584: typedef typename base_type::real_section_type real_section_type;
1585: typedef typename base_type::numbering_type numbering_type;
1586: typedef typename base_type::order_type order_type;
1587: typedef typename base_type::send_overlap_type send_overlap_type;
1588: typedef typename base_type::recv_overlap_type recv_overlap_type;
1589: typedef std::set<std::string> names_type;
1590: typedef std::vector<typename PETSc::Point<3> > holes_type;
1591: typedef std::map<std::string, Obj<Discretization> > discretizations_type;
1592: protected:
1593: int _dim;
1594: bool _calculatedOverlap;
1595: Obj<send_overlap_type> _sendOverlap;
1596: Obj<recv_overlap_type> _recvOverlap;
1597: std::map<int,double> _periodicity;
1598: holes_type _holes;
1599: discretizations_type _discretizations;
1600: int _maxDof;
1601: public:
1602: IMesh(MPI_Comm comm, int dim, int debug = 0) : base_type(comm, debug), _dim(dim) {
1603: this->_calculatedOverlap = false;
1604: this->_sendOverlap = new send_overlap_type(comm, debug);
1605: this->_recvOverlap = new recv_overlap_type(comm, debug);
1606: this->_maxDof = -1;
1607: };
1608: public: // Accessors
1609: int getDimension() const {return this->_dim;};
1610: void setDimension(const int dim) {this->_dim = dim;};
1611: bool getCalculatedOverlap() const {return this->_calculatedOverlap;};
1612: void setCalculatedOverlap(const bool calc) {this->_calculatedOverlap = calc;};
1613: const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
1614: void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
1615: const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
1616: void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
1617: bool getPeriodicity(const int d) {return this->_periodicity[d];};
1618: void setPeriodicity(const int d, const double length) {this->_periodicity[d] = length;};
1619: const holes_type& getHoles() const {return this->_holes;};
1620: void addHole(const double hole[]) {
1621: this->_holes.push_back(hole);
1622: };
1623: void copyHoles(const Obj<IMesh>& m) {
1624: const holes_type& holes = m->getHoles();
1626: for(holes_type::const_iterator h_iter = holes.begin(); h_iter != holes.end(); ++h_iter) {
1627: this->_holes.push_back(*h_iter);
1628: }
1629: };
1630: const Obj<Discretization>& getDiscretization() {return this->getDiscretization("default");};
1631: const Obj<Discretization>& getDiscretization(const std::string& name) {return this->_discretizations[name];};
1632: void setDiscretization(const Obj<Discretization>& disc) {this->setDiscretization("default", disc);};
1633: void setDiscretization(const std::string& name, const Obj<Discretization>& disc) {this->_discretizations[name] = disc;};
1634: Obj<names_type> getDiscretizations() const {
1635: Obj<names_type> names = names_type();
1637: for(discretizations_type::const_iterator d_iter = this->_discretizations.begin(); d_iter != this->_discretizations.end(); ++d_iter) {
1638: names->insert(d_iter->first);
1639: }
1640: return names;
1641: };
1642: int getMaxDof() const {return this->_maxDof;};
1643: void setMaxDof(const int maxDof) {this->_maxDof = maxDof;};
1644: public: // Sizes
1645: template<typename Section>
1646: int size(const Obj<Section>& section, const point_type& p) {
1647: typedef ISieveVisitor::SizeVisitor<sieve_type,Section> size_visitor_type;
1648: typedef ISieveVisitor::TransitiveClosureVisitor<sieve_type,size_visitor_type> closure_visitor_type;
1649: size_visitor_type sV(*section);
1650: closure_visitor_type cV(*this->getSieve(), sV);
1652: this->getSieve()->cone(p, cV);
1653: if (!sV.getSize()) sV.visitPoint(p);
1654: return sV.getSize();
1655: };
1656: template<typename Section>
1657: int sizeWithBC(const Obj<Section>& section, const point_type& p) {
1658: typedef ISieveVisitor::SizeWithBCVisitor<sieve_type,Section> size_visitor_type;
1659: typedef ISieveVisitor::TransitiveClosureVisitor<sieve_type,size_visitor_type> closure_visitor_type;
1660: size_visitor_type sV(*section);
1661: closure_visitor_type cV(*this->getSieve(), sV);
1663: this->getSieve()->cone(p, cV);
1664: if (!sV.getSize()) sV.visitPoint(p);
1665: return sV.getSize();
1666: };
1667: template<typename Section>
1668: void allocate(const Obj<Section>& section) {
1669: section->allocatePoint();
1670: };
1671: public: // Restrict/Update closures
1672: template<typename Sieve, typename Visitor>
1673: void closure1(const Sieve& sieve, const point_type& p, Visitor& v)
1674: {
1675: v.visitPoint(p, 0);
1676: // Cone is guarateed to be ordered correctly
1677: sieve.orientedCone(p, v);
1678: };
1679: // Return the values for the closure of this point
1680: template<typename Section>
1681: const typename Section::value_type *restrictClosure(const Obj<Section>& section, const point_type& p) {
1682: const int size = this->sizeWithBC(section, p);
1683: ISieveVisitor::RestrictVisitor<Section> rV(*section, size, section->getRawArray(size));
1685: if (this->depth() == 1) {
1686: closure1(*this->getSieve(), p, rV);
1687: } else {
1688: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::RestrictVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, rV, true);
1690: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1691: }
1692: return rV.getValues();
1693: };
1694: template<typename Section>
1695: const typename Section::value_type *restrictClosure(const Obj<Section>& section, const point_type& p, typename Section::value_type *values, const int valuesSize) {
1696: const int size = this->sizeWithBC(section, p);
1697: if (valuesSize < size) {throw ALE::Exception("Input array to small for restrictClosure()");}
1698: ISieveVisitor::RestrictVisitor<Section> rV(*section, size, values);
1700: if (this->depth() == 1) {
1701: closure1(*this->getSieve(), p, rV);
1702: } else {
1703: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::RestrictVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, rV, true);
1705: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1706: }
1707: return rV.getValues();
1708: };
1709: template<typename Visitor>
1710: void restrictClosure(const point_type& p, Visitor& v) {
1711: if (this->depth() == 1) {
1712: closure1(*this->getSieve(), p, v);
1713: } else {
1714: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, v);
1715: }
1716: };
1717: // Replace the values for the closure of this point
1718: template<typename Section>
1719: void update(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1720: ISieveVisitor::UpdateVisitor<Section> uV(*section, v);
1722: if (this->depth() == 1) {
1723: closure1(*this->getSieve(), p, uV);
1724: } else {
1725: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);
1727: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1728: }
1729: };
1730: // Replace the values for the closure of this point, including points constrained by BC
1731: template<typename Section>
1732: void updateAll(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1733: ISieveVisitor::UpdateAllVisitor<Section> uV(*section, v);
1735: if (this->depth() == 1) {
1736: closure1(*this->getSieve(), p, uV);
1737: } else {
1738: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateAllVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);
1740: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1741: }
1742: };
1743: // Augment the values for the closure of this point
1744: template<typename Section>
1745: void updateAdd(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1746: ISieveVisitor::UpdateAddVisitor<Section> uV(*section, v);
1748: if (this->depth() == 1) {
1749: closure1(*this->getSieve(), p, uV);
1750: } else {
1751: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateAddVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);
1753: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1754: }
1755: };
1756: // Augment the values for the closure of this point
1757: template<typename Visitor>
1758: void updateClosure(const point_type& p, Visitor& v) {
1759: if (this->depth() == 1) {
1760: closure1(*this->getSieve(), p, v);
1761: } else {
1762: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, v);
1763: }
1764: };
1765: public: // Overlap
1766: void constructOverlap() {
1767: if (!this->_calculatedOverlap && (this->commSize() > 1)) {throw ALE::Exception("Must calculate overlap during distribution");}
1768: };
1769: public: // Cell topology and geometry
1770: int getNumCellCorners(const point_type& p, const int depth = -1) const {
1771: const int d = depth < 0 ? this->depth() : depth;
1773: if (d == 1) {
1774: return this->_sieve->getConeSize(p);
1775: } else if (d <= 0) {
1776: return 0;
1777: }
1778: // Warning: this is slow
1779: ISieveVisitor::NConeRetriever<sieve_type> ncV(*this->_sieve, (int) pow((double) this->_sieve->getMaxConeSize(), this->depth()));
1780: ALE::ISieveTraversal<sieve_type>::orientedClosure(*this->_sieve, p, ncV);
1781: return ncV.getOrientedSize();
1782: };
1783: int getNumCellCorners() {
1784: return getNumCellCorners(*this->heightStratum(0)->begin());
1785: };
1786: void setupCoordinates(const Obj<real_section_type>& coordinates) {
1787: const Obj<label_sequence>& vertices = this->depthStratum(0);
1789: coordinates->setChart(typename real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()),
1790: *std::max_element(vertices->begin(), vertices->end())+1));
1791: };
1792: point_type locatePoint(const typename real_section_type::value_type point[], point_type guess = -1) {
1793: //guess overrides this by saying that we already know the relation of this point to this mesh. We will need to make it a more robust "guess" later for more than P1
1794: if (guess != -1) {
1795: return guess;
1796: } else {
1797: throw ALE::Exception("No point location for mesh dimension");
1798: }
1799: };
1800: void computeTriangleGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1801: const double *coords = this->restrictClosure(coordinates, e);
1802: const int dim = 2;
1803: double invDet;
1805: if (v0) {
1806: for(int d = 0; d < dim; d++) {
1807: v0[d] = coords[d];
1808: }
1809: }
1810: if (J) {
1811: for(int d = 0; d < dim; d++) {
1812: for(int f = 0; f < dim; f++) {
1813: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1814: }
1815: }
1816: detJ = J[0]*J[3] - J[1]*J[2];
1817: if (detJ < 0.0) {
1818: const double xLength = this->_periodicity[0];
1820: if (xLength != 0.0) {
1821: double v0x = coords[0*dim+0];
1823: if (v0x == 0.0) {
1824: v0x = v0[0] = xLength;
1825: }
1826: for(int f = 0; f < dim; f++) {
1827: const double px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
1829: J[0*dim+f] = 0.5*(px - v0x);
1830: }
1831: }
1832: detJ = J[0]*J[3] - J[1]*J[2];
1833: }
1834: PetscLogFlopsNoError(8.0 + 3.0);
1835: }
1836: if (invJ) {
1837: invDet = 1.0/detJ;
1838: invJ[0] = invDet*J[3];
1839: invJ[1] = -invDet*J[1];
1840: invJ[2] = -invDet*J[2];
1841: invJ[3] = invDet*J[0];
1842: PetscLogFlopsNoError(5.0);
1843: }
1844: };
1845: void computeTetrahedronGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1846: const double *coords = this->restrictClosure(coordinates, e);
1847: const int dim = 3;
1848: double invDet;
1850: if (v0) {
1851: for(int d = 0; d < dim; d++) {
1852: v0[d] = coords[d];
1853: }
1854: }
1855: if (J) {
1856: for(int d = 0; d < dim; d++) {
1857: for(int f = 0; f < dim; f++) {
1858: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1859: }
1860: }
1861: // The minus sign is here since I orient the first face to get the outward normal
1862: detJ = -(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
1863: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
1864: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
1865: PetscLogFlopsNoError(18.0 + 12.0);
1866: }
1867: if (invJ) {
1868: invDet = -1.0/detJ;
1869: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
1870: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
1871: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
1872: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
1873: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
1874: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
1875: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1876: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
1877: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
1878: PetscLogFlopsNoError(37.0);
1879: }
1880: };
1881: void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1882: if (this->_dim == 2) {
1883: computeTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
1884: } else if (this->_dim == 3) {
1885: computeTetrahedronGeometry(coordinates, e, v0, J, invJ, detJ);
1886: } else {
1887: throw ALE::Exception("Unsupported dimension for element geometry computation");
1888: }
1889: };
1890: void computeBdSegmentGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1891: const double *coords = this->restrictClosure(coordinates, e);
1892: const int dim = 2;
1893: double invDet;
1895: if (v0) {
1896: for(int d = 0; d < dim; d++) {
1897: v0[d] = coords[d];
1898: }
1899: }
1900: if (J) {
1901: //r2 = coords[1*dim+0]*coords[1*dim+0] + coords[1*dim+1]*coords[1*dim+1];
1902: J[0] = (coords[1*dim+0] - coords[0*dim+0])*0.5; J[1] = (-coords[1*dim+1] + coords[0*dim+1])*0.5;
1903: J[2] = (coords[1*dim+1] - coords[0*dim+1])*0.5; J[3] = ( coords[1*dim+0] - coords[0*dim+0])*0.5;
1904: detJ = J[0]*J[3] - J[1]*J[2];
1905: if (detJ < 0.0) {
1906: const double xLength = this->_periodicity[0];
1908: if (xLength != 0.0) {
1909: double v0x = coords[0*dim+0];
1911: if (v0x == 0.0) {
1912: v0x = v0[0] = xLength;
1913: }
1914: for(int f = 0; f < dim; f++) {
1915: const double px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
1917: J[0*dim+f] = 0.5*(px - v0x);
1918: }
1919: }
1920: detJ = J[0]*J[3] - J[1]*J[2];
1921: }
1922: PetscLogFlopsNoError(8.0 + 3.0);
1923: }
1924: if (invJ) {
1925: invDet = 1.0/detJ;
1926: invJ[0] = invDet*J[3];
1927: invJ[1] = -invDet*J[1];
1928: invJ[2] = -invDet*J[2];
1929: invJ[3] = invDet*J[0];
1930: PetscLogFlopsNoError(5.0);
1931: }
1932: };
1933: void computeBdElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1934: if (this->_dim == 1) {
1935: computeBdSegmentGeometry(coordinates, e, v0, J, invJ, detJ);
1936: // } else if (this->_dim == 2) {
1937: // computeBdTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
1938: } else {
1939: throw ALE::Exception("Unsupported dimension for element geometry computation");
1940: }
1941: };
1942: double getMaxVolume() {
1943: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
1944: const Obj<label_sequence>& cells = this->heightStratum(0);
1945: const int dim = this->getDimension();
1946: double v0[3], J[9], invJ[9], detJ, refVolume = 0.0, maxVolume = 0.0;
1948: if (dim == 1) refVolume = 2.0;
1949: if (dim == 2) refVolume = 2.0;
1950: if (dim == 3) refVolume = 4.0/3.0;
1951: for(typename label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1952: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
1953: maxVolume = std::max(maxVolume, detJ*refVolume);
1954: }
1955: return maxVolume;
1956: };
1957: public:
1958: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1959: if (comm == MPI_COMM_NULL) {
1960: comm = this->comm();
1961: }
1962: if (name == "") {
1963: PetscPrintf(comm, "viewing a Mesh\n");
1964: } else {
1965: PetscPrintf(comm, "viewing Mesh '%s'\n", name.c_str());
1966: }
1967: this->getSieve()->view("mesh sieve", comm);
1968: Obj<names_type> sections = this->getRealSections();
1970: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1971: this->getRealSection(*name)->view(*name);
1972: }
1973: sections = this->getIntSections();
1974: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1975: this->getIntSection(*name)->view(*name);
1976: }
1977: sections = this->getArrowSections();
1978: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1979: this->getArrowSection(*name)->view(*name);
1980: }
1981: };
1982: public: // Discretization
1983: void markBoundaryCells(const std::string& name, const int marker = 1, const int newMarker = 2, const bool onlyVertices = false) {
1984: const Obj<label_type>& label = this->getLabel(name);
1985: const Obj<label_sequence>& boundary = this->getLabelStratum(name, marker);
1986: const typename label_sequence::iterator end = boundary->end();
1987: const Obj<sieve_type>& sieve = this->getSieve();
1989: if (!onlyVertices) {
1990: typename ISieveVisitor::MarkVisitor<sieve_type,label_type> mV(*label, newMarker);
1992: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != end; ++e_iter) {
1993: if (this->height(*e_iter) == 1) {
1994: sieve->support(*e_iter, mV);
1995: }
1996: }
1997: } else {
1998: #if 1
1999: throw ALE::Exception("Rewrite this to first mark bounadry edges/faces.");
2000: #else
2001: const int depth = this->depth();
2003: for(typename label_sequence::iterator v_iter = boundary->begin(); v_iter != end; ++v_iter) {
2004: const Obj<supportArray> support = sieve->nSupport(*v_iter, depth);
2005: const typename supportArray::iterator sEnd = support->end();
2007: for(typename supportArray::iterator c_iter = support->begin(); c_iter != sEnd; ++c_iter) {
2008: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*c_iter);
2009: const typename sieve_type::traits::coneSequence::iterator cEnd = cone->end();
2011: for(typename sieve_type::traits::coneSequence::iterator e_iter = cone->begin(); e_iter != cEnd; ++e_iter) {
2012: if (sieve->support(*e_iter)->size() == 1) {
2013: this->setValue(label, *c_iter, newMarker);
2014: break;
2015: }
2016: }
2017: }
2018: }
2019: #endif
2020: }
2021: };
2022: int setFiberDimensions(const Obj<real_section_type>& s, const Obj<names_type>& discs, names_type& bcLabels) {
2023: const int debug = this->debug();
2024: int maxDof = 0;
2026: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
2027: s->addSpace();
2028: }
2029: for(int d = 0; d <= this->_dim; ++d) {
2030: int numDof = 0;
2031: int f = 0;
2033: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2034: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2035: const int sDof = disc->getNumDof(d);
2037: numDof += sDof;
2038: if (sDof) s->setFiberDimension(this->depthStratum(d), sDof, f);
2039: }
2040: if (numDof) s->setFiberDimension(this->depthStratum(d), numDof);
2041: maxDof = std::max(maxDof, numDof);
2042: }
2043: // Process exclusions
2044: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2045: int f = 0;
2047: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2048: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2049: std::string labelName = "exclude-"+*f_iter;
2050: std::set<point_type> seen;
2051: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth()), true);
2053: if (this->hasLabel(labelName)) {
2054: const Obj<label_type>& label = this->getLabel(labelName);
2055: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
2056: const typename label_sequence::iterator end = exclusion->end();
2057: if (debug > 1) {label->view(labelName.c_str());}
2059: for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2060: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *e_iter, pV);
2061: const typename Visitor::point_type *oPoints = pV.getPoints();
2062: const int oSize = pV.getSize();
2064: for(int cl = 0; cl < oSize; ++cl) {
2065: if (seen.find(oPoints[cl]) != seen.end()) continue;
2066: if (this->getValue(label, oPoints[cl]) == 1) {
2067: seen.insert(oPoints[cl]);
2068: s->setFiberDimension(oPoints[cl], 0, f);
2069: s->addFiberDimension(oPoints[cl], -disc->getNumDof(this->depth(oPoints[cl])));
2070: if (debug > 1) {std::cout << " point: " << oPoints[cl] << " dim: " << disc->getNumDof(this->depth(oPoints[cl])) << std::endl;}
2071: }
2072: }
2073: pV.clear();
2074: }
2075: }
2076: }
2077: // Process constraints
2078: f = 0;
2079: for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2080: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2081: const Obj<std::set<std::string> > bcs = disc->getBoundaryConditions();
2082: std::string excludeName = "exclude-"+*f_iter;
2084: for(std::set<std::string>::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter) {
2085: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
2086: const Obj<label_sequence>& boundary = this->getLabelStratum(bc->getLabelName(), bc->getMarker());
2088: bcLabels.insert(bc->getLabelName());
2089: if (this->hasLabel(excludeName)) {
2090: const Obj<label_type>& label = this->getLabel(excludeName);
2092: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2093: if (!this->getValue(label, *e_iter)) {
2094: const int numDof = disc->getNumDof(this->depth(*e_iter));
2096: if (numDof) s->addConstraintDimension(*e_iter, numDof);
2097: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2098: }
2099: }
2100: } else {
2101: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2102: const int numDof = disc->getNumDof(this->depth(*e_iter));
2104: if (numDof) s->addConstraintDimension(*e_iter, numDof);
2105: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2106: }
2107: }
2108: }
2109: }
2110: return maxDof;
2111: };
2112: void calculateIndices() {
2113: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2114: // Should have an iterator over the whole tree
2115: Obj<names_type> discs = this->getDiscretizations();
2116: const int debug = this->debug();
2117: std::map<std::string, std::pair<int, int*> > indices;
2119: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2120: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2122: indices[*d_iter] = std::pair<int, int*>(0, new int[disc->sizeV(*this)]);
2123: disc->setIndices(indices[*d_iter].second);
2124: }
2125: const Obj<label_sequence>& cells = this->heightStratum(0);
2126: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);
2127: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *cells->begin(), pV);
2128: const typename Visitor::point_type *oPoints = pV.getPoints();
2129: const int oSize = pV.getSize();
2130: int offset = 0;
2132: if (debug > 1) {std::cout << "Closure for first element" << std::endl;}
2133: for(int cl = 0; cl < oSize; ++cl) {
2134: const int dim = this->depth(oPoints[cl]);
2136: if (debug > 1) {std::cout << " point " << oPoints[cl] << " depth " << dim << std::endl;}
2137: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2138: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2139: const int num = disc->getNumDof(dim);
2141: if (debug > 1) {std::cout << " disc " << disc->getName() << " numDof " << num << std::endl;}
2142: for(int o = 0; o < num; ++o) {
2143: indices[*d_iter].second[indices[*d_iter].first++] = offset++;
2144: }
2145: }
2146: }
2147: pV.clear();
2148: if (debug > 1) {
2149: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2150: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2152: std::cout << "Discretization " << disc->getName() << " indices:";
2153: for(int i = 0; i < indices[*d_iter].first; ++i) {
2154: std::cout << " " << indices[*d_iter].second[i];
2155: }
2156: std::cout << std::endl;
2157: }
2158: }
2159: };
2160: void calculateIndicesExcluded(const Obj<real_section_type>& s, const Obj<names_type>& discs) {
2161: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2162: typedef std::map<std::string, std::pair<int, indexSet> > indices_type;
2163: const Obj<label_type>& indexLabel = this->createLabel("cellExclusion");
2164: const int debug = this->debug();
2165: int marker = 0;
2166: std::map<indices_type, int> indexMap;
2167: indices_type indices;
2168: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);
2170: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2171: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2172: const int size = disc->sizeV(*this);
2174: indices[*d_iter].second.resize(size);
2175: }
2176: const typename names_type::const_iterator dBegin = discs->begin();
2177: const typename names_type::const_iterator dEnd = discs->end();
2178: std::set<point_type> seen;
2179: int f = 0;
2181: for(typename names_type::const_iterator f_iter = dBegin; f_iter != dEnd; ++f_iter, ++f) {
2182: std::string labelName = "exclude-"+*f_iter;
2184: if (this->hasLabel(labelName)) {
2185: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
2186: const typename label_sequence::iterator end = exclusion->end();
2188: if (debug > 1) {std::cout << "Processing exclusion " << labelName << std::endl;}
2189: for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2190: if (this->height(*e_iter)) continue;
2191: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *e_iter, pV);
2192: const typename Visitor::point_type *oPoints = pV.getPoints();
2193: const int oSize = pV.getSize();
2194: int offset = 0;
2196: if (debug > 1) {std::cout << " Closure for cell " << *e_iter << std::endl;}
2197: for(int cl = 0; cl < oSize; ++cl) {
2198: int g = 0;
2200: if (debug > 1) {std::cout << " point " << oPoints[cl] << std::endl;}
2201: for(typename names_type::const_iterator g_iter = dBegin; g_iter != dEnd; ++g_iter, ++g) {
2202: const int fDim = s->getFiberDimension(oPoints[cl], g);
2204: if (debug > 1) {std::cout << " disc " << *g_iter << " numDof " << fDim << std::endl;}
2205: for(int d = 0; d < fDim; ++d) {
2206: indices[*g_iter].second[indices[*g_iter].first++] = offset++;
2207: }
2208: }
2209: }
2210: pV.clear();
2211: const std::map<indices_type, int>::iterator entry = indexMap.find(indices);
2213: if (debug > 1) {
2214: for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
2215: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2216: std::cout << "Discretization (" << i_iter->second << ") " << *g_iter << " indices:";
2217: for(int i = 0; i < ((indices_type) i_iter->first)[*g_iter].first; ++i) {
2218: std::cout << " " << ((indices_type) i_iter->first)[*g_iter].second[i];
2219: }
2220: std::cout << std::endl;
2221: }
2222: std::cout << "Comparison: " << (indices == i_iter->first) << std::endl;
2223: }
2224: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2225: std::cout << "Discretization " << *g_iter << " indices:";
2226: for(int i = 0; i < indices[*g_iter].first; ++i) {
2227: std::cout << " " << indices[*g_iter].second[i];
2228: }
2229: std::cout << std::endl;
2230: }
2231: }
2232: if (entry != indexMap.end()) {
2233: this->setValue(indexLabel, *e_iter, entry->second);
2234: if (debug > 1) {std::cout << " Found existing indices with marker " << entry->second << std::endl;}
2235: } else {
2236: indexMap[indices] = ++marker;
2237: this->setValue(indexLabel, *e_iter, marker);
2238: if (debug > 1) {std::cout << " Created new indices with marker " << marker << std::endl;}
2239: }
2240: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2241: indices[*g_iter].first = 0;
2242: for(unsigned int i = 0; i < indices[*g_iter].second.size(); ++i) indices[*g_iter].second[i] = 0;
2243: }
2244: }
2245: }
2246: }
2247: if (debug > 1) {indexLabel->view("cellExclusion");}
2248: for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
2249: if (debug > 1) {std::cout << "Setting indices for marker " << i_iter->second << std::endl;}
2250: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2251: const Obj<Discretization>& disc = this->getDiscretization(*g_iter);
2252: const indexSet indSet = ((indices_type) i_iter->first)[*g_iter].second;
2253: const int size = indSet.size();
2254: int *_indices = new int[size];
2256: if (debug > 1) {std::cout << " field " << *g_iter << std::endl;}
2257: for(int i = 0; i < size; ++i) {
2258: _indices[i] = indSet[i];
2259: if (debug > 1) {std::cout << " indices["<<i<<"] = " << _indices[i] << std::endl;}
2260: }
2261: disc->setIndices(_indices, i_iter->second);
2262: }
2263: }
2264: };
2265: void setupField(const Obj<real_section_type>& s, const int cellMarker = 2, const bool noUpdate = false) {
2266: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2267: const Obj<names_type>& discs = this->getDiscretizations();
2268: const int debug = s->debug();
2269: names_type bcLabels;
2271: s->setChart(this->getSieve()->getChart());
2272: this->_maxDof = this->setFiberDimensions(s, discs, bcLabels);
2273: this->calculateIndices();
2274: this->calculateIndicesExcluded(s, discs);
2275: this->allocate(s);
2276: s->defaultConstraintDof();
2277: const Obj<label_type>& cellExclusion = this->getLabel("cellExclusion");
2279: if (debug > 1) {std::cout << "Setting boundary values" << std::endl;}
2280: for(typename names_type::const_iterator n_iter = bcLabels.begin(); n_iter != bcLabels.end(); ++n_iter) {
2281: const Obj<label_sequence>& boundaryCells = this->getLabelStratum(*n_iter, cellMarker);
2282: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2283: const Obj<names_type>& discs = this->getDiscretizations();
2284: const point_type firstCell = *boundaryCells->begin();
2285: const int numFields = discs->size();
2286: typename real_section_type::value_type *values = new typename real_section_type::value_type[this->sizeWithBC(s, firstCell)];
2287: int *dofs = new int[this->_maxDof];
2288: int *v = new int[numFields];
2289: double *v0 = new double[this->getDimension()];
2290: double *J = new double[this->getDimension()*this->getDimension()];
2291: double detJ;
2292: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);
2294: for(typename label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
2295: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *c_iter, pV);
2296: const typename Visitor::point_type *oPoints = pV.getPoints();
2297: const int oSize = pV.getSize();
2299: if (debug > 1) {std::cout << " Boundary cell " << *c_iter << std::endl;}
2300: this->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
2301: for(int f = 0; f < numFields; ++f) v[f] = 0;
2302: for(int cl = 0; cl < oSize; ++cl) {
2303: const int cDim = s->getConstraintDimension(oPoints[cl]);
2304: int off = 0;
2305: int f = 0;
2306: int i = -1;
2308: if (debug > 1) {std::cout << " point " << oPoints[cl] << std::endl;}
2309: if (cDim) {
2310: if (debug > 1) {std::cout << " constrained excMarker: " << this->getValue(cellExclusion, *c_iter) << std::endl;}
2311: for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2312: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2313: const Obj<names_type> bcs = disc->getBoundaryConditions();
2314: const int fDim = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
2315: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
2316: int b = 0;
2318: for(typename names_type::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter, ++b) {
2319: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
2320: const int value = this->getValue(this->getLabel(bc->getLabelName()), oPoints[cl]);
2322: if (b > 0) v[f] -= fDim;
2323: if (value == bc->getMarker()) {
2324: if (debug > 1) {std::cout << " field " << *f_iter << " marker " << value << std::endl;}
2325: for(int d = 0; d < fDim; ++d, ++v[f]) {
2326: dofs[++i] = off+d;
2327: if (!noUpdate) values[indices[v[f]]] = (*bc->getDualIntegrator())(v0, J, v[f], bc->getFunction());
2328: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2329: }
2330: // Allow only one condition per point
2331: ++b;
2332: break;
2333: } else {
2334: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
2335: for(int d = 0; d < fDim; ++d, ++v[f]) {
2336: values[indices[v[f]]] = 0.0;
2337: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2338: }
2339: }
2340: }
2341: if (b == 0) {
2342: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
2343: for(int d = 0; d < fDim; ++d, ++v[f]) {
2344: values[indices[v[f]]] = 0.0;
2345: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2346: }
2347: }
2348: off += fDim;
2349: }
2350: if (i != cDim-1) {throw ALE::Exception("Invalid constraint initialization");}
2351: s->setConstraintDof(oPoints[cl], dofs);
2352: } else {
2353: if (debug > 1) {std::cout << " unconstrained" << std::endl;}
2354: for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2355: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2356: const int fDim = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
2357: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
2359: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
2360: for(int d = 0; d < fDim; ++d, ++v[f]) {
2361: values[indices[v[f]]] = 0.0;
2362: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2363: }
2364: }
2365: }
2366: }
2367: if (debug > 1) {
2368: for(int f = 0; f < numFields; ++f) v[f] = 0;
2369: for(int cl = 0; cl < oSize; ++cl) {
2370: int f = 0;
2371: for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2372: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2373: const int fDim = s->getFiberDimension(oPoints[cl], f);
2374: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
2376: for(int d = 0; d < fDim; ++d, ++v[f]) {
2377: std::cout << " "<<*f_iter<<"-value["<<indices[v[f]]<<"] " << values[indices[v[f]]] << std::endl;
2378: }
2379: }
2380: }
2381: }
2382: if (!noUpdate) {
2383: this->updateAll(s, *c_iter, values);
2384: }
2385: pV.clear();
2386: }
2387: delete [] dofs;
2388: delete [] values;
2389: delete [] v0;
2390: delete [] J;
2391: }
2392: if (debug > 1) {s->view("");}
2393: };
2394: };
2395: }
2397: namespace ALE {
2398: class Mesh : public Bundle<ALE::Sieve<int,int,int>, GeneralSection<int, double> > {
2399: public:
2400: typedef Bundle<ALE::Sieve<int,int,int>, GeneralSection<int, double> > base_type;
2401: typedef base_type::sieve_type sieve_type;
2402: typedef sieve_type::point_type point_type;
2403: typedef base_type::alloc_type alloc_type;
2404: typedef base_type::label_sequence label_sequence;
2405: typedef base_type::real_section_type real_section_type;
2406: typedef base_type::numbering_type numbering_type;
2407: typedef base_type::order_type order_type;
2408: typedef base_type::send_overlap_type send_overlap_type;
2409: typedef base_type::recv_overlap_type recv_overlap_type;
2410: typedef base_type::sieve_alg_type sieve_alg_type;
2411: typedef std::set<std::string> names_type;
2412: typedef std::map<std::string, Obj<Discretization> > discretizations_type;
2413: typedef std::vector<PETSc::Point<3> > holes_type;
2414: protected:
2415: int _dim;
2416: discretizations_type _discretizations;
2417: std::map<int,double> _periodicity;
2418: holes_type _holes;
2419: public:
2420: Mesh(MPI_Comm comm, int dim, int debug = 0) : base_type(comm, debug), _dim(dim) {
2421: ///this->_factory = MeshNumberingFactory::singleton(debug);
2422: //std::cout << "["<<this->commRank()<<"]: Creating an ALE::Mesh" << std::endl;
2423: };
2424: ~Mesh() {
2425: //std::cout << "["<<this->commRank()<<"]: Destroying an ALE::Mesh" << std::endl;
2426: };
2427: public: // Accessors
2428: int getDimension() const {return this->_dim;};
2429: void setDimension(const int dim) {this->_dim = dim;};
2430: const Obj<Discretization>& getDiscretization() {return this->getDiscretization("default");};
2431: const Obj<Discretization>& getDiscretization(const std::string& name) {return this->_discretizations[name];};
2432: void setDiscretization(const Obj<Discretization>& disc) {this->setDiscretization("default", disc);};
2433: void setDiscretization(const std::string& name, const Obj<Discretization>& disc) {this->_discretizations[name] = disc;};
2434: Obj<names_type> getDiscretizations() const {
2435: Obj<names_type> names = names_type();
2437: for(discretizations_type::const_iterator d_iter = this->_discretizations.begin(); d_iter != this->_discretizations.end(); ++d_iter) {
2438: names->insert(d_iter->first);
2439: }
2440: return names;
2441: };
2442: bool getPeriodicity(const int d) {return this->_periodicity[d];};
2443: void setPeriodicity(const int d, const double length) {this->_periodicity[d] = length;};
2444: const holes_type& getHoles() const {return this->_holes;};
2445: void addHole(const double hole[]) {
2446: this->_holes.push_back(hole);
2447: };
2448: void copyHoles(const Obj<Mesh>& m) {
2449: const holes_type& holes = m->getHoles();
2451: for(holes_type::const_iterator h_iter = holes.begin(); h_iter != holes.end(); ++h_iter) {
2452: this->_holes.push_back(*h_iter);
2453: }
2454: };
2455: void copy(const Obj<Mesh>& m) {
2456: this->setSieve(m->getSieve());
2457: this->setLabel("height", m->getLabel("height"));
2458: this->_maxHeight = m->height();
2459: this->setLabel("depth", m->getLabel("depth"));
2460: this->_maxDepth = m->depth();
2461: this->setLabel("marker", m->getLabel("marker"));
2462: this->setRealSection("coordinates", m->getRealSection("coordinates"));
2463: this->setArrowSection("orientation", m->getArrowSection("orientation"));
2464: };
2465: public: // Cell topology and geometry
2466: int getNumCellCorners(const point_type& p, const int depth = -1) const {
2467: return (this->getDimension() > 0) ? this->_sieve->nCone(p, depth < 0 ? this->depth() : depth)->size() : 1;
2468: };
2469: int getNumCellCorners() {
2470: return getNumCellCorners(*(this->heightStratum(0)->begin()));
2471: };
2472: void setupCoordinates(const Obj<real_section_type>& coordinates) {};
2473: void computeTriangleGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2474: const double *coords = this->restrictClosure(coordinates, e);
2475: const int dim = 2;
2476: double invDet;
2478: if (v0) {
2479: for(int d = 0; d < dim; d++) {
2480: v0[d] = coords[d];
2481: }
2482: }
2483: if (J) {
2484: for(int d = 0; d < dim; d++) {
2485: for(int f = 0; f < dim; f++) {
2486: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
2487: }
2488: }
2489: detJ = J[0]*J[3] - J[1]*J[2];
2490: if (detJ < 0.0) {
2491: const double xLength = this->_periodicity[0];
2493: if (xLength != 0.0) {
2494: double v0x = coords[0*dim+0];
2496: if (v0x == 0.0) {
2497: v0x = v0[0] = xLength;
2498: }
2499: for(int f = 0; f < dim; f++) {
2500: const double px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
2502: J[0*dim+f] = 0.5*(px - v0x);
2503: }
2504: }
2505: detJ = J[0]*J[3] - J[1]*J[2];
2506: }
2507: PetscLogFlopsNoError(8.0 + 3.0);
2508: }
2509: if (invJ) {
2510: invDet = 1.0/detJ;
2511: invJ[0] = invDet*J[3];
2512: invJ[1] = -invDet*J[1];
2513: invJ[2] = -invDet*J[2];
2514: invJ[3] = invDet*J[0];
2515: PetscLogFlopsNoError(5.0);
2516: }
2517: };
2518: void computeQuadrilateralGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double point[], double v0[], double J[], double invJ[], double& detJ) {
2519: const double *coords = this->restrictClosure(coordinates, e);
2520: const int dim = 2;
2521: double invDet;
2523: if (v0) {
2524: for(int d = 0; d < dim; d++) {
2525: v0[d] = coords[d];
2526: }
2527: }
2528: if (J) {
2529: double x_1 = coords[2] - coords[0];
2530: double y_1 = coords[3] - coords[1];
2531: double x_2 = coords[6] - coords[0];
2532: double y_2 = coords[7] - coords[1];
2533: double x_3 = coords[4] - coords[0];
2534: double y_3 = coords[5] - coords[1];
2536: J[0] = x_1 + (x_3 - x_1 - x_2)*point[1];
2537: J[1] = x_2 + (x_3 - x_1 - x_2)*point[0];
2538: J[2] = y_1 + (y_3 - y_1 - y_2)*point[1];
2539: J[3] = y_1 + (y_3 - y_1 - y_2)*point[0];
2540: detJ = J[0]*J[3] - J[1]*J[2];
2541: PetscLogFlopsNoError(6.0 + 16.0 + 3.0);
2542: }
2543: if (invJ) {
2544: invDet = 1.0/detJ;
2545: invJ[0] = invDet*J[3];
2546: invJ[1] = -invDet*J[1];
2547: invJ[2] = -invDet*J[2];
2548: invJ[3] = invDet*J[0];
2549: PetscLogFlopsNoError(5.0);
2550: }
2551: };
2552: void computeTetrahedronGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2553: const double *coords = this->restrictClosure(coordinates, e);
2554: const int dim = 3;
2555: double invDet;
2557: if (v0) {
2558: for(int d = 0; d < dim; d++) {
2559: v0[d] = coords[d];
2560: }
2561: }
2562: if (J) {
2563: for(int d = 0; d < dim; d++) {
2564: for(int f = 0; f < dim; f++) {
2565: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
2566: }
2567: }
2568: // The minus sign is here since I orient the first face to get the outward normal
2569: detJ = -(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
2570: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
2571: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
2572: PetscLogFlopsNoError(18.0 + 12.0);
2573: }
2574: if (invJ) {
2575: invDet = -1.0/detJ;
2576: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
2577: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
2578: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
2579: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
2580: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
2581: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
2582: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
2583: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
2584: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
2585: PetscLogFlopsNoError(37.0);
2586: }
2587: };
2588: void computeHexahedralGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double point[], double v0[], double J[], double invJ[], double& detJ) {
2589: const double *coords = this->restrictClosure(coordinates, e);
2590: const int dim = 3;
2591: double invDet;
2593: if (v0) {
2594: for(int d = 0; d < dim; d++) {
2595: v0[d] = coords[d];
2596: }
2597: }
2598: if (J) {
2599: double x_1 = coords[3] - coords[0];
2600: double y_1 = coords[4] - coords[1];
2601: double z_1 = coords[5] - coords[2];
2602: double x_2 = coords[6] - coords[0];
2603: double y_2 = coords[7] - coords[1];
2604: double z_2 = coords[8] - coords[2];
2605: double x_3 = coords[9] - coords[0];
2606: double y_3 = coords[10] - coords[1];
2607: double z_3 = coords[11] - coords[2];
2608: double x_4 = coords[12] - coords[0];
2609: double y_4 = coords[13] - coords[1];
2610: double z_4 = coords[14] - coords[2];
2611: double x_5 = coords[15] - coords[0];
2612: double y_5 = coords[16] - coords[1];
2613: double z_5 = coords[17] - coords[2];
2614: double x_6 = coords[18] - coords[0];
2615: double y_6 = coords[19] - coords[1];
2616: double z_6 = coords[20] - coords[2];
2617: double x_7 = coords[21] - coords[0];
2618: double y_7 = coords[22] - coords[1];
2619: double z_7 = coords[23] - coords[2];
2620: double g_x = x_1 - x_2 + x_3 + x_4 - x_5 + x_6 - x_7;
2621: double g_y = y_1 - y_2 + y_3 + y_4 - y_5 + y_6 - y_7;
2622: double g_z = z_1 - z_2 + z_3 + z_4 - z_5 + z_6 - z_7;
2624: J[0] = x_1 + (x_2 - x_1 - x_3)*point[1] + (x_5 - x_1 - x_4)*point[2] + g_x*point[1]*point[2];
2625: J[1] = x_3 + (x_2 - x_1 - x_3)*point[0] + (x_7 - x_3 - x_4)*point[2] + g_x*point[2]*point[0];
2626: J[2] = x_4 + (x_7 - x_3 - x_4)*point[1] + (x_5 - x_1 - x_4)*point[0] + g_x*point[0]*point[1];
2627: J[3] = y_1 + (y_2 - y_1 - y_3)*point[1] + (y_5 - y_1 - y_4)*point[2] + g_y*point[1]*point[2];
2628: J[4] = y_3 + (y_2 - y_1 - y_3)*point[0] + (y_7 - y_3 - y_4)*point[2] + g_y*point[2]*point[0];
2629: J[5] = y_4 + (y_7 - y_3 - y_4)*point[1] + (y_5 - y_1 - y_4)*point[0] + g_y*point[0]*point[1];
2630: J[6] = z_1 + (z_2 - z_1 - z_3)*point[1] + (z_5 - z_1 - z_4)*point[2] + g_z*point[1]*point[2];
2631: J[7] = z_3 + (z_2 - z_1 - z_3)*point[0] + (z_7 - z_3 - z_4)*point[2] + g_z*point[2]*point[0];
2632: J[8] = z_4 + (z_7 - z_3 - z_4)*point[1] + (z_5 - z_1 - z_4)*point[0] + g_z*point[0]*point[1];
2633: detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
2634: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
2635: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
2636: PetscLogFlopsNoError(39.0 + 81.0 + 12.0);
2637: }
2638: if (invJ) {
2639: invDet = 1.0/detJ;
2640: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
2641: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
2642: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
2643: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
2644: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
2645: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
2646: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
2647: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
2648: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
2649: PetscLogFlopsNoError(37.0);
2650: }
2651: };
2652: void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2653: if (this->_dim == 2) {
2654: computeTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
2655: } else if (this->_dim == 3) {
2656: computeTetrahedronGeometry(coordinates, e, v0, J, invJ, detJ);
2657: } else {
2658: throw ALE::Exception("Unsupported dimension for element geometry computation");
2659: }
2660: };
2661: void computeLineFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2662: const arrow_section_type::point_type arrow(cell, face);
2663: const bool reversed = (this->getArrowSection("orientation")->restrictPoint(arrow)[0] == -2);
2664: const int dim = this->getDimension();
2665: double norm = 0.0;
2666: double *vec = tangent;
2668: if (f == 0) {
2669: vec[0] = 0.0; vec[1] = -1.0;
2670: } else if (f == 1) {
2671: vec[0] = 0.70710678; vec[1] = 0.70710678;
2672: } else if (f == 2) {
2673: vec[0] = -1.0; vec[1] = 0.0;
2674: }
2675: for(int d = 0; d < dim; ++d) {
2676: normal[d] = 0.0;
2677: for(int e = 0; e < dim; ++e) normal[d] += cellInvJ[e*dim+d]*vec[e];
2678: if (reversed) normal[d] = -normal[d];
2679: norm += normal[d]*normal[d];
2680: }
2681: norm = std::sqrt(norm);
2682: for(int d = 0; d < dim; ++d) {
2683: normal[d] /= norm;
2684: }
2685: tangent[0] = normal[1];
2686: tangent[1] = -normal[0];
2687: if (this->debug()) {
2688: std::cout << "Cell: " << cell << " Face: " << face << "("<<f<<")" << std::endl;
2689: for(int d = 0; d < dim; ++d) {
2690: std::cout << "Normal["<<d<<"]: " << normal[d] << " Tangent["<<d<<"]: " << tangent[d] << std::endl;
2691: }
2692: }
2693: // Now get 1D Jacobian info
2694: // Should be a way to get this directly
2695: const double *coords = this->restrictClosure(this->getRealSection("coordinates"), face);
2696: detJ = std::sqrt(PetscSqr(coords[1*2+0] - coords[0*2+0]) + PetscSqr(coords[1*2+1] - coords[0*2+1]))/2.0;
2697: invJ[0] = 1.0/detJ;
2698: };
2699: void computeTriangleFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2700: const arrow_section_type::point_type arrow(cell, face);
2701: const bool reversed = this->getArrowSection("orientation")->restrictPoint(arrow)[0] < 0;
2702: const int dim = this->getDimension();
2703: const int faceDim = dim-1;
2704: double norm = 0.0;
2705: double *vec = tangent;
2707: if (f == 0) {
2708: vec[0] = 0.0; vec[1] = 0.0; vec[2] = -1.0;
2709: } else if (f == 1) {
2710: vec[0] = 0.0; vec[1] = -1.0; vec[2] = 0.0;
2711: } else if (f == 2) {
2712: vec[0] = 0.57735027; vec[1] = 0.57735027; vec[2] = 0.57735027;
2713: } else if (f == 3) {
2714: vec[0] = -1.0; vec[1] = 0.0; vec[2] = 0.0;
2715: }
2716: for(int d = 0; d < dim; ++d) {
2717: normal[d] = 0.0;
2718: for(int e = 0; e < dim; ++e) normal[d] += cellInvJ[e*dim+d]*vec[e];
2719: if (reversed) normal[d] = -normal[d];
2720: norm += normal[d]*normal[d];
2721: }
2722: norm = std::sqrt(norm);
2723: for(int d = 0; d < dim; ++d) {
2724: normal[d] /= norm;
2725: }
2726: // Get tangents
2727: tangent[0] = normal[1] - normal[2];
2728: tangent[1] = normal[2] - normal[0];
2729: tangent[2] = normal[0] - normal[1];
2730: norm = 0.0;
2731: for(int d = 0; d < dim; ++d) {
2732: norm += tangent[d]*tangent[d];
2733: }
2734: norm = std::sqrt(norm);
2735: for(int d = 0; d < dim; ++d) {
2736: tangent[d] /= norm;
2737: }
2738: tangent[3] = normal[1]*tangent[2] - normal[2]*tangent[1];
2739: tangent[4] = normal[2]*tangent[0] - normal[0]*tangent[2];
2740: tangent[5] = normal[0]*tangent[1] - normal[1]*tangent[0];
2741: if (this->debug()) {
2742: std::cout << "Cell: " << cell << " Face: " << face << "("<<f<<")" << std::endl;
2743: for(int d = 0; d < dim; ++d) {
2744: std::cout << "Normal["<<d<<"]: " << normal[d] << " TangentA["<<d<<"]: " << tangent[d] << " TangentB["<<d<<"]: " << tangent[dim+d] << std::endl;
2745: }
2746: }
2747: // Now get 2D Jacobian info
2748: // Should be a way to get this directly
2749: const double *coords = this->restrictClosure(this->getRealSection("coordinates"), face);
2750: // Rotate so that normal in z
2751: double invR[9], R[9];
2752: double detR, invDetR;
2753: for(int d = 0; d < dim; d++) {
2754: invR[d*dim+0] = tangent[d];
2755: invR[d*dim+1] = tangent[dim+d];
2756: invR[d*dim+2] = normal[d];
2757: }
2758: invDetR = (invR[0*3+0]*(invR[1*3+1]*invR[2*3+2] - invR[1*3+2]*invR[2*3+1]) +
2759: invR[0*3+1]*(invR[1*3+2]*invR[2*3+0] - invR[1*3+0]*invR[2*3+2]) +
2760: invR[0*3+2]*(invR[1*3+0]*invR[2*3+1] - invR[1*3+1]*invR[2*3+0]));
2761: detR = 1.0/invDetR;
2762: R[0*3+0] = detR*(invR[1*3+1]*invR[2*3+2] - invR[1*3+2]*invR[2*3+1]);
2763: R[0*3+1] = detR*(invR[0*3+2]*invR[2*3+1] - invR[0*3+1]*invR[2*3+2]);
2764: R[0*3+2] = detR*(invR[0*3+1]*invR[1*3+2] - invR[0*3+2]*invR[1*3+1]);
2765: R[1*3+0] = detR*(invR[1*3+2]*invR[2*3+0] - invR[1*3+0]*invR[2*3+2]);
2766: R[1*3+1] = detR*(invR[0*3+0]*invR[2*3+2] - invR[0*3+2]*invR[2*3+0]);
2767: R[1*3+2] = detR*(invR[0*3+2]*invR[1*3+0] - invR[0*3+0]*invR[1*3+2]);
2768: R[2*3+0] = detR*(invR[1*3+0]*invR[2*3+1] - invR[1*3+1]*invR[2*3+0]);
2769: R[2*3+1] = detR*(invR[0*3+1]*invR[2*3+0] - invR[0*3+0]*invR[2*3+1]);
2770: R[2*3+2] = detR*(invR[0*3+0]*invR[1*3+1] - invR[0*3+1]*invR[1*3+0]);
2771: for(int d = 0; d < dim; d++) {
2772: for(int e = 0; e < dim; e++) {
2773: invR[d*dim+e] = 0.0;
2774: for(int g = 0; g < dim; g++) {
2775: invR[d*dim+e] += R[e*dim+g]*coords[d*dim+g];
2776: }
2777: }
2778: }
2779: for(int d = dim-1; d >= 0; --d) {
2780: invR[d*dim+2] -= invR[0*dim+2];
2781: if (this->debug() && (d == dim-1)) {
2782: double ref[9];
2783: for(int q = 0; q < dim; q++) {
2784: for(int e = 0; e < dim; e++) {
2785: ref[q*dim+e] = 0.0;
2786: for(int g = 0; g < dim; g++) {
2787: ref[q*dim+e] += cellInvJ[e*dim+g]*coords[q*dim+g];
2788: }
2789: }
2790: }
2791: std::cout << "f: " << f << std::endl;
2792: std::cout << this->printMatrix(std::string("coords"), dim, dim, coords) << std::endl;
2793: std::cout << this->printMatrix(std::string("ref coords"), dim, dim, ref) << std::endl;
2794: std::cout << this->printMatrix(std::string("R"), dim, dim, R) << std::endl;
2795: std::cout << this->printMatrix(std::string("invR"), dim, dim, invR) << std::endl;
2796: }
2797: if (fabs(invR[d*dim+2]) > 1.0e-8) {
2798: throw ALE::Exception("Invalid rotation");
2799: }
2800: }
2801: double J[4];
2802: for(int d = 0; d < faceDim; d++) {
2803: for(int e = 0; e < faceDim; e++) {
2804: J[d*faceDim+e] = 0.5*(invR[(e+1)*dim+d] - invR[0*dim+d]);
2805: }
2806: }
2807: detJ = fabs(J[0]*J[3] - J[1]*J[2]);
2808: // Probably need something here if detJ < 0
2809: const double invDet = 1.0/detJ;
2810: invJ[0] = invDet*J[3];
2811: invJ[1] = -invDet*J[1];
2812: invJ[2] = -invDet*J[2];
2813: invJ[3] = invDet*J[0];
2814: };
2815: void computeFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2816: if (this->_dim == 2) {
2817: computeLineFaceGeometry(cell, face, f, cellInvJ, invJ, detJ, normal, tangent);
2818: } else if (this->_dim == 3) {
2819: computeTriangleFaceGeometry(cell, face, f, cellInvJ, invJ, detJ, normal, tangent);
2820: } else {
2821: throw ALE::Exception("Unsupported dimension for element geometry computation");
2822: }
2823: };
2824: double getMaxVolume() {
2825: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2826: const Obj<label_sequence>& cells = this->heightStratum(0);
2827: const int dim = this->getDimension();
2828: double v0[3], J[9], invJ[9], detJ, refVolume = 0.0, maxVolume = 0.0;
2830: if (dim == 1) refVolume = 2.0;
2831: if (dim == 2) refVolume = 2.0;
2832: if (dim == 3) refVolume = 4.0/3.0;
2833: for(label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2834: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2835: maxVolume = std::max(maxVolume, detJ*refVolume);
2836: }
2837: return maxVolume;
2838: };
2839: // Find the cell in which this point lies (stupid algorithm)
2840: point_type locatePoint_2D(const real_section_type::value_type point[]) {
2841: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2842: const Obj<label_sequence>& cells = this->heightStratum(0);
2843: const int embedDim = 2;
2844: double v0[2], J[4], invJ[4], detJ;
2846: for(label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2847: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2848: double xi = invJ[0*embedDim+0]*(point[0] - v0[0]) + invJ[0*embedDim+1]*(point[1] - v0[1]);
2849: double eta = invJ[1*embedDim+0]*(point[0] - v0[0]) + invJ[1*embedDim+1]*(point[1] - v0[1]);
2851: if ((xi >= 0.0) && (eta >= 0.0) && (xi + eta <= 2.0)) {
2852: return *c_iter;
2853: }
2854: }
2855: throw ALE::Exception("Could not locate point");
2856: };
2857: // Assume a simplex and 3D
2858: point_type locatePoint_3D(const real_section_type::value_type point[]) {
2859: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2860: const Obj<label_sequence>& cells = this->heightStratum(0);
2861: const int embedDim = 3;
2862: double v0[3], J[9], invJ[9], detJ;
2864: for(label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2865: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2866: double xi = invJ[0*embedDim+0]*(point[0] - v0[0]) + invJ[0*embedDim+1]*(point[1] - v0[1]) + invJ[0*embedDim+2]*(point[2] - v0[2]);
2867: double eta = invJ[1*embedDim+0]*(point[0] - v0[0]) + invJ[1*embedDim+1]*(point[1] - v0[1]) + invJ[1*embedDim+2]*(point[2] - v0[2]);
2868: double zeta = invJ[2*embedDim+0]*(point[0] - v0[0]) + invJ[2*embedDim+1]*(point[1] - v0[1]) + invJ[2*embedDim+2]*(point[2] - v0[2]);
2870: if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) {
2871: return *c_iter;
2872: }
2873: }
2874: throw ALE::Exception("Could not locate point");
2875: };
2876: point_type locatePoint(const real_section_type::value_type point[], point_type guess = -1) {
2877: //guess overrides this by saying that we already know the relation of this point to this mesh. We will need to make it a more robust "guess" later for more than P1
2878: if (guess != -1) {
2879: return guess;
2880: }else if (this->_dim == 2) {
2881: return locatePoint_2D(point);
2882: } else if (this->_dim == 3) {
2883: return locatePoint_3D(point);
2884: } else {
2885: throw ALE::Exception("No point location for mesh dimension");
2886: }
2887: };
2888: public: // Discretization
2889: void markBoundaryCells(const std::string& name, const int marker = 1, const int newMarker = 2, const bool onlyVertices = false) {
2890: const Obj<label_type>& label = this->getLabel(name);
2891: const Obj<label_sequence>& boundary = this->getLabelStratum(name, marker);
2892: const Obj<sieve_type>& sieve = this->getSieve();
2894: if (!onlyVertices) {
2895: const label_sequence::iterator end = boundary->end();
2897: for(label_sequence::iterator e_iter = boundary->begin(); e_iter != end; ++e_iter) {
2898: if (this->height(*e_iter) == 1) {
2899: const point_type cell = *sieve->support(*e_iter)->begin();
2901: this->setValue(label, cell, newMarker);
2902: }
2903: }
2904: } else {
2905: const label_sequence::iterator end = boundary->end();
2906: const int depth = this->depth();
2908: for(label_sequence::iterator v_iter = boundary->begin(); v_iter != end; ++v_iter) {
2909: const Obj<supportArray> support = sieve->nSupport(*v_iter, depth);
2910: const supportArray::iterator sEnd = support->end();
2912: for(supportArray::iterator c_iter = support->begin(); c_iter != sEnd; ++c_iter) {
2913: const Obj<sieve_type::traits::coneSequence>& cone = sieve->cone(*c_iter);
2914: const sieve_type::traits::coneSequence::iterator cEnd = cone->end();
2916: for(sieve_type::traits::coneSequence::iterator e_iter = cone->begin(); e_iter != cEnd; ++e_iter) {
2917: if (sieve->support(*e_iter)->size() == 1) {
2918: this->setValue(label, *c_iter, newMarker);
2919: break;
2920: }
2921: }
2922: }
2923: }
2924: }
2925: };
2926: int setFiberDimensions(const Obj<real_section_type>& s, const Obj<names_type>& discs, names_type& bcLabels) {
2927: const int debug = this->debug();
2928: int maxDof = 0;
2930: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
2931: s->addSpace();
2932: }
2933: for(int d = 0; d <= this->_dim; ++d) {
2934: int numDof = 0;
2935: int f = 0;
2937: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2938: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2939: const int sDof = disc->getNumDof(d);
2941: numDof += sDof;
2942: if (sDof) s->setFiberDimension(this->depthStratum(d), sDof, f);
2943: }
2944: if (numDof) s->setFiberDimension(this->depthStratum(d), numDof);
2945: maxDof = std::max(maxDof, numDof);
2946: }
2947: // Process exclusions
2948: int f = 0;
2950: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2951: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2952: std::string labelName = "exclude-"+*f_iter;
2953: std::set<point_type> seen;
2955: if (this->hasLabel(labelName)) {
2956: const Obj<label_type>& label = this->getLabel(labelName);
2957: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
2958: const label_sequence::iterator end = exclusion->end();
2959: if (debug > 1) {label->view(labelName.c_str());}
2961: for(label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2962: const Obj<coneArray> closure = ALE::SieveAlg<ALE::Mesh>::closure(this, this->getArrowSection("orientation"), *e_iter);
2963: const coneArray::iterator cEnd = closure->end();
2965: for(coneArray::iterator c_iter = closure->begin(); c_iter != cEnd; ++c_iter) {
2966: if (seen.find(*c_iter) != seen.end()) continue;
2967: if (this->getValue(label, *c_iter) == 1) {
2968: seen.insert(*c_iter);
2969: s->setFiberDimension(*c_iter, 0, f);
2970: s->addFiberDimension(*c_iter, -disc->getNumDof(this->depth(*c_iter)));
2971: if (debug > 1) {std::cout << " cell: " << *c_iter << " dim: " << disc->getNumDof(this->depth(*c_iter)) << std::endl;}
2972: }
2973: }
2974: }
2975: }
2976: }
2977: // Process constraints
2978: f = 0;
2979: for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2980: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2981: const Obj<std::set<std::string> > bcs = disc->getBoundaryConditions();
2982: std::string excludeName = "exclude-"+*f_iter;
2984: for(std::set<std::string>::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter) {
2985: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
2986: const Obj<label_sequence>& boundary = this->getLabelStratum(bc->getLabelName(), bc->getMarker());
2988: bcLabels.insert(bc->getLabelName());
2989: if (this->hasLabel(excludeName)) {
2990: const Obj<label_type>& label = this->getLabel(excludeName);
2992: for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2993: if (!this->getValue(label, *e_iter)) {
2994: const int numDof = disc->getNumDof(this->depth(*e_iter));
2996: if (numDof) s->addConstraintDimension(*e_iter, numDof);
2997: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2998: }
2999: }
3000: } else {
3001: for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
3002: const int numDof = disc->getNumDof(this->depth(*e_iter));
3004: if (numDof) s->addConstraintDimension(*e_iter, numDof);
3005: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
3006: }
3007: }
3008: }
3009: }
3010: return maxDof;
3011: };
3012: void calculateIndices() {
3013: // Should have an iterator over the whole tree
3014: Obj<names_type> discs = this->getDiscretizations();
3015: Obj<Mesh> mesh = this;
3016: const int debug = this->debug();
3017: std::map<std::string, std::pair<int, int*> > indices;
3019: mesh.addRef();
3020: for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3021: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3023: indices[*d_iter] = std::pair<int, int*>(0, new int[disc->size(mesh)]);
3024: disc->setIndices(indices[*d_iter].second);
3025: }
3026: const Obj<label_sequence>& cells = this->heightStratum(0);
3027: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *cells->begin());
3028: const coneArray::iterator end = closure->end();
3029: int offset = 0;
3031: if (debug > 1) {std::cout << "Closure for first element" << std::endl;}
3032: for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3033: const int dim = this->depth(*cl_iter);
3035: if (debug > 1) {std::cout << " point " << *cl_iter << " depth " << dim << std::endl;}
3036: for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3037: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3038: const int num = disc->getNumDof(dim);
3040: if (debug > 1) {std::cout << " disc " << disc->getName() << " numDof " << num << std::endl;}
3041: for(int o = 0; o < num; ++o) {
3042: indices[*d_iter].second[indices[*d_iter].first++] = offset++;
3043: }
3044: }
3045: }
3046: if (debug > 1) {
3047: for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3048: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3050: std::cout << "Discretization " << disc->getName() << " indices:";
3051: for(int i = 0; i < indices[*d_iter].first; ++i) {
3052: std::cout << " " << indices[*d_iter].second[i];
3053: }
3054: std::cout << std::endl;
3055: }
3056: }
3057: };
3058: void calculateIndicesExcluded(const Obj<real_section_type>& s, const Obj<names_type>& discs) {
3059: typedef std::map<std::string, std::pair<int, indexSet> > indices_type;
3060: const Obj<label_type>& indexLabel = this->createLabel("cellExclusion");
3061: const int debug = this->debug();
3062: int marker = 0;
3063: Obj<Mesh> mesh = this;
3064: std::map<indices_type, int> indexMap;
3065: indices_type indices;
3067: mesh.addRef();
3068: for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3069: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3070: const int size = disc->size(mesh);
3072: indices[*d_iter].second.resize(size);
3073: }
3074: const names_type::const_iterator dBegin = discs->begin();
3075: const names_type::const_iterator dEnd = discs->end();
3076: std::set<point_type> seen;
3077: int f = 0;
3079: for(names_type::const_iterator f_iter = dBegin; f_iter != dEnd; ++f_iter, ++f) {
3080: std::string labelName = "exclude-"+*f_iter;
3082: if (this->hasLabel(labelName)) {
3083: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
3084: const label_sequence::iterator end = exclusion->end();
3086: if (debug > 1) {std::cout << "Processing exclusion " << labelName << std::endl;}
3087: for(label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
3088: if (this->height(*e_iter)) continue;
3089: const Obj<coneArray> closure = ALE::SieveAlg<ALE::Mesh>::closure(this, this->getArrowSection("orientation"), *e_iter);
3090: const coneArray::iterator clEnd = closure->end();
3091: int offset = 0;
3093: if (debug > 1) {std::cout << " Closure for cell " << *e_iter << std::endl;}
3094: for(coneArray::iterator cl_iter = closure->begin(); cl_iter != clEnd; ++cl_iter) {
3095: int g = 0;
3097: if (debug > 1) {std::cout << " point " << *cl_iter << std::endl;}
3098: for(names_type::const_iterator g_iter = dBegin; g_iter != dEnd; ++g_iter, ++g) {
3099: const int fDim = s->getFiberDimension(*cl_iter, g);
3101: if (debug > 1) {std::cout << " disc " << *g_iter << " numDof " << fDim << std::endl;}
3102: for(int d = 0; d < fDim; ++d) {
3103: indices[*g_iter].second[indices[*g_iter].first++] = offset++;
3104: }
3105: }
3106: }
3107: const std::map<indices_type, int>::iterator entry = indexMap.find(indices);
3109: if (debug > 1) {
3110: for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
3111: for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3112: std::cout << "Discretization (" << i_iter->second << ") " << *g_iter << " indices:";
3113: for(int i = 0; i < ((indices_type) i_iter->first)[*g_iter].first; ++i) {
3114: std::cout << " " << ((indices_type) i_iter->first)[*g_iter].second[i];
3115: }
3116: std::cout << std::endl;
3117: }
3118: std::cout << "Comparison: " << (indices == i_iter->first) << std::endl;
3119: }
3120: for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3121: std::cout << "Discretization " << *g_iter << " indices:";
3122: for(int i = 0; i < indices[*g_iter].first; ++i) {
3123: std::cout << " " << indices[*g_iter].second[i];
3124: }
3125: std::cout << std::endl;
3126: }
3127: }
3128: if (entry != indexMap.end()) {
3129: this->setValue(indexLabel, *e_iter, entry->second);
3130: if (debug > 1) {std::cout << " Found existing indices with marker " << entry->second << std::endl;}
3131: } else {
3132: indexMap[indices] = ++marker;
3133: this->setValue(indexLabel, *e_iter, marker);
3134: if (debug > 1) {std::cout << " Created new indices with marker " << marker << std::endl;}
3135: }
3136: for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3137: indices[*g_iter].first = 0;
3138: for(unsigned int i = 0; i < indices[*g_iter].second.size(); ++i) indices[*g_iter].second[i] = 0;
3139: }
3140: }
3141: }
3142: }
3143: if (debug > 1) {indexLabel->view("cellExclusion");}
3144: for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
3145: if (debug > 1) {std::cout << "Setting indices for marker " << i_iter->second << std::endl;}
3146: for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3147: const Obj<Discretization>& disc = this->getDiscretization(*g_iter);
3148: const indexSet indSet = ((indices_type) i_iter->first)[*g_iter].second;
3149: const int size = indSet.size();
3150: int *_indices = new int[size];
3152: if (debug > 1) {std::cout << " field " << *g_iter << std::endl;}
3153: for(int i = 0; i < size; ++i) {
3154: _indices[i] = indSet[i];
3155: if (debug > 1) {std::cout << " indices["<<i<<"] = " << _indices[i] << std::endl;}
3156: }
3157: disc->setIndices(_indices, i_iter->second);
3158: }
3159: }
3160: };
3161: void setupField(const Obj<real_section_type>& s, const int cellMarker = 2, const bool noUpdate = false) {
3162: const Obj<names_type>& discs = this->getDiscretizations();
3163: const int debug = s->debug();
3164: names_type bcLabels;
3165: int maxDof;
3167: maxDof = this->setFiberDimensions(s, discs, bcLabels);
3168: this->calculateIndices();
3169: this->calculateIndicesExcluded(s, discs);
3170: this->allocate(s);
3171: s->defaultConstraintDof();
3172: const Obj<label_type>& cellExclusion = this->getLabel("cellExclusion");
3174: if (debug > 1) {std::cout << "Setting boundary values" << std::endl;}
3175: for(names_type::const_iterator n_iter = bcLabels.begin(); n_iter != bcLabels.end(); ++n_iter) {
3176: const Obj<label_sequence>& boundaryCells = this->getLabelStratum(*n_iter, cellMarker);
3177: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
3178: const Obj<names_type>& discs = this->getDiscretizations();
3179: const point_type firstCell = *boundaryCells->begin();
3180: const int numFields = discs->size();
3181: real_section_type::value_type *values = new real_section_type::value_type[this->sizeWithBC(s, firstCell)];
3182: int *dofs = new int[maxDof];
3183: int *v = new int[numFields];
3184: double *v0 = new double[this->getDimension()];
3185: double *J = new double[this->getDimension()*this->getDimension()];
3186: double detJ;
3188: for(label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
3189: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
3190: const coneArray::iterator end = closure->end();
3192: if (debug > 1) {std::cout << " Boundary cell " << *c_iter << std::endl;}
3193: this->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
3194: for(int f = 0; f < numFields; ++f) v[f] = 0;
3195: for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3196: const int cDim = s->getConstraintDimension(*cl_iter);
3197: int off = 0;
3198: int f = 0;
3199: int i = -1;
3201: if (debug > 1) {std::cout << " point " << *cl_iter << std::endl;}
3202: if (cDim) {
3203: if (debug > 1) {std::cout << " constrained excMarker: " << this->getValue(cellExclusion, *c_iter) << std::endl;}
3204: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3205: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3206: const Obj<names_type> bcs = disc->getBoundaryConditions();
3207: const int fDim = s->getFiberDimension(*cl_iter, f);//disc->getNumDof(this->depth(*cl_iter));
3208: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
3209: int b = 0;
3211: for(names_type::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter, ++b) {
3212: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
3213: const int value = this->getValue(this->getLabel(bc->getLabelName()), *cl_iter);
3215: if (b > 0) v[f] -= fDim;
3216: if (value == bc->getMarker()) {
3217: if (debug > 1) {std::cout << " field " << *f_iter << " marker " << value << std::endl;}
3218: for(int d = 0; d < fDim; ++d, ++v[f]) {
3219: dofs[++i] = off+d;
3220: if (!noUpdate) values[indices[v[f]]] = (*bc->getDualIntegrator())(v0, J, v[f], bc->getFunction());
3221: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3222: }
3223: // Allow only one condition per point
3224: ++b;
3225: break;
3226: } else {
3227: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
3228: for(int d = 0; d < fDim; ++d, ++v[f]) {
3229: values[indices[v[f]]] = 0.0;
3230: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3231: }
3232: }
3233: }
3234: if (b == 0) {
3235: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
3236: for(int d = 0; d < fDim; ++d, ++v[f]) {
3237: values[indices[v[f]]] = 0.0;
3238: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3239: }
3240: }
3241: off += fDim;
3242: }
3243: if (i != cDim-1) {throw ALE::Exception("Invalid constraint initialization");}
3244: s->setConstraintDof(*cl_iter, dofs);
3245: } else {
3246: if (debug > 1) {std::cout << " unconstrained" << std::endl;}
3247: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3248: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3249: const int fDim = s->getFiberDimension(*cl_iter, f);//disc->getNumDof(this->depth(*cl_iter));
3250: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
3252: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
3253: for(int d = 0; d < fDim; ++d, ++v[f]) {
3254: values[indices[v[f]]] = 0.0;
3255: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3256: }
3257: }
3258: }
3259: }
3260: if (debug > 1) {
3261: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
3262: const coneArray::iterator end = closure->end();
3264: for(int f = 0; f < numFields; ++f) v[f] = 0;
3265: for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3266: int f = 0;
3267: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3268: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3269: const int fDim = s->getFiberDimension(*cl_iter, f);
3270: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
3272: for(int d = 0; d < fDim; ++d, ++v[f]) {
3273: std::cout << " "<<*f_iter<<"-value["<<indices[v[f]]<<"] " << values[indices[v[f]]] << std::endl;
3274: }
3275: }
3276: }
3277: }
3278: if (!noUpdate) {
3279: this->updateAll(s, *c_iter, values);
3280: }
3281: }
3282: delete [] dofs;
3283: delete [] values;
3284: delete [] v0;
3285: delete [] J;
3286: }
3287: if (debug > 1) {s->view("");}
3288: };
3289: public:
3290: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
3291: if (comm == MPI_COMM_NULL) {
3292: comm = this->comm();
3293: }
3294: if (name == "") {
3295: PetscPrintf(comm, "viewing a Mesh\n");
3296: } else {
3297: PetscPrintf(comm, "viewing Mesh '%s'\n", name.c_str());
3298: }
3299: this->getSieve()->view("mesh sieve", comm);
3300: Obj<names_type> sections = this->getRealSections();
3302: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3303: this->getRealSection(*name)->view(*name);
3304: }
3305: sections = this->getIntSections();
3306: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3307: this->getIntSection(*name)->view(*name);
3308: }
3309: sections = this->getArrowSections();
3310: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3311: this->getArrowSection(*name)->view(*name);
3312: }
3313: };
3314: template<typename value_type>
3315: static std::string printMatrix(const std::string& name, const int rows, const int cols, const value_type matrix[], const int rank = -1)
3316: {
3317: ostringstream output;
3318: ostringstream rankStr;
3320: if (rank >= 0) {
3321: rankStr << "[" << rank << "]";
3322: }
3323: output << rankStr.str() << name << " = " << std::endl;
3324: for(int r = 0; r < rows; r++) {
3325: if (r == 0) {
3326: output << rankStr.str() << " /";
3327: } else if (r == rows-1) {
3328: output << rankStr.str() << " \\";
3329: } else {
3330: output << rankStr.str() << " |";
3331: }
3332: for(int c = 0; c < cols; c++) {
3333: output << " " << matrix[r*cols+c];
3334: }
3335: if (r == 0) {
3336: output << " \\" << std::endl;
3337: } else if (r == rows-1) {
3338: output << " /" << std::endl;
3339: } else {
3340: output << " |" << std::endl;
3341: }
3342: }
3343: return output.str();
3344: };
3345: };
3346: template<typename Mesh>
3347: class MeshBuilder {
3348: public:
3351: /*
3352: Simple square boundary:
3354: 18--5-17--4--16
3355: | | |
3356: 6 10 3
3357: | | |
3358: 19-11-20--9--15
3359: | | |
3360: 7 8 2
3361: | | |
3362: 12--0-13--1--14
3363: */
3364: static Obj<Mesh> createSquareBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int edges[], const int debug = 0) {
3365: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3366: int numVertices = (edges[0]+1)*(edges[1]+1);
3367: int numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
3368: double *coords = new double[numVertices*2];
3369: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3370: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3371: int order = 0;
3373: mesh->setSieve(sieve);
3374: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3375: if (mesh->commRank() == 0) {
3376: /* Create sieve and ordering */
3377: for(int v = numEdges; v < numEdges+numVertices; v++) {
3378: vertices[v-numEdges] = typename Mesh::point_type(v);
3379: }
3380: for(int vy = 0; vy <= edges[1]; vy++) {
3381: for(int ex = 0; ex < edges[0]; ex++) {
3382: typename Mesh::point_type edge(vy*edges[0] + ex);
3383: int vertex = vy*(edges[0]+1) + ex;
3385: sieve->addArrow(vertices[vertex+0], edge, order++);
3386: sieve->addArrow(vertices[vertex+1], edge, order++);
3387: if ((vy == 0) || (vy == edges[1])) {
3388: mesh->setValue(markers, edge, 1);
3389: mesh->setValue(markers, vertices[vertex], 1);
3390: if (ex == edges[0]-1) {
3391: mesh->setValue(markers, vertices[vertex+1], 1);
3392: }
3393: }
3394: }
3395: }
3396: for(int vx = 0; vx <= edges[0]; vx++) {
3397: for(int ey = 0; ey < edges[1]; ey++) {
3398: typename Mesh::point_type edge(vx*edges[1] + ey + edges[0]*(edges[1]+1));
3399: int vertex = ey*(edges[0]+1) + vx;
3401: sieve->addArrow(vertices[vertex], edge, order++);
3402: sieve->addArrow(vertices[vertex+edges[0]+1], edge, order++);
3403: if ((vx == 0) || (vx == edges[0])) {
3404: mesh->setValue(markers, edge, 1);
3405: mesh->setValue(markers, vertices[vertex], 1);
3406: if (ey == edges[1]-1) {
3407: mesh->setValue(markers, vertices[vertex+edges[0]+1], 1);
3408: }
3409: }
3410: }
3411: }
3412: }
3413: mesh->stratify();
3414: for(int vy = 0; vy <= edges[1]; ++vy) {
3415: for(int vx = 0; vx <= edges[0]; ++vx) {
3416: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
3417: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
3418: }
3419: }
3420: delete [] vertices;
3421: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3422: return mesh;
3423: };
3426: /*
3427: Simple square boundary:
3429: 14--5-13--4--12
3430: | |
3431: 6 3
3432: | |
3433: 15 11
3434: | |
3435: 7 2
3436: | |
3437: 8--0--9--1--10
3438: */
3439: static Obj<Mesh> createSquareBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int edges, const int debug = 0) {
3440: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3441: int numVertices = edges*4;
3442: int numEdges = edges*4;
3443: double *coords = new double[numVertices*2];
3444: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3445: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3447: mesh->setSieve(sieve);
3448: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3449: if (mesh->commRank() == 0) {
3450: /* Create sieve and ordering */
3451: for(int v = numEdges; v < numEdges+numVertices; v++) {
3452: vertices[v-numEdges] = typename Mesh::point_type(v);
3453: }
3454: for(int e = 0; e < numEdges; ++e) {
3455: typename Mesh::point_type edge(e);
3456: int order = 0;
3458: sieve->addArrow(vertices[e], edge, order++);
3459: sieve->addArrow(vertices[(e+1)%numVertices], edge, order++);
3460: mesh->setValue(markers, edge, 2);
3461: mesh->setValue(markers, vertices[e], 1);
3462: mesh->setValue(markers, vertices[(e+1)%numVertices], 1);
3463: }
3464: }
3465: mesh->stratify();
3466: for(int v = 0; v < edges; ++v) {
3467: coords[(v+edges*0)*2+0] = lower[0] + ((upper[0] - lower[0])/edges)*v;
3468: coords[(v+edges*0)*2+1] = lower[1];
3469: }
3470: for(int v = 0; v < edges; ++v) {
3471: coords[(v+edges*1)*2+0] = upper[0];
3472: coords[(v+edges*1)*2+1] = lower[1] + ((upper[1] - lower[1])/edges)*v;
3473: }
3474: for(int v = 0; v < edges; ++v) {
3475: coords[(v+edges*2)*2+0] = upper[0] - ((upper[0] - lower[0])/edges)*v;
3476: coords[(v+edges*2)*2+1] = upper[1];
3477: }
3478: for(int v = 0; v < edges; ++v) {
3479: coords[(v+edges*3)*2+0] = lower[0];
3480: coords[(v+edges*3)*2+1] = upper[1] - ((upper[1] - lower[1])/edges)*v;
3481: }
3482: delete [] vertices;
3483: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3484: // Build normals for cells
3485: const Obj<typename Mesh::real_section_type>& normals = mesh->getRealSection("normals");
3486: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
3488: //normals->setChart(typename Mesh::real_section_type::chart_type(*std::min_element(cells->begin(), cells->end()),
3489: // *std::max_element(cells->begin(), cells->end())+1));
3490: normals->setFiberDimension(cells, mesh->getDimension()+1);
3491: mesh->allocate(normals);
3492: for(int e = 0; e < edges; ++e) {
3493: double normal[2] = {0.0, -1.0};
3494: normals->updatePoint(e+edges*0, normal);
3495: }
3496: for(int e = 0; e < edges; ++e) {
3497: double normal[2] = {1.0, 0.0};
3498: normals->updatePoint(e+edges*1, normal);
3499: }
3500: for(int e = 0; e < edges; ++e) {
3501: double normal[2] = {0.0, 1.0};
3502: normals->updatePoint(e+edges*2, normal);
3503: }
3504: for(int e = 0; e < edges; ++e) {
3505: double normal[2] = {-1.0, 0.0};
3506: normals->updatePoint(e+edges*3, normal);
3507: }
3508: return mesh;
3509: };
3512: /*
3513: Simple square boundary:
3515: 18--5-17--4--16
3516: | | |
3517: 6 10 3
3518: | | |
3519: 19-11-20--9--15
3520: | | |
3521: 7 8 2
3522: | | |
3523: 12--0-13--1--14
3524: */
3525: static Obj<Mesh> createParticleInSquareBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int edges[], const double radius, const int partEdges, const int debug = 0) {
3526: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3527: const int numSquareVertices = (edges[0]+1)*(edges[1]+1);
3528: const int numVertices = numSquareVertices + partEdges;
3529: const int numSquareEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
3530: const int numEdges = numSquareEdges + partEdges;
3531: double *coords = new double[numVertices*2];
3532: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3533: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3534: int order = 0;
3536: mesh->setSieve(sieve);
3537: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3538: if (mesh->commRank() == 0) {
3539: /* Create sieve and ordering */
3540: for(int v = numEdges; v < numEdges+numVertices; v++) {
3541: vertices[v-numEdges] = typename Mesh::point_type(v);
3542: }
3543: // Make square
3544: for(int vy = 0; vy <= edges[1]; vy++) {
3545: for(int ex = 0; ex < edges[0]; ex++) {
3546: typename Mesh::point_type edge(vy*edges[0] + ex);
3547: int vertex = vy*(edges[0]+1) + ex;
3549: sieve->addArrow(vertices[vertex+0], edge, order++);
3550: sieve->addArrow(vertices[vertex+1], edge, order++);
3551: if ((vy == 0) || (vy == edges[1])) {
3552: mesh->setValue(markers, edge, 1);
3553: mesh->setValue(markers, vertices[vertex], 1);
3554: if (ex == edges[0]-1) {
3555: mesh->setValue(markers, vertices[vertex+1], 1);
3556: }
3557: }
3558: }
3559: }
3560: for(int vx = 0; vx <= edges[0]; vx++) {
3561: for(int ey = 0; ey < edges[1]; ey++) {
3562: typename Mesh::point_type edge(vx*edges[1] + ey + edges[0]*(edges[1]+1));
3563: int vertex = ey*(edges[0]+1) + vx;
3565: sieve->addArrow(vertices[vertex], edge, order++);
3566: sieve->addArrow(vertices[vertex+edges[0]+1], edge, order++);
3567: if ((vx == 0) || (vx == edges[0])) {
3568: mesh->setValue(markers, edge, 1);
3569: mesh->setValue(markers, vertices[vertex], 1);
3570: if (ey == edges[1]-1) {
3571: mesh->setValue(markers, vertices[vertex+edges[0]+1], 1);
3572: }
3573: }
3574: }
3575: }
3576: // Make particle
3577: for(int ep = 0; ep < partEdges; ++ep) {
3578: typename Mesh::point_type edge(numSquareEdges + ep);
3579: const int vertexA = numSquareVertices + ep;
3580: const int vertexB = numSquareVertices + (ep+1)%partEdges;
3582: sieve->addArrow(vertices[vertexA], edge, order++);
3583: sieve->addArrow(vertices[vertexB], edge, order++);
3584: mesh->setValue(markers, edge, 2);
3585: mesh->setValue(markers, vertices[vertexA], 2);
3586: mesh->setValue(markers, vertices[vertexB], 2);
3587: }
3588: }
3589: mesh->stratify();
3590: for(int vy = 0; vy <= edges[1]; ++vy) {
3591: for(int vx = 0; vx <= edges[0]; ++vx) {
3592: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
3593: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
3594: }
3595: }
3596: const double centroidX = 0.5*(upper[0] + lower[0]);
3597: const double centroidY = 0.5*(upper[1] + lower[1]);
3598: for(int vp = 0; vp < partEdges; ++vp) {
3599: const double rad = 2.0*PETSC_PI*vp/partEdges;
3600: coords[(numSquareVertices+vp)*2+0] = centroidX + radius*cos(rad);
3601: coords[(numSquareVertices+vp)*2+1] = centroidY + radius*sin(rad);
3602: }
3603: delete [] vertices;
3604: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3605: return mesh;
3606: };
3609: /*
3610: Simple boundary with reentrant singularity:
3612: 12--5-11
3613: | |
3614: | 4
3615: | |
3616: 6 10--3--9
3617: | |
3618: | 2
3619: | |
3620: 7-----1-----8
3621: */
3622: static Obj<Mesh> createReentrantBoundary(const MPI_Comm comm, const double lower[], const double upper[], double notchpercent[], const int debug = 0) {
3623: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3624: int numVertices = 6;
3625: int numEdges = numVertices;
3626: double *coords = new double[numVertices*2];
3627: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3628: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3630: mesh->setSieve(sieve);
3631: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3632: if (mesh->commRank() == 0) {
3633: /* Create sieve and ordering */
3634: for (int b = 0; b < numVertices; b++) {
3635: sieve->addArrow(numEdges+b, b);
3636: sieve->addArrow(numEdges+b, (b+1)%numVertices);
3637: mesh->setValue(markers, b, 1);
3638: mesh->setValue(markers, b+numVertices, 1);
3639: }
3640: coords[0] = upper[0];
3641: coords[1] = lower[1];
3643: coords[2] = lower[0];
3644: coords[3] = lower[1];
3645:
3646: coords[4] = lower[0];
3647: coords[5] = notchpercent[1]*lower[1] + (1 - notchpercent[1])*upper[1];
3648:
3649: coords[6] = notchpercent[0]*upper[0] + (1 - notchpercent[0])*lower[0];
3650: coords[7] = notchpercent[1]*lower[1] + (1 - notchpercent[1])*upper[1];
3651:
3652:
3653: coords[8] = notchpercent[0]*upper[0] + (1 - notchpercent[0])*lower[0];
3654: coords[9] = upper[1];
3656: coords[10] = upper[0];
3657: coords[11] = upper[1];
3658: mesh->stratify();
3659: }
3660: delete [] vertices;
3661: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3662: return mesh;
3663: }
3667: /*
3668: Circular boundary with reentrant singularity:
3670: ---1
3671: -- |
3672: - |
3673: | |
3674: | 0-----n
3675: | |
3676: - -
3677: -- --
3678: -------
3679: */
3680: static Obj<Mesh> createCircularReentrantBoundary(const MPI_Comm comm, const int segments, const double radius, const double arc_percent, const int debug = 0) {
3681: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3682: int numVertices = segments+2;
3683: int numEdges = numVertices;
3684: double *coords = new double[numVertices*2];
3685: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3686: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3688: mesh->setSieve(sieve);
3689: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3690: if (mesh->commRank() == 0) {
3691: /* Create sieve and ordering */
3693: int startvertex = 1;
3694: if (arc_percent < 1.) {
3695: coords[0] = 0.;
3696: coords[1] = 0.;
3697: } else {
3698: numVertices = segments;
3699: numEdges = numVertices;
3700: startvertex = 0;
3701: }
3703: for (int b = 0; b < numVertices; b++) {
3704: sieve->addArrow(numEdges+b, b);
3705: sieve->addArrow(numEdges+b, (b+1)%numVertices);
3706: mesh->setValue(markers, b, 1);
3707: mesh->setValue(markers, b+numVertices, 1);
3708: }
3710: double anglestep = arc_percent*2.*3.14159265/((float)segments);
3712: for (int i = startvertex; i < numVertices; i++) {
3713: coords[2*i] = radius * sin(anglestep*(i-startvertex));
3714: coords[2*i+1] = radius*cos(anglestep*(i-startvertex));
3715: }
3716: mesh->stratify();
3717: }
3718: delete [] vertices;
3719: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3720: return mesh;
3721: };
3724: static Obj<Mesh> createAnnularBoundary(const MPI_Comm comm, const int segments, const double centers[4], const double radii[2], const int debug = 0) {
3725: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3726: int numVertices = segments*2;
3727: int numEdges = numVertices;
3728: double *coords = new double[numVertices*2];
3729: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3730: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3732: mesh->setSieve(sieve);
3733: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3734: if (mesh->commRank() == 0) {
3735: for (int e = 0; e < segments; ++e) {
3736: sieve->addArrow(numEdges+e, e);
3737: sieve->addArrow(numEdges+(e+1)%segments, e);
3738: sieve->addArrow(numEdges+segments+e, e+segments);
3739: sieve->addArrow(numEdges+segments+(e+1)%segments, e+segments);
3740: mesh->setValue(markers, e, 1);
3741: mesh->setValue(markers, e+segments, 1);
3742: mesh->setValue(markers, e+numEdges, 1);
3743: mesh->setValue(markers, e+numEdges+segments, 1);
3744: }
3745: const double anglestep = 2.0*M_PI/segments;
3747: for (int v = 0; v < segments; ++v) {
3748: coords[v*2] = centers[0] + radii[0]*cos(anglestep*v);
3749: coords[v*2+1] = centers[1] + radii[0]*sin(anglestep*v);
3750: coords[(v+segments)*2] = centers[2] + radii[1]*cos(anglestep*v);
3751: coords[(v+segments)*2+1] = centers[3] + radii[1]*sin(anglestep*v);
3752: }
3753: mesh->addHole(¢ers[2]);
3754: }
3755: mesh->stratify();
3756: delete [] vertices;
3757: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3758: return mesh;
3759: };
3762: /*
3763: Simple cubic boundary:
3765: 30----31-----32
3766: | | |
3767: | 3 | 2 |
3768: | | |
3769: 27----28-----29
3770: | | |
3771: | 0 | 1 |
3772: | | |
3773: 24----25-----26
3774: */
3775: static Obj<Mesh> createCubeBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int faces[], const int debug = 0) {
3776: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
3777: int numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
3778: int numFaces = 6;
3779: double *coords = new double[numVertices*3];
3780: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3781: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3782: int order = 0;
3784: mesh->setSieve(sieve);
3785: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3786: if (mesh->commRank() == 0) {
3787: /* Create sieve and ordering */
3788: for(int v = numFaces; v < numFaces+numVertices; v++) {
3789: vertices[v-numFaces] = typename Mesh::point_type(v);
3790: mesh->setValue(markers, vertices[v-numFaces], 1);
3791: }
3792: {
3793: // Side 0 (Front)
3794: typename Mesh::point_type face(0);
3795: sieve->addArrow(vertices[0], face, order++);
3796: sieve->addArrow(vertices[1], face, order++);
3797: sieve->addArrow(vertices[2], face, order++);
3798: sieve->addArrow(vertices[3], face, order++);
3799: mesh->setValue(markers, face, 1);
3800: }
3801: {
3802: // Side 1 (Back)
3803: typename Mesh::point_type face(1);
3804: sieve->addArrow(vertices[5], face, order++);
3805: sieve->addArrow(vertices[4], face, order++);
3806: sieve->addArrow(vertices[7], face, order++);
3807: sieve->addArrow(vertices[6], face, order++);
3808: mesh->setValue(markers, face, 1);
3809: }
3810: {
3811: // Side 2 (Bottom)
3812: typename Mesh::point_type face(2);
3813: sieve->addArrow(vertices[4], face, order++);
3814: sieve->addArrow(vertices[5], face, order++);
3815: sieve->addArrow(vertices[1], face, order++);
3816: sieve->addArrow(vertices[0], face, order++);
3817: mesh->setValue(markers, face, 1);
3818: }
3819: {
3820: // Side 3 (Top)
3821: typename Mesh::point_type face(3);
3822: sieve->addArrow(vertices[3], face, order++);
3823: sieve->addArrow(vertices[2], face, order++);
3824: sieve->addArrow(vertices[6], face, order++);
3825: sieve->addArrow(vertices[7], face, order++);
3826: mesh->setValue(markers, face, 1);
3827: }
3828: {
3829: // Side 4 (Left)
3830: typename Mesh::point_type face(4);
3831: sieve->addArrow(vertices[4], face, order++);
3832: sieve->addArrow(vertices[0], face, order++);
3833: sieve->addArrow(vertices[3], face, order++);
3834: sieve->addArrow(vertices[7], face, order++);
3835: mesh->setValue(markers, face, 1);
3836: }
3837: {
3838: // Side 5 (Right)
3839: typename Mesh::point_type face(5);
3840: sieve->addArrow(vertices[1], face, order++);
3841: sieve->addArrow(vertices[5], face, order++);
3842: sieve->addArrow(vertices[6], face, order++);
3843: sieve->addArrow(vertices[2], face, order++);
3844: mesh->setValue(markers, face, 1);
3845: }
3846: }
3847: mesh->stratify();
3848: #if 0
3849: for(int vz = 0; vz <= edges[2]; ++vz) {
3850: for(int vy = 0; vy <= edges[1]; ++vy) {
3851: for(int vx = 0; vx <= edges[0]; ++vx) {
3852: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
3853: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
3854: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
3855: }
3856: }
3857: }
3858: #else
3859: coords[0*3+0] = lower[0];
3860: coords[0*3+1] = lower[1];
3861: coords[0*3+2] = upper[2];
3862: coords[1*3+0] = upper[0];
3863: coords[1*3+1] = lower[1];
3864: coords[1*3+2] = upper[2];
3865: coords[2*3+0] = upper[0];
3866: coords[2*3+1] = upper[1];
3867: coords[2*3+2] = upper[2];
3868: coords[3*3+0] = lower[0];
3869: coords[3*3+1] = upper[1];
3870: coords[3*3+2] = upper[2];
3871: coords[4*3+0] = lower[0];
3872: coords[4*3+1] = lower[1];
3873: coords[4*3+2] = lower[2];
3874: coords[5*3+0] = upper[0];
3875: coords[5*3+1] = lower[1];
3876: coords[5*3+2] = lower[2];
3877: coords[6*3+0] = upper[0];
3878: coords[6*3+1] = upper[1];
3879: coords[6*3+2] = lower[2];
3880: coords[7*3+0] = lower[0];
3881: coords[7*3+1] = upper[1];
3882: coords[7*3+2] = lower[2];
3883: #endif
3884: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3885: return mesh;
3886: };
3888: // Creates a triangular prism boundary
3889: static Obj<Mesh> createPrismBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int debug = 0) {
3890: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
3891: int numVertices = 6;
3892: int numFaces = 5;
3893: double *coords = new double[numVertices*3];
3894: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3895: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3896: int order = 0;
3898: mesh->setSieve(sieve);
3899: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3900: if (mesh->commRank() == 0) {
3901: /* Create sieve and ordering */
3902: for(int v = numFaces; v < numFaces+numVertices; v++) {
3903: vertices[v-numFaces] = typename Mesh::point_type(v);
3904: mesh->setValue(markers, vertices[v-numFaces], 1);
3905: }
3906: {
3907: // Side 0 (Top)
3908: typename Mesh::point_type face(0);
3909: sieve->addArrow(vertices[0], face, order++);
3910: sieve->addArrow(vertices[1], face, order++);
3911: sieve->addArrow(vertices[2], face, order++);
3912: mesh->setValue(markers, face, 1);
3913: }
3914: {
3915: // Side 1 (Bottom)
3916: typename Mesh::point_type face(1);
3917: sieve->addArrow(vertices[3], face, order++);
3918: sieve->addArrow(vertices[5], face, order++);
3919: sieve->addArrow(vertices[4], face, order++);
3920: mesh->setValue(markers, face, 1);
3921: }
3922: {
3923: // Side 2 (Front)
3924: typename Mesh::point_type face(2);
3925: sieve->addArrow(vertices[0], face, order++);
3926: sieve->addArrow(vertices[3], face, order++);
3927: sieve->addArrow(vertices[4], face, order++);
3928: sieve->addArrow(vertices[1], face, order++);
3929: mesh->setValue(markers, face, 1);
3930: }
3931: {
3932: // Side 3 (Left)
3933: typename Mesh::point_type face(3);
3934: sieve->addArrow(vertices[1], face, order++);
3935: sieve->addArrow(vertices[4], face, order++);
3936: sieve->addArrow(vertices[5], face, order++);
3937: sieve->addArrow(vertices[2], face, order++);
3938: mesh->setValue(markers, face, 1);
3939: }
3940: {
3941: // Side 4 (Right)
3942: typename Mesh::point_type face(4);
3943: sieve->addArrow(vertices[2], face, order++);
3944: sieve->addArrow(vertices[5], face, order++);
3945: sieve->addArrow(vertices[3], face, order++);
3946: sieve->addArrow(vertices[0], face, order++);
3947: mesh->setValue(markers, face, 1);
3948: }
3949: }
3950: mesh->stratify();
3951: coords[0*3+0] = lower[0];
3952: coords[0*3+1] = lower[1];
3953: coords[0*3+2] = upper[2];
3954: coords[1*3+0] = upper[0];
3955: coords[1*3+1] = lower[1];
3956: coords[1*3+2] = upper[2];
3957: coords[2*3+0] = upper[0];
3958: coords[2*3+1] = upper[1];
3959: coords[2*3+2] = upper[2];
3960: coords[3*3+0] = lower[0];
3961: coords[3*3+1] = upper[1];
3962: coords[3*3+2] = upper[2];
3963: coords[4*3+0] = lower[0];
3964: coords[4*3+1] = lower[1];
3965: coords[4*3+2] = lower[2];
3966: coords[5*3+0] = upper[0];
3967: coords[5*3+1] = lower[1];
3968: coords[5*3+2] = lower[2];
3969: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3970: return mesh;
3971: };
3975: /* v0
3976: / \
3977: / \
3978: 2 / 4 \ 1
3979: / \
3980: / v12 \
3981: v6|\ /|\ /|v5
3982: | \v8 | v7/ | z
3983: | |7 |8 | | |
3984: | |v13\ | | <-v4 / \
3985: | v9/ 9 \v10| x y
3986: v1 | 5 \ / 6 |v2
3987: \ \ / /
3988: \ v11 /
3989: \ | /
3990: 3 \ | /
3991: \|/
3992: v3
3993: */
3994: static Obj<Mesh> createFicheraCornerBoundary(const MPI_Comm comm, const double lower[], const double upper[], const double offset[], const int debug = 0) {
3995: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
3996: const int nVertices = 14;
3997: const int nFaces = 12;
3998: double ilower[3];
3999: ilower[0] = lower[0]*(1. - offset[0]) + upper[0]*offset[0];
4000: ilower[1] = lower[1]*(1. - offset[1]) + upper[1]*offset[1];
4001: ilower[2] = lower[2]*(1. - offset[2]) + upper[2]*offset[2];
4002: double coords[nVertices*3];
4003: //outer square-triplet
4004: coords[0*3+0] = lower[0];
4005: coords[0*3+1] = lower[1];
4006: coords[0*3+2] = upper[2];
4007: coords[1*3+0] = upper[0];
4008: coords[1*3+1] = lower[1];
4009: coords[1*3+2] = lower[2];
4010: coords[2*3+0] = lower[0];
4011: coords[2*3+1] = upper[1];
4012: coords[2*3+2] = lower[2];
4013: coords[3*3+0] = upper[0];
4014: coords[3*3+1] = upper[1];
4015: coords[3*3+2] = lower[2];
4016: coords[4*3+0] = lower[0];
4017: coords[4*3+1] = lower[1];
4018: coords[4*3+2] = lower[2];
4019: coords[5*3+0] = lower[0];
4020: coords[5*3+1] = upper[1];
4021: coords[5*3+2] = upper[2];
4022: coords[6*3+0] = upper[0];
4023: coords[6*3+1] = lower[1];
4024: coords[6*3+2] = upper[2];
4026: //inner square-triplet
4027: coords[7*3+0] = ilower[0];
4028: coords[7*3+1] = upper[1];
4029: coords[7*3+2] = upper[2];
4030: coords[8*3+0] = upper[0];
4031: coords[8*3+1] = ilower[1];
4032: coords[8*3+2] = upper[2];
4033: coords[9*3+0] = upper[0];
4034: coords[9*3+1] = ilower[1];
4035: coords[9*3+2] = ilower[2];
4036: coords[10*3+0] = ilower[0];
4037: coords[10*3+1] = upper[1];
4038: coords[10*3+2] = ilower[2];
4039: coords[11*3+0] = upper[0];
4040: coords[11*3+1] = upper[1];
4041: coords[11*3+2] = ilower[2];
4042: coords[12*3+0] = ilower[0];
4043: coords[12*3+1] = ilower[1];
4044: coords[12*3+2] = upper[2];
4045: coords[13*3+0] = ilower[0];
4046: coords[13*3+1] = ilower[1];
4047: coords[13*3+2] = ilower[2];
4049:
4050: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
4051: mesh->setSieve(sieve);
4052: typename Mesh::point_type p[nVertices];
4053: typename Mesh::point_type f[nFaces];
4054: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
4055: for (int i = 0; i < nVertices; i++) {
4056: p[i] = typename Mesh::point_type(i+nFaces);
4057: mesh->setValue(markers, p[i], 1);
4058: }
4059: for (int i = 0; i < nFaces; i++) {
4060: f[i] = typename Mesh::point_type(i);
4061: }
4062: int order = 0;
4063: //assemble the larger square sides
4064: sieve->addArrow(p[0], f[0], order++);
4065: sieve->addArrow(p[5], f[0], order++);
4066: sieve->addArrow(p[2], f[0], order++);
4067: sieve->addArrow(p[4], f[0], order++);
4068: mesh->setValue(markers, f[0], 1);
4070: sieve->addArrow(p[0], f[1], order++);
4071: sieve->addArrow(p[4], f[1], order++);
4072: sieve->addArrow(p[1], f[1], order++);
4073: sieve->addArrow(p[6], f[1], order++);
4074: mesh->setValue(markers, f[1], 1);
4076: sieve->addArrow(p[4], f[2], order++);
4077: sieve->addArrow(p[1], f[2], order++);
4078: sieve->addArrow(p[3], f[2], order++);
4079: sieve->addArrow(p[2], f[2], order++);
4080: mesh->setValue(markers, f[2], 1);
4081:
4082: //assemble the L-shaped sides
4084: sieve->addArrow(p[0], f[3], order++);
4085: sieve->addArrow(p[12], f[3], order++);
4086: sieve->addArrow(p[7], f[3], order++);
4087: sieve->addArrow(p[5], f[3], order++);
4088: mesh->setValue(markers, f[3], 1);
4090: sieve->addArrow(p[0], f[4], order++);
4091: sieve->addArrow(p[12],f[4], order++);
4092: sieve->addArrow(p[8], f[4], order++);
4093: sieve->addArrow(p[6], f[4], order++);
4094: mesh->setValue(markers, f[4], 1);
4096: sieve->addArrow(p[9], f[5], order++);
4097: sieve->addArrow(p[1], f[5], order++);
4098: sieve->addArrow(p[3], f[5], order++);
4099: sieve->addArrow(p[11], f[5], order++);
4100: mesh->setValue(markers, f[5], 1);
4102: sieve->addArrow(p[9], f[6], order++);
4103: sieve->addArrow(p[1], f[6], order++);
4104: sieve->addArrow(p[6], f[6], order++);
4105: sieve->addArrow(p[8], f[6], order++);
4106: mesh->setValue(markers, f[6], 1);
4108: sieve->addArrow(p[10], f[7], order++);
4109: sieve->addArrow(p[2], f[7], order++);
4110: sieve->addArrow(p[5], f[7], order++);
4111: sieve->addArrow(p[7], f[7], order++);
4112: mesh->setValue(markers, f[7], 1);
4114: sieve->addArrow(p[10], f[8], order++);
4115: sieve->addArrow(p[2], f[8], order++);
4116: sieve->addArrow(p[3], f[8], order++);
4117: sieve->addArrow(p[11], f[8], order++);
4118: mesh->setValue(markers, f[8], 1);
4120: //assemble the smaller square sides
4122: sieve->addArrow(p[13], f[9], order++);
4123: sieve->addArrow(p[10], f[9], order++);
4124: sieve->addArrow(p[11], f[9], order++);
4125: sieve->addArrow(p[9], f[9], order++);
4126: mesh->setValue(markers, f[9], 1);
4128: sieve->addArrow(p[12], f[10], order++);
4129: sieve->addArrow(p[7], f[10], order++);
4130: sieve->addArrow(p[10], f[10], order++);
4131: sieve->addArrow(p[13], f[10], order++);
4132: mesh->setValue(markers, f[10], 1);
4134: sieve->addArrow(p[8], f[11], order++);
4135: sieve->addArrow(p[12], f[11], order++);
4136: sieve->addArrow(p[13], f[11], order++);
4137: sieve->addArrow(p[9], f[11], order++);
4138: mesh->setValue(markers, f[11], 1);
4140: mesh->stratify();
4141: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
4142: Obj<typename Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
4143: //coordinates->view("coordinates");
4144: return mesh;
4146: }
4150: /*
4151: //"sphere" out a cube
4153: */
4154: #if 0
4155: static Obj<Mesh> createSphereBoundary(const MPI_Comm comm, const double radius, const int refinement, const int debug = 0) {
4156: Obj<Mesh> m = new Mesh(comm, 2, debug);
4157: Obj<Mesh::sieve_type> s = new Mesh::sieve_type(comm, debug);
4158: m->setSieve(s);
4159: Mesh::point_type p = 0;
4160: int nVertices = 8+12*(refinement)+6*(refinement)*(refinement);
4161: Mesh::point_type vertices[nVertices];
4162: double coords[3*nVertices];
4163: int nCells = 6*2*(refinement+1)*(refinement+1);
4164: double delta = 2./((double)(refinement+1));
4165: Mesh::point_type cells[nCells];
4166: for (int i = 0; i < nCells; i++) {
4167: cells[i] = p;
4168: p++;
4169: }
4170: for (int i = 0; i < nVertices; i++) {
4171: vertices[i] = p;
4172: p++;
4173: }
4174: //set up the corners;
4175: //lll
4176: coords[0*3+0] = -1.;
4177: coords[0*3+1] = -1.;
4178: coords[0*3+2] = -1.;
4179: //llh
4180: coords[1*3+0] = -1.;
4181: coords[1*3+1] = -1.;
4182: coords[1*3+2] = 1.;
4183: //lhh
4184: coords[2*3+0] = -1.;
4185: coords[2*3+1] = 1.;
4186: coords[2*3+2] = 1.;
4187: //lhl
4188: coords[3*3+0] = -1.;
4189: coords[3*3+1] = 1.;
4190: coords[3*3+2] = -1.;
4191: //hhl
4192: coords[4*3+0] = 1.;
4193: coords[4*3+1] = 1.;
4194: coords[4*3+2] = -1.;
4195: //hhh
4196: coords[5*3+0] = 1.;
4197: coords[5*3+1] = 1.;
4198: coords[5*3+2] = 1.;
4199: //hlh
4200: coords[6*3+0] = 1.;
4201: coords[6*3+1] = -1.;
4202: coords[6*3+2] = 1.;
4203: //hll
4204: coords[7*3+0] = 1.;
4205: coords[7*3+1] = -1.;
4206: coords[7*3+2] = -1.;
4207: //set up the edges (always go low to high)
4208: //xll
4209: for (int i = 0; i < refinement; i++) {
4210: coords[3*(8+0*refinement+i)+0] = -1. + delta*i;
4211: coords[3*(8+0*refinement+i)+1] = -1.;
4212: coords[3*(8+0*refinement+i)+2] = -1.;
4213: }
4214: //xlh
4215: for (int i = 0; i < refinement; i++) {
4216: coords[3*(8+1*refinement+i)+0] = -1. + delta*i;
4217: coords[3*(8+1*refinement+i)+1] = -1.;
4218: coords[3*(8+1*refinement+i)+2] = 1.;
4219: }
4220: //xhh
4221: for (int i = 0; i < refinement; i++) {
4222: coords[3*(8+2*refinement+i)+0] = -1. + delta*i;
4223: coords[3*(8+2*refinement+i)+1] = 1.;
4224: coords[3*(8+2*refinement+i)+2] = 1.;
4225: }
4226: //xhl
4227: for (int i = 0; i < refinement; i++) {
4228: coords[3*(8+3*refinement+i)+0] = -1. + delta*i;
4229: coords[3*(8+3*refinement+i)+1] = 1.;
4230: coords[3*(8+3*refinement+i)+2] = -1.;
4231: }
4232: //lxl
4233: for (int i = 0; i < refinement; i++) {
4234: coords[3*(8+4*refinement+i)+0] = -1.;
4235: coords[3*(8+4*refinement+i)+1] = -1. + delta*i;
4236: coords[3*(8+4*refinement+i)+2] = -1.;
4237: }
4238: //lxh
4239: for (int i = 0; i < refinement; i++) {
4240: coords[3*(8+5*refinement+i)+0] = -1.;
4241: coords[3*(8+5*refinement+i)+1] = -1. + delta*i;
4242: coords[3*(8+5*refinement+i)+2] = 1.;
4243: }
4244: //hxh
4245: for (int i = 0; i < refinement; i++) {
4246: coords[3*(8+6*refinement+i)+0] = 1.;
4247: coords[3*(8+6*refinement+i)+1] = -1. + delta*i;
4248: coords[3*(8+6*refinement+i)+2] = 1.;
4249: }
4250: //hxl
4251: for (int i = 0; i < refinement; i++) {
4252: coords[3*(8+7*refinement+i)+0] = 1.;
4253: coords[3*(8+7*refinement+i)+1] = -1. + delta*i;
4254: coords[3*(8+7*refinement+i)+2] = -1.;
4255: }
4256: //llx
4257: for (int i = 0; i < refinement; i++) {
4258: coords[3*(8+8*refinement+i)+0] = -1.;
4259: coords[3*(8+8*refinement+i)+1] = -1.;
4260: coords[3*(8+8*refinement+i)+2] = -1. + delta*i;
4261: }
4262: //lhx
4263: for (int i = 0; i < refinement; i++) {
4264: coords[3*(8+9*refinement+i)+0] = -1.;
4265: coords[3*(8+9*refinement+i)+1] = 1.;
4266: coords[3*(8+9*refinement+i)+2] = -1. + delta*i;
4267: }
4268: //hhx
4269: for (int i = 0; i < refinement; i++) {
4270: coords[3*(8+10*refinement+i)+0] = 1.;
4271: coords[3*(8+10*refinement+i)+1] = 1.;
4272: coords[3*(8+10*refinement+i)+2] = -1. + delta*i;
4273: }
4274: //hlx
4275: for (int i = 0; i < refinement; i++) {
4276: coords[3*(8+11*refinement+i)+0] = 1.;
4277: coords[3*(8+11*refinement+i)+1] = -1.;
4278: coords[3*(8+11*refinement+i)+2] = -1. + delta*i;
4279: }
4280: //set up the faces
4281: //lxx
4282: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4283: coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+0] = -1.;
4284: coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4285: coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4286: }
4287: //hxx
4288: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4289: coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+0] = 1.;
4290: coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4291: coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4292: }
4293: //xlx
4294: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4295: coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+0] = -1. + delta*j;
4296: coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+1] = -1.;
4297: coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4298: }
4299: //xhx
4300: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4301: coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+0] = -1. + delta*j;
4302: coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+1] = 1.;
4303: coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4304: }
4305: //xxl
4306: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4307: coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+0] = -1.;
4308: coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4309: coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4310: }
4311: //xxh
4312: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4313: coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+0] = 1.;
4314: coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4315: coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4316: }
4317: //stitch the corners up with the edges and the faces
4318:
4319: //stitch the edges to the faces
4320: //fill in the faces
4321: int face_offset = 8 + 12*refinement;
4322: for (int i = 0; i < 6; i++) for (int j = 0; j < refinement; j++) for (int k = 0; k < refinement; k++) {
4323: //build each square doublet
4324: }
4325: }
4327: #endif
4331: /*
4332: Simple cubic boundary:
4334: 30----31-----32
4335: | | |
4336: | 3 | 2 |
4337: | | |
4338: 27----28-----29
4339: | | |
4340: | 0 | 1 |
4341: | | |
4342: 24----25-----26
4343: */
4344: static Obj<Mesh> createParticleInCubeBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int faces[], const double radius, const int thetaEdges, const int phiSlices, const int debug = 0) {
4345: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
4346: const int numCubeVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
4347: const int numPartVertices = (thetaEdges - 1)*phiSlices + 2;
4348: const int numVertices = numCubeVertices + numPartVertices;
4349: const int numCubeFaces = 6;
4350: const int numFaces = numCubeFaces + thetaEdges*phiSlices;
4351: double *coords = new double[numVertices*3];
4352: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
4353: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
4354: int order = 0;
4356: mesh->setSieve(sieve);
4357: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
4358: if (mesh->commRank() == 0) {
4359: // Make cube
4360: for(int v = numFaces; v < numFaces+numVertices; v++) {
4361: vertices[v-numFaces] = typename Mesh::point_type(v);
4362: mesh->setValue(markers, vertices[v-numFaces], 1);
4363: }
4364: {
4365: // Side 0 (Front)
4366: typename Mesh::point_type face(0);
4367: sieve->addArrow(vertices[0], face, order++);
4368: sieve->addArrow(vertices[1], face, order++);
4369: sieve->addArrow(vertices[2], face, order++);
4370: sieve->addArrow(vertices[3], face, order++);
4371: mesh->setValue(markers, face, 1);
4372: }
4373: {
4374: // Side 1 (Back)
4375: typename Mesh::point_type face(1);
4376: sieve->addArrow(vertices[5], face, order++);
4377: sieve->addArrow(vertices[4], face, order++);
4378: sieve->addArrow(vertices[7], face, order++);
4379: sieve->addArrow(vertices[6], face, order++);
4380: mesh->setValue(markers, face, 1);
4381: }
4382: {
4383: // Side 2 (Bottom)
4384: typename Mesh::point_type face(2);
4385: sieve->addArrow(vertices[4], face, order++);
4386: sieve->addArrow(vertices[5], face, order++);
4387: sieve->addArrow(vertices[1], face, order++);
4388: sieve->addArrow(vertices[0], face, order++);
4389: mesh->setValue(markers, face, 1);
4390: }
4391: {
4392: // Side 3 (Top)
4393: typename Mesh::point_type face(3);
4394: sieve->addArrow(vertices[3], face, order++);
4395: sieve->addArrow(vertices[2], face, order++);
4396: sieve->addArrow(vertices[6], face, order++);
4397: sieve->addArrow(vertices[7], face, order++);
4398: mesh->setValue(markers, face, 1);
4399: }
4400: {
4401: // Side 4 (Left)
4402: typename Mesh::point_type face(4);
4403: sieve->addArrow(vertices[4], face, order++);
4404: sieve->addArrow(vertices[0], face, order++);
4405: sieve->addArrow(vertices[3], face, order++);
4406: sieve->addArrow(vertices[7], face, order++);
4407: mesh->setValue(markers, face, 1);
4408: }
4409: {
4410: // Side 5 (Right)
4411: typename Mesh::point_type face(5);
4412: sieve->addArrow(vertices[1], face, order++);
4413: sieve->addArrow(vertices[5], face, order++);
4414: sieve->addArrow(vertices[6], face, order++);
4415: sieve->addArrow(vertices[2], face, order++);
4416: mesh->setValue(markers, face, 1);
4417: }
4418: // Make particle
4419: for(int s = 0; s < phiSlices; ++s) {
4420: for(int ep = 0; ep < thetaEdges; ++ep) {
4421: // Vertices on each slice are 0..thetaEdges
4422: typename Mesh::point_type face(numCubeFaces + s*thetaEdges + ep);
4423: int vertexA = numCubeVertices + ep + 0 + s*(thetaEdges+1);
4424: int vertexB = numCubeVertices + ep + 1 + s*(thetaEdges+1);
4425: int vertexC = numCubeVertices + (ep + 1 + (s+1)*(thetaEdges+1))%((thetaEdges+1)*phiSlices);
4426: int vertexD = numCubeVertices + (ep + 0 + (s+1)*(thetaEdges+1))%((thetaEdges+1)*phiSlices);
4427: const int correction1 = (s > 0)*((s-1)*2 + 1);
4428: const int correction2 = (s < phiSlices-1)*(s*2 + 1);
4430: if ((vertexA - numCubeVertices)%(thetaEdges+1) == 0) {
4431: vertexA = vertexD = numCubeVertices;
4432: vertexB -= correction1;
4433: vertexC -= correction2;
4434: } else if ((vertexB - numCubeVertices)%(thetaEdges+1) == thetaEdges) {
4435: vertexA -= correction1;
4436: vertexD -= correction2;
4437: vertexB = vertexC = numCubeVertices + thetaEdges;
4438: } else {
4439: vertexA -= correction1;
4440: vertexB -= correction1;
4441: vertexC -= correction2;
4442: vertexD -= correction2;
4443: }
4444: if ((vertexA >= numVertices) || (vertexB >= numVertices) || (vertexC >= numVertices) || (vertexD >= numVertices)) {
4445: throw ALE::Exception("Bad vertex");
4446: }
4447: sieve->addArrow(vertices[vertexA], face, order++);
4448: sieve->addArrow(vertices[vertexB], face, order++);
4449: if (vertexB != vertexC) sieve->addArrow(vertices[vertexC], face, order++);
4450: if (vertexA != vertexD) sieve->addArrow(vertices[vertexD], face, order++);
4451: mesh->setValue(markers, face, 2);
4452: mesh->setValue(markers, vertices[vertexA], 2);
4453: mesh->setValue(markers, vertices[vertexB], 2);
4454: if (vertexB != vertexC) mesh->setValue(markers, vertices[vertexC], 2);
4455: if (vertexA != vertexD) mesh->setValue(markers, vertices[vertexD], 2);
4456: }
4457: }
4458: }
4459: mesh->stratify();
4460: #if 0
4461: for(int vz = 0; vz <= edges[2]; ++vz) {
4462: for(int vy = 0; vy <= edges[1]; ++vy) {
4463: for(int vx = 0; vx <= edges[0]; ++vx) {
4464: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
4465: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
4466: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
4467: }
4468: }
4469: }
4470: #else
4471: coords[0*3+0] = lower[0];
4472: coords[0*3+1] = lower[1];
4473: coords[0*3+2] = upper[2];
4474: coords[1*3+0] = upper[0];
4475: coords[1*3+1] = lower[1];
4476: coords[1*3+2] = upper[2];
4477: coords[2*3+0] = upper[0];
4478: coords[2*3+1] = upper[1];
4479: coords[2*3+2] = upper[2];
4480: coords[3*3+0] = lower[0];
4481: coords[3*3+1] = upper[1];
4482: coords[3*3+2] = upper[2];
4483: coords[4*3+0] = lower[0];
4484: coords[4*3+1] = lower[1];
4485: coords[4*3+2] = lower[2];
4486: coords[5*3+0] = upper[0];
4487: coords[5*3+1] = lower[1];
4488: coords[5*3+2] = lower[2];
4489: coords[6*3+0] = upper[0];
4490: coords[6*3+1] = upper[1];
4491: coords[6*3+2] = lower[2];
4492: coords[7*3+0] = lower[0];
4493: coords[7*3+1] = upper[1];
4494: coords[7*3+2] = lower[2];
4495: #endif
4496: const double centroidX = 0.5*(upper[0] + lower[0]);
4497: const double centroidY = 0.5*(upper[1] + lower[1]);
4498: const double centroidZ = 0.5*(upper[2] + lower[2]);
4499: for(int s = 0; s < phiSlices; ++s) {
4500: for(int v = 0; v <= thetaEdges; ++v) {
4501: int vertex = numCubeVertices + v + s*(thetaEdges+1);
4502: const double theta = v*(PETSC_PI/thetaEdges);
4503: const double phi = s*(2.0*PETSC_PI/phiSlices);
4504: const int correction = (s > 0)*((s-1)*2 + 1);
4506: if ((vertex- numCubeVertices)%(thetaEdges+1) == 0) {
4507: vertex = numCubeVertices;
4508: } else if ((vertex - numCubeVertices)%(thetaEdges+1) == thetaEdges) {
4509: vertex = numCubeVertices + thetaEdges;
4510: } else {
4511: vertex -= correction;
4512: }
4513: coords[vertex*3+0] = centroidX + radius*sin(theta)*cos(phi);
4514: coords[vertex*3+1] = centroidY + radius*sin(theta)*sin(phi);
4515: coords[vertex*3+2] = centroidZ + radius*cos(theta);
4516: }
4517: }
4518: delete [] vertices;
4519: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
4520: return mesh;
4521: };
4522: // This method takes a tetrahedral mesh and performs a 1 --> 8 refinement of each cell
4523: // It does this by adding a new vertex at the midpoint of each edge
4524: template<typename MeshType, typename MapType>
4525: static void refineTetrahedra(MeshType& mesh, MeshType& newMesh, MapType& edge2vertex) {
4526: typedef typename MeshType::sieve_type sieve_type;
4527: typedef typename MeshType::point_type point_type;
4528: typedef typename MapType::key_type edge_type;
4530: const int numCells = mesh.heightStratum(0)->size();
4531: const int numVertices = mesh.depthStratum(0)->size();
4532: // Calculate number of new cells
4533: const int numNewCells = numCells * 8;
4534: // Bound number of new vertices
4535: //const int maxNewVertices = numCells * 6;
4536: int curNewVertex = numNewCells + numVertices;
4538: // Loop over cells
4539: const Obj<sieve_type>& sieve = mesh.getSieve();
4540: const Obj<sieve_type>& newSieve = newMesh.getSieve();
4541: ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
4543: // First compute map from edges to new vertices
4544: for(int c = 0; c < numCells; ++c) {
4545: sieve->cone(c, cV);
4546: assert(cV.getSize() == 4);
4547: const point_type *cone = cV.getPoints();
4549: // As per Brad's diagram
4550: edge_type edges[6] = {edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1])),
4551: edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2])),
4552: edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0])),
4553: edge_type(std::min(cone[0], cone[3]), std::max(cone[0], cone[3])),
4554: edge_type(std::min(cone[1], cone[3]), std::max(cone[1], cone[3])),
4555: edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]))};
4556: // Check that vertex does not yet exist
4557: for(int v = 0; v < 6; ++v) {
4558: if (edge2vertex.find(edges[v]) == edge2vertex.end()) {
4559: edge2vertex[edges[v]] = curNewVertex++;
4560: }
4561: }
4562: cV.clear();
4563: }
4564: // Reallocate the sieve chart
4565: newSieve->setChart(typename sieve_type::chart_type(0, curNewVertex));
4566: // Create new sieve with correct sizes for refined cells
4567: for(int c = 0; c < numCells; ++c) {
4568: sieve->cone(c, cV);
4569: assert(cV.getSize() == 4);
4570: const point_type *cone = cV.getPoints();
4572: // As per Brad's diagram
4573: edge_type edges[6] = {edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1])),
4574: edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2])),
4575: edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0])),
4576: edge_type(std::min(cone[0], cone[3]), std::max(cone[0], cone[3])),
4577: edge_type(std::min(cone[1], cone[3]), std::max(cone[1], cone[3])),
4578: edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]))};
4579: // Check that vertex does not yet exist
4580: point_type newVertices[6];
4582: for(int v = 0; v < 6; ++v) {
4583: newVertices[v] = edge2vertex[edges[v]];
4584: }
4585: // Set new sizes
4586: for(int nc = 0; nc < 8; ++nc) {newSieve->setConeSize(c*8+nc, 4);}
4587: const point_type offset = numNewCells - numCells;
4589: point_type cell0[4] = {cone[0]+offset, newVertices[3], newVertices[0], newVertices[2]};
4590: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell0[v], 1);}
4591: point_type cell1[4] = {cone[1]+offset, newVertices[4], newVertices[1], newVertices[0]};
4592: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell1[v], 1);}
4593: point_type cell2[4] = {cone[2]+offset, newVertices[5], newVertices[2], newVertices[1]};
4594: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell2[v], 1);}
4595: point_type cell3[4] = {cone[3]+offset, newVertices[3], newVertices[5], newVertices[4]};
4596: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell3[v], 1);}
4597: point_type cell4[4] = {newVertices[0], newVertices[3], newVertices[4], newVertices[2]};
4598: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell4[v], 1);}
4599: point_type cell5[4] = {newVertices[1], newVertices[4], newVertices[5], newVertices[3]};
4600: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell5[v], 1);}
4601: point_type cell6[4] = {newVertices[2], newVertices[5], newVertices[3], newVertices[1]};
4602: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell6[v], 1);}
4603: point_type cell7[4] = {newVertices[0], newVertices[1], newVertices[2], newVertices[3]};
4604: for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell7[v], 1);}
4605: cV.clear();
4606: }
4607: newSieve->allocate();
4608: const int numNewVertices = curNewVertex - numNewCells;
4609: point_type *vertex2edge = new point_type[numNewVertices*2];
4611: // Create refined cells in new sieve
4612: for(int c = 0; c < numCells; ++c) {
4613: sieve->cone(c, cV);
4614: assert(cV.getSize() == 4);
4615: const point_type *cone = cV.getPoints();
4617: // As per Brad's diagram
4618: edge_type edges[6] = {edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1])),
4619: edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2])),
4620: edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0])),
4621: edge_type(std::min(cone[0], cone[3]), std::max(cone[0], cone[3])),
4622: edge_type(std::min(cone[1], cone[3]), std::max(cone[1], cone[3])),
4623: edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]))};
4624: // Check that vertex does not yet exist
4625: point_type newVertices[6];
4627: for(int v = 0; v < 6; ++v) {
4628: if (edge2vertex.find(edges[v]) == edge2vertex.end()) {
4629: throw ALE::Exception("Missing edge in refined mesh");
4630: }
4631: newVertices[v] = edge2vertex[edges[v]];
4632: vertex2edge[(newVertices[v]-numNewCells-numVertices)*2+0] = edges[v].first;
4633: vertex2edge[(newVertices[v]-numNewCells-numVertices)*2+1] = edges[v].second;
4634: }
4635: // Create new cells
4636: const point_type offset = numNewCells - numCells;
4638: point_type cell0[4] = {cone[0]+offset, newVertices[3], newVertices[0], newVertices[2]};
4639: newSieve->setCone(cell0, c*8+0);
4640: point_type cell1[4] = {cone[1]+offset, newVertices[4], newVertices[1], newVertices[0]};
4641: newSieve->setCone(cell1, c*8+1);
4642: point_type cell2[4] = {cone[2]+offset, newVertices[5], newVertices[2], newVertices[1]};
4643: newSieve->setCone(cell2, c*8+2);
4644: point_type cell3[4] = {cone[3]+offset, newVertices[3], newVertices[5], newVertices[4]};
4645: newSieve->setCone(cell3, c*8+3);
4646: point_type cell4[4] = {newVertices[0], newVertices[3], newVertices[4], newVertices[2]};
4647: newSieve->setCone(cell4, c*8+4);
4648: point_type cell5[4] = {newVertices[1], newVertices[4], newVertices[5], newVertices[3]};
4649: newSieve->setCone(cell5, c*8+5);
4650: point_type cell6[4] = {newVertices[2], newVertices[5], newVertices[3], newVertices[1]};
4651: newSieve->setCone(cell6, c*8+6);
4652: point_type cell7[4] = {newVertices[0], newVertices[1], newVertices[2], newVertices[3]};
4653: newSieve->setCone(cell7, c*8+7);
4654: cV.clear();
4655: }
4656: newSieve->symmetrize();
4657: // Create new coordinates
4658: const Obj<typename MeshType::real_section_type>& coordinates = mesh.getRealSection("coordinates");
4659: const Obj<typename MeshType::real_section_type>& newCoordinates = newMesh.getRealSection("coordinates");
4661: newCoordinates->setChart(typename sieve_type::chart_type(numNewCells, curNewVertex));
4662: for(int v = numNewCells; v < curNewVertex; ++v) {
4663: newCoordinates->setFiberDimension(v, 3);
4664: }
4665: newCoordinates->allocatePoint();
4666: for(int v = 0; v < numVertices; ++v) {
4667: newCoordinates->updatePoint(v+numNewCells, coordinates->restrictPoint(v+numCells));
4668: }
4669: for(int v = numNewCells+numVertices; v < curNewVertex; ++v) {
4670: const int endpointA = vertex2edge[(v-numNewCells-numVertices)*2+0];
4671: const int endpointB = vertex2edge[(v-numNewCells-numVertices)*2+1];
4672: const double *coordsA = coordinates->restrictPoint(endpointA);
4673: double coords[3];
4675: for(int d = 0; d < 3; ++d) {
4676: coords[d] = coordsA[d];
4677: }
4678: const double *coordsB = coordinates->restrictPoint(endpointB);
4679: for(int d = 0; d < 3; ++d) {
4680: coords[d] += coordsB[d];
4681: coords[d] *= 0.5;
4682: }
4683: newCoordinates->updatePoint(v, coords);
4684: }
4685: delete [] vertex2edge;
4686: // Fast stratification
4687: const Obj<typename MeshType::label_type>& height = newMesh.createLabel("height");
4688: const Obj<typename MeshType::label_type>& depth = newMesh.createLabel("depth");
4689: for(int c = 0; c < numNewCells; ++c) {
4690: height->setCone(0, c);
4691: depth->setCone(1, c);
4692: }
4693: for(int v = numNewCells; v < numNewCells+numNewVertices; ++v) {
4694: height->setCone(1, v);
4695: depth->setCone(0, v);
4696: }
4697: newMesh.setHeight(1);
4698: newMesh.setDepth(1);
4699: // Exchange new boundary vertices
4700: // We can convert endpoints, and then just match to new vertex on this side
4701: // 1) Create the overlap of edges which are vertex pairs (do not need for interpolated meshes)
4702: // 2) Create a section of overlap edge --> new vertex (this will generalize to other split points in interpolated meshes)
4703: // 3) Copy across new overlap
4704: // 4) Fuse matches new vertex pairs and inserts them into the old overlap
4707: // Create the parallel overlap
4708: int *numCellsP = new int[mesh.commSize()];
4709: int *numNewCellsP = new int[newMesh.commSize()];
4710: int ierr;
4712: MPI_Allgather((void *) &numCells, 1, MPI_INT, numCellsP, 1, MPI_INT, mesh.comm());CHKERRXX(ierr);
4713: MPI_Allgather((void *) &numNewCells, 1, MPI_INT, numNewCellsP, 1, MPI_INT, newMesh.comm());CHKERRXX(ierr);
4714: Obj<typename MeshType::send_overlap_type> newSendOverlap = newMesh.getSendOverlap();
4715: Obj<typename MeshType::recv_overlap_type> newRecvOverlap = newMesh.getRecvOverlap();
4716: const Obj<typename MeshType::send_overlap_type>& sendOverlap = mesh.getSendOverlap();
4717: const Obj<typename MeshType::recv_overlap_type>& recvOverlap = mesh.getRecvOverlap();
4718: Obj<typename MeshType::send_overlap_type::traits::capSequence> sendPoints = sendOverlap->cap();
4719: const typename MeshType::send_overlap_type::source_type localOffset = numNewCellsP[newMesh.commRank()] - numCellsP[mesh.commRank()];
4721: for(typename MeshType::send_overlap_type::traits::capSequence::iterator p_iter = sendPoints->begin(); p_iter != sendPoints->end(); ++p_iter) {
4722: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& ranks = sendOverlap->support(*p_iter);
4723: const typename MeshType::send_overlap_type::source_type& localPoint = *p_iter;
4725: for(typename MeshType::send_overlap_type::traits::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
4726: const int rank = *r_iter;
4727: const typename MeshType::send_overlap_type::source_type& remotePoint = r_iter.color();
4728: const typename MeshType::send_overlap_type::source_type remoteOffset = numNewCellsP[rank] - numCellsP[rank];
4730: newSendOverlap->addArrow(localPoint+localOffset, rank, remotePoint+remoteOffset);
4731: }
4732: }
4733: Obj<typename MeshType::recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();
4735: for(typename MeshType::recv_overlap_type::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
4736: const Obj<typename MeshType::recv_overlap_type::traits::coneSequence>& ranks = recvOverlap->cone(*p_iter);
4737: const typename MeshType::recv_overlap_type::target_type& localPoint = *p_iter;
4739: for(typename MeshType::recv_overlap_type::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
4740: const int rank = *r_iter;
4741: const typename MeshType::recv_overlap_type::target_type& remotePoint = r_iter.color();
4742: const typename MeshType::recv_overlap_type::target_type remoteOffset = numNewCellsP[rank] - numCellsP[rank];
4744: newRecvOverlap->addArrow(rank, localPoint+localOffset, remotePoint+remoteOffset);
4745: }
4746: }
4747: newMesh.setCalculatedOverlap(true);
4748: delete [] numCellsP;
4749: delete [] numNewCellsP;
4750: // Check edges in edge2vertex for both endpoints sent to same process
4751: // Put it in section with point being the lowest numbered vertex and value (other endpoint, new vertex)
4752: Obj<ALE::Section<point_type, edge_type> > newVertices = new ALE::Section<point_type, edge_type>(mesh.comm());
4753: std::map<edge_type, std::vector<int> > bdedge2rank;
4755: for(typename std::map<edge_type, point_type>::const_iterator e_iter = edge2vertex.begin(); e_iter != edge2vertex.end(); ++e_iter) {
4756: const point_type left = e_iter->first.first;
4757: const point_type right = e_iter->first.second;
4759: if (sendOverlap->capContains(left) && sendOverlap->capContains(right)) {
4760: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& leftRanksSeq = sendOverlap->support(left);
4761: std::list<int> leftRanks(leftRanksSeq->begin(), leftRanksSeq->end());
4762: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& rightRanks = sendOverlap->support(right);
4763: std::list<int> ranks;
4764: std::set_intersection(leftRanks.begin(), leftRanks.end(), rightRanks->begin(), rightRanks->end(),
4765: std::insert_iterator<std::list<int> >(ranks, ranks.begin()));
4767: if(ranks.size()) {
4768: newVertices->addFiberDimension(std::min(e_iter->first.first, e_iter->first.second)+localOffset, 1);
4769: for(typename std::list<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
4770: bdedge2rank[e_iter->first].push_back(*r_iter);
4771: }
4772: }
4773: }
4774: }
4775: newVertices->allocatePoint();
4776: const typename ALE::Section<point_type, edge_type>::chart_type& chart = newVertices->getChart();
4778: for(typename ALE::Section<point_type, edge_type>::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
4779: typedef typename ALE::Section<point_type, edge_type>::value_type value_type;
4780: const point_type p = *c_iter;
4781: const int dim = newVertices->getFiberDimension(p);
4782: int v = 0;
4783: value_type *values = new value_type[dim];
4785: for(typename std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end() && v < dim; ++e_iter) {
4786: if (std::min(e_iter->first.first, e_iter->first.second)+localOffset == p) {
4787: values[v++] = edge_type(std::max(e_iter->first.first, e_iter->first.second)+localOffset, edge2vertex[e_iter->first]);
4788: }
4789: }
4790: newVertices->updatePoint(p, values);
4791: delete [] values;
4792: }
4793: // Copy across overlap
4794: typedef ALE::Pair<int, point_type> overlap_point_type;
4795: Obj<ALE::Section<overlap_point_type, edge_type> > overlapVertices = new ALE::Section<overlap_point_type, edge_type>(mesh.comm());
4797: ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVertices, overlapVertices);
4798: // Merge by translating edge to local points, finding edge in edge2vertex, and adding (local new vetex, remote new vertex) to overlap
4799: for(typename std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end(); ++e_iter) {
4800: const point_type localPoint = edge2vertex[e_iter->first];
4802: for(typename std::vector<int>::const_iterator r_iter = e_iter->second.begin(); r_iter != e_iter->second.end(); ++r_iter) {
4803: point_type remoteLeft = -1, remoteRight = -1;
4804: const int rank = *r_iter;
4806: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& leftRanks = newSendOverlap->support(e_iter->first.first+localOffset);
4807: for(typename MeshType::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks->begin(); lr_iter != leftRanks->end(); ++lr_iter) {
4808: if (rank == *lr_iter) {
4809: remoteLeft = lr_iter.color();
4810: break;
4811: }
4812: }
4813: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& rightRanks = newSendOverlap->support(e_iter->first.second+localOffset);
4814: for(typename MeshType::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks->begin(); rr_iter != rightRanks->end(); ++rr_iter) {
4815: if (rank == *rr_iter) {
4816: remoteRight = rr_iter.color();
4817: break;
4818: }
4819: }
4820: const point_type remoteMin = std::min(remoteLeft, remoteRight);
4821: const point_type remoteMax = std::max(remoteLeft, remoteRight);
4822: const int remoteSize = overlapVertices->getFiberDimension(overlap_point_type(rank, remoteMin));
4823: const edge_type *remoteVals = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
4824: point_type remotePoint = -1;
4826: for(int d = 0; d < remoteSize; ++d) {
4827: if (remoteVals[d].first == remoteMax) {
4828: remotePoint = remoteVals[d].second;
4829: break;
4830: }
4831: }
4832: newSendOverlap->addArrow(localPoint, rank, remotePoint);
4833: newRecvOverlap->addArrow(rank, localPoint, remotePoint);
4834: }
4835: }
4836: };
4837: };
4838: class MeshSerializer {
4839: public:
4840: template<typename Mesh>
4841: static void writeMesh(const std::string& filename, Mesh& mesh) {
4842: std::ofstream fs;
4844: if (mesh.commRank() == 0) {
4845: fs.open(filename.c_str());
4846: }
4847: writeMesh(fs, mesh);
4848: if (mesh.commRank() == 0) {
4849: fs.close();
4850: }
4851: };
4852: template<typename Mesh>
4853: static void writeMesh(std::ofstream& fs, Mesh& mesh) {
4854: ISieveSerializer::writeSieve(fs, *mesh.getSieve());
4855: // Write labels
4856: const typename Mesh::labels_type& labels = mesh.getLabels();
4858: if (!mesh.commRank()) {fs << labels.size() << std::endl;}
4859: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
4860: if (!mesh.commRank()) {fs << l_iter->first << std::endl;}
4861: LabelSifterSerializer::writeLabel(fs, *l_iter->second);
4862: }
4863: // Write sections
4864: Obj<std::set<std::string> > realNames = mesh.getRealSections();
4866: if (!mesh.commRank()) {fs << realNames->size() << std::endl;}
4867: for(std::set<std::string>::const_iterator n_iter = realNames->begin(); n_iter != realNames->end(); ++n_iter) {
4868: if (!mesh.commRank()) {fs << *n_iter << std::endl;}
4869: SectionSerializer::writeSection(fs, *mesh.getRealSection(*n_iter));
4870: }
4871: Obj<std::set<std::string> > intNames = mesh.getIntSections();
4873: if (!mesh.commRank()) {fs << intNames->size() << std::endl;}
4874: for(std::set<std::string>::const_iterator n_iter = intNames->begin(); n_iter != intNames->end(); ++n_iter) {
4875: if (!mesh.commRank()) {fs << *n_iter << std::endl;}
4876: SectionSerializer::writeSection(fs, *mesh.getIntSection(*n_iter));
4877: }
4878: // Write overlap
4879: SifterSerializer::writeSifter(fs, *mesh.getSendOverlap());
4880: SifterSerializer::writeSifter(fs, *mesh.getRecvOverlap());
4881: // Write distribution overlap
4882: // Write renumbering
4883: };
4884: template<typename Mesh>
4885: static void loadMesh(const std::string& filename, Mesh& mesh) {
4886: std::ifstream fs;
4888: if (mesh.commRank() == 0) {
4889: fs.open(filename.c_str());
4890: }
4891: loadMesh(fs, mesh);
4892: if (mesh.commRank() == 0) {
4893: fs.close();
4894: }
4895: };
4896: template<typename Mesh>
4897: static void loadMesh(std::ifstream& fs, Mesh& mesh) {
4898: ALE::Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh.comm(), mesh.debug());
4899: PetscErrorCode ierr;
4901: ISieveSerializer::loadSieve(fs, *sieve);
4902: mesh.setSieve(sieve);
4903: // Load labels
4904: int numLabels;
4906: if (!mesh.commRank()) {fs >> numLabels;}
4907: MPI_Bcast(&numLabels, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4908: for(int l = 0; l < numLabels; ++l) {
4909: ALE::Obj<typename Mesh::label_type> label = new typename Mesh::label_type(mesh.comm(), mesh.debug());
4910: std::string name;
4911: int len;
4913: if (!mesh.commRank()) {
4914: fs >> name;
4915: len = name.size();
4916: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4917: MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4918: } else {
4919: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4920: char *n = new char[len+1];
4921: MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4922: name = n;
4923: delete [] n;
4924: }
4925: LabelSifterSerializer::loadLabel(fs, *label);
4926: mesh.setLabel(name, label);
4927: }
4928: // Load sections
4929: int numRealSections;
4931: if (!mesh.commRank()) {fs >> numRealSections;}
4932: MPI_Bcast(&numRealSections, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4933: for(int s = 0; s < numRealSections; ++s) {
4934: ALE::Obj<typename Mesh::real_section_type> section = new typename Mesh::real_section_type(mesh.comm(), mesh.debug());
4935: std::string name;
4936: int len;
4938: if (!mesh.commRank()) {
4939: fs >> name;
4940: len = name.size();
4941: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4942: MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4943: } else {
4944: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4945: char *n = new char[len+1];
4946: MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4947: name = n;
4948: delete [] n;
4949: }
4950: SectionSerializer::loadSection(fs, *section);
4951: mesh.setRealSection(name, section);
4952: }
4953: int numIntSections;
4955: if (!mesh.commRank()) {fs >> numIntSections;}
4956: MPI_Bcast(&numIntSections, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4957: for(int s = 0; s < numIntSections; ++s) {
4958: ALE::Obj<typename Mesh::int_section_type> section = new typename Mesh::int_section_type(mesh.comm(), mesh.debug());
4959: std::string name;
4960: int len;
4962: if (!mesh.commRank()) {
4963: fs >> name;
4964: len = name.size();
4965: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4966: MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4967: } else {
4968: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4969: char *n = new char[len+1];
4970: MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4971: name = n;
4972: delete [] n;
4973: }
4974: SectionSerializer::loadSection(fs, *section);
4975: mesh.setIntSection(name, section);
4976: }
4977: // Load overlap
4978: SifterSerializer::loadSifter(fs, *mesh.getSendOverlap());
4979: SifterSerializer::loadSifter(fs, *mesh.getRecvOverlap());
4980: // Load distribution overlap
4981: // Load renumbering
4982: };
4983: };
4984: } // namespace ALE
4985: #endif