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
det.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_Det_H
30 #define __FEELPP_VF_Det_H 1
31 
32 namespace Feel
33 {
34 namespace vf
35 {
37 
44 template<typename ExprT>
45 class Det
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 Det<ExprT> this_type;
75 
76 
78 
82 
83  explicit Det( expression_type const & __expr )
84  :
85  M_expr( __expr )
86  {}
87  Det( Det const & te )
88  :
89  M_expr( te.M_expr )
90  {}
91  ~Det()
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_ptr<vf::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,Scalar,false,false> shape;
138 
139  typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> vector_type;
140 
141  template <class Args> struct sig
142  {
143  typedef value_type type;
144  };
145 
146  struct is_zero
147  {
148  static const bool value = tensor_expr_type::is_zero::value;
149  };
150 
151  tensor( this_type const& expr,
152  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
153  :
154  M_tensor_expr( expr.expression(), geom, fev, feu ),
155  M_det( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
156  {
157  }
158 
159  tensor( this_type const& expr,
160  Geo_t const& geom, Basis_i_t const& fev )
161  :
162  M_tensor_expr( expr.expression(), geom, fev ),
163  M_det( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
164  {
165  }
166 
167  tensor( this_type const& expr, Geo_t const& geom )
168  :
169  M_tensor_expr( expr.expression(), geom ),
170  M_det( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
171  {
172  }
173  template<typename IM>
174  void init( IM const& im )
175  {
176  M_tensor_expr.init( im );
177  }
178  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
179  {
180  update(geom);
181  }
182  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
183  {
184  update(geom);
185  }
186  void update( Geo_t const& geom )
187  {
188  M_tensor_expr.update( geom );
189  computeDet( mpl::int_<expr_shape::nDim>() );
190  }
191  void update( Geo_t const& geom, uint16_type face )
192  {
193  M_tensor_expr.update( geom, face );
194  computeDet( mpl::int_<expr_shape::nDim>() );
195  }
196 
197 
198  value_type
199  evalij( uint16_type i, uint16_type j ) const
200  {
201  return M_tensor_expr.evalij( i, j );
202  }
203 
204 
205  value_type
206  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type c2, uint16_type q ) const
207  {
208  return evalq( c1, c2, q );
209  }
210 
211  value_type
212  evaliq( uint16_type /*i*/, uint16_type c1, uint16_type c2, uint16_type q ) const
213  {
214  return evalq( c1, c2, q );
215  }
216 
217  value_type
218  evalq( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q ) const
219  {
220  return M_det(q);
221  }
222 
223  private:
224  void
225  computeDet( mpl::int_<1> )
226  {
227  for( int q = 0; q < M_det.rows(); ++q )
228  {
229  M_det(q) = M_tensor_expr.evalq( 0, 0, q );
230  }
231  }
232  void
233  computeDet( mpl::int_<2> )
234  {
235  for( int q = 0; q < M_det.rows(); ++q )
236  {
237  double a = M_tensor_expr.evalq( 0, 0, q );
238  double b = M_tensor_expr.evalq( 0, 1, q );
239  double c = M_tensor_expr.evalq( 1, 0, q );
240  double d = M_tensor_expr.evalq( 1, 1, q );
241  M_det(q) = a*d-b*c;
242  }
243  }
244  void
245  computeDet( mpl::int_<3> )
246  {
247 
248  for( int q = 0; q < M_det.rows(); ++q )
249  {
250  double a = M_tensor_expr.evalq( 0, 0, q );
251  double b = M_tensor_expr.evalq( 0, 1, q );
252  double c = M_tensor_expr.evalq( 0, 2, q );
253  double d = M_tensor_expr.evalq( 1, 0, q );
254  double e = M_tensor_expr.evalq( 1, 1, q );
255  double f = M_tensor_expr.evalq( 1, 2, q );
256  double g = M_tensor_expr.evalq( 2, 0, q );
257  double h = M_tensor_expr.evalq( 2, 1, q );
258  double i = M_tensor_expr.evalq( 2, 2, q );
259  M_det(q) = a*e*i+b*f*g+c*d*h - (c*e*g+b*d*i+a*f*h);
260  }
261  }
262  private:
263  tensor_expr_type M_tensor_expr;
264  vector_type M_det;
265  };
266 
267 private:
268  mutable expression_type M_expr;
269 };
271 
275 template<typename ExprT>
276 inline
277 Expr< Det<ExprT> >
278 det( ExprT v )
279 {
280  typedef Det<ExprT> det_t;
281  return Expr< det_t >( det_t( v ) );
282 }
283 
284 }
285 }
286 #endif /* __Det_H */

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