AFEPack
HGeometry.h
浏览该文件的文档。
00001 
00011 #ifndef _HGeometry_h_
00012 #define _HGeometry_h_
00013 
00014 #include <iostream>
00015 #include <fstream>
00016 #include <string>
00017 #include <vector>
00018 #include <list>
00019 #include <iterator>
00020 #include <algorithm>
00021 
00022 #include <base/exceptions.h>
00023 
00024 #include "DerefIterator.h"
00025 #include "Geometry.h"
00026 #include "TemplateElement.h"
00027 
00028 template <int DIM, int DOW> class HGeometry;
00029 template <int DOW> class HGeometry<0,DOW>;
00030 template <int DOW> class HGeometry<1,DOW>;
00031 template <int DOW> class HGeometry<2,DOW>;
00032 template <int DOW> class HGeometry<3,DOW>;
00033 template <int DIM, int DOW> class HGeometryTree;
00034 template <int DIM, int DOW> class RegularMesh;
00035 template <int DIM, int DOW> class IrregularMesh;
00036 template <int DIM, int DOW> class HElement;
00037 
00038 template <int DIM, int DOW> class ElementIterator;
00039 template <int DIM, int DOW> class RootFirstElementIterator;
00040 template <int DIM, int DOW> class ActiveElementIterator;
00041 template <int DIM, int DOW> class IrregularMeshPair;
00042 template <int DIM, int DOW> class ActiveElementPairIterator;
00043 
00044 template <int DIM, int DOW> std::ostream& operator<<(std::ostream&, const HGeometry<DIM,DOW>&);
00045 template <int DIM, int DOW> std::ostream& operator<<(std::ostream&, const HElement<DIM, DOW>&);
00046 
00047 template <int DIM, int DOW> std::ostream& operator<<(std::ostream&, IrregularMesh<DIM, DOW>&);
00048 template <int DOW> std::ostream& operator<<(std::ostream&, const HGeometry<0,DOW>&);
00049 
00050 template <int DIM, int DOW> bool operator==(const ElementIterator<DIM, DOW>&,
00051                                             const ElementIterator<DIM, DOW>&);
00052 template <int DIM, int DOW> bool operator!=(const ElementIterator<DIM, DOW>&,
00053                                             const ElementIterator<DIM, DOW>&);
00054 template <int DIM, int DOW> bool operator==(const ElementIterator<DIM, DOW>&,
00055                                             ElementIterator<DIM, DOW>&);
00056 template <int DIM, int DOW> bool operator!=(const ElementIterator<DIM, DOW>&,
00057                                             ElementIterator<DIM, DOW>&);
00058 
00059 template <int DIM, int DOW> bool operator==(const ActiveElementPairIterator<DIM, DOW>&,
00060                                             const ActiveElementPairIterator<DIM, DOW>&);
00061 template <int DIM, int DOW> bool operator!=(const ActiveElementPairIterator<DIM, DOW>&,
00062                                             const ActiveElementPairIterator<DIM, DOW>&);
00063 template <int DIM, int DOW> bool operator==(const ActiveElementPairIterator<DIM, DOW>&,
00064                                             ActiveElementPairIterator<DIM, DOW>&);
00065 template <int DIM, int DOW> bool operator!=(const ActiveElementPairIterator<DIM, DOW>&,
00066                                             ActiveElementPairIterator<DIM, DOW>&);
00067 
00068 template <int DIM>
00069 struct HGeometryInfo {
00070   enum {
00071     dimension = DIM,
00072     n_vertex = DIM + 1,
00073     n_boundary = DIM + 1,
00074     n_child = (1<<DIM)
00075   };
00076 };
00077 
00078 template <> 
00079 struct HGeometryInfo<0> {
00080   enum {
00081     dimension = 0,
00082     n_vertex = 0,
00083     n_boundary = 0,
00084     n_child = 0
00085   };
00086 };
00087 
00088 template <> 
00089 struct HGeometryInfo<1> {
00090   enum {
00091     dimension = 1,
00092     n_vertex = 2,
00093     n_boundary = 0,
00094     n_child = 2
00095   };
00096 };
00097 
00098 typedef int bmark_t;
00099 
00100 #include "PropertyTable.h"
00101 
00102 #ifdef __SERIALIZATION__
00103 #include "Migration.h"
00104 struct HGeometryBase : public PropertyTable, public Migration::HBuffer
00105 #else
00106 struct HGeometryBase : public PropertyTable
00107 #endif // __SERIALIZATION__
00108 {
00109   virtual ~HGeometryBase() {};
00110 };
00111 
00121 struct HTools {
00122 
00123   enum { USED = -8, UNUSED = -7,
00124          INACTIVE = -9, ACTIVE = USED };
00125 
00127 
00135   template <class GEO> bool
00136   isGeometryUsed(const GEO& geo) const {
00137     return (geo.index == USED);
00138   }
00139 
00143   template <class GEO> bool
00144   isGeometryActive(const GEO& geo) const {
00145     return (geo.index == ACTIVE);
00146   }
00147 
00151   template <class GEO> bool
00152   isGeometryInactive(const GEO& geo) const {
00153     return (geo.index == INACTIVE);
00154   }
00155 
00159   template <class GEO> bool
00160   isGeometryIndexed(const GEO& geo) const {
00161     return (geo.index >= 0);
00162   }
00164 
00166 
00174   template <class GEO> void 
00175   setGeometryUsed(GEO& geo) const {
00176     geo.index = USED;
00177     for (u_int i = 0;i < GEO::n_boundary;++ i) {
00178       this->setGeometryUsed(*geo.boundary[i]);
00179     }
00180   }
00181   template <int DOW> void 
00182   setGeometryUsed(HGeometry<0,DOW>& geo) const {
00183     geo.index = USED;
00184   }
00185   template <int DOW> void 
00186   setGeometryUsed(HGeometry<1,DOW>& geo) const {
00187     geo.index = USED;
00188   }
00189   template <class GEO> void 
00190   setGeometryActive(GEO& geo) const { 
00191     setGeometryUsed(geo);
00192   }
00194 
00196 
00199   template <class GEO> void 
00200   setGeometryUnused(GEO& geo) const {
00201     geo.index = UNUSED;
00202     for (u_int i = 0;i < GEO::n_boundary;++ i) {
00203       this->setGeometryUnused(*geo.boundary[i]);
00204     }
00205   }
00206   template <int DOW> void 
00207   setGeometryUnused(HGeometry<1,DOW>& geo) const {
00208     geo.index = UNUSED;
00209   }
00211 
00213 
00216   template <class GEO> void 
00217   setGeometryInactive(GEO& geo) const {
00218     geo.index = INACTIVE;
00219     for (u_int i = 0;i < GEO::n_boundary;++ i) {
00220       this->setGeometryInactive(*geo.boundary[i]);
00221     }
00222   }
00223   template <int DOW> void 
00224   setGeometryInactive(HGeometry<1,DOW>& geo) const {
00225     geo.index = INACTIVE;
00226   }
00228 
00232   template <class GEO> void 
00233   setGeometryUnusedRecursively(GEO& geo) const {
00234     setGeometryUnused(geo);
00235 
00236     if (geo.isRefined()) {
00237       for (u_int i = 0;i < GEO::n_child;++ i) {
00238         this->setGeometryUnusedRecursively(*geo.child[i]);
00239       }
00240     }
00241   }
00243 
00247   template <class GEO> bool 
00248   isRefined(const GEO& geo) const {
00249     if (geo.isRefined()) {
00250       for (int n = 0;n < geo.n_child;++ n) {
00251         if (this->isGeometryUsed(*geo.child[n])) {
00252           return true;
00253         }
00254       }
00255     }
00256 
00257     return false;
00258   }
00259 
00266   template <int DOW> bool 
00267   isSemiregular(const HGeometry<1,DOW>& geo) const {
00268     assert (this->isGeometryUsed(geo));
00269 
00270     bool result = true;
00271     if (geo.isRefined()) {
00272       if (this->isRefined(*geo.child[0]) ||
00273           this->isRefined(*geo.child[1])) {
00274         result = false;
00275       }
00276     }
00277     return result;
00278   }
00279 
00288   template <int DOW> bool 
00289   isSemiregular(const HGeometry<2,DOW>& geo) const {
00290     assert (this->isGeometryUsed(geo));
00291 
00292     u_int n_refined_edge = 0;
00293     for (u_int i = 0;i < geo.n_boundary;++ i) {
00294       const HGeometry<1,DOW>& edge = *geo.boundary[i];
00295       if(! this->isSemiregular(edge)) return false;
00296       if (this->isRefined(edge)) {
00297         n_refined_edge += 1;
00298       }
00299     }
00300     bool result = (n_refined_edge <= 1);
00301     return result;
00302   }
00303 
00304   template <class HGEO> bool
00305   semiregularizeBoundary(HGEO& geo) const {
00306     return false;
00307   }
00308 
00309   template <int DOW> bool
00310   semiregularizeBoundary(HGeometry<3,DOW>& geo) const {
00311     assert (this->isGeometryUsed(geo));
00312 
00313     bool result = false;
00314     for (int i = 0;i < geo.n_boundary;++ i) {
00315       HGeometry<2,DOW>& face = *geo.boundary[i];
00316       u_int n_refined_face_edge = 0;
00317       for (int j = 0;j < face.n_boundary;++ j) {
00318         const HGeometry<1,DOW> * edge = face.boundary[j];
00319         if (this->isRefined(*edge)) {
00320           n_refined_face_edge += 1;
00321         }
00322       }
00323 
00324       bool is_operated = false;
00325       if (n_refined_face_edge >= 2) {
00326         if (! face.isRefined()) face.refine();
00327         
00328         for (int j = 0;j < face.n_child;++ j) {
00329           HGeometry<2,DOW>& chd = *(face.child[j]);
00330           if (this->isGeometryUsed(chd)) continue;
00331 
00332           this->setGeometryUsed(chd);
00333           is_operated = true;
00334         }
00335       }
00336       if (is_operated) {
00337         i = -1;
00338         result = true;
00339       }
00340     }
00341     return result;
00342   }
00343 
00351   template <int DOW> bool 
00352   isSemiregular(const HGeometry<3,DOW>& geo) const {
00353     assert (this->isGeometryUsed(geo));
00354 
00355     u_int n_refined_edge = 0;
00356     u_int n_refined_face = 0;
00357     for (int i = 0;i < geo.n_boundary;++ i) {
00358       const HGeometry<2,DOW>& face = *geo.boundary[i];
00359       u_int n_refined_face_edge = 0;
00360       for (int j = 0;j < face.n_boundary;++ j) {
00361         const HGeometry<1,DOW> * edge = face.boundary[j];
00362         if (! this->isSemiregular(*edge)) return false;
00363         if (this->isRefined(*edge)) {
00364           n_refined_face_edge += 1;
00365         }
00366       }
00367 
00368       if (n_refined_face_edge == 3) { 
00369         n_refined_face += 1;
00371         const HGeometry<2,DOW> * chd3 = face.child[3];
00372         for (int j = 0;j < chd3->n_boundary;++ j) {
00373           const HGeometry<1,DOW> * edge = chd3->boundary[j];
00374           if (this->isRefined(*edge)) return false;
00375         }
00376       }
00377 
00378       n_refined_edge += n_refined_face_edge;
00379     }
00380     if (n_refined_edge%2 != 0) abort();
00381     n_refined_edge /= 2; // 每条边都被算了两遍
00382     bool result = (n_refined_edge <= 1 || 
00383                    (n_refined_edge == 3 && 
00384                     n_refined_face == 1));
00385     return result;
00386   }
00387 
00391   template <class GEO> void
00392   setParentInactive(GEO& geo) const {
00393     if (geo.parent == NULL) return;
00394     setGeometryInactive(*geo.parent);
00395   }
00396 
00398 
00411   template <class GEO> void
00412   clearIndex(GEO& geo) const {
00413     geo.index = 0;
00414     for (u_int i = 0;i < geo.n_boundary;++ i) {
00415       clearIndex(*geo.boundary[i]);
00416     }
00417     if (geo.isRefined()) {
00418       for (u_int i = 0;i < geo.n_child;++ i) {
00419         clearIndex(*geo.child[i]);
00420       }
00421     }
00422   }
00423   template <int DOW> void
00424   clearIndex(HGeometry<0,DOW>& geo) const {
00425     geo.index = 0;
00426   }
00427   template <int DOW> void
00428   clearIndex(HGeometry<1,DOW>& geo) const {
00429     geo.index = 0;
00430     for (u_int i = 0;i < geo.n_vertex;++ i) {
00431       clearIndex(*geo.vertex[i]);
00432     }
00433     if (geo.isRefined()) {
00434       for (u_int i = 0;i < geo.n_child;++ i) {
00435         clearIndex(*geo.child[i]);
00436       }
00437     }
00438   }
00440 
00442 
00445   template <class GEO> void
00446   incrIndex(GEO& geo) const {
00447     geo.index += 1;
00448     for (u_int i = 0;i < geo.n_boundary;++ i) {
00449       incrIndex(*geo.boundary[i]);
00450     }
00451     if (geo.isRefined()) {
00452       for (u_int i = 0;i < geo.n_child;++ i) {
00453         incrIndex(*geo.child[i]);
00454       }
00455     }
00456   }
00457   template <int DOW> void
00458   incrIndex(HGeometry<0,DOW>& geo) const {
00459     geo.index += 1;
00460   }
00461   template <int DOW> void
00462   incrIndex(HGeometry<1,DOW>& geo) const {
00463     geo.index += 1;
00464     for (u_int i = 0;i < geo.n_vertex;++ i) {
00465       incrIndex(*geo.vertex[i]);
00466     }
00467     if (geo.isRefined()) {
00468       for (u_int i = 0;i < geo.n_child;++ i) {
00469         incrIndex(*geo.child[i]);
00470       }
00471     }
00472   }
00474 
00476 
00479   template <class GEO> void
00480   decrIndex(GEO& geo) const {
00481     geo.index -= 1;
00482     if (geo.isRefined()) {
00483       for (u_int i = 0;i < geo.n_child;++ i) {
00484         decrIndex(*geo.child[i]);
00485       }
00486     }
00487     for (u_int i = 0;i < geo.n_boundary;++ i) {
00488       decrIndex(*geo.boundary[i]);
00489     }
00490     if (geo.index == 0) delete &geo;
00491   }
00492   template <int DOW> void
00493   decrIndex(HGeometry<0,DOW>& geo) const {
00494     geo.index -= 1;
00495     if (geo.index == 0) delete &geo;
00496   }
00497   template <int DOW> void
00498   decrIndex(HGeometry<1,DOW>& geo) const {
00499     geo.index -= 1;
00500     if (geo.isRefined()) {
00501       for (u_int i = 0;i < geo.n_child;++ i) {
00502         decrIndex(*geo.child[i]);
00503       }
00504     }
00505     for (u_int i = 0;i < geo.n_vertex;++ i) {
00506       decrIndex(*geo.vertex[i]);
00507     }
00508     if (geo.index == 0) delete &geo;
00509   }
00511 
00512 
00514   template <class HGEO, class VTX> void
00515   regularize_add_node(HGEO& hgeo, 
00516                       GeometryBM& geo,
00517                       VTX& vtx) const {
00518     int idx = hgeo.index;
00519     if (geo.index() == -1) {
00520       vtx = hgeo; 
00521 
00522       geo.index() = idx;
00523       geo.vertex().resize(1, idx);
00524       geo.boundary().resize(1, idx);
00525       geo.boundaryMark() = hgeo.bmark;
00526     }
00527   }
00528 
00529   template <class HGEO> void
00530   regularize_add_side(HGEO& hgeo, 
00531                       GeometryBM& geo) const {
00532     int idx = hgeo.index;
00533     if (geo.index() == -1) {
00534       geo.index() = idx;
00535 
00536       geo.vertex().resize(2);
00537       geo.vertex(0) = hgeo.vertex[0]->index;
00538       geo.vertex(1) = hgeo.vertex[1]->index;
00539 
00540       geo.boundary().resize(2);
00541       geo.boundary(0) = hgeo.vertex[0]->index;
00542       geo.boundary(1) = hgeo.vertex[1]->index;
00543 
00544       geo.boundaryMark() = hgeo.bmark;
00545     }
00546   }
00547 
00548   template <class HGEO> void
00549   regularize_add_triangle(HGEO& hgeo, 
00550                           GeometryBM& geo) const {
00551     int idx = hgeo.index;
00552     if (geo.index() == -1) {
00553       geo.index() = idx;
00554 
00555       geo.vertex().resize(3);
00556       geo.vertex(0) = hgeo.vertex[0]->index;
00557       geo.vertex(1) = hgeo.vertex[1]->index;
00558       geo.vertex(2) = hgeo.vertex[2]->index;
00559 
00560       geo.boundary().resize(3);
00561       geo.boundary(0) = hgeo.boundary[0]->index;
00562       geo.boundary(1) = hgeo.boundary[1]->index;
00563       geo.boundary(2) = hgeo.boundary[2]->index;
00564 
00565       geo.boundaryMark() = hgeo.bmark;
00566     }
00567   }
00568 
00569   template <class HGEO> void
00570   regularize_add_twin_triangle(HGEO& hgeo, 
00571                                GeometryBM& geo, 
00572                                int k) const {
00573     int idx = hgeo.index;
00574     int ii[] = {0,1,2,0,1,2,0,1,2};
00575     if (geo.index() == -1) {
00576       geo.index() = idx;
00577 
00578       geo.vertex().resize(4);
00579       geo.vertex(0) = hgeo.vertex[k]->index;
00580       geo.vertex(1) = hgeo.vertex[ii[k + 1]]->index;
00581       geo.vertex(2) = hgeo.boundary[k]->child[0]->vertex[1]->index;
00582       geo.vertex(3) = hgeo.vertex[ii[k + 2]]->index;
00583 
00584       geo.boundary().resize(4);
00585       geo.boundary(0) = hgeo.boundary[ii[k + 2]]->index;
00586       if (hgeo.boundary[k]->child[0]->vertex[0] == hgeo.vertex[ii[k + 1]]) {
00587         geo.boundary(1) = hgeo.boundary[k]->child[0]->index;
00588         geo.boundary(2) = hgeo.boundary[k]->child[1]->index;
00589       } else {
00590         geo.boundary(1) = hgeo.boundary[k]->child[1]->index;
00591         geo.boundary(2) = hgeo.boundary[k]->child[0]->index;
00592       }
00593       geo.boundary(3) = hgeo.boundary[ii[k + 1]]->index;
00594 
00595       geo.boundaryMark() = hgeo.bmark;
00596     }
00597   }
00599 };
00600 
00606 template <int DIM, int DOW=DIM>
00607   class HGeometry : public HGeometryInfo<DIM>, public HGeometryBase
00608   {
00609   public:
00610   enum { dim = DIM, dow = DOW };
00611   typedef HGeometry<0,dow> vertex_t;
00612   typedef HGeometry<dim-1,dow> bound_t;
00613   typedef HGeometry<dim,dow> this_t;
00614   typedef this_t child_t;
00615   typedef this_t parent_t;
00616 
00617   int index;
00618   std::vector<vertex_t *> vertex;
00619   std::vector<bound_t *> boundary;
00620   parent_t * parent;
00621   std::vector<child_t *> child;
00622   bmark_t bmark;
00623 
00624   public:
00625   HGeometry();
00626   virtual ~HGeometry() {}
00627 
00628   public:
00629   bool isRefined() const;
00630   bool isIncludePoint(const afepack::Point<DOW>&) const;
00631   void refine();
00632   void checkIntegrity() const;
00633 
00634   friend std::ostream& operator<< <>(std::ostream&, const HGeometry<DIM,DOW>&);
00635   };
00636 
00642 template <int DOW>
00643 class HGeometry<0,DOW> : public afepack::Point<DOW>, public HGeometryInfo<0>, public HGeometryBase
00644 {
00645  public:
00646   enum { dim = 0, dow = DOW };
00647   typedef HGeometry<dim,dow> this_t;
00648   typedef this_t vertex_t;
00649   typedef this_t bound_t;
00650   typedef this_t child_t;
00651   typedef this_t parent_t;
00652 
00653   int index;
00654   bmark_t bmark;
00655 
00656   static parent_t * parent;
00657   static std::vector<vertex_t *> vertex;
00658   static std::vector<bound_t *> boundary;
00659   static std::vector<child_t *> child;
00660  public:
00661   HGeometry();
00662   virtual ~HGeometry() {}
00663 
00664   bool isRefined() const { return false; }
00665   void refine() {}
00666 };
00667 
00668 
00669 
00670 template <int DOW=1>
00671   class HGeometry<1,DOW> : public HGeometryInfo<1>, public HGeometryBase
00672   {
00673   public:
00674   enum { dim = 1, dow = DOW };
00675   typedef HGeometry<0,dow> vertex_t;
00676   typedef HGeometry<dim-1,dow> bound_t;
00677   typedef HGeometry<dim,dow> this_t;
00678   typedef this_t child_t;
00679   typedef this_t parent_t;
00680 
00681   static void (*mid_point)(const afepack::Point<DOW>&, 
00682                            const afepack::Point<DOW>&,
00683                            bmark_t,
00684                            afepack::Point<DOW>&);
00685   public:
00686   int index;
00687   std::vector<vertex_t *> vertex;
00688   std::vector<bound_t *> boundary;
00689   parent_t * parent;
00690   std::vector<HGeometry<1,DOW> *> child;
00691   bmark_t bmark;
00692   public:
00693   HGeometry();
00694   virtual ~HGeometry() {}
00695   public:
00696   bool isRefined() const;
00697   bool isIncludePoint(const afepack::Point<DOW>&) const;
00698   void refine();
00699   void checkIntegrity() const;
00700 
00701   friend std::ostream& operator<< <>(std::ostream&, const HGeometry<1,DOW>&);
00702   };
00703 
00704 template <int DOW=2>
00705   class HGeometry<2,DOW> : public HGeometryInfo<2>, public HGeometryBase
00706   {
00707   public:
00708   enum { dim = 2, dow = DOW };
00709   typedef HGeometry<0,dow> vertex_t;
00710   typedef HGeometry<dim-1,dow> bound_t;
00711   typedef HGeometry<dim,dow> this_t;
00712   typedef this_t child_t;
00713   typedef this_t parent_t;
00714 
00715   int index;
00716   std::vector<vertex_t *> vertex;
00717   std::vector<bound_t *> boundary;
00718   parent_t * parent;
00719   std::vector<child_t *> child;
00720   bmark_t bmark;
00721   public:
00722   HGeometry();
00723   virtual ~HGeometry() {}
00724   public:
00725   bool isRefined() const;
00726   bool isIncludePoint(const afepack::Point<DOW>&) const;
00727   void refine();
00728   void checkIntegrity() const;
00729 
00733   static double 
00734   triangle_area(const afepack::Point<DOW>& v0,
00735                 const afepack::Point<DOW>& v1,
00736                 const afepack::Point<DOW>& v2) {
00737     return 0.5*((v1[0] - v0[0])*(v2[1] - v0[1]) - 
00738                 (v1[1] - v0[1])*(v2[0] - v0[0]));
00739   }
00740 
00741   friend std::ostream& operator<< <>(std::ostream&, const HGeometry<2,DOW>&);
00742   };
00743 
00744 template <int DOW=3>
00745   class HGeometry<3,DOW> : public HGeometryInfo<3>, public HGeometryBase
00746   {
00747   public:
00748   enum { dim = 3, dow = DOW };
00749   typedef HGeometry<0,dow> vertex_t;
00750   typedef HGeometry<dim-1,dow> bound_t;
00751   typedef HGeometry<dim,dow> this_t;
00752   typedef this_t child_t;
00753   typedef this_t parent_t;
00754 
00755   static const int      REFINE_MODEL_03 = 0;
00756   static const int      REFINE_MODEL_14 = 1;
00757   static const int      REFINE_MODEL_25 = 2;
00758   int                   refine_model;
00759   public:
00760   int index;
00761   std::vector<vertex_t *> vertex;
00762   std::vector<bound_t *> boundary;
00763   parent_t * parent;
00764   std::vector<child_t *> child;
00765   bmark_t bmark;
00766   public:
00767   HGeometry();
00768   virtual ~HGeometry() {}
00769   public:
00770   bool isRefined() const;
00771   bool isIncludePoint(const afepack::Point<DOW>&) const;
00772   void refine();
00773   void checkIntegrity() const;
00774 
00778   static double
00779   tetrahedron_volume(const afepack::Point<DOW>& v0,
00780                      const afepack::Point<DOW>& v1,
00781                      const afepack::Point<DOW>& v2,
00782                      const afepack::Point<DOW>& v3) {
00783     return ((v1[0] - v0[0])*(v2[1] - v0[1])*(v3[2] - v0[2]) +
00784             (v1[1] - v0[1])*(v2[2] - v0[2])*(v3[0] - v0[0]) +
00785             (v1[2] - v0[2])*(v2[0] - v0[0])*(v3[1] - v0[1]) -
00786             (v1[0] - v0[0])*(v2[2] - v0[2])*(v3[1] - v0[1]) -
00787             (v1[1] - v0[1])*(v2[0] - v0[0])*(v3[2] - v0[2]) -
00788             (v1[2] - v0[2])*(v2[1] - v0[1])*(v3[0] - v0[0]))/6.;
00789   }
00790 
00791   friend std::ostream& operator<< <>(std::ostream&, const HGeometry<3,DOW>&);
00792   };
00793 
00794 
00799 template <int DIM, int DOW=DIM>
00800   class HGeometryTree
00801   {
00802   public:
00803   enum { dim = DIM, dow = DOW };
00804 
00805   private:
00806   typedef HGeometry<DIM,DOW> entry_t;
00807   public:
00808   typedef std::list<entry_t *> container_t;
00809   private:
00810   container_t root_element;
00811 
00812   bool _is_locked; 
00819   public:
00820   typedef _Deref_iterator<typename container_t::iterator, entry_t> RootIterator;
00821   typedef _Deref_iterator<typename container_t::const_iterator, const entry_t> ConstRootIterator;
00822 
00823   typedef HTools Tools;
00824   public:
00825   HGeometryTree() : _is_locked(false) {};
00826   virtual ~HGeometryTree() {clear();};
00827 
00828   protected:
00829   bool lock() {
00830     if (_is_locked) return false;
00831     else {
00832       _is_locked = true;
00833       return true;
00834     }
00835   }
00836   void unlock() {
00837     _is_locked = false;
00838   }
00839 
00840   public:
00841   container_t& rootElement() { return root_element; }
00842   const container_t& rootElement() const { return root_element; }
00843 
00844   unsigned int n_rootElement() const {return root_element.size();}
00845 
00846   RootIterator beginRootElement() {return RootIterator(root_element.begin());}
00847   RootIterator endRootElement() {return RootIterator(root_element.end());}
00848 
00849   ConstRootIterator beginRootElement() const {return ConstRootIterator(root_element.begin());}
00850   ConstRootIterator endRootElement() const {return ConstRootIterator(root_element.end());}
00851 
00852   void clear(); // this method need very complex implementation, left for future
00853   void checkIntegrity();
00854 
00855   bool is_locked() const { return _is_locked; }
00856   bool& is_locked() { return _is_locked; }
00857         
00862   void readEasyMesh(const std::string&);
00863         
00867   void readMesh(const std::string&);
00868 
00869   friend class IrregularMesh<DIM, DOW>;
00870   };
00871 
00872 
00873 
00881 template <int DIM, int DOW=DIM>
00882   class HElement : public HGeometryInfo<DIM>, public HGeometryBase {
00883  public:
00884  enum { dim = DIM, dow = DOW };
00885  typedef HGeometry<dim,dow> h_element_t;
00886  typedef HElement<dim,dow> element_t;
00887  typedef element_t parent_t;
00888  typedef element_t child_t;
00889 
00890  public:
00891  typedef int                            ElementType;
00892  static const ElementType               NOT_ACTIVE = -1;
00893         
00894  // for 2 dimensional case
00895  static const ElementType               TRIANGLE = 0;
00896  static const ElementType               QUADRILATERAL = 1;
00897         
00898  // for 3 dimensional case
00899  static const ElementType               TETRAHEDRON = 0;
00900  static const ElementType               TWIN_TETRAHEDRON = 1;
00901  static const ElementType               FOUR_TETRAHEDRON = 2;
00902  public:
00903  int                                    index;
00904  double                                 indicator;
00905  int                                    value; // default: -1
00906  h_element_t *                          h_element; // default: NULL
00907  element_t *                            parent; // default: NULL
00908  std::vector<element_t *>               child; // default: NULL
00909  public:
00910  HElement();
00911  HElement(const element_t&);
00912  virtual ~HElement();
00913  public:
00914  element_t& operator=(const element_t&);
00915  /* bool is_dummy() const { return h_element->is_dummy(); } */
00916  bool isRefined() const;
00917  bool isIncludePoint(const afepack::Point<DOW>&) const;
00918  void refine();
00919  void checkIntegrity() const;
00920 
00921  friend std::ostream& operator<< <>(std::ostream&, const HElement<DIM, DOW>&);
00922 };
00923 
00924 
00925 
00937 template <int DIM, int DOW=DIM>
00938   class IrregularMesh
00939   {
00940   public:
00941   enum { dim = DIM, dow = DOW };
00942   typedef RegularMesh<DIM,DOW> mesh_t;
00943   typedef HGeometryTree<DIM,DOW> tree_t;
00944   typedef IrregularMesh<DIM,DOW> ir_mesh_t;
00945 
00946   protected:
00947   typedef HGeometry<DIM,DOW> h_element_t;
00948   typedef HElement<DIM,DOW> element_t;
00949   typedef std::list<element_t *> container_t;
00950   typedef HTools Tools;
00951 
00952   private:
00953   tree_t * geometry_tree; // default: NULL
00954   container_t root_element; // default: NULL
00955   mesh_t * regular_mesh; // default: NULL
00956 
00957   public:
00958   typedef _Deref_iterator<typename container_t::iterator, element_t> RootIterator;
00959   typedef _Deref_iterator<typename container_t::const_iterator, const element_t> ConstRootIterator;
00960 
00961   typedef RootFirstElementIterator<DIM, DOW> RootFirstIterator;
00962   typedef ActiveElementIterator<DIM, DOW> ActiveIterator;
00963 
00964   public:
00965   IrregularMesh();
00966   explicit IrregularMesh(tree_t&);
00967   IrregularMesh(const ir_mesh_t&);
00968   virtual ~IrregularMesh();
00969 
00970   public:
00971   void clear();
00972   ir_mesh_t& operator=(const ir_mesh_t&);
00973 
00974   public:
00975   RootIterator beginRootElement() {return root_element.begin();}
00976   RootIterator endRootElement() {return root_element.end();}
00977 
00978   ConstRootIterator beginRootElement() const {return root_element.begin();};
00979   ConstRootIterator endRootElement() const {return root_element.end();};
00980 
00981   RootFirstIterator beginRootFirstElement();
00982   RootFirstIterator endRootFirstElement();
00983 
00984   ActiveIterator beginActiveElement();
00985   ActiveIterator endActiveElement();
00986 
00987   public:
00999   void reinit(tree_t& htree, bool is_bare = false);
01000   tree_t& geometryTree() const {return *geometry_tree;};
01001   container_t& rootElement() {return root_element;};
01002   mesh_t& regularMesh() {return *regular_mesh;};
01003   const mesh_t& regularMesh() const {return *regular_mesh;};
01004 
01005   virtual void semiregularize();
01006   void regularize(bool renumerate = true);
01007   void globalRefine(unsigned int i = 1);
01008   void randomRefine(double percent = 50.0);
01009   void writeFormatted(const std::string&);
01010 
01011   void copyTree(const ir_mesh_t&);
01012   void copyNonnegtiveSubtree(const ir_mesh_t&);
01013 
01014   void copyTree(const element_t *, element_t *);
01015   void copyNonnegtiveSubtree(const element_t *, element_t *);
01016 
01017   void deleteTree(element_t *);
01018 
01019   friend std::ostream& operator<< <>(std::ostream&, IrregularMesh<DIM, DOW>&);
01020   friend class IrregularMeshPair<DIM, DOW>;
01021 
01022   protected:
01023   void checkIntegrity();
01024   void setGeometryTree(tree_t *);
01025 
01026   void semiregularizeHelper(bool&, int&);
01027   void semiregularizeHelper(bool&, element_t&, int&);
01028 
01029   void prepareSemiregularize();
01030   void prepareSemiregularizeHelper(h_element_t *);
01031 
01032   void renumerateElement();
01033 
01034   void refineElement(element_t& h_ele);
01035 
01036   public:
01037   friend class RegularMesh<DIM, DOW>;
01038   };
01039 
01040 
01041 
01047 template <int DIM, int DOW=DIM>
01048   class RegularMesh : public Mesh<DIM,DOW>
01049   {
01050   public:
01051   enum { dim = DIM, dow = DOW };
01052   typedef IrregularMesh<DIM,DOW> ir_mesh_t;
01053   typedef RegularMesh<DIM,DOW> mesh_t;
01054 
01055   private:
01056   ir_mesh_t * irregular_mesh;
01057 
01058 #ifdef __SERIALIZATION__
01059   std::vector<std::vector<HGeometryBase*> > h_geometry_ptr; 
01060 #endif
01061 
01062   private: 
01063   RegularMesh() {}
01064   RegularMesh(const mesh_t& m) {}
01065   mesh_t& operator=(const mesh_t& m) {}
01066   explicit RegularMesh(ir_mesh_t * m) : irregular_mesh(m) {};
01067   public:
01068   ir_mesh_t& irregularMesh() const {return *irregular_mesh;};
01069   void renumerateElement() {irregular_mesh->renumerateElement();};
01070 
01071 #ifdef __SERIALIZATION__
01072   std::vector<std::vector<HGeometryBase*> >& h_geometry() { return h_geometry_ptr;}
01073   const std::vector<std::vector<HGeometryBase*> >& h_geometry() const { return h_geometry_ptr;}
01074   HGeometryBase * h_geometry(int dim, u_int idx) const { return h_geometry_ptr[dim][idx]; }
01075 
01079   template <int GDIM> HGeometry<GDIM,DOW> *
01080   h_geometry(u_int idx) const { 
01081     return (HGeometry<GDIM,DOW> *)h_geometry_ptr[GDIM][idx]; 
01082   }
01083 
01093   template <class T, int GDIM> T * 
01094   new_property(u_int idx, const property_id_t<T>& pid) const {
01095     return this->template h_geometry<GDIM>(idx)->new_property(pid);
01096   }
01106   template <class T, int GDIM> T * 
01107   get_property(int idx, const property_id_t<T>& pid) const {
01108     return this->template h_geometry<GDIM>(idx)->get_property(pid);
01109   }
01110   template <class T, int GDIM> void
01111   free_property(int idx, const property_id_t<T>& pid) const {
01112     this->template h_geometry<GDIM>(idx)->free_property(pid);
01113   }
01114 
01118   template <class T, int GDIM> T * 
01119   new_property(const GeometryBM& geo, const property_id_t<T>& pid) const {
01120     return this->template h_geometry<GDIM>(geo.index())->new_property(pid);
01121   }
01125   template <class T, int GDIM> T * 
01126   get_property(const GeometryBM& geo, const property_id_t<T>& pid) const {
01127     return this->template h_geometry<GDIM>(geo.index())->get_property(pid);
01128   }
01129   template <class T, int GDIM> void
01130   free_property(const GeometryBM& geo, const property_id_t<T>& pid) const {
01131     this->template h_geometry<GDIM>(geo.index())->free_property(pid);
01132   }
01133 
01137   template <class T> T *
01138   new_property(int dim, 
01139                const GeometryBM& geo, 
01140                const property_id_t<T>& pid) const {
01141     switch(dim) {
01142     case 0: return this->template new_property<T,0>(geo, pid);
01143     case 1: return this->template new_property<T,1>(geo, pid);
01144     case 2: return this->template new_property<T,2>(geo, pid);
01145     case 3: return this->template new_property<T,3>(geo, pid);
01146     }
01147   }
01157   template <class T> T *
01158   get_property(int dim, 
01159                int idx, 
01160                const property_id_t<T>& pid) const {
01161     switch(dim) {
01162     case 0: return this->template get_property<T,0>(idx, pid);
01163     case 1: return this->template get_property<T,1>(idx, pid);
01164     case 2: return this->template get_property<T,2>(idx, pid);
01165     case 3: return this->template get_property<T,3>(idx, pid);
01166     }
01167   }
01168   template <class T> void
01169   free_property(int dim, 
01170                          int idx, 
01171                          const property_id_t<T>& pid) const {
01172     switch(dim) {
01173     case 0: this->template free_property<T,0>(idx, pid); break;
01174     case 1: this->template free_property<T,1>(idx, pid); break;
01175     case 2: this->template free_property<T,2>(idx, pid); break;
01176     case 3: this->template free_property<T,3>(idx, pid); break;
01177     }
01178   }
01179 
01189   template <class T> T *
01190   new_property(int dim, 
01191                int idx, 
01192                const property_id_t<T>& pid) const {
01193     switch(dim) {
01194     case 0: return this->template new_property<T,0>(idx, pid);
01195     case 1: return this->template new_property<T,1>(idx, pid);
01196     case 2: return this->template new_property<T,2>(idx, pid);
01197     case 3: return this->template new_property<T,3>(idx, pid);
01198     }
01199   }
01203   template <class T> T *
01204   get_property(int dim, 
01205                const GeometryBM& geo, 
01206                const property_id_t<T>& pid) const {
01207     switch(dim) {
01208     case 0: return this->template get_property<T,0>(geo, pid);
01209     case 1: return this->template get_property<T,1>(geo, pid);
01210     case 2: return this->template get_property<T,2>(geo, pid);
01211     case 3: return this->template get_property<T,3>(geo, pid);
01212     }
01213   }
01214   template <class T> void
01215   free_property(int dim, 
01216                 const GeometryBM& geo, 
01217                 const property_id_t<T>& pid) const {
01218     switch(dim) {
01219     case 0: this->template free_property<T,0>(geo, pid); break;
01220     case 1: this->template free_property<T,1>(geo, pid); break;
01221     case 2: this->template free_property<T,2>(geo, pid); break;
01222     case 3: this->template free_property<T,3>(geo, pid); break;
01223     }
01224   }
01225 
01226 #endif // __SERIALIZATION__
01227 
01228   template <class LABEL>
01229   void refineLabelled(LABEL& flag) {
01230     typename ir_mesh_t::ActiveIterator
01231     the_ele = irregular_mesh->beginActiveElement(),
01232     end_ele = irregular_mesh->endActiveElement();
01233     for (;the_ele != end_ele;) {
01234       typename ir_mesh_t::ActiveIterator it = the_ele;
01235       ++ the_ele;
01236       if (flag[it->index]) {
01237         irregular_mesh->refineElement(*it);
01238       }
01239     }
01240   }
01241   template <class LABEL>
01242   void coarsenLabelled(LABEL& flag) {
01243     typedef HElement<DIM,DOW> h_element_t;
01244     std::list<h_element_t *> coarsen_list;
01245     property_id_t<int> pid;
01246     new_property_id(pid);
01247     typename ir_mesh_t::ActiveIterator
01248     the_ele = irregular_mesh->beginActiveElement(),
01249     end_ele = irregular_mesh->endActiveElement();
01250     for (;the_ele != end_ele;++ the_ele) {
01251       if (flag[the_ele->index]) {
01252         HElement<DIM,DOW> * p_ele = the_ele->parent;
01253         if (p_ele != NULL) {
01254           int * p_prp = p_ele->get_property(pid);
01255           if (p_prp == NULL) {
01256             p_prp = p_ele->new_property(pid);
01257             *p_prp = 1;
01258           } else {
01259             *p_prp += 1;
01260           }
01261           if (*p_prp == p_ele->n_child) {
01262             coarsen_list.push_back(p_ele);
01263             p_ele->free_property(pid);
01264           }
01265         }
01266       }
01267     }
01268     free_property_id(pid);
01269 
01270     typename std::list<h_element_t *>::iterator
01271     the_ele_ptr = coarsen_list.begin(),
01272     end_ele_ptr = coarsen_list.end();
01273     for (;the_ele_ptr != end_ele_ptr;++ the_ele_ptr) {
01274       (*the_ele_ptr)->value = 0;
01275     }
01276   } 
01277 
01278   void renumerateElementHSFC(void (*f)(double,double,double,double&,double&,double&)=NULL);
01282   virtual void writeEasyMesh(const std::string&) const;
01286   virtual void writeTecplotData(const std::string&) const;
01290   virtual void writeOpenDXData(const std::string&) const;
01297   virtual void writeSimplestSimplexMesh(const std::string&) const;
01301   virtual void writeSimplexMesh(const std::string&) const;
01302   public:
01303   friend class IrregularMesh<DIM,DOW>;
01304   };
01305 
01313 template <int DIM, int DOW=DIM>
01314   class ElementIterator
01315   {
01316   public:
01317   enum { n_child = HGeometry<DIM,DOW>::n_child };
01318 
01319   typedef HElement<DIM,DOW> value_t;
01320   typedef value_t * pointer_t;
01321   typedef value_t& reference_t;
01322 
01323   typedef typename std::list<pointer_t>::iterator root_t;
01324   typedef ElementIterator<DIM,DOW> this_t;
01325   typedef IrregularMesh<DIM, DOW> ir_mesh_t;
01326 
01327   protected:
01328   ir_mesh_t * mesh;
01329   root_t root_element;
01330   pointer_t element;
01331 
01332   public:
01333   ElementIterator();
01334   ElementIterator(ir_mesh_t *, root_t&, pointer_t);
01335   ElementIterator(const this_t&);
01336   virtual ~ElementIterator();
01337 
01338   public:
01339   this_t& operator=(const this_t&);
01340 
01341   const reference_t operator*() const {return *element;};
01342   reference_t operator*() {return *element;};
01343 
01344   operator const pointer_t() const {return element;};
01345   operator pointer_t() {return element;};
01346 
01347   const pointer_t operator->() const {return element;};
01348   pointer_t operator->() {return element;};
01349 
01350   virtual this_t& operator++() = 0;
01351 
01352   friend bool operator== <>(const this_t&, this_t&);
01353   friend bool operator!= <>(const this_t&, this_t&);
01354 
01355   public:
01356   friend class IrregularMesh<DIM, DOW>;
01357   friend class HElement<DIM, DOW>;
01358   friend class ActiveElementPairIterator<DIM, DOW>;
01359   };
01360 
01361 
01362 
01367 template <int DIM, int DOW=DIM>
01368   class RootFirstElementIterator : public ElementIterator<DIM, DOW>
01369   {
01370   public:
01371   enum { n_child = HGeometry<DIM,DOW>::n_child };
01372 
01373   typedef ElementIterator<DIM,DOW> base_t;
01374   typedef typename base_t::root_t root_t;
01375   typedef RootFirstElementIterator<DIM,DOW> this_t;
01376 
01377   using base_t::mesh;
01378   using base_t::root_element;
01379   using base_t::element;
01380 
01381   public:
01382   RootFirstElementIterator() {};
01383   RootFirstElementIterator(IrregularMesh<DIM, DOW> * m, 
01384                            root_t& i, 
01385                            HElement<DIM, DOW> * e) :
01386   base_t::ElementIterator(m, i, e) {};
01387   public:
01388   virtual this_t& operator++();
01389   };
01390 
01391 
01392 
01397 template <int DIM, int DOW=DIM>
01398   class ActiveElementIterator : public RootFirstElementIterator<DIM, DOW>
01399   {
01400   public:
01401   enum { n_child = HGeometry<DIM,DOW>::n_child };
01402 
01403   typedef RootFirstElementIterator<DIM,DOW> base_t;
01404   typedef typename base_t::root_t root_t;
01405   typedef ActiveElementIterator<DIM,DOW> this_t;
01406   public:
01407   ActiveElementIterator() {};
01408   ActiveElementIterator(IrregularMesh<DIM, DOW> * m, 
01409                         root_t& i,
01410                         HElement<DIM, DOW> * e) : base_t(m, i, e) {};
01411   ActiveElementIterator(const base_t& it) : base_t(it) {}
01412   public:
01413   virtual this_t& operator++();
01414   };
01415 
01416 
01417 
01423 template <int DIM, int DOW=DIM>
01424   class IrregularMeshPair
01425   {
01426   public:
01427   enum { dim = DIM, dow = DOW };
01428 
01429   typedef IrregularMesh<DIM, DOW> ir_mesh_t;
01430   typedef ActiveElementPairIterator<DIM, DOW> iterator;
01431 
01432   ir_mesh_t *   mesh0;
01433   ir_mesh_t *   mesh1;
01434   public:
01435   IrregularMeshPair(ir_mesh_t&, ir_mesh_t&);
01436   IrregularMeshPair(ir_mesh_t *, ir_mesh_t *);
01437   ~IrregularMeshPair();
01438   public:
01439   iterator beginActiveElementPair();
01440   iterator endActiveElementPair();
01441   };
01442 
01443 
01444 
01452 template <int DIM, int DOW=DIM>
01453   class ActiveElementPairIterator
01454   {
01455   public:
01456   typedef IrregularMeshPair<DIM, DOW> ir_mesh_pair_t;
01457 
01458   public:
01459   typedef int State;
01460   static const State GREAT_THAN                 = -1;   // 0 is the ancestor of 1
01461   static const State EQUAL                      = 0;    // equal
01462   static const State LESS_THAN                  = 1;    // 0 is a descendant of 1
01463 
01464   private:
01465   typedef RootFirstElementIterator<DIM,DOW> iterator;
01466   typedef ActiveElementPairIterator<DIM,DOW> this_t;
01467 
01468   public:
01469   typedef HElement<DIM,DOW> value_t;
01470   typedef value_t& reference_t;
01471   typedef value_t * pointer_t;
01472 
01473   ir_mesh_pair_t * mesh_pair;
01474   State                                         st;
01475   iterator              iterator0;
01476   iterator              iterator1;
01477 
01478   public:
01479   ActiveElementPairIterator() : mesh_pair(NULL) {};
01480   ActiveElementPairIterator(ir_mesh_pair_t * mp,
01481                             State s,
01482                             const iterator& it0,
01483                             const iterator& it1) :
01484   mesh_pair(mp), st(s), iterator0(it0), iterator1(it1) {};
01485   ActiveElementPairIterator(const this_t&);
01486   ~ActiveElementPairIterator() {};
01487 
01488   public:
01489   const reference_t operator()(u_int i) const {
01490     if (i == 0) return *iterator0;
01491     else if (i == 1) return *iterator1;
01492     else Assert (false, ExcInternalError()); // something must be wrong
01493   }
01494   reference_t operator()(u_int i) {
01495     if (i == 0) return *iterator0;
01496     else if (i == 1) return *iterator1;
01497     else {
01498       Assert (false, ExcInternalError()); // something must be wrong
01499       return *((HElement<DIM, DOW> *)NULL);
01500     }
01501   };
01502 
01503   const State& state() const {return st;};
01504 
01505   this_t& operator=(const this_t&);
01506   this_t& operator++();
01507 
01508   friend bool operator== <>(const this_t&, this_t&);
01509   friend bool operator!= <>(const this_t&, this_t&);
01510   public:
01511   friend class IrregularMeshPair<DIM, DOW>;
01512   };
01513 
01514 
01515 
01522 template <int DIM, int DOW=DIM>
01523   class Indicator : public std::vector<double>
01524   {
01525   public:
01526   enum { dim = DIM, dow = DOW };
01527   typedef RegularMesh<DIM, DOW> mesh_t;
01528 
01529   public:
01530   mesh_t *                              msh;
01531   public:
01532   Indicator() : msh(NULL) {};
01533   explicit Indicator(mesh_t& m) : msh(&m) {
01534     resize(msh->n_geometry(DIM));
01535   };
01536   ~Indicator() {};
01537   public:
01538   const mesh_t& mesh() const {return *msh;}
01539   void reinit(mesh_t& m, bool is_bare = false) {
01540     msh = &m;
01541     if (! is_bare) {
01542       resize(msh->n_geometry(DIM));
01543       std::fill(begin(), end(), 0.0);
01544     }
01545   }
01546   };
01547 
01548 
01549 
01563 template <int DIM, int DOW=DIM>
01564   class MeshAdaptor
01565   {
01566   public:
01567   enum { dim = DIM, dow = DOW };
01568   typedef IrregularMesh<DIM,DOW> ir_mesh_t;
01569   typedef Indicator<DIM,DOW> indicator_t;
01570 
01571   private:
01572   ir_mesh_t *                   from_mesh;
01573   ir_mesh_t *                   to_mesh;
01574   indicator_t *                 ind;
01575   double                        tol;
01576   double                        convergence_order;
01577   int                           refine_step;
01578   double                        refine_threshold;
01579   double                        coarse_threshold;
01580 
01581   bool                          _is_refine_only;
01582 
01583   public:
01584   MeshAdaptor();
01585   explicit MeshAdaptor(ir_mesh_t& f);
01586   MeshAdaptor(ir_mesh_t& f, ir_mesh_t& t);
01587   ~MeshAdaptor();
01588 
01589   public:
01590   void reinit(ir_mesh_t& f) {from_mesh = &f; to_mesh = &f;};
01591   void reinit(ir_mesh_t& f, ir_mesh_t& t)  {from_mesh = &f; to_mesh = &t;};
01592   const ir_mesh_t& fromMesh() const {return *from_mesh;};
01593   void setFromMesh(ir_mesh_t& f) {from_mesh = &f;};
01594   const ir_mesh_t& toMesh() const {return *to_mesh;};
01595   void  setToMesh(ir_mesh_t& t) {to_mesh = &t;};
01596   const indicator_t& indicator() const {return *ind;};
01597   const double& indicator(const int& i) const {return (*ind)[i];};
01598   double& indicator(const int& i) {return (*ind)[i];};
01599   void setIndicator(indicator_t& i) {
01600     ind = &i;
01601     Assert (&(ind->mesh()) == &(from_mesh->regularMesh()), ExcInternalError());
01602   };
01603   const double& tolerence() const {return tol;};
01604   double& tolerence() {return tol;};
01605   const double& convergenceOrder() const {return convergence_order;};
01606   double& convergenceOrder() {return convergence_order;};
01607   const int& refineStep() const {return refine_step;};
01608   int& refineStep() {return refine_step;};
01609   const double& refineThreshold() const {return refine_threshold;};
01610   double& refineThreshold() {return refine_threshold;};
01611   const double& coarseThreshold() const {return coarse_threshold;};
01612   double& coarseThreshold() {return coarse_threshold;};
01613 
01614   bool is_refine_only() const { return _is_refine_only; }
01615   bool& is_refine_only() { return _is_refine_only; }
01616 
01617   bool is_indicator_underflow(double ind) const {
01618     return (ind < coarse_threshold*tolerence());
01619   }
01620   bool is_indicator_overflow(double ind) const {
01621     double convergence_coefficient = pow(2.0, DIM + convergenceOrder());
01622     return (ind > refine_threshold*convergence_coefficient*tolerence());
01623   }
01624   bool is_indicator_overflow(double ind, double convergence_coefficient) const {
01625     return (ind > refine_threshold*convergence_coefficient*tolerence());
01626   }
01627 
01628   public:
01629   void globalRefine(unsigned int i = 1);
01630   void randomRefine(double percent = 50.0);
01631   void adapt();
01632   private:
01633   void collectIndicator(HElement<DIM,DOW>&, double);
01634   void collectIndicator();
01635   void prepareToMesh();
01636   void implementAdaption();
01637   void adaptElement(HElement<DIM,DOW>&, double, int);
01638   };
01639 
01640 #endif // _HGeometry_h_
01641