NGSolve  4.9
comp/meshaccess.hpp
00001 #ifndef FILE_MESHACCESS
00002 #define FILE_MESHACCESS
00003 
00004 /*********************************************************************/
00005 /* File:   meshaccess.hpp                                            */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   25. Mar. 2000                                             */
00008 /*********************************************************************/
00009 
00010 
00011 #include <nginterface.h>
00012 #include <nginterface_v2.hpp>
00013 
00014 namespace ngfem
00015 {
00016   class ElementTransformation;
00017   class IntegrationPoint;
00018 }
00019 
00020 
00021 
00022 namespace ngcomp
00023 {
00024   using netgen::Ng_Element;
00025   using netgen::Ng_Point;
00026   using netgen::Ng_Node;
00027 
00028   using netgen::Ng_GetPoint;
00029   using netgen::Ng_GetElement;
00030   using netgen::Ng_GetNode;
00031 
00032 
00037   inline ELEMENT_TYPE ConvertElementType (NG_ELEMENT_TYPE type)
00038   {
00039     switch (type)
00040       {
00041       case NG_PNT: 
00042         return ET_POINT;
00043 
00044       case NG_SEGM: case NG_SEGM3:
00045         return ET_SEGM;
00046 
00047       case NG_TRIG: case NG_TRIG6: 
00048         return ET_TRIG;
00049         
00050       case NG_QUAD: case NG_QUAD6:
00051         return ET_QUAD;
00052         
00053       case NG_TET: case NG_TET10:
00054         return ET_TET;
00055 
00056       case NG_PRISM: case NG_PRISM12:
00057         return ET_PRISM;
00058 
00059       case NG_PYRAMID:
00060         return ET_PYRAMID;
00061 
00062       case NG_HEX:
00063         return ET_HEX;
00064       }
00065     throw Exception ("Netgen2NgS type conversion: Unhandled element type");
00066   }
00067 
00068 
00069 
00070 
00071 
00084   class NGS_DLL_HEADER MeshAccess : public BaseStatusHandler
00085   {
00086 
00089     int dim;
00090   
00092     int nnodes[4];
00093 
00094     // number of nodes of co-dimension i 
00095     // these are NC, NF, NE, NV  in 3D,
00096     // and NF, NE, NV, undef, in 2D
00097     int nnodes_cd[4];
00098 
00099 
00101     int nelements[4];  
00103     int nelements_cd[4];
00104 
00106     int nlevels;
00107 
00109     int ndomains;
00110 
00112     int nboundaries;
00113 
00114   public:
00116     MeshAccess ();
00118     virtual ~MeshAccess ();
00119 
00121     int GetDimension() const { return dim; }  
00122 
00124     int GetNP() const  { return Ng_GetNP(); }
00125 
00127     int GetNV() const  { return nnodes[0]; }  
00128 
00130     int GetNE() const  { return nelements_cd[0]; }  
00131 
00133     int GetNSE() const { return nelements_cd[1]; }  
00134 
00136     int GetNEdges() const { return nnodes[1]; }     
00137 
00139     int GetNFaces() const { return nnodes[2]; }    
00140 
00141 
00143     int GetNDomains () const  { return ndomains; }
00144 
00146     int GetNBoundaries () const { return nboundaries; }
00147 
00149     template <int D>
00150     void GetPoint (int pi, Vec<D> & p) const
00151     { 
00152       Ng_Point pt = Ng_GetPoint (pi);
00153       for (int j = 0; j < D; j++) p(j) = pt[j];
00154     }
00155 
00157     template <int D>
00158     Vec<D> GetPoint (int pi) const
00159     { 
00160       Vec<D> p;
00161       Ng_Point pt = Ng_GetPoint (pi);
00162       for (int j = 0; j < D; j++) p(j) = pt[j];
00163       return p;
00164     }
00165 
00167     ELEMENT_TYPE GetElType (int elnr) const
00168     { return ConvertElementType (GetElement(elnr).GetType()); }
00169 
00171     ELEMENT_TYPE GetSElType (int elnr) const
00172     { return ConvertElementType (GetSElement(elnr).GetType()); }
00173 
00175     int GetElIndex (int elnr) const
00176     { return Ng_GetElementIndex (elnr+1)-1; }
00177 
00179     int GetSElIndex (const int elnr) const
00180     { return Ng_GetSurfaceElementIndex (elnr+1)-1; }
00181 
00183     void SetElIndex (int elnr, int index) const
00184     { Ng_SetElementIndex (elnr+1,index+1); }
00185 
00187     string GetElMaterial (int elnr) const
00188     { return Ng_GetElementMaterial (elnr+1); }
00189 
00191     string GetDomainMaterial (int domnr) const
00192     { return Ng_GetDomainMaterial(domnr+1); }
00193       
00194 
00196     string GetSElBCName (int selnr) const
00197     { return Ng_GetSurfaceElementBCName (selnr+1); }
00198 
00200     string GetBCNumBCName (const int bcnr) const
00201     { return Ng_GetBCNumBCName(bcnr); }
00202 
00203 
00204 
00206     int GetSElSurfaceNumber (const int elnr) const
00207     { return Ng_GetSurfaceElementSurfaceNumber (elnr+1)-1; }
00208 
00210     int GetSElFDNumber (const int elnr) const
00211     { return Ng_GetSurfaceElementFDNumber (elnr+1)-1; }
00212 
00213 
00214 
00217     void GetSElNeighbouringDomains(const int elnr, int & in, int & out) const
00218     { 
00219       Ng_GetSurfaceElementNeighbouringDomains(elnr+1,in,out);
00220       //in--; out--;
00221     }
00222   
00223 
00226     void UpdateBuffers();
00227 
00228 
00229 
00230     /*
00231       Nodes are an abstraction for vertices, edges, faces, and cells
00232     */
00233 
00235     int GetNElements (int dim) const { return nelements[dim]; }
00236 
00238     int GetNNodes (NODE_TYPE nt) const { return nnodes[nt]; }  
00239 
00240 
00243     void GetTopologicElement (int elnr, TopologicElement & topel) const;
00244 
00245   
00246 
00247     /*
00248       So moechte ich es gern (JS):
00249 
00250       Anstelle von GetFaceEdges (fanr, edges) etc macht man dann 
00251       GetClosureNodes (NT_FACE, fanr, NodeSet (NT_EDGE), nodes);
00252 
00253       und GetTopologicElement in 3D ruft
00254       GetClosureNodes (NT_CELL, cell, NodeSet (NT_VERTEX, NT_EDGE, NT_FACE, NT_CELL), nodes);
00255     */
00256     // void GetClosureNodes (NODE_TYPE nt, int nodenr, NODE_SET ns, Array<Node> & nodes);
00257 
00264     Ng_Element GetElement (int elnr) const
00265     {
00266       switch (dim)
00267         {
00268         case 1: return Ng_GetElement<1> (elnr);
00269         case 2: return Ng_GetElement<2> (elnr);
00270         case 3: return Ng_GetElement<3> (elnr);
00271         default:
00272           throw Exception ("GetElement, illegal dimension");
00273         }
00274     }
00275 
00282     Ng_Element GetSElement (int elnr) const
00283     {
00284       switch (dim)
00285         {
00286         case 1: return Ng_GetElement<0> (elnr);
00287         case 2: return Ng_GetElement<1> (elnr);
00288         case 3: return Ng_GetElement<2> (elnr);
00289         default:
00290           throw Exception ("GetSElement, illegal dimension");
00291         }
00292     }
00293 
00297     template <int DIM>
00298     Ng_Element GetElement (int elnr) const
00299     {
00300       return Ng_GetElement<DIM> (elnr);
00301     }
00302 
00303     
00309     template <int DIM>
00310     Ng_Node<DIM> GetNode (int nr) const
00311     {
00312       return Ng_GetNode<DIM> (nr);
00313     }
00314 
00315     // von astrid
00316     int GetElNV (int elnr) const;
00317 
00320     void GetElPNums (int elnr, Array<int> & pnums) const;
00321 
00324     void GetSElPNums (int selnr, Array<int> & pnums) const;
00325 
00327     void GetElVertices (int elnr, Array<int> & vnums) const
00328     { vnums = ArrayObject (GetElement(elnr).vertices); }
00329 
00331     void GetSElVertices (int selnr, Array<int> & vnums) const
00332     { vnums = ArrayObject (GetSElement(selnr).vertices); }
00333 
00335     void GetElEdges (int elnr, Array<int> & ednums) const
00336     { ednums = ArrayObject (GetElement(elnr).edges); }
00337 
00338     // returns edge numbers and edge orientation of an element. (old style function)
00339     void GetElEdges (int elnr, Array<int> & ednums, Array<int> & orient) const;
00340 
00342     void GetSElEdges (int selnr, Array<int> & ednums) const
00343     { ednums = ArrayObject (GetSElement(selnr).edges); }
00344 
00345     // returns edge numbers and edge orientation of an element. (old style function)
00346     void GetSElEdges (int selnr, Array<int> & ednums, Array<int> & orient) const;
00347 
00349     void GetElFaces (int elnr, Array<int> & fnums) const
00350     { fnums = ArrayObject (GetElement(elnr).faces); }
00351 
00352     // returns face numbers and face orientation of an element. (old style function)
00353     void GetElFaces (int elnr, Array<int> & fnums, Array<int> & orient) const;
00354 
00356     int GetSElFace (int selnr) const;
00357 
00358     // returns face number and orientation of surface element
00359     void GetSElFace (int selnr, int & fnum, int & orient) const;
00360 
00362     void GetFacePNums (int fnr, Array<int> & pnums) const;
00364     void GetEdgePNums (int enr, int & pn1, int & pn2) const;
00366     void GetEdgePNums (int enr, Array<int> & pnums) const;
00368     void GetEdgeElements (int enr, Array<int> & elnums) const;
00370     void GetFaceEdges (int fnr, Array<int> & edges) const;
00372     void GetFaceElements (int fnr, Array<int> & elnums) const;
00374     void GetSegmentPNums (int snr, Array<int> & pnums) const;
00376     int GetSegmentIndex (int snr) const;
00377     
00380     int GetNFacets() const { return nnodes_cd[1]; } 
00382     void GetElFacets (int elnr, Array<int> & fnums) const;
00384     void GetSElFacets (int selnr, Array<int> & fnums) const;
00386     void GetFacetPNums (int fnr, Array<int> & pnums) const;
00388     ELEMENT_TYPE GetFacetType (int fnr) const;
00390     void GetFacetElements (int fnr, Array<int> & elnums) const
00391     { (dim == 2) ? GetEdgeElements(fnr, elnums) : GetFaceElements(fnr, elnums); }    
00392 
00393     // void GetVertexElements (int vnr, Array<int> & elnrs) const;
00395     int GetElOrder (int enr) const
00396     { return Ng_GetElementOrder (enr+1); } 
00398     INT<3> GetElOrders (int enr) const
00399     { 
00400       INT<3> eo; 
00401       Ng_GetElementOrders(enr+1,&eo[0],&eo[1],&eo[2]); 
00402       return eo; 
00403     } 
00405     void SetElOrder (int enr, int order) const
00406     { Ng_SetElementOrder (enr+1,order); }
00408     void SetElOrders (int enr, int ox, int oy, int oz) const
00409     { Ng_SetElementOrders (enr+1, ox,oy,oz); }
00411     int GetSElOrder (int enr) const
00412     { return Ng_GetSurfaceElementOrder (enr+1); } 
00414     INT<2> GetSElOrders (int enr) const
00415     { 
00416       INT<2> eo; 
00417       Ng_GetSurfaceElementOrders(enr+1,&eo[0],&eo[1]); 
00418       return eo; 
00419     }
00421     void SetSElOrder (int enr, int order) const
00422     { Ng_SetSurfaceElementOrder (enr+1,order); }
00424     void SetSElOrders (int enr, int ox, int oy) const
00425     { Ng_SetSurfaceElementOrders (enr+1, ox,oy); }
00426     
00427 
00429     double ElementVolume (int elnr) const;
00431     double SurfaceElementVolume (int selnr) const;
00432       
00433 
00435     int GetNLevels() const { return nlevels; }  
00437     void GetParentNodes (int pi, int * parents) const
00438     { 
00439       Ng_GetParentNodes (pi+1, parents);
00440       parents[0]--; parents[1]--; 
00441     }
00443     int GetParentElement (int ei) const
00444     { return Ng_GetParentElement (ei+1)-1; }
00446     int GetParentSElement (int ei) const
00447     { return Ng_GetParentSElement (ei+1)-1; }
00448   
00450     int GetClusterRepVertex (int pi) const
00451     { return Ng_GetClusterRepVertex (pi+1)-1; }
00453     int GetClusterRepEdge (int pi) const
00454     { return Ng_GetClusterRepEdge (pi+1)-1; }
00456     int GetClusterRepFace (int pi) const
00457     { return Ng_GetClusterRepFace (pi+1)-1; }
00459     int GetClusterRepElement (int pi) const
00460     { return Ng_GetClusterRepElement (pi+1)-1; }
00461     
00462     
00466     ngfem::ElementTransformation & GetTrafo (int elnr, bool boundary, LocalHeap & lh) const;
00467 
00468 
00469     // (old style optimization)
00470     void SetPointSearchStartElement(const int el) const;
00471 
00472 
00473     
00474     int FindElementOfPoint (FlatVector<double> point,
00475                             IntegrationPoint & ip, 
00476                             bool build_searchtree,
00477                             const Array<int> * const indices = NULL) const;
00478     int FindElementOfPoint (FlatVector<double> point,
00479                             ngfem::IntegrationPoint & ip, 
00480                             bool build_searchtree,
00481                             const int index) const;
00482     int FindSurfaceElementOfPoint (FlatVector<double> point,
00483                                    ngfem::IntegrationPoint & ip, 
00484                                    bool build_searchtree,
00485                                    const Array<int> * const indices = NULL) const;
00486     int FindSurfaceElementOfPoint (FlatVector<double> point,
00487                                    ngfem::IntegrationPoint & ip, 
00488                                    bool build_searchtree,
00489                                    const int index) const;
00490 
00492     bool IsElementCurved (int elnr) const
00493     { return bool (Ng_IsElementCurved (elnr+1)); }
00494     
00495     
00496     void GetPeriodicVertices ( Array<ngstd::INT<2> > & pairs) const;
00497     int GetNPairsPeriodicVertices () const;
00498     void GetPeriodicVertices (int idnr, Array<ngstd::INT<2> > & pairs) const;
00499     int GetNPairsPeriodicVertices (int idnr) const;  
00500 
00501     void GetPeriodicEdges ( Array<ngstd::INT<2> > & pairs) const;
00502     int GetNPairsPeriodicEdges () const;
00503     void GetPeriodicEdges (int idnr, Array<ngstd::INT<2> > & pairs) const;
00504     int GetNPairsPeriodicEdges (int idnr) const;  
00505 
00506 
00507     virtual void PushStatus (const char * str) const;
00508     virtual void PopStatus () const;
00509     virtual void SetThreadPercentage (double percent) const;
00510     virtual void GetStatus (string & str, double & percent) const;
00511 
00512     virtual void SetTerminate(void) const;
00513     virtual void UnSetTerminate(void) const;
00514     virtual bool ShouldTerminate(void) const;
00515   
00517     void GetVertexElements( int vnr, Array<int>& elems) const;
00518 
00519     void GetVertexSurfaceElements( int vnr, Array<int>& elems) const;
00520 
00521 
00522   private:
00523     Array<bool> higher_integration_order;
00524   public:
00525     void SetHigherIntegrationOrder(int elnr);
00526     void UnSetHigherIntegrationOrder(int elnr);
00527 
00528 
00529 
00530     void LoadMeshFromString(const string & str)
00531     {
00532       Ng_LoadMeshFromString(const_cast<char*>(str.c_str()));
00533     }
00534 
00535 
00536 
00537     // void PrecomputeGeometryData(int intorder);
00538 
00539     void InitPointCurve(double red = 1, double green = 0, double blue = 0) const;
00540     void AddPointCurvePoint(const Vec<3> & point) const;
00541 
00542 
00543     MPI_Comm GetCommunicator () const { return ngs_comm; }
00544 
00549     void GetDistantProcs (Node node, Array<int> & procs) const;
00550 
00555     int GetGlobalNodeNum (Node node) const;
00556 
00557 
00559     template <typename T>
00560     void AllReduceNodalData (NODE_TYPE nt, Array<T> & data, MPI_Op op) const;
00561   };
00562 
00563 
00564 
00565 
00571   class ProgressOutput
00572   {
00573     const MeshAccess & ma;
00574     string task;
00575     int total;
00576     double prevtime;
00577     bool is_root;
00578 
00579     int cnt;
00580     bool done_called;
00581   public:
00582     ProgressOutput (const MeshAccess & ama,
00583                     string atask, int atotal);
00584     ~ProgressOutput ();
00585 
00586     void Update ();
00587     void Update (int nr);
00588     void Done();
00589   };
00590 
00591 
00592 
00594   template <typename T>
00595   void AllReduceNodalData (NODE_TYPE nt, Array<T> & data, MPI_Op op, const MeshAccess & ma)
00596   {
00597     ma.AllReduceNodalData (nt, data, op);
00598   }
00599 
00600 
00601 #ifdef PARALLEL
00602 
00603   template <typename T>
00604   void MeshAccess::
00605   AllReduceNodalData (NODE_TYPE nt, Array<T> & data, MPI_Op op) const
00606   {
00607     MPI_Comm comm = GetCommunicator();
00608     MPI_Datatype type = MyGetMPIType<T>();
00609 
00610     int ntasks = MyMPI_GetNTasks (comm);
00611     if (ntasks <= 1) return;
00612 
00613     Array<int> cnt(ntasks), distp;
00614     cnt = 0;
00615     for (int i = 0; i < GetNNodes(nt); i++)
00616       {
00617         GetDistantProcs (Node (nt, i), distp);
00618         for (int j = 0; j < distp.Size(); j++)
00619           cnt[distp[j]]++;
00620       }
00621 
00622     Table<T> dist_data(cnt), recv_data(cnt);
00623 
00624     cnt = 0;
00625     for (int i = 0; i < GetNNodes(nt); i++)
00626       {
00627         GetDistantProcs (Node (nt, i), distp);
00628         for (int j = 0; j < distp.Size(); j++)
00629           dist_data[distp[j]][cnt[distp[j]]++] = data[i];
00630       }
00631 
00632     Array<MPI_Request> requests;
00633     for (int i = 0; i < ntasks; i++)
00634       if (cnt[i])
00635         {
00636           requests.Append (MyMPI_ISend (dist_data[i], i, MPI_TAG_SOLVE, comm));
00637           requests.Append (MyMPI_IRecv (recv_data[i], i, MPI_TAG_SOLVE, comm));
00638         }
00639     MyMPI_WaitAll (requests);
00640     
00641     cnt = 0;
00642     for (int i = 0; i < data.Size(); i++)
00643       {
00644         GetDistantProcs (Node (nt, i), distp);
00645         for (int j = 0; j < distp.Size(); j++)
00646           MPI_Reduce_local (&recv_data[distp[j]][cnt[distp[j]]++],
00647                             &data[i], 1, type, op);
00648       }
00649   }
00650 
00651 
00652 #else
00653 
00654   template <typename T>
00655   void MeshAccess::
00656   AllReduceNodalData (NODE_TYPE nt, Array<T> & data, MPI_Op op) const { ; }
00657 
00658 #endif
00659 
00660 }
00661 
00662 #endif