Actual source code: Mesh.hh

  1: #ifndef included_ALE_Mesh_hh
  2: #define included_ALE_Mesh_hh

  4: #include <list>
  5: #include <valarray>

  7: #ifndef  included_ALE_Numbering_hh
  8: #include <Numbering.hh>
  9: #endif

 11: #ifndef  included_ALE_INumbering_hh
 12: #include <INumbering.hh>
 13: #endif

 15: #ifndef  included_ALE_Field_hh
 16: #include <Field.hh>
 17: #endif

 19: #ifndef  included_ALE_IField_hh
 20: #include <IField.hh>
 21: #endif

 23: #ifndef  included_ALE_SieveBuilder_hh
 24: #include <SieveBuilder.hh>
 25: #endif

 27: #ifndef  included_ALE_LabelSifter_hh
 28: #include <LabelSifter.hh>
 29: #endif

 31: #ifndef  included_ALE_Partitioner_hh
 32: #include <Partitioner.hh>
 33: #endif

 35: namespace ALE {
 36:   class indexSet : public std::valarray<int> {
 37:   public:
 38:     inline bool
 39:     operator<(const indexSet& __x) {
 40:       if (__x.size() != this->size()) return __x.size() < this->size();
 41:       for(unsigned int i = 0; i < __x.size(); ++i) {
 42:         if (__x[i] == (*this)[i]) continue;
 43:           return __x[i] < (*this)[i];
 44:       }
 45:       return false;
 46:     }
 47:   };
 48:   inline bool
 49:   operator<(const indexSet& __x, const indexSet& __y) {
 50:     if (__x.size() != __y.size()) return __x.size() < __y.size();
 51:     for(unsigned int i = 0; i < __x.size(); ++i) {
 52:       if (__x[i] == __y[i]) continue;
 53:       return __x[i] < __y[i];
 54:     }
 55:     return false;
 56:   };
 57:   inline bool
 58:   operator<=(const indexSet& __x, const indexSet& __y) {
 59:     if (__x.size() != __y.size()) return __x.size() < __y.size();
 60:     for(unsigned int i = 0; i < __x.size(); ++i) {
 61:       if (__x[i] == __y[i]) continue;
 62:       return __x[i] < __y[i];
 63:     }
 64:     return true;
 65:   };
 66:   inline bool
 67:   operator==(const indexSet& __x, const indexSet& __y) {
 68:     if (__x.size() != __y.size()) return false;
 69:     for(unsigned int i = 0; i < __x.size(); ++i) {
 70:       if (__x[i] != __y[i]) return false;
 71:     }
 72:     return true;
 73:   };
 74:   inline bool
 75:   operator!=(const indexSet& __x, const indexSet& __y) {
 76:     if (__x.size() != __y.size()) return true;
 77:     for(unsigned int i = 0; i < __x.size(); ++i) {
 78:       if (__x[i] != __y[i]) return true;
 79:     }
 80:     return false;
 81:   };

 83:   template<typename Sieve_,
 84:            typename RealSection_  = Section<typename Sieve_::point_type, double>,
 85:            typename IntSection_   = Section<typename Sieve_::point_type, int>,
 86:            typename ArrowSection_ = UniformSection<MinimalArrow<typename Sieve_::point_type, typename Sieve_::point_type>, int> >
 87:   class Bundle : public ALE::ParallelObject {
 88:   public:
 89:     typedef Sieve_                                                    sieve_type;
 90:     typedef RealSection_                                              real_section_type;
 91:     typedef IntSection_                                               int_section_type;
 92:     typedef ArrowSection_                                             arrow_section_type;
 93:     typedef Bundle<Sieve_,RealSection_,IntSection_,ArrowSection_>     this_type;
 94:     typedef typename sieve_type::point_type                           point_type;
 95:     typedef malloc_allocator<point_type>                              alloc_type;
 96:     typedef typename ALE::LabelSifter<int, point_type>                label_type;
 97:     typedef typename std::map<const std::string, Obj<label_type> >    labels_type;
 98:     typedef typename label_type::supportSequence                      label_sequence;
 99:     typedef std::map<std::string, Obj<arrow_section_type> >           arrow_sections_type;
100:     typedef std::map<std::string, Obj<real_section_type> >            real_sections_type;
101:     typedef std::map<std::string, Obj<int_section_type> >             int_sections_type;
102:     typedef ALE::Point                                                index_type;
103:     typedef std::pair<index_type, int>                                oIndex_type;
104:     typedef std::vector<oIndex_type>                                  oIndexArray;
105:     typedef std::pair<int *, int>                                     indices_type;
106:     typedef NumberingFactory<this_type>                               MeshNumberingFactory;
107:     typedef typename ALE::Partitioner<>::part_type                    rank_type;
108:     typedef typename ALE::Sifter<point_type,rank_type,point_type>     send_overlap_type;
109:     typedef typename ALE::Sifter<rank_type,point_type,point_type>     recv_overlap_type;
110:     typedef typename MeshNumberingFactory::numbering_type             numbering_type;
111:     typedef typename MeshNumberingFactory::order_type                 order_type;
112:     typedef std::map<point_type, point_type>                          renumbering_type;
113:     typedef typename ALE::SieveAlg<this_type>                         sieve_alg_type;
114:     typedef typename sieve_alg_type::coneArray                        coneArray;
115:     typedef typename sieve_alg_type::orientedConeArray                oConeArray;
116:     typedef typename sieve_alg_type::supportArray                     supportArray;
117:   protected:
118:     Obj<sieve_type>       _sieve;
119:     labels_type           _labels;
120:     int                   _maxHeight;
121:     int                   _maxDepth;
122:     arrow_sections_type   _arrowSections;
123:     real_sections_type    _realSections;
124:     int_sections_type     _intSections;
125:     Obj<oIndexArray>      _indexArray;
126:     Obj<MeshNumberingFactory> _factory;
127:     bool                   _calculatedOverlap;
128:     Obj<send_overlap_type> _sendOverlap;
129:     Obj<recv_overlap_type> _recvOverlap;
130:     Obj<send_overlap_type> _distSendOverlap;
131:     Obj<recv_overlap_type> _distRecvOverlap;
132:     renumbering_type       _renumbering; // Maps global points to local points
133:     // Work space
134:     Obj<std::set<point_type> > _modifiedPoints;
135:   public:
136:     Bundle(MPI_Comm comm, int debug = 0) : ALE::ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1) {
137:       this->_indexArray        = new oIndexArray();
138:       this->_modifiedPoints    = new std::set<point_type>();
139:       this->_factory           = MeshNumberingFactory::singleton(this->comm(), this->debug());
140:       this->_calculatedOverlap = false;
141:       this->_sendOverlap       = new send_overlap_type(comm, debug);
142:       this->_recvOverlap       = new recv_overlap_type(comm, debug);
143:     };
144:     Bundle(const Obj<sieve_type>& sieve) : ALE::ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _maxHeight(-1), _maxDepth(-1) {
145:       this->_indexArray        = new oIndexArray();
146:       this->_modifiedPoints    = new std::set<point_type>();
147:       this->_factory           = MeshNumberingFactory::singleton(this->comm(), this->debug());
148:       this->_calculatedOverlap = false;
149:       this->_sendOverlap       = new send_overlap_type(comm, debug);
150:       this->_recvOverlap       = new recv_overlap_type(comm, debug);
151:     };
152:     virtual ~Bundle() {};
153:   public: // Verifiers
154:     bool hasLabel(const std::string& name) {
155:       if (this->_labels.find(name) != this->_labels.end()) {
156:         return true;
157:       }
158:       return false;
159:     };
160:     void checkLabel(const std::string& name) {
161:       if (!this->hasLabel(name)) {
162:         ostringstream msg;
163:         msg << "Invalid label name: " << name << std::endl;
164:         throw ALE::Exception(msg.str().c_str());
165:       }
166:     };
167:   public: // Accessors
168:     const Obj<sieve_type>& getSieve() const {return this->_sieve;};
169:     void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
170:     bool hasArrowSection(const std::string& name) const {
171:       return this->_arrowSections.find(name) != this->_arrowSections.end();
172:     };
173:     const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
174:       if (!this->hasArrowSection(name)) {
175:         Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());

177:         section->setName(name);
178:         if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
179:         this->_arrowSections[name] = section;
180:       }
181:       return this->_arrowSections[name];
182:     };
183:     void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
184:       this->_arrowSections[name] = section;
185:     };
186:     Obj<std::set<std::string> > getArrowSections() const {
187:       Obj<std::set<std::string> > names = std::set<std::string>();

189:       for(typename arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
190:         names->insert(s_iter->first);
191:       }
192:       return names;
193:     };
194:     bool hasRealSection(const std::string& name) const {
195:       return this->_realSections.find(name) != this->_realSections.end();
196:     };
197:     const Obj<real_section_type>& getRealSection(const std::string& name) {
198:       if (!this->hasRealSection(name)) {
199:         Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());

201:         section->setName(name);
202:         if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
203:         this->_realSections[name] = section;
204:       }
205:       return this->_realSections[name];
206:     };
207:     void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
208:       this->_realSections[name] = section;
209:     };
210:     Obj<std::set<std::string> > getRealSections() const {
211:       Obj<std::set<std::string> > names = std::set<std::string>();

213:       for(typename real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
214:         names->insert(s_iter->first);
215:       }
216:       return names;
217:     };
218:     bool hasIntSection(const std::string& name) const {
219:       return this->_intSections.find(name) != this->_intSections.end();
220:     };
221:     const Obj<int_section_type>& getIntSection(const std::string& name) {
222:       if (!this->hasIntSection(name)) {
223:         Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());

225:         section->setName(name);
226:         if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
227:         this->_intSections[name] = section;
228:       }
229:       return this->_intSections[name];
230:     };
231:     void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
232:       this->_intSections[name] = section;
233:     };
234:     Obj<std::set<std::string> > getIntSections() const {
235:       Obj<std::set<std::string> > names = std::set<std::string>();

237:       for(typename int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
238:         names->insert(s_iter->first);
239:       }
240:       return names;
241:     };
242:     const Obj<MeshNumberingFactory>& getFactory() const {return this->_factory;};
243:     bool getCalculatedOverlap() const {return this->_calculatedOverlap;};
244:     void setCalculatedOverlap(const bool calc) {this->_calculatedOverlap = calc;};
245:     const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
246:     void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
247:     const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
248:     void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
249:     const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
250:     void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
251:     const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
252:     void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
253:     renumbering_type& getRenumbering() {return this->_renumbering;};
254:   public: // Labels
255:     int getValue (const Obj<label_type>& label, const point_type& point, const int defValue = 0) {
256:       const Obj<typename label_type::coneSequence>& cone = label->cone(point);

258:       if (cone->size() == 0) return defValue;
259:       return *cone->begin();
260:     };
261:     void setValue(const Obj<label_type>& label, const point_type& point, const int value) {
262:       label->setCone(value, point);
263:     };
264:     template<typename InputPoints>
265:     int getMaxValue (const Obj<label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
266:       int maxValue = defValue;

268:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
269:         maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
270:       }
271:       return maxValue;
272:     };
273:     const Obj<label_type>& createLabel(const std::string& name) {
274:       this->_labels[name] = new label_type(this->comm(), this->debug());
275:       return this->_labels[name];
276:     };
277:     const Obj<label_type>& getLabel(const std::string& name) {
278:       this->checkLabel(name);
279:       return this->_labels[name];
280:     };
281:     void setLabel(const std::string& name, const Obj<label_type>& label) {
282:       this->_labels[name] = label;
283:     };
284:     const labels_type& getLabels() {
285:       return this->_labels;
286:     };
287:     virtual const Obj<label_sequence>& getLabelStratum(const std::string& name, int value) {
288:       this->checkLabel(name);
289:       return this->_labels[name]->support(value);
290:     };
291:   public: // Stratification
292:     template<class InputPoints>
293:     void computeHeight(const Obj<label_type>& height, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxHeight) {
294:       this->_modifiedPoints->clear();

296:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
297:         // Compute the max height of the points in the support of p, and add 1
298:         int h0 = this->getValue(height, *p_iter, -1);
299:         int h1 = this->getMaxValue(height, sieve->support(*p_iter), -1) + 1;

301:         if(h1 != h0) {
302:           this->setValue(height, *p_iter, h1);
303:           if (h1 > maxHeight) maxHeight = h1;
304:           this->_modifiedPoints->insert(*p_iter);
305:         }
306:       }
307:       // FIX: We would like to avoid the copy here with cone()
308:       if(this->_modifiedPoints->size() > 0) {
309:         this->computeHeight(height, sieve, sieve->cone(this->_modifiedPoints), maxHeight);
310:       }
311:     };
312:     void computeHeights() {
313:       const Obj<label_type>& label = this->createLabel(std::string("height"));

315:       this->_maxHeight = -1;
316:       this->computeHeight(label, this->_sieve, this->_sieve->leaves(), this->_maxHeight);
317:     };
318:     virtual int height() const {return this->_maxHeight;};
319:     virtual int height(const point_type& point) {
320:       return this->getValue(this->_labels["height"], point, -1);
321:     };
322:     virtual const Obj<label_sequence>& heightStratum(int height) {
323:       return this->getLabelStratum("height", height);
324:     };
325:     void setHeight(const Obj<label_type>& label) {
326:       this->_labels["height"] = label;
327:       const Obj<typename label_type::traits::capSequence> cap = label->cap();

329:       for(typename label_type::traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
330:         this->_maxHeight = std::max(this->_maxHeight, *c_iter);
331:       }
332:     };
333:     template<class InputPoints>
334:     void computeDepth(const Obj<label_type>& depth, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxDepth) {
335:       this->_modifiedPoints->clear();

337:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
338:         // Compute the max depth of the points in the cone of p, and add 1
339:         int d0 = this->getValue(depth, *p_iter, -1);
340:         int d1 = this->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;

342:         if(d1 != d0) {
343:           this->setValue(depth, *p_iter, d1);
344:           if (d1 > maxDepth) maxDepth = d1;
345:           this->_modifiedPoints->insert(*p_iter);
346:         }
347:       }
348:       // FIX: We would like to avoid the copy here with support()
349:       if(this->_modifiedPoints->size() > 0) {
350:         this->computeDepth(depth, sieve, sieve->support(this->_modifiedPoints), maxDepth);
351:       }
352:     };
353:     void computeDepths() {
354:       const Obj<label_type>& label = this->createLabel(std::string("depth"));

356:       this->_maxDepth = -1;
357:       this->computeDepth(label, this->_sieve, this->_sieve->roots(), this->_maxDepth);
358:     };
359:     virtual int depth() const {return this->_maxDepth;};
360:     virtual int depth(const point_type& point) {
361:       return this->getValue(this->_labels["depth"], point, -1);
362:     };
363:     virtual const Obj<label_sequence>& depthStratum(int depth) {
364:       return this->getLabelStratum("depth", depth);
365:     };
368:     virtual void stratify() {
369:       ALE_LOG_EVENT_BEGIN;
370:       this->computeHeights();
371:       this->computeDepths();
372:       ALE_LOG_EVENT_END;
373:     };
374:   public: // Size traversal
375:     template<typename Section_>
376:     int size(const Obj<Section_>& section, const point_type& p) {
377:       const typename Section_::chart_type& chart = section->getChart();
378:       int                                  size  = 0;

380:       if (this->height() < 2) {
381:         const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
382:         typename sieve_type::coneSequence::iterator   end  = cone->end();

384:         if (chart.count(p)) {
385:           size += section->getConstrainedFiberDimension(p);
386:         }
387:         for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
388:           if (chart.count(*c_iter)) {
389:             size += section->getConstrainedFiberDimension(*c_iter);
390:           }
391:         }
392:       } else {
393:         const Obj<coneArray>         closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), p);
394:         typename coneArray::iterator end     = closure->end();

396:         for(typename coneArray::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
397:           if (chart.count(*c_iter)) {
398:             size += section->getConstrainedFiberDimension(*c_iter);
399:           }
400:         }
401:       }
402:       return size;
403:     };
404:     template<typename Section_>
405:     int sizeWithBC(const Obj<Section_>& section, const point_type& p) {
406:       const typename Section_::chart_type& chart = section->getChart();
407:       int                                  size  = 0;

409:       if (this->height() < 2) {
410:         const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
411:         typename sieve_type::coneSequence::iterator   end  = cone->end();

413:         if (chart.count(p)) {
414:           size += section->getFiberDimension(p);
415:         }
416:         for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
417:           if (chart.count(*c_iter)) {
418:             size += section->getFiberDimension(*c_iter);
419:           }
420:         }
421:       } else {
422:         const Obj<coneArray>         closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), p);
423:         typename coneArray::iterator end     = closure->end();

425:         for(typename coneArray::iterator c_iter = closure->begin(); c_iter != end; ++c_iter) {
426:           if (chart.count(*c_iter)) {
427:             size += section->getFiberDimension(*c_iter);
428:           }
429:         }
430:       }
431:       return size;
432:     };
433:   protected:
434:     int *getIndexArray(const int size) {
435:       static int *array   = NULL;
436:       static int  maxSize = 0;

438:       if (size > maxSize) {
439:         maxSize = size;
440:         if (array) delete [] array;
441:         array = new int[maxSize];
442:       };
443:       return array;
444:     };
445:   public: // Index traversal
446:     void expandInterval(const index_type& interval, PetscInt indices[], PetscInt *indx) {
447:       const int end = interval.prefix + interval.index;

449:       for(int i = interval.index; i < end; ++i) {
450:         indices[(*indx)++] = i;
451:       }
452:     };
453:     void expandInterval(const index_type& interval, const int orientation, PetscInt indices[], PetscInt *indx) {
454:       if (orientation >= 0) {
455:         for(int i = 0; i < interval.prefix; ++i) {
456:           indices[(*indx)++] = interval.index + i;
457:         }
458:       } else {
459:         for(int i = interval.prefix-1; i >= 0; --i) {
460:           indices[(*indx)++] = interval.index + i;
461:         }
462:       }
463:       for(int i = 0; i < -interval.prefix; ++i) {
464:         indices[(*indx)++] = -1;
465:       }
466:     };
467:     void expandIntervals(Obj<oIndexArray> intervals, PetscInt *indices) {
468:       int k = 0;

470:       for(typename oIndexArray::iterator i_iter = intervals->begin(); i_iter != intervals->end(); ++i_iter) {
471:         this->expandInterval(i_iter->first, i_iter->second, indices, &k);
472:       }
473:     }
474:     template<typename Section_>
475:     const indices_type getIndicesRaw(const Obj<Section_>& section, const point_type& p) {
476:       int *indexArray = NULL;
477:       int  size       = 0;

479:       const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
480:       typename oConeArray::iterator begin   = closure->begin();
481:       typename oConeArray::iterator end     = closure->end();

483:       for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
484:         size    += section->getFiberDimension(p_iter->first);
485:       }
486:       indexArray = this->getIndexArray(size);
487:       int  k     = 0;
488:       for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
489:         section->getIndicesRaw(p_iter->first, section->getIndex(p_iter->first), indexArray, &k, p_iter->second);
490:       }
491:       return indices_type(indexArray, size);
492:     };
493:     template<typename Section_>
494:     const indices_type getIndices(const Obj<Section_>& section, const point_type& p, const int level = -1) {
495:       int *indexArray = NULL;
496:       int  size       = 0;

498:       if (level == 0) {
499:         size      += section->getFiberDimension(p);
500:         indexArray = this->getIndexArray(size);
501:         int  k     = 0;

503:         section->getIndices(p, indexArray, &k);
504:       } else if ((level == 1) || (this->height() == 1)) {
505:         const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
506:         typename sieve_type::coneSequence::iterator   end  = cone->end();

508:         size      += section->getFiberDimension(p);
509:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
510:           size    += section->getFiberDimension(*p_iter);
511:         }
512:         indexArray = this->getIndexArray(size);
513:         int  k     = 0;

515:         section->getIndices(p, indexArray, &k);
516:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
517:           section->getIndices(*p_iter, indexArray, &k);
518:         }
519:       } else if (level == -1) {
520:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
521:         typename oConeArray::iterator begin   = closure->begin();
522:         typename oConeArray::iterator end     = closure->end();

524:         for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
525:           size    += section->getFiberDimension(p_iter->first);
526:         }
527:         indexArray = this->getIndexArray(size);
528:         int  k     = 0;
529:         for(typename oConeArray::iterator p_iter = begin; p_iter != end; ++p_iter) {
530:           section->getIndices(p_iter->first, indexArray, &k, p_iter->second);
531:         }
532:       } else {
533:         throw ALE::Exception("Bundle has not yet implemented getIndices() for an arbitrary level");
534:       }
535:       if (this->debug()) {
536:         for(int i = 0; i < size; ++i) {
537:           printf("[%d]index %d: %d\n", this->commRank(), i, indexArray[i]);
538:         }
539:       }
540:       return indices_type(indexArray, size);
541:     };
542:     template<typename Section_, typename Numbering>
543:     const indices_type getIndices(const Obj<Section_>& section, const point_type& p, const Obj<Numbering>& numbering, const int level = -1) {
544:       int *indexArray = NULL;
545:       int  size       = 0;

547:       if (level == 0) {
548:         size      += section->getFiberDimension(p);
549:         indexArray = this->getIndexArray(size);
550:         int  k     = 0;

552:         section->getIndices(p, numbering, indexArray, &k);
553:       } else if ((level == 1) || (this->height() == 1)) {
554:         const Obj<typename sieve_type::coneSequence>& cone = this->_sieve->cone(p);
555:         typename sieve_type::coneSequence::iterator   end  = cone->end();

557:         size      += section->getFiberDimension(p);
558:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
559:           size    += section->getFiberDimension(*p_iter);
560:         }
561:         indexArray = this->getIndexArray(size);
562:         int  k     = 0;

564:         section->getIndices(p, numbering, indexArray, &k);
565:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
566:           section->getIndices(*p_iter, numbering, indexArray, &k);
567:         }
568:       } else if (level == -1) {
569:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
570:         typename oConeArray::iterator end     = closure->end();

572:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
573:           size    += section->getFiberDimension(p_iter->first);
574:         }
575:         indexArray = this->getIndexArray(size);
576:         int  k     = 0;
577:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
578:           section->getIndices(p_iter->first, numbering, indexArray, &k, p_iter->second);
579:         }
580:       } else {
581:         throw ALE::Exception("Bundle has not yet implemented getIndices() for an arbitrary level");
582:       }
583:       return indices_type(indexArray, size);
584:     };
585:   public: // Retrieval traversal
586:     // Return the values for the closure of this point
587:     //   use a smart pointer?
588:     template<typename Section_>
589:     const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p) {
590:       const int size = this->sizeWithBC(section, p);
591:       return this->restrictClosure(section, p, section->getRawArray(size), size);
592:     };
593:     template<typename Section_>
594:     const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const point_type& p, typename Section_::value_type  *values, const int valuesSize) {
595:       const int size = this->sizeWithBC(section, p);
596:       int       j    = -1;
597:       if (valuesSize < size) throw ALE::Exception("Input array too small");

599:       // We could actually ask for the height of the individual point
600:       if (this->height() < 2) {
601:         const int& dim = section->getFiberDimension(p);
602:         const typename Section_::value_type *array = section->restrictPoint(p);

604:         for(int i = 0; i < dim; ++i) {
605:           values[++j] = array[i];
606:         }
607:         const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);
608:         typename sieve_type::coneSequence::iterator   end  = cone->end();

610:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
611:           const int& dim = section->getFiberDimension(*p_iter);

613:           array = section->restrictPoint(*p_iter);
614:           for(int i = 0; i < dim; ++i) {
615:             values[++j] = array[i];
616:           }
617:         }
618:       } else {
619:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
620:         typename oConeArray::iterator end     = closure->end();

622:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
623:           const typename Section_::value_type *array = section->restrictPoint(p_iter->first);
624:           const int& dim = section->getFiberDimension(p_iter->first);

626:           if (p_iter->second >= 0) {
627:             for(int i = 0; i < dim; ++i) {
628:               values[++j] = array[i];
629:             }
630:           } else {
631:             for(int i = dim-1; i >= 0; --i) {
632:               values[++j] = array[i];
633:             }
634:           }
635:         }
636:       }
637:       if (j != size-1) {
638:         ostringstream txt;

640:         txt << "Invalid restrict to point " << p << std::endl;
641:         txt << "  j " << j << " should be " << (size-1) << std::endl;
642:         std::cout << txt.str();
643:         throw ALE::Exception(txt.str().c_str());
644:       }
645:       return values;
646:     };
647:     template<typename Section_>
648:     const typename Section_::value_type *restrictNew(const Obj<Section_>& section, const point_type& p) {
649:       const int size = this->sizeWithBC(section, p);
650:       return this->restrictNew(section, p, section->getRawArray(size), size);
651:     };
652:     template<typename Section_>
653:     const typename Section_::value_type *restrictNew(const Obj<Section_>& section, const point_type& p, typename Section_::value_type  *values, const int valuesSize) {
654:       const int                     size    = this->sizeWithBC(section, p);
655:       const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
656:       typename oConeArray::iterator end     = closure->end();
657:       int                           j       = -1;
658:       if (valuesSize < size) throw ALE::Exception("Input array too small");

660:       for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
661:         const typename Section_::value_type *array = section->restrictPoint(p_iter->first);

663:         if (p_iter->second >= 0) {
664:           const int& dim = section->getFiberDimension(p_iter->first);

666:           for(int i = 0; i < dim; ++i) {
667:             values[++j] = array[i];
668:           }
669:         } else {
670:           int offset = 0;

672:           for(int space = 0; space < section->getNumSpaces(); ++space) {
673:             const int& dim = section->getFiberDimension(p_iter->first, space);

675:             for(int i = dim-1; i >= 0; --i) {
676:               values[++j] = array[i+offset];
677:             }
678:             offset += dim;
679:           }
680:         }
681:       }
682:       if (j != size-1) {
683:         ostringstream txt;

685:         txt << "Invalid restrict to point " << p << std::endl;
686:         txt << "  j " << j << " should be " << (size-1) << std::endl;
687:         std::cout << txt.str();
688:         throw ALE::Exception(txt.str().c_str());
689:       }
690:       return values;
691:     };
692:     template<typename Section_>
693:     void update(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
694:       int j = 0;

696:       if (this->height() < 2) {
697:         section->updatePoint(p, &v[j]);
698:         j += section->getFiberDimension(p);
699:         const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

701:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
702:           section->updatePoint(*p_iter, &v[j]);
703:           j += section->getFiberDimension(*p_iter);
704:         }
705:       } else {
706:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
707:         typename oConeArray::iterator end     = closure->end();

709:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
710:           section->updatePoint(p_iter->first, &v[j], p_iter->second);
711:           j += section->getFiberDimension(p_iter->first);
712:         }
713:       }
714:     };
715:     template<typename Section_>
716:     void updateAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
717:       int j = 0;

719:       if (this->height() < 2) {
720:         section->updateAddPoint(p, &v[j]);
721:         j += section->getFiberDimension(p);
722:         const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

724:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
725:           section->updateAddPoint(*p_iter, &v[j]);
726:           j += section->getFiberDimension(*p_iter);
727:         }
728:       } else {
729:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
730:         typename oConeArray::iterator end     = closure->end();

732:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
733:           section->updateAddPoint(p_iter->first, &v[j], p_iter->second);
734:           j += section->getFiberDimension(p_iter->first);
735:         }
736:       }
737:     };
738:     template<typename Section_>
739:     void updateBC(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
740:       int j = 0;

742:       if (this->height() < 2) {
743:         section->updatePointBC(p, &v[j]);
744:         j += section->getFiberDimension(p);
745:         const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

747:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
748:           section->updatePointBC(*p_iter, &v[j]);
749:           j += section->getFiberDimension(*p_iter);
750:         }
751:       } else {
752:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
753:         typename oConeArray::iterator end     = closure->end();

755:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
756:           section->updatePointBC(p_iter->first, &v[j], p_iter->second);
757:           j += section->getFiberDimension(p_iter->first);
758:         }
759:       }
760:     };
761:     template<typename Section_>
762:     void updateAll(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
763:       int j = 0;

765:       if (this->height() < 2) {
766:         section->updatePointAll(p, &v[j]);
767:         j += section->getFiberDimension(p);
768:         const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

770:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
771:           section->updatePointAll(*p_iter, &v[j]);
772:           j += section->getFiberDimension(*p_iter);
773:         }
774:       } else {
775:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
776:         typename oConeArray::iterator end     = closure->end();

778:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
779:           section->updatePointAll(p_iter->first, &v[j], p_iter->second);
780:           j += section->getFiberDimension(p_iter->first);
781:         }
782:       }
783:     };
784:     template<typename Section_>
785:     void updateAllAdd(const Obj<Section_>& section, const point_type& p, const typename Section_::value_type v[]) {
786:       int j = 0;

788:       if (this->height() < 2) {
789:         section->updatePointAllAdd(p, &v[j]);
790:         j += section->getFiberDimension(p);
791:         const Obj<typename sieve_type::coneSequence>& cone = this->getSieve()->cone(p);

793:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
794:           section->updatePointAllAdd(*p_iter, &v[j]);
795:           j += section->getFiberDimension(*p_iter);
796:         }
797:       } else {
798:         const Obj<oConeArray>         closure = sieve_alg_type::orientedClosure(this, this->getArrowSection("orientation"), p);
799:         typename oConeArray::iterator end     = closure->end();

801:         for(typename oConeArray::iterator p_iter = closure->begin(); p_iter != end; ++p_iter) {
802:           section->updatePointAllAdd(p_iter->first, &v[j], p_iter->second);
803:           j += section->getFiberDimension(p_iter->first);
804:         }
805:       }
806:     };
807:   public: // Optimization
808:     // Calculate a custom atlas for the given traversal
809:     //   This returns the tag value assigned to the traversal
810:     template<typename Section_, typename Sequence_>
811:     int calculateCustomAtlas(const Obj<Section_>& section, const Obj<Sequence_>& points) {
812:       const typename Sequence_::iterator begin    = points->begin();
813:       const typename Sequence_::iterator end      = points->end();
814:       const int                          num      = points->size();
815:       int                               *rOffsets = new int[num+1];
816:       int                               *rIndices;
817:       int                               *uOffsets = new int[num+1];
818:       int                               *uIndices;
819:       int                                p;

821:       p = 0;
822:       rOffsets[p] = 0;
823:       uOffsets[p] = 0;
824:       for(typename Sequence_::iterator p_iter = begin; p_iter != end; ++p_iter, ++p) {
825:         rOffsets[p+1] = rOffsets[p] + this->sizeWithBC(section, *p_iter);
826:         uOffsets[p+1] = rOffsets[p+1];
827:         //uOffsets[p+1] = uOffsets[p] + this->size(section, *p_iter);
828:       }
829:       rIndices = new int[rOffsets[p]];
830:       uIndices = new int[uOffsets[p]];
831:       p = 0;
832:       for(typename Sequence_::iterator p_iter = begin; p_iter != end; ++p_iter, ++p) {
833:         const indices_type rIdx = this->getIndicesRaw(section, *p_iter);
834:         for(int i = 0, k = rOffsets[p]; k < rOffsets[p+1]; ++i, ++k) rIndices[k] = rIdx.first[i];

836:         const indices_type uIdx = this->getIndices(section, *p_iter);
837:         for(int i = 0, k = uOffsets[p]; k < uOffsets[p+1]; ++i, ++k) uIndices[k] = uIdx.first[i];
838:       }
839:       return section->setCustomAtlas(rOffsets, rIndices, uOffsets, uIndices);
840:     };
841:     template<typename Section_>
842:     const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p) {
843:       const int *offsets, *indices;

845:       section->getCustomRestrictAtlas(tag, &offsets, &indices);
846:       const int size = offsets[p+1] - offsets[p];
847:       return this->restrictClosure(section, tag, p, section->getRawArray(size), offsets, indices);
848:     };
849:     template<typename Section_>
850:     const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p, typename Section_::value_type  *values, const int valuesSize) {
851:       const int *offsets, *indices;

853:       section->getCustomRestrictAtlas(tag, &offsets, &indices);
854:       const int size = offsets[p+1] - offsets[p];
855:       if (valuesSize < size) {throw ALE::Exception("Input array too small");}
856:       return this->restrictClosure(section, tag, p, values, offsets, indices);
857:     };
858:     template<typename Section_>
859:     const typename Section_::value_type *restrictClosure(const Obj<Section_>& section, const int tag, const int p, typename Section_::value_type  *values, const int offsets[], const int indices[]) {
860:       const typename Section_::value_type *array = section->restrictSpace();

862:       const int size = offsets[p+1] - offsets[p];
863:       for(int j = 0, k = offsets[p]; j < size; ++j, ++k) {
864:         values[j] = array[indices[k]];
865:       }
866:       return values;
867:     };
868:     template<typename Section_>
869:     void updateAdd(const Obj<Section_>& section, const int tag, const int p, const typename Section_::value_type values[]) {
870:       typename Section_::value_type *array = (typename Section_::value_type *) section->restrictSpace();
871:       const int *offsets, *indices;

873:       section->getCustomUpdateAtlas(tag, &offsets, &indices);
874:       const int size = offsets[p+1] - offsets[p];
875:       for(int j = 0, k = offsets[p]; j < size; ++j, ++k) {
876:         if (indices[k] < 0) continue;
877:         array[indices[k]] += values[j];
878:       }
879:     };
880:   public: // Allocation
881:     template<typename Section_>
882:     void allocate(const Obj<Section_>& section, const Obj<send_overlap_type>& sendOverlap = NULL) {
883:       bool doGhosts = !sendOverlap.isNull();

885:       this->_factory->orderPatch(section, this->getSieve(), sendOverlap);
886:       if (doGhosts) {
887:         if (this->_debug > 1) {std::cout << "Ordering patch for ghosts" << std::endl;}
888:         const typename Section_::chart_type& points = section->getChart();
889:         typename Section_::index_type::index_type offset = 0;

891:         for(typename Section_::chart_type::const_iterator point = points.begin(); point != points.end(); ++point) {
892:           const typename Section_::index_type& idx = section->getIndex(*point);

894:           offset = std::max(offset, idx.index + std::abs(idx.prefix));
895:         }
896:         this->_factory->orderPatch(section, this->getSieve(), NULL, offset);
897:         if (offset != section->sizeWithBC()) throw ALE::Exception("Inconsistent array sizes in section");
898:       }
899:       section->allocateStorage();
900:     };
901:     template<typename Section_>
902:     void reallocate(const Obj<Section_>& section) {
903:       if (section->getNewAtlas().isNull()) return;
904:       // Since copy() preserves offsets, we must reinitialize them before ordering
905:       const Obj<typename Section_::atlas_type>         atlas    = section->getAtlas();
906:       const Obj<typename Section_::atlas_type>&        newAtlas = section->getNewAtlas();
907:       const typename Section_::atlas_type::chart_type& chart    = newAtlas->getChart();
908:       const typename Section_::atlas_type::chart_type& oldChart = atlas->getChart();
909:       int                                              newSize  = 0;
910:       typename Section_::index_type                    defaultIdx(0, -1);

912:       for(typename Section_::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
913:         defaultIdx.prefix = newAtlas->restrictPoint(*c_iter)[0].prefix;
914:         newAtlas->updatePoint(*c_iter, &defaultIdx);
915:         newSize += defaultIdx.prefix;
916:       }
917:       section->setAtlas(newAtlas);
918:       this->_factory->orderPatch(section, this->getSieve());
919:       // Copy over existing values
920:       typedef typename alloc_type::template rebind<typename Section_::value_type>::other value_alloc_type;
921:       value_alloc_type value_allocator;
922:       typename Section_::value_type                   *newArray = value_allocator.allocate(newSize);
923:       for(int i = 0; i < newSize; ++i) {value_allocator.construct(newArray+i, typename Section_::value_type());}
924:       ///typename Section_::value_type                   *newArray = new typename Section_::value_type[newSize];
925:       const typename Section_::value_type             *array    = section->restrictSpace();

927:       for(typename Section_::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
928:         const int& dim       = section->getFiberDimension(*c_iter);
929:         const int& offset    = atlas->restrictPoint(*c_iter)->index;
930:         const int& newOffset = newAtlas->restrictPoint(*c_iter)->index;

932:         for(int i = 0; i < dim; ++i) {
933:           newArray[newOffset+i] = array[offset+i];
934:         }
935:       }
936:       section->replaceStorage(newArray);
937:     };
938:   public: // Overlap
939:     template<typename Sequence>
940:     void constructOverlap(const Obj<Sequence>& points, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
941:       point_type *sendBuf = new point_type[points->size()];
942:       int         size    = 0;
943:       for(typename Sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
944:         sendBuf[size++] = *l_iter;
945:       }
946:       int *sizes   = new int[this->commSize()];   // The number of points coming from each process
947:       int *offsets = new int[this->commSize()+1]; // Prefix sums for sizes
948:       int *oldOffs = new int[this->commSize()+1]; // Temporary storage
949:       point_type *remotePoints = NULL;            // The points from each process
950:       int        *remoteRanks  = NULL;            // The rank and number of overlap points of each process that overlaps another

952:       // Change to Allgather() for the correct binning algorithm
953:       MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, this->comm());
954:       if (this->commRank() == 0) {
955:         offsets[0] = 0;
956:         for(int p = 1; p <= this->commSize(); p++) {
957:           offsets[p] = offsets[p-1] + sizes[p-1];
958:         }
959:         remotePoints = new point_type[offsets[this->commSize()]];
960:       }
961:       MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, this->comm());
962:       delete [] sendBuf;
963:       std::map<int, std::map<int, std::set<point_type> > > overlapInfo; // Maps (p,q) to their set of overlap points

965:       if (this->commRank() == 0) {
966:         for(int p = 0; p < this->commSize(); p++) {
967:           std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]]);
968:         }
969:         for(int p = 0; p <= this->commSize(); p++) {
970:           oldOffs[p] = offsets[p];
971:         }
972:         for(int p = 0; p < this->commSize(); p++) {
973:           for(int q = p+1; q < this->commSize(); q++) {
974:             std::set_intersection(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
975:                                   &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
976:                                   std::insert_iterator<std::set<point_type> >(overlapInfo[p][q], overlapInfo[p][q].begin()));
977:             overlapInfo[q][p] = overlapInfo[p][q];
978:           }
979:           sizes[p]     = overlapInfo[p].size()*2;
980:           offsets[p+1] = offsets[p] + sizes[p];
981:         }
982:         remoteRanks = new int[offsets[this->commSize()]];
983:         int       k = 0;
984:         for(int p = 0; p < this->commSize(); p++) {
985:           for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
986:             remoteRanks[k*2]   = r_iter->first;
987:             remoteRanks[k*2+1] = r_iter->second.size();
988:             k++;
989:           }
990:         }
991:       }
992:       int numOverlaps;                          // The number of processes overlapping this process
993:       MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, this->comm());
994:       int *overlapRanks = new int[numOverlaps]; // The rank and overlap size for each overlapping process
995:       MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, this->comm());
996:       point_type *sendPoints = NULL;            // The points to send to each process
997:       if (this->commRank() == 0) {
998:         for(int p = 0, k = 0; p < this->commSize(); p++) {
999:           sizes[p] = 0;
1000:           for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
1001:             sizes[p] += remoteRanks[k*2+1];
1002:             k++;
1003:           }
1004:           offsets[p+1] = offsets[p] + sizes[p];
1005:         }
1006:         sendPoints = new point_type[offsets[this->commSize()]];
1007:         for(int p = 0, k = 0; p < this->commSize(); p++) {
1008:           for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
1009:             int rank = r_iter->first;
1010:             for(typename std::set<point_type>::iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
1011:               sendPoints[k++] = *p_iter;
1012:             }
1013:           }
1014:         }
1015:       }
1016:       int numOverlapPoints = 0;
1017:       for(int r = 0; r < numOverlaps/2; r++) {
1018:         numOverlapPoints += overlapRanks[r*2+1];
1019:       }
1020:       point_type *overlapPoints = new point_type[numOverlapPoints];
1021:       MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints, MPI_INT, 0, this->comm());

1023:       for(int r = 0, k = 0; r < numOverlaps/2; r++) {
1024:         int rank = overlapRanks[r*2];

1026:         for(int p = 0; p < overlapRanks[r*2+1]; p++) {
1027:           point_type point = overlapPoints[k++];

1029:           sendOverlap->addArrow(point, rank, point);
1030:           recvOverlap->addArrow(rank, point, point);
1031:         }
1032:       }

1034:       delete [] overlapPoints;
1035:       delete [] overlapRanks;
1036:       delete [] sizes;
1037:       delete [] offsets;
1038:       delete [] oldOffs;
1039:       if (this->commRank() == 0) {
1040:         delete [] remoteRanks;
1041:         delete [] remotePoints;
1042:         delete [] sendPoints;
1043:       }
1044:     };
1045:     void constructOverlap() {
1046:       if (this->_calculatedOverlap) return;
1047:       this->constructOverlap(this->getSieve()->base(), this->getSendOverlap(), this->getRecvOverlap());
1048:       this->constructOverlap(this->getSieve()->cap(),  this->getSendOverlap(), this->getRecvOverlap());
1049:       if (this->debug()) {
1050:         this->_sendOverlap->view("Send overlap");
1051:         this->_recvOverlap->view("Receive overlap");
1052:       }
1053:       this->_calculatedOverlap = true;
1054:     };
1055:   };
1056:   class BoundaryCondition : public ALE::ParallelObject {
1057:   public:
1058:     typedef double (*function_type)(const double []);
1059:     typedef double (*integrator_type)(const double [], const double [], const int, function_type);
1060:   protected:
1061:     std::string     _labelName;
1062:     int             _marker;
1063:     function_type   _func;
1064:     integrator_type _integrator;
1065:   public:
1066:     BoundaryCondition(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {};
1067:     ~BoundaryCondition() {};
1068:   public:
1069:     const std::string& getLabelName() const {return this->_labelName;};
1070:     void setLabelName(const std::string& name) {this->_labelName = name;};
1071:     int getMarker() const {return this->_marker;};
1072:     void setMarker(const int marker) {this->_marker = marker;};
1073:     function_type getFunction() const {return this->_func;};
1074:     void setFunction(function_type func) {this->_func = func;};
1075:     integrator_type getDualIntegrator() const {return this->_integrator;};
1076:     void setDualIntegrator(integrator_type integrator) {this->_integrator = integrator;};
1077:   public:
1078:     double evaluate(const double coords[]) const {return this->_func(coords);};
1079:     double integrateDual(const double v0[], const double J[], const int dualIndex) const {return this->_integrator(v0, J, dualIndex, this->_func);};
1080:   };
1081:   class Discretization : public ALE::ParallelObject {
1082:     typedef std::map<std::string, Obj<BoundaryCondition> > boundaryConditions_type;
1083:   protected:
1084:     boundaryConditions_type _boundaryConditions;
1085:     Obj<BoundaryCondition> _exactSolution;
1086:     std::map<int,int> _dim2dof;
1087:     std::map<int,int> _dim2class;
1088:     int           _quadSize;
1089:     const double *_points;
1090:     const double *_weights;
1091:     int           _basisSize;
1092:     const double *_basis;
1093:     const double *_basisDer;
1094:     const int    *_indices;
1095:     std::map<int, const int *> _exclusionIndices;
1096:   public:
1097:     Discretization(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _quadSize(0), _points(NULL), _weights(NULL), _basisSize(0), _basis(NULL), _basisDer(NULL), _indices(NULL) {};
1098:     virtual ~Discretization() {
1099:       if (this->_indices) {delete [] this->_indices;}
1100:       for(std::map<int, const int *>::iterator i_iter = _exclusionIndices.begin(); i_iter != _exclusionIndices.end(); ++i_iter) {
1101:         delete [] i_iter->second;
1102:       }
1103:     };
1104:   public:
1105:     const bool hasBoundaryCondition() {return (this->_boundaryConditions.find("default") != this->_boundaryConditions.end());};
1106:     const Obj<BoundaryCondition>& getBoundaryCondition() {return this->getBoundaryCondition("default");};
1107:     void setBoundaryCondition(const Obj<BoundaryCondition>& boundaryCondition) {this->setBoundaryCondition("default", boundaryCondition);};
1108:     const Obj<BoundaryCondition>& getBoundaryCondition(const std::string& name) {return this->_boundaryConditions[name];};
1109:     void setBoundaryCondition(const std::string& name, const Obj<BoundaryCondition>& boundaryCondition) {this->_boundaryConditions[name] = boundaryCondition;};
1110:     Obj<std::set<std::string> > getBoundaryConditions() const {
1111:       Obj<std::set<std::string> > names = std::set<std::string>();

1113:       for(boundaryConditions_type::const_iterator d_iter = this->_boundaryConditions.begin(); d_iter != this->_boundaryConditions.end(); ++d_iter) {
1114:         names->insert(d_iter->first);
1115:       }
1116:       return names;
1117:     };
1118:     const Obj<BoundaryCondition>& getExactSolution() {return this->_exactSolution;};
1119:     void setExactSolution(const Obj<BoundaryCondition>& exactSolution) {this->_exactSolution = exactSolution;};
1120:     const int     getQuadratureSize() {return this->_quadSize;};
1121:     void          setQuadratureSize(const int size) {this->_quadSize = size;};
1122:     const double *getQuadraturePoints() {return this->_points;};
1123:     void          setQuadraturePoints(const double *points) {this->_points = points;};
1124:     const double *getQuadratureWeights() {return this->_weights;};
1125:     void          setQuadratureWeights(const double *weights) {this->_weights = weights;};
1126:     const int     getBasisSize() {return this->_basisSize;};
1127:     void          setBasisSize(const int size) {this->_basisSize = size;};
1128:     const double *getBasis() {return this->_basis;};
1129:     void          setBasis(const double *basis) {this->_basis = basis;};
1130:     const double *getBasisDerivatives() {return this->_basisDer;};
1131:     void          setBasisDerivatives(const double *basisDer) {this->_basisDer = basisDer;};
1132:     int  getNumDof(const int dim) {return this->_dim2dof[dim];};
1133:     void setNumDof(const int dim, const int numDof) {this->_dim2dof[dim] = numDof;};
1134:     int  getDofClass(const int dim) {return this->_dim2class[dim];};
1135:     void setDofClass(const int dim, const int dofClass) {this->_dim2class[dim] = dofClass;};
1136:   public:
1137:     const int *getIndices() {return this->_indices;};
1138:     const int *getIndices(const int marker) {
1139:       if (!marker) return this->getIndices();
1140:       return this->_exclusionIndices[marker];
1141:     };
1142:     void       setIndices(const int *indices) {this->_indices = indices;};
1143:     void       setIndices(const int *indices, const int marker) {
1144:       if (!marker) this->_indices = indices;
1145:       this->_exclusionIndices[marker] = indices;
1146:     };
1147:     template<typename Bundle>
1148:     int sizeV(Bundle& mesh) {
1149:       typedef typename ISieveVisitor::PointRetriever<typename Bundle::sieve_type> Visitor;
1150:       Visitor pV((int) pow((double) mesh.getSieve()->getMaxConeSize(), mesh.depth())+1, true);
1151:       ISieveTraversal<typename Bundle::sieve_type>::orientedClosure(*mesh.getSieve(), *mesh.heightStratum(0)->begin(), pV);
1152:       const typename Visitor::point_type *oPoints = pV.getPoints();
1153:       const int                           oSize   = pV.getSize();
1154:       int                                 size    = 0;

1156:       for(int cl = 0; cl < oSize; ++cl) {
1157:         size += this->_dim2dof[mesh.depth(oPoints[cl])];
1158:       }
1159:       return size;
1160:     };
1161:     template<typename Bundle>
1162:     int size(const Obj<Bundle>& mesh) {
1163:       const Obj<typename Bundle::label_sequence>& cells   = mesh->heightStratum(0);
1164:       const Obj<typename Bundle::coneArray>       closure = ALE::SieveAlg<Bundle>::closure(mesh, *cells->begin());
1165:       const typename Bundle::coneArray::iterator  end     = closure->end();
1166:       int                                         size    = 0;

1168:       for(typename Bundle::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1169:         size += this->_dim2dof[mesh->depth(*cl_iter)];
1170:       }
1171:       return size;
1172:     };
1173:   };
1174: }

1176: namespace ALE {
1177:   template<typename Sieve_,
1178:            typename RealSection_  = Section<typename Sieve_::point_type, double>,
1179:            typename IntSection_   = Section<typename Sieve_::point_type, int>,
1180:            typename Label_        = LabelSifter<int, typename Sieve_::point_type>,
1181:            typename ArrowSection_ = UniformSection<MinimalArrow<typename Sieve_::point_type, typename Sieve_::point_type>, int> >
1182:   class IBundle : public ALE::ParallelObject {
1183:   public:
1184:     typedef Sieve_                                                    sieve_type;
1185:     typedef RealSection_                                              real_section_type;
1186:     typedef IntSection_                                               int_section_type;
1187:     typedef ArrowSection_                                             arrow_section_type;
1188:     typedef IBundle<Sieve_,RealSection_,IntSection_,Label_,ArrowSection_> this_type;
1189:     typedef typename sieve_type::point_type                           point_type;
1190:     typedef malloc_allocator<point_type>                              alloc_type;
1191:     typedef Label_                                                    label_type;
1192:     typedef typename std::map<const std::string, Obj<label_type> >    labels_type;
1193:     typedef typename label_type::supportSequence                      label_sequence;
1194:     typedef std::map<std::string, Obj<arrow_section_type> >           arrow_sections_type;
1195:     typedef std::map<std::string, Obj<real_section_type> >            real_sections_type;
1196:     typedef std::map<std::string, Obj<int_section_type> >             int_sections_type;
1197:     typedef ALE::Point                                                index_type;
1198:     typedef std::pair<index_type, int>                                oIndex_type;
1199:     typedef std::vector<oIndex_type>                                  oIndexArray;
1200:     typedef std::pair<int *, int>                                     indices_type;
1201:     typedef NumberingFactory<this_type>                               MeshNumberingFactory;
1202:     typedef typename ALE::Partitioner<>::part_type                    rank_type;
1203:     typedef typename ALE::Sifter<point_type,rank_type,point_type>     send_overlap_type;
1204:     typedef typename ALE::Sifter<rank_type,point_type,point_type>     recv_overlap_type;
1205:     typedef typename MeshNumberingFactory::numbering_type             numbering_type;
1206:     typedef typename MeshNumberingFactory::order_type                 order_type;
1207:     typedef std::map<point_type, point_type>                          renumbering_type;
1208:     // These should go away
1209:     typedef typename ALE::SieveAlg<this_type>                         sieve_alg_type;
1210:     typedef typename sieve_alg_type::coneArray                        coneArray;
1211:     typedef typename sieve_alg_type::orientedConeArray                oConeArray;
1212:     typedef typename sieve_alg_type::supportArray                     supportArray;
1213:   public:
1214:     class LabelVisitor {
1215:     protected:
1216:       label_type& label;
1217:       int         defaultValue;
1218:       int         value;
1219:     public:
1220:       LabelVisitor(label_type& l, const int defValue) : label(l), defaultValue(defValue), value(defValue) {};
1221:       int getLabelValue(const point_type& point) const {
1222:         const Obj<typename label_type::coneSequence>& cone = this->label.cone(point);

1224:         if (cone->size() == 0) return this->defaultValue;
1225:         return *cone->begin();
1226:       };
1227:       void setLabelValue(const point_type& point, const int value) {
1228:         this->label.setCone(value, point);
1229:       };
1230:       int getValue() const {return this->value;};
1231:     };
1232:     class MaxConeVisitor : public LabelVisitor {
1233:     public:
1234:       MaxConeVisitor(label_type& l, const int defValue) : LabelVisitor(l, defValue) {};
1235:       void visitPoint(const typename sieve_type::point_type& point) {};
1236:       void visitArrow(const typename sieve_type::arrow_type& arrow) {
1237:         this->value = std::max(this->value, this->getLabelValue(arrow.source));
1238:       };
1239:     };
1240:     class MaxSupportVisitor : public LabelVisitor {
1241:     public:
1242:       MaxSupportVisitor(label_type& l, const int defValue) : LabelVisitor(l, defValue) {};
1243:       void visitPoint(const typename sieve_type::point_type& point) {};
1244:       void visitArrow(const typename sieve_type::arrow_type& arrow) {
1245:         this->value = std::max(this->value, this->getLabelValue(arrow.target));
1246:       };
1247:     };
1248:     class HeightVisitor {
1249:     protected:
1250:       const sieve_type& sieve;
1251:       label_type&       height;
1252:       int               maxHeight;
1253:       std::set<typename sieve_type::point_type> modifiedPoints;
1254:     public:
1255:       HeightVisitor(const sieve_type& s, label_type& h) : sieve(s), height(h), maxHeight(-1) {};
1256:       void visitPoint(const typename sieve_type::point_type& point) {
1257:         MaxSupportVisitor v(height, -1);

1259:         // Compute the max height of the points in the support of p, and add 1
1260:         this->sieve.support(point, v);
1261:         const int h0 = v.getLabelValue(point);
1262:         const int h1 = v.getValue() + 1;

1264:         if(h1 != h0) {
1265:           v.setLabelValue(point, h1);
1266:           if (h1 > this->maxHeight) this->maxHeight = h1;
1267:           this->modifiedPoints.insert(point);
1268:         }
1269:       };
1270:       void visitArrow(const typename sieve_type::arrow_type& arrow) {
1271:         this->visitPoint(arrow.source);
1272:       };
1273:       int getMaxHeight() const {return this->maxHeight;};
1274:       bool isModified() const {return this->modifiedPoints.size() > 0;};
1275:       const std::set<typename sieve_type::point_type>& getModifiedPoints() const {return this->modifiedPoints;};
1276:       void clear() {this->modifiedPoints.clear();};
1277:     };
1278:     class DepthVisitor {
1279:     public:
1280:       typedef typename sieve_type::point_type point_type;
1281:     protected:
1282:       const sieve_type& sieve;
1283:       label_type&       depth;
1284:       int               maxDepth;
1285:       const point_type  limitPoint;
1286:       std::set<point_type> modifiedPoints;
1287:     public:
1288:       DepthVisitor(const sieve_type& s, label_type& d) : sieve(s), depth(d), maxDepth(-1), limitPoint(sieve.getChart().max()+1) {};
1289:       DepthVisitor(const sieve_type& s, const point_type& limit, label_type& d) : sieve(s), depth(d), maxDepth(-1), limitPoint(limit) {};
1290:       void visitPoint(const point_type& point) {
1291:         if (point >= this->limitPoint) return;
1292:         MaxConeVisitor v(depth, -1);

1294:         // Compute the max height of the points in the support of p, and add 1
1295:         this->sieve.cone(point, v);
1296:         const int d0 = v.getLabelValue(point);
1297:         const int d1 = v.getValue() + 1;

1299:         if(d1 != d0) {
1300:           v.setLabelValue(point, d1);
1301:           if (d1 > this->maxDepth) this->maxDepth = d1;
1302:           this->modifiedPoints.insert(point);
1303:         }
1304:       };
1305:       void visitArrow(const typename sieve_type::arrow_type& arrow) {
1306:         this->visitPoint(arrow.target);
1307:       };
1308:       int getMaxDepth() const {return this->maxDepth;};
1309:       bool isModified() const {return this->modifiedPoints.size() > 0;};
1310:       const std::set<typename sieve_type::point_type>& getModifiedPoints() const {return this->modifiedPoints;};
1311:       void clear() {this->modifiedPoints.clear();};
1312:     };
1313:   protected:
1314:     Obj<sieve_type>       _sieve;
1315:     labels_type           _labels;
1316:     int                   _maxHeight;
1317:     int                   _maxDepth;
1318:     arrow_sections_type   _arrowSections;
1319:     real_sections_type    _realSections;
1320:     int_sections_type     _intSections;
1321:     Obj<oIndexArray>      _indexArray;
1322:     Obj<MeshNumberingFactory> _factory;
1323:     bool                   _calculatedOverlap;
1324:     Obj<send_overlap_type> _sendOverlap;
1325:     Obj<recv_overlap_type> _recvOverlap;
1326:     Obj<send_overlap_type> _distSendOverlap;
1327:     Obj<recv_overlap_type> _distRecvOverlap;
1328:     renumbering_type       _renumbering;
1329:     // Work space
1330:     Obj<std::set<point_type> > _modifiedPoints;
1331:   public:
1332:     IBundle(MPI_Comm comm, int debug = 0) : ALE::ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1) {
1333:       this->_indexArray        = new oIndexArray();
1334:       this->_modifiedPoints    = new std::set<point_type>();
1335:       this->_factory           = MeshNumberingFactory::singleton(this->comm(), this->debug());
1336:       this->_calculatedOverlap = false;
1337:       this->_sendOverlap       = new send_overlap_type(comm, debug);
1338:       this->_recvOverlap       = new recv_overlap_type(comm, debug);
1339:     };
1340:     IBundle(const Obj<sieve_type>& sieve) : ALE::ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _maxHeight(-1), _maxDepth(-1) {
1341:       this->_indexArray        = new oIndexArray();
1342:       this->_modifiedPoints    = new std::set<point_type>();
1343:       this->_factory           = MeshNumberingFactory::singleton(this->comm(), this->debug());
1344:       this->_calculatedOverlap = false;
1345:       this->_sendOverlap       = new send_overlap_type(comm, debug);
1346:       this->_recvOverlap       = new recv_overlap_type(comm, debug);
1347:     };
1348:     virtual ~IBundle() {};
1349:   public: // Verifiers
1350:     bool hasLabel(const std::string& name) {
1351:       if (this->_labels.find(name) != this->_labels.end()) {
1352:         return true;
1353:       }
1354:       return false;
1355:     };
1356:     void checkLabel(const std::string& name) {
1357:       if (!this->hasLabel(name)) {
1358:         ostringstream msg;
1359:         msg << "Invalid label name: " << name << std::endl;
1360:         throw ALE::Exception(msg.str().c_str());
1361:       }
1362:     };
1363:   public: // Accessors
1364:     const Obj<sieve_type>& getSieve() const {return this->_sieve;};
1365:     void setSieve(const Obj<sieve_type>& sieve) {this->_sieve = sieve;};
1366:     bool hasArrowSection(const std::string& name) const {
1367:       return this->_arrowSections.find(name) != this->_arrowSections.end();
1368:     };
1369:     const Obj<arrow_section_type>& getArrowSection(const std::string& name) {
1370:       if (!this->hasArrowSection(name)) {
1371:         Obj<arrow_section_type> section = new arrow_section_type(this->comm(), this->debug());

1373:         section->setName(name);
1374:         if (this->_debug) {std::cout << "Creating new arrow section: " << name << std::endl;}
1375:         this->_arrowSections[name] = section;
1376:       }
1377:       return this->_arrowSections[name];
1378:     };
1379:     void setArrowSection(const std::string& name, const Obj<arrow_section_type>& section) {
1380:       this->_arrowSections[name] = section;
1381:     };
1382:     Obj<std::set<std::string> > getArrowSections() const {
1383:       Obj<std::set<std::string> > names = std::set<std::string>();

1385:       for(typename arrow_sections_type::const_iterator s_iter = this->_arrowSections.begin(); s_iter != this->_arrowSections.end(); ++s_iter) {
1386:         names->insert(s_iter->first);
1387:       }
1388:       return names;
1389:     };
1390:     bool hasRealSection(const std::string& name) const {
1391:       return this->_realSections.find(name) != this->_realSections.end();
1392:     };
1393:     const Obj<real_section_type>& getRealSection(const std::string& name) {
1394:       if (!this->hasRealSection(name)) {
1395:         Obj<real_section_type> section = new real_section_type(this->comm(), this->debug());

1397:         section->setName(name);
1398:         if (this->_debug) {std::cout << "Creating new real section: " << name << std::endl;}
1399:         this->_realSections[name] = section;
1400:       }
1401:       return this->_realSections[name];
1402:     };
1403:     void setRealSection(const std::string& name, const Obj<real_section_type>& section) {
1404:       this->_realSections[name] = section;
1405:     };
1406:     Obj<std::set<std::string> > getRealSections() const {
1407:       Obj<std::set<std::string> > names = std::set<std::string>();

1409:       for(typename real_sections_type::const_iterator s_iter = this->_realSections.begin(); s_iter != this->_realSections.end(); ++s_iter) {
1410:         names->insert(s_iter->first);
1411:       }
1412:       return names;
1413:     };
1414:     bool hasIntSection(const std::string& name) const {
1415:       return this->_intSections.find(name) != this->_intSections.end();
1416:     };
1417:     const Obj<int_section_type>& getIntSection(const std::string& name) {
1418:       if (!this->hasIntSection(name)) {
1419:         Obj<int_section_type> section = new int_section_type(this->comm(), this->debug());

1421:         section->setName(name);
1422:         if (this->_debug) {std::cout << "Creating new int section: " << name << std::endl;}
1423:         this->_intSections[name] = section;
1424:       }
1425:       return this->_intSections[name];
1426:     };
1427:     void setIntSection(const std::string& name, const Obj<int_section_type>& section) {
1428:       this->_intSections[name] = section;
1429:     };
1430:     Obj<std::set<std::string> > getIntSections() const {
1431:       Obj<std::set<std::string> > names = std::set<std::string>();

1433:       for(typename int_sections_type::const_iterator s_iter = this->_intSections.begin(); s_iter != this->_intSections.end(); ++s_iter) {
1434:         names->insert(s_iter->first);
1435:       }
1436:       return names;
1437:     };
1438:     const Obj<MeshNumberingFactory>& getFactory() const {return this->_factory;};
1439:     renumbering_type& getRenumbering() {return this->_renumbering;};
1440:   public: // Labels
1441:     int getValue (const Obj<label_type>& label, const point_type& point, const int defValue = 0) {
1442:       const Obj<typename label_type::coneSequence>& cone = label->cone(point);

1444:       if (cone->size() == 0) return defValue;
1445:       return *cone->begin();
1446:     };
1447:     void setValue(const Obj<label_type>& label, const point_type& point, const int value) {
1448:       label->setCone(value, point);
1449:     };
1450:     void addValue(const Obj<label_type>& label, const point_type& point, const int value) {
1451:       label->addCone(value, point);
1452:     };
1453:     template<typename InputPoints>
1454:     int getMaxValue (const Obj<label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
1455:       int maxValue = defValue;

1457:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1458:         maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
1459:       }
1460:       return maxValue;
1461:     };
1462:     const Obj<label_type>& createLabel(const std::string& name) {
1463:       this->_labels[name] = new label_type(this->comm(), this->debug());
1464:       return this->_labels[name];
1465:     };
1466:     const Obj<label_type>& getLabel(const std::string& name) {
1467:       this->checkLabel(name);
1468:       return this->_labels[name];
1469:     };
1470:     void setLabel(const std::string& name, const Obj<label_type>& label) {
1471:       this->_labels[name] = label;
1472:     };
1473:     const labels_type& getLabels() {
1474:       return this->_labels;
1475:     };
1476:     virtual const Obj<label_sequence>& getLabelStratum(const std::string& name, int value) {
1477:       this->checkLabel(name);
1478:       return this->_labels[name]->support(value);
1479:     };
1480:   public: // Stratification
1481:     void computeHeights() {
1482:       const Obj<label_type>& label = this->createLabel(std::string("height"));
1483:       HeightVisitor          h(*this->_sieve, *label);

1485: #ifdef IMESH_NEW_LABELS
1486:       label->setChart(this->getSieve()->getChart());
1487:       for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setConeSize(p, 1);}
1488:       if (label->getChart().size()) {label->setSupportSize(0, label->getChart().size());}
1489:       label->allocate();
1490:       for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setCone(-1, p);}
1491: #endif
1492:       this->_sieve->leaves(h);
1493:       while(h.isModified()) {
1494:         // FIX: Avoid the copy here somehow by fixing the traversal
1495:         std::vector<point_type> modifiedPoints(h.getModifiedPoints().begin(), h.getModifiedPoints().end());

1497:         h.clear();
1498:         this->_sieve->cone(modifiedPoints, h);
1499:       }
1500: #ifdef IMESH_NEW_LABELS
1501:       // Recalculate supportOffsets and populate support
1502:       label->recalculateLabel();
1503: #endif
1504:       this->_maxHeight = h.getMaxHeight();
1505:     };
1506:     virtual int height() const {return this->_maxHeight;};
1507:     virtual void setHeight(const int height) {this->_maxHeight = height;};
1508:     virtual int height(const point_type& point) {
1509:       return this->getValue(this->_labels["height"], point, -1);
1510:     };
1511:     virtual void setHeight(const point_type& point, const int height) {
1512:       return this->setValue(this->_labels["height"], point, height);
1513:     };
1514:     virtual const Obj<label_sequence>& heightStratum(int height) {
1515:       return this->getLabelStratum("height", height);
1516:     };
1517:     void computeDepths() {
1518:       const Obj<label_type>& label = this->createLabel(std::string("depth"));
1519:       DepthVisitor           d(*this->_sieve, *label);

1521: #ifdef IMESH_NEW_LABELS
1522:       label->setChart(this->getSieve()->getChart());
1523:       for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setConeSize(p, 1);}
1524:       if (label->getChart().size()) {label->setSupportSize(0, label->getChart().size());}
1525:       label->allocate();
1526:       for(point_type p = label->getChart().min(); p < label->getChart().max(); ++p) {label->setCone(-1, p);}
1527: #endif
1528:       this->_sieve->roots(d);
1529:       while(d.isModified()) {
1530:         // FIX: Avoid the copy here somehow by fixing the traversal
1531:         std::vector<point_type> modifiedPoints(d.getModifiedPoints().begin(), d.getModifiedPoints().end());

1533:         d.clear();
1534:         this->_sieve->support(modifiedPoints, d);
1535:       }
1536: #ifdef IMESH_NEW_LABELS
1537:       // Recalculate supportOffsets and populate support
1538:       label->recalculateLabel();
1539: #endif
1540:       this->_maxDepth = d.getMaxDepth();
1541:     };
1542:     virtual int depth() const {return this->_maxDepth;};
1543:     virtual void setDepth(const int depth) {this->_maxDepth = depth;};
1544:     virtual int depth(const point_type& point) {
1545:       return this->getValue(this->_labels["depth"], point, -1);
1546:     };
1547:     virtual void setDepth(const point_type& point, const int depth) {
1548:       return this->setValue(this->_labels["depth"], point, depth);
1549:     };
1550:     virtual const Obj<label_sequence>& depthStratum(int depth) {
1551:       return this->getLabelStratum("depth", depth);
1552:     };
1553:     void stratify() {
1554:       this->computeHeights();
1555:       this->computeDepths();
1556:     };
1557:   protected:
1558:     template<typename Value>
1559:     static bool lt1(const Value& a, const Value& b) {
1560:       return a.first < b.first;
1561:     };
1562:   public: // Allocation
1563:     template<typename Section_>
1564:     void reallocate(const Obj<Section_>& section) {
1565:       if (!section->hasNewPoints()) return;
1566:       typename Section_::chart_type newChart(std::min(std::min_element(section->getNewPoints().begin(), section->getNewPoints().end(), lt1<typename Section_::newpoint_type>)->first, section->getChart().min()),
1567:                                              std::max(std::max_element(section->getNewPoints().begin(), section->getNewPoints().end(), lt1<typename Section_::newpoint_type>)->first, section->getChart().max()-1)+1);
1568:       section->reallocatePoint(newChart);
1569:     };
1570:   };
1571: #ifdef IMESH_NEW_LABELS
1572:   template<typename Label_ = IFSieve<int> >
1573: #else
1574:   template<typename Label_ = LabelSifter<int, int> >
1575: #endif
1576:   class IMesh : public IBundle<IFSieve<int>, IGeneralSection<int, double>, IGeneralSection<int, int>,  Label_> {
1577:   public:
1578:     typedef IBundle<IFSieve<int>, IGeneralSection<int, double>, IGeneralSection<int, int>, Label_> base_type;
1579:     typedef typename base_type::sieve_type            sieve_type;
1580:     typedef typename sieve_type::point_type           point_type;
1581:     typedef typename base_type::alloc_type            alloc_type;
1582:     typedef typename base_type::label_type            label_type;
1583:     typedef typename base_type::label_sequence        label_sequence;
1584:     typedef typename base_type::real_section_type     real_section_type;
1585:     typedef typename base_type::numbering_type        numbering_type;
1586:     typedef typename base_type::order_type            order_type;
1587:     typedef typename base_type::send_overlap_type     send_overlap_type;
1588:     typedef typename base_type::recv_overlap_type     recv_overlap_type;
1589:     typedef std::set<std::string>                     names_type;
1590:     typedef std::vector<typename PETSc::Point<3> >    holes_type;
1591:     typedef std::map<std::string, Obj<Discretization> > discretizations_type;
1592:   protected:
1593:     int _dim;
1594:     bool                   _calculatedOverlap;
1595:     Obj<send_overlap_type> _sendOverlap;
1596:     Obj<recv_overlap_type> _recvOverlap;
1597:     std::map<int,double>   _periodicity;
1598:     holes_type             _holes;
1599:     discretizations_type   _discretizations;
1600:     int                    _maxDof;
1601:   public:
1602:     IMesh(MPI_Comm comm, int dim, int debug = 0) : base_type(comm, debug), _dim(dim) {
1603:       this->_calculatedOverlap = false;
1604:       this->_sendOverlap       = new send_overlap_type(comm, debug);
1605:       this->_recvOverlap       = new recv_overlap_type(comm, debug);
1606:       this->_maxDof            = -1;
1607:     };
1608:   public: // Accessors
1609:     int getDimension() const {return this->_dim;};
1610:     void setDimension(const int dim) {this->_dim = dim;};
1611:     bool getCalculatedOverlap() const {return this->_calculatedOverlap;};
1612:     void setCalculatedOverlap(const bool calc) {this->_calculatedOverlap = calc;};
1613:     const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
1614:     void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
1615:     const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
1616:     void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
1617:     bool getPeriodicity(const int d) {return this->_periodicity[d];};
1618:     void setPeriodicity(const int d, const double length) {this->_periodicity[d] = length;};
1619:     const holes_type& getHoles() const {return this->_holes;};
1620:     void addHole(const double hole[]) {
1621:       this->_holes.push_back(hole);
1622:     };
1623:     void copyHoles(const Obj<IMesh>& m) {
1624:       const holes_type& holes = m->getHoles();

1626:       for(holes_type::const_iterator h_iter = holes.begin(); h_iter != holes.end(); ++h_iter) {
1627:         this->_holes.push_back(*h_iter);
1628:       }
1629:     };
1630:     const Obj<Discretization>& getDiscretization() {return this->getDiscretization("default");};
1631:     const Obj<Discretization>& getDiscretization(const std::string& name) {return this->_discretizations[name];};
1632:     void setDiscretization(const Obj<Discretization>& disc) {this->setDiscretization("default", disc);};
1633:     void setDiscretization(const std::string& name, const Obj<Discretization>& disc) {this->_discretizations[name] = disc;};
1634:     Obj<names_type> getDiscretizations() const {
1635:       Obj<names_type> names = names_type();

1637:       for(discretizations_type::const_iterator d_iter = this->_discretizations.begin(); d_iter != this->_discretizations.end(); ++d_iter) {
1638:         names->insert(d_iter->first);
1639:       }
1640:       return names;
1641:     };
1642:     int getMaxDof() const {return this->_maxDof;};
1643:     void setMaxDof(const int maxDof) {this->_maxDof = maxDof;};
1644:   public: // Sizes
1645:     template<typename Section>
1646:     int size(const Obj<Section>& section, const point_type& p) {
1647:       typedef ISieveVisitor::SizeVisitor<sieve_type,Section>                        size_visitor_type;
1648:       typedef ISieveVisitor::TransitiveClosureVisitor<sieve_type,size_visitor_type> closure_visitor_type;
1649:       size_visitor_type    sV(*section);
1650:       closure_visitor_type cV(*this->getSieve(), sV);

1652:       this->getSieve()->cone(p, cV);
1653:       if (!sV.getSize()) sV.visitPoint(p);
1654:       return sV.getSize();
1655:     };
1656:     template<typename Section>
1657:     int sizeWithBC(const Obj<Section>& section, const point_type& p) {
1658:       typedef ISieveVisitor::SizeWithBCVisitor<sieve_type,Section>                  size_visitor_type;
1659:       typedef ISieveVisitor::TransitiveClosureVisitor<sieve_type,size_visitor_type> closure_visitor_type;
1660:       size_visitor_type    sV(*section);
1661:       closure_visitor_type cV(*this->getSieve(), sV);

1663:       this->getSieve()->cone(p, cV);
1664:       if (!sV.getSize()) sV.visitPoint(p);
1665:       return sV.getSize();
1666:     };
1667:     template<typename Section>
1668:     void allocate(const Obj<Section>& section) {
1669:       section->allocatePoint();
1670:     };
1671:   public: // Restrict/Update closures
1672:     template<typename Sieve, typename Visitor>
1673:     void closure1(const Sieve& sieve, const point_type& p, Visitor& v)
1674:     {
1675:       v.visitPoint(p, 0);
1676:       // Cone is guarateed to be ordered correctly
1677:       sieve.orientedCone(p, v);
1678:     };
1679:     // Return the values for the closure of this point
1680:     template<typename Section>
1681:     const typename Section::value_type *restrictClosure(const Obj<Section>& section, const point_type& p) {
1682:       const int size = this->sizeWithBC(section, p);
1683:       ISieveVisitor::RestrictVisitor<Section> rV(*section, size, section->getRawArray(size));

1685:       if (this->depth() == 1) {
1686:         closure1(*this->getSieve(), p, rV);
1687:       } else {
1688:         ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::RestrictVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, rV, true);

1690:         ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1691:       }
1692:       return rV.getValues();
1693:     };
1694:     template<typename Section>
1695:     const typename Section::value_type *restrictClosure(const Obj<Section>& section, const point_type& p, typename Section::value_type *values, const int valuesSize) {
1696:       const int size = this->sizeWithBC(section, p);
1697:       if (valuesSize < size) {throw ALE::Exception("Input array to small for restrictClosure()");}
1698:       ISieveVisitor::RestrictVisitor<Section> rV(*section, size, values);

1700:       if (this->depth() == 1) {
1701:         closure1(*this->getSieve(), p, rV);
1702:       } else {
1703:         ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::RestrictVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, rV, true);

1705:         ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1706:       }
1707:       return rV.getValues();
1708:     };
1709:     template<typename Visitor>
1710:     void restrictClosure(const point_type& p, Visitor& v) {
1711:       if (this->depth() == 1) {
1712:         closure1(*this->getSieve(), p, v);
1713:       } else {
1714:         ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, v);
1715:       }
1716:     };
1717:     // Replace the values for the closure of this point
1718:     template<typename Section>
1719:     void update(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1720:       ISieveVisitor::UpdateVisitor<Section> uV(*section, v);

1722:       if (this->depth() == 1) {
1723:         closure1(*this->getSieve(), p, uV);
1724:       } else {
1725:         ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);

1727:         ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1728:       }
1729:     };
1730:     // Replace the values for the closure of this point, including points constrained by BC
1731:     template<typename Section>
1732:     void updateAll(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1733:       ISieveVisitor::UpdateAllVisitor<Section> uV(*section, v);

1735:       if (this->depth() == 1) {
1736:         closure1(*this->getSieve(), p, uV);
1737:       } else {
1738:         ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateAllVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);

1740:         ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1741:       }
1742:     };
1743:     // Augment the values for the closure of this point
1744:     template<typename Section>
1745:     void updateAdd(const Obj<Section>& section, const point_type& p, const typename Section::value_type *v) {
1746:       ISieveVisitor::UpdateAddVisitor<Section> uV(*section, v);

1748:       if (this->depth() == 1) {
1749:         closure1(*this->getSieve(), p, uV);
1750:       } else {
1751:         ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::UpdateAddVisitor<Section> > pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, uV, true);

1753:         ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, pV);
1754:       }
1755:     };
1756:     // Augment the values for the closure of this point
1757:     template<typename Visitor>
1758:     void updateClosure(const point_type& p, Visitor& v) {
1759:       if (this->depth() == 1) {
1760:         closure1(*this->getSieve(), p, v);
1761:       } else {
1762:         ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), p, v);
1763:       }
1764:     };
1765:   public: // Overlap
1766:     void constructOverlap() {
1767:       if (!this->_calculatedOverlap && (this->commSize() > 1)) {throw ALE::Exception("Must calculate overlap during distribution");}
1768:     };
1769:   public: // Cell topology and geometry
1770:     int getNumCellCorners(const point_type& p, const int depth = -1) const {
1771:       const int d = depth < 0 ? this->depth() : depth;

1773:       if (d == 1) {
1774:         return this->_sieve->getConeSize(p);
1775:       } else if (d <= 0) {
1776:         return 0;
1777:       }
1778:       // Warning: this is slow
1779:       ISieveVisitor::NConeRetriever<sieve_type> ncV(*this->_sieve, (int) pow((double) this->_sieve->getMaxConeSize(), this->depth()));
1780:       ALE::ISieveTraversal<sieve_type>::orientedClosure(*this->_sieve, p, ncV);
1781:       return ncV.getOrientedSize();
1782:     };
1783:     int getNumCellCorners() {
1784:       return getNumCellCorners(*this->heightStratum(0)->begin());
1785:     };
1786:     void setupCoordinates(const Obj<real_section_type>& coordinates) {
1787:       const Obj<label_sequence>& vertices = this->depthStratum(0);

1789:       coordinates->setChart(typename real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()),
1790:                                                                    *std::max_element(vertices->begin(), vertices->end())+1));
1791:     };
1792:     point_type locatePoint(const typename real_section_type::value_type point[], point_type guess = -1) {
1793:       //guess overrides this by saying that we already know the relation of this point to this mesh.  We will need to make it a more robust "guess" later for more than P1
1794:       if (guess != -1) {
1795:         return guess;
1796:       } else {
1797:         throw ALE::Exception("No point location for mesh dimension");
1798:       }
1799:     };
1800:     void computeTriangleGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1801:       const double *coords = this->restrictClosure(coordinates, e);
1802:       const int     dim    = 2;
1803:       double        invDet;

1805:       if (v0) {
1806:         for(int d = 0; d < dim; d++) {
1807:           v0[d] = coords[d];
1808:         }
1809:       }
1810:       if (J) {
1811:         for(int d = 0; d < dim; d++) {
1812:           for(int f = 0; f < dim; f++) {
1813:             J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1814:           }
1815:         }
1816:         detJ = J[0]*J[3] - J[1]*J[2];
1817:         if (detJ < 0.0) {
1818:           const double  xLength = this->_periodicity[0];

1820:           if (xLength != 0.0) {
1821:             double v0x = coords[0*dim+0];

1823:             if (v0x == 0.0) {
1824:               v0x = v0[0] = xLength;
1825:             }
1826:             for(int f = 0; f < dim; f++) {
1827:               const double px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];

1829:               J[0*dim+f] = 0.5*(px - v0x);
1830:             }
1831:           }
1832:           detJ = J[0]*J[3] - J[1]*J[2];
1833:         }
1834:         PetscLogFlopsNoError(8.0 + 3.0);
1835:       }
1836:       if (invJ) {
1837:         invDet  = 1.0/detJ;
1838:         invJ[0] =  invDet*J[3];
1839:         invJ[1] = -invDet*J[1];
1840:         invJ[2] = -invDet*J[2];
1841:         invJ[3] =  invDet*J[0];
1842:         PetscLogFlopsNoError(5.0);
1843:       }
1844:     };
1845:     void computeTetrahedronGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1846:       const double *coords = this->restrictClosure(coordinates, e);
1847:       const int     dim    = 3;
1848:       double        invDet;

1850:       if (v0) {
1851:         for(int d = 0; d < dim; d++) {
1852:           v0[d] = coords[d];
1853:         }
1854:       }
1855:       if (J) {
1856:         for(int d = 0; d < dim; d++) {
1857:           for(int f = 0; f < dim; f++) {
1858:             J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1859:           }
1860:         }
1861:         // The minus sign is here since I orient the first face to get the outward normal
1862:         detJ = -(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
1863:                  J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
1864:                  J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
1865:         PetscLogFlopsNoError(18.0 + 12.0);
1866:       }
1867:       if (invJ) {
1868:         invDet  = -1.0/detJ;
1869:         invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
1870:         invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
1871:         invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
1872:         invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
1873:         invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
1874:         invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
1875:         invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1876:         invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
1877:         invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
1878:         PetscLogFlopsNoError(37.0);
1879:       }
1880:     };
1881:     void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1882:       if (this->_dim == 2) {
1883:         computeTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
1884:       } else if (this->_dim == 3) {
1885:         computeTetrahedronGeometry(coordinates, e, v0, J, invJ, detJ);
1886:       } else {
1887:         throw ALE::Exception("Unsupported dimension for element geometry computation");
1888:       }
1889:     };
1890:     void computeBdSegmentGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1891:       const double *coords = this->restrictClosure(coordinates, e);
1892:       const int     dim    = 2;
1893:       double        invDet;

1895:       if (v0) {
1896:         for(int d = 0; d < dim; d++) {
1897:           v0[d] = coords[d];
1898:         }
1899:       }
1900:       if (J) {
1901:         //r2   = coords[1*dim+0]*coords[1*dim+0] + coords[1*dim+1]*coords[1*dim+1];
1902:         J[0] =  (coords[1*dim+0] - coords[0*dim+0])*0.5; J[1] = (-coords[1*dim+1] + coords[0*dim+1])*0.5;
1903:         J[2] =  (coords[1*dim+1] - coords[0*dim+1])*0.5; J[3] = ( coords[1*dim+0] - coords[0*dim+0])*0.5;
1904:         detJ = J[0]*J[3] - J[1]*J[2];
1905:         if (detJ < 0.0) {
1906:           const double  xLength = this->_periodicity[0];

1908:           if (xLength != 0.0) {
1909:             double v0x = coords[0*dim+0];

1911:             if (v0x == 0.0) {
1912:               v0x = v0[0] = xLength;
1913:             }
1914:             for(int f = 0; f < dim; f++) {
1915:               const double px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];

1917:               J[0*dim+f] = 0.5*(px - v0x);
1918:             }
1919:           }
1920:           detJ = J[0]*J[3] - J[1]*J[2];
1921:         }
1922:         PetscLogFlopsNoError(8.0 + 3.0);
1923:       }
1924:       if (invJ) {
1925:         invDet  = 1.0/detJ;
1926:         invJ[0] =  invDet*J[3];
1927:         invJ[1] = -invDet*J[1];
1928:         invJ[2] = -invDet*J[2];
1929:         invJ[3] =  invDet*J[0];
1930:         PetscLogFlopsNoError(5.0);
1931:       }
1932:     };
1933:     void computeBdElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
1934:       if (this->_dim == 1) {
1935:         computeBdSegmentGeometry(coordinates, e, v0, J, invJ, detJ);
1936:         //      } else if (this->_dim == 2) {
1937:         //        computeBdTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
1938:       } else {
1939:         throw ALE::Exception("Unsupported dimension for element geometry computation");
1940:       }
1941:     };
1942:     double getMaxVolume() {
1943:       const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
1944:       const Obj<label_sequence>&    cells       = this->heightStratum(0);
1945:       const int                     dim         = this->getDimension();
1946:       double v0[3], J[9], invJ[9], detJ, refVolume = 0.0, maxVolume = 0.0;

1948:       if (dim == 1) refVolume = 2.0;
1949:       if (dim == 2) refVolume = 2.0;
1950:       if (dim == 3) refVolume = 4.0/3.0;
1951:       for(typename label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1952:         this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
1953:         maxVolume = std::max(maxVolume, detJ*refVolume);
1954:       }
1955:       return maxVolume;
1956:     };
1957:   public:
1958:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1959:       if (comm == MPI_COMM_NULL) {
1960:         comm = this->comm();
1961:       }
1962:       if (name == "") {
1963:         PetscPrintf(comm, "viewing a Mesh\n");
1964:       } else {
1965:         PetscPrintf(comm, "viewing Mesh '%s'\n", name.c_str());
1966:       }
1967:       this->getSieve()->view("mesh sieve", comm);
1968:       Obj<names_type> sections = this->getRealSections();

1970:       for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1971:         this->getRealSection(*name)->view(*name);
1972:       }
1973:       sections = this->getIntSections();
1974:       for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1975:         this->getIntSection(*name)->view(*name);
1976:       }
1977:       sections = this->getArrowSections();
1978:       for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
1979:         this->getArrowSection(*name)->view(*name);
1980:       }
1981:     };
1982:   public: // Discretization
1983:     void markBoundaryCells(const std::string& name, const int marker = 1, const int newMarker = 2, const bool onlyVertices = false) {
1984:       const Obj<label_type>&                  label    = this->getLabel(name);
1985:       const Obj<label_sequence>&              boundary = this->getLabelStratum(name, marker);
1986:       const typename label_sequence::iterator end      = boundary->end();
1987:       const Obj<sieve_type>&                  sieve    = this->getSieve();

1989:       if (!onlyVertices) {
1990:         typename ISieveVisitor::MarkVisitor<sieve_type,label_type> mV(*label, newMarker);

1992:         for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != end; ++e_iter) {
1993:           if (this->height(*e_iter) == 1) {
1994:             sieve->support(*e_iter, mV);
1995:           }
1996:         }
1997:       } else {
1998: #if 1
1999:         throw ALE::Exception("Rewrite this to first mark bounadry edges/faces.");
2000: #else
2001:         const int depth = this->depth();

2003:         for(typename label_sequence::iterator v_iter = boundary->begin(); v_iter != end; ++v_iter) {
2004:           const Obj<supportArray>               support = sieve->nSupport(*v_iter, depth);
2005:           const typename supportArray::iterator sEnd    = support->end();

2007:           for(typename supportArray::iterator c_iter = support->begin(); c_iter != sEnd; ++c_iter) {
2008:             const Obj<typename sieve_type::traits::coneSequence>&     cone = sieve->cone(*c_iter);
2009:             const typename sieve_type::traits::coneSequence::iterator cEnd = cone->end();

2011:             for(typename sieve_type::traits::coneSequence::iterator e_iter = cone->begin(); e_iter != cEnd; ++e_iter) {
2012:               if (sieve->support(*e_iter)->size() == 1) {
2013:                 this->setValue(label, *c_iter, newMarker);
2014:                 break;
2015:               }
2016:             }
2017:           }
2018:         }
2019: #endif
2020:       }
2021:     };
2022:     int setFiberDimensions(const Obj<real_section_type>& s, const Obj<names_type>& discs, names_type& bcLabels) {
2023:       const int debug  = this->debug();
2024:       int       maxDof = 0;

2026:       for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
2027:         s->addSpace();
2028:       }
2029:       for(int d = 0; d <= this->_dim; ++d) {
2030:         int numDof = 0;
2031:         int f      = 0;

2033:         for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2034:           const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2035:           const int                       sDof = disc->getNumDof(d);

2037:           numDof += sDof;
2038:           if (sDof) s->setFiberDimension(this->depthStratum(d), sDof, f);
2039:         }
2040:         if (numDof) s->setFiberDimension(this->depthStratum(d), numDof);
2041:         maxDof = std::max(maxDof, numDof);
2042:       }
2043:       // Process exclusions
2044:       typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2045:       int f = 0;

2047:       for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2048:         const Obj<ALE::Discretization>& disc      = this->getDiscretization(*f_iter);
2049:         std::string                     labelName = "exclude-"+*f_iter;
2050:         std::set<point_type>            seen;
2051:         Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth()), true);

2053:         if (this->hasLabel(labelName)) {
2054:           const Obj<label_type>&                  label     = this->getLabel(labelName);
2055:           const Obj<label_sequence>&              exclusion = this->getLabelStratum(labelName, 1);
2056:           const typename label_sequence::iterator end       = exclusion->end();
2057:           if (debug > 1) {label->view(labelName.c_str());}

2059:           for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2060:             ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *e_iter, pV);
2061:             const typename Visitor::point_type *oPoints = pV.getPoints();
2062:             const int                           oSize   = pV.getSize();

2064:             for(int cl = 0; cl < oSize; ++cl) {
2065:               if (seen.find(oPoints[cl]) != seen.end()) continue;
2066:               if (this->getValue(label, oPoints[cl]) == 1) {
2067:                 seen.insert(oPoints[cl]);
2068:                 s->setFiberDimension(oPoints[cl], 0, f);
2069:                 s->addFiberDimension(oPoints[cl], -disc->getNumDof(this->depth(oPoints[cl])));
2070:                 if (debug > 1) {std::cout << "  point: " << oPoints[cl] << " dim: " << disc->getNumDof(this->depth(oPoints[cl])) << std::endl;}
2071:               }
2072:             }
2073:             pV.clear();
2074:           }
2075:         }
2076:       }
2077:       // Process constraints
2078:       f = 0;
2079:       for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2080:         const Obj<ALE::Discretization>&   disc        = this->getDiscretization(*f_iter);
2081:         const Obj<std::set<std::string> > bcs         = disc->getBoundaryConditions();
2082:         std::string                       excludeName = "exclude-"+*f_iter;

2084:         for(std::set<std::string>::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter) {
2085:           const Obj<ALE::BoundaryCondition>& bc       = disc->getBoundaryCondition(*bc_iter);
2086:           const Obj<label_sequence>&         boundary = this->getLabelStratum(bc->getLabelName(), bc->getMarker());

2088:           bcLabels.insert(bc->getLabelName());
2089:           if (this->hasLabel(excludeName)) {
2090:             const Obj<label_type>& label = this->getLabel(excludeName);

2092:             for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2093:               if (!this->getValue(label, *e_iter)) {
2094:                 const int numDof = disc->getNumDof(this->depth(*e_iter));

2096:                 if (numDof) s->addConstraintDimension(*e_iter, numDof);
2097:                 if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2098:               }
2099:             }
2100:           } else {
2101:             for(typename label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2102:               const int numDof = disc->getNumDof(this->depth(*e_iter));

2104:               if (numDof) s->addConstraintDimension(*e_iter, numDof);
2105:               if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2106:             }
2107:           }
2108:         }
2109:       }
2110:       return maxDof;
2111:     };
2112:     void calculateIndices() {
2113:       typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2114:       // Should have an iterator over the whole tree
2115:       Obj<names_type> discs = this->getDiscretizations();
2116:       const int       debug = this->debug();
2117:       std::map<std::string, std::pair<int, int*> > indices;

2119:       for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2120:         const Obj<Discretization>& disc = this->getDiscretization(*d_iter);

2122:         indices[*d_iter] = std::pair<int, int*>(0, new int[disc->sizeV(*this)]);
2123:         disc->setIndices(indices[*d_iter].second);
2124:       }
2125:       const Obj<label_sequence>& cells   = this->heightStratum(0);
2126:       Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);
2127:       ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *cells->begin(), pV);
2128:       const typename Visitor::point_type *oPoints = pV.getPoints();
2129:       const int                           oSize   = pV.getSize();
2130:       int                                 offset  = 0;

2132:       if (debug > 1) {std::cout << "Closure for first element" << std::endl;}
2133:       for(int cl = 0; cl < oSize; ++cl) {
2134:         const int dim = this->depth(oPoints[cl]);

2136:         if (debug > 1) {std::cout << "  point " << oPoints[cl] << " depth " << dim << std::endl;}
2137:         for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2138:           const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2139:           const int                  num  = disc->getNumDof(dim);

2141:           if (debug > 1) {std::cout << "    disc " << disc->getName() << " numDof " << num << std::endl;}
2142:           for(int o = 0; o < num; ++o) {
2143:             indices[*d_iter].second[indices[*d_iter].first++] = offset++;
2144:           }
2145:         }
2146:       }
2147:       pV.clear();
2148:       if (debug > 1) {
2149:         for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2150:           const Obj<Discretization>& disc = this->getDiscretization(*d_iter);

2152:           std::cout << "Discretization " << disc->getName() << " indices:";
2153:           for(int i = 0; i < indices[*d_iter].first; ++i) {
2154:             std::cout << " " << indices[*d_iter].second[i];
2155:           }
2156:           std::cout << std::endl;
2157:         }
2158:       }
2159:     };
2160:     void calculateIndicesExcluded(const Obj<real_section_type>& s, const Obj<names_type>& discs) {
2161:       typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2162:       typedef std::map<std::string, std::pair<int, indexSet> > indices_type;
2163:       const Obj<label_type>& indexLabel = this->createLabel("cellExclusion");
2164:       const int debug  = this->debug();
2165:       int       marker = 0;
2166:       std::map<indices_type, int> indexMap;
2167:       indices_type                indices;
2168:       Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);

2170:       for(typename names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
2171:         const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
2172:         const int                  size = disc->sizeV(*this);

2174:         indices[*d_iter].second.resize(size);
2175:       }
2176:       const typename names_type::const_iterator dBegin = discs->begin();
2177:       const typename names_type::const_iterator dEnd   = discs->end();
2178:       std::set<point_type> seen;
2179:       int f = 0;

2181:       for(typename names_type::const_iterator f_iter = dBegin; f_iter != dEnd; ++f_iter, ++f) {
2182:         std::string labelName = "exclude-"+*f_iter;

2184:         if (this->hasLabel(labelName)) {
2185:           const Obj<label_sequence>&              exclusion = this->getLabelStratum(labelName, 1);
2186:           const typename label_sequence::iterator end       = exclusion->end();

2188:           if (debug > 1) {std::cout << "Processing exclusion " << labelName << std::endl;}
2189:           for(typename label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2190:             if (this->height(*e_iter)) continue;
2191:             ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *e_iter, pV);
2192:             const typename Visitor::point_type *oPoints = pV.getPoints();
2193:             const int                           oSize   = pV.getSize();
2194:             int                                 offset  = 0;

2196:             if (debug > 1) {std::cout << "  Closure for cell " << *e_iter << std::endl;}
2197:             for(int cl = 0; cl < oSize; ++cl) {
2198:               int g = 0;

2200:               if (debug > 1) {std::cout << "    point " << oPoints[cl] << std::endl;}
2201:               for(typename names_type::const_iterator g_iter = dBegin; g_iter != dEnd; ++g_iter, ++g) {
2202:                 const int fDim = s->getFiberDimension(oPoints[cl], g);

2204:                 if (debug > 1) {std::cout << "      disc " << *g_iter << " numDof " << fDim << std::endl;}
2205:                 for(int d = 0; d < fDim; ++d) {
2206:                   indices[*g_iter].second[indices[*g_iter].first++] = offset++;
2207:                 }
2208:               }
2209:             }
2210:             pV.clear();
2211:             const std::map<indices_type, int>::iterator entry = indexMap.find(indices);

2213:             if (debug > 1) {
2214:               for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
2215:                 for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2216:                   std::cout << "Discretization (" << i_iter->second << ") " << *g_iter << " indices:";
2217:                   for(int i = 0; i < ((indices_type) i_iter->first)[*g_iter].first; ++i) {
2218:                     std::cout << " " << ((indices_type) i_iter->first)[*g_iter].second[i];
2219:                   }
2220:                   std::cout << std::endl;
2221:                 }
2222:                 std::cout << "Comparison: " << (indices == i_iter->first) << std::endl;
2223:               }
2224:               for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2225:                 std::cout << "Discretization " << *g_iter << " indices:";
2226:                 for(int i = 0; i < indices[*g_iter].first; ++i) {
2227:                   std::cout << " " << indices[*g_iter].second[i];
2228:                 }
2229:                 std::cout << std::endl;
2230:               }
2231:             }
2232:             if (entry != indexMap.end()) {
2233:               this->setValue(indexLabel, *e_iter, entry->second);
2234:               if (debug > 1) {std::cout << "  Found existing indices with marker " << entry->second << std::endl;}
2235:             } else {
2236:               indexMap[indices] = ++marker;
2237:               this->setValue(indexLabel, *e_iter, marker);
2238:               if (debug > 1) {std::cout << "  Created new indices with marker " << marker << std::endl;}
2239:             }
2240:             for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2241:               indices[*g_iter].first  = 0;
2242:               for(unsigned int i = 0; i < indices[*g_iter].second.size(); ++i) indices[*g_iter].second[i] = 0;
2243:             }
2244:           }
2245:         }
2246:       }
2247:       if (debug > 1) {indexLabel->view("cellExclusion");}
2248:       for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
2249:         if (debug > 1) {std::cout << "Setting indices for marker " << i_iter->second << std::endl;}
2250:         for(typename names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
2251:           const Obj<Discretization>& disc = this->getDiscretization(*g_iter);
2252:           const indexSet  indSet   = ((indices_type) i_iter->first)[*g_iter].second;
2253:           const int       size     = indSet.size();
2254:           int            *_indices = new int[size];

2256:           if (debug > 1) {std::cout << "  field " << *g_iter << std::endl;}
2257:           for(int i = 0; i < size; ++i) {
2258:             _indices[i] = indSet[i];
2259:             if (debug > 1) {std::cout << "    indices["<<i<<"] = " << _indices[i] << std::endl;}
2260:           }
2261:           disc->setIndices(_indices, i_iter->second);
2262:         }
2263:       }
2264:     };
2265:     void setupField(const Obj<real_section_type>& s, const int cellMarker = 2, const bool noUpdate = false) {
2266:       typedef ISieveVisitor::PointRetriever<sieve_type> Visitor;
2267:       const Obj<names_type>& discs  = this->getDiscretizations();
2268:       const int              debug  = s->debug();
2269:       names_type             bcLabels;

2271:       s->setChart(this->getSieve()->getChart());
2272:       this->_maxDof = this->setFiberDimensions(s, discs, bcLabels);
2273:       this->calculateIndices();
2274:       this->calculateIndicesExcluded(s, discs);
2275:       this->allocate(s);
2276:       s->defaultConstraintDof();
2277:       const Obj<label_type>& cellExclusion = this->getLabel("cellExclusion");

2279:       if (debug > 1) {std::cout << "Setting boundary values" << std::endl;}
2280:       for(typename names_type::const_iterator n_iter = bcLabels.begin(); n_iter != bcLabels.end(); ++n_iter) {
2281:         const Obj<label_sequence>&              boundaryCells = this->getLabelStratum(*n_iter, cellMarker);
2282:         const Obj<real_section_type>&           coordinates   = this->getRealSection("coordinates");
2283:         const Obj<names_type>&                  discs         = this->getDiscretizations();
2284:         const point_type                        firstCell     = *boundaryCells->begin();
2285:         const int                               numFields     = discs->size();
2286:         typename real_section_type::value_type *values        = new typename real_section_type::value_type[this->sizeWithBC(s, firstCell)];
2287:         int                                    *dofs          = new int[this->_maxDof];
2288:         int                                    *v             = new int[numFields];
2289:         double                                 *v0            = new double[this->getDimension()];
2290:         double                                 *J             = new double[this->getDimension()*this->getDimension()];
2291:         double                                  detJ;
2292:         Visitor pV((int) pow((double) this->getSieve()->getMaxConeSize(), this->depth())+1, true);

2294:         for(typename label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
2295:           ISieveTraversal<sieve_type>::orientedClosure(*this->getSieve(), *c_iter, pV);
2296:           const typename Visitor::point_type *oPoints = pV.getPoints();
2297:           const int                           oSize   = pV.getSize();

2299:           if (debug > 1) {std::cout << "  Boundary cell " << *c_iter << std::endl;}
2300:           this->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
2301:           for(int f = 0; f < numFields; ++f) v[f] = 0;
2302:           for(int cl = 0; cl < oSize; ++cl) {
2303:             const int cDim = s->getConstraintDimension(oPoints[cl]);
2304:             int       off  = 0;
2305:             int       f    = 0;
2306:             int       i    = -1;

2308:             if (debug > 1) {std::cout << "    point " << oPoints[cl] << std::endl;}
2309:             if (cDim) {
2310:               if (debug > 1) {std::cout << "      constrained excMarker: " << this->getValue(cellExclusion, *c_iter) << std::endl;}
2311:               for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2312:                 const Obj<ALE::Discretization>& disc    = this->getDiscretization(*f_iter);
2313:                 const Obj<names_type>           bcs     = disc->getBoundaryConditions();
2314:                 const int                       fDim    = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
2315:                 const int                      *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
2316:                 int                             b       = 0;

2318:                 for(typename names_type::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter, ++b) {
2319:                   const Obj<ALE::BoundaryCondition>& bc    = disc->getBoundaryCondition(*bc_iter);
2320:                   const int                          value = this->getValue(this->getLabel(bc->getLabelName()), oPoints[cl]);

2322:                   if (b > 0) v[f] -= fDim;
2323:                   if (value == bc->getMarker()) {
2324:                     if (debug > 1) {std::cout << "      field " << *f_iter << " marker " << value << std::endl;}
2325:                     for(int d = 0; d < fDim; ++d, ++v[f]) {
2326:                       dofs[++i] = off+d;
2327:                       if (!noUpdate) values[indices[v[f]]] = (*bc->getDualIntegrator())(v0, J, v[f], bc->getFunction());
2328:                       if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2329:                     }
2330:                     // Allow only one condition per point
2331:                     ++b;
2332:                     break;
2333:                   } else {
2334:                     if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
2335:                     for(int d = 0; d < fDim; ++d, ++v[f]) {
2336:                       values[indices[v[f]]] = 0.0;
2337:                       if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2338:                     }
2339:                   }
2340:                 }
2341:                 if (b == 0) {
2342:                   if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
2343:                   for(int d = 0; d < fDim; ++d, ++v[f]) {
2344:                     values[indices[v[f]]] = 0.0;
2345:                     if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2346:                   }
2347:                 }
2348:                 off += fDim;
2349:               }
2350:               if (i != cDim-1) {throw ALE::Exception("Invalid constraint initialization");}
2351:               s->setConstraintDof(oPoints[cl], dofs);
2352:             } else {
2353:               if (debug > 1) {std::cout << "      unconstrained" << std::endl;}
2354:               for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2355:                 const Obj<ALE::Discretization>& disc    = this->getDiscretization(*f_iter);
2356:                 const int                       fDim    = s->getFiberDimension(oPoints[cl], f);//disc->getNumDof(this->depth(oPoints[cl]));
2357:                 const int                      *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));

2359:                 if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
2360:                 for(int d = 0; d < fDim; ++d, ++v[f]) {
2361:                   values[indices[v[f]]] = 0.0;
2362:                   if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
2363:                 }
2364:               }
2365:             }
2366:           }
2367:           if (debug > 1) {
2368:             for(int f = 0; f < numFields; ++f) v[f] = 0;
2369:             for(int cl = 0; cl < oSize; ++cl) {
2370:               int f = 0;
2371:               for(typename names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2372:                 const Obj<ALE::Discretization>& disc    = this->getDiscretization(*f_iter);
2373:                 const int                       fDim    = s->getFiberDimension(oPoints[cl], f);
2374:                 const int                      *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));

2376:                 for(int d = 0; d < fDim; ++d, ++v[f]) {
2377:                   std::cout << "    "<<*f_iter<<"-value["<<indices[v[f]]<<"] " << values[indices[v[f]]] << std::endl;
2378:                 }
2379:               }
2380:             }
2381:           }
2382:           if (!noUpdate) {
2383:             this->updateAll(s, *c_iter, values);
2384:           }
2385:           pV.clear();
2386:         }
2387:         delete [] dofs;
2388:         delete [] values;
2389:         delete [] v0;
2390:         delete [] J;
2391:       }
2392:       if (debug > 1) {s->view("");}
2393:     };
2394:   };
2395: }

2397: namespace ALE {
2398:   class Mesh : public Bundle<ALE::Sieve<int,int,int>, GeneralSection<int, double> > {
2399:   public:
2400:     typedef Bundle<ALE::Sieve<int,int,int>, GeneralSection<int, double> > base_type;
2401:     typedef base_type::sieve_type            sieve_type;
2402:     typedef sieve_type::point_type           point_type;
2403:     typedef base_type::alloc_type            alloc_type;
2404:     typedef base_type::label_sequence        label_sequence;
2405:     typedef base_type::real_section_type     real_section_type;
2406:     typedef base_type::numbering_type        numbering_type;
2407:     typedef base_type::order_type            order_type;
2408:     typedef base_type::send_overlap_type     send_overlap_type;
2409:     typedef base_type::recv_overlap_type     recv_overlap_type;
2410:     typedef base_type::sieve_alg_type        sieve_alg_type;
2411:     typedef std::set<std::string>            names_type;
2412:     typedef std::map<std::string, Obj<Discretization> > discretizations_type;
2413:     typedef std::vector<PETSc::Point<3> >    holes_type;
2414:   protected:
2415:     int                  _dim;
2416:     discretizations_type _discretizations;
2417:     std::map<int,double> _periodicity;
2418:     holes_type           _holes;
2419:   public:
2420:     Mesh(MPI_Comm comm, int dim, int debug = 0) : base_type(comm, debug), _dim(dim) {
2421:       ///this->_factory = MeshNumberingFactory::singleton(debug);
2422:       //std::cout << "["<<this->commRank()<<"]: Creating an ALE::Mesh" << std::endl;
2423:     };
2424:     ~Mesh() {
2425:       //std::cout << "["<<this->commRank()<<"]: Destroying an ALE::Mesh" << std::endl;
2426:     };
2427:   public: // Accessors
2428:     int getDimension() const {return this->_dim;};
2429:     void setDimension(const int dim) {this->_dim = dim;};
2430:     const Obj<Discretization>& getDiscretization() {return this->getDiscretization("default");};
2431:     const Obj<Discretization>& getDiscretization(const std::string& name) {return this->_discretizations[name];};
2432:     void setDiscretization(const Obj<Discretization>& disc) {this->setDiscretization("default", disc);};
2433:     void setDiscretization(const std::string& name, const Obj<Discretization>& disc) {this->_discretizations[name] = disc;};
2434:     Obj<names_type> getDiscretizations() const {
2435:       Obj<names_type> names = names_type();

2437:       for(discretizations_type::const_iterator d_iter = this->_discretizations.begin(); d_iter != this->_discretizations.end(); ++d_iter) {
2438:         names->insert(d_iter->first);
2439:       }
2440:       return names;
2441:     };
2442:     bool getPeriodicity(const int d) {return this->_periodicity[d];};
2443:     void setPeriodicity(const int d, const double length) {this->_periodicity[d] = length;};
2444:     const holes_type& getHoles() const {return this->_holes;};
2445:     void addHole(const double hole[]) {
2446:       this->_holes.push_back(hole);
2447:     };
2448:     void copyHoles(const Obj<Mesh>& m) {
2449:       const holes_type& holes = m->getHoles();

2451:       for(holes_type::const_iterator h_iter = holes.begin(); h_iter != holes.end(); ++h_iter) {
2452:         this->_holes.push_back(*h_iter);
2453:       }
2454:     };
2455:     void copy(const Obj<Mesh>& m) {
2456:       this->setSieve(m->getSieve());
2457:       this->setLabel("height", m->getLabel("height"));
2458:       this->_maxHeight = m->height();
2459:       this->setLabel("depth", m->getLabel("depth"));
2460:       this->_maxDepth  = m->depth();
2461:       this->setLabel("marker", m->getLabel("marker"));
2462:       this->setRealSection("coordinates", m->getRealSection("coordinates"));
2463:       this->setArrowSection("orientation", m->getArrowSection("orientation"));
2464:     };
2465:   public: // Cell topology and geometry
2466:     int getNumCellCorners(const point_type& p, const int depth = -1) const {
2467:       return (this->getDimension() > 0) ? this->_sieve->nCone(p, depth < 0 ? this->depth() : depth)->size() : 1;
2468:     };
2469:     int getNumCellCorners() {
2470:       return getNumCellCorners(*(this->heightStratum(0)->begin()));
2471:     };
2472:     void setupCoordinates(const Obj<real_section_type>& coordinates) {};
2473:     void computeTriangleGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2474:       const double *coords = this->restrictClosure(coordinates, e);
2475:       const int     dim    = 2;
2476:       double        invDet;

2478:       if (v0) {
2479:         for(int d = 0; d < dim; d++) {
2480:           v0[d] = coords[d];
2481:         }
2482:       }
2483:       if (J) {
2484:         for(int d = 0; d < dim; d++) {
2485:           for(int f = 0; f < dim; f++) {
2486:             J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
2487:           }
2488:         }
2489:         detJ = J[0]*J[3] - J[1]*J[2];
2490:         if (detJ < 0.0) {
2491:           const double  xLength = this->_periodicity[0];

2493:           if (xLength != 0.0) {
2494:             double v0x = coords[0*dim+0];

2496:             if (v0x == 0.0) {
2497:               v0x = v0[0] = xLength;
2498:             }
2499:             for(int f = 0; f < dim; f++) {
2500:               const double px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];

2502:               J[0*dim+f] = 0.5*(px - v0x);
2503:             }
2504:           }
2505:           detJ = J[0]*J[3] - J[1]*J[2];
2506:         }
2507:         PetscLogFlopsNoError(8.0 + 3.0);
2508:       }
2509:       if (invJ) {
2510:         invDet  = 1.0/detJ;
2511:         invJ[0] =  invDet*J[3];
2512:         invJ[1] = -invDet*J[1];
2513:         invJ[2] = -invDet*J[2];
2514:         invJ[3] =  invDet*J[0];
2515:         PetscLogFlopsNoError(5.0);
2516:       }
2517:     };
2518:     void computeQuadrilateralGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double point[], double v0[], double J[], double invJ[], double& detJ) {
2519:       const double *coords = this->restrictClosure(coordinates, e);
2520:       const int     dim    = 2;
2521:       double        invDet;

2523:       if (v0) {
2524:         for(int d = 0; d < dim; d++) {
2525:           v0[d] = coords[d];
2526:         }
2527:       }
2528:       if (J) {
2529:         double x_1 = coords[2] - coords[0];
2530:         double y_1 = coords[3] - coords[1];
2531:         double x_2 = coords[6] - coords[0];
2532:         double y_2 = coords[7] - coords[1];
2533:         double x_3 = coords[4] - coords[0];
2534:         double y_3 = coords[5] - coords[1];

2536:         J[0] = x_1 + (x_3 - x_1 - x_2)*point[1];
2537:         J[1] = x_2 + (x_3 - x_1 - x_2)*point[0];
2538:         J[2] = y_1 + (y_3 - y_1 - y_2)*point[1];
2539:         J[3] = y_1 + (y_3 - y_1 - y_2)*point[0];
2540:         detJ = J[0]*J[3] - J[1]*J[2];
2541:         PetscLogFlopsNoError(6.0 + 16.0 + 3.0);
2542:       }
2543:       if (invJ) {
2544:         invDet  = 1.0/detJ;
2545:         invJ[0] =  invDet*J[3];
2546:         invJ[1] = -invDet*J[1];
2547:         invJ[2] = -invDet*J[2];
2548:         invJ[3] =  invDet*J[0];
2549:         PetscLogFlopsNoError(5.0);
2550:       }
2551:     };
2552:     void computeTetrahedronGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2553:       const double *coords = this->restrictClosure(coordinates, e);
2554:       const int     dim    = 3;
2555:       double        invDet;

2557:       if (v0) {
2558:         for(int d = 0; d < dim; d++) {
2559:           v0[d] = coords[d];
2560:         }
2561:       }
2562:       if (J) {
2563:         for(int d = 0; d < dim; d++) {
2564:           for(int f = 0; f < dim; f++) {
2565:             J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
2566:           }
2567:         }
2568:         // The minus sign is here since I orient the first face to get the outward normal
2569:         detJ = -(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
2570:                  J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
2571:                  J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
2572:         PetscLogFlopsNoError(18.0 + 12.0);
2573:       }
2574:       if (invJ) {
2575:         invDet  = -1.0/detJ;
2576:         invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
2577:         invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
2578:         invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
2579:         invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
2580:         invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
2581:         invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
2582:         invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
2583:         invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
2584:         invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
2585:         PetscLogFlopsNoError(37.0);
2586:       }
2587:     };
2588:     void computeHexahedralGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double point[], double v0[], double J[], double invJ[], double& detJ) {
2589:       const double *coords = this->restrictClosure(coordinates, e);
2590:       const int     dim    = 3;
2591:       double        invDet;

2593:       if (v0) {
2594:         for(int d = 0; d < dim; d++) {
2595:           v0[d] = coords[d];
2596:         }
2597:       }
2598:       if (J) {
2599:         double x_1 = coords[3]  - coords[0];
2600:         double y_1 = coords[4]  - coords[1];
2601:         double z_1 = coords[5]  - coords[2];
2602:         double x_2 = coords[6]  - coords[0];
2603:         double y_2 = coords[7]  - coords[1];
2604:         double z_2 = coords[8]  - coords[2];
2605:         double x_3 = coords[9]  - coords[0];
2606:         double y_3 = coords[10] - coords[1];
2607:         double z_3 = coords[11] - coords[2];
2608:         double x_4 = coords[12] - coords[0];
2609:         double y_4 = coords[13] - coords[1];
2610:         double z_4 = coords[14] - coords[2];
2611:         double x_5 = coords[15] - coords[0];
2612:         double y_5 = coords[16] - coords[1];
2613:         double z_5 = coords[17] - coords[2];
2614:         double x_6 = coords[18] - coords[0];
2615:         double y_6 = coords[19] - coords[1];
2616:         double z_6 = coords[20] - coords[2];
2617:         double x_7 = coords[21] - coords[0];
2618:         double y_7 = coords[22] - coords[1];
2619:         double z_7 = coords[23] - coords[2];
2620:         double g_x = x_1 - x_2 + x_3 + x_4 - x_5 + x_6 - x_7;
2621:         double g_y = y_1 - y_2 + y_3 + y_4 - y_5 + y_6 - y_7;
2622:         double g_z = z_1 - z_2 + z_3 + z_4 - z_5 + z_6 - z_7;

2624:         J[0] = x_1 + (x_2 - x_1 - x_3)*point[1] + (x_5 - x_1 - x_4)*point[2] + g_x*point[1]*point[2];
2625:         J[1] = x_3 + (x_2 - x_1 - x_3)*point[0] + (x_7 - x_3 - x_4)*point[2] + g_x*point[2]*point[0];
2626:         J[2] = x_4 + (x_7 - x_3 - x_4)*point[1] + (x_5 - x_1 - x_4)*point[0] + g_x*point[0]*point[1];
2627:         J[3] = y_1 + (y_2 - y_1 - y_3)*point[1] + (y_5 - y_1 - y_4)*point[2] + g_y*point[1]*point[2];
2628:         J[4] = y_3 + (y_2 - y_1 - y_3)*point[0] + (y_7 - y_3 - y_4)*point[2] + g_y*point[2]*point[0];
2629:         J[5] = y_4 + (y_7 - y_3 - y_4)*point[1] + (y_5 - y_1 - y_4)*point[0] + g_y*point[0]*point[1];
2630:         J[6] = z_1 + (z_2 - z_1 - z_3)*point[1] + (z_5 - z_1 - z_4)*point[2] + g_z*point[1]*point[2];
2631:         J[7] = z_3 + (z_2 - z_1 - z_3)*point[0] + (z_7 - z_3 - z_4)*point[2] + g_z*point[2]*point[0];
2632:         J[8] = z_4 + (z_7 - z_3 - z_4)*point[1] + (z_5 - z_1 - z_4)*point[0] + g_z*point[0]*point[1];
2633:         detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
2634:                 J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
2635:                 J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
2636:         PetscLogFlopsNoError(39.0 + 81.0 + 12.0);
2637:       }
2638:       if (invJ) {
2639:         invDet  = 1.0/detJ;
2640:         invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
2641:         invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
2642:         invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
2643:         invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
2644:         invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
2645:         invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
2646:         invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
2647:         invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
2648:         invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
2649:         PetscLogFlopsNoError(37.0);
2650:       }
2651:     };
2652:     void computeElementGeometry(const Obj<real_section_type>& coordinates, const point_type& e, double v0[], double J[], double invJ[], double& detJ) {
2653:       if (this->_dim == 2) {
2654:         computeTriangleGeometry(coordinates, e, v0, J, invJ, detJ);
2655:       } else if (this->_dim == 3) {
2656:         computeTetrahedronGeometry(coordinates, e, v0, J, invJ, detJ);
2657:       } else {
2658:         throw ALE::Exception("Unsupported dimension for element geometry computation");
2659:       }
2660:     };
2661:     void computeLineFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2662:       const arrow_section_type::point_type arrow(cell, face);
2663:       const bool reversed = (this->getArrowSection("orientation")->restrictPoint(arrow)[0] == -2);
2664:       const int  dim      = this->getDimension();
2665:       double     norm     = 0.0;
2666:       double    *vec      = tangent;

2668:       if (f == 0) {
2669:         vec[0] = 0.0;        vec[1] = -1.0;
2670:       } else if (f == 1) {
2671:         vec[0] = 0.70710678; vec[1] = 0.70710678;
2672:       } else if (f == 2) {
2673:         vec[0] = -1.0;       vec[1] = 0.0;
2674:       }
2675:       for(int d = 0; d < dim; ++d) {
2676:         normal[d] = 0.0;
2677:         for(int e = 0; e < dim; ++e) normal[d] += cellInvJ[e*dim+d]*vec[e];
2678:         if (reversed) normal[d] = -normal[d];
2679:         norm += normal[d]*normal[d];
2680:       }
2681:       norm = std::sqrt(norm);
2682:       for(int d = 0; d < dim; ++d) {
2683:         normal[d] /= norm;
2684:       }
2685:       tangent[0] =  normal[1];
2686:       tangent[1] = -normal[0];
2687:       if (this->debug()) {
2688:         std::cout << "Cell: " << cell << " Face: " << face << "("<<f<<")" << std::endl;
2689:         for(int d = 0; d < dim; ++d) {
2690:           std::cout << "Normal["<<d<<"]: " << normal[d] << " Tangent["<<d<<"]: " << tangent[d] << std::endl;
2691:         }
2692:       }
2693:       // Now get 1D Jacobian info
2694:       //   Should be a way to get this directly
2695:       const double *coords = this->restrictClosure(this->getRealSection("coordinates"), face);
2696:       detJ    = std::sqrt(PetscSqr(coords[1*2+0] - coords[0*2+0]) + PetscSqr(coords[1*2+1] - coords[0*2+1]))/2.0;
2697:       invJ[0] = 1.0/detJ;
2698:     };
2699:     void computeTriangleFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2700:       const arrow_section_type::point_type arrow(cell, face);
2701:       const bool reversed = this->getArrowSection("orientation")->restrictPoint(arrow)[0] < 0;
2702:       const int  dim      = this->getDimension();
2703:       const int  faceDim  = dim-1;
2704:       double     norm     = 0.0;
2705:       double    *vec      = tangent;

2707:       if (f == 0) {
2708:         vec[0] = 0.0;        vec[1] = 0.0;        vec[2] = -1.0;
2709:       } else if (f == 1) {
2710:         vec[0] = 0.0;        vec[1] = -1.0;       vec[2] = 0.0;
2711:       } else if (f == 2) {
2712:         vec[0] = 0.57735027; vec[1] = 0.57735027; vec[2] = 0.57735027;
2713:       } else if (f == 3) {
2714:         vec[0] = -1.0;       vec[1] = 0.0;        vec[2] = 0.0;
2715:       }
2716:       for(int d = 0; d < dim; ++d) {
2717:         normal[d] = 0.0;
2718:         for(int e = 0; e < dim; ++e) normal[d] += cellInvJ[e*dim+d]*vec[e];
2719:         if (reversed) normal[d] = -normal[d];
2720:         norm += normal[d]*normal[d];
2721:       }
2722:       norm = std::sqrt(norm);
2723:       for(int d = 0; d < dim; ++d) {
2724:         normal[d] /= norm;
2725:       }
2726:       // Get tangents
2727:       tangent[0] = normal[1] - normal[2];
2728:       tangent[1] = normal[2] - normal[0];
2729:       tangent[2] = normal[0] - normal[1];
2730:       norm = 0.0;
2731:       for(int d = 0; d < dim; ++d) {
2732:         norm += tangent[d]*tangent[d];
2733:       }
2734:       norm = std::sqrt(norm);
2735:       for(int d = 0; d < dim; ++d) {
2736:         tangent[d] /= norm;
2737:       }
2738:       tangent[3] = normal[1]*tangent[2] - normal[2]*tangent[1];
2739:       tangent[4] = normal[2]*tangent[0] - normal[0]*tangent[2];
2740:       tangent[5] = normal[0]*tangent[1] - normal[1]*tangent[0];
2741:       if (this->debug()) {
2742:         std::cout << "Cell: " << cell << " Face: " << face << "("<<f<<")" << std::endl;
2743:         for(int d = 0; d < dim; ++d) {
2744:           std::cout << "Normal["<<d<<"]: " << normal[d] << " TangentA["<<d<<"]: " << tangent[d] << " TangentB["<<d<<"]: " << tangent[dim+d] << std::endl;
2745:         }
2746:       }
2747:       // Now get 2D Jacobian info
2748:       //   Should be a way to get this directly
2749:       const double *coords = this->restrictClosure(this->getRealSection("coordinates"), face);
2750:       // Rotate so that normal in z
2751:       double invR[9], R[9];
2752:       double detR, invDetR;
2753:       for(int d = 0; d < dim; d++) {
2754:         invR[d*dim+0] = tangent[d];
2755:         invR[d*dim+1] = tangent[dim+d];
2756:         invR[d*dim+2] = normal[d];
2757:       }
2758:       invDetR = (invR[0*3+0]*(invR[1*3+1]*invR[2*3+2] - invR[1*3+2]*invR[2*3+1]) +
2759:                  invR[0*3+1]*(invR[1*3+2]*invR[2*3+0] - invR[1*3+0]*invR[2*3+2]) +
2760:                  invR[0*3+2]*(invR[1*3+0]*invR[2*3+1] - invR[1*3+1]*invR[2*3+0]));
2761:       detR  = 1.0/invDetR;
2762:       R[0*3+0] = detR*(invR[1*3+1]*invR[2*3+2] - invR[1*3+2]*invR[2*3+1]);
2763:       R[0*3+1] = detR*(invR[0*3+2]*invR[2*3+1] - invR[0*3+1]*invR[2*3+2]);
2764:       R[0*3+2] = detR*(invR[0*3+1]*invR[1*3+2] - invR[0*3+2]*invR[1*3+1]);
2765:       R[1*3+0] = detR*(invR[1*3+2]*invR[2*3+0] - invR[1*3+0]*invR[2*3+2]);
2766:       R[1*3+1] = detR*(invR[0*3+0]*invR[2*3+2] - invR[0*3+2]*invR[2*3+0]);
2767:       R[1*3+2] = detR*(invR[0*3+2]*invR[1*3+0] - invR[0*3+0]*invR[1*3+2]);
2768:       R[2*3+0] = detR*(invR[1*3+0]*invR[2*3+1] - invR[1*3+1]*invR[2*3+0]);
2769:       R[2*3+1] = detR*(invR[0*3+1]*invR[2*3+0] - invR[0*3+0]*invR[2*3+1]);
2770:       R[2*3+2] = detR*(invR[0*3+0]*invR[1*3+1] - invR[0*3+1]*invR[1*3+0]);
2771:       for(int d = 0; d < dim; d++) {
2772:         for(int e = 0; e < dim; e++) {
2773:           invR[d*dim+e] = 0.0;
2774:           for(int g = 0; g < dim; g++) {
2775:             invR[d*dim+e] += R[e*dim+g]*coords[d*dim+g];
2776:           }
2777:         }
2778:       }
2779:       for(int d = dim-1; d >= 0; --d) {
2780:         invR[d*dim+2] -= invR[0*dim+2];
2781:         if (this->debug() && (d == dim-1)) {
2782:           double ref[9];
2783:           for(int q = 0; q < dim; q++) {
2784:             for(int e = 0; e < dim; e++) {
2785:               ref[q*dim+e] = 0.0;
2786:               for(int g = 0; g < dim; g++) {
2787:                 ref[q*dim+e] += cellInvJ[e*dim+g]*coords[q*dim+g];
2788:               }
2789:             }
2790:           }
2791:           std::cout << "f: " << f << std::endl;
2792:           std::cout << this->printMatrix(std::string("coords"), dim, dim, coords) << std::endl;
2793:           std::cout << this->printMatrix(std::string("ref coords"), dim, dim, ref) << std::endl;
2794:           std::cout << this->printMatrix(std::string("R"), dim, dim, R) << std::endl;
2795:           std::cout << this->printMatrix(std::string("invR"), dim, dim, invR) << std::endl;
2796:         }
2797:         if (fabs(invR[d*dim+2]) > 1.0e-8) {
2798:           throw ALE::Exception("Invalid rotation");
2799:         }
2800:       }
2801:       double J[4];
2802:       for(int d = 0; d < faceDim; d++) {
2803:         for(int e = 0; e < faceDim; e++) {
2804:           J[d*faceDim+e] = 0.5*(invR[(e+1)*dim+d] - invR[0*dim+d]);
2805:         }
2806:       }
2807:       detJ = fabs(J[0]*J[3] - J[1]*J[2]);
2808:       // Probably need something here if detJ < 0
2809:       const double invDet = 1.0/detJ;
2810:       invJ[0] =  invDet*J[3];
2811:       invJ[1] = -invDet*J[1];
2812:       invJ[2] = -invDet*J[2];
2813:       invJ[3] =  invDet*J[0];
2814:     };
2815:     void computeFaceGeometry(const point_type& cell, const point_type& face, const int f, const double cellInvJ[], double invJ[], double& detJ, double normal[], double tangent[]) {
2816:       if (this->_dim == 2) {
2817:         computeLineFaceGeometry(cell, face, f, cellInvJ, invJ, detJ, normal, tangent);
2818:       } else if (this->_dim == 3) {
2819:         computeTriangleFaceGeometry(cell, face, f, cellInvJ, invJ, detJ, normal, tangent);
2820:       } else {
2821:         throw ALE::Exception("Unsupported dimension for element geometry computation");
2822:       }
2823:     };
2824:     double getMaxVolume() {
2825:       const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2826:       const Obj<label_sequence>&    cells       = this->heightStratum(0);
2827:       const int                     dim         = this->getDimension();
2828:       double v0[3], J[9], invJ[9], detJ, refVolume = 0.0, maxVolume = 0.0;

2830:       if (dim == 1) refVolume = 2.0;
2831:       if (dim == 2) refVolume = 2.0;
2832:       if (dim == 3) refVolume = 4.0/3.0;
2833:       for(label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2834:         this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2835:         maxVolume = std::max(maxVolume, detJ*refVolume);
2836:       }
2837:       return maxVolume;
2838:     };
2839:     // Find the cell in which this point lies (stupid algorithm)
2840:     point_type locatePoint_2D(const real_section_type::value_type point[]) {
2841:       const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2842:       const Obj<label_sequence>&    cells       = this->heightStratum(0);
2843:       const int                     embedDim    = 2;
2844:       double v0[2], J[4], invJ[4], detJ;

2846:       for(label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2847:         this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2848:         double xi   = invJ[0*embedDim+0]*(point[0] - v0[0]) + invJ[0*embedDim+1]*(point[1] - v0[1]);
2849:         double eta  = invJ[1*embedDim+0]*(point[0] - v0[0]) + invJ[1*embedDim+1]*(point[1] - v0[1]);

2851:         if ((xi >= 0.0) && (eta >= 0.0) && (xi + eta <= 2.0)) {
2852:           return *c_iter;
2853:         }
2854:       }
2855:       throw ALE::Exception("Could not locate point");
2856:     };
2857:     //   Assume a simplex and 3D
2858:     point_type locatePoint_3D(const real_section_type::value_type point[]) {
2859:       const Obj<real_section_type>& coordinates = this->getRealSection("coordinates");
2860:       const Obj<label_sequence>&    cells       = this->heightStratum(0);
2861:       const int                     embedDim    = 3;
2862:       double v0[3], J[9], invJ[9], detJ;

2864:       for(label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
2865:         this->computeElementGeometry(coordinates, *c_iter, v0, J, invJ, detJ);
2866:         double xi   = invJ[0*embedDim+0]*(point[0] - v0[0]) + invJ[0*embedDim+1]*(point[1] - v0[1]) + invJ[0*embedDim+2]*(point[2] - v0[2]);
2867:         double eta  = invJ[1*embedDim+0]*(point[0] - v0[0]) + invJ[1*embedDim+1]*(point[1] - v0[1]) + invJ[1*embedDim+2]*(point[2] - v0[2]);
2868:         double zeta = invJ[2*embedDim+0]*(point[0] - v0[0]) + invJ[2*embedDim+1]*(point[1] - v0[1]) + invJ[2*embedDim+2]*(point[2] - v0[2]);

2870:         if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) {
2871:           return *c_iter;
2872:         }
2873:       }
2874:       throw ALE::Exception("Could not locate point");
2875:     };
2876:     point_type locatePoint(const real_section_type::value_type point[], point_type guess = -1) {
2877:       //guess overrides this by saying that we already know the relation of this point to this mesh.  We will need to make it a more robust "guess" later for more than P1
2878:       if (guess != -1) {
2879:         return guess;
2880:       }else if (this->_dim == 2) {
2881:         return locatePoint_2D(point);
2882:       } else if (this->_dim == 3) {
2883:         return locatePoint_3D(point);
2884:       } else {
2885:         throw ALE::Exception("No point location for mesh dimension");
2886:       }
2887:     };
2888:   public: // Discretization
2889:     void markBoundaryCells(const std::string& name, const int marker = 1, const int newMarker = 2, const bool onlyVertices = false) {
2890:       const Obj<label_type>&     label    = this->getLabel(name);
2891:       const Obj<label_sequence>& boundary = this->getLabelStratum(name, marker);
2892:       const Obj<sieve_type>&     sieve    = this->getSieve();

2894:       if (!onlyVertices) {
2895:         const label_sequence::iterator end = boundary->end();

2897:         for(label_sequence::iterator e_iter = boundary->begin(); e_iter != end; ++e_iter) {
2898:           if (this->height(*e_iter) == 1) {
2899:             const point_type cell = *sieve->support(*e_iter)->begin();

2901:             this->setValue(label, cell, newMarker);
2902:           }
2903:         }
2904:       } else {
2905:         const label_sequence::iterator end   = boundary->end();
2906:         const int                      depth = this->depth();

2908:         for(label_sequence::iterator v_iter = boundary->begin(); v_iter != end; ++v_iter) {
2909:           const Obj<supportArray>      support = sieve->nSupport(*v_iter, depth);
2910:           const supportArray::iterator sEnd    = support->end();

2912:           for(supportArray::iterator c_iter = support->begin(); c_iter != sEnd; ++c_iter) {
2913:             const Obj<sieve_type::traits::coneSequence>&     cone = sieve->cone(*c_iter);
2914:             const sieve_type::traits::coneSequence::iterator cEnd = cone->end();

2916:             for(sieve_type::traits::coneSequence::iterator e_iter = cone->begin(); e_iter != cEnd; ++e_iter) {
2917:               if (sieve->support(*e_iter)->size() == 1) {
2918:                 this->setValue(label, *c_iter, newMarker);
2919:                 break;
2920:               }
2921:             }
2922:           }
2923:         }
2924:       }
2925:     };
2926:     int setFiberDimensions(const Obj<real_section_type>& s, const Obj<names_type>& discs, names_type& bcLabels) {
2927:       const int debug  = this->debug();
2928:       int       maxDof = 0;

2930:       for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
2931:         s->addSpace();
2932:       }
2933:       for(int d = 0; d <= this->_dim; ++d) {
2934:         int numDof = 0;
2935:         int f      = 0;

2937:         for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2938:           const Obj<ALE::Discretization>& disc = this->getDiscretization(*f_iter);
2939:           const int                       sDof = disc->getNumDof(d);

2941:           numDof += sDof;
2942:           if (sDof) s->setFiberDimension(this->depthStratum(d), sDof, f);
2943:         }
2944:         if (numDof) s->setFiberDimension(this->depthStratum(d), numDof);
2945:         maxDof = std::max(maxDof, numDof);
2946:       }
2947:       // Process exclusions
2948:       int f = 0;

2950:       for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2951:         const Obj<ALE::Discretization>& disc      = this->getDiscretization(*f_iter);
2952:         std::string                     labelName = "exclude-"+*f_iter;
2953:         std::set<point_type>            seen;

2955:         if (this->hasLabel(labelName)) {
2956:           const Obj<label_type>&         label     = this->getLabel(labelName);
2957:           const Obj<label_sequence>&     exclusion = this->getLabelStratum(labelName, 1);
2958:           const label_sequence::iterator end       = exclusion->end();
2959:           if (debug > 1) {label->view(labelName.c_str());}

2961:           for(label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
2962:             const Obj<coneArray>      closure = ALE::SieveAlg<ALE::Mesh>::closure(this, this->getArrowSection("orientation"), *e_iter);
2963:             const coneArray::iterator cEnd    = closure->end();

2965:             for(coneArray::iterator c_iter = closure->begin(); c_iter != cEnd; ++c_iter) {
2966:               if (seen.find(*c_iter) != seen.end()) continue;
2967:               if (this->getValue(label, *c_iter) == 1) {
2968:                 seen.insert(*c_iter);
2969:                 s->setFiberDimension(*c_iter, 0, f);
2970:                 s->addFiberDimension(*c_iter, -disc->getNumDof(this->depth(*c_iter)));
2971:                 if (debug > 1) {std::cout << "  cell: " << *c_iter << " dim: " << disc->getNumDof(this->depth(*c_iter)) << std::endl;}
2972:               }
2973:             }
2974:           }
2975:         }
2976:       }
2977:       // Process constraints
2978:       f = 0;
2979:       for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
2980:         const Obj<ALE::Discretization>&    disc        = this->getDiscretization(*f_iter);
2981:         const Obj<std::set<std::string> >  bcs         = disc->getBoundaryConditions();
2982:         std::string                        excludeName = "exclude-"+*f_iter;

2984:         for(std::set<std::string>::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter) {
2985:           const Obj<ALE::BoundaryCondition>& bc       = disc->getBoundaryCondition(*bc_iter);
2986:           const Obj<label_sequence>&         boundary = this->getLabelStratum(bc->getLabelName(), bc->getMarker());

2988:           bcLabels.insert(bc->getLabelName());
2989:           if (this->hasLabel(excludeName)) {
2990:             const Obj<label_type>& label = this->getLabel(excludeName);

2992:             for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
2993:               if (!this->getValue(label, *e_iter)) {
2994:                 const int numDof = disc->getNumDof(this->depth(*e_iter));

2996:                 if (numDof) s->addConstraintDimension(*e_iter, numDof);
2997:                 if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
2998:               }
2999:             }
3000:           } else {
3001:             for(label_sequence::iterator e_iter = boundary->begin(); e_iter != boundary->end(); ++e_iter) {
3002:               const int numDof = disc->getNumDof(this->depth(*e_iter));

3004:               if (numDof) s->addConstraintDimension(*e_iter, numDof);
3005:               if (numDof) s->setConstraintDimension(*e_iter, numDof, f);
3006:             }
3007:           }
3008:         }
3009:       }
3010:       return maxDof;
3011:     };
3012:     void calculateIndices() {
3013:       // Should have an iterator over the whole tree
3014:       Obj<names_type> discs = this->getDiscretizations();
3015:       Obj<Mesh>       mesh  = this;
3016:       const int       debug = this->debug();
3017:       std::map<std::string, std::pair<int, int*> > indices;

3019:       mesh.addRef();
3020:       for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3021:         const Obj<Discretization>& disc = this->getDiscretization(*d_iter);

3023:         indices[*d_iter] = std::pair<int, int*>(0, new int[disc->size(mesh)]);
3024:         disc->setIndices(indices[*d_iter].second);
3025:       }
3026:       const Obj<label_sequence>& cells   = this->heightStratum(0);
3027:       const Obj<coneArray>       closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *cells->begin());
3028:       const coneArray::iterator  end     = closure->end();
3029:       int                        offset  = 0;

3031:       if (debug > 1) {std::cout << "Closure for first element" << std::endl;}
3032:       for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3033:         const int dim = this->depth(*cl_iter);

3035:         if (debug > 1) {std::cout << "  point " << *cl_iter << " depth " << dim << std::endl;}
3036:         for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3037:           const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3038:           const int                  num  = disc->getNumDof(dim);

3040:           if (debug > 1) {std::cout << "    disc " << disc->getName() << " numDof " << num << std::endl;}
3041:           for(int o = 0; o < num; ++o) {
3042:             indices[*d_iter].second[indices[*d_iter].first++] = offset++;
3043:           }
3044:         }
3045:       }
3046:       if (debug > 1) {
3047:         for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3048:           const Obj<Discretization>& disc = this->getDiscretization(*d_iter);

3050:           std::cout << "Discretization " << disc->getName() << " indices:";
3051:           for(int i = 0; i < indices[*d_iter].first; ++i) {
3052:             std::cout << " " << indices[*d_iter].second[i];
3053:           }
3054:           std::cout << std::endl;
3055:         }
3056:       }
3057:     };
3058:     void calculateIndicesExcluded(const Obj<real_section_type>& s, const Obj<names_type>& discs) {
3059:       typedef std::map<std::string, std::pair<int, indexSet> > indices_type;
3060:       const Obj<label_type>& indexLabel = this->createLabel("cellExclusion");
3061:       const int debug  = this->debug();
3062:       int       marker = 0;
3063:       Obj<Mesh> mesh   = this;
3064:       std::map<indices_type, int> indexMap;
3065:       indices_type                indices;

3067:       mesh.addRef();
3068:       for(names_type::const_iterator d_iter = discs->begin(); d_iter != discs->end(); ++d_iter) {
3069:         const Obj<Discretization>& disc = this->getDiscretization(*d_iter);
3070:         const int                  size = disc->size(mesh);

3072:         indices[*d_iter].second.resize(size);
3073:       }
3074:       const names_type::const_iterator dBegin = discs->begin();
3075:       const names_type::const_iterator dEnd   = discs->end();
3076:       std::set<point_type> seen;
3077:       int f = 0;

3079:       for(names_type::const_iterator f_iter = dBegin; f_iter != dEnd; ++f_iter, ++f) {
3080:         std::string labelName = "exclude-"+*f_iter;

3082:         if (this->hasLabel(labelName)) {
3083:           const Obj<label_sequence>&     exclusion = this->getLabelStratum(labelName, 1);
3084:           const label_sequence::iterator end       = exclusion->end();

3086:           if (debug > 1) {std::cout << "Processing exclusion " << labelName << std::endl;}
3087:           for(label_sequence::iterator e_iter = exclusion->begin(); e_iter != end; ++e_iter) {
3088:             if (this->height(*e_iter)) continue;
3089:             const Obj<coneArray>      closure = ALE::SieveAlg<ALE::Mesh>::closure(this, this->getArrowSection("orientation"), *e_iter);
3090:             const coneArray::iterator clEnd   = closure->end();
3091:             int                       offset  = 0;

3093:             if (debug > 1) {std::cout << "  Closure for cell " << *e_iter << std::endl;}
3094:             for(coneArray::iterator cl_iter = closure->begin(); cl_iter != clEnd; ++cl_iter) {
3095:               int g = 0;

3097:               if (debug > 1) {std::cout << "    point " << *cl_iter << std::endl;}
3098:               for(names_type::const_iterator g_iter = dBegin; g_iter != dEnd; ++g_iter, ++g) {
3099:                 const int fDim = s->getFiberDimension(*cl_iter, g);

3101:                 if (debug > 1) {std::cout << "      disc " << *g_iter << " numDof " << fDim << std::endl;}
3102:                 for(int d = 0; d < fDim; ++d) {
3103:                   indices[*g_iter].second[indices[*g_iter].first++] = offset++;
3104:                 }
3105:               }
3106:             }
3107:             const std::map<indices_type, int>::iterator entry = indexMap.find(indices);

3109:             if (debug > 1) {
3110:               for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
3111:                 for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3112:                   std::cout << "Discretization (" << i_iter->second << ") " << *g_iter << " indices:";
3113:                   for(int i = 0; i < ((indices_type) i_iter->first)[*g_iter].first; ++i) {
3114:                     std::cout << " " << ((indices_type) i_iter->first)[*g_iter].second[i];
3115:                   }
3116:                   std::cout << std::endl;
3117:                 }
3118:                 std::cout << "Comparison: " << (indices == i_iter->first) << std::endl;
3119:               }
3120:               for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3121:                 std::cout << "Discretization " << *g_iter << " indices:";
3122:                 for(int i = 0; i < indices[*g_iter].first; ++i) {
3123:                   std::cout << " " << indices[*g_iter].second[i];
3124:                 }
3125:                 std::cout << std::endl;
3126:               }
3127:             }
3128:             if (entry != indexMap.end()) {
3129:               this->setValue(indexLabel, *e_iter, entry->second);
3130:               if (debug > 1) {std::cout << "  Found existing indices with marker " << entry->second << std::endl;}
3131:             } else {
3132:               indexMap[indices] = ++marker;
3133:               this->setValue(indexLabel, *e_iter, marker);
3134:               if (debug > 1) {std::cout << "  Created new indices with marker " << marker << std::endl;}
3135:             }
3136:             for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3137:               indices[*g_iter].first  = 0;
3138:               for(unsigned int i = 0; i < indices[*g_iter].second.size(); ++i) indices[*g_iter].second[i] = 0;
3139:             }
3140:           }
3141:         }
3142:       }
3143:       if (debug > 1) {indexLabel->view("cellExclusion");}
3144:       for(std::map<indices_type, int>::iterator i_iter = indexMap.begin(); i_iter != indexMap.end(); ++i_iter) {
3145:         if (debug > 1) {std::cout << "Setting indices for marker " << i_iter->second << std::endl;}
3146:         for(names_type::const_iterator g_iter = discs->begin(); g_iter != discs->end(); ++g_iter) {
3147:           const Obj<Discretization>& disc = this->getDiscretization(*g_iter);
3148:           const indexSet  indSet   = ((indices_type) i_iter->first)[*g_iter].second;
3149:           const int       size     = indSet.size();
3150:           int            *_indices = new int[size];

3152:           if (debug > 1) {std::cout << "  field " << *g_iter << std::endl;}
3153:           for(int i = 0; i < size; ++i) {
3154:             _indices[i] = indSet[i];
3155:             if (debug > 1) {std::cout << "    indices["<<i<<"] = " << _indices[i] << std::endl;}
3156:           }
3157:           disc->setIndices(_indices, i_iter->second);
3158:         }
3159:       }
3160:     };
3161:     void setupField(const Obj<real_section_type>& s, const int cellMarker = 2, const bool noUpdate = false) {
3162:       const Obj<names_type>& discs  = this->getDiscretizations();
3163:       const int              debug  = s->debug();
3164:       names_type             bcLabels;
3165:       int                    maxDof;

3167:       maxDof = this->setFiberDimensions(s, discs, bcLabels);
3168:       this->calculateIndices();
3169:       this->calculateIndicesExcluded(s, discs);
3170:       this->allocate(s);
3171:       s->defaultConstraintDof();
3172:       const Obj<label_type>& cellExclusion = this->getLabel("cellExclusion");

3174:       if (debug > 1) {std::cout << "Setting boundary values" << std::endl;}
3175:       for(names_type::const_iterator n_iter = bcLabels.begin(); n_iter != bcLabels.end(); ++n_iter) {
3176:         const Obj<label_sequence>&     boundaryCells = this->getLabelStratum(*n_iter, cellMarker);
3177:         const Obj<real_section_type>&  coordinates   = this->getRealSection("coordinates");
3178:         const Obj<names_type>&         discs         = this->getDiscretizations();
3179:         const point_type               firstCell     = *boundaryCells->begin();
3180:         const int                      numFields     = discs->size();
3181:         real_section_type::value_type *values        = new real_section_type::value_type[this->sizeWithBC(s, firstCell)];
3182:         int                           *dofs          = new int[maxDof];
3183:         int                           *v             = new int[numFields];
3184:         double                        *v0            = new double[this->getDimension()];
3185:         double                        *J             = new double[this->getDimension()*this->getDimension()];
3186:         double                         detJ;

3188:         for(label_sequence::iterator c_iter = boundaryCells->begin(); c_iter != boundaryCells->end(); ++c_iter) {
3189:           const Obj<coneArray>      closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
3190:           const coneArray::iterator end     = closure->end();

3192:           if (debug > 1) {std::cout << "  Boundary cell " << *c_iter << std::endl;}
3193:           this->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
3194:           for(int f = 0; f < numFields; ++f) v[f] = 0;
3195:           for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3196:             const int cDim = s->getConstraintDimension(*cl_iter);
3197:             int       off  = 0;
3198:             int       f    = 0;
3199:             int       i    = -1;

3201:             if (debug > 1) {std::cout << "    point " << *cl_iter << std::endl;}
3202:             if (cDim) {
3203:               if (debug > 1) {std::cout << "      constrained excMarker: " << this->getValue(cellExclusion, *c_iter) << std::endl;}
3204:               for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3205:                 const Obj<ALE::Discretization>& disc    = this->getDiscretization(*f_iter);
3206:                 const Obj<names_type>           bcs     = disc->getBoundaryConditions();
3207:                 const int                       fDim    = s->getFiberDimension(*cl_iter, f);//disc->getNumDof(this->depth(*cl_iter));
3208:                 const int                      *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));
3209:                 int                             b       = 0;

3211:                 for(names_type::const_iterator bc_iter = bcs->begin(); bc_iter != bcs->end(); ++bc_iter, ++b) {
3212:                   const Obj<ALE::BoundaryCondition>& bc    = disc->getBoundaryCondition(*bc_iter);
3213:                   const int                          value = this->getValue(this->getLabel(bc->getLabelName()), *cl_iter);

3215:                   if (b > 0) v[f] -= fDim;
3216:                   if (value == bc->getMarker()) {
3217:                     if (debug > 1) {std::cout << "      field " << *f_iter << " marker " << value << std::endl;}
3218:                     for(int d = 0; d < fDim; ++d, ++v[f]) {
3219:                       dofs[++i] = off+d;
3220:                       if (!noUpdate) values[indices[v[f]]] = (*bc->getDualIntegrator())(v0, J, v[f], bc->getFunction());
3221:                       if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3222:                     }
3223:                     // Allow only one condition per point
3224:                     ++b;
3225:                     break;
3226:                   } else {
3227:                     if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
3228:                     for(int d = 0; d < fDim; ++d, ++v[f]) {
3229:                       values[indices[v[f]]] = 0.0;
3230:                       if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3231:                     }
3232:                   }
3233:                 }
3234:                 if (b == 0) {
3235:                   if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
3236:                   for(int d = 0; d < fDim; ++d, ++v[f]) {
3237:                     values[indices[v[f]]] = 0.0;
3238:                     if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3239:                   }
3240:                 }
3241:                 off += fDim;
3242:               }
3243:               if (i != cDim-1) {throw ALE::Exception("Invalid constraint initialization");}
3244:               s->setConstraintDof(*cl_iter, dofs);
3245:             } else {
3246:               if (debug > 1) {std::cout << "      unconstrained" << std::endl;}
3247:               for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3248:                 const Obj<ALE::Discretization>& disc    = this->getDiscretization(*f_iter);
3249:                 const int                       fDim    = s->getFiberDimension(*cl_iter, f);//disc->getNumDof(this->depth(*cl_iter));
3250:                 const int                      *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));

3252:                 if (debug > 1) {std::cout << "      field " << *f_iter << std::endl;}
3253:                 for(int d = 0; d < fDim; ++d, ++v[f]) {
3254:                   values[indices[v[f]]] = 0.0;
3255:                   if (debug > 1) {std::cout << "      setting values["<<indices[v[f]]<<"] = " << values[indices[v[f]]] << std::endl;}
3256:                 }
3257:               }
3258:             }
3259:           }
3260:           if (debug > 1) {
3261:             const Obj<coneArray>      closure = sieve_alg_type::closure(this, this->getArrowSection("orientation"), *c_iter);
3262:             const coneArray::iterator end     = closure->end();

3264:             for(int f = 0; f < numFields; ++f) v[f] = 0;
3265:             for(coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
3266:               int f = 0;
3267:               for(names_type::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter, ++f) {
3268:                 const Obj<ALE::Discretization>& disc    = this->getDiscretization(*f_iter);
3269:                 const int                       fDim    = s->getFiberDimension(*cl_iter, f);
3270:                 const int                      *indices = disc->getIndices(this->getValue(cellExclusion, *c_iter));

3272:                 for(int d = 0; d < fDim; ++d, ++v[f]) {
3273:                   std::cout << "    "<<*f_iter<<"-value["<<indices[v[f]]<<"] " << values[indices[v[f]]] << std::endl;
3274:                 }
3275:               }
3276:             }
3277:           }
3278:           if (!noUpdate) {
3279:             this->updateAll(s, *c_iter, values);
3280:           }
3281:         }
3282:         delete [] dofs;
3283:         delete [] values;
3284:         delete [] v0;
3285:         delete [] J;
3286:       }
3287:       if (debug > 1) {s->view("");}
3288:     };
3289:   public:
3290:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
3291:       if (comm == MPI_COMM_NULL) {
3292:         comm = this->comm();
3293:       }
3294:       if (name == "") {
3295:         PetscPrintf(comm, "viewing a Mesh\n");
3296:       } else {
3297:         PetscPrintf(comm, "viewing Mesh '%s'\n", name.c_str());
3298:       }
3299:       this->getSieve()->view("mesh sieve", comm);
3300:       Obj<names_type> sections = this->getRealSections();

3302:       for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3303:         this->getRealSection(*name)->view(*name);
3304:       }
3305:       sections = this->getIntSections();
3306:       for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3307:         this->getIntSection(*name)->view(*name);
3308:       }
3309:       sections = this->getArrowSections();
3310:       for(names_type::iterator name = sections->begin(); name != sections->end(); ++name) {
3311:         this->getArrowSection(*name)->view(*name);
3312:       }
3313:     };
3314:     template<typename value_type>
3315:     static std::string printMatrix(const std::string& name, const int rows, const int cols, const value_type matrix[], const int rank = -1)
3316:     {
3317:       ostringstream output;
3318:       ostringstream rankStr;

3320:       if (rank >= 0) {
3321:         rankStr << "[" << rank << "]";
3322:       }
3323:       output << rankStr.str() << name << " = " << std::endl;
3324:       for(int r = 0; r < rows; r++) {
3325:         if (r == 0) {
3326:           output << rankStr.str() << " /";
3327:         } else if (r == rows-1) {
3328:           output << rankStr.str() << " \\";
3329:         } else {
3330:           output << rankStr.str() << " |";
3331:         }
3332:         for(int c = 0; c < cols; c++) {
3333:           output << " " << matrix[r*cols+c];
3334:         }
3335:         if (r == 0) {
3336:           output << " \\" << std::endl;
3337:         } else if (r == rows-1) {
3338:           output << " /" << std::endl;
3339:         } else {
3340:           output << " |" << std::endl;
3341:         }
3342:       }
3343:       return output.str();
3344:     };
3345:   };
3346:   template<typename Mesh>
3347:   class MeshBuilder {
3348:   public:
3351:     /*
3352:       Simple square boundary:

3354:      18--5-17--4--16
3355:       |     |     |
3356:       6    10     3
3357:       |     |     |
3358:      19-11-20--9--15
3359:       |     |     |
3360:       7     8     2
3361:       |     |     |
3362:      12--0-13--1--14
3363:     */
3364:     static Obj<Mesh> createSquareBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int edges[], const int debug = 0) {
3365:       Obj<Mesh> mesh        = new Mesh(comm, 1, debug);
3366:       int       numVertices = (edges[0]+1)*(edges[1]+1);
3367:       int       numEdges    = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
3368:       double   *coords      = new double[numVertices*2];
3369:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3370:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];
3371:       int                         order    = 0;

3373:       mesh->setSieve(sieve);
3374:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3375:       if (mesh->commRank() == 0) {
3376:         /* Create sieve and ordering */
3377:         for(int v = numEdges; v < numEdges+numVertices; v++) {
3378:           vertices[v-numEdges] = typename Mesh::point_type(v);
3379:         }
3380:         for(int vy = 0; vy <= edges[1]; vy++) {
3381:           for(int ex = 0; ex < edges[0]; ex++) {
3382:             typename Mesh::point_type edge(vy*edges[0] + ex);
3383:             int vertex = vy*(edges[0]+1) + ex;

3385:             sieve->addArrow(vertices[vertex+0], edge, order++);
3386:             sieve->addArrow(vertices[vertex+1], edge, order++);
3387:             if ((vy == 0) || (vy == edges[1])) {
3388:               mesh->setValue(markers, edge, 1);
3389:               mesh->setValue(markers, vertices[vertex], 1);
3390:               if (ex == edges[0]-1) {
3391:                 mesh->setValue(markers, vertices[vertex+1], 1);
3392:               }
3393:             }
3394:           }
3395:         }
3396:         for(int vx = 0; vx <= edges[0]; vx++) {
3397:           for(int ey = 0; ey < edges[1]; ey++) {
3398:             typename Mesh::point_type edge(vx*edges[1] + ey + edges[0]*(edges[1]+1));
3399:             int vertex = ey*(edges[0]+1) + vx;

3401:             sieve->addArrow(vertices[vertex],            edge, order++);
3402:             sieve->addArrow(vertices[vertex+edges[0]+1], edge, order++);
3403:             if ((vx == 0) || (vx == edges[0])) {
3404:               mesh->setValue(markers, edge, 1);
3405:               mesh->setValue(markers, vertices[vertex], 1);
3406:               if (ey == edges[1]-1) {
3407:                 mesh->setValue(markers, vertices[vertex+edges[0]+1], 1);
3408:               }
3409:             }
3410:           }
3411:         }
3412:       }
3413:       mesh->stratify();
3414:       for(int vy = 0; vy <= edges[1]; ++vy) {
3415:         for(int vx = 0; vx <= edges[0]; ++vx) {
3416:           coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
3417:           coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
3418:         }
3419:       }
3420:       delete [] vertices;
3421:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3422:       return mesh;
3423:     };
3426:     /*
3427:       Simple square boundary:

3429:      14--5-13--4--12
3430:       |           |
3431:       6           3
3432:       |           |
3433:      15           11
3434:       |           |
3435:       7           2
3436:       |           |
3437:       8--0--9--1--10
3438:     */
3439:     static Obj<Mesh> createSquareBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int edges, const int debug = 0) {
3440:       Obj<Mesh> mesh        = new Mesh(comm, 1, debug);
3441:       int       numVertices = edges*4;
3442:       int       numEdges    = edges*4;
3443:       double   *coords      = new double[numVertices*2];
3444:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3445:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];

3447:       mesh->setSieve(sieve);
3448:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3449:       if (mesh->commRank() == 0) {
3450:         /* Create sieve and ordering */
3451:         for(int v = numEdges; v < numEdges+numVertices; v++) {
3452:           vertices[v-numEdges] = typename Mesh::point_type(v);
3453:         }
3454:         for(int e = 0; e < numEdges; ++e) {
3455:           typename Mesh::point_type edge(e);
3456:           int order = 0;

3458:           sieve->addArrow(vertices[e],                 edge, order++);
3459:           sieve->addArrow(vertices[(e+1)%numVertices], edge, order++);
3460:           mesh->setValue(markers, edge, 2);
3461:           mesh->setValue(markers, vertices[e], 1);
3462:           mesh->setValue(markers, vertices[(e+1)%numVertices], 1);
3463:         }
3464:       }
3465:       mesh->stratify();
3466:       for(int v = 0; v < edges; ++v) {
3467:         coords[(v+edges*0)*2+0] = lower[0] + ((upper[0] - lower[0])/edges)*v;
3468:         coords[(v+edges*0)*2+1] = lower[1];
3469:       }
3470:       for(int v = 0; v < edges; ++v) {
3471:         coords[(v+edges*1)*2+0] = upper[0];
3472:         coords[(v+edges*1)*2+1] = lower[1] + ((upper[1] - lower[1])/edges)*v;
3473:       }
3474:       for(int v = 0; v < edges; ++v) {
3475:         coords[(v+edges*2)*2+0] = upper[0] - ((upper[0] - lower[0])/edges)*v;
3476:         coords[(v+edges*2)*2+1] = upper[1];
3477:       }
3478:       for(int v = 0; v < edges; ++v) {
3479:         coords[(v+edges*3)*2+0] = lower[0];
3480:         coords[(v+edges*3)*2+1] = upper[1] - ((upper[1] - lower[1])/edges)*v;
3481:       }
3482:       delete [] vertices;
3483:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3484:       // Build normals for cells
3485:       const Obj<typename Mesh::real_section_type>& normals = mesh->getRealSection("normals");
3486:       const Obj<typename Mesh::label_sequence>&    cells   = mesh->heightStratum(0);

3488:       //normals->setChart(typename Mesh::real_section_type::chart_type(*std::min_element(cells->begin(), cells->end()),
3489:       //                                                               *std::max_element(cells->begin(), cells->end())+1));
3490:       normals->setFiberDimension(cells, mesh->getDimension()+1);
3491:       mesh->allocate(normals);
3492:       for(int e = 0; e < edges; ++e) {
3493:         double normal[2] = {0.0, -1.0};
3494:         normals->updatePoint(e+edges*0, normal);
3495:       }
3496:       for(int e = 0; e < edges; ++e) {
3497:         double normal[2] = {1.0, 0.0};
3498:         normals->updatePoint(e+edges*1, normal);
3499:       }
3500:       for(int e = 0; e < edges; ++e) {
3501:         double normal[2] = {0.0, 1.0};
3502:         normals->updatePoint(e+edges*2, normal);
3503:       }
3504:       for(int e = 0; e < edges; ++e) {
3505:         double normal[2] = {-1.0, 0.0};
3506:         normals->updatePoint(e+edges*3, normal);
3507:       }
3508:       return mesh;
3509:     };
3512:     /*
3513:       Simple square boundary:

3515:      18--5-17--4--16
3516:       |     |     |
3517:       6    10     3
3518:       |     |     |
3519:      19-11-20--9--15
3520:       |     |     |
3521:       7     8     2
3522:       |     |     |
3523:      12--0-13--1--14
3524:     */
3525:     static Obj<Mesh> createParticleInSquareBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int edges[], const double radius, const int partEdges, const int debug = 0) {
3526:       Obj<Mesh> mesh              = new Mesh(comm, 1, debug);
3527:       const int numSquareVertices = (edges[0]+1)*(edges[1]+1);
3528:       const int numVertices       = numSquareVertices + partEdges;
3529:       const int numSquareEdges    = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
3530:       const int numEdges          = numSquareEdges + partEdges;
3531:       double   *coords            = new double[numVertices*2];
3532:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3533:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];
3534:       int                         order    = 0;

3536:       mesh->setSieve(sieve);
3537:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3538:       if (mesh->commRank() == 0) {
3539:         /* Create sieve and ordering */
3540:         for(int v = numEdges; v < numEdges+numVertices; v++) {
3541:           vertices[v-numEdges] = typename Mesh::point_type(v);
3542:         }
3543:         // Make square
3544:         for(int vy = 0; vy <= edges[1]; vy++) {
3545:           for(int ex = 0; ex < edges[0]; ex++) {
3546:             typename Mesh::point_type edge(vy*edges[0] + ex);
3547:             int vertex = vy*(edges[0]+1) + ex;

3549:             sieve->addArrow(vertices[vertex+0], edge, order++);
3550:             sieve->addArrow(vertices[vertex+1], edge, order++);
3551:             if ((vy == 0) || (vy == edges[1])) {
3552:               mesh->setValue(markers, edge, 1);
3553:               mesh->setValue(markers, vertices[vertex], 1);
3554:               if (ex == edges[0]-1) {
3555:                 mesh->setValue(markers, vertices[vertex+1], 1);
3556:               }
3557:             }
3558:           }
3559:         }
3560:         for(int vx = 0; vx <= edges[0]; vx++) {
3561:           for(int ey = 0; ey < edges[1]; ey++) {
3562:             typename Mesh::point_type edge(vx*edges[1] + ey + edges[0]*(edges[1]+1));
3563:             int vertex = ey*(edges[0]+1) + vx;

3565:             sieve->addArrow(vertices[vertex],            edge, order++);
3566:             sieve->addArrow(vertices[vertex+edges[0]+1], edge, order++);
3567:             if ((vx == 0) || (vx == edges[0])) {
3568:               mesh->setValue(markers, edge, 1);
3569:               mesh->setValue(markers, vertices[vertex], 1);
3570:               if (ey == edges[1]-1) {
3571:                 mesh->setValue(markers, vertices[vertex+edges[0]+1], 1);
3572:               }
3573:             }
3574:           }
3575:         }
3576:         // Make particle
3577:         for(int ep = 0; ep < partEdges; ++ep) {
3578:           typename Mesh::point_type edge(numSquareEdges + ep);
3579:           const int vertexA = numSquareVertices + ep;
3580:           const int vertexB = numSquareVertices + (ep+1)%partEdges;

3582:           sieve->addArrow(vertices[vertexA], edge, order++);
3583:           sieve->addArrow(vertices[vertexB], edge, order++);
3584:           mesh->setValue(markers, edge, 2);
3585:           mesh->setValue(markers, vertices[vertexA], 2);
3586:           mesh->setValue(markers, vertices[vertexB], 2);
3587:         }
3588:       }
3589:       mesh->stratify();
3590:       for(int vy = 0; vy <= edges[1]; ++vy) {
3591:         for(int vx = 0; vx <= edges[0]; ++vx) {
3592:           coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
3593:           coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
3594:         }
3595:       }
3596:       const double centroidX = 0.5*(upper[0] + lower[0]);
3597:       const double centroidY = 0.5*(upper[1] + lower[1]);
3598:       for(int vp = 0; vp < partEdges; ++vp) {
3599:         const double rad = 2.0*PETSC_PI*vp/partEdges;
3600:         coords[(numSquareVertices+vp)*2+0] = centroidX + radius*cos(rad);
3601:         coords[(numSquareVertices+vp)*2+1] = centroidY + radius*sin(rad);
3602:       }
3603:       delete [] vertices;
3604:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3605:       return mesh;
3606:     };
3609:     /*
3610:       Simple boundary with reentrant singularity:

3612:      12--5-11
3613:       |     |
3614:       |     4
3615:       |     |
3616:       6    10--3--9
3617:       |           |
3618:       |           2
3619:       |           |
3620:       7-----1-----8
3621:     */
3622:     static Obj<Mesh> createReentrantBoundary(const MPI_Comm comm, const double lower[], const double upper[], double notchpercent[], const int debug = 0) {
3623:       Obj<Mesh> mesh        = new Mesh(comm, 1, debug);
3624:       int       numVertices = 6;
3625:       int       numEdges    = numVertices;
3626:       double   *coords      = new double[numVertices*2];
3627:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3628:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];

3630:       mesh->setSieve(sieve);
3631:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3632:       if (mesh->commRank() == 0) {
3633:         /* Create sieve and ordering */
3634:         for (int b = 0; b < numVertices; b++) {
3635:           sieve->addArrow(numEdges+b, b);
3636:           sieve->addArrow(numEdges+b, (b+1)%numVertices);
3637:           mesh->setValue(markers, b, 1);
3638:           mesh->setValue(markers, b+numVertices, 1);
3639:         }
3640:         coords[0] = upper[0];
3641:         coords[1] = lower[1];

3643:         coords[2] = lower[0];
3644:         coords[3] = lower[1];
3645: 
3646:         coords[4] = lower[0];
3647:         coords[5] = notchpercent[1]*lower[1] + (1 - notchpercent[1])*upper[1];
3648: 
3649:         coords[6] = notchpercent[0]*upper[0] + (1 - notchpercent[0])*lower[0];
3650:         coords[7] = notchpercent[1]*lower[1] + (1 - notchpercent[1])*upper[1];
3651: 
3652: 
3653:         coords[8] = notchpercent[0]*upper[0] + (1 - notchpercent[0])*lower[0];
3654:         coords[9] = upper[1];

3656:         coords[10] = upper[0];
3657:         coords[11] = upper[1];
3658:         mesh->stratify();
3659:       }
3660:       delete [] vertices;
3661:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3662:       return mesh;
3663:     }

3667:     /*
3668:       Circular boundary with reentrant singularity:

3670:          ---1
3671:       --    |
3672:      -      |
3673:      |      |
3674:      |      0-----n
3675:      |            |
3676:      -           -
3677:       --       --
3678:         -------
3679:     */
3680:     static Obj<Mesh> createCircularReentrantBoundary(const MPI_Comm comm, const int segments, const double radius, const double arc_percent, const int debug = 0) {
3681:       Obj<Mesh> mesh        = new Mesh(comm, 1, debug);
3682:       int       numVertices = segments+2;
3683:       int       numEdges    = numVertices;
3684:       double   *coords      = new double[numVertices*2];
3685:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3686:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];

3688:       mesh->setSieve(sieve);
3689:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3690:       if (mesh->commRank() == 0) {
3691:         /* Create sieve and ordering */

3693:         int startvertex = 1;
3694:         if (arc_percent < 1.) {
3695:           coords[0] = 0.;
3696:           coords[1] = 0.;
3697:         } else {
3698:           numVertices = segments;
3699:           numEdges = numVertices;
3700:           startvertex = 0;
3701:         }

3703:         for (int b = 0; b < numVertices; b++) {
3704:           sieve->addArrow(numEdges+b, b);
3705:           sieve->addArrow(numEdges+b, (b+1)%numVertices);
3706:           mesh->setValue(markers, b, 1);
3707:           mesh->setValue(markers, b+numVertices, 1);
3708:         }

3710:         double anglestep = arc_percent*2.*3.14159265/((float)segments);

3712:         for (int i = startvertex; i < numVertices; i++) {
3713:           coords[2*i] = radius * sin(anglestep*(i-startvertex));
3714:           coords[2*i+1] = radius*cos(anglestep*(i-startvertex));
3715:         }
3716:         mesh->stratify();
3717:       }
3718:       delete [] vertices;
3719:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3720:       return mesh;
3721:     };
3724:     static Obj<Mesh> createAnnularBoundary(const MPI_Comm comm, const int segments, const double centers[4], const double radii[2], const int debug = 0) {
3725:       Obj<Mesh> mesh        = new Mesh(comm, 1, debug);
3726:       int       numVertices = segments*2;
3727:       int       numEdges    = numVertices;
3728:       double   *coords      = new double[numVertices*2];
3729:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3730:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];

3732:       mesh->setSieve(sieve);
3733:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3734:       if (mesh->commRank() == 0) {
3735:         for (int e = 0; e < segments; ++e) {
3736:           sieve->addArrow(numEdges+e,              e);
3737:           sieve->addArrow(numEdges+(e+1)%segments, e);
3738:           sieve->addArrow(numEdges+segments+e,              e+segments);
3739:           sieve->addArrow(numEdges+segments+(e+1)%segments, e+segments);
3740:           mesh->setValue(markers, e,          1);
3741:           mesh->setValue(markers, e+segments, 1);
3742:           mesh->setValue(markers, e+numEdges,          1);
3743:           mesh->setValue(markers, e+numEdges+segments, 1);
3744:         }
3745:         const double anglestep = 2.0*M_PI/segments;

3747:         for (int v = 0; v < segments; ++v) {
3748:           coords[v*2]              = centers[0] + radii[0]*cos(anglestep*v);
3749:           coords[v*2+1]            = centers[1] + radii[0]*sin(anglestep*v);
3750:           coords[(v+segments)*2]   = centers[2] + radii[1]*cos(anglestep*v);
3751:           coords[(v+segments)*2+1] = centers[3] + radii[1]*sin(anglestep*v);
3752:         }
3753:         mesh->addHole(&centers[2]);
3754:       }
3755:       mesh->stratify();
3756:       delete [] vertices;
3757:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3758:       return mesh;
3759:     };
3762:     /*
3763:       Simple cubic boundary:

3765:      30----31-----32
3766:       |     |     |
3767:       |  3  |  2  |
3768:       |     |     |
3769:      27----28-----29
3770:       |     |     |
3771:       |  0  |  1  |
3772:       |     |     |
3773:      24----25-----26
3774:     */
3775:     static Obj<Mesh> createCubeBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int faces[], const int debug = 0) {
3776:       Obj<Mesh> mesh        = new Mesh(comm, 2, debug);
3777:       int       numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
3778:       int       numFaces    = 6;
3779:       double   *coords      = new double[numVertices*3];
3780:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3781:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];
3782:       int                         order    = 0;

3784:       mesh->setSieve(sieve);
3785:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3786:       if (mesh->commRank() == 0) {
3787:         /* Create sieve and ordering */
3788:         for(int v = numFaces; v < numFaces+numVertices; v++) {
3789:           vertices[v-numFaces] = typename Mesh::point_type(v);
3790:           mesh->setValue(markers, vertices[v-numFaces], 1);
3791:         }
3792:         {
3793:           // Side 0 (Front)
3794:           typename Mesh::point_type face(0);
3795:           sieve->addArrow(vertices[0], face, order++);
3796:           sieve->addArrow(vertices[1], face, order++);
3797:           sieve->addArrow(vertices[2], face, order++);
3798:           sieve->addArrow(vertices[3], face, order++);
3799:           mesh->setValue(markers, face, 1);
3800:         }
3801:         {
3802:           // Side 1 (Back)
3803:           typename Mesh::point_type face(1);
3804:           sieve->addArrow(vertices[5], face, order++);
3805:           sieve->addArrow(vertices[4], face, order++);
3806:           sieve->addArrow(vertices[7], face, order++);
3807:           sieve->addArrow(vertices[6], face, order++);
3808:           mesh->setValue(markers, face, 1);
3809:         }
3810:         {
3811:           // Side 2 (Bottom)
3812:           typename Mesh::point_type face(2);
3813:           sieve->addArrow(vertices[4], face, order++);
3814:           sieve->addArrow(vertices[5], face, order++);
3815:           sieve->addArrow(vertices[1], face, order++);
3816:           sieve->addArrow(vertices[0], face, order++);
3817:           mesh->setValue(markers, face, 1);
3818:         }
3819:         {
3820:           // Side 3 (Top)
3821:           typename Mesh::point_type face(3);
3822:           sieve->addArrow(vertices[3], face, order++);
3823:           sieve->addArrow(vertices[2], face, order++);
3824:           sieve->addArrow(vertices[6], face, order++);
3825:           sieve->addArrow(vertices[7], face, order++);
3826:           mesh->setValue(markers, face, 1);
3827:         }
3828:         {
3829:           // Side 4 (Left)
3830:           typename Mesh::point_type face(4);
3831:           sieve->addArrow(vertices[4], face, order++);
3832:           sieve->addArrow(vertices[0], face, order++);
3833:           sieve->addArrow(vertices[3], face, order++);
3834:           sieve->addArrow(vertices[7], face, order++);
3835:           mesh->setValue(markers, face, 1);
3836:         }
3837:         {
3838:           // Side 5 (Right)
3839:           typename Mesh::point_type face(5);
3840:           sieve->addArrow(vertices[1], face, order++);
3841:           sieve->addArrow(vertices[5], face, order++);
3842:           sieve->addArrow(vertices[6], face, order++);
3843:           sieve->addArrow(vertices[2], face, order++);
3844:           mesh->setValue(markers, face, 1);
3845:         }
3846:       }
3847:       mesh->stratify();
3848: #if 0
3849:       for(int vz = 0; vz <= edges[2]; ++vz) {
3850:         for(int vy = 0; vy <= edges[1]; ++vy) {
3851:           for(int vx = 0; vx <= edges[0]; ++vx) {
3852:             coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
3853:             coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
3854:             coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
3855:           }
3856:         }
3857:       }
3858: #else
3859:       coords[0*3+0] = lower[0];
3860:       coords[0*3+1] = lower[1];
3861:       coords[0*3+2] = upper[2];
3862:       coords[1*3+0] = upper[0];
3863:       coords[1*3+1] = lower[1];
3864:       coords[1*3+2] = upper[2];
3865:       coords[2*3+0] = upper[0];
3866:       coords[2*3+1] = upper[1];
3867:       coords[2*3+2] = upper[2];
3868:       coords[3*3+0] = lower[0];
3869:       coords[3*3+1] = upper[1];
3870:       coords[3*3+2] = upper[2];
3871:       coords[4*3+0] = lower[0];
3872:       coords[4*3+1] = lower[1];
3873:       coords[4*3+2] = lower[2];
3874:       coords[5*3+0] = upper[0];
3875:       coords[5*3+1] = lower[1];
3876:       coords[5*3+2] = lower[2];
3877:       coords[6*3+0] = upper[0];
3878:       coords[6*3+1] = upper[1];
3879:       coords[6*3+2] = lower[2];
3880:       coords[7*3+0] = lower[0];
3881:       coords[7*3+1] = upper[1];
3882:       coords[7*3+2] = lower[2];
3883: #endif
3884:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3885:       return mesh;
3886:     };

3888:     // Creates a triangular prism boundary
3889:     static Obj<Mesh> createPrismBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int debug = 0) {
3890:       Obj<Mesh> mesh        = new Mesh(comm, 2, debug);
3891:       int       numVertices = 6;
3892:       int       numFaces    = 5;
3893:       double   *coords      = new double[numVertices*3];
3894:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
3895:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];
3896:       int                         order    = 0;

3898:       mesh->setSieve(sieve);
3899:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
3900:       if (mesh->commRank() == 0) {
3901:         /* Create sieve and ordering */
3902:         for(int v = numFaces; v < numFaces+numVertices; v++) {
3903:           vertices[v-numFaces] = typename Mesh::point_type(v);
3904:           mesh->setValue(markers, vertices[v-numFaces], 1);
3905:         }
3906:         {
3907:           // Side 0 (Top)
3908:           typename Mesh::point_type face(0);
3909:           sieve->addArrow(vertices[0], face, order++);
3910:           sieve->addArrow(vertices[1], face, order++);
3911:           sieve->addArrow(vertices[2], face, order++);
3912:           mesh->setValue(markers, face, 1);
3913:         }
3914:         {
3915:           // Side 1 (Bottom)
3916:           typename Mesh::point_type face(1);
3917:           sieve->addArrow(vertices[3], face, order++);
3918:           sieve->addArrow(vertices[5], face, order++);
3919:           sieve->addArrow(vertices[4], face, order++);
3920:           mesh->setValue(markers, face, 1);
3921:         }
3922:         {
3923:           // Side 2 (Front)
3924:           typename Mesh::point_type face(2);
3925:           sieve->addArrow(vertices[0], face, order++);
3926:           sieve->addArrow(vertices[3], face, order++);
3927:           sieve->addArrow(vertices[4], face, order++);
3928:           sieve->addArrow(vertices[1], face, order++);
3929:           mesh->setValue(markers, face, 1);
3930:         }
3931:         {
3932:           // Side 3 (Left)
3933:           typename Mesh::point_type face(3);
3934:           sieve->addArrow(vertices[1], face, order++);
3935:           sieve->addArrow(vertices[4], face, order++);
3936:           sieve->addArrow(vertices[5], face, order++);
3937:           sieve->addArrow(vertices[2], face, order++);
3938:           mesh->setValue(markers, face, 1);
3939:         }
3940:         {
3941:           // Side 4 (Right)
3942:           typename Mesh::point_type face(4);
3943:           sieve->addArrow(vertices[2], face, order++);
3944:           sieve->addArrow(vertices[5], face, order++);
3945:           sieve->addArrow(vertices[3], face, order++);
3946:           sieve->addArrow(vertices[0], face, order++);
3947:           mesh->setValue(markers, face, 1);
3948:         }
3949:       }
3950:       mesh->stratify();
3951:       coords[0*3+0] = lower[0];
3952:       coords[0*3+1] = lower[1];
3953:       coords[0*3+2] = upper[2];
3954:       coords[1*3+0] = upper[0];
3955:       coords[1*3+1] = lower[1];
3956:       coords[1*3+2] = upper[2];
3957:       coords[2*3+0] = upper[0];
3958:       coords[2*3+1] = upper[1];
3959:       coords[2*3+2] = upper[2];
3960:       coords[3*3+0] = lower[0];
3961:       coords[3*3+1] = upper[1];
3962:       coords[3*3+2] = upper[2];
3963:       coords[4*3+0] = lower[0];
3964:       coords[4*3+1] = lower[1];
3965:       coords[4*3+2] = lower[2];
3966:       coords[5*3+0] = upper[0];
3967:       coords[5*3+1] = lower[1];
3968:       coords[5*3+2] = lower[2];
3969:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
3970:       return mesh;
3971:     };

3975:     /*    v0
3976:          / \
3977:         /   \
3978:     2  /  4  \  1
3979:       /       \
3980:      /   v12   \
3981:   v6|\   /|\   /|v5
3982:     | \v8 | v7/ |          z  
3983:     |  |7 |8 |  |          | 
3984:     |  |v13\ |  |  <-v4   / \
3985:     | v9/ 9 \v10|        x   y
3986:  v1 | 5 \   / 6 |v2
3987:      \   \ /   /
3988:       \  v11  /
3989:        \  |  /
3990:      3  \ | /
3991:          \|/
3992:           v3
3993:     */
3994:     static Obj<Mesh> createFicheraCornerBoundary(const MPI_Comm comm, const double lower[], const double upper[], const double offset[], const int debug = 0) {
3995:       Obj<Mesh> mesh            = new Mesh(comm, 2, debug);
3996:       const int nVertices = 14;
3997:       const int nFaces = 12;
3998:       double ilower[3];
3999:       ilower[0] = lower[0]*(1. - offset[0]) + upper[0]*offset[0];
4000:       ilower[1] = lower[1]*(1. - offset[1]) + upper[1]*offset[1];
4001:       ilower[2] = lower[2]*(1. - offset[2]) + upper[2]*offset[2];
4002:       double coords[nVertices*3];
4003:       //outer square-triplet
4004:       coords[0*3+0] = lower[0];
4005:       coords[0*3+1] = lower[1];
4006:       coords[0*3+2] = upper[2];
4007:       coords[1*3+0] = upper[0];
4008:       coords[1*3+1] = lower[1];
4009:       coords[1*3+2] = lower[2];
4010:       coords[2*3+0] = lower[0];
4011:       coords[2*3+1] = upper[1];
4012:       coords[2*3+2] = lower[2];
4013:       coords[3*3+0] = upper[0];
4014:       coords[3*3+1] = upper[1];
4015:       coords[3*3+2] = lower[2];
4016:       coords[4*3+0] = lower[0];
4017:       coords[4*3+1] = lower[1];
4018:       coords[4*3+2] = lower[2];
4019:       coords[5*3+0] = lower[0];
4020:       coords[5*3+1] = upper[1];
4021:       coords[5*3+2] = upper[2];
4022:       coords[6*3+0] = upper[0];
4023:       coords[6*3+1] = lower[1];
4024:       coords[6*3+2] = upper[2];

4026:       //inner square-triplet
4027:       coords[7*3+0] = ilower[0];
4028:       coords[7*3+1] = upper[1];
4029:       coords[7*3+2] = upper[2];
4030:       coords[8*3+0] = upper[0];
4031:       coords[8*3+1] = ilower[1];
4032:       coords[8*3+2] = upper[2];
4033:       coords[9*3+0] = upper[0];
4034:       coords[9*3+1] = ilower[1];
4035:       coords[9*3+2] = ilower[2];
4036:       coords[10*3+0] = ilower[0];
4037:       coords[10*3+1] = upper[1];
4038:       coords[10*3+2] = ilower[2];
4039:       coords[11*3+0] = upper[0];
4040:       coords[11*3+1] = upper[1];
4041:       coords[11*3+2] = ilower[2];
4042:       coords[12*3+0] = ilower[0];
4043:       coords[12*3+1] = ilower[1];
4044:       coords[12*3+2] = upper[2];
4045:       coords[13*3+0] = ilower[0];
4046:       coords[13*3+1] = ilower[1];
4047:       coords[13*3+2] = ilower[2];

4049: 
4050:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
4051:       mesh->setSieve(sieve);
4052:       typename Mesh::point_type p[nVertices];
4053:       typename Mesh::point_type f[nFaces];
4054:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
4055:       for (int i = 0; i < nVertices; i++) {
4056:         p[i] = typename Mesh::point_type(i+nFaces);
4057:         mesh->setValue(markers, p[i], 1);
4058:       }
4059:       for (int i = 0; i < nFaces; i++) {
4060:         f[i] = typename Mesh::point_type(i);
4061:       }
4062:       int order = 0;
4063:      //assemble the larger square sides
4064:       sieve->addArrow(p[0], f[0], order++);
4065:       sieve->addArrow(p[5], f[0], order++);
4066:       sieve->addArrow(p[2], f[0], order++);
4067:       sieve->addArrow(p[4], f[0], order++);
4068:       mesh->setValue(markers, f[0], 1);

4070:       sieve->addArrow(p[0], f[1], order++);
4071:       sieve->addArrow(p[4], f[1], order++);
4072:       sieve->addArrow(p[1], f[1], order++);
4073:       sieve->addArrow(p[6], f[1], order++);
4074:       mesh->setValue(markers, f[1], 1);

4076:       sieve->addArrow(p[4], f[2], order++);
4077:       sieve->addArrow(p[1], f[2], order++);
4078:       sieve->addArrow(p[3], f[2], order++);
4079:       sieve->addArrow(p[2], f[2], order++);
4080:       mesh->setValue(markers, f[2], 1);
4081: 
4082:       //assemble the L-shaped sides

4084:       sieve->addArrow(p[0], f[3], order++);
4085:       sieve->addArrow(p[12], f[3], order++);
4086:       sieve->addArrow(p[7], f[3], order++);
4087:       sieve->addArrow(p[5], f[3], order++);
4088:       mesh->setValue(markers, f[3], 1);

4090:       sieve->addArrow(p[0], f[4], order++);
4091:       sieve->addArrow(p[12],f[4], order++);
4092:       sieve->addArrow(p[8], f[4], order++);
4093:       sieve->addArrow(p[6], f[4], order++);
4094:       mesh->setValue(markers, f[4], 1);

4096:       sieve->addArrow(p[9], f[5], order++);
4097:       sieve->addArrow(p[1], f[5], order++);
4098:       sieve->addArrow(p[3], f[5], order++);
4099:       sieve->addArrow(p[11], f[5], order++);
4100:       mesh->setValue(markers, f[5], 1);

4102:       sieve->addArrow(p[9], f[6], order++);
4103:       sieve->addArrow(p[1], f[6], order++);
4104:       sieve->addArrow(p[6], f[6], order++);
4105:       sieve->addArrow(p[8], f[6], order++);
4106:       mesh->setValue(markers, f[6], 1);

4108:       sieve->addArrow(p[10], f[7], order++);
4109:       sieve->addArrow(p[2], f[7], order++);
4110:       sieve->addArrow(p[5], f[7], order++);
4111:       sieve->addArrow(p[7], f[7], order++);
4112:       mesh->setValue(markers, f[7], 1);

4114:       sieve->addArrow(p[10], f[8], order++);
4115:       sieve->addArrow(p[2], f[8], order++);
4116:       sieve->addArrow(p[3], f[8], order++);
4117:       sieve->addArrow(p[11], f[8], order++);
4118:       mesh->setValue(markers, f[8], 1);

4120:       //assemble the smaller square sides

4122:       sieve->addArrow(p[13], f[9], order++);
4123:       sieve->addArrow(p[10], f[9], order++);
4124:       sieve->addArrow(p[11], f[9], order++);
4125:       sieve->addArrow(p[9], f[9], order++);
4126:       mesh->setValue(markers, f[9], 1);

4128:       sieve->addArrow(p[12], f[10], order++);
4129:       sieve->addArrow(p[7], f[10], order++);
4130:       sieve->addArrow(p[10], f[10], order++);
4131:       sieve->addArrow(p[13], f[10], order++);
4132:       mesh->setValue(markers, f[10], 1);

4134:       sieve->addArrow(p[8], f[11], order++);
4135:       sieve->addArrow(p[12], f[11], order++);
4136:       sieve->addArrow(p[13], f[11], order++);
4137:       sieve->addArrow(p[9], f[11], order++);
4138:       mesh->setValue(markers, f[11], 1);

4140:       mesh->stratify();
4141:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
4142:       Obj<typename Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
4143:       //coordinates->view("coordinates");
4144:       return mesh;

4146:     }

4150:     /*
4151:       //"sphere" out a cube 

4153:     */
4154: #if 0
4155:     static Obj<Mesh> createSphereBoundary(const MPI_Comm comm, const double radius, const int refinement, const int debug = 0) {
4156:       Obj<Mesh> m = new Mesh(comm, 2, debug);
4157:       Obj<Mesh::sieve_type> s = new Mesh::sieve_type(comm, debug);
4158:       m->setSieve(s);
4159:       Mesh::point_type p = 0;
4160:       int nVertices = 8+12*(refinement)+6*(refinement)*(refinement);
4161:       Mesh::point_type vertices[nVertices];
4162:       double coords[3*nVertices];
4163:       int nCells = 6*2*(refinement+1)*(refinement+1);
4164:       double delta = 2./((double)(refinement+1));
4165:       Mesh::point_type cells[nCells];
4166:       for (int i = 0; i < nCells; i++) {
4167:         cells[i] = p;
4168:         p++;
4169:       }
4170:       for (int i = 0; i < nVertices; i++) {
4171:         vertices[i] = p;
4172:         p++;
4173:       }
4174:       //set up the corners;
4175:       //lll
4176:       coords[0*3+0] = -1.;
4177:       coords[0*3+1] = -1.;
4178:       coords[0*3+2] = -1.;
4179:       //llh
4180:       coords[1*3+0] = -1.;
4181:       coords[1*3+1] = -1.;
4182:       coords[1*3+2] = 1.;
4183:       //lhh
4184:       coords[2*3+0] = -1.;
4185:       coords[2*3+1] = 1.;
4186:       coords[2*3+2] = 1.;
4187:       //lhl
4188:       coords[3*3+0] = -1.;
4189:       coords[3*3+1] = 1.;
4190:       coords[3*3+2] = -1.;
4191:       //hhl
4192:       coords[4*3+0] = 1.;
4193:       coords[4*3+1] = 1.;
4194:       coords[4*3+2] = -1.;
4195:       //hhh
4196:       coords[5*3+0] = 1.;
4197:       coords[5*3+1] = 1.;
4198:       coords[5*3+2] = 1.;
4199:       //hlh
4200:       coords[6*3+0] = 1.;
4201:       coords[6*3+1] = -1.;
4202:       coords[6*3+2] = 1.;
4203:       //hll
4204:       coords[7*3+0] = 1.;
4205:       coords[7*3+1] = -1.;
4206:       coords[7*3+2] = -1.;
4207:       //set up the edges (always go low to high)
4208:       //xll
4209:       for (int i = 0; i < refinement; i++) {
4210:         coords[3*(8+0*refinement+i)+0] = -1. + delta*i;
4211:         coords[3*(8+0*refinement+i)+1] = -1.;
4212:         coords[3*(8+0*refinement+i)+2] = -1.;
4213:       }
4214:       //xlh
4215:       for (int i = 0; i < refinement; i++) {
4216:         coords[3*(8+1*refinement+i)+0] = -1. + delta*i;
4217:         coords[3*(8+1*refinement+i)+1] = -1.;
4218:         coords[3*(8+1*refinement+i)+2] = 1.;
4219:       }
4220:       //xhh
4221:       for (int i = 0; i < refinement; i++) {
4222:         coords[3*(8+2*refinement+i)+0] = -1. + delta*i;
4223:         coords[3*(8+2*refinement+i)+1] = 1.;
4224:         coords[3*(8+2*refinement+i)+2] = 1.;
4225:       }
4226:       //xhl
4227:       for (int i = 0; i < refinement; i++) {
4228:         coords[3*(8+3*refinement+i)+0] = -1. + delta*i;
4229:         coords[3*(8+3*refinement+i)+1] = 1.;
4230:         coords[3*(8+3*refinement+i)+2] = -1.;
4231:       }
4232:       //lxl
4233:       for (int i = 0; i < refinement; i++) {
4234:         coords[3*(8+4*refinement+i)+0] = -1.;
4235:         coords[3*(8+4*refinement+i)+1] = -1. + delta*i;
4236:         coords[3*(8+4*refinement+i)+2] = -1.;
4237:       }
4238:       //lxh
4239:       for (int i = 0; i < refinement; i++) {
4240:         coords[3*(8+5*refinement+i)+0] = -1.;
4241:         coords[3*(8+5*refinement+i)+1] = -1. + delta*i;
4242:         coords[3*(8+5*refinement+i)+2] = 1.;
4243:       }
4244:       //hxh
4245:       for (int i = 0; i < refinement; i++) {
4246:         coords[3*(8+6*refinement+i)+0] = 1.;
4247:         coords[3*(8+6*refinement+i)+1] = -1. + delta*i;
4248:         coords[3*(8+6*refinement+i)+2] = 1.;
4249:       }
4250:       //hxl
4251:       for (int i = 0; i < refinement; i++) {
4252:         coords[3*(8+7*refinement+i)+0] = 1.;
4253:         coords[3*(8+7*refinement+i)+1] = -1. + delta*i;
4254:         coords[3*(8+7*refinement+i)+2] = -1.;
4255:       }
4256:       //llx
4257:       for (int i = 0; i < refinement; i++) {
4258:         coords[3*(8+8*refinement+i)+0] = -1.;
4259:         coords[3*(8+8*refinement+i)+1] = -1.;
4260:         coords[3*(8+8*refinement+i)+2] = -1. + delta*i;
4261:       }
4262:       //lhx
4263:       for (int i = 0; i < refinement; i++) {
4264:         coords[3*(8+9*refinement+i)+0] = -1.;
4265:         coords[3*(8+9*refinement+i)+1] = 1.;
4266:         coords[3*(8+9*refinement+i)+2] = -1. + delta*i;
4267:       }
4268:       //hhx
4269:       for (int i = 0; i < refinement; i++) {
4270:         coords[3*(8+10*refinement+i)+0] = 1.;
4271:         coords[3*(8+10*refinement+i)+1] = 1.;
4272:         coords[3*(8+10*refinement+i)+2] = -1. + delta*i;
4273:       }
4274:       //hlx
4275:       for (int i = 0; i < refinement; i++) {
4276:         coords[3*(8+11*refinement+i)+0] = 1.;
4277:         coords[3*(8+11*refinement+i)+1] = -1.;
4278:         coords[3*(8+11*refinement+i)+2] = -1. + delta*i;
4279:       }
4280:       //set up the faces
4281:       //lxx
4282:       for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4283:         coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+0] = -1.;
4284:         coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4285:         coords[3*(8+12*refinement+0*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4286:       }
4287:       //hxx
4288:       for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4289:         coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+0] = 1.;
4290:         coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4291:         coords[3*(8+12*refinement+1*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4292:       }
4293:       //xlx
4294:       for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4295:         coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+0] = -1. + delta*j;
4296:         coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+1] = -1.;
4297:         coords[3*(8+12*refinement+2*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4298:       }
4299:       //xhx
4300:       for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4301:         coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+0] = -1. + delta*j;
4302:         coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+1] = 1.;
4303:         coords[3*(8+12*refinement+3*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4304:       }
4305:       //xxl
4306:       for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4307:         coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+0] = -1.;
4308:         coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4309:         coords[3*(8+12*refinement+4*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4310:       }
4311:       //xxh
4312:       for (int i = 0; i < refinement; i++) for (int j = 0; j < refinement; j++) {
4313:         coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+0] = 1.;
4314:         coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+1] = -1. + delta*j;
4315:         coords[3*(8+12*refinement+5*refinement*refinement+i*refinement+j)+2] = -1. + delta*i;
4316:       }
4317:       //stitch the corners up with the edges and the faces
4318: 
4319:       //stitch the edges to the faces
4320:       //fill in the faces
4321:       int face_offset = 8 + 12*refinement;
4322:       for (int i = 0; i < 6; i++) for (int j = 0; j < refinement; j++) for (int k = 0; k < refinement; k++) {
4323:         //build each square doublet
4324:       }
4325:     }

4327: #endif

4331:     /*
4332:       Simple cubic boundary:

4334:      30----31-----32
4335:       |     |     |
4336:       |  3  |  2  |
4337:       |     |     |
4338:      27----28-----29
4339:       |     |     |
4340:       |  0  |  1  |
4341:       |     |     |
4342:      24----25-----26
4343:     */
4344:     static Obj<Mesh> createParticleInCubeBoundary(const MPI_Comm comm, const double lower[], const double upper[], const int faces[], const double radius, const int thetaEdges, const int phiSlices, const int debug = 0) {
4345:       Obj<Mesh> mesh            = new Mesh(comm, 2, debug);
4346:       const int numCubeVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
4347:       const int numPartVertices = (thetaEdges - 1)*phiSlices + 2;
4348:       const int numVertices     = numCubeVertices + numPartVertices;
4349:       const int numCubeFaces    = 6;
4350:       const int numFaces        = numCubeFaces + thetaEdges*phiSlices;
4351:       double   *coords          = new double[numVertices*3];
4352:       const Obj<typename Mesh::sieve_type> sieve    = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
4353:       typename Mesh::point_type           *vertices = new typename Mesh::point_type[numVertices];
4354:       int                         order    = 0;

4356:       mesh->setSieve(sieve);
4357:       const Obj<typename Mesh::label_type>& markers = mesh->createLabel("marker");
4358:       if (mesh->commRank() == 0) {
4359:         // Make cube
4360:         for(int v = numFaces; v < numFaces+numVertices; v++) {
4361:           vertices[v-numFaces] = typename Mesh::point_type(v);
4362:           mesh->setValue(markers, vertices[v-numFaces], 1);
4363:         }
4364:         {
4365:           // Side 0 (Front)
4366:           typename Mesh::point_type face(0);
4367:           sieve->addArrow(vertices[0], face, order++);
4368:           sieve->addArrow(vertices[1], face, order++);
4369:           sieve->addArrow(vertices[2], face, order++);
4370:           sieve->addArrow(vertices[3], face, order++);
4371:           mesh->setValue(markers, face, 1);
4372:         }
4373:         {
4374:           // Side 1 (Back)
4375:           typename Mesh::point_type face(1);
4376:           sieve->addArrow(vertices[5], face, order++);
4377:           sieve->addArrow(vertices[4], face, order++);
4378:           sieve->addArrow(vertices[7], face, order++);
4379:           sieve->addArrow(vertices[6], face, order++);
4380:           mesh->setValue(markers, face, 1);
4381:         }
4382:         {
4383:           // Side 2 (Bottom)
4384:           typename Mesh::point_type face(2);
4385:           sieve->addArrow(vertices[4], face, order++);
4386:           sieve->addArrow(vertices[5], face, order++);
4387:           sieve->addArrow(vertices[1], face, order++);
4388:           sieve->addArrow(vertices[0], face, order++);
4389:           mesh->setValue(markers, face, 1);
4390:         }
4391:         {
4392:           // Side 3 (Top)
4393:           typename Mesh::point_type face(3);
4394:           sieve->addArrow(vertices[3], face, order++);
4395:           sieve->addArrow(vertices[2], face, order++);
4396:           sieve->addArrow(vertices[6], face, order++);
4397:           sieve->addArrow(vertices[7], face, order++);
4398:           mesh->setValue(markers, face, 1);
4399:         }
4400:         {
4401:           // Side 4 (Left)
4402:           typename Mesh::point_type face(4);
4403:           sieve->addArrow(vertices[4], face, order++);
4404:           sieve->addArrow(vertices[0], face, order++);
4405:           sieve->addArrow(vertices[3], face, order++);
4406:           sieve->addArrow(vertices[7], face, order++);
4407:           mesh->setValue(markers, face, 1);
4408:         }
4409:         {
4410:           // Side 5 (Right)
4411:           typename Mesh::point_type face(5);
4412:           sieve->addArrow(vertices[1], face, order++);
4413:           sieve->addArrow(vertices[5], face, order++);
4414:           sieve->addArrow(vertices[6], face, order++);
4415:           sieve->addArrow(vertices[2], face, order++);
4416:           mesh->setValue(markers, face, 1);
4417:         }
4418:         // Make particle
4419:         for(int s = 0; s < phiSlices; ++s) {
4420:           for(int ep = 0; ep < thetaEdges; ++ep) {
4421:             // Vertices on each slice are 0..thetaEdges
4422:             typename Mesh::point_type face(numCubeFaces + s*thetaEdges + ep);
4423:             int vertexA = numCubeVertices + ep + 0 +     s*(thetaEdges+1);
4424:             int vertexB = numCubeVertices + ep + 1 +     s*(thetaEdges+1);
4425:             int vertexC = numCubeVertices + (ep + 1 + (s+1)*(thetaEdges+1))%((thetaEdges+1)*phiSlices);
4426:             int vertexD = numCubeVertices + (ep + 0 + (s+1)*(thetaEdges+1))%((thetaEdges+1)*phiSlices);
4427:             const int correction1 = (s > 0)*((s-1)*2 + 1);
4428:             const int correction2 = (s < phiSlices-1)*(s*2 + 1);

4430:             if ((vertexA - numCubeVertices)%(thetaEdges+1) == 0) {
4431:               vertexA = vertexD = numCubeVertices;
4432:               vertexB -= correction1;
4433:               vertexC -= correction2;
4434:             } else if ((vertexB - numCubeVertices)%(thetaEdges+1) == thetaEdges) {
4435:               vertexA -= correction1;
4436:               vertexD -= correction2;
4437:               vertexB = vertexC = numCubeVertices + thetaEdges;
4438:             } else {
4439:               vertexA -= correction1;
4440:               vertexB -= correction1;
4441:               vertexC -= correction2;
4442:               vertexD -= correction2;
4443:             }
4444:             if ((vertexA >= numVertices) || (vertexB >= numVertices) || (vertexC >= numVertices) || (vertexD >= numVertices)) {
4445:               throw ALE::Exception("Bad vertex");
4446:             }
4447:             sieve->addArrow(vertices[vertexA], face, order++);
4448:             sieve->addArrow(vertices[vertexB], face, order++);
4449:             if (vertexB != vertexC) sieve->addArrow(vertices[vertexC], face, order++);
4450:             if (vertexA != vertexD) sieve->addArrow(vertices[vertexD], face, order++);
4451:             mesh->setValue(markers, face, 2);
4452:             mesh->setValue(markers, vertices[vertexA], 2);
4453:             mesh->setValue(markers, vertices[vertexB], 2);
4454:             if (vertexB != vertexC) mesh->setValue(markers, vertices[vertexC], 2);
4455:             if (vertexA != vertexD) mesh->setValue(markers, vertices[vertexD], 2);
4456:           }
4457:         }
4458:       }
4459:       mesh->stratify();
4460: #if 0
4461:       for(int vz = 0; vz <= edges[2]; ++vz) {
4462:         for(int vy = 0; vy <= edges[1]; ++vy) {
4463:           for(int vx = 0; vx <= edges[0]; ++vx) {
4464:             coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
4465:             coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
4466:             coords[((vz*(edges[1]+1)+vy)*(edges[0]+1)+vx)*2+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
4467:           }
4468:         }
4469:       }
4470: #else
4471:       coords[0*3+0] = lower[0];
4472:       coords[0*3+1] = lower[1];
4473:       coords[0*3+2] = upper[2];
4474:       coords[1*3+0] = upper[0];
4475:       coords[1*3+1] = lower[1];
4476:       coords[1*3+2] = upper[2];
4477:       coords[2*3+0] = upper[0];
4478:       coords[2*3+1] = upper[1];
4479:       coords[2*3+2] = upper[2];
4480:       coords[3*3+0] = lower[0];
4481:       coords[3*3+1] = upper[1];
4482:       coords[3*3+2] = upper[2];
4483:       coords[4*3+0] = lower[0];
4484:       coords[4*3+1] = lower[1];
4485:       coords[4*3+2] = lower[2];
4486:       coords[5*3+0] = upper[0];
4487:       coords[5*3+1] = lower[1];
4488:       coords[5*3+2] = lower[2];
4489:       coords[6*3+0] = upper[0];
4490:       coords[6*3+1] = upper[1];
4491:       coords[6*3+2] = lower[2];
4492:       coords[7*3+0] = lower[0];
4493:       coords[7*3+1] = upper[1];
4494:       coords[7*3+2] = lower[2];
4495: #endif
4496:       const double centroidX = 0.5*(upper[0] + lower[0]);
4497:       const double centroidY = 0.5*(upper[1] + lower[1]);
4498:       const double centroidZ = 0.5*(upper[2] + lower[2]);
4499:       for(int s = 0; s < phiSlices; ++s) {
4500:         for(int v = 0; v <= thetaEdges; ++v) {
4501:           int          vertex  = numCubeVertices + v + s*(thetaEdges+1);
4502:           const double theta   = v*(PETSC_PI/thetaEdges);
4503:           const double phi     = s*(2.0*PETSC_PI/phiSlices);
4504:           const int correction = (s > 0)*((s-1)*2 + 1);

4506:           if ((vertex- numCubeVertices)%(thetaEdges+1) == 0) {
4507:             vertex = numCubeVertices;
4508:           } else if ((vertex - numCubeVertices)%(thetaEdges+1) == thetaEdges) {
4509:             vertex = numCubeVertices + thetaEdges;
4510:           } else {
4511:             vertex -= correction;
4512:           }
4513:           coords[vertex*3+0] = centroidX + radius*sin(theta)*cos(phi);
4514:           coords[vertex*3+1] = centroidY + radius*sin(theta)*sin(phi);
4515:           coords[vertex*3+2] = centroidZ + radius*cos(theta);
4516:         }
4517:       }
4518:       delete [] vertices;
4519:       ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, mesh->getDimension()+1, coords);
4520:       return mesh;
4521:     };
4522:     // This method takes a tetrahedral mesh and performs a 1 --> 8 refinement of each cell
4523:     //   It does this by adding a new vertex at the midpoint of each edge
4524:     template<typename MeshType, typename MapType>
4525:     static void refineTetrahedra(MeshType& mesh, MeshType& newMesh, MapType& edge2vertex) {
4526:       typedef typename MeshType::sieve_type sieve_type;
4527:       typedef typename MeshType::point_type point_type;
4528:       typedef typename MapType::key_type    edge_type;

4530:       const int numCells       = mesh.heightStratum(0)->size();
4531:       const int numVertices    = mesh.depthStratum(0)->size();
4532:       // Calculate number of new cells
4533:       const int numNewCells    = numCells * 8;
4534:       // Bound number of new vertices
4535:       //const int maxNewVertices = numCells * 6;
4536:       int       curNewVertex   = numNewCells + numVertices;

4538:       // Loop over cells
4539:       const Obj<sieve_type>&                         sieve    = mesh.getSieve();
4540:       const Obj<sieve_type>&                         newSieve = newMesh.getSieve();
4541:       ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, sieve->getMaxConeSize()));

4543:       // First compute map from edges to new vertices
4544:       for(int c = 0; c < numCells; ++c) {
4545:         sieve->cone(c, cV);
4546:         assert(cV.getSize() == 4);
4547:         const point_type *cone = cV.getPoints();

4549:         //   As per Brad's diagram
4550:         edge_type edges[6] = {edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1])),
4551:                               edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2])),
4552:                               edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0])),
4553:                               edge_type(std::min(cone[0], cone[3]), std::max(cone[0], cone[3])),
4554:                               edge_type(std::min(cone[1], cone[3]), std::max(cone[1], cone[3])),
4555:                               edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]))};
4556:         //   Check that vertex does not yet exist
4557:         for(int v = 0; v < 6; ++v) {
4558:           if (edge2vertex.find(edges[v]) == edge2vertex.end()) {
4559:             edge2vertex[edges[v]] = curNewVertex++;
4560:           }
4561:         }
4562:         cV.clear();
4563:       }
4564:       // Reallocate the sieve chart
4565:       newSieve->setChart(typename sieve_type::chart_type(0, curNewVertex));
4566:       // Create new sieve with correct sizes for refined cells
4567:       for(int c = 0; c < numCells; ++c) {
4568:         sieve->cone(c, cV);
4569:         assert(cV.getSize() == 4);
4570:         const point_type *cone = cV.getPoints();

4572:         // As per Brad's diagram
4573:         edge_type edges[6] = {edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1])),
4574:                               edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2])),
4575:                               edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0])),
4576:                               edge_type(std::min(cone[0], cone[3]), std::max(cone[0], cone[3])),
4577:                               edge_type(std::min(cone[1], cone[3]), std::max(cone[1], cone[3])),
4578:                               edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]))};
4579:         //   Check that vertex does not yet exist
4580:         point_type newVertices[6];

4582:         for(int v = 0; v < 6; ++v) {
4583:           newVertices[v] = edge2vertex[edges[v]];
4584:         }
4585:         // Set new sizes
4586:         for(int nc = 0; nc < 8; ++nc) {newSieve->setConeSize(c*8+nc, 4);}
4587:         const point_type offset = numNewCells - numCells;

4589:         point_type cell0[4] = {cone[0]+offset, newVertices[3], newVertices[0], newVertices[2]};
4590:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell0[v], 1);}
4591:         point_type cell1[4] = {cone[1]+offset, newVertices[4], newVertices[1], newVertices[0]};
4592:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell1[v], 1);}
4593:         point_type cell2[4] = {cone[2]+offset, newVertices[5], newVertices[2], newVertices[1]};
4594:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell2[v], 1);}
4595:         point_type cell3[4] = {cone[3]+offset, newVertices[3], newVertices[5], newVertices[4]};
4596:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell3[v], 1);}
4597:         point_type cell4[4] = {newVertices[0], newVertices[3], newVertices[4], newVertices[2]};
4598:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell4[v], 1);}
4599:         point_type cell5[4] = {newVertices[1], newVertices[4], newVertices[5], newVertices[3]};
4600:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell5[v], 1);}
4601:         point_type cell6[4] = {newVertices[2], newVertices[5], newVertices[3], newVertices[1]};
4602:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell6[v], 1);}
4603:         point_type cell7[4] = {newVertices[0], newVertices[1], newVertices[2], newVertices[3]};
4604:         for(int v = 0; v < 4; ++v) {newSieve->addSupportSize(cell7[v], 1);}
4605:         cV.clear();
4606:       }
4607:       newSieve->allocate();
4608:       const int   numNewVertices = curNewVertex - numNewCells;
4609:       point_type *vertex2edge    = new point_type[numNewVertices*2];

4611:       // Create refined cells in new sieve
4612:       for(int c = 0; c < numCells; ++c) {
4613:         sieve->cone(c, cV);
4614:         assert(cV.getSize() == 4);
4615:         const point_type *cone = cV.getPoints();

4617:         // As per Brad's diagram
4618:         edge_type edges[6] = {edge_type(std::min(cone[0], cone[1]), std::max(cone[0], cone[1])),
4619:                               edge_type(std::min(cone[1], cone[2]), std::max(cone[1], cone[2])),
4620:                               edge_type(std::min(cone[2], cone[0]), std::max(cone[2], cone[0])),
4621:                               edge_type(std::min(cone[0], cone[3]), std::max(cone[0], cone[3])),
4622:                               edge_type(std::min(cone[1], cone[3]), std::max(cone[1], cone[3])),
4623:                               edge_type(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]))};
4624:         //   Check that vertex does not yet exist
4625:         point_type newVertices[6];

4627:         for(int v = 0; v < 6; ++v) {
4628:           if (edge2vertex.find(edges[v]) == edge2vertex.end()) {
4629:             throw ALE::Exception("Missing edge in refined mesh");
4630:           }
4631:           newVertices[v] = edge2vertex[edges[v]];
4632:           vertex2edge[(newVertices[v]-numNewCells-numVertices)*2+0] = edges[v].first;
4633:           vertex2edge[(newVertices[v]-numNewCells-numVertices)*2+1] = edges[v].second;
4634:         }
4635:         // Create new cells
4636:         const point_type offset = numNewCells - numCells;

4638:         point_type cell0[4] = {cone[0]+offset, newVertices[3], newVertices[0], newVertices[2]};
4639:         newSieve->setCone(cell0, c*8+0);
4640:         point_type cell1[4] = {cone[1]+offset, newVertices[4], newVertices[1], newVertices[0]};
4641:         newSieve->setCone(cell1, c*8+1);
4642:         point_type cell2[4] = {cone[2]+offset, newVertices[5], newVertices[2], newVertices[1]};
4643:         newSieve->setCone(cell2, c*8+2);
4644:         point_type cell3[4] = {cone[3]+offset, newVertices[3], newVertices[5], newVertices[4]};
4645:         newSieve->setCone(cell3, c*8+3);
4646:         point_type cell4[4] = {newVertices[0], newVertices[3], newVertices[4], newVertices[2]};
4647:         newSieve->setCone(cell4, c*8+4);
4648:         point_type cell5[4] = {newVertices[1], newVertices[4], newVertices[5], newVertices[3]};
4649:         newSieve->setCone(cell5, c*8+5);
4650:         point_type cell6[4] = {newVertices[2], newVertices[5], newVertices[3], newVertices[1]};
4651:         newSieve->setCone(cell6, c*8+6);
4652:         point_type cell7[4] = {newVertices[0], newVertices[1], newVertices[2], newVertices[3]};
4653:         newSieve->setCone(cell7, c*8+7);
4654:         cV.clear();
4655:       }
4656:       newSieve->symmetrize();
4657:       // Create new coordinates
4658:       const Obj<typename MeshType::real_section_type>& coordinates    = mesh.getRealSection("coordinates");
4659:       const Obj<typename MeshType::real_section_type>& newCoordinates = newMesh.getRealSection("coordinates");

4661:       newCoordinates->setChart(typename sieve_type::chart_type(numNewCells, curNewVertex));
4662:       for(int v = numNewCells; v < curNewVertex; ++v) {
4663:         newCoordinates->setFiberDimension(v, 3);
4664:       }
4665:       newCoordinates->allocatePoint();
4666:       for(int v = 0; v < numVertices; ++v) {
4667:         newCoordinates->updatePoint(v+numNewCells, coordinates->restrictPoint(v+numCells));
4668:       }
4669:       for(int v = numNewCells+numVertices; v < curNewVertex; ++v) {
4670:         const int     endpointA = vertex2edge[(v-numNewCells-numVertices)*2+0];
4671:         const int     endpointB = vertex2edge[(v-numNewCells-numVertices)*2+1];
4672:         const double *coordsA   = coordinates->restrictPoint(endpointA);
4673:         double        coords[3];

4675:         for(int d = 0; d < 3; ++d) {
4676:           coords[d]  = coordsA[d];
4677:         }
4678:         const double *coordsB = coordinates->restrictPoint(endpointB);
4679:         for(int d = 0; d < 3; ++d) {
4680:           coords[d] += coordsB[d];
4681:           coords[d] *= 0.5;
4682:         }
4683:         newCoordinates->updatePoint(v, coords);
4684:       }
4685:       delete [] vertex2edge;
4686:       // Fast stratification
4687:       const Obj<typename MeshType::label_type>& height = newMesh.createLabel("height");
4688:       const Obj<typename MeshType::label_type>& depth  = newMesh.createLabel("depth");
4689:       for(int c = 0; c < numNewCells; ++c) {
4690:         height->setCone(0, c);
4691:         depth->setCone(1, c);
4692:       }
4693:       for(int v = numNewCells; v < numNewCells+numNewVertices; ++v) {
4694:         height->setCone(1, v);
4695:         depth->setCone(0, v);
4696:       }
4697:       newMesh.setHeight(1);
4698:       newMesh.setDepth(1);
4699:       // Exchange new boundary vertices
4700:       //   We can convert endpoints, and then just match to new vertex on this side
4701:       //   1) Create the overlap of edges which are vertex pairs (do not need for interpolated meshes)
4702:       //   2) Create a section of overlap edge --> new vertex (this will generalize to other split points in interpolated meshes)
4703:       //   3) Copy across new overlap
4704:       //   4) Fuse matches new vertex pairs and inserts them into the old overlap


4707:       // Create the parallel overlap
4708:       int *numCellsP    = new int[mesh.commSize()];
4709:       int *numNewCellsP = new int[newMesh.commSize()];
4710:       int  ierr;

4712:       MPI_Allgather((void *) &numCells, 1, MPI_INT, numCellsP, 1, MPI_INT, mesh.comm());CHKERRXX(ierr);
4713:       MPI_Allgather((void *) &numNewCells, 1, MPI_INT, numNewCellsP, 1, MPI_INT, newMesh.comm());CHKERRXX(ierr);
4714:       Obj<typename MeshType::send_overlap_type> newSendOverlap = newMesh.getSendOverlap();
4715:       Obj<typename MeshType::recv_overlap_type> newRecvOverlap = newMesh.getRecvOverlap();
4716:       const Obj<typename MeshType::send_overlap_type>& sendOverlap = mesh.getSendOverlap();
4717:       const Obj<typename MeshType::recv_overlap_type>& recvOverlap = mesh.getRecvOverlap();
4718:       Obj<typename MeshType::send_overlap_type::traits::capSequence> sendPoints  = sendOverlap->cap();
4719:       const typename MeshType::send_overlap_type::source_type        localOffset = numNewCellsP[newMesh.commRank()] - numCellsP[mesh.commRank()];

4721:       for(typename MeshType::send_overlap_type::traits::capSequence::iterator p_iter = sendPoints->begin(); p_iter != sendPoints->end(); ++p_iter) {
4722:         const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& ranks      = sendOverlap->support(*p_iter);
4723:         const typename MeshType::send_overlap_type::source_type&                  localPoint = *p_iter;

4725:         for(typename MeshType::send_overlap_type::traits::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
4726:           const int                                   rank         = *r_iter;
4727:           const typename MeshType::send_overlap_type::source_type& remotePoint  = r_iter.color();
4728:           const typename MeshType::send_overlap_type::source_type  remoteOffset = numNewCellsP[rank] - numCellsP[rank];

4730:           newSendOverlap->addArrow(localPoint+localOffset, rank, remotePoint+remoteOffset);
4731:         }
4732:       }
4733:       Obj<typename MeshType::recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();

4735:       for(typename MeshType::recv_overlap_type::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
4736:         const Obj<typename MeshType::recv_overlap_type::traits::coneSequence>& ranks      = recvOverlap->cone(*p_iter);
4737:         const typename MeshType::recv_overlap_type::target_type&               localPoint = *p_iter;

4739:         for(typename MeshType::recv_overlap_type::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
4740:           const int                                        rank         = *r_iter;
4741:           const typename MeshType::recv_overlap_type::target_type& remotePoint  = r_iter.color();
4742:           const typename MeshType::recv_overlap_type::target_type  remoteOffset = numNewCellsP[rank] - numCellsP[rank];

4744:           newRecvOverlap->addArrow(rank, localPoint+localOffset, remotePoint+remoteOffset);
4745:         }
4746:       }
4747:       newMesh.setCalculatedOverlap(true);
4748:       delete [] numCellsP;
4749:       delete [] numNewCellsP;
4750:       // Check edges in edge2vertex for both endpoints sent to same process
4751:       //   Put it in section with point being the lowest numbered vertex and value (other endpoint, new vertex)
4752:       Obj<ALE::Section<point_type, edge_type> > newVertices = new ALE::Section<point_type, edge_type>(mesh.comm());
4753:       std::map<edge_type, std::vector<int> > bdedge2rank;

4755:       for(typename std::map<edge_type, point_type>::const_iterator e_iter = edge2vertex.begin(); e_iter != edge2vertex.end(); ++e_iter) {
4756:         const point_type left  = e_iter->first.first;
4757:         const point_type right = e_iter->first.second;

4759:         if (sendOverlap->capContains(left) && sendOverlap->capContains(right)) {
4760:           const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& leftRanksSeq = sendOverlap->support(left);
4761:           std::list<int> leftRanks(leftRanksSeq->begin(), leftRanksSeq->end());
4762:           const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& rightRanks   = sendOverlap->support(right);
4763:           std::list<int> ranks;
4764:           std::set_intersection(leftRanks.begin(), leftRanks.end(), rightRanks->begin(), rightRanks->end(),
4765:                                 std::insert_iterator<std::list<int> >(ranks, ranks.begin()));

4767:           if(ranks.size()) {
4768:             newVertices->addFiberDimension(std::min(e_iter->first.first, e_iter->first.second)+localOffset, 1);
4769:             for(typename std::list<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
4770:               bdedge2rank[e_iter->first].push_back(*r_iter);
4771:             }
4772:           }
4773:         }
4774:       }
4775:       newVertices->allocatePoint();
4776:       const typename ALE::Section<point_type, edge_type>::chart_type& chart = newVertices->getChart();

4778:       for(typename ALE::Section<point_type, edge_type>::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
4779:         typedef typename ALE::Section<point_type, edge_type>::value_type value_type;
4780:         const point_type p      = *c_iter;
4781:         const int        dim    = newVertices->getFiberDimension(p);
4782:         int              v      = 0;
4783:         value_type      *values = new value_type[dim];

4785:         for(typename std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end() && v < dim; ++e_iter) {
4786:           if (std::min(e_iter->first.first, e_iter->first.second)+localOffset == p) {
4787:             values[v++] = edge_type(std::max(e_iter->first.first, e_iter->first.second)+localOffset, edge2vertex[e_iter->first]);
4788:           }
4789:         }
4790:         newVertices->updatePoint(p, values);
4791:         delete [] values;
4792:       }
4793:       // Copy across overlap
4794:       typedef ALE::Pair<int, point_type> overlap_point_type;
4795:       Obj<ALE::Section<overlap_point_type, edge_type> > overlapVertices = new ALE::Section<overlap_point_type, edge_type>(mesh.comm());

4797:       ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVertices, overlapVertices);
4798:       // Merge by translating edge to local points, finding edge in edge2vertex, and adding (local new vetex, remote new vertex) to overlap
4799:       for(typename std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end(); ++e_iter) {
4800:         const point_type localPoint = edge2vertex[e_iter->first];

4802:         for(typename std::vector<int>::const_iterator r_iter = e_iter->second.begin(); r_iter != e_iter->second.end(); ++r_iter) {
4803:           point_type remoteLeft = -1, remoteRight = -1;
4804:           const int  rank       = *r_iter;

4806:           const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& leftRanks = newSendOverlap->support(e_iter->first.first+localOffset);
4807:           for(typename MeshType::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks->begin(); lr_iter != leftRanks->end(); ++lr_iter) {
4808:             if (rank == *lr_iter) {
4809:               remoteLeft = lr_iter.color();
4810:               break;
4811:             }
4812:           }
4813:           const Obj<typename MeshType::send_overlap_type::traits::supportSequence>& rightRanks = newSendOverlap->support(e_iter->first.second+localOffset);
4814:           for(typename MeshType::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks->begin(); rr_iter != rightRanks->end(); ++rr_iter) {
4815:             if (rank == *rr_iter) {
4816:               remoteRight = rr_iter.color();
4817:               break;
4818:             }
4819:           }
4820:           const point_type remoteMin   = std::min(remoteLeft, remoteRight);
4821:           const point_type remoteMax   = std::max(remoteLeft, remoteRight);
4822:           const int        remoteSize  = overlapVertices->getFiberDimension(overlap_point_type(rank, remoteMin));
4823:           const edge_type *remoteVals  = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
4824:           point_type       remotePoint = -1;

4826:           for(int d = 0; d < remoteSize; ++d) {
4827:             if (remoteVals[d].first == remoteMax) {
4828:               remotePoint = remoteVals[d].second;
4829:               break;
4830:             }
4831:           }
4832:           newSendOverlap->addArrow(localPoint, rank, remotePoint);
4833:           newRecvOverlap->addArrow(rank, localPoint, remotePoint);
4834:         }
4835:       }
4836:     };
4837:   };
4838:   class MeshSerializer {
4839:   public:
4840:     template<typename Mesh>
4841:     static void writeMesh(const std::string& filename, Mesh& mesh) {
4842:       std::ofstream fs;

4844:       if (mesh.commRank() == 0) {
4845:         fs.open(filename.c_str());
4846:       }
4847:       writeMesh(fs, mesh);
4848:       if (mesh.commRank() == 0) {
4849:         fs.close();
4850:       }
4851:     };
4852:     template<typename Mesh>
4853:     static void writeMesh(std::ofstream& fs, Mesh& mesh) {
4854:       ISieveSerializer::writeSieve(fs, *mesh.getSieve());
4855:       // Write labels
4856:       const typename Mesh::labels_type& labels = mesh.getLabels();

4858:       if (!mesh.commRank()) {fs << labels.size() << std::endl;}
4859:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
4860:         if (!mesh.commRank()) {fs << l_iter->first << std::endl;}
4861:         LabelSifterSerializer::writeLabel(fs, *l_iter->second);
4862:       }
4863:       // Write sections
4864:       Obj<std::set<std::string> > realNames = mesh.getRealSections();

4866:       if (!mesh.commRank()) {fs << realNames->size() << std::endl;}
4867:       for(std::set<std::string>::const_iterator n_iter = realNames->begin(); n_iter != realNames->end(); ++n_iter) {
4868:         if (!mesh.commRank()) {fs << *n_iter << std::endl;}
4869:         SectionSerializer::writeSection(fs, *mesh.getRealSection(*n_iter));
4870:       }
4871:       Obj<std::set<std::string> > intNames = mesh.getIntSections();

4873:       if (!mesh.commRank()) {fs << intNames->size() << std::endl;}
4874:       for(std::set<std::string>::const_iterator n_iter = intNames->begin(); n_iter != intNames->end(); ++n_iter) {
4875:         if (!mesh.commRank()) {fs << *n_iter << std::endl;}
4876:         SectionSerializer::writeSection(fs, *mesh.getIntSection(*n_iter));
4877:       }
4878:       // Write overlap
4879:       SifterSerializer::writeSifter(fs, *mesh.getSendOverlap());
4880:       SifterSerializer::writeSifter(fs, *mesh.getRecvOverlap());
4881:       // Write distribution overlap
4882:       // Write renumbering
4883:     };
4884:     template<typename Mesh>
4885:     static void loadMesh(const std::string& filename, Mesh& mesh) {
4886:       std::ifstream fs;

4888:       if (mesh.commRank() == 0) {
4889:         fs.open(filename.c_str());
4890:       }
4891:       loadMesh(fs, mesh);
4892:       if (mesh.commRank() == 0) {
4893:         fs.close();
4894:       }
4895:     };
4896:     template<typename Mesh>
4897:     static void loadMesh(std::ifstream& fs, Mesh& mesh) {
4898:       ALE::Obj<typename Mesh::sieve_type> sieve = new typename Mesh::sieve_type(mesh.comm(), mesh.debug());
4899:       PetscErrorCode                      ierr;

4901:       ISieveSerializer::loadSieve(fs, *sieve);
4902:       mesh.setSieve(sieve);
4903:       // Load labels
4904:       int numLabels;

4906:       if (!mesh.commRank()) {fs >> numLabels;}
4907:       MPI_Bcast(&numLabels, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4908:       for(int l = 0; l < numLabels; ++l) {
4909:         ALE::Obj<typename Mesh::label_type> label = new typename Mesh::label_type(mesh.comm(), mesh.debug());
4910:         std::string                         name;
4911:         int                                 len;

4913:         if (!mesh.commRank()) {
4914:           fs >> name;
4915:           len = name.size();
4916:           MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4917:           MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4918:         } else {
4919:           MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4920:           char *n = new char[len+1];
4921:           MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4922:           name = n;
4923:           delete [] n;
4924:         }
4925:         LabelSifterSerializer::loadLabel(fs, *label);
4926:         mesh.setLabel(name, label);
4927:       }
4928:       // Load sections
4929:       int numRealSections;

4931:       if (!mesh.commRank()) {fs >> numRealSections;}
4932:       MPI_Bcast(&numRealSections, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4933:       for(int s = 0; s < numRealSections; ++s) {
4934:         ALE::Obj<typename Mesh::real_section_type> section = new typename Mesh::real_section_type(mesh.comm(), mesh.debug());
4935:         std::string                                name;
4936:         int                                        len;

4938:         if (!mesh.commRank()) {
4939:           fs >> name;
4940:           len = name.size();
4941:           MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4942:           MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4943:         } else {
4944:           MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4945:           char *n = new char[len+1];
4946:           MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4947:           name = n;
4948:           delete [] n;
4949:         }
4950:         SectionSerializer::loadSection(fs, *section);
4951:         mesh.setRealSection(name, section);
4952:       }
4953:       int numIntSections;

4955:       if (!mesh.commRank()) {fs >> numIntSections;}
4956:       MPI_Bcast(&numIntSections, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4957:       for(int s = 0; s < numIntSections; ++s) {
4958:         ALE::Obj<typename Mesh::int_section_type> section = new typename Mesh::int_section_type(mesh.comm(), mesh.debug());
4959:         std::string                               name;
4960:         int                                       len;

4962:         if (!mesh.commRank()) {
4963:           fs >> name;
4964:           len = name.size();
4965:           MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4966:           MPI_Bcast((void *) name.c_str(), len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4967:         } else {
4968:           MPI_Bcast(&len, 1, MPI_INT, 0, mesh.comm());CHKERRXX(ierr);
4969:           char *n = new char[len+1];
4970:           MPI_Bcast(n, len+1, MPI_CHAR, 0, mesh.comm());CHKERRXX(ierr);
4971:           name = n;
4972:           delete [] n;
4973:         }
4974:         SectionSerializer::loadSection(fs, *section);
4975:         mesh.setIntSection(name, section);
4976:       }
4977:       // Load overlap
4978:       SifterSerializer::loadSifter(fs, *mesh.getSendOverlap());
4979:       SifterSerializer::loadSifter(fs, *mesh.getRecvOverlap());
4980:       // Load distribution overlap
4981:       // Load renumbering
4982:     };
4983:   };
4984: } // namespace ALE
4985: #endif