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
integratoron.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: 2005-03-15
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2006-2011 Universite Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __INTEGRATORON_HPP
31 #define __INTEGRATORON_HPP 1
32 
33 #include <boost/timer.hpp>
34 #include <boost/foreach.hpp>
35 
36 #include <feel/feelalg/enums.hpp>
37 
38 namespace Feel
39 {
40 namespace vf
41 {
43 template<typename T>
44 struct access_value
45 {
46 };
47 #if defined( FEELPP_HAS_QD_REAL )
48 template<>
49 struct access_value<dd_real>
50 {
51  access_value( dd_real val, int /*n*/ )
52  {
53  v = val;
54  }
55  dd_real operator()() const
56  {
57  return v;
58  }
59  dd_real v;
60 };
61 template<>
62 struct access_value<qd_real>
63 {
64  access_value( qd_real val, int /*n*/ )
65  {
66  v = val;
67  }
68  qd_real operator()() const
69  {
70  return v;
71  }
72  qd_real v;
73 };
74 #endif /*FEELPP_HAS_QD_REAL*/
75 
76 #if defined(FEELPP_HAS_MPFR)
77 template<>
78 struct access_value<mp_type>
79 {
80  access_value( mp_type val, int /*n*/ )
81  {
82  v = val;
83  }
84  mp_type operator()() const
85  {
86  return v;
87  }
88  mp_type v;
89 };
90 #endif /* FEELPP_HAS_MPFR */
91 
92 template<>
93 struct access_value<double>
94 {
95  access_value( double val, int /*n*/ )
96  {
97  v = val;
98  }
99  double operator()() const
100  {
101  return v;
102  }
103  double v;
104 };
105 template<>
106 struct access_value<int>
107 {
108  access_value( int val, int /*n*/ )
109  {
110  v = val;
111  }
112  double operator()() const
113  {
114  return v;
115  }
116  double v;
117 };
118 template<>
119 struct access_value<node_type>
120 {
121  access_value( node_type vec, int n )
122  {
123  v = vec[n];
124  }
125  double operator()() const
126  {
127  return v;
128  }
129  double v;
130 };
140 template<typename ElementRange, typename Elem, typename RhsElem, typename OnExpr >
141 class IntegratorOnExpr
142 {
143 public:
144 
145 
149  static const size_type context = OnExpr::context|vm::POINT;
150  static const size_type is_terminal = false;
151 
152  static const uint16_type imorder = OnExpr::imorder;
153  static const bool imIsPoly = OnExpr::imIsPoly;
154 
155  typedef typename boost::tuples::template element<1, ElementRange>::type element_iterator;
156 
157  typedef Elem element_type;
158  typedef RhsElem rhs_element_type;
159  typedef typename element_type::value_type value_type;
160  typedef typename element_type::return_type return_type;
161  typedef boost::function<return_type ( node_type const& )> bc_type;
162  typedef OnExpr expression_type;
163 
164  template<typename Func>
165  struct HasTestFunction
166  {
167  static const bool result = true;
168  };
169 
170  template<typename Func>
171  struct HasTrialFunction
172  {
173  static const bool result = boost::is_same<Func,typename element_type::functionspace_type::basis_type>::value;
174  };
175 
176  static const uint16_type nComponents = element_type::nComponents;
177 
179 
183 
184  IntegratorOnExpr( ElementRange const& __elts,
185  element_type const& __u,
186  rhs_element_type const& __rhs,
187  expression_type const& __expr,
188  size_type __on )
189  :
190  M_eltbegin( __elts.template get<1>() ),
191  M_eltend( __elts.template get<2>() ),
192  M_u( __u ),
193  M_rhs( __rhs ),
194  M_expr( __expr ),
195  M_on_strategy( __on )
196  {
197  }
198  IntegratorOnExpr( IntegratorOnExpr const& ioe )
199  :
200  M_eltbegin( ioe.M_eltbegin ),
201  M_eltend( ioe.M_eltend ),
202  M_u( ioe.M_u ),
203  M_rhs( ioe.M_rhs ),
204  M_expr( ioe.M_expr ),
205  M_on_strategy( ioe.M_on_strategy )
206  {
207  }
208 
209  ~IntegratorOnExpr() {}
210 
212 
216 
217 
222  element_iterator beginElement() const
223  {
224  return M_eltbegin;
225  }
226 
231  element_iterator endElement() const
232  {
233  return M_eltend;
234  }
235 
236 
238 
241 
246  template<typename Elem1, typename Elem2, typename FormType>
247  void assemble( boost::shared_ptr<Elem1> const& __u,
248  boost::shared_ptr<Elem2> const& __v,
249  FormType& __f ) const
250  {
251  typedef typename Elem::functionspace_type functionspace_type;
252  DVLOG(2) << "[IntegratorOn::assemble()] is_same: "
253  << mpl::bool_<boost::is_same<functionspace_type,Elem1>::value>::value << "\n";
254  assemble( __u, __v, __f, mpl::bool_<boost::is_same<functionspace_type,Elem1>::value>() );
255  }
257 private:
258 
259  template<typename Elem1, typename Elem2, typename FormType>
260  void assemble( boost::shared_ptr<Elem1> const& /*__u*/,
261  boost::shared_ptr<Elem2> const& /*__v*/,
262  FormType& /*__f*/, mpl::bool_<false> ) const {}
263 
264  template<typename Elem1, typename Elem2, typename FormType>
265  void assemble( boost::shared_ptr<Elem1> const& __u,
266  boost::shared_ptr<Elem2> const& __v,
267  FormType& __f, mpl::bool_<true> ) const;
268 
269 private:
270 
271  element_iterator M_eltbegin;
272  element_iterator M_eltend;
273 
274  element_type const& M_u;
275  mutable rhs_element_type M_rhs;
276  expression_type M_expr;
277  Context M_on_strategy;
278 };
279 
280 template<typename ElementRange, typename Elem, typename RhsElem, typename OnExpr>
281 template<typename Elem1, typename Elem2, typename FormType>
282 void
283 IntegratorOnExpr<ElementRange, Elem, RhsElem, OnExpr>::assemble( boost::shared_ptr<Elem1> const& /*__u*/,
284  boost::shared_ptr<Elem2> const& /*__v*/,
285  FormType& __form,
286  mpl::bool_<true> ) const
287 {
288 #if 0
289 
290  if ( !boost::is_same<Elem1, typename Elem::functionspace_type>::value ||
291  !boost::is_same<Elem2, typename Elem::functionspace_type>::value )
292  return;
293 
294 #endif
295  DVLOG(2) << "call on::assemble() " << "\n";
296  //
297  // a few typedefs
298  //
299 
300  // mesh element
301  typedef typename element_type::functionspace_type::mesh_type::element_type geoelement_type;
302  typedef typename geoelement_type::face_type face_type;
303 
304  // geometric mapping context
305  typedef typename element_type::functionspace_type::mesh_type::gm_type gm_type;
306  typedef boost::shared_ptr<gm_type> gm_ptrtype;
307  typedef typename gm_type::template Context<context, geoelement_type> gmc_type;
308  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
309  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
310 
311 
312  // dof
313  typedef typename element_type::functionspace_type::dof_type dof_type;
314 
315  // basis
316  typedef typename element_type::functionspace_type::fe_type fe_type;
317  typedef typename fe_type::template Context< context, fe_type, gm_type, geoelement_type> fecontext_type;
318  typedef boost::shared_ptr<fecontext_type> fecontext_ptrtype;
319  //typedef fusion::map<fusion::pair<vf::detail::gmc<0>, fecontext_ptrtype> > map_gmc_type;
320 
321  // expression
322  //typedef typename expression_type::template tensor<map_gmc_type,fecontext_type> t_expr_type;
323  typedef typename expression_type::template tensor<map_gmc_type> t_expr_type;
324  typedef typename t_expr_type::shape shape;
325 
326  // make sure that the form is close, ie the associated matrix is assembled
327  __form.matrix().close();
328  // make sure that the right hand side is closed, ie the associated vector is assembled
329  M_rhs->close();
330 
331  //
332  // start
333  //
334  DVLOG(2) << "assembling Dirichlet conditions\n";
335  boost::timer __timer;
336 
337  std::vector<int> dofs;
338  std::vector<value_type> values;
339  element_iterator __face_it = this->beginElement();
340  element_iterator __face_en = this->endElement();
341  if ( __face_it != __face_en )
342  {
343  // get the first face properly connected
344  for( ; __face_it != __face_en; ++__face_it )
345  if ( __face_it->isConnectedTo0() )
346  break;
347 
348 
349  dof_type const* __dof = M_u.functionSpace()->dof().get();
350 
351  fe_type const* __fe = M_u.functionSpace()->fe().get();
352 
353  gm_ptrtype __gm( new gm_type );
354 
355 
356  //
357  // Precompute some data in the reference element for
358  // geometric mapping and reference finite element
359  //
360  typedef typename geoelement_type::permutation_type permutation_type;
361  typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
362  typedef typename gm_type::precompute_type geopc_type;
363  DVLOG(2) << "[integratoron] numTopologicalFaces = " << geoelement_type::numTopologicalFaces << "\n";
364  std::vector<std::map<permutation_type, geopc_ptrtype> > __geopc( geoelement_type::numTopologicalFaces );
365 
366  for ( uint16_type __f = 0; __f < geoelement_type::numTopologicalFaces; ++__f )
367  {
368  for ( permutation_type __p( permutation_type::IDENTITY );
369  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
370  {
371  __geopc[__f][__p] = geopc_ptrtype( new geopc_type( __gm, __fe->points( __f ) ) );
372  //DVLOG(2) << "[geopc] FACE_ID = " << __f << " ref pts=" << __fe->dual().points( __f ) << "\n";
373  FEELPP_ASSERT( __geopc[__f][__p]->nPoints() ).error( "invalid number of points" );
374  }
375  }
376 
377  uint16_type __face_id = __face_it->pos_first();
378  gmc_ptrtype __c( new gmc_type( __gm, __face_it->element( 0 ), __geopc, __face_id ) );
379 
380  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
381  //t_expr_type expr( M_expr, mapgmc );
382 
383 
384  DVLOG(2) << "face_type::numVertices = " << face_type::numVertices << ", fe_type::nDofPerVertex = " << fe_type::nDofPerVertex << "\n"
385  << "face_type::numEdges = " << face_type::numEdges << ", fe_type::nDofPerEdge = " << fe_type::nDofPerEdge << "\n"
386  << "face_type::numFaces = " << face_type::numFaces << ", fe_type::nDofPerFace = " << fe_type::nDofPerFace << "\n";
387 
389 
390  if ( !fe_type::is_modal )
391  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
392  face_type::numEdges * fe_type::nDofPerEdge +
393  face_type::numFaces * fe_type::nDofPerFace );
394 
395  else
396  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
397 
398  DVLOG(2) << "nbFaceDof = " << nbFaceDof << "\n";
399  //const size_type nbFaceDof = __fe->boundaryFE()->points().size2();
400 
401  for ( ;
402  __face_it != this->endElement();
403  ++__face_it )
404  {
405  if ( !__face_it->isConnectedTo0() )
406  {
407  LOG( WARNING ) << "face not connected" << *__face_it;
408 
409  continue;
410  }
411  // do not process the face if it is a ghost face: belonging to two
412  // processes and being in a process id greater than the one
413  // corresponding face
414  if ( __face_it->isGhostFace() )
415  {
416  LOG(WARNING) << "face id : " << __face_it->id() << " is a ghost face";
417  continue;
418  }
419 
420  DVLOG(2) << "FACE_ID = " << __face_it->id()
421  << " element id= " << __face_it->ad_first()
422  << " pos in elt= " << __face_it->pos_first()
423  << " marker: " << __face_it->marker() << "\n";
424  DVLOG(2) << "FACE_ID = " << __face_it->id() << " face pts=" << __face_it->G() << "\n";
425 
426  uint16_type __face_id = __face_it->pos_first();
427  __c->update( __face_it->element( 0 ), __face_id );
428 
429  DVLOG(2) << "FACE_ID = " << __face_it->id() << " ref pts=" << __c->xRefs() << "\n";
430  DVLOG(2) << "FACE_ID = " << __face_it->id() << " real pts=" << __c->xReal() << "\n";
431 
432  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
433 
434  t_expr_type expr( M_expr, mapgmc );
435  expr.update( mapgmc );
436 
437  std::pair<size_type,size_type> range_dof( std::make_pair( M_u.start(),
438  M_u.functionSpace()->nDof() ) );
439  DVLOG(2) << "[integratoron] dof start = " << range_dof.first << "\n";
440  DVLOG(2) << "[integratoron] dof range = " << range_dof.second << "\n";
441 
442  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
443  for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
444  {
445  for ( uint16_type l = 0; l < nbFaceDof; ++l )
446  {
447  DVLOG(2) << "[integratoronexpr] local dof=" << l
448  << " |comp1=" << c1 << " comp 2= " << c2 << " | pt = " << __c->xReal( l ) << "\n";
449  typename expression_type::value_type __value = expr.evalq( c1, c2, l );
450  DVLOG(2) << "[integratoronexpr] value=" << __value << "\n";
451 
452  // global Dof
453  size_type thedof = M_u.start() +
454  boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
455 
456  //size_type thedof_nproc = __dof->dofNProc( thedof );
457  if ( std::find( dofs.begin(),
458  dofs.end(),
459  thedof ) != dofs.end() )
460  continue;
461 
462  if ( M_on_strategy.test( ON_ELIMINATION|ON_ELIMINATION_SYMMETRIC ) )
463  {
464  DVLOG(2) << "Eliminating row " << thedof << " using value : " << __value << "\n";
465 
466  // this can be quite expensive depending on the
467  // matrix storage format.
468  //__form.diagonalize( thedof, range_dof, M_rhs, __value, thedof_nproc );
469 
470  // only the real dof ( not the ghosts )
471  //if ( __form.testSpace()->mapOn().dofGlobalClusterIsOnProc( __form.testSpace()->mapOn().mapGlobalProcessToGlobalCluster( thedof ) ) )
472  {
473  dofs.push_back( thedof );
474  values.push_back( __value );
475  }
476 
477  //M_rhs.set( thedof, __value );
478  }
479 
480  else if ( M_on_strategy.test( ON_PENALISATION ) &&
481  !M_on_strategy.test( ON_ELIMINATION | ON_ELIMINATION_SYMMETRIC ) )
482  {
483  __form.set( thedof, thedof, 1.0*1e30 );
484  M_rhs->set( thedof, __value*1e30 );
485  }
486  } // loop on space components
487 
488  } // loop on face dof
489  }
490 
491  } // __face_it != __face_en
492 
493 
494  if ( __form.rowStartInMatrix()!=0)
495  {
496  auto const thedofshift = __form.rowStartInMatrix();
497  for (auto itd=dofs.begin(),end=dofs.end() ; itd!=end ; ++itd)
498  *itd+=thedofshift;
499  }
500  auto x = M_rhs->clone();
501  //CHECK( dofs.size() > 0 ) << "Invalid number of Dirichlet dof, should be > 0 ";
502  CHECK( values.size() == dofs.size() ) << "Invalid dofs/values size: " << dofs.size() << "/" << values.size();
503  //x->zero();
504  x->addVector( dofs.data(), dofs.size(), values.data() );
505  //values->zero();
506 
507  __form.zeroRows( dofs, *x, *M_rhs, M_on_strategy );
508  x.reset();
509 }
510 
511 
512 namespace detail
513 {
514 template<typename T >
515 struct v_ptr1
516 {
517  typedef T type;
518 };
519 template<typename T >
520 struct v_ptr2
521 {
522  typedef typename T::vector_ptrtype type;
523 };
524 
525 template<typename Args>
526 struct integratoron_type
527 {
528  typedef typename clean_type<Args,tag::range>::type _range_type;
529  typedef typename clean_type<Args,tag::rhs>::type _rhs_type;
530  typedef typename clean_type<Args,tag::element>::type _element_type;
531  typedef typename clean_type<Args,tag::expr>::type _expr_type;
532 #if 1
533  typedef typename mpl::if_<Feel::detail::is_vector_ptr<_rhs_type>,
534  mpl::identity<v_ptr1<_rhs_type> >,
535  mpl::identity<v_ptr2<_rhs_type> > >::type::type::type the_rhs_type;
536 #else
537  typedef _rhs_type the_rhs_type;
538 #endif
539 typedef IntegratorOnExpr<_range_type, _element_type, the_rhs_type,
540  typename mpl::if_<boost::is_arithmetic<_expr_type>,
541  mpl::identity<Expr<Cst<_expr_type> > >,
542  mpl::identity<_expr_type> >::type::type> type;
543  typedef Expr<type> expr_type;
544 };
545 
546 template<typename V>
547 typename V::vector_ptrtype
548 getRhsVector( V const& v, mpl::false_ )
549 {
550  return v.vectorPtr();
551 }
552 
553 template<typename V>
554 V
555 getRhsVector( V const& v, mpl::true_ )
556 {
557  return v;
558 }
559 
560 template<typename V>
561 typename mpl::if_<Feel::detail::is_vector_ptr<V>,
562  mpl::identity<v_ptr1<V> >,
563  mpl::identity<v_ptr2<V> > >::type::type::type
564 getRhsVector( V const& v )
565 {
566  return getRhsVector( v, Feel::detail::is_vector_ptr<V>() );
567 }
568 
569 
570 }
581 BOOST_PARAMETER_FUNCTION(
582  ( typename vf::detail::integratoron_type<Args>::expr_type ), // return type
583  on, // 2. function name
584 
585  tag, // 3. namespace of tag types
586 
587  ( required
588  ( range, * )
589  ( element, * )
590  ( rhs, * )
591  ( expr, * )
592  ) // 4. one required parameter, and
593 
594  ( optional
595  ( prefix, ( std::string ), "" )
596  ( type, ( int ), option(_prefix=prefix,_name="on.type").template as<int>() )
597  ( verbose, ( bool ), option(_prefix=prefix,_name="on.verbose").template as<bool>() )
598  )
599  )
600 {
601  typename vf::detail::integratoron_type<Args>::type ion( range,
602  element,
603  Feel::vf::detail::getRhsVector(rhs),
604  expr,
605  on_context_type(type) );
606  if ( verbose )
607  {
608  LOG(INFO) << "Dirichlet condition over : "<< nelements(range) << " faces";
609  switch( type )
610  {
611  case ON_ELIMINATION:
612  LOG(INFO) << "treatment of Dirichlet condition: " << type << " (elimination, unsymmetric)";
613  break;
615  LOG(INFO) << "treatment of Dirichlet condition: " << type << " (elimination and keep diagonal, unsymmetric)";
616  break;
618  LOG(INFO) << "treatment of Dirichlet condition: " << type << " (elimination, symmetric, more expensive than unsymmetric treatment)";
619  break;
621  LOG(INFO) << "treatment of Dirichlet condition: " << type << " (elimination and keep diagonal, symmetric, more expensive than unsymmetric treatment)";
622  break;
623  case ON_PENALISATION:
624  LOG(INFO) << "treatment of Dirichlet condition: " << type << " (penalisation, symmetric, very big value on diagonal)";
625  break;
626  default:
627  break;
628  }
629  }
630  //typename vf::detail::integratoron_type<Args>::type ion( range, element, rhs, expr, type );
631  return typename vf::detail::integratoron_type<Args>::expr_type( ion );
632 }
633 
634 
635 } // vf
636 } // feel
637 #endif

Generated on Sun Oct 20 2013 08:25:00 for Feel++ by doxygen 1.8.4