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 : unsigned char { 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, PRISM15 = 27, PYRAMID13 = 28, 25 HEX = 25, HEX20 = 26 26 }; 27 28 /* 29 typedef int ELEMENT_EDGE[2]; // initial point, end point 30 typedef int ELEMENT_FACE[4]; // points, last one is -1 for trig 31 */ 32 33 struct ELEMENT_EDGE 34 { 35 int vals[2]; operator []netgen::ELEMENT_EDGE36 int & operator[] (size_t i) { return vals[i]; } operator []netgen::ELEMENT_EDGE37 int operator[] (size_t i) const { return vals[i]; } 38 }; 39 40 struct ELEMENT_FACE 41 { 42 int vals[4]; operator []netgen::ELEMENT_FACE43 int & operator[] (size_t i) { return vals[i]; } operator []netgen::ELEMENT_FACE44 int operator[] (size_t i) const { return vals[i]; } 45 }; 46 47 48 #define ELEMENT_MAXPOINTS 20 49 #define ELEMENT2D_MAXPOINTS 8 50 51 52 enum POINTTYPE : unsigned char { FIXEDPOINT = 1, EDGEPOINT = 2, SURFACEPOINT = 3, INNERPOINT = 4 }; 53 enum ELEMENTTYPE { FREEELEMENT, FIXEDELEMENT }; 54 enum OPTIMIZEGOAL { OPT_QUALITY, OPT_CONFORM, OPT_REST, OPT_WORSTCASE, OPT_LEGAL }; 55 56 57 extern DLL_HEADER size_t timestamp; GetTimeStamp()58 inline size_t GetTimeStamp() 59 { 60 return timestamp; 61 } 62 NextTimeStamp()63 inline size_t NextTimeStamp() 64 { 65 timestamp++; 66 return timestamp; 67 } 68 69 /* 70 extern DLL_HEADER int GetTimeStamp(); 71 extern DLL_HEADER int NextTimeStamp(); 72 */ 73 class PointGeomInfo 74 { 75 public: 76 int trignum; // for STL Meshing 77 double u, v; // for OCC Meshing 78 79 PointGeomInfo () = default; 80 // : trignum(-1), u(0), v(0) { ; } 81 PointGeomInfo (const PointGeomInfo&) = default; 82 PointGeomInfo (PointGeomInfo &&) = default; 83 PointGeomInfo & operator= (const PointGeomInfo&) = default; 84 PointGeomInfo & operator= (PointGeomInfo&&) = default; 85 86 }; 87 operator <<(ostream & ost,const PointGeomInfo & gi)88 inline ostream & operator<< (ostream & ost, const PointGeomInfo & gi) 89 { 90 return (ost << gi.trignum << " " << gi.u << " " << gi.v); 91 } 92 operator >>(istream & ist,PointGeomInfo & gi)93 inline istream & operator>> (istream & ist, PointGeomInfo & gi) 94 { 95 return (ist >> gi.trignum >> gi.u >> gi.v); 96 } 97 98 99 100 class MultiPointGeomInfo 101 { 102 ArrayMem<PointGeomInfo, 100> mgi; 103 public: 104 int AddPointGeomInfo (const PointGeomInfo & gi); Init()105 void Init () { mgi.SetSize(0); } DeleteAll()106 void DeleteAll () { mgi.SetSize(0); } 107 GetNPGI() const108 int GetNPGI () const { return mgi.Size(); } GetPGI(int i) const109 const PointGeomInfo & GetPGI (int i) const { return mgi[i-1]; } 110 111 MultiPointGeomInfo () = default; 112 MultiPointGeomInfo (const MultiPointGeomInfo&) = default; 113 MultiPointGeomInfo (MultiPointGeomInfo &&) = default; 114 MultiPointGeomInfo & operator= (const MultiPointGeomInfo&) = delete; 115 MultiPointGeomInfo & operator= (MultiPointGeomInfo&&) = default; 116 }; 117 118 119 class EdgePointGeomInfo 120 { 121 public: 122 int edgenr; 123 int body; // for ACIS 124 double dist; // for 2d meshing 125 double u, v; // for OCC Meshing 126 127 public: EdgePointGeomInfo()128 EdgePointGeomInfo () 129 : edgenr(0), body(0), dist(0.0), u(0.0), v(0.0) { ; } 130 131 operator =(const EdgePointGeomInfo & gi2)132 EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) 133 { 134 edgenr = gi2.edgenr; 135 body = gi2.body; 136 dist = gi2.dist; 137 u = gi2.u; v = gi2.v; 138 return *this; 139 } 140 }; 141 operator <<(ostream & ost,const EdgePointGeomInfo & gi)142 inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi) 143 { 144 ost << "epgi: edgnr=" << gi.edgenr << ", dist=" << gi.dist; 145 return ost; 146 } 147 148 149 150 151 152 class PointIndex 153 { 154 int i; 155 public: 156 class t_invalid { public: constexpr t_invalid() = default; }; 157 static constexpr t_invalid INVALID{}; 158 159 PointIndex () = default; 160 PointIndex (const PointIndex&) = default; 161 PointIndex (PointIndex &&) = default; 162 PointIndex & operator= (const PointIndex&) = default; 163 PointIndex & operator= (PointIndex&&) = default; 164 PointIndex(int ai)165 constexpr PointIndex (int ai) : i(ai) 166 { 167 #ifdef DEBUG 168 if (ai < PointIndex::BASE) 169 cout << "illegal PointIndex, use PointIndex::INVALID instead" << endl; 170 // throw Exception("illegal PointIndex, use PointIndex::INVALID instead"); 171 #endif 172 } PointIndex(t_invalid inv)173 constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; } 174 // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } operator int() const175 constexpr operator int () const { return i; } operator ++(int)176 PointIndex operator++ (int) { PointIndex hi(*this); i++; return hi; } operator --(int)177 PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; } operator ++()178 PointIndex & operator++ () { i++; return *this; } operator --()179 PointIndex operator-- () { i--; return *this; } operator +=(int add)180 PointIndex operator+= (int add) { i += add; return *this; } Invalidate()181 void Invalidate() { i = PointIndex::BASE-1; } IsValid() const182 bool IsValid() const { return i != PointIndex::BASE-1; } 183 #ifdef BASE0 184 static constexpr size_t BASE = 0; 185 #else 186 static constexpr size_t BASE = 1; 187 #endif 188 DoArchive(Archive & ar)189 void DoArchive (Archive & ar) { ar & i; } 190 }; 191 192 } 193 194 namespace ngcore 195 { 196 template<> IndexBASE()197 constexpr netgen::PointIndex IndexBASE<netgen::PointIndex> () { return netgen::PointIndex(netgen::PointIndex::BASE); } 198 } 199 200 namespace netgen 201 { 202 203 operator >>(istream & ist,PointIndex & pi)204 inline istream & operator>> (istream & ist, PointIndex & pi) 205 { 206 int i; ist >> i; pi = PointIndex(i); return ist; 207 } 208 operator <<(ostream & ost,const PointIndex & pi)209 inline ostream & operator<< (ostream & ost, const PointIndex & pi) 210 { 211 return (ost << int(pi)); 212 } 213 214 template <int N> class PointIndices; 215 template <> class PointIndices<2> : public INDEX_2 216 { 217 public: 218 PointIndices () = default; PointIndices(INDEX_2 i2)219 PointIndices (INDEX_2 i2) : INDEX_2(i2) { ; } PointIndices(PointIndex i1,PointIndex i2)220 PointIndices (PointIndex i1, PointIndex i2) : INDEX_2(i1,i2) { ; } operator [](int i) const221 PointIndex operator[] (int i) const { return PointIndex(INDEX_2::operator[](i)); } operator [](int i)222 PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_2::operator[](i)); } Sort(PointIndex i1,PointIndex i2)223 static PointIndices Sort(PointIndex i1, PointIndex i2) { return INDEX_2::Sort(i1, i2); } 224 template <size_t J> get() const225 PointIndex get() const { return PointIndex(INDEX_2::operator[](J)); } 226 }; 227 } 228 229 namespace std 230 { 231 // structured binding support 232 template <auto N> 233 struct tuple_size<netgen::PointIndices<N>> : std::integral_constant<std::size_t, N> {}; 234 template<size_t N, auto M> struct tuple_element<N,netgen::PointIndices<M>> { using type = netgen::PointIndex; }; 235 } 236 237 namespace netgen 238 { 239 240 class ElementIndex 241 { 242 int i; 243 public: 244 ElementIndex () = default; ElementIndex(int ai)245 constexpr ElementIndex (int ai) : i(ai) { ; } operator =(const ElementIndex & ai)246 ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; } operator =(int ai)247 ElementIndex & operator= (int ai) { i = ai; return *this; } operator int() const248 operator int () const { return i; } operator ++(int)249 ElementIndex operator++ (int) { return ElementIndex(i++); } operator --(int)250 ElementIndex operator-- (int) { return ElementIndex(i--); } operator ++()251 ElementIndex & operator++ () { ++i; return *this; } operator --()252 ElementIndex & operator-- () { --i; return *this; } 253 }; 254 operator >>(istream & ist,ElementIndex & pi)255 inline istream & operator>> (istream & ist, ElementIndex & pi) 256 { 257 int i; ist >> i; pi = i; return ist; 258 } 259 operator <<(ostream & ost,const ElementIndex & pi)260 inline ostream & operator<< (ostream & ost, const ElementIndex & pi) 261 { 262 return (ost << int(pi)); 263 } 264 265 266 class SurfaceElementIndex 267 { 268 int i; 269 public: 270 SurfaceElementIndex () = default; SurfaceElementIndex(int ai)271 constexpr SurfaceElementIndex (int ai) : i(ai) { ; } 272 /* 273 SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) 274 { i = ai.i; return *this; } 275 */ 276 SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) = default; operator =(int ai)277 SurfaceElementIndex & operator= (int ai) { i = ai; return *this; } operator int() const278 constexpr operator int () const { return i; } operator ++(int)279 SurfaceElementIndex operator++ (int) { SurfaceElementIndex hi(*this); i++; return hi; } operator --(int)280 SurfaceElementIndex operator-- (int) { SurfaceElementIndex hi(*this); i--; return hi; } operator ++()281 SurfaceElementIndex & operator++ () { ++i; return *this; } operator --()282 SurfaceElementIndex & operator-- () { --i; return *this; } operator +=(int inc)283 SurfaceElementIndex & operator+= (int inc) { i+=inc; return *this; } DoArchive(Archive & ar)284 void DoArchive (Archive & ar) { ar & i; } 285 }; 286 SetInvalid(SurfaceElementIndex & id)287 inline void SetInvalid (SurfaceElementIndex & id) { id = -1; } IsInvalid(SurfaceElementIndex & id)288 inline bool IsInvalid (SurfaceElementIndex & id) { return id == -1; } 289 operator >>(istream & ist,SurfaceElementIndex & pi)290 inline istream & operator>> (istream & ist, SurfaceElementIndex & pi) 291 { 292 int i; ist >> i; pi = i; return ist; 293 } 294 operator <<(ostream & ost,const SurfaceElementIndex & pi)295 inline ostream & operator<< (ostream & ost, const SurfaceElementIndex & pi) 296 { 297 return (ost << int(pi)); 298 } 299 300 class SegmentIndex 301 { 302 int i; 303 public: 304 SegmentIndex () = default; SegmentIndex(int ai)305 constexpr SegmentIndex (int ai) : i(ai) { ; } operator =(const SegmentIndex & ai)306 SegmentIndex & operator= (const SegmentIndex & ai) 307 { i = ai.i; return *this; } operator =(int ai)308 SegmentIndex & operator= (int ai) { i = ai; return *this; } operator int() const309 constexpr operator int () const { return i; } operator ++()310 SegmentIndex& operator++ () { ++i; return *this; } operator --()311 SegmentIndex& operator-- () { --i; return *this; } operator ++(int)312 SegmentIndex operator++ (int) { return i++; } operator --(int)313 SegmentIndex operator-- (int) { return i--; } 314 }; 315 SetInvalid(SegmentIndex & id)316 inline void SetInvalid (SegmentIndex & id) { id = -1; } IsInvalid(SegmentIndex & id)317 inline bool IsInvalid (SegmentIndex & id) { return id == -1; } 318 319 operator >>(istream & ist,SegmentIndex & pi)320 inline istream & operator>> (istream & ist, SegmentIndex & pi) 321 { 322 int i; ist >> i; pi = i; return ist; 323 } 324 operator <<(ostream & ost,const SegmentIndex & pi)325 inline ostream & operator<< (ostream & ost, const SegmentIndex & pi) 326 { 327 return (ost << int(pi)); 328 } 329 330 331 332 333 /** 334 Point in the mesh. 335 Contains layer (a new feature in 4.3 for overlapping meshes. 336 */ 337 class MeshPoint : public Point<3> 338 { 339 double singular; // singular factor for hp-refinement 340 int layer; 341 POINTTYPE type; 342 343 344 public: MeshPoint()345 MeshPoint () 346 { 347 ; 348 } 349 MeshPoint(const Point<3> & ap,int alayer=1,POINTTYPE apt=INNERPOINT)350 MeshPoint (const Point<3> & ap, int alayer = 1, POINTTYPE apt = INNERPOINT) 351 : Point<3> (ap), layer(alayer), singular(0.),type(apt) 352 { 353 ; 354 } 355 SetPoint(const Point<3> & ap)356 void SetPoint (const Point<3> & ap) 357 { 358 Point<3>::operator= (ap); 359 layer = 0; 360 singular = 0; 361 } 362 Scale(double factor)363 void Scale(double factor) { *testout << "before: " << x[0] << endl; x[0] *= factor; x[1] *= factor; x[2] *= factor; *testout << "after: " << x[0] << endl;} 364 GetLayer() const365 int GetLayer() const { return layer; } 366 Type() const367 POINTTYPE Type() const { return type; } SetType(POINTTYPE at)368 void SetType(POINTTYPE at) { type = at; } 369 Singularity() const370 double Singularity() const { return singular; } Singularity(double s)371 void Singularity(double s) { singular = s; } IsSingular() const372 bool IsSingular() const { return (singular != 0.0); } 373 374 #ifdef PARALLEL 375 static MPI_Datatype MyGetMPIType ( ); 376 #endif 377 DoArchive(Archive & ar)378 void DoArchive (Archive & ar) 379 { 380 // ar & x[0] & x[1] & x[2] & layer & singular; 381 // ar.Do(&x[0], 3); 382 // ar & layer & singular; 383 // ar & (unsigned char&)(type); 384 ar.DoPacked (x[0], x[1], x[2], layer, singular, (unsigned char&)(type)); 385 } 386 }; 387 operator <<(ostream & s,const MeshPoint & pt)388 inline ostream & operator<<(ostream & s, const MeshPoint & pt) 389 { 390 return (s << Point<3> (pt)); 391 } 392 393 394 395 396 // typedef NgArray<MeshPoint, PointIndex::BASE, PointIndex> T_POINTS; 397 typedef Array<MeshPoint, PointIndex> T_POINTS; 398 399 400 401 /** 402 Triangle element for surface mesh generation. 403 */ 404 class Element2d 405 { 406 /// point numbers 407 PointIndex pnum[ELEMENT2D_MAXPOINTS]; 408 /// geom info of points 409 PointGeomInfo geominfo[ELEMENT2D_MAXPOINTS]; 410 411 /// surface nr 412 int index; 413 /// 414 ELEMENT_TYPE typ; 415 /// number of points 416 int8_t np; 417 bool badel:1; 418 bool refflag:1; // marked for refinement 419 bool strongrefflag:1; 420 bool deleted:1; // element is deleted 421 422 // Philippose - 08 August 2010 423 // Set a new property for each element, to 424 // control whether it is visible or not 425 bool visible:1; // element visible 426 bool is_curved:1; // element is (high order) curved 427 /// order for hp-FEM 428 unsigned int orderx:6; 429 unsigned int ordery:6; 430 431 // #ifdef PARALLEL 432 // int partitionNumber; 433 // #endif 434 435 /// a linked list for all segments in the same face 436 SurfaceElementIndex next; 437 438 public: 439 /// 440 Element2d (); 441 Element2d (const Element2d &) = default; 442 Element2d (Element2d &&) = default; 443 Element2d & operator= (const Element2d &) = default; 444 Element2d & operator= (Element2d &&) = default; operator =(initializer_list<PointIndex> list)445 Element2d & operator= (initializer_list<PointIndex> list) 446 { 447 size_t cnt = 0; 448 for (auto val : list) 449 pnum[cnt++] = val; 450 return *this; 451 } operator =(initializer_list<std::tuple<PointIndex,PointGeomInfo>> list)452 Element2d & operator= (initializer_list<std::tuple<PointIndex,PointGeomInfo>> list) 453 { 454 size_t cnt = 0; 455 for (auto val : list) 456 { 457 pnum[cnt] = get<0>(val); 458 geominfo[cnt++] = get<1>(val); 459 } 460 return *this; 461 } 462 /// 463 DLL_HEADER Element2d (int anp); 464 /// 465 DLL_HEADER Element2d (ELEMENT_TYPE type); 466 /// 467 DLL_HEADER Element2d (int pi1, int pi2, int pi3); 468 /// 469 DLL_HEADER Element2d (int pi1, int pi2, int pi3, int pi4); 470 /// GetType() const471 ELEMENT_TYPE GetType () const { return typ; } 472 /// SetType(ELEMENT_TYPE atyp)473 void SetType (ELEMENT_TYPE atyp) 474 { 475 typ = atyp; 476 switch (typ) 477 { 478 case TRIG: np = 3; break; 479 case QUAD: np = 4; break; 480 case TRIG6: np = 6; break; 481 case QUAD6: np = 6; break; 482 case QUAD8: np = 8; break; 483 default: 484 PrintSysError ("Element2d::SetType, illegal type ", int(typ)); 485 } 486 is_curved = (np >= 4); 487 } 488 /// GetNP() const489 int GetNP() const { return np; } 490 /// GetNV() const491 int GetNV() const 492 { 493 if (typ == TRIG || typ == TRIG6) 494 return 3; 495 else 496 { 497 #ifdef DEBUG 498 if (typ != QUAD && typ != QUAD6 && typ != QUAD8) 499 PrintSysError ("element2d::GetNV not implemented for typ", int(typ)); 500 #endif 501 return 4; 502 } 503 /* 504 switch (typ) 505 { 506 case TRIG: 507 case TRIG6: return 3; 508 509 case QUAD: 510 case QUAD8: 511 case QUAD6: return 4; 512 default: 513 #ifdef DEBUG 514 PrintSysError ("element2d::GetNV not implemented for typ", typ) 515 #endif 516 ; 517 } 518 return np; 519 */ 520 } 521 522 /// operator [](int i)523 PointIndex & operator[] (int i) { return pnum[i]; } 524 /// operator [](int i) const525 const PointIndex & operator[] (int i) const { return pnum[i]; } 526 PNums() const527 auto PNums () const { return FlatArray<const PointIndex> (np, &pnum[0]); } PNums()528 auto PNums () { return FlatArray<PointIndex> (np, &pnum[0]); } 529 template <int NP> PNums() const530 auto PNums() const { return FlatArray<const PointIndex> (NP, &pnum[0]); } Vertices() const531 auto Vertices() const { return FlatArray<const PointIndex> (GetNV(), &pnum[0]); } 532 GeomInfo() const533 auto GeomInfo() const { return FlatArray<const PointGeomInfo> (np, &geominfo[0]); } GeomInfo()534 auto GeomInfo() { return FlatArray<PointGeomInfo> (np, &geominfo[0]); } 535 536 /// PNum(int i)537 PointIndex & PNum (int i) { return pnum[i-1]; } 538 /// PNum(int i) const539 const PointIndex & PNum (int i) const { return pnum[i-1]; } 540 /// PNumMod(int i)541 PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } 542 /// PNumMod(int i) const543 const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } 544 /// 545 546 /// GeomInfoPi(int i)547 PointGeomInfo & GeomInfoPi (int i) { return geominfo[i-1]; } 548 /// GeomInfoPi(int i) const549 const PointGeomInfo & GeomInfoPi (int i) const { return geominfo[i-1]; } 550 /// GeomInfoPiMod(int i)551 PointGeomInfo & GeomInfoPiMod (int i) { return geominfo[(i-1) % np]; } 552 /// GeomInfoPiMod(int i) const553 const PointGeomInfo & GeomInfoPiMod (int i) const { return geominfo[(i-1) % np]; } 554 DoArchive(Archive & ar)555 void DoArchive (Archive & ar) 556 { 557 short _np, _typ; 558 bool _curved, _vis, _deleted; 559 if (ar.Output()) 560 { _np = np; _typ = typ; _curved = is_curved; 561 _vis = visible; _deleted = deleted; } 562 // ar & _np & _typ & index & _curved & _vis & _deleted; 563 ar.DoPacked (_np, _typ, index, _curved, _vis, _deleted); 564 // ar & next; don't need 565 if (ar.Input()) 566 { np = _np; typ = ELEMENT_TYPE(_typ); is_curved = _curved; 567 visible = _vis; deleted = _deleted; } 568 /* 569 for (size_t i = 0; i < np; i++) 570 ar & pnum[i]; 571 */ 572 static_assert(sizeof(int) == sizeof (PointIndex)); 573 ar.Do( (int*)&pnum[0], np); 574 } 575 576 #ifdef PARALLEL 577 static MPI_Datatype MyGetMPIType(); 578 #endif 579 580 SetIndex(int si)581 void SetIndex (int si) { index = si; } 582 /// GetIndex() const583 int GetIndex () const { return index; } 584 GetOrder() const585 int GetOrder () const { return orderx; } SetOrder(int aorder)586 void SetOrder (int aorder) { orderx = ordery = aorder; } 587 588 GetOrder(int & ox,int & oy) const589 void GetOrder (int & ox, int & oy) const { ox = orderx, oy =ordery;}; GetOrder(int & ox,int & oy,int & oz) const590 void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz=0; } SetOrder(int ox,int oy,int)591 void SetOrder (int ox, int oy, int /* oz */) { orderx = ox; ordery = oy;} SetOrder(int ox,int oy)592 void SetOrder (int ox, int oy) { orderx = ox; ordery = oy;} 593 594 595 /// 596 void GetBox (const T_POINTS & points, Box3d & box) const; 597 /// invert orientation 598 inline void Invert (); 599 /// 600 void Invert2 (); 601 /// first point number is smallest 602 inline void NormalizeNumbering (); 603 /// 604 void NormalizeNumbering2 (); 605 BadElement() const606 bool BadElement() const { return badel; } 607 608 // friend ostream & operator<<(ostream & s, const Element2d & el); 609 friend class Mesh; 610 611 612 /// get number of 'integration points' 613 int GetNIP () const; 614 void GetIntegrationPoint (int ip, Point<2> & p, double & weight) const; 615 616 void GetTransformation (int ip, const NgArray<Point<2>> & points, 617 class DenseMatrix & trans) const; 618 void GetTransformation (int ip, class DenseMatrix & pmat, 619 class DenseMatrix & trans) const; 620 621 void GetShape (const Point<2> & p, class Vector & shape) const; 622 void GetShapeNew (const Point<2> & p, class FlatVector & shape) const; 623 template <typename T> 624 void GetShapeNew (const Point<2,T> & p, TFlatVector<T> shape) const; 625 /// matrix 2 * np 626 void GetDShape (const Point<2> & p, class DenseMatrix & dshape) const; 627 template <typename T> 628 void GetDShapeNew (const Point<2,T> & p, class MatrixFixWidth<2,T> & dshape) const; 629 630 /// matrix 2 * np 631 void GetPointMatrix (const NgArray<Point<2>> & points, 632 class DenseMatrix & pmat) const; 633 634 void ComputeIntegrationPointData () const; 635 636 637 double CalcJacobianBadness (const NgArray<Point<2>> & points) const; 638 double CalcJacobianBadness (const T_POINTS & points, 639 const Vec<3> & n) const; 640 double CalcJacobianBadnessDirDeriv (const NgArray<Point<2>> & points, 641 int pi, Vec<2> & dir, double & dd) const; 642 643 644 Delete()645 void Delete () 646 { 647 deleted = 1; 648 // for (PointIndex & p : pnum) p.Invalidate(); 649 } 650 IsDeleted() const651 bool IsDeleted () const 652 { 653 #ifdef DEBUG 654 if (pnum[0] < PointIndex::BASE && !deleted) 655 cerr << "Surfelement has illegal pnum, but not marked as deleted" << endl; 656 #endif 657 return deleted; 658 } 659 660 // Philippose - 08 August 2010 661 // Access functions for the new property: visible Visible(bool vis=true)662 void Visible(bool vis = true) 663 { visible = vis; } IsVisible() const664 bool IsVisible () const 665 { return visible; } 666 SetRefinementFlag(bool rflag=true)667 void SetRefinementFlag (bool rflag = true) 668 { refflag = rflag; } TestRefinementFlag() const669 bool TestRefinementFlag () const 670 { return refflag; } 671 SetStrongRefinementFlag(bool rflag=true)672 void SetStrongRefinementFlag (bool rflag = true) 673 { strongrefflag = rflag; } TestStrongRefinementFlag() const674 bool TestStrongRefinementFlag () const 675 { return strongrefflag; } 676 677 IsCurved() const678 bool IsCurved () const { return is_curved; } SetCurved(bool acurved)679 void SetCurved (bool acurved) { is_curved = acurved; } 680 NextElement()681 SurfaceElementIndex NextElement() { return next; } 682 683 bool operator==(const Element2d & el2) const; 684 685 int HasFace(const Element2d& el) const; 686 /// 687 int meshdocval; 688 /// 689 int hp_elnr; 690 691 /* 692 #ifdef PARALLEL 693 int GetPartition () const { return partitionNumber; } 694 void SetPartition (int nr) { partitionNumber = nr; }; 695 #endif 696 */ 697 }; 698 699 ostream & operator<<(ostream & s, const Element2d & el); 700 701 702 703 704 705 class IntegrationPointData 706 { 707 public: 708 Point<3> p; 709 double weight; 710 Vector shape; 711 DenseMatrix dshape; 712 }; 713 714 715 716 717 718 719 720 721 /** 722 Volume element 723 */ 724 class Element 725 { 726 private: 727 /// point numbers 728 PointIndex pnum[ELEMENT_MAXPOINTS]; 729 /// 730 ELEMENT_TYPE typ; 731 /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) 732 int8_t np; 733 /// 734 class flagstruct { 735 public: 736 bool marked:1; // marked for refinement 737 bool badel:1; // angles worse then limit 738 bool reverse:1; // for refinement a la Bey 739 bool illegal:1; // illegal, will be split or swapped 740 bool illegal_valid:1; // is illegal-flag valid ? 741 bool badness_valid:1; // is badness valid ? 742 bool refflag:1; // mark element for refinement 743 bool strongrefflag:1; 744 bool deleted:1; // element is deleted, will be removed from array 745 bool fixed:1; // don't change element in optimization 746 }; 747 748 /// sub-domain index 749 int index; 750 /// order for hp-FEM 751 unsigned int orderx:6; 752 unsigned int ordery:6; 753 unsigned int orderz:6; 754 /* unsigned int levelx:6; 755 unsigned int levely:6; 756 unsigned int levelz:6; */ 757 /// stored shape-badness of element 758 float badness; 759 bool is_curved:1; // element is (high order) curved 760 761 // #ifdef PARALLEL 762 /// number of partition for parallel computation 763 // int partitionNumber; 764 765 // #endif 766 767 public: 768 flagstruct flags; 769 770 /// 771 DLL_HEADER Element () = default; 772 Element (const Element &) = default; 773 Element (Element &&) = default; 774 Element & operator= (const Element &) = default; 775 Element & operator= (Element &&) = default; 776 777 /// 778 Element (int anp); 779 /// 780 Element (ELEMENT_TYPE type); 781 /// 782 // Element & operator= (const Element & el2); 783 784 /// 785 void SetNP (int anp); 786 /// 787 void SetType (ELEMENT_TYPE atyp); 788 /// GetNP() const789 int GetNP () const { return np; } 790 /// GetNV() const791 uint8_t GetNV() const 792 { 793 __assume(typ >= TET && typ <= PYRAMID13); 794 switch (typ) 795 { 796 case TET: 797 case TET10: 798 return 4; 799 case PRISM12: 800 case PRISM15: 801 case PRISM: 802 return 6; 803 case PYRAMID: 804 case PYRAMID13: 805 return 5; 806 case HEX: 807 case HEX20: 808 return 8; 809 default: // not a 3D element 810 #ifdef DEBUG 811 PrintSysError ("Element3d::GetNV not implemented for typ ", int(typ)); 812 #endif 813 __assume(false); 814 return -1; 815 } 816 } 817 818 bool operator==(const Element & el2) const; 819 820 // old style: NP() const821 int NP () const { return np; } 822 823 /// GetType() const824 ELEMENT_TYPE GetType () const { return typ; } 825 826 /// operator [](int i)827 PointIndex & operator[] (int i) { return pnum[i]; } 828 /// operator [](int i) const829 const PointIndex & operator[] (int i) const { return pnum[i]; } 830 PNums() const831 auto PNums () const { return FlatArray<const PointIndex> (np, &pnum[0]); } PNums()832 auto PNums () { return FlatArray<PointIndex> (np, &pnum[0]); } 833 template <int NP> PNums() const834 auto PNums() const { return FlatArray<const PointIndex> (NP, &pnum[0]); } 835 Vertices() const836 FlatArray<const PointIndex> Vertices() const { return { GetNV(), &pnum[0] }; } 837 838 /// PNum(int i)839 PointIndex & PNum (int i) { return pnum[i-1]; } 840 /// PNum(int i) const841 const PointIndex & PNum (int i) const { return pnum[i-1]; } 842 /// PNumMod(int i)843 PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } 844 /// PNumMod(int i) const845 const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } 846 DoArchive(Archive & ar)847 void DoArchive (Archive & ar) 848 { 849 short _np, _typ; 850 bool _curved; 851 if (ar.Output()) 852 { _np = np; _typ = typ; _curved = is_curved; } 853 // ar & _np & _typ & index & _curved; 854 ar.DoPacked (_np, _typ, index, _curved); 855 856 if (ar.Input()) 857 { 858 np = _np; 859 typ = ELEMENT_TYPE(_typ); 860 is_curved = _curved; 861 flags.marked = 1; 862 flags.badel = 0; 863 flags.reverse = 0; 864 flags.illegal = 0; 865 flags.illegal_valid = 0; 866 flags.badness_valid = 0; 867 flags.refflag = 1; 868 flags.strongrefflag = false; 869 flags.deleted = 0; 870 flags.fixed = 0; 871 } 872 873 static_assert(sizeof(int) == sizeof (PointIndex)); 874 ar.Do( (int*)&pnum[0], np); 875 } 876 877 #ifdef PARALLEL 878 static MPI_Datatype MyGetMPIType(); 879 #endif 880 881 /// SetIndex(int si)882 void SetIndex (int si) { index = si; } 883 /// GetIndex() const884 int GetIndex () const { return index; } 885 GetOrder() const886 int GetOrder () const { return orderx; } 887 void SetOrder (const int aorder) ; 888 GetOrder(int & ox,int & oy,int & oz) const889 void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz = orderz; } 890 void SetOrder (const int ox, const int oy, const int oz); 891 // void GetLevel (int & ox, int & oy, int & oz) const { ox = levelx; oy = levely; oz = levelz; } 892 // void SetLevel (int ox, int oy, int oz) { levelx = ox; levely = oy; levelz = oz; } 893 894 895 /// 896 void GetBox (const T_POINTS & points, Box3d & box) const; 897 /// Calculates Volume of element 898 double Volume (const T_POINTS & points) const; 899 /// 900 DLL_HEADER void Print (ostream & ost) const; 901 /// GetNFaces() const902 int GetNFaces () const 903 { 904 switch (typ) 905 { 906 case TET: 907 case TET10: return 4; 908 case PYRAMID: case PYRAMID13: return 5; 909 case PRISM: 910 case PRISM15: 911 case PRISM12: return 5; 912 case HEX: case HEX20: 913 return 6; 914 default: 915 #ifdef DEBUG 916 PrintSysError ("element3d::GetNFaces not implemented for typ", int(typ)) 917 #endif 918 ; 919 } 920 return 0; 921 } 922 /// 923 inline void GetFace (int i, Element2d & face) const; 924 /// 925 void GetFace2 (int i, Element2d & face) const; 926 /// 927 void Invert (); 928 929 /// split into 4 node tets 930 void GetTets (NgArray<Element> & locels) const; 931 /// split into 4 node tets, local point nrs 932 void GetTetsLocal (NgArray<Element> & locels) const; 933 /// returns coordinates of nodes 934 // void GetNodesLocal (NgArray<Point<3> > & points) const; 935 void GetNodesLocalNew (NgArray<Point<3> > & points) const; 936 937 /// split surface into 3 node trigs 938 void GetSurfaceTriangles (NgArray<Element2d> & surftrigs) const; 939 940 941 /// get number of 'integration points' 942 int GetNIP () const; 943 void GetIntegrationPoint (int ip, Point<3> & p, double & weight) const; 944 945 void GetTransformation (int ip, const T_POINTS & points, 946 class DenseMatrix & trans) const; 947 void GetTransformation (int ip, class DenseMatrix & pmat, 948 class DenseMatrix & trans) const; 949 950 void GetShape (const Point<3> & p, class Vector & shape) const; 951 // void GetShapeNew (const Point<3> & p, class FlatVector & shape) const; 952 template <typename T> 953 void GetShapeNew (const Point<3,T> & p, TFlatVector<T> shape) const; 954 /// matrix 2 * np 955 void GetDShape (const Point<3> & p, class DenseMatrix & dshape) const; 956 template <typename T> 957 void GetDShapeNew (const Point<3,T> & p, class MatrixFixWidth<3,T> & dshape) const; 958 /// matrix 3 * np 959 void GetPointMatrix (const T_POINTS & points, 960 class DenseMatrix & pmat) const; 961 962 void ComputeIntegrationPointData () const; 963 964 965 double CalcJacobianBadness (const T_POINTS & points) const; 966 double CalcJacobianBadnessDirDeriv (const T_POINTS & points, 967 int pi, Vec<3> & dir, double & dd) const; 968 double CalcJacobianBadnessGradient (const T_POINTS & points, 969 int pi, Vec<3> & grad) const; 970 971 /// 972 // friend ostream & operator<<(ostream & s, const Element & el); 973 SetRefinementFlag(bool rflag=1)974 void SetRefinementFlag (bool rflag = 1) 975 { flags.refflag = rflag; } TestRefinementFlag() const976 int TestRefinementFlag () const 977 { return flags.refflag; } 978 SetStrongRefinementFlag(bool rflag=1)979 void SetStrongRefinementFlag (bool rflag = 1) 980 { flags.strongrefflag = rflag; } TestStrongRefinementFlag() const981 int TestStrongRefinementFlag () const 982 { return flags.strongrefflag; } 983 Illegal() const984 int Illegal () const 985 { return flags.illegal; } IllegalValid() const986 int IllegalValid () const 987 { return flags.illegal_valid; } SetIllegal(int aillegal)988 void SetIllegal (int aillegal) 989 { 990 flags.illegal = aillegal ? 1 : 0; 991 flags.illegal_valid = 1; 992 } SetLegal(int alegal)993 void SetLegal (int alegal) 994 { 995 flags.illegal = alegal ? 0 : 1; 996 flags.illegal_valid = 1; 997 } 998 Delete()999 void Delete () { flags.deleted = 1; } IsDeleted() const1000 bool IsDeleted () const 1001 { 1002 #ifdef DEBUG 1003 if (pnum[0] < PointIndex::BASE && !flags.deleted) 1004 cerr << "Volelement has illegal pnum, but not marked as deleted" << endl; 1005 #endif 1006 1007 return flags.deleted; 1008 } 1009 IsCurved() const1010 bool IsCurved () const { return is_curved; } SetCurved(bool acurved)1011 void SetCurved (bool acurved) { is_curved = acurved; } 1012 1013 /* 1014 #ifdef PARALLEL 1015 int GetPartition () const { return partitionNumber; } 1016 void SetPartition (int nr) { partitionNumber = nr; }; 1017 #else 1018 int GetPartition () const { return 0; } 1019 #endif 1020 */ 1021 1022 int hp_elnr; 1023 }; 1024 1025 ostream & operator<<(ostream & s, const Element & el); 1026 1027 1028 1029 1030 1031 1032 /** 1033 Edge segment. 1034 */ 1035 class Segment 1036 { 1037 public: 1038 /// 1039 DLL_HEADER Segment(); 1040 DLL_HEADER Segment (const Segment& other); 1041 ~Segment()1042 ~Segment() 1043 { ; } 1044 1045 // friend ostream & operator<<(ostream & s, const Segment & seg); 1046 1047 PointIndex pnums[3]; // p1, p2, pmid 1048 1049 int edgenr; 1050 /// 1051 double singedge_left; 1052 double singedge_right; 1053 1054 /// 0.. not first segment of segs, 1..first of class, 2..first of class, inverse 1055 unsigned int seginfo:2; 1056 1057 /// surface decoding index 1058 int si; 1059 /// co dim 2 deconding index 1060 int cd2i; 1061 /// domain number inner side 1062 int domin; 1063 /// domain number outer side 1064 int domout; 1065 /// top-level object number of surface 1066 int tlosurf; 1067 /// 1068 PointGeomInfo geominfo[2]; 1069 1070 /// surfaces describing edge 1071 int surfnr1, surfnr2; 1072 /// 1073 EdgePointGeomInfo epgeominfo[2]; 1074 /// 1075 // int pmid; // for second order 1076 /// 1077 int meshdocval; 1078 1079 // #ifdef PARALLEL 1080 /// number of partition for parallel computation 1081 // int partitionNumber; 1082 // #endif 1083 1084 private: 1085 bool is_curved; 1086 1087 public: 1088 /* 1089 PointIndex operator[] (int i) const 1090 { return (i == 0) ? p1 : p2; } 1091 1092 PointIndex & operator[] (int i) 1093 { return (i == 0) ? p1 : p2; } 1094 */ 1095 1096 Segment& operator=(const Segment & other); 1097 1098 1099 int hp_elnr; 1100 GetNP() const1101 int GetNP() const 1102 { 1103 return pnums[2].IsValid() ? 3 : 2; 1104 } 1105 PNums() const1106 auto PNums() const { return FlatArray<const PointIndex> (GetNP(), &pnums[0]); } PNums()1107 auto PNums() { return FlatArray<PointIndex> (GetNP(), &pnums[0]); } 1108 1109 GetType() const1110 ELEMENT_TYPE GetType() const 1111 { 1112 return pnums[2].IsValid() ? SEGMENT3 : SEGMENT; 1113 } 1114 operator [](int i)1115 PointIndex & operator[] (int i) { return pnums[i]; } operator [](int i) const1116 const PointIndex & operator[] (int i) const { return pnums[i]; } 1117 1118 IsCurved() const1119 bool IsCurved () const { return is_curved; } SetCurved(bool acurved)1120 void SetCurved (bool acurved) { is_curved = acurved; } 1121 1122 /* 1123 #ifdef PARALLEL 1124 int GetPartition () const { return partitionNumber; } 1125 void SetPartition (int nr) { partitionNumber = nr; }; 1126 #else 1127 int GetPartition () const { return 0; } 1128 #endif 1129 */ 1130 1131 void DoArchive (Archive & ar); 1132 #ifdef PARALLEL 1133 static MPI_Datatype MyGetMPIType(); 1134 #endif 1135 1136 }; 1137 1138 ostream & operator<<(ostream & s, const Segment & seg); 1139 1140 1141 class Element0d 1142 { 1143 public: 1144 PointIndex pnum; 1145 string name; 1146 int index; 1147 Element0d () = default; Element0d(PointIndex _pnum,int _index)1148 Element0d (PointIndex _pnum, int _index) 1149 : pnum(_pnum), index(_index) { ; } 1150 }; 1151 1152 ostream & operator<<(ostream & s, const Element0d & el); 1153 1154 // class Surface; 1155 // class FaceDescriptor; 1156 1157 /// 1158 class FaceDescriptor 1159 { 1160 /// which surface, 0 if not available 1161 int surfnr; 1162 /// domain nr inside 1163 int domin; 1164 /// domain nr outside 1165 int domout; 1166 /// top level object number of surface 1167 int tlosurf; 1168 /// boundary condition property 1169 int bcprop; 1170 // Philippose - 06/07/2009 1171 // Add capability to store surface colours along with 1172 // other face data 1173 /// surface colour (Default: R=0.0 ; G=1.0 ; B=0.0) 1174 Vec<3> surfcolour; 1175 1176 /// 1177 static string default_bcname; 1178 string * bcname = &default_bcname; 1179 /// root of linked list 1180 SurfaceElementIndex firstelement; 1181 1182 double domin_singular; 1183 double domout_singular; 1184 1185 public: 1186 DLL_HEADER FaceDescriptor(); 1187 DLL_HEADER FaceDescriptor(int surfnri, int domini, int domouti, int tlosurfi); 1188 DLL_HEADER FaceDescriptor(const Segment & seg); 1189 DLL_HEADER FaceDescriptor(const FaceDescriptor& other); ~FaceDescriptor()1190 DLL_HEADER ~FaceDescriptor() { ; } 1191 1192 DLL_HEADER int SegmentFits (const Segment & seg); 1193 SurfNr() const1194 int SurfNr () const { return surfnr; } DomainIn() const1195 int DomainIn () const { return domin; } DomainOut() const1196 int DomainOut () const { return domout; } TLOSurface() const1197 int TLOSurface () const { return tlosurf; } BCProperty() const1198 int BCProperty () const { return bcprop; } 1199 1200 DomainInSingular() const1201 double DomainInSingular() const { return domin_singular; } DomainOutSingular() const1202 double DomainOutSingular() const { return domout_singular; } 1203 1204 // Philippose - 06/07/2009 1205 // Get Surface colour SurfColour() const1206 Vec<3> SurfColour () const { return surfcolour; } GetBCName() const1207 DLL_HEADER const string & GetBCName () const { return *bcname; } 1208 // string * BCNamePtr () { return bcname; } 1209 // const string * BCNamePtr () const { return bcname; } SetSurfNr(int sn)1210 void SetSurfNr (int sn) { surfnr = sn; } SetDomainIn(int di)1211 void SetDomainIn (int di) { domin = di; } SetDomainOut(int dom)1212 void SetDomainOut (int dom) { domout = dom; } SetBCProperty(int bc)1213 void SetBCProperty (int bc) { bcprop = bc; } 1214 void SetBCName (string * bcn); // { bcname = bcn; } 1215 // Philippose - 06/07/2009 1216 // Set the surface colour SetSurfColour(Vec<3> colour)1217 void SetSurfColour (Vec<3> colour) { surfcolour = colour; } 1218 SetDomainInSingular(double v)1219 void SetDomainInSingular (double v) { domin_singular = v; } SetDomainOutSingular(double v)1220 void SetDomainOutSingular (double v) { domout_singular = v; } 1221 FirstElement()1222 SurfaceElementIndex FirstElement() { return firstelement; } 1223 // friend ostream & operator<<(ostream & s, const FaceDescriptor & fd); 1224 friend class Mesh; 1225 1226 void DoArchive (Archive & ar); 1227 }; 1228 1229 ostream & operator<< (ostream & s, const FaceDescriptor & fd); 1230 1231 1232 class EdgeDescriptor 1233 { 1234 int tlosurf; 1235 int surfnr[2]; 1236 public: EdgeDescriptor()1237 EdgeDescriptor () 1238 : tlosurf(-1) 1239 { surfnr[0] = surfnr[1] = -1; } 1240 SurfNr(int i) const1241 int SurfNr (int i) const { return surfnr[i]; } SetSurfNr(int i,int nr)1242 void SetSurfNr (int i, int nr) { surfnr[i] = nr; } 1243 TLOSurface() const1244 int TLOSurface() const { return tlosurf; } SetTLOSurface(int nr)1245 void SetTLOSurface (int nr) { tlosurf = nr; } 1246 }; 1247 1248 1249 1250 class DLL_HEADER MeshingParameters 1251 { 1252 public: 1253 /** 1254 3d optimization strategy: 1255 // m .. move nodes 1256 // M .. move nodes, cheap functional 1257 // s .. swap faces 1258 // c .. combine elements 1259 // d .. divide elements 1260 // D .. divide and join opposite edges, remove element 1261 // p .. plot, no pause 1262 // P .. plot, Pause 1263 // h .. Histogramm, no pause 1264 // H .. Histogramm, pause 1265 */ 1266 string optimize3d = "cmdDmustm"; 1267 /// number of 3d optimization steps 1268 int optsteps3d = 3; 1269 /** 1270 2d optimization strategy: 1271 // s .. swap, opt 6 lines/node 1272 // S .. swap, optimal elements 1273 // m .. move nodes 1274 // p .. plot, no pause 1275 // P .. plot, pause 1276 // c .. combine 1277 **/ 1278 string optimize2d = "smcmSmcmSmcm"; 1279 /// number of 2d optimization steps 1280 int optsteps2d = 3; 1281 /// power of error (to approximate max err optimization) 1282 double opterrpow = 2; 1283 /// do block filling ? 1284 bool blockfill = true; 1285 /// block filling up to distance 1286 double filldist = 0.1; 1287 /// radius of local environment (times h) 1288 double safety = 5; 1289 /// radius of active environment (times h) 1290 double relinnersafety = 3; 1291 /// use local h ? 1292 bool uselocalh = true; 1293 /// grading for local h 1294 double grading = 0.3; 1295 /// use delaunay for 3d meshing 1296 bool delaunay = true; 1297 /// use delaunay for 2d meshing 1298 bool delaunay2d = false; 1299 /// maximal mesh size 1300 double maxh = 1e10; 1301 /// minimal mesh size 1302 double minh = 0.0; 1303 /// file for meshsize 1304 string meshsizefilename = ""; 1305 /// restrict h based on close edges 1306 optional<double> closeedgefac = nullopt; 1307 /// start surfacemeshing from everywhere in surface 1308 bool startinsurface = false; 1309 /// check overlapping surfaces (debug) 1310 bool checkoverlap = true; 1311 /// check overlapping surface mesh before volume meshing 1312 bool checkoverlappingboundary = true; 1313 /// check chart boundary (sometimes too restrictive) 1314 bool checkchartboundary = true; 1315 /// safety factor for curvatures (elements per radius) 1316 double curvaturesafety = 2; 1317 /// minimal number of segments per edge 1318 double segmentsperedge = 1; 1319 /// use parallel threads 1320 bool parthread = 0; 1321 /// weight of element size w.r.t element shape 1322 double elsizeweight = 0.2; 1323 /// init with default values 1324 1325 /// start at step 1326 int perfstepsstart = 0; 1327 /// end at step 1328 int perfstepsend = 6; 1329 1330 1331 /// from mp3: 1332 /// give up quality class, 2d meshing 1333 int giveuptol2d = 200; 1334 /// give up quality class, 3d meshing 1335 int giveuptol = 10; 1336 /// maximal outer steps 1337 int maxoutersteps = 10; 1338 /// class starting star-shape filling 1339 int starshapeclass = 5; 1340 /// if non-zero, baseelement must have baseelnp points 1341 int baseelnp = 0; 1342 /// quality tolerances are handled less careful 1343 int sloppy = 1; 1344 1345 /// limit for max element angle (150-180) 1346 double badellimit = 175; 1347 1348 bool check_impossible = false; 1349 1350 int only3D_domain_nr = 0; 1351 1352 /// 1353 bool secondorder = false; 1354 /// high order element curvature 1355 int elementorder = 1; 1356 /// quad-dominated surface meshing 1357 bool quad = false; 1358 /// 1359 bool try_hexes = false; 1360 /// 1361 bool inverttets = false; 1362 /// 1363 bool inverttrigs = false; 1364 /// 1365 bool autozrefine = false; 1366 1367 bool parallel_meshing = true; 1368 int nthreads = 4; 1369 1370 Flags geometrySpecificParameters; 1371 /// 1372 MeshingParameters (); 1373 /// 1374 MeshingParameters (const MeshingParameters & mp2) = default; 1375 MeshingParameters (MeshingParameters && mp2) = default; 1376 MeshingParameters & operator= (const MeshingParameters & mp2) = default; 1377 MeshingParameters & operator= (MeshingParameters && mp2) = default; 1378 /// 1379 void Print (ostream & ost) const; 1380 /// 1381 // void CopyFrom(const MeshingParameters & other); 1382 1383 class MeshSizePoint 1384 { 1385 public: 1386 Point<3> pnt; 1387 double h; MeshSizePoint(Point<3> _pnt,double _h)1388 MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; } 1389 MeshSizePoint () = default; 1390 MeshSizePoint (const MeshSizePoint &) = default; 1391 MeshSizePoint (MeshSizePoint &&) = default; 1392 MeshSizePoint & operator= (const MeshSizePoint &) = default; 1393 MeshSizePoint & operator= (MeshSizePoint &&) = default; 1394 }; 1395 NgArray<MeshSizePoint> meshsize_points; 1396 1397 void (*render_function)(bool) = NULL; Render(bool blocking=false) const1398 void Render(bool blocking = false) const 1399 { 1400 if (render_function) 1401 (*render_function)(blocking); 1402 } 1403 }; 1404 operator <<(ostream & ost,const MeshingParameters & mp)1405 inline ostream & operator<< (ostream & ost, const MeshingParameters & mp) 1406 { 1407 mp.Print (ost); 1408 return ost; 1409 } 1410 1411 class DebugParameters 1412 { 1413 public: 1414 /// 1415 int debugoutput; 1416 /// use slow checks 1417 int slowchecks; 1418 /// 1419 int haltsuccess; 1420 /// 1421 int haltnosuccess; 1422 /// 1423 int haltlargequalclass; 1424 /// 1425 int haltsegment; 1426 /// 1427 int haltnode; 1428 /// 1429 int haltsegmentp1; 1430 /// 1431 int haltsegmentp2; 1432 /// 1433 int haltexistingline; 1434 /// 1435 int haltoverlap; 1436 /// 1437 int haltface; 1438 /// 1439 int haltfacenr; 1440 /// 1441 DebugParameters (); 1442 }; 1443 1444 1445 1446 Invert()1447 inline void Element2d :: Invert() 1448 { 1449 if (typ == TRIG) 1450 Swap (PNum(2), PNum(3)); 1451 else 1452 Invert2(); 1453 } 1454 1455 1456 1457 NormalizeNumbering()1458 inline void Element2d :: NormalizeNumbering () 1459 { 1460 if (GetNP() == 3) 1461 { 1462 if (PNum(1) < PNum(2) && PNum(1) < PNum(3)) 1463 return; 1464 else 1465 { 1466 if (PNum(2) < PNum(3)) 1467 { 1468 PointIndex pi1 = PNum(2); 1469 PNum(2) = PNum(3); 1470 PNum(3) = PNum(1); 1471 PNum(1) = pi1; 1472 } 1473 else 1474 { 1475 PointIndex pi1 = PNum(3); 1476 PNum(3) = PNum(2); 1477 PNum(2) = PNum(1); 1478 PNum(1) = pi1; 1479 } 1480 } 1481 } 1482 else 1483 NormalizeNumbering2(); 1484 } 1485 1486 1487 1488 static const int gftetfacesa[4][3] = 1489 { { 1, 2, 3 }, 1490 { 2, 0, 3 }, 1491 { 0, 1, 3 }, 1492 { 1, 0, 2 } }; 1493 GetFace(int i,Element2d & face) const1494 inline void Element :: GetFace (int i, Element2d & face) const 1495 { 1496 if (typ == TET) 1497 { 1498 face.SetType(TRIG); 1499 face[0] = pnum[gftetfacesa[i-1][0]]; 1500 face[1] = pnum[gftetfacesa[i-1][1]]; 1501 face[2] = pnum[gftetfacesa[i-1][2]]; 1502 } 1503 else 1504 GetFace2 (i, face); 1505 } 1506 1507 1508 1509 1510 1511 1512 1513 /** 1514 Identification of periodic surfaces, close surfaces, etc. 1515 */ 1516 class Identifications 1517 { 1518 public: 1519 enum ID_TYPE : unsigned char { UNDEFINED = 1, PERIODIC = 2, CLOSESURFACES = 3, CLOSEEDGES = 4}; 1520 1521 1522 private: 1523 class Mesh & mesh; 1524 1525 /// identify points (thin layers, periodic b.c.) 1526 INDEX_2_HASHTABLE<int> identifiedpoints; 1527 1528 /// the same, with info about the id-nr 1529 INDEX_3_HASHTABLE<int> identifiedpoints_nr; 1530 1531 /// sorted by identification nr 1532 TABLE<INDEX_2> idpoints_table; 1533 1534 NgArray<ID_TYPE> type; 1535 1536 /// number of identifications (or, actually used identifications ?) 1537 int maxidentnr; 1538 1539 public: 1540 /// 1541 DLL_HEADER Identifications (class Mesh & amesh); 1542 /// 1543 DLL_HEADER ~Identifications (); 1544 1545 DLL_HEADER void Delete (); 1546 1547 /* 1548 Identify points pi1 and pi2, due to 1549 identification nr identnr 1550 */ 1551 DLL_HEADER void Add (PointIndex pi1, PointIndex pi2, int identnr); 1552 1553 1554 int Get (PointIndex pi1, PointIndex pi2) const; 1555 int GetSymmetric (PointIndex pi1, PointIndex pi2) const; 1556 1557 bool Get (PointIndex pi1, PointIndex pi2, int identnr) const; 1558 bool GetSymmetric (PointIndex pi1, PointIndex pi2, int identnr) const; 1559 1560 // bool HasIdentifiedPoints() const { return identifiedpoints != nullptr; } 1561 /// GetIdentifiedPoints()1562 INDEX_2_HASHTABLE<int> & GetIdentifiedPoints () 1563 { 1564 return identifiedpoints; 1565 } 1566 Used(PointIndex pi1,PointIndex pi2)1567 bool Used (PointIndex pi1, PointIndex pi2) 1568 { 1569 return identifiedpoints.Used (INDEX_2 (pi1, pi2)); 1570 } 1571 UsedSymmetric(PointIndex pi1,PointIndex pi2)1572 bool UsedSymmetric (PointIndex pi1, PointIndex pi2) 1573 { 1574 return 1575 identifiedpoints.Used (INDEX_2 (pi1, pi2)) || 1576 identifiedpoints.Used (INDEX_2 (pi2, pi1)); 1577 } 1578 1579 /// 1580 void GetMap (int identnr, NgArray<int,PointIndex::BASE> & identmap, bool symmetric = false) const; 1581 /// GetType(int identnr) const1582 ID_TYPE GetType(int identnr) const 1583 { 1584 if(identnr <= type.Size()) 1585 return type[identnr-1]; 1586 else 1587 return UNDEFINED; 1588 } SetType(int identnr,ID_TYPE t)1589 void SetType(int identnr, ID_TYPE t) 1590 { 1591 while(type.Size() < identnr) 1592 type.Append(UNDEFINED); 1593 type[identnr-1] = t; 1594 } 1595 1596 /// 1597 void GetPairs (int identnr, NgArray<INDEX_2> & identpairs) const; 1598 /// GetMaxNr() const1599 int GetMaxNr () const { return maxidentnr; } 1600 1601 /// remove secondorder 1602 void SetMaxPointNr (int maxpnum); 1603 1604 DLL_HEADER void Print (ostream & ost) const; 1605 1606 void DoArchive (Archive & ar); 1607 }; 1608 } 1609 1610 1611 #ifdef PARALLEL 1612 namespace ngcore 1613 { 1614 template <> struct MPI_typetrait<netgen::PointIndex> { MPITypengcore::MPI_typetrait1615 static MPI_Datatype MPIType () { return MPI_INT; } 1616 }; 1617 1618 template <> struct MPI_typetrait<netgen::ELEMENT_TYPE> { MPITypengcore::MPI_typetrait1619 static MPI_Datatype MPIType () { return MPI_CHAR; } 1620 }; 1621 1622 template <> struct MPI_typetrait<netgen::MeshPoint> { MPITypengcore::MPI_typetrait1623 static MPI_Datatype MPIType () { return netgen::MeshPoint::MyGetMPIType(); } 1624 }; 1625 1626 template <> struct MPI_typetrait<netgen::Element> { MPITypengcore::MPI_typetrait1627 static MPI_Datatype MPIType () { return netgen::Element::MyGetMPIType(); } 1628 }; 1629 template <> struct MPI_typetrait<netgen::Element2d> { MPITypengcore::MPI_typetrait1630 static MPI_Datatype MPIType () { return netgen::Element2d::MyGetMPIType(); } 1631 }; 1632 template <> struct MPI_typetrait<netgen::Segment> { MPITypengcore::MPI_typetrait1633 static MPI_Datatype MPIType () { return netgen::Segment::MyGetMPIType(); } 1634 }; 1635 1636 } 1637 #endif 1638 1639 1640 #endif 1641 1642