29 #ifndef __INTEGRATORDIRAC_HPP
30 #define __INTEGRATORDIRAC_HPP 1
32 #include <boost/timer.hpp>
33 #include <boost/foreach.hpp>
35 #include <feel/feelalg/enums.hpp>
50 template<
typename ElementRange,
typename Pts,
typename DiracExpr >
59 static const size_type context = DiracExpr::context;
61 static const uint16_type imorder = DiracExpr::imorder;
62 static const bool imIsPoly = DiracExpr::imIsPoly;
64 template<
typename Func>
65 struct HasTestFunction
67 static const bool result = DiracExpr::template HasTestFunction<Func>::result;
70 template<
typename Func>
71 struct HasTrialFunction
73 static const bool result = DiracExpr::template HasTrialFunction<Func>::result;
76 static const size_type iDim = boost::tuples::template element<0, ElementRange>::type::value;
78 typedef IntegratorDirac<ElementRange, Pts, DiracExpr> self_type;
80 typedef typename boost::tuples::template element<1, ElementRange>::type element_iterator;
82 typedef Pts points_type;
83 typedef DiracExpr expression_type;
84 typedef typename expression_type::value_type value_type;
91 typedef typename boost::remove_reference<typename element_iterator::reference>::type const_t;
92 typedef typename boost::remove_const<const_t>::type the_face_element_type;
93 typedef typename the_face_element_type::super2::template Element<the_face_element_type>::type the_element_type;
94 typedef the_element_type element_type;
95 typedef typename the_element_type::gm_type gm_type;
96 typedef boost::shared_ptr<gm_type> gm_ptrtype;
97 typedef typename gm_type::template Context<expression_type::context, the_element_type> gmc_type;
98 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
99 typedef typename gm_type::precompute_ptrtype gmcpc_ptrtype;
102 typedef typename expression_type::value_type value_type;
108 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
109 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
110 typedef typename eval_expr_type::shape shape;
116 typedef ublas::matrix<value_type> ret_type;
119 typedef ublas::matrix<value_type> ret_type;
126 IntegratorDirac( ElementRange
const& __elts,
128 expression_type
const& __expr )
130 M_eltbegin( __elts.template get<1>() ),
131 M_eltend( __elts.template get<2>() ),
136 IntegratorDirac( IntegratorDirac
const& ioe )
138 M_eltbegin( ioe.M_eltbegin ),
139 M_eltend( ioe.M_eltend ),
145 ~IntegratorDirac() {}
158 element_iterator beginElement()
const
167 element_iterator endElement()
const
182 template<
typename Elem1,
typename Elem2,
typename FormType>
183 void assemble( boost::shared_ptr<Elem1>
const& __u,
184 boost::shared_ptr<Elem2>
const& __v,
185 FormType& __f )
const
187 typedef typename Elem1::functionspace_type functionspace_type;
188 DVLOG(2) <<
"[IntegratorDirac::assemble()] is_same: "
189 << mpl::bool_<boost::is_same<functionspace_type,Elem1>::value>::value <<
"\n";
190 assemble( __u, __v, __f, mpl::bool_<boost::is_same<functionspace_type,Elem1>::value>() );
197 return evaluate( mpl::int_<iDim>() );
203 template<
typename Elem1,
typename Elem2,
typename FormType>
204 void assemble( boost::shared_ptr<Elem1>
const& ,
205 boost::shared_ptr<Elem2>
const& ,
206 FormType& , mpl::bool_<false> )
const {}
208 template<
typename Elem1,
typename Elem2,
typename FormType>
209 void assemble( boost::shared_ptr<Elem1>
const& __u,
210 boost::shared_ptr<Elem2>
const& __v,
211 FormType& __f, mpl::bool_<true> )
const;
213 ret_type evaluate( mpl::int_<MESH_ELEMENTS> )
const;
214 ret_type evaluate( mpl::int_<MESH_FACES> )
const;
218 element_iterator M_eltbegin;
219 element_iterator M_eltend;
222 expression_type M_expr;
225 template<
typename ElementRange,
typename Pts,
typename DiracExpr>
226 template<
typename Elem1,
typename Elem2,
typename FormType>
228 IntegratorDirac<ElementRange, Pts, DiracExpr>::assemble( boost::shared_ptr<Elem1>
const& ,
229 boost::shared_ptr<Elem2>
const& ,
231 mpl::bool_<true> )
const
234 std::map<std::string,std::pair<boost::timer,value_type> > timer;
235 timer[
"init intvrho"].first.restart();
237 element_type v( M_Space,
"v" );
238 std::fill( v.begin(), v.end(), .0 );
240 molecule_type::atoms_const_iterator_type atom( molecule.begin() );
242 if ( atom == molecule.end() )
return;
244 mesh_type::Inverse meshinv( M_mesh );
245 std::vector<value_type> atomcharges( molecule.size() );
248 for (
size_type atomid = 0; atom != molecule.end(); ++atom, ++atomid )
250 meshinv.addPointWithId( element_prod( M_invStretch , atom->center() - M_translation ),
252 atomcharges[atomid] = atom->charge();
255 meshinv.distribute();
257 std::vector<bool> dof_done( molecule.size() );
258 std::fill( dof_done.begin(), dof_done.end(), false );
259 std::vector<size_type> itab;
260 matrix_node<value_type>::type pts( mesh_type::nDim, 1 );
261 typedef mesh_type::element_type geoelement_type;
262 typedef mesh_type::element_iterator mesh_element_iterator;
264 mesh_element_iterator it = M_mesh->beginElementWithProcessId( M_mesh->comm().rank() );
265 mesh_element_iterator en = M_mesh->endElementWithProcessId( M_mesh->comm().rank() );
268 typedef mesh_type::gm_type gm_type;
269 typedef boost::shared_ptr<gm_type> gm_ptrtype;
270 typedef gm_type::Context<vm::POINT, geoelement_type> gmc_type;
271 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
274 typedef space_type::basis_type basis_type;
275 basis_type
const* __basis = M_Space->basis().get();
276 gm_ptrtype __gm = M_Space->gm();
277 typedef gm_type::precompute_ptrtype geopc_ptrtype;
278 typedef gm_type::precompute_type geopc_type;
279 geopc_ptrtype __geopc(
new geopc_type( __gm, __basis->dual().points() ) );
280 gmc_ptrtype __c(
new gmc_type( __gm, *it, __geopc ) );
284 if ( !v.areGlobalValuesUpdated() )
285 v.updateGlobalValues();
290 timer[
"init intvrho"].second = timer[
"init intvrho"].first.elapsed();
291 LOG(INFO) <<
"[timer] init intvrho(): " << timer[
"init intvrho"].second <<
"\n";
293 timer[
"intvrho"].first.restart();
296 for ( ; it != en; ++ it )
298 __c->update( *it, __geopc );
299 meshinv.pointsInConvex( it->id(), itab );
301 if ( itab.size() == 0 )
304 for (
size_type i = 0; i < itab.size(); ++i )
309 if ( !dof_done[dof] )
312 ublas::column( pts, 0 ) = meshinv.referenceCoords()[dof];
313 element_type::pc_type pc( v.functionSpace()->fe(), pts );
314 __geopc->update( pts );
315 __c->update( *it, __geopc );
317 for ( uint16_type loc_ind=0; loc_ind < basis_type::nLocalDof ; loc_ind++ )
319 for ( uint16_type comp = 0; comp < basis_type::nComponents; ++comp )
321 size_type globaldof = boost::get<0>( M_Space->dof()->localToGlobal( it->id(), loc_ind, comp ) );
324 if ( globaldof >= v.firstLocalIndex() &&
325 globaldof < v.lastLocalIndex() )
327 v.setGlobalValue( globaldof, 1 );
329 element_type::id_type interpfunc( v.id( *__c, pc ) );
332 rhs->add( globaldof, atomcharges[dof] * interpfunc( comp, 0, 0 ) );
335 v.setGlobalValue( globaldof, 0 );
345 timer[
"intvrho"].second = timer[
"intvrho"].first.elapsed();
346 LOG(INFO) <<
"[timer] intvrho(): " << timer[
"intvrho"].second <<
"\n";
349 template<
typename Elements,
typename Pts,
typename DiracExpr>
350 typename IntegratorDirac<Elements, Pts, DiracExpr>::ret_type
351 IntegratorDirac<Elements, Pts, DiracExpr>::evaluate( mpl::int_<MESH_ELEMENTS> )
const
353 mesh_type::Inverse meshinv( M_mesh );
355 pts_iterator ptit = M_pts.begin();
356 pts_iterator pten = M_pts.end();
359 for (
size_type ptid = 0; ptit != pten; ++ptit, ++ptid )
361 meshinv.addPointWithId( ptit->node(), ptid );
364 meshinv.distribute();
366 std::vector<bool> dof_done( M_pts.size() );
367 std::fill( dof_done.begin(), dof_done.end(), false );
368 std::vector<size_type> itab;
369 matrix_node<value_type>::type pts( mesh_type::nDim, 1 );
370 typedef mesh_type::element_type geoelement_type;
371 typedef mesh_type::element_iterator mesh_element_iterator;
373 mesh_element_iterator it = M_mesh->beginElementWithProcessId( M_mesh->comm().size() );
374 mesh_element_iterator en = M_mesh->endElementWithProcessId( M_mesh->comm().size() );
377 typedef mesh_type::gm_type gm_type;
378 typedef boost::shared_ptr<gm_type> gm_ptrtype;
379 typedef gm_type::Context<vm::POINT, geoelement_type> gmc_type;
380 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
386 gm_ptrtype gm(
new gm_type );
387 typename gm_type::precompute_ptrtype __geopc(
new typename gm_type::precompute_type( gm, pts ) );
392 return typename eval::ret_type( eval::shape::M, eval::shape::N );;
394 gmc_ptrtype __c(
new gmc_type( gm, *it, __geopc ) );
395 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
396 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
398 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
399 eval_expr_type expr( expression(), mapgmc );
400 typedef typename eval_expr_type::shape shape;
402 typename eval::ret_type res( eval::shape::M, eval::shape::N );
406 for ( ; it != en; ++ it )
408 __c->update( *it, __geopc );
409 meshinv.pointsInConvex( it->id(), itab );
411 if ( itab.size() == 0 )
414 for (
size_type i = 0; i < itab.size(); ++i )
419 if ( !dof_done[dof] )
422 ublas::column( pts, 0 ) = meshinv.referenceCoords()[dof];
423 element_type::pc_type pc( v.functionSpace()->fe(), pts );
424 __geopc->update( pts );
425 __c->update( *it, __geopc );
427 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
429 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
453 template<
typename ElementRange,
typename Po
intRange,
typename DiracExpr>
454 Expr<IntegratorDirac<ElementRange, PointRange, DiracExpr> >
455 dirac( ElementRange
const& elem_range, PointRange
const& pt_range, DiracExpr
const& expr )
457 typedef IntegratorDirac<ElementRange, PointRange, DiracExpr> expr_t;
458 return Expr<expr_t>( expr_t( elem_range, pt_range, expr ) );