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
linearformcontext.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: 2010-04-27
7 
8  Copyright (C) 2010 Université 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 2.1 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 __LinearFormContext_H
30 #define __LinearFormContext_H 1
31 
32 
33 namespace Feel
34 {
35 namespace vf
36 {
37 namespace detail
38 {
39 //
40 // Context class for linear forms
41 //
42 template<typename SpaceType, typename VectorType, typename ElemContType>
43 template<typename GeomapContext,typename ExprT, typename IM,typename GeomapExprContext>
44 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::Context( form_type& __form,
45  map_test_geometric_mapping_context_type const& _gmcTest,
46  map_trial_geometric_mapping_context_type const & _gmcTrial,
47  map_geometric_mapping_expr_context_type const& _gmcExpr,
48  ExprT const& expr,
49  IM const& im )
50  :
51  //super(),
52  M_form( __form ),
53  M_test_dof( __form.functionSpace()->dof().get() ),
54  M_lb( __form.blockList() ),
55 
56  M_test_pc( new test_precompute_type( M_form.testSpace()->fe(), fusion::at_key<gmc<0> >( _gmcTest )->pc()->nodes() ) ),
57  M_test_pc_face( precomputeTestBasisAtPoints( im ) ),
58 
59  M_gmc( _gmcTest ),
60  M_gmc_left( fusion::at_key<gmc<0> >( _gmcTest ) ),
61  M_left_map( fusion::make_map<gmc<0> >( M_gmc_left ) ),
62  M_test_fec( fusion::transform( M_gmc,vf::detail::FEContextInit<0,form_context_type>( __form.functionSpace()->fe(), *this ) ) ),
63  M_test_fec0( fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) ) ),
64  M_rep(),
65  M_rep_2(),
66  M_eval0_expr( new eval0_expr_type( expr, _gmcExpr, M_test_fec0 ) ),
67  M_eval1_expr(),
68  M_integrator( im )
69 {
70  M_eval0_expr->init( im );
71 }
72 
73 template<typename SpaceType, typename VectorType, typename ElemContType>
74 template<typename GeomapContext,typename ExprT, typename IM,typename GeomapExprContext>
75 template<typename IM2>
76 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::Context( form_type& __form,
77  map_test_geometric_mapping_context_type const& _gmcTest,
78  map_trial_geometric_mapping_context_type const & _gmcTrial,
79  map_geometric_mapping_expr_context_type const& _gmcExpr,
80  ExprT const& expr,
81  IM const& im,
82  IM2 const& im2 )
83  :
84  //super(),
85  M_form( __form ),
86  M_test_dof( __form.functionSpace()->dof().get() ),
87  M_lb( __form.blockList() ),
88 
89  M_test_pc( new test_precompute_type( M_form.testSpace()->fe(), im2.points() ) ),
90  M_test_pc_face( precomputeTestBasisAtPoints( im2 ) ),
91 
92  M_gmc( _gmcTest ),
93  M_gmc_left( fusion::at_key<gmc<0> >( _gmcTest ) ),
94  M_left_map( fusion::make_map<gmc<0> >( M_gmc_left ) ),
95  M_test_fec( fusion::transform( M_gmc,vf::detail::FEContextInit<0,form_context_type>( __form.functionSpace()->fe(), *this ) ) ),
96  M_test_fec0( fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) ) ),
97  M_rep(),
98  M_rep_2(),
99  M_eval0_expr( new eval0_expr_type( expr, _gmcExpr, M_test_fec0 ) ),
100  M_eval1_expr(),
101  M_integrator( im )
102 {
103  M_eval0_expr->init( im2 );
104 }
105 template<typename SpaceType, typename VectorType, typename ElemContType>
106 template<typename GeomapContext,typename ExprT, typename IM,typename GeomapExprContext>
107 template<typename IM2>
108 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::Context( form_type& __form,
109  map_test_geometric_mapping_context_type const& _gmcTest,
110  map_trial_geometric_mapping_context_type const & _gmcTrial,
111  map_geometric_mapping_expr_context_type const& _gmcExpr,
112  ExprT const& expr,
113  IM const& im,
114  IM2 const& im2,
115  mpl::int_<2> )
116  :
117  //super(),
118  M_form( __form ),
119  M_test_dof( __form.functionSpace()->dof().get() ),
120  M_lb( __form.blockList() ),
121 
122  M_test_pc( new test_precompute_type( M_form.testSpace()->fe(), im2.points() ) ),
123  M_test_pc_face( precomputeTestBasisAtPoints( im2 ) ),
124 
125  M_gmc( _gmcTest ),
126  M_gmc_left( fusion::at_key<gmc<0> >( _gmcTest ) ),
127  M_gmc_right( fusion::at_key<gmc1 >( _gmcTest ) ),
128  M_left_map( fusion::make_map<gmc<0> >( M_gmc_left ) ),
129  M_right_map( fusion::make_map<gmc<0> >( M_gmc_right ) ),
130  M_test_fec( fusion::transform( M_gmc,vf::detail::FEContextInit<0,form_context_type>( __form.functionSpace()->fe(), *this ) ) ),
131  M_test_fec0( fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) ) ),
132  M_test_fec1( fusion::make_pair<gmc1 >( fusion::at_key<gmc1 >( M_test_fec ) ) ),
133  M_rep(),
134  M_rep_2(),
135  M_eval0_expr( new eval0_expr_type( expr, _gmcExpr, M_test_fec0 ) ),
136  M_eval1_expr( new eval1_expr_type( expr, _gmcExpr, M_test_fec1 ) ),
137  M_integrator( im )
138 
139 {
140  M_eval0_expr->init( im2 );
141  M_eval1_expr->init( im2 );
142 }
143 template<typename SpaceType, typename VectorType, typename ElemContType>
144 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
145 void
146 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
147  map_trial_geometric_mapping_context_type const & _gmcTrial,
148  map_geometric_mapping_expr_context_type const& _gmcExpr )
149 {
150  M_gmc = _gmcTest;
151  M_gmc_left = fusion::at_key<gmc<0> >( _gmcTest );
152  M_left_map = fusion::make_map<gmc<0> >( M_gmc_left );
153  fusion::for_each( M_test_fec,vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
154  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
155  M_eval0_expr->update( _gmcExpr, M_test_fec0 );
156 
157  M_integrator.update( *fusion::at_key<gmc<0> >( _gmcTest ) );
158 }
159 template<typename SpaceType, typename VectorType, typename ElemContType>
160 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
161 void
162 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& _gmcTest,
163  map_trial_geometric_mapping_context_type const & _gmcTrial,
164  map_geometric_mapping_expr_context_type const& _gmcExpr,
165  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad )
166 {
167  M_gmc = _gmcTest;
168  M_gmc_left = fusion::at_key<gmc<0> >( _gmcTest );
169  M_left_map = fusion::make_map<gmc<0> >( M_gmc_left );
170  precomputeBasisAtPoints( fusion::at_key<gmc<0> >( _gmcTest )->xRefs() );
171  fusion::for_each( M_test_fec,vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
172  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
173  M_eval0_expr->update( _gmcExpr, M_test_fec0 );
174 
175  M_integrator.update( *fusion::at_key<gmc<0> >( _gmcExpr ),indexLocalToQuad );
176 }
177 
178 template<typename SpaceType, typename VectorType, typename ElemContType>
179 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
180 void
181 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
182  map_trial_geometric_mapping_context_type const & _gmcTrial,
183  map_geometric_mapping_expr_context_type const& _gmcExpr,
184  mpl::int_<2> )
185 {
186  // typedef mpl::int_<fusion::result_of::template size<map_geometric_mapping_context_type>::type::value> map_size;
187  /*BOOST_MPL_ASSERT_MSG( (mpl::equal_to<mpl::int_<map_size::value>,mpl::int_<2> >::value),
188  INVALID_GEOMAP,
189  (map_size,map_geometric_mapping_context_type ));*/
190  //M_gmc = _gmc;
191 #if 0
192  M_gmc_left = fusion::at_key<gmc<0> >( _gmcTest );
193  M_gmc_right = fusion::at_key<gmc1 >( _gmcTest );
194  M_left_map = fusion::make_map<gmc<0> >( M_gmc_left );
195  M_right_map = fusion::make_map<gmc<0> >( M_gmc_right );
196 #endif
197  fusion::for_each( M_test_fec,vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
198  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
199  M_test_fec1 = fusion::make_map<gmc1 >( fusion::at_key<gmc1 >( M_test_fec ) );
200  M_eval0_expr->update( _gmcExpr, M_test_fec0 );
201  M_eval1_expr->update( _gmcExpr, M_test_fec1 );
202 
203  M_integrator.update( *fusion::at_key<gmc<0> >( _gmcTest ) );
204 }
205 template<typename SpaceType, typename VectorType, typename ElemContType>
206 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
207 void
208 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
209  map_trial_geometric_mapping_context_type const & _gmcTrial,
210  map_geometric_mapping_expr_context_type const& _gmcExpr,
211  IM const& im )
212 {
213  M_integrator = im;
214  this->update( _gmcTest, _gmcTrial, _gmcExpr );
215 }
216 template<typename SpaceType, typename VectorType, typename ElemContType>
217 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
218 void
219 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
220  map_trial_geometric_mapping_context_type const & _gmcTrial,
221  map_geometric_mapping_expr_context_type const& _gmcExpr,
222  IM const& im, mpl::int_<2> )
223 {
224  M_integrator = im;
225  this->update( _gmcTest, _gmcTrial, _gmcExpr, mpl::int_<2>() );
226 }
227 
228 
229 template<typename SpaceType, typename VectorType, typename ElemContType>
230 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
231 void
233 {
234  typedef typename eval0_expr_type::shape shape;
235  BOOST_MPL_ASSERT_MSG( ( shape::M == 1 && shape::N == 1 ),
236  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_0,
237  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
238 
239  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
240  {
241  M_rep( i ) = M_integrator( *M_eval0_expr, i, 0, 0 );
242  }
243 }
244 template<typename SpaceType, typename VectorType, typename ElemContType>
245 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
246 void
248 {
249  typedef mpl::int_<fusion::result_of::template size<map_test_geometric_mapping_context_type>::type::value> map_size;
250  BOOST_MPL_ASSERT_MSG( map_size::value == 2, INVALID_GEOMAP, ( map_size,map_test_geometric_mapping_context_type ) );
251 
252  typedef typename eval0_expr_type::shape shape;
253  BOOST_MPL_ASSERT_MSG( ( shape::M == 1 && shape::N == 1 ),
254  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_0,
255  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
256 
257  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
258  {
259  uint16_type ii = i;
260  // test dof element 0
261  M_rep_2( ii ) = M_integrator( *M_eval0_expr, i, 0, 0 );
262 
263  ii = i + test_dof_type::nDofPerElement;
264  // test dof element 1
265  M_rep_2( ii ) = M_integrator( *M_eval1_expr, i, 0, 0 );
266  }
267 }
268 template<typename SpaceType, typename VectorType, typename ElemContType>
269 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
270 void
271 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::integrateInCaseOfInterpolate( mpl::int_<1>,
272  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
273  bool isFirstExperience )
274 {
275  typedef typename eval0_expr_type::shape shape;
276  BOOST_MPL_ASSERT_MSG( ( shape::M == 1 && shape::N == 1 ),
277  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_0,
278  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
279 
280  if ( isFirstExperience )
281  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
282  {
283  M_rep( i ) = M_integrator( *M_eval0_expr, i, 0, 0, indexLocalToQuad );
284  }
285 
286  else
287  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
288  {
289  M_rep( i ) += M_integrator( *M_eval0_expr, i, 0, 0, indexLocalToQuad );
290  }
291 
292 
293 }
294 template<typename SpaceType, typename VectorType, typename ElemContType>
295 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
296 void
297 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::assemble( size_type elt_0 )
298 {
299  size_type row_start = M_lb.front().globalRowStart();
300  M_local_rows = M_test_dof->localToGlobalIndices( elt_0 ).array() + row_start;
301 
302  if ( test_dof_type::is_modal )
303  {
304  M_local_rowsigns = M_test_dof->localToGlobalSigns( elt_0 );
305  M_rep.array() *= M_local_rowsigns.array().template cast<value_type>();
306  }
307 
308  M_form.addVector( M_local_rows.data(), M_local_rows.size(),
309  M_rep.data() );
310 }
311 template<typename SpaceType, typename VectorType, typename ElemContType>
312 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext>
313 void
314 LinearForm<SpaceType, VectorType, ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext>::assemble( size_type elt_0, size_type elt_1 )
315 {
316  size_type row_start = M_lb.front().globalRowStart();
317  auto local_rows_0 = M_test_dof->localToGlobalIndices( elt_0 ).array() + row_start;
318  auto local_rows_1 = M_test_dof->localToGlobalIndices( elt_1 ).array() + row_start;
319  M_local_rows_2.template head<test_dof_type::nDofPerElement>() = local_rows_0;
320  M_local_rows_2.template tail<test_dof_type::nDofPerElement>() = local_rows_1;
321 
322  if ( test_dof_type::is_modal )
323  {
324  auto local_rowsigns_0 = M_test_dof->localToGlobalSigns( elt_0 );
325  auto local_rowsigns_1 = M_test_dof->localToGlobalSigns( elt_1 );
326  M_local_rowsigns_2.template head<test_dof_type::nDofPerElement>() = local_rowsigns_0;
327  M_local_rowsigns_2.template tail<test_dof_type::nDofPerElement>() = local_rowsigns_1;
328 
329  M_rep_2.array() *= M_local_rowsigns_2.array().template cast<value_type>();
330  }
331 
332  M_form.addVector( M_local_rows_2.data(), M_local_rows_2.size(),
333  M_rep_2.data() );
334 }
335 
336 }
337 }
338 }
339 #endif /* __LinearFormContext_H */

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