30 #ifndef __INTEGRATORON_HPP
31 #define __INTEGRATORON_HPP 1
33 #include <boost/timer.hpp>
34 #include <boost/foreach.hpp>
36 #include <feel/feelalg/enums.hpp>
47 #if defined( FEELPP_HAS_QD_REAL )
49 struct access_value<dd_real>
51 access_value( dd_real
val,
int )
55 dd_real operator()()
const
62 struct access_value<qd_real>
64 access_value( qd_real
val,
int )
68 qd_real operator()()
const
76 #if defined(FEELPP_HAS_MPFR)
78 struct access_value<mp_type>
80 access_value( mp_type
val,
int )
84 mp_type operator()()
const
93 struct access_value<double>
95 access_value(
double val,
int )
99 double operator()()
const
106 struct access_value<int>
108 access_value(
int val,
int )
112 double operator()()
const
119 struct access_value<node_type>
121 access_value( node_type vec,
int n )
125 double operator()()
const
140 template<
typename ElementRange,
typename Elem,
typename RhsElem,
typename OnExpr >
141 class IntegratorOnExpr
149 static const size_type context = OnExpr::context|vm::POINT;
150 static const size_type is_terminal =
false;
152 static const uint16_type imorder = OnExpr::imorder;
153 static const bool imIsPoly = OnExpr::imIsPoly;
155 typedef typename boost::tuples::template element<1, ElementRange>::type element_iterator;
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;
164 template<
typename Func>
165 struct HasTestFunction
167 static const bool result =
true;
170 template<
typename Func>
171 struct HasTrialFunction
173 static const bool result = boost::is_same<Func,typename element_type::functionspace_type::basis_type>::value;
176 static const uint16_type nComponents = element_type::nComponents;
184 IntegratorOnExpr( ElementRange
const& __elts,
185 element_type
const& __u,
186 rhs_element_type
const& __rhs,
187 expression_type
const& __expr,
190 M_eltbegin( __elts.template get<1>() ),
191 M_eltend( __elts.template get<2>() ),
195 M_on_strategy( __on )
198 IntegratorOnExpr( IntegratorOnExpr
const& ioe )
200 M_eltbegin( ioe.M_eltbegin ),
201 M_eltend( ioe.M_eltend ),
204 M_expr( ioe.M_expr ),
205 M_on_strategy( ioe.M_on_strategy )
209 ~IntegratorOnExpr() {}
222 element_iterator beginElement()
const
231 element_iterator endElement()
const
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
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>() );
259 template<
typename Elem1,
typename Elem2,
typename FormType>
260 void assemble( boost::shared_ptr<Elem1>
const& ,
261 boost::shared_ptr<Elem2>
const& ,
262 FormType& , mpl::bool_<false> )
const {}
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;
271 element_iterator M_eltbegin;
272 element_iterator M_eltend;
274 element_type
const& M_u;
275 mutable rhs_element_type M_rhs;
276 expression_type M_expr;
277 Context M_on_strategy;
280 template<
typename ElementRange,
typename Elem,
typename RhsElem,
typename OnExpr>
281 template<
typename Elem1,
typename Elem2,
typename FormType>
283 IntegratorOnExpr<ElementRange, Elem, RhsElem, OnExpr>::assemble( boost::shared_ptr<Elem1>
const& ,
284 boost::shared_ptr<Elem2>
const& ,
286 mpl::bool_<true> )
const
290 if ( !boost::is_same<Elem1, typename Elem::functionspace_type>::value ||
291 !boost::is_same<Elem2, typename Elem::functionspace_type>::value )
295 DVLOG(2) <<
"call on::assemble() " <<
"\n";
301 typedef typename element_type::functionspace_type::mesh_type::element_type geoelement_type;
302 typedef typename geoelement_type::face_type face_type;
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;
313 typedef typename element_type::functionspace_type::dof_type dof_type;
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;
323 typedef typename expression_type::template tensor<map_gmc_type> t_expr_type;
324 typedef typename t_expr_type::shape shape;
327 __form.matrix().close();
334 DVLOG(2) <<
"assembling Dirichlet conditions\n";
335 boost::timer __timer;
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 )
344 for( ; __face_it != __face_en; ++__face_it )
345 if ( __face_it->isConnectedTo0() )
349 dof_type
const* __dof = M_u.functionSpace()->dof().get();
351 fe_type
const* __fe = M_u.functionSpace()->fe().get();
353 gm_ptrtype __gm(
new gm_type );
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 );
366 for ( uint16_type __f = 0; __f < geoelement_type::numTopologicalFaces; ++__f )
368 for ( permutation_type __p( permutation_type::IDENTITY );
369 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
371 __geopc[__f][__p] = geopc_ptrtype(
new geopc_type( __gm, __fe->points( __f ) ) );
373 FEELPP_ASSERT( __geopc[__f][__p]->nPoints() ).error(
"invalid number of points" );
377 uint16_type __face_id = __face_it->pos_first();
378 gmc_ptrtype __c(
new gmc_type( __gm, __face_it->element( 0 ), __geopc, __face_id ) );
380 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
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";
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 );
396 nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
398 DVLOG(2) <<
"nbFaceDof = " << nbFaceDof <<
"\n";
402 __face_it != this->endElement();
405 if ( !__face_it->isConnectedTo0() )
407 LOG( WARNING ) <<
"face not connected" << *__face_it;
414 if ( __face_it->isGhostFace() )
416 LOG(WARNING) <<
"face id : " << __face_it->id() <<
" is a ghost face";
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";
426 uint16_type __face_id = __face_it->pos_first();
427 __c->update( __face_it->element( 0 ), __face_id );
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";
432 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
434 t_expr_type expr( M_expr, mapgmc );
435 expr.update( mapgmc );
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";
442 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
443 for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
445 for ( uint16_type l = 0; l < nbFaceDof; ++l )
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";
454 boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
457 if ( std::find( dofs.begin(),
459 thedof ) != dofs.end() )
464 DVLOG(2) <<
"Eliminating row " << thedof <<
" using value : " << __value <<
"\n";
473 dofs.push_back( thedof );
474 values.push_back( __value );
483 __form.set( thedof, thedof, 1.0*1e30 );
484 M_rhs->set( thedof, __value*1e30 );
494 if ( __form.rowStartInMatrix()!=0)
496 auto const thedofshift = __form.rowStartInMatrix();
497 for (
auto itd=dofs.begin(),end=dofs.end() ; itd!=end ; ++itd)
500 auto x = M_rhs->clone();
502 CHECK( values.size() == dofs.size() ) <<
"Invalid dofs/values size: " << dofs.size() <<
"/" << values.size();
504 x->addVector( dofs.data(), dofs.size(), values.data() );
507 __form.zeroRows( dofs, *x, *M_rhs, M_on_strategy );
514 template<
typename T >
519 template<
typename T >
522 typedef typename T::vector_ptrtype type;
525 template<
typename Args>
526 struct integratoron_type
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;
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;
537 typedef _rhs_type the_rhs_type;
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;
547 typename V::vector_ptrtype
548 getRhsVector( V
const& v, mpl::false_ )
550 return v.vectorPtr();
555 getRhsVector( V
const& v, mpl::true_ )
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 )
566 return getRhsVector( v, Feel::detail::is_vector_ptr<V>() );
581 BOOST_PARAMETER_FUNCTION(
582 (
typename vf::detail::integratoron_type<Args>::expr_type ),
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>() )
601 typename vf::detail::integratoron_type<Args>::type ion( range,
603 Feel::vf::detail::getRhsVector(rhs),
608 LOG(INFO) <<
"Dirichlet condition over : "<<
nelements(range) <<
" faces";
612 LOG(INFO) <<
"treatment of Dirichlet condition: " << type <<
" (elimination, unsymmetric)";
615 LOG(INFO) <<
"treatment of Dirichlet condition: " << type <<
" (elimination and keep diagonal, unsymmetric)";
618 LOG(INFO) <<
"treatment of Dirichlet condition: " << type <<
" (elimination, symmetric, more expensive than unsymmetric treatment)";
621 LOG(INFO) <<
"treatment of Dirichlet condition: " << type <<
" (elimination and keep diagonal, symmetric, more expensive than unsymmetric treatment)";
624 LOG(INFO) <<
"treatment of Dirichlet condition: " << type <<
" (penalisation, symmetric, very big value on diagonal)";
631 return typename vf::detail::integratoron_type<Args>::expr_type( ion );