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
integratordirac.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: 2008-02-10
7 
8  Copyright (C) 2008 Universite 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 3.0 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 __INTEGRATORDIRAC_HPP
30 #define __INTEGRATORDIRAC_HPP 1
31 
32 #include <boost/timer.hpp>
33 #include <boost/foreach.hpp>
34 
35 #include <feel/feelalg/enums.hpp>
36 
37 namespace Feel
38 {
39 namespace vf
40 {
42 
50 template<typename ElementRange, typename Pts, typename DiracExpr >
51 class IntegratorDirac
52 {
53 public:
54 
55 
59  static const size_type context = DiracExpr::context;
60 
61  static const uint16_type imorder = DiracExpr::imorder;
62  static const bool imIsPoly = DiracExpr::imIsPoly;
63 
64  template<typename Func>
65  struct HasTestFunction
66  {
67  static const bool result = DiracExpr::template HasTestFunction<Func>::result;
68  };
69 
70  template<typename Func>
71  struct HasTrialFunction
72  {
73  static const bool result = DiracExpr::template HasTrialFunction<Func>::result;
74  };
75 
76  static const size_type iDim = boost::tuples::template element<0, ElementRange>::type::value;
77 
78  typedef IntegratorDirac<ElementRange, Pts, DiracExpr> self_type;
79 
80  typedef typename boost::tuples::template element<1, ElementRange>::type element_iterator;
81 
82  typedef Pts points_type;
83  typedef DiracExpr expression_type;
84  typedef typename expression_type::value_type value_type;
85 
86  struct eval
87  {
88  //
89  // some typedefs
90  //
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;
100  //typedef typename eval_expr_type::value_type value_type;
101  //typedef typename strongest_numeric_type<typename Pts::value_type, typename expression_type::value_type>::type value_type;
102  typedef typename expression_type::value_type value_type;
103 
104  //
105  // Precompute some data in the reference element for
106  // geometric mapping and reference finite element
107  //
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;
111  //typedef typename shape_type::storage<value_type> storage_type;
112  /*
113  typedef mpl::if_<mpl::bool_<shape_type::is_scalar>,
114  mpl::identity<value_type>,
115  mpl::identity<ublas::vector<value_type,storage_type> > >::type::type ret_type;*/
116  typedef ublas::matrix<value_type> ret_type;
117  };
118 
119  typedef ublas::matrix<value_type> ret_type;
121 
125 
126  IntegratorDirac( ElementRange const& __elts,
127  Pts const& pts,
128  expression_type const& __expr )
129  :
130  M_eltbegin( __elts.template get<1>() ),
131  M_eltend( __elts.template get<2>() ),
132  M_pts( pts ),
133  M_expr( __expr )
134  {
135  }
136  IntegratorDirac( IntegratorDirac const& ioe )
137  :
138  M_eltbegin( ioe.M_eltbegin ),
139  M_eltend( ioe.M_eltend ),
140  M_pts( ioe.M_pts ),
141  M_expr( ioe.M_expr )
142  {
143  }
144 
145  ~IntegratorDirac() {}
146 
148 
152 
153 
158  element_iterator beginElement() const
159  {
160  return M_eltbegin;
161  }
162 
167  element_iterator endElement() const
168  {
169  return M_eltend;
170  }
171 
172 
174 
177 
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
186  {
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>() );
191  }
192 
193  //typename expression_type::template tensor<Geo_t>::value_type
194  ret_type
195  evaluate() const
196  {
197  return evaluate( mpl::int_<iDim>() );
198  }
199 
201 private:
202 
203  template<typename Elem1, typename Elem2, typename FormType>
204  void assemble( boost::shared_ptr<Elem1> const& /*__u*/,
205  boost::shared_ptr<Elem2> const& /*__v*/,
206  FormType& /*__f*/, mpl::bool_<false> ) const {}
207 
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;
212 
213  ret_type evaluate( mpl::int_<MESH_ELEMENTS> ) const;
214  ret_type evaluate( mpl::int_<MESH_FACES> ) const;
215 
216 private:
217 
218  element_iterator M_eltbegin;
219  element_iterator M_eltend;
220 
221  Pts const& M_pts;
222  expression_type M_expr;
223 };
224 
225 template<typename ElementRange, typename Pts, typename DiracExpr>
226 template<typename Elem1, typename Elem2, typename FormType>
227 void
228 IntegratorDirac<ElementRange, Pts, DiracExpr>::assemble( boost::shared_ptr<Elem1> const& /*__u*/,
229  boost::shared_ptr<Elem2> const& /*__v*/,
230  FormType& __form,
231  mpl::bool_<true> ) const
232 {
233 #if 0
234  std::map<std::string,std::pair<boost::timer,value_type> > timer;
235  timer["init intvrho"].first.restart();
236 
237  element_type v( M_Space, "v" );
238  std::fill( v.begin(), v.end(), .0 );
239 
240  molecule_type::atoms_const_iterator_type atom( molecule.begin() );
241 
242  if ( atom == molecule.end() ) return;
243 
244  mesh_type::Inverse meshinv( M_mesh );
245  std::vector<value_type> atomcharges( molecule.size() );
246 
247  /* initialisation of the mesh::inverse data structure */
248  for ( size_type atomid = 0; atom != molecule.end(); ++atom, ++atomid )
249  {
250  meshinv.addPointWithId( element_prod( M_invStretch , atom->center() - M_translation ),
251  atomid );
252  atomcharges[atomid] = atom->charge();
253  }
254 
255  meshinv.distribute();
256 
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;
263 
264  mesh_element_iterator it = M_mesh->beginElementWithProcessId( M_mesh->comm().rank() );
265  mesh_element_iterator en = M_mesh->endElementWithProcessId( M_mesh->comm().rank() );
266 
267  // geometric mapping context
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;
272 
273  // basis
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 ) );
281 
282 
283  /* note: meshmover.hpp has this two lines: */
284  if ( !v.areGlobalValuesUpdated() )
285  v.updateGlobalValues();
286 
287  /*shall we use this ?
288  */
289 
290  timer["init intvrho"].second = timer["init intvrho"].first.elapsed();
291  LOG(INFO) << "[timer] init intvrho(): " << timer["init intvrho"].second << "\n";
292 
293  timer["intvrho"].first.restart();
294 
295  // size_type first_dof = M_Space->dof()->firstDof();
296  for ( ; it != en; ++ it )
297  {
298  __c->update( *it, __geopc );
299  meshinv.pointsInConvex( it->id(), itab );
300 
301  if ( itab.size() == 0 )
302  continue;
303 
304  for ( size_type i = 0; i < itab.size(); ++i )
305  {
306  // get dof id in target dof table
307  size_type dof = itab[i];
308 
309  if ( !dof_done[dof] )
310  {
311  dof_done[dof]=true;
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 );
316 
317  for ( uint16_type loc_ind=0; loc_ind < basis_type::nLocalDof ; loc_ind++ )
318  {
319  for ( uint16_type comp = 0; comp < basis_type::nComponents; ++comp )
320  {
321  size_type globaldof = boost::get<0>( M_Space->dof()->localToGlobal( it->id(), loc_ind, comp ) );
322 
323  // update only values on the processor
324  if ( globaldof >= v.firstLocalIndex() &&
325  globaldof < v.lastLocalIndex() )
326  {
327  v.setGlobalValue( globaldof, 1 );
328 
329  element_type::id_type interpfunc( v.id( *__c, pc ) );
330  //std::cout << "interpfunc : " << interpfunc << "\n";
331 
332  rhs->add( globaldof, atomcharges[dof] * interpfunc( comp, 0, 0 ) );
333  // DVLOG(2) << "rhs( " << globaldof << ")=" << (*rhs)( globaldof )
334  // << " (just added " << atomcharges[dof] * interpfunc( comp, 0, 0 ) << " )" << "\n";
335  v.setGlobalValue( globaldof, 0 );
336 
337  }
338  }
339  }
340 
341  }
342  }
343  } // element
344 
345  timer["intvrho"].second = timer["intvrho"].first.elapsed();
346  LOG(INFO) << "[timer] intvrho(): " << timer["intvrho"].second << "\n";
347 #endif
348 }
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
352 {
353  mesh_type::Inverse meshinv( M_mesh );
354 
355  pts_iterator ptit = M_pts.begin();
356  pts_iterator pten = M_pts.end();
357 
358  /* initialisation of the mesh::inverse data structure */
359  for ( size_type ptid = 0; ptit != pten; ++ptit, ++ptid )
360  {
361  meshinv.addPointWithId( ptit->node(), ptid );
362  }
363 
364  meshinv.distribute();
365 
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;
372 
373  mesh_element_iterator it = M_mesh->beginElementWithProcessId( M_mesh->comm().size() );
374  mesh_element_iterator en = M_mesh->endElementWithProcessId( M_mesh->comm().size() );
375 
376  // geometric mapping context
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;
381 
382  //
383  // Precompute some data in the reference element for
384  // geometric mapping and reference finite element
385  //
386  gm_ptrtype gm( new gm_type );
387  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( gm, pts ) );
388 
389  // make sure that we have elements to iterate over (return 0
390  // otherwise)
391  if ( it == en )
392  return typename eval::ret_type( eval::shape::M, eval::shape::N );;
393 
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 ) );
397 
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;
401 
402  typename eval::ret_type res( eval::shape::M, eval::shape::N );
403  res.clear();
404 
405  // size_type first_dof = M_Space->dof()->firstDof();
406  for ( ; it != en; ++ it )
407  {
408  __c->update( *it, __geopc );
409  meshinv.pointsInConvex( it->id(), itab );
410 
411  if ( itab.size() == 0 )
412  continue;
413 
414  for ( size_type i = 0; i < itab.size(); ++i )
415  {
416  // get dof id in target dof table
417  size_type dof = itab[i];
418 
419  if ( !dof_done[dof] )
420  {
421  dof_done[dof]=true;
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 );
426 
427  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
428  {
429  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
430  {
431  //res( c1, c2 ) += ;
432  }
433  }
434 
435  }
436  }
437  } // element
438 
439  return res;
440 }
441 
442 
444 
453 template<typename ElementRange, typename PointRange, typename DiracExpr>
454 Expr<IntegratorDirac<ElementRange, PointRange, DiracExpr> >
455 dirac( ElementRange const& elem_range, PointRange const& pt_range, DiracExpr const& expr )
456 {
457  typedef IntegratorDirac<ElementRange, PointRange, DiracExpr> expr_t;
458  return Expr<expr_t>( expr_t( elem_range, pt_range, expr ) );
459 
460 }
461 
462 } // vf
463 } // feel
464 #endif

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