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
bilinearformcontext.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 __BilinearFormContext_H
30 #define __BilinearFormContext_H 1
31 
32 namespace Feel
33 {
34 namespace vf
35 {
36 namespace detail
37 {
38 //
39 // Context
40 //
41 template<typename FE1, typename FE2, typename ElemContType>
42 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
43 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::Context( form_type& __form,
44  map_test_geometric_mapping_context_type const& _gmcTest,
45  map_trial_geometric_mapping_context_type const& _gmcTrial,
46  map_geometric_mapping_expr_context_type const & gmcExpr,
47  ExprT const& expr,
48  IM const& im )
49  :
50  M_form( __form ),
51  M_lb( __form.blockList() ),
52  M_test_dof( __form.testSpace()->dof().get() ),
53  M_trial_dof( __form.trialSpace()->dof().get() ),
54 
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  M_trial_pc( new trial_precompute_type( M_form.trialSpace()->fe(), fusion::at_key<gmc<0> >( _gmcTrial )->pc()->nodes() ) ),
59  M_trial_pc_face( precomputeTrialBasisAtPoints( im ) ),
60 
61  M_test_gmc( _gmcTest ),
62  M_trial_gmc( _gmcTrial ),
63 
64  M_test_fec( fusion::transform( _gmcTest,
65  vf::detail::FEContextInit<0,form_context_type>( __form.testSpace()->fe(),
66  *this ) ) ),
67  M_test_fec0( fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) ) ),
68  M_trial_fec( getMap( M_test_fec, fusion::transform( _gmcTrial,
69  vf::detail::FEContextInit<1,form_context_type>( __form.trialSpace()->fe(),
70  *this ) ) ) ),
71  M_trial_fec0( getMapL( M_test_fec0, fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_trial_fec ) ) ) ),
72 
73  M_rep(),
74  M_rep_2(),
75  M_eval_expr00( new eval00_expr_type( expr, gmcExpr, M_test_fec0, M_trial_fec0 ) ),
76 
77  M_eval_expr01(),
78  M_eval_expr10(),
79  M_eval_expr11(),
80 
81  M_integrator( im )
82 {
83  M_eval_expr00->init( im );
84 }
85 
86 template<typename FE1, typename FE2, typename ElemContType>
87 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
88 template<typename IM2>
89 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::Context( form_type& __form,
90  map_test_geometric_mapping_context_type const& _gmcTest,
91  map_trial_geometric_mapping_context_type const& _gmcTrial,
92  map_geometric_mapping_expr_context_type const & _gmcExpr,
93  ExprT const& expr,
94  IM const& im,
95  IM2 const& im2 )
96  :
97  M_form( __form ),
98  M_lb( __form.blockList() ),
99  M_test_dof( __form.testSpace()->dof().get() ),
100  M_trial_dof( __form.trialSpace()->dof().get() ),
101 
102  M_test_pc( new test_precompute_type( M_form.testSpace()->fe(), im2.points() ) ),
103  M_test_pc_face( precomputeTestBasisAtPoints( im2 ) ),
104  M_trial_pc( new trial_precompute_type( M_form.trialSpace()->fe(), im2.points() ) ),
105  M_trial_pc_face( precomputeTrialBasisAtPoints( im2 ) ),
106 
107  M_test_gmc( _gmcTest ),
108  M_trial_gmc( _gmcTrial ),
109 
110  M_test_fec( fusion::transform( _gmcTest, vf::detail::FEContextInit<0,form_context_type>( __form.testSpace()->fe(), *this ) ) ),
111  M_test_fec0( fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) ) ),
112  M_trial_fec( getMap( M_test_fec, fusion::transform( _gmcTrial, vf::detail::FEContextInit<1,form_context_type>( __form.trialSpace()->fe(), *this ) ) ) ),
113  M_trial_fec0( getMapL( M_test_fec0, fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_trial_fec ) ) ) ),
114  M_rep(),
115  M_rep_2(),
116  M_eval_expr00( new eval00_expr_type( expr, _gmcExpr, M_test_fec0, M_trial_fec0 ) ),
117  M_eval_expr01(),
118  M_eval_expr10(),
119  M_eval_expr11(),
120  M_integrator( im )
121 {
122  // faces
123  M_eval_expr00->init( im2 );
124 }
125 
126 template<typename FE1, typename FE2, typename ElemContType>
127 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
128 template</*typename IM2, */typename IMTest,typename IMTrial>
129 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::Context( form_type& __form,
130  map_test_geometric_mapping_context_type const& _gmcTest,
131  map_trial_geometric_mapping_context_type const& _gmcTrial,
132  map_geometric_mapping_expr_context_type const & _gmcExpr,
133  ExprT const& expr,
134  IM const& im,
135  //IM2 const& im2,
136  IMTest const& imTest, IMTrial const& imTrial )
137  :
138  M_form( __form ),
139  M_lb( __form.blockList() ),
140  M_test_dof( __form.testSpace()->dof().get() ),
141  M_trial_dof( __form.trialSpace()->dof().get() ),
142 
143  M_test_pc( new test_precompute_type( M_form.testSpace()->fe(), fusion::at_key<gmc<0> >( _gmcTest )->pc()->nodes() ) ),
144  M_test_pc_face( precomputeTestBasisAtPoints( imTest ) ),
145  M_trial_pc( new trial_precompute_type( M_form.trialSpace()->fe(), fusion::at_key<gmc<0> >( _gmcTrial )->pc()->nodes() ) ),
146  M_trial_pc_face( precomputeTrialBasisAtPoints( imTrial ) ),
147 
148  /*M_test_pc( new test_precompute_type( M_form.testSpace()->fe(), im2.points() ) ),
149  M_test_pc_face( precomputeTestBasisAtPoints( im2 ) ),
150  M_trial_pc( new trial_precompute_type( M_form.trialSpace()->fe(), im2.points() ) ),
151  M_trial_pc_face( precomputeTrialBasisAtPoints( im2 ) ),*/
152 
153  M_test_gmc( _gmcTest ),
154  M_trial_gmc( _gmcTrial ),
155 
156  M_test_fec( fusion::transform( _gmcTest, vf::detail::FEContextInit<0,form_context_type>( __form.testSpace()->fe(), *this ) ) ),
157  M_test_fec0( fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) ) ),
158  M_trial_fec( getMap( M_test_fec, fusion::transform( _gmcTrial, vf::detail::FEContextInit<1,form_context_type>( __form.trialSpace()->fe(), *this ) ) ) ),
159  M_trial_fec0( getMapL( M_test_fec0, fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_trial_fec ) ) ) ),
160  M_rep(),
161  M_rep_2(),
162  M_eval_expr00( new eval00_expr_type( expr, _gmcExpr, M_test_fec0, M_trial_fec0 ) ),
163  M_eval_expr01(),
164  M_eval_expr10(),
165  M_eval_expr11(),
166  M_integrator( im )
167 {
168  // faces
169  M_eval_expr00->init( im );
170 }
171 
172 template<typename FE1, typename FE2, typename ElemContType>
173 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
174 template<typename IM2>
175 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::Context( form_type& __form,
176  map_test_geometric_mapping_context_type const& _gmcTest,
177  map_trial_geometric_mapping_context_type const& _gmcTrial,
178  map_geometric_mapping_expr_context_type const & _gmcExpr,
179  ExprT const& expr,
180  IM const& im,
181  IM2 const& im2,
182  mpl::int_<2> )
183  :
184  M_form( __form ),
185  M_lb( __form.blockList() ),
186  M_test_dof( __form.testSpace()->dof().get() ),
187  M_trial_dof( __form.trialSpace()->dof().get() ),
188 
189  M_test_pc( new test_precompute_type( M_form.testSpace()->fe(), im2.points() ) ),
190  M_test_pc_face( precomputeTestBasisAtPoints( im2 ) ),
191  M_trial_pc( new trial_precompute_type( M_form.trialSpace()->fe(), im2.points() ) ),
192  M_trial_pc_face( precomputeTrialBasisAtPoints( im2 ) ),
193 
194  M_test_gmc( _gmcTest ),
195  M_trial_gmc( _gmcTrial ),
196 
197  M_test_fec( fusion::transform( _gmcTest, vf::detail::FEContextInit<0,form_context_type>( __form.testSpace()->fe(), *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<test_gmc1 >( fusion::at_key<test_gmc1 >( M_test_fec ) ) ),
200  M_trial_fec( fusion::transform( _gmcTrial, vf::detail::FEContextInit<1,form_context_type>( __form.trialSpace()->fe(), *this ) ) ),
201  M_trial_fec0( fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_trial_fec ) ) ),
202  M_trial_fec1( fusion::make_map<trial_gmc1 >( fusion::at_key<trial_gmc1 >( M_trial_fec ) ) ),
203  M_rep(),
204  M_rep_2(),
205  M_eval_expr00( new eval00_expr_type( expr, _gmcExpr, M_test_fec0, M_trial_fec0 ) ),
206  M_eval_expr01( new eval01_expr_type( expr, _gmcExpr, M_test_fec0, M_trial_fec1 ) ),
207  M_eval_expr10( new eval10_expr_type( expr, _gmcExpr, M_test_fec1, M_trial_fec0 ) ),
208  M_eval_expr11( new eval11_expr_type( expr, _gmcExpr, M_test_fec1, M_trial_fec1 ) ),
209  M_integrator( im )
210 {
211  FEELPP_ASSERT( fusion::at_key<gmc<0> >( M_test_fec0 ).get() != 0 ).error( "invalid test_fec" );
212  FEELPP_ASSERT( fusion::at_key<test_gmc1 >( M_test_fec1 ).get() != 0 ).error( "invalid test_fec" );
213  FEELPP_ASSERT( fusion::at_key<gmc<0> >( M_trial_fec0 ).get() != 0 ).error( "invalid trial_fec" );
214  FEELPP_ASSERT( fusion::at_key<trial_gmc1 >( M_trial_fec1 ).get() != 0 ).error( "invalid trial_fec" );
215 
216  M_eval_expr00->init( im2 );
217  M_eval_expr01->init( im2 );
218  M_eval_expr10->init( im2 );
219  M_eval_expr11->init( im2 );
220 }
221 
222 template<typename FE1, typename FE2, typename ElemContType>
223 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
224 void
225 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
226  map_trial_geometric_mapping_context_type const& _gmcTrial,
227  map_geometric_mapping_expr_context_type const& _gmcExpr )
228 {
229  update( _gmcTest, _gmcTrial, _gmcExpr, boost::is_same<map_test_fecontext_type, map_trial_fecontext_type>() );
230  // if we know that the result will be zero, don't update the integrator and return immediately
231  M_integrator.update( *fusion::at_key<gmc<0> >( _gmcExpr ) );
232 }
233 template<typename FE1, typename FE2, typename ElemContType>
234 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
235 void
236 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
237  map_trial_geometric_mapping_context_type const& _gmcTrial,
238  map_geometric_mapping_expr_context_type const& _gmcExpr,
239  mpl::bool_<false> )
240 {
241  fusion::for_each( M_test_fec, vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
242  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
243  fusion::for_each( M_trial_fec, vf::detail::FEContextUpdate<1,form_context_type>( _gmcTrial, *this ) );
244  M_trial_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_trial_fec ) );
245  M_eval_expr00->update( _gmcExpr, M_test_fec0, M_trial_fec0 );
246 }
247 template<typename FE1, typename FE2, typename ElemContType>
248 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
249 void
250 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
251  map_trial_geometric_mapping_context_type const& _gmcTrial,
252  map_geometric_mapping_expr_context_type const& _gmcExpr,
253  mpl::bool_<true> )
254 {
255  fusion::for_each( M_test_fec, vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
256  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
257  M_eval_expr00->update( _gmcExpr, M_test_fec0, M_test_fec0 );
258 }
259 template<typename FE1, typename FE2, typename ElemContType>
260 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
261 void
262 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::update( map_test_geometric_mapping_context_type const& _gmcTest,
263  map_trial_geometric_mapping_context_type const& _gmcTrial,
264  map_geometric_mapping_expr_context_type const& _gmcExpr,
265  mpl::int_<2> )
266 {
267  fusion::for_each( M_test_fec, vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
268  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
269  M_test_fec1 = fusion::make_map<test_gmc1 >( fusion::at_key<test_gmc1 >( M_test_fec ) );
270  fusion::for_each( M_trial_fec, vf::detail::FEContextUpdate<1,form_context_type>( _gmcTrial, *this ) );
271  M_trial_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_trial_fec ) );
272  M_trial_fec1 = fusion::make_map<trial_gmc1 >( fusion::at_key<trial_gmc1 >( M_trial_fec ) );
273 
274  FEELPP_ASSERT( fusion::at_key<gmc<0> >( M_test_fec0 ).get() != 0 )
275  ( 0 ).error( "invalid test_fec0" );
276  FEELPP_ASSERT( fusion::at_key<gmc<1> >( M_test_fec1 ).get() != 0 )
277  ( 1 ).error( "invalid test_fec1" );
278  FEELPP_ASSERT( fusion::at_key<gmc<0> >( M_trial_fec0 ).get() != 0 )
279  ( 0 ).error( "invalid trial_fec0" );
280  FEELPP_ASSERT( fusion::at_key<gmc<1> >( M_trial_fec1 ).get() != 0 )
281  ( 0 ).error( "invalid trial_fec1" );
282 
283  M_eval_expr00->update( _gmcExpr, M_test_fec0, M_trial_fec0 );
284  M_eval_expr01->update( _gmcExpr, M_test_fec0, M_trial_fec1 );
285  M_eval_expr10->update( _gmcExpr, M_test_fec1, M_trial_fec0 );
286  M_eval_expr11->update( _gmcExpr, M_test_fec1, M_trial_fec1 );
287 
288  M_integrator.update( *fusion::at_key<gmc<0> >( _gmcTest ) );
289 }
290 
291 template<typename FE1, typename FE2, typename ElemContType>
292 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
293 void
294 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& _gmcTest,
295  map_trial_geometric_mapping_context_type const& _gmcTrial,
296  map_geometric_mapping_expr_context_type const& _gmcExpr,
297  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad )
298 {
299  M_test_gmc = _gmcTest;
300  M_trial_gmc = _gmcTrial;
301  precomputeBasisAtPoints( fusion::at_key<gmc<0> >( _gmcTest )->xRefs(),
302  fusion::at_key<gmc<0> >( _gmcTrial )->xRefs() );
303  //updateInCaseOfInterpolate( _gmcTest, _gmcTrial, _gmcExpr, boost::is_same<map_test_fecontext_type, map_trial_fecontext_type>() );
304  updateInCaseOfInterpolate( _gmcTest, _gmcTrial, _gmcExpr, mpl::bool_<false>() ); // forcage!
305  M_integrator.update( *fusion::at_key<gmc<0> >( _gmcExpr ), indexLocalToQuad );
306 }
307 template<typename FE1, typename FE2, typename ElemContType>
308 template<typename GeomapContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
309 void
310 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& _gmcTest,
311  map_trial_geometric_mapping_context_type const& _gmcTrial,
312  map_geometric_mapping_expr_context_type const& _gmcExpr,
313  mpl::bool_<false> )
314 {
315  fusion::for_each( M_test_fec, vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
316  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
317  fusion::for_each( M_trial_fec, vf::detail::FEContextUpdate<1,form_context_type>( _gmcTrial, *this ) );
318  M_trial_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_trial_fec ) );
319  M_eval_expr00->update( _gmcExpr, M_test_fec0, M_trial_fec0 );
320 }
321 template<typename FE1, typename FE2, typename ElemContType>
322 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
323 void
324 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& _gmcTest,
325  map_trial_geometric_mapping_context_type const& _gmcTrial,
326  map_geometric_mapping_expr_context_type const& _gmcExpr,
327  mpl::bool_<true> )
328 {
329  fusion::for_each( M_test_fec, vf::detail::FEContextUpdate<0,form_context_type>( _gmcTest, *this ) );
330  M_test_fec0 = fusion::make_map<gmc<0> >( fusion::at_key<gmc<0> >( M_test_fec ) );
331  M_eval_expr00->update( _gmcExpr, M_test_fec0, M_test_fec0 );
332 }
333 
334 
335 
336 template<typename FE1, typename FE2, typename ElemContType>
337 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
338 void
340 {
341 
342  typedef test_geometric_mapping_context_type gmc_type;
343  typedef typename eval00_expr_type::shape shape;
344  static const bool cond = ( shape::M == 1 && shape::N == 1 );
345  BOOST_MPL_ASSERT_MSG( cond,
346  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_0,
347  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
348 
349 #if !defined(NDEBUG)
350  test_geometric_mapping_context_type const& _gmc = *fusion::at_key<gmc<0> >( M_test_gmc );
351  DVLOG(2) << "[BilinearForm::integrate] local assembly in element " << _gmc.id() << "\n";
352 #endif /* NDEBUG */
353 
354  if ( M_form.isPatternDefault() && boost::is_same<trial_dof_type,test_dof_type>::value &&
355  trial_dof_type::is_product )
356  {
357  M_rep = local_matrix_type::Zero();
358 
359  if ( M_form.isPatternSymmetric() )
360  {
361  for ( uint16_type c = 0; c < trial_dof_type::nComponents1; ++c )
362  for ( uint16_type j = 0; j < trial_dof_type::fe_type::nLocalDof; ++j )
363  for ( uint16_type i = 0; i <= j; ++i )
364  {
365  uint16_type testLocalDofIndex = i+c*test_dof_type::fe_type::nLocalDof;
366  uint16_type trialLocalDofIndex = j+c*trial_dof_type::fe_type::nLocalDof;
367  M_rep( testLocalDofIndex, trialLocalDofIndex ) = M_integrator( *M_eval_expr00, testLocalDofIndex, trialLocalDofIndex, 0, 0 );
368  M_rep( trialLocalDofIndex, testLocalDofIndex ) = M_rep( testLocalDofIndex, trialLocalDofIndex );
369  }
370  }
371 
372  else
373  {
374  for ( uint16_type c = 0; c < trial_dof_type::nComponents1; ++c )
375  for ( uint16_type j = 0; j < trial_dof_type::fe_type::nLocalDof; ++j )
376  for ( uint16_type i = 0; i < test_dof_type::fe_type::nLocalDof; ++i )
377  {
378  uint16_type testLocalDofIndex = i+c*test_dof_type::fe_type::nLocalDof;
379  uint16_type trialLocalDofIndex = j+c*trial_dof_type::fe_type::nLocalDof;
380  M_rep( testLocalDofIndex, trialLocalDofIndex ) = M_integrator( *M_eval_expr00, testLocalDofIndex, trialLocalDofIndex, 0, 0 );
381  }
382  }
383  }
384 
385  else
386  {
387  if ( boost::is_same<trial_dof_type,test_dof_type>::value && M_form.isPatternSymmetric() )
388  {
389  for ( uint16_type j = 0; j < trial_dof_type::nDofPerElement; ++j )
390  for ( uint16_type i = 0; i <= j; ++i )
391  {
392  M_rep( i, j ) = M_integrator( *M_eval_expr00, i, j, 0, 0 );
393  M_rep( j,i )=M_rep( i,j );
394  }
395  }
396 
397  else
398  {
399  for ( uint16_type j = 0; j < trial_dof_type::nDofPerElement; ++j )
400  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
401  {
402  M_rep( i, j ) = M_integrator( *M_eval_expr00, i, j, 0, 0 );
403  }
404  }
405  }
406 }
407 template<typename FE1, typename FE2, typename ElemContType>
408 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
409 void
411 {
412  //geometric_mapping_context_type const& _gmc = *fusion::at_key<gmc<0> >( M_gmc );
413  typedef test_geometric_mapping_context_type gmc_type;
414  typedef typename eval00_expr_type::shape shape;
415  BOOST_MPL_ASSERT_MSG( ( mpl::and_<mpl::equal_to<mpl::int_<shape::M>,mpl::int_<1> >,
416  mpl::equal_to<mpl::int_<shape::N>,mpl::int_<1> > >::value ),
417  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_0,
418  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
419 
420 
421  for ( uint16_type j = 0; j < trial_dof_type::nDofPerElement; ++j )
422  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
423  {
424  uint16_type ii = i;
425  uint16_type jj = j;
426  // test dof element 0 - trial dof element 0
427  M_rep_2( i, j ) = M_integrator( *M_eval_expr00, i, j, 0, 0 );
428 
429  ii = i;
430  jj = j + trial_dof_type::nDofPerElement;
431  // test dof element 0 - trial dof element 1
432  M_rep_2( ii,jj ) = M_integrator( *M_eval_expr01, i, j, 0, 0 );
433 
434  ii = i + test_dof_type::nDofPerElement;
435  jj = j;
436  // test dof element 1 - trial dof element 0
437  M_rep_2( ii,jj ) = M_integrator( *M_eval_expr10, i, j, 0, 0 );
438 
439  ii = i + test_dof_type::nDofPerElement;
440  jj = j + trial_dof_type::nDofPerElement;
441  // test dof element 1 - trial dof element 1
442  M_rep_2( ii,jj ) = M_integrator( *M_eval_expr11, i, j, 0, 0 );
443  }
444 }
445 template<typename FE1, typename FE2, typename ElemContType>
446 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
447 void
448 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::integrateInCaseOfInterpolate( mpl::int_<1>,
449  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
450  bool isFirstExperience )
451 {
452 
453  typedef test_geometric_mapping_context_type gmc_type;
454  typedef typename eval00_expr_type::shape shape;
455  static const bool cond = ( shape::M == 1 && shape::N == 1 );
456  BOOST_MPL_ASSERT_MSG( cond,
457  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_0,
458  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
459 
460 #if !defined(NDEBUG)
461  test_geometric_mapping_context_type const& _gmc = *fusion::at_key<gmc<0> >( M_test_gmc );
462  DVLOG(2) << "[BilinearForm::integrate] local assembly in element " << _gmc.id() << "\n";
463 #endif /* NDEBUG */
464 
465  if ( isFirstExperience )
466  for ( uint16_type j = 0; j < trial_dof_type::nDofPerElement; ++j )
467  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
468  {
469  M_rep( i, j ) = M_integrator( *M_eval_expr00, i, j, 0, 0, indexLocalToQuad );
470  }
471 
472  else
473  for ( uint16_type j = 0; j < trial_dof_type::nDofPerElement; ++j )
474  for ( uint16_type i = 0; i < test_dof_type::nDofPerElement; ++i )
475  {
476  M_rep( i, j ) += M_integrator( *M_eval_expr00, i, j, 0, 0, indexLocalToQuad );
477  }
478 }
479 template<typename FE1, typename FE2, typename ElemContType>
480 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
481 void
482 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::assemble( size_type elt_0 )
483 {
484  size_type row_start = M_lb.front().globalRowStart();
485  size_type col_start = M_lb.front().globalColumnStart();
486 
487 #if !defined(NDEBUG)
488  DVLOG(2) << "[BilinearForm::assemble] global assembly in element " << elt_0 << "\n";
489  DVLOG(2) << "[BilinearForm::assemble] row start " << row_start << "\n";
490  DVLOG(2) << "[BilinearForm::assemble] col start " << col_start << "\n";
491 #endif /* NDEBUG */
492  bool do_less = ( ( M_form.isPatternDefault() &&
493  ( M_test_dof->nComponents == M_trial_dof->nComponents ) ) &&
494  !M_form.isPatternCoupled() );
495 
496  if ( do_less )
497  {
498  for ( uint16_type c = 0; c < trial_dof_type::nComponents1; ++c )
499  {
500  M_c_rep = M_rep.block( c*test_dof_type::fe_type::nLocalDof, c*trial_dof_type::fe_type::nLocalDof,
501  test_dof_type::fe_type::nLocalDof, trial_dof_type::fe_type::nLocalDof );
502  M_c_local_rows.array() = M_test_dof->localToGlobalIndices( elt_0 ).array().segment( c*test_dof_type::fe_type::nLocalDof,
503  test_dof_type::fe_type::nLocalDof ) + row_start;
504  M_c_local_cols.array() = M_trial_dof->localToGlobalIndices( elt_0 ).array().segment( c*trial_dof_type::fe_type::nLocalDof,
505  trial_dof_type::fe_type::nLocalDof ) + col_start;
506 
507  if ( test_dof_type::is_modal || trial_dof_type::is_modal )
508  {
509  M_c_local_rowsigns = M_test_dof->localToGlobalSigns( elt_0 ).segment( c*test_dof_type::fe_type::nLocalDof,test_dof_type::fe_type::nLocalDof );
510  M_c_local_colsigns = M_trial_dof->localToGlobalSigns( elt_0 ).segment( c*trial_dof_type::fe_type::nLocalDof,trial_dof_type::fe_type::nLocalDof );
511  M_c_rep.array() *= ( M_c_local_rowsigns*M_c_local_colsigns.transpose() ).array().template cast<value_type>();
512  }
513 
514  M_form.addMatrix( M_c_local_rows.data(), M_c_local_rows.size(),
515  M_c_local_cols.data(), M_c_local_cols.size(),
516  M_c_rep.data() );
517  }
518  }
519 
520  else
521  {
522  size_type trial_eid= this->trialElementId( elt_0 );
523  //
524  DCHECK( trial_eid != invalid_size_type_value )
525  << "this case should have been taken care of earlier before the assembly process\n";
526 
527  M_local_rows.array() = M_test_dof->localToGlobalIndices( elt_0 ).array() + row_start;
528  M_local_cols.array() = M_trial_dof->localToGlobalIndices( trial_eid ).array() + col_start;
529 
530 
531  if ( test_dof_type::is_modal || trial_dof_type::is_modal )
532  {
533  M_local_rowsigns = M_test_dof->localToGlobalSigns( elt_0 );
534  M_local_colsigns = M_trial_dof->localToGlobalSigns( trial_eid );
535  M_rep.array() *= ( M_local_rowsigns*M_local_colsigns.transpose() ).array().template cast<value_type>();
536  }
537 
538  M_form.addMatrix( M_local_rows.data(), M_local_rows.size(),
539  M_local_cols.data(), M_local_cols.size(),
540  M_rep.data() );
541  }
542 }
543 
544 template<typename FE1, typename FE2, typename ElemContType>
545 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
546 void
547 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::assemble( size_type elt_0, size_type elt_1 )
548 {
549  size_type row_start = M_lb.front().globalRowStart();
550  size_type col_start = M_lb.front().globalColumnStart();
551 
552  auto local_rows_0 = M_test_dof->localToGlobalIndices( elt_0 ).array() + row_start;
553  auto local_rows_1 = M_test_dof->localToGlobalIndices( elt_1 ).array() + row_start;
554  M_local_rows_2.template head<test_dof_type::nDofPerElement>().array() = local_rows_0;
555  M_local_rows_2.template tail<test_dof_type::nDofPerElement>().array() = local_rows_1;
556 
557  auto local_cols_0 = M_trial_dof->localToGlobalIndices( elt_0 ).array() + col_start;
558  auto local_cols_1 = M_trial_dof->localToGlobalIndices( elt_1 ).array() + col_start;
559  M_local_cols_2.template head<trial_dof_type::nDofPerElement>().array() = local_cols_0;
560  M_local_cols_2.template tail<trial_dof_type::nDofPerElement>().array() = local_cols_1;
561 
562 
563  if ( test_dof_type::is_modal || trial_dof_type::is_modal )
564  {
565  auto local_rowsigns_0 = M_test_dof->localToGlobalSigns( elt_0 );
566  auto local_rowsigns_1 = M_test_dof->localToGlobalSigns( elt_1 );
567  M_local_rowsigns_2.template head<test_dof_type::nDofPerElement>() = local_rowsigns_0;
568  M_local_rowsigns_2.template tail<test_dof_type::nDofPerElement>() = local_rowsigns_1;
569 
570  auto local_colsigns_0 = M_trial_dof->localToGlobalSigns( elt_0 );
571  auto local_colsigns_1 = M_trial_dof->localToGlobalSigns( elt_1 );
572  M_local_colsigns_2.template head<trial_dof_type::nDofPerElement>() = local_colsigns_0;
573  M_local_colsigns_2.template tail<trial_dof_type::nDofPerElement>() = local_colsigns_1;
574 
575  M_rep_2.array() *= ( M_local_rowsigns_2*M_local_colsigns_2.transpose() ).array().template cast<value_type>();
576  }
577 
578  M_form.addMatrix( M_local_rows_2.data(), M_local_rows_2.size(),
579  M_local_cols_2.data(), M_local_cols_2.size(),
580  M_rep_2.data() );
581 }
582 
583 template<typename FE1, typename FE2, typename ElemContType>
584 template<typename GeomapTestContext,typename ExprT,typename IM,typename GeomapExprContext,typename GeomapTrialContext>
585 void
586 BilinearForm<FE1,FE2,ElemContType>::Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext>::assembleInCaseOfInterpolate()
587 {
588  size_type row_start = M_lb.front().globalRowStart();
589  size_type col_start = M_lb.front().globalColumnStart();
590 
591  auto eltTest = fusion::at_key<gmc<0> >( M_test_gmc )->id();
592  auto eltTrial = fusion::at_key<gmc<0> >( M_trial_gmc )->id();
593 
594  M_local_rows.array() = M_test_dof->localToGlobalIndices( eltTest ).array() + row_start;
595  M_local_cols.array() = M_trial_dof->localToGlobalIndices( eltTrial ).array() + col_start;
596 
597 #if 0
598  bool do_less = ( ( M_form.isPatternDefault() &&
599  ( M_test_dof->nComponents == M_trial_dof->nComponents ) ) &&
600  !M_form.isPatternCoupled() );
601 #endif
602 
603  if ( test_dof_type::is_modal || trial_dof_type::is_modal )
604  {
605  M_local_rowsigns = M_test_dof->localToGlobalSigns( eltTest );
606  M_local_colsigns = M_trial_dof->localToGlobalSigns( eltTrial );
607  M_rep.array() *= ( M_local_rowsigns*M_local_colsigns.transpose() ).array().template cast<value_type>();
608  }
609 
610  M_form.addMatrix( M_local_rows.data(), M_local_rows.size(),
611  M_local_cols.data(), M_local_cols.size(),
612  M_rep.data() );
613 
614 }
615 
616 
617 }
618 }
619 }
620 #endif /* __BilinearFormContext_H */

Generated on Sun Oct 20 2013 08:24:54 for Feel++ by doxygen 1.8.4