NGSolve
4.9
|
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