0.95.0-final
Finite Element Embedded Library and Language in C++
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
34
#include <
feel/feelpoly/quadmapped.hpp
>
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
1.8.4