29 #ifndef __FEELPP_VF_Det_H
30 #define __FEELPP_VF_Det_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 Det<ExprT> this_type;
83 explicit Det( 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,Scalar,false,false> shape;
139 typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> vector_type;
141 template <
class Args>
struct sig
143 typedef value_type type;
148 static const bool value = tensor_expr_type::is_zero::value;
151 tensor( this_type
const& expr,
152 Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
154 M_tensor_expr( expr.expression(), geom, fev, feu ),
155 M_det( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
159 tensor( this_type
const& expr,
160 Geo_t
const& geom, Basis_i_t
const& fev )
162 M_tensor_expr( expr.expression(), geom, fev ),
163 M_det( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
167 tensor( this_type
const& expr, Geo_t
const& geom )
169 M_tensor_expr( expr.expression(), geom ),
170 M_det( vf::detail::ExtractGm<Geo_t>::get( geom )->nPoints() )
173 template<
typename IM>
174 void init( IM
const& im )
176 M_tensor_expr.init( im );
178 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
182 void update( Geo_t
const& geom, Basis_i_t
const& )
186 void update( Geo_t
const& geom )
188 M_tensor_expr.update( geom );
189 computeDet( mpl::int_<expr_shape::nDim>() );
191 void update( Geo_t
const& geom, uint16_type face )
193 M_tensor_expr.update( geom, face );
194 computeDet( mpl::int_<expr_shape::nDim>() );
199 evalij( uint16_type i, uint16_type j )
const
201 return M_tensor_expr.evalij( i, j );
206 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
208 return evalq( c1, c2, q );
212 evaliq( uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
214 return evalq( c1, c2, q );
218 evalq( uint16_type , uint16_type , uint16_type q )
const
225 computeDet( mpl::int_<1> )
227 for(
int q = 0; q < M_det.rows(); ++q )
229 M_det(q) = M_tensor_expr.evalq( 0, 0, q );
233 computeDet( mpl::int_<2> )
235 for(
int q = 0; q < M_det.rows(); ++q )
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 );
245 computeDet( mpl::int_<3> )
248 for(
int q = 0; q < M_det.rows(); ++q )
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);
263 tensor_expr_type M_tensor_expr;
268 mutable expression_type M_expr;
275 template<
typename ExprT>
280 typedef Det<ExprT> det_t;
281 return Expr< det_t >( det_t( v ) );