00001 #ifndef DUNE_ORTHONORMALSHAPEFUNCTIONS_HH
00002 #define DUNE_ORTHONORMALSHAPEFUNCTIONS_HH
00003
00004 #include <iostream>
00005 #include <vector>
00006 #include <map>
00007 #include <dune/common/fvector.hh>
00008 #include <dune/common/exceptions.hh>
00009 #include <dune/common/stdstreams.hh>
00010 #include <dune/common/misc.hh>
00011 #include <dune/grid/common/grid.hh>
00012
00013 #include <dune/disc/shapefunctions/shapefunctions.hh>
00014
00021 #define PMAX 8
00022
00023 namespace Dune
00024 {
00025
00031 template<typename C,typename T, int dim>
00032 class OrthonormalShapeFunction;
00033 template<typename C,typename T, int dim>
00034 class OrthonormalShapeFunctionSet;
00035 template<typename C,typename T, int dim,
00036 int maxOrder=PMAX>
00037 class OrthonormalShapeFunctionSetContainer;
00038
00044 template <int dim, int order>
00045 struct OrthonormalShapeFunctionSetSize
00046 {
00047 enum
00048 {
00049 maxSize =
00050 OrthonormalShapeFunctionSetSize<dim-1,order>::maxSize
00051 * (order+dim) / dim
00052 };
00053 };
00054
00060 template <int order>
00061 struct OrthonormalShapeFunctionSetSize<1,order>
00062 {
00063 enum { maxSize = order + 1 };
00064 };
00065
00071 template<typename C,typename T, int d>
00072 class OrthonormalShapeFunction : public ShapeFunction<C,T,d,1>
00073 {
00074 public:
00075 enum { dim=d };
00076 typedef C CoordType;
00077 typedef T ResultType;
00078 typedef FieldVector<C,dim> VectorType;
00079 typedef FieldVector<T,dim> GradientType;
00080 typedef OrthonormalShapeFunctionSet<C,T,dim> ShapeFunctionSet;
00081
00089 OrthonormalShapeFunction(int n, const ShapeFunctionSet & s):
00090 n_(n), s_(s) {}
00091
00093 virtual ResultType
00094 evaluateFunction (int, const FieldVector<CoordType,dim>& x) const
00095 {
00096 return s_.evaluateFunction(n_, x);
00097 }
00098
00100 virtual ResultType
00101 evaluateDerivative(int, int dir, const FieldVector<CoordType,dim>& x) const
00102 {
00103 GradientType g;
00104 s_.evaluateGradient(n_, x, g);
00105 return g[dir];
00106 }
00107
00109 virtual int localindex (int comp) const { return n_; };
00110
00112 virtual int codim () const { return 0; };
00113
00115 virtual int entity () const { return 0; };
00116
00118 virtual int entityindex () const { return 0; };
00119
00120 private:
00122 int n_;
00124 const ShapeFunctionSet & s_;
00125 };
00126
00142 template<typename C,typename T, int d>
00143 class OrthonormalShapeFunctionSet : public ShapeFunctionSet<C, T, d, 1>
00144 {
00145 public:
00146 enum { dim=d };
00147 typedef C CoordType;
00148 typedef T ResultType;
00149 typedef FieldVector<C,dim> VectorType;
00150 typedef FieldVector<T,dim> GradientType;
00151 typedef OrthonormalShapeFunction<C,T,dim> ShapeFunction;
00152
00156 OrthonormalShapeFunctionSet(const GeometryType & type, int order):
00157 type_(type), order_(order), n_(sz(order)), shapeFunctions(n_)
00158 {
00159 for (int i=0; i < n_; i++)
00160 {
00161 shapeFunctions[i] = new ShapeFunction(i,*this);
00162 }
00163 }
00164
00165 virtual ~OrthonormalShapeFunctionSet()
00166 {
00167 for (int i=0; i < n_; i++)
00168 {
00169 if (shapeFunctions[i]) delete shapeFunctions[i];
00170 }
00171 }
00172
00174 virtual int size() const
00175 {
00176 assert ((size_t)n_ == shapeFunctions.size());
00177 return n_;
00178 }
00179
00181 virtual int size (int entity, int codim) const
00182 {
00183 if (codim == 0) return n_;
00184 return 0;
00185 };
00186
00188 virtual const ShapeFunction& operator[](int i) const
00189 {
00190 assert (i < n_);
00191 assert ((size_t)n_ == shapeFunctions.size());
00192 return *shapeFunctions[i];
00193 }
00194
00196 virtual int order() const
00197 {
00198 return order_;
00199 }
00200
00202 virtual GeometryType type () const
00203 {
00204 return type_;
00205 };
00206
00208 virtual ResultType
00209 evaluateFunction (int i, const FieldVector<CoordType,dim>& x) const
00210 {
00211 if (type_.isTriangle())
00212 return eval_triangle_2d (i,x);
00213 if (type_.isQuadrilateral())
00214 return eval_quadrilateral_2d (i,x);
00215 if (type_.isTetrahedron())
00216 return eval_tetrahedron_3d (i,x);
00217 if (type_.isPyramid())
00218 return eval_pyramid_3d (i,x);
00219 if (type_.isPrism())
00220 return eval_prism_3d (i,x);
00221 if (type_.isHexahedron())
00222 return eval_hexahedron_3d (i,x);
00223 DUNE_THROW(Exception, "unhandled geometry type: " << type_);
00224 }
00225
00227 virtual void
00228 evaluateGradient(int i, const FieldVector<CoordType,dim>& x, GradientType & g) const
00229 {
00230 if (type_.isTriangle())
00231 return grad_triangle_2d (i,x,g);
00232 if (type_.isQuadrilateral())
00233 return grad_quadrilateral_2d (i,x,g);
00234 if (type_.isTetrahedron())
00235 return grad_tetrahedron_3d (i,x,g);
00236 if (type_.isPyramid())
00237 return grad_pyramid_3d (i,x,g);
00238 if (type_.isPrism())
00239 return grad_prism_3d (i,x,g);
00240 if (type_.isHexahedron())
00241 return grad_hexahedron_3d (i,x,g);
00242 DUNE_THROW(Exception, "unhandled geometry type: " << type_);
00243 }
00244
00245 private:
00246 const GeometryType type_;
00248 const int order_;
00250 const int n_;
00252 static int sz(int order)
00253 {
00254 int s;
00255 s = 1;
00256 for (int i=1; i<=dim; i++) {
00257 s = s * (order + i) / i;
00258 }
00259 return s;
00260 }
00262 std::vector< ShapeFunction* > shapeFunctions;
00263
00264 double eval_triangle_2d ( int i, const VectorType & xi ) const;
00265 double eval_quadrilateral_2d ( int i, const VectorType & xi ) const;
00266 double eval_tetrahedron_3d ( int i, const VectorType & xi ) const;
00267 double eval_pyramid_3d ( int i, const VectorType & xi ) const;
00268 double eval_prism_3d ( int i, const VectorType & xi ) const;
00269 double eval_hexahedron_3d ( int i, const VectorType & xi ) const;
00270 void grad_triangle_2d ( int i, const VectorType & xi, GradientType & grad ) const;
00271 void grad_quadrilateral_2d ( int i, const VectorType & xi, GradientType & grad ) const;
00272 void grad_tetrahedron_3d ( int i, const VectorType & xi, GradientType & grad ) const;
00273 void grad_pyramid_3d ( int i, const VectorType & xi, GradientType & grad ) const;
00274 void grad_prism_3d ( int i, const VectorType & xi, GradientType & grad ) const;
00275 void grad_hexahedron_3d ( int i, const VectorType & xi, GradientType & grad ) const;
00276
00277 };
00278
00293 template<typename C,typename T, int d, int maxOrder>
00294 class OrthonormalShapeFunctionSetContainer :
00295 public ShapeFunctionSetContainer<C,T,d,1,OrthonormalShapeFunctionSetSize<d,maxOrder>::maxSize >
00296 {
00297 public:
00298 enum { dim = d };
00299 typedef C CoordType;
00300 typedef T ResultType;
00301 typedef FieldVector<C,dim> VectorType;
00302 typedef FieldVector<T,dim> GradientType;
00303 typedef OrthonormalShapeFunctionSet<C,T,dim> ShapeFunctionSet;
00304
00308 OrthonormalShapeFunctionSetContainer()
00309 {
00310 assert (maxOrder <= PMAX);
00311 }
00312
00313 virtual ~OrthonormalShapeFunctionSetContainer()
00314 {
00315 typename
00316 std::map< std::pair<GeometryType,int>, ShapeFunctionSet* >::iterator
00317 it = shapeFunctionSets.begin();
00318 for (; it != shapeFunctionSets.end(); ++it)
00319 {
00320 delete it->second;
00321 }
00322 }
00323
00325 virtual const ShapeFunctionSet& operator() (GeometryType type, int order) const
00326 {
00327 assert(order <= maxOrder);
00328 std::pair<GeometryType, int> key(type, order);
00329 typename
00330 std::map< std::pair<GeometryType,int>, ShapeFunctionSet* >::iterator
00331 it = shapeFunctionSets.find(key);
00332 if (it == shapeFunctionSets.end())
00333 {
00334 ShapeFunctionSet *sf = new ShapeFunctionSet(type, order);
00335 shapeFunctionSets[key] = sf;
00336 return *sf;
00337 }
00338 return *(it->second);
00339 }
00340
00341 private:
00343 mutable std::map< std::pair<GeometryType,int>, ShapeFunctionSet* > shapeFunctionSets;
00344 };
00345
00348 }
00349
00350 #include "orthonormalshapefunctionsset.cc"
00351
00352 #undef PMAX
00353
00354 #endif // DUNE_ORTHONORMALSHAPEFUNCTIONS_HH