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
expansions.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: 2005-07-29
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 __expansions_H
30 #define __expansions_H 1
31 
32 #include <feel/feelalg/glas.hpp>
34 #include <feel/feelpoly/jacobi.hpp>
35 
36 namespace Feel
37 {
38 enum Shape { LINE = 1, TRIANGLE = 2, TETRAHEDRON = 3 };
39 template<Shape sh>
40 struct DimFromShape
41 {
42  static const uint16_type value = mpl::if_<mpl::equal_to<mpl::int_<sh>, mpl::int_<LINE> >,
43  mpl::int_<1>,
44  typename mpl::if_<mpl::equal_to<mpl::int_<sh>, mpl::int_<TRIANGLE> >,
45  mpl::int_<2>,
46  mpl::int_<3> >::type>::type::value;
47 };
48 
50 
53 namespace details
54 {
55 
59 template<Shape sh, typename T = double>
60 struct xi
61 {
62 };
63 
64 template<typename T>
65 struct xi<TRIANGLE, T>
66 {
67  typedef T value_type;
68  typedef typename node<value_type>::type node_type;
69 
70  xi()
71  :
72  M_xi( 2 )
73  {}
74 
75  xi( node_type const& eta )
76  :
77  M_xi( 2 )
78  {
79  M_xi[0] = 0.5*( 1.0 + eta[0] )*( 1.0 - eta[1] ) - 1.0;
80  M_xi[1] = eta[1];
81  }
82 
83  node_type const& operator()( node_type const& eta )
84  {
85  M_xi[0] = 0.5*( 1.0 + eta[0] )*( 1.0 - eta[1] ) - 1.0;
86  M_xi[1] = eta[1];
87  return M_xi;
88  }
89  node_type const& operator()() const
90  {
91  return M_xi;
92  }
93  node_type M_xi;
94 };
95 
96 template<typename T>
97 struct xi<TETRAHEDRON, T>
98 {
99  typedef T value_type;
100  typedef typename node<value_type>::type node_type;
101 
102  xi()
103  :
104  M_xi( 3 )
105  {}
106 
107  xi( node_type const& eta )
108  :
109  M_xi( 3 )
110  {
111  M_xi[0] = 0.25*( 1.0 + eta[0] )*( 1.0 - eta[1] )*( 1.0 - eta[2] ) - 1.0;
112  M_xi[1] = 0.5*( 1.0 + eta[1] )*( 1.0 - eta[2] ) - 1.0;
113  M_xi[2] = eta[2];
114  }
115 
116  node_type const& operator()( node_type const& eta )
117  {
118  M_xi[0] = 0.25*( 1.0 + eta[0] )*( 1.0 - eta[1] )*( 1.0 - eta[2] ) - 1.0;
119  M_xi[1] = 0.5*( 1.0 + eta[1] )*( 1.0 - eta[2] ) - 1.0;
120  M_xi[2] = eta[2];
121  return M_xi;
122  }
123  node_type const& operator()() const
124  {
125  return M_xi;
126  }
127  node_type M_xi;
128 };
129 
130 
134 template<Shape sh, typename T = double>
135 struct eta
136 {
137 };
138 template<Shape sh, typename T = double>
139 struct etas
140 {
141 };
142 
143 template<typename T>
144 struct eta<TRIANGLE, T>
145 {
146  typedef T value_type;
147  typedef typename node<value_type>::type node_type;
148 
149  eta()
150  :
151  M_eta( 2 )
152  {}
153 
154  eta( node_type const& xi )
155  :
156  M_eta( 2 )
157  {
158  if ( xi[1] == 1.0 )
159  M_eta[0] = -1.0;
160 
161  else
162  M_eta[0] = 2.0 * ( 1.0 + xi[0] ) / ( 1.0 - xi[1] ) - 1.0;
163 
164  M_eta[1] = xi[1];
165  }
166 
167  node_type const& operator()( node_type const& xi )
168  {
169  if ( xi[1] == 1.0 )
170  M_eta[0] = -1.0;
171 
172  else
173  M_eta[0] = 2.0 * ( 1.0 + xi[0] ) / ( 1.0 - xi[1] ) - 1.0;
174 
175  M_eta[1] = xi[1];
176  return M_eta;
177  }
178  node_type const& operator()() const
179  {
180  return M_eta;
181  }
182  node_type M_eta;
183 };
184 template<typename T>
185 struct etas<TRIANGLE, T>
186 {
187  typedef T value_type;
188  typedef typename matrix_node<value_type>::type matrix_node_type;
189 
190  etas()
191  :
192  M_eta()
193  {}
194 
195  etas( matrix_node_type const& xi )
196  :
197  M_eta( xi.size1(), xi.size2() )
198  {
199  for ( size_type i = 0; i < xi.size2(); ++i )
200  {
201  if ( xi( 1, i ) == 1.0 )
202  M_eta( 0, i ) = -1.0;
203 
204  else
205  M_eta( 0, i ) = 2.0 * ( 1.0 + xi( 0, i ) ) / ( 1.0 - xi( 1, i ) ) - 1.0;
206 
207  M_eta( 1, i ) = xi( 1, i );
208  }
209 
210  }
211 
212  matrix_node_type const& operator()( matrix_node_type const& xi )
213  {
214  M_eta.resize( xi.size1(), xi.size2() );
215 
216  for ( size_type i = 0; i < xi.size2(); ++i )
217  {
218  if ( xi( 1, i ) == 1.0 )
219  M_eta( 0, i ) = -1.0;
220 
221  else
222  M_eta( 0, i ) = 2.0 * ( 1.0 + xi( 0, i ) ) / ( 1.0 - xi( 1, i ) ) - 1.0;
223 
224  M_eta( 1, i ) = xi( 1, i );
225  }
226 
227  return M_eta;
228  }
229  matrix_node_type const& operator()() const
230  {
231  return M_eta;
232  }
233  matrix_node_type M_eta;
234 };
235 
236 template<typename T>
237 struct etas<TETRAHEDRON, T>
238 {
239  typedef T value_type;
240  typedef typename matrix_node<value_type>::type matrix_node_type;
241 
242  etas()
243  :
244  M_eta()
245  {}
246 
247  etas( matrix_node_type const& xi )
248  :
249  M_eta( xi.size1(), xi.size2() )
250  {
251  for ( size_type i = 0; i < xi.size2(); ++i )
252  {
253  if ( xi( 1, i ) + xi( 2, i ) == 0. )
254  M_eta( 0, i ) = 1.;
255 
256  else
257  M_eta( 0, i ) = -2. * ( 1. + xi( 0, i ) ) / ( xi( 1, i ) + xi( 2, i ) ) - 1.;
258 
259  if ( xi( 2, i ) == 1. )
260  M_eta( 1, i ) = -1.;
261 
262  else
263  M_eta( 1, i ) = 2. * ( 1. + xi( 1, i ) ) / ( 1. - xi( 2, i ) ) - 1.;
264 
265  M_eta( 2, i ) = xi( 2, i );
266  }
267  }
268 
269  matrix_node_type const& operator()( matrix_node_type const& xi )
270  {
271  M_eta.resize( xi.size1(), xi.size2() );
272 
273  for ( size_type i = 0; i < xi.size2(); ++i )
274  {
275  if ( xi( 1, i ) + xi( 2, i ) == 0. )
276  M_eta( 0, i ) = 1.;
277 
278  else
279  M_eta( 0, i ) = -2. * ( 1. + xi( 0, i ) ) / ( xi( 1, i ) + xi( 2, i ) ) - 1.;
280 
281  if ( xi( 2, i ) == 1. )
282  M_eta( 1, i ) = -1.;
283 
284  else
285  M_eta( 1, i ) = 2. * ( 1. + xi( 1, i ) ) / ( 1. - xi( 2, i ) ) - 1.;
286 
287  M_eta( 2, i ) = xi( 2, i );
288  }
289 
290  return M_eta;
291  }
292  matrix_node_type const& operator()() const
293  {
294  return M_eta;
295  }
296  matrix_node_type M_eta;
297 };
298 
311 template<typename T>
312 class psitilde
313 {
314 public:
315  typedef T value_type;
316  typedef Feel::dyna::Jacobi<T> P;
317  psitilde()
318  :
319  p( 0, 0.0, 0.0 ),
320  M_b( 0.0 )
321  {}
322  psitilde( int i )
323  :
324  p( i, 0.0, 0.0 ),
325  M_b( i )
326  {
327  }
328  psitilde( int i, int j )
329  :
330  p( j, 2*i+1, 0.0 ),
331  M_b( i )
332  {
333  }
334  psitilde( int i, int j, int k )
335  :
336  p( k, 2*( i+j+1 ), 0.0 ),
337  M_b( i+j )
338  {
339  }
345  value_type a( value_type const& x ) const
346  {
347  return p( x );
348  }
354  value_type b( value_type const& x ) const
355  {
356  return math::pow( 0.5 *( 1.0-x ), M_b )*p( x );
357  }
363  value_type c( value_type const& x ) const
364  {
365  return math::pow( 0.5 *( 1.0-x ), M_b )*p( x );
366  }
367  P p;
368  value_type M_b;
369 };
370 
377 template<uint16_type N, typename T = double>
378 struct scalings
379 {
380  typedef T value_type;
381  typedef ublas::matrix<value_type> matrix_type;
382  typedef typename matrix_node<value_type>::type matrix_node_type;
383 
389  scalings( ublas::vector<value_type> const& pts )
390  :
391  M_s( N+1, pts.size() )
392  {
393 #if 0
394  ublas::vector<value_type> one( ublas::scalar_vector<value_type>( pts.size(), 1.0 ) );
395  ublas::row( M_s, 0 ) = one;
396 #else
397  ublas::row( M_s, 0 ) = ublas::scalar_vector<value_type>( pts.size(), 1.0 );
398 #endif
399 
400  if ( N > 0 )
401  {
402  ublas::row( M_s, 1 ) = 0.5 * ( ublas::row( M_s, 0 ) - pts );
403 
404  for ( uint16_type k = 2; k < N+1; ++k )
405  {
406  ublas::row( M_s, k ) = ublas::element_prod( ublas::row( M_s, k-1 ),
407  ublas::row( M_s, 1 ) );
408  }
409  }
410  }
411  matrix_type const& operator()() const
412  {
413  return M_s;
414  }
415 
416  matrix_type M_s;
417 };
418 } // details
420 
421 
422 }
423 
424 #endif /* __expansions_H */

Generated on Fri Oct 25 2013 14:24:08 for Feel++ by doxygen 1.8.4