1 #ifndef MESHTYPE 2 #define MESHTYPE 3 4 5 /**************************************************************************/ 6 /* File: meshtype.hpp */ 7 /* Author: Joachim Schoeberl */ 8 /* Date: 01. Okt. 95 */ 9 /**************************************************************************/ 10 11 namespace netgen 12 { 13 14 /* 15 Classes for NETGEN 16 */ 17 18 19 20 enum ELEMENT_TYPE { 21 SEGMENT = 1, SEGMENT3 = 2, 22 TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14, 23 TET = 20, TET10 = 21, 24 PYRAMID = 22, PRISM = 23, PRISM12 = 24, 25 HEX = 25 26 }; 27 28 typedef int ELEMENT_EDGE[2]; // initial point, end point 29 typedef int ELEMENT_FACE[4]; // points, last one is -1 for trig 30 31 32 #define ELEMENT_MAXPOINTS 12 33 #define ELEMENT2D_MAXPOINTS 8 34 35 36 enum POINTTYPE { FIXEDPOINT = 1, EDGEPOINT = 2, SURFACEPOINT = 3, INNERPOINT = 4 }; 37 enum ELEMENTTYPE { FREEELEMENT, FIXEDELEMENT }; 38 enum OPTIMIZEGOAL { OPT_QUALITY, OPT_CONFORM, OPT_REST, OPT_WORSTCASE, OPT_LEGAL }; 39 40 41 42 extern DLL_HEADER int GetTimeStamp(); 43 extern DLL_HEADER int NextTimeStamp(); 44 45 class PointGeomInfo 46 { 47 public: 48 int trignum; // for STL Meshing 49 double u, v; // for OCC Meshing 50 PointGeomInfo()51 PointGeomInfo () 52 : trignum(-1), u(0), v(0) { ; } 53 }; 54 operator <<(ostream & ost,const PointGeomInfo & gi)55 inline ostream & operator<< (ostream & ost, const PointGeomInfo & gi) 56 { 57 return (ost << gi.trignum << " " << gi.u << " " << gi.v); 58 } 59 operator >>(istream & ist,PointGeomInfo & gi)60 inline istream & operator>> (istream & ist, PointGeomInfo & gi) 61 { 62 return (ist >> gi.trignum >> gi.u >> gi.v); 63 } 64 65 66 67 #define MULTIPOINTGEOMINFO_MAX 100 68 class MultiPointGeomInfo 69 { 70 int cnt; 71 PointGeomInfo mgi[MULTIPOINTGEOMINFO_MAX]; 72 public: MultiPointGeomInfo()73 MultiPointGeomInfo () { cnt = 0; } 74 int AddPointGeomInfo (const PointGeomInfo & gi); Init()75 void Init () { cnt = 0; } DeleteAll()76 void DeleteAll () { cnt = 0; } 77 GetNPGI() const78 int GetNPGI () const { return cnt; } GetPGI(int i) const79 const PointGeomInfo & GetPGI (int i) const { return mgi[i-1]; } 80 }; 81 82 83 class EdgePointGeomInfo 84 { 85 public: 86 int edgenr; 87 int body; // for ACIS 88 double dist; // for 2d meshing 89 double u, v; // for OCC Meshing 90 91 public: EdgePointGeomInfo()92 EdgePointGeomInfo () 93 : edgenr(0), body(0), dist(0.0), u(0.0), v(0.0) { ; } 94 95 operator =(const EdgePointGeomInfo & gi2)96 EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) 97 { 98 edgenr = gi2.edgenr; 99 body = gi2.body; 100 dist = gi2.dist; 101 u = gi2.u; v = gi2.v; 102 return *this; 103 } 104 }; 105 operator <<(ostream & ost,const EdgePointGeomInfo & gi)106 inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi) 107 { 108 ost << "epgi: edgnr=" << gi.edgenr << ", dist=" << gi.dist; 109 return ost; 110 } 111 112 113 114 115 116 class PointIndex 117 { 118 int i; 119 public: PointIndex()120 PointIndex () { ; } PointIndex(int ai)121 PointIndex (int ai) : i(ai) { ; } operator =(const PointIndex & ai)122 PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } operator =(int ai)123 PointIndex & operator= (int ai) { i = ai; return *this; } operator int() const124 operator int () const { return i; } GetInt() const125 int GetInt () const { return i; } operator ++(int)126 PointIndex operator++ (int) { int hi = i; i++; return PointIndex(hi); } operator --(int)127 PointIndex operator-- (int) { int hi = i; i--; return PointIndex(hi); } 128 129 #ifdef BASE0 130 enum { BASE = 0 }; 131 #else 132 enum { BASE = 1 }; 133 #endif 134 }; 135 operator >>(istream & ist,PointIndex & pi)136 inline istream & operator>> (istream & ist, PointIndex & pi) 137 { 138 int i; ist >> i; pi = i; return ist; 139 } 140 operator <<(ostream & ost,const PointIndex & pi)141 inline ostream & operator<< (ostream & ost, const PointIndex & pi) 142 { 143 return (ost << pi.GetInt()); 144 } 145 146 147 148 149 class ElementIndex 150 { 151 int i; 152 public: ElementIndex()153 ElementIndex () { ; } ElementIndex(int ai)154 ElementIndex (int ai) : i(ai) { ; } operator =(const ElementIndex & ai)155 ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; } operator =(int ai)156 ElementIndex & operator= (int ai) { i = ai; return *this; } operator int() const157 operator int () const { return i; } operator ++(int)158 ElementIndex & operator++ (int) { i++; return *this; } operator --(int)159 ElementIndex & operator-- (int) { i--; return *this; } 160 }; 161 operator >>(istream & ist,ElementIndex & pi)162 inline istream & operator>> (istream & ist, ElementIndex & pi) 163 { 164 int i; ist >> i; pi = i; return ist; 165 } 166 operator <<(ostream & ost,const ElementIndex & pi)167 inline ostream & operator<< (ostream & ost, const ElementIndex & pi) 168 { 169 return (ost << int(pi)); 170 } 171 172 173 class SurfaceElementIndex 174 { 175 int i; 176 public: SurfaceElementIndex()177 SurfaceElementIndex () { ; } SurfaceElementIndex(int ai)178 SurfaceElementIndex (int ai) : i(ai) { ; } operator =(const SurfaceElementIndex & ai)179 SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) 180 { i = ai.i; return *this; } operator =(int ai)181 SurfaceElementIndex & operator= (int ai) { i = ai; return *this; } operator int() const182 operator int () const { return i; } operator ++(int)183 SurfaceElementIndex & operator++ (int) { i++; return *this; } operator --(int)184 SurfaceElementIndex & operator-- (int) { i--; return *this; } 185 }; 186 operator >>(istream & ist,SurfaceElementIndex & pi)187 inline istream & operator>> (istream & ist, SurfaceElementIndex & pi) 188 { 189 int i; ist >> i; pi = i; return ist; 190 } 191 operator <<(ostream & ost,const SurfaceElementIndex & pi)192 inline ostream & operator<< (ostream & ost, const SurfaceElementIndex & pi) 193 { 194 return (ost << int(pi)); 195 } 196 197 class SegmentIndex 198 { 199 int i; 200 public: SegmentIndex()201 SegmentIndex () { ; } SegmentIndex(int ai)202 SegmentIndex (int ai) : i(ai) { ; } operator =(const SegmentIndex & ai)203 SegmentIndex & operator= (const SegmentIndex & ai) 204 { i = ai.i; return *this; } operator =(int ai)205 SegmentIndex & operator= (int ai) { i = ai; return *this; } operator int() const206 operator int () const { return i; } operator ++(int)207 SegmentIndex & operator++ (int) { i++; return *this; } operator --(int)208 SegmentIndex & operator-- (int) { i--; return *this; } 209 }; 210 operator >>(istream & ist,SegmentIndex & pi)211 inline istream & operator>> (istream & ist, SegmentIndex & pi) 212 { 213 int i; ist >> i; pi = i; return ist; 214 } 215 operator <<(ostream & ost,const SegmentIndex & pi)216 inline ostream & operator<< (ostream & ost, const SegmentIndex & pi) 217 { 218 return (ost << int(pi)); 219 } 220 221 222 223 224 /** 225 Point in the mesh. 226 Contains layer (a new feature in 4.3 for overlapping meshes. 227 */ 228 class MeshPoint : public Point<3> 229 { 230 int layer; 231 double singular; // singular factor for hp-refinement 232 POINTTYPE type; 233 234 #ifdef PARALLEL 235 bool isghost; 236 #endif 237 238 public: MeshPoint()239 MeshPoint () 240 { 241 #ifdef PARALLEL 242 isghost = 0; 243 #endif 244 ; 245 } 246 MeshPoint(const Point<3> & ap,int alayer=1,POINTTYPE apt=INNERPOINT)247 MeshPoint (const Point<3> & ap, int alayer = 1, POINTTYPE apt = INNERPOINT) 248 : Point<3> (ap), layer(alayer), singular(0.),type(apt) 249 { 250 #ifdef PARALLEL 251 isghost = 0; 252 #endif 253 ; 254 } 255 SetPoint(const Point<3> & ap)256 void SetPoint (const Point<3> & ap) 257 { 258 Point<3>::operator= (ap); 259 layer = 0; 260 singular = 0; 261 #ifdef PARALLEL 262 isghost = 0; 263 #endif 264 } 265 GetLayer() const266 int GetLayer() const { return layer; } 267 Type() const268 POINTTYPE Type() const { return type; } SetType(POINTTYPE at)269 void SetType(POINTTYPE at) { type = at; } 270 Singularity() const271 double Singularity() const { return singular; } Singularity(double s)272 void Singularity(double s) { singular = s; } IsSingular() const273 bool IsSingular() const { return (singular != 0.0); } 274 275 #ifdef PARALLEL IsGhost() const276 bool IsGhost () const { return isghost; } SetGhost(bool aisghost)277 void SetGhost ( bool aisghost ) { isghost = aisghost; } 278 279 static MPI_Datatype MyGetMPIType ( ); 280 #endif 281 282 }; 283 operator <<(ostream & s,const MeshPoint & pt)284 inline ostream & operator<<(ostream & s, const MeshPoint & pt) 285 { 286 return (s << Point<3> (pt)); 287 } 288 289 290 291 292 typedef Array<MeshPoint, PointIndex::BASE> T_POINTS; 293 294 295 296 /** 297 Triangle element for surface mesh generation. 298 */ 299 class Element2d 300 { 301 /// point numbers 302 PointIndex pnum[ELEMENT2D_MAXPOINTS]; 303 /// geom info of points 304 PointGeomInfo geominfo[ELEMENT2D_MAXPOINTS]; 305 306 /// surface nr 307 int index:16; 308 /// 309 ELEMENT_TYPE typ:6; 310 /// number of points 311 unsigned int np:4; 312 bool badel:1; 313 bool refflag:1; // marked for refinement 314 bool strongrefflag:1; 315 bool deleted:1; // element is deleted 316 317 // Philippose - 08 August 2010 318 // Set a new property for each element, to 319 // control whether it is visible or not 320 bool visible:1; // element visible 321 322 /// order for hp-FEM 323 unsigned int orderx:6; 324 unsigned int ordery:6; 325 326 #ifdef PARALLEL 327 bool isghost; 328 int partitionNumber; 329 #endif 330 331 /// a linked list for all segments in the same face 332 SurfaceElementIndex next; 333 334 public: 335 /// 336 Element2d (); 337 /// 338 Element2d (int anp); 339 /// 340 DLL_HEADER Element2d (ELEMENT_TYPE type); 341 /// 342 Element2d (int pi1, int pi2, int pi3); 343 /// 344 Element2d (int pi1, int pi2, int pi3, int pi4); 345 /// GetType() const346 ELEMENT_TYPE GetType () const { return typ; } 347 /// SetType(ELEMENT_TYPE atyp)348 void SetType (ELEMENT_TYPE atyp) 349 { 350 typ = atyp; 351 switch (typ) 352 { 353 case TRIG: np = 3; break; 354 case QUAD: np = 4; break; 355 case TRIG6: np = 6; break; 356 case QUAD6: np = 6; break; 357 case QUAD8: np = 8; break; 358 default: 359 PrintSysError ("Element2d::SetType, illegal type ", typ); 360 } 361 } 362 /// GetNP() const363 int GetNP() const { return np; } 364 /// GetNV() const365 int GetNV() const 366 { 367 switch (typ) 368 { 369 case TRIG: 370 case TRIG6: return 3; 371 372 case QUAD: 373 case QUAD8: 374 case QUAD6: return 4; 375 default: 376 #ifdef DEBUG 377 PrintSysError ("element2d::GetNV not implemented for typ", typ) 378 #endif 379 ; 380 } 381 return np; 382 } 383 384 /// operator [](int i)385 PointIndex & operator[] (int i) { return pnum[i]; } 386 /// operator [](int i) const387 const PointIndex & operator[] (int i) const { return pnum[i]; } 388 389 /// PNum(int i)390 PointIndex & PNum (int i) { return pnum[i-1]; } 391 /// PNum(int i) const392 const PointIndex & PNum (int i) const { return pnum[i-1]; } 393 /// PNumMod(int i)394 PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } 395 /// PNumMod(int i) const396 const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } 397 /// 398 399 /// GeomInfoPi(int i)400 PointGeomInfo & GeomInfoPi (int i) { return geominfo[i-1]; } 401 /// GeomInfoPi(int i) const402 const PointGeomInfo & GeomInfoPi (int i) const { return geominfo[i-1]; } 403 /// GeomInfoPiMod(int i)404 PointGeomInfo & GeomInfoPiMod (int i) { return geominfo[(i-1) % np]; } 405 /// GeomInfoPiMod(int i) const406 const PointGeomInfo & GeomInfoPiMod (int i) const { return geominfo[(i-1) % np]; } 407 408 SetIndex(int si)409 void SetIndex (int si) { index = si; } 410 /// GetIndex() const411 int GetIndex () const { return index; } 412 GetOrder() const413 int GetOrder () const { return orderx; } SetOrder(int aorder)414 void SetOrder (int aorder) { orderx = ordery = aorder; } 415 416 GetOrder(int & ox,int & oy) const417 void GetOrder (int & ox, int & oy) const { ox = orderx, oy =ordery;}; GetOrder(int & ox,int & oy,int & oz) const418 void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz=0; } SetOrder(int ox,int oy,int)419 void SetOrder (int ox, int oy, int /* oz */) { orderx = ox; ordery = oy;} SetOrder(int ox,int oy)420 void SetOrder (int ox, int oy) { orderx = ox; ordery = oy;} 421 422 423 /// 424 void GetBox (const T_POINTS & points, Box3d & box) const; 425 /// invert orientation 426 inline void Invert (); 427 /// 428 void Invert2 (); 429 /// first point number is smallest 430 inline void NormalizeNumbering (); 431 /// 432 void NormalizeNumbering2 (); 433 BadElement() const434 bool BadElement() const { return badel; } 435 436 // friend ostream & operator<<(ostream & s, const Element2d & el); 437 friend class Mesh; 438 439 440 /// get number of 'integration points' 441 int GetNIP () const; 442 void GetIntegrationPoint (int ip, Point2d & p, double & weight) const; 443 444 void GetTransformation (int ip, const Array<Point2d> & points, 445 class DenseMatrix & trans) const; 446 void GetTransformation (int ip, class DenseMatrix & pmat, 447 class DenseMatrix & trans) const; 448 449 void GetShape (const Point2d & p, class Vector & shape) const; 450 void GetShapeNew (const Point<2> & p, class FlatVector & shape) const; 451 /// matrix 2 * np 452 void GetDShape (const Point2d & p, class DenseMatrix & dshape) const; 453 void GetDShapeNew (const Point<2> & p, class MatrixFixWidth<2> & dshape) const; 454 /// matrix 2 * np 455 void GetPointMatrix (const Array<Point2d> & points, 456 class DenseMatrix & pmat) const; 457 458 void ComputeIntegrationPointData () const; 459 460 461 double CalcJacobianBadness (const Array<Point2d> & points) const; 462 double CalcJacobianBadness (const T_POINTS & points, 463 const Vec<3> & n) const; 464 double CalcJacobianBadnessDirDeriv (const Array<Point2d> & points, 465 int pi, Vec2d & dir, double & dd) const; 466 467 468 Delete()469 void Delete () { deleted = 1; pnum[0] = pnum[1] = pnum[2] = pnum[3] = PointIndex::BASE-1; } IsDeleted() const470 bool IsDeleted () const 471 { 472 #ifdef DEBUG 473 if (pnum[0] < PointIndex::BASE && !deleted) 474 cerr << "Surfelement has illegal pnum, but not marked as deleted" << endl; 475 #endif 476 return deleted; 477 } 478 479 // Philippose - 08 August 2010 480 // Access functions for the new property: visible Visible(bool vis=1)481 void Visible(bool vis = 1) 482 { visible = vis; } IsVisible() const483 bool IsVisible () const 484 { return visible; } 485 SetRefinementFlag(bool rflag=1)486 void SetRefinementFlag (bool rflag = 1) 487 { refflag = rflag; } TestRefinementFlag() const488 bool TestRefinementFlag () const 489 { return refflag; } 490 SetStrongRefinementFlag(bool rflag=1)491 void SetStrongRefinementFlag (bool rflag = 1) 492 { strongrefflag = rflag; } TestStrongRefinementFlag() const493 bool TestStrongRefinementFlag () const 494 { return strongrefflag; } 495 496 NextElement()497 SurfaceElementIndex NextElement() { return next; } 498 499 bool operator==(const Element2d & el2) const; 500 501 int HasFace(const Element2d& el) const; 502 /// 503 int meshdocval; 504 /// 505 int hp_elnr; 506 507 #ifdef PARALLEL IsGhost() const508 bool IsGhost () const { return isghost; } SetGhost(bool aisghost)509 void SetGhost ( bool aisghost ) { isghost = aisghost; } 510 511 // by JS, only for 2D meshes ???? GetPartition() const512 int GetPartition () const { return partitionNumber; } SetPartition(int nr)513 void SetPartition (int nr) { partitionNumber = nr; }; 514 515 #else IsGhost() const516 bool IsGhost () const { return false; } 517 #endif 518 }; 519 520 521 ostream & operator<<(ostream & s, const Element2d & el); 522 523 524 525 526 527 class IntegrationPointData 528 { 529 public: 530 Point<3> p; 531 double weight; 532 Vector shape; 533 DenseMatrix dshape; 534 }; 535 536 537 538 539 540 541 542 543 /** 544 Volume element 545 */ 546 class Element 547 { 548 private: 549 /// point numbers 550 PointIndex pnum[ELEMENT_MAXPOINTS]; 551 /// 552 ELEMENT_TYPE typ:6; 553 /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) 554 int np:5; 555 /// 556 class flagstruct { 557 public: 558 bool marked:1; // marked for refinement 559 bool badel:1; // angles worse then limit 560 bool reverse:1; // for refinement a la Bey 561 bool illegal:1; // illegal, will be split or swaped 562 bool illegal_valid:1; // is illegal-flag valid ? 563 bool badness_valid:1; // is badness valid ? 564 bool refflag:1; // mark element for refinement 565 bool strongrefflag:1; 566 bool deleted:1; // element is deleted, will be removed from array 567 bool fixed:1; // don't change element in optimization 568 }; 569 570 /// sub-domain index 571 short int index; 572 /// order for hp-FEM 573 unsigned int orderx:6; 574 unsigned int ordery:6; 575 unsigned int orderz:6; 576 /* unsigned int levelx:6; 577 unsigned int levely:6; 578 unsigned int levelz:6; */ 579 /// stored shape-badness of element 580 float badness; 581 582 #ifdef PARALLEL 583 /// number of partition for parallel computation 584 int partitionNumber; 585 bool isghost; 586 #endif 587 588 public: 589 flagstruct flags; 590 591 /// 592 Element (); 593 /// 594 Element (int anp); 595 /// 596 Element (ELEMENT_TYPE type); 597 /// 598 Element & operator= (const Element & el2); 599 600 /// 601 void SetNP (int anp); 602 /// 603 void SetType (ELEMENT_TYPE atyp); 604 /// GetNP() const605 int GetNP () const { return np; } 606 /// GetNV() const607 int GetNV() const 608 { 609 switch (typ) 610 { 611 case TET: 612 case TET10: 613 return 4; 614 case PRISM12: 615 case PRISM: 616 return 6; 617 case PYRAMID: 618 return 5; 619 case HEX: 620 return 8; 621 default: 622 #ifdef DEBUG 623 PrintSysError ("Element3d::GetNV not implemented for typ ", typ) 624 #endif 625 ; 626 } 627 return np; 628 } 629 630 bool operator==(const Element & el2) const; 631 632 // old style: NP() const633 int NP () const { return np; } 634 635 /// GetType() const636 ELEMENT_TYPE GetType () const { return typ; } 637 638 /// operator [](int i)639 PointIndex & operator[] (int i) { return pnum[i]; } 640 /// operator [](int i) const641 const PointIndex & operator[] (int i) const { return pnum[i]; } 642 643 /// PNum(int i)644 PointIndex & PNum (int i) { return pnum[i-1]; } 645 /// PNum(int i) const646 const PointIndex & PNum (int i) const { return pnum[i-1]; } 647 /// PNumMod(int i)648 PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } 649 /// PNumMod(int i) const650 const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } 651 652 /// SetIndex(int si)653 void SetIndex (int si) { index = si; } 654 /// GetIndex() const655 int GetIndex () const { return index; } 656 GetOrder() const657 int GetOrder () const { return orderx; } 658 void SetOrder (const int aorder) ; 659 GetOrder(int & ox,int & oy,int & oz) const660 void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz = orderz; } 661 void SetOrder (const int ox, const int oy, const int oz); 662 // void GetLevel (int & ox, int & oy, int & oz) const { ox = levelx; oy = levely; oz = levelz; } 663 // void SetLevel (int ox, int oy, int oz) { levelx = ox; levely = oy; levelz = oz; } 664 665 666 /// 667 void GetBox (const T_POINTS & points, Box3d & box) const; 668 /// Calculates Volume of elemenet 669 double Volume (const T_POINTS & points) const; 670 /// 671 virtual void Print (ostream & ost) const; 672 /// GetNFaces() const673 int GetNFaces () const 674 { 675 switch (typ) 676 { 677 case TET: 678 case TET10: return 4; 679 case PYRAMID: return 5; 680 case PRISM: 681 case PRISM12: return 5; 682 default: 683 #ifdef DEBUG 684 PrintSysError ("element3d::GetNFaces not implemented for typ", typ) 685 #endif 686 ; 687 } 688 return 0; 689 } 690 /// 691 inline void GetFace (int i, Element2d & face) const; 692 /// 693 void GetFace2 (int i, Element2d & face) const; 694 /// 695 void Invert (); 696 697 /// split into 4 node tets 698 void GetTets (Array<Element> & locels) const; 699 /// split into 4 node tets, local point nrs 700 void GetTetsLocal (Array<Element> & locels) const; 701 /// returns coordinates of nodes 702 // void GetNodesLocal (Array<Point<3> > & points) const; 703 void GetNodesLocalNew (Array<Point<3> > & points) const; 704 705 /// split surface into 3 node trigs 706 void GetSurfaceTriangles (Array<Element2d> & surftrigs) const; 707 708 709 /// get number of 'integration points' 710 int GetNIP () const; 711 void GetIntegrationPoint (int ip, Point<3> & p, double & weight) const; 712 713 void GetTransformation (int ip, const T_POINTS & points, 714 class DenseMatrix & trans) const; 715 void GetTransformation (int ip, class DenseMatrix & pmat, 716 class DenseMatrix & trans) const; 717 718 void GetShape (const Point<3> & p, class Vector & shape) const; 719 void GetShapeNew (const Point<3> & p, class FlatVector & shape) const; 720 /// matrix 2 * np 721 void GetDShape (const Point<3> & p, class DenseMatrix & dshape) const; 722 void GetDShapeNew (const Point<3> & p, class MatrixFixWidth<3> & dshape) const; 723 /// matrix 3 * np 724 void GetPointMatrix (const T_POINTS & points, 725 class DenseMatrix & pmat) const; 726 727 void ComputeIntegrationPointData () const; 728 729 730 double CalcJacobianBadness (const T_POINTS & points) const; 731 double CalcJacobianBadnessDirDeriv (const T_POINTS & points, 732 int pi, Vec<3> & dir, double & dd) const; 733 double CalcJacobianBadnessGradient (const T_POINTS & points, 734 int pi, Vec<3> & grad) const; 735 736 /// 737 // friend ostream & operator<<(ostream & s, const Element & el); 738 SetRefinementFlag(bool rflag=1)739 void SetRefinementFlag (bool rflag = 1) 740 { flags.refflag = rflag; } TestRefinementFlag() const741 int TestRefinementFlag () const 742 { return flags.refflag; } 743 SetStrongRefinementFlag(bool rflag=1)744 void SetStrongRefinementFlag (bool rflag = 1) 745 { flags.strongrefflag = rflag; } TestStrongRefinementFlag() const746 int TestStrongRefinementFlag () const 747 { return flags.strongrefflag; } 748 Illegal() const749 int Illegal () const 750 { return flags.illegal; } IllegalValid() const751 int IllegalValid () const 752 { return flags.illegal_valid; } SetIllegal(int aillegal)753 void SetIllegal (int aillegal) 754 { 755 flags.illegal = aillegal ? 1 : 0; 756 flags.illegal_valid = 1; 757 } SetLegal(int alegal)758 void SetLegal (int alegal) 759 { 760 flags.illegal = alegal ? 0 : 1; 761 flags.illegal_valid = 1; 762 } 763 Delete()764 void Delete () { flags.deleted = 1; } IsDeleted() const765 bool IsDeleted () const 766 { 767 #ifdef DEBUG 768 if (pnum[0] < PointIndex::BASE && !flags.deleted) 769 cerr << "Volelement has illegal pnum, but not marked as deleted" << endl; 770 #endif 771 772 return flags.deleted; 773 } 774 775 776 777 #ifdef PARALLEL GetPartition() const778 int GetPartition () const { return partitionNumber; } SetPartition(int nr)779 void SetPartition (int nr) { partitionNumber = nr; }; 780 IsGhost() const781 bool IsGhost () const { return isghost; } SetGhost(const bool aisghost)782 void SetGhost ( const bool aisghost ) { isghost = aisghost; } 783 #else IsGhost() const784 bool IsGhost () const { return false; } GetPartition() const785 int GetPartition () const { return 0; } 786 #endif 787 788 int hp_elnr; 789 }; 790 791 ostream & operator<<(ostream & s, const Element & el); 792 793 794 795 796 797 798 /** 799 Edge segment. 800 */ 801 class Segment 802 { 803 public: 804 /// 805 DLL_HEADER Segment(); 806 DLL_HEADER Segment (const Segment& other); 807 ~Segment()808 ~Segment() 809 { ; } 810 811 // friend ostream & operator<<(ostream & s, const Segment & seg); 812 813 PointIndex pnums[3]; // p1, p2, pmid 814 815 int edgenr; 816 /// 817 double singedge_left; 818 double singedge_right; 819 820 /// 0.. not first segment of segs, 1..first of class, 2..first of class, inverse 821 unsigned int seginfo:2; 822 823 /// surface decoding index 824 int si; 825 /// domain number inner side 826 int domin; 827 /// domain number outer side 828 int domout; 829 /// top-level object number of surface 830 int tlosurf; 831 /// 832 PointGeomInfo geominfo[2]; 833 834 /// surfaces describing edge 835 int surfnr1, surfnr2; 836 /// 837 EdgePointGeomInfo epgeominfo[2]; 838 /// 839 // int pmid; // for second order 840 /// 841 int meshdocval; 842 843 private: 844 string* bcname; 845 846 public: 847 /* 848 PointIndex operator[] (int i) const 849 { return (i == 0) ? p1 : p2; } 850 851 PointIndex & operator[] (int i) 852 { return (i == 0) ? p1 : p2; } 853 */ 854 855 Segment& operator=(const Segment & other); 856 857 858 int hp_elnr; 859 SetBCName(string * abcname)860 void SetBCName ( string * abcname ) 861 { 862 bcname = abcname; 863 } 864 BCNamePtr()865 string * BCNamePtr () 866 { return bcname; } 867 BCNamePtr() const868 const string * BCNamePtr () const 869 { return bcname; } 870 GetBCName() const871 string GetBCName () const 872 { 873 if (! bcname ) 874 { 875 return "default"; 876 } 877 return *bcname; 878 } 879 GetNP() const880 int GetNP() const 881 { 882 return (pnums[2] < 0) ? 2 : 3; 883 } 884 GetType() const885 ELEMENT_TYPE GetType() const 886 { 887 return (pnums[2] < 0) ? SEGMENT : SEGMENT3; 888 } 889 operator [](int i)890 PointIndex & operator[] (int i) { return pnums[i]; } operator [](int i) const891 const PointIndex & operator[] (int i) const { return pnums[i]; } 892 }; 893 894 ostream & operator<<(ostream & s, const Segment & seg); 895 896 897 898 899 // class Surface; 900 // class FaceDescriptor; 901 902 /// 903 class FaceDescriptor 904 { 905 /// which surface, 0 if not available 906 int surfnr; 907 /// domain nr inside 908 int domin; 909 /// domain nr outside 910 int domout; 911 /// top level object number of surface 912 int tlosurf; 913 /// boundary condition property 914 int bcprop; 915 // Philippose - 06/07/2009 916 // Add capability to store surface colours along with 917 // other face data 918 /// surface colour (Default: R=0.0 ; G=1.0 ; B=0.0) 919 Vec3d surfcolour; 920 921 /// 922 string * bcname; 923 /// root of linked list 924 SurfaceElementIndex firstelement; 925 926 double domin_singular; 927 double domout_singular; 928 929 public: 930 DLL_HEADER FaceDescriptor(); 931 DLL_HEADER FaceDescriptor(int surfnri, int domini, int domouti, int tlosurfi); 932 DLL_HEADER FaceDescriptor(const Segment & seg); 933 DLL_HEADER FaceDescriptor(const FaceDescriptor& other); ~FaceDescriptor()934 DLL_HEADER ~FaceDescriptor() { ; } 935 936 DLL_HEADER int SegmentFits (const Segment & seg); 937 SurfNr() const938 int SurfNr () const { return surfnr; } DomainIn() const939 int DomainIn () const { return domin; } DomainOut() const940 int DomainOut () const { return domout; } TLOSurface() const941 int TLOSurface () const { return tlosurf; } BCProperty() const942 int BCProperty () const { return bcprop; } 943 944 DomainInSingular() const945 double DomainInSingular() const { return domin_singular; } DomainOutSingular() const946 double DomainOutSingular() const { return domout_singular; } 947 948 // Philippose - 06/07/2009 949 // Get Surface colour SurfColour() const950 Vec3d SurfColour () const { return surfcolour; } 951 string GetBCName () const; 952 // string * BCNamePtr () { return bcname; } 953 // const string * BCNamePtr () const { return bcname; } SetSurfNr(int sn)954 void SetSurfNr (int sn) { surfnr = sn; } SetDomainIn(int di)955 void SetDomainIn (int di) { domin = di; } SetDomainOut(int dom)956 void SetDomainOut (int dom) { domout = dom; } SetBCProperty(int bc)957 void SetBCProperty (int bc) { bcprop = bc; } SetBCName(string * bcn)958 void SetBCName (string * bcn) { bcname = bcn; } 959 // Philippose - 06/07/2009 960 // Set the surface colour SetSurfColour(Vec3d colour)961 void SetSurfColour (Vec3d colour) { surfcolour = colour; } 962 SetDomainInSingular(double v)963 void SetDomainInSingular (double v) { domin_singular = v; } SetDomainOutSingular(double v)964 void SetDomainOutSingular (double v) { domout_singular = v; } 965 FirstElement()966 SurfaceElementIndex FirstElement() { return firstelement; } 967 // friend ostream & operator<<(ostream & s, const FaceDescriptor & fd); 968 friend class Mesh; 969 }; 970 971 ostream & operator<< (ostream & s, const FaceDescriptor & fd); 972 973 974 class EdgeDescriptor 975 { 976 int tlosurf; 977 int surfnr[2]; 978 public: EdgeDescriptor()979 EdgeDescriptor () 980 : tlosurf(-1) 981 { surfnr[0] = surfnr[1] = -1; } 982 SurfNr(int i) const983 int SurfNr (int i) const { return surfnr[i]; } SetSurfNr(int i,int nr)984 void SetSurfNr (int i, int nr) { surfnr[i] = nr; } 985 TLOSurface() const986 int TLOSurface() const { return tlosurf; } SetTLOSurface(int nr)987 void SetTLOSurface (int nr) { tlosurf = nr; } 988 }; 989 990 991 992 class DLL_HEADER MeshingParameters 993 { 994 public: 995 /** 996 3d optimization strategy: 997 // m .. move nodes 998 // M .. move nodes, cheap functional 999 // s .. swap faces 1000 // c .. combine elements 1001 // d .. divide elements 1002 // p .. plot, no pause 1003 // P .. plot, Pause 1004 // h .. Histogramm, no pause 1005 // H .. Histogramm, pause 1006 */ 1007 const char * optimize3d; 1008 /// number of 3d optimization steps 1009 int optsteps3d; 1010 /** 1011 2d optimization strategy: 1012 // s .. swap, opt 6 lines/node 1013 // S .. swap, optimal elements 1014 // m .. move nodes 1015 // p .. plot, no pause 1016 // P .. plot, pause 1017 // c .. combine 1018 **/ 1019 const char * optimize2d; 1020 /// number of 2d optimization steps 1021 int optsteps2d; 1022 /// power of error (to approximate max err optimization) 1023 double opterrpow; 1024 /// do block filling ? 1025 int blockfill; 1026 /// block filling up to distance 1027 double filldist; 1028 /// radius of local environment (times h) 1029 double safety; 1030 /// radius of active environment (times h) 1031 double relinnersafety; 1032 /// use local h ? 1033 int uselocalh; 1034 /// grading for local h 1035 double grading; 1036 /// use delaunay meshing 1037 int delaunay; 1038 /// maximal mesh size 1039 double maxh; 1040 /// minimal mesh size 1041 double minh; 1042 /// file for meshsize 1043 const char * meshsizefilename; 1044 /// start surfacemeshing from everywhere in surface 1045 int startinsurface; 1046 /// check overlapping surfaces (debug) 1047 int checkoverlap; 1048 /// check overlapping surface mesh before volume meshing 1049 int checkoverlappingboundary; 1050 /// check chart boundary (sometimes too restrictive) 1051 int checkchartboundary; 1052 /// safty factor for curvatures (elemetns per radius) 1053 double curvaturesafety; 1054 /// minimal number of segments per edge 1055 double segmentsperedge; 1056 /// use parallel threads 1057 int parthread; 1058 /// weight of element size w.r.t element shape 1059 double elsizeweight; 1060 /// init with default values 1061 1062 1063 /// from mp3: 1064 /// give up quality class, 2d meshing 1065 int giveuptol2d; 1066 /// give up quality class, 3d meshing 1067 int giveuptol; 1068 /// maximal outer steps 1069 int maxoutersteps; 1070 /// class starting star-shape filling 1071 int starshapeclass; 1072 /// if non-zero, baseelement must have baseelnp points 1073 int baseelnp; 1074 /// quality tolerances are handled less careful 1075 int sloppy; 1076 1077 /// limit for max element angle (150-180) 1078 double badellimit; 1079 1080 bool check_impossible; 1081 1082 /// 1083 int secondorder; 1084 /// high order element curvature 1085 int elementorder; 1086 /// quad-dominated surface meshing 1087 int quad; 1088 /// 1089 int inverttets; 1090 /// 1091 int inverttrigs; 1092 /// 1093 int autozrefine; 1094 /// 1095 MeshingParameters (); 1096 /// 1097 void Print (ostream & ost) const; 1098 1099 void CopyFrom(const MeshingParameters & other); 1100 }; 1101 1102 1103 1104 class DebugParameters 1105 { 1106 public: 1107 /// 1108 int debugoutput; 1109 /// use slow checks 1110 int slowchecks; 1111 /// 1112 int haltsuccess; 1113 /// 1114 int haltnosuccess; 1115 /// 1116 int haltlargequalclass; 1117 /// 1118 int haltsegment; 1119 /// 1120 int haltnode; 1121 /// 1122 int haltsegmentp1; 1123 /// 1124 int haltsegmentp2; 1125 /// 1126 int haltexistingline; 1127 /// 1128 int haltoverlap; 1129 /// 1130 int haltface; 1131 /// 1132 int haltfacenr; 1133 /// 1134 DebugParameters (); 1135 }; 1136 1137 1138 1139 Invert()1140 inline void Element2d :: Invert() 1141 { 1142 if (typ == TRIG) 1143 Swap (PNum(2), PNum(3)); 1144 else 1145 Invert2(); 1146 } 1147 1148 1149 1150 NormalizeNumbering()1151 inline void Element2d :: NormalizeNumbering () 1152 { 1153 if (GetNP() == 3) 1154 { 1155 if (PNum(1) < PNum(2) && PNum(1) < PNum(3)) 1156 return; 1157 else 1158 { 1159 if (PNum(2) < PNum(3)) 1160 { 1161 PointIndex pi1 = PNum(2); 1162 PNum(2) = PNum(3); 1163 PNum(3) = PNum(1); 1164 PNum(1) = pi1; 1165 } 1166 else 1167 { 1168 PointIndex pi1 = PNum(3); 1169 PNum(3) = PNum(2); 1170 PNum(2) = PNum(1); 1171 PNum(1) = pi1; 1172 } 1173 } 1174 } 1175 else 1176 NormalizeNumbering2(); 1177 } 1178 1179 1180 1181 static const int gftetfacesa[4][3] = 1182 { { 1, 2, 3 }, 1183 { 2, 0, 3 }, 1184 { 0, 1, 3 }, 1185 { 1, 0, 2 } }; 1186 GetFace(int i,Element2d & face) const1187 inline void Element :: GetFace (int i, Element2d & face) const 1188 { 1189 if (typ == TET) 1190 { 1191 face.SetType(TRIG); 1192 face[0] = pnum[gftetfacesa[i-1][0]]; 1193 face[1] = pnum[gftetfacesa[i-1][1]]; 1194 face[2] = pnum[gftetfacesa[i-1][2]]; 1195 } 1196 else 1197 GetFace2 (i, face); 1198 } 1199 1200 1201 1202 1203 1204 1205 1206 /** 1207 Identification of periodic surfaces, close surfaces, etc. 1208 */ 1209 class Identifications 1210 { 1211 public: 1212 enum ID_TYPE { UNDEFINED = 1, PERIODIC = 2, CLOSESURFACES = 3, CLOSEEDGES = 4}; 1213 1214 1215 private: 1216 class Mesh & mesh; 1217 1218 /// identify points (thin layers, periodic b.c.) 1219 INDEX_2_HASHTABLE<int> * identifiedpoints; 1220 1221 /// the same, with info about the id-nr 1222 INDEX_3_HASHTABLE<int> * identifiedpoints_nr; 1223 1224 /// sorted by identification nr 1225 TABLE<INDEX_2> idpoints_table; 1226 1227 Array<ID_TYPE> type; 1228 1229 /// number of identifications (or, actually used identifications ?) 1230 int maxidentnr; 1231 1232 public: 1233 /// 1234 DLL_HEADER Identifications (class Mesh & amesh); 1235 /// 1236 DLL_HEADER ~Identifications (); 1237 1238 DLL_HEADER void Delete (); 1239 1240 /* 1241 Identify points pi1 and pi2, due to 1242 identification nr identnr 1243 */ 1244 DLL_HEADER void Add (PointIndex pi1, PointIndex pi2, int identnr); 1245 1246 1247 int Get (PointIndex pi1, PointIndex pi2) const; 1248 int GetSymmetric (PointIndex pi1, PointIndex pi2) const; 1249 1250 bool Get (PointIndex pi1, PointIndex pi2, int identnr) const; 1251 bool GetSymmetric (PointIndex pi1, PointIndex pi2, int identnr) const; 1252 1253 /// GetIdentifiedPoints()1254 INDEX_2_HASHTABLE<int> & GetIdentifiedPoints () 1255 { 1256 return *identifiedpoints; 1257 } 1258 Used(PointIndex pi1,PointIndex pi2)1259 bool Used (PointIndex pi1, PointIndex pi2) 1260 { 1261 return identifiedpoints->Used (INDEX_2 (pi1, pi2)); 1262 } 1263 UsedSymmetric(PointIndex pi1,PointIndex pi2)1264 bool UsedSymmetric (PointIndex pi1, PointIndex pi2) 1265 { 1266 return 1267 identifiedpoints->Used (INDEX_2 (pi1, pi2)) || 1268 identifiedpoints->Used (INDEX_2 (pi2, pi1)); 1269 } 1270 1271 /// 1272 void GetMap (int identnr, Array<int,PointIndex::BASE> & identmap, bool symmetric = false) const; 1273 /// GetType(int identnr) const1274 ID_TYPE GetType(int identnr) const 1275 { 1276 if(identnr <= type.Size()) 1277 return type[identnr-1]; 1278 else 1279 return UNDEFINED; 1280 } SetType(int identnr,ID_TYPE t)1281 void SetType(int identnr, ID_TYPE t) 1282 { 1283 while(type.Size() < identnr) 1284 type.Append(UNDEFINED); 1285 type[identnr-1] = t; 1286 } 1287 1288 /// 1289 void GetPairs (int identnr, Array<INDEX_2> & identpairs) const; 1290 /// GetMaxNr() const1291 int GetMaxNr () const { return maxidentnr; } 1292 1293 /// remove secondorder 1294 void SetMaxPointNr (int maxpnum); 1295 1296 void Print (ostream & ost) const; 1297 }; 1298 1299 1300 } 1301 1302 1303 1304 1305 #endif 1306 1307