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