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
feelpoly/policy.hpp
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: 2005-12-03
7 
8  Copyright (C) 2005,2006 EPFL
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 __policy_H
30 #define __policy_H 1
31 
32 
33 
34 #include <boost/mpl/vector.hpp>
35 #include <boost/numeric/ublas/matrix.hpp>
36 #include <boost/numeric/ublas/matrix_proxy.hpp>
37 #include <boost/numeric/ublas/vector.hpp>
38 
39 #include <boost/mpl/if.hpp>
40 #include <boost/mpl/find.hpp>
41 
42 #include <feel/feelcore/feel.hpp>
43 #include <feel/feelalg/glas.hpp>
44 namespace Feel
45 {
46 namespace ublas = boost::numeric::ublas;
47 
48 namespace fem
49 {
50 enum transformation_type { LINEAR, BILINEAR, NONLINEAR };
51 }
52 
58 template<uint16_type Dim>
59 struct Scalar
60 {
61  static const uint16_type rank = 0;
62  static const uint16_type nDim = Dim;
63 
64  static const bool is_scalar = true;
65  static const bool is_vectorial = false;
66  static const bool is_tensor2 = false;
67  static const bool is_tensor3 = false;
68 
69  static const uint16_type nComponents = 1;
70  static const uint16_type nComponents1 = 1;
71  static const uint16_type nComponents2 = 1;
72  static const uint16_type nComponents3 = 1;
73  static const uint16_type nComponentsLast = 1;
74 
75  template<typename T>
76  static
77  inline ublas::matrix<T> const&
78  toMatrix( ublas::matrix<T> const& __c )
79  {
80  return __c;
81  }
82  template<typename T>
83  static
84  inline ublas::matrix<T> const&
85  toType( ublas::matrix<T> const& __c )
86  {
87  return __c;
88  }
89 };
90 
91 
97 template<uint16_type Dim>
98 struct Vectorial
99 {
100  static const uint16_type rank = 1;
101  static const uint16_type nDim = Dim;
102 
103  static const bool is_scalar = false;
104  static const bool is_vectorial = true;
105  static const bool is_tensor2 = false;
106  static const bool is_tensor3 = false;
107 
108  static const uint16_type nComponents = nDim;
109  static const uint16_type nComponents1 = nDim;
110  static const uint16_type nComponents2 = 1;
111  static const uint16_type nComponents3 = 1;
112  static const uint16_type nComponentsLast = nComponents1;
113  template<typename T>
114  static
115  ublas::matrix<T>
116  toMatrix( ublas::matrix<T> const& __c )
117  {
118  typedef T value_type;
119  // reshape the coefficients in the vectorial case
120  const size_type nRows = __c.size1()/nComponents;
121  const size_type nCols = __c.size2();
122  ublas::matrix<T> __c_reshaped( ublas::zero_matrix<value_type>( nRows, nComponents*nCols ) );
123 
124  for ( int c = 0; c < nComponents; ++c )
125  {
126  ublas::project( __c_reshaped,
127  ublas::range( 0, nRows ),
128  ublas::range( c*nCols, ( c+1 )*nCols ) ) = ublas::project( __c,
129  ublas::slice( c, nComponents, nRows ),
130  ublas::slice( 0, 1, nCols ) );
131  }
132 
133  return __c_reshaped;
134  }
135 
136  template<typename AE>
137  static
138  ublas::matrix<typename ublas::matrix_expression<AE>::value_type>
139  toType( ublas::matrix_expression<AE> const& __c )
140  {
141  typedef typename ublas::matrix_expression<AE>::value_type value_type;
142 
143  // reshape the coefficients in the vectorial case
144  const size_type nRows = __c().size1()*nComponents;
145  const size_type nCols = __c().size2()/nComponents;
146  ublas::matrix<value_type> __c_reshaped( ublas::zero_matrix<value_type>( nRows, nCols ) );
147 
148  for ( int c = 0; c < nComponents; ++c )
149  {
150  ublas::project( __c_reshaped,
151  ublas::slice( c, nComponents, nRows/nComponents ),
152  ublas::slice( 0, 1, nCols ) ) = ublas::project( __c(),
153  ublas::range( 0, nRows/nComponents ),
154  ublas::range( c*nCols, ( c+1 )*nCols ) );
155  }
156 
157  return __c_reshaped;
158  }
159 
160  template<typename T>
161  static
162  ublas::matrix<T>
163  toType( ublas::matrix<T> const& __c )
164  {
165  typedef T value_type;
166 
167  // reshape the coefficients in the vectorial case
168  const size_type nRows = __c.size1()*nComponents;
169  const size_type nCols = __c.size2()/nComponents;
170  ublas::matrix<value_type> __c_reshaped( ublas::zero_matrix<value_type>( nRows, nCols ) );
171 
172  for ( int c = 0; c < nComponents; ++c )
173  {
174  ublas::project( __c_reshaped,
175  ublas::slice( c, nComponents, nRows/nComponents ),
176  ublas::slice( 0, 1, nCols ) ) = ublas::project( __c,
177  ublas::range( 0, nRows/nComponents ),
178  ublas::range( c*nCols, ( c+1 )*nCols ) );
179  }
180 
181  return __c_reshaped;
182  }
183 };
184 
189 namespace detail
190 {
191 template<uint16_type N, uint16_type M>
192 struct Field
193 {
194  static const uint16_type rank = ( M > 1 );
195  static const uint16_type nDim = N;
196  static const uint16_type nVariables = N;
197 
198  static const bool is_scalar = ( M==1 );
199  static const bool is_vectorial = ( N==M );
200  static const bool is_tensor2 = false;
201  static const bool is_tensor3 = false;
202 
203  static const uint16_type nComponents = M;
204  static const uint16_type nComponents1 = M;
205  static const uint16_type nComponents2 = 1;
206  static const uint16_type nComponents3 = 1;
207  static const uint16_type nComponentsLast = 1;
208 
209  template<typename T>
210  static
211  inline ublas::matrix<T> const&
212  toMatrix( ublas::matrix<T> const& __c )
213  {
214  return __c;
215  }
216  template<typename T>
217  static
218  inline ublas::matrix<T> const&
219  toType( ublas::matrix<T> const& __c )
220  {
221  return __c;
222  }
223 };
224 }
225 
226 template<uint16_type M>
227 struct Field
228 {
229  template <uint16_type Nvar>
230  struct apply
231  {
232  typedef detail::Field<Nvar,M> type;
233  };
234 };
235 
241 template<uint16_type Dim>
242 struct Tensor2
243 {
244  static const uint16_type rank = 2;
245  static const uint16_type nDim = Dim;
246 
247  static const bool is_scalar = false;
248  static const bool is_vectorial = false;
249  static const bool is_tensor2 = true;
250  static const bool is_tensor3 = false;
251 
252  static const uint16_type nComponents = nDim*nDim;
253  static const uint16_type nComponents1 = nDim;
254  static const uint16_type nComponents2 = nDim;
255  static const uint16_type nComponents3 = 1;
256  static const uint16_type nComponentsLast = nComponents2;
257 
258  template<typename T>
259  static
260  ublas::matrix<T>
261  toMatrix( ublas::matrix<T> const& __c )
262  {
263  typedef T value_type;
264  // reshape the coefficients in the vectorial case
265  const size_type nRows = __c.size1();
266  const size_type nRows1= __c.size2()*nComponents;
267  const size_type nCols = __c.size2();
268  ublas::matrix<T> __c_reshaped( nRows/nComponents, nCols*nComponents );
269 
270  //__c_reshaped = ublas::scalar_matrix<value_type>( nRows/nComponents, nCols*nComponents, -1 );
271  for ( int c1 = 0; c1 < nComponents; ++c1 )
272  {
273  uint16_type i1 = nRows1*c1;
274 
275  for ( int c2 = 0; c2 < nComponents; ++c2 )
276  {
277  ublas::project( __c_reshaped,
278  ublas::range( i1/nComponents, ( i1+nRows1 )/nComponents ),
279  ublas::range( c2*nCols, ( c2+1 )*nCols ) ) = ublas::project( __c,
280  ublas::slice( i1+c2, nComponents, nRows1/nComponents ),
281  ublas::slice( 0, 1, nCols ) );
282  }
283  }
284 
285  return __c_reshaped;
286  }
287 
288  template<typename AE>
289  static
290  ublas::matrix<typename ublas::matrix_expression<AE>::value_type>
291  toType( ublas::matrix_expression<AE> const& __c )
292  {
293 
294  }
295 
296  template<typename T>
297  static
298  ublas::matrix<T>
299  toType( ublas::matrix<T> const& __c )
300  {
301  typedef T value_type;
302  // reshape the coefficients in the vectorial case
303  const size_type nRows = __c.size1()*nComponents;
304  const size_type nRows1= __c.size2();
305  const size_type nCols = __c.size2()/nComponents;
306  ublas::matrix<T> __c_reshaped( nRows, nCols );
307 
308  //__c_reshaped = ublas::scalar_matrix<value_type>( nRows, nCols, -1 );
309  for ( int c1 = 0; c1 < nComponents; ++c1 )
310  {
311  uint16_type i1 = nRows1*c1;
312 
313  for ( int c2 = 0; c2 < nComponents; ++c2 )
314  {
315  ublas::project( __c_reshaped,
316  ublas::slice( i1+c2, nComponents, nRows1/nComponents ),
317  ublas::slice( 0, 1, nCols ) ) = ublas::project( __c,
318  ublas::range( i1/nComponents, ( i1+nRows1 )/nComponents ),
319  ublas::range( c2*nCols, ( c2+1 )*nCols ) );
320  }
321  }
322 
323  return __c_reshaped;
324  }
325 };
326 
332 template<uint16_type Dim>
333 struct Tensor3
334 {
335  static const uint16_type rank = 3;
336  static const uint16_type nDim = Dim;
337 
338  static const bool is_scalar = false;
339  static const bool is_vectorial = false;
340  static const bool is_tensor2 = false;
341  static const bool is_tensor3 = true;
342 
343  static const uint16_type nComponents = nDim*nDim*nDim;
344  static const uint16_type nComponents1 = nDim;
345  static const uint16_type nComponents2 = nDim;
346  static const uint16_type nComponents3 = nDim;
347  static const uint16_type nComponentsLast = nComponents3;
348 };
349 
350 template<int Dim>
351 struct ListReturnTypes
352 {
353  typedef mpl::vector<Scalar<Dim>, Vectorial<Dim>, Tensor2<Dim>, Tensor3<Dim> > return_types;
354 };
355 template<typename T1, typename T2>
356 struct ReturnSelect
357 {
358  typedef typename mpl::if_<boost::is_same<T1, T2>,
359  mpl::identity<T1>,
360  typename mpl::if_<mpl::greater<mpl::int_<T1::rank>, mpl::int_<T2::rank> >,
361  mpl::identity<T1>,
362  mpl::identity<T2> >::type>::type::type type;
363 };
364 
365 enum EnumIndex
366 {
367  GLOBAL_COMPONENT,
368  COMPONENT_IN_COMPONENT,
369  GLOBAL_FUNCTION_INDEX,
370  PER_COMPONENT_FUNCTION_INDEX,
371  COMPONENT_IN_COMPONENT_FUNCTION_INDEX,
372  FUNCTION_INDEX
373 };
374 
375 const mpl::int_<GLOBAL_COMPONENT> INDEX_GLOBAL_COMPONENT = mpl::int_<GLOBAL_COMPONENT>();
376 const mpl::int_<COMPONENT_IN_COMPONENT> INDEX_COMPONENT_IN_COMPONENT = mpl::int_<COMPONENT_IN_COMPONENT>() ;
377 const mpl::int_<GLOBAL_FUNCTION_INDEX> INDEX_GLOBAL_FUNCTION_INDEX = mpl::int_<GLOBAL_FUNCTION_INDEX>();
378 const mpl::int_<PER_COMPONENT_FUNCTION_INDEX> INDEX_PER_COMPONENT_FUNCTION_INDEX = mpl::int_<PER_COMPONENT_FUNCTION_INDEX>();
379 const mpl::int_<COMPONENT_IN_COMPONENT_FUNCTION_INDEX> INDEX_COMPONENT_IN_COMPONENT_FUNCTION_INDEX = mpl::int_<COMPONENT_IN_COMPONENT_FUNCTION_INDEX>();
380 const mpl::int_<FUNCTION_INDEX> INDEX_FUNCTION_INDEX = mpl::int_<FUNCTION_INDEX>();
381 
389 template<typename T>
390 struct Component
391 {
392  static const uint16_type nDim = T::nDim;
393  typedef mpl::vector<Scalar<nDim>, Vectorial<nDim>, Tensor2<nDim> > types;
394 #if 0
395  typedef typename mpl::find<types, T>::type iter;
396  typedef typename mpl::if_<boost::is_same<T, Scalar<nDim> >,
397  mpl::identity<Scalar<nDim> >,
398  mpl::identity<typename mpl::deref<typename mpl::prior<iter>::type>::type> >::type::type type;
399 #else
400  typedef typename mpl::if_<boost::is_same<T, Scalar<nDim> >,
401  mpl::identity<Scalar<nDim> >,
402  typename mpl::if_<boost::is_same<T, Vectorial<nDim> >,
403  mpl::identity<Vectorial<nDim> >,
404  typename mpl::if_<boost::is_same<T, Tensor2<nDim> >,
405  mpl::identity<Tensor2<nDim> > >::type>::type>::type::type type;
406 #endif
407 };
408 
416 template<typename T>
417 struct RankUp
418 {
419  static const uint16_type nDim = T::nDim;
420  typedef mpl::vector<Scalar<nDim>, Vectorial<nDim>, Tensor2<nDim>, Tensor3<nDim> > types;
421  typedef typename mpl::find<types, T>::type iter;
422  typedef typename mpl::deref<typename mpl::next<iter>::type>::type type;
423 };
424 
425 template<typename T>
426 struct RankUp2
427 {
428  static const uint16_type nDim = T::nDim;
429  typedef mpl::vector<Scalar<nDim>, Vectorial<nDim>, Tensor2<nDim>, Tensor3<nDim> > types;
430  typedef typename mpl::find<types, T>::type iter;
431  typedef typename mpl::deref<typename mpl::next<typename mpl::next<iter>::type>::type>::type type;
432 };
433 
434 template<typename T>
435 struct RankSame
436 {
437  static const uint16_type nDim = T::nDim;
438  typedef T type;
439 };
440 
441 template<typename T>
442 struct RankDown
443 {
444  static const uint16_type nDim = T::nDim;
445  typedef mpl::vector<Scalar<nDim>, Vectorial<nDim>, Tensor2<nDim>, Tensor3<nDim> > types;
446  typedef typename mpl::find<types, T>::type iter;
447  typedef typename mpl::deref<typename mpl::prior<iter>::type>::type _type;
448  typedef typename mpl::if_<boost::is_same<_type,mpl::void_>,
449  mpl::identity<Scalar<nDim> >,
450  mpl::identity<_type> >::type::type type;
451 };
457 template<bool normalized>
459 {
460  static const bool is_normalized = normalized;
461 };
462 
468 template<typename T>
470 {
471  typedef T value_type;
472  typedef ublas::vector<value_type> vector_type;
473  typedef ublas::matrix<value_type> matrix_type;
474  typedef ublas::vector<matrix_type> vector_matrix_type;
475  typedef ublas::vector<vector_matrix_type> vector_vector_matrix_type;
476  typedef typename matrix_node<value_type>::type matrix_node_type;
477  typedef typename matrix_node<value_type>::type points_type;
478  typedef typename node<value_type>::type node_type;
479 
480 };
481 
482 } // Feel
483 
484 #endif /* __policy_H */
Definition: feelpoly/policy.hpp:458
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
Definition: feelpoly/policy.hpp:98
Definition: feelpoly/policy.hpp:390
Definition: glas.hpp:157
Definition: feelpoly/policy.hpp:417
Definition: glas.hpp:125
Definition: feelpoly/policy.hpp:469
Definition: feelpoly/policy.hpp:59
Definition: feelpoly/policy.hpp:333
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Definition: feelpoly/policy.hpp:242

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