AFEPack
|
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