diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 750d5f0f..80a40e1e 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -7,41 +7,42 @@ /* Date: 23. Aug. 09 */ /**************************************************************************/ -class Tcl_Interp; + +struct Tcl_Interp; namespace netgen { -class NetgenGeometry -{ -public: - virtual ~NetgenGeometry () { ; } + class NetgenGeometry + { + public: + virtual ~NetgenGeometry () { ; } - virtual int GenerateMesh (Mesh*& mesh, MeshingParameters & mparam, - int perfstepsstart, int perfstepsend); + virtual int GenerateMesh (Mesh*& mesh, MeshingParameters & mparam, + int perfstepsstart, int perfstepsend); - virtual const Refinement & GetRefinement () const; + virtual const Refinement & GetRefinement () const; - virtual void Save (string filename) const; - virtual void SaveToMeshFile (ostream & ost) const { ; } -}; + virtual void Save (string filename) const; + virtual void SaveToMeshFile (ostream & ost) const { ; } + }; -class GeometryRegister -{ -public: - virtual ~GeometryRegister(); - virtual NetgenGeometry * Load (string filename) const = 0; - virtual NetgenGeometry * LoadFromMeshFile (istream & ist) const { return NULL; } - virtual class VisualScene * GetVisualScene (const NetgenGeometry * geom) const - { return NULL; } - virtual void SetParameters (Tcl_Interp * interp) { ; } -}; + class GeometryRegister + { + public: + virtual ~GeometryRegister(); + virtual NetgenGeometry * Load (string filename) const = 0; + virtual NetgenGeometry * LoadFromMeshFile (istream & ist) const { return NULL; } + virtual class VisualScene * GetVisualScene (const NetgenGeometry * geom) const + { return NULL; } + virtual void SetParameters (Tcl_Interp * interp) { ; } + }; -extern Array geometryregister; + extern Array geometryregister; } diff --git a/libsrc/meshing/global.hpp b/libsrc/meshing/global.hpp index 2f2d3011..9805b1a2 100644 --- a/libsrc/meshing/global.hpp +++ b/libsrc/meshing/global.hpp @@ -12,39 +12,43 @@ global functions and variables */ -/// -DLL_HEADER extern double GetTime (); -extern void ResetTime (); - -/// -extern int testmode; - -/// calling parameters -// extern Flags parameters; - -extern MeshingParameters mparam; - -extern Array tets_in_qualclass; - -class multithreadt +namespace netgen { -public: - int pause; - int testmode; - int redraw; - int drawing; - int terminate; - int running; - double percent; - const char * task; - bool demorunning; - multithreadt(); -}; -extern volatile multithreadt multithread; + /// + DLL_HEADER extern double GetTime (); + extern void ResetTime (); -extern string ngdir; -extern DebugParameters debugparam; -extern bool verbose; + /// + extern int testmode; + + /// calling parameters + // extern Flags parameters; + + extern MeshingParameters mparam; + + extern Array tets_in_qualclass; + + class multithreadt + { + public: + int pause; + int testmode; + int redraw; + int drawing; + int terminate; + int running; + double percent; + const char * task; + bool demorunning; + multithreadt(); + }; + + extern volatile multithreadt multithread; + + extern string ngdir; + extern DebugParameters debugparam; + extern bool verbose; +} #endif diff --git a/libsrc/meshing/localh.hpp b/libsrc/meshing/localh.hpp index 583b6e0a..ee279c86 100644 --- a/libsrc/meshing/localh.hpp +++ b/libsrc/meshing/localh.hpp @@ -8,176 +8,180 @@ /**************************************************************************/ - - -/// box for grading -class GradingBox +namespace netgen { - /// xmid - float xmid[3]; - /// half edgelength - float h2; - /// - GradingBox * childs[8]; - /// - GradingBox * father; - /// - double hopt; - /// -public: - struct + + /// box for grading + class GradingBox { - unsigned int cutboundary:1; - unsigned int isinner:1; - unsigned int oldcell:1; - unsigned int pinner:1; - } flags; + /// xmid + float xmid[3]; + /// half edgelength + float h2; + /// + GradingBox * childs[8]; + /// + GradingBox * father; + /// + double hopt; + /// + public: - /// - GradingBox (const double * ax1, const double * ax2); - /// - void DeleteChilds(); - /// + struct + { + unsigned int cutboundary:1; + unsigned int isinner:1; + unsigned int oldcell:1; + unsigned int pinner:1; + } flags; - Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); } - double H2() const { return h2; } + /// + GradingBox (const double * ax1, const double * ax2); + /// + void DeleteChilds(); + /// - friend class LocalH; + Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); } + double H2() const { return h2; } - static BlockAllocator ball; - void * operator new(size_t); - void operator delete (void *); -}; + friend class LocalH; + + static BlockAllocator ball; + void * operator new(size_t); + void operator delete (void *); + }; -/** - Control of 3D mesh grading - */ -class LocalH -{ - /// - GradingBox * root; - /// - double grading; - /// - Array boxes; - /// - Box3d boundingbox; -public: - /// - LocalH (const Point3d & pmin, const Point3d & pmax, double grading); - /// - LocalH (const Box<3> & box, double grading); - /// - ~LocalH(); - /// - void Delete(); - /// - void SetGrading (double agrading) { grading = agrading; } - /// - void SetH (const Point3d & x, double h); - /// - double GetH (const Point3d & x) const; - /// minimal h in box (pmin, pmax) - double GetMinH (const Point3d & pmin, const Point3d & pmax) const; + /** + Control of 3D mesh grading + */ + class LocalH + { + /// + GradingBox * root; + /// + double grading; + /// + Array boxes; + /// + Box3d boundingbox; + public: + /// + LocalH (const Point3d & pmin, const Point3d & pmax, double grading); + /// + LocalH (const Box<3> & box, double grading); + /// + ~LocalH(); + /// + void Delete(); + /// + void SetGrading (double agrading) { grading = agrading; } + /// + void SetH (const Point3d & x, double h); + /// + double GetH (const Point3d & x) const; + /// minimal h in box (pmin, pmax) + double GetMinH (const Point3d & pmin, const Point3d & pmax) const; - /// mark boxes intersecting with boundary-box - // void CutBoundary (const Point3d & pmin, const Point3d & pmax) - // { CutBoundaryRec (pmin, pmax, root); } - void CutBoundary (const Box<3> & box) - { CutBoundaryRec (box.PMin(), box.PMax(), root); } + /// mark boxes intersecting with boundary-box + // void CutBoundary (const Point3d & pmin, const Point3d & pmax) + // { CutBoundaryRec (pmin, pmax, root); } + void CutBoundary (const Box<3> & box) + { CutBoundaryRec (box.PMin(), box.PMax(), root); } - /// find inner boxes - void FindInnerBoxes (class AdFront3 * adfront, - int (*testinner)(const Point3d & p1)); + /// find inner boxes + void FindInnerBoxes (class AdFront3 * adfront, + int (*testinner)(const Point3d & p1)); - void FindInnerBoxes (class AdFront2 * adfront, - int (*testinner)(const Point<2> & p1)); + void FindInnerBoxes (class AdFront2 * adfront, + int (*testinner)(const Point<2> & p1)); - /// clears all flags - void ClearFlags () + /// clears all flags + void ClearFlags () { ClearFlagsRec(root); } - /// widen refinement zone - void WidenRefinement (); + /// widen refinement zone + void WidenRefinement (); - /// get points in inner elements - void GetInnerPoints (Array > & points); + /// get points in inner elements + void GetInnerPoints (Array > & points); - /// get points in outer closure - void GetOuterPoints (Array > & points); + /// get points in outer closure + void GetOuterPoints (Array > & points); - /// - void Convexify (); - /// - int GetNBoxes () { return boxes.Size(); } - const Box3d & GetBoundingBox () const - { return boundingbox; } - /// - void PrintMemInfo (ostream & ost) const; -private: - /// - double GetMinHRec (const Point3d & pmin, const Point3d & pmax, - const GradingBox * box) const; - /// - void CutBoundaryRec (const Point3d & pmin, const Point3d & pmax, - GradingBox * box); + /// + void Convexify (); + /// + int GetNBoxes () { return boxes.Size(); } + const Box3d & GetBoundingBox () const + { return boundingbox; } + /// + void PrintMemInfo (ostream & ost) const; + private: + /// + double GetMinHRec (const Point3d & pmin, const Point3d & pmax, + const GradingBox * box) const; + /// + void CutBoundaryRec (const Point3d & pmin, const Point3d & pmax, + GradingBox * box); - /// - void FindInnerBoxesRec ( int (*inner)(const Point3d & p), - GradingBox * box); + /// + void FindInnerBoxesRec ( int (*inner)(const Point3d & p), + GradingBox * box); - /// - void FindInnerBoxesRec2 (GradingBox * box, - class AdFront3 * adfront, - Array & faceboxes, - Array & finds, int nfinbox); + /// + void FindInnerBoxesRec2 (GradingBox * box, + class AdFront3 * adfront, + Array & faceboxes, + Array & finds, int nfinbox); - void FindInnerBoxesRec ( int (*inner)(const Point<2> & p), - GradingBox * box); + void FindInnerBoxesRec ( int (*inner)(const Point<2> & p), + GradingBox * box); - /// - void FindInnerBoxesRec2 (GradingBox * box, - class AdFront2 * adfront, - Array > & faceboxes, - Array & finds, int nfinbox); + /// + void FindInnerBoxesRec2 (GradingBox * box, + class AdFront2 * adfront, + Array > & faceboxes, + Array & finds, int nfinbox); - /// - void SetInnerBoxesRec (GradingBox * box); + /// + void SetInnerBoxesRec (GradingBox * box); - /// - void ClearFlagsRec (GradingBox * box); + /// + void ClearFlagsRec (GradingBox * box); - /// - void ConvexifyRec (GradingBox * box); + /// + void ConvexifyRec (GradingBox * box); - friend ostream & operator<< (ostream & ost, const LocalH & loch); -}; + friend ostream & operator<< (ostream & ost, const LocalH & loch); + }; -inline ostream & operator<< (ostream & ost, const GradingBox & box) -{ - ost << "gradbox, pmid = " << box.PMid() << ", h2 = " << box.H2() - << " cutbound = " << box.flags.cutboundary << " isinner = " << box.flags.isinner - << endl; - return ost; -} + inline ostream & operator<< (ostream & ost, const GradingBox & box) + { + ost << "gradbox, pmid = " << box.PMid() << ", h2 = " << box.H2() + << " cutbound = " << box.flags.cutboundary << " isinner = " << box.flags.isinner + << endl; + return ost; + } + + inline ostream & operator<< (ostream & ost, const LocalH & loch) + { + for (int i = 0; i < loch.boxes.Size(); i++) + ost << "box[" << i << "] = " << *(loch.boxes[i]); + return ost; + } -inline ostream & operator<< (ostream & ost, const LocalH & loch) -{ - for (int i = 0; i < loch.boxes.Size(); i++) - ost << "box[" << i << "] = " << *(loch.boxes[i]); - return ost; } #endif diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 361cf385..1d3cbbc1 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -11,752 +11,754 @@ The mesh class */ - - -enum resthtype { RESTRICTH_FACE, RESTRICTH_EDGE, - RESTRICTH_SURFACEELEMENT, RESTRICTH_POINT, RESTRICTH_SEGMENT }; - -class HPRefElement; - - -/// 2d/3d mesh -class Mesh +namespace netgen { -public: - typedef ::netgen::T_POINTS T_POINTS; - // typedef MoveableArray T_POINTS; - // typedef MoveableArray T_VOLELEMENTS; - // typedef MoveableArray T_SURFELEMENTS; + enum resthtype { RESTRICTH_FACE, RESTRICTH_EDGE, + RESTRICTH_SURFACEELEMENT, RESTRICTH_POINT, RESTRICTH_SEGMENT }; - // typedef Array T_POINTS; - typedef Array T_VOLELEMENTS; - typedef Array T_SURFELEMENTS; + class HPRefElement; -private: - /// point coordinates - T_POINTS points; + /// 2d/3d mesh + class Mesh + { + public: + typedef ::netgen::T_POINTS T_POINTS; - /// line-segments at edges - Array segments; - /// surface elements, 2d-inner elements - T_SURFELEMENTS surfelements; - /// volume elements - T_VOLELEMENTS volelements; - /// points will be fixed forever - Array lockedpoints; + // typedef MoveableArray T_POINTS; + // typedef MoveableArray T_VOLELEMENTS; + // typedef MoveableArray T_SURFELEMENTS; + + // typedef Array T_POINTS; + typedef Array T_VOLELEMENTS; + typedef Array T_SURFELEMENTS; - /// surface indices at boundary nodes - TABLE surfacesonnode; - /// boundary edges (1..normal bedge, 2..segment) - INDEX_2_CLOSED_HASHTABLE * boundaryedges; - /// - INDEX_2_CLOSED_HASHTABLE * segmentht; - /// - INDEX_3_CLOSED_HASHTABLE * surfelementht; + private: + /// point coordinates + T_POINTS points; - /// faces of rest-solid - Array openelements; - /// open segmenets for surface meshing - Array opensegments; + /// line-segments at edges + Array segments; + /// surface elements, 2d-inner elements + T_SURFELEMENTS surfelements; + /// volume elements + T_VOLELEMENTS volelements; + /// points will be fixed forever + Array lockedpoints; + + + /// surface indices at boundary nodes + TABLE surfacesonnode; + /// boundary edges (1..normal bedge, 2..segment) + INDEX_2_CLOSED_HASHTABLE * boundaryedges; + /// + INDEX_2_CLOSED_HASHTABLE * segmentht; + /// + INDEX_3_CLOSED_HASHTABLE * surfelementht; + + /// faces of rest-solid + Array openelements; + /// open segmenets for surface meshing + Array opensegments; - /** - Representation of local mesh-size h - */ - LocalH * lochfunc; - /// - double hglob; - /// - double hmin; - /// - Array maxhdomain; + /** + Representation of local mesh-size h + */ + LocalH * lochfunc; + /// + double hglob; + /// + double hmin; + /// + Array maxhdomain; - /** - the face-index of the surface element maps into - this table. - */ - Array facedecoding; + /** + the face-index of the surface element maps into + this table. + */ + Array facedecoding; - /** - the edge-index of the line element maps into - this table. - */ - Array edgedecoding; + /** + the edge-index of the line element maps into + this table. + */ + Array edgedecoding; - /// sub-domain materials - Array materials; + /// sub-domain materials + Array materials; - /// labels for boundary conditions - Array bcnames; + /// labels for boundary conditions + Array bcnames; - /// Periodic surface, close surface, etc. identifications - Identifications * ident; + /// Periodic surface, close surface, etc. identifications + Identifications * ident; - /// number of vertices (if < 0, use np) - int numvertices; + /// number of vertices (if < 0, use np) + int numvertices; - /// geometric search tree for interval intersection search - Box3dTree * elementsearchtree; - /// time stamp for tree - int elementsearchtreets; + /// geometric search tree for interval intersection search + Box3dTree * elementsearchtree; + /// time stamp for tree + int elementsearchtreets; - /// element -> face, element -> edge etc ... - class MeshTopology * topology; - /// methods for high order elements - class CurvedElements * curvedelems; + /// element -> face, element -> edge etc ... + class MeshTopology * topology; + /// methods for high order elements + class CurvedElements * curvedelems; - /// nodes identified by close points - class AnisotropicClusters * clusters; + /// nodes identified by close points + class AnisotropicClusters * clusters; - /// space dimension (2 or 3) - int dimension; + /// space dimension (2 or 3) + int dimension; - /// changed by every minor modification (addpoint, ...) - int timestamp; - /// changed after finishing global algorithm (improve, ...) - int majortimestamp; + /// changed by every minor modification (addpoint, ...) + int timestamp; + /// changed after finishing global algorithm (improve, ...) + int majortimestamp; - /// mesh access semaphors. - NgMutex mutex; - /// mesh access semaphors. - NgMutex majormutex; + /// mesh access semaphors. + NgMutex mutex; + /// mesh access semaphors. + NgMutex majormutex; - SYMBOLTABLE< Array* > userdata_int; - SYMBOLTABLE< Array* > userdata_double; + SYMBOLTABLE< Array* > userdata_int; + SYMBOLTABLE< Array* > userdata_double; - mutable Array< Point3d > pointcurves; - mutable Array pointcurves_startpoint; - mutable Array pointcurves_red,pointcurves_green,pointcurves_blue; + mutable Array< Point3d > pointcurves; + mutable Array pointcurves_startpoint; + mutable Array pointcurves_red,pointcurves_green,pointcurves_blue; - /// start element for point search (GetElementOfPoint) - mutable int ps_startelement; + /// start element for point search (GetElementOfPoint) + mutable int ps_startelement; #ifdef PARALLEL - /// connection to parallel meshes - class ParallelMeshTopology * paralleltop; + /// connection to parallel meshes + class ParallelMeshTopology * paralleltop; #endif -private: - void BuildBoundaryEdges(void); + private: + void BuildBoundaryEdges(void); -public: - bool PointContainedIn2DElement(const Point3d & p, - double lami[3], - const int element, - bool consider3D = false) const; - bool PointContainedIn3DElement(const Point3d & p, - double lami[3], - const int element) const; - bool PointContainedIn3DElementOld(const Point3d & p, - double lami[3], - const int element) const; + public: + bool PointContainedIn2DElement(const Point3d & p, + double lami[3], + const int element, + bool consider3D = false) const; + bool PointContainedIn3DElement(const Point3d & p, + double lami[3], + const int element) const; + bool PointContainedIn3DElementOld(const Point3d & p, + double lami[3], + const int element) const; -public: + public: - // store coarse mesh before hp-refinement - Array * hpelements; - Mesh * coarsemesh; + // store coarse mesh before hp-refinement + Array * hpelements; + Mesh * coarsemesh; - /// number of refinement levels - int mglevels; - /// refinement hierarchy - Array mlbetweennodes; - /// parent element of volume element - Array mlparentelement; - /// parent element of surface element - Array mlparentsurfaceelement; + /// number of refinement levels + int mglevels; + /// refinement hierarchy + Array mlbetweennodes; + /// parent element of volume element + Array mlparentelement; + /// parent element of surface element + Array mlparentsurfaceelement; - /// - Mesh(); - /// - ~Mesh(); + /// + Mesh(); + /// + ~Mesh(); - Mesh & operator= (const Mesh & mesh2); + Mesh & operator= (const Mesh & mesh2); - /// - void DeleteMesh(); + /// + void DeleteMesh(); - /// - void ClearSurfaceElements(); + /// + void ClearSurfaceElements(); - /// - void ClearVolumeElements() - { - volelements.SetSize(0); - timestamp = NextTimeStamp(); - } + /// + void ClearVolumeElements() + { + volelements.SetSize(0); + timestamp = NextTimeStamp(); + } - /// - void ClearSegments() - { - segments.SetSize(0); - timestamp = NextTimeStamp(); - } + /// + void ClearSegments() + { + segments.SetSize(0); + timestamp = NextTimeStamp(); + } - /// - bool TestOk () const; + /// + bool TestOk () const; - void SetAllocSize(int nnodes, int nsegs, int nsel, int nel); + void SetAllocSize(int nnodes, int nsegs, int nsel, int nel); - PointIndex AddPoint (const Point3d & p, int layer = 1); - PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type); + PointIndex AddPoint (const Point3d & p, int layer = 1); + PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type); #ifdef PARALLEL - PointIndex AddPoint (const Point3d & p, bool aisghost, int layer = 1); - PointIndex AddPoint (const Point3d & p, bool aisghost, int layer, POINTTYPE type); + PointIndex AddPoint (const Point3d & p, bool aisghost, int layer = 1); + PointIndex AddPoint (const Point3d & p, bool aisghost, int layer, POINTTYPE type); #endif - int GetNP () const { return points.Size(); } + int GetNP () const { return points.Size(); } - MeshPoint & Point(int i) { return points.Elem(i); } - MeshPoint & Point(PointIndex pi) { return points[pi]; } - const MeshPoint & Point(int i) const { return points.Get(i); } - const MeshPoint & Point(PointIndex pi) const { return points[pi]; } + MeshPoint & Point(int i) { return points.Elem(i); } + MeshPoint & Point(PointIndex pi) { return points[pi]; } + const MeshPoint & Point(int i) const { return points.Get(i); } + const MeshPoint & Point(PointIndex pi) const { return points[pi]; } - const MeshPoint & operator[] (PointIndex pi) const { return points[pi]; } - MeshPoint & operator[] (PointIndex pi) { return points[pi]; } + const MeshPoint & operator[] (PointIndex pi) const { return points[pi]; } + MeshPoint & operator[] (PointIndex pi) { return points[pi]; } - const T_POINTS & Points() const { return points; } - T_POINTS & Points() { return points; } + const T_POINTS & Points() const { return points; } + T_POINTS & Points() { return points; } - SegmentIndex AddSegment (const Segment & s); - void DeleteSegment (int segnr) - { - segments.Elem(segnr)[0] = PointIndex::BASE-1; - segments.Elem(segnr)[1] = PointIndex::BASE-1; - } - void FullDeleteSegment (int segnr) // von wem ist das ??? - { - segments.Delete(segnr-PointIndex::BASE); - } + SegmentIndex AddSegment (const Segment & s); + void DeleteSegment (int segnr) + { + segments.Elem(segnr)[0] = PointIndex::BASE-1; + segments.Elem(segnr)[1] = PointIndex::BASE-1; + } + void FullDeleteSegment (int segnr) // von wem ist das ??? + { + segments.Delete(segnr-PointIndex::BASE); + } - int GetNSeg () const { return segments.Size(); } - Segment & LineSegment(int i) { return segments.Elem(i); } - const Segment & LineSegment(int i) const { return segments.Get(i); } + int GetNSeg () const { return segments.Size(); } + Segment & LineSegment(int i) { return segments.Elem(i); } + const Segment & LineSegment(int i) const { return segments.Get(i); } - Segment & LineSegment(SegmentIndex si) { return segments[si]; } - const Segment & LineSegment(SegmentIndex si) const { return segments[si]; } - const Segment & operator[] (SegmentIndex si) const { return segments[si]; } - Segment & operator[] (SegmentIndex si) { return segments[si]; } + Segment & LineSegment(SegmentIndex si) { return segments[si]; } + const Segment & LineSegment(SegmentIndex si) const { return segments[si]; } + const Segment & operator[] (SegmentIndex si) const { return segments[si]; } + Segment & operator[] (SegmentIndex si) { return segments[si]; } - SurfaceElementIndex AddSurfaceElement (const Element2d & el); - void DeleteSurfaceElement (int eli) - { - surfelements.Elem(eli).Delete(); - surfelements.Elem(eli).PNum(1) = -1; - surfelements.Elem(eli).PNum(2) = -1; - surfelements.Elem(eli).PNum(3) = -1; - timestamp = NextTimeStamp(); - } + SurfaceElementIndex AddSurfaceElement (const Element2d & el); + void DeleteSurfaceElement (int eli) + { + surfelements.Elem(eli).Delete(); + surfelements.Elem(eli).PNum(1) = -1; + surfelements.Elem(eli).PNum(2) = -1; + surfelements.Elem(eli).PNum(3) = -1; + timestamp = NextTimeStamp(); + } - void DeleteSurfaceElement (SurfaceElementIndex eli) - { - DeleteSurfaceElement (int(eli)+1); - } + void DeleteSurfaceElement (SurfaceElementIndex eli) + { + DeleteSurfaceElement (int(eli)+1); + } - int GetNSE () const { return surfelements.Size(); } - Element2d & SurfaceElement(int i) { return surfelements.Elem(i); } - const Element2d & SurfaceElement(int i) const { return surfelements.Get(i); } - Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; } - const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; } + int GetNSE () const { return surfelements.Size(); } + Element2d & SurfaceElement(int i) { return surfelements.Elem(i); } + const Element2d & SurfaceElement(int i) const { return surfelements.Get(i); } + Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; } + const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; } - const Element2d & operator[] (SurfaceElementIndex ei) const - { return surfelements[ei]; } - Element2d & operator[] (SurfaceElementIndex ei) - { return surfelements[ei]; } + const Element2d & operator[] (SurfaceElementIndex ei) const + { return surfelements[ei]; } + Element2d & operator[] (SurfaceElementIndex ei) + { return surfelements[ei]; } - void RebuildSurfaceElementLists (); - void GetSurfaceElementsOfFace (int facenr, Array & sei) const; + void RebuildSurfaceElementLists (); + void GetSurfaceElementsOfFace (int facenr, Array & sei) const; - ElementIndex AddVolumeElement (const Element & el); + ElementIndex AddVolumeElement (const Element & el); - int GetNE () const { return volelements.Size(); } + int GetNE () const { return volelements.Size(); } - Element & VolumeElement(int i) { return volelements.Elem(i); } - const Element & VolumeElement(int i) const { return volelements.Get(i); } - Element & VolumeElement(ElementIndex i) { return volelements[i]; } - const Element & VolumeElement(ElementIndex i) const { return volelements[i]; } + Element & VolumeElement(int i) { return volelements.Elem(i); } + const Element & VolumeElement(int i) const { return volelements.Get(i); } + Element & VolumeElement(ElementIndex i) { return volelements[i]; } + const Element & VolumeElement(ElementIndex i) const { return volelements[i]; } - const Element & operator[] (ElementIndex ei) const - { return volelements[ei]; } - Element & operator[] (ElementIndex ei) - { return volelements[ei]; } + const Element & operator[] (ElementIndex ei) const + { return volelements[ei]; } + Element & operator[] (ElementIndex ei) + { return volelements[ei]; } - ELEMENTTYPE ElementType (ElementIndex i) const - { return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; } + ELEMENTTYPE ElementType (ElementIndex i) const + { return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; } - const T_VOLELEMENTS & VolumeElements() const { return volelements; } - T_VOLELEMENTS & VolumeElements() { return volelements; } + const T_VOLELEMENTS & VolumeElements() const { return volelements; } + T_VOLELEMENTS & VolumeElements() { return volelements; } - /// - double ElementError (int eli) const; + /// + double ElementError (int eli) const; - /// - void AddLockedPoint (PointIndex pi); - /// - void ClearLockedPoints (); + /// + void AddLockedPoint (PointIndex pi); + /// + void ClearLockedPoints (); - const Array & LockedPoints() const - { return lockedpoints; } + const Array & LockedPoints() const + { return lockedpoints; } - /// Returns number of domains - int GetNDomains() const; + /// Returns number of domains + int GetNDomains() const; - /// - int GetDimension() const - { return dimension; } - void SetDimension(int dim) - { dimension = dim; } + /// + int GetDimension() const + { return dimension; } + void SetDimension(int dim) + { dimension = dim; } - /// sets internal tables - void CalcSurfacesOfNode (); + /// sets internal tables + void CalcSurfacesOfNode (); - /// additional (temporarily) fix points - void FixPoints (const BitArray & fixpoints); + /// additional (temporarily) fix points + void FixPoints (const BitArray & fixpoints); - /** - finds elements without neighbour and - boundary elements without inner element. - Results are stored in openelements. - if dom == 0, all sub-domains, else subdomain dom */ - void FindOpenElements (int dom = 0); + /** + finds elements without neighbour and + boundary elements without inner element. + Results are stored in openelements. + if dom == 0, all sub-domains, else subdomain dom */ + void FindOpenElements (int dom = 0); - /** - finds segments without surface element, - and surface elements without neighbours. - store in opensegmentsy - */ - void FindOpenSegments (int surfnr = 0); - /** - remove one layer of surface elements - */ - void RemoveOneLayerSurfaceElements (); + /** + finds segments without surface element, + and surface elements without neighbours. + store in opensegmentsy + */ + void FindOpenSegments (int surfnr = 0); + /** + remove one layer of surface elements + */ + void RemoveOneLayerSurfaceElements (); - int GetNOpenSegments () { return opensegments.Size(); } - const Segment & GetOpenSegment (int nr) { return opensegments.Get(nr); } + int GetNOpenSegments () { return opensegments.Size(); } + const Segment & GetOpenSegment (int nr) { return opensegments.Get(nr); } - /** - Checks overlap of boundary - return == 1, iff overlap - */ - int CheckOverlappingBoundary (); - /** - Checks consistent boundary - return == 0, everything ok - */ - int CheckConsistentBoundary () const; + /** + Checks overlap of boundary + return == 1, iff overlap + */ + int CheckOverlappingBoundary (); + /** + Checks consistent boundary + return == 0, everything ok + */ + int CheckConsistentBoundary () const; - /* - checks element orientation - */ - int CheckVolumeMesh () const; + /* + checks element orientation + */ + int CheckVolumeMesh () const; - /** - finds average h of surface surfnr if surfnr > 0, - else of all surfaces. - */ - double AverageH (int surfnr = 0) const; - /// Calculates localh - void CalcLocalH (); - /// - void SetLocalH (const Point3d & pmin, const Point3d & pmax, double grading); - /// - void RestrictLocalH (const Point3d & p, double hloc); - /// - void RestrictLocalHLine (const Point3d & p1, const Point3d & p2, - double hloc); - /// number of elements per radius - void CalcLocalHFromSurfaceCurvature(double elperr); - /// - void CalcLocalHFromPointDistances(void); - /// - void RestrictLocalH (resthtype rht, int nr, double loch); - /// - void LoadLocalMeshSize (const char * meshsizefilename); - /// - void SetGlobalH (double h); - /// - void SetMinimalH (double h); - /// - double MaxHDomain (int dom) const; - /// - void SetMaxHDomain (const Array & mhd); - /// - double GetH (const Point3d & p) const; - /// - double GetMinH (const Point3d & pmin, const Point3d & pmax); - /// - LocalH & LocalHFunction () { return * lochfunc; } - /// - bool LocalHFunctionGenerated(void) const { return (lochfunc != NULL); } + /** + finds average h of surface surfnr if surfnr > 0, + else of all surfaces. + */ + double AverageH (int surfnr = 0) const; + /// Calculates localh + void CalcLocalH (); + /// + void SetLocalH (const Point3d & pmin, const Point3d & pmax, double grading); + /// + void RestrictLocalH (const Point3d & p, double hloc); + /// + void RestrictLocalHLine (const Point3d & p1, const Point3d & p2, + double hloc); + /// number of elements per radius + void CalcLocalHFromSurfaceCurvature(double elperr); + /// + void CalcLocalHFromPointDistances(void); + /// + void RestrictLocalH (resthtype rht, int nr, double loch); + /// + void LoadLocalMeshSize (const char * meshsizefilename); + /// + void SetGlobalH (double h); + /// + void SetMinimalH (double h); + /// + double MaxHDomain (int dom) const; + /// + void SetMaxHDomain (const Array & mhd); + /// + double GetH (const Point3d & p) const; + /// + double GetMinH (const Point3d & pmin, const Point3d & pmax); + /// + LocalH & LocalHFunction () { return * lochfunc; } + /// + bool LocalHFunctionGenerated(void) const { return (lochfunc != NULL); } - /// Find bounding box - void GetBox (Point3d & pmin, Point3d & pmax, int dom = -1) const; + /// Find bounding box + void GetBox (Point3d & pmin, Point3d & pmax, int dom = -1) const; - /// Find bounding box of points of typ ptyp or less - void GetBox (Point3d & pmin, Point3d & pmax, POINTTYPE ptyp ) const; + /// Find bounding box of points of typ ptyp or less + void GetBox (Point3d & pmin, Point3d & pmax, POINTTYPE ptyp ) const; - /// - int GetNOpenElements() const - { return openelements.Size(); } - /// - const Element2d & OpenElement(int i) const - { return openelements.Get(i); } + /// + int GetNOpenElements() const + { return openelements.Size(); } + /// + const Element2d & OpenElement(int i) const + { return openelements.Get(i); } - /// are also quads open elements - bool HasOpenQuads () const; + /// are also quads open elements + bool HasOpenQuads () const; - /// split into connected pieces - void SplitIntoParts (); + /// split into connected pieces + void SplitIntoParts (); - /// - void SplitSeparatedFaces (); + /// + void SplitSeparatedFaces (); - /// Refines mesh and projects points to true surface - // void Refine (int levels, const CSGeometry * geom); + /// Refines mesh and projects points to true surface + // void Refine (int levels, const CSGeometry * geom); - bool BoundaryEdge (PointIndex pi1, PointIndex pi2) const - { - if(!boundaryedges) - const_cast(this)->BuildBoundaryEdges(); + bool BoundaryEdge (PointIndex pi1, PointIndex pi2) const + { + if(!boundaryedges) + const_cast(this)->BuildBoundaryEdges(); - INDEX_2 i2 (pi1, pi2); - i2.Sort(); - return boundaryedges->Used (i2); - } + INDEX_2 i2 (pi1, pi2); + i2.Sort(); + return boundaryedges->Used (i2); + } - bool IsSegment (PointIndex pi1, PointIndex pi2) const - { - INDEX_2 i2 (pi1, pi2); - i2.Sort(); - return segmentht->Used (i2); - } + bool IsSegment (PointIndex pi1, PointIndex pi2) const + { + INDEX_2 i2 (pi1, pi2); + i2.Sort(); + return segmentht->Used (i2); + } - SegmentIndex SegmentNr (PointIndex pi1, PointIndex pi2) const - { - INDEX_2 i2 (pi1, pi2); - i2.Sort(); - return segmentht->Get (i2); - } + SegmentIndex SegmentNr (PointIndex pi1, PointIndex pi2) const + { + INDEX_2 i2 (pi1, pi2); + i2.Sort(); + return segmentht->Get (i2); + } - /** - Remove unused points. etc. - */ - void Compress (); + /** + Remove unused points. etc. + */ + void Compress (); - /// - void Save (ostream & outfile) const; - /// - void Load (istream & infile); - /// - void Merge (istream & infile, const int surfindex_offset = 0); - /// - void Save (const string & filename) const; - /// - void Load (const string & filename); - /// - void Merge (const string & filename, const int surfindex_offset = 0); + /// + void Save (ostream & outfile) const; + /// + void Load (istream & infile); + /// + void Merge (istream & infile, const int surfindex_offset = 0); + /// + void Save (const string & filename) const; + /// + void Load (const string & filename); + /// + void Merge (const string & filename, const int surfindex_offset = 0); - /// - void ImproveMesh (OPTIMIZEGOAL goal = OPT_QUALITY); + /// + void ImproveMesh (OPTIMIZEGOAL goal = OPT_QUALITY); - /// - void ImproveMeshJacobian (OPTIMIZEGOAL goal = OPT_QUALITY, const BitArray * usepoint = NULL); - /// - void ImproveMeshJacobianOnSurface (const BitArray & usepoint, - const Array< Vec<3>* > & nv, - OPTIMIZEGOAL goal = OPT_QUALITY, - const Array< Array* > * idmaps = NULL); - /** - free nodes in environment of openelements - for optimiztion - */ - void FreeOpenElementsEnvironment (int layers); + /// + void ImproveMeshJacobian (OPTIMIZEGOAL goal = OPT_QUALITY, const BitArray * usepoint = NULL); + /// + void ImproveMeshJacobianOnSurface (const BitArray & usepoint, + const Array< Vec<3>* > & nv, + OPTIMIZEGOAL goal = OPT_QUALITY, + const Array< Array* > * idmaps = NULL); + /** + free nodes in environment of openelements + for optimiztion + */ + void FreeOpenElementsEnvironment (int layers); - /// - bool LegalTet (Element & el) const - { - if (el.IllegalValid()) - return !el.Illegal(); - return LegalTet2 (el); - } - /// - bool LegalTet2 (Element & el) const; + /// + bool LegalTet (Element & el) const + { + if (el.IllegalValid()) + return !el.Illegal(); + return LegalTet2 (el); + } + /// + bool LegalTet2 (Element & el) const; - /// - bool LegalTrig (const Element2d & el) const; - /** - if values non-null, return values in 4-double array: - triangle angles min/max, tetangles min/max - if null, output results on cout - */ - void CalcMinMaxAngle (double badellimit, double * retvalues = NULL); + /// + bool LegalTrig (const Element2d & el) const; + /** + if values non-null, return values in 4-double array: + triangle angles min/max, tetangles min/max + if null, output results on cout + */ + void CalcMinMaxAngle (double badellimit, double * retvalues = NULL); - /* - Marks elements which are dangerous to refine - return: number of illegal elements - */ - int MarkIllegalElements (); + /* + Marks elements which are dangerous to refine + return: number of illegal elements + */ + int MarkIllegalElements (); - /// orient surface mesh, for one sub-domain only - void SurfaceMeshOrientation (); + /// orient surface mesh, for one sub-domain only + void SurfaceMeshOrientation (); - /// convert mixed element mesh to tet-mesh - void Split2Tets(); + /// convert mixed element mesh to tet-mesh + void Split2Tets(); - /// build box-search tree - void BuildElementSearchTree (); + /// build box-search tree + void BuildElementSearchTree (); - void SetPointSearchStartElement(const int el) const {ps_startelement = el;} + void SetPointSearchStartElement(const int el) const {ps_startelement = el;} - /// gives element of point, barycentric coordinates - int GetElementOfPoint (const Point3d & p, - double * lami, - bool build_searchtree = 0, - const int index = -1, - const bool allowindex = true) const; - int GetElementOfPoint (const Point3d & p, - double * lami, - const Array * const indices, - bool build_searchtree = 0, - const bool allowindex = true) const; - int GetSurfaceElementOfPoint (const Point3d & p, - double * lami, - bool build_searchtree = 0, - const int index = -1, - const bool allowindex = true) const; - int GetSurfaceElementOfPoint (const Point3d & p, - double * lami, - const Array * const indices, - bool build_searchtree = 0, - const bool allowindex = true) const; + /// gives element of point, barycentric coordinates + int GetElementOfPoint (const Point3d & p, + double * lami, + bool build_searchtree = 0, + const int index = -1, + const bool allowindex = true) const; + int GetElementOfPoint (const Point3d & p, + double * lami, + const Array * const indices, + bool build_searchtree = 0, + const bool allowindex = true) const; + int GetSurfaceElementOfPoint (const Point3d & p, + double * lami, + bool build_searchtree = 0, + const int index = -1, + const bool allowindex = true) const; + int GetSurfaceElementOfPoint (const Point3d & p, + double * lami, + const Array * const indices, + bool build_searchtree = 0, + const bool allowindex = true) const; - /// give list of vol elements which are int the box(p1,p2) - void GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, - Array & locels) const; + /// give list of vol elements which are int the box(p1,p2) + void GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, + Array & locels) const; - /// - int AddFaceDescriptor(const FaceDescriptor& fd) - { return facedecoding.Append(fd); } + /// + int AddFaceDescriptor(const FaceDescriptor& fd) + { return facedecoding.Append(fd); } - int AddEdgeDescriptor(const EdgeDescriptor & fd) - { return edgedecoding.Append(fd) - 1; } + int AddEdgeDescriptor(const EdgeDescriptor & fd) + { return edgedecoding.Append(fd) - 1; } - /// - void SetMaterial (int domnr, const char * mat); - /// - const char * GetMaterial (int domnr) const; + /// + void SetMaterial (int domnr, const char * mat); + /// + const char * GetMaterial (int domnr) const; - void SetNBCNames ( int nbcn ); + void SetNBCNames ( int nbcn ); - void SetBCName ( int bcnr, const string & abcname ); + void SetBCName ( int bcnr, const string & abcname ); - string GetBCName ( int bcnr ) const; + string GetBCName ( int bcnr ) const; - string * GetBCNamePtr ( int bcnr ) - { return bcnames[bcnr]; } + string * GetBCNamePtr ( int bcnr ) + { return bcnames[bcnr]; } - /// - void ClearFaceDescriptors() - { facedecoding.SetSize(0); } + /// + void ClearFaceDescriptors() + { facedecoding.SetSize(0); } - /// - int GetNFD () const - { return facedecoding.Size(); } + /// + int GetNFD () const + { return facedecoding.Size(); } - const FaceDescriptor & GetFaceDescriptor (int i) const - { return facedecoding.Get(i); } + const FaceDescriptor & GetFaceDescriptor (int i) const + { return facedecoding.Get(i); } - const EdgeDescriptor & GetEdgeDescriptor (int i) const - { return edgedecoding[i]; } + const EdgeDescriptor & GetEdgeDescriptor (int i) const + { return edgedecoding[i]; } - /// - FaceDescriptor & GetFaceDescriptor (int i) - { return facedecoding.Elem(i); } + /// + FaceDescriptor & GetFaceDescriptor (int i) + { return facedecoding.Elem(i); } -// #ifdef NONE -// /* -// Identify points pi1 and pi2, due to -// identification nr identnr -// */ -// void AddIdentification (int pi1, int pi2, int identnr); + // #ifdef NONE + // /* + // Identify points pi1 and pi2, due to + // identification nr identnr + // */ + // void AddIdentification (int pi1, int pi2, int identnr); -// int GetIdentification (int pi1, int pi2) const; -// int GetIdentificationSym (int pi1, int pi2) const; -// /// -// INDEX_2_HASHTABLE & GetIdentifiedPoints () -// { -// return *identifiedpoints; -// } + // int GetIdentification (int pi1, int pi2) const; + // int GetIdentificationSym (int pi1, int pi2) const; + // /// + // INDEX_2_HASHTABLE & GetIdentifiedPoints () + // { + // return *identifiedpoints; + // } -// /// -// void GetIdentificationMap (int identnr, Array & identmap) const; -// /// -// void GetIdentificationPairs (int identnr, Array & identpairs) const; -// /// -// int GetMaxIdentificationNr () const -// { -// return maxidentnr; -// } -// #endif + // /// + // void GetIdentificationMap (int identnr, Array & identmap) const; + // /// + // void GetIdentificationPairs (int identnr, Array & identpairs) const; + // /// + // int GetMaxIdentificationNr () const + // { + // return maxidentnr; + // } + // #endif - /// return periodic, close surface etc. identifications - Identifications & GetIdentifications () { return *ident; } - /// return periodic, close surface etc. identifications - const Identifications & GetIdentifications () const { return *ident; } + /// return periodic, close surface etc. identifications + Identifications & GetIdentifications () { return *ident; } + /// return periodic, close surface etc. identifications + const Identifications & GetIdentifications () const { return *ident; } - void InitPointCurve(double red = 1, double green = 0, double blue = 0) const; - void AddPointCurvePoint(const Point3d & pt) const; - int GetNumPointCurves(void) const; - int GetNumPointsOfPointCurve(int curve) const; - Point3d & GetPointCurvePoint(int curve, int n) const; - void GetPointCurveColor(int curve, double & red, double & green, double & blue) const; + void InitPointCurve(double red = 1, double green = 0, double blue = 0) const; + void AddPointCurvePoint(const Point3d & pt) const; + int GetNumPointCurves(void) const; + int GetNumPointsOfPointCurve(int curve) const; + Point3d & GetPointCurvePoint(int curve, int n) const; + void GetPointCurveColor(int curve, double & red, double & green, double & blue) const; - /// find number of vertices - void ComputeNVertices (); - /// number of vertices (no edge-midpoints) - int GetNV () const; - /// remove edge points - void SetNP (int np); + /// find number of vertices + void ComputeNVertices (); + /// number of vertices (no edge-midpoints) + int GetNV () const; + /// remove edge points + void SetNP (int np); - /* - /// build connected nodes along prism stack - void BuildConnectedNodes (); - void ConnectToNodeRec (int node, int tonode, - const TABLE & conto); - */ + /* + /// build connected nodes along prism stack + void BuildConnectedNodes (); + void ConnectToNodeRec (int node, int tonode, + const TABLE & conto); + */ - bool PureTrigMesh (int faceindex = 0) const; - bool PureTetMesh () const; + bool PureTrigMesh (int faceindex = 0) const; + bool PureTetMesh () const; - const class MeshTopology & GetTopology () const - { return *topology; } + const class MeshTopology & GetTopology () const + { return *topology; } - void UpdateTopology(); + void UpdateTopology(); - class CurvedElements & GetCurvedElements () const - { return *curvedelems; } + class CurvedElements & GetCurvedElements () const + { return *curvedelems; } - const class AnisotropicClusters & GetClusters () const - { return *clusters; } + const class AnisotropicClusters & GetClusters () const + { return *clusters; } - int GetTimeStamp() const { return timestamp; } - void SetNextTimeStamp() - { timestamp = NextTimeStamp(); } + int GetTimeStamp() const { return timestamp; } + void SetNextTimeStamp() + { timestamp = NextTimeStamp(); } - int GetMajorTimeStamp() const { return majortimestamp; } - void SetNextMajorTimeStamp() - { majortimestamp = timestamp = NextTimeStamp(); } + int GetMajorTimeStamp() const { return majortimestamp; } + void SetNextMajorTimeStamp() + { majortimestamp = timestamp = NextTimeStamp(); } - /// return mutex - NgMutex & Mutex () { return mutex; } - NgMutex & MajorMutex () { return majormutex; } + /// return mutex + NgMutex & Mutex () { return mutex; } + NgMutex & MajorMutex () { return majormutex; } - /// - void SetUserData(const char * id, Array & data); - /// - bool GetUserData(const char * id, Array & data, int shift = 0) const; - /// - void SetUserData(const char * id, Array & data); - /// - bool GetUserData(const char * id, Array & data, int shift = 0) const; + /// + void SetUserData(const char * id, Array & data); + /// + bool GetUserData(const char * id, Array & data, int shift = 0) const; + /// + void SetUserData(const char * id, Array & data); + /// + bool GetUserData(const char * id, Array & data, int shift = 0) const; - /// - friend void OptimizeRestart (Mesh & mesh3d); - /// - void PrintMemInfo (ostream & ost) const; - /// - friend class Meshing3; + /// + friend void OptimizeRestart (Mesh & mesh3d); + /// + void PrintMemInfo (ostream & ost) const; + /// + friend class Meshing3; - enum GEOM_TYPE { NO_GEOM = 0, GEOM_2D = 1, GEOM_CSG = 10, GEOM_STL = 11, GEOM_OCC = 12, GEOM_ACIS = 13 }; - GEOM_TYPE geomtype; + enum GEOM_TYPE { NO_GEOM = 0, GEOM_2D = 1, GEOM_CSG = 10, GEOM_STL = 11, GEOM_OCC = 12, GEOM_ACIS = 13 }; + GEOM_TYPE geomtype; #ifdef PARALLEL - /// returns parallel topology - class ParallelMeshTopology & GetParallelTopology () const - { return *paralleltop; } + /// returns parallel topology + class ParallelMeshTopology & GetParallelTopology () const + { return *paralleltop; } - /// distributes the master-mesh to local meshes - void Distribute (); + /// distributes the master-mesh to local meshes + void Distribute (); - /// loads a mesh sent from master processor - void ReceiveParallelMesh (); + /// loads a mesh sent from master processor + void ReceiveParallelMesh (); - /// find connection to parallel meshes -// void FindExchangePoints () ; + /// find connection to parallel meshes + // void FindExchangePoints () ; -// void FindExchangeEdges (); -// void FindExchangeFaces (); + // void FindExchangeEdges (); + // void FindExchangeFaces (); - /// use metis to decompose master mesh - void ParallelMetis (); // Array & neloc ); - void PartHybridMesh (); // Array & neloc ); - void PartDualHybridMesh (); // Array & neloc ); - void PartDualHybridMesh2D (); // ( Array & neloc ); + /// use metis to decompose master mesh + void ParallelMetis (); // Array & neloc ); + void PartHybridMesh (); // Array & neloc ); + void PartDualHybridMesh (); // Array & neloc ); + void PartDualHybridMesh2D (); // ( Array & neloc ); - /// send mesh to parallel machine, keep global mesh at master - void SendMesh ( ) const; // Mesh * mastermesh, Array & neloc) const; + /// send mesh to parallel machine, keep global mesh at master + void SendMesh ( ) const; // Mesh * mastermesh, Array & neloc) const; - void UpdateOverlap (); + void UpdateOverlap (); #endif -}; + }; + + inline ostream& operator<<(ostream& ost, const Mesh& mesh) + { + ost << "mesh: " << endl; + mesh.Save(ost); + return ost; + } -inline ostream& operator<<(ostream& ost, const Mesh& mesh) -{ - ost << "mesh: " << endl; - mesh.Save(ost); - return ost; } - #endif diff --git a/libsrc/meshing/meshing.hpp b/libsrc/meshing/meshing.hpp index 81daa843..7fd9c196 100644 --- a/libsrc/meshing/meshing.hpp +++ b/libsrc/meshing/meshing.hpp @@ -12,20 +12,22 @@ namespace netgen { - - extern int printmessage_importance; + // extern int printmessage_importance; class CSGeometry; + class NetgenGeometry; +} #include "msghandler.hpp" - #include "meshtype.hpp" #include "localh.hpp" #include "meshclass.hpp" #include "global.hpp" +namespace netgen +{ #include "meshtool.hpp" #include "ruler2.hpp" #include "adfront2.hpp" @@ -37,17 +39,12 @@ namespace netgen #include "adfront3.hpp" #include "ruler3.hpp" -#ifndef SMALLLIB #define _INCLUDE_MORE -#endif -#ifdef LINUX -#define _INCLUDE_MORE -#endif -#ifdef _INCLUDE_MORE + #include "meshing3.hpp" #include "improve3.hpp" -#endif + #include "findip.hpp" #include "findip2.hpp" @@ -55,23 +52,20 @@ namespace netgen #include "curvedelems.hpp" #include "clusters.hpp" -#ifdef _INCLUDE_MORE #include "meshfunc.hpp" -#endif + #include "bisect.hpp" #include "hprefinement.hpp" #include "boundarylayer.hpp" #include "specials.hpp" +} #include "validate.hpp" +#include "basegeom.hpp" #ifdef PARALLEL #include "paralleltop.hpp" -// #include "../parallel/parallelmesh.hpp" #endif -} -#include "basegeom.hpp" - #endif diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 7d3e695b..ed7b6c78 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -8,663 +8,667 @@ /* Date: 01. Okt. 95 */ /**************************************************************************/ -/* +namespace netgen +{ + + /* Classes for NETGEN -*/ + */ -enum ELEMENT_TYPE { - SEGMENT = 1, SEGMENT3 = 2, - TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14, - TET = 20, TET10 = 21, - PYRAMID = 22, PRISM = 23, PRISM12 = 24, - HEX = 25 -}; -typedef int ELEMENT_EDGE[2]; // initial point, end point -typedef int ELEMENT_FACE[4]; // points, last one is -1 for trig + enum ELEMENT_TYPE { + SEGMENT = 1, SEGMENT3 = 2, + TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14, + TET = 20, TET10 = 21, + PYRAMID = 22, PRISM = 23, PRISM12 = 24, + HEX = 25 + }; + + typedef int ELEMENT_EDGE[2]; // initial point, end point + typedef int ELEMENT_FACE[4]; // points, last one is -1 for trig #define ELEMENT_MAXPOINTS 12 #define ELEMENT2D_MAXPOINTS 8 -enum POINTTYPE { FIXEDPOINT = 1, EDGEPOINT = 2, SURFACEPOINT = 3, INNERPOINT = 4 }; -enum ELEMENTTYPE { FREEELEMENT, FIXEDELEMENT }; -enum OPTIMIZEGOAL { OPT_QUALITY, OPT_CONFORM, OPT_REST, OPT_WORSTCASE, OPT_LEGAL }; + enum POINTTYPE { FIXEDPOINT = 1, EDGEPOINT = 2, SURFACEPOINT = 3, INNERPOINT = 4 }; + enum ELEMENTTYPE { FREEELEMENT, FIXEDELEMENT }; + enum OPTIMIZEGOAL { OPT_QUALITY, OPT_CONFORM, OPT_REST, OPT_WORSTCASE, OPT_LEGAL }; -extern int GetTimeStamp(); -extern int NextTimeStamp(); + extern int GetTimeStamp(); + extern int NextTimeStamp(); -class PointGeomInfo -{ -public: - int trignum; // for STL Meshing - double u, v; // for OCC Meshing + class PointGeomInfo + { + public: + int trignum; // for STL Meshing + double u, v; // for OCC Meshing - PointGeomInfo () - : trignum(-1), u(0), v(0) { ; } -}; + PointGeomInfo () + : trignum(-1), u(0), v(0) { ; } + }; -inline ostream & operator<< (ostream & ost, const PointGeomInfo & gi) -{ - return (ost << gi.trignum << " " << gi.u << " " << gi.v); -} + inline ostream & operator<< (ostream & ost, const PointGeomInfo & gi) + { + return (ost << gi.trignum << " " << gi.u << " " << gi.v); + } -inline istream & operator>> (istream & ist, PointGeomInfo & gi) -{ - return (ist >> gi.trignum >> gi.u >> gi.v); -} + inline istream & operator>> (istream & ist, PointGeomInfo & gi) + { + return (ist >> gi.trignum >> gi.u >> gi.v); + } #define MULTIPOINTGEOMINFO_MAX 100 -class MultiPointGeomInfo -{ - int cnt; - PointGeomInfo mgi[MULTIPOINTGEOMINFO_MAX]; -public: - MultiPointGeomInfo () { cnt = 0; } - int AddPointGeomInfo (const PointGeomInfo & gi); - void Init () { cnt = 0; } - void DeleteAll () { cnt = 0; } - - int GetNPGI () const { return cnt; } - const PointGeomInfo & GetPGI (int i) const { return mgi[i-1]; } -}; - - -class EdgePointGeomInfo -{ -public: - int edgenr; - int body; // for ACIS - double dist; // for 2d meshing - double u, v; // for OCC Meshing - -public: - EdgePointGeomInfo () - : edgenr(0), body(0), dist(0.0), u(0.0), v(0.0) { ; } - - - EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) + class MultiPointGeomInfo { - edgenr = gi2.edgenr; - body = gi2.body; - dist = gi2.dist; - u = gi2.u; v = gi2.v; - return *this; - } -}; - -inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi) -{ - ost << "epgi: edgnr=" << gi.edgenr << ", dist=" << gi.dist; - return ost; -} - - - - - -class PointIndex -{ - int i; -public: - PointIndex () { ; } - PointIndex (int ai) : i(ai) { ; } - PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } - PointIndex & operator= (int ai) { i = ai; return *this; } - operator int () const { return i; } - int GetInt () const { return i; } - PointIndex operator++ (int) { int hi = i; i++; return PointIndex(hi); } - PointIndex operator-- (int) { int hi = i; i--; return PointIndex(hi); } - -#ifdef BASE0 - enum { BASE = 0 }; -#else - enum { BASE = 1 }; -#endif -}; - -inline istream & operator>> (istream & ist, PointIndex & pi) -{ - int i; ist >> i; pi = i; return ist; -} - -inline ostream & operator<< (ostream & ost, const PointIndex & pi) -{ - return (ost << pi.GetInt()); -} - - - - -class ElementIndex -{ - int i; -public: - ElementIndex () { ; } - ElementIndex (int ai) : i(ai) { ; } - ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; } - ElementIndex & operator= (int ai) { i = ai; return *this; } - operator int () const { return i; } - ElementIndex & operator++ (int) { i++; return *this; } - ElementIndex & operator-- (int) { i--; return *this; } -}; - -inline istream & operator>> (istream & ist, ElementIndex & pi) -{ - int i; ist >> i; pi = i; return ist; -} - -inline ostream & operator<< (ostream & ost, const ElementIndex & pi) -{ - return (ost << int(pi)); -} - - -class SurfaceElementIndex -{ - int i; -public: - SurfaceElementIndex () { ; } - SurfaceElementIndex (int ai) : i(ai) { ; } - SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) - { i = ai.i; return *this; } - SurfaceElementIndex & operator= (int ai) { i = ai; return *this; } - operator int () const { return i; } - SurfaceElementIndex & operator++ (int) { i++; return *this; } - SurfaceElementIndex & operator-- (int) { i--; return *this; } -}; - -inline istream & operator>> (istream & ist, SurfaceElementIndex & pi) -{ - int i; ist >> i; pi = i; return ist; -} - -inline ostream & operator<< (ostream & ost, const SurfaceElementIndex & pi) -{ - return (ost << int(pi)); -} - -class SegmentIndex -{ - int i; -public: - SegmentIndex () { ; } - SegmentIndex (int ai) : i(ai) { ; } - SegmentIndex & operator= (const SegmentIndex & ai) - { i = ai.i; return *this; } - SegmentIndex & operator= (int ai) { i = ai; return *this; } - operator int () const { return i; } - SegmentIndex & operator++ (int) { i++; return *this; } - SegmentIndex & operator-- (int) { i--; return *this; } -}; - -inline istream & operator>> (istream & ist, SegmentIndex & pi) -{ - int i; ist >> i; pi = i; return ist; -} - -inline ostream & operator<< (ostream & ost, const SegmentIndex & pi) -{ - return (ost << int(pi)); -} - - - - -/** - Point in the mesh. - Contains layer (a new feature in 4.3 for overlapping meshes. - */ -class MeshPoint : public Point<3> -{ - int layer; - double singular; // singular factor for hp-refinement - POINTTYPE type; - -#ifdef PARALLEL - bool isghost; -#endif - -public: - MeshPoint () - { -#ifdef PARALLEL - isghost = 0; -#endif - ; -} - - MeshPoint (const Point<3> & ap, int alayer = 1, POINTTYPE apt = INNERPOINT) - : Point<3> (ap), layer(alayer), singular(0.),type(apt) - { -#ifdef PARALLEL - isghost = 0; -#endif - ; - } - - void SetPoint (const Point<3> & ap) - { - Point<3>::operator= (ap); - layer = 0; - singular = 0; -#ifdef PARALLEL - isghost = 0; -#endif - } - - int GetLayer() const { return layer; } - - POINTTYPE Type() const { return type; } - void SetType(POINTTYPE at) { type = at; } - - double Singularity() const { return singular; } - void Singularity(double s) { singular = s; } - bool IsSingular() const { return (singular != 0.0); } - -#ifdef PARALLEL - bool IsGhost () const { return isghost; } - void SetGhost ( bool aisghost ) { isghost = aisghost; } -#endif - -}; - -inline ostream & operator<<(ostream & s, const MeshPoint & pt) -{ - return (s << Point<3> (pt)); -} - - - - -typedef Array T_POINTS; - - - -/** - Triangle element for surface mesh generation. - */ -class Element2d -{ - /// point numbers - PointIndex pnum[ELEMENT2D_MAXPOINTS]; - /// geom info of points - PointGeomInfo geominfo[ELEMENT2D_MAXPOINTS]; - - /// surface nr - int index:16; - /// - ELEMENT_TYPE typ:6; - /// number of points - unsigned int np:4; - bool badel:1; - bool refflag:1; // marked for refinement - bool strongrefflag:1; - bool deleted:1; // element is deleted - - // Philippose - 08 August 2010 - // Set a new property for each element, to - // control whether it is visible or not - bool visible:1; // element visible - - /// order for hp-FEM - unsigned int orderx:6; - unsigned int ordery:6; - -#ifdef PARALLEL - bool isghost; - int partitionNumber; -#endif - - /// a linked list for all segments in the same face - SurfaceElementIndex next; - -public: - /// - Element2d (); - /// - Element2d (int anp); - /// - Element2d (ELEMENT_TYPE type); - /// - Element2d (int pi1, int pi2, int pi3); - /// - Element2d (int pi1, int pi2, int pi3, int pi4); - /// - ELEMENT_TYPE GetType () const { return typ; } - /// - void SetType (ELEMENT_TYPE atyp) - { - typ = atyp; - switch (typ) - { - case TRIG: np = 3; break; - case QUAD: np = 4; break; - case TRIG6: np = 6; break; - case QUAD6: np = 6; break; - case QUAD8: np = 8; break; - default: - PrintSysError ("Element2d::SetType, illegal type ", typ); - } - } - /// - int GetNP() const { return np; } - /// - int GetNV() const - { - switch (typ) - { - case TRIG: - case TRIG6: return 3; - - case QUAD: - case QUAD8: - case QUAD6: return 4; - default: -#ifdef DEBUG - PrintSysError ("element2d::GetNV not implemented for typ", typ) -#endif - ; - } - return np; - } - - /// - PointIndex & operator[] (int i) { return pnum[i]; } - /// - const PointIndex & operator[] (int i) const { return pnum[i]; } - - /// - PointIndex & PNum (int i) { return pnum[i-1]; } - /// - const PointIndex & PNum (int i) const { return pnum[i-1]; } - /// - PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } - /// - const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } - /// - - /// - PointGeomInfo & GeomInfoPi (int i) { return geominfo[i-1]; } - /// - const PointGeomInfo & GeomInfoPi (int i) const { return geominfo[i-1]; } - /// - PointGeomInfo & GeomInfoPiMod (int i) { return geominfo[(i-1) % np]; } - /// - const PointGeomInfo & GeomInfoPiMod (int i) const { return geominfo[(i-1) % np]; } - - - void SetIndex (int si) { index = si; } - /// - int GetIndex () const { return index; } - - int GetOrder () const { return orderx; } - void SetOrder (int aorder) { orderx = ordery = aorder; } - - - void GetOrder (int & ox, int & oy) const { ox = orderx, oy =ordery;}; - void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz=0; } - void SetOrder (int ox, int oy, int /* oz */) { orderx = ox; ordery = oy;} - void SetOrder (int ox, int oy) { orderx = ox; ordery = oy;} - - - /// - void GetBox (const T_POINTS & points, Box3d & box) const; - /// invert orientation - inline void Invert (); - /// - void Invert2 (); - /// first point number is smallest - inline void NormalizeNumbering (); - /// - void NormalizeNumbering2 (); - - bool BadElement() const { return badel; } - - // friend ostream & operator<<(ostream & s, const Element2d & el); - friend class Mesh; - - - /// get number of 'integration points' - int GetNIP () const; - void GetIntegrationPoint (int ip, Point2d & p, double & weight) const; - - void GetTransformation (int ip, const Array & points, - class DenseMatrix & trans) const; - void GetTransformation (int ip, class DenseMatrix & pmat, - class DenseMatrix & trans) const; - - void GetShape (const Point2d & p, class Vector & shape) const; - void GetShapeNew (const Point<2> & p, class FlatVector & shape) const; - /// matrix 2 * np - void GetDShape (const Point2d & p, class DenseMatrix & dshape) const; - void GetDShapeNew (const Point<2> & p, class MatrixFixWidth<2> & dshape) const; - /// matrix 2 * np - void GetPointMatrix (const Array & points, - class DenseMatrix & pmat) const; - - void ComputeIntegrationPointData () const; - - - double CalcJacobianBadness (const Array & points) const; - double CalcJacobianBadness (const T_POINTS & points, - const Vec<3> & n) const; - double CalcJacobianBadnessDirDeriv (const Array & points, - int pi, Vec2d & dir, double & dd) const; - - - - void Delete () { deleted = 1; pnum[0] = pnum[1] = pnum[2] = pnum[3] = PointIndex::BASE-1; } - bool IsDeleted () const - { -#ifdef DEBUG - if (pnum[0] < PointIndex::BASE && !deleted) - cerr << "Surfelement has illegal pnum, but not marked as deleted" << endl; -#endif - return deleted; - } - - // Philippose - 08 August 2010 - // Access functions for the new property: visible - void Visible(bool vis = 1) - { visible = vis; } - bool IsVisible () const - { return visible; } - - void SetRefinementFlag (bool rflag = 1) - { refflag = rflag; } - bool TestRefinementFlag () const - { return refflag; } - - void SetStrongRefinementFlag (bool rflag = 1) - { strongrefflag = rflag; } - bool TestStrongRefinementFlag () const - { return strongrefflag; } - - - SurfaceElementIndex NextElement() { return next; } - - bool operator==(const Element2d & el2) const; - - int HasFace(const Element2d& el) const; - /// - int meshdocval; - /// - int hp_elnr; - -#ifdef PARALLEL - bool IsGhost () const { return isghost; } - void SetGhost ( bool aisghost ) { isghost = aisghost; } - - // by JS, only for 2D meshes ???? - int GetPartition () const { return partitionNumber; } - void SetPartition (int nr) { partitionNumber = nr; }; - -#else - bool IsGhost () const { return false; } -#endif -}; - - -ostream & operator<<(ostream & s, const Element2d & el); - - - - - -class IntegrationPointData -{ -public: - Point<3> p; - double weight; - Vector shape; - DenseMatrix dshape; -}; - - - - - - - - -/** - Volume element - */ -class Element -{ -private: - /// point numbers - PointIndex pnum[ELEMENT_MAXPOINTS]; - /// - ELEMENT_TYPE typ:6; - /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) - int np:5; - /// - class flagstruct { + int cnt; + PointGeomInfo mgi[MULTIPOINTGEOMINFO_MAX]; public: - bool marked:1; // marked for refinement - bool badel:1; // angles worse then limit - bool reverse:1; // for refinement a la Bey - bool illegal:1; // illegal, will be split or swaped - bool illegal_valid:1; // is illegal-flag valid ? - bool badness_valid:1; // is badness valid ? - bool refflag:1; // mark element for refinement - bool strongrefflag:1; - bool deleted:1; // element is deleted, will be removed from array - bool fixed:1; // don't change element in optimization + MultiPointGeomInfo () { cnt = 0; } + int AddPointGeomInfo (const PointGeomInfo & gi); + void Init () { cnt = 0; } + void DeleteAll () { cnt = 0; } + + int GetNPGI () const { return cnt; } + const PointGeomInfo & GetPGI (int i) const { return mgi[i-1]; } }; - /// sub-domain index - short int index; - /// order for hp-FEM - unsigned int orderx:6; - unsigned int ordery:6; - unsigned int orderz:6; - /* unsigned int levelx:6; - unsigned int levely:6; - unsigned int levelz:6; */ - /// stored shape-badness of element - float badness; - -#ifdef PARALLEL - /// number of partition for parallel computation - int partitionNumber; - bool isghost; -#endif -public: - flagstruct flags; - - /// - Element (); - /// - Element (int anp); - /// - Element (ELEMENT_TYPE type); - /// - Element & operator= (const Element & el2); - - /// - void SetNP (int anp); - /// - void SetType (ELEMENT_TYPE atyp); - /// - int GetNP () const { return np; } - /// - int GetNV() const + class EdgePointGeomInfo { - switch (typ) - { - case TET: - case TET10: - return 4; - case PRISM12: - case PRISM: - return 6; - case PYRAMID: - return 5; - case HEX: - return 8; - default: -#ifdef DEBUG - PrintSysError ("Element3d::GetNV not implemented for typ ", typ) -#endif - ; - } - return np; + public: + int edgenr; + int body; // for ACIS + double dist; // for 2d meshing + double u, v; // for OCC Meshing + + public: + EdgePointGeomInfo () + : edgenr(0), body(0), dist(0.0), u(0.0), v(0.0) { ; } + + + EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) + { + edgenr = gi2.edgenr; + body = gi2.body; + dist = gi2.dist; + u = gi2.u; v = gi2.v; + return *this; + } + }; + + inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi) + { + ost << "epgi: edgnr=" << gi.edgenr << ", dist=" << gi.dist; + return ost; } - bool operator==(const Element & el2) const; - // old style: - int NP () const { return np; } - /// - ELEMENT_TYPE GetType () const { return typ; } - /// - PointIndex & operator[] (int i) { return pnum[i]; } - /// - const PointIndex & operator[] (int i) const { return pnum[i]; } - /// - PointIndex & PNum (int i) { return pnum[i-1]; } - /// - const PointIndex & PNum (int i) const { return pnum[i-1]; } - /// - PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } - /// - const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } + class PointIndex + { + int i; + public: + PointIndex () { ; } + PointIndex (int ai) : i(ai) { ; } + PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } + PointIndex & operator= (int ai) { i = ai; return *this; } + operator int () const { return i; } + int GetInt () const { return i; } + PointIndex operator++ (int) { int hi = i; i++; return PointIndex(hi); } + PointIndex operator-- (int) { int hi = i; i--; return PointIndex(hi); } + +#ifdef BASE0 + enum { BASE = 0 }; +#else + enum { BASE = 1 }; +#endif + }; + + inline istream & operator>> (istream & ist, PointIndex & pi) + { + int i; ist >> i; pi = i; return ist; + } + + inline ostream & operator<< (ostream & ost, const PointIndex & pi) + { + return (ost << pi.GetInt()); + } + + + + + class ElementIndex + { + int i; + public: + ElementIndex () { ; } + ElementIndex (int ai) : i(ai) { ; } + ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; } + ElementIndex & operator= (int ai) { i = ai; return *this; } + operator int () const { return i; } + ElementIndex & operator++ (int) { i++; return *this; } + ElementIndex & operator-- (int) { i--; return *this; } + }; + + inline istream & operator>> (istream & ist, ElementIndex & pi) + { + int i; ist >> i; pi = i; return ist; + } + + inline ostream & operator<< (ostream & ost, const ElementIndex & pi) + { + return (ost << int(pi)); + } + + + class SurfaceElementIndex + { + int i; + public: + SurfaceElementIndex () { ; } + SurfaceElementIndex (int ai) : i(ai) { ; } + SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) + { i = ai.i; return *this; } + SurfaceElementIndex & operator= (int ai) { i = ai; return *this; } + operator int () const { return i; } + SurfaceElementIndex & operator++ (int) { i++; return *this; } + SurfaceElementIndex & operator-- (int) { i--; return *this; } + }; + + inline istream & operator>> (istream & ist, SurfaceElementIndex & pi) + { + int i; ist >> i; pi = i; return ist; + } + + inline ostream & operator<< (ostream & ost, const SurfaceElementIndex & pi) + { + return (ost << int(pi)); + } + + class SegmentIndex + { + int i; + public: + SegmentIndex () { ; } + SegmentIndex (int ai) : i(ai) { ; } + SegmentIndex & operator= (const SegmentIndex & ai) + { i = ai.i; return *this; } + SegmentIndex & operator= (int ai) { i = ai; return *this; } + operator int () const { return i; } + SegmentIndex & operator++ (int) { i++; return *this; } + SegmentIndex & operator-- (int) { i--; return *this; } + }; + + inline istream & operator>> (istream & ist, SegmentIndex & pi) + { + int i; ist >> i; pi = i; return ist; + } + + inline ostream & operator<< (ostream & ost, const SegmentIndex & pi) + { + return (ost << int(pi)); + } + + + + + /** + Point in the mesh. + Contains layer (a new feature in 4.3 for overlapping meshes. + */ + class MeshPoint : public Point<3> + { + int layer; + double singular; // singular factor for hp-refinement + POINTTYPE type; + +#ifdef PARALLEL + bool isghost; +#endif + + public: + MeshPoint () + { +#ifdef PARALLEL + isghost = 0; +#endif + ; + } + + MeshPoint (const Point<3> & ap, int alayer = 1, POINTTYPE apt = INNERPOINT) + : Point<3> (ap), layer(alayer), singular(0.),type(apt) + { +#ifdef PARALLEL + isghost = 0; +#endif + ; + } - /// - void SetIndex (int si) { index = si; } - /// - int GetIndex () const { return index; } + void SetPoint (const Point<3> & ap) + { + Point<3>::operator= (ap); + layer = 0; + singular = 0; +#ifdef PARALLEL + isghost = 0; +#endif + } - int GetOrder () const { return orderx; } - void SetOrder (const int aorder) ; + int GetLayer() const { return layer; } - void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz = orderz; } - void SetOrder (const int ox, const int oy, const int oz); - // void GetLevel (int & ox, int & oy, int & oz) const { ox = levelx; oy = levely; oz = levelz; } - // void SetLevel (int ox, int oy, int oz) { levelx = ox; levely = oy; levelz = oz; } + POINTTYPE Type() const { return type; } + void SetType(POINTTYPE at) { type = at; } + + double Singularity() const { return singular; } + void Singularity(double s) { singular = s; } + bool IsSingular() const { return (singular != 0.0); } + +#ifdef PARALLEL + bool IsGhost () const { return isghost; } + void SetGhost ( bool aisghost ) { isghost = aisghost; } +#endif + + }; + + inline ostream & operator<<(ostream & s, const MeshPoint & pt) + { + return (s << Point<3> (pt)); + } - /// - void GetBox (const T_POINTS & points, Box3d & box) const; - /// Calculates Volume of elemenet - double Volume (const T_POINTS & points) const; - /// - virtual void Print (ostream & ost) const; - /// - int GetNFaces () const + + + typedef Array T_POINTS; + + + + /** + Triangle element for surface mesh generation. + */ + class Element2d + { + /// point numbers + PointIndex pnum[ELEMENT2D_MAXPOINTS]; + /// geom info of points + PointGeomInfo geominfo[ELEMENT2D_MAXPOINTS]; + + /// surface nr + int index:16; + /// + ELEMENT_TYPE typ:6; + /// number of points + unsigned int np:4; + bool badel:1; + bool refflag:1; // marked for refinement + bool strongrefflag:1; + bool deleted:1; // element is deleted + + // Philippose - 08 August 2010 + // Set a new property for each element, to + // control whether it is visible or not + bool visible:1; // element visible + + /// order for hp-FEM + unsigned int orderx:6; + unsigned int ordery:6; + +#ifdef PARALLEL + bool isghost; + int partitionNumber; +#endif + + /// a linked list for all segments in the same face + SurfaceElementIndex next; + + public: + /// + Element2d (); + /// + Element2d (int anp); + /// + Element2d (ELEMENT_TYPE type); + /// + Element2d (int pi1, int pi2, int pi3); + /// + Element2d (int pi1, int pi2, int pi3, int pi4); + /// + ELEMENT_TYPE GetType () const { return typ; } + /// + void SetType (ELEMENT_TYPE atyp) + { + typ = atyp; + switch (typ) + { + case TRIG: np = 3; break; + case QUAD: np = 4; break; + case TRIG6: np = 6; break; + case QUAD6: np = 6; break; + case QUAD8: np = 8; break; + default: + PrintSysError ("Element2d::SetType, illegal type ", typ); + } + } + /// + int GetNP() const { return np; } + /// + int GetNV() const + { + switch (typ) + { + case TRIG: + case TRIG6: return 3; + + case QUAD: + case QUAD8: + case QUAD6: return 4; + default: +#ifdef DEBUG + PrintSysError ("element2d::GetNV not implemented for typ", typ) +#endif + ; + } + return np; + } + + /// + PointIndex & operator[] (int i) { return pnum[i]; } + /// + const PointIndex & operator[] (int i) const { return pnum[i]; } + + /// + PointIndex & PNum (int i) { return pnum[i-1]; } + /// + const PointIndex & PNum (int i) const { return pnum[i-1]; } + /// + PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } + /// + const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } + /// + + /// + PointGeomInfo & GeomInfoPi (int i) { return geominfo[i-1]; } + /// + const PointGeomInfo & GeomInfoPi (int i) const { return geominfo[i-1]; } + /// + PointGeomInfo & GeomInfoPiMod (int i) { return geominfo[(i-1) % np]; } + /// + const PointGeomInfo & GeomInfoPiMod (int i) const { return geominfo[(i-1) % np]; } + + + void SetIndex (int si) { index = si; } + /// + int GetIndex () const { return index; } + + int GetOrder () const { return orderx; } + void SetOrder (int aorder) { orderx = ordery = aorder; } + + + void GetOrder (int & ox, int & oy) const { ox = orderx, oy =ordery;}; + void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz=0; } + void SetOrder (int ox, int oy, int /* oz */) { orderx = ox; ordery = oy;} + void SetOrder (int ox, int oy) { orderx = ox; ordery = oy;} + + + /// + void GetBox (const T_POINTS & points, Box3d & box) const; + /// invert orientation + inline void Invert (); + /// + void Invert2 (); + /// first point number is smallest + inline void NormalizeNumbering (); + /// + void NormalizeNumbering2 (); + + bool BadElement() const { return badel; } + + // friend ostream & operator<<(ostream & s, const Element2d & el); + friend class Mesh; + + + /// get number of 'integration points' + int GetNIP () const; + void GetIntegrationPoint (int ip, Point2d & p, double & weight) const; + + void GetTransformation (int ip, const Array & points, + class DenseMatrix & trans) const; + void GetTransformation (int ip, class DenseMatrix & pmat, + class DenseMatrix & trans) const; + + void GetShape (const Point2d & p, class Vector & shape) const; + void GetShapeNew (const Point<2> & p, class FlatVector & shape) const; + /// matrix 2 * np + void GetDShape (const Point2d & p, class DenseMatrix & dshape) const; + void GetDShapeNew (const Point<2> & p, class MatrixFixWidth<2> & dshape) const; + /// matrix 2 * np + void GetPointMatrix (const Array & points, + class DenseMatrix & pmat) const; + + void ComputeIntegrationPointData () const; + + + double CalcJacobianBadness (const Array & points) const; + double CalcJacobianBadness (const T_POINTS & points, + const Vec<3> & n) const; + double CalcJacobianBadnessDirDeriv (const Array & points, + int pi, Vec2d & dir, double & dd) const; + + + + void Delete () { deleted = 1; pnum[0] = pnum[1] = pnum[2] = pnum[3] = PointIndex::BASE-1; } + bool IsDeleted () const + { +#ifdef DEBUG + if (pnum[0] < PointIndex::BASE && !deleted) + cerr << "Surfelement has illegal pnum, but not marked as deleted" << endl; +#endif + return deleted; + } + + // Philippose - 08 August 2010 + // Access functions for the new property: visible + void Visible(bool vis = 1) + { visible = vis; } + bool IsVisible () const + { return visible; } + + void SetRefinementFlag (bool rflag = 1) + { refflag = rflag; } + bool TestRefinementFlag () const + { return refflag; } + + void SetStrongRefinementFlag (bool rflag = 1) + { strongrefflag = rflag; } + bool TestStrongRefinementFlag () const + { return strongrefflag; } + + + SurfaceElementIndex NextElement() { return next; } + + bool operator==(const Element2d & el2) const; + + int HasFace(const Element2d& el) const; + /// + int meshdocval; + /// + int hp_elnr; + +#ifdef PARALLEL + bool IsGhost () const { return isghost; } + void SetGhost ( bool aisghost ) { isghost = aisghost; } + + // by JS, only for 2D meshes ???? + int GetPartition () const { return partitionNumber; } + void SetPartition (int nr) { partitionNumber = nr; }; + +#else + bool IsGhost () const { return false; } +#endif + }; + + + ostream & operator<<(ostream & s, const Element2d & el); + + + + + + class IntegrationPointData + { + public: + Point<3> p; + double weight; + Vector shape; + DenseMatrix dshape; + }; + + + + + + + + + /** + Volume element + */ + class Element + { + private: + /// point numbers + PointIndex pnum[ELEMENT_MAXPOINTS]; + /// + ELEMENT_TYPE typ:6; + /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) + int np:5; + /// + class flagstruct { + public: + bool marked:1; // marked for refinement + bool badel:1; // angles worse then limit + bool reverse:1; // for refinement a la Bey + bool illegal:1; // illegal, will be split or swaped + bool illegal_valid:1; // is illegal-flag valid ? + bool badness_valid:1; // is badness valid ? + bool refflag:1; // mark element for refinement + bool strongrefflag:1; + bool deleted:1; // element is deleted, will be removed from array + bool fixed:1; // don't change element in optimization + }; + + /// sub-domain index + short int index; + /// order for hp-FEM + unsigned int orderx:6; + unsigned int ordery:6; + unsigned int orderz:6; + /* unsigned int levelx:6; + unsigned int levely:6; + unsigned int levelz:6; */ + /// stored shape-badness of element + float badness; + +#ifdef PARALLEL + /// number of partition for parallel computation + int partitionNumber; + bool isghost; +#endif + + public: + flagstruct flags; + + /// + Element (); + /// + Element (int anp); + /// + Element (ELEMENT_TYPE type); + /// + Element & operator= (const Element & el2); + + /// + void SetNP (int anp); + /// + void SetType (ELEMENT_TYPE atyp); + /// + int GetNP () const { return np; } + /// + int GetNV() const + { + switch (typ) + { + case TET: + case TET10: + return 4; + case PRISM12: + case PRISM: + return 6; + case PYRAMID: + return 5; + case HEX: + return 8; + default: +#ifdef DEBUG + PrintSysError ("Element3d::GetNV not implemented for typ ", typ) +#endif + ; + } + return np; + } + + bool operator==(const Element & el2) const; + + // old style: + int NP () const { return np; } + + /// + ELEMENT_TYPE GetType () const { return typ; } + + /// + PointIndex & operator[] (int i) { return pnum[i]; } + /// + const PointIndex & operator[] (int i) const { return pnum[i]; } + + /// + PointIndex & PNum (int i) { return pnum[i-1]; } + /// + const PointIndex & PNum (int i) const { return pnum[i-1]; } + /// + PointIndex & PNumMod (int i) { return pnum[(i-1) % np]; } + /// + const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; } + + /// + void SetIndex (int si) { index = si; } + /// + int GetIndex () const { return index; } + + int GetOrder () const { return orderx; } + void SetOrder (const int aorder) ; + + void GetOrder (int & ox, int & oy, int & oz) const { ox = orderx; oy = ordery; oz = orderz; } + void SetOrder (const int ox, const int oy, const int oz); + // void GetLevel (int & ox, int & oy, int & oz) const { ox = levelx; oy = levely; oz = levelz; } + // void SetLevel (int ox, int oy, int oz) { levelx = ox; levely = oy; levelz = oz; } + + + /// + void GetBox (const T_POINTS & points, Box3d & box) const; + /// Calculates Volume of elemenet + double Volume (const T_POINTS & points) const; + /// + virtual void Print (ostream & ost) const; + /// + int GetNFaces () const { switch (typ) { @@ -681,619 +685,619 @@ public: } return 0; } - /// - inline void GetFace (int i, Element2d & face) const; - /// - void GetFace2 (int i, Element2d & face) const; - /// - void Invert (); + /// + inline void GetFace (int i, Element2d & face) const; + /// + void GetFace2 (int i, Element2d & face) const; + /// + void Invert (); - /// split into 4 node tets - void GetTets (Array & locels) const; - /// split into 4 node tets, local point nrs - void GetTetsLocal (Array & locels) const; - /// returns coordinates of nodes - // void GetNodesLocal (Array > & points) const; - void GetNodesLocalNew (Array > & points) const; + /// split into 4 node tets + void GetTets (Array & locels) const; + /// split into 4 node tets, local point nrs + void GetTetsLocal (Array & locels) const; + /// returns coordinates of nodes + // void GetNodesLocal (Array > & points) const; + void GetNodesLocalNew (Array > & points) const; - /// split surface into 3 node trigs - void GetSurfaceTriangles (Array & surftrigs) const; + /// split surface into 3 node trigs + void GetSurfaceTriangles (Array & surftrigs) const; - /// get number of 'integration points' - int GetNIP () const; - void GetIntegrationPoint (int ip, Point<3> & p, double & weight) const; + /// get number of 'integration points' + int GetNIP () const; + void GetIntegrationPoint (int ip, Point<3> & p, double & weight) const; - void GetTransformation (int ip, const T_POINTS & points, - class DenseMatrix & trans) const; - void GetTransformation (int ip, class DenseMatrix & pmat, - class DenseMatrix & trans) const; + void GetTransformation (int ip, const T_POINTS & points, + class DenseMatrix & trans) const; + void GetTransformation (int ip, class DenseMatrix & pmat, + class DenseMatrix & trans) const; - void GetShape (const Point<3> & p, class Vector & shape) const; - void GetShapeNew (const Point<3> & p, class FlatVector & shape) const; - /// matrix 2 * np - void GetDShape (const Point<3> & p, class DenseMatrix & dshape) const; - void GetDShapeNew (const Point<3> & p, class MatrixFixWidth<3> & dshape) const; - /// matrix 3 * np - void GetPointMatrix (const T_POINTS & points, - class DenseMatrix & pmat) const; + void GetShape (const Point<3> & p, class Vector & shape) const; + void GetShapeNew (const Point<3> & p, class FlatVector & shape) const; + /// matrix 2 * np + void GetDShape (const Point<3> & p, class DenseMatrix & dshape) const; + void GetDShapeNew (const Point<3> & p, class MatrixFixWidth<3> & dshape) const; + /// matrix 3 * np + void GetPointMatrix (const T_POINTS & points, + class DenseMatrix & pmat) const; - void ComputeIntegrationPointData () const; + void ComputeIntegrationPointData () const; - double CalcJacobianBadness (const T_POINTS & points) const; - double CalcJacobianBadnessDirDeriv (const T_POINTS & points, - int pi, Vec<3> & dir, double & dd) const; - double CalcJacobianBadnessGradient (const T_POINTS & points, - int pi, Vec<3> & grad) const; + double CalcJacobianBadness (const T_POINTS & points) const; + double CalcJacobianBadnessDirDeriv (const T_POINTS & points, + int pi, Vec<3> & dir, double & dd) const; + double CalcJacobianBadnessGradient (const T_POINTS & points, + int pi, Vec<3> & grad) const; - /// - // friend ostream & operator<<(ostream & s, const Element & el); + /// + // friend ostream & operator<<(ostream & s, const Element & el); - void SetRefinementFlag (bool rflag = 1) - { flags.refflag = rflag; } - int TestRefinementFlag () const - { return flags.refflag; } + void SetRefinementFlag (bool rflag = 1) + { flags.refflag = rflag; } + int TestRefinementFlag () const + { return flags.refflag; } - void SetStrongRefinementFlag (bool rflag = 1) - { flags.strongrefflag = rflag; } - int TestStrongRefinementFlag () const - { return flags.strongrefflag; } + void SetStrongRefinementFlag (bool rflag = 1) + { flags.strongrefflag = rflag; } + int TestStrongRefinementFlag () const + { return flags.strongrefflag; } - int Illegal () const + int Illegal () const { return flags.illegal; } - int IllegalValid () const + int IllegalValid () const { return flags.illegal_valid; } - void SetIllegal (int aillegal) - { - flags.illegal = aillegal ? 1 : 0; - flags.illegal_valid = 1; - } - void SetLegal (int alegal) - { - flags.illegal = alegal ? 0 : 1; - flags.illegal_valid = 1; - } + void SetIllegal (int aillegal) + { + flags.illegal = aillegal ? 1 : 0; + flags.illegal_valid = 1; + } + void SetLegal (int alegal) + { + flags.illegal = alegal ? 0 : 1; + flags.illegal_valid = 1; + } - void Delete () { flags.deleted = 1; } - bool IsDeleted () const - { + void Delete () { flags.deleted = 1; } + bool IsDeleted () const + { #ifdef DEBUG - if (pnum[0] < PointIndex::BASE && !flags.deleted) - cerr << "Volelement has illegal pnum, but not marked as deleted" << endl; + if (pnum[0] < PointIndex::BASE && !flags.deleted) + cerr << "Volelement has illegal pnum, but not marked as deleted" << endl; #endif - return flags.deleted; - } + return flags.deleted; + } #ifdef PARALLEL - int GetPartition () const { return partitionNumber; } - void SetPartition (int nr) { partitionNumber = nr; }; + int GetPartition () const { return partitionNumber; } + void SetPartition (int nr) { partitionNumber = nr; }; #endif - int hp_elnr; + int hp_elnr; #ifdef PARALLEL - bool IsGhost () const { return isghost; } + bool IsGhost () const { return isghost; } - void SetGhost ( const bool aisghost ) { isghost = aisghost; } + void SetGhost ( const bool aisghost ) { isghost = aisghost; } #else - bool IsGhost () const { return false; } + bool IsGhost () const { return false; } #endif - // friend class Mesh; -}; + // friend class Mesh; + }; -ostream & operator<<(ostream & s, const Element & el); + ostream & operator<<(ostream & s, const Element & el); -/** - Edge segment. + /** + Edge segment. */ -class Segment -{ -public: - /// - Segment(); - Segment (const Segment& other); + class Segment + { + public: + /// + Segment(); + Segment (const Segment& other); - ~Segment() - { ; } + ~Segment() + { ; } - // friend ostream & operator<<(ostream & s, const Segment & seg); + // friend ostream & operator<<(ostream & s, const Segment & seg); - PointIndex pnums[3]; // p1, p2, pmid + PointIndex pnums[3]; // p1, p2, pmid - int edgenr; - /// - double singedge_left; - double singedge_right; + int edgenr; + /// + double singedge_left; + double singedge_right; - /// 0.. not first segment of segs, 1..first of class, 2..first of class, inverse - unsigned int seginfo:2; + /// 0.. not first segment of segs, 1..first of class, 2..first of class, inverse + unsigned int seginfo:2; - /// surface decoding index - int si; - /// domain number inner side - int domin; - /// domain number outer side - int domout; - /// top-level object number of surface - int tlosurf; - /// - PointGeomInfo geominfo[2]; + /// surface decoding index + int si; + /// domain number inner side + int domin; + /// domain number outer side + int domout; + /// top-level object number of surface + int tlosurf; + /// + PointGeomInfo geominfo[2]; - /// surfaces describing edge - int surfnr1, surfnr2; - /// - EdgePointGeomInfo epgeominfo[2]; - /// - // int pmid; // for second order - /// - int meshdocval; + /// surfaces describing edge + int surfnr1, surfnr2; + /// + EdgePointGeomInfo epgeominfo[2]; + /// + // int pmid; // for second order + /// + int meshdocval; -private: - string* bcname; + private: + string* bcname; -public: - /* - PointIndex operator[] (int i) const - { return (i == 0) ? p1 : p2; } + public: + /* + PointIndex operator[] (int i) const + { return (i == 0) ? p1 : p2; } - PointIndex & operator[] (int i) - { return (i == 0) ? p1 : p2; } - */ + PointIndex & operator[] (int i) + { return (i == 0) ? p1 : p2; } + */ - Segment& operator=(const Segment & other); + Segment& operator=(const Segment & other); - int hp_elnr; + int hp_elnr; - void SetBCName ( string * abcname ) - { - bcname = abcname; - } + void SetBCName ( string * abcname ) + { + bcname = abcname; + } - string * BCNamePtr () - { return bcname; } + string * BCNamePtr () + { return bcname; } - const string * BCNamePtr () const - { return bcname; } + const string * BCNamePtr () const + { return bcname; } - string GetBCName () const - { - if (! bcname ) - { - return "default"; - } - return *bcname; - } + string GetBCName () const + { + if (! bcname ) + { + return "default"; + } + return *bcname; + } - int GetNP() const - { - return (pnums[2] < 0) ? 2 : 3; - } + int GetNP() const + { + return (pnums[2] < 0) ? 2 : 3; + } - ELEMENT_TYPE GetType() const - { - return (pnums[2] < 0) ? SEGMENT : SEGMENT3; - } + ELEMENT_TYPE GetType() const + { + return (pnums[2] < 0) ? SEGMENT : SEGMENT3; + } - PointIndex & operator[] (int i) { return pnums[i]; } - const PointIndex & operator[] (int i) const { return pnums[i]; } -}; + PointIndex & operator[] (int i) { return pnums[i]; } + const PointIndex & operator[] (int i) const { return pnums[i]; } + }; -ostream & operator<<(ostream & s, const Segment & seg); + ostream & operator<<(ostream & s, const Segment & seg); -// class Surface; -// class FaceDescriptor; - -/// -class FaceDescriptor -{ - /// which surface, 0 if not available - int surfnr; - /// domain nr inside - int domin; - /// domain nr outside - int domout; - /// top level object number of surface - int tlosurf; - /// boundary condition property - int bcprop; - // Philippose - 06/07/2009 - // Add capability to store surface colours along with - // other face data - /// surface colour (Default: R=0.0 ; G=1.0 ; B=0.0) - Vec3d surfcolour; + // class Surface; + // class FaceDescriptor; /// - string * bcname; - /// root of linked list - SurfaceElementIndex firstelement; + class FaceDescriptor + { + /// which surface, 0 if not available + int surfnr; + /// domain nr inside + int domin; + /// domain nr outside + int domout; + /// top level object number of surface + int tlosurf; + /// boundary condition property + int bcprop; + // Philippose - 06/07/2009 + // Add capability to store surface colours along with + // other face data + /// surface colour (Default: R=0.0 ; G=1.0 ; B=0.0) + Vec3d surfcolour; + + /// + string * bcname; + /// root of linked list + SurfaceElementIndex firstelement; - double domin_singular; - double domout_singular; + double domin_singular; + double domout_singular; -public: - FaceDescriptor(); - FaceDescriptor(int surfnri, int domini, int domouti, int tlosurfi); - FaceDescriptor(const Segment & seg); - FaceDescriptor(const FaceDescriptor& other); - ~FaceDescriptor() { ; } + public: + FaceDescriptor(); + FaceDescriptor(int surfnri, int domini, int domouti, int tlosurfi); + FaceDescriptor(const Segment & seg); + FaceDescriptor(const FaceDescriptor& other); + ~FaceDescriptor() { ; } - int SegmentFits (const Segment & seg); + int SegmentFits (const Segment & seg); - int SurfNr () const { return surfnr; } - int DomainIn () const { return domin; } - int DomainOut () const { return domout; } - int TLOSurface () const { return tlosurf; } - int BCProperty () const { return bcprop; } + int SurfNr () const { return surfnr; } + int DomainIn () const { return domin; } + int DomainOut () const { return domout; } + int TLOSurface () const { return tlosurf; } + int BCProperty () const { return bcprop; } - double DomainInSingular() const { return domin_singular; } - double DomainOutSingular() const { return domout_singular; } + double DomainInSingular() const { return domin_singular; } + double DomainOutSingular() const { return domout_singular; } - // Philippose - 06/07/2009 - // Get Surface colour - Vec3d SurfColour () const { return surfcolour; } - string GetBCName () const; - // string * BCNamePtr () { return bcname; } - // const string * BCNamePtr () const { return bcname; } - void SetSurfNr (int sn) { surfnr = sn; } - void SetDomainIn (int di) { domin = di; } - void SetDomainOut (int dom) { domout = dom; } - void SetBCProperty (int bc) { bcprop = bc; } - void SetBCName (string * bcn) { bcname = bcn; } - // Philippose - 06/07/2009 - // Set the surface colour - void SetSurfColour (Vec3d colour) { surfcolour = colour; } + // Philippose - 06/07/2009 + // Get Surface colour + Vec3d SurfColour () const { return surfcolour; } + string GetBCName () const; + // string * BCNamePtr () { return bcname; } + // const string * BCNamePtr () const { return bcname; } + void SetSurfNr (int sn) { surfnr = sn; } + void SetDomainIn (int di) { domin = di; } + void SetDomainOut (int dom) { domout = dom; } + void SetBCProperty (int bc) { bcprop = bc; } + void SetBCName (string * bcn) { bcname = bcn; } + // Philippose - 06/07/2009 + // Set the surface colour + void SetSurfColour (Vec3d colour) { surfcolour = colour; } - void SetDomainInSingular (double v) { domin_singular = v; } - void SetDomainOutSingular (double v) { domout_singular = v; } + void SetDomainInSingular (double v) { domin_singular = v; } + void SetDomainOutSingular (double v) { domout_singular = v; } - SurfaceElementIndex FirstElement() { return firstelement; } - // friend ostream & operator<<(ostream & s, const FaceDescriptor & fd); - friend class Mesh; -}; + SurfaceElementIndex FirstElement() { return firstelement; } + // friend ostream & operator<<(ostream & s, const FaceDescriptor & fd); + friend class Mesh; + }; -ostream & operator<< (ostream & s, const FaceDescriptor & fd); + ostream & operator<< (ostream & s, const FaceDescriptor & fd); -class EdgeDescriptor -{ - int tlosurf; - int surfnr[2]; -public: - EdgeDescriptor () - : tlosurf(-1) - { surfnr[0] = surfnr[1] = -1; } - - int SurfNr (int i) const { return surfnr[i]; } - void SetSurfNr (int i, int nr) { surfnr[i] = nr; } - - int TLOSurface() const { return tlosurf; } - void SetTLOSurface (int nr) { tlosurf = nr; } -}; - - - -class MeshingParameters -{ -public: - /** - 3d optimization strategy: - // m .. move nodes - // M .. move nodes, cheap functional - // s .. swap faces - // c .. combine elements - // d .. divide elements - // p .. plot, no pause - // P .. plot, Pause - // h .. Histogramm, no pause - // H .. Histogramm, pause - */ - const char * optimize3d; - /// number of 3d optimization steps - int optsteps3d; - /** - 2d optimization strategy: - // s .. swap, opt 6 lines/node - // S .. swap, optimal elements - // m .. move nodes - // p .. plot, no pause - // P .. plot, pause - // c .. combine - **/ - const char * optimize2d; - /// number of 2d optimization steps - int optsteps2d; - /// power of error (to approximate max err optimization) - double opterrpow; - /// do block filling ? - int blockfill; - /// block filling up to distance - double filldist; - /// radius of local environment (times h) - double safety; - /// radius of active environment (times h) - double relinnersafety; - /// use local h ? - int uselocalh; - /// grading for local h - double grading; - /// use delaunay meshing - int delaunay; - /// maximal mesh size - double maxh; - /// minimal mesh size - double minh; - /// file for meshsize - const char * meshsizefilename; - /// start surfacemeshing from everywhere in surface - int startinsurface; - /// check overlapping surfaces (debug) - int checkoverlap; - /// check overlapping surface mesh before volume meshing - int checkoverlappingboundary; - /// check chart boundary (sometimes too restrictive) - int checkchartboundary; - /// safty factor for curvatures (elemetns per radius) - double curvaturesafety; - /// minimal number of segments per edge - double segmentsperedge; - /// use parallel threads - int parthread; - /// weight of element size w.r.t element shape - double elsizeweight; - /// init with default values - - - /// from mp3: - /// give up quality class, 2d meshing - int giveuptol2d; - /// give up quality class, 3d meshing - int giveuptol; - /// maximal outer steps - int maxoutersteps; - /// class starting star-shape filling - int starshapeclass; - /// if non-zero, baseelement must have baseelnp points - int baseelnp; - /// quality tolerances are handled less careful - int sloppy; - - /// limit for max element angle (150-180) - double badellimit; - - bool check_impossible; - - /// - int secondorder; - /// high order element curvature - int elementorder; - /// quad-dominated surface meshing - int quad; - /// - int inverttets; - /// - int inverttrigs; - /// - int autozrefine; - /// - MeshingParameters (); - /// - void Print (ostream & ost) const; - - void CopyFrom(const MeshingParameters & other); -}; - - - -class DebugParameters -{ -public: - /// - int debugoutput; - /// use slow checks - int slowchecks; - /// - int haltsuccess; - /// - int haltnosuccess; - /// - int haltlargequalclass; - /// - int haltsegment; - /// - int haltnode; - /// - int haltsegmentp1; - /// - int haltsegmentp2; - /// - int haltexistingline; - /// - int haltoverlap; - /// - int haltface; - /// - int haltfacenr; - /// - DebugParameters (); -}; - - - - -inline void Element2d :: Invert() -{ - if (typ == TRIG) - Swap (PNum(2), PNum(3)); - else - Invert2(); -} - - - - -inline void Element2d :: NormalizeNumbering () -{ - if (GetNP() == 3) - { - if (PNum(1) < PNum(2) && PNum(1) < PNum(3)) - return; - else - { - if (PNum(2) < PNum(3)) - { - PointIndex pi1 = PNum(2); - PNum(2) = PNum(3); - PNum(3) = PNum(1); - PNum(1) = pi1; - } - else - { - PointIndex pi1 = PNum(3); - PNum(3) = PNum(2); - PNum(2) = PNum(1); - PNum(1) = pi1; - } - } - } - else - NormalizeNumbering2(); -} - - - -static const int gftetfacesa[4][3] = - { { 1, 2, 3 }, - { 2, 0, 3 }, - { 0, 1, 3 }, - { 1, 0, 2 } }; - -inline void Element :: GetFace (int i, Element2d & face) const -{ - if (typ == TET) - { - face.SetType(TRIG); - face[0] = pnum[gftetfacesa[i-1][0]]; - face[1] = pnum[gftetfacesa[i-1][1]]; - face[2] = pnum[gftetfacesa[i-1][2]]; - } - else - GetFace2 (i, face); -} - - - - - - - -/** - Identification of periodic surfaces, close surfaces, etc. - */ -class Identifications -{ -public: - enum ID_TYPE { UNDEFINED = 1, PERIODIC = 2, CLOSESURFACES = 3, CLOSEEDGES = 4}; - - -private: - class Mesh & mesh; - - /// identify points (thin layers, periodic b.c.) - INDEX_2_HASHTABLE * identifiedpoints; - - /// the same, with info about the id-nr - INDEX_3_HASHTABLE * identifiedpoints_nr; - - /// sorted by identification nr - TABLE idpoints_table; - - Array type; - - /// number of identifications (or, actually used identifications ?) - int maxidentnr; - -public: - /// - Identifications (class Mesh & amesh); - /// - ~Identifications (); - - void Delete (); - - /* - Identify points pi1 and pi2, due to - identification nr identnr - */ - void Add (PointIndex pi1, PointIndex pi2, int identnr); - - - int Get (PointIndex pi1, PointIndex pi2) const; - int GetSymmetric (PointIndex pi1, PointIndex pi2) const; - - bool Get (PointIndex pi1, PointIndex pi2, int identnr) const; - bool GetSymmetric (PointIndex pi1, PointIndex pi2, int identnr) const; - - /// - INDEX_2_HASHTABLE & GetIdentifiedPoints () - { - return *identifiedpoints; - } - - bool Used (PointIndex pi1, PointIndex pi2) + class EdgeDescriptor { - return identifiedpoints->Used (INDEX_2 (pi1, pi2)); - } + int tlosurf; + int surfnr[2]; + public: + EdgeDescriptor () + : tlosurf(-1) + { surfnr[0] = surfnr[1] = -1; } - bool UsedSymmetric (PointIndex pi1, PointIndex pi2) - { - return - identifiedpoints->Used (INDEX_2 (pi1, pi2)) || - identifiedpoints->Used (INDEX_2 (pi2, pi1)); - } + int SurfNr (int i) const { return surfnr[i]; } + void SetSurfNr (int i, int nr) { surfnr[i] = nr; } - /// - void GetMap (int identnr, Array & identmap, bool symmetric = false) const; - /// - ID_TYPE GetType(int identnr) const + int TLOSurface() const { return tlosurf; } + void SetTLOSurface (int nr) { tlosurf = nr; } + }; + + + + class MeshingParameters { - if(identnr <= type.Size()) - return type[identnr-1]; + public: + /** + 3d optimization strategy: + // m .. move nodes + // M .. move nodes, cheap functional + // s .. swap faces + // c .. combine elements + // d .. divide elements + // p .. plot, no pause + // P .. plot, Pause + // h .. Histogramm, no pause + // H .. Histogramm, pause + */ + const char * optimize3d; + /// number of 3d optimization steps + int optsteps3d; + /** + 2d optimization strategy: + // s .. swap, opt 6 lines/node + // S .. swap, optimal elements + // m .. move nodes + // p .. plot, no pause + // P .. plot, pause + // c .. combine + **/ + const char * optimize2d; + /// number of 2d optimization steps + int optsteps2d; + /// power of error (to approximate max err optimization) + double opterrpow; + /// do block filling ? + int blockfill; + /// block filling up to distance + double filldist; + /// radius of local environment (times h) + double safety; + /// radius of active environment (times h) + double relinnersafety; + /// use local h ? + int uselocalh; + /// grading for local h + double grading; + /// use delaunay meshing + int delaunay; + /// maximal mesh size + double maxh; + /// minimal mesh size + double minh; + /// file for meshsize + const char * meshsizefilename; + /// start surfacemeshing from everywhere in surface + int startinsurface; + /// check overlapping surfaces (debug) + int checkoverlap; + /// check overlapping surface mesh before volume meshing + int checkoverlappingboundary; + /// check chart boundary (sometimes too restrictive) + int checkchartboundary; + /// safty factor for curvatures (elemetns per radius) + double curvaturesafety; + /// minimal number of segments per edge + double segmentsperedge; + /// use parallel threads + int parthread; + /// weight of element size w.r.t element shape + double elsizeweight; + /// init with default values + + + /// from mp3: + /// give up quality class, 2d meshing + int giveuptol2d; + /// give up quality class, 3d meshing + int giveuptol; + /// maximal outer steps + int maxoutersteps; + /// class starting star-shape filling + int starshapeclass; + /// if non-zero, baseelement must have baseelnp points + int baseelnp; + /// quality tolerances are handled less careful + int sloppy; + + /// limit for max element angle (150-180) + double badellimit; + + bool check_impossible; + + /// + int secondorder; + /// high order element curvature + int elementorder; + /// quad-dominated surface meshing + int quad; + /// + int inverttets; + /// + int inverttrigs; + /// + int autozrefine; + /// + MeshingParameters (); + /// + void Print (ostream & ost) const; + + void CopyFrom(const MeshingParameters & other); + }; + + + + class DebugParameters + { + public: + /// + int debugoutput; + /// use slow checks + int slowchecks; + /// + int haltsuccess; + /// + int haltnosuccess; + /// + int haltlargequalclass; + /// + int haltsegment; + /// + int haltnode; + /// + int haltsegmentp1; + /// + int haltsegmentp2; + /// + int haltexistingline; + /// + int haltoverlap; + /// + int haltface; + /// + int haltfacenr; + /// + DebugParameters (); + }; + + + + + inline void Element2d :: Invert() + { + if (typ == TRIG) + Swap (PNum(2), PNum(3)); else - return UNDEFINED; + Invert2(); } - void SetType(int identnr, ID_TYPE t) + + + + + inline void Element2d :: NormalizeNumbering () { - while(type.Size() < identnr) - type.Append(UNDEFINED); - type[identnr-1] = t; + if (GetNP() == 3) + { + if (PNum(1) < PNum(2) && PNum(1) < PNum(3)) + return; + else + { + if (PNum(2) < PNum(3)) + { + PointIndex pi1 = PNum(2); + PNum(2) = PNum(3); + PNum(3) = PNum(1); + PNum(1) = pi1; + } + else + { + PointIndex pi1 = PNum(3); + PNum(3) = PNum(2); + PNum(2) = PNum(1); + PNum(1) = pi1; + } + } + } + else + NormalizeNumbering2(); } + + + + static const int gftetfacesa[4][3] = + { { 1, 2, 3 }, + { 2, 0, 3 }, + { 0, 1, 3 }, + { 1, 0, 2 } }; + + inline void Element :: GetFace (int i, Element2d & face) const + { + if (typ == TET) + { + face.SetType(TRIG); + face[0] = pnum[gftetfacesa[i-1][0]]; + face[1] = pnum[gftetfacesa[i-1][1]]; + face[2] = pnum[gftetfacesa[i-1][2]]; + } + else + GetFace2 (i, face); + } + + + + + + + + /** + Identification of periodic surfaces, close surfaces, etc. + */ + class Identifications + { + public: + enum ID_TYPE { UNDEFINED = 1, PERIODIC = 2, CLOSESURFACES = 3, CLOSEEDGES = 4}; + + + private: + class Mesh & mesh; + + /// identify points (thin layers, periodic b.c.) + INDEX_2_HASHTABLE * identifiedpoints; + + /// the same, with info about the id-nr + INDEX_3_HASHTABLE * identifiedpoints_nr; + + /// sorted by identification nr + TABLE idpoints_table; + + Array type; + + /// number of identifications (or, actually used identifications ?) + int maxidentnr; + + public: + /// + Identifications (class Mesh & amesh); + /// + ~Identifications (); + + void Delete (); + + /* + Identify points pi1 and pi2, due to + identification nr identnr + */ + void Add (PointIndex pi1, PointIndex pi2, int identnr); + + + int Get (PointIndex pi1, PointIndex pi2) const; + int GetSymmetric (PointIndex pi1, PointIndex pi2) const; + + bool Get (PointIndex pi1, PointIndex pi2, int identnr) const; + bool GetSymmetric (PointIndex pi1, PointIndex pi2, int identnr) const; + + /// + INDEX_2_HASHTABLE & GetIdentifiedPoints () + { + return *identifiedpoints; + } + + bool Used (PointIndex pi1, PointIndex pi2) + { + return identifiedpoints->Used (INDEX_2 (pi1, pi2)); + } + + bool UsedSymmetric (PointIndex pi1, PointIndex pi2) + { + return + identifiedpoints->Used (INDEX_2 (pi1, pi2)) || + identifiedpoints->Used (INDEX_2 (pi2, pi1)); + } + + /// + void GetMap (int identnr, Array & identmap, bool symmetric = false) const; + /// + ID_TYPE GetType(int identnr) const + { + if(identnr <= type.Size()) + return type[identnr-1]; + else + return UNDEFINED; + } + void SetType(int identnr, ID_TYPE t) + { + while(type.Size() < identnr) + type.Append(UNDEFINED); + type[identnr-1] = t; + } - /// - void GetPairs (int identnr, Array & identpairs) const; - /// - int GetMaxNr () const { return maxidentnr; } + /// + void GetPairs (int identnr, Array & identpairs) const; + /// + int GetMaxNr () const { return maxidentnr; } - /// remove secondorder - void SetMaxPointNr (int maxpnum); - - void Print (ostream & ost) const; -}; + /// remove secondorder + void SetMaxPointNr (int maxpnum); + void Print (ostream & ost) const; + }; +} diff --git a/libsrc/meshing/msghandler.hpp b/libsrc/meshing/msghandler.hpp index 843cd186..e7632873 100644 --- a/libsrc/meshing/msghandler.hpp +++ b/libsrc/meshing/msghandler.hpp @@ -8,45 +8,49 @@ /**************************************************************************/ -extern void PrintDot(char ch = '.'); +namespace netgen +{ + + extern void PrintDot(char ch = '.'); -//Message Pipeline: + //Message Pipeline: -//importance: importance of message: 1=very important, 3=middle, 5=low, 7=unimportant -extern void PrintMessage(int importance, - const MyStr& s1, const MyStr& s2=MyStr()); -extern void PrintMessage(int importance, - const MyStr& s1, const MyStr& s2, const MyStr& s3, const MyStr& s4=MyStr()); -extern void PrintMessage(int importance, - const MyStr& s1, const MyStr& s2, const MyStr& s3, const MyStr& s4, - const MyStr& s5, const MyStr& s6=MyStr(), const MyStr& s7=MyStr(), const MyStr& s8=MyStr()); + //importance: importance of message: 1=very important, 3=middle, 5=low, 7=unimportant + extern void PrintMessage(int importance, + const MyStr& s1, const MyStr& s2=MyStr()); + extern void PrintMessage(int importance, + const MyStr& s1, const MyStr& s2, const MyStr& s3, const MyStr& s4=MyStr()); + extern void PrintMessage(int importance, + const MyStr& s1, const MyStr& s2, const MyStr& s3, const MyStr& s4, + const MyStr& s5, const MyStr& s6=MyStr(), const MyStr& s7=MyStr(), const MyStr& s8=MyStr()); -// CR without line-feed -extern void PrintMessageCR(int importance, - const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + // CR without line-feed + extern void PrintMessageCR(int importance, + const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); + extern void PrintFnStart(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void PrintFnStart(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + extern void PrintWarning(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); + extern void PrintError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void PrintWarning(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", - const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void PrintError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", - const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void PrintFileError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", - const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void PrintSysError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", - const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void PrintUserError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", - const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void PrintTime(const MyStr& s1="", const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", - const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); -extern void SetStatMsg(const MyStr& s); + extern void PrintFileError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); + extern void PrintSysError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); + extern void PrintUserError(const MyStr& s1, const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); + extern void PrintTime(const MyStr& s1="", const MyStr& s2="", const MyStr& s3="", const MyStr& s4="", + const MyStr& s5="", const MyStr& s6="", const MyStr& s7="", const MyStr& s8=""); + extern void SetStatMsg(const MyStr& s); -extern void PushStatus(const MyStr& s); -extern void PushStatusF(const MyStr& s); -extern void PopStatus(); -extern void SetThreadPercent(double percent); -extern void GetStatus(MyStr & s, double & percentage); + extern void PushStatus(const MyStr& s); + extern void PushStatusF(const MyStr& s); + extern void PopStatus(); + extern void SetThreadPercent(double percent); + extern void GetStatus(MyStr & s, double & percentage); +} #endif diff --git a/libsrc/meshing/paralleltop.hpp b/libsrc/meshing/paralleltop.hpp index 971394ca..f8429d98 100644 --- a/libsrc/meshing/paralleltop.hpp +++ b/libsrc/meshing/paralleltop.hpp @@ -1,267 +1,268 @@ #ifndef FILE_PARALLELTOP #define FILE_PARALLELTOP -#include - -extern int ntasks; - -class ParallelMeshTopology +namespace netgen { - const Mesh & mesh; - // number of local elements, vertices, points (?), edges, faces - int ne, nv, np, ned, nfa; + extern int ntasks; - // number of local segments and surface elements - int nseg, nsurfel; - - // number of global elements, vertices, ???, faces - int neglob, nvglob; - int nparel; - int nfaglob; - - /** - mapping from local to distant vertex number - each row of the table corresponds to one vertex - each row contains a list of pairs (procnr, dist_vnum) - */ - TABLE loc2distvert; - TABLE loc2distedge, loc2distface, loc2distel; - TABLE loc2distsegm, loc2distsurfel; - - BitArray * isexchangeface, * isexchangevert, * isexchangeedge, * isexchangeel; - - BitArray isghostedge, isghostface; - - bool coarseupdate; - int overlap; - -public: - - ParallelMeshTopology (const Mesh & amesh); - ~ParallelMeshTopology (); - - /// set number of local vertices, reset sizes of loc2dist_vert, isexchangevert... - void SetNV ( int anv ); - void SetNE ( int ane ); - void SetNSE ( int anse ); - void SetNSegm ( int anseg ); - - void Reset (); - - void SetLoc2Glob_Vert ( int locnum, int globnum ) { loc2distvert[locnum][0] = globnum; } - void SetLoc2Glob_VolEl ( int locnum, int globnum ) { loc2distel[locnum-1][0] = globnum; } - void SetLoc2Glob_SurfEl ( int locnum, int globnum ) { loc2distsurfel[locnum-1][0] = globnum; } - void SetLoc2Glob_Segm ( int locnum, int globnum ) { loc2distsegm[locnum-1][0] = globnum; } - - int GetLoc2Glob_Vert ( int locnum ) const { return loc2distvert[locnum][0]; } - int GetLoc2Glob_VolEl ( int locnum ) const { return loc2distel[locnum-1][0]; } - - void GetVertNeighbours ( int vnum, Array & dests ) const; - - int Glob2Loc_SurfEl ( int globnum ); - int Glob2Loc_VolEl ( int globnum ); - int Glob2Loc_Segm ( int globnum ); - int Glob2Loc_Vert ( int globnum ); - - int GetNDistantPNums ( int locpnum ) const { return loc2distvert[locpnum].Size() / 2 + 1; } - int GetNDistantFaceNums ( int locfacenum ) const { return loc2distface[locfacenum-1].Size() / 2 + 1; } - int GetNDistantEdgeNums ( int locedgenum ) const { return loc2distedge[locedgenum-1].Size() / 2 + 1; } - int GetNDistantElNums ( int locelnum ) const { return loc2distel[locelnum-1].Size() / 2 + 1; } - - int GetDistantPNum ( int proc, int locpnum ) const; - int GetDistantEdgeNum ( int proc, int locedgenum ) const; - int GetDistantFaceNum ( int proc, int locedgenum ) const; - int GetDistantElNum ( int proc, int locelnum ) const; - - int GetDistantPNums ( int locpnum, int * distpnums ) const; - int GetDistantEdgeNums ( int locedgenum, int * distedgenums ) const; - int GetDistantFaceNums ( int locedgenum, int * distfacenums ) const; - int GetDistantElNums ( int locelnum, int * distfacenums ) const; - - void Print() const; - - void SetExchangeFace ( int fnr ) { isexchangeface->Set( (fnr-1)*(ntasks+1) ); } - void SetExchangeVert ( int vnum ) { isexchangevert->Set( (vnum-1)*(ntasks+1) ); } - void SetExchangeElement ( int elnum ) { isexchangeel->Set((elnum-1)*(ntasks+1) ); } - void SetExchangeEdge ( int ednum ) { isexchangeedge->Set ((ednum-1)*(ntasks+1)); } - - - // fuer diese Fkts kann man ja O(N) - bitarrays lassen - bool IsExchangeVert ( PointIndex vnum ) const + class ParallelMeshTopology { - return loc2distvert[vnum].Size() > 1; - // return isexchangevert->Test((vnum-1)*(ntasks+1)); - } + const Mesh & mesh; - bool IsExchangeEdge ( int ednum ) const - { - /* - if ( isexchangeedge->Test( (ednum-1)*(ntasks+1) ) != - ( loc2distedge[ednum-1].Size() > 1) ) - { + // number of local elements, vertices, points (?), edges, faces + int ne, nv, np, ned, nfa; + + // number of local segments and surface elements + int nseg, nsurfel; + + // number of global elements, vertices, ???, faces + int neglob, nvglob; + int nparel; + int nfaglob; + + /** + mapping from local to distant vertex number + each row of the table corresponds to one vertex + each row contains a list of pairs (procnr, dist_vnum) + */ + TABLE loc2distvert; + TABLE loc2distedge, loc2distface, loc2distel; + TABLE loc2distsegm, loc2distsurfel; + + BitArray * isexchangeface, * isexchangevert, * isexchangeedge, * isexchangeel; + + BitArray isghostedge, isghostface; + + bool coarseupdate; + int overlap; + + public: + + ParallelMeshTopology (const Mesh & amesh); + ~ParallelMeshTopology (); + + /// set number of local vertices, reset sizes of loc2dist_vert, isexchangevert... + void SetNV ( int anv ); + void SetNE ( int ane ); + void SetNSE ( int anse ); + void SetNSegm ( int anseg ); + + void Reset (); + + void SetLoc2Glob_Vert ( int locnum, int globnum ) { loc2distvert[locnum][0] = globnum; } + void SetLoc2Glob_VolEl ( int locnum, int globnum ) { loc2distel[locnum-1][0] = globnum; } + void SetLoc2Glob_SurfEl ( int locnum, int globnum ) { loc2distsurfel[locnum-1][0] = globnum; } + void SetLoc2Glob_Segm ( int locnum, int globnum ) { loc2distsegm[locnum-1][0] = globnum; } + + int GetLoc2Glob_Vert ( int locnum ) const { return loc2distvert[locnum][0]; } + int GetLoc2Glob_VolEl ( int locnum ) const { return loc2distel[locnum-1][0]; } + + void GetVertNeighbours ( int vnum, Array & dests ) const; + + int Glob2Loc_SurfEl ( int globnum ); + int Glob2Loc_VolEl ( int globnum ); + int Glob2Loc_Segm ( int globnum ); + int Glob2Loc_Vert ( int globnum ); + + int GetNDistantPNums ( int locpnum ) const { return loc2distvert[locpnum].Size() / 2 + 1; } + int GetNDistantFaceNums ( int locfacenum ) const { return loc2distface[locfacenum-1].Size() / 2 + 1; } + int GetNDistantEdgeNums ( int locedgenum ) const { return loc2distedge[locedgenum-1].Size() / 2 + 1; } + int GetNDistantElNums ( int locelnum ) const { return loc2distel[locelnum-1].Size() / 2 + 1; } + + int GetDistantPNum ( int proc, int locpnum ) const; + int GetDistantEdgeNum ( int proc, int locedgenum ) const; + int GetDistantFaceNum ( int proc, int locedgenum ) const; + int GetDistantElNum ( int proc, int locelnum ) const; + + int GetDistantPNums ( int locpnum, int * distpnums ) const; + int GetDistantEdgeNums ( int locedgenum, int * distedgenums ) const; + int GetDistantFaceNums ( int locedgenum, int * distfacenums ) const; + int GetDistantElNums ( int locelnum, int * distfacenums ) const; + + void Print() const; + + void SetExchangeFace ( int fnr ) { isexchangeface->Set( (fnr-1)*(ntasks+1) ); } + void SetExchangeVert ( int vnum ) { isexchangevert->Set( (vnum-1)*(ntasks+1) ); } + void SetExchangeElement ( int elnum ) { isexchangeel->Set((elnum-1)*(ntasks+1) ); } + void SetExchangeEdge ( int ednum ) { isexchangeedge->Set ((ednum-1)*(ntasks+1)); } + + + // fuer diese Fkts kann man ja O(N) - bitarrays lassen + bool IsExchangeVert ( PointIndex vnum ) const + { + return loc2distvert[vnum].Size() > 1; + // return isexchangevert->Test((vnum-1)*(ntasks+1)); + } + + bool IsExchangeEdge ( int ednum ) const + { + /* + if ( isexchangeedge->Test( (ednum-1)*(ntasks+1) ) != + ( loc2distedge[ednum-1].Size() > 1) ) + { cerr << "isexedge differs, id = " << id << ", ednum = " << ednum << endl; - } - */ - // return loc2distedge[ednum-1].Size() > 1; - int i = (ednum-1)*(ntasks+1); - return isexchangeedge->Test(i); - } + } + */ + // return loc2distedge[ednum-1].Size() > 1; + int i = (ednum-1)*(ntasks+1); + return isexchangeedge->Test(i); + } - bool IsExchangeFace ( int fnum ) const - { - /* - if ( isexchangeface->Test( (fnum-1)*(ntasks+1) ) != - ( loc2distface[fnum-1].Size() > 1) ) - { + bool IsExchangeFace ( int fnum ) const + { + /* + if ( isexchangeface->Test( (fnum-1)*(ntasks+1) ) != + ( loc2distface[fnum-1].Size() > 1) ) + { cerr << "it differs, id = " << id << ", fnum = " << fnum << endl; - } - */ - // nur hier funktioniert's so nicht ???? (JS) - // return loc2distface[fnum-1].Size() > 1; - return isexchangeface->Test( (fnum-1)*(ntasks+1) ); - } + } + */ + // nur hier funktioniert's so nicht ???? (JS) + // return loc2distface[fnum-1].Size() > 1; + return isexchangeface->Test( (fnum-1)*(ntasks+1) ); + } - bool IsExchangeElement ( int elnum ) const - { - // return loc2distel[elnum-1].Size() > 1; - return isexchangeel->Test((elnum-1)*(ntasks+1)); - } + bool IsExchangeElement ( int elnum ) const + { + // return loc2distel[elnum-1].Size() > 1; + return isexchangeel->Test((elnum-1)*(ntasks+1)); + } - bool IsExchangeSEl ( int selnum ) const - { - return loc2distsurfel[selnum-1].Size() > 1; - } + bool IsExchangeSEl ( int selnum ) const + { + return loc2distsurfel[selnum-1].Size() > 1; + } - void SetExchangeFace (int dest, int fnr ) - { - // isexchangeface->Set((fnr-1)*(ntasks+1) + (dest+1)); - } + void SetExchangeFace (int dest, int fnr ) + { + // isexchangeface->Set((fnr-1)*(ntasks+1) + (dest+1)); + } - void SetExchangeVert (int dest, int vnum ) - { - // isexchangevert->Set((vnum-1)*(ntasks+1) + (dest+1) ); - } + void SetExchangeVert (int dest, int vnum ) + { + // isexchangevert->Set((vnum-1)*(ntasks+1) + (dest+1) ); + } - void SetExchangeElement (int dest, int vnum ) - { - isexchangeel->Set( (vnum-1)*(ntasks+1) + (dest+1) ); - } + void SetExchangeElement (int dest, int vnum ) + { + isexchangeel->Set( (vnum-1)*(ntasks+1) + (dest+1) ); + } - void SetExchangeEdge (int dest, int ednum ) - { - // isexchangeedge->Set ( (ednum-1)*(ntasks+1) + (dest+1) ); - } + void SetExchangeEdge (int dest, int ednum ) + { + // isexchangeedge->Set ( (ednum-1)*(ntasks+1) + (dest+1) ); + } - bool IsExchangeVert (int dest, int vnum ) const - { - FlatArray exchange = loc2distvert[vnum]; - for (int i = 1; i < exchange.Size(); i += 2) - if (exchange[i] == dest) return true; - return false; - // return isexchangevert->Test((vnum-1)*(ntasks+1) + (dest+1) ); - } + bool IsExchangeVert (int dest, int vnum ) const + { + FlatArray exchange = loc2distvert[vnum]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; + // return isexchangevert->Test((vnum-1)*(ntasks+1) + (dest+1) ); + } - bool IsExchangeEdge (int dest, int ednum ) const - { - FlatArray exchange = loc2distedge[ednum-1]; - for (int i = 1; i < exchange.Size(); i += 2) - if (exchange[i] == dest) return true; - return false; + bool IsExchangeEdge (int dest, int ednum ) const + { + FlatArray exchange = loc2distedge[ednum-1]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; - // return isexchangeedge->Test((ednum-1)*(ntasks+1) + (dest+1)); - } + // return isexchangeedge->Test((ednum-1)*(ntasks+1) + (dest+1)); + } - bool IsExchangeFace (int dest, int fnum ) const - { - FlatArray exchange = loc2distface[fnum-1]; - for (int i = 1; i < exchange.Size(); i += 2) - if (exchange[i] == dest) return true; - return false; + bool IsExchangeFace (int dest, int fnum ) const + { + FlatArray exchange = loc2distface[fnum-1]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; - // return isexchangeface->Test((fnum-1)*(ntasks+1) + (dest+1) ); - } + // return isexchangeface->Test((fnum-1)*(ntasks+1) + (dest+1) ); + } - bool IsExchangeElement (int dest, int elnum ) const - { - // das geht auch nicht - /* - FlatArray exchange = loc2distel[elnum-1]; - for (int i = 1; i < exchange.Size(); i += 2) - if (exchange[i] == dest) return true; - return false; - */ - return isexchangeel->Test((elnum-1)*(ntasks+1) + (dest+1)); - } + bool IsExchangeElement (int dest, int elnum ) const + { + // das geht auch nicht + /* + FlatArray exchange = loc2distel[elnum-1]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; + */ + return isexchangeel->Test((elnum-1)*(ntasks+1) + (dest+1)); + } - int Overlap() const - { return overlap; } + int Overlap() const + { return overlap; } - int IncreaseOverlap () - { overlap++; return overlap; } + int IncreaseOverlap () + { overlap++; return overlap; } - void Update(); + void Update(); - void UpdateCoarseGrid(); - void UpdateCoarseGridOverlap(); - void UpdateRefinement (); - void UpdateTopology (); - void UpdateExchangeElements(); + void UpdateCoarseGrid(); + void UpdateCoarseGridOverlap(); + void UpdateRefinement (); + void UpdateTopology (); + void UpdateExchangeElements(); - void UpdateCoarseGridGlobal(); + void UpdateCoarseGridGlobal(); - bool DoCoarseUpdate() const - { return !coarseupdate; } + bool DoCoarseUpdate() const + { return !coarseupdate; } - void SetDistantFaceNum ( int dest, int locnum, int distnum ); - void SetDistantPNum ( int dest, int locnum, int distnum ); - void SetDistantEdgeNum ( int dest, int locnum, int distnum ); - void SetDistantEl ( int dest, int locnum, int distnum ); - void SetDistantSurfEl ( int dest, int locnum, int distnum ); - void SetDistantSegm ( int dest, int locnum, int distnum ); + void SetDistantFaceNum ( int dest, int locnum, int distnum ); + void SetDistantPNum ( int dest, int locnum, int distnum ); + void SetDistantEdgeNum ( int dest, int locnum, int distnum ); + void SetDistantEl ( int dest, int locnum, int distnum ); + void SetDistantSurfEl ( int dest, int locnum, int distnum ); + void SetDistantSegm ( int dest, int locnum, int distnum ); - bool IsGhostEl ( int elnum ) const { return mesh.VolumeElement(elnum).IsGhost(); } - bool IsGhostFace ( int fanum ) const { return isghostface.Test(fanum-1); } - bool IsGhostEdge ( int ednum ) const { return isghostedge.Test(ednum-1); } - bool IsGhostVert ( int pnum ) const { return mesh.Point(pnum).IsGhost(); } + bool IsGhostEl ( int elnum ) const { return mesh.VolumeElement(elnum).IsGhost(); } + bool IsGhostFace ( int fanum ) const { return isghostface.Test(fanum-1); } + bool IsGhostEdge ( int ednum ) const { return isghostedge.Test(ednum-1); } + bool IsGhostVert ( int pnum ) const { return mesh.Point(pnum).IsGhost(); } -// inline void SetGhostVert ( const int pnum ) -// { isghostvert.Set(pnum-1); } + // inline void SetGhostVert ( const int pnum ) + // { isghostvert.Set(pnum-1); } - void SetGhostEdge ( int ednum ) - { isghostedge.Set(ednum-1); } + void SetGhostEdge ( int ednum ) + { isghostedge.Set(ednum-1); } - void ClearGhostEdge ( int ednum ) - { isghostedge.Clear(ednum-1); } + void ClearGhostEdge ( int ednum ) + { isghostedge.Clear(ednum-1); } - void SetGhostFace ( int fanum ) - { isghostface.Set(fanum-1); } + void SetGhostFace ( int fanum ) + { isghostface.Set(fanum-1); } - void ClearGhostFace ( int fanum ) - { isghostface.Clear(fanum-1); } + void ClearGhostFace ( int fanum ) + { isghostface.Clear(fanum-1); } - bool IsElementInPartition ( int elnum, int dest ) const - { - return IsExchangeElement ( dest, elnum); - } + bool IsElementInPartition ( int elnum, int dest ) const + { + return IsExchangeElement ( dest, elnum); + } - void SetElementInPartition ( int elnum, int dest ) - { - SetExchangeElement ( dest, elnum ); - } + void SetElementInPartition ( int elnum, int dest ) + { + SetExchangeElement ( dest, elnum ); + } - void SetNVGlob ( int anvglob ) { nvglob = anvglob; } - void SetNEGlob ( int aneglob ) { neglob = aneglob; } + void SetNVGlob ( int anvglob ) { nvglob = anvglob; } + void SetNEGlob ( int aneglob ) { neglob = aneglob; } - int GetNVGlob () { return nvglob; } - int GetNEGlob () { return neglob; } -}; + int GetNVGlob () { return nvglob; } + int GetNEGlob () { return neglob; } + }; - +} diff --git a/libsrc/meshing/validate.hpp b/libsrc/meshing/validate.hpp index 71e2f2e4..1cc01bec 100644 --- a/libsrc/meshing/validate.hpp +++ b/libsrc/meshing/validate.hpp @@ -1,17 +1,21 @@ #ifndef VALIDATE_HPP #define VALIDATE_HPP +namespace netgen +{ + + void GetPureBadness(Mesh & mesh, Array & pure_badness, + const BitArray & isnewpoint); + double Validate(const Mesh & mesh, Array & bad_elements, + const Array & pure_badness, + double max_worsening, const bool uselocalworsening, + Array * quality_loss = NULL); + void RepairBisection(Mesh & mesh, Array & bad_elements, + const BitArray & isnewpoint, const Refinement & refinement, + const Array & pure_badness, + double max_worsening, const bool uselocalworsening, + const Array< Array* > & idmaps); - -void GetPureBadness(Mesh & mesh, Array & pure_badness, - const BitArray & isnewpoint); -double Validate(const Mesh & mesh, Array & bad_elements, - const Array & pure_badness, - double max_worsening, const bool uselocalworsening, - Array * quality_loss = NULL); -void RepairBisection(Mesh & mesh, Array & bad_elements, const BitArray & isnewpoint, const Refinement & refinement, - const Array & pure_badness, - double max_worsening, const bool uselocalworsening, - const Array< Array* > & idmaps); +} #endif // VALIDATE_HPP diff --git a/tutorials/Makefile.am b/tutorials/Makefile.am index 17ea0eab..1e8db8e8 100644 --- a/tutorials/Makefile.am +++ b/tutorials/Makefile.am @@ -5,7 +5,7 @@ dist_pkgdata_DATA = boxcyl.geo circle_on_cube.geo cone.geo cube.geo \ cylsphere.geo ellipsoid.geo ellipticcyl.geo extrusion.geo fichera.geo lshape3d.geo \ manyholes.geo manyholes2.geo matrix.geo ortho.geo period.geo revolution.geo \ sculpture.geo shaft.geo shell.geo sphere.geo sphereincube.geo torus.geo trafo.geo \ - twobricks.geo twocubes.geo twocyl.geo \ + twobricks.geo twocubes.geo twocyl.geo boundarycondition.geo \ hinge.stl part1.stl frame.step screw.step \ squarehole.in2d squarecircle.in2d square.in2d diff --git a/tutorials/boundarycondition.geo b/tutorials/boundarycondition.geo new file mode 100644 index 00000000..2e3b7fad --- /dev/null +++ b/tutorials/boundarycondition.geo @@ -0,0 +1,16 @@ +algebraic3d + +solid p1 = plane (0.5, 0, 0; 1, 0, 0); + +# since surfaces of both bricks are identic they get the same bc id: +solid brick1 = orthobrick (0,0,0; 1,1,1) and p1 -bc=1; +solid brick2 = orthobrick (0,0,-1; 1,1,0) and p1 -bc=2; + + +tlo brick1; +tlo brick2; + +# override bc number: +# all faces of solid p1 belonging to the boundary of tlo brick1 get bc=3 + +boundarycondition p1 brick1 3;