Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
orthonormalpolynomialset.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2009-04-30
7 
8  Copyright (C) 2009 Universite Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __OrthonormalPolynomialSet_H
30 #define __OrthonormalPolynomialSet_H 1
31 namespace Feel
32 {
34 namespace detail
35 {
44 template<uint16_type Dim,
45  uint16_type Order,
46  uint16_type RealDim,
47  template<uint16_type> class PolySetType = Scalar,
48  typename T = double,
49  uint16_type TheTAG = 0,
50  template<uint16_type,uint16_type,uint16_type> class Convex = Simplex>
51 class OrthonormalPolynomialSet
52 {};
53 
54 template<uint16_type Dim,
55  uint16_type Order,
56  uint16_type RealDim,
57  template<uint16_type> class PolySetType,
58  typename T,
59  uint16_type TheTAG>
60 class OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType, T, TheTAG, Simplex>
61  :
62 public PolynomialSet<Dubiner<Dim, RealDim, Order, Normalized<true>, T, StorageUBlas>, PolySetType >
63 {
64  typedef PolynomialSet<Dubiner<Dim, RealDim, Order, Normalized<true>, T, StorageUBlas>, PolySetType > super;
65 public:
66 
67  static const uint16_type TAG = TheTAG;
68  static const uint16_type nDim = Dim;
69  static const uint16_type nOrder = Order;
70  static const uint16_type nRealDim = RealDim;
71  static const bool isTransformationEquivalent = true;
72  typedef OrthonormalPolynomialSet<Dim, Order,RealDim, PolySetType, T, TheTAG, Simplex> self_type;
73  typedef self_type component_basis_type;
74 
75  typedef typename super::polyset_type polyset_type;
76  static const bool is_tensor2 = polyset_type::is_tensor2;
77  static const bool is_vectorial = polyset_type::is_vectorial;
78  static const bool is_scalar = polyset_type::is_scalar;
79  static const bool is_continuous = false;
80  static const bool is_modal = true;
81  static const uint16_type nComponents = polyset_type::nComponents;
82  static const bool is_product = true;
83  static const bool isContinuous = false;
84  typedef Discontinuous continuity_type;
85 
86  typedef typename super::component_type component_type;
87 
88  typedef T value_type;
89  typedef Dubiner<Dim, RealDim, Order, Normalized<true>, T, StorageUBlas> basis_type;
90  typedef Simplex<Dim, Order, /*RealDim*/Dim> convex_type;
91  template<int O>
92  struct convex
93  {
94  typedef Simplex<Dim, O, /*RealDim*/Dim> type;
95  };
96  typedef Reference<convex_type, nDim, nOrder, nDim/*nRealDim*/, value_type> reference_convex_type;
97 
98  typedef typename super::polynomial_type polynomial_type;
99 
101  static const uint16_type nDofPerVertex = convex_type::nbPtsPerVertex;
103  static const uint16_type nDofPerEdge = convex_type::nbPtsPerEdge;
105  static const uint16_type nDofPerFace = convex_type::nbPtsPerFace;
106 
108  static const uint16_type nDofPerVolume = convex_type::nbPtsPerVolume;
109 
110  static const uint16_type nLocalDof = convex_type::numPoints;
111 
112  static const uint16_type nDof = nLocalDof;
113  static const uint16_type nNodes = nDof;
114  static const uint16_type nDofGrad = super::nDim*nDof;
115  static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
116  typedef typename matrix_node<value_type>::type points_type;
117 
118  OrthonormalPolynomialSet()
119  :
120  super( basis_type() )
121 
122  {
123  ublas::matrix<value_type> m( ublas::identity_matrix<value_type>( nComponents*convex_type::polyDims( nOrder ) ) );
124 
125  if ( is_tensor2 )
126  std::cout << "[orthonormalpolynomialset] m = " << m << "\n";
127 
128  if ( !( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
129  m ) < 1e-10 ) )
130  std::cout << "m1=" << m << "\n"
131  << "m2=" << polyset_type::toMatrix( polyset_type::toType( m ) ) << "\n"
132  << ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) - m ) << "\n";
133 
134  FEELPP_ASSERT( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
135  m ) < 1e-10 )( m ).warn ( "invalid transformation" );
136  this->setCoefficient( polyset_type::toType( m ), true );
137  }
138 
139  OrthonormalPolynomialSet<Dim, Order, RealDim, Scalar,T, TheTAG, Simplex > toScalar() const
140  {
141  return OrthonormalPolynomialSet<Dim, Order, RealDim, Scalar,T, TheTAG, Simplex >();
142  }
143 
147  std::string familyName() const
148  {
149  return "dubiner";
150  }
151 
152 
153  points_type points() const
154  {
155  return points_type();
156  }
157  points_type points( int f ) const
158  {
159  return points_type();
160  }
161 };
162 
163 template<uint16_type Dim,
164  uint16_type Order,
165  uint16_type RealDim,
166  template<uint16_type> class PolySetType,
167  typename T,
168  uint16_type TheTAG>
169 const uint16_type OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType,T, TheTAG, Simplex>::nLocalDof;
170 
171 
172 template<uint16_type Dim,
173  uint16_type Order,
174  uint16_type RealDim,
175  template<uint16_type> class PolySetType,
176  typename T,
177  uint16_type TheTAG>
178 class OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType, T, TheTAG, Hypercube>
179  :
180 public PolynomialSet<Legendre<Dim, RealDim, Order, Normalized<true>, T>, PolySetType >
181 {
182  typedef PolynomialSet<Legendre<Dim, RealDim, Order, Normalized<true>, T>, PolySetType > super;
183 public:
184 
185  static const uint16_type nDim = Dim;
186  static const uint16_type nOrder = Order;
187  static const uint16_type nRealDim = RealDim;
188  static const bool isTransformationEquivalent = true;
189 
190  typedef OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType, T, TheTAG, Hypercube> self_type;
191  typedef self_type component_basis_type;
192 
193  typedef typename super::polyset_type polyset_type;
194  static const bool is_tensor2 = polyset_type::is_tensor2;
195  static const bool is_vectorial = polyset_type::is_vectorial;
196  static const bool is_scalar = polyset_type::is_scalar;
197  static const bool is_continuous = false;
198  static const bool is_modal = true;
199  static const uint16_type nComponents = polyset_type::nComponents;
200  static const bool is_product = true;
201  static const bool isContinuous = false;
202  typedef Discontinuous continuity_type;
203 
204  typedef typename super::component_type component_type;
205  typedef T value_type;
206  typedef Legendre<Dim, RealDim, Order, Normalized<true>, T> basis_type;
207  typedef Hypercube<Dim, Order, /*RealDim*/Dim> convex_type;
208  typedef typename matrix_node<value_type>::type points_type;
209 
210  template<int O>
211  struct convex
212  {
213  typedef Hypercube<Dim, O, nDim/*RealDim*/> type;
214  };
215  typedef Reference<convex_type, nDim, nOrder, nDim/*nRealDim*/, value_type> reference_convex_type;
216 
217  typedef typename super::polynomial_type polynomial_type;
218 
220  static const uint16_type nDofPerVertex = convex_type::nbPtsPerVertex;
222  static const uint16_type nDofPerEdge = convex_type::nbPtsPerEdge;
224  static const uint16_type nDofPerFace = convex_type::nbPtsPerFace;
225 
227  static const uint16_type nDofPerVolume = convex_type::nbPtsPerVolume;
228 
229  static const uint16_type nLocalDof = convex_type::numPoints;
230 
231  static const uint16_type nDof = nLocalDof;
232  static const uint16_type nNodes = nDof;
233  static const uint16_type nDofGrad = super::nDim*nDof;
234  static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
235 
236  OrthonormalPolynomialSet()
237  :
238  super( basis_type() )
239 
240  {
241  ublas::matrix<value_type> m( ublas::identity_matrix<value_type>( nComponents*convex_type::polyDims( nOrder ) ) );
242 
243  if ( is_tensor2 )
244  std::cout << "[orthonormalpolynomialset] m = " << m << "\n";
245 
246  FEELPP_ASSERT( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
247  m ) < 1e-10 )( m ).warn ( "invalid transformation" );
248  this->setCoefficient( polyset_type::toType( m ), true );
249  }
250 
251  OrthonormalPolynomialSet<Dim, Order, RealDim,Scalar,T, TheTAG, Hypercube > toScalar() const
252  {
253  return OrthonormalPolynomialSet<Dim, Order, RealDim, Scalar,T, TheTAG, Hypercube >();
254  }
255  std::string familyName() const
256  {
257  return "legendre";
258  }
259  points_type points() const
260  {
261  return points_type();
262  }
263  points_type points( int f ) const
264  {
265  return points_type();
266  }
267 };
268 
269 template<uint16_type Dim,
270  uint16_type Order,
271  uint16_type RealDim,
272  template<uint16_type> class PolySetType,
273  typename T,
274  uint16_type TheTAG>
275 const uint16_type OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType,T, TheTAG, Hypercube>::nLocalDof;
276 } // detail
278 
279 template<uint16_type Order,
280  template<uint16_type Dim> class PolySetType = Scalar,
281  uint16_type TheTAG=0 >
282 class OrthonormalPolynomialSet
283 {
284 public:
285  template<uint16_type N,
286  uint16_type RealDim,
287  typename T = double,
288  typename Convex = Simplex<N> >
289  struct apply
290  {
291  typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
292  mpl::identity<detail::OrthonormalPolynomialSet<N,Order,RealDim,PolySetType,T,TheTAG,Simplex> >,
293  mpl::identity<detail::OrthonormalPolynomialSet<N,Order,RealDim,PolySetType,T,TheTAG,Hypercube> > >::type::type result_type;
294  typedef result_type type;
295  };
296 
297  template<uint16_type TheNewTAG>
298  struct ChangeTag
299  {
300  typedef OrthonormalPolynomialSet<Order,PolySetType,TheNewTAG> type;
301  };
302 
303  typedef OrthonormalPolynomialSet<Order,Scalar,TheTAG> component_basis_type;
304 
305  static const uint16_type nOrder = Order;
306  static const uint16_type TAG = TheTAG;
307 };
308 
309 } // Feel
310 #endif /* __OrthonormalPolynomialSet_H */
boost::tuple< mpl::size_t< MESH_POINTS >, typename MeshTraits< MeshType >::point_const_iterator, typename MeshTraits< MeshType >::point_const_iterator > points(MeshType const &mesh)
Definition: filters.hpp:1296

Generated on Sun Dec 22 2013 13:11:10 for Feel++ by doxygen 1.8.5