NGSolve  4.9
comp/fespace.hpp
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