29 #ifndef __FEELPP_VF_Inv_H
30 #define __FEELPP_VF_Inv_H 1
44 template<
typename ExprT>
49 static const size_type context = ExprT::context;
50 static const bool is_terminal =
false;
52 static const uint16_type imorder = ExprT::imorder;
53 static const bool imIsPoly = ExprT::imIsPoly;
55 template<
typename Func>
56 struct HasTestFunction
58 static const bool result = ExprT::template HasTestFunction<Func>::result;
61 template<
typename Func>
62 struct HasTrialFunction
64 static const bool result = ExprT::template HasTrialFunction<Func>::result;
72 typedef ExprT expression_type;
73 typedef typename expression_type::value_type value_type;
74 typedef Inv<ExprT> this_type;
83 explicit Inv( expression_type
const & __expr )
121 expression_type
const& expression()
const
129 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t>
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;
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;
140 typedef Eigen::Matrix<value_type,shape::M,shape::N> matrix_type;
141 typedef std::vector<matrix_type> inv_matrix_type;
144 template <
class Args>
struct sig
146 typedef value_type type;
151 static const bool value = tensor_expr_type::is_zero::value;
154 tensor( this_type
const& expr,
155 Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
157 M_tensor_expr( expr.expression(), geom, fev, feu ),
158 M_inv( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
162 tensor( this_type
const& expr,
163 Geo_t
const& geom, Basis_i_t
const& fev )
165 M_tensor_expr( expr.expression(), geom, fev ),
166 M_inv( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
170 tensor( this_type
const& expr, Geo_t
const& geom )
172 M_tensor_expr( expr.expression(), geom ),
173 M_inv( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
176 template<
typename IM>
177 void init( IM
const& im )
179 M_tensor_expr.init( im );
181 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
185 void update( Geo_t
const& geom, Basis_i_t
const& )
189 void update( Geo_t
const& geom )
191 M_tensor_expr.update( geom );
192 computeInv( mpl::int_<shape::N>() );
194 void update( Geo_t
const& geom, uint16_type face )
196 M_tensor_expr.update( geom, face );
197 computeInv( mpl::int_<shape::N>() );
202 evalij( uint16_type i, uint16_type j )
const
204 return M_tensor_expr.evalij( i, j );
209 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
211 return evalq( c1, c2, q );
215 evaliq( uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
217 return evalq( c1, c2, q );
221 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
223 return M_inv[q](c1,c2);
228 computeInv( mpl::int_<1> )
230 for(
int q = 0; q < M_inv.size(); ++q )
232 M_inv[q](0,0) = 1./M_tensor_expr.evalq( 0, 0, q );
236 computeInv( mpl::int_<2> )
238 for(
int q = 0; q < M_inv.size(); ++q )
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 );
245 double determinant = a*d-c*b;
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;
254 computeInv( mpl::int_<3> )
256 for(
int q = 0; q < M_inv.size(); ++q )
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 );
268 double determinant = a*(e*l-f*h)-b*(d*l-g*h)+c*(d*h-g*e);
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;
282 tensor_expr_type M_tensor_expr;
283 inv_matrix_type M_inv;
287 mutable expression_type M_expr;
294 template<
typename ExprT>
299 typedef Inv<ExprT> inv_t;
300 return Expr< inv_t >( inv_t( v ) );
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319