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