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
quadptlocalization.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): Vincent Chabannes <vincent.chabannes@imag.fr>
6  Date: 2008-01-03
7 
8  Copyright (C) 2011 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 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 */
30 #ifndef __quadptlocalization_H
31 #define __quadptlocalization_H 1
32 
33 
34 #include <feel/feeldiscr/mesh.hpp>
35 
36 #define FEELPP_EXPORT_QUADLOCALIZATION 0
37 #if FEELPP_EXPORT_QUADLOCALIZATION
38 #include <feel/feelfilters/exporter.hpp>
39 #endif
40 
41 
42 namespace Feel
43 {
44 
45 
46 template <typename IteratorRange,typename Im,typename Expr>
47 class QuadPtLocalization
48 {
49 public :
50 
51  static const size_type iDim = boost::tuples::template element<0, IteratorRange>::type::value;
52 
53  //static const size_type context = Expr::context|vm::POINT;
54  static const size_type context = mpl::if_< boost::is_same< mpl::int_<iDim>, mpl::int_<MESH_FACES> >,
55  mpl::int_<Expr::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT>,
56  mpl::int_<Expr::context|vm::POINT> >::type::value;
57 
58 
59  //expression_type::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT
60 
61  typedef IteratorRange range_iterator;
62 
63 
64  typedef typename boost::tuples::template element<1, IteratorRange>::type element_iterator_type;
65  typedef typename element_iterator_type::value_type::GeoShape GeoShape;
66  //typedef Mesh<GeoShape> mesh_type;
67 #if 0
68  typedef typename mesh_type::gm_type gm_type;
69  typedef typename mesh_type::element_type geoelement_type;
70 #else
71 
72  //typedef typename element_iterator_type::value_type geoelement_type;
73  //typedef typename geoelement_type::gm_type gm_type;
74 
75  typedef typename boost::remove_reference<typename element_iterator_type::reference>::type const_t;
76  typedef typename boost::remove_const<const_t>::type the_face_element_type;
77  typedef typename the_face_element_type::super2::template Element<the_face_element_type>::type the_element_type;
78  typedef the_element_type geoelement_type;
79  typedef typename geoelement_type::gm_type gm_type;
80 
81 
82 
83 #endif
84 
85 
86  typedef typename gm_type::template Context<context, geoelement_type> gmc_type;
87  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
88  typedef typename gm_type::precompute_type pc_type;
89  typedef typename gm_type::precompute_ptrtype pc_ptrtype;
90 #if 0
91  typedef typename mesh_type::value_type value_type;
92  typedef typename mesh_type::node_type node_type;
93  typedef typename matrix_node<value_type>::type matrix_node_type;
94 #else
95  typedef typename geoelement_type::value_type value_type;
96  typedef typename geoelement_type::node_type node_type;
97  typedef typename geoelement_type::matrix_node_type matrix_node_type;
98 #endif
99  //--------------------------------------------------------------------------------------//
100 
101  typedef typename mpl::if_<mpl::bool_<geoelement_type::is_simplex>,
102  mpl::identity<typename Im::template apply<geoelement_type::nDim, value_type, Simplex>::type >,
103  mpl::identity<typename Im::template apply<geoelement_type::nDim, value_type, Hypercube>::type >
104  >::type::type im_type;
105 
106  typedef typename im_type::face_quadrature_type im_face_type;
107 
108  typedef typename QuadMapped<im_type>::permutation_type permutation_type;
109  typedef typename QuadMapped<im_type>::permutation_points_type permutation_points_type;
110 
111  //--------------------------------------------------------------------------------------//
112 
113  // temporary container
114  typedef std::map< size_type,std::list<boost::tuple<size_type,size_type,node_type,node_type> > > result_temp1_type;//idTrial->list( q, idq, ptRefTest, ptRefTrial )
115  // sub container
116  typedef std::map< size_type,boost::tuple< std::vector< boost::tuple<size_type,size_type > >,matrix_node_type,matrix_node_type> > result_temp2_type;//idTrial->( mapFor<q,idq>, ptRefsTest, ptsRefTrial)
117  // result container
118  typedef std::list< boost::tuple< size_type, result_temp2_type> > result_container_type;// list(idTest, idTrial->( mapForQ, ptRefsTest, ptsRefTrial) )
119 
120  // result container for linearform
121  typedef std::list< boost::tuple<size_type,std::vector< boost::tuple<size_type,size_type > >,matrix_node_type > > result_container_linear_type;// list(idTest, mapForQ, ptRefsTest)
122 
123  //--------------------------------------------------------------------------------------//
124 
125 
126  QuadPtLocalization( IteratorRange const& elts /*, im_type const& *//*__im*/ )
127  :
128  M_listRange(),
129  M_im( ),
130  M_qm( ),
131  M_ppts( M_qm( this->im() ) ),
132  M_hasPrecompute(false)
133  {
134  M_listRange.push_back( elts );
135  }
136 
137  QuadPtLocalization( std::list<IteratorRange> const& elts )
138  :
139  M_listRange( elts ),
140  M_im( ),
141  M_qm( ),
142  M_ppts( M_qm( this->im() ) ),
143  M_hasPrecompute(false)
144  {}
145 
146 
150  im_type const& im() const
151  {
152  return M_im;
153  }
154 
158  im_face_type im( uint16_type f ) const
159  {
160  return M_im.face( f );
161  }
162 
166  element_iterator_type beginElement() const
167  {
168  return M_listRange.front().template get<1>();// M_eltbegin;
169  }
170 #if 0
171 
174  element_iterator_type endElement() const
175  {
176  return M_listRange.front().template get<2>();//M_eltend;
177  }
178 #endif
179 
182  result_container_linear_type const & resultLinear() const
183  {
184  return M_resLinear;
185  }
186 
190  result_container_type const & result() const
191  {
192  return M_resBilinear;
193  }
194 
195 
196  bool hasPrecompute() const { return M_hasPrecompute; }
197 
198  //--------------------------------------------------------------------------------------//
199 
200 
201  size_type eltForThisQuadPt( size_type theIdq, mpl::int_<MESH_ELEMENTS> )
202  {
203  return theIdq/this->im().nPoints();
204  }
205  size_type eltForThisQuadPt( size_type theIdq, mpl::int_<MESH_FACES> )
206  {
207  return theIdq/this->im().nPointsOnFace();
208  }
209 
210  size_type qIndexForThisGlobalQuadPt( size_type theIdq,mpl::int_<MESH_ELEMENTS> )
211  {
212  return theIdq%this->im().nPoints();
213  }
214  size_type qIndexForThisGlobalQuadPt( size_type theIdq,mpl::int_<MESH_FACES> )
215  {
216  return theIdq%this->im().nPointsOnFace();
217  }
218 
219  //--------------------------------------------------------------------------------------//
220 
221  gmc_ptrtype gmcForThisElt( size_type theIdElt,
222  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
223  mpl::int_<MESH_ELEMENTS> )
224  {
225 #if 1
226  element_iterator_type elt_it = this->beginElement(), elt_en = this->beginElement();
227  auto itListRange = M_listRange.begin();
228  auto const enListRange = M_listRange.end();
229  bool findElt=false;
230  for ( size_type ide = 0 ; itListRange!=enListRange && !findElt; ++itListRange)
231  {
232  boost::tie( boost::tuples::ignore, elt_it, elt_en ) = *itListRange;
233  const size_type distRange = std::distance(elt_it, elt_en);
234  if ( (ide+distRange-1) < theIdElt) { ide+=distRange; continue; }
235 
236  for ( ; ide<theIdElt ; ++ide ) ++elt_it;
237  findElt = true;
238  }
239 #else
240  //search element
241  auto elt_it = this->beginElement();
242  for ( size_type i=0; i<theIdElt; ++i ) ++elt_it;
243 #endif
244  // get only usefull quad point and reoder
245  uint16_type nContextPt = indexLocalToQuad.size();
246  matrix_node_type const& quadPtsRef = this->im().points();
247  matrix_node_type newquadPtsRef( quadPtsRef.size1() , nContextPt );
248 
249  for ( uint16_type i=0; i<nContextPt; ++i )
250  {
251  ublas::column( newquadPtsRef, i ) = ublas::column( quadPtsRef, indexLocalToQuad[i].template get<0>() );
252  }
253 
254  // create context
255  pc_ptrtype geopc( new pc_type( elt_it->gm(), newquadPtsRef ) );
256  gmc_ptrtype gmc( new gmc_type( elt_it->gm(),*elt_it, geopc ) );
257 
258  return gmc;
259  }
260 
261  //--------------------------------------------------------------------------------------//
262 
263  gmc_ptrtype gmcForThisElt( size_type theIdElt,
264  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
265  mpl::int_<MESH_FACES> )
266  {
267 #if 1
268  element_iterator_type elt_it=this->beginElement(), elt_en=this->beginElement();
269  auto itListRange = M_listRange.begin();
270  auto const enListRange = M_listRange.end();
271  bool findElt=false;
272  for ( size_type ide = 0 ; itListRange!=enListRange && !findElt; ++itListRange)
273  {
274  boost::tie( boost::tuples::ignore, elt_it, elt_en ) = *itListRange;
275  const size_type distRange = std::distance(elt_it, elt_en);
276  if ( (ide+distRange-1) < theIdElt) { ide+=distRange; continue; }
277 
278  for ( ; ide<theIdElt ; ++ide ) ++elt_it;
279  findElt = true;
280  }
281 #else
282  //search element
283  auto elt_it = this->beginElement();
284  for ( size_type i=0; i<theIdElt; ++i ) ++elt_it;
285 #endif
286  // get only usefull quad point and reoder
287  const uint16_type nContextPt = indexLocalToQuad.size();
288  const uint16_type __face_id_in_elt_0 = elt_it->pos_first();
289 
290  auto const __perm = elt_it->element( 0 ).permutation( __face_id_in_elt_0 );
291  matrix_node_type const& quadPtsRef = M_ppts[ __face_id_in_elt_0].find( __perm )->second;
292  //matrix_node_type quadPtsRef = this->im().points();
293  matrix_node_type newquadPtsRef( quadPtsRef.size1() , nContextPt );
294 
295  for ( uint16_type i=0; i<nContextPt; ++i )
296  {
297  ublas::column( newquadPtsRef, i ) = ublas::column( quadPtsRef, indexLocalToQuad[i].template get<0>() );
298  }
299 
300  // create context
301  std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( this->im().nFaces() );
302 
303  for ( uint16_type __f = 0; __f < this->im().nFaces(); ++__f )
304  {
305  for ( permutation_type __p( permutation_type::IDENTITY );
306  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
307  {
308  __geopc[__f][__p] = pc_ptrtype( new pc_type( elt_it->element( 0 ).gm(),
309  newquadPtsRef ) );
310  }
311  }
312 
313  gmc_ptrtype gmc( new gmc_type( elt_it->element( 0 ).gm(),
314  elt_it->element( 0 ),
315  __geopc,
316  __face_id_in_elt_0 ) );
317 
318  return gmc;
319  }
320 
321  //--------------------------------------------------------------------------------------//
322 
323  std::vector< boost::tuple< std::vector<boost::tuple<size_type,size_type> >,gmc_ptrtype,matrix_node_type > >
324  getUsableDataInFormContext( std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
325  matrix_node_type const & ptsRefTest )
326  {
327  // compute the number of several elements
328  uint16_type nContextPt = indexLocalToQuad.size();
329  std::map<size_type,std::list<size_type> > mapEltId;
330 
331  for ( uint16_type i=0; i<nContextPt; ++i )
332  {
333  size_type eltId = this->eltForThisQuadPt( indexLocalToQuad[i].template get<1>(),mpl::int_<iDim>() );
334  mapEltId[eltId].push_back( i );
335  }
336 
337  // number of elements
338  auto nEltInContext = mapEltId.size();
339 
340  // the vector result
341  std::vector< boost::tuple< std::vector<boost::tuple<size_type,size_type> >, gmc_ptrtype, matrix_node_type > > vec_res( nEltInContext );
342 
343  auto map_it = mapEltId.begin();
344  auto map_en = mapEltId.end();
345  size_type cptId = 0;
346 
347  for ( ; map_it!= map_en ; ++map_it, ++cptId )
348  {
349  // subdivide into elements the map indexLocalToQuad
350  std::vector<boost::tuple<size_type,size_type> > newindexLocalToQuad( map_it->second.size() );
351  // subdivide into elements the ptsRef
352  matrix_node_type newptsRefTest( ptsRefTest.size1(), map_it->second.size() );
353  auto sublist_it = map_it->second.begin();
354  auto sublist_en = map_it->second.end();
355  size_type index = 0;
356 
357  for ( ; sublist_it != sublist_en ; ++sublist_it,++index )
358  {
359  newindexLocalToQuad[index] = indexLocalToQuad[*sublist_it];
360  ublas::column( newptsRefTest, index ) = ublas::column( ptsRefTest,*sublist_it );
361  }
362  // get the corresponding gmc
363  auto thegmc = gmcForThisElt( map_it->first,newindexLocalToQuad,mpl::int_<iDim>() );
364  // add to result
365  vec_res[cptId] = boost::make_tuple( newindexLocalToQuad, thegmc, newptsRefTest );
366  }
367 
368  return vec_res;
369  }
370 
371  //--------------------------------------------------------------------------------------//
372 
373  std::vector< boost::tuple< std::vector<boost::tuple<size_type,size_type> >,gmc_ptrtype,matrix_node_type, matrix_node_type > >
374  getUsableDataInFormContext( std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
375  matrix_node_type const & ptsRefTest,
376  matrix_node_type const & ptsRefTrial )
377  {
378  // compute the number of several elements
379  uint16_type nContextPt = indexLocalToQuad.size();
380  std::map<size_type,std::list<size_type> > mapEltId;
381 
382  for ( uint16_type i=0; i<nContextPt; ++i )
383  {
384  size_type eltId = this->eltForThisQuadPt( indexLocalToQuad[i].template get<1>(),mpl::int_<iDim>() );
385  mapEltId[eltId].push_back( i );
386  }
387 
388  // number of elements
389  auto nEltInContext = mapEltId.size();
390 
391  // the vector result
392  std::vector< boost::tuple< std::vector<boost::tuple<size_type,size_type> >, gmc_ptrtype, matrix_node_type, matrix_node_type > > vec_res( nEltInContext );
393 
394  auto map_it = mapEltId.begin();
395  auto map_en = mapEltId.end();
396  size_type cptId = 0;
397 
398  for ( ; map_it!= map_en ; ++map_it, ++cptId )
399  {
400  // subdivide into elements the map indexLocalToQuad
401  std::vector<boost::tuple<size_type,size_type> > newindexLocalToQuad( map_it->second.size() );
402  // subdivide into elements the ptsRef
403  matrix_node_type newptsRefTest( ptsRefTest.size1(), map_it->second.size() );
404  matrix_node_type newptsRefTrial( ptsRefTrial.size1(), map_it->second.size() );
405 
406  auto sublist_it = map_it->second.begin();
407  auto sublist_en = map_it->second.end();
408  size_type index = 0;
409 
410  for ( ; sublist_it != sublist_en ; ++sublist_it,++index )
411  {
412  newindexLocalToQuad[index] = indexLocalToQuad[*sublist_it];
413  ublas::column( newptsRefTest, index ) = ublas::column( ptsRefTest,*sublist_it );
414  ublas::column( newptsRefTrial, index ) = ublas::column( ptsRefTrial,*sublist_it );
415  }
416 
417  // get the corresponding gmc
418  auto thegmc = gmcForThisElt( map_it->first,newindexLocalToQuad,mpl::int_<iDim>() );
419  // add to result
420  vec_res[cptId] = boost::make_tuple( newindexLocalToQuad, thegmc, newptsRefTest,newptsRefTrial );
421  }
422 
423  return vec_res;
424  }
425 
426  //--------------------------------------------------------------------------------------//
427 
428  template <typename Mesh1Type>
429  void
430  localization( mpl::int_<MESH_FACES> ,
431  boost::shared_ptr<Mesh1Type> meshTest,
432  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > & testEltToPtsQuad )
433  {
434 
435  auto begin_elt_it = this->beginElement();
436 
437  std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( this->im().nFaces() );
438  //typedef typename im_type::face_quadrature_type face_im_type;
439 
440  //std::vector<face_im_type> face_ims( this->im().nFaces() );
441 
442  for ( uint16_type __f = 0; __f < this->im().nFaces(); ++__f )
443  {
444  //face_ims[__f] = this->im( __f );
445 
446  for ( permutation_type __p( permutation_type::IDENTITY );
447  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
448  {
449  //FEELPP_ASSERT( ppts[__f].find(__p)->second.size2() != 0 ).warn( "invalid quadrature type" );
450  __geopc[__f][__p] = pc_ptrtype( new pc_type( begin_elt_it->element( 0 ).gm(), M_ppts[__f].find( __p )->second ) );
451  }
452  }
453 
454 
455  uint16_type __face_id_in_elt_0 = begin_elt_it->pos_first();
456 
457  gmc_ptrtype gmc( new gmc_type( begin_elt_it->element( 0 ).gm(),
458  begin_elt_it->element( 0 ),
459  __geopc,
460  __face_id_in_elt_0 ) );
461 
462  auto meshTestLocalization = meshTest->tool_localization();
463  meshTestLocalization->updateForUse();
464  const auto nbNearNeighborAtStartTest = meshTestLocalization->kdtree()->nPtMaxNearNeighbor();
465  bool notUseOptLocTest = Mesh1Type::nDim!=Mesh1Type::nRealDim;
466  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(Mesh1Type::element_type::numPoints);
467 
468  matrix_node_type ptsReal( begin_elt_it->vertices().size1(), 1 );
469  size_type testIdElt=0;
470  node_type testNodeRef;
471 
472  element_iterator_type elt_it, elt_en;
473 
474  auto itListRange = M_listRange.begin();
475  auto const enListRange = M_listRange.end();
476  for ( size_type ide = 0 ; itListRange!=enListRange ; ++itListRange)
477  {
478  boost::tie( boost::tuples::ignore, elt_it, elt_en ) = *itListRange;
479  for ( ; elt_it != elt_en; ++elt_it, ++ide )
480  {
481  __face_id_in_elt_0 = elt_it->pos_first();
482  //if ( elt_it->isConnectedTo1()) std::cout << "\n AIEEEEEE!!!!!!!!!!!\n";
483 
484  gmc->update( elt_it->element( 0 ), __face_id_in_elt_0 );
485 
486  //std::cout << "\n quad gmc "<< gmc->xReal();
487  for ( int q = 0; q < gmc->nPoints(); ++ q )
488  {
489  size_type idq = gmc->nPoints()*ide+q;
490 #if 0
491  auto testAnalysis = meshTestLocalization->searchElement( gmc->xReal( q ) );
492  testIdElt = testAnalysis.template get<1>();
493  testNodeRef = testAnalysis.template get<2>();
494 #else
495  ublas::column(ptsReal,0 ) = gmc->xReal( q );
496  if (notUseOptLocTest) testIdElt=invalid_size_type_value;
497  auto resLocalisationTest = meshTestLocalization->run_analysis(ptsReal,testIdElt,elt_it->vertices(),mpl::int_<0>());
498  testIdElt = resLocalisationTest.template get<1>();
499  testNodeRef = meshTestLocalization->result_analysis().begin()->second.begin()->template get<1>();
500 #endif
501  testEltToPtsQuad[testIdElt].push_back( boost::make_tuple( idq,q,testNodeRef ) );
502  }
503  }
504  }
505 
506  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTest);
507 
508  }
509 
510  //--------------------------------------------------------------------------------------//
511 
512  template <typename Mesh1Type>
513  void
514  localization( mpl::int_<MESH_ELEMENTS> ,
515  boost::shared_ptr<Mesh1Type> meshTest,
516  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > & testEltToPtsQuad )
517  {
518 
519  auto begin_elt_it = this->beginElement();
520  //auto elt_en = this->endElement();
521 
522  pc_ptrtype geopc( new pc_type( begin_elt_it->gm(), this->im().points() ) );
523  gmc_ptrtype gmc( new gmc_type( begin_elt_it->gm(),*begin_elt_it, geopc ) );
524 
525  auto meshTestLocalization = meshTest->tool_localization();
526  meshTestLocalization->updateForUse();
527  const auto nbNearNeighborAtStartTest = meshTestLocalization->kdtree()->nPtMaxNearNeighbor();
528  bool notUseOptLocTest = Mesh1Type::nDim!=Mesh1Type::nRealDim;
529  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(Mesh1Type::element_type::numPoints);
530 
531  matrix_node_type ptsReal( begin_elt_it->vertices().size1(), 1 );
532  size_type testIdElt=0;
533  node_type testNodeRef;
534 
535  element_iterator_type elt_it, elt_en;
536 
537  auto itListRange = M_listRange.begin();
538  auto const enListRange = M_listRange.end();
539  for ( size_type ide = 0 ; itListRange!=enListRange ; ++itListRange)
540  {
541  boost::tie( boost::tuples::ignore, elt_it, elt_en ) = *itListRange;
542  for ( ; elt_it != elt_en; ++elt_it, ++ide )
543  {
544  gmc->update( *elt_it );
545 
546  for ( int q = 0; q < gmc->nPoints(); ++ q )
547  {
548  // cpt of quad pt
549  size_type idq = gmc->nPoints()*ide+q;
550  // search in test mesh
551 #if 0
552  auto testAnalysis = meshTestLocalization->searchElement( gmc->xReal( q ) );
553  testIdElt = testAnalysis.template get<1>();
554  testNodeRef = testAnalysis.template get<2>();
555 #else
556  ublas::column(ptsReal,0 ) = gmc->xReal( q );
557  if (notUseOptLocTest) testIdElt=invalid_size_type_value;
558  auto resLocalisationTest = meshTestLocalization->run_analysis(ptsReal,testIdElt,elt_it->vertices(),mpl::int_<0>());
559  testIdElt = resLocalisationTest.template get<1>();
560  testNodeRef = meshTestLocalization->result_analysis().begin()->second.begin()->template get<1>();
561 #endif
562 
563  testEltToPtsQuad[testIdElt].push_back( boost::make_tuple( idq,q,testNodeRef ) );
564  }
565  }
566  }
567 
568  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTest);
569 
570  }
571 
572  //--------------------------------------------------------------------------------------//
573 
574  template <typename Mesh1Type,typename Mesh2Type>
575  void
576  localization( mpl::int_<MESH_FACES> ,
577  boost::shared_ptr<Mesh1Type> meshTest,
578  boost::shared_ptr<Mesh2Type> meshTrial,
579  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > & testEltToPtsQuad,
580  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > & trialEltToPtsQuad,
581  std::vector<std::list<size_type> > & EltCoupled )
582  {
583  //std::cout << "[QuadPtLocalization] : localization<MESH_FACES>(bilinear form start" << std::endl;
584 
585  auto begin_elt_it = this->beginElement();
586 
587  std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( this->im().nFaces() );
588 
589  for ( uint16_type __f = 0; __f < this->im().nFaces(); ++__f )
590  {
591  for ( permutation_type __p( permutation_type::IDENTITY );
592  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
593  {
594  //FEELPP_ASSERT( ppts[__f].find(__p)->second.size2() != 0 ).warn( "invalid quadrature type" );
595  __geopc[__f][__p] = pc_ptrtype( new pc_type( begin_elt_it->element( 0 ).gm(), M_ppts[__f].find( __p )->second ) );
596  }
597  }
598 
599  uint16_type __face_id_in_elt_0 = begin_elt_it->pos_first();
600 
601  gmc_ptrtype gmc( new gmc_type( begin_elt_it->element( 0 ).gm(),
602  begin_elt_it->element( 0 ),
603  __geopc,
604  __face_id_in_elt_0 ) );
605 
606 
607  auto meshTrialLocalization = meshTrial->tool_localization();
608  meshTrialLocalization->updateForUse();
609  const auto nbNearNeighborAtStartTrial = meshTrialLocalization->kdtree()->nPtMaxNearNeighbor();
610  bool notUseOptLocTrial = Mesh2Type::nDim!=Mesh2Type::nRealDim;
611  if (notUseOptLocTrial) meshTrialLocalization->kdtree()->nbNearNeighbor(Mesh2Type::element_type::numPoints);
612 
613  auto meshTestLocalization = meshTest->tool_localization();
614  meshTestLocalization->updateForUse();
615  const auto nbNearNeighborAtStartTest = meshTestLocalization->kdtree()->nPtMaxNearNeighbor();
616  bool notUseOptLocTest = Mesh1Type::nDim!=Mesh1Type::nRealDim;
617  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(Mesh1Type::element_type::numPoints);
618 
619 
620  matrix_node_type ptsReal( begin_elt_it->vertices().size1(), 1 );
621  size_type trialIdElt = 0, testIdElt=0;
622  node_type trialNodeRef,testNodeRef;
623 
624  bool quadMeshIsSameThatTrialMesh=false,quadMeshIsSameThatTestMesh=false;
625  if ( dynamic_cast<void*>( const_cast<MeshBase*>( begin_elt_it->mesh() ) ) == dynamic_cast<void*>( meshTrial.get() ) )
626  quadMeshIsSameThatTrialMesh=true;
627  if ( dynamic_cast<void*>( const_cast<MeshBase*>( begin_elt_it->mesh() ) ) == dynamic_cast<void*>( meshTest.get() ) )
628  quadMeshIsSameThatTestMesh=true;
629 
630 #if FEELPP_EXPORT_QUADLOCALIZATION
631  std::map<size_type,std::list<size_type> > mapBetweenMeshes_test;
632  std::map<size_type,std::list<size_type> > mapBetweenMeshes_trial;
633 #endif
634 
635  element_iterator_type elt_it, elt_en;
636 
637  auto itListRange = M_listRange.begin();
638  auto const enListRange = M_listRange.end();
639  for ( size_type ide = 0 ; itListRange!=enListRange ; ++itListRange)
640  {
641  boost::tie( boost::tuples::ignore, elt_it, elt_en ) = *itListRange;
642  for ( ; elt_it != elt_en; ++elt_it, ++ide )
643  {
644  __face_id_in_elt_0 = elt_it->pos_first();
645  //if ( elt_it->isConnectedTo1()) std::cout << "\n AIEEEEEE!!!!!!!!!!!\n";
646 
647  gmc->update( elt_it->element( 0 ), __face_id_in_elt_0 );
648 
649  //std::cout << "\n quad gmc "<< gmc->xReal();
650 
651  for ( int q = 0; q < gmc->nPoints(); ++ q )
652  {
653  // cpt of quad pt
654  size_type idq = gmc->nPoints()*ide+q;
655 
656  // search in trial mesh
657  ublas::column(ptsReal,0 ) = gmc->xReal( q );
658  if (!quadMeshIsSameThatTrialMesh)
659  {
660  if (notUseOptLocTrial) trialIdElt=invalid_size_type_value;
661  auto resLocalisationTrial = meshTrialLocalization->run_analysis(ptsReal,trialIdElt,elt_it->vertices(),mpl::int_<0>());
662  trialIdElt = resLocalisationTrial.template get<1>();
663  trialNodeRef = meshTrialLocalization->result_analysis().begin()->second.begin()->template get<1>();
664  }
665  else
666  {
667  trialIdElt = gmc->id();
668  trialNodeRef = gmc->xRef(q);
669  }
670  trialEltToPtsQuad[trialIdElt].push_back( boost::make_tuple( idq,q,trialNodeRef ) );
671 
672  // search in test mesh
673  if (!quadMeshIsSameThatTestMesh)
674  {
675  if (notUseOptLocTest) testIdElt=invalid_size_type_value;
676  auto resLocalisationTest = meshTestLocalization->run_analysis(ptsReal,testIdElt,elt_it->vertices(),mpl::int_<0>());
677  testIdElt = resLocalisationTest.template get<1>();
678  testNodeRef = meshTestLocalization->result_analysis().begin()->second.begin()->template get<1>();
679  }
680  else
681  {
682  testIdElt = gmc->id();
683  testNodeRef = gmc->xRef(q);
684  }
685  testEltToPtsQuad[testIdElt].push_back( boost::make_tuple( idq,q,testNodeRef ) );
686 
687  // relation between test and trial
688  if ( std::find( EltCoupled[testIdElt].begin(),EltCoupled[testIdElt].end(),trialIdElt )==EltCoupled[testIdElt].end() )
689  {
690  EltCoupled[testIdElt].push_back( trialIdElt );
691  /*std::cout << "gmc->xRef(q)" << gmc->xReal(q)
692  << " testIdElt " << testIdElt << " meshTest->element().G()" << meshTest->element(testIdElt).G()
693  << " trialIdElt " << trialIdElt << " meshTrial->element().G()" << meshTrial->element(trialIdElt).G()
694  << std::endl;*/
695 #if FEELPP_EXPORT_QUADLOCALIZATION
696  mapBetweenMeshes_test[testIdElt].push_back( trialIdElt );
697  mapBetweenMeshes_trial[trialIdElt].push_back( testIdElt );
698 #endif
699  }
700  } // for ( int q = 0; q < gmc->nPoints(); ++ q )
701  } // for ( ; elt_it ...)
702  } // end for( size_type ide ... )
703 #if FEELPP_EXPORT_QUADLOCALIZATION
704  this->exportQuadLocalization(meshTest,meshTrial,mapBetweenMeshes_test,mapBetweenMeshes_trial);
705 #endif
706 
707  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTest);
708  if (notUseOptLocTrial) meshTrialLocalization->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTrial);
709 
710  //std::cout << "[QuadPtLocalization] : localization<MESH_FACES>(bilinear form finish" << std::endl;
711 
712  } // localization
713 
714  //--------------------------------------------------------------------------------------//
715 
716  template <typename Mesh1Type,typename Mesh2Type>
717  void
718  localization( mpl::int_<MESH_ELEMENTS> ,
719  boost::shared_ptr<Mesh1Type> meshTest,
720  boost::shared_ptr<Mesh2Type> meshTrial,
721  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > & testEltToPtsQuad,
722  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > & trialEltToPtsQuad,
723  std::vector<std::list<size_type> > & EltCoupled )
724  {
725  //std::cout << "[QuadPtLocalization] : localization<MESH_ELEMENTS>(bilinear form start" << std::endl;
726 
727  auto begin_elt_it = this->beginElement();
728  //auto elt_it = this->beginElement();
729  //auto elt_en = this->endElement();
730 
731  pc_ptrtype geopc( new pc_type( begin_elt_it->gm(), this->im().points() ) );
732  gmc_ptrtype gmc( new gmc_type( begin_elt_it->gm(),*begin_elt_it, geopc ) );
733 
734  auto meshTrialLocalization = meshTrial->tool_localization();
735  meshTrialLocalization->updateForUse();
736  const auto nbNearNeighborAtStartTrial = meshTrialLocalization->kdtree()->nPtMaxNearNeighbor();
737  bool notUseOptLocTrial = Mesh2Type::nDim!=Mesh2Type::nRealDim;
738  if (notUseOptLocTrial) meshTrialLocalization->kdtree()->nbNearNeighbor(Mesh2Type::element_type::numPoints);
739 
740  auto meshTestLocalization = meshTest->tool_localization();
741  meshTestLocalization->updateForUse();
742  const auto nbNearNeighborAtStartTest = meshTestLocalization->kdtree()->nPtMaxNearNeighbor();
743  bool notUseOptLocTest = Mesh1Type::nDim!=Mesh1Type::nRealDim;
744  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(Mesh1Type::element_type::numPoints);
745 
746 #if FEELPP_EXPORT_QUADLOCALIZATION
747  std::map<size_type,std::list<size_type> > mapBetweenMeshes_test;
748  std::map<size_type,std::list<size_type> > mapBetweenMeshes_trial;
749 #endif
750 
751 
752  //auto nQuadPtsInElt = this->im().nPoints();
753  //auto nElts = std::distance(elt_it,elt_en);
754  //auto nQuadPts = nElts*nQuadPtsInElt;
755  // auto nEltTrial= meshTrial->numElements();
756  //auto nEltTest= meshTest->numElements();
757 
758  matrix_node_type ptsReal( begin_elt_it->vertices().size1(), 1 );
759  size_type trialIdElt = 0, testIdElt=0;
760  node_type trialNodeRef,testNodeRef;
761 
762  bool quadMeshIsSameThatTrialMesh=false,quadMeshIsSameThatTestMesh=false;
763  if ( dynamic_cast<void*>( const_cast<MeshBase*>( begin_elt_it->mesh() ) ) == dynamic_cast<void*>( meshTrial.get() ) )
764  quadMeshIsSameThatTrialMesh=true;
765  if ( dynamic_cast<void*>( const_cast<MeshBase*>( begin_elt_it->mesh() ) ) == dynamic_cast<void*>( meshTest.get() ) )
766  quadMeshIsSameThatTestMesh=true;
767 
768  element_iterator_type elt_it, elt_en;
769 
770  auto itListRange = M_listRange.begin();
771  auto const enListRange = M_listRange.end();
772  for ( size_type ide = 0 ; itListRange!=enListRange ; ++itListRange)
773  {
774  boost::tie( boost::tuples::ignore, elt_it, elt_en ) = *itListRange;
775  for ( ; elt_it != elt_en; ++elt_it, ++ide )
776  {
777  gmc->update( *elt_it );
778 
779  for ( int q = 0; q < gmc->nPoints(); ++ q )
780  {
781  // cpt of quad pt
782  size_type idq = gmc->nPoints()*ide+q;
783 
784  // search in trial mesh
785  ublas::column(ptsReal,0 ) = gmc->xReal( q );
786  if (!quadMeshIsSameThatTrialMesh)
787  {
788  if (notUseOptLocTrial) trialIdElt=invalid_size_type_value;
789  auto resLocalisationTrial = meshTrialLocalization->run_analysis(ptsReal,trialIdElt,elt_it->vertices(),mpl::int_<0>());
790  trialIdElt = resLocalisationTrial.template get<1>();
791  trialNodeRef = meshTrialLocalization->result_analysis().begin()->second.begin()->template get<1>();
792  }
793  else
794  {
795  trialIdElt = gmc->id();
796  trialNodeRef = gmc->xRef(q);
797  }
798  //std::cout << "\n trialNodeRef " << trialNodeRef << std::endl;
799  trialEltToPtsQuad[trialIdElt].push_back( boost::make_tuple( idq,q,trialNodeRef ) );
800 
801  // search in test mesh
802  if (!quadMeshIsSameThatTestMesh)
803  {
804  if (notUseOptLocTest) testIdElt=invalid_size_type_value;
805  auto resLocalisationTest = meshTestLocalization->run_analysis(ptsReal,testIdElt,elt_it->vertices(),mpl::int_<0>());
806  testIdElt = resLocalisationTest.template get<1>();
807  testNodeRef = meshTestLocalization->result_analysis().begin()->second.begin()->template get<1>();
808  }
809  else
810  {
811  testIdElt = gmc->id();
812  testNodeRef = gmc->xRef(q);
813  }
814  testEltToPtsQuad[testIdElt].push_back( boost::make_tuple( idq,q,testNodeRef ) );
815 
816  //std::cout << " AAtestIdElt " << testIdElt << " AAtrialIdElt " << trialIdElt << std::endl;
817 
818  // relation between test and trial
819  if ( std::find( EltCoupled[testIdElt].begin(),EltCoupled[testIdElt].end(),trialIdElt )==EltCoupled[testIdElt].end() )
820  {
821  EltCoupled[testIdElt].push_back( trialIdElt );
822  //std::cout << " AAtestIdElt " << testIdElt << " AAtrialIdElt " << trialIdElt << std::endl;
823 #if FEELPP_EXPORT_QUADLOCALIZATION
824  mapBetweenMeshes_test[testIdElt].push_back( trialIdElt );
825  mapBetweenMeshes_trial[trialIdElt].push_back( testIdElt );
826 #endif
827  }
828 
829 
830 #if 0 // try with neighbors
831  auto const& geoeltTrial = meshTrial->element( trialIdElt );
832  for ( uint16_type ms=0; ms < geoeltTrial.nNeighbors(); ms++ )
833  {
834  const size_type neighborTrial_id = geoeltTrial.neighbor( ms ).first;
835  if ( neighborTrial_id==invalid_size_type_value ) continue;
836 
837  auto resNeighboor = meshTrialLocalization->isIn( neighborTrial_id,gmc->xReal( q ),elt_it->vertices(),mpl::int_<1>() );
838  if (resNeighboor.get<0>())
839  {
840  trialEltToPtsQuad[neighborTrial_id].push_back( boost::make_tuple( idq,q,resNeighboor.get<1>() ) );
841 
842  if ( std::find( EltCoupled[testIdElt].begin(),EltCoupled[testIdElt].end(),neighborTrial_id )==EltCoupled[testIdElt].end() )
843  {
844  EltCoupled[testIdElt].push_back( neighborTrial_id );
845  }
846 #if FEELPP_EXPORT_QUADLOCALIZATION
847  mapBetweenMeshes_trial[neighborTrial_id].push_back( testIdElt );
848 #endif
849  }
850  } // for ( uint16_type ms=0;...
851 
852  auto const& geoeltTest = meshTest->element( testIdElt );
853  for ( uint16_type ms=0; ms < geoeltTest.nNeighbors(); ms++ )
854  {
855  const size_type neighborTest_id = geoeltTest.neighbor( ms ).first;
856  if ( neighborTest_id==invalid_size_type_value ) continue;
857 
858  auto resNeighboor = meshTestLocalization->isIn( neighborTest_id,gmc->xReal( q ),elt_it->vertices(),mpl::int_<1>() );
859  if (resNeighboor.get<0>())
860  {
861  testEltToPtsQuad[neighborTest_id].push_back( boost::make_tuple( idq,q,resNeighboor.get<1>() ) );
862 
863  if ( std::find( EltCoupled[neighborTest_id].begin(),EltCoupled[neighborTest_id].end(),trialIdElt )==EltCoupled[neighborTest_id].end() )
864  {
865  EltCoupled[neighborTest_id].push_back( trialIdElt );
866  }
867 #if FEELPP_EXPORT_QUADLOCALIZATION
868  mapBetweenMeshes_test[neighborTest_id].push_back( trialIdElt );
869 #endif
870  }
871  } // for ( uint16_type ms=0;...
872 #endif
873  } // for ( int q = 0; q < gmc->nPoints(); ++ q )
874  } // for ( ; elt_it != elt_en; ++elt_it, ++ide )
875  } // end for( size_type ide ... )
876 #if FEELPP_EXPORT_QUADLOCALIZATION
877  this->exportQuadLocalization(meshTest,meshTrial,mapBetweenMeshes_test,mapBetweenMeshes_trial);
878 #endif
879 
880  if (notUseOptLocTest) meshTestLocalization->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTest);
881  if (notUseOptLocTrial) meshTrialLocalization->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTrial);
882 
883  //std::cout << "[QuadPtLocalization] : localization<MESH_ELEMENTS>(bilinear form finish" << std::endl;
884  }
885 
886  //--------------------------------------------------------------------------------------//
887 
888  template <typename Mesh1Type>
889  void
890  update( boost::shared_ptr<Mesh1Type> meshTest )
891  {
892  typedef Mesh1Type mesh_test_type;
893 
894  //clean result
895  M_resLinear.clear();
896 
897  // localize quad pts on meshTest -> storage in testEltToPtsQuad
898  auto nEltTest= meshTest->numElements();
899  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > testEltToPtsQuad( nEltTest );
900  this->localization( mpl::int_<iDim>(),meshTest,testEltToPtsQuad );
901 
902  //build a efficent container : M_resLinear
903  size_type theIdEltTest = 0;
904  auto eltTest_it = testEltToPtsQuad.begin();
905  auto eltTest_en = testEltToPtsQuad.end();
906 
907  for ( ; eltTest_it != eltTest_en ; ++eltTest_it, ++theIdEltTest )
908  {
909  if ( eltTest_it->size()>0 )
910  {
911  auto nPtsRef = eltTest_it->size();
912  matrix_node_type ptsRefTest( mesh_test_type::nDim,nPtsRef );
913  std::vector<boost::tuple<size_type,size_type> > indexLocalToQuad( nPtsRef );
914  // build
915  auto idq_it = eltTest_it->begin();
916  auto idq_en = eltTest_it->end();
917  uint16_type cptIdq = 0;
918 
919  for ( ; idq_it != idq_en ; ++idq_it, ++cptIdq )
920  {
921  indexLocalToQuad[cptIdq] = boost::make_tuple( idq_it->template get<1>(),idq_it->template get<0>() );
922  ublas::column( ptsRefTest, cptIdq ) = idq_it->template get<2>();
923  }
924 
925  M_resLinear.push_back( boost::make_tuple( theIdEltTest, indexLocalToQuad, ptsRefTest ) );
926  }
927  }
928 
929  }
930 
931  //--------------------------------------------------------------------------------------//
932 
933  template <typename Mesh1Type,typename Mesh2Type>
934  void
935  update( boost::shared_ptr<Mesh1Type> meshTest, boost::shared_ptr<Mesh2Type> meshTrial )
936  {
937  typedef Mesh1Type mesh_test_type;
938  typedef Mesh2Type mesh_trial_type;
939 
940  //clean result
941  M_resBilinear.clear();
942 
943  auto nEltTrial= meshTrial->numElements();
944  auto nEltTest= meshTest->numElements();
945  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > trialEltToPtsQuad( nEltTrial ) ;
946  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > testEltToPtsQuad( nEltTest );
947  std::vector<std::list<size_type> > EltCoupled( nEltTest );
948 
949  // localize quad pts on meshTest and meshTrial -> storage in testEltToPtsQuad, trialEltToPtsQuad and EltCoupled
950  this->localization( mpl::int_<iDim>(),
951  meshTest,
952  meshTrial,
953  testEltToPtsQuad,
954  trialEltToPtsQuad,
955  EltCoupled );
956 
957  //build a efficent container : M_resBilinear
958  auto eltCoupled_it = EltCoupled.begin();
959  auto eltCoupled_en = EltCoupled.end();
960  size_type theIdEltTest = 0;
961 
962  for ( ; eltCoupled_it != eltCoupled_en ; ++eltCoupled_it, ++theIdEltTest )
963  {
964  if (eltCoupled_it->size()==0) continue;
965 
966  //init map
967  result_temp1_type mapTrial2idq;
968 
969  // search
970  auto quadPtTest_it = testEltToPtsQuad[theIdEltTest].begin();
971  auto const quadPtTest_en = testEltToPtsQuad[theIdEltTest].end();
972  for ( ; quadPtTest_it != quadPtTest_en ; ++quadPtTest_it )
973  {
974  // get test data
975  auto theIdq = quadPtTest_it->template get<0>();
976  auto theq = quadPtTest_it->template get<1>();
977  //std::cout << "\nidq modulo : " << theIdq%nQuadPtsInElt << " q " << theq << " elt " << theIdq/nQuadPtsInElt;
978  auto theNodeRefTest = quadPtTest_it->template get<2>();
979  // search in trial data
980  auto theRes = this->findQuadPt( *eltCoupled_it,trialEltToPtsQuad,theIdq );
981  // get trial data
982  auto const theIdEltTrial = theRes.template get<0>();
983  auto const theNodeRefTrial = theRes.template get<1>();
984  // add into the temporary map
985 
986  //std::cout << " theIdEltTrial " << theIdEltTrial <<"("<<meshTrial->numElements()<<")" << " theIdEltTest " << theIdEltTest <<"("<<meshTest->numElements()<<")"<<std::endl;
987  mapTrial2idq[theIdEltTrial].push_back( boost::make_tuple( theq,theIdq,theNodeRefTest,theNodeRefTrial ) );
988  }
989 
990  result_temp2_type mapiIdTrial2qAndPtRef;
991 
992  // build matrix_node of ptsRef (test and trial)
993  auto map_it = mapTrial2idq.begin();
994  auto map_en = mapTrial2idq.end();
995  for ( ; map_it != map_en ; ++map_it )
996  {
997  // init matrix node
998  auto nPtsRef = map_it->second.size();
999  //std::cout << "nPtsRef " << nPtsRef << std::endl;
1000  matrix_node_type ptsRefTest( mesh_test_type::nDim,nPtsRef );
1001  matrix_node_type ptsRefTrial( mesh_trial_type::nDim,nPtsRef );
1002  // init map between index node and index quad point
1003  std::vector<boost::tuple<size_type,size_type> > indexLocalToQuad( nPtsRef );
1004  // build
1005  auto idq_it = map_it->second.begin();
1006  auto idq_en = map_it->second.end();
1007  uint16_type cptIdq = 0;
1008 
1009  for ( ; idq_it != idq_en ; ++idq_it, ++cptIdq )
1010  {
1011  indexLocalToQuad[cptIdq] = boost::make_tuple( idq_it->template get<0>(),idq_it->template get<1>() );
1012  ublas::column( ptsRefTest, cptIdq ) = idq_it->template get<2>();
1013  ublas::column( ptsRefTrial, cptIdq ) = idq_it->template get<3>();
1014  }
1015 
1016  mapiIdTrial2qAndPtRef[map_it->first] = boost::make_tuple( indexLocalToQuad,ptsRefTest,ptsRefTrial );
1017  }
1018 
1019  // add to result container
1020  M_resBilinear.push_back( boost::make_tuple( theIdEltTest,mapiIdTrial2qAndPtRef ) );
1021  }
1022 
1023  } // update
1024 
1025 private :
1026 
1027  boost::tuple<size_type,node_type >
1028  findQuadPt( std::list<size_type> const & eltCoupled_it,
1029  std::vector<std::list<boost::tuple< size_type,size_type,node_type> > > const & trialPtQuadToElt,
1030  size_type idq )
1031  {
1032  size_type resId=0;
1033  ublas::vector<double> resNodeRef;
1034  bool find=false;
1035 
1036  auto eltTrial_it = eltCoupled_it.begin();
1037  auto eltTrial_en = eltCoupled_it.end();
1038 
1039  while ( !find && eltTrial_it != eltTrial_en )
1040  {
1041  auto quadPtTrial_it = trialPtQuadToElt[*eltTrial_it].begin();
1042  auto quadPtTrial_en = trialPtQuadToElt[*eltTrial_it].end();
1043 
1044  while ( !find && quadPtTrial_it != quadPtTrial_en )
1045  {
1046  if ( quadPtTrial_it->template get<0>() == idq )
1047  {
1048  find = true;
1049  resNodeRef=quadPtTrial_it->template get<2>();
1050  }
1051 
1052  else ++quadPtTrial_it;
1053  }
1054 
1055  if ( !find ) ++eltTrial_it;
1056 
1057  else resId=*eltTrial_it;
1058  }
1059 
1060  if ( !find ) std::cout << "BUG in findQuadPt " << std::endl;
1061 
1062  return boost::make_tuple( resId,resNodeRef );
1063  }
1064 
1065 
1066 #if FEELPP_EXPORT_QUADLOCALIZATION
1067  template <typename Mesh1Type,typename Mesh2Type>
1068  void
1069  exportQuadLocalization( boost::shared_ptr<Mesh1Type> meshTest,
1070  boost::shared_ptr<Mesh2Type> meshTrial,
1071  std::map<size_type,std::list<size_type> > const& mapBetweenMeshes_test,
1072  std::map<size_type,std::list<size_type> > const& mapBetweenMeshes_trial ) const
1073  {
1074  typedef FunctionSpace<Mesh1Type, bases<Lagrange<0, Scalar,Discontinuous> > > space_test_disc_type;
1075  auto spaceGraphProjTest = space_test_disc_type::New( meshTest );
1076  auto elemTest_itt = meshTest->beginElementWithProcessId( 0 );
1077  auto elemTest_ent = meshTest->endElementWithProcessId( 0 );
1078  auto graphProjTest = spaceGraphProjTest->element();graphProjTest.zero();
1079  for ( ; elemTest_itt != elemTest_ent; ++elemTest_itt )
1080  {
1081  if ( mapBetweenMeshes_test.find( elemTest_itt->id() ) != mapBetweenMeshes_test.end() )
1082  {
1083  if ( mapBetweenMeshes_test.find( elemTest_itt->id() )->second.size()>0 )
1084  {
1085  auto element_dof1 = spaceGraphProjTest->dof()->getIndices( elemTest_itt->id() );
1086  graphProjTest( element_dof1[0] ) = 1;
1087  }
1088  }
1089  }
1090  //auto exporterTest = boost::shared_ptr<Feel::Exporter<Mesh1Type> >( Feel::Exporter<Mesh1Type>::New( "ensight", "ExportQuadLocalizationTest" ) );
1091  auto exporterTest = Feel::Exporter<Mesh1Type>::New( "ensight", "ExportQuadLocalizationTest" );
1092  exporterTest->step( 0 )->setMesh( graphProjTest.mesh() );
1093  exporterTest->step( 0 )->add( "quadLocalizationProjTest", graphProjTest );
1094  exporterTest->save();
1095 
1096  typedef FunctionSpace<Mesh2Type, bases<Lagrange<0, Scalar,Discontinuous> > > space_trial_disc_type;
1097  auto spaceGraphProjTrial = space_trial_disc_type::New( meshTrial );
1098  auto elemTrial_itt = meshTrial->beginElementWithProcessId( 0 );
1099  auto elemTrial_ent = meshTrial->endElementWithProcessId( 0 );
1100  auto graphProjTrial = spaceGraphProjTrial->element();graphProjTrial.zero();
1101  for ( ; elemTrial_itt != elemTrial_ent; ++elemTrial_itt )
1102  {
1103  if ( mapBetweenMeshes_trial.find( elemTrial_itt->id() ) != mapBetweenMeshes_trial.end() )
1104  {
1105  if ( mapBetweenMeshes_trial.find( elemTrial_itt->id() )->second.size()>0 )
1106  {
1107  auto element_dof2 = spaceGraphProjTrial->dof()->getIndices( elemTrial_itt->id() );
1108  graphProjTrial( element_dof2[0] ) = 1;
1109  }
1110  }
1111  }
1112  //auto exporterTrial = boost::shared_ptr<Feel::Exporter<Mesh2Type> >( Feel::Exporter<Mesh2Type>::New( "ensight", "ExportQuadLocalizationTrial" ) );
1113  auto exporterTrial = Feel::Exporter<Mesh2Type>::New( "ensight", "ExportQuadLocalizationTrial" );
1114  exporterTrial->step( 0 )->setMesh( graphProjTrial.mesh() );
1115  exporterTrial->step( 0 )->add( "quadLocalizationProjTrial", graphProjTrial );
1116  exporterTrial->save();
1117 
1118  }
1119 #endif
1120 
1121 public :
1122 
1123 
1124 
1125 template<typename FE1,typename FE2,typename ElemContType>
1126 struct bilinearformContext
1127 {
1128  typedef Expr expression_type;
1129 
1130  // typedef on form (trial and test):
1131  typedef vf::detail::BilinearForm<FE1,FE2,ElemContType> FormType;
1132  // test
1133  typedef typename FormType::gm_1_type gm_formTest_type;
1134  typedef typename FormType::mesh_element_1_type geoelement_formTest_type;
1135  typedef typename gm_formTest_type::template Context<expression_type::context|vm::POINT,geoelement_formTest_type> gmc_formTest_type;
1136  typedef boost::shared_ptr<gmc_formTest_type> gmc_formTest_ptrtype;
1137  typedef typename gm_formTest_type::precompute_type pc_formTest_type;
1138  typedef typename gm_formTest_type::precompute_ptrtype pc_formTest_ptrtype;
1139  //typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTest_ptrtype> > map_gmc_formTest_type;
1140  // trial
1141  typedef typename FormType::gm_2_type gm_formTrial_type;
1142  typedef typename FormType::mesh_element_2_type geoelement_formTrial_type;
1143  typedef typename gm_formTrial_type::template Context<expression_type::context|vm::POINT,geoelement_formTrial_type> gmc_formTrial_type;
1144  typedef boost::shared_ptr<gmc_formTrial_type> gmc_formTrial_ptrtype;
1145  typedef typename gm_formTrial_type::precompute_type pc_formTrial_type;
1146  typedef typename gm_formTrial_type::precompute_ptrtype pc_formTrial_ptrtype;
1147  //typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTrial_ptrtype> > map_gmc_formTrial_type;
1148 
1149  typedef std::list<boost::tuple< std::vector<boost::tuple<size_type,size_type> >,
1150  gmc_ptrtype,
1151  gmc_formTest_ptrtype,
1152  gmc_formTrial_ptrtype > > return_loc_type;
1153  typedef std::list<return_loc_type> return_type;
1154 
1155 };
1156 
1157 
1158 template<typename FE1,typename FE2,typename ElemContType>
1159 void
1160 precompute(vf::detail::BilinearForm<FE1,FE2,ElemContType>const& __form)
1161 {
1162  typedef typename bilinearformContext<FE1,FE2,ElemContType>::gmc_formTest_type gmc_formTest_type;
1163  typedef typename bilinearformContext<FE1,FE2,ElemContType>::gmc_formTest_ptrtype gmc_formTest_ptrtype;
1164  typedef typename bilinearformContext<FE1,FE2,ElemContType>::pc_formTest_type pc_formTest_type;
1165  typedef typename bilinearformContext<FE1,FE2,ElemContType>::pc_formTest_ptrtype pc_formTest_ptrtype;
1166 
1167  typedef typename bilinearformContext<FE1,FE2,ElemContType>::gmc_formTrial_type gmc_formTrial_type;
1168  typedef typename bilinearformContext<FE1,FE2,ElemContType>::gmc_formTrial_ptrtype gmc_formTrial_ptrtype;
1169  typedef typename bilinearformContext<FE1,FE2,ElemContType>::pc_formTrial_type pc_formTrial_type;
1170  typedef typename bilinearformContext<FE1,FE2,ElemContType>::pc_formTrial_ptrtype pc_formTrial_ptrtype;
1171 
1172  auto meshTrial = __form.trialSpace()->mesh();
1173  auto meshTest = __form.testSpace()->mesh();
1174 
1175  typename bilinearformContext<FE1,FE2,ElemContType>::return_type theres;
1176 
1177  auto res_it = this->result().begin();
1178  auto const res_en = this->result().end();
1179  for ( ; res_it != res_en ; ++res_it )
1180  {
1181  auto const idEltTest = res_it->template get<0>();
1182  auto const& map = res_it->template get<1>();
1183  auto map_it = map.begin();
1184  auto const map_en = map.end();
1185  for ( ; map_it != map_en ; ++map_it )
1186  {
1187  auto const idEltTrial = map_it->first;
1188  auto const& eltTrial = meshTrial->element( idEltTrial );
1189  auto const& eltTest = meshTest->element( idEltTest );
1190 
1191  auto const& ptRefTest = map_it->second.template get<1>();
1192  auto const& ptRefTrial = map_it->second.template get<2>();
1193  auto const& themapQuad = map_it->second.template get<0>();
1194  auto vec_gmcExpr = this->getUsableDataInFormContext( themapQuad,ptRefTest,ptRefTrial );
1195  auto gmcExpr_it = vec_gmcExpr.begin();
1196  auto const gmcExpr_en = vec_gmcExpr.end();
1197  typename bilinearformContext<FE1,FE2,ElemContType>::return_loc_type theresloc;
1198  for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
1199  {
1200  pc_formTest_ptrtype geopcFormTest( new pc_formTest_type( __form.gm(), gmcExpr_it->template get<2>()/*__form.testSpace()->fe()->points()*/ ) );
1201  gmc_formTest_ptrtype gmcFormTest( new gmc_formTest_type( __form.gm(), eltTest /*__form.testSpace()->mesh()->element( 0 )*/, geopcFormTest ) );
1202 
1203  pc_formTrial_ptrtype geopcFormTrial( new pc_formTrial_type( __form.gmTrial(), gmcExpr_it->template get<3>() /* __form.trialSpace()->fe()->points()*/ ) );
1204  gmc_formTrial_ptrtype gmcFormTrial( new gmc_formTrial_type( __form.gmTrial(), eltTrial/*__form.trialSpace()->mesh()->element( 0 )*/, geopcFormTrial ) );
1205 
1206  theresloc.push_back(boost::make_tuple(gmcExpr_it->template get<0>(),gmcExpr_it->template get<1>(),gmcFormTest,gmcFormTrial));
1207  }
1208  theres.push_back(theresloc);
1209  }
1210  }
1211 
1212  M_precompute = theres;
1213  M_hasPrecompute=true;
1214 
1215 }
1216 
1217 
1218 template<typename FE,typename VectorType,typename ElemContType>
1219 struct linearformContext
1220 {
1221  typedef Expr expression_type;
1222 
1223  //typedef on form test :
1224  typedef vf::detail::LinearForm<FE,VectorType,ElemContType> FormType;
1225  typedef typename FormType::gm_type gm_form_type;
1226  typedef typename FormType::mesh_test_element_type geoelement_form_type;
1227  typedef typename gm_form_type::template Context<expression_type::context|vm::POINT,geoelement_form_type> gmc_form_type;
1228  typedef boost::shared_ptr<gmc_form_type> gmc_form_ptrtype;
1229  typedef typename gm_form_type::precompute_type pc_form_type;
1230  typedef typename gm_form_type::precompute_ptrtype pc_form_ptrtype;
1231 
1232  typedef std::list<boost::tuple< std::vector<boost::tuple<size_type,size_type> >,
1233  gmc_ptrtype,
1234  gmc_form_ptrtype > > return_loc_type;
1235  typedef std::list<return_loc_type> return_type;
1236 
1237 };
1238 
1239 template<typename FE,typename VectorType,typename ElemContType>
1240 void
1241 precompute(vf::detail::LinearForm<FE,VectorType,ElemContType> const& __form)
1242 {
1243 
1244  typedef typename linearformContext<FE,VectorType,ElemContType>::gmc_form_type gmc_form_type;
1245  typedef typename linearformContext<FE,VectorType,ElemContType>::gmc_form_ptrtype gmc_form_ptrtype;
1246  typedef typename linearformContext<FE,VectorType,ElemContType>::pc_form_type pc_form_type;
1247  typedef typename linearformContext<FE,VectorType,ElemContType>::pc_form_ptrtype pc_form_ptrtype;
1248 
1249  auto meshTest = __form.testSpace()->mesh();
1250 
1251  typename linearformContext<FE,VectorType,ElemContType>::return_type theres;
1252 
1253  auto res_it = this->resultLinear().begin();
1254  auto const res_en = this->resultLinear().end();
1255  for ( ; res_it != res_en ; ++res_it )
1256  {
1257  auto const idEltTest = res_it->template get<0>();
1258  auto const& eltTest = meshTest->element( idEltTest );
1259 
1260  auto ptRefTest = res_it->template get<2>();
1261  auto themapQuad = res_it->template get<1>();
1262 
1263  auto vec_gmcExpr = this->getUsableDataInFormContext( themapQuad,ptRefTest );
1264  auto gmcExpr_it = vec_gmcExpr.begin();
1265  auto const gmcExpr_en = vec_gmcExpr.end();
1266  typename linearformContext<FE,VectorType,ElemContType>::return_loc_type theresloc;
1267  for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
1268  {
1269 
1270  pc_form_ptrtype geopcForm( new pc_form_type( __form.gm(), gmcExpr_it->template get<2>() /*this->im().points()*/ ) );
1271  gmc_form_ptrtype gmcForm( new gmc_form_type( __form.gm(), eltTest, geopcForm ) );
1272  theresloc.push_back(boost::make_tuple(gmcExpr_it->template get<0>(),gmcExpr_it->template get<1>(),gmcForm));
1273  }
1274  theres.push_back(theresloc);
1275 
1276  }
1277 
1278  M_precompute = theres;
1279  M_hasPrecompute=true;
1280 
1281 }
1282 
1283 
1284 template<typename FE1,typename FE2,typename ElemContType>
1285 typename bilinearformContext<FE1,FE2,ElemContType>::return_type const&
1286 getPrecompute(vf::detail::BilinearForm<FE1,FE2,ElemContType>const& __form) const
1287 {
1288  return boost::any_cast<typename bilinearformContext<FE1,FE2,ElemContType>::return_type const&>( M_precompute);
1289 }
1290 
1291 template<typename FE,typename VectorType,typename ElemContType>
1292 typename linearformContext<FE,VectorType,ElemContType>::return_type const&
1293 getPrecompute(vf::detail::LinearForm<FE,VectorType,ElemContType> const& __form) const
1294 {
1295  return boost::any_cast<typename linearformContext<FE,VectorType,ElemContType>::return_type const&>( M_precompute);
1296 }
1297 
1298 private :
1299 
1300  std::list<range_iterator> M_listRange;
1301 
1302  element_iterator_type M_eltbegin;
1303  element_iterator_type M_eltend;
1304  mutable im_type M_im;
1305  QuadMapped<im_type> M_qm;
1306  permutation_points_type M_ppts;
1307 
1308  result_container_type M_resBilinear;
1309 
1310  result_container_linear_type M_resLinear;
1311 
1312  bool M_hasPrecompute;
1313  boost::any M_precompute;
1314 
1315 
1316 }; // QuadPtLocalization
1317 
1318 
1319 
1320 
1321 namespace detail
1322 {
1323 template <typename RangeType>
1324 struct quadptlocrangetype
1325 {
1326  typedef typename mpl::if_< boost::is_std_list<RangeType>,
1327  mpl::identity<RangeType>,
1328  mpl::identity<std::list<RangeType> > >::type::type::value_type type;
1329 };
1330 }
1331 
1332 
1333 template<typename MeshTestType, typename MeshTrialType, typename IteratorRange,typename Im,typename Expr>
1334 boost::shared_ptr<QuadPtLocalization<typename Feel::detail::quadptlocrangetype<IteratorRange>::type,Im,Expr> >
1335 quadPtLocPtr( boost::shared_ptr<MeshTestType> meshTest, boost::shared_ptr<MeshTrialType> meshTrial, IteratorRange const& elts,Im const& im,Expr const& expr )
1336 {
1337  typedef QuadPtLocalization<typename Feel::detail::quadptlocrangetype<IteratorRange>::type,Im,Expr> quadptloc_type;
1338  boost::shared_ptr<quadptloc_type> res(new quadptloc_type(elts) );
1339  res->update( meshTest,meshTrial );
1340  return res;
1341 }
1342 
1343 template<typename MeshTestType, typename IteratorRange,typename Im,typename Expr>
1344 boost::shared_ptr<QuadPtLocalization<typename Feel::detail::quadptlocrangetype<IteratorRange>::type,Im,Expr> >
1345 quadPtLocPtr( boost::shared_ptr<MeshTestType> meshTest, IteratorRange const& elts,Im const& im,Expr const& expr )
1346 {
1347  typedef QuadPtLocalization<typename Feel::detail::quadptlocrangetype<IteratorRange>::type,Im,Expr> quadptloc_type;
1348  boost::shared_ptr<quadptloc_type> res(new quadptloc_type(elts) );
1349  res->update( meshTest );
1350  return res;
1351 }
1352 
1353 
1354 
1355 } // Feel
1356 #endif
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
static boost::shared_ptr< Exporter< MeshType, N > > New(std::string const &exportername, std::string prefix=Environment::about().appName(), WorldComm const &worldComm=Environment::worldComm())
Definition: exporterimpl.hpp:127
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319

Generated on Sun Dec 22 2013 13:11:12 for Feel++ by doxygen 1.8.5