NGSolve
4.9
|
00001 #ifndef FILE_FESPACE 00002 #define FILE_FESPACE 00003 00004 /*********************************************************************/ 00005 /* File: fespace.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 25. Mar. 2000 */ 00008 /*********************************************************************/ 00009 00010 00011 namespace ngcomp 00012 { 00013 00014 /* 00015 Finite Element Space 00016 */ 00017 00018 00023 enum TRANSFORM_TYPE { TRANSFORM_MAT_LEFT = 1, 00024 TRANSFORM_MAT_RIGHT = 2, 00025 TRANSFORM_MAT_LEFT_RIGHT = 3, 00026 TRANSFORM_RHS = 4, 00027 TRANSFORM_SOL = 8 }; 00035 enum COUPLING_TYPE { UNUSED_DOF = 0, 00036 LOCAL_DOF = 1, 00037 INTERFACE_DOF = 2, 00038 NONWIREBASKET_DOF = 3, 00039 WIREBASKET_DOF = 4, 00040 EXTERNAL_DOF = 6, 00041 ANY_DOF = 7 00042 }; 00043 00044 00045 00046 00052 class NGS_DLL_HEADER FESpace : public NGS_Object 00053 { 00054 protected: 00056 int order; 00058 int dimension; 00060 bool iscomplex; 00061 00063 bool dgjumps; 00064 00066 bool print; 00067 00069 ngmg::Prolongation *prol; 00071 int level_updated; 00072 00074 Array<int> definedon; 00076 Array<int> definedonbound; 00077 00079 BitArray dirichlet_boundaries; 00080 00082 BitArray dirichlet_dofs; 00083 BitArray free_dofs; 00084 BitArray external_free_dofs; 00085 00086 00087 Array<bool> dirichlet_vertex; 00088 Array<bool> dirichlet_edge; 00089 Array<bool> dirichlet_face; 00090 00091 00093 FiniteElement * tet; 00095 FiniteElement * prism; 00097 FiniteElement * pyramid; 00099 FiniteElement * hex; 00101 FiniteElement * trig; 00103 FiniteElement * quad; 00105 FiniteElement * segm; 00107 FiniteElement * point; 00108 00110 DifferentialOperator * evaluator; 00112 DifferentialOperator * boundary_evaluator; 00113 00115 BilinearFormIntegrator * integrator; 00117 BilinearFormIntegrator * boundary_integrator; 00118 00119 00120 00122 FESpace * low_order_space; 00124 bool is_low_order_space; 00125 00127 Array<bool> directsolverclustered; 00128 00129 Array<string> directsolvermaterials; 00130 00131 mutable Array<int> adddirectsolverdofs; 00132 00133 Array<int> directvertexclusters; 00134 Array<int> directedgeclusters; 00135 Array<int> directfaceclusters; 00136 Array<int> directelementclusters; 00137 00138 Table<int> * element_coloring; 00139 Array<COUPLING_TYPE> ctofdof; 00140 00141 ParallelDofs * paralleldofs; 00142 00143 public: 00152 FESpace (const MeshAccess & ama, const Flags & flags, 00153 bool checkflags = false); 00155 virtual ~FESpace (); 00156 00158 // virtual void Update(); 00159 00161 virtual void Update(LocalHeap & lh); 00162 00164 virtual void FinalizeUpdate(LocalHeap & lh); 00165 00167 int GetLevelUpdated() const { return level_updated; } 00168 00169 const Table<int> & ElementColoring() const { return *element_coloring; } 00170 00172 virtual void PrintReport (ostream & ost); 00173 00175 int GetOrder () const { return order; } 00176 00178 int GetDimension () const { return dimension; } 00179 00181 bool IsComplex () const { return iscomplex; } 00182 00183 00185 virtual int GetNDof () const = 0; 00187 virtual int GetNDofLevel (int level) const; 00188 00190 const FiniteElement & GetFE (int elnr, bool boundary, LocalHeap & lh) const 00191 { 00192 return boundary ? GetSFE(elnr, lh) : GetFE(elnr, lh); 00193 } 00194 00196 virtual const FiniteElement & GetFE (int elnr, LocalHeap & lh) const; 00197 00199 virtual void GetDofNrs (int elnr, Array<int> & dnums) const = 0; 00200 00202 void GetDofNrs (int elnr, bool boundary, Array<int> & dnums) const 00203 { 00204 if (boundary) 00205 GetSDofNrs (elnr, dnums); 00206 else 00207 GetDofNrs (elnr, dnums); 00208 } 00209 00210 00212 virtual void GetDofCouplingTypes (int elnr, Array<COUPLING_TYPE> & dnums) const; 00213 00215 virtual COUPLING_TYPE GetDofCouplingType (int dof) const; 00216 00218 void GetDofNrs (int elnr, Array<int> & dnums, COUPLING_TYPE ctype) const; 00219 00221 virtual void GetNodeDofNrs (NODE_TYPE nt, int nr, Array<int> & dnums) const; 00223 // virtual int GetNLowOrderNodeDofs ( NODE_TYPE nt ) const; 00224 // { return lodofs_per_node[nt]; } 00225 00227 virtual void GetVertexDofNrs (int vnr, Array<int> & dnums) const; 00229 virtual void GetEdgeDofNrs (int ednr, Array<int> & dnums) const; 00231 virtual void GetFaceDofNrs (int fanr, Array<int> & dnums) const; 00233 virtual void GetInnerDofNrs (int elnr, Array<int> & dnums) const; 00234 00235 virtual bool UsesDGCoupling () const {return dgjumps;}; 00236 00238 virtual const FiniteElement & GetSFE (int selnr, LocalHeap & lh) const; 00240 virtual void GetSDofNrs (int selnr, Array<int> & dnums) const = 0; 00241 00242 00244 bool DefinedOn (int domnr) const 00245 { return !definedon.Size() || definedon[domnr]; } 00247 bool DefinedOnBoundary (int bnr) const 00248 {return !definedonbound.Size() || definedonbound[bnr]; } 00249 00251 void SetDefinedOn (const BitArray & defon); 00253 void SetDefinedOnBoundary (const BitArray & defon); 00254 00256 void SetDirichletBoundaries (const BitArray & dirbnds); 00258 const FiniteElement & GetFE (ELEMENT_TYPE type) const; 00259 00261 FESpace & LowOrderFESpace () { return *low_order_space; } 00263 const FESpace & LowOrderFESpace () const { return *low_order_space; } 00265 void SetLowOrderSpace (bool los) { is_low_order_space = los; } 00267 bool IsLowOrderSpace () const { return is_low_order_space; } 00268 00270 virtual const BitArray * GetFreeDofs (bool external = false) const; 00272 bool IsDirichletDof (int i) const 00273 { return dirichlet_dofs.Size() && dirichlet_dofs[i]; } 00274 00275 bool IsDirichletBoundary (int i) const 00276 { return dirichlet_boundaries.Size() && dirichlet_boundaries[i]; } 00277 00279 bool IsDirichletVertex (int i) const { return dirichlet_vertex.Size() && dirichlet_vertex[i]; } 00281 bool IsDirichletEdge (int i) const { return dirichlet_edge.Size() && dirichlet_edge[i]; } 00283 bool IsDirichletFace (int i) const { return dirichlet_face.Size() && dirichlet_face[i]; } 00284 00285 void GetFilteredDofs(COUPLING_TYPE doffilter, BitArray & output, bool freedofsonly=true) const; 00287 virtual Table<int> * CreateSmoothingBlocks (const Flags & flags) const; 00289 virtual Array<int> * CreateDirectSolverClusters (const Flags & flags) const 00290 { return 0; } 00291 00292 virtual void AddDirectSolverClusterDof(int dn) const 00293 { adddirectsolverdofs.Append(dn); } 00294 00295 virtual Array<int> & DirectVertexClusters(void) 00296 { return directvertexclusters; } 00297 virtual Array<int> & DirectEdgeClusters(void) 00298 { return directedgeclusters; } 00299 virtual Array<int> & DirectFaceClusters(void) 00300 { return directfaceclusters; } 00301 virtual Array<int> & DirectElementClusters(void) 00302 { return directelementclusters; } 00303 00304 00305 void TransformMat (int elnr, bool boundary, 00306 const SliceMatrix<double> & mat, TRANSFORM_TYPE type) const 00307 { 00308 VTransformMR (elnr, boundary, mat, type); 00309 } 00310 00311 void TransformMat (int elnr, bool boundary, 00312 const SliceMatrix<Complex> & mat, TRANSFORM_TYPE type) const 00313 { 00314 VTransformMC (elnr, boundary, mat, type); 00315 } 00316 00317 00318 00319 void TransformVec (int elnr, bool boundary, 00320 const FlatVector<double> & vec, TRANSFORM_TYPE type) const 00321 { 00322 VTransformVR (elnr, boundary, vec, type); 00323 } 00324 00325 void TransformVec (int elnr, bool boundary, 00326 const FlatVector<Complex> & vec, TRANSFORM_TYPE type) const 00327 { 00328 VTransformVC (elnr, boundary, vec, type); 00329 } 00330 00331 00332 template < int S, class T > 00333 void TransformVec (int elnr, bool boundary, 00334 const FlatVector< Vec<S,T> >& vec, TRANSFORM_TYPE type) const; 00335 00336 00337 00338 virtual void VTransformMR (int elnr, bool boundary, 00339 const SliceMatrix<double> & mat, TRANSFORM_TYPE type) const 00340 { ; } 00341 virtual void VTransformMC (int elnr, bool boundary, 00342 const SliceMatrix<Complex> & mat, TRANSFORM_TYPE type) const 00343 { ; } 00344 00345 00346 virtual void VTransformVR (int elnr, bool boundary, 00347 const FlatVector<double> & vec, TRANSFORM_TYPE type) const 00348 { ; } 00349 virtual void VTransformVC (int elnr, bool boundary, 00350 const FlatVector<Complex> & vec, TRANSFORM_TYPE type) const 00351 { ; } 00352 00353 00355 virtual const ngmg::Prolongation * GetProlongation () const 00356 { return prol; } 00358 // void SetProlongation (ngmg::Prolongation * aprol) 00359 // { prol = aprol; } 00360 00361 00363 const DifferentialOperator * GetEvaluator (bool boundary = false) const 00364 { 00365 if (boundary) 00366 return boundary_evaluator; 00367 else 00368 return evaluator; 00369 } 00370 00371 00373 const BilinearFormIntegrator * GetIntegrator (bool boundary = false) const 00374 { 00375 if (boundary) 00376 return boundary_integrator; 00377 else 00378 return integrator; 00379 } 00380 00381 00383 const DifferentialOperator * GetBoundaryEvaluator () const 00384 { return boundary_evaluator; } 00385 const BilinearFormIntegrator * GetBoundaryIntegrator () const 00386 { return boundary_integrator; } 00387 00388 00390 Array<SpecialElement*> specialelements; 00391 00392 void AppendSpecialElement (SpecialElement * spel) 00393 { specialelements.Append (spel); } 00394 00395 const Array<SpecialElement*> & GetSpecialElements() const {return specialelements;} 00396 00397 00398 ParallelDofs & GetParallelDofs () const { return *paralleldofs; } 00399 virtual void UpdateParallelDofs (); 00400 00402 bool IsParallel() const; 00403 00405 int GetNDofGlobal() const; 00406 00407 virtual int GetRelOrder() const 00408 { 00409 cout << "virtual GetRelOrder called for FiniteElementSpace, not available ! " << endl; 00410 return 0; 00411 } 00412 00413 virtual bool VarOrder() const { return 0; } 00414 00415 bool timing; 00416 void Timing () const; 00417 }; 00418 00419 00420 00421 00422 00423 00428 class NGS_DLL_HEADER NodalFESpace : public FESpace 00429 { 00431 Array<int> ndlevel; 00432 00433 public: 00434 00436 NodalFESpace (const MeshAccess & ama, const Flags & flags, bool parseflags=false); 00438 virtual ~NodalFESpace (); 00439 00441 virtual string GetClassName () const 00442 { 00443 return "NodalFESpace"; 00444 } 00445 00447 virtual void Update (LocalHeap & lh); 00449 virtual int GetNDof () const; 00451 virtual int GetNDofLevel (int level) const; 00453 virtual void GetDofNrs (int elnr, Array<int> & dnums) const; 00455 virtual void GetSDofNrs (int selnr, Array<int> & dnums) const; 00456 00457 00458 virtual void GetVertexDofNrs (int vnr, Array<int> & dnums) const; 00459 virtual void GetEdgeDofNrs (int ednr, Array<int> & dnums) const; 00460 virtual void GetFaceDofNrs (int fanr, Array<int> & dnums) const; 00461 virtual void GetInnerDofNrs (int elnr, Array<int> & dnums) const; 00462 00463 virtual Array<int> * CreateDirectSolverClusters (const Flags & flags) const; 00464 }; 00465 00466 00467 00468 00469 00470 00472 class NGS_DLL_HEADER NonconformingFESpace : public FESpace 00473 { 00475 Array<int> ndlevel; 00476 00477 public: 00478 NonconformingFESpace (const MeshAccess & ama, const Flags & flags, bool parseflags=false); 00479 virtual ~NonconformingFESpace (); 00480 00481 virtual string GetClassName () const 00482 { return "Nonconforming FESpace"; } 00483 00485 virtual void Update(LocalHeap & lh); 00487 virtual int GetNDof () const; 00489 virtual void GetDofNrs (int elnr, Array<int> & dnums) const; 00491 virtual void GetSDofNrs (int selnr, Array<int> & dnums) const; 00492 }; 00493 00494 00495 00496 00497 00498 00499 00501 class NGS_DLL_HEADER ElementFESpace : public FESpace 00502 { 00504 Array<int> ndlevel; 00505 int n_el_dofs; 00506 public: 00508 ElementFESpace (const MeshAccess & ama, const Flags& flags, bool parseflags=false); 00509 00511 ~ElementFESpace (); 00512 00513 virtual string GetClassName () const 00514 { 00515 return "ElementFESpace"; 00516 } 00517 00519 virtual void Update(LocalHeap & lh); 00520 00522 virtual int GetNDof () const { return ndlevel.Last(); } 00523 00525 virtual void GetDofNrs (int elnr, Array<int> & dnums) const; 00526 00528 virtual int GetNDofLevel (int level) const; 00529 00531 virtual void GetSDofNrs (int selnr, Array<int> & dnums) const; 00532 00533 00534 virtual void GetVertexDofNrs (int vnr, Array<int> & dnums) const 00535 { dnums.SetSize (0); } 00536 virtual void GetEdgeDofNrs (int ednr, Array<int> & dnums) const 00537 { dnums.SetSize (0); } 00538 virtual void GetFaceDofNrs (int fanr, Array<int> & dnums) const 00539 { dnums.SetSize (0); } 00540 virtual void GetInnerDofNrs (int elnr, Array<int> & dnums) const 00541 { GetDofNrs (elnr, dnums); } 00542 }; 00543 00544 00545 00546 00547 00549 class NGS_DLL_HEADER SurfaceElementFESpace : public FESpace 00550 { 00552 Array<int> ndlevel; 00553 int n_el_dofs; 00554 public: 00556 SurfaceElementFESpace (const MeshAccess & ama, const Flags& flags, 00557 bool checkflags = false); 00558 00560 ~SurfaceElementFESpace (); 00561 00563 virtual string GetClassName() const 00564 { return "SurfaceElement"; } 00565 00567 virtual void Update(LocalHeap & lh); 00568 00570 virtual int GetNDof () const { return ndlevel.Last(); } 00571 00573 virtual const FiniteElement & GetFE (int elnr, LocalHeap & lh) const; 00574 00576 virtual void GetDofNrs (int elnr, Array<int> & dnums) const; 00577 00579 virtual int GetNDofLevel (int level) const; 00580 00582 virtual void GetSDofNrs (int selnr, Array<int> & dnums) const; 00583 }; 00584 00585 00586 00587 00588 00589 00590 00592 class NGS_DLL_HEADER CompoundFESpace : public FESpace 00593 { 00594 protected: 00596 Array<FESpace*> spaces; 00598 Array<int> cummulative_nd; 00600 Array<int> ndlevel; 00601 public: 00604 CompoundFESpace (const MeshAccess & ama, 00605 const Flags & flags, bool parseflags = false); 00608 CompoundFESpace (const MeshAccess & ama, 00609 const Array<FESpace*> & aspaces, 00610 const Flags & flags, bool parseflags = false); 00611 00614 virtual ~CompoundFESpace (); 00615 00617 void AddSpace (FESpace * fes); 00618 00620 virtual string GetClassName () const 00621 { 00622 return "CompoundFESpace"; 00623 } 00624 00626 virtual void Update(LocalHeap & lh); 00628 virtual void FinalizeUpdate(LocalHeap & lh); 00629 00631 virtual void UpdateCouplingDofArray(); 00632 00634 virtual int GetNDof () const { return ndlevel.Last(); } 00636 virtual int GetNDofLevel (int level) const { return ndlevel[level]; } 00637 00638 /* 00639 // returns start and end points of dofs corresponding to space "spacenr" 00640 // first space: spacenr = 0 00641 int GetStorageStart(int spacenr) const 00642 { return cummulative_nd[spacenr]; } 00644 int GetStorageEnd(int spacenr) const 00645 { return cummulative_nd[spacenr+1]; } 00646 */ 00647 00648 IntRange GetRange (int spacenr) const 00649 { 00650 return IntRange(cummulative_nd[spacenr], cummulative_nd[spacenr+1]); 00651 } 00652 00654 const FESpace * operator[] (int i) const { return spaces[i]; } 00656 FESpace * operator[] (int i) { return const_cast<FESpace*> (spaces[i]); } 00657 00659 virtual const FiniteElement & GetFE (int elnr, LocalHeap & lh) const; 00661 virtual const FiniteElement & GetSFE (int selnr, LocalHeap & lh) const; 00663 virtual void GetDofNrs (int elnr, Array<int> & dnums) const; 00665 virtual void GetSDofNrs (int selnr, Array<int> & dnums) const; 00667 virtual void GetVertexDofNrs (int vnr, Array<int> & dnums) const; 00668 virtual void GetEdgeDofNrs (int ednr, Array<int> & dnums) const; 00669 virtual void GetFaceDofNrs (int fanr, Array<int> & dnums) const; 00670 virtual void GetInnerDofNrs (int elnr, Array<int> & dnums) const; 00671 00672 00673 template <class MAT> NGS_DLL_HEADER 00674 void TransformMat (int elnr, bool boundary, 00675 MAT & mat, TRANSFORM_TYPE tt) const; 00676 00677 template <class VEC> NGS_DLL_HEADER 00678 void TransformVec (int elnr, bool boundary, 00679 VEC & vec, TRANSFORM_TYPE tt) const; 00680 00681 virtual void VTransformMR (int elnr, bool boundary, 00682 const SliceMatrix<double> & mat, TRANSFORM_TYPE tt) const; 00683 virtual void VTransformMC (int elnr, bool boundary, 00684 const SliceMatrix<Complex> & mat, TRANSFORM_TYPE tt) const; 00685 virtual void VTransformVR (int elnr, bool boundary, 00686 const FlatVector<double> & vec, TRANSFORM_TYPE tt) const; 00687 virtual void VTransformVC (int elnr, bool boundary, 00688 const FlatVector<Complex> & vec, TRANSFORM_TYPE tt) const; 00689 00691 inline int GetNSpaces () const { return spaces.Size(); } 00692 }; 00693 00694 00695 00696 00697 00699 class NGS_DLL_HEADER FESpaceClasses 00700 { 00701 public: 00704 struct FESpaceInfo 00705 { 00707 string name; 00709 FESpace* (*creator)(const MeshAccess & ma, const Flags & flags); 00711 FESpaceInfo (const string & aname, 00712 FESpace* (*acreator)(const MeshAccess & ma, const Flags & flags)) 00713 : name(aname), creator(acreator) {;} 00714 }; 00715 private: 00716 Array<FESpaceInfo*> fesa; 00717 00718 public: 00720 FESpaceClasses() { ; } 00722 ~FESpaceClasses(); 00723 00725 void AddFESpace (const string & aname, 00726 FESpace* (*acreator)(const MeshAccess & ma, const Flags & flags)); 00727 00729 const Array<FESpaceInfo*> & GetFESpaces() { return fesa; } 00730 00732 const FESpaceInfo * GetFESpace(const string & name); 00733 00735 void Print (ostream & ost) const; 00736 }; 00737 00739 extern NGS_DLL_HEADER FESpaceClasses & GetFESpaceClasses (); 00740 00742 extern NGS_DLL_HEADER FESpace * CreateFESpace (const string & type, 00743 const MeshAccess & ma, 00744 const Flags & flags); 00745 00746 00751 template <typename FES> 00752 class RegisterFESpace 00753 { 00754 public: 00756 RegisterFESpace (string label) 00757 { 00758 GetFESpaceClasses().AddFESpace (label, Create); 00759 // cout << "register fespace '" << label << "'" << endl; 00760 } 00761 00763 static FESpace * Create (const MeshAccess & ma, const Flags & flags) 00764 { 00765 return new FES (ma, flags); 00766 } 00767 }; 00768 00769 00770 00771 00772 00773 00774 00775 00776 00777 00778 00779 00780 00781 00782 00783 #ifdef PARALLEL 00784 00785 class ParallelMeshDofs : public ParallelDofs 00786 { 00787 const MeshAccess & ma; 00788 Array<Node> dofnodes; 00789 public: 00790 ParallelMeshDofs (const MeshAccess & ama, const Array<Node> & adofnodes, 00791 int dim = 1, bool iscomplex = false); 00792 00793 const MeshAccess & GetMeshAccess() const { return ma; } 00794 const Array<Node> & GetDofNodes() const { return dofnodes; } 00795 }; 00796 00797 #else 00798 00799 00800 class ParallelMeshDofs : public ParallelDofs 00801 { 00802 public: 00803 ParallelMeshDofs (const MeshAccess & ama, const Array<Node> & adofnodes, 00804 int dim = 1, bool iscomplex = false) 00805 { ndof = adofnodes.Size(); } 00806 }; 00807 00808 #endif 00809 00810 } 00811 00812 00813 00814 00815 00816 #ifdef PARALLEL 00817 namespace ngstd 00818 { 00819 template<> 00820 class MPI_Traits<ngcomp::COUPLING_TYPE> 00821 { 00822 public: 00824 static MPI_Datatype MPIType () 00825 { 00826 if (sizeof(ngcomp::COUPLING_TYPE) == sizeof(int)) return MPI_INT; 00827 cout << "please provide MPI_Datatype for COUPLING_TYPE" << endl; 00828 exit(1); 00829 } 00830 }; 00831 } 00832 #endif 00833 00834 00835 #endif