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
inv.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: 2012-04-26
7 
8  Copyright (C) 2012 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 __FEELPP_VF_Inv_H
30 #define __FEELPP_VF_Inv_H 1
31 
32 namespace Feel
33 {
34 namespace vf
35 {
37 
44 template<typename ExprT>
45 class Inv
46 {
47 public:
48 
49  static const size_type context = ExprT::context;
50  static const bool is_terminal = false;
51 
52  static const uint16_type imorder = ExprT::imorder;
53  static const bool imIsPoly = ExprT::imIsPoly;
54 
55  template<typename Func>
56  struct HasTestFunction
57  {
58  static const bool result = ExprT::template HasTestFunction<Func>::result;
59  };
60 
61  template<typename Func>
62  struct HasTrialFunction
63  {
64  static const bool result = ExprT::template HasTrialFunction<Func>::result;
65  };
66 
67 
71 
72  typedef ExprT expression_type;
73  typedef typename expression_type::value_type value_type;
74  typedef Inv<ExprT> this_type;
75 
76 
78 
82 
83  explicit Inv( expression_type const & __expr )
84  :
85  M_expr( __expr )
86  {}
87  Inv( Inv const & te )
88  :
89  M_expr( te.M_expr )
90  {}
91  ~Inv()
92  {}
93 
95 
99 
100 
102 
106 
107 
109 
113 
114 
116 
120 
121  expression_type const& expression() const
122  {
123  return M_expr;
124  }
125 
127 
128  //template<typename Geo_t, typename Basis_i_t = fusion::map<fusion::pair<vf::detail::gmc<0>,boost::shared_ptrvf::detail::gmc<0> > > >, typename Basis_j_t = Basis_i_t>
129  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
130  struct tensor
131  {
132  typedef typename expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> tensor_expr_type;
133  typedef typename tensor_expr_type::value_type value_type;
134 
135  typedef typename tensor_expr_type::shape expr_shape;
136  BOOST_MPL_ASSERT_MSG( ( boost::is_same<mpl::int_<expr_shape::M>,mpl::int_<expr_shape::N> >::value ), INVALID_TENSOR_SHOULD_BE_RANK_2_OR_0, ( mpl::int_<expr_shape::M>, mpl::int_<expr_shape::N> ) );
137  typedef Shape<expr_shape::nDim,Tensor2,false,false> shape;
138 
139 
140  typedef Eigen::Matrix<value_type,shape::M,shape::N> matrix_type;
141  typedef std::vector<matrix_type> inv_matrix_type;
142 
143 
144  template <class Args> struct sig
145  {
146  typedef value_type type;
147  };
148 
149  struct is_zero
150  {
151  static const bool value = tensor_expr_type::is_zero::value;
152  };
153 
154  tensor( this_type const& expr,
155  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
156  :
157  M_tensor_expr( expr.expression(), geom, fev, feu ),
158  M_inv( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
159  {
160  }
161 
162  tensor( this_type const& expr,
163  Geo_t const& geom, Basis_i_t const& fev )
164  :
165  M_tensor_expr( expr.expression(), geom, fev ),
166  M_inv( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
167  {
168  }
169 
170  tensor( this_type const& expr, Geo_t const& geom )
171  :
172  M_tensor_expr( expr.expression(), geom ),
173  M_inv( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
174  {
175  }
176  template<typename IM>
177  void init( IM const& im )
178  {
179  M_tensor_expr.init( im );
180  }
181  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
182  {
183  update( geom );
184  }
185  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
186  {
187  update( geom );
188  }
189  void update( Geo_t const& geom )
190  {
191  M_tensor_expr.update( geom );
192  computeInv( mpl::int_<shape::N>() );
193  }
194  void update( Geo_t const& geom, uint16_type face )
195  {
196  M_tensor_expr.update( geom, face );
197  computeInv( mpl::int_<shape::N>() );
198  }
199 
200 
201  value_type
202  evalij( uint16_type i, uint16_type j ) const
203  {
204  return M_tensor_expr.evalij( i, j );
205  }
206 
207 
208  value_type
209  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type c2, uint16_type q ) const
210  {
211  return evalq( c1, c2, q );
212  }
213 
214  value_type
215  evaliq( uint16_type /*i*/, uint16_type c1, uint16_type c2, uint16_type q ) const
216  {
217  return evalq( c1, c2, q );
218  }
219 
220  value_type
221  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
222  {
223  return M_inv[q](c1,c2);
224  }
225 
226  private:
227  void
228  computeInv( mpl::int_<1> )
229  {
230  for( int q = 0; q < M_inv.size(); ++q )
231  {
232  M_inv[q](0,0) = 1./M_tensor_expr.evalq( 0, 0, q );
233  }
234  }
235  void
236  computeInv( mpl::int_<2> )
237  {
238  for( int q = 0; q < M_inv.size(); ++q )
239  {
240  double a = M_tensor_expr.evalq( 0, 0, q );
241  double b = M_tensor_expr.evalq( 0, 1, q );
242  double c = M_tensor_expr.evalq( 1, 0, q );
243  double d = M_tensor_expr.evalq( 1, 1, q );
244 
245  double determinant = a*d-c*b;
246 
247  M_inv[q](0,0) = d/determinant;
248  M_inv[q](0,1) = -b/determinant;
249  M_inv[q](1,0) = -c/determinant;
250  M_inv[q](1,1) = a/determinant;
251  }
252  }
253  void
254  computeInv( mpl::int_<3> )
255  {
256  for( int q = 0; q < M_inv.size(); ++q )
257  {
258  double a = M_tensor_expr.evalq( 0, 0, q );
259  double b = M_tensor_expr.evalq( 0, 1, q );
260  double c = M_tensor_expr.evalq( 0, 2, q );
261  double d = M_tensor_expr.evalq( 1, 0, q );
262  double e = M_tensor_expr.evalq( 1, 1, q );
263  double f = M_tensor_expr.evalq( 1, 2, q );
264  double g = M_tensor_expr.evalq( 2, 0, q );
265  double h = M_tensor_expr.evalq( 2, 1, q );
266  double l = M_tensor_expr.evalq( 2, 2, q );
267 
268  double determinant = a*(e*l-f*h)-b*(d*l-g*h)+c*(d*h-g*e);
269 
270  M_inv[q](0,0) = (e*l-f*h)/determinant;
271  M_inv[q](0,1) = (c*h-b*l)/determinant;
272  M_inv[q](0,2) = (b*f-c*e)/determinant;
273  M_inv[q](1,0) = (f*g-d*l)/determinant;
274  M_inv[q](1,1) = (a*l-c*g)/determinant;
275  M_inv[q](1,2) = (c*d-a*f)/determinant;
276  M_inv[q](2,0) = (d*h-e*g)/determinant;
277  M_inv[q](2,1) = (b*g-a*h)/determinant;
278  M_inv[q](2,2) = (a*e-b*d)/determinant;
279  }
280  }
281  private:
282  tensor_expr_type M_tensor_expr;
283  inv_matrix_type M_inv;
284  };
285 
286 private:
287  mutable expression_type M_expr;
288 };
290 
294 template<typename ExprT>
295 inline
296 Expr< Inv<ExprT> >
297 inv( ExprT v )
298 {
299  typedef Inv<ExprT> inv_t;
300  return Expr< inv_t >( inv_t( v ) );
301 }
302 
303 }
304 }
305 #endif /* __Inv_H */

Generated on Sun Oct 20 2013 08:25:00 for Feel++ by doxygen 1.8.4