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 <sieve/Numbering.hh>
9: #endif
11: #ifndef included_ALE_INumbering_hh
12: #include <sieve/INumbering.hh>
13: #endif
15: #ifndef included_ALE_Field_hh
16: #include <sieve/Field.hh>
17: #endif
19: #ifndef included_ALE_IField_hh
20: #include <sieve/IField.hh>
21: #endif
23: #ifndef included_ALE_SieveBuilder_hh
24: #include <sieve/SieveBuilder.hh>
25: #endif
27: #ifndef included_ALE_LabelSifter_hh
28: #include <sieve/LabelSifter.hh>
29: #endif
31: #ifndef included_ALE_Partitioner_hh
32: #include <sieve/Partitioner.hh>
33: #endif
35: #ifndef included_ALE_Ordering_hh
36: #include <sieve/Ordering.hh>
37: #endif
39: #ifndef included_PETSc_Overlap_hh
40: #include <sieve/Overlap.hh>
41: #endif
43: namespace ALE {
44: class indexSet : public std::valarray<int> {
45: public:
46: inline bool
47: operator<(const indexSet& __x) {
48: if (__x.size() != this->size()) return __x.size() < this->size();
49: for(unsigned int i = 0; i < __x.size(); ++i) {
50: if (__x[i] == (*this)[i]) continue;
51: return __x[i] < (*this)[i];
52: }
53: return false;
54: }
55: };
56: inline bool
57: operator<(const indexSet& __x, const indexSet& __y) {
58: if (__x.size() != __y.size()) return __x.size() < __y.size();
59: for(unsigned int i = 0; i < __x.size(); ++i) {
60: if (__x[i] == __y[i]) continue;
61: return __x[i] < __y[i];
62: }
63: return false;
64: };
65: inline bool
66: operator<=(const indexSet& __x, const indexSet& __y) {
67: if (__x.size() != __y.size()) return __x.size() < __y.size();
68: for(unsigned int i = 0; i < __x.size(); ++i) {
69: if (__x[i] == __y[i]) continue;
70: return __x[i] < __y[i];
71: }
72: return true;
73: };
74: inline bool
75: operator==(const indexSet& __x, const indexSet& __y) {
76: if (__x.size() != __y.size()) return false;
77: for(unsigned int i = 0; i < __x.size(); ++i) {
78: if (__x[i] != __y[i]) return false;
79: }
80: return true;
81: };
82: inline bool
83: operator!=(const indexSet& __x, const indexSet& __y) {
84: if (__x.size() != __y.size()) return true;
85: for(unsigned int i = 0; i < __x.size(); ++i) {
86: if (__x[i] != __y[i]) return true;
87: }
88: return false;
89: };
91: template<typename Sieve_,
92: typename RealSection_ = Section<typename Sieve_::point_type, double>,
93: typename IntSection_ = Section<typename Sieve_::point_type, int>,
94: typename ArrowSection_ = UniformSection<MinimalArrow<typename Sieve_::point_type, typename Sieve_::point_type>, int> >
95: class Bundle : public ALE::ParallelObject {
96: public:
97: typedef Sieve_ sieve_type;
98: typedef RealSection_ real_section_type;
99: typedef IntSection_ int_section_type;
100: typedef ArrowSection_ arrow_section_type;
101: typedef Bundle<Sieve_,RealSection_,IntSection_,ArrowSection_> this_type;
102: typedef typename sieve_type::point_type point_type;
103: typedef malloc_allocator<point_type> alloc_type;
104: typedef typename ALE::LabelSifter<int, point_type> label_type;
105: typedef typename std::map<const std::string, Obj<label_type> > labels_type;
106: typedef typename label_type::supportSequence label_sequence;
107: typedef std::map<std::string, Obj<arrow_section_type> > arrow_sections_type;
108: typedef std::map<std::string, Obj<real_section_type> > real_sections_type;
109: typedef std::map<std::string, Obj<int_section_type> > int_sections_type;
110: typedef ALE::Point index_type;
111: typedef std::pair<index_type, int> oIndex_type;
112: typedef std::vector<oIndex_type> oIndexArray;
113: typedef std::pair<int *, int> indices_type;
114: typedef NumberingFactory<this_type> MeshNumberingFactory;
115: typedef typename ALE::Partitioner<>::part_type rank_type;
116: typedef typename ALE::Sifter<point_type,rank_type,point_type> send_overlap_type;
117: typedef typename ALE::Sifter<rank_type,point_type,point_type> recv_overlap_type;
118: typedef typename MeshNumberingFactory::numbering_type numbering_type;
119: typedef typename MeshNumberingFactory::order_type order_type;
120: typedef std::map<point_type, point_type> renumbering_type;
121: typedef typename ALE::SieveAlg<this_type> sieve_alg_type;
122: typedef typename sieve_alg_type::coneArray coneArray;
123: typedef typename sieve_alg_type::orientedConeArray oConeArray;
124: typedef typename sieve_alg_type::supportArray supportArray;
125: protected:
126: Obj<sieve_type> _sieve;
127: labels_type _labels;
128: int _maxHeight;
129: int _maxDepth;
130: arrow_sections_type _arrowSections;
131: real_sections_type _realSections;
132: int_sections_type _intSections;
133: Obj<oIndexArray> _indexArray;
134: Obj<MeshNumberingFactory> _factory;
135: bool _calculatedOverlap;
136: Obj<send_overlap_type> _sendOverlap;
137: Obj<recv_overlap_type> _recvOverlap;
138: Obj<send_overlap_type> _distSendOverlap;
139: Obj<recv_overlap_type> _distRecvOverlap;
140: renumbering_type _renumbering; // Maps global points to local points
141: // Work space
142: Obj<std::set<point_type> > _modifiedPoints;
143: public:
144: Bundle(MPI_Comm comm, int debug = 0) : ALE::ParallelObject(comm, debug), _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: Bundle(const Obj<sieve_type>& sieve) : ALE::ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _maxHeight(-1), _maxDepth(-1) {
153: this->_indexArray = new oIndexArray();
154: this->_modifiedPoints = new std::set<point_type>();
155: this->_factory = MeshNumberingFactory::singleton(this->comm(), this->debug());
156: this->_calculatedOverlap = false;
157: this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
158: this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
159: };
160: virtual ~Bundle() {};
161: public: // Verifiers
162: bool hasLabel(const std::string& name) {
163: if (this->_labels.find(name) != this->_labels.end()) {
164: return true;
165: }
166: return false;
167: };
168: void checkLabel(const std::string& name) {
169: if (!this->hasLabel(name)) {
170: ostringstream msg;
171: msg << "Invalid label name: " << name << std::endl;
172: throw ALE::Exception(msg.str().c_str());
173: }
174: };
175: public: // Accessors
176: const Obj<sieve_type>& getSieve() const {return this->_sieve;};
177: void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
178: bool hasArrowSection(const std::string& name) const {
179: return this->_arrowSections.find(name) != this->_arrowSections.end();
180: };
181: const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
182: if (!this->hasArrowSection(name)) {
183: Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());
185: section->setName(name);
186: if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
187: this->_arrowSections[name] = section;
188: }
189: return this->_arrowSections[name];
190: };
191: void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
192: this->_arrowSections[name] = section;
193: };
194: Obj<std::set<std::string> > getArrowSections() const {
195: Obj<std::set<std::string> > names = std::set<std::string>();
197: for(typename arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
198: names->insert(s_iter->first);
199: }
200: return names;
201: };
202: bool hasRealSection(const std::string& name) const {
203: return this->_realSections.find(name) != this->_realSections.end();
204: };
205: const Obj<real_section_type>& getRealSection(const std::string& name) {
206: if (!this->hasRealSection(name)) {
207: Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());
209: section->setName(name);
210: if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
211: this->_realSections[name] = section;
212: }
213: return this->_realSections[name];
214: };
215: void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
216: this->_realSections[name] = section;
217: };
218: Obj<std::set<std::string> > getRealSections() const {
219: Obj<std::set<std::string> > names = std::set<std::string>();
221: for(typename real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
222: names->insert(s_iter->first);
223: }
224: return names;
225: };
226: bool hasIntSection(const std::string& name) const {
227: return this->_intSections.find(name) != this->_intSections.end();
228: };
229: const Obj<int_section_type>& getIntSection(const std::string& name) {
230: if (!this->hasIntSection(name)) {
231: Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());
233: section->setName(name);
234: if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
235: this->_intSections[name] = section;
236: }
237: return this->_intSections[name];
238: };
239: void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
240: this->_intSections[name] = section;
241: };
242: Obj<std::set<std::string> > getIntSections() const {
243: Obj<std::set<std::string> > names = std::set<std::string>();
245: for(typename int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
246: names->insert(s_iter->first);
247: }
248: return names;
249: };
250: const Obj<MeshNumberingFactory>& getFactory() const {return this->_factory;};
251: bool getCalculatedOverlap() const {return this->_calculatedOverlap;};
252: void setCalculatedOverlap(const bool calc) {this->_calculatedOverlap = calc;};
253: const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
254: void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
255: const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
256: void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
257: const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
258: void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
259: const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
260: void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
261: renumbering_type& getRenumbering() {return this->_renumbering;};
262: public: // Labels
263: int getValue (const Obj<label_type>& label, const point_type& point, const int defValue = 0) {
264: const Obj<typename label_type::coneSequence>& cone = label->cone(point);
266: if (cone->size() == 0) return defValue;
267: return *cone->begin();
268: };
269: void setValue(const Obj<label_type>& label, const point_type& point, const int value) {
270: label->setCone(value, point);
271: };
272: template<typename InputPoints>
273: int getMaxValue (const Obj<label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
274: int maxValue = defValue;
276: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
277: maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
278: }
279: return maxValue;
280: }
281: const Obj<label_type>& createLabel(const std::string& name) {
282: this->_labels[name] = new label_type(this->comm(), this->debug());
283: return this->_labels[name];
284: };
285: const Obj<label_type>& getLabel(const std::string& name) {
286: this->checkLabel(name);
287: return this->_labels[name];
288: };
289: void setLabel(const std::string& name, const Obj<label_type>& label) {
290: this->_labels[name] = label;
291: };
292: const labels_type& getLabels() {
293: return this->_labels;
294: };
295: virtual const Obj<label_sequence>& getLabelStratum(const std::string& name, int value) {
296: this->checkLabel(name);
297: return this->_labels[name]->support(value);
298: };
299: public: // Stratification
300: template<class InputPoints>
301: void computeHeight(const Obj<label_type>& height, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxHeight) {
302: this->_modifiedPoints->clear();
304: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
305: // Compute the max height of the points in the support of p, and add 1
306: int h0 = this->getValue(height, *p_iter, -1);
307: int h1 = this->getMaxValue(height, sieve->support(*p_iter), -1) + 1;
309: if(h1 != h0) {
310: this->setValue(height, *p_iter, h1);
311: if (h1 > maxHeight) maxHeight = h1;
312: this->_modifiedPoints->insert(*p_iter);
313: }
314: }
315: // FIX: We would like to avoid the copy here with cone()
316: if(this->_modifiedPoints->size() > 0) {
317: this->computeHeight(height, sieve, sieve->cone(this->_modifiedPoints), maxHeight);
318: }
319: }
320: void computeHeights() {
321: const Obj<label_type>& label = this->createLabel(std::string("height"));
323: this->_maxHeight = -1;
324: this->computeHeight(label, this->_sieve, this->_sieve->leaves(), this->_maxHeight);
325: };
326: virtual int height() const {return this->_maxHeight;};
327: virtual int height(const point_type& point) {
328: return this->getValue(this->_labels["height"], point, -1);
329: };
330: virtual const Obj<label_sequence>& heightStratum(int height) {
331: return this->getLabelStratum("height", height);
332: };
333: void setHeight(const Obj<label_type>& label) {
334: this->_labels["height"] = label;
335: const Obj<typename label_type::traits::capSequence> cap = label->cap();
337: for(typename label_type::traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
338: this->_maxHeight = std::max(this->_maxHeight, *c_iter);
339: }
340: };
341: template<class InputPoints>
342: void computeDepth(const Obj<label_type>& depth, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxDepth) {
343: this->_modifiedPoints->clear();
345: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
346: // Compute the max depth of the points in the cone of p, and add 1
347: int d0 = this->getValue(depth, *p_iter, -1);
348: int d1 = this->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;
350: if(d1 != d0) {
351: this->setValue(depth, *p_iter, d1);
352: if (d1 > maxDepth) maxDepth = d1;
353: this->_modifiedPoints->insert(*p_iter);
354: }
355: }
356: // FIX: We would like to avoid the copy here with support()
357: if(this->_modifiedPoints->size() > 0) {
358: this->computeDepth(depth, sieve, sieve->support(this->_modifiedPoints), maxDepth);
359: }
360: }
361: void computeDepths() {
362: const Obj<label_type>& label = this->createLabel(std::string("depth"));
364: this->_maxDepth = -1;
365: this->computeDepth(label, this->_sieve, this->_sieve->roots(), this->_maxDepth);
366: };
367: virtual int depth() const {return this->_maxDepth;};
368: virtual int depth(const point_type& point) {
369: return this->getValue(this->_labels["depth"], point, -1);
370: };
371: virtual const Obj<label_sequence>& depthStratum(int depth) {
372: return this->getLabelStratum("depth", depth);
373: };
376: virtual void stratify() {
377: ALE_LOG_EVENT_BEGIN;
378: this->computeHeights();
379: this->computeDepths();
380: ALE_LOG_EVENT_END;
381: };
382: public: // Size traversal
383: template<typename Section_>
384: int size(const Obj<Section_>& section, const point_type& p) {
385: const typename Section_::chart_type& chart = section->getChart();
386: int size = 0;
388: if (this->height() < 2) {
389: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
390: typename sieve_type::coneSequence::iterator end = cone->end();
392: if (chart.count(p)) {
393: size += section->getConstrainedFiberDimension(p);
394: }
395: for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
396: if (chart.count(*c_iter)) {
397: size += section->getConstrainedFiberDimension(*c_iter);
398: }
399: }
400: } else {
401: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), p);
402: typename coneArray::iterator end = closure->end();
404: for(typename coneArray::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
405: if (chart.count(*c_iter)) {
406: size += section->getConstrainedFiberDimension(*c_iter);
407: }
408: }
409: }
410: return size;
411: }
412: template<typename Section_>
413: int sizeWithBC(const Obj<Section_>& section, const point_type& p) {
414: const typename Section_::chart_type& chart = section->getChart();
415: int size = 0;
417: if (this->height() < 2) {
418: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
419: typename sieve_type::coneSequence::iterator end = cone->end();
421: if (chart.count(p)) {
422: size += section->getFiberDimension(p);
423: }
424: for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
425: if (chart.count(*c_iter)) {
426: size += section->getFiberDimension(*c_iter);
427: }
428: }
429: } else {
430: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), p);
431: typename coneArray::iterator end = closure->end();
433: for(typename coneArray::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
434: if (chart.count(*c_iter)) {
435: size += section->getFiberDimension(*c_iter);
436: }
437: }
438: }
439: return size;
440: }
441: protected:
442: int *getIndexArray(const int size) {
443: static int *array = NULL;
444: static int maxSize = 0;
446: if (size > maxSize) {
447: maxSize = size;
448: if (array) delete [] array;
449: array = new int[maxSize];
450: };
451: return array;
452: };
453: public: // Index traversal
454: void expandInterval(const index_type& interval, PetscInt indices[], PetscInt *indx) {
455: const int end = interval.prefix + interval.index;
457: for(int i = interval.index; i < end; ++i) {
458: indices[(*indx)++] = i;
459: }
460: };
461: void expandInterval(const index_type& interval, const int orientation, PetscInt indices[], PetscInt *indx) {
462: if (orientation >= 0) {
463: for(int i = 0; i < interval.prefix; ++i) {
464: indices[(*indx)++] = interval.index + i;
465: }
466: } else {
467: for(int i = interval.prefix-1; i >= 0; --i) {
468: indices[(*indx)++] = interval.index + i;
469: }
470: }
471: for(int i = 0; i < -interval.prefix; ++i) {
472: indices[(*indx)++] = -1;
473: }
474: };
475: void expandIntervals(Obj<oIndexArray> intervals, PetscInt *indices) {
476: int k = 0;
478: for(typename oIndexArray::iterator i_iter = intervals->begin(); i_iter != intervals->end(); ++i_iter) {
479: this->expandInterval(i_iter->first, i_iter->second, indices, &k);
480: }
481: }
482: template<typename Section_>
483: const indices_type getIndicesRaw(const Obj<Section_>& section, const point_type& p) {
484: int *indexArray = NULL;
485: int size = 0;
487: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
488: typename oConeArray::iterator begin = closure->begin();
489: typename oConeArray::iterator end = closure->end();
491: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
492: size += section->getFiberDimension(p_iter->first);
493: }
494: indexArray = this->getIndexArray(size);
495: int k = 0;
496: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
497: section->getIndicesRaw(p_iter->first, section->getIndex(p_iter->first), indexArray, &k, p_iter->second);
498: }
499: return indices_type(indexArray, size);
500: }
501: template<typename Section_>
502: const indices_type getIndices(const Obj<Section_>& section, const point_type& p, const int level = -1) {
503: int *indexArray = NULL;
504: int size = 0;
506: if (level == 0) {
507: size += section->getFiberDimension(p);
508: indexArray = this->getIndexArray(size);
509: int k = 0;
511: section->getIndices(p, indexArray, &k);
512: } else if ((level == 1) || (this->height() == 1)) {
513: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
514: typename sieve_type::coneSequence::iterator end = cone->end();
516: size += section->getFiberDimension(p);
517: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
518: size += section->getFiberDimension(*p_iter);
519: }
520: indexArray = this->getIndexArray(size);
521: int k = 0;
523: section->getIndices(p, indexArray, &k);
524: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
525: section->getIndices(*p_iter, indexArray, &k);
526: }
527: } else if (level == -1) {
528: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
529: typename oConeArray::iterator begin = closure->begin();
530: typename oConeArray::iterator end = closure->end();
532: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
533: size += section->getFiberDimension(p_iter->first);
534: }
535: indexArray = this->getIndexArray(size);
536: int k = 0;
537: for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
538: section->getIndices(p_iter->first, indexArray, &k, p_iter->second);
539: }
540: } else {
541: throw ALE::Exception("Bundle has not yet implemented getIndices() for an arbitrary level");
542: }
543: if (this->debug()) {
544: for(int i = 0; i < size; ++i) {
545: printf("[%d]index %d: %d\n", this->commRank(), i, indexArray[i]);
546: }
547: }
548: return indices_type(indexArray, size);
549: }
550: template<typename Section_, typename Numbering>
551: const indices_type getIndices(const Obj<Section_>& section, const point_type& p, const Obj<Numbering>& numbering, const int level = -1) {
552: int *indexArray = NULL;
553: int size = 0;
555: if (level == 0) {
556: size += section->getFiberDimension(p);
557: indexArray = this->getIndexArray(size);
558: int k = 0;
560: section->getIndices(p, numbering, indexArray, &k);
561: } else if ((level == 1) || (this->height() == 1)) {
562: const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
563: typename sieve_type::coneSequence::iterator end = cone->end();
565: size += section->getFiberDimension(p);
566: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
567: size += section->getFiberDimension(*p_iter);
568: }
569: indexArray = this->getIndexArray(size);
570: int k = 0;
572: section->getIndices(p, numbering, indexArray, &k);
573: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
574: section->getIndices(*p_iter, numbering, indexArray, &k);
575: }
576: } else if (level == -1) {
577: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
578: typename oConeArray::iterator end = closure->end();
580: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
581: size += section->getFiberDimension(p_iter->first);
582: }
583: indexArray = this->getIndexArray(size);
584: int k = 0;
585: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
586: section->getIndices(p_iter->first, numbering, indexArray, &k, p_iter->second);
587: }
588: } else {
589: throw ALE::Exception("Bundle has not yet implemented getIndices() for an arbitrary level");
590: }
591: return indices_type(indexArray, size);
592: }
593: public: // Retrieval traversal
594: // Return the values for the closure of this point
595: // use a smart pointer?
596: template<typename Section_>
597: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p) {
598: const int size = this->sizeWithBC(section, p);
599: return this->restrictClosure(section, p, section->getRawArray(size), size);
600: }
601: template<typename Section_>
602: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p, typename Section_::value_type *values, const int valuesSize) {
603: const int size = this->sizeWithBC(section, p);
604: int j = -1;
605: if (valuesSize < size) throw ALE::Exception("Input array too small");
607: // We could actually ask for the height of the individual point
608: if (this->height() < 2) {
609: const int& dim = section->getFiberDimension(p);
610: const typename Section_::value_type *array = section->restrictPoint(p);
612: for(int i = 0; i < dim; ++i) {
613: values[++j] = array[i];
614: }
615: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
616: typename sieve_type::coneSequence::iterator end = cone->end();
618: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
619: const int& dim = section->getFiberDimension(*p_iter);
621: array = section->restrictPoint(*p_iter);
622: for(int i = 0; i < dim; ++i) {
623: values[++j] = array[i];
624: }
625: }
626: } else {
627: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
628: typename oConeArray::iterator end = closure->end();
630: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
631: const typename Section_::value_type *array = section->restrictPoint(p_iter->first);
632: const int& dim = section->getFiberDimension(p_iter->first);
634: if (p_iter->second >= 0) {
635: for(int i = 0; i < dim; ++i) {
636: values[++j] = array[i];
637: }
638: } else {
639: for(int i = dim-1; i >= 0; --i) {
640: values[++j] = array[i];
641: }
642: }
643: }
644: }
645: if (j != size-1) {
646: ostringstream txt;
648: txt << "Invalid restrict to point " << p << std::endl;
649: txt << " j " << j << " should be " << (size-1) << std::endl;
650: std::cout << txt.str();
651: throw ALE::Exception(txt.str().c_str());
652: }
653: return values;
654: }
655: template<typename Section_>
656: const typename Section_::value_type *restrictNew(const Obj<Section_>& section, const point_type& p) {
657: const int size = this->sizeWithBC(section, p);
658: return this->restrictNew(section, p, section->getRawArray(size), size);
659: }
660: template<typename Section_>
661: const typename Section_::value_type *restrictNew(const Obj<Section_>& section, const point_type& p, typename Section_::value_type *values, const int valuesSize) {
662: const int size = this->sizeWithBC(section, p);
663: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
664: typename oConeArray::iterator end = closure->end();
665: int j = -1;
666: if (valuesSize < size) throw ALE::Exception("Input array too small");
668: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
669: const typename Section_::value_type *array = section->restrictPoint(p_iter->first);
671: if (p_iter->second >= 0) {
672: const int& dim = section->getFiberDimension(p_iter->first);
674: for(int i = 0; i < dim; ++i) {
675: values[++j] = array[i];
676: }
677: } else {
678: int offset = 0;
680: for(int space = 0; space < section->getNumSpaces(); ++space) {
681: const int& dim = section->getFiberDimension(p_iter->first, space);
683: for(int i = dim-1; i >= 0; --i) {
684: values[++j] = array[i+offset];
685: }
686: offset += dim;
687: }
688: }
689: }
690: if (j != size-1) {
691: ostringstream txt;
693: txt << "Invalid restrict to point " << p << std::endl;
694: txt << " j " << j << " should be " << (size-1) << std::endl;
695: std::cout << txt.str();
696: throw ALE::Exception(txt.str().c_str());
697: }
698: return values;
699: }
700: template<typename Section_>
701: void update(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
702: int j = 0;
704: if (this->height() < 2) {
705: section->updatePoint(p, &v[j]);
706: j += section->getFiberDimension(p);
707: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
709: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
710: section->updatePoint(*p_iter, &v[j]);
711: j += section->getFiberDimension(*p_iter);
712: }
713: } else {
714: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
715: typename oConeArray::iterator end = closure->end();
717: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
718: section->updatePoint(p_iter->first, &v[j], p_iter->second);
719: j += section->getFiberDimension(p_iter->first);
720: }
721: }
722: }
723: template<typename Section_>
724: void updateAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
725: int j = 0;
727: if (this->height() < 2) {
728: section->updateAddPoint(p, &v[j]);
729: j += section->getFiberDimension(p);
730: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
732: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
733: section->updateAddPoint(*p_iter, &v[j]);
734: j += section->getFiberDimension(*p_iter);
735: }
736: } else {
737: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
738: typename oConeArray::iterator end = closure->end();
740: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
741: section->updateAddPoint(p_iter->first, &v[j], p_iter->second);
742: j += section->getFiberDimension(p_iter->first);
743: }
744: }
745: }
746: template<typename Section_>
747: void updateBC(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
748: int j = 0;
750: if (this->height() < 2) {
751: section->updatePointBC(p, &v[j]);
752: j += section->getFiberDimension(p);
753: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
755: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
756: section->updatePointBC(*p_iter, &v[j]);
757: j += section->getFiberDimension(*p_iter);
758: }
759: } else {
760: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
761: typename oConeArray::iterator end = closure->end();
763: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
764: section->updatePointBC(p_iter->first, &v[j], p_iter->second);
765: j += section->getFiberDimension(p_iter->first);
766: }
767: }
768: }
769: template<typename Section_>
770: void updateAll(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
771: int j = 0;
773: if (this->height() < 2) {
774: section->updatePointAll(p, &v[j]);
775: j += section->getFiberDimension(p);
776: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
778: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
779: section->updatePointAll(*p_iter, &v[j]);
780: j += section->getFiberDimension(*p_iter);
781: }
782: } else {
783: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
784: typename oConeArray::iterator end = closure->end();
786: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
787: section->updatePointAll(p_iter->first, &v[j], p_iter->second);
788: j += section->getFiberDimension(p_iter->first);
789: }
790: }
791: }
792: template<typename Section_>
793: void updateAllAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
794: int j = 0;
796: if (this->height() < 2) {
797: section->updatePointAllAdd(p, &v[j]);
798: j += section->getFiberDimension(p);
799: const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
801: for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
802: section->updatePointAllAdd(*p_iter, &v[j]);
803: j += section->getFiberDimension(*p_iter);
804: }
805: } else {
806: const Obj<oConeArray> closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
807: typename oConeArray::iterator end = closure->end();
809: for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
810: section->updatePointAllAdd(p_iter->first, &v[j], p_iter->second);
811: j += section->getFiberDimension(p_iter->first);
812: }
813: }
814: }
815: public: // Optimization
816: // Calculate a custom atlas for the given traversal
817: // This returns the tag value assigned to the traversal
818: template<typename Section_, typename Sequence_>
819: int calculateCustomAtlas(const Obj<Section_>& section, const Obj<Sequence_>& points) {
820: const typename Sequence_::iterator begin = points->begin();
821: const typename Sequence_::iterator end = points->end();
822: const int num = points->size();
823: int *rOffsets = new int[num+1];
824: int *rIndices;
825: int *uOffsets = new int[num+1];
826: int *uIndices;
827: int p;
829: p = 0;
830: rOffsets[p] = 0;
831: uOffsets[p] = 0;
832: for(typename Sequence_::iterator p_iter = begin; p_iter != end; ++p_iter, ++p) {
833: rOffsets[p+1] = rOffsets[p] + this->sizeWithBC(section, *p_iter);
834: uOffsets[p+1] = rOffsets[p+1];
835: //uOffsets[p+1] = uOffsets[p] + this->size(section, *p_iter);
836: }
837: rIndices = new int[rOffsets[p]];
838: uIndices = new int[uOffsets[p]];
839: p = 0;
840: for(typename Sequence_::iterator p_iter = begin; p_iter != end; ++p_iter, ++p) {
841: const indices_type rIdx = this->getIndicesRaw(section, *p_iter);
842: for(int i = 0, k = rOffsets[p]; k < rOffsets[p+1]; ++i, ++k) rIndices[k] = rIdx.first[i];
844: const indices_type uIdx = this->getIndices(section, *p_iter);
845: for(int i = 0, k = uOffsets[p]; k < uOffsets[p+1]; ++i, ++k) uIndices[k] = uIdx.first[i];
846: }
847: return section->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices);
848: }
849: template<typename Section_>
850: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p) {
851: const int *offsets, *indices;
853: section->getCustomRestrictAtlas(tag, &offsets, &indices);
854: const int size = offsets[p+1] - offsets[p];
855: return this->restrictClosure(section, tag, p, section->getRawArray(size), offsets, indices);
856: }
857: template<typename Section_>
858: const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p, typename Section_::value_type *values, const int valuesSize) {
859: const int *offsets, *indices;
861: section->getCustomRestrictAtlas(tag, &offsets, &indices);
862: const int size = offsets[p+1] - offsets[p];
863: if (valuesSize < size) {throw ALE::Exception("Input array too small");}
864: return this->restrictClosure(section, tag, p, values, offsets, indices);
865: }
866: template<typename Section_>
867: 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[]) {
868: const typename Section_::value_type *array = section->restrictSpace();
870: const int size = offsets[p+1] - offsets[p];
871: for(int j = 0, k = offsets[p]; j < size; ++j, ++k) {
872: values[j] = array[indices[k]];
873: }
874: return values;
875: }
876: template<typename Section_>
877: void updateAdd(const Obj<Section_>& section, const int tag, const int p, const typename Section_::value_type values[]) {
878: typename Section_::value_type *array = (typename Section_::value_type *) section->restrictSpace();
879: const int *offsets, *indices;
881: section->getCustomUpdateAtlas(tag, &offsets, &indices);
882: const int size = offsets[p+1] - offsets[p];
883: for(int j = 0, k = offsets[p]; j < size; ++j, ++k) {
884: if (indices[k] < 0) continue;
885: array[indices[k]] += values[j];
886: }
887: }
888: public: // Allocation
889: template<typename Section_>
890: void allocate(const Obj<Section_>& section, const Obj<send_overlap_type>& sendOverlap = NULL) {
891: bool doGhosts = !sendOverlap.isNull();
893: this->_factory->orderPatch(section, this->getSieve(), sendOverlap);
894: if (doGhosts) {
895: if (this->_debug > 1) {std::cout << "Ordering patch for ghosts" << std::endl;}
896: const typename Section_::chart_type& points = section->getChart();
897: typename Section_::index_type::index_type offset = 0;
899: for(typename Section_::chart_type::const_iterator point = points.begin(); point != points.end(); ++point) {
900: const typename Section_::index_type& idx = section->getIndex(*point);
902: offset = std::max(offset, idx.index + std::abs(idx.prefix));
903: }
904: this->_factory->orderPatch(section, this->getSieve(), NULL, offset);
905: if (offset != section->sizeWithBC()) throw ALE::Exception("Inconsistent array sizes in section");
906: }
907: section->allocateStorage();
908: }
909: template<typename Section_>
910: void reallocate(const Obj<Section_>& section) {
911: if (section->getNewAtlas().isNull()) return;
912: // Since copy() preserves offsets, we must reinitialize them before ordering
913: const Obj<typename Section_::atlas_type> atlas = section->getAtlas();
914: const Obj<typename Section_::atlas_type>& newAtlas = section->getNewAtlas();
915: const typename Section_::atlas_type::chart_type& chart = newAtlas->getChart();
916: const typename Section_::atlas_type::chart_type& oldChart = atlas->getChart();
917: int newSize = 0;
918: typename Section_::index_type defaultIdx(0, -1);
920: for(typename Section_::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
921: defaultIdx.prefix = newAtlas->restrictPoint(*c_iter)[0].prefix;
922: newAtlas->updatePoint(*c_iter, &defaultIdx);
923: newSize += defaultIdx.prefix;
924: }
925: section->setAtlas(newAtlas);
926: this->_factory->orderPatch(section, this->getSieve());
927: // Copy over existing values
928: typedef typename alloc_type::template rebind<typename Section_::value_type>::other value_alloc_type;
929: value_alloc_type value_allocator;
930: typename Section_::value_type *newArray = value_allocator.allocate(newSize);
931: for(int i = 0; i < newSize; ++i) {value_allocator.construct(newArray+i, typename Section_::value_type());}
932: ///typename Section_::value_type *newArray = new typename Section_::value_type[newSize];
933: const typename Section_::value_type *array = section->restrictSpace();
935: for(typename Section_::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
936: const int& dim = section->getFiberDimension(*c_iter);
937: const int& offset = atlas->restrictPoint(*c_iter)->index;
938: const int& newOffset = newAtlas->restrictPoint(*c_iter)->index;
940: for(int i = 0; i < dim; ++i) {
941: newArray[newOffset+i] = array[offset+i];
942: }
943: }
944: section->replaceStorage(newArray);
945: }
946: public: // Overlap
947: template<typename Sequence>
948: void constructOverlap(const Obj<Sequence>& points, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
949: point_type *sendBuf = new point_type[points->size()];
950: int size = 0;
951: for(typename Sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
952: sendBuf[size++] = *l_iter;
953: }
954: int *sizes = new int[this->commSize()]; // The number of points coming from each process
955: int *offsets = new int[this->commSize()+1]; // Prefix sums for sizes
956: int *oldOffs = new int[this->commSize()+1]; // Temporary storage
957: point_type *remotePoints = NULL; // The points from each process
958: int *remoteRanks = NULL; // The rank and number of overlap points of each process that overlaps another
960: // Change to Allgather() for the correct binning algorithm
961: MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, this->comm());
962: if (this->commRank() == 0) {
963: offsets[0] = 0;
964: for(int p = 1; p <= this->commSize(); p++) {
965: offsets[p] = offsets[p-1] + sizes[p-1];
966: }
967: remotePoints = new point_type[offsets[this->commSize()]];
968: }
969: MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, this->comm());
970: delete [] sendBuf;
971: std::map<int, std::map<int, std::set<point_type> > > overlapInfo; // Maps (p,q) to their set of overlap points
973: if (this->commRank() == 0) {
974: for(int p = 0; p < this->commSize(); p++) {
975: std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]]);
976: }
977: for(int p = 0; p <= this->commSize(); p++) {
978: oldOffs[p] = offsets[p];
979: }
980: for(int p = 0; p < this->commSize(); p++) {
981: for(int q = p+1; q < this->commSize(); q++) {
982: std::set_intersection(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
983: &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
984: std::insert_iterator<std::set<point_type> >(overlapInfo[p][q], overlapInfo[p][q].begin()));
985: overlapInfo[q][p] = overlapInfo[p][q];
986: }
987: sizes[p] = overlapInfo[p].size()*2;
988: offsets[p+1] = offsets[p] + sizes[p];
989: }
990: remoteRanks = new int[offsets[this->commSize()]];
991: int k = 0;
992: for(int p = 0; p < this->commSize(); p++) {
993: for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
994: remoteRanks[k*2] = r_iter->first;
995: remoteRanks[k*2+1] = r_iter->second.size();
996: k++;
997: }
998: }
999: }
1000: int numOverlaps; // The number of processes overlapping this process
1001: MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, this->comm());
1002: int *overlapRanks = new int[numOverlaps]; // The rank and overlap size for each overlapping process
1003: MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, this->comm());
1004: point_type *sendPoints = NULL; // The points to send to each process
1005: if (this->commRank() == 0) {
1006: for(int p = 0, k = 0; p < this->commSize(); p++) {
1007: sizes[p] = 0;
1008: for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
1009: sizes[p] += remoteRanks[k*2+1];
1010: k++;
1011: }
1012: offsets[p+1] = offsets[p] + sizes[p];
1013: }
1014: sendPoints = new point_type[offsets[this->commSize()]];
1015: for(int p = 0, k = 0; p < this->commSize(); p++) {
1016: for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
1017: int rank = r_iter->first;
1018: for(typename std::set<point_type>::iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
1019: sendPoints[k++] = *p_iter;
1020: }
1021: }
1022: }
1023: }
1024: int numOverlapPoints = 0;
1025: for(int r = 0; r < numOverlaps/2; r++) {
1026: numOverlapPoints += overlapRanks[r*2+1];
1027: }
1028: point_type *overlapPoints = new point_type[numOverlapPoints];
1029: MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints, MPI_INT, 0, this->comm());
1031: for(int r = 0, k = 0; r < numOverlaps/2; r++) {
1032: int rank = overlapRanks[r*2];
1034: for(int p = 0; p < overlapRanks[r*2+1]; p++) {
1035: point_type point = overlapPoints[k++];
1037: sendOverlap->addArrow(point, rank, point);
1038: recvOverlap->addArrow(rank, point, point);
1039: }
1040: }
1042: delete [] overlapPoints;
1043: delete [] overlapRanks;
1044: delete [] sizes;
1045: delete [] offsets;
1046: delete [] oldOffs;
1047: if (this->commRank() == 0) {
1048: delete [] remoteRanks;
1049: delete [] remotePoints;
1050: delete [] sendPoints;
1051: }
1052: }
1053: void constructOverlap() {
1054: if (this->_calculatedOverlap) return;
1055: this->constructOverlap(this->getSieve()->base(), this->getSendOverlap(), this->getRecvOverlap());
1056: this->constructOverlap(this->getSieve()->cap(), this->getSendOverlap(), this->getRecvOverlap());
1057: if (this->debug()) {
1058: this->_sendOverlap->view("Send overlap");
1059: this->_recvOverlap->view("Receive overlap");
1060: }
1061: this->_calculatedOverlap = true;
1062: }
1063: };
1064: class BoundaryCondition : public ALE::ParallelObject {
1065: public:
1066: typedef double (*function_type)(const PetscReal []);
1067: typedef double (*integrator_type)(const PetscReal [], const PetscReal [], const int, function_type);
1068: protected:
1069: std::string _labelName;
1070: int _marker;
1071: function_type _func;
1072: integrator_type _integrator;
1073: public:
1074: BoundaryCondition(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
1075: ~BoundaryCondition() {};
1076: public:
1077: const std::string& getLabelName() const {return this->_labelName;};
1078: void setLabelName(const std::string& name) {this->_labelName = name;};
1079: int getMarker() const {return this->_marker;};
1080: void setMarker(const int marker) {this->_marker = marker;};
1081: function_type getFunction() const {return this->_func;};
1082: void setFunction(function_type func) {this->_func = func;};
1083: integrator_type getDualIntegrator() const {return this->_integrator;};
1084: void setDualIntegrator(integrator_type integrator) {this->_integrator = integrator;};
1085: public:
1086: PetscReal evaluate(const PetscReal coords[]) const {return this->_func(coords);};
1087: PetscReal integrateDual(const PetscReal v0[], const PetscReal J[], const int dualIndex) const {return this->_integrator(v0, J, dualIndex, this->_func);};
1088: };
1089: class Discretization : public ALE::ParallelObject {
1090: typedef std::map<std::string, Obj<BoundaryCondition> > boundaryConditions_type;
1091: protected:
1092: boundaryConditions_type _boundaryConditions;
1093: Obj<BoundaryCondition> _exactSolution;
1094: std::map<int,int> _dim2dof;
1095: std::map<int,int> _dim2class;
1096: int _quadSize;
1097: const double *_points;
1098: const double *_weights;
1099: int _basisSize;
1100: const double *_basis;
1101: const double *_basisDer;
1102: const int *_indices;
1103: std::map<int, const int *> _exclusionIndices;
1104: public:
1105: 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) {};
1106: virtual ~Discretization() {
1107: if (this->_indices) {delete [] this->_indices;}
1108: for(std::map<int, const int *>::iterator i_iter = _exclusionIndices.begin(); i_iter != _exclusionIndices.end(); ++i_iter) {
1109: delete [] i_iter->second;
1110: }
1111: };
1112: public:
1113: bool hasBoundaryCondition() {return (this->_boundaryConditions.find("default") != this->_boundaryConditions.end());};
1114: const Obj<BoundaryCondition>& getBoundaryCondition() {return this->getBoundaryCondition("default");};
1115: void setBoundaryCondition(const Obj<BoundaryCondition>& boundaryCondition) {this->setBoundaryCondition("default", boundaryCondition);};
1116: const Obj<BoundaryCondition>& getBoundaryCondition(const std::string& name) {return this->_boundaryConditions[name];};
1117: void setBoundaryCondition(const std::string& name, const Obj<BoundaryCondition>& boundaryCondition) {this->_boundaryConditions[name] = boundaryCondition;};
1118: Obj<std::set<std::string> > getBoundaryConditions() const {
1119: Obj<std::set<std::string> > names = std::set<std::string>();
1121: for(boundaryConditions_type::const_iterator d_iter = this->_boundaryConditions.begin(); d_iter != this->_boundaryConditions.end(); ++d_iter) {
1122: names->insert(d_iter->first);
1123: }
1124: return names;
1125: };
1126: const Obj<BoundaryCondition>& getExactSolution() {return this->_exactSolution;};
1127: void setExactSolution(const Obj<BoundaryCondition>& exactSolution) {this->_exactSolution = exactSolution;};
1128: int getQuadratureSize() {return this->_quadSize;};
1129: void setQuadratureSize(const int size) {this->_quadSize = size;};
1130: const double *getQuadraturePoints() {return this->_points;};
1131: void setQuadraturePoints(const double *points) {this->_points = points;};
1132: const double *getQuadratureWeights() {return this->_weights;};
1133: void setQuadratureWeights(const double *weights) {this->_weights = weights;};
1134: int getBasisSize() {return this->_basisSize;};
1135: void setBasisSize(const int size) {this->_basisSize = size;};
1136: const double *getBasis() {return this->_basis;};
1137: void setBasis(const double *basis) {this->_basis = basis;};
1138: const double *getBasisDerivatives() {return this->_basisDer;};
1139: void setBasisDerivatives(const double *basisDer) {this->_basisDer = basisDer;};
1140: int getNumDof(const int dim) {return this->_dim2dof[dim];};
1141: void setNumDof(const int dim, const int numDof) {this->_dim2dof[dim] = numDof;};
1142: int getDofClass(const int dim) {return this->_dim2class[dim];};
1143: void setDofClass(const int dim, const int dofClass) {this->_dim2class[dim] = dofClass;};
1144: public:
1145: const int *getIndices() {return this->_indices;};
1146: const int *getIndices(const int marker) {
1147: if (!marker) return this->getIndices();
1148: return this->_exclusionIndices[marker];
1149: };
1150: void setIndices(const int *indices) {this->_indices = indices;};
1151: void setIndices(const int *indices, const int marker) {
1152: if (!marker) this->_indices = indices;
1153: this->_exclusionIndices[marker] = indices;
1154: };
1155: template<typename Bundle>
1156: int sizeV(Bundle& mesh) {
1157: typedef typename ISieveVisitor::PointRetriever<typename Bundle::sieve_type> Visitor;
1158: Visitor pV((int) pow((double) mesh.getSieve()->getMaxConeSize(), mesh.depth())+1, true);
1159: ISieveTraversal<typename Bundle::sieve_type>::orientedClosure(*mesh.getSieve(), *mesh.heightStratum(0)->begin(), pV);
1160: const typename Visitor::point_type *oPoints = pV.getPoints();
1161: const int oSize = pV.getSize();
1162: int size = 0;
1164: for(int cl = 0; cl < oSize; ++cl) {
1165: size += this->_dim2dof[mesh.depth(oPoints[cl])];
1166: }
1167: return size;
1168: }
1169: template<typename Bundle>
1170: int size(const Obj<Bundle>& mesh) {
1171: const Obj<typename Bundle::label_sequence>& cells = mesh->heightStratum(0);
1172: const Obj<typename Bundle::coneArray> closure = ALE::SieveAlg<Bundle>::closure(mesh, *cells->begin());
1173: const typename Bundle::coneArray::iterator end = closure->end();
1174: int size = 0;
1176: for(typename Bundle::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1177: size += this->_dim2dof[mesh->depth(*cl_iter)];
1178: }
1179: return size;
1180: }
1181: };
1182: }
1184: namespace ALE {
1185: template<typename Sieve_,
1186: typename RealSection_ = Section<typename Sieve_::point_type, double>,
1187: typename IntSection_ = Section<typename Sieve_::point_type, int>,
1188: typename Label_ = LabelSifter<int, typename Sieve_::point_type>,
1189: typename ArrowSection_ = UniformSection<MinimalArrow<typename Sieve_::point_type, typename Sieve_::point_type>, int> >
1190: class IBundle : public ALE::ParallelObject {
1191: public:
1192: typedef Sieve_ sieve_type;
1193: typedef RealSection_ real_section_type;
1194: typedef IntSection_ int_section_type;
1195: typedef ArrowSection_ arrow_section_type;
1196: typedef IBundle<Sieve_,RealSection_,IntSection_,Label_,ArrowSection_> this_type;
1197: typedef typename sieve_type::point_type point_type;
1198: typedef malloc_allocator<point_type> alloc_type;
1199: typedef Label_ label_type;
1200: typedef typename std::map<const std::string, Obj<label_type> > labels_type;
1201: typedef typename label_type::supportSequence label_sequence;
1202: typedef std::map<std::string, Obj<arrow_section_type> > arrow_sections_type;
1203: typedef std::map<std::string, Obj<real_section_type> > real_sections_type;
1204: typedef std::map<std::string, Obj<int_section_type> > int_sections_type;
1205: typedef ALE::Point index_type;
1206: typedef std::pair<index_type, int> oIndex_type;
1207: typedef std::vector<oIndex_type> oIndexArray;
1208: typedef std::pair<int *, int> indices_type;
1209: typedef NumberingFactory<this_type> MeshNumberingFactory;
1210: typedef typename ALE::Partitioner<>::part_type rank_type;
1211: #ifdef USE_NEW_OVERLAP
1212: typedef typename PETSc::SendOverlap<point_type,rank_type> send_overlap_type;
1213: typedef typename PETSc::RecvOverlap<point_type,rank_type> recv_overlap_type;
1214: #else
1215: typedef typename ALE::Sifter<point_type,rank_type,point_type> send_overlap_type;
1216: typedef typename ALE::Sifter<rank_type,point_type,point_type> recv_overlap_type;
1217: #endif
1218: typedef typename MeshNumberingFactory::numbering_type numbering_type;
1219: typedef typename MeshNumberingFactory::order_type order_type;
1220: typedef std::map<point_type, point_type> renumbering_type;
1221: // These should go away
1222: typedef typename ALE::SieveAlg<this_type> sieve_alg_type;
1223: typedef typename sieve_alg_type::coneArray coneArray;
1224: typedef typename sieve_alg_type::orientedConeArray oConeArray;
1225: typedef typename sieve_alg_type::supportArray supportArray;
1226: public:
1227: class LabelVisitor {
1228: protected:
1229: label_type& label;
1230: int defaultValue;
1231: int value;
1232: public:
1233: LabelVisitor(label_type& l, const int defValue) : label(l), defaultValue(defValue), value(defValue) {};
1234: int getLabelValue(const point_type& point) const {
1235: const Obj<typename label_type::coneSequence>& cone = this->label.cone(point);
1237: if (cone->size() == 0) return this->defaultValue;
1238: return *cone->begin();
1239: };
1240: void setLabelValue(const point_type& point, const int value) {
1241: this->label.setCone(value, point);
1242: };
1243: int getValue() const {return this->value;};
1244: };
1245: class MaxConeVisitor : public LabelVisitor {
1246: public:
1247: MaxConeVisitor(label_type& l, const int defValue) : LabelVisitor(l, defValue) {};
1248: void visitPoint(const typename sieve_type::point_type& point) {};
1249: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1250: this->value = std::max(this->value, this->getLabelValue(arrow.source));
1251: };
1252: };
1253: class MaxSupportVisitor : public LabelVisitor {
1254: public:
1255: MaxSupportVisitor(label_type& l, const int defValue) : LabelVisitor(l, defValue) {};
1256: void visitPoint(const typename sieve_type::point_type& point) {};
1257: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1258: this->value = std::max(this->value, this->getLabelValue(arrow.target));
1259: };
1260: };
1261: class HeightVisitor {
1262: protected:
1263: const sieve_type& sieve;
1264: label_type& height;
1265: int maxHeight;
1266: std::set<typename sieve_type::point_type> modifiedPoints;
1267: public:
1268: HeightVisitor(const sieve_type& s, label_type& h) : sieve(s), height(h), maxHeight(-1) {};
1269: void visitPoint(const typename sieve_type::point_type& point) {
1270: MaxSupportVisitor v(height, -1);
1272: // Compute the max height of the points in the support of p, and add 1
1273: this->sieve.support(point, v);
1274: const int h0 = v.getLabelValue(point);
1275: const int h1 = v.getValue() + 1;
1277: if(h1 != h0) {
1278: v.setLabelValue(point, h1);
1279: if (h1 > this->maxHeight) this->maxHeight = h1;
1280: this->modifiedPoints.insert(point);
1281: }
1282: };
1283: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1284: this->visitPoint(arrow.source);
1285: };
1286: int getMaxHeight() const {return this->maxHeight;};
1287: bool isModified() const {return this->modifiedPoints.size() > 0;};
1288: const std::set<typename sieve_type::point_type>& getModifiedPoints() const {return this->modifiedPoints;};
1289: void clear() {this->modifiedPoints.clear();};
1290: };
1291: class DepthVisitor {
1292: public:
1293: typedef typename sieve_type::point_type point_type;
1294: protected:
1295: const sieve_type& sieve;
1296: label_type& depth;
1297: int maxDepth;
1298: const point_type limitPoint;
1299: std::set<point_type> modifiedPoints;
1300: public:
1301: DepthVisitor(const sieve_type& s, label_type& d) : sieve(s), depth(d), maxDepth(-1), limitPoint(sieve.getChart().max()+1) {};
1302: DepthVisitor(const sieve_type& s, const point_type& limit, label_type& d) : sieve(s), depth(d), maxDepth(-1), limitPoint(limit) {};
1303: void visitPoint(const point_type& point) {
1304: if (point >= this->limitPoint) return;
1305: MaxConeVisitor v(depth, -1);
1307: // Compute the max height of the points in the support of p, and add 1
1308: this->sieve.cone(point, v);
1309: const int d0 = v.getLabelValue(point);
1310: const int d1 = v.getValue() + 1;
1312: if(d1 != d0) {
1313: v.setLabelValue(point, d1);
1314: if (d1 > this->maxDepth) this->maxDepth = d1;
1315: this->modifiedPoints.insert(point);
1316: }
1317: };
1318: void visitArrow(const typename sieve_type::arrow_type& arrow) {
1319: this->visitPoint(arrow.target);
1320: };
1321: int getMaxDepth() const {return this->maxDepth;};
1322: bool isModified() const {return this->modifiedPoints.size() > 0;};
1323: const std::set<typename sieve_type::point_type>& getModifiedPoints() const {return this->modifiedPoints;};
1324: void clear() {this->modifiedPoints.clear();};
1325: };
1326: protected:
1327: Obj<sieve_type> _sieve;
1328: labels_type _labels;
1329: int _maxHeight;
1330: int _maxDepth;
1331: arrow_sections_type _arrowSections;
1332: real_sections_type _realSections;
1333: int_sections_type _intSections;
1334: Obj<oIndexArray> _indexArray;
1335: Obj<MeshNumberingFactory> _factory;
1336: bool _calculatedOverlap;
1337: Obj<send_overlap_type> _sendOverlap;
1338: Obj<recv_overlap_type> _recvOverlap;
1339: Obj<send_overlap_type> _distSendOverlap;
1340: Obj<recv_overlap_type> _distRecvOverlap;
1341: renumbering_type _renumbering;
1342: // Work space
1343: Obj<std::set<point_type> > _modifiedPoints;
1344: public:
1345: IBundle(MPI_Comm comm, int debug = 0) : ALE::ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1) {
1346: this->_indexArray = new oIndexArray();
1347: this->_modifiedPoints = new std::set<point_type>();
1348: this->_factory = MeshNumberingFactory::singleton(this->comm(), this->debug());
1349: this->_calculatedOverlap = false;
1350: this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
1351: this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
1352: };
1353: IBundle(const Obj<sieve_type>& sieve) : ALE::ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _maxHeight(-1), _maxDepth(-1) {
1354: this->_indexArray = new oIndexArray();
1355: this->_modifiedPoints = new std::set<point_type>();
1356: this->_factory = MeshNumberingFactory::singleton(this->comm(), this->debug());
1357: this->_calculatedOverlap = false;
1358: this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
1359: this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
1360: };
1361: virtual ~IBundle() {};
1362: public: // Verifiers
1363: bool hasLabel(const std::string& name) {
1364: if (this->_labels.find(name) != this->_labels.end()) {
1365: return true;
1366: }
1367: return false;
1368: };
1369: void checkLabel(const std::string& name) {
1370: if (!this->hasLabel(name)) {
1371: ostringstream msg;
1372: msg << "Invalid label name: " << name << std::endl;
1373: throw ALE::Exception(msg.str().c_str());
1374: }
1375: };
1376: public: // Accessors
1377: const Obj<sieve_type>& getSieve() const {return this->_sieve;};
1378: void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
1379: bool hasArrowSection(const std::string& name) const {
1380: return this->_arrowSections.find(name) != this->_arrowSections.end();
1381: };
1382: const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
1383: if (!this->hasArrowSection(name)) {
1384: Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());
1386: section->setName(name);
1387: if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
1388: this->_arrowSections[name] = section;
1389: }
1390: return this->_arrowSections[name];
1391: };
1392: void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
1393: this->_arrowSections[name] = section;
1394: };
1395: Obj<std::set<std::string> > getArrowSections() const {
1396: Obj<std::set<std::string> > names = std::set<std::string>();
1398: for(typename arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
1399: names->insert(s_iter->first);
1400: }
1401: return names;
1402: };
1403: bool hasRealSection(const std::string& name) const {
1404: return this->_realSections.find(name) != this->_realSections.end();
1405: };
1406: const Obj<real_section_type>& getRealSection(const std::string& name) {
1407: if (!this->hasRealSection(name)) {
1408: Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());
1410: section->setName(name);
1411: if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
1412: this->_realSections[name] = section;
1413: }
1414: return this->_realSections[name];
1415: };
1416: void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
1417: this->_realSections[name] = section;
1418: };
1419: Obj<std::set<std::string> > getRealSections() const {
1420: Obj<std::set<std::string> > names = std::set<std::string>();
1422: for(typename real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
1423: names->insert(s_iter->first);
1424: }
1425: return names;
1426: };
1427: bool hasIntSection(const std::string& name) const {
1428: return this->_intSections.find(name) != this->_intSections.end();
1429: };
1430: const Obj<int_section_type>& getIntSection(const std::string& name) {
1431: if (!this->hasIntSection(name)) {
1432: Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());
1434: section->setName(name);
1435: if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
1436: this->_intSections[name] = section;
1437: }
1438: return this->_intSections[name];
1439: };
1440: void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
1441: this->_intSections[name] = section;
1442: };
1443: Obj<std::set<std::string> > getIntSections() const {
1444: Obj<std::set<std::string> > names = std::set<std::string>();
1446: for(typename int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
1447: names->insert(s_iter->first);
1448: }
1449: return names;
1450: };
1451: const Obj<MeshNumberingFactory>& getFactory() const {return this->_factory;};
1452: renumbering_type& getRenumbering() {return this->_renumbering;};
1453: public: // Labels
1454: int getValue (const Obj<label_type>& label, const point_type& point, const int defValue = 0) {
1455: const Obj<typename label_type::coneSequence>& cone = label->cone(point);
1457: if (cone->size() == 0) return defValue;
1458: return *cone->begin();
1459: };
1460: void setValue(const Obj<label_type>& label, const point_type& point, const int value) {
1461: label->setCone(value, point);
1462: };
1463: void addValue(const Obj<label_type>& label, const point_type& point, const int value) {
1464: label->addCone(value, point);
1465: };
1466: template<typename InputPoints>
1467: int getMaxValue (const Obj<label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
1468: int maxValue = defValue;
1470: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1471: maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
1472: }
1473: return maxValue;
1474: }
1475: const Obj<label_type>& createLabel(const std::string& name) {
1476: this->_labels[name] = new label_type(this->comm(), this->debug());
1477: return this->_labels[name];
1478: };
1479: const Obj<label_type>& getLabel(const std::string& name) {
1480: this->checkLabel(name);
1481: return this->_labels[name];
1482: };
1483: void setLabel(const std::string& name, const Obj<label_type>& label) {
1484: this->_labels[name] = label;
1485: };
1486: const labels_type& getLabels() {
1487: return this->_labels;
1488: };
1489: virtual const Obj<label_sequence>& getLabelStratum(const std::string& name, int value) {
1490: this->checkLabel(name);
1491: return this->_labels[name]->support(value);
1492: };
1493: public: // Stratification
1494: void computeHeights() {
1495: const Obj<label_type>& label = this->createLabel(std::string("height"));
1496: HeightVisitor h(*this->_sieve, *label);
1498: #ifdef IMESH_NEW_LABELS
1499: label->setChart(this->getSieve()->getChart());
1500: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setConeSize(p, 1);}
1501: if (label->getChart().size()) {label->setSupportSize(0, label->getChart().size());}
1502: label->allocate();
1503: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setCone(-1, p);}
1504: #endif
1505: this->_sieve->leaves(h);
1506: while(h.isModified()) {
1507: // FIX: Avoid the copy here somehow by fixing the traversal
1508: std::vector<point_type> modifiedPoints(h.getModifiedPoints().begin(), h.getModifiedPoints().end());
1510: h.clear();
1511: this->_sieve->cone(modifiedPoints, h);
1512: }
1513: #ifdef IMESH_NEW_LABELS
1514: // Recalculate supportOffsets and populate support
1515: label->recalculateLabel();
1516: #endif
1517: this->_maxHeight = h.getMaxHeight();
1518: };
1519: virtual int height() const {return this->_maxHeight;};
1520: virtual void setHeight(const int height) {this->_maxHeight = height;};
1521: virtual int height(const point_type& point) {
1522: return this->getValue(this->_labels["height"], point, -1);
1523: };
1524: virtual void setHeight(const point_type& point, const int height) {
1525: return this->setValue(this->_labels["height"], point, height);
1526: };
1527: virtual const Obj<label_sequence>& heightStratum(int height) {
1528: return this->getLabelStratum("height", height);
1529: };
1530: void computeDepths() {
1531: const Obj<label_type>& label = this->createLabel(std::string("depth"));
1532: DepthVisitor d(*this->_sieve, *label);
1534: #ifdef IMESH_NEW_LABELS
1535: label->setChart(this->getSieve()->getChart());
1536: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setConeSize(p, 1);}
1537: if (label->getChart().size()) {label->setSupportSize(0, label->getChart().size());}
1538: label->allocate();
1539: for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setCone(-1, p);}
1540: #endif
1541: this->_sieve->roots(d);
1542: while(d.isModified()) {
1543: // FIX: Avoid the copy here somehow by fixing the traversal
1544: std::vector<point_type> modifiedPoints(d.getModifiedPoints().begin(), d.getModifiedPoints().end());
1546: d.clear();
1547: this->_sieve->support(modifiedPoints, d);
1548: }
1549: #ifdef IMESH_NEW_LABELS
1550: // Recalculate supportOffsets and populate support
1551: label->recalculateLabel();
1552: #endif
1553: this->_maxDepth = d.getMaxDepth();
1554: };
1555: virtual int depth() const {return this->_maxDepth;};
1556: virtual void setDepth(const int depth) {this->_maxDepth = depth;};
1557: virtual int depth(const point_type& point) {
1558: return this->getValue(this->_labels["depth"], point, -1);
1559: };
1560: virtual void setDepth(const point_type& point, const int depth) {
1561: return this->setValue(this->_labels["depth"], point, depth);
1562: };
1563: virtual const Obj<label_sequence>& depthStratum(int depth) {
1564: return this->getLabelStratum("depth", depth);
1565: };
1566: void stratify() {
1567: this->computeHeights();
1568: this->computeDepths();
1569: };
1570: protected:
1571: template<typename Value>
1572: static bool lt1(const Value& a, const Value& b) {
1573: return a.first < b.first;
1574: }
1575: public: // Allocation
1576: template<typename Section_>
1577: void reallocate(const Obj<Section_>& section) {
1578: if (!section->hasNewPoints()) return;
1579: 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()),
1580: std::max(std::max_element(section->getNewPoints().begin(), section->getNewPoints().end(), lt1<typename Section_::newpoint_type>)->first, section->getChart().max()-1)+1);
1581: section->reallocatePoint(newChart);
1582: }
1583: };
1584: #ifdef IMESH_NEW_LABELS
1585: template<typename Label_ = IFSieve<int> >
1586: #else
1587: template<typename IndexType, typename ScalarType, typename Label_ = LabelSifter<IndexType, IndexType> >
1588: #endif
1589: class IMesh : public IBundle<IFSieve<IndexType>, IGeneralSection<IndexType, ScalarType>, IGeneralSection<IndexType, IndexType>, Label_> {
1590: public:
1591: typedef IBundle<IFSieve<IndexType>, IGeneralSection<IndexType, ScalarType>, IGeneralSection<IndexType, IndexType>, Label_> base_type;
1592: typedef typename base_type::sieve_type sieve_type;
1593: typedef typename sieve_type::point_type point_type;
1594: typedef typename base_type::alloc_type alloc_type;
1595: typedef typename base_type::label_type label_type;
1596: typedef typename base_type::labels_type labels_type;
1597: typedef typename base_type::label_sequence label_sequence;
1598: typedef typename base_type::real_section_type real_section_type;
1599: typedef typename base_type::int_section_type int_section_type;
1600: typedef typename base_type::numbering_type numbering_type;
1601: typedef typename base_type::order_type order_type;
1602: typedef typename base_type::send_overlap_type send_overlap_type;
1603: typedef typename base_type::recv_overlap_type recv_overlap_type;
1604: typedef std::set<std::string> names_type;
1605: typedef std::vector<typename PETSc::Point<3> > holes_type;
1606: typedef std::map<std::string, Obj<Discretization> > discretizations_type;
1607: protected:
1608: int _dim;
1609: bool _calculatedOverlap;
1610: Obj<send_overlap_type> _sendOverlap;
1611: Obj<recv_overlap_type> _recvOverlap;
1612: std::map<int,double> _periodicity;
1613: holes_type _holes;
1614: discretizations_type _discretizations;
1615: int _maxDof;
1616: public:
1617: IMesh(MPI_Comm comm, int dim, int debug = 0) : base_type(comm, debug), _dim(dim) {
1618: this->_calculatedOverlap = false;
1619: this->_sendOverlap = new send_overlap_type(comm, debug);
1620: this->_recvOverlap = new recv_overlap_type(comm, debug);
1621: this->_maxDof = -1;
1622: };
1623: public: // Accessors
1624: int getDimension() const {return this->_dim;};
1625: void setDimension(const int dim) {this->_dim = dim;};
1626: bool getCalculatedOverlap() const {return this->_calculatedOverlap;};
1627: void setCalculatedOverlap(const bool calc) {this->_calculatedOverlap = calc;};
1628: const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
1629: void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
1630: const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
1631: void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
1632: bool getPeriodicity(const int d) {return this->_periodicity[d];};
1633: void setPeriodicity(const int d, const double length) {this->_periodicity[d] = length;};
1634: const holes_type& getHoles() const {return this->_holes;};
1635: void addHole(const double hole[]) {
1636: this->_holes.push_back(hole);
1637: };
1638: void copyHoles(const Obj<IMesh>& m) {
1639: const holes_type& holes = m->getHoles();
1641: for(holes_type::const_iterator h_iter = holes.begin(); h_iter != holes.end(); ++h_iter) {
1642: this->_holes.push_back(*h_iter);
1643: }
1644: };
1645: const Obj<Discretization>& getDiscretization() {return this->getDiscretization("default");};
1646: const Obj<Discretization>& getDiscretization(const std::string& name) {return this->_discretizations[name];};
1647: void setDiscretization(const Obj<Discretization>& disc) {this->setDiscretization("default", disc);};
1648: void setDiscretization(const std::string& name, const Obj<Discretization>& disc) {this->_discretizations[name] = disc;};
1649: Obj<names_type> getDiscretizations() const {
1650: Obj<names_type> names = names_type();
1652: for(discretizations_type::const_iterator d_iter = this->_discretizations.begin(); d_iter != this->_discretizations.end(); ++d_iter) {
1653: names->insert(d_iter->first);
1654: }
1655: return names;
1656: };
1657: int getMaxDof() const {return this->_maxDof;};
1658: void setMaxDof(const int maxDof) {this->_maxDof = maxDof;};
1659: public: // Sizes
1660: template<typename Section>
1661: int size(const Obj<Section>& section, const point_type& p) {
1662: typedef ISieveVisitor::SizeVisitor<sieve_type,Section> size_visitor_type;
1663: typedef ISieveVisitor::TransitiveClosureVisitor<sieve_type,size_visitor_type> closure_visitor_type;
1664: size_visitor_type sV(*section);
1665: closure_visitor_type cV(*this->getSieve(), sV);
1667: this->getSieve()->cone(p, cV);
1668: if (!sV.getSize()) sV.visitPoint(p);
1669: return sV.getSize();
1670: }
1671: template<typename Section>
1672: int sizeWithBC(const Obj<Section>& section, const point_type& p) {
1673: typedef ISieveVisitor::SizeWithBCVisitor<sieve_type,Section> size_visitor_type;
1674: typedef ISieveVisitor::TransitiveClosureVisitor<sieve_type,size_visitor_type> closure_visitor_type;
1675: size_visitor_type sV(*section);
1676: closure_visitor_type cV(*this->getSieve(), sV);
1678: this->getSieve()->cone(p, cV);
1679: if (!sV.getSize()) sV.visitPoint(p);
1680: return sV.getSize();
1681: }
1682: template<typename Section>
1683: void allocate(const Obj<Section>& section) {
1684: section->allocatePoint();
1685: }
1686: public: // Restrict/Update closures
1687: template<typename Sieve, typename Visitor>
1688: void closure1(const Sieve& sieve, const point_type& p, Visitor& v)
1689: {
1690: v.visitPoint(p, 0);
1691: // Cone is guarateed to be ordered correctly
1692: sieve.orientedCone(p, v);
1693: }
1694: // Return the values for the closure of this point
1695: template<typename Section>
1696: const typename Section::value_type *restrictClosure(const Obj<Section>& section, const point_type& p) {
1697: const int size = this->sizeWithBC(section, p);
1698: ISieveVisitor::RestrictVisitor<Section> rV(*section, size, section->getRawArray(size));
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 Section>
1710: const typename Section::value_type *restrictClosure(const Obj<Section>& section, const point_type& p, typename Section::value_type *values, const int valuesSize) {
1711: const int size = this->sizeWithBC(section, p);
1712: if (valuesSize < size) {throw ALE::Exception("Input array to small for restrictClosure()");}
1713: ISieveVisitor::RestrictVisitor<Section> rV(*section, size, values);
1715: if (this->depth() == 1) {
1716: closure1(*this->getSieve(), p, rV);
1717: } else {
1718: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::RestrictVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, rV, true);
1720: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1721: }
1722: return rV.getValues();
1723: }
1724: template<typename Visitor>
1725: void restrictClosure(const point_type& p, Visitor& v) {
1726: if (this->depth() == 1) {
1727: closure1(*this->getSieve(), p, v);
1728: } else {
1729: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, v);
1730: }
1731: }
1732: // Replace the values for the closure of this point
1733: template<typename Section>
1734: void update(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1735: ISieveVisitor::UpdateVisitor<Section> uV(*section, v);
1737: if (this->depth() == 1) {
1738: closure1(*this->getSieve(), p, uV);
1739: } else {
1740: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);
1742: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1743: }
1744: }
1745: // Replace the values for the closure of this point, including points constrained by BC
1746: template<typename Section>
1747: void updateAll(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1748: ISieveVisitor::UpdateAllVisitor<Section> uV(*section, v);
1750: if (this->depth() == 1) {
1751: closure1(*this->getSieve(), p, uV);
1752: } else {
1753: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateAllVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);
1755: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1756: }
1757: }
1758: // Augment the values for the closure of this point
1759: template<typename Section>
1760: void updateAdd(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1761: ISieveVisitor::UpdateAddVisitor<Section> uV(*section, v);
1763: if (this->depth() == 1) {
1764: closure1(*this->getSieve(), p, uV);
1765: } else {
1766: ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateAddVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);
1768: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1769: }
1770: }
1771: // Augment the values for the closure of this point
1772: template<typename Visitor>
1773: void updateClosure(const point_type& p, Visitor& v) {
1774: if (this->depth() == 1) {
1775: closure1(*this->getSieve(), p, v);
1776: } else {
1777: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, v);
1778: }
1779: }
1780: public: // Overlap
1781: void constructOverlap() {
1782: if (!this->_calculatedOverlap && (this->commSize() > 1)) {throw ALE::Exception("Must calculate overlap during distribution");}
1783: };
1784: public: // Cell topology and geometry
1785: int getNumCellCorners(const point_type& p, const int depth = -1) const {
1786: const int d = depth < 0 ? this->depth() : depth;
1788: if (d == 1) {
1789: return this->_sieve->getConeSize(p);
1790: } else if (d <= 0) {
1791: return 0;
1792: }
1793: // Warning: this is slow
1794: ISieveVisitor::NConeRetriever<sieve_type> ncV(*this->_sieve, (int) pow((double) this->_sieve->getMaxConeSize(), this->depth()));
1795: ALE::ISieveTraversal<sieve_type>::orientedClosure(*this->_sieve, p, ncV);
1796: return ncV.getOrientedSize();
1797: };
1798: int getNumCellCorners() {
1799: return getNumCellCorners(*this->heightStratum(0)->begin());
1800: };
1801: void setupCoordinates(const Obj<real_section_type>& coordinates) {
1802: const Obj<label_sequence>& vertices = this->depthStratum(0);
1804: coordinates->setChart(typename real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()),
1805: *std::max_element(vertices->begin(), vertices->end())+1));
1806: };
1807: point_type locatePoint(const typename real_section_type::value_type point[], point_type guess = -1) {
1808: //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
1809: if (guess != -1) {
1810: return guess;
1811: } else {
1812: throw ALE::Exception("No point location for mesh dimension");
1813: }
1814: };
1815: void computeTriangleGeometry(const Obj<real_section_type>& coordinates, const point_type& e, typename real_section_type::value_type v0[], typename real_section_type::value_type J[], typename real_section_type::value_type invJ[], typename real_section_type::value_type& detJ) {
1816: const PetscReal *coords = this->restrictClosure(coordinates, e);
1817: const int dim = 2;
1818: typename real_section_type::value_type invDet;
1820: if (v0) {
1821: for(int d = 0; d < dim; d++) {
1822: v0[d] = coords[d];
1823: }
1824: }
1825: if (J) {
1826: for(int d = 0; d < dim; d++) {
1827: for(int f = 0; f < dim; f++) {
1828: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1829: }
1830: }
1831: detJ = J[0]*J[3] - J[1]*J[2];
1832: if (detJ < 0.0) {
1833: const typename real_section_type::value_type xLength = this->_periodicity[0];
1835: if (xLength != 0.0) {
1836: typename real_section_type::value_type v0x = coords[0*dim+0];
1838: if (v0x == 0.0) {
1839: v0x = v0[0] = xLength;
1840: }
1841: for(int f = 0; f < dim; f++) {
1842: const typename real_section_type::value_type px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
1844: J[0*dim+f] = 0.5*(px - v0x);
1845: }
1846: }
1847: detJ = J[0]*J[3] - J[1]*J[2];
1848: }
1849: PetscLogFlopsNoError(8.0 + 3.0);
1850: }
1851: if (invJ) {
1852: invDet = 1.0/detJ;
1853: invJ[0] = invDet*J[3];
1854: invJ[1] = -invDet*J[1];
1855: invJ[2] = -invDet*J[2];
1856: invJ[3] = invDet*J[0];
1857: PetscLogFlopsNoError(5.0);
1858: }
1859: };
1860: void computeTetrahedronGeometry(const Obj<real_section_type>& coordinates, const point_type& e, typename real_section_type::value_type v0[], typename real_section_type::value_type J[], typename real_section_type::value_type invJ[], typename real_section_type::value_type& detJ) {
1861: const PetscReal *coords = this->restrictClosure(coordinates, e);
1862: const int dim = 3;
1863: typename real_section_type::value_type invDet;
1865: if (v0) {
1866: for(int d = 0; d < dim; d++) {
1867: v0[d] = coords[d];
1868: }
1869: }
1870: if (J) {
1871: for(int d = 0; d < dim; d++) {
1872: for(int f = 0; f < dim; f++) {
1873: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1874: }
1875: }
1876: // The minus sign is here since I orient the first face to get the outward normal
1877: detJ = -(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
1878: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
1879: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
1880: PetscLogFlopsNoError(18.0 + 12.0);
1881: }
1882: if (invJ) {
1883: invDet = -1.0/detJ;
1884: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
1885: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
1886: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
1887: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
1888: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
1889: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
1890: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1891: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
1892: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
1893: PetscLogFlopsNoError(37.0);
1894: }
1895: };
1896: void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, typename real_section_type::value_type v0[], typename real_section_type::value_type J[], typename real_section_type::value_type invJ[], typename real_section_type::value_type& detJ) {
1897: if (this->_dim == 2) {
1898: computeTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
1899: } else if (this->_dim == 3) {
1900: computeTetrahedronGeometry(coordinates, e, v0, J, invJ, detJ);
1901: } else {
1902: throw ALE::Exception("Unsupported dimension for element geometry computation");
1903: }
1904: };
1905: void computeBdSegmentGeometry(const Obj<real_section_type>& coordinates, const point_type& e, typename real_section_type::value_type v0[], typename real_section_type::value_type J[], typename real_section_type::value_type invJ[], typename real_section_type::value_type& detJ) {
1906: const typename real_section_type::value_type *coords = this->restrictClosure(coordinates, e);
1907: const int dim = 2;
1908: typename real_section_type::value_type invDet;
1910: if (v0) {
1911: for(int d = 0; d < dim; d++) {
1912: v0[d] = coords[d];
1913: }
1914: }
1915: if (J) {
1916: //r2 = coords[1*dim+0]*coords[1*dim+0] + coords[1*dim+1]*coords[1*dim+1];
1917: J[0] = (coords[1*dim+0] - coords[0*dim+0])*0.5; J[1] = (-coords[1*dim+1] + coords[0*dim+1])*0.5;
1918: J[2] = (coords[1*dim+1] - coords[0*dim+1])*0.5; J[3] = ( coords[1*dim+0] - coords[0*dim+0])*0.5;
1919: detJ = J[0]*J[3] - J[1]*J[2];
1920: if (detJ < 0.0) {
1921: const typename real_section_type::value_type xLength = this->_periodicity[0];
1923: if (xLength != 0.0) {
1924: typename real_section_type::value_type v0x = coords[0*dim+0];
1926: if (v0x == 0.0) {
1927: v0x = v0[0] = xLength;
1928: }
1929: for(int f = 0; f < dim; f++) {
1930: const typename real_section_type::value_type px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
1932: J[0*dim+f] = 0.5*(px - v0x);
1933: }
1934: }
1935: detJ = J[0]*J[3] - J[1]*J[2];
1936: }
1937: PetscLogFlopsNoError(8.0 + 3.0);
1938: }
1939: if (invJ) {
1940: invDet = 1.0/detJ;
1941: invJ[0] = invDet*J[3];
1942: invJ[1] = -invDet*J[1];
1943: invJ[2] = -invDet*J[2];
1944: invJ[3] = invDet*J[0];
1945: PetscLogFlopsNoError(5.0);
1946: }
1947: };
1948: void computeBdElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, typename real_section_type::value_type v0[], typename real_section_type::value_type J[], typename real_section_type::value_type invJ[], typename real_section_type::value_type& detJ) {
1949: if (this->_dim == 1) {
1950: computeBdSegmentGeometry(coordinates, e, v0, J, invJ, detJ);
1951: // } else if (this->_dim == 2) {
1952: // computeBdTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
1953: } else {
1954: throw ALE::Exception("Unsupported dimension for element geometry computation");
1955: }
1956: };
1957: typename real_section_type::value_type getMaxVolume() {
1958: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
1959: const Obj<label_sequence>& cells = this->heightStratum(0);
1960: const int dim = this->getDimension();
1961: typename real_section_type::value_type v0[3], J[9], invJ[9], detJ, refVolume = 0.0, maxVolume = 0.0;
1963: if (dim == 1) refVolume = 2.0;
1964: if (dim == 2) refVolume = 2.0;
1965: if (dim == 3) refVolume = 4.0/3.0;
1966: for(typename label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1967: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
1968: maxVolume = std::max(maxVolume, detJ*refVolume);
1969: }
1970: return maxVolume;
1971: };
1972: public:
1973: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1974: if (comm == MPI_COMM_NULL) {
1975: comm = this->comm();
1976: }
1977: if (name == "") {
1978: PetscPrintf(comm, "viewing a Mesh\n");
1979: } else {
1980: PetscPrintf(comm, "viewing Mesh '%s'\n", name.c_str());
1981: }
1982: this->getSieve()->view("mesh sieve", comm);
1983: Obj<names_type> sections = this->getRealSections();
1985: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1986: this->getRealSection(*name)->view(*name);
1987: }
1988: sections = this->getIntSections();
1989: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1990: this->getIntSection(*name)->view(*name);
1991: }
1992: sections = this->getArrowSections();
1993: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1994: this->getArrowSection(*name)->view(*name);
1995: }
1996: };
1997: public: // Discretization
1998: void markBoundaryCells(const std::string& name, const int marker = 1, const int newMarker = 2, const bool onlyVertices = false) {
1999: const Obj<label_type>& label = this->getLabel(name);
2000: const Obj<label_sequence>& boundary = this->getLabelStratum(name, marker);
2001: const typename label_sequence::iterator end = boundary->end();
2002: const Obj<sieve_type>& sieve = this->getSieve();
2004: if (!onlyVertices) {
2005: typename ISieveVisitor::MarkVisitor<sieve_type,label_type> mV(*label, newMarker);
2007: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != end; ++e_iter) {
2008: if (this->height(*e_iter) == 1) {
2009: sieve->support(*e_iter, mV);
2010: }
2011: }
2012: } else {
2013: #if 1
2014: throw ALE::Exception("Rewrite this to first mark boundary edges/faces.");
2015: #else
2016: const int depth = this->depth();
2018: for(typename label_sequence::iterator v_iter = boundary->begin(); v_iter != end; ++v_iter) {
2019: const Obj<supportArray> support = sieve->nSupport(*v_iter, depth);
2020: const typename supportArray::iterator sEnd = support->end();
2022: for(typename supportArray::iterator c_iter = support->begin(); c_iter != sEnd; ++c_iter) {
2023: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*c_iter);
2024: const typename sieve_type::traits::coneSequence::iterator cEnd = cone->end();
2026: for(typename sieve_type::traits::coneSequence::iterator e_iter = cone->begin(); e_iter != cEnd; ++e_iter) {
2027: if (sieve->support(*e_iter)->size() == 1) {
2028: this->setValue(label, *c_iter, newMarker);
2029: break;
2030: }
2031: }
2032: }
2033: }
2034: #endif
2035: }
2036: };
2037: int setFiberDimensions(const Obj<real_section_type>& s, const Obj<names_type>& discs, names_type& bcLabels) {
2038: const int debug = this->debug();
2039: int maxDof = 0;
2041: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
2042: s->addSpace();
2043: }
2044: for(int d = 0; d <= this->_dim; ++d) {
2045: int numDof = 0;
2046: int f = 0;
2048: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2049: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2050: const int sDof = disc->getNumDof(d);
2052: numDof += sDof;
2053: if (sDof) s->setFiberDimension(this->depthStratum(d), sDof, f);
2054: }
2055: if (numDof) s->setFiberDimension(this->depthStratum(d), numDof);
2056: maxDof = std::max(maxDof, numDof);
2057: }
2058: // Process exclusions
2059: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2060: int f = 0;
2062: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2063: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2064: std::string labelName = "exclude-"+*f_iter;
2065: std::set<point_type> seen;
2066: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth()), true);
2068: if (this->hasLabel(labelName)) {
2069: const Obj<label_type>& label = this->getLabel(labelName);
2070: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
2071: const typename label_sequence::iterator end = exclusion->end();
2072: if (debug > 1) {label->view(labelName.c_str());}
2074: for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2075: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *e_iter, pV);
2076: const typename Visitor::point_type *oPoints = pV.getPoints();
2077: const int oSize = pV.getSize();
2079: for(int cl = 0; cl < oSize; ++cl) {
2080: if (seen.find(oPoints[cl]) != seen.end()) continue;
2081: if (this->getValue(label, oPoints[cl]) == 1) {
2082: seen.insert(oPoints[cl]);
2083: s->setFiberDimension(oPoints[cl], 0, f);
2084: s->addFiberDimension(oPoints[cl], -disc->getNumDof(this->depth(oPoints[cl])));
2085: if (debug > 1) {std::cout << " point: " << oPoints[cl] << " dim: " << disc->getNumDof(this->depth(oPoints[cl])) << std::endl;}
2086: }
2087: }
2088: pV.clear();
2089: }
2090: }
2091: }
2092: // Process constraints
2093: f = 0;
2094: for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2095: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2096: const Obj<std::set<std::string> > bcs = disc->getBoundaryConditions();
2097: std::string excludeName = "exclude-"+*f_iter;
2099: for(std::set<std::string>::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter) {
2100: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
2101: const Obj<label_sequence>& boundary = this->getLabelStratum(bc->getLabelName(), bc->getMarker());
2103: bcLabels.insert(bc->getLabelName());
2104: if (this->hasLabel(excludeName)) {
2105: const Obj<label_type>& label = this->getLabel(excludeName);
2107: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2108: if (!this->getValue(label, *e_iter)) {
2109: const int numDof = disc->getNumDof(this->depth(*e_iter));
2111: if (numDof) s->addConstraintDimension(*e_iter, numDof);
2112: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2113: }
2114: }
2115: } else {
2116: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2117: const int numDof = disc->getNumDof(this->depth(*e_iter));
2119: if (numDof) s->addConstraintDimension(*e_iter, numDof);
2120: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2121: }
2122: }
2123: }
2124: }
2125: return maxDof;
2126: };
2127: void calculateIndices() {
2128: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2129: // Should have an iterator over the whole tree
2130: Obj<names_type> discs = this->getDiscretizations();
2131: const int debug = this->debug();
2132: std::map<std::string, std::pair<int, int*> > indices;
2134: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2135: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2137: indices[*d_iter] = std::pair<int, int*>(0, new int[disc->sizeV(*this)]);
2138: disc->setIndices(indices[*d_iter].second);
2139: }
2140: const Obj<label_sequence>& cells = this->heightStratum(0);
2141: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);
2142: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *cells->begin(), pV);
2143: const typename Visitor::point_type *oPoints = pV.getPoints();
2144: const int oSize = pV.getSize();
2145: int offset = 0;
2147: if (debug > 1) {std::cout << "Closure for first element" << std::endl;}
2148: for(int cl = 0; cl < oSize; ++cl) {
2149: const int dim = this->depth(oPoints[cl]);
2151: if (debug > 1) {std::cout << " point " << oPoints[cl] << " depth " << dim << std::endl;}
2152: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2153: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2154: const int num = disc->getNumDof(dim);
2156: if (debug > 1) {std::cout << " disc " << disc->getName() << " numDof " << num << std::endl;}
2157: for(int o = 0; o < num; ++o) {
2158: indices[*d_iter].second[indices[*d_iter].first++] = offset++;
2159: }
2160: }
2161: }
2162: pV.clear();
2163: if (debug > 1) {
2164: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2165: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2167: std::cout << "Discretization " << disc->getName() << " indices:";
2168: for(int i = 0; i < indices[*d_iter].first; ++i) {
2169: std::cout << " " << indices[*d_iter].second[i];
2170: }
2171: std::cout << std::endl;
2172: }
2173: }
2174: };
2175: void calculateIndicesExcluded(const Obj<real_section_type>& s, const Obj<names_type>& discs) {
2176: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2177: typedef std::map<std::string, std::pair<int, indexSet> > indices_type;
2178: const Obj<label_type>& indexLabel = this->createLabel("cellExclusion");
2179: const int debug = this->debug();
2180: int marker = 0;
2181: std::map<indices_type, int> indexMap;
2182: indices_type indices;
2183: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);
2185: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2186: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2187: const int size = disc->sizeV(*this);
2189: indices[*d_iter].second.resize(size);
2190: }
2191: const typename names_type::const_iterator dBegin = discs->begin();
2192: const typename names_type::const_iterator dEnd = discs->end();
2193: std::set<point_type> seen;
2194: int f = 0;
2196: for(typename names_type::const_iterator f_iter = dBegin; f_iter != dEnd; ++f_iter, ++f) {
2197: std::string labelName = "exclude-"+*f_iter;
2199: if (this->hasLabel(labelName)) {
2200: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
2201: const typename label_sequence::iterator end = exclusion->end();
2203: if (debug > 1) {std::cout << "Processing exclusion " << labelName << std::endl;}
2204: for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2205: if (this->height(*e_iter)) continue;
2206: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *e_iter, pV);
2207: const typename Visitor::point_type *oPoints = pV.getPoints();
2208: const int oSize = pV.getSize();
2209: int offset = 0;
2211: if (debug > 1) {std::cout << " Closure for cell " << *e_iter << std::endl;}
2212: for(int cl = 0; cl < oSize; ++cl) {
2213: int g = 0;
2215: if (debug > 1) {std::cout << " point " << oPoints[cl] << std::endl;}
2216: for(typename names_type::const_iterator g_iter = dBegin; g_iter != dEnd; ++g_iter, ++g) {
2217: const int fDim = s->getFiberDimension(oPoints[cl], g);
2219: if (debug > 1) {std::cout << " disc " << *g_iter << " numDof " << fDim << std::endl;}
2220: for(int d = 0; d < fDim; ++d) {
2221: indices[*g_iter].second[indices[*g_iter].first++] = offset++;
2222: }
2223: }
2224: }
2225: pV.clear();
2226: const std::map<indices_type, int>::iterator entry = indexMap.find(indices);
2228: if (debug > 1) {
2229: for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
2230: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2231: std::cout << "Discretization (" << i_iter->second << ") " << *g_iter << " indices:";
2232: for(int i = 0; i < ((indices_type) i_iter->first)[*g_iter].first; ++i) {
2233: std::cout << " " << ((indices_type) i_iter->first)[*g_iter].second[i];
2234: }
2235: std::cout << std::endl;
2236: }
2237: std::cout << "Comparison: " << (indices == i_iter->first) << std::endl;
2238: }
2239: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2240: std::cout << "Discretization " << *g_iter << " indices:";
2241: for(int i = 0; i < indices[*g_iter].first; ++i) {
2242: std::cout << " " << indices[*g_iter].second[i];
2243: }
2244: std::cout << std::endl;
2245: }
2246: }
2247: if (entry != indexMap.end()) {
2248: this->setValue(indexLabel, *e_iter, entry->second);
2249: if (debug > 1) {std::cout << " Found existing indices with marker " << entry->second << std::endl;}
2250: } else {
2251: indexMap[indices] = ++marker;
2252: this->setValue(indexLabel, *e_iter, marker);
2253: if (debug > 1) {std::cout << " Created new indices with marker " << marker << std::endl;}
2254: }
2255: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2256: indices[*g_iter].first = 0;
2257: for(unsigned int i = 0; i < indices[*g_iter].second.size(); ++i) indices[*g_iter].second[i] = 0;
2258: }
2259: }
2260: }
2261: }
2262: if (debug > 1) {indexLabel->view("cellExclusion");}
2263: for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
2264: if (debug > 1) {std::cout << "Setting indices for marker " << i_iter->second << std::endl;}
2265: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2266: const Obj<Discretization>& disc = this->getDiscretization(*g_iter);
2267: const indexSet indSet = ((indices_type) i_iter->first)[*g_iter].second;
2268: const int size = indSet.size();
2269: int *_indices = new int[size];
2271: if (debug > 1) {std::cout << " field " << *g_iter << std::endl;}
2272: for(int i = 0; i < size; ++i) {
2273: _indices[i] = indSet[i];
2274: if (debug > 1) {std::cout << " indices["<<i<<"] = " << _indices[i] << std::endl;}
2275: }
2276: disc->setIndices(_indices, i_iter->second);
2277: }
2278: }
2279: };
2280: void setupField(const Obj<real_section_type>& s, const int cellMarker = 2, const bool noUpdate = false) {
2281: typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2282: const Obj<names_type>& discs = this->getDiscretizations();
2283: const int debug = s->debug();
2284: names_type bcLabels;
2286: s->setChart(this->getSieve()->getChart());
2287: this->_maxDof = this->setFiberDimensions(s, discs, bcLabels);
2288: this->calculateIndices();
2289: this->calculateIndicesExcluded(s, discs);
2290: this->allocate(s);
2291: s->defaultConstraintDof();
2292: const Obj<label_type>& cellExclusion = this->getLabel("cellExclusion");
2294: if (debug > 1) {std::cout << "Setting boundary values" << std::endl;}
2295: for(typename names_type::const_iterator n_iter = bcLabels.begin(); n_iter != bcLabels.end(); ++n_iter) {
2296: const Obj<label_sequence>& boundaryCells = this->getLabelStratum(*n_iter, cellMarker);
2297: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2298: const Obj<names_type>& discs = this->getDiscretizations();
2299: const point_type firstCell = *boundaryCells->begin();
2300: const int numFields = discs->size();
2301: typename real_section_type::value_type *values = new typename real_section_type::value_type[this->sizeWithBC(s, firstCell)];
2302: int *dofs = new int[this->_maxDof];
2303: int *v = new int[numFields];
2304: typename real_section_type::value_type *v0 = new typename real_section_type::value_type[this->getDimension()];
2305: typename real_section_type::value_type *J = new typename real_section_type::value_type[this->getDimension()*this->getDimension()];
2306: typename real_section_type::value_type detJ;
2307: Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);
2309: for(typename label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
2310: ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *c_iter, pV);
2311: const typename Visitor::point_type *oPoints = pV.getPoints();
2312: const int oSize = pV.getSize();
2314: if (debug > 1) {std::cout << " Boundary cell " << *c_iter << std::endl;}
2315: this->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
2316: for(int f = 0; f < numFields; ++f) v[f] = 0;
2317: for(int cl = 0; cl < oSize; ++cl) {
2318: const int cDim = s->getConstraintDimension(oPoints[cl]);
2319: int off = 0;
2320: int f = 0;
2321: int i = -1;
2323: if (debug > 1) {std::cout << " point " << oPoints[cl] << std::endl;}
2324: if (cDim) {
2325: if (debug > 1) {std::cout << " constrained excMarker: " << this->getValue(cellExclusion, *c_iter) << std::endl;}
2326: for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2327: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2328: const Obj<names_type> bcs = disc->getBoundaryConditions();
2329: const int fDim = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
2330: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
2331: int b = 0;
2333: for(typename names_type::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter, ++b) {
2334: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
2335: const int value = this->getValue(this->getLabel(bc->getLabelName()), oPoints[cl]);
2337: if (b > 0) v[f] -= fDim;
2338: if (value == bc->getMarker()) {
2339: if (debug > 1) {std::cout << " field " << *f_iter << " marker " << value << std::endl;}
2340: for(int d = 0; d < fDim; ++d, ++v[f]) {
2341: dofs[++i] = off+d;
2342: if (!noUpdate) values[indices[v[f]]] = (*bc->getDualIntegrator())(v0, J, v[f], bc->getFunction());
2343: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2344: }
2345: // Allow only one condition per point
2346: ++b;
2347: break;
2348: } else {
2349: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
2350: for(int d = 0; d < fDim; ++d, ++v[f]) {
2351: values[indices[v[f]]] = 0.0;
2352: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2353: }
2354: }
2355: }
2356: if (b == 0) {
2357: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
2358: for(int d = 0; d < fDim; ++d, ++v[f]) {
2359: values[indices[v[f]]] = 0.0;
2360: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2361: }
2362: }
2363: off += fDim;
2364: }
2365: if (i != cDim-1) {throw ALE::Exception("Invalid constraint initialization");}
2366: s->setConstraintDof(oPoints[cl], dofs);
2367: } else {
2368: if (debug > 1) {std::cout << " unconstrained" << std::endl;}
2369: for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2370: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2371: const int fDim = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
2372: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
2374: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
2375: for(int d = 0; d < fDim; ++d, ++v[f]) {
2376: values[indices[v[f]]] = 0.0;
2377: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2378: }
2379: }
2380: }
2381: }
2382: if (debug > 1) {
2383: for(int f = 0; f < numFields; ++f) v[f] = 0;
2384: for(int cl = 0; cl < oSize; ++cl) {
2385: int f = 0;
2386: for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2387: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2388: const int fDim = s->getFiberDimension(oPoints[cl], f);
2389: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
2391: for(int d = 0; d < fDim; ++d, ++v[f]) {
2392: std::cout << " "<<*f_iter<<"-value["<<indices[v[f]]<<"] " << values[indices[v[f]]] << std::endl;
2393: }
2394: }
2395: }
2396: }
2397: if (!noUpdate) {
2398: this->updateAll(s, *c_iter, values);
2399: }
2400: pV.clear();
2401: }
2402: delete [] dofs;
2403: delete [] values;
2404: delete [] v0;
2405: delete [] J;
2406: }
2407: if (debug > 1) {s->view("");}
2408: };
2409: public:
2410: // Take in a map for the cells labels
2411: template<typename Section_>
2412: void relabel(Section_& labeling) {
2413: this->getSieve()->relabel(labeling);
2414: // Relabel sections
2415: Obj<std::set<std::string> > realNames = this->getRealSections();
2417: for(std::set<std::string>::const_iterator n_iter = realNames->begin(); n_iter != realNames->end(); ++n_iter) {
2418: Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());
2420: section->setName(*n_iter);
2421: ALE::Ordering<>::relabelSection(*this->getRealSection(*n_iter), labeling, *section);
2422: this->setRealSection(*n_iter, section);
2423: }
2424: Obj<std::set<std::string> > intNames = this->getIntSections();
2426: for(std::set<std::string>::const_iterator n_iter = intNames->begin(); n_iter != intNames->end(); ++n_iter) {
2427: Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());
2429: section->setName(*n_iter);
2430: ALE::Ordering<>::relabelSection(*this->getIntSection(*n_iter), labeling, *section);
2431: this->setIntSection(*n_iter, section);
2432: }
2433: // Relabel labels
2434: const labels_type& labels = this->getLabels();
2436: for(typename labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
2437: Obj<label_type> label = new label_type(this->comm(), this->debug());
2439: l_iter->second->relabel(labeling, *label);
2440: this->setLabel(l_iter->first, label);
2441: }
2442: // Relabel overlap
2443: Obj<send_overlap_type> sendOverlap = new send_overlap_type(this->comm(), this->debug());
2444: Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(this->comm(), this->debug());
2446: this->getSendOverlap()->relabel(labeling, *sendOverlap);
2447: this->setSendOverlap(sendOverlap);
2448: this->getRecvOverlap()->relabel(labeling, *recvOverlap);
2449: this->setRecvOverlap(recvOverlap);
2450: // Relabel distribution overlap ???
2451: // Relabel renumbering
2452: }
2453: };
2454: }
2456: namespace ALE {
2457: template<typename IndexType, typename ScalarType>
2458: class Mesh : public Bundle<ALE::Sieve<IndexType,IndexType,IndexType>, GeneralSection<IndexType,ScalarType> > {
2459: public:
2460: typedef Bundle<ALE::Sieve<IndexType,IndexType,IndexType>, GeneralSection<IndexType,ScalarType> > base_type;
2461: typedef typename base_type::sieve_type sieve_type;
2462: typedef typename sieve_type::point_type point_type;
2463: typedef typename base_type::alloc_type alloc_type;
2464: typedef typename base_type::label_type label_type;
2465: typedef typename base_type::label_sequence label_sequence;
2466: typedef typename base_type::coneArray coneArray;
2467: typedef typename base_type::supportArray supportArray;
2468: typedef typename base_type::arrow_section_type arrow_section_type;
2469: typedef typename base_type::real_section_type real_section_type;
2470: typedef typename base_type::numbering_type numbering_type;
2471: typedef typename base_type::order_type order_type;
2472: typedef typename base_type::send_overlap_type send_overlap_type;
2473: typedef typename base_type::recv_overlap_type recv_overlap_type;
2474: typedef typename base_type::sieve_alg_type sieve_alg_type;
2475: typedef std::set<std::string> names_type;
2476: typedef std::map<std::string, Obj<Discretization> > discretizations_type;
2477: typedef std::vector<PETSc::Point<3> > holes_type;
2478: protected:
2479: int _dim;
2480: discretizations_type _discretizations;
2481: std::map<int,double> _periodicity;
2482: holes_type _holes;
2483: public:
2484: Mesh(MPI_Comm comm, int dim, int debug = 0) : base_type(comm, debug), _dim(dim) {
2485: ///this->_factory = MeshNumberingFactory::singleton(debug);
2486: //std::cout << "["<<this->commRank()<<"]: Creating an ALE::Mesh" << std::endl;
2487: };
2488: ~Mesh() {
2489: //std::cout << "["<<this->commRank()<<"]: Destroying an ALE::Mesh" << std::endl;
2490: };
2491: public: // Accessors
2492: int getDimension() const {return this->_dim;};
2493: void setDimension(const int dim) {this->_dim = dim;};
2494: const Obj<Discretization>& getDiscretization() {return this->getDiscretization("default");};
2495: const Obj<Discretization>& getDiscretization(const std::string& name) {return this->_discretizations[name];};
2496: void setDiscretization(const Obj<Discretization>& disc) {this->setDiscretization("default", disc);};
2497: void setDiscretization(const std::string& name, const Obj<Discretization>& disc) {this->_discretizations[name] = disc;};
2498: Obj<names_type> getDiscretizations() const {
2499: Obj<names_type> names = names_type();
2501: for(discretizations_type::const_iterator d_iter = this->_discretizations.begin(); d_iter != this->_discretizations.end(); ++d_iter) {
2502: names->insert(d_iter->first);
2503: }
2504: return names;
2505: };
2506: bool getPeriodicity(const int d) {return this->_periodicity[d];};
2507: void setPeriodicity(const int d, const double length) {this->_periodicity[d] = length;};
2508: const holes_type& getHoles() const {return this->_holes;};
2509: void addHole(const double hole[]) {
2510: this->_holes.push_back(hole);
2511: };
2512: void copyHoles(const Obj<Mesh>& m) {
2513: const holes_type& holes = m->getHoles();
2515: for(holes_type::const_iterator h_iter = holes.begin(); h_iter != holes.end(); ++h_iter) {
2516: this->_holes.push_back(*h_iter);
2517: }
2518: };
2519: void copy(const Obj<Mesh>& m) {
2520: this->setSieve(m->getSieve());
2521: this->setLabel("height", m->getLabel("height"));
2522: this->_maxHeight = m->height();
2523: this->setLabel("depth", m->getLabel("depth"));
2524: this->_maxDepth = m->depth();
2525: this->setLabel("marker", m->getLabel("marker"));
2526: this->setRealSection("coordinates", m->getRealSection("coordinates"));
2527: this->setArrowSection("orientation", m->getArrowSection("orientation"));
2528: };
2529: public: // Cell topology and geometry
2530: int getNumCellCorners(const point_type& p, const int depth = -1) const {
2531: return (this->getDimension() > 0) ? this->_sieve->nCone(p, depth < 0 ? this->depth() : depth)->size() : 1;
2532: };
2533: int getNumCellCorners() {
2534: return getNumCellCorners(*(this->heightStratum(0)->begin()));
2535: };
2536: void setupCoordinates(const Obj<real_section_type>& coordinates) {};
2537: void computeTriangleGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2538: const double *coords = this->restrictClosure(coordinates, e);
2539: const int dim = 2;
2540: double invDet;
2542: if (v0) {
2543: for(int d = 0; d < dim; d++) {
2544: v0[d] = coords[d];
2545: }
2546: }
2547: if (J) {
2548: for(int d = 0; d < dim; d++) {
2549: for(int f = 0; f < dim; f++) {
2550: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
2551: }
2552: }
2553: detJ = J[0]*J[3] - J[1]*J[2];
2554: if (detJ < 0.0) {
2555: const double xLength = this->_periodicity[0];
2557: if (xLength != 0.0) {
2558: double v0x = coords[0*dim+0];
2560: if (v0x == 0.0) {
2561: v0x = v0[0] = xLength;
2562: }
2563: for(int f = 0; f < dim; f++) {
2564: const double px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
2566: J[0*dim+f] = 0.5*(px - v0x);
2567: }
2568: }
2569: detJ = J[0]*J[3] - J[1]*J[2];
2570: }
2571: PetscLogFlopsNoError(8.0 + 3.0);
2572: }
2573: if (invJ) {
2574: invDet = 1.0/detJ;
2575: invJ[0] = invDet*J[3];
2576: invJ[1] = -invDet*J[1];
2577: invJ[2] = -invDet*J[2];
2578: invJ[3] = invDet*J[0];
2579: PetscLogFlopsNoError(5.0);
2580: }
2581: };
2582: void computeQuadrilateralGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double point[], double v0[], double J[], double invJ[], double& detJ) {
2583: const double *coords = this->restrictClosure(coordinates, e);
2584: const int dim = 2;
2585: double invDet;
2587: if (v0) {
2588: for(int d = 0; d < dim; d++) {
2589: v0[d] = coords[d];
2590: }
2591: }
2592: if (J) {
2593: double x_1 = coords[2] - coords[0];
2594: double y_1 = coords[3] - coords[1];
2595: double x_2 = coords[6] - coords[0];
2596: double y_2 = coords[7] - coords[1];
2597: double x_3 = coords[4] - coords[0];
2598: double y_3 = coords[5] - coords[1];
2600: J[0] = x_1 + (x_3 - x_1 - x_2)*point[1];
2601: J[1] = x_2 + (x_3 - x_1 - x_2)*point[0];
2602: J[2] = y_1 + (y_3 - y_1 - y_2)*point[1];
2603: J[3] = y_1 + (y_3 - y_1 - y_2)*point[0];
2604: detJ = J[0]*J[3] - J[1]*J[2];
2605: PetscLogFlopsNoError(6.0 + 16.0 + 3.0);
2606: }
2607: if (invJ) {
2608: invDet = 1.0/detJ;
2609: invJ[0] = invDet*J[3];
2610: invJ[1] = -invDet*J[1];
2611: invJ[2] = -invDet*J[2];
2612: invJ[3] = invDet*J[0];
2613: PetscLogFlopsNoError(5.0);
2614: }
2615: };
2616: void computeTetrahedronGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2617: const double *coords = this->restrictClosure(coordinates, e);
2618: const int dim = 3;
2619: double invDet;
2621: if (v0) {
2622: for(int d = 0; d < dim; d++) {
2623: v0[d] = coords[d];
2624: }
2625: }
2626: if (J) {
2627: for(int d = 0; d < dim; d++) {
2628: for(int f = 0; f < dim; f++) {
2629: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
2630: }
2631: }
2632: // The minus sign is here since I orient the first face to get the outward normal
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(18.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 computeHexahedralGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double point[], double v0[], double J[], double invJ[], double& detJ) {
2653: const double *coords = this->restrictClosure(coordinates, e);
2654: const int dim = 3;
2655: double invDet;
2657: if (v0) {
2658: for(int d = 0; d < dim; d++) {
2659: v0[d] = coords[d];
2660: }
2661: }
2662: if (J) {
2663: double x_1 = coords[3] - coords[0];
2664: double y_1 = coords[4] - coords[1];
2665: double z_1 = coords[5] - coords[2];
2666: double x_2 = coords[6] - coords[0];
2667: double y_2 = coords[7] - coords[1];
2668: double z_2 = coords[8] - coords[2];
2669: double x_3 = coords[9] - coords[0];
2670: double y_3 = coords[10] - coords[1];
2671: double z_3 = coords[11] - coords[2];
2672: double x_4 = coords[12] - coords[0];
2673: double y_4 = coords[13] - coords[1];
2674: double z_4 = coords[14] - coords[2];
2675: double x_5 = coords[15] - coords[0];
2676: double y_5 = coords[16] - coords[1];
2677: double z_5 = coords[17] - coords[2];
2678: double x_6 = coords[18] - coords[0];
2679: double y_6 = coords[19] - coords[1];
2680: double z_6 = coords[20] - coords[2];
2681: double x_7 = coords[21] - coords[0];
2682: double y_7 = coords[22] - coords[1];
2683: double z_7 = coords[23] - coords[2];
2684: double g_x = x_1 - x_2 + x_3 + x_4 - x_5 + x_6 - x_7;
2685: double g_y = y_1 - y_2 + y_3 + y_4 - y_5 + y_6 - y_7;
2686: double g_z = z_1 - z_2 + z_3 + z_4 - z_5 + z_6 - z_7;
2688: 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];
2689: 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];
2690: 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];
2691: 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];
2692: 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];
2693: 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];
2694: 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];
2695: 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];
2696: 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];
2697: detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
2698: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
2699: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
2700: PetscLogFlopsNoError(39.0 + 81.0 + 12.0);
2701: }
2702: if (invJ) {
2703: invDet = 1.0/detJ;
2704: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
2705: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
2706: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
2707: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
2708: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
2709: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
2710: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
2711: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
2712: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
2713: PetscLogFlopsNoError(37.0);
2714: }
2715: };
2716: void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2717: if (this->_dim == 2) {
2718: computeTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
2719: } else if (this->_dim == 3) {
2720: computeTetrahedronGeometry(coordinates, e, v0, J, invJ, detJ);
2721: } else {
2722: throw ALE::Exception("Unsupported dimension for element geometry computation");
2723: }
2724: };
2725: void computeLineFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2726: const typename arrow_section_type::point_type arrow(cell, face);
2727: const bool reversed = (this->getArrowSection("orientation")->restrictPoint(arrow)[0] == -2);
2728: const int dim = this->getDimension();
2729: double norm = 0.0;
2730: double *vec = tangent;
2732: if (f == 0) {
2733: vec[0] = 0.0; vec[1] = -1.0;
2734: } else if (f == 1) {
2735: vec[0] = 0.70710678; vec[1] = 0.70710678;
2736: } else if (f == 2) {
2737: vec[0] = -1.0; vec[1] = 0.0;
2738: }
2739: for(int d = 0; d < dim; ++d) {
2740: normal[d] = 0.0;
2741: for(int e = 0; e < dim; ++e) normal[d] += cellInvJ[e*dim+d]*vec[e];
2742: if (reversed) normal[d] = -normal[d];
2743: norm += normal[d]*normal[d];
2744: }
2745: norm = std::sqrt(norm);
2746: for(int d = 0; d < dim; ++d) {
2747: normal[d] /= norm;
2748: }
2749: tangent[0] = normal[1];
2750: tangent[1] = -normal[0];
2751: if (this->debug()) {
2752: std::cout << "Cell: " << cell << " Face: " << face << "("<<f<<")" << std::endl;
2753: for(int d = 0; d < dim; ++d) {
2754: std::cout << "Normal["<<d<<"]: " << normal[d] << " Tangent["<<d<<"]: " << tangent[d] << std::endl;
2755: }
2756: }
2757: // Now get 1D Jacobian info
2758: // Should be a way to get this directly
2759: const double *coords = this->restrictClosure(this->getRealSection("coordinates"), face);
2760: detJ = std::sqrt(PetscSqr(coords[1*2+0] - coords[0*2+0]) + PetscSqr(coords[1*2+1] - coords[0*2+1]))/2.0;
2761: invJ[0] = 1.0/detJ;
2762: };
2763: void computeTriangleFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2764: const typename arrow_section_type::point_type arrow(cell, face);
2765: const bool reversed = this->getArrowSection("orientation")->restrictPoint(arrow)[0] < 0;
2766: const int dim = this->getDimension();
2767: const int faceDim = dim-1;
2768: double norm = 0.0;
2769: double *vec = tangent;
2771: if (f == 0) {
2772: vec[0] = 0.0; vec[1] = 0.0; vec[2] = -1.0;
2773: } else if (f == 1) {
2774: vec[0] = 0.0; vec[1] = -1.0; vec[2] = 0.0;
2775: } else if (f == 2) {
2776: vec[0] = 0.57735027; vec[1] = 0.57735027; vec[2] = 0.57735027;
2777: } else if (f == 3) {
2778: vec[0] = -1.0; vec[1] = 0.0; vec[2] = 0.0;
2779: }
2780: for(int d = 0; d < dim; ++d) {
2781: normal[d] = 0.0;
2782: for(int e = 0; e < dim; ++e) normal[d] += cellInvJ[e*dim+d]*vec[e];
2783: if (reversed) normal[d] = -normal[d];
2784: norm += normal[d]*normal[d];
2785: }
2786: norm = std::sqrt(norm);
2787: for(int d = 0; d < dim; ++d) {
2788: normal[d] /= norm;
2789: }
2790: // Get tangents
2791: tangent[0] = normal[1] - normal[2];
2792: tangent[1] = normal[2] - normal[0];
2793: tangent[2] = normal[0] - normal[1];
2794: norm = 0.0;
2795: for(int d = 0; d < dim; ++d) {
2796: norm += tangent[d]*tangent[d];
2797: }
2798: norm = std::sqrt(norm);
2799: for(int d = 0; d < dim; ++d) {
2800: tangent[d] /= norm;
2801: }
2802: tangent[3] = normal[1]*tangent[2] - normal[2]*tangent[1];
2803: tangent[4] = normal[2]*tangent[0] - normal[0]*tangent[2];
2804: tangent[5] = normal[0]*tangent[1] - normal[1]*tangent[0];
2805: if (this->debug()) {
2806: std::cout << "Cell: " << cell << " Face: " << face << "("<<f<<")" << std::endl;
2807: for(int d = 0; d < dim; ++d) {
2808: std::cout << "Normal["<<d<<"]: " << normal[d] << " TangentA["<<d<<"]: " << tangent[d] << " TangentB["<<d<<"]: " << tangent[dim+d] << std::endl;
2809: }
2810: }
2811: // Now get 2D Jacobian info
2812: // Should be a way to get this directly
2813: const double *coords = this->restrictClosure(this->getRealSection("coordinates"), face);
2814: // Rotate so that normal in z
2815: double invR[9], R[9];
2816: double detR, invDetR;
2817: for(int d = 0; d < dim; d++) {
2818: invR[d*dim+0] = tangent[d];
2819: invR[d*dim+1] = tangent[dim+d];
2820: invR[d*dim+2] = normal[d];
2821: }
2822: invDetR = (invR[0*3+0]*(invR[1*3+1]*invR[2*3+2] - invR[1*3+2]*invR[2*3+1]) +
2823: invR[0*3+1]*(invR[1*3+2]*invR[2*3+0] - invR[1*3+0]*invR[2*3+2]) +
2824: invR[0*3+2]*(invR[1*3+0]*invR[2*3+1] - invR[1*3+1]*invR[2*3+0]));
2825: detR = 1.0/invDetR;
2826: R[0*3+0] = detR*(invR[1*3+1]*invR[2*3+2] - invR[1*3+2]*invR[2*3+1]);
2827: R[0*3+1] = detR*(invR[0*3+2]*invR[2*3+1] - invR[0*3+1]*invR[2*3+2]);
2828: R[0*3+2] = detR*(invR[0*3+1]*invR[1*3+2] - invR[0*3+2]*invR[1*3+1]);
2829: R[1*3+0] = detR*(invR[1*3+2]*invR[2*3+0] - invR[1*3+0]*invR[2*3+2]);
2830: R[1*3+1] = detR*(invR[0*3+0]*invR[2*3+2] - invR[0*3+2]*invR[2*3+0]);
2831: R[1*3+2] = detR*(invR[0*3+2]*invR[1*3+0] - invR[0*3+0]*invR[1*3+2]);
2832: R[2*3+0] = detR*(invR[1*3+0]*invR[2*3+1] - invR[1*3+1]*invR[2*3+0]);
2833: R[2*3+1] = detR*(invR[0*3+1]*invR[2*3+0] - invR[0*3+0]*invR[2*3+1]);
2834: R[2*3+2] = detR*(invR[0*3+0]*invR[1*3+1] - invR[0*3+1]*invR[1*3+0]);
2835: for(int d = 0; d < dim; d++) {
2836: for(int e = 0; e < dim; e++) {
2837: invR[d*dim+e] = 0.0;
2838: for(int g = 0; g < dim; g++) {
2839: invR[d*dim+e] += R[e*dim+g]*coords[d*dim+g];
2840: }
2841: }
2842: }
2843: for(int d = dim-1; d >= 0; --d) {
2844: invR[d*dim+2] -= invR[0*dim+2];
2845: if (this->debug() && (d == dim-1)) {
2846: double ref[9];
2847: for(int q = 0; q < dim; q++) {
2848: for(int e = 0; e < dim; e++) {
2849: ref[q*dim+e] = 0.0;
2850: for(int g = 0; g < dim; g++) {
2851: ref[q*dim+e] += cellInvJ[e*dim+g]*coords[q*dim+g];
2852: }
2853: }
2854: }
2855: std::cout << "f: " << f << std::endl;
2856: std::cout << this->printMatrix(std::string("coords"), dim, dim, coords) << std::endl;
2857: std::cout << this->printMatrix(std::string("ref coords"), dim, dim, ref) << std::endl;
2858: std::cout << this->printMatrix(std::string("R"), dim, dim, R) << std::endl;
2859: std::cout << this->printMatrix(std::string("invR"), dim, dim, invR) << std::endl;
2860: }
2861: if (fabs(invR[d*dim+2]) > 1.0e-8) {
2862: throw ALE::Exception("Invalid rotation");
2863: }
2864: }
2865: double J[4];
2866: for(int d = 0; d < faceDim; d++) {
2867: for(int e = 0; e < faceDim; e++) {
2868: J[d*faceDim+e] = 0.5*(invR[(e+1)*dim+d] - invR[0*dim+d]);
2869: }
2870: }
2871: detJ = fabs(J[0]*J[3] - J[1]*J[2]);
2872: // Probably need something here if detJ < 0
2873: const double invDet = 1.0/detJ;
2874: invJ[0] = invDet*J[3];
2875: invJ[1] = -invDet*J[1];
2876: invJ[2] = -invDet*J[2];
2877: invJ[3] = invDet*J[0];
2878: };
2879: void computeFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2880: if (this->_dim == 2) {
2881: computeLineFaceGeometry(cell, face, f, cellInvJ, invJ, detJ, normal, tangent);
2882: } else if (this->_dim == 3) {
2883: computeTriangleFaceGeometry(cell, face, f, cellInvJ, invJ, detJ, normal, tangent);
2884: } else {
2885: throw ALE::Exception("Unsupported dimension for element geometry computation");
2886: }
2887: };
2888: double getMaxVolume() {
2889: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2890: const Obj<label_sequence>& cells = this->heightStratum(0);
2891: const int dim = this->getDimension();
2892: double v0[3], J[9], invJ[9], detJ, refVolume = 0.0, maxVolume = 0.0;
2894: if (dim == 1) refVolume = 2.0;
2895: if (dim == 2) refVolume = 2.0;
2896: if (dim == 3) refVolume = 4.0/3.0;
2897: for(typename label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2898: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2899: maxVolume = std::max(maxVolume, detJ*refVolume);
2900: }
2901: return maxVolume;
2902: };
2903: // Find the cell in which this point lies (stupid algorithm)
2904: point_type locatePoint_2D(const typename real_section_type::value_type point[]) {
2905: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2906: const Obj<label_sequence>& cells = this->heightStratum(0);
2907: const int embedDim = 2;
2908: double v0[2], J[4], invJ[4], detJ;
2910: for(typename label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2911: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2912: double xi = invJ[0*embedDim+0]*(point[0] - v0[0]) + invJ[0*embedDim+1]*(point[1] - v0[1]);
2913: double eta = invJ[1*embedDim+0]*(point[0] - v0[0]) + invJ[1*embedDim+1]*(point[1] - v0[1]);
2915: if ((xi >= 0.0) && (eta >= 0.0) && (xi + eta <= 2.0)) {
2916: return *c_iter;
2917: }
2918: }
2919: throw ALE::Exception("Could not locate point");
2920: };
2921: // Assume a simplex and 3D
2922: point_type locatePoint_3D(const typename real_section_type::value_type point[]) {
2923: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2924: const Obj<label_sequence>& cells = this->heightStratum(0);
2925: const int embedDim = 3;
2926: double v0[3], J[9], invJ[9], detJ;
2928: for(typename label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2929: this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2930: 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]);
2931: 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]);
2932: 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]);
2934: if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) {
2935: return *c_iter;
2936: }
2937: }
2938: throw ALE::Exception("Could not locate point");
2939: };
2940: point_type locatePoint(const typename real_section_type::value_type point[], point_type guess = -1) {
2941: //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
2942: if (guess != -1) {
2943: return guess;
2944: }else if (this->_dim == 2) {
2945: return locatePoint_2D(point);
2946: } else if (this->_dim == 3) {
2947: return locatePoint_3D(point);
2948: } else {
2949: throw ALE::Exception("No point location for mesh dimension");
2950: }
2951: };
2952: public: // Discretization
2953: void markBoundaryCells(const std::string& name, const int marker = 1, const int newMarker = 2, const bool onlyVertices = false) {
2954: const Obj<label_type>& label = this->getLabel(name);
2955: const Obj<label_sequence>& boundary = this->getLabelStratum(name, marker);
2956: const Obj<sieve_type>& sieve = this->getSieve();
2958: if (!onlyVertices) {
2959: const typename label_sequence::iterator end = boundary->end();
2961: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != end; ++e_iter) {
2962: if (this->height(*e_iter) == 1) {
2963: const point_type cell = *sieve->support(*e_iter)->begin();
2965: this->setValue(label, cell, newMarker);
2966: }
2967: }
2968: } else {
2969: const typename label_sequence::iterator end = boundary->end();
2970: const int depth = this->depth();
2972: for(typename label_sequence::iterator v_iter = boundary->begin(); v_iter != end; ++v_iter) {
2973: const Obj<supportArray> support = sieve->nSupport(*v_iter, depth);
2974: const typename supportArray::iterator sEnd = support->end();
2976: for(typename supportArray::iterator c_iter = support->begin(); c_iter != sEnd; ++c_iter) {
2977: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*c_iter);
2978: const typename sieve_type::traits::coneSequence::iterator cEnd = cone->end();
2980: for(typename sieve_type::traits::coneSequence::iterator e_iter = cone->begin(); e_iter != cEnd; ++e_iter) {
2981: if (sieve->support(*e_iter)->size() == 1) {
2982: this->setValue(label, *c_iter, newMarker);
2983: break;
2984: }
2985: }
2986: }
2987: }
2988: }
2989: };
2990: int setFiberDimensions(const Obj<real_section_type>& s, const Obj<names_type>& discs, names_type& bcLabels) {
2991: const int debug = this->debug();
2992: int maxDof = 0;
2994: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
2995: s->addSpace();
2996: }
2997: for(int d = 0; d <= this->_dim; ++d) {
2998: int numDof = 0;
2999: int f = 0;
3001: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3002: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3003: const int sDof = disc->getNumDof(d);
3005: numDof += sDof;
3006: if (sDof) s->setFiberDimension(this->depthStratum(d), sDof, f);
3007: }
3008: if (numDof) s->setFiberDimension(this->depthStratum(d), numDof);
3009: maxDof = std::max(maxDof, numDof);
3010: }
3011: // Process exclusions
3012: int f = 0;
3014: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3015: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3016: std::string labelName = "exclude-"+*f_iter;
3017: std::set<point_type> seen;
3019: if (this->hasLabel(labelName)) {
3020: const Obj<label_type>& label = this->getLabel(labelName);
3021: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
3022: const typename label_sequence::iterator end = exclusion->end();
3023: if (debug > 1) {label->view(labelName.c_str());}
3025: for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
3026: const Obj<coneArray> closure = ALE::SieveAlg<ALE::Mesh<IndexType,ScalarType> >::closure(this, this->getArrowSection("orientation"), *e_iter);
3027: const typename coneArray::iterator cEnd = closure->end();
3029: for(typename coneArray::iterator c_iter = closure->begin(); c_iter != cEnd; ++c_iter) {
3030: if (seen.find(*c_iter) != seen.end()) continue;
3031: if (this->getValue(label, *c_iter) == 1) {
3032: seen.insert(*c_iter);
3033: s->setFiberDimension(*c_iter, 0, f);
3034: s->addFiberDimension(*c_iter, -disc->getNumDof(this->depth(*c_iter)));
3035: if (debug > 1) {std::cout << " cell: " << *c_iter << " dim: " << disc->getNumDof(this->depth(*c_iter)) << std::endl;}
3036: }
3037: }
3038: }
3039: }
3040: }
3041: // Process constraints
3042: f = 0;
3043: for(typename std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3044: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3045: const Obj<std::set<std::string> > bcs = disc->getBoundaryConditions();
3046: std::string excludeName = "exclude-"+*f_iter;
3048: for(typename std::set<std::string>::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter) {
3049: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
3050: const Obj<label_sequence>& boundary = this->getLabelStratum(bc->getLabelName(), bc->getMarker());
3052: bcLabels.insert(bc->getLabelName());
3053: if (this->hasLabel(excludeName)) {
3054: const Obj<label_type>& label = this->getLabel(excludeName);
3056: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
3057: if (!this->getValue(label, *e_iter)) {
3058: const int numDof = disc->getNumDof(this->depth(*e_iter));
3060: if (numDof) s->addConstraintDimension(*e_iter, numDof);
3061: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
3062: }
3063: }
3064: } else {
3065: for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
3066: const int numDof = disc->getNumDof(this->depth(*e_iter));
3068: if (numDof) s->addConstraintDimension(*e_iter, numDof);
3069: if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
3070: }
3071: }
3072: }
3073: }
3074: return maxDof;
3075: };
3076: void calculateIndices() {
3077: // Should have an iterator over the whole tree
3078: Obj<names_type> discs = this->getDiscretizations();
3079: Obj<Mesh> mesh = this;
3080: const int debug = this->debug();
3081: std::map<std::string, std::pair<int, int*> > indices;
3083: mesh.addRef();
3084: for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3085: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3087: indices[*d_iter] = std::pair<int, int*>(0, new int[disc->size(mesh)]);
3088: disc->setIndices(indices[*d_iter].second);
3089: }
3090: const Obj<label_sequence>& cells = this->heightStratum(0);
3091: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *cells->begin());
3092: const typename coneArray::iterator end = closure->end();
3093: int offset = 0;
3095: if (debug > 1) {std::cout << "Closure for first element" << std::endl;}
3096: for(typename coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3097: const int dim = this->depth(*cl_iter);
3099: if (debug > 1) {std::cout << " point " << *cl_iter << " depth " << dim << std::endl;}
3100: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3101: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3102: const int num = disc->getNumDof(dim);
3104: if (debug > 1) {std::cout << " disc " << disc->getName() << " numDof " << num << std::endl;}
3105: for(int o = 0; o < num; ++o) {
3106: indices[*d_iter].second[indices[*d_iter].first++] = offset++;
3107: }
3108: }
3109: }
3110: if (debug > 1) {
3111: for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3112: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3114: std::cout << "Discretization " << disc->getName() << " indices:";
3115: for(int i = 0; i < indices[*d_iter].first; ++i) {
3116: std::cout << " " << indices[*d_iter].second[i];
3117: }
3118: std::cout << std::endl;
3119: }
3120: }
3121: };
3122: void calculateIndicesExcluded(const Obj<real_section_type>& s, const Obj<names_type>& discs) {
3123: typedef std::map<std::string, std::pair<int, indexSet> > indices_type;
3124: const Obj<label_type>& indexLabel = this->createLabel("cellExclusion");
3125: const int debug = this->debug();
3126: int marker = 0;
3127: Obj<Mesh> mesh = this;
3128: std::map<indices_type, int> indexMap;
3129: indices_type indices;
3131: mesh.addRef();
3132: for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3133: const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3134: const int size = disc->size(mesh);
3136: indices[*d_iter].second.resize(size);
3137: }
3138: const names_type::const_iterator dBegin = discs->begin();
3139: const names_type::const_iterator dEnd = discs->end();
3140: std::set<point_type> seen;
3141: int f = 0;
3143: for(names_type::const_iterator f_iter = dBegin; f_iter != dEnd; ++f_iter, ++f) {
3144: std::string labelName = "exclude-"+*f_iter;
3146: if (this->hasLabel(labelName)) {
3147: const Obj<label_sequence>& exclusion = this->getLabelStratum(labelName, 1);
3148: const typename label_sequence::iterator end = exclusion->end();
3150: if (debug > 1) {std::cout << "Processing exclusion " << labelName << std::endl;}
3151: for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
3152: if (this->height(*e_iter)) continue;
3153: const Obj<coneArray> closure = ALE::SieveAlg<ALE::Mesh<IndexType,ScalarType> >::closure(this, this->getArrowSection("orientation"), *e_iter);
3154: const typename coneArray::iterator clEnd = closure->end();
3155: int offset = 0;
3157: if (debug > 1) {std::cout << " Closure for cell " << *e_iter << std::endl;}
3158: for(typename coneArray::iterator cl_iter = closure->begin(); cl_iter != clEnd; ++cl_iter) {
3159: int g = 0;
3161: if (debug > 1) {std::cout << " point " << *cl_iter << std::endl;}
3162: for(typename names_type::const_iterator g_iter = dBegin; g_iter != dEnd; ++g_iter, ++g) {
3163: const int fDim = s->getFiberDimension(*cl_iter, g);
3165: if (debug > 1) {std::cout << " disc " << *g_iter << " numDof " << fDim << std::endl;}
3166: for(int d = 0; d < fDim; ++d) {
3167: indices[*g_iter].second[indices[*g_iter].first++] = offset++;
3168: }
3169: }
3170: }
3171: const typename std::map<indices_type, int>::iterator entry = indexMap.find(indices);
3173: if (debug > 1) {
3174: for(typename std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
3175: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3176: std::cout << "Discretization (" << i_iter->second << ") " << *g_iter << " indices:";
3177: for(int i = 0; i < ((indices_type) i_iter->first)[*g_iter].first; ++i) {
3178: std::cout << " " << ((indices_type) i_iter->first)[*g_iter].second[i];
3179: }
3180: std::cout << std::endl;
3181: }
3182: std::cout << "Comparison: " << (indices == i_iter->first) << std::endl;
3183: }
3184: for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3185: std::cout << "Discretization " << *g_iter << " indices:";
3186: for(int i = 0; i < indices[*g_iter].first; ++i) {
3187: std::cout << " " << indices[*g_iter].second[i];
3188: }
3189: std::cout << std::endl;
3190: }
3191: }
3192: if (entry != indexMap.end()) {
3193: this->setValue(indexLabel, *e_iter, entry->second);
3194: if (debug > 1) {std::cout << " Found existing indices with marker " << entry->second << std::endl;}
3195: } else {
3196: indexMap[indices] = ++marker;
3197: this->setValue(indexLabel, *e_iter, marker);
3198: if (debug > 1) {std::cout << " Created new indices with marker " << marker << std::endl;}
3199: }
3200: for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3201: indices[*g_iter].first = 0;
3202: for(unsigned int i = 0; i < indices[*g_iter].second.size(); ++i) indices[*g_iter].second[i] = 0;
3203: }
3204: }
3205: }
3206: }
3207: if (debug > 1) {indexLabel->view("cellExclusion");}
3208: for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
3209: if (debug > 1) {std::cout << "Setting indices for marker " << i_iter->second << std::endl;}
3210: for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3211: const Obj<Discretization>& disc = this->getDiscretization(*g_iter);
3212: const indexSet indSet = ((indices_type) i_iter->first)[*g_iter].second;
3213: const int size = indSet.size();
3214: int *_indices = new int[size];
3216: if (debug > 1) {std::cout << " field " << *g_iter << std::endl;}
3217: for(int i = 0; i < size; ++i) {
3218: _indices[i] = indSet[i];
3219: if (debug > 1) {std::cout << " indices["<<i<<"] = " << _indices[i] << std::endl;}
3220: }
3221: disc->setIndices(_indices, i_iter->second);
3222: }
3223: }
3224: };
3225: void setupField(const Obj<real_section_type>& s, const int cellMarker = 2, const bool noUpdate = false) {
3226: const Obj<names_type>& discs = this->getDiscretizations();
3227: const int debug = s->debug();
3228: names_type bcLabels;
3229: int maxDof;
3231: maxDof = this->setFiberDimensions(s, discs, bcLabels);
3232: this->calculateIndices();
3233: this->calculateIndicesExcluded(s, discs);
3234: this->allocate(s);
3235: s->defaultConstraintDof();
3236: const Obj<label_type>& cellExclusion = this->getLabel("cellExclusion");
3238: if (debug > 1) {std::cout << "Setting boundary values" << std::endl;}
3239: for(names_type::const_iterator n_iter = bcLabels.begin(); n_iter != bcLabels.end(); ++n_iter) {
3240: const Obj<label_sequence>& boundaryCells = this->getLabelStratum(*n_iter, cellMarker);
3241: const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
3242: const Obj<names_type>& discs = this->getDiscretizations();
3243: const point_type firstCell = *boundaryCells->begin();
3244: const int numFields = discs->size();
3245: typename real_section_type::value_type *values = new typename real_section_type::value_type[this->sizeWithBC(s, firstCell)];
3246: int *dofs = new int[maxDof];
3247: int *v = new int[numFields];
3248: typename real_section_type::value_type *v0 = new typename real_section_type::value_type[this->getDimension()];
3249: typename real_section_type::value_type *J = new typename real_section_type::value_type[this->getDimension()*this->getDimension()];
3250: typename real_section_type::value_type detJ;
3252: for(typename label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
3253: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
3254: const typename coneArray::iterator end = closure->end();
3256: if (debug > 1) {std::cout << " Boundary cell " << *c_iter << std::endl;}
3257: this->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
3258: for(int f = 0; f < numFields; ++f) v[f] = 0;
3259: for(typename coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3260: const int cDim = s->getConstraintDimension(*cl_iter);
3261: int off = 0;
3262: int f = 0;
3263: int i = -1;
3265: if (debug > 1) {std::cout << " point " << *cl_iter << std::endl;}
3266: if (cDim) {
3267: if (debug > 1) {std::cout << " constrained excMarker: " << this->getValue(cellExclusion, *c_iter) << std::endl;}
3268: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3269: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3270: const Obj<names_type> bcs = disc->getBoundaryConditions();
3271: const int fDim = s->getFiberDimension(*cl_iter, f);//disc->getNumDof(this->depth(*cl_iter));
3272: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
3273: int b = 0;
3275: for(names_type::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter, ++b) {
3276: const Obj<ALE::BoundaryCondition>& bc = disc->getBoundaryCondition(*bc_iter);
3277: const int value = this->getValue(this->getLabel(bc->getLabelName()), *cl_iter);
3279: if (b > 0) v[f] -= fDim;
3280: if (value == bc->getMarker()) {
3281: if (debug > 1) {std::cout << " field " << *f_iter << " marker " << value << std::endl;}
3282: for(int d = 0; d < fDim; ++d, ++v[f]) {
3283: dofs[++i] = off+d;
3284: if (!noUpdate) values[indices[v[f]]] = (*bc->getDualIntegrator())(v0, J, v[f], bc->getFunction());
3285: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3286: }
3287: // Allow only one condition per point
3288: ++b;
3289: break;
3290: } else {
3291: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
3292: for(int d = 0; d < fDim; ++d, ++v[f]) {
3293: values[indices[v[f]]] = 0.0;
3294: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3295: }
3296: }
3297: }
3298: if (b == 0) {
3299: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
3300: for(int d = 0; d < fDim; ++d, ++v[f]) {
3301: values[indices[v[f]]] = 0.0;
3302: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3303: }
3304: }
3305: off += fDim;
3306: }
3307: if (i != cDim-1) {throw ALE::Exception("Invalid constraint initialization");}
3308: s->setConstraintDof(*cl_iter, dofs);
3309: } else {
3310: if (debug > 1) {std::cout << " unconstrained" << std::endl;}
3311: for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3312: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3313: const int fDim = s->getFiberDimension(*cl_iter, f);//disc->getNumDof(this->depth(*cl_iter));
3314: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
3316: if (debug > 1) {std::cout << " field " << *f_iter << std::endl;}
3317: for(int d = 0; d < fDim; ++d, ++v[f]) {
3318: values[indices[v[f]]] = 0.0;
3319: if (debug > 1) {std::cout << " setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3320: }
3321: }
3322: }
3323: }
3324: if (debug > 1) {
3325: const Obj<coneArray> closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
3326: const typename coneArray::iterator end = closure->end();
3328: for(int f = 0; f < numFields; ++f) v[f] = 0;
3329: for(typename coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3330: int f = 0;
3331: for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3332: const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
3333: const int fDim = s->getFiberDimension(*cl_iter, f);
3334: const int *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
3336: for(int d = 0; d < fDim; ++d, ++v[f]) {
3337: std::cout << " "<<*f_iter<<"-value["<<indices[v[f]]<<"] " << values[indices[v[f]]] << std::endl;
3338: }
3339: }
3340: }
3341: }
3342: if (!noUpdate) {
3343: this->updateAll(s, *c_iter, values);
3344: }
3345: }
3346: delete [] dofs;
3347: delete [] values;
3348: delete [] v0;
3349: delete [] J;
3350: }
3351: if (debug > 1) {s->view("");}
3352: };
3353: public:
3354: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
3355: if (comm == MPI_COMM_NULL) {
3356: comm = this->comm();
3357: }
3358: if (name == "") {
3359: PetscPrintf(comm, "viewing a Mesh\n");
3360: } else {
3361: PetscPrintf(comm, "viewing Mesh '%s'\n", name.c_str());
3362: }
3363: this->getSieve()->view("mesh sieve", comm);
3364: Obj<names_type> sections = this->getRealSections();
3366: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3367: this->getRealSection(*name)->view(*name);
3368: }
3369: sections = this->getIntSections();
3370: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3371: this->getIntSection(*name)->view(*name);
3372: }
3373: sections = this->getArrowSections();
3374: for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3375: this->getArrowSection(*name)->view(*name);
3376: }
3377: };
3378: template<typename value_type>
3379: static std::string printMatrix(const std::string& name, const int rows, const int cols, const value_type matrix[], const int rank = -1)
3380: {
3381: ostringstream output;
3382: ostringstream rankStr;
3384: if (rank >= 0) {
3385: rankStr << "[" << rank << "]";
3386: }
3387: output << rankStr.str() << name << " = " << std::endl;
3388: for(int r = 0; r < rows; r++) {
3389: if (r == 0) {
3390: output << rankStr.str() << " /";
3391: } else if (r == rows-1) {
3392: output << rankStr.str() << " \\";
3393: } else {
3394: output << rankStr.str() << " |";
3395: }
3396: for(int c = 0; c < cols; c++) {
3397: output << " " << matrix[r*cols+c];
3398: }
3399: if (r == 0) {
3400: output << " \\" << std::endl;
3401: } else if (r == rows-1) {
3402: output << " /" << std::endl;
3403: } else {
3404: output << " |" << std::endl;
3405: }
3406: }
3407: return output.str();
3408: }
3409: };
3410: template<typename Mesh>
3411: class MeshBuilder {
3412: public:
3413: typedef typename Mesh::real_section_type::value_type real;
3414: public:
3417: /*
3418: Simple square boundary:
3420: 18--5-17--4--16
3421: | | |
3422: 6 10 3
3423: | | |
3424: 19-11-20--9--15
3425: | | |
3426: 7 8 2
3427: | | |
3428: 12--0-13--1--14
3429: */
3430: static Obj<Mesh> createSquareBoundary(const MPI_Comm comm, const real lower[], const real upper[], const int edges[], const int debug = 0) {
3431: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3432: int numVertices = (edges[0]+1)*(edges[1]+1);
3433: int numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
3434: real *coords = new real[numVertices*2];
3435: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3436: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3437: int order = 0;
3439: mesh->setSieve(sieve);
3440: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3441: if (mesh->commRank() == 0) {
3442: /* Create sieve and ordering */
3443: for(int v = numEdges; v < numEdges+numVertices; v++) {
3444: vertices[v-numEdges] = typename Mesh::point_type(v);
3445: }
3446: for(int vy = 0; vy <= edges[1]; vy++) {
3447: for(int ex = 0; ex < edges[0]; ex++) {
3448: typename Mesh::point_type edge(vy*edges[0] + ex);
3449: int vertex = vy*(edges[0]+1) + ex;
3451: sieve->addArrow(vertices[vertex+0], edge, order++);
3452: sieve->addArrow(vertices[vertex+1], edge, order++);
3453: if ((vy == 0) || (vy == edges[1])) {
3454: mesh->setValue(markers, edge, 1);
3455: mesh->setValue(markers, vertices[vertex], 1);
3456: if (ex == edges[0]-1) {
3457: mesh->setValue(markers, vertices[vertex+1], 1);
3458: }
3459: }
3460: }
3461: }
3462: for(int vx = 0; vx <= edges[0]; vx++) {
3463: for(int ey = 0; ey < edges[1]; ey++) {
3464: typename Mesh::point_type edge(vx*edges[1] + ey + edges[0]*(edges[1]+1));
3465: int vertex = ey*(edges[0]+1) + vx;
3467: sieve->addArrow(vertices[vertex], edge, order++);
3468: sieve->addArrow(vertices[vertex+edges[0]+1], edge, order++);
3469: if ((vx == 0) || (vx == edges[0])) {
3470: mesh->setValue(markers, edge, 1);
3471: mesh->setValue(markers, vertices[vertex], 1);
3472: if (ey == edges[1]-1) {
3473: mesh->setValue(markers, vertices[vertex+edges[0]+1], 1);
3474: }
3475: }
3476: }
3477: }
3478: }
3479: mesh->stratify();
3480: for(int vy = 0; vy <= edges[1]; ++vy) {
3481: for(int vx = 0; vx <= edges[0]; ++vx) {
3482: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
3483: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
3484: }
3485: }
3486: delete [] vertices;
3487: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3488: return mesh;
3489: };
3492: /*
3493: Simple square boundary:
3495: 14--5-13--4--12
3496: | |
3497: 6 3
3498: | |
3499: 15 11
3500: | |
3501: 7 2
3502: | |
3503: 8--0--9--1--10
3504: */
3505: static Obj<Mesh> createSquareBoundary(const MPI_Comm comm, const real lower[], const real upper[], const int edges, const int debug = 0) {
3506: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3507: int numVertices = edges*4;
3508: int numEdges = edges*4;
3509: real *coords = new real[numVertices*2];
3510: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3511: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3513: mesh->setSieve(sieve);
3514: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3515: if (mesh->commRank() == 0) {
3516: /* Create sieve and ordering */
3517: for(int v = numEdges; v < numEdges+numVertices; v++) {
3518: vertices[v-numEdges] = typename Mesh::point_type(v);
3519: }
3520: for(int e = 0; e < numEdges; ++e) {
3521: typename Mesh::point_type edge(e);
3522: int order = 0;
3524: sieve->addArrow(vertices[e], edge, order++);
3525: sieve->addArrow(vertices[(e+1)%numVertices], edge, order++);
3526: mesh->setValue(markers, edge, 2);
3527: mesh->setValue(markers, vertices[e], 1);
3528: mesh->setValue(markers, vertices[(e+1)%numVertices], 1);
3529: }
3530: }
3531: mesh->stratify();
3532: for(int v = 0; v < edges; ++v) {
3533: coords[(v+edges*0)*2+0] = lower[0] + ((upper[0] - lower[0])/edges)*v;
3534: coords[(v+edges*0)*2+1] = lower[1];
3535: }
3536: for(int v = 0; v < edges; ++v) {
3537: coords[(v+edges*1)*2+0] = upper[0];
3538: coords[(v+edges*1)*2+1] = lower[1] + ((upper[1] - lower[1])/edges)*v;
3539: }
3540: for(int v = 0; v < edges; ++v) {
3541: coords[(v+edges*2)*2+0] = upper[0] - ((upper[0] - lower[0])/edges)*v;
3542: coords[(v+edges*2)*2+1] = upper[1];
3543: }
3544: for(int v = 0; v < edges; ++v) {
3545: coords[(v+edges*3)*2+0] = lower[0];
3546: coords[(v+edges*3)*2+1] = upper[1] - ((upper[1] - lower[1])/edges)*v;
3547: }
3548: delete [] vertices;
3549: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3550: // Build normals for cells
3551: const Obj<typename Mesh::real_section_type>& normals = mesh->getRealSection("normals");
3552: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
3554: //normals->setChart(typename Mesh::real_section_type::chart_type(*std::min_element(cells->begin(), cells->end()),
3555: // *std::max_element(cells->begin(), cells->end())+1));
3556: normals->setFiberDimension(cells, mesh->getDimension()+1);
3557: mesh->allocate(normals);
3558: for(int e = 0; e < edges; ++e) {
3559: real normal[2] = {0.0, -1.0};
3560: normals->updatePoint(e+edges*0, normal);
3561: }
3562: for(int e = 0; e < edges; ++e) {
3563: real normal[2] = {1.0, 0.0};
3564: normals->updatePoint(e+edges*1, normal);
3565: }
3566: for(int e = 0; e < edges; ++e) {
3567: real normal[2] = {0.0, 1.0};
3568: normals->updatePoint(e+edges*2, normal);
3569: }
3570: for(int e = 0; e < edges; ++e) {
3571: real normal[2] = {-1.0, 0.0};
3572: normals->updatePoint(e+edges*3, normal);
3573: }
3574: return mesh;
3575: };
3578: /*
3579: Simple square boundary:
3581: 18--5-17--4--16
3582: | | |
3583: 6 10 3
3584: | | |
3585: 19-11-20--9--15
3586: | | |
3587: 7 8 2
3588: | | |
3589: 12--0-13--1--14
3590: */
3591: static Obj<Mesh> createParticleInSquareBoundary(const MPI_Comm comm, const real lower[], const real upper[], const int edges[], const real radius, const int partEdges, const int debug = 0) {
3592: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3593: const int numSquareVertices = (edges[0]+1)*(edges[1]+1);
3594: const int numVertices = numSquareVertices + partEdges;
3595: const int numSquareEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
3596: const int numEdges = numSquareEdges + partEdges;
3597: real *coords = new real[numVertices*2];
3598: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3599: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3600: int order = 0;
3602: mesh->setSieve(sieve);
3603: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3604: if (mesh->commRank() == 0) {
3605: /* Create sieve and ordering */
3606: for(int v = numEdges; v < numEdges+numVertices; v++) {
3607: vertices[v-numEdges] = typename Mesh::point_type(v);
3608: }
3609: // Make square
3610: for(int vy = 0; vy <= edges[1]; vy++) {
3611: for(int ex = 0; ex < edges[0]; ex++) {
3612: typename Mesh::point_type edge(vy*edges[0] + ex);
3613: int vertex = vy*(edges[0]+1) + ex;
3615: sieve->addArrow(vertices[vertex+0], edge, order++);
3616: sieve->addArrow(vertices[vertex+1], edge, order++);
3617: if ((vy == 0) || (vy == edges[1])) {
3618: mesh->setValue(markers, edge, 1);
3619: mesh->setValue(markers, vertices[vertex], 1);
3620: if (ex == edges[0]-1) {
3621: mesh->setValue(markers, vertices[vertex+1], 1);
3622: }
3623: }
3624: }
3625: }
3626: for(int vx = 0; vx <= edges[0]; vx++) {
3627: for(int ey = 0; ey < edges[1]; ey++) {
3628: typename Mesh::point_type edge(vx*edges[1] + ey + edges[0]*(edges[1]+1));
3629: int vertex = ey*(edges[0]+1) + vx;
3631: sieve->addArrow(vertices[vertex], edge, order++);
3632: sieve->addArrow(vertices[vertex+edges[0]+1], edge, order++);
3633: if ((vx == 0) || (vx == edges[0])) {
3634: mesh->setValue(markers, edge, 1);
3635: mesh->setValue(markers, vertices[vertex], 1);
3636: if (ey == edges[1]-1) {
3637: mesh->setValue(markers, vertices[vertex+edges[0]+1], 1);
3638: }
3639: }
3640: }
3641: }
3642: // Make particle
3643: for(int ep = 0; ep < partEdges; ++ep) {
3644: typename Mesh::point_type edge(numSquareEdges + ep);
3645: const int vertexA = numSquareVertices + ep;
3646: const int vertexB = numSquareVertices + (ep+1)%partEdges;
3648: sieve->addArrow(vertices[vertexA], edge, order++);
3649: sieve->addArrow(vertices[vertexB], edge, order++);
3650: mesh->setValue(markers, edge, 2);
3651: mesh->setValue(markers, vertices[vertexA], 2);
3652: mesh->setValue(markers, vertices[vertexB], 2);
3653: }
3654: }
3655: mesh->stratify();
3656: for(int vy = 0; vy <= edges[1]; ++vy) {
3657: for(int vx = 0; vx <= edges[0]; ++vx) {
3658: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
3659: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
3660: }
3661: }
3662: const real centroidX = 0.5*(upper[0] + lower[0]);
3663: const real centroidY = 0.5*(upper[1] + lower[1]);
3664: for(int vp = 0; vp < partEdges; ++vp) {
3665: const real rad = 2.0*PETSC_PI*vp/partEdges;
3666: coords[(numSquareVertices+vp)*2+0] = centroidX + radius*cos(rad);
3667: coords[(numSquareVertices+vp)*2+1] = centroidY + radius*sin(rad);
3668: }
3669: delete [] vertices;
3670: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3671: return mesh;
3672: };
3675: /*
3676: Simple boundary with reentrant singularity:
3678: 12--5-11
3679: | |
3680: | 4
3681: | |
3682: 6 10--3--9
3683: | |
3684: | 2
3685: | |
3686: 7-----1-----8
3687: */
3688: static Obj<Mesh> createReentrantBoundary(const MPI_Comm comm, const real lower[], const real upper[], real notchpercent[], const int debug = 0) {
3689: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3690: int numVertices = 6;
3691: int numEdges = numVertices;
3692: real *coords = new real[numVertices*2];
3693: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3694: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3696: mesh->setSieve(sieve);
3697: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3698: if (mesh->commRank() == 0) {
3699: /* Create sieve and ordering */
3700: for (int b = 0; b < numVertices; b++) {
3701: sieve->addArrow(numEdges+b, b);
3702: sieve->addArrow(numEdges+b, (b+1)%numVertices);
3703: mesh->setValue(markers, b, 1);
3704: mesh->setValue(markers, b+numVertices, 1);
3705: }
3706: coords[0] = upper[0];
3707: coords[1] = lower[1];
3709: coords[2] = lower[0];
3710: coords[3] = lower[1];
3711:
3712: coords[4] = lower[0];
3713: coords[5] = notchpercent[1]*lower[1] + (1 - notchpercent[1])*upper[1];
3714:
3715: coords[6] = notchpercent[0]*upper[0] + (1 - notchpercent[0])*lower[0];
3716: coords[7] = notchpercent[1]*lower[1] + (1 - notchpercent[1])*upper[1];
3717:
3718:
3719: coords[8] = notchpercent[0]*upper[0] + (1 - notchpercent[0])*lower[0];
3720: coords[9] = upper[1];
3722: coords[10] = upper[0];
3723: coords[11] = upper[1];
3724: mesh->stratify();
3725: }
3726: delete [] vertices;
3727: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3728: return mesh;
3729: }
3733: /*
3734: Circular boundary with reentrant singularity:
3736: ---1
3737: -- |
3738: - |
3739: | |
3740: | 0-----n
3741: | |
3742: - -
3743: -- --
3744: -------
3745: */
3746: static Obj<Mesh> createCircularReentrantBoundary(const MPI_Comm comm, const int segments, const real radius, const real arc_percent, const int debug = 0) {
3747: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3748: int numVertices = segments+2;
3749: int numEdges = numVertices;
3750: real *coords = new real[numVertices*2];
3751: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3752: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3754: mesh->setSieve(sieve);
3755: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3756: if (mesh->commRank() == 0) {
3757: /* Create sieve and ordering */
3759: int startvertex = 1;
3760: if (arc_percent < 1.) {
3761: coords[0] = 0.;
3762: coords[1] = 0.;
3763: } else {
3764: numVertices = segments;
3765: numEdges = numVertices;
3766: startvertex = 0;
3767: }
3769: for (int b = 0; b < numVertices; b++) {
3770: sieve->addArrow(numEdges+b, b);
3771: sieve->addArrow(numEdges+b, (b+1)%numVertices);
3772: mesh->setValue(markers, b, 1);
3773: mesh->setValue(markers, b+numVertices, 1);
3774: }
3776: real anglestep = arc_percent*2.*3.14159265/((float)segments);
3778: for (int i = startvertex; i < numVertices; i++) {
3779: coords[2*i] = radius * sin(anglestep*(i-startvertex));
3780: coords[2*i+1] = radius*cos(anglestep*(i-startvertex));
3781: }
3782: mesh->stratify();
3783: }
3784: delete [] vertices;
3785: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3786: return mesh;
3787: };
3790: static Obj<Mesh> createAnnularBoundary(const MPI_Comm comm, const int segments, const real centers[4], const real radii[2], const int debug = 0) {
3791: Obj<Mesh> mesh = new Mesh(comm, 1, debug);
3792: int numVertices = segments*2;
3793: int numEdges = numVertices;
3794: real *coords = new real[numVertices*2];
3795: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3796: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3798: mesh->setSieve(sieve);
3799: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3800: if (mesh->commRank() == 0) {
3801: for (int e = 0; e < segments; ++e) {
3802: sieve->addArrow(numEdges+e, e);
3803: sieve->addArrow(numEdges+(e+1)%segments, e);
3804: sieve->addArrow(numEdges+segments+e, e+segments);
3805: sieve->addArrow(numEdges+segments+(e+1)%segments, e+segments);
3806: mesh->setValue(markers, e, 1);
3807: mesh->setValue(markers, e+segments, 1);
3808: mesh->setValue(markers, e+numEdges, 1);
3809: mesh->setValue(markers, e+numEdges+segments, 1);
3810: }
3811: const real anglestep = 2.0*M_PI/segments;
3813: for (int v = 0; v < segments; ++v) {
3814: coords[v*2] = centers[0] + radii[0]*cos(anglestep*v);
3815: coords[v*2+1] = centers[1] + radii[0]*sin(anglestep*v);
3816: coords[(v+segments)*2] = centers[2] + radii[1]*cos(anglestep*v);
3817: coords[(v+segments)*2+1] = centers[3] + radii[1]*sin(anglestep*v);
3818: }
3819: mesh->addHole(¢ers[2]);
3820: }
3821: mesh->stratify();
3822: delete [] vertices;
3823: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3824: return mesh;
3825: };
3828: /*
3829: Simple cubic boundary:
3831: 30----31-----32
3832: | | |
3833: | 3 | 2 |
3834: | | |
3835: 27----28-----29
3836: | | |
3837: | 0 | 1 |
3838: | | |
3839: 24----25-----26
3840: */
3841: static Obj<Mesh> createCubeBoundary(const MPI_Comm comm, const real lower[], const real upper[], const int faces[], const int debug = 0) {
3842: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
3843: int numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
3844: int numFaces = 6;
3845: real *coords = new real[numVertices*3];
3846: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3847: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3848: int order = 0;
3850: mesh->setSieve(sieve);
3851: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3852: if (mesh->commRank() == 0) {
3853: /* Create sieve and ordering */
3854: for(int v = numFaces; v < numFaces+numVertices; v++) {
3855: vertices[v-numFaces] = typename Mesh::point_type(v);
3856: mesh->setValue(markers, vertices[v-numFaces], 1);
3857: }
3858: {
3859: // Side 0 (Front)
3860: typename Mesh::point_type face(0);
3861: sieve->addArrow(vertices[0], face, order++);
3862: sieve->addArrow(vertices[1], face, order++);
3863: sieve->addArrow(vertices[2], face, order++);
3864: sieve->addArrow(vertices[3], face, order++);
3865: mesh->setValue(markers, face, 1);
3866: }
3867: {
3868: // Side 1 (Back)
3869: typename Mesh::point_type face(1);
3870: sieve->addArrow(vertices[5], face, order++);
3871: sieve->addArrow(vertices[4], face, order++);
3872: sieve->addArrow(vertices[7], face, order++);
3873: sieve->addArrow(vertices[6], face, order++);
3874: mesh->setValue(markers, face, 1);
3875: }
3876: {
3877: // Side 2 (Bottom)
3878: typename Mesh::point_type face(2);
3879: sieve->addArrow(vertices[4], face, order++);
3880: sieve->addArrow(vertices[5], face, order++);
3881: sieve->addArrow(vertices[1], face, order++);
3882: sieve->addArrow(vertices[0], face, order++);
3883: mesh->setValue(markers, face, 1);
3884: }
3885: {
3886: // Side 3 (Top)
3887: typename Mesh::point_type face(3);
3888: sieve->addArrow(vertices[3], face, order++);
3889: sieve->addArrow(vertices[2], face, order++);
3890: sieve->addArrow(vertices[6], face, order++);
3891: sieve->addArrow(vertices[7], face, order++);
3892: mesh->setValue(markers, face, 1);
3893: }
3894: {
3895: // Side 4 (Left)
3896: typename Mesh::point_type face(4);
3897: sieve->addArrow(vertices[4], face, order++);
3898: sieve->addArrow(vertices[0], face, order++);
3899: sieve->addArrow(vertices[3], face, order++);
3900: sieve->addArrow(vertices[7], face, order++);
3901: mesh->setValue(markers, face, 1);
3902: }
3903: {
3904: // Side 5 (Right)
3905: typename Mesh::point_type face(5);
3906: sieve->addArrow(vertices[1], face, order++);
3907: sieve->addArrow(vertices[5], face, order++);
3908: sieve->addArrow(vertices[6], face, order++);
3909: sieve->addArrow(vertices[2], face, order++);
3910: mesh->setValue(markers, face, 1);
3911: }
3912: }
3913: mesh->stratify();
3914: #if 0
3915: for(int vz = 0; vz <= edges[2]; ++vz) {
3916: for(int vy = 0; vy <= edges[1]; ++vy) {
3917: for(int vx = 0; vx <= edges[0]; ++vx) {
3918: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
3919: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
3920: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
3921: }
3922: }
3923: }
3924: #else
3925: coords[0*3+0] = lower[0];
3926: coords[0*3+1] = lower[1];
3927: coords[0*3+2] = upper[2];
3928: coords[1*3+0] = upper[0];
3929: coords[1*3+1] = lower[1];
3930: coords[1*3+2] = upper[2];
3931: coords[2*3+0] = upper[0];
3932: coords[2*3+1] = upper[1];
3933: coords[2*3+2] = upper[2];
3934: coords[3*3+0] = lower[0];
3935: coords[3*3+1] = upper[1];
3936: coords[3*3+2] = upper[2];
3937: coords[4*3+0] = lower[0];
3938: coords[4*3+1] = lower[1];
3939: coords[4*3+2] = lower[2];
3940: coords[5*3+0] = upper[0];
3941: coords[5*3+1] = lower[1];
3942: coords[5*3+2] = lower[2];
3943: coords[6*3+0] = upper[0];
3944: coords[6*3+1] = upper[1];
3945: coords[6*3+2] = lower[2];
3946: coords[7*3+0] = lower[0];
3947: coords[7*3+1] = upper[1];
3948: coords[7*3+2] = lower[2];
3949: #endif
3950: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3951: return mesh;
3952: };
3954: // Creates a triangular prism boundary
3955: static Obj<Mesh> createPrismBoundary(const MPI_Comm comm, const real lower[], const real upper[], const int debug = 0) {
3956: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
3957: int numVertices = 6;
3958: int numFaces = 5;
3959: real *coords = new real[numVertices*3];
3960: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3961: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
3962: int order = 0;
3964: mesh->setSieve(sieve);
3965: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3966: if (mesh->commRank() == 0) {
3967: /* Create sieve and ordering */
3968: for(int v = numFaces; v < numFaces+numVertices; v++) {
3969: vertices[v-numFaces] = typename Mesh::point_type(v);
3970: mesh->setValue(markers, vertices[v-numFaces], 1);
3971: }
3972: {
3973: // Side 0 (Top)
3974: typename Mesh::point_type face(0);
3975: sieve->addArrow(vertices[0], face, order++);
3976: sieve->addArrow(vertices[1], face, order++);
3977: sieve->addArrow(vertices[2], face, order++);
3978: mesh->setValue(markers, face, 1);
3979: }
3980: {
3981: // Side 1 (Bottom)
3982: typename Mesh::point_type face(1);
3983: sieve->addArrow(vertices[3], face, order++);
3984: sieve->addArrow(vertices[5], face, order++);
3985: sieve->addArrow(vertices[4], face, order++);
3986: mesh->setValue(markers, face, 1);
3987: }
3988: {
3989: // Side 2 (Front)
3990: typename Mesh::point_type face(2);
3991: sieve->addArrow(vertices[0], face, order++);
3992: sieve->addArrow(vertices[3], face, order++);
3993: sieve->addArrow(vertices[4], face, order++);
3994: sieve->addArrow(vertices[1], face, order++);
3995: mesh->setValue(markers, face, 1);
3996: }
3997: {
3998: // Side 3 (Left)
3999: typename Mesh::point_type face(3);
4000: sieve->addArrow(vertices[1], face, order++);
4001: sieve->addArrow(vertices[4], face, order++);
4002: sieve->addArrow(vertices[5], face, order++);
4003: sieve->addArrow(vertices[2], face, order++);
4004: mesh->setValue(markers, face, 1);
4005: }
4006: {
4007: // Side 4 (Right)
4008: typename Mesh::point_type face(4);
4009: sieve->addArrow(vertices[2], face, order++);
4010: sieve->addArrow(vertices[5], face, order++);
4011: sieve->addArrow(vertices[3], face, order++);
4012: sieve->addArrow(vertices[0], face, order++);
4013: mesh->setValue(markers, face, 1);
4014: }
4015: }
4016: mesh->stratify();
4017: coords[0*3+0] = lower[0];
4018: coords[0*3+1] = lower[1];
4019: coords[0*3+2] = upper[2];
4020: coords[1*3+0] = upper[0];
4021: coords[1*3+1] = lower[1];
4022: coords[1*3+2] = upper[2];
4023: coords[2*3+0] = upper[0];
4024: coords[2*3+1] = upper[1];
4025: coords[2*3+2] = upper[2];
4026: coords[3*3+0] = lower[0];
4027: coords[3*3+1] = upper[1];
4028: coords[3*3+2] = upper[2];
4029: coords[4*3+0] = lower[0];
4030: coords[4*3+1] = lower[1];
4031: coords[4*3+2] = lower[2];
4032: coords[5*3+0] = upper[0];
4033: coords[5*3+1] = lower[1];
4034: coords[5*3+2] = lower[2];
4035: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
4036: return mesh;
4037: };
4041: /* v0
4042: / \
4043: / \
4044: 2 / 4 \ 1
4045: / \
4046: / v12 \
4047: v6|\ /|\ /|v5
4048: | \v8 | v7/ | z
4049: | |7 |8 | | |
4050: | |v13\ | | <-v4 / \
4051: | v9/ 9 \v10| x y
4052: v1 | 5 \ / 6 |v2
4053: \ \ / /
4054: \ v11 /
4055: \ | /
4056: 3 \ | /
4057: \|/
4058: v3
4059: */
4060: static Obj<Mesh> createFicheraCornerBoundary(const MPI_Comm comm, const real lower[], const real upper[], const real offset[], const int debug = 0) {
4061: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
4062: const int nVertices = 14;
4063: const int nFaces = 12;
4064: real ilower[3];
4065: ilower[0] = lower[0]*(1. - offset[0]) + upper[0]*offset[0];
4066: ilower[1] = lower[1]*(1. - offset[1]) + upper[1]*offset[1];
4067: ilower[2] = lower[2]*(1. - offset[2]) + upper[2]*offset[2];
4068: real coords[nVertices*3];
4069: //outer square-triplet
4070: coords[0*3+0] = lower[0];
4071: coords[0*3+1] = lower[1];
4072: coords[0*3+2] = upper[2];
4073: coords[1*3+0] = upper[0];
4074: coords[1*3+1] = lower[1];
4075: coords[1*3+2] = lower[2];
4076: coords[2*3+0] = lower[0];
4077: coords[2*3+1] = upper[1];
4078: coords[2*3+2] = lower[2];
4079: coords[3*3+0] = upper[0];
4080: coords[3*3+1] = upper[1];
4081: coords[3*3+2] = lower[2];
4082: coords[4*3+0] = lower[0];
4083: coords[4*3+1] = lower[1];
4084: coords[4*3+2] = lower[2];
4085: coords[5*3+0] = lower[0];
4086: coords[5*3+1] = upper[1];
4087: coords[5*3+2] = upper[2];
4088: coords[6*3+0] = upper[0];
4089: coords[6*3+1] = lower[1];
4090: coords[6*3+2] = upper[2];
4092: //inner square-triplet
4093: coords[7*3+0] = ilower[0];
4094: coords[7*3+1] = upper[1];
4095: coords[7*3+2] = upper[2];
4096: coords[8*3+0] = upper[0];
4097: coords[8*3+1] = ilower[1];
4098: coords[8*3+2] = upper[2];
4099: coords[9*3+0] = upper[0];
4100: coords[9*3+1] = ilower[1];
4101: coords[9*3+2] = ilower[2];
4102: coords[10*3+0] = ilower[0];
4103: coords[10*3+1] = upper[1];
4104: coords[10*3+2] = ilower[2];
4105: coords[11*3+0] = upper[0];
4106: coords[11*3+1] = upper[1];
4107: coords[11*3+2] = ilower[2];
4108: coords[12*3+0] = ilower[0];
4109: coords[12*3+1] = ilower[1];
4110: coords[12*3+2] = upper[2];
4111: coords[13*3+0] = ilower[0];
4112: coords[13*3+1] = ilower[1];
4113: coords[13*3+2] = ilower[2];
4115:
4116: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
4117: mesh->setSieve(sieve);
4118: typename Mesh::point_type p[nVertices];
4119: typename Mesh::point_type f[nFaces];
4120: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
4121: for (int i = 0; i < nVertices; i++) {
4122: p[i] = typename Mesh::point_type(i+nFaces);
4123: mesh->setValue(markers, p[i], 1);
4124: }
4125: for (int i = 0; i < nFaces; i++) {
4126: f[i] = typename Mesh::point_type(i);
4127: }
4128: int order = 0;
4129: //assemble the larger square sides
4130: sieve->addArrow(p[0], f[0], order++);
4131: sieve->addArrow(p[5], f[0], order++);
4132: sieve->addArrow(p[2], f[0], order++);
4133: sieve->addArrow(p[4], f[0], order++);
4134: mesh->setValue(markers, f[0], 1);
4136: sieve->addArrow(p[0], f[1], order++);
4137: sieve->addArrow(p[4], f[1], order++);
4138: sieve->addArrow(p[1], f[1], order++);
4139: sieve->addArrow(p[6], f[1], order++);
4140: mesh->setValue(markers, f[1], 1);
4142: sieve->addArrow(p[4], f[2], order++);
4143: sieve->addArrow(p[1], f[2], order++);
4144: sieve->addArrow(p[3], f[2], order++);
4145: sieve->addArrow(p[2], f[2], order++);
4146: mesh->setValue(markers, f[2], 1);
4147:
4148: //assemble the L-shaped sides
4150: sieve->addArrow(p[0], f[3], order++);
4151: sieve->addArrow(p[12], f[3], order++);
4152: sieve->addArrow(p[7], f[3], order++);
4153: sieve->addArrow(p[5], f[3], order++);
4154: mesh->setValue(markers, f[3], 1);
4156: sieve->addArrow(p[0], f[4], order++);
4157: sieve->addArrow(p[12],f[4], order++);
4158: sieve->addArrow(p[8], f[4], order++);
4159: sieve->addArrow(p[6], f[4], order++);
4160: mesh->setValue(markers, f[4], 1);
4162: sieve->addArrow(p[9], f[5], order++);
4163: sieve->addArrow(p[1], f[5], order++);
4164: sieve->addArrow(p[3], f[5], order++);
4165: sieve->addArrow(p[11], f[5], order++);
4166: mesh->setValue(markers, f[5], 1);
4168: sieve->addArrow(p[9], f[6], order++);
4169: sieve->addArrow(p[1], f[6], order++);
4170: sieve->addArrow(p[6], f[6], order++);
4171: sieve->addArrow(p[8], f[6], order++);
4172: mesh->setValue(markers, f[6], 1);
4174: sieve->addArrow(p[10], f[7], order++);
4175: sieve->addArrow(p[2], f[7], order++);
4176: sieve->addArrow(p[5], f[7], order++);
4177: sieve->addArrow(p[7], f[7], order++);
4178: mesh->setValue(markers, f[7], 1);
4180: sieve->addArrow(p[10], f[8], order++);
4181: sieve->addArrow(p[2], f[8], order++);
4182: sieve->addArrow(p[3], f[8], order++);
4183: sieve->addArrow(p[11], f[8], order++);
4184: mesh->setValue(markers, f[8], 1);
4186: //assemble the smaller square sides
4188: sieve->addArrow(p[13], f[9], order++);
4189: sieve->addArrow(p[10], f[9], order++);
4190: sieve->addArrow(p[11], f[9], order++);
4191: sieve->addArrow(p[9], f[9], order++);
4192: mesh->setValue(markers, f[9], 1);
4194: sieve->addArrow(p[12], f[10], order++);
4195: sieve->addArrow(p[7], f[10], order++);
4196: sieve->addArrow(p[10], f[10], order++);
4197: sieve->addArrow(p[13], f[10], order++);
4198: mesh->setValue(markers, f[10], 1);
4200: sieve->addArrow(p[8], f[11], order++);
4201: sieve->addArrow(p[12], f[11], order++);
4202: sieve->addArrow(p[13], f[11], order++);
4203: sieve->addArrow(p[9], f[11], order++);
4204: mesh->setValue(markers, f[11], 1);
4206: mesh->stratify();
4207: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
4208: Obj<typename Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
4209: //coordinates->view("coordinates");
4210: return mesh;
4212: }
4216: /*
4217: //"sphere" out a cube
4219: */
4220: #if 0
4221: static Obj<Mesh> createSphereBoundary(const MPI_Comm comm, const real radius, const int refinement, const int debug = 0) {
4222: Obj<Mesh> m = new Mesh(comm, 2, debug);
4223: Obj<Mesh::sieve_type> s = new Mesh::sieve_type(comm, debug);
4224: m->setSieve(s);
4225: Mesh::point_type p = 0;
4226: int nVertices = 8+12*(refinement)+6*(refinement)*(refinement);
4227: Mesh::point_type vertices[nVertices];
4228: real coords[3*nVertices];
4229: int nCells = 6*2*(refinement+1)*(refinement+1);
4230: real delta = 2./((real)(refinement+1));
4231: Mesh::point_type cells[nCells];
4232: for (int i = 0; i < nCells; i++) {
4233: cells[i] = p;
4234: p++;
4235: }
4236: for (int i = 0; i < nVertices; i++) {
4237: vertices[i] = p;
4238: p++;
4239: }
4240: //set up the corners;
4241: //lll
4242: coords[0*3+0] = -1.;
4243: coords[0*3+1] = -1.;
4244: coords[0*3+2] = -1.;
4245: //llh
4246: coords[1*3+0] = -1.;
4247: coords[1*3+1] = -1.;
4248: coords[1*3+2] = 1.;
4249: //lhh
4250: coords[2*3+0] = -1.;
4251: coords[2*3+1] = 1.;
4252: coords[2*3+2] = 1.;
4253: //lhl
4254: coords[3*3+0] = -1.;
4255: coords[3*3+1] = 1.;
4256: coords[3*3+2] = -1.;
4257: //hhl
4258: coords[4*3+0] = 1.;
4259: coords[4*3+1] = 1.;
4260: coords[4*3+2] = -1.;
4261: //hhh
4262: coords[5*3+0] = 1.;
4263: coords[5*3+1] = 1.;
4264: coords[5*3+2] = 1.;
4265: //hlh
4266: coords[6*3+0] = 1.;
4267: coords[6*3+1] = -1.;
4268: coords[6*3+2] = 1.;
4269: //hll
4270: coords[7*3+0] = 1.;
4271: coords[7*3+1] = -1.;
4272: coords[7*3+2] = -1.;
4273: //set up the edges (always go low to high)
4274: //xll
4275: for (int i = 0; i < refinement; i++) {
4276: coords[3*(8+0*refinement+i)+0] = -1. + delta*i;
4277: coords[3*(8+0*refinement+i)+1] = -1.;
4278: coords[3*(8+0*refinement+i)+2] = -1.;
4279: }
4280: //xlh
4281: for (int i = 0; i < refinement; i++) {
4282: coords[3*(8+1*refinement+i)+0] = -1. + delta*i;
4283: coords[3*(8+1*refinement+i)+1] = -1.;
4284: coords[3*(8+1*refinement+i)+2] = 1.;
4285: }
4286: //xhh
4287: for (int i = 0; i < refinement; i++) {
4288: coords[3*(8+2*refinement+i)+0] = -1. + delta*i;
4289: coords[3*(8+2*refinement+i)+1] = 1.;
4290: coords[3*(8+2*refinement+i)+2] = 1.;
4291: }
4292: //xhl
4293: for (int i = 0; i < refinement; i++) {
4294: coords[3*(8+3*refinement+i)+0] = -1. + delta*i;
4295: coords[3*(8+3*refinement+i)+1] = 1.;
4296: coords[3*(8+3*refinement+i)+2] = -1.;
4297: }
4298: //lxl
4299: for (int i = 0; i < refinement; i++) {
4300: coords[3*(8+4*refinement+i)+0] = -1.;
4301: coords[3*(8+4*refinement+i)+1] = -1. + delta*i;
4302: coords[3*(8+4*refinement+i)+2] = -1.;
4303: }
4304: //lxh
4305: for (int i = 0; i < refinement; i++) {
4306: coords[3*(8+5*refinement+i)+0] = -1.;
4307: coords[3*(8+5*refinement+i)+1] = -1. + delta*i;
4308: coords[3*(8+5*refinement+i)+2] = 1.;
4309: }
4310: //hxh
4311: for (int i = 0; i < refinement; i++) {
4312: coords[3*(8+6*refinement+i)+0] = 1.;
4313: coords[3*(8+6*refinement+i)+1] = -1. + delta*i;
4314: coords[3*(8+6*refinement+i)+2] = 1.;
4315: }
4316: //hxl
4317: for (int i = 0; i < refinement; i++) {
4318: coords[3*(8+7*refinement+i)+0] = 1.;
4319: coords[3*(8+7*refinement+i)+1] = -1. + delta*i;
4320: coords[3*(8+7*refinement+i)+2] = -1.;
4321: }
4322: //llx
4323: for (int i = 0; i < refinement; i++) {
4324: coords[3*(8+8*refinement+i)+0] = -1.;
4325: coords[3*(8+8*refinement+i)+1] = -1.;
4326: coords[3*(8+8*refinement+i)+2] = -1. + delta*i;
4327: }
4328: //lhx
4329: for (int i = 0; i < refinement; i++) {
4330: coords[3*(8+9*refinement+i)+0] = -1.;
4331: coords[3*(8+9*refinement+i)+1] = 1.;
4332: coords[3*(8+9*refinement+i)+2] = -1. + delta*i;
4333: }
4334: //hhx
4335: for (int i = 0; i < refinement; i++) {
4336: coords[3*(8+10*refinement+i)+0] = 1.;
4337: coords[3*(8+10*refinement+i)+1] = 1.;
4338: coords[3*(8+10*refinement+i)+2] = -1. + delta*i;
4339: }
4340: //hlx
4341: for (int i = 0; i < refinement; i++) {
4342: coords[3*(8+11*refinement+i)+0] = 1.;
4343: coords[3*(8+11*refinement+i)+1] = -1.;
4344: coords[3*(8+11*refinement+i)+2] = -1. + delta*i;
4345: }
4346: //set up the faces
4347: //lxx
4348: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4349: coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+0] = -1.;
4350: coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4351: coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4352: }
4353: //hxx
4354: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4355: coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+0] = 1.;
4356: coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4357: coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4358: }
4359: //xlx
4360: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4361: coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+0] = -1. + delta*j;
4362: coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+1] = -1.;
4363: coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4364: }
4365: //xhx
4366: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4367: coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+0] = -1. + delta*j;
4368: coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+1] = 1.;
4369: coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4370: }
4371: //xxl
4372: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4373: coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+0] = -1.;
4374: coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4375: coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4376: }
4377: //xxh
4378: for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4379: coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+0] = 1.;
4380: coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4381: coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4382: }
4383: //stitch the corners up with the edges and the faces
4384:
4385: //stitch the edges to the faces
4386: //fill in the faces
4387: int face_offset = 8 + 12*refinement;
4388: for (int i = 0; i < 6; i++) for (int j = 0; j < refinement; j++) for (int k = 0; k < refinement; k++) {
4389: //build each square doublet
4390: }
4391: }
4393: #endif
4397: /*
4398: Simple cubic boundary:
4400: 30----31-----32
4401: | | |
4402: | 3 | 2 |
4403: | | |
4404: 27----28-----29
4405: | | |
4406: | 0 | 1 |
4407: | | |
4408: 24----25-----26
4409: */
4410: static Obj<Mesh> createParticleInCubeBoundary(const MPI_Comm comm, const real lower[], const real upper[], const int faces[], const real radius, const int thetaEdges, const int phiSlices, const int debug = 0) {
4411: Obj<Mesh> mesh = new Mesh(comm, 2, debug);
4412: const int numCubeVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
4413: const int numPartVertices = (thetaEdges - 1)*phiSlices + 2;
4414: const int numVertices = numCubeVertices + numPartVertices;
4415: const int numCubeFaces = 6;
4416: const int numFaces = numCubeFaces + thetaEdges*phiSlices;
4417: real *coords = new real[numVertices*3];
4418: const Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
4419: typename Mesh::point_type *vertices = new typename Mesh::point_type[numVertices];
4420: int order = 0;
4422: mesh->setSieve(sieve);
4423: const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
4424: if (mesh->commRank() == 0) {
4425: // Make cube
4426: for(int v = numFaces; v < numFaces+numVertices; v++) {
4427: vertices[v-numFaces] = typename Mesh::point_type(v);
4428: mesh->setValue(markers, vertices[v-numFaces], 1);
4429: }
4430: {
4431: // Side 0 (Front)
4432: typename Mesh::point_type face(0);
4433: sieve->addArrow(vertices[0], face, order++);
4434: sieve->addArrow(vertices[1], face, order++);
4435: sieve->addArrow(vertices[2], face, order++);
4436: sieve->addArrow(vertices[3], face, order++);
4437: mesh->setValue(markers, face, 1);
4438: }
4439: {
4440: // Side 1 (Back)
4441: typename Mesh::point_type face(1);
4442: sieve->addArrow(vertices[5], face, order++);
4443: sieve->addArrow(vertices[4], face, order++);
4444: sieve->addArrow(vertices[7], face, order++);
4445: sieve->addArrow(vertices[6], face, order++);
4446: mesh->setValue(markers, face, 1);
4447: }
4448: {
4449: // Side 2 (Bottom)
4450: typename Mesh::point_type face(2);
4451: sieve->addArrow(vertices[4], face, order++);
4452: sieve->addArrow(vertices[5], face, order++);
4453: sieve->addArrow(vertices[1], face, order++);
4454: sieve->addArrow(vertices[0], face, order++);
4455: mesh->setValue(markers, face, 1);
4456: }
4457: {
4458: // Side 3 (Top)
4459: typename Mesh::point_type face(3);
4460: sieve->addArrow(vertices[3], face, order++);
4461: sieve->addArrow(vertices[2], face, order++);
4462: sieve->addArrow(vertices[6], face, order++);
4463: sieve->addArrow(vertices[7], face, order++);
4464: mesh->setValue(markers, face, 1);
4465: }
4466: {
4467: // Side 4 (Left)
4468: typename Mesh::point_type face(4);
4469: sieve->addArrow(vertices[4], face, order++);
4470: sieve->addArrow(vertices[0], face, order++);
4471: sieve->addArrow(vertices[3], face, order++);
4472: sieve->addArrow(vertices[7], face, order++);
4473: mesh->setValue(markers, face, 1);
4474: }
4475: {
4476: // Side 5 (Right)
4477: typename Mesh::point_type face(5);
4478: sieve->addArrow(vertices[1], face, order++);
4479: sieve->addArrow(vertices[5], face, order++);
4480: sieve->addArrow(vertices[6], face, order++);
4481: sieve->addArrow(vertices[2], face, order++);
4482: mesh->setValue(markers, face, 1);
4483: }
4484: // Make particle
4485: for(int s = 0; s < phiSlices; ++s) {
4486: for(int ep = 0; ep < thetaEdges; ++ep) {
4487: // Vertices on each slice are 0..thetaEdges
4488: typename Mesh::point_type face(numCubeFaces + s*thetaEdges + ep);
4489: int vertexA = numCubeVertices + ep + 0 + s*(thetaEdges+1);
4490: int vertexB = numCubeVertices + ep + 1 + s*(thetaEdges+1);
4491: int vertexC = numCubeVertices + (ep + 1 + (s+1)*(thetaEdges+1))%((thetaEdges+1)*phiSlices);
4492: int vertexD = numCubeVertices + (ep + 0 + (s+1)*(thetaEdges+1))%((thetaEdges+1)*phiSlices);
4493: const int correction1 = (s > 0)*((s-1)*2 + 1);
4494: const int correction2 = (s < phiSlices-1)*(s*2 + 1);
4496: if ((vertexA - numCubeVertices)%(thetaEdges+1) == 0) {
4497: vertexA = vertexD = numCubeVertices;
4498: vertexB -= correction1;
4499: vertexC -= correction2;
4500: } else if ((vertexB - numCubeVertices)%(thetaEdges+1) == thetaEdges) {
4501: vertexA -= correction1;
4502: vertexD -= correction2;
4503: vertexB = vertexC = numCubeVertices + thetaEdges;
4504: } else {
4505: vertexA -= correction1;
4506: vertexB -= correction1;
4507: vertexC -= correction2;
4508: vertexD -= correction2;
4509: }
4510: if ((vertexA >= numVertices) || (vertexB >= numVertices) || (vertexC >= numVertices) || (vertexD >= numVertices)) {
4511: throw ALE::Exception("Bad vertex");
4512: }
4513: sieve->addArrow(vertices[vertexA], face, order++);
4514: sieve->addArrow(vertices[vertexB], face, order++);
4515: if (vertexB != vertexC) sieve->addArrow(vertices[vertexC], face, order++);
4516: if (vertexA != vertexD) sieve->addArrow(vertices[vertexD], face, order++);
4517: mesh->setValue(markers, face, 2);
4518: mesh->setValue(markers, vertices[vertexA], 2);
4519: mesh->setValue(markers, vertices[vertexB], 2);
4520: if (vertexB != vertexC) mesh->setValue(markers, vertices[vertexC], 2);
4521: if (vertexA != vertexD) mesh->setValue(markers, vertices[vertexD], 2);
4522: }
4523: }
4524: }
4525: mesh->stratify();
4526: #if 0
4527: for(int vz = 0; vz <= edges[2]; ++vz) {
4528: for(int vy = 0; vy <= edges[1]; ++vy) {
4529: for(int vx = 0; vx <= edges[0]; ++vx) {
4530: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
4531: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
4532: coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
4533: }
4534: }
4535: }
4536: #else
4537: coords[0*3+0] = lower[0];
4538: coords[0*3+1] = lower[1];
4539: coords[0*3+2] = upper[2];
4540: coords[1*3+0] = upper[0];
4541: coords[1*3+1] = lower[1];
4542: coords[1*3+2] = upper[2];
4543: coords[2*3+0] = upper[0];
4544: coords[2*3+1] = upper[1];
4545: coords[2*3+2] = upper[2];
4546: coords[3*3+0] = lower[0];
4547: coords[3*3+1] = upper[1];
4548: coords[3*3+2] = upper[2];
4549: coords[4*3+0] = lower[0];
4550: coords[4*3+1] = lower[1];
4551: coords[4*3+2] = lower[2];
4552: coords[5*3+0] = upper[0];
4553: coords[5*3+1] = lower[1];
4554: coords[5*3+2] = lower[2];
4555: coords[6*3+0] = upper[0];
4556: coords[6*3+1] = upper[1];
4557: coords[6*3+2] = lower[2];
4558: coords[7*3+0] = lower[0];
4559: coords[7*3+1] = upper[1];
4560: coords[7*3+2] = lower[2];
4561: #endif
4562: const real centroidX = 0.5*(upper[0] + lower[0]);
4563: const real centroidY = 0.5*(upper[1] + lower[1]);
4564: const real centroidZ = 0.5*(upper[2] + lower[2]);
4565: for(int s = 0; s < phiSlices; ++s) {
4566: for(int v = 0; v <= thetaEdges; ++v) {
4567: int vertex = numCubeVertices + v + s*(thetaEdges+1);
4568: const real theta = v*(PETSC_PI/thetaEdges);
4569: const real phi = s*(2.0*PETSC_PI/phiSlices);
4570: const int correction = (s > 0)*((s-1)*2 + 1);
4572: if ((vertex- numCubeVertices)%(thetaEdges+1) == 0) {
4573: vertex = numCubeVertices;
4574: } else if ((vertex - numCubeVertices)%(thetaEdges+1) == thetaEdges) {
4575: vertex = numCubeVertices + thetaEdges;
4576: } else {
4577: vertex -= correction;
4578: }
4579: coords[vertex*3+0] = centroidX + radius*sin(theta)*cos(phi);
4580: coords[vertex*3+1] = centroidY + radius*sin(theta)*sin(phi);
4581: coords[vertex*3+2] = centroidZ + radius*cos(theta);
4582: }
4583: }
4584: delete [] vertices;
4585: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
4586: return mesh;
4587: };
4589: #if 0 // WORKING ON REDESIGN WITHIN PYLITH SOURCE TREE
4590: template<typename MeshType, typename EdgeType>
4591: class CellRefiner {
4592: public:
4593: typedef typename MeshType::point_type point_type;
4594: typedef EdgeType edge_type;
4595: typedef typename std::map<edge_type, point_type> edge_map_type;
4596: typedef enum {LINE, LINE_LAGRANGE, TRIANGLE, QUADRILATERAL, TETRAHEDRON, HEXAHEDRON, TRIANGULAR_PRISM, TRIANGULAR_PRISM_LAGRANGE, HEXAHEDRON_LAGRANGE} CellType;
4597: protected:
4598: const MeshType& _mesh;
4599: int _dim;
4600: point_type _vertexOffset;
4601: edge_map_type _edge2vertex;
4602: public:
4603: CellRefiner(MeshType& mesh) : _mesh(mesh) {
4604: _dim = _mesh.getDimension();
4605: };
4606: ~CellRefiner() {};
4607: protected:
4608: CellType getCellType(const point_type cell) {
4609: const int corners = _mesh.getSieve()->getConeSize(cell);
4610: switch (_dim) {
4611: return LINE;
4612: case 2:
4613: switch (corners) {
4614: case 3:
4615: return TRIANGLE;
4616: case 4:
4617: throw ALE::Exception("Not implemented.");
4618: return QUADRILATERAL;
4619: case 6:
4620: return LINE_LAGRANGE;
4621: case 0: {
4622: std::ostringstream msg;
4623: std::cerr << "Internal error. Cone size for mesh point " << cell << " is zero. May be a vertex.";
4624: assert(0);
4625: throw ALE::Exception("Could not determine 2-D cell type.");
4626: } // case 0
4627: default : {
4628: std::ostringstream msg;
4629: std::cerr << "Internal error. Unknown cone size for mesh point " << cell << ". Unknown cell type.";
4630: assert(0);
4631: throw ALE::Exception("Could not determine 2-D cell type.");
4632: } // default
4633: } // switch
4634: case 3:
4635: switch (corners) {
4636: case 4:
4637: return TETRAHEDRON;
4638: case 6:
4639: return TRIANGULAR_PRISM;
4640: case 9:
4641: return TRIANGULAR_PRISM_LAGRANGE;
4642: case 12:
4643: throw ALE::Exception("Not implemented.");
4644: return HEXAHEDRON_LAGRANGE;
4645: case 0: {
4646: std::ostringstream msg;
4647: std::cerr << "Internal error. Cone size for mesh point " << cell << " is zero. May be a vertex.";
4648: assert(0);
4649: throw ALE::Exception("Could not determine 3-D cell type.");
4650: } // case 0
4651: default : {
4652: std::ostringstream msg;
4653: std::cerr << "Internal error. Unknown cone size for mesh point " << cell << ". Unknown cell type.";
4654: assert(0);
4655: throw ALE::Exception("Could not determine 3-D cell type.");
4656: } // default
4657: } // switch
4658: } // switch
4659: };
4660: void getEdges_TRIANGLE(const int coneSize, const point_type cone[], int *numEdges, const edge_type **edges) {
4661: static edge_type triEdges[3];
4663: assert(coneSize == 3);
4664: triEdges[0] = edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
4665: triEdges[1] = edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
4666: triEdges[2] = edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0]));
4667: *numEdges = 3;
4668: *edges = triEdges;
4669: };
4670: void getEdges_LINE_LAGRANGE(const int coneSize, const point_type cone[], int *numEdges, const edge_type **edges) {
4671: static edge_type lineEdges[6];
4673: assert(coneSize == 6);
4674: lineEdges[0] = edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
4675: lineEdges[1] = edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]));
4676: lineEdges[2] = edge_type(std::min(cone[4], cone[5]), std::max(cone[4], cone[5]));
4677: *numEdges = 3;
4678: *edges = lineEdges;
4679: };
4680: void getEdges_TETRAHEDRON(const int coneSize, const point_type cone[], int *numEdges, const edge_type **edges) {
4681: static edge_type tetEdges[6];
4683: assert(coneSize == 4);
4684: // As per Brad's diagram
4685: tetEdges[0] = edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
4686: tetEdges[1] = edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
4687: tetEdges[2] = edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0]));
4688: tetEdges[3] = edge_type(std::min(cone[0], cone[3]), std::max(cone[0], cone[3]));
4689: tetEdges[4] = edge_type(std::min(cone[1], cone[3]), std::max(cone[1], cone[3]));
4690: tetEdges[5] = edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]));
4691: *numEdges = 6;
4692: *edges = tetEdges;
4693: };
4694: void getEdges_TRIANGULAR_PRISM(const int coneSize, const point_type cone[], int *numEdges, const edge_type **edges) {
4695: static edge_type triPrismEdges[6];
4697: assert(coneSize == 6);
4698: triPrismEdges[0] = edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
4699: triPrismEdges[1] = edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
4700: triPrismEdges[2] = edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0]));
4701: triPrismEdges[3] = edge_type(std::min(cone[3], cone[4]), std::max(cone[3], cone[4]));
4702: triPrismEdges[4] = edge_type(std::min(cone[4], cone[5]), std::max(cone[4], cone[5]));
4703: triPrismEdges[5] = edge_type(std::min(cone[5], cone[3]), std::max(cone[5], cone[3]));
4704: *numEdges = 6;
4705: *edges = triPrismEdges;
4706: };
4707: void getEdges_TRIANGULAR_PRISM_LAGRANGE(const int coneSize, const point_type cone[], int *numEdges, const edge_type **edges) {
4708: static edge_type triPrismLEdges[9];
4710: assert(coneSize == 9);
4711: triPrismLEdges[0] = edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
4712: triPrismLEdges[1] = edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
4713: triPrismLEdges[2] = edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0]));
4714: triPrismLEdges[3] = edge_type(std::min(cone[3], cone[4]), std::max(cone[3], cone[4]));
4715: triPrismLEdges[4] = edge_type(std::min(cone[4], cone[5]), std::max(cone[4], cone[5]));
4716: triPrismLEdges[5] = edge_type(std::min(cone[5], cone[3]), std::max(cone[5], cone[3]));
4717: triPrismLEdges[6] = edge_type(std::min(cone[6], cone[7]), std::max(cone[6], cone[7]));
4718: triPrismLEdges[7] = edge_type(std::min(cone[7], cone[8]), std::max(cone[7], cone[8]));
4719: triPrismLEdges[8] = edge_type(std::min(cone[8], cone[6]), std::max(cone[8], cone[6]));
4720: *numEdges = 9;
4721: *edges = triPrismLEdges;
4722: };
4723: void getNewCells_TRIANGLE(const int coneSize, const point_type cone[], int *numCells, const point_type **cells) {
4724: int numEdges;
4725: const edge_type *edges;
4726: static point_type triCells[4*3];
4727: point_type newVertices[3];
4729: getEdges_TRIANGLE(coneSize, cone, &numEdges, &edges);
4730: assert(numEdges == 3);
4731: for(int e = 0; e < numEdges; ++e) {
4732: if (_edge2vertex.find(edges[e]) == _edge2vertex.end()) {
4733: throw ALE::Exception("Missing edge in refined mesh");
4734: }
4735: newVertices[e] = _edge2vertex[edges[e]];
4736: }
4737: triCells[0*3+0] = cone[0]+_vertexOffset; triCells[0*3+1] = newVertices[0]; triCells[0*3+2] = newVertices[2];
4738: triCells[1*3+0] = newVertices[0]; triCells[1*3+1] = newVertices[1]; triCells[1*3+2] = newVertices[2];
4739: triCells[2*3+0] = cone[1]+_vertexOffset; triCells[2*3+1] = newVertices[1]; triCells[2*3+2] = newVertices[0];
4740: triCells[3*3+0] = cone[2]+_vertexOffset; triCells[3*3+1] = newVertices[2]; triCells[3*3+2] = newVertices[1];
4741: *numCells = 4;
4742: *cells = triCells;
4743: };
4744: void getNewCells_LINE_LAGRANGE(const int coneSize, const point_type cone[], int *numCells, const point_type **cells) {
4745: int numEdges;
4746: const edge_type *edges;
4747: static point_type lineCells[2*6];
4748: point_type newVertices[3];
4750: getEdges_LINE_LAGRANGE(coneSize, cone, &numEdges, &edges);
4751: assert(numEdges == 3);
4752: for(int e = 0; e < numEdges; ++e) {
4753: if (_edge2vertex.find(edges[e]) == _edge2vertex.end()) {
4754: throw ALE::Exception("Missing edge in refined mesh");
4755: }
4756: newVertices[e] = _edge2vertex[edges[e]];
4757: }
4758: lineCells[0*6+0] = cone[0]+_vertexOffset; // new cell 0
4759: lineCells[0*6+1] = newVertices[0];
4760: lineCells[0*6+2] = cone[2]+_vertexOffset;
4761: lineCells[0*6+3] = newVertices[1];
4762: lineCells[0*6+4] = cone[4]+_vertexOffset;
4763: lineCells[0*6+5] = newVertices[2];
4765: lineCells[1*6+0] = newVertices[0]; // new cell 1
4766: lineCells[1*6+1] = cone[1]+_vertexOffset;
4767: lineCells[1*6+2] = newVertices[1];
4768: lineCells[1*6+3] = cone[3]+_vertexOffset;
4769: lineCells[1*6+4] = newVertices[2];
4770: lineCells[1*6+5] = cone[5]+_vertexOffset;
4772: *numCells = 2;
4773: *cells = lineCells;
4774: };
4775: void getNewCells_TETRAHEDRON(const int coneSize, const point_type cone[], int *numCells, const point_type **cells) {
4776: int numEdges;
4777: const edge_type *edges;
4778: static point_type tetCells[8*4];
4779: point_type newVertices[6];
4781: getEdges_TETRAHEDRON(coneSize, cone, &numEdges, &edges);
4782: assert(numEdges == 6);
4783: for(int e = 0; e < numEdges; ++e) {
4784: if (_edge2vertex.find(edges[e]) == _edge2vertex.end()) {
4785: throw ALE::Exception("Missing edge in refined mesh");
4786: }
4787: newVertices[e] = _edge2vertex[edges[e]];
4788: }
4789: tetCells[0*4+0] = cone[0]+_vertexOffset; tetCells[0*4+1] = newVertices[3]; tetCells[0*4+2] = newVertices[0]; tetCells[0*4+3] = newVertices[2];
4790: tetCells[1*4+0] = newVertices[0]; tetCells[1*4+1] = newVertices[1]; tetCells[1*4+2] = newVertices[2]; tetCells[1*4+3] = newVertices[3];
4791: tetCells[2*4+0] = newVertices[0]; tetCells[2*4+1] = newVertices[3]; tetCells[2*4+2] = newVertices[4]; tetCells[2*4+3] = newVertices[1];
4792: tetCells[3*4+0] = cone[1]+_vertexOffset; tetCells[3*4+1] = newVertices[4]; tetCells[3*4+2] = newVertices[1]; tetCells[3*4+3] = newVertices[0];
4793: tetCells[4*4+0] = newVertices[2]; tetCells[4*4+1] = newVertices[5]; tetCells[4*4+2] = newVertices[3]; tetCells[4*4+3] = newVertices[1];
4794: tetCells[5*4+0] = cone[2]+_vertexOffset; tetCells[5*4+1] = newVertices[5]; tetCells[5*4+2] = newVertices[2]; tetCells[5*4+3] = newVertices[1];
4795: tetCells[6*4+0] = newVertices[1]; tetCells[6*4+1] = newVertices[4]; tetCells[6*4+2] = newVertices[5]; tetCells[6*4+3] = newVertices[3];
4796: tetCells[7*4+0] = cone[3]+_vertexOffset; tetCells[7*4+1] = newVertices[3]; tetCells[7*4+2] = newVertices[5]; tetCells[7*4+3] = newVertices[4];
4797: *numCells = 8;
4798: *cells = tetCells;
4799: };
4800: void getNewCells_TRIANGULAR_PRISM_LAGRANGE(const int coneSize, const point_type cone[], int *numCells, const point_type **cells) {
4801: int numEdges;
4802: const edge_type *edges;
4803: static point_type tcells[4*9];
4804: point_type newVertices[9];
4806: getEdges_TRIANGULAR_PRISM_LAGRANGE(coneSize, cone, &numEdges, &edges);
4807: assert(numEdges == 9);
4808: for(int e = 0; e < numEdges; ++e) {
4809: if (_edge2vertex.find(edges[e]) == _edge2vertex.end()) {
4810: throw ALE::Exception("Missing edge in refined mesh");
4811: }
4812: newVertices[e] = _edge2vertex[edges[e]];
4813: }
4814: tcells[0*9+0] = cone[0]+_vertexOffset; // New cell 0
4815: tcells[0*9+1] = newVertices[0];
4816: tcells[0*9+2] = newVertices[2];
4817: tcells[0*9+3] = cone[3]+_vertexOffset;
4818: tcells[0*9+4] = newVertices[3];
4819: tcells[0*9+5] = newVertices[5];
4820: tcells[0*9+6] = cone[6]+_vertexOffset;
4821: tcells[0*9+7] = newVertices[6];
4822: tcells[0*9+8] = newVertices[8];
4824: tcells[1*9+0] = newVertices[0]; // New cell 1
4825: tcells[1*9+1] = newVertices[1];
4826: tcells[1*9+2] = newVertices[2];
4827: tcells[1*9+3] = newVertices[3];
4828: tcells[1*9+4] = newVertices[4];
4829: tcells[1*9+5] = newVertices[5];
4830: tcells[1*9+6] = newVertices[6];
4831: tcells[1*9+7] = newVertices[7];
4832: tcells[1*9+8] = newVertices[8];
4834: tcells[2*9+0] = cone[1]+_vertexOffset; // New cell 2
4835: tcells[2*9+1] = newVertices[1];
4836: tcells[2*9+2] = newVertices[0];
4837: tcells[2*9+3] = cone[4]+_vertexOffset;
4838: tcells[2*9+4] = newVertices[4];
4839: tcells[2*9+5] = newVertices[3];
4840: tcells[2*9+6] = cone[7]+_vertexOffset;
4841: tcells[2*9+7] = newVertices[7];
4842: tcells[2*9+8] = newVertices[6];
4844: tcells[3*9+0] = cone[2]+_vertexOffset; // New cell 3
4845: tcells[3*9+1] = newVertices[2];
4846: tcells[3*9+2] = newVertices[1];
4847: tcells[3*9+3] = cone[5]+_vertexOffset;
4848: tcells[3*9+4] = newVertices[5];
4849: tcells[3*9+5] = newVertices[4];
4850: tcells[3*9+6] = cone[8]+_vertexOffset;
4851: tcells[3*9+7] = newVertices[8];
4852: tcells[3*9+8] = newVertices[7];
4854: *numCells = 4;
4855: *cells = tcells;
4856: };
4857: public:
4858: point_type getVertexRelativeOffset() {return _vertexOffset;};
4859: void setVertexRelativeOffset(const point_type offset) {_vertexOffset = offset;};
4860: edge_map_type& getEdgeToVertex() {return _edge2vertex;};
4861: int numNewCells(const point_type cell) {
4862: switch(this->getCellType(cell)) {
4863: case TRIANGLE:
4864: return 4;
4865: case LINE_LAGRANGE:
4866: return 2;
4867: case TETRAHEDRON:
4868: return 8;
4869: case TRIANGULAR_PRISM:
4870: case TRIANGULAR_PRISM_LAGRANGE:
4871: return 4;
4872: }
4873: throw ALE::Exception("Could not determine number of new cells for this cell type");
4874: };
4875: void splitEdge(const point_type cell, const int coneSize, const point_type cone[], point_type& curNewVertex) {
4876: const CellType t = this->getCellType(cell);
4877: int numEdges;
4878: const edge_type *edges;
4880: switch(t) {
4881: case TRIANGLE:
4882: getEdges_TRIANGLE(coneSize, cone, &numEdges, &edges);
4883: break;
4884: case LINE_LAGRANGE:
4885: getEdges_LINE_LAGRANGE(coneSize, cone, &numEdges, &edges);
4886: break;
4887: case TETRAHEDRON:
4888: getEdges_TETRAHEDRON(coneSize, cone, &numEdges, &edges);
4889: break;
4890: case TRIANGULAR_PRISM:
4891: getEdges_TRIANGULAR_PRISM(coneSize, cone, &numEdges, &edges);
4892: break;
4893: case TRIANGULAR_PRISM_LAGRANGE:
4894: getEdges_TRIANGULAR_PRISM_LAGRANGE(coneSize, cone, &numEdges, &edges);
4895: break;
4896: default:
4897: throw ALE::Exception("Could not determine number of new cells for this cell type");
4898: }
4899: // Check that vertex does not yet exist
4900: for(int v = 0; v < numEdges; ++v) {
4901: if (_edge2vertex.find(edges[v]) == _edge2vertex.end()) {
4902: std::cout << "Edge: " << edges[v] << ", new vertex: " << curNewVertex << std::endl;
4903: _edge2vertex[edges[v]] = curNewVertex++;
4904: }
4905: }
4906: };
4907: void getNewCell(const point_type cell, const int coneSize, const point_type cone[], int newCellNumber, int *newConeSize, const point_type **newCone) {
4908: const CellType t = this->getCellType(cell);
4909: int numCells;
4910: const point_type *cells;
4912: switch(t) {
4913: case TRIANGLE:
4914: getNewCells_TRIANGLE(coneSize, cone, &numCells, &cells);
4915: *newConeSize = 3;
4916: *newCone = &cells[newCellNumber*3];
4917: break;
4918: case LINE_LAGRANGE:
4919: getNewCells_LINE_LAGRANGE(coneSize, cone, &numCells, &cells);
4920: *newConeSize = 6;
4921: *newCone = &cells[newCellNumber*6];
4922: break;
4923: case TETRAHEDRON:
4924: getNewCells_TETRAHEDRON(coneSize, cone, &numCells, &cells);
4925: *newConeSize = 4;
4926: *newCone = &cells[newCellNumber*4];
4927: break;
4928: case TRIANGULAR_PRISM_LAGRANGE:
4929: getNewCells_TRIANGULAR_PRISM_LAGRANGE(coneSize, cone, &numCells, &cells);
4930: *newConeSize = 9;
4931: *newCone = &cells[newCellNumber*9];
4932: break;
4933: default:
4934: throw ALE::Exception("Could not create new cell for this cell type");
4935: }
4936: };
4937: void getNeighboringVertices(const point_type cell, const int coneSize, const point_type cone[], const point_type firstNewVertex, point_type vertex2edge[]) {
4938: const CellType t = this->getCellType(cell);
4939: int numEdges;
4940: const edge_type *edges;
4942: switch(t) {
4943: case TRIANGLE:
4944: getEdges_TRIANGLE(coneSize, cone, &numEdges, &edges);
4945: break;
4946: case LINE_LAGRANGE:
4947: getEdges_LINE_LAGRANGE(coneSize, cone, &numEdges, &edges);
4948: break;
4949: case TETRAHEDRON:
4950: getEdges_TETRAHEDRON(coneSize, cone, &numEdges, &edges);
4951: break;
4952: case TRIANGULAR_PRISM:
4953: getEdges_TRIANGULAR_PRISM(coneSize, cone, &numEdges, &edges);
4954: break;
4955: case TRIANGULAR_PRISM_LAGRANGE:
4956: getEdges_TRIANGULAR_PRISM_LAGRANGE(coneSize, cone, &numEdges, &edges);
4957: break;
4958: default:
4959: throw ALE::Exception("Could not determine number of new cells for this cell type");
4960: }
4961: for(int v = 0; v < numEdges; ++v) {
4962: point_type newVertex = _edge2vertex[edges[v]];
4964: std::cout << "VERTEX2EDGE index: " << newVertex-firstNewVertex << ", first: " << edges[v].first << ", second: " << edges[v].second << std::endl;
4965: vertex2edge[(newVertex-firstNewVertex)*2+0] = edges[v].first;
4966: vertex2edge[(newVertex-firstNewVertex)*2+1] = edges[v].second;
4967: }
4968: };
4969: };
4970: // This method takes a mesh and performs a refinement of each cell
4971: //
4972: // triangle: 1 --> 4 refinement, adding a new vertex at the midpoint of each edge
4973: // tetrahedra: 1 --> 8 refinement, adding a new vertex at the midpoint of each edge
4974: //
4975: // :WARNING: This method currently only works for uninterpolated meshes with tri and tet cells.
4976: template<typename MeshType, typename Refiner>
4977: static void refineGeneral(const Obj<MeshType>& mesh, const Obj<MeshType>& newMesh, Refiner& refiner) {
4978: typedef typename MeshType::sieve_type sieve_type;
4979: typedef typename MeshType::point_type point_type;
4980: typedef typename Refiner::edge_type edge_type;
4982: // :WARNING: Assumed order of mesh points (cells and vertices):
4983: //
4984: // normal cells (in censored depth)
4985: // normal vertices (in censored depth)
4986: // other vertices
4987: // other cells
4988: //
4989: // This permits omitting in output the other vertices (e.g.,
4990: // Lagrange multipliers) and other cells (e.g., cohesive cells)
4991: // which have a custom reference cell that is not recognized.
4993: assert(!mesh.isNull());
4994: assert(!newMesh.isNull());
4996: // Get original mesh stuff.
4997: const Obj<typename MeshType::label_sequence>& cells = mesh->heightStratum(0);
4998: assert(!cells.isNull());
4999: const typename MeshType::label_sequence::iterator cellsEnd = cells->end();
5001: const Obj<typename MeshType::label_sequence>& vertices = mesh->depthStratum(0);
5002: assert(!vertices.isNull());
5004: const Obj<sieve_type>& sieve = mesh->getSieve();
5005: assert(!sieve.isNull());
5006: ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
5008: if (mesh->hasLabel("censored depth")) {
5009: // :WARNING: Assume all cells in the censored depth come before
5010: // any other cells. This guarantees that we add vertices in the
5011: // censored depth before adding other vertices.
5013: int counterBegin = 0;
5015: int oldNumCellsNormal = 0;
5016: int oldNumCellsOther = 0;
5017: int oldNumVerticesNormal = 0;
5018: int oldNumVerticesOther = 0;
5020: int newNumCellsNormal = 0;
5021: int newNumCellsOther = 0;
5022: int newNumVerticesNormal = 0;
5023: int newNumVerticesOther = 0;
5025: // Count number of cells in censored depth (normal cells).
5026: const Obj<typename MeshType::label_sequence>& cellsNormal = mesh->getLabelStratum("censored depth", mesh->depth());
5027: assert(!cellsNormal.isNull());
5028: const typename MeshType::label_sequence::iterator cellsNormalEnd = cellsNormal->end();
5029: oldNumCellsNormal = cellsNormal->size();
5030: for(typename MeshType::label_sequence::iterator c_iter = cellsNormal->begin(); c_iter != cellsNormalEnd; ++c_iter)
5031: newNumCellsNormal += refiner.numNewCells(*c_iter);
5033: // Count number of remaining cells (other cells).
5034: const int numSkip = oldNumCellsNormal;
5035: oldNumCellsOther = cells->size() - oldNumCellsNormal;
5036: typename MeshType::label_sequence::iterator c_iter = cells->begin();
5037: for (int i=0; i < numSkip; ++i)
5038: ++c_iter;
5039: for (; c_iter != cellsEnd; ++c_iter)
5040: newNumCellsOther += refiner.numNewCells(*c_iter);
5042: // Get number of old normal vertices.
5043: assert(!mesh->getFactory.isNull());
5044: Obj<typename Mesh::numbering_type> vNumbering = mesh->getFactory()->getNumbering(mesh, "censored depth", 0);
5045: assert(!vNumbering.isNull());
5046: oldNumVerticesNormal = vNumbering->size();
5048: // Count number of new normal vertices.
5049: int counterBegin = newNumCellsNormal + vertices->size();
5050: const point_type curNewVertex = counterBegin;
5051: for(typename MeshType::label_sequence::iterator c_iter = cellsNormal->begin(); c_iter != cellsNormalEnd; ++c_iter) {
5052: cV.clear();
5053: sieve->cone(*c_iter, cV);
5054: refiner.splitEdge(*c_iter, cV.getSize(), cV.getPoints(), curNewVertex);
5055: } // for
5056: newNumVerticesNormal = curNewVertex - counterBegin;
5058: // Count number of remaining vertices (other vertices).
5059: oldNumVerticesOther = vertices->size() - oldNumVerticessNormal;
5060: counterBegin = curNewVertex + oldNumVerticesOther;
5061: curNewVertex = counterBegin;
5062: c_iter = cells->begin();
5063: for (int i=0; i < numSkip; ++i)
5064: ++c_iter;
5065: for (; c_iter != cellsEnd; ++c_iter) {
5066: cV.clear();
5067: sieve->cone(*c_iter, cV);
5068: refiner.splitEdge(*c_iter, cV.getSize(), cV.getPoints(), curNewVertex);
5069: } // for
5070: newNumVerticesOther = curNewVertex - counterBegin;
5072: Interval<point_type> oldCellsNormalRange(0, oldNumCellsNormal);
5073: Interval<point_type> newCellsNormalRange(0, newNumCellsNormal);
5075: Interval<point_type> oldVerticesNormalRange(oldNumCellsNormal, oldNumCellsNormal+oldNumVerticesNormal);
5076: Interval<point_type> newVerticesNormalRange(newNumCellsNormal, newNumCellsNormal+newNumVerticesNormal);
5078: Interval<point_type> oldVerticesOtherRange(oldNumCellsNormal+oldNumVerticesNormal ,
5079: oldNumCellsNormal+oldNumVerticesNormal+oldNumVerticesOther);
5080: Interval<point_type> newVerticesOtherRange(newNumCellsNormal+newNumVerticesNormal ,
5081: newNumCellsNormal+newNumVerticesNormal+newNumVerticesOther);
5083: Interval<point_type> oldCellssOtherRange(oldNumCellsNormal+oldNumVerticesNormal+oldNumVerticesOther,
5084: oldNumCellsNormal+oldNumVerticesNormal+oldNumVerticesOther+oldNumCellsOther);
5085: Interval<point_type> newCellssOtherRange(newNumCellsNormal+newNumVerticesNormal+newNumVerticesOther,
5086: newNumCellsNormal+newNumVerticesNormal+newNumVerticesOther+newNumCellsOther);
5089: // Allocate chart for new sieve.
5090: const Obj<sieve_type>& newSieve = newMesh->getSieve();
5091: assert(!newSieve.isNull());
5092: newSieve->setChart(typename sieve_type::chart_type(0, newCellsOtherRange.end()));
5093: refiner.setVertexRelativeOffset(newNumCellsNormal-oldNumCellsNormal); // THIS DOES NOT WORK FOR COHESIVE CELLS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5095: // Create new sieve with correct sizes for refined cells
5097: // Start with normal cells.
5098: point_type curNewCell = newCellsNormalRange.begin();
5099: const typename Interval<point_type>::const_iterator oldCellsNormalRangeEnd = oldCellsNormalRange.end();
5100: for (typename Interval<point_type>::const_iterator c_iter=oldCellsNormalRange.begin();
5101: c_iter != oldCellsNormalRangeEnd;
5102: ++c_iter) {
5103: // Set new cone and support sizes
5104: cV.clear();
5105: sieve->cone(*c_iter, cV);
5106: const point_type *cone = cV.getPoints();
5107: const int coneSize = cV.getSize();
5108: const int newCells = refiner.numNewCells(*c_iter);
5110: for(int iCell=0; iCell < newCells; ++iCell, ++curNewCell) {
5111: const point_type *newCone;
5112: int newConeSize;
5114: newSieve->setConeSize(curNewCell, sieve->getConeSize(*c_iter));
5115: // OPTIMIZE THIS
5116: refiner.getNewCell(*c_iter, coneSize, cone, iCell, &newConeSize, &newCone);
5117: for(int v = 0; v < newConeSize; ++v)
5118: newSieve->addSupportSize(newCone[v], 1);
5119: } // for
5120: } // for
5122: // Continue with other cells.
5123: point_type curNewCell = newCellsOtherRange.begin();
5124: const typename Interval<point_type>::const_iterator oldCellsOtherRangeEnd = oldCellsOtherRange.end();
5125: for (typename Interval<point_type>::const_iterator c_iter=oldCellOtherRange.begin();
5126: c_iter != oldCellsOtherRangeEnd;
5127: ++c_iter) {
5128: // Set new cone and support sizes
5129: cV.clear();
5130: sieve->cone(*c_iter, cV);
5131: const point_type *cone = cV.getPoints();
5132: const int coneSize = cV.getSize();
5133: const int newCells = refiner.numNewCells(*c_iter);
5135: for(int iCell=0; iCell < newCells; ++iCell, ++curNewCell) {
5136: const point_type *newCone;
5137: int newConeSize;
5139: newSieve->setConeSize(curNewCell, sieve->getConeSize(*c_iter));
5140: // OPTIMIZE THIS
5141: refiner.getNewCell(*c_iter, coneSize, cone, iCell, &newConeSize, &newCone);
5142: for(int v = 0; v < newConeSize; ++v)
5143: newSieve->addSupportSize(newCone[v], 1);
5144: } // for
5145: } // for
5146: newSieve->allocate();
5147: point_type* vertex2edge = new point_type[(newNumVerticesNormal+newNumVerticesOther)*2];
5148: typename Refiner::edge_map_type& edge2vertex = refiner.getEdgeToVertex();
5151: // Create refined normal cells in new sieve.
5152: curNewCell = newCellsNormalRange.begin();
5153: for (typename Interval<point_type>::const_iterator c_iter=oldCellsNormalRange.begin();
5154: c_iter != oldCellsNormalRangeEnd;
5155: ++c_iter) {
5156: cV.clear();
5157: sieve->cone(*c_iter, cV);
5158: const point_type *cone = cV.getPoints();
5159: const int coneSize = cV.getSize();
5160: const int newCells = refiner.numNewCells(*c_iter);
5162: for (int iCell=0; iCell < newCells; ++iCell, ++curNewCell) {
5163: const point_type *newCone;
5164: int newConeSize;
5166: refiner.getNewCell(*c_iter, coneSize, cone, iCell, &newConeSize, &newCone);
5167: newSieve->setCone(newCone, curNewCell);
5168: } // for
5170: refiner.getNeighboringVertices(*c_iter, coneSize, cone, firstNewVertex, vertex2edge);
5171: } // for
5173: // Create refined other cells in new sieve.
5174: curNewCell = newCellsOtherRange.begin();
5175: for (typename Interval<point_type>::const_iterator c_iter=oldCellsOtherRange.begin();
5176: c_iter != oldCellsOtherRangeEnd;
5177: ++c_iter) {
5178: cV.clear();
5179: sieve->cone(*c_iter, cV);
5180: const point_type *cone = cV.getPoints();
5181: const int coneSize = cV.getSize();
5182: const int newCells = refiner.numNewCells(*c_iter);
5184: for (int iCell=0; iCell < newCells; ++iCell, ++curNewCell) {
5185: const point_type *newCone;
5186: int newConeSize;
5188: refiner.getNewCell(*c_iter, coneSize, cone, iCell, &newConeSize, &newCone);
5189: newSieve->setCone(newCone, curNewCell);
5190: } // for
5192: // FIX THIS!!! LAGRANGE VERTICES MUST BE AFTER ALL OTHER VERTICES (INCLUDING OLD LAGRANGE VERTICES)
5193: refiner.getNeighboringVertices(*c_iter, coneSize, cone, firstNewVertex, vertex2edge);
5194: } // for
5195: newSieve->symmetrize();
5197: // Set coordinates in refined mesh.
5198: const Obj<typename MeshType::real_section_type>& coordinates = mesh->getRealSection("coordinates");
5199: assert(!coordinates.isNull());
5200: const Obj<typename MeshType::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");
5201: assert(!newCoordinates.isNull());
5203: const typename MeshType::label_sequence::const_iterator verticesEnd = vertices->end();
5204: assert(vertices->size() > 0);
5205: const int spaceDim = coordinates->getFiberDimension(*vertices->begin());
5206: assert(spaceDim > 0);
5207: newCoordinates->setChart(typename sieve_type::chart_type(newNumCellsNormal, newNumCellsNormal+newNumVertices));
5209: for (int iVertex=0, offset=newNumCellsNormal; iVertex < newNumVertices; ++iVertex) {
5210: const point_type vNew = iVertex + offset;
5211: newCoordinates->setFiberDimension(vNew, spaceDim);
5212: } // for
5213: newCoordinates->allocatePoint();
5215: for (int iVertex=0, oldOffset=oldNumCellsNormal, newOffset=newNumCellsNormal; iVertex < oldNumVertices; ++iVertex) {
5216: const point_type vOld = iVertex + oldOffset;
5217: const point_type vNew = iVertex + newOffset;
5218: newCoordinates->updatePoint(vNew, coordinates->restrictPoint(vOld));
5219: } // for
5220: for(int v=0, iVertex=oldNumVertices, newOffset=newNumCellsNormal; iVertex < newNumVertices; ++v, ++iVertex) {
5221: const point_type vNew = iVertex + newOffset;
5222: const point_type endpointA = vertex2edge[v*2+0];
5223: const point_type endpointB = vertex2edge[v*2+1];
5224: std::cout << "Setting coordinates of vertex " << vNew << " between vertices "
5225: << endpointA << " and " << endpointB << std::endl;
5226: const real *coordsA = coordinates->restrictPoint(endpointA);
5227: real coords[3];
5229: for(int d = 0; d < 3; ++d)
5230: coords[d] = coordsA[d];
5231: const real *coordsB = coordinates->restrictPoint(endpointB);
5232: for(int d = 0; d < 3; ++d) {
5233: coords[d] += coordsB[d];
5234: coords[d] *= 0.5;
5235: } // for
5236: newCoordinates->updatePoint(vNew, coords);
5237: } // for
5238: delete [] vertex2edge;
5239: // Fast stratification
5240: const Obj<typename MeshType::label_type>& height = newMesh->createLabel("height");
5241: const Obj<typename MeshType::label_type>& depth = newMesh->createLabel("depth");
5243: for (int iCell=0; iCell < newNumCellsNormal; ++iCell) {
5244: const point_type cNew = iCell;
5245: height->setCone(0, cNew);
5246: depth->setCone(1, cNew);
5247: } // for
5248: for (int iCell=newNumCellsNormal, offset=newNumVertices; iCell < newNumCells; ++iCell) {
5249: const point_type cNew = iCell + offset;
5250: height->setCone(0, cNew);
5251: depth->setCone(1, cNew);
5252: } // for
5253: for (int iVertex=0, newOffset=newNumCellsNormal; iVertex < newNumVertices; ++iVertex) {
5254: const point_type vNew = iVertex + newOffset;
5255: height->setCone(1, vNew);
5256: depth->setCone(0, vNew);
5257: } // for
5258: newMesh->setHeight(1);
5259: newMesh->setDepth(1);
5261: } else {
5262: int counterBegin = 0;
5264: int oldNumCells = 0;
5265: int oldNumVertices = 0;
5267: int newNumCells = 0;
5268: int newNumVertices = 0;
5270: // Count number of cells.
5271: oldNumCells = cells->size();
5272: for (typename MeshType::label_sequence::iterator c_iter = cells->begin(); c_iter != cellsEnd; ++c_iter)
5273: newNumCells += refiner.numNewCells(*c_iter);
5275: // Count number of vertices (normal vertices).
5276: oldNumVertices = vertices->size();
5277: int counterBegin = newNumCells + oldNumVertices;
5278: const point_type curNewVertex = counterBegin;
5279: for(typename MeshType::label_sequence::iterator c_iter = cells->begin(); c_iter != cellsEnd; ++c_iter) {
5280: cV.clear();
5281: sieve->cone(*c_iter, cV);
5282: refiner.splitEdge(*c_iter, cV.getSize(), cV.getPoints(), curNewVertex);
5283: } // for
5284: newNumVertices = curNewVertex - counterBegin;
5286: Interval<point_type> oldCellsRange(0, oldNumCells);
5287: Interval<point_type> newCellsRange(0, newNumCells);
5289: Interval<point_type> oldVerticesRange(oldNumCells, oldNumCells+oldNumVertices);
5290: Interval<point_type> newVerticesRange(newNumCells, newNumCells+newNumVertices);
5292: // Allocate chart for new sieve.
5293: const Obj<sieve_type>& newSieve = newMesh->getSieve();
5294: assert(!newSieve.isNull());
5295: newSieve->setChart(typename sieve_type::chart_type(0, newCellsOtherRange.end()));
5296: refiner.setVertexRelativeOffset(newNumCells-oldNumCells);
5298: // Create new sieve with correct sizes for refined cells
5300: // Start with normal cells.
5301: point_type curNewCell = newCellsRange.begin();
5302: const typename Interval<point_type>::const_iterator oldCellsRangeEnd = oldCellsRange.end();
5303: for (typename Interval<point_type>::const_iterator c_iter=oldCellsRange.begin();
5304: c_iter != oldCellsRangeEnd;
5305: ++c_iter) {
5306: // Set new cone and support sizes
5307: cV.clear();
5308: sieve->cone(*c_iter, cV);
5309: const point_type *cone = cV.getPoints();
5310: const int coneSize = cV.getSize();
5311: const int newCells = refiner.numNewCells(*c_iter);
5313: for(int iCell=0; iCell < newCells; ++iCell, ++curNewCell) {
5314: const point_type *newCone;
5315: int newConeSize;
5317: newSieve->setConeSize(curNewCell, sieve->getConeSize(*c_iter));
5318: // OPTIMIZE THIS
5319: refiner.getNewCell(*c_iter, coneSize, cone, iCell, &newConeSize, &newCone);
5320: for(int v = 0; v < newConeSize; ++v)
5321: newSieve->addSupportSize(newCone[v], 1);
5322: } // for
5323: } // for
5325: // Create refined normal cells in new sieve.
5326: curNewCell = newCellsRange.begin();
5327: for (typename Interval<point_type>::const_iterator c_iter=oldCellsRange.begin();
5328: c_iter != oldCellsRangeEnd;
5329: ++c_iter) {
5330: cV.clear();
5331: sieve->cone(*c_iter, cV);
5332: const point_type *cone = cV.getPoints();
5333: const int coneSize = cV.getSize();
5334: const int newCells = refiner.numNewCells(*c_iter);
5336: for (int iCell=0; iCell < newCells; ++iCell, ++curNewCell) {
5337: const point_type *newCone;
5338: int newConeSize;
5340: refiner.getNewCell(*c_iter, coneSize, cone, iCell, &newConeSize, &newCone);
5341: newSieve->setCone(newCone, curNewCell);
5342: } // for
5344: refiner.getNeighboringVertices(*c_iter, coneSize, cone, firstNewVertex, vertex2edge);
5345: } // for
5346: newSieve->symmetrize();
5348: // Set coordinates in refined mesh.
5349: const Obj<typename MeshType::real_section_type>& coordinates = mesh->getRealSection("coordinates");
5350: assert(!coordinates.isNull());
5351: const Obj<typename MeshType::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");
5352: assert(!newCoordinates.isNull());
5354: const typename MeshType::label_sequence::const_iterator verticesEnd = vertices->end();
5355: assert(vertices->size() > 0);
5356: const int spaceDim = coordinates->getFiberDimension(*vertices->begin());
5357: assert(spaceDim > 0);
5358: newCoordinates->setChart(typename sieve_type::chart_type(newVerticesRange.begin(), newVerticesRange.end()));
5360: const typename Interval<point_type>::const_iterator newVerticesRangeEnd = newVerticesRange.end();
5361: for (typename Interval<point_type>::const_iterator v_iter=newVerticesRange.begin(); v_iter != newVerticesRangeEnd; ++v_iter)
5362: newCoordinates->setFiberDimension(v_iter, spaceDim);
5363: newCoordinates->allocatePoint();
5365: const typename Interval<point_type>::const_iterator oldVerticesRangeEnd = oldVerticesRange.end();
5366: for (typename Interval<point_type>::const_iterator vOld_iter=oldVerticesRange.begin(), vNew_iter=newVerticesRange.begin(); vOld_iter != oldVerticesRangeEnd; ++vOld_iter)
5367: newCoordinates->updatePoint(*vNew_iter, coordinates->restrictPoint(*vOld_iter));
5368: for(int v=0, iVertex=oldNumVertices; iVertex < newNumVertices; ++v, ++iVertex) {
5369: const point_type vNew = newVerticesRange.begin() + iVertex;
5370: const point_type endpointA = vertex2edge[v*2+0];
5371: const point_type endpointB = vertex2edge[v*2+1];
5372: std::cout << "Setting coordinates of vertex " << vNew << " between vertices "
5373: << endpointA << " and " << endpointB << std::endl;
5374: const real *coordsA = coordinates->restrictPoint(endpointA);
5375: real coords[3];
5377: for(int d = 0; d < 3; ++d)
5378: coords[d] = coordsA[d];
5379: const real *coordsB = coordinates->restrictPoint(endpointB);
5380: for(int d = 0; d < 3; ++d) {
5381: coords[d] += coordsB[d];
5382: coords[d] *= 0.5;
5383: } // for
5384: newCoordinates->updatePoint(vNew, coords);
5385: } // for
5386: delete [] vertex2edge;
5388: // Fast stratification
5389: const Obj<typename MeshType::label_type>& height = newMesh->createLabel("height");
5390: const Obj<typename MeshType::label_type>& depth = newMesh->createLabel("depth");
5391: for (int iCell=0; iCell < newNumCells; ++iCell) {
5392: const point_type cNew = iCell;
5393: height->setCone(0, cNew);
5394: depth->setCone(1, cNew);
5395: } // for
5396: for (int iVertex=0, newOffset=newNumCellsNormal; iVertex < newNumVertices; ++iVertex) {
5397: const point_type vNew = iVertex + newOffset;
5398: height->setCone(1, vNew);
5399: depth->setCone(0, vNew);
5400: } // for
5401: newMesh->setHeight(1);
5402: newMesh->setDepth(1);
5403: } // if/else
5405: // Exchange new boundary vertices
5406: // We can convert endpoints, and then just match to new vertex on this side
5407: // 1) Create the overlap of edges which are vertex pairs (do not need for interpolated meshes)
5408: // 2) Create a section of overlap edge --> new vertex (this will generalize to other split points in interpolated meshes)
5409: // 3) Copy across new overlap
5410: // 4) Fuse matches new vertex pairs and inserts them into the old overlap
5413: // Create the parallel overlap
5414: int *oldNumCellsP = new int[mesh->commSize()];
5415: int *newNumCellsP = new int[newMesh->commSize()];
5416: int ierr;
5418: MPI_Allgather((void *) &oldNumCells, 1, MPI_INT, oldNumCellsP, 1, MPI_INT, mesh->comm());CHKERRXX(ierr);
5419: MPI_Allgather((void *) &newNumCells, 1, MPI_INT, newNumCellsP, 1, MPI_INT, newMesh->comm());CHKERRXX(ierr);
5420: Obj<typename MeshType::send_overlap_type> newSendOverlap = newMesh->getSendOverlap();
5421: Obj<typename MeshType::recv_overlap_type> newRecvOverlap = newMesh->getRecvOverlap();
5422: const Obj<typename MeshType::send_overlap_type>& sendOverlap = mesh->getSendOverlap();
5423: const Obj<typename MeshType::recv_overlap_type>& recvOverlap = mesh->getRecvOverlap();
5424: Obj<typename MeshType::send_overlap_type::traits::capSequence> sendPoints = sendOverlap->cap();
5425: const typename MeshType::send_overlap_type::source_type localOffset = newNumCellsP[newMesh->commRank()] - oldNumCellsP[mesh->commRank()];
5427: for(typename MeshType::send_overlap_type::traits::capSequence::iterator p_iter = sendPoints->begin(); p_iter != sendPoints->end(); ++p_iter) {
5428: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& ranks = sendOverlap->support(*p_iter);
5429: const typename MeshType::send_overlap_type::source_type& localPoint = *p_iter;
5431: for(typename MeshType::send_overlap_type::traits::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
5432: const int rank = *r_iter;
5433: const typename MeshType::send_overlap_type::source_type& remotePoint = r_iter.color();
5434: const typename MeshType::send_overlap_type::source_type remoteOffset = newNumCellsP[rank] - oldNumCellsP[rank];
5436: newSendOverlap->addArrow(localPoint+localOffset, rank, remotePoint+remoteOffset);
5437: }
5438: }
5439: Obj<typename MeshType::recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();
5441: for(typename MeshType::recv_overlap_type::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
5442: const Obj<typename MeshType::recv_overlap_type::traits::coneSequence>& ranks = recvOverlap->cone(*p_iter);
5443: const typename MeshType::recv_overlap_type::target_type& localPoint = *p_iter;
5445: for(typename MeshType::recv_overlap_type::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
5446: const int rank = *r_iter;
5447: const typename MeshType::recv_overlap_type::target_type& remotePoint = r_iter.color();
5448: const typename MeshType::recv_overlap_type::target_type remoteOffset = newNumCellsP[rank] - oldNumCellsP[rank];
5450: newRecvOverlap->addArrow(rank, localPoint+localOffset, remotePoint+remoteOffset);
5451: }
5452: }
5453: newMesh->setCalculatedOverlap(true);
5454: delete [] oldNumCellsP;
5455: delete [] newNumCellsP;
5456: // Check edges in edge2vertex for both endpoints sent to same process
5457: // Put it in section with point being the lowest numbered vertex and value (other endpoint, new vertex)
5458: Obj<ALE::Section<point_type, edge_type> > newVerticesSection = new ALE::Section<point_type, edge_type>(mesh->comm());
5459: std::map<edge_type, std::vector<int> > bdedge2rank;
5461: for(typename std::map<edge_type, point_type>::const_iterator e_iter = edge2vertex.begin(); e_iter != edge2vertex.end(); ++e_iter) {
5462: const point_type left = e_iter->first.first;
5463: const point_type right = e_iter->first.second;
5465: if (sendOverlap->capContains(left) && sendOverlap->capContains(right)) {
5466: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& leftRanksSeq = sendOverlap->support(left);
5467: std::set<int> leftRanks(leftRanksSeq->begin(), leftRanksSeq->end());
5468: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& rightRanksSeq = sendOverlap->support(right);
5469: std::set<int> rightRanks(rightRanksSeq->begin(), rightRanksSeq->end());
5470: std::set<int> ranks;
5471: std::set_intersection(leftRanks.begin(), leftRanks.end(), rightRanks->begin(), rightRanks->end(),
5472: std::insert_iterator<std::list<int> >(ranks, ranks.begin()));
5474: if(ranks.size()) {
5475: newVerticesSection->addFiberDimension(std::min(e_iter->first.first, e_iter->first.second)+localOffset, 1);
5476: for(typename std::list<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
5477: bdedge2rank[e_iter->first].push_back(*r_iter);
5478: }
5479: }
5480: }
5481: }
5482: newVerticesSection->allocatePoint();
5483: const typename ALE::Section<point_type, edge_type>::chart_type& chart = newVerticesSection->getChart();
5485: for(typename ALE::Section<point_type, edge_type>::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
5486: typedef typename ALE::Section<point_type, edge_type>::value_type value_type;
5487: const point_type p = *c_iter;
5488: const int dim = newVerticesSection->getFiberDimension(p);
5489: int v = 0;
5490: value_type *values = new value_type[dim];
5492: for(typename std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end() && v < dim; ++e_iter) {
5493: if (std::min(e_iter->first.first, e_iter->first.second)+localOffset == p) {
5494: values[v++] = edge_type(std::max(e_iter->first.first, e_iter->first.second)+localOffset, edge2vertex[e_iter->first]);
5495: }
5496: }
5497: newVerticesSection->updatePoint(p, values);
5498: delete [] values;
5499: }
5500: // Copy across overlap
5501: typedef ALE::Pair<int, point_type> overlap_point_type;
5502: Obj<ALE::Section<overlap_point_type, edge_type> > overlapVertices = new ALE::Section<overlap_point_type, edge_type>(mesh->comm());
5504: ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
5505: // Merge by translating edge to local points, finding edge in edge2vertex, and adding (local new vetex, remote new vertex) to overlap
5506: for(typename std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end(); ++e_iter) {
5507: const point_type localPoint = edge2vertex[e_iter->first];
5509: for(typename std::vector<int>::const_iterator r_iter = e_iter->second.begin(); r_iter != e_iter->second.end(); ++r_iter) {
5510: point_type remoteLeft = -1, remoteRight = -1;
5511: const int rank = *r_iter;
5513: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& leftRanks = newSendOverlap->support(e_iter->first.first+localOffset);
5514: for(typename MeshType::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks->begin(); lr_iter != leftRanks->end(); ++lr_iter) {
5515: if (rank == *lr_iter) {
5516: remoteLeft = lr_iter.color();
5517: break;
5518: }
5519: }
5520: const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& rightRanks = newSendOverlap->support(e_iter->first.second+localOffset);
5521: for(typename MeshType::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks->begin(); rr_iter != rightRanks->end(); ++rr_iter) {
5522: if (rank == *rr_iter) {
5523: remoteRight = rr_iter.color();
5524: break;
5525: }
5526: }
5527: const point_type remoteMin = std::min(remoteLeft, remoteRight);
5528: const point_type remoteMax = std::max(remoteLeft, remoteRight);
5529: const int remoteSize = overlapVertices->getFiberDimension(overlap_point_type(rank, remoteMin));
5530: const edge_type *remoteVals = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
5531: point_type remotePoint = -1;
5533: for(int d = 0; d < remoteSize; ++d) {
5534: if (remoteVals[d].first == remoteMax) {
5535: remotePoint = remoteVals[d].second;
5536: break;
5537: }
5538: }
5539: newSendOverlap->addArrow(localPoint, rank, remotePoint);
5540: newRecvOverlap->addArrow(rank, localPoint, remotePoint);
5541: }
5542: }
5543: };
5544: #endif
5545: };
5547: class MeshSerializer {
5548: public:
5549: template<typename Mesh>
5550: static void writeMesh(const std::string& filename, Mesh& mesh) {
5551: std::ofstream fs;
5553: if (mesh.commRank() == 0) {
5554: fs.open(filename.c_str());
5555: }
5556: writeMesh(fs, mesh);
5557: if (mesh.commRank() == 0) {
5558: fs.close();
5559: }
5560: }
5561: template<typename Mesh>
5562: static void writeMesh(std::ofstream& fs, Mesh& mesh) {
5563: ISieveSerializer::writeSieve(fs, *mesh.getSieve());
5564: // Write labels
5565: const typename Mesh::labels_type& labels = mesh.getLabels();
5567: if (!mesh.commRank()) {fs << labels.size() << std::endl;}
5568: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
5569: if (!mesh.commRank()) {fs << l_iter->first << std::endl;}
5570: LabelSifterSerializer::writeLabel(fs, *l_iter->second);
5571: }
5572: // Write sections
5573: Obj<std::set<std::string> > realNames = mesh.getRealSections();
5575: if (!mesh.commRank()) {fs << realNames->size() << std::endl;}
5576: for(std::set<std::string>::const_iterator n_iter = realNames->begin(); n_iter != realNames->end(); ++n_iter) {
5577: if (!mesh.commRank()) {fs << *n_iter << std::endl;}
5578: SectionSerializer::writeSection(fs, *mesh.getRealSection(*n_iter));
5579: }
5580: Obj<std::set<std::string> > intNames = mesh.getIntSections();
5582: if (!mesh.commRank()) {fs << intNames->size() << std::endl;}
5583: for(std::set<std::string>::const_iterator n_iter = intNames->begin(); n_iter != intNames->end(); ++n_iter) {
5584: if (!mesh.commRank()) {fs << *n_iter << std::endl;}
5585: SectionSerializer::writeSection(fs, *mesh.getIntSection(*n_iter));
5586: }
5587: // Write overlap
5588: #ifdef USE_NEW_OVERLAP
5589: PETSc::OverlapSerializer::writeOverlap(fs, *mesh.getSendOverlap());
5590: PETSc::OverlapSerializer::writeOverlap(fs, *mesh.getRecvOverlap());
5591: #else
5592: SifterSerializer::writeSifter(fs, *mesh.getSendOverlap());
5593: SifterSerializer::writeSifter(fs, *mesh.getRecvOverlap());
5594: #endif
5595: // Write distribution overlap
5596: // Write renumbering
5597: }
5598: template<typename Mesh>
5599: static void loadMesh(const std::string& filename, Mesh& mesh) {
5600: std::ifstream fs;
5602: if (mesh.commRank() == 0) {
5603: fs.open(filename.c_str());
5604: }
5605: loadMesh(fs, mesh);
5606: if (mesh.commRank() == 0) {
5607: fs.close();
5608: }
5609: }
5610: template<typename Mesh>
5611: static void loadMesh(std::ifstream& fs, Mesh& mesh) {
5612: ALE::Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh.comm(), mesh.debug());
5613: PetscErrorCode ierr;
5615: ISieveSerializer::loadSieve(fs, *sieve);
5616: mesh.setSieve(sieve);
5617: // Load labels
5618: int numLabels;
5620: if (!mesh.commRank()) {fs >> numLabels;}
5621: MPI_Bcast(&numLabels, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5622: for(int l = 0; l < numLabels; ++l) {
5623: ALE::Obj<typename Mesh::label_type> label = new typename Mesh::label_type(mesh.comm(), mesh.debug());
5624: std::string name;
5625: int len;
5627: if (!mesh.commRank()) {
5628: fs >> name;
5629: len = name.size();
5630: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5631: MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
5632: } else {
5633: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5634: char *n = new char[len+1];
5635: MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
5636: name = n;
5637: delete [] n;
5638: }
5639: LabelSifterSerializer::loadLabel(fs, *label);
5640: mesh.setLabel(name, label);
5641: }
5642: // Load sections
5643: int numRealSections;
5645: if (!mesh.commRank()) {fs >> numRealSections;}
5646: MPI_Bcast(&numRealSections, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5647: for(int s = 0; s < numRealSections; ++s) {
5648: ALE::Obj<typename Mesh::real_section_type> section = new typename Mesh::real_section_type(mesh.comm(), mesh.debug());
5649: std::string name;
5650: int len;
5652: if (!mesh.commRank()) {
5653: fs >> name;
5654: len = name.size();
5655: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5656: MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
5657: } else {
5658: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5659: char *n = new char[len+1];
5660: MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
5661: name = n;
5662: delete [] n;
5663: }
5664: SectionSerializer::loadSection(fs, *section);
5665: mesh.setRealSection(name, section);
5666: }
5667: int numIntSections;
5669: if (!mesh.commRank()) {fs >> numIntSections;}
5670: MPI_Bcast(&numIntSections, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5671: for(int s = 0; s < numIntSections; ++s) {
5672: ALE::Obj<typename Mesh::int_section_type> section = new typename Mesh::int_section_type(mesh.comm(), mesh.debug());
5673: std::string name;
5674: int len;
5676: if (!mesh.commRank()) {
5677: fs >> name;
5678: len = name.size();
5679: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5680: MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
5681: } else {
5682: MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
5683: char *n = new char[len+1];
5684: MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
5685: name = n;
5686: delete [] n;
5687: }
5688: SectionSerializer::loadSection(fs, *section);
5689: mesh.setIntSection(name, section);
5690: }
5691: // Load overlap
5692: #ifdef USE_NEW_OVERLAP
5693: PETSc::OverlapSerializer::loadOverlap(fs, *mesh.getSendOverlap());
5694: PETSc::OverlapSerializer::loadOverlap(fs, *mesh.getRecvOverlap());
5695: #else
5696: SifterSerializer::loadSifter(fs, *mesh.getSendOverlap());
5697: SifterSerializer::loadSifter(fs, *mesh.getRecvOverlap());
5698: #endif
5699: // Load distribution overlap
5700: // Load renumbering
5701: }
5702: };
5703: } // namespace ALE
5704: #endif