AFEPack
FEMSpace.h
浏览该文件的文档。
00001 
00011 #ifndef _FEMSpace_h_
00012 #define _FEMSpace_h_
00013 
00014 #include <iostream>
00015 #include <string>
00016 #include <vector>
00017 #include <iterator>
00018 #include <algorithm>
00019 
00020 #include <base/exceptions.h>
00021 #include <lac/vector.h>
00022 #include <lac/sparse_matrix.h>
00023 
00024 #include "Geometry.h"
00025 #include "HGeometry.h"
00026 #include "TemplateElement.h"
00027 #ifdef MULTITHREAD
00028 #include "Thread.h"
00029 #endif
00030 
00031 template <class value_type, int DIM, int DOW, int TDIM> class Element;
00032 template <class value_type, int DIM, int DOW, int TDIM> class FEMSpace;
00033 template <class value_type, int DIM, int DOW, int TDIM, typename Number> class FEMFunction;
00034 
00036 
00043 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM>
00044   class Element
00045   {
00046   public:
00047   enum { dim = DIM, dow = DOW, tdim = TDIM };
00048   typedef value_type value_t;
00049   typedef FEMSpace<value_t,DIM,DOW,TDIM> fe_space_t;
00050   typedef Element<value_t,DIM,DOW,TDIM> element_t;
00051   typedef TemplateElement<value_t,DOW,TDIM> template_t;
00052 
00053   private:
00054   fe_space_t * sp; 
00055   int geometry_index; 
00056   int template_element_index; 
00057   std::vector<int> dof_index; 
00058   std::vector<std::vector<int> > geo_img; 
00059   public:
00060   explicit Element(fe_space_t& = *((fe_space_t *)NULL)); 
00061   Element(const element_t&);  
00062   ~Element(); 
00063   public:
00064   void reinit(fe_space_t&, int, int); 
00065   void reinit(int, int); 
00066   void reinit(fe_space_t&, int, int, const std::vector<int>&); 
00067   element_t& operator=(const element_t&); 
00068   const fe_space_t& femSpace() const; 
00069   fe_space_t& femSpace(); 
00070   Mesh<DIM,DOW>& mesh() const; 
00071   int index() const; 
00072   int& index(); 
00073   GeometryBM& geometry() const; 
00074   void geometry(const GeometryBM&); 
00075   void geometry(int); 
00076   int templateElementIndex() const; 
00077   template_t& templateElement() const; 
00078   void templateElement(int); 
00079   TemplateGeometry<TDIM>& templateGeometry() const; 
00080   const std::vector<int>& dof() const; 
00081   std::vector<int>& dof(); 
00082   unsigned int n_dof() const; 
00083   const std::vector<std::vector<int> >& geometryImage() const; 
00084   std::vector<std::vector<int> >& geometryImage(); 
00085   void buildGeometryImage(); 
00086   void lazyBuildGeometryImage(); 
00087   void buildVertexArray(std::vector<afepack::Point<DOW> >&) const; 
00088   const double ** buildVertexArray() const; 
00090   const typename BasisFunction<value_t,DIM,TDIM>::Identity& basis_function_identity(int i) const; 
00091   value_t basis_function_value(int i, const afepack::Point<DOW>&) const; 
00092   std::vector<value_t> basis_function_value(const afepack::Point<DOW>&) const; 
00093   std::vector<value_t> basis_function_value(int i, const std::vector<afepack::Point<DOW> >&) const; 
00094   std::vector<std::vector<value_t> > basis_function_value(const std::vector<afepack::Point<DOW> >&) const; 
00095   std::vector<value_t> basis_function_gradient(int i, const afepack::Point<DOW>&) const; 
00096   std::vector<std::vector<value_t> > basis_function_gradient(const afepack::Point<DOW>&) const; 
00097   std::vector<std::vector<value_t> > basis_function_gradient(int i, const std::vector<afepack::Point<DOW> >&) const; 
00098   std::vector<std::vector<std::vector<value_t> > > basis_function_gradient(const std::vector<afepack::Point<DOW> >&) const; 
00100 #ifdef QUADRATIC_ELEMENT_SUPPORT
00101   std::vector<std::vector<value_t> > 
00102   basis_function_hesse(int i, 
00103                        const afepack::Point<DOW>&) const; 
00104   std::vector<std::vector<std::vector<value_t> > > 
00105   basis_function_hesse(const afepack::Point<DOW>&) const;
00106   std::vector<std::vector<std::vector<value_t> > >
00107   basis_function_hesse(int i, 
00108                        const std::vector<afepack::Point<DOW> >&) const;
00109   std::vector<std::vector<std::vector<std::vector<value_t> > > >
00110   basis_function_gradient(const std::vector<afepack::Point<DOW> >&) const;
00111 #endif // QUADRATIC_ELEMENT_SUPPORT
00112 
00113   afepack::Point<DOW> local_to_global(const afepack::Point<TDIM>&) const; 
00114   afepack::Point<TDIM> global_to_local(const afepack::Point<DOW>&) const; 
00115   double local_to_global_jacobian(const afepack::Point<TDIM>&) const; 
00116   double global_to_local_jacobian(const afepack::Point<DOW>&) const; 
00117   std::vector<afepack::Point<DOW> > local_to_global(const std::vector<afepack::Point<TDIM> >&) const; 
00118   std::vector<afepack::Point<TDIM> > global_to_local(const std::vector<afepack::Point<DOW> >&) const; 
00119   std::vector<double> local_to_global_jacobian(const std::vector<afepack::Point<TDIM> >&) const; 
00120   std::vector<double> global_to_local_jacobian(const std::vector<afepack::Point<DOW> >&) const; 
00122   std::vector<double> unitOutNormal(const afepack::Point<DOW>&, int) const;
00123   std::vector<std::vector<double> > unitOutNormal(const std::vector<afepack::Point<DOW> >&, int) const;
00124         
00125   const QuadratureInfo<TDIM>& findQuadratureInfo(int) const; 
00127   public:
00128   DeclException1(ExcMeshData, std::string, << "Mesh data uncompatible: " << arg1);
00129 #ifdef __SERIALIZATION__ 
00130   template <class T> T *
00131   new_property(const property_id_t<T>& pid) const {
00132     return femSpace().new_property(*this, pid);
00133   }
00134   template <class T> T *
00135   get_property(const property_id_t<T>& pid) const {
00136     return femSpace().get_property(*this, pid);
00137   }
00138   template <class T> void
00139   free_property(const property_id_t<T>& pid) const {
00140     femSpace().free_property(*this, pid);
00141   }
00142 #endif // __SERIALIZATION__
00143   };
00144 
00153 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM>
00154   class FEMSpace
00155   {
00156   public:
00157   enum { dim = DIM, dow = DOW, tdim = TDIM };
00158 
00159   typedef value_type value_t;
00160   typedef Element<value_t,DIM,DOW,TDIM> element_t;
00161   typedef FEMSpace<value_t,DIM,DOW,TDIM> fe_space_t;
00162   typedef DOFInfo<value_t, DIM, DOW,TDIM> dof_info_t;
00163   typedef typename element_t::template_t template_t;
00164   typedef Mesh<DIM,DOW> mesh_t;
00165   typedef typename mesh_t::bmark_t bmark_t;
00166   typedef typename std::vector<element_t>::iterator ElementIterator;
00167   typedef typename std::vector<element_t>::const_iterator ConstElementIterator;
00168 
00169   private:
00170   mesh_t * msh; 
00171   std::vector<template_t> * tmp_ele; 
00172   std::vector<element_t> ele; 
00173   DegreeOfFreedom df; 
00174   std::vector<dof_info_t> df_in; 
00176   u_int effi_flag;              
00178   public:
00179   explicit FEMSpace(mesh_t& = *((mesh_t *)NULL),
00180                     std::vector<template_t>& 
00181                     = *((std::vector<template_t> *)NULL)); 
00182   FEMSpace(const fe_space_t&); 
00183   virtual ~FEMSpace(); 
00184   public:
00185   fe_space_t& operator=(const fe_space_t&); 
00186   void reinit(mesh_t& = *((mesh_t *)NULL),
00187               std::vector<template_t>& 
00188               = *((std::vector<template_t> *)NULL)); 
00190   const mesh_t& mesh() const {return *msh;}
00191   mesh_t& mesh() {return *msh;}
00192 
00193   const std::vector<template_t>& templateElement() const {return *tmp_ele;}
00194   std::vector<template_t>& templateElement() {return *tmp_ele;}
00195   const template_t& templateElement(int i) const {return (*tmp_ele)[i];}
00196   template_t& templateElement(int i) {return (*tmp_ele)[i];}
00197 
00198   int n_element() {return ele.size();};
00199   const std::vector<element_t>& element() const {return ele;}
00200   std::vector<element_t>& element() {return ele;}
00201   const element_t& element(int i) const {return ele[i];}
00202   element_t& element(int i) {return ele[i];}
00203 
00204   const DegreeOfFreedom& dof() const {return df;}
00205   DegreeOfFreedom& dof() {return df;}
00206 
00207   const std::vector<dof_info_t>& dofInfo() const {return df_in;}
00208   std::vector<dof_info_t>& dofInfo() {return df_in;}
00209   const dof_info_t& dofInfo(int i) const {return df_in[i];}
00210   dof_info_t& dofInfo(int i) {return df_in[i];}
00211 
00212   unsigned int n_dof() const {return dof().n_dof;}
00213   const std::vector<int>& geometryDof(int i, int j) const {return dof().geometry_dof[i][j];}
00214   const std::vector<int>& geometryDof(int i, const Geometry& g) const {return dof().geometry_dof[i][g.index()];}
00215 
00216   const std::vector<DOFIndex>& dofIndex() const {return dof().dof_index;}
00217   const DOFIndex& dofIndex(int i) const {return dof().dof_index[i];}
00218 
00219   const bmark_t& dofBoundaryMark(int i) const {return dofInfo(i).boundary_mark;}
00220   bmark_t& dofBoundaryMark(int i) {return dofInfo(i).boundary_mark;}
00221 
00222   virtual void buildElement(bool lazy = false); 
00224   virtual void buildDofBoundaryMark(); 
00226   void buildDof(); 
00227   void updateDofInterpPoint(); 
00229   ElementIterator beginElement() {return element().begin();};
00230   ConstElementIterator beginElement() const {return element().begin();};
00232   ElementIterator endElement() {return element().end();};
00233   ConstElementIterator endElement() const {return element().end();};
00234         
00235   DeclException1(ExcDOFData, std::string, << "DOF data uncompatible: " << arg1);
00236 
00237   const unsigned int& efficiencyFlag() const {return effi_flag;};
00238   unsigned int& efficiencyFlag() {return effi_flag;}
00239 
00240   public:
00241 #ifdef __SERIALIZATION__ 
00242   template <class T> T *
00243   new_property(const element_t& e,
00244                const property_id_t<T>& pid) const {
00245     typedef RegularMesh<DIM,DOW> reg_mesh_t;
00246     const reg_mesh_t * mesh = dynamic_cast<const reg_mesh_t *>(msh);
00247     return mesh->template new_property<T,DIM>(e.geometry().index(), pid);
00248   }
00249   template <class T> T *
00250   get_property(const element_t& e,
00251                const property_id_t<T>& pid) const {
00252     typedef RegularMesh<DIM,DOW> reg_mesh_t;
00253     const reg_mesh_t * mesh = dynamic_cast<const reg_mesh_t *>(msh);
00254     return mesh->template get_property<T,DIM>(e.geometry().index(), pid);
00255   }
00256   template <class T> void
00257   free_property(const element_t& e,
00258                const property_id_t<T>& pid) const {
00259     typedef RegularMesh<DIM,DOW> reg_mesh_t;
00260     const reg_mesh_t * mesh = dynamic_cast<const reg_mesh_t *>(msh);
00261     mesh->template free_property<T,DIM>(e.geometry().index(), pid);
00262   }
00263 #endif // __SERIALIZATION__
00264 
00265   friend class Element<value_t,DIM,DOW,TDIM>;
00266 
00267 #ifdef MULTITHREAD
00268   void threadBuildElement(bool lazy, int n_thread = 1, int rank = 0);
00269   void threadBuildDof0(std::vector<std::vector<bool> >&, 
00270                        pthread_mutex_t&, int n_thread = 1, int rank = 0);
00271   void threadBuildDof1(std::vector<std::vector<bool> >&, 
00272                        pthread_mutex_t&, int n_thread = 1, int rank = 0);
00273 #endif
00274   };
00275 
00281 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double>
00282   class FEMFunction : public Vector<Number>
00283   {
00284   public:
00285   enum { dim = DIM, dow = DOW, tdim = TDIM };
00286 
00287   typedef value_type value_t;
00288   typedef Number number_t;
00289   typedef Element<value_t,DIM,DOW,TDIM> element_t;
00290   typedef FEMSpace<value_t,DIM,DOW,TDIM> fe_space_t;
00291   typedef FEMFunction<value_t,DIM,DOW,TDIM,number_t> fe_func_t;
00292 
00293   private:
00294   fe_space_t * sp; 
00295   public:
00296   FEMFunction(fe_space_t & = *((fe_space_t *)NULL)); 
00297   virtual ~FEMFunction(); 
00298   public:
00299   fe_space_t& femSpace(); 
00300   const fe_space_t& femSpace() const { return *sp; };
00301   void reinit(fe_space_t&, bool = false); 
00303   value_type value(const afepack::Point<DOW>&, const element_t&) const; 
00304   std::vector<value_type> gradient(const afepack::Point<DOW>&, const element_t&) const; 
00305   std::vector<value_type> value(const std::vector<afepack::Point<DOW> >&, const element_t&) const; 
00306   std::vector<std::vector<value_type> > gradient(const std::vector<afepack::Point<DOW> >&, const element_t&) const; 
00313   value_type value(const std::vector<value_type>&, const element_t&) const; 
00314   std::vector<value_type> gradient(const std::vector<std::vector<value_type> >&, const element_t&) const; 
00315   std::vector<value_type> value(const std::vector<std::vector<value_type> >&, const element_t&) const; 
00316   std::vector<std::vector<value_type> > gradient(const std::vector<std::vector<std::vector<value_type> > >&, const element_t&) const; 
00318   public:
00319   void loadData(const std::string& filename); 
00320   void writeData(const std::string& filename); 
00322   void writeEasyMeshData(const std::string& filename); 
00332   void writeTecplotData(const std::string& filename);
00344   void writeOpenDXData(const std::string& filename,
00345                        int flag = 0) const;
00346   };
00347 
00348 
00349 
00350 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double>
00351   class LocalFEMFunction : public Vector<Number>
00352   {
00353   private:
00354   Element<value_type,DIM,DOW,TDIM> * ele;
00355   public:
00356   LocalFEMFunction(Element<value_type,DIM,DOW,TDIM> & = *((Element<value_type,DIM,DOW,TDIM> *)NULL));
00357   virtual ~LocalFEMFunction() {};
00358   public:
00359   Element<value_type,DIM,DOW,TDIM>& element();
00360   void reinit(Element<value_type,DIM,DOW,TDIM> &);
00361   value_type value(const afepack::Point<DOW>&) const;
00362   std::vector<value_type> gradient(const afepack::Point<DOW>&) const;
00363   std::vector<value_type> value(const std::vector<afepack::Point<DOW> >&) const;
00364   std::vector<std::vector<value_type> > gradient(const std::vector<afepack::Point<DOW> >&) const;       
00365   };
00366 
00367 
00368 
00376 struct BoundaryConditionInfo
00377 {
00378 public:
00379   typedef int   Type;
00380   static Type   DIRICHLET;
00381   static Type   NEUMANN;
00382   static Type   ROBIN;
00383 };
00384 
00385 
00386 
00387 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double>
00388   class BoundaryCondition : public BoundaryConditionInfo
00389   {
00390   public:
00391   typedef typename FEMSpace<value_type,DIM,DOW,TDIM>::bmark_t bmark_t;
00392   typedef typename BoundaryConditionInfo::Type Type;
00393   private:
00394   Type                                  boundary_type;
00395   bmark_t                       boundary_mark;
00396   public:
00397   BoundaryCondition() {};
00398   BoundaryCondition(const Type& t, const bmark_t& m) :
00399   boundary_type(t), boundary_mark(m) {};
00400   BoundaryCondition(const BoundaryCondition<value_type,DIM>& b) :
00401   boundary_type(b.boundaryType()), boundary_mark(b.boundaryMark()) {};
00402   virtual ~BoundaryCondition() {};
00403   public:
00404   const Type& boundaryType() const {return boundary_type;};
00405   Type& boundaryType() {return boundary_type;};
00406   const bmark_t& boundaryMark() const {return boundary_mark;};
00407   bmark_t& boundaryMark() {return boundary_mark;};
00408   void reinit(const Type& t, const bmark_t& m) {
00409     boundary_type = t;
00410     boundary_mark = m;
00411   }
00412   BoundaryCondition<value_type,DIM,DOW,TDIM>& operator=(const BoundaryCondition<value_type,DIM,DOW,TDIM>& b) {
00413     boundary_type = b.boundaryType();
00414     boundary_mark = b.boundaryMark();
00415     return *this;
00416   };
00417   public:
00418   virtual Number value(const afepack::Point<DOW>&) const {return 0;};
00419   };
00420 
00421 
00422 
00427 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double>
00428   class BoundaryFunction : public BoundaryCondition<value_type,DIM,DOW,TDIM,Number>
00429   {
00430   public:
00431   typedef typename BoundaryCondition<value_type,DIM,DOW,TDIM,Number>::bmark_t bmark_t;
00432   typedef typename BoundaryCondition<value_type,DIM,DOW,TDIM,Number>::Type Type;
00433   private:
00434   bool is_newed;
00435   const Function<Number> * f;
00436   public:
00437   BoundaryFunction() : is_newed(false), f(NULL) {};
00438   BoundaryFunction(const Type& t,
00439                    const bmark_t& m,
00440                    const Function<Number>& fun) :
00441   BoundaryCondition<value_type,DIM>(t,m),
00442   is_newed(false),
00443   f(&fun)       {};
00444   BoundaryFunction(const Type& t,
00445                    const bmark_t& m, 
00446                    value_type (*fun)(const double *),
00447                    std::vector<value_type> (*grad)(const double *) = NULL) :
00448   BoundaryCondition<value_type,DIM,DOW,TDIM,Number>(t,m) 
00449   {
00450     is_newed = true;
00451     f = new FunctionFunction<Number>(fun, grad);
00452   };
00453   BoundaryFunction(const BoundaryFunction<value_type,DIM,DOW,TDIM,Number>& b) :
00454   BoundaryCondition<value_type,DIM,DOW,TDIM,Number>(b), 
00455   is_newed(false), 
00456   f(b.f) {};
00457   virtual ~BoundaryFunction() {
00458     if (is_newed && f != NULL) delete f;
00459   };
00460   public:
00461   void reinit(const Type& t,
00462               const bmark_t& m, 
00463               value_type (*fun)(const double *),
00464               std::vector<value_type> (*grad)(const double *) = NULL) {
00465     BoundaryCondition<value_type,DIM,DOW,TDIM,Number>::reinit(t, m);
00466     if (is_newed && f != NULL) delete f;
00467     is_newed = true;
00468     f = new FunctionFunction<Number>(fun, grad);
00469   }
00470   virtual Number value(const afepack::Point<DOW>& p) const {return f->value(p);};
00471   virtual std::vector<Number> gradient(const afepack::Point<DOW>& p) const {return f->gradient(p);};
00472   };
00473 
00474 
00475 
00482 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double>
00483   class BoundaryConditionAdmin : std::vector<const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *>
00484   {
00485   public:
00486   typedef typename FEMSpace<value_type,DIM,DOW,TDIM>::bmark_t bmark_t;
00487   private:
00488   std::vector<int>                                      index_map;
00489   const FEMSpace<value_type,DIM,DOW,TDIM> *             fem_space;
00490   public:
00491   BoundaryConditionAdmin(const FEMSpace<value_type,DIM,DOW,TDIM>& sp = 
00492                          *((const FEMSpace<value_type,DIM,DOW,TDIM> *)(NULL))) {fem_space = &sp;};
00493   ~BoundaryConditionAdmin() {};
00494   public:
00495   void reinit(const FEMSpace<value_type,DIM,DOW,TDIM>& sp) {fem_space = &sp; index_map.clear();};
00496   void setFemSpace(const FEMSpace<value_type,DIM,DOW,TDIM>& sp) {fem_space = &sp;};
00497   const FEMSpace<value_type,DIM,DOW,TDIM>& femSpace() const {return *fem_space;};
00498   void apply(SparseMatrix<double>& A, Vector<double>& u, Vector<double>& f, 
00499              bool preserve_symmetry = true); 
00501   void clearEntry(Vector<double>& f);
00502   void add(const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& b);
00503 
00504   bool isValid(const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& bc) const {
00505     return (&bc != NULL);
00506   };
00507   const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& find(const bmark_t& bm) const {
00508     if ((unsigned int)bm >= index_map.size() || index_map[bm] == -1)
00509       return *((const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *)NULL);
00510     else
00511       return *((*this)[index_map[bm]]);
00512   };
00513   };
00514 
00516 
00525 template <typename value_type, int DIM>
00526   struct ElementAdditionalData : public GeometryAdditionalData<DIM>
00527 {
00528  public:
00529   std::vector<std::vector<value_type> > basis_value;
00530   std::vector<std::vector<std::vector<value_type> > > basis_gradient;
00531 };
00532 
00533 #endif //_FEMSpace_h_
00534