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
operators.hpp
Go to the documentation of this file.
1 
2 /* -*- 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
3 
4  This file is part of the Feel library
5 
6  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
7  Date: 2005-01-17
8 
9  Copyright (C) 2005,2006 EPFL
10  Copyright (C) 2006,2007 Université Joseph Fourier (Grenoble I)
11 
12  This library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU Lesser General Public
14  License as published by the Free Software Foundation; either
15  version 3.0 of the License, or (at your option) any later version.
16 
17  This library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  Lesser General Public License for more details.
21 
22  You should have received a copy of the GNU Lesser General Public
23  License along with this library; if not, write to the Free Software
24  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 */
31 #if !defined( __FEELPP_VF_OPERATORS_HPP )
32 #define __FEELPP_VF_OPERATORS_HPP 1
33 
35 
36 # include <boost/preprocessor/comparison/less.hpp>
37 # include <boost/preprocessor/logical/and.hpp>
38 # include <boost/preprocessor/control/if.hpp>
39 # include <boost/preprocessor/list/at.hpp>
40 # include <boost/preprocessor/list/cat.hpp>
41 # include <boost/preprocessor/list/for_each_product.hpp>
42 # include <boost/preprocessor/logical/or.hpp>
43 # include <boost/preprocessor/tuple/to_list.hpp>
44 # include <boost/preprocessor/tuple/eat.hpp>
45 # include <boost/preprocessor/facilities/empty.hpp>
46 # include <boost/preprocessor/punctuation/comma.hpp>
47 # include <boost/preprocessor/facilities/identity.hpp>
48 # include <boost/preprocessor/stringize.hpp>
49 
50 
51 namespace Feel
52 {
53 namespace vf
54 {
56 # /* Accessors for the operator datatype. */
57 # define VF_OPERATOR_NAME(O) BOOST_PP_TUPLE_ELEM(11, 0, O)
58 # define VF_OPERATOR_SYMBOL(O) BOOST_PP_TUPLE_ELEM(11, 1, O)
59 # define VF_OPERATOR_TERM(O) BOOST_PP_TUPLE_ELEM(11, 2, O)
60 # // 0 means any rank, 1 or more is the rank for which the operator is defined
61 # define VF_OPERATOR_TYPE_RANK_DEF(O) BOOST_PP_TUPLE_ELEM(11, 3, O)
62 # define VF_OPERATOR_TYPE_IS_COMP(O) BOOST_PP_TUPLE_ELEM(11, 4, O)
63 # define VF_OPERATOR_TYPE_COMP(O) BOOST_PP_TUPLE_ELEM(11, 5, O)
64 # define VF_OPERATOR_CONTEXT(O) BOOST_PP_TUPLE_ELEM(11, 6, O)
65 # define VF_OPERATOR_RANK(O) BOOST_PP_TUPLE_ELEM(11, 7, O)
66 # define VF_OPERATOR_TRANSPOSE(O) BOOST_PP_TUPLE_ELEM(11, 8, O)
67 # define VF_OPERATOR_DIFFORDERIM(O) BOOST_PP_TUPLE_ELEM(11, 9, O)
68 # define VF_OPERATOR_TERMINAL(O) BOOST_PP_TUPLE_ELEM(11,10, O)
69 
70 
71 
72 # /* List of applicative operators. */
73 # define VF_OPERATORS \
74  BOOST_PP_TUPLE_TO_LIST( \
75  12, \
76  ( \
77  ( OpId , id , id , 0, 0, 0, vm::JACOBIAN , RankSame,false, 0, 1 ), \
78  ( OpDx , dx , dx , 0, 1, 0, vm::JACOBIAN|vm::KB|vm::GRAD , RankSame,false,-1,1 ), \
79  ( OpDy , dy , dy , 0, 1, 1, vm::JACOBIAN|vm::KB|vm::GRAD , RankSame,false,-1,1 ), \
80  ( OpDz , dz , dz , 0, 1, 2, vm::JACOBIAN|vm::KB|vm::GRAD , RankSame,false,-1,1 ), \
81  ( OpDn , dn , dn , 0, 0, 0, vm::JACOBIAN|vm::KB|vm::NORMAL|vm::FIRST_DERIVATIVE|vm::FIRST_DERIVATIVE_NORMAL , RankSame,false,-1,1 ), \
82  ( OpGrad , grad , grad , 0, 0, 0, vm::JACOBIAN|vm::KB|vm::GRAD , RankUp,true,-1,1 ), \
83  ( OpDiv , div , div , 1, 0, 0, vm::DIV|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE , RankDown,false,-1,1 ), \
84  ( OpCurl , curl , curl , 1, 0, 0, vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE , RankSame,false,-1,1 ), \
85  ( OpCurlX, curlx, curlx, 1, 1, 0, vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE , RankDown,false,-1,1 ), \
86  ( OpCurlY, curly, curly, 1, 1, 1, vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE , RankDown,false,-1,1 ), \
87  ( OpCurlZ, curlz, curlz, 1, 1, 2, vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE , RankDown,false,-1,1 ), \
88  ( OpHess , hess , hess, 0, 0, 0, vm::JACOBIAN|vm::KB|vm::HESSIAN|vm::FIRST_DERIVATIVE , RankUp2,false,-2,1 ) \
89  ) \
90  ) \
91 
92 #
93 enum OperatorType { __TEST, __TRIAL, __VALUE };
94 # define VF_OP_TYPE_TYPE(T) BOOST_PP_TUPLE_ELEM(6, 0, T)
95 # define VF_OP_TYPE_IS_GENERIC(T) BOOST_PP_TUPLE_ELEM(6, 1, T)
96 # define VF_OP_TYPE_IS_TEST(T) BOOST_PP_TUPLE_ELEM(6, 2, T)
97 # define VF_OP_TYPE_IS_TRIAL(T) BOOST_PP_TUPLE_ELEM(6, 3, T)
98 # define VF_OP_TYPE_IS_VALUE(T) BOOST_PP_TUPLE_ELEM(6, 4, T)
99 # define VF_OP_TYPE_SUFFIX(T) BOOST_PP_TUPLE_ELEM(6, 5, T)
100 # define VF_OPERATORS_TYPE \
101  BOOST_PP_TUPLE_TO_LIST( \
102  4, \
103  ( \
104  ( OperatorType, 1, 0, 0, 0, g ),\
105  ( __TEST , 0, 1, 0, 0, ),\
106  ( __TRIAL , 0, 0, 1, 0, t ),\
107  ( __VALUE , 0, 0, 0, 1, v ) \
108  ) \
109  ) \
110 
111 #
112 
113 # define VF_TRIAL_NAME(T) BOOST_PP_TUPLE_ELEM(3, 0, T)
114 # define VF_TRIAL_IS_GENERIC(T) BOOST_PP_TUPLE_ELEM(3, 1, T)
115 # define VF_TRIAL_VALUE(T) BOOST_PP_TUPLE_ELEM(3, 2, T)
116 #
117 
118 # define VF_TRIAL \
119  BOOST_PP_TUPLE_TO_LIST( \
120  3, \
121  ( \
122  ( bool, 1, 0 ),\
123  ( false, 0, 0 ),\
124  ( true , 0, 1 ) \
125  ) \
126  ) \
127 
128 #
129 # define VF_OP_SPECIALIZATION_IF_NOT_GENERIC( E, T ) \
130  BOOST_PP_IF ( BOOST_PP_NOT( VF_OP_TYPE_IS_GENERIC( T ) ),\
131  BOOST_PP_IDENTITY( <E ), \
132  BOOST_PP_EMPTY )() \
133  BOOST_PP_IF ( BOOST_PP_NOT( VF_OP_TYPE_IS_GENERIC( T ) ),\
134  BOOST_PP_COMMA, \
135  BOOST_PP_EMPTY )() \
136  BOOST_PP_IF ( BOOST_PP_NOT( VF_OP_TYPE_IS_GENERIC( T ) ),\
137  BOOST_PP_IDENTITY( VF_OP_TYPE_TYPE( T )> ),\
138  BOOST_PP_EMPTY )() \
139 
140 #
141 # define VF_OP_SWITCH( T, A, B ) \
142  BOOST_PP_IF ( T, \
143  BOOST_PP_IDENTITY( A ), \
144  BOOST_PP_IDENTITY( B ) \
145  )() \
146 
147 # define VF_OP_SWITCH_ELSE_EMPTY( T, A ) \
148  BOOST_PP_IF ( T, \
149  BOOST_PP_IDENTITY( A ), \
150  BOOST_PP_EMPTY \
151  )() \
152 
153 #
154 # define VF_OP_ADD_COMP( O ) \
155  BOOST_PP_IF ( VF_OPERATOR_TYPE_IS_COMP( O ), \
156  BOOST_PP_IDENTITY( ( VF_OPERATOR_TYPE_COMP( O ) ) ), \
157  BOOST_PP_EMPTY \
158  )() \
159 
160 #
161 #
162 #define VF_OP_TYPE_OBJECT(T) \
163  BOOST_PP_IF ( VF_OP_TYPE_IS_GENERIC( T ), \
164  BOOST_PP_IDENTITY( sw ), \
165  BOOST_PP_IDENTITY( VF_OP_TYPE_TYPE( T ) ) )() \
166 
167 #
168 #
169 # /* Generates code for all binary operators and integral type pairs. */
170 # define VF_ARRAY_OPERATOR(_, OT) \
171  VF_ARRAY_OPERATOR_CODE OT \
172 
173 #define VF_FUSION1(T) BOOST_PP_IF( BOOST_PP_NOT(VF_OP_TYPE_IS_VALUE( T )), \
174  BOOST_PP_IDENTITY(typedef typename fusion::result_of::value_at_key<map_basis_context_type), \
175  BOOST_PP_EMPTY )() \
176  BOOST_PP_IF( BOOST_PP_NOT(VF_OP_TYPE_IS_VALUE( T )), \
177  BOOST_PP_COMMA, \
178  BOOST_PP_EMPTY )() \
179  BOOST_PP_IF( BOOST_PP_NOT(VF_OP_TYPE_IS_VALUE( T )), \
180  BOOST_PP_IDENTITY(key_type>::type::element_type basis_context_type), \
181  BOOST_PP_IDENTITY( typedef boost::none_t basis_context_type ))() \
182 
183 
184 #define VF_FUSION2(T) BOOST_PP_IF( BOOST_PP_NOT(VF_OP_TYPE_IS_VALUE( T )), \
185  BOOST_PP_IDENTITY(typedef typename fusion::result_of::value_at_key<map_basis_context_type), \
186  BOOST_PP_EMPTY )() \
187  BOOST_PP_IF( BOOST_PP_NOT(VF_OP_TYPE_IS_VALUE( T )), \
188  BOOST_PP_COMMA, \
189  BOOST_PP_EMPTY )() \
190  BOOST_PP_IF( BOOST_PP_NOT(VF_OP_TYPE_IS_VALUE( T )), \
191  BOOST_PP_IDENTITY(key_type>::type::element_type* basis_context_ptrtype), \
192  BOOST_PP_IDENTITY( typedef boost::none_t basis_context_ptrtype ) )() \
193 
194 
195 #define VF_ARRAY_OPERATOR_CODE(O,T) \
196  template <class Element \
197  BOOST_PP_IF( VF_OP_TYPE_IS_GENERIC( T ), BOOST_PP_COMMA, BOOST_PP_EMPTY )() \
198  BOOST_PP_IF( VF_OP_TYPE_IS_GENERIC( T ), BOOST_PP_IDENTITY( VF_OP_TYPE_TYPE( T ) sw ), BOOST_PP_EMPTY )() > \
199  class VF_OPERATOR_NAME( O ) VF_OP_SPECIALIZATION_IF_NOT_GENERIC( Element, T ) \
200  { \
201  public: \
202  \
203  static const size_type context = VF_OPERATOR_CONTEXT( O ); \
204  \
205  typedef Element element_type; \
206  typedef boost::shared_ptr<element_type> element_ptrtype; \
207  typedef VF_OPERATOR_NAME( O )<element_type, VF_OP_TYPE_OBJECT(T)> this_type; \
208  typedef this_type self_type; \
209  \
210  typedef typename element_type::functionspace_type functionspace_type; \
211  typedef typename functionspace_type::reference_element_type* fe_ptrtype; \
212  typedef typename functionspace_type::reference_element_type fe_type; \
213  typedef typename functionspace_type::geoelement_type geoelement_type; \
214  typedef typename functionspace_type::gm_type gm_type; \
215  typedef typename functionspace_type::value_type value_type; \
216  static const uint16_type rank = fe_type::rank; \
217  static const uint16_type nComponents1 = fe_type::nComponents1; \
218  static const uint16_type nComponents2 = fe_type::nComponents2; \
219  static const bool is_terminal = VF_OPERATOR_TERMINAL(O); \
220  \
221  static const uint16_type imorder_test=element_type::functionspace_type::basis_type::nOrder + VF_OPERATOR_DIFFORDERIM(O); \
222  static const uint16_type imorder = (imorder_test==invalid_uint16_type_value)?0:imorder_test; \
223  static const bool imIsPoly = true; \
224  \
225  template<typename Func> \
226  struct HasTestFunction \
227  { \
228  static const bool result = VF_OP_SWITCH( BOOST_PP_OR( VF_OP_TYPE_IS_TRIAL( T ), VF_OP_TYPE_IS_VALUE( T ) ), false , (boost::is_same<Func,fe_type>::value) ); \
229  }; \
230  \
231  template<typename Func> \
232  struct HasTrialFunction \
233  { \
234  static const bool result = VF_OP_SWITCH( VF_OP_TYPE_IS_TRIAL( T ), (boost::is_same<Func,fe_type>::value), false ); \
235  }; \
236  \
237  \
238  VF_OPERATOR_NAME( O ) ( element_type const& v ) \
239  : M_v ( boost::cref(v) ) \
240  { \
241  if ( VF_OP_TYPE_IS_VALUE( T ) ) \
242  v.updateGlobalValues(); \
243  DVLOG(2) << "[" BOOST_PP_STRINGIZE(VF_OPERATOR_NAME( O )) "] default constructor\n"; \
244  } \
245  VF_OPERATOR_NAME( O )( VF_OPERATOR_NAME( O ) const& op ) \
246  : M_v ( op.M_v ) \
247  { \
248  DVLOG(2) << "[" BOOST_PP_STRINGIZE(VF_OPERATOR_NAME( O )) "] copy constructor\n"; \
249  } \
250  template<typename TheExpr> \
251  struct Lambda \
252  { \
253  typedef this_type type; \
254  }; \
255  template<typename TheExpr> \
256  typename Lambda<TheExpr>::type \
257  operator()( TheExpr const& e ) { return *this; } \
258  \
259  element_type const& e() const { return M_v; } \
260  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t> \
261  struct tensor \
262  { \
263  typedef this_type expression_type; \
264  typedef BOOST_PP_CAT( Basis_,BOOST_PP_CAT(VF_OP_SWITCH( BOOST_PP_NOT( VF_OP_TYPE_IS_TRIAL( T ) ), i ,j ), _t)) map_basis_context_type; \
265  typedef typename mpl::if_<mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>, \
266  typename mpl::if_<fusion::result_of::has_key<Geo_t,vf::detail::gmc<0> >, \
267  mpl::identity<vf::detail::gmc<0> >, \
268  mpl::identity<vf::detail::gmc<1> > >::type, \
269  typename mpl::if_<fusion::result_of::has_key<map_basis_context_type,vf::detail::gmc<0> >, \
270  mpl::identity<vf::detail::gmc<0> >, \
271  mpl::identity<vf::detail::gmc<1> > >::type>::type::type key_type; \
272  typedef typename mpl::if_<fusion::result_of::has_key<map_basis_context_type,vf::detail::gmc<0> >, \
273  mpl::identity<vf::detail::gmc<0> >, \
274  mpl::identity<vf::detail::gmc<1> > >::type::type basis_context_key_type; \
275  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type; \
276  typedef boost::shared_ptr<gmc_type> gmc_ptrtype; \
277  typedef typename gmc_type::gm_type gm_type; \
278  BOOST_MPL_ASSERT_MSG( ( fusion::result_of::has_key<map_basis_context_type, basis_context_key_type >::value ), INVALID_BASISMAP_OP, (map_basis_context_type, key_type, basis_context_key_type, mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>, mpl::bool_<VF_OP_TYPE_IS_TRIAL( T )> )); \
279  typedef typename mpl::if_<mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>, \
280  mpl::identity<mpl::int_<0> >, \
281  mpl::identity<typename fusion::result_of::value_at_key<map_basis_context_type,basis_context_key_type>::type::element_type > >::type::type basis_context_type; \
282  typedef typename mpl::if_<mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>, \
283  mpl::identity<mpl::int_<0> >, \
284  mpl::identity<typename fusion::result_of::value_at_key<map_basis_context_type,basis_context_key_type>::type::element_type* > >::type::type basis_context_ptrtype; \
285  typedef typename element_type::value_type value_type; \
286  typedef typename matrix_node<value_type>::type matrix_node_type; \
287  \
288  typedef value_type result_type; \
289  typedef typename element_type::polyset_type function_rank_type; \
290  typedef typename VF_OPERATOR_RANK( O )<function_rank_type>::type return_value_type; \
291  \
292  typedef typename mpl::if_<mpl::equal_to<mpl::int_<return_value_type::rank>, \
293  mpl::int_<0> >, \
294  mpl::identity<Shape<gmc_type::NDim, Scalar, false> >, \
295  typename mpl::if_<mpl::equal_to<mpl::int_<return_value_type::rank>, \
296  mpl::int_<1> >, \
297  mpl::identity<Shape<gmc_type::NDim, Vectorial, VF_OPERATOR_TRANSPOSE(O)> >, \
298  mpl::identity<Shape<gmc_type::NDim, Tensor2, false> > >::type>::type::type shape; \
299  typedef typename fe_type::PreCompute pc_type; \
300  typedef boost::shared_ptr<pc_type> pc_ptrtype; \
301  typedef typename fe_type::template Context<context, fe_type, gm_type,geoelement_type,gmc_type::context> ctx_type; \
302  typedef boost::shared_ptr<ctx_type> ctx_ptrtype; \
303  typedef Eigen::Matrix<value_type,shape::M,shape::N> loc_type; \
304  typedef Eigen::Matrix<value_type,shape::M,shape::N> ret_type; \
305  typedef boost::multi_array<loc_type,1> array_type; \
306  \
307  \
308  template<typename E>\
309  struct ttt { \
310  typedef typename mpl::if_<boost::is_same<E,mpl::int_<0> >, \
311  mpl::identity<functionspace_type>, \
312  mpl::identity<basis_context_type> >::type::type type; \
313  }; \
314  static const bool dim_ok = (VF_OPERATOR_TYPE_COMP(O) < gmc_type::NDim); \
315  static const bool fe_ok = mpl::if_<mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>, \
316  mpl::bool_<true>, \
317  boost::is_same<typename ttt<basis_context_type>::type::reference_element_type, fe_type> >::type::value; \
318  struct is_zero { \
319  /*static const bool value = !(dim_ok && fe_ok);*/ \
320  static const bool value = false; \
321  }; \
322  \
323  static const uint16_type rank = return_value_type::rank+1; \
324  static const uint16_type nComponents = return_value_type::nComponents; \
325  \
326  static const bool isSameGeo = boost::is_same<typename gmc_type::element_type,geoelement_type>::value; \
327  \
328  tensor( tensor const& t ) \
329  : \
330  M_expr( t.M_expr ), \
331  M_geot( new gmc_type( *t.M_geot ) ), \
332  M_fec( VF_OP_SWITCH( VF_OP_TYPE_IS_VALUE( T ), , t.M_fec/*new basis_context_type( *t.M_fec )*/ ) ), \
333  M_np( M_geot->nPoints() ), \
334  M_pc( new pc_type( M_expr.e().functionSpace()->fe(), M_geot->xRefs() )), \
335  /*M_pcf(),*/ \
336  M_ctx( VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), (new ctx_type( M_expr.e().functionSpace()->fe(), M_geot, (pc_ptrtype const&)M_pc ) ) ) ), \
337  M_loc(VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), M_expr.e().BOOST_PP_CAT(VF_OPERATOR_TERM( O ),Extents)(*M_geot) ) ), \
338  M_zero( ret_type::Zero() ), \
339  M_did_init( t.M_did_init ), \
340  M_same_mesh( t.M_same_mesh ) \
341  { \
342  } \
343  \
344  tensor( this_type const& expr, \
345  Geo_t const& geom, \
346  Basis_i_t const& VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TEST( T ), fev ), \
347  Basis_j_t const& VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TRIAL( T ), feu ) ) \
348  : \
349  M_expr( expr ), \
350  M_geot( fusion::at_key<key_type>( geom ) ), \
351  M_fec( VF_OP_SWITCH( VF_OP_TYPE_IS_TEST( T ), \
352  fusion::at_key<basis_context_key_type>( fev ).get() , \
353  VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TRIAL( T ), \
354  fusion::at_key<basis_context_key_type>( feu ).get() ) ) ), \
355  M_np( fusion::at_key<key_type>( geom )->nPoints() ), \
356  M_pc( this->createPcIfSameGeom(expr,geom, mpl::bool_<isSameGeo>() ) ) /*new pc_type( expr.e().functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ))*/, \
357  /*M_pcf(),*/ \
358  M_ctx( VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), \
359  ( this->createCtxIfSameGeom(expr,geom, mpl::bool_<isSameGeo>() )) ) ), \
360  M_loc(VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), expr.e().BOOST_PP_CAT(VF_OPERATOR_TERM( O ),Extents)(*fusion::at_key<key_type>( geom )) ) ), \
361  M_zero( ret_type::Zero() ), \
362  M_did_init( false ), \
363  M_same_mesh( fusion::at_key<key_type>( geom )->element().mesh()->isRelatedTo( expr.e().functionSpace()->mesh()) && isSameGeo ) \
364  { \
365  if(!M_same_mesh) \
366  expr.e().functionSpace()->mesh()->tool_localization()->updateForUse(); \
367  /*update( geom );*/ \
368  } \
369  tensor( this_type const& expr, \
370  Geo_t const& geom, \
371  Basis_i_t const& VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TEST( T ), fev ) ) \
372  : \
373  M_expr( expr ), \
374  M_geot( fusion::at_key<key_type>( geom ) ), \
375  M_fec( VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TEST( T ), \
376  fusion::at_key<basis_context_key_type>( fev ).get() ) ), \
377  M_np( fusion::at_key<key_type>( geom )->nPoints() ), \
378  M_pc( this->createPcIfSameGeom(expr,geom, mpl::bool_<isSameGeo>() ) ) /*new pc_type( expr.e().functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ) )*/, \
379  /*M_pcf(),*/ \
380  M_ctx( VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), \
381  ( this->createCtxIfSameGeom(expr,geom, mpl::bool_<isSameGeo>() )) ) ), \
382  M_loc(VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), expr.e().BOOST_PP_CAT(VF_OPERATOR_TERM( O ),Extents)(*fusion::at_key<key_type>( geom )) ) ), \
383  M_zero( ret_type::Zero() ), \
384  M_did_init( false ), \
385  M_same_mesh( fusion::at_key<key_type>( geom )->element().mesh()->isRelatedTo( expr.e().functionSpace()->mesh()) && isSameGeo ) \
386  { \
387  if(!M_same_mesh) \
388  expr.e().functionSpace()->mesh()->tool_localization()->updateForUse(); \
389  /*update( geom );*/ \
390  } \
391  tensor( this_type const& expr, \
392  Geo_t const& geom ) \
393  : \
394  M_expr( expr ), \
395  M_geot( fusion::at_key<key_type>( geom ) ), \
396  M_np( fusion::at_key<key_type>( geom )->nPoints() ), \
397  M_pc( this->createPcIfSameGeom(expr,geom, mpl::bool_<isSameGeo>() ) ) /* new pc_type( expr.e().functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ) )*/, \
398  /*M_pcf(),*/ \
399  M_ctx( VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), \
400  ( this->createCtxIfSameGeom(expr,geom, mpl::bool_<isSameGeo>() )) ) ), \
401  M_loc(VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_VALUE( T ), expr.e().BOOST_PP_CAT(VF_OPERATOR_TERM( O ),Extents)(*fusion::at_key<key_type>( geom )) ) ), \
402  M_zero( ret_type::Zero() ), \
403  M_did_init( false ), \
404  M_same_mesh( fusion::at_key<key_type>( geom )->element().mesh()->isRelatedTo( expr.e().functionSpace()->mesh()) && isSameGeo ) \
405  { \
406  if(!M_same_mesh) \
407  expr.e().functionSpace()->mesh()->tool_localization()->updateForUse(); \
408  /*update( geom ); */ \
409  BOOST_MPL_ASSERT_MSG( VF_OP_TYPE_IS_VALUE( T ), INVALID_CALL_TO_CONSTRUCTOR, ()); \
410  } \
411  template<typename IM> \
412  void init( IM const& im ) \
413  { \
414  M_did_init = true; \
415  /* \
416  QuadMapped<IM> qm; \
417  typedef typename QuadMapped<IM>::permutation_type permutation_type; \
418  typename QuadMapped<IM>::permutation_points_type ppts( qm( im ) ); \
419  \
420  M_pcf.resize( im.nFaces() ); \
421  for ( uint16_type __f = 0; __f < im.nFaces(); ++__f ) \
422  { \
423  for( permutation_type __p( permutation_type::IDENTITY ); \
424  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p ) \
425  { \
426  M_pcf[__f][__p.value()] = pc_ptrtype( new pc_type( M_expr.e().functionSpace()->fe(), \
427  ppts[__f].find(__p)->second ) ); \
428  } \
429  } \
430  */ \
431  /* expr.e().functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ), */ \
432  } \
433  \
434  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu ) \
435  { \
436  update( geom, fev, feu, mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>() ); \
437  } \
438  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu , mpl::bool_<true> ) \
439  { \
440  update( geom, mpl::bool_<true>() ); \
441  } \
442  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu , mpl::bool_<false> ) \
443  { \
444  if (M_same_mesh) \
445  updateInCaseOfInterpolate( geom, fev, feu, mpl::bool_<false>() ); \
446  else \
447  updateInCaseOfInterpolate( geom, fev, feu, mpl::bool_<true>() ); \
448  } \
449  void updateInCaseOfInterpolate( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu , mpl::bool_<false> ) \
450  { \
451  /*nothing : always same context*/ \
452  } \
453  void updateInCaseOfInterpolate( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu , mpl::bool_<true> ) \
454  { /*with interp*/ \
455  VF_OP_SWITCH( VF_OP_TYPE_IS_TEST( T ), \
456  M_fec =fusion::at_key<basis_context_key_type>( fev ).get() , \
457  VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TRIAL( T ), \
458  M_fec = fusion::at_key<basis_context_key_type>( feu ).get() ) ) ; \
459  } \
460  void update( Geo_t const& geom, Basis_i_t const& fev ) \
461  { \
462  update( geom, fev, mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>() ); \
463  } \
464  void update( Geo_t const& geom, Basis_i_t const& fev, mpl::bool_<true> ) \
465  { \
466  update( geom, mpl::bool_<true>() ); \
467  } \
468  void update( Geo_t const& geom, Basis_i_t const& fev , mpl::bool_<false>) \
469  { \
470  if (M_same_mesh) \
471  updateInCaseOfInterpolate( geom, fev, mpl::bool_<false>() ); \
472  else \
473  updateInCaseOfInterpolate( geom, fev, mpl::bool_<true>() ); \
474  } \
475  void updateInCaseOfInterpolate( Geo_t const& geom, Basis_i_t const& fev, mpl::bool_<false> ) \
476  { \
477  /*no interp*/ \
478  } \
479  void updateInCaseOfInterpolate( Geo_t const& geom, Basis_i_t const& fev, mpl::bool_<true> ) \
480  { /*with interp*/ \
481  VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TEST( T ), \
482  M_fec = fusion::at_key<basis_context_key_type>( fev ).get() ) ; \
483  } \
484  void update( Geo_t const& geom ) \
485  { \
486  /*BOOST_STATIC_ASSERT( dim_ok );*/ \
487  update( geom, mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>() ); \
488  } \
489  void update( Geo_t const& geom, uint16_type face ) \
490  { \
491  /*BOOST_STATIC_ASSERT( dim_ok );*/ \
492  update( geom, face, mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>() ); \
493  } \
494  void update( Geo_t const& geom, uint16_type face1, mpl::bool_<true> ) \
495  { \
496  std::fill( M_loc.data(), M_loc.data()+M_loc.num_elements(), loc_type::Zero() ); \
497  this->updateCtxFaceIfSameGeom(geom,mpl::bool_<isSameGeo>() ); \
498  if (M_same_mesh) \
499  M_expr.e().VF_OPERATOR_SYMBOL( O )( *M_ctx, M_loc ); \
500  else { \
501  matrix_node_type __ptsreal = M_expr.e().ptsInContext(*fusion::at_key<key_type>( geom ), mpl::int_<2>()); \
502  M_expr.e().BOOST_PP_CAT(VF_OPERATOR_SYMBOL( O ),Interpolate)( __ptsreal, M_loc ); \
503  } \
504  } \
505  void update( Geo_t const& geom, mpl::bool_<true> ) \
506  { \
507  std::fill( M_loc.data(), M_loc.data()+M_loc.num_elements(), loc_type::Zero() ); \
508  this->updateCtxIfSameGeom(geom,mpl::bool_<isSameGeo>() ); \
509  if (M_same_mesh) { \
510  /*std::cout << "\n idv no interp \n";*/ \
511  M_expr.e().VF_OPERATOR_SYMBOL( O )( *M_ctx, M_loc ); \
512  } \
513  else { \
514  /*std::cout << "\n idv with interp \n";*/ \
515  matrix_node_type __ptsreal = M_expr.e().ptsInContext(*fusion::at_key<key_type>( geom ), mpl::int_<1>()); \
516  M_expr.e().BOOST_PP_CAT(VF_OPERATOR_SYMBOL( O ),Interpolate)( __ptsreal, M_loc ); \
517  } \
518  } \
519  void update( Geo_t const& geom, mpl::bool_<false> ) \
520  { \
521  /*std::cout << "\n idv no interp \n";*/ \
522  Feel::detail::ignore_unused_variable_warning(geom); \
523  } \
524  \
525  ret_type const& \
526  evalijq( uint16_type i, \
527  uint16_type VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TRIAL( T ), j ), \
528  uint16_type q ) const \
529  { \
530  Feel::detail::ignore_unused_variable_warning(i); \
531  return evaliq( VF_OP_SWITCH( BOOST_PP_NOT( VF_OP_TYPE_IS_TRIAL( T ) ),i,j), q ); \
532  } \
533  result_type \
534  evalijq( uint16_type i, \
535  uint16_type VF_OP_SWITCH_ELSE_EMPTY( VF_OP_TYPE_IS_TRIAL( T ), j ), \
536  uint16_type c1, uint16_type c2, uint16_type q ) const \
537  { \
538  Feel::detail::ignore_unused_variable_warning(i); \
539  return evaliq( VF_OP_SWITCH( BOOST_PP_NOT( VF_OP_TYPE_IS_TRIAL( T ) ),i,j), c1, c2, q ); \
540  } \
541  \
542  template<int PatternContext> \
543  result_type \
544  evalijq( uint16_type i, \
545  uint16_type j, \
546  uint16_type c1, uint16_type c2, uint16_type q, \
547  mpl::int_<PatternContext> ) const \
548  { \
549  return evalijq( i, j, c1, c2, q ); \
550  } \
551  \
552  ret_type const& \
553  evaliq( uint16_type i, uint16_type q ) const \
554  { \
555  return evaliq_( i, q, mpl::bool_<dim_ok && fe_ok>() ); \
556  } \
557  result_type \
558  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const \
559  { \
560  return evaliq_( i, c1, c2, q, mpl::bool_<dim_ok && fe_ok>() ); \
561  } \
562  result_type \
563  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const \
564  { \
565  BOOST_MPL_ASSERT_MSG( VF_OP_TYPE_IS_VALUE( T ), INVALID_CALL_TO_EVALQ, ()); \
566  return evalq( c1, c2, q, mpl::int_<shape::rank>() ); \
567  } \
568  ret_type const& \
569  evalq( uint16_type q ) const \
570  { \
571  BOOST_MPL_ASSERT_MSG( VF_OP_TYPE_IS_VALUE( T ), INVALID_CALL_TO_EVALQ, ()); \
572  return M_loc[q]; \
573  } \
574  private: \
575  \
576  result_type \
577  evaliq_( uint16_type /*i*/, \
578  uint16_type /*c1*/, uint16_type /*c2*/, \
579  int /*q*/, \
580  mpl::bool_<false> ) const \
581  { \
582  return 0; \
583  } \
584  ret_type const& \
585  evaliq_( uint16_type /*i*/, \
586  int /*q*/, \
587  mpl::bool_<false> ) const \
588  { \
589  return M_zero; \
590  } \
591  \
592  ret_type const& \
593  evaliq_( uint16_type i, uint16_type q, mpl::bool_<true> ) const \
594  { \
595  return evaliq__( i, q, mpl::bool_<true>(), mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>() ); \
596  } \
597  result_type \
598  evaliq_( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q, mpl::bool_<true> ) const \
599  { \
600  return evaliq__( i, c1, c2, q, mpl::bool_<true>(), mpl::bool_<VF_OP_TYPE_IS_VALUE( T )>() ); \
601  } \
602  \
603  ret_type const& \
604  evaliq__( uint16_type /*i*/, uint16_type q, \
605  mpl::bool_<true>, mpl::bool_<true> ) const \
606  { \
607  return M_loc[q]; \
608  } \
609  result_type \
610  evaliq__( uint16_type /*i*/, uint16_type c1, uint16_type c2, uint16_type q, \
611  mpl::bool_<true>, mpl::bool_<true> ) const \
612  { \
613  return evalq( c1, c2, q, mpl::int_<shape::rank>() ); \
614  } \
615  \
616  ret_type const& \
617  evaliq__( uint16_type i, uint16_type q, mpl::bool_<true>, mpl::bool_<false> ) const \
618  { \
619  return M_fec->VF_OPERATOR_TERM( O )( i, q ); \
620  } \
621  result_type \
622  evaliq__( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q, mpl::bool_<true>, mpl::bool_<false> ) const \
623  { \
624  return M_fec->VF_OPERATOR_TERM( O )( i, c1, c2, q ); \
625  } \
626  \
627  result_type \
628  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<0> ) const \
629  { \
630  Feel::detail::ignore_unused_variable_warning(c1); \
631  Feel::detail::ignore_unused_variable_warning(c2); \
632  return M_loc[q](0,0); \
633  } \
634  result_type \
635  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<1> ) const \
636  { \
637  return evalq( c1, c2, q, mpl::int_<1>(), mpl::bool_<shape::is_transposed>() ); \
638  } \
639  result_type \
640  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<1>, mpl::bool_<false> ) const \
641  { \
642  Feel::detail::ignore_unused_variable_warning(c1); \
643  Feel::detail::ignore_unused_variable_warning(c2); \
644  return M_loc[q](c1,0); \
645  } \
646  result_type \
647  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<1>, mpl::bool_<true> ) const \
648  { \
649  Feel::detail::ignore_unused_variable_warning(c1); \
650  Feel::detail::ignore_unused_variable_warning(c2); \
651  return M_loc[q](0,c2); \
652  } \
653  result_type \
654  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<2> ) const \
655  { \
656  return M_loc[q](c1,c2); \
657  } \
658  pc_ptrtype createPcIfSameGeom(this_type const& expr, Geo_t const& geom,mpl::bool_<true>) \
659  { \
660  return pc_ptrtype (new pc_type( expr.e().functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ) ); \
661  } \
662  pc_ptrtype createPcIfSameGeom(this_type const& expr, Geo_t const& geom,mpl::bool_<false>) \
663  { \
664  return pc_ptrtype(); \
665  } \
666  ctx_ptrtype createCtxIfSameGeom(this_type const& expr, Geo_t const& geom,mpl::bool_<true>) \
667  { \
668  return ctx_ptrtype( new ctx_type( expr.e().functionSpace()->fe(),fusion::at_key<key_type>( geom ),(pc_ptrtype const&)M_pc ) ); \
669  } \
670  ctx_ptrtype createCtxIfSameGeom(this_type const& expr, Geo_t const& geom,mpl::bool_<false>) \
671  { \
672  return ctx_ptrtype( /*new ctx_type( )*/ ); \
673  } \
674  void updateCtxIfSameGeom(Geo_t const& geom, mpl::bool_<true> ) \
675  { if (fusion::at_key<key_type>( geom )->faceId() != invalid_uint16_type_value ) /*face case*/ \
676  M_pc->update(fusion::at_key<key_type>( geom )->pc()->nodes() ); \
677  M_ctx->update( fusion::at_key<key_type>( geom ), (pc_ptrtype const&) M_pc ); \
678  } \
679  void updateCtxIfSameGeom(Geo_t const& geom, mpl::bool_<false> ) \
680  { \
681  } \
682  void updateCtxFaceIfSameGeom(Geo_t const& geom, mpl::bool_<true> ) \
683  { \
684  /*uint16_type face = fusion::at_key<key_type>( geom )->faceId(); \
685  uint16_type perm = fusion::at_key<key_type>( geom )->permutation().value(); \
686  M_ctx->update( fusion::at_key<key_type>( geom ), (pc_ptrtype const&) M_pcf[face][perm] );*/ \
687  M_pc->update(fusion::at_key<key_type>( geom )->pc()->nodes() ); \
688  M_ctx->update( fusion::at_key<key_type>( geom ), (pc_ptrtype const&) M_pc ); \
689  } \
690  void updateCtxFaceIfSameGeom(Geo_t const& geom, mpl::bool_<false> ) \
691  { \
692  } \
693  this_type const& M_expr; \
694  gmc_ptrtype M_geot; \
695  basis_context_ptrtype M_fec; \
696  const uint16_type M_np; \
697  pc_ptrtype M_pc; \
698  /*std::vector<std::map<uint16_type, pc_ptrtype> > M_pcf;*/ \
699  ctx_ptrtype M_ctx; \
700  array_type M_loc; \
701  ret_type M_zero; \
702  /*typename element_type::BOOST_PP_CAT( VF_OPERATOR_TERM( O ), _type) M_loc;*/ \
703  bool M_did_init; \
704  const bool M_same_mesh; \
705  }; \
706  \
707  protected: \
708  VF_OPERATOR_NAME( O ) () {} \
709  boost::reference_wrapper<const element_type> M_v; \
710  }; \
711  template <class ELEM \
712  BOOST_PP_IF( VF_OP_TYPE_IS_GENERIC( T ), BOOST_PP_COMMA, BOOST_PP_EMPTY )() \
713  BOOST_PP_IF( VF_OP_TYPE_IS_GENERIC( T ), BOOST_PP_IDENTITY( VF_OP_TYPE_TYPE( T ) sw ), BOOST_PP_EMPTY )() > \
714  inline Expr< VF_OPERATOR_NAME( O )< ELEM, VF_OP_TYPE_OBJECT(T)> > \
715  BOOST_PP_CAT( VF_OPERATOR_SYMBOL(O), VF_OP_TYPE_SUFFIX(T) )( ELEM const& expr ) \
716  { \
717  typedef VF_OPERATOR_NAME( O )< ELEM, VF_OP_TYPE_OBJECT(T)> expr_t; \
718  return Expr< expr_t >( expr_t(expr) ); \
719  } \
720  \
721  template <class ELEM \
722  BOOST_PP_IF( VF_OP_TYPE_IS_GENERIC( T ), BOOST_PP_COMMA, BOOST_PP_EMPTY )() \
723  BOOST_PP_IF( VF_OP_TYPE_IS_GENERIC( T ), BOOST_PP_IDENTITY( VF_OP_TYPE_TYPE( T ) sw ), BOOST_PP_EMPTY )() > \
724  inline Expr< VF_OPERATOR_NAME( O )< ELEM, VF_OP_TYPE_OBJECT(T)> > \
725  BOOST_PP_CAT( VF_OPERATOR_SYMBOL(O), VF_OP_TYPE_SUFFIX(T) )( boost::shared_ptr<ELEM> expr ) \
726  { \
727  typedef VF_OPERATOR_NAME( O )< ELEM, VF_OP_TYPE_OBJECT(T)> expr_t; \
728  return Expr< expr_t >( expr_t(*expr) ); \
729  } \
730 
731 
732 #
733 //
734 // Generate the code
735 //
736 BOOST_PP_LIST_FOR_EACH_PRODUCT( VF_ARRAY_OPERATOR, 2, ( VF_OPERATORS, VF_OPERATORS_TYPE ) )
738 }
739 }
740 
741 #endif /* __FEELPP_VF_OPERATORS_HPP */

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