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
meshimpl.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: 2007-08-07
7 
8  Copyright (C) 2007-2010 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 */
29 #ifndef __MESHIMPL_HPP
30 #define __MESHIMPL_HPP 1
31 
32 #include <boost/preprocessor/comparison/greater_equal.hpp>
37 #include <feel/feeldiscr/mesh.hpp>
38 #include <feel/feelalg/glas.hpp>
41 
42 namespace Feel
43 {
44 template<typename Shape, typename T, int Tag>
45 Mesh<Shape, T, Tag>::Mesh( WorldComm const& worldComm )
46  :
47  super(worldComm),
48  M_gm( new gm_type ),
49  M_gm1( new gm1_type ),
50  M_meas( 0 ),
51  M_measbdy( 0 ),
52  //M_part(),
53  M_tool_localization( new Localization() )
54 {
55  VLOG(2) << "[Mesh] constructor called\n";
56 }
57 template<typename Shape, typename T, int Tag>
58 void
59 Mesh<Shape, T, Tag>::partition ( const uint16_type n_parts )
60 {
61  //M_part->partition( *this, n_parts );
62 }
63 
64 template<typename Shape, typename T, int Tag>
65 flag_type
66 Mesh<Shape, T, Tag>::markerId ( boost::any const& __marker )
67 {
68  flag_type theflag = -1;
69  if ( boost::any_cast<flag_type>( &__marker ) )
70  {
71  theflag = boost::any_cast<flag_type>( __marker);
72  }
73  else if ( boost::any_cast<int>( &__marker ) )
74  {
75  theflag = boost::any_cast<int>( __marker);
76  }
77  else if ( boost::any_cast<size_type>( &__marker ) )
78  {
79  theflag = boost::any_cast<size_type>( __marker);
80  }
81  else if ( boost::any_cast<uint16_type>( &__marker ) )
82  {
83  theflag = boost::any_cast<uint16_type>( __marker);
84  }
85  else if ( boost::any_cast<std::string>( &__marker ) )
86  {
87  theflag = this->markerName( boost::any_cast<std::string>( __marker) );
88  }
89  else
90  CHECK( theflag != -1 ) << "invalid flag type\n";
91  return theflag;
92 }
93 
94 template<typename Shape, typename T, int Tag>
95 void
97 {
98  VLOG(2) << "component MESH_RENUMBER: " << this->components().test( MESH_RENUMBER ) << "\n";
99  VLOG(2) << "component MESH_UPDATE_EDGES: " << this->components().test( MESH_UPDATE_EDGES ) << "\n";
100  VLOG(2) << "component MESH_UPDATE_FACES: " << this->components().test( MESH_UPDATE_FACES ) << "\n";
101  VLOG(2) << "component MESH_PARTITION: " << this->components().test( MESH_PARTITION ) << "\n";
102 
103  if ( this->numElements() == 0 )
104  {
105  VLOG(2) << "No elements in Mesh? (with process rank " << this->worldComm().rank() <<")\n";
106  if ( this->worldComm().localSize()==1 ) return;
107  }
108 
109  boost::timer ti;
110 
111  VLOG(2) << "is already updated? : " << this->isUpdatedForUse() << "\n";
112  if ( !this->isUpdatedForUse() )
113  {
114  if ( this->components().test( MESH_RENUMBER ) )
115  {
116 
117  //this->renumber();
118  VLOG(2) << "[Mesh::updateForUse] renumber : " << ti.elapsed() << "\n";
119  }
120 
121  //
122  // compute the Adjacency graph
123  //
124  ti.restart();
125  VLOG(2) << "Compute adjacency graph\n";
126  element_iterator iv, en;
127  std::map<size_type,boost::tuple<size_type, uint16_type, size_type> > f2e;
128 
129  std::map<std::set<int>, size_type > _faces;
130  typename std::map<std::set<int>, size_type >::iterator _faceit;
131  int next_face = 0;
132 
133  boost::tie( iv, en ) = this->elementsRange();
134 
135  for ( ; iv != en; ++iv )
136  {
137  for ( size_type j = 0; j < this->numLocalFaces(); j++ )
138  {
139  std::set<int> s;
140 
141  for ( int f = 0; f < face_type::numVertices; ++f )
142  {
143  if ( nDim == 1 )
144  s.insert( iv->point( j ).id() );
145 
146  else
147  s.insert( iv->point( iv->fToP( j, f ) ).id() );
148  }
149 
150  bool faceinserted = false;
151  boost::tie( _faceit, faceinserted ) = _faces.insert( std::make_pair( s, next_face ) );
152 
153  if ( faceinserted )
154  ++next_face;
155 
156  DVLOG(2) << "------------------------------------------------------------\n";
157  DVLOG(2) << "Element id: " << iv->id() << " local face id: " << j << " process id:" << iv->processId() << "\n";
158 
159  if ( faceinserted )
160  {
161  DVLOG(2) << " new face " << _faceit->second << " is now in store with elt "
162  << f2e[_faceit->second].template get<0>() << " and local face id "
163  << f2e[_faceit->second].template get<1>() << "\n";
164 
165 
166  f2e[_faceit->second].template get<0>() = iv->id();
167  f2e[_faceit->second].template get<1>() = j;
168  M_e2f[std::make_pair(iv->id(),j)]=boost::make_tuple( _faceit->second, invalid_size_type_value );
169  }
170 
171  else // already stored
172  {
173  DVLOG(2) << "old face " << _faceit->second << " was already in store with elt " << f2e[_faceit->second].template get<0>() << " and local face id " << f2e[_faceit->second].template get<1>() << "\n";
174 
175  f2e[_faceit->second].template get<2>() = iv->id();
176  M_e2f[std::make_pair(iv->id(),j)]=boost::make_tuple( _faceit->second, f2e[_faceit->second].template get<0>() );
177  M_e2f[std::make_pair(f2e[_faceit->second].template get<0>(), f2e[_faceit->second].template get<1>())] =
178  boost::make_tuple( _faceit->second, iv->id() );
179  }
180 
181 
182  } // local face
183  } // element loop
184 
185  VLOG(2) << "Compute adjacency graph done in " << ti.elapsed() << "\n";
186 #if 0
187 
188  // partition mesh
189  if ( this->components().test( MESH_PARTITION ) || ( this->worldComm().localSize() > 1 ) )
190  {
191  boost::timer ti1;
192  this->partition();
193  VLOG(2) << "[Mesh::updateForUse] partition time : " << ti1.elapsed() << "\n";
194  }
195 
196 #endif
197 
198  if ( this->components().test( MESH_UPDATE_FACES ) )
199  {
200  ti.restart();
201  // update connectivities of entities of co dimension 1
202  this->updateEntitiesCoDimensionOne();
203  // update permutation of entities of co-dimension 1
204  this->updateEntitiesCoDimensionOnePermutation();
205  // update in ghost cells of entities of co-dimension 1
206  //if (this->worldComm().localSize()>1)
207  // this->updateEntitiesCoDimensionOneGhostCell();
208 
209  VLOG(2) << "[Mesh::updateForUse] update entities of codimension 1 : " << ti.elapsed() << "\n";
210 
211  }
212 
213  if ( this->components().test( MESH_UPDATE_EDGES ) )
214  {
215  ti.restart();
216  // update connectivities of entities of co dimension 2
217  // (edges in 3D)
218  this->updateEntitiesCoDimensionTwo();
219  VLOG(2) << "[Mesh::updateForUse] update edges : " << ti.elapsed() << "\n";
220  }
221 
222  if ( this->components().test( MESH_UPDATE_FACES ) ||
223  this->components().test( MESH_UPDATE_EDGES )
224  )
225  {
226  updateOnBoundary();
227  }
228  this->setUpdatedForUse( true );
229  }
230 
231  {
232  element_iterator iv, en;
233  if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
234  {
235  boost::tie( iv, en ) = this->elementsRange();
236  for ( ; iv != en; ++iv )
237  {
238  this->elements().modify( iv, typename super_elements::ElementConnectPointToElement() );
239  }
240  }
241 
242  boost::tie( iv, en ) = this->elementsRange();
243  auto pc = M_gm->preCompute( M_gm, M_gm->referenceConvex().vertices() );
244  auto pcf = M_gm->preComputeOnFaces( M_gm, M_gm->referenceConvex().barycenterFaces() );
245  M_meas = 0;
246  M_measbdy = 0;
247 
248  for ( ; iv != en; ++iv )
249  {
250  this->elements().modify( iv,
251  [=,&pc,&pcf]( element_type& e )
252  {
253  for ( int i = 0; i < e.numPoints; ++i )
254  e.point( i ).addElement( e.id() );
255  e.setMeshAndGm( this, M_gm, M_gm1 );
256  e.updateWithPc(pc, boost::ref( pcf) );
257  } );
258 #if 0
259  lambda::bind( &element_type::setMeshAndGm,
260  lambda::_1,
261  this, M_gm, M_gm1 ) );
262 
263  this->elements().modify( iv,
264  lambda::bind( &element_type::updateWithPc,
265  lambda::_1, pc, boost::ref( pcf ) ) );
266 #endif
267 
268  M_meas += iv->measure();
269  auto _faces = iv->faces();
270 
271  if ( nDim == 1 )
272  M_measbdy = 0;
273  else
274  for ( ; _faces.first != _faces.second; ++_faces.first )
275  if ( ( *_faces.first ) && ( *_faces.first )->isOnBoundary() )
276  M_measbdy += ( *_faces.first )->measure();
277  }
278 
279  // now that all elements have been updated, build inter element
280  // data such as the measure of point element neighbors
281  boost::tie( iv, en ) = this->elementsRange();
282 
283  for ( ; iv != en; ++iv )
284  {
285  value_type meas = 0;
286  for( auto _elt: iv->pointElementNeighborIds() )
287  {
288  if ( this->hasElement( _elt ) )
289  meas += this->element( _elt ).measure();
290  }
291  this->elements().modify( iv, [meas]( element_type& e ){ e.setMeasurePointElementNeighbors( meas ); } );
292  }
293 
294  for ( auto itf = this->beginFace(), ite = this->endFace(); itf != ite; ++ itf )
295  {
296  this->faces().modify( itf,[this]( face_type& f ) { f.setMesh( this ); } );
297  }
298  }
299  //std::cout<<"this->worldComm().localSize()= "<< this->worldComm().localSize() << std::endl;
300 #if defined(FEELPP_ENABLE_MPI_MODE)
301 
302  if ( this->components().test( MESH_UPDATE_FACES ) && this->worldComm().localSize()>1 )
303  {
304  this->updateEntitiesCoDimensionOneGhostCell();
305  }
306 #endif
307 
308  if ( this->components().test( MESH_PROPAGATE_MARKERS ) )
309  propagateMarkers(mpl::int_<nDim>() );
310 
311  // check mesh connectivity
312  this->check();
313  //std::cout<<"pass hier\n";
314 
315  M_gm->initCache( this );
316  M_gm1->initCache( this );
317 
318  M_tool_localization->setMesh( this->shared_from_this(),false );
319 
320 
321 
322  VLOG(2) << "[Mesh::updateForUse] total time : " << ti.elapsed() << "\n";
323 }
324 template<typename Shape, typename T, int Tag>
326 Mesh<Shape, T, Tag>::operator+=( self_type const& m )
327 {
328  std::map<std::vector<double>, size_type> mapDel;
329  for( auto it = m.beginPoint(), en = m.endPoint();
330  it != en;
331  ++it )
332  {
333  bool found = false;
334  if( it->isOnBoundary() )
335  {
336  for( auto pit = this->beginPoint(), pen = this->endPoint();
337  pit != pen;
338  ++pit )
339  {
340  if ( pit->isOnBoundary() )
341  {
342  if ( ublas::norm_2( it->node() - pit->node() ) < 1e-10 )
343  {
344  found = true;
345  std::vector<double> ublasCopy(nDim);
346  std::copy(it->node().begin(), it->node().end(), ublasCopy.begin());
347  // std::cout << "Attention ! Le point " << it->id() << " doit etre remplace par le point " << pit->id() << " (" << it->node() << " ~= " << pit->node() << ")" << std::endl;
348  mapDel.insert(std::pair<std::vector<double>, size_type>(ublasCopy, pit->id()));
349  break;
350  }
351  }
352  }
353  }
354  if(!found)
355  this->addPoint( *it );
356  }
357  for( auto it = m.beginFace(), en = m.endFace();it != en;++it )
358  {
359  // eventually need work
360  if ( !it->isOnBoundary() )
361  {
362  face_type f = *it;
363  f.disconnect();
364  this->addFace( f );
365  }
366  }
367  for( auto it = m.beginElement(), en = m.endElement();it != en;++it )
368  {
369  element_type e = *it;
370  for( uint16_type p = 0; p < e.numPoints; ++p ) {
371  if ( e.point(p).isOnBoundary() ) {
372  std::vector<double> ublasCopy(nDim);
373  std::copy(e.point(p).node().begin(), e.point(p).node().end(), ublasCopy.begin());
374  auto itMap = mapDel.find(ublasCopy);
375  if (itMap != mapDel.end() )
376  {
377  // std::cout << "Attention ! Le point " << it->point(p).id() << " va etre remplace par le point " << this->point(itMap->second).id() << std::endl;
378  e.setPoint(p, this->point(itMap->second));
379  }
380  }
381  }
382  this->addElement( e );
383  }
384  this->setUpdatedForUse( false );
385  this->updateForUse();
386  return *this;
387 }
388 template<typename Shape, typename T, int Tag>
389 void
390 Mesh<Shape, T, Tag>::propagateMarkers( mpl::int_<2> )
391 {
392  // propagate top-down marker from face if edge has not been marked
393  std::for_each( this->beginFace(), this->endFace(),
394  [this]( face_type const& f )
395  {
396  if ( f.marker().isOff() )
397  return;
398 
399  // update points
400  for( int i = 0; i < face_type::numPoints; ++i )
401  {
402  if ( f.point( i ).marker().isOff() )
403  {
404  // inherit marker from edge
405  this->points().modify( this->points().iterator_to( f.point(i) ),
406  [&f] ( point_type& p )
407  {
408  p.setMarker( f.marker().value() );
409  } );
410 
411  }
412  }
413  } );
414 }
415 
416 template<typename Shape, typename T, int Tag>
417 void
418 Mesh<Shape, T, Tag>::propagateMarkers( mpl::int_<3> )
419 {
420  // first propagate top-down marker from edges if points have not been marked
421  std::for_each( this->beginEdge(), this->endEdge(),
422  [this]( edge_type const& e )
423  {
424  if ( e.marker().isOff() )
425  return;
426 
427  for( int i = 0; i < edge_type::numPoints; ++i )
428  {
429  if ( e.point( i ).marker().isOff() )
430  {
431  // inherit marker from edge
432  this->points().modify( this->points().iterator_to( e.point(i) ),
433  [&e] ( point_type& p )
434  {
435  p.setMarker( e.marker().value() );
436  } );
437 
438  }
439  } } );
440 
441  // then propagate top-down marker from face if edge has not been marked
442  std::for_each( this->beginFace(), this->endFace(),
443  [this]( face_type const& f )
444  {
445  if ( f.marker().isOff() )
446  return;
447 
448  // update points
449  for( int i = 0; i < face_type::numPoints; ++i )
450  {
451  if ( f.point( i ).marker().isOff() )
452  {
453  // inherit marker from edge
454  this->points().modify( this->points().iterator_to( f.point(i) ),
455  [&f] ( point_type& p )
456  {
457  p.setMarker( f.marker().value() );
458  } );
459 
460  }
461  }
462  // update edges
463  for( int i = 0; i < face_type::numEdges; ++i )
464  {
465  if ( f.edge( i ).marker().isOff() )
466  {
467  // inherit marker from edge
468  this->edges().modify( this->edges().iterator_to( f.edge(i) ),
469  [&f] ( edge_type& e )
470  {
471  e.setMarker( f.marker().value() );
472  } );
473 
474  }
475  }
476  } );
477 }
478 
479 template<typename Shape, typename T, int Tag>
480 void
481 Mesh<Shape, T, Tag>::renumber( mpl::bool_<true> )
482 {
483  size_type next_free_node = 0;
484 
485  // map old/new ids
486  std::vector<size_type> node_map( this->numPoints(), invalid_size_type_value );
487  //std::map<size_type> node_map( this->numPoints(), invalid_size_type_value );
488  //std::vector<point_type> pt_map( this->numPoints() );
489  typedef std::map<size_type,point_type> ptmap_type;
490  ptmap_type pt_map;
491 
492  typedef typename ptmap_type::iterator ptmap_iterator;
493  typedef std::vector<size_type>::iterator nm_iterator;
494 
495  // first collect all point id that will be swapped in a dictionary
496  for ( element_const_iterator elt = this->beginElement();
497  elt != this->endElement(); ++elt )
498  {
499  element_type const& __element = *elt;
500 
501  DVLOG(2) << "mesh::renumber] element id " << __element.id() << " proc " << __element.processId() << "\n";
502 #if !defined(NDEBUG)
503  for ( int i = 0; i < __element.nPoints(); ++i )
504  {
505  DVLOG(2) << "point id = " << __element.point( i ).id() << "\n";
506  }
507 #endif
508 
509  // renumber the nodes of the element
510  for ( int i = 0; i < __element.nPoints(); ++i )
511  {
512  point_type __pt = __element.point( i );
513  size_type __true_id = __pt.id();
514  size_type __id = __true_id;
515 
516  // did we already change the id? if yes get the
517  // id it was changed to
518  if ( node_map[ __true_id ] != invalid_size_type_value )
519  {
520  __id = node_map[__true_id];
521 
522  DVLOG_IF(2,(__id >= next_free_node )) << "next_free_node = " << next_free_node
523  << " point id = " << __id << "\n";
524  }
525 
526  // don't renumber if already done
527  if ( __id == next_free_node )
528  {
529  node_map[ __true_id ] = next_free_node;
530  pt_map[next_free_node] = __pt;
531  pt_map[next_free_node].setId( next_free_node );
532  ++next_free_node;
533  }
534 
535  else if ( __id > next_free_node )
536  {
537  // first check if next_free_node has been
538  // already taken by another point, if so
539  // then affect this point id to __id
540  nm_iterator nm_it = std::find( node_map.begin(),
541  node_map.end(),
542  next_free_node );
543 
544  if ( nm_it != node_map.end() )
545  *nm_it = __id;
546 
547  else
548  node_map[ next_free_node ] = __id;
549 
550  // must use \p __true_id here as it is the
551  // original point id we started with
552  node_map[ __true_id ] = next_free_node;
553  pt_map[next_free_node] = __pt;
554  pt_map[next_free_node].setId( next_free_node );
555 
556  DVLOG(2) << "next_free_node = " << next_free_node
557  << " swapping point " << __true_id
558  << " and point " << next_free_node << "\n";
559 
560  // next free node
561  ++next_free_node;
562  } // if
563  } // for
564 
565  //__element.setId(
566  }
567 
568  VLOG(2) << "[mesh::renumber] done collecting ids\n";
569 #if !defined(NDEBUG)
570  std::vector<size_type> check_id( node_map );
571  std::unique( check_id.begin(), check_id.end() );
572 
573  FEELPP_ASSERT( check_id.size() == node_map.size() )( node_map.size() )( check_id.size() ).error( "all ids must be unique" );
574 
575  // in parallel this will generate a lot of warning since all the points are
576  //loaded into the mesh data structure including the ones not belonging to
577  //the current processor. Do the test only in sequential
578  if ( Environment::numberOfProcessors() == 1 )
579  {
580  FEELPP_ASSERT( std::find( node_map.begin(),node_map.end(), invalid_size_type_value ) == node_map.end() ).warn("invalid size_type value found as id " );
581  }
582 
583 #endif /* NDEBUG */
584 
585  // now we can replace the points id
586  //for( size_type p = 0; p < pt_map.size(); ++p )
587  ptmap_iterator ptmapit = pt_map.begin();
588  ptmap_iterator ptmapen = pt_map.end();
589 
590  for ( ; ptmapit != ptmapen; ++ptmapit )
591  {
592  point_iterator ptit = this->points().find( ptmapit->first );
593 #if !defined( NDEBUG )
594  VLOG(2) << "[mesh::replace] replacing point " << ptit->id()
595  << " with " << ptmapit->second.id() << "\n";
596 #endif
597  bool __rep1 = this->points().replace( ptit, ptmapit->second );
598  FEELPP_ASSERT( __rep1 )( __rep1 )( ptit->id() )( ptmapit->second.id() ) .warn( "invalid point replacement" );
599 
600  }
601 
602  VLOG(2) << "[mesh::renumber] done replace point ids\n";
603 
604  for ( element_iterator elt = this->beginElement();
605  elt != this->endElement(); ++elt )
606  {
607  element_type __element = *elt;
608  DVLOG(2) << "mesh::renumber] element id " << __element.id() << " proc " << __element.processId() << "\n";
609 
610  // renumber the nodes of the element
611  for ( int i = 0; i < __element.nPoints(); ++i )
612  {
613 
614  size_type __true_id =__element.point( i ).id();
615  auto& new_pt = this->point( node_map[__true_id] );
616  this->elements().modify( elt,[i,&new_pt]( element_type& e ) { e.setPoint( i, new_pt ); } );
617 
618  DVLOG(2) << "point id = " << __true_id << " node map " << node_map[__true_id] << " "
619  << "new point id = " << elt->point( i ).id() << "\n";
620  }
621  }
622 
623  VLOG(2) << "[mesh::renumber] done replace point ids in elements\n";
624 
625  for ( face_iterator elt = this->beginFace();
626  elt != this->endFace(); ++elt )
627  {
628  face_type __face = *elt;
629  DVLOG(2) << "face id: " << __face.id()
630  << " marker: " << __face.marker() << "\n";
631 
632  // renumber the nodes of the face
633  for ( int i = 0; i < __face.nPoints(); ++i )
634  {
635  size_type __true_id =__face.point( i ).id();
636  this->faces().modify( elt,
637  lambda::bind( &face_type::setPoint,
638  lambda::_1,
639  lambda::constant( i ),
640  boost::cref( this->point( node_map[__true_id] ) ) ) );
641  face_type __face2 = *elt;
642 
643  DVLOG(2) << "id1= " << __face2.point( 0 ).id() << " id2= " << __face2.point( 1 ).id()<< "\n";
644  DVLOG(2) << "point lid = " << i << " id = " << __true_id
645  << " nid = " << this->point( node_map[__true_id] ).id()
646  << " new point id = " << elt->point( i ).id() << "\n";
647  }
648  }
649  renumber( node_map, mpl::int_<nDim>() );
650 
651  if ( Shape == SHAPE_TETRA && nOrder==1 )
652  {
653  localrenumber();
654 
655 #if !defined(NDEBUG)
656  checkLocalPermutation( mpl::bool_< ( Shape == SHAPE_TETRA ) >() );
657 #endif
658  }
659 
660  // should we renumber also the faces and elements ?
661 
662 
663 }
664 
665 template<typename Shape, typename T, int Tag>
666 void
667 Mesh<Shape, T, Tag>::renumber( std::vector<size_type> const& node_map, mpl::int_<1> )
668 {
669 }
670 template<typename Shape, typename T, int Tag>
671 void
672 Mesh<Shape, T, Tag>::renumber( std::vector<size_type> const& node_map, mpl::int_<2> )
673 {
674 }
675 
676 template<typename Shape, typename T, int Tag>
677 void
678 Mesh<Shape, T, Tag>::renumber( std::vector<size_type> const& node_map, mpl::int_<3> )
679 {
680  for ( auto elt = this->beginEdge();
681  elt != this->endEdge(); ++elt )
682  {
683  edge_type __edge = *elt;
684 
685  DVLOG(2) << "edge id: " << __edge.id()
686  << " marker: " << __edge.marker() << "\n";
687 
688  // renumber the nodes of the face
689  for ( int i = 0; i < __edge.nPoints(); ++i )
690  {
691  size_type __true_id =__edge.point( i ).id();
692  this->edges().modify( elt,
693  lambda::bind( &edge_type::setPoint,
694  lambda::_1,
695  lambda::constant( i ),
696  boost::cref( this->point( node_map[__true_id] ) ) ) );
697  edge_type __edge2 = *elt;
698 
699  DVLOG(2) << "renumber edge: id1= " << __edge2.point( 0 ).id() << " id2= " << __edge2.point( 1 ).id()<< "\n";
700  DVLOG(2) << "renumber edge: point lid = " << i << " id = " << __true_id
701  << " nid = " << this->point( node_map[__true_id] ).id()
702  << " new point id = " << elt->point( i ).id() << "\n";
703 
704  }
705  }
706 }
707 template<typename Shape, typename T, int Tag>
708 void
710 {
711  for ( element_const_iterator elt = this->beginElement();
712  elt != this->endElement(); ++elt )
713  {
714  element_type __element = *elt;
715 
716  int16_type n_swap = 0;
717 
718  for ( uint16_type i=Shape::numVertices-1; i > 1; --i )
719  {
720  uint16_type min_index = 0;
721 
722  for ( uint16_type j=1; j < __element.nPoints()+i-( Shape::numVertices-1 ); ++j )
723  {
724  if ( __element.point( j ).id() < __element.point( min_index ).id() )
725  min_index = j;
726  }
727 
728  if ( i != min_index )
729  {
730  __element.swapPoints( i, min_index );
731  n_swap += int16_type( 1 );
732  }
733  }
734 
735  // Adds consistent permutation anti-clockwise for the two last points
736  if ( n_swap == 1 )
737  __element.swapPoints( 0,1 );
738 
739  // Modify the points ids
740  if ( n_swap > 0 )
741  {
742  for ( uint16_type i=0; i < __element.nPoints(); ++i )
743  {
744  auto& new_pt = this->point( __element.point( i ).id() );
745  this->elements().modify( elt, [i,&new_pt]( element_type& e ) { e.setPoint( i, new_pt ); } );
746  }
747  }
748  }
749 }
751 template<typename Shape, typename T, int Tag>
752 void
754 {
755  //updateEntitiesCoDimensionOne( mpl::bool_<nDim==nRealDim>() );
756  updateEntitiesCoDimensionOne( mpl::bool_<true>() );
757 }
758 template<typename Shape, typename T, int Tag>
759 void
761 {}
762 template<typename Shape, typename T, int Tag>
763 void
765 {
766  boost::timer ti;
767  face_type face;
768  //face.setWorldComm(this->worldComm());
769  face.setProcessIdInPartition( this->worldComm().localRank() );
770 
771  std::map<std::set<int>, size_type > _faces;
772  typename std::map<std::set<int>, size_type >::iterator _faceit;
773  int next_face = 0;
774  element_type ele;
775 
776  // First We check if we have already Faces stored
777  if ( ! this->faces().empty() )
778  {
779  // dump all faces in the container, to maintain the correct numbering
780  // if everything is correct the numbering in the bareface structure
781  // will reflect the actual face numbering However, if I want to create
782  // the internal faces I need to make sure that I am processing only the
783  // boundary ones. So I resize the container!
784  //if ( cf ) this->faces().resize( _numBFaces );
785 
786 
787  face_iterator __it = this->beginFace();
788  face_iterator __en = this->endFace();
789  LOG(INFO) << "We have " << std::distance( __it, __en ) << " faces in the database";
790  auto itb = this->beginFaceOnBoundary();
791  auto enb = this->beginFaceOnBoundary();
792  LOG(INFO) << "We have " << std::distance( itb, enb ) << " faces on the boundary in the database";
793  for ( ; __it!=__en; )
794  {
795  std::set<int> s;
796 
797  for ( int f = 0; f < face_type::numVertices; ++f )
798  {
799  s.insert( __it->point( f ).id() );
800  }
801 
802  bool faceinserted = false;
803  boost::tie( _faceit, faceinserted ) = _faces.insert( std::make_pair( s, next_face ) );
804 
805  if ( faceinserted )
806  ++next_face;
807 
808  DVLOG_IF(2,faceinserted) << "added face with id " << __it->id () << "\n";
809  DVLOG_IF(2,!faceinserted) << "not added face with id " << __it->id ()
810  << " was already face with id = " << _faceit->second << "\n";
811 
812 
813  if ( faceinserted == false )
814  {
815  // here we get the next face or \c end()
816  size_type theid = __it->id();
817 
818  // find the other face
819  face_iterator __other = this->faces().find( face_type( _faceit->second ) );
820  FEELPP_ASSERT( __other->id() != theid )
821  ( __other->id() )
822  ( theid ).error( "faces should have different ids " );
823 
824  if ( __it->marker() != __other->marker() )
825  {
826  this->faces().modify( __other, [__it]( face_type& f ) { f.setMarker2( __it->marker().value() ); } );
827  }
828 
829  __it = this->eraseFace( __it );
830 
831 
832 
833  }
834 
835  else
836  {
837  // ensure that item handler ids are in sync with
838  // faces ids
839  face_type __f = *__it;
840  __f.setId( _faceit->second );
841 
842  // set the id for this partition
843  //__f.setIdInPartition( this->worldComm().localRank(),__f.id() );
844 
845  DVLOG(2) << "set face id " << __f.id()
846  << " iterator id = " << __it->id()
847  << " check id = " << _faceit->second << "\n";
848 
849  this->faces().replace( __it, __f );
850  ++__it;
851  }
852  }
853  }
854 
855  VLOG(2) << "[Mesh::updateFaces] adding faces : " << ti.elapsed() << "\n";
856  ti.restart();
857 
858  VLOG(2) << "[Mesh::updateFaces] numLocalFaces : " << this->numLocalFaces() << "\n";
859  VLOG(2) << "[Mesh::updateFaces] face_type::numVertices : " << face_type::numVertices << "\n";
860  element_iterator iv, en;
861  boost::tie( iv, en ) = this->elementsRange();
862 
863  for ( ; iv != en; ++iv )
864  {
865  element_type const& __element = *iv;
866  size_type __element_id = __element.id();
867 
868  //MakeBareEntity<element_type,nDim> baremaker( __element );
869  for ( size_type j = 0; j < this->numLocalFaces(); j++ )
870  {
871  DVLOG(2) << "------------------------------------------------------------\n";
872  DVLOG(2) << "Element id: " << iv->id() << " local face id: " << j << "\n";
873 
874  std::set<int> s;
875 
876  for ( int f = 0; f < face_type::numVertices; ++f )
877  {
878  uint16_type pt_localid = ( nDim==1 )?j:iv->fToP( j, f );
879  s.insert( iv->point( pt_localid ).id() );
880  VLOG(3) << "add point local id " << f << " to face " << j << " " << iv->fToP( j, f )
881  << " global id " << iv->point( pt_localid ).id() << "\n";
882  }
883 
884  bool faceinserted = false;
885  boost::tie( _faceit, faceinserted ) = _faces.insert( std::make_pair( s, next_face ) );
886 
887  if ( faceinserted )
888  ++next_face;
889 
890  DVLOG(2) << "------------------------------------------------------------\n";
891  DVLOG(2) << "Element id: " << iv->id() << " local face id: " << j << "\n";
892 
893  if ( faceinserted )
894  {
895  DVLOG(2) << "creating the face:" << _faceit->second << "\n";
896 
897  // set face id
898  face.setId( _faceit->second );
899  face.disconnect();
900  if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
901  face.addElement( __element_id );
902 
903  // set the process id from element
904  face.setProcessId( __element.processId() );
905 
906  // set the id for this partition
907  //face.setIdInPartition( this->worldComm().localRank(),face.id() );
908 
909  // set the vertices of the face
910  for ( size_type k = 0; k < face_type::numPoints; ++k )
911  face.setPoint( k, iv->point( ele.fToP( j, k ) ) );
912 
913  // set the connection with the element
914  face.setConnection0( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) );
915  face.setOnBoundary( true, face_type::nDim );
916 
917  // adding the face
918  bool inserted = false;
919  face_iterator __fit;
920  boost::tie( __fit, inserted ) = this->addFace( face );
921  FEELPP_ASSERT( inserted && __fit != this->endFace() )
922  ( _faceit->second )
923  ( iv->id() )
924  ( __fit->id() )
925  ( face.id() ).error( "invalid face iterator" );
926  this->elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
927 
928 
929  DVLOG(2) << "Adding [new] face info : \n";
930  DVLOG(2) << "element id: " << __element_id << "\n";
931  DVLOG(2) << "process id: " << __fit->processId() << "\n";
932  DVLOG(2) << "id: " << __fit->id() << "\n";
933  DVLOG(2) << "bdy: " << __fit->isOnBoundary() << "\n";
934  DVLOG(2) << "marker: " << __fit->marker() << "\n";
935  DVLOG(2) << "ad_first: " << __fit->ad_first() << "\n";
936  DVLOG(2) << "pos_first: " << __fit->pos_first() << "\n";
937  DVLOG(2) << "proc_first: " << __fit->proc_first() << "\n";
938  DVLOG(2) << "ad_second: " << __fit->ad_second() << "\n";
939  DVLOG(2) << "pos_second: " << __fit->pos_second() << "\n";
940  DVLOG(2) << "proc_second: " << __fit->proc_second() << "\n";
941 
942  }
943 
944  else
945  {
946  DVLOG(2) << "found the face:" << _faceit->second << " in element " << __element_id << " and local face: " << j << "\n";
947 
948  // look in the face table for the face
949  face_iterator __fit = this->faces().find( face_type( _faceit->second ) );
950  FEELPP_ASSERT( __fit != this->endFace() )( _faceit->second ).error( "face is not in face container" );
951 
952 
953  face_type face = *__fit;
954  DVLOG(2) << "the face id :" << __fit->id() << "\n";
955 
956  if ( this->components().test( MESH_ADD_ELEMENTS_INFO ) )
957  face.addElement( __element_id );
958 
959  // the three conditions below typically arise after reading a serialized mesh
960  if ( __fit->isConnectedTo0() && __fit->connection0().template get<0>() == 0 && ( __element.id() == __fit->ad_first() ) )
961  {
962  DVLOG(2) << "fixing connection 0 in face\n";
963  // reconnect the elements to the face
964  auto connect0 = __fit->connection0();
965  connect0.template get<0>() = boost::addressof( __element );
966  face.setConnection0( connect0 );
967 
968  // need to reconstruct the neighbors
969  this->elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
970 
971  }
972  if ( __fit->isConnectedTo1() && __fit->connection1().template get<0>() == 0 && ( __element.id() == __fit->ad_second() ) )
973  {
974  DVLOG(2) << "fixing connection 1 in face\n";
975 
976  // reconnect the elements to the face
977  auto connect1 = __fit->connection1();
978  connect1.template get<0>() = boost::addressof( __element );
979  face.setConnection1( connect1 );
980 
981  // need to reconstruct the neighbors
982  element_iterator elt1 = this->elementIterator( __fit->ad_first(), __fit->proc_first() );
983  this->elements().modify( elt1, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
984  this->elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
985 
986  }
987  if ( __fit->isConnectedTo0() && __fit->isConnectedTo1() )
988  {
989  DVLOG(2) << "internal face, fixing process id if necessary\n";
990  if ( !iv->facePtr( j ) )
991  this->elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
992  FEELPP_ASSERT( iv->facePtr( j ) )( j )( iv->id() ).warn( "invalid element face error" );
993  FEELPP_ASSERT( face.isConnectedTo0() && face.isConnectedTo1() )
994  ( face.isConnectedTo0() )( face.isConnectedTo1() ).error ("inconsistent data structure" );
995  if ( face.processId()!=this->worldComm().localRank() )
996  {
997  if ( ( face.element0().processId()==this->worldComm().localRank() ) ||
998  ( face.element1().processId()==this->worldComm().localRank() ) )
999  face.setProcessId( this->worldComm().localRank() );
1000  }
1001  }
1002 
1003 
1004  // the face could have been entered apriori given by
1005  // the mesh generator, so just set the connection0
1006  // properly .
1007  if ( !__fit->isConnectedTo0() )
1008  {
1009  DVLOG(2) << "[updateFaces][boundary] element: " << __element_id
1010  << " face: " << j << " id: " << _faceit->second << "\n";
1011 
1012  // set the connection with the element
1013  face.setConnection0( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) );
1014  // set the process id from element
1015  face.setProcessId( __element.processId() );
1016  face.setOnBoundary( true, face_type::nDim );
1017  //this->faces().modify( __fit,
1018  //detail::UpdateFaceConnection0<typename face_type::element_connectivity_type>( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) ) );
1019 
1020  this->faces().replace( __fit, face );
1021 
1022  this->elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
1023 
1024  DVLOG(2) << "adding [!isConnectedTo0] face info : \n";
1025  DVLOG(2) << "id: " << __fit->id() << "\n";
1026  DVLOG(2) << "process id: " << __fit->processId() << "\n";
1027  DVLOG(2) << "bdy: " << __fit->isOnBoundary() << "\n";
1028  DVLOG(2) << "marker: " << __fit->marker() << "\n";
1029  DVLOG(2) << "ad_first: " << __fit->ad_first() << "\n";
1030  DVLOG(2) << "pos_first: " << __fit->pos_first() << "\n";
1031  DVLOG(2) << "proc_first: " << __fit->proc_first() << "\n";
1032  DVLOG(2) << "ad_second: " << __fit->ad_second() << "\n";
1033  DVLOG(2) << "pos_second: " << __fit->pos_second() << "\n";
1034  DVLOG(2) << "proc_second: " << __fit->proc_second() << "\n";
1035  DVLOG(2) << "element process id: " << iv->processId() << "\n";
1036  }
1037 
1038  // we found an internal face
1039  else if ( !__fit->isConnectedTo1() )
1040  {
1041 
1042 
1043 #if 0
1044  this->faces().modify( __fit,
1045  detail::UpdateFaceConnection1<typename face_type::element_connectivity_type>( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) ) );
1046 
1047 
1048  FEELPP_ASSERT( __fit->isConnectedTo0() && __fit->isConnectedTo1() )
1049  ( __fit->isConnectedTo0() )( __fit->isConnectedTo1() ).error( "invalid face connection" );
1050 #endif
1051  detail::UpdateFaceConnection1<typename face_type::element_connectivity_type> update1( boost::make_tuple( boost::addressof( __element ), __element_id, j, __element.processId() ) );
1052  update1( face );
1053  face.setOnBoundary( false );
1054 
1055  // force processId equal to M_comm.rank() if face on interprocessfaces
1056  if ( face.processId()!=this->worldComm().localRank() )
1057  {
1058  if ( ( face.element0().processId()==this->worldComm().localRank() ) || ( face.element1().processId()==this->worldComm().localRank() ) )
1059  face.setProcessId( this->worldComm().localRank() );
1060  }
1061 
1062  this->faces().replace( __fit, face );
1063 #if 0
1064  // update neighbors for each element and replace in element container
1065  element_iterator elt1 = this->elementIterator( __fit->ad_first(), __fit->proc_first() );
1066  this->elements().modify( elt1, update_element_neighbor_type( __fit->pos_first(),
1067  __fit->ad_second() ) );
1068  this->elements().modify( iv, update_element_neighbor_type( __fit->pos_second(),
1069  __fit->ad_first() ) );
1070 #endif
1071  // be careful: this step is crucial to set proper neighbor
1072  // connectivity
1073  element_iterator elt1 = this->elementIterator( __fit->ad_first(), __fit->proc_first() );
1074  this->elements().modify( elt1, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
1075  this->elements().modify( iv, detail::UpdateFace<face_type>( boost::cref( *__fit ) ) );
1076 
1077 
1078  DVLOG(2) << "adding face info : \n";
1079  DVLOG(2) << "id: " << __fit->id() << "\n";
1080  DVLOG(2) << "process id: " << __fit->processId() << "\n";
1081  DVLOG(2) << "bdy: " << __fit->isOnBoundary() << "\n";
1082  DVLOG(2) << "marker: " << __fit->marker() << "\n";
1083  DVLOG(2) << "ad_first: " << __fit->ad_first() << "\n";
1084  DVLOG(2) << "pos_first: " << __fit->pos_first() << "\n";
1085  DVLOG(2) << "proc_first: " << __fit->proc_first() << "\n";
1086  DVLOG(2) << "ad_second: " << __fit->ad_second() << "\n";
1087  DVLOG(2) << "pos_second: " << __fit->pos_second() << "\n";
1088  DVLOG(2) << "proc_second: " << __fit->proc_second() << "\n";
1089  DVLOG(2) << "element1 process id: " << elt1->processId() << "\n";
1090  DVLOG(2) << "element2 process id: " << iv->processId() << "\n";
1091  }
1092 
1093  CHECK( (__fit->processId() == __fit->proc_first()) ||
1094  (__fit->processId() == __fit->proc_second()) )
1095  << "invalid process id " << __fit->processId() << " with element proc first = " << __fit->proc_first()
1096  << " and element proc second " << __fit->proc_second();
1097  }
1098 
1099  LOG_IF( WARNING, !iv->facePtr( j ) ) << "invalid element " << iv->id() << " with local face id " << j;
1100  } // face loop
1101  } // element loop
1102 
1103 
1104  face_iterator f_it = this->beginFace();
1105  face_iterator f_en = this->endFace();
1106 
1107  for ( ; f_it!=f_en; )
1108  {
1109  // cleanup the face data structure :
1110 
1111  if ( !f_it->isConnectedTo0() )
1112  {
1113  DLOG(INFO) << "removing face id : " << f_it->id()
1114  << " marker : " << f_it->marker();
1115  // remove all faces that are not connected to any elements
1116  f_it = this->faces().erase( f_it );
1117  //++f_it;
1118  }
1119  else
1120  {
1121  ++f_it;
1122  }
1123 
1124  }
1125 
1126  LOG(INFO) << "We have now " << nelements(boundaryfaces(this)) << " faces on the boundary in the database";
1127  VLOG(2) << "element/face connectivity : " << ti.elapsed() << "\n";
1128  ti.restart();
1129 }
1130 
1131 template<typename Shape, typename T, int Tag>
1132 void
1133 Mesh<Shape, T, Tag>::modifyEdgesOnBoundary( face_iterator& it , mpl::bool_<true> )
1134 {
1135  // loop over face edges
1136  for ( int f = 0; f < face_type::numEdges; ++f )
1137  {
1138  if ( it->edge( f ).isOnBoundary() == false )
1139  {
1140  auto eit = this->edgeIterator( it->edge(f).id() );
1141  this->edges().modify( eit,
1142  []( edge_type& e )
1143  {e.setOnBoundary(true, 1 );} );
1144  }
1145 
1146  }
1147  DVLOG(3) << "We have " << nelements(boundaryedges(this)) << " boundary edges";
1148 
1149 }
1150 template<typename Shape, typename T, int Tag>
1151 void
1152 Mesh<Shape, T, Tag>::modifyEdgesOnBoundary( face_iterator& f, mpl::bool_<false> )
1153 {
1154 }
1155 
1156 template<typename Shape, typename T, int Tag>
1157 bool
1158 Mesh<Shape, T, Tag>::modifyElementOnBoundaryFromEdge( element_iterator& e, mpl::bool_<false> )
1159 {
1160  return false;
1161 }
1162 template<typename Shape, typename T, int Tag>
1163 bool
1164 Mesh<Shape, T, Tag>::modifyElementOnBoundaryFromEdge( element_iterator& iv, mpl::bool_<true> )
1165 {
1166  // in 3D check if the edges of the element touch the boundary
1167  bool isOnBoundary = false;
1168  for ( size_type j = 0; j < iv->nEdges(); j++ )
1169  {
1170  isOnBoundary |= iv->edge( j ).isOnBoundary();
1171  }
1172  if ( isOnBoundary )
1173  {
1174  // the element touches the boundary with just an edge
1175  this->elements().modify( iv, [isOnBoundary]( element_type& e ) { e.setOnBoundary( isOnBoundary, 1 ); } );
1176  return isOnBoundary;
1177  }
1178  return isOnBoundary;
1179 }
1180 
1181 template<typename Shape, typename T, int Tag>
1182 void
1183 Mesh<Shape, T, Tag>::updateOnBoundary()
1184 {
1185  // first go through all the faces and set the points of the boundary
1186  // faces to be on the boundary
1187  LOG(INFO) << "update boundary points...";
1188  for( auto it = this->beginFace(), en = this->endFace(); it != en; ++ it )
1189  {
1190  if ( it->isOnBoundary() == true )
1191  {
1192  modifyEdgesOnBoundary( it, mpl::bool_<(nDim >= 3)>() );
1193  // loop over face points
1194  for ( int f = 0; f < face_type::numPoints; ++f )
1195  {
1196  if ( it->point( f ).isOnBoundary() == false )
1197  {
1198  auto pit = this->pointIterator( it->point(f).id() );
1199  this->points().modify( pit,
1200  []( point_type& p )
1201  {p.setOnBoundary(true, 0 );} );
1202  }
1203 
1204  }
1205  }
1206  }
1207  LOG(INFO) << "We have " << nelements(boundarypoints(this)) << " boundary points";
1208  LOG(INFO) << "update boundary elements...";
1209  // loop through faces to set the elements having a face on the boundary
1210  for ( auto iv = this->beginElement(), en = this->endElement();
1211  iv != en; ++iv )
1212  {
1213  bool isOnBoundary = false;
1214 
1215  // first check if a face is on the boundary
1216  for ( size_type j = 0; j < iv->nTopologicalFaces(); j++ )
1217  {
1218  isOnBoundary |= iv->face( j ).isOnBoundary();
1219  }
1220 
1221  if ( isOnBoundary )
1222  {
1223  VLOG(3) << "checking " << iv->nTopologicalFaces() << " faces, isOnBoundary: " << isOnBoundary << " face_type::nDim: " << face_type::nDim;
1224  this->elements().modify( iv, [isOnBoundary]( element_type& e ) { e.setOnBoundary( isOnBoundary, face_type::nDim ); } );
1225  // go to the next element, no need to look further
1226  continue;
1227  }
1228  bool e_modified = modifyElementOnBoundaryFromEdge( iv, mpl::bool_<(nDim>=3)>() );
1229  // go to next element if element is on boundary
1230  if ( e_modified )
1231  {
1232  continue;
1233  }
1234 
1235  // finally check if a point of the element touches the boundary
1236  for ( size_type j = 0; j < iv->nPoints(); j++ )
1237  {
1238  isOnBoundary |= iv->point( j ).isOnBoundary();
1239  }
1240 
1241  if ( isOnBoundary )
1242  {
1243  VLOG(3) << "checking " << iv->nPoints() << " points, isOnBoundary: " << isOnBoundary;
1244  this->elements().modify( iv, [isOnBoundary]( element_type& e ) { e.setOnBoundary( isOnBoundary, 0 ); } );
1245  }
1246  } // loop over the elements
1247  LOG(INFO) << "[updateOnBoundary] We have " << nelements(boundaryelements(this))
1248  << " elements sharing a point, a edge or a face with the boundary in the database";
1249  BOOST_FOREACH( auto e, this->boundaryElements( 0, 2, 0 ) )
1250  {
1251  VLOG(3) << "boundary element : " << e.id()
1252  << " entity on boundary max dim " << e.boundaryEntityDimension()
1253  << " process id : " << e.processId();
1254  }
1255 
1256 }
1257 template<typename Shape, typename T, int Tag>
1258 void
1259 Mesh<Shape, T, Tag>::removeFacesFromBoundary( std::initializer_list<uint16_type> markers )
1260 {
1261  std::for_each( markers.begin(), markers.end(),
1262  [=]( uint16_type marker )
1263  {
1264  auto range=markedfaces( this, marker );
1265  LOG(INFO) << "removing " << nelements(range) << " faces marked " << marker << " from boundary faces\n";
1266 
1267  for( auto it = range.template get<1>(), en = range.template get<2>(); it != en; ++it )
1268  {
1269  if ( it->isOnBoundary() )
1270  {
1271  DVLOG(3) << "removing face " << it->id() << "\n";
1272  auto it2 = this->faces().template project<0>( it );
1273  this->faces().modify( it2, []( face_type & f ) { f.setOnBoundary( false ); } );
1274 
1275 
1276  //face_type f = *it2;
1277  //f.setOnBoundary( false );
1278  //this->faces().replace( it2, f );
1279  CHECK( it2->isOnBoundary() ==false ) << " face should not be on the boundary anymore\n";
1280  CHECK( it->isOnBoundary() ==false ) << " face should not be on the boundary anymore\n";
1281  }
1282 
1283  }
1284  } );
1285 }
1286 
1287 #if defined(FEELPP_ENABLE_MPI_MODE)
1288 template<typename Shape, typename T, int Tag>
1289 void
1291 {
1292  typedef std::vector< boost::tuple<size_type, std::vector<double> > > resultghost_type;
1293 
1294  VLOG(2) << "[Mesh::updateEntitiesCoDimensionOneGhostCell] start on god rank "<< this->worldComm().godRank() << "\n";
1295 
1296  std::vector<int> nbMsgToSend( this->worldComm().localSize() );
1297  std::fill( nbMsgToSend.begin(),nbMsgToSend.end(),0 );
1298 
1299  std::vector< std::map<int,int> > mapMsg( this->worldComm().localSize() );
1300 
1301 
1302  auto iv = this->beginGhostElement();
1303  auto en = this->endGhostElement();
1304 
1305  for ( ; iv != en; ++iv )
1306  {
1307  element_type const& __element = *iv;
1308  const int IdProcessOfGhost = __element.processId();
1309  const size_type idInPartition = __element.idInOthersPartitions( IdProcessOfGhost );
1310 
1311  this->elements().modify( this->elementIterator(iv->id(),IdProcessOfGhost) , typename super_elements::ElementGhostConnectPointToElement() );
1312  if ( nDim==3 )
1313  this->elements().modify( this->elementIterator(iv->id(),IdProcessOfGhost) , typename super_elements::ElementGhostConnectEdgeToElement() );
1314 
1315  // send
1316  this->worldComm().localComm().send( IdProcessOfGhost , nbMsgToSend[IdProcessOfGhost], idInPartition );
1317 #if 0
1318  std::cout<< "I am the proc" << this->worldComm().localRank()<<" , I send to proc " << IdProcessOfGhost
1319  <<" with tag "<< nbMsgToSend[IdProcessOfGhost]
1320  << " idSend " << idInPartition
1321  << " it_ghost->G() " << __element.G()
1322  << std::endl;
1323 #endif
1324  // save tag of request
1325  mapMsg[IdProcessOfGhost].insert( std::make_pair( nbMsgToSend[IdProcessOfGhost],__element.id() ) );
1326  // update nb send
1327  ++nbMsgToSend[IdProcessOfGhost];
1328  }
1329 
1330  //------------------------------------------------------------------------------------------------//
1331  // counter of msg received for each process
1332  std::vector<int> nbMsgToRecv;
1333  mpi::all_to_all( this->worldComm().localComm(),
1334  nbMsgToSend,
1335  nbMsgToRecv );
1336 
1337  //------------------------------------------------------------------------------------------------//
1338  // recv id asked and re-send set of face id
1339  for ( int proc=0; proc<this->worldComm().localSize(); ++proc )
1340  {
1341  for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
1342  {
1343  //recv
1344  size_type idRecv;
1345  this->worldComm().localComm().recv( proc, cpt, idRecv );
1346 #if 0
1347  std::cout<< "I am the proc" << this->worldComm().localRank()<<" I receive to proc " << proc
1348  <<" with tag "<< cpt << " idRecv " << idRecv
1349  << " it_ghost->G() " << this->element( idRecv ).G()
1350  << std::endl;
1351 #endif
1352  auto const& theelt = this->element( idRecv );
1353 
1354 
1355  // get faces id and bary
1356  resultghost_type idFacesWithBary(this->numLocalFaces(),boost::make_tuple(0, std::vector<double>(nRealDim) ));
1357  for ( size_type j = 0; j < this->numLocalFaces(); j++ )
1358  {
1359  auto const& theface = theelt.face( j );
1360  idFacesWithBary[j].template get<0>() = theface.id();
1361  //auto const& theGj = theface.G();
1362  auto const& theGj = theface.vertices();
1363 #if 1
1364  //compute face barycenter
1365  typename Localization::matrix_node_type v( theGj.size1(), 1 );
1366  ublas::scalar_vector<T> avg( theGj.size2(), T( 1 ) );
1367  T n_val = int( theGj.size2() );
1368 
1369  for ( size_type i = 0; i < theGj.size1(); ++i )
1370  v( i, 0 ) = ublas::inner_prod( ublas::row( theGj, i ), avg )/n_val;
1371 
1372  auto baryFace = ublas::column( v,0 );
1373 #else // doesn't work
1374  /*auto GLASbaryFace =*/ //Feel::glas::average(jkj);
1375  auto baryFace = ublas::column( glas::average( theface.G() ),0 );
1376 #endif
1377 
1378  // save facebarycenter by components
1379  for ( uint16_type comp = 0; comp<nRealDim; ++comp )
1380  {
1381  idFacesWithBary[j].template get<1>()[comp]=baryFace[comp];
1382  }
1383  }
1384 
1385  // get points id and nodes
1386  resultghost_type idPointsWithNode(element_type::numLocalVertices,boost::make_tuple(0, std::vector<double>(nRealDim) ));
1387  for ( size_type j = 0; j < element_type::numLocalVertices; j++ )
1388  {
1389  auto const& thepoint = theelt.point( j );
1390  idPointsWithNode[j].template get<0>() = thepoint.id();
1391  for ( uint16_type comp = 0; comp<nRealDim; ++comp )
1392  {
1393  idPointsWithNode[j].template get<1>()[comp]=thepoint(comp);
1394  }
1395  }
1396 
1397  std::vector<resultghost_type> theresponse(2);
1398  theresponse[0] = idPointsWithNode;
1399  theresponse[1] = idFacesWithBary;
1400  //auto theresponse = boost::make_tuple( idPointsWithNode, idFacesWithBary );
1401  // send response
1402  //this->worldComm().localComm().send( proc, cpt, idFacesWithBary );
1403  this->worldComm().localComm().send( proc, cpt, theresponse );
1404  }
1405  }
1406 
1407  //------------------------------------------------------------------------------------------------//
1408  // get response to initial request and update Feel::Mesh::Faces data
1409  for ( int proc=0; proc<this->worldComm().localSize(); ++proc )
1410  {
1411  for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
1412  {
1413  //recv
1414 #if 0
1415  std::vector< boost::tuple<size_type, std::vector<double> > > idFacesWithBaryRecv(this->numLocalFaces());
1416  this->worldComm().localComm().recv( proc, cpt, idFacesWithBaryRecv );
1417 #elif 0
1418  typedef std::vector< boost::tuple<size_type, std::vector<double> > > resultghost_face_type;
1419  typedef std::vector< boost::tuple<size_type, std::vector<double> > > resultghost_point_type;
1420  boost::tuple<resultghost_point_type,resultghost_face_type> requestRecv;
1421  this->worldComm().localComm().recv( proc, cpt, requestRecv );
1422  auto const& idPointsWithNodeRecv = requestRecv.template get<0>();
1423  auto const& idFacesWithBaryRecv = requestRecv.template get<1>();
1424 #else
1425  std::vector<resultghost_type> requestRecv;
1426  this->worldComm().localComm().recv( proc, cpt, requestRecv );
1427  auto const& idPointsWithNodeRecv = requestRecv[0];
1428  auto const& idFacesWithBaryRecv = requestRecv[1];
1429 #endif
1430 
1431 #if 0
1432  std::cout<< "I am the proc " << this->worldComm().localRank()<<" I receive to proc " << proc
1433  <<" with tag "<< cpt << std::endl;
1434 #endif
1435  auto const& theelt = this->element( mapMsg[proc][cpt],proc );
1436 
1437  //update faces data
1438  for ( size_type j = 0; j < this->numLocalFaces(); j++ )
1439  {
1440  auto const& idFaceRecv = idFacesWithBaryRecv[j].template get<0>();
1441  auto const& baryFaceRecv = idFacesWithBaryRecv[j].template get<1>();
1442 
1443  //objective : find face_it (hence jBis in theelt ) (permutations would be necessary)
1444  uint16_type jBis = invalid_uint16_type_value;
1445  bool hasFind=false;
1446  for ( uint16_type j2 = 0; j2 < this->numLocalFaces() && !hasFind; j2++ )
1447  {
1448 
1449  auto const& thefacej2 = theelt.face( j2 );
1450  //auto const& theGj2 = thefacej2.G();
1451  auto const& theGj2 = thefacej2.vertices();
1452 #if 1
1453  //compute face barycenter
1454  typename Localization::matrix_node_type v( theGj2.size1(), 1 );
1455  ublas::scalar_vector<T> avg( theGj2.size2(), T( 1 ) );
1456  T n_val = int( theGj2.size2() );
1457 
1458  for ( size_type i = 0; i < theGj2.size1(); ++i )
1459  v( i, 0 ) = ublas::inner_prod( ublas::row( theGj2, i ), avg )/n_val;
1460 
1461  auto baryFace = ublas::column( v,0 );
1462 #else // doesn't compile (I don't now why)
1463  auto const& thefacej2 =theelt.face( j2 );
1464  auto baryFace = ublas::column( glas::average( thefacej2.G() ),0 );
1465 #endif
1466  // compare barycenters
1467  bool find2=true;
1468  for (uint16_type d=0;d<nRealDim;++d)
1469  {
1470  find2 = find2 && ( std::abs( baryFace[d]-baryFaceRecv[d] )<1e-9 );
1471  }
1472  if (find2) { hasFind=true;jBis=j2; }
1473 
1474  } //for ( uint16_type j2 = 0; j2 < this->numLocalFaces() && !hasFind; j2++ )
1475 
1476  CHECK ( hasFind ) << "[mesh::updateEntitiesCoDimensionOneGhostCell] : invalid partitioning data, ghost face cells are not available\n";
1477 
1478  // get the good face
1479  auto face_it = this->faceIterator( theelt.face( jBis ).id() );
1480  //update the face
1481  this->faces().modify( face_it, detail::updateIdInOthersPartitions( proc, idFaceRecv ) );
1482 
1483  } // for ( size_type j = 0; j < this->numLocalFaces(); j++ )
1484 
1485 
1486  for ( size_type j = 0; j < element_type::numLocalVertices; j++ )
1487  {
1488  auto const& idPointRecv = idPointsWithNodeRecv[j].template get<0>();
1489  auto const& nodePointRecv = idPointsWithNodeRecv[j].template get<1>();
1490 
1491  uint16_type jBis = invalid_uint16_type_value;
1492  bool hasFind=false;
1493  for ( uint16_type j2 = 0; j2 < element_type::numLocalVertices && !hasFind; j2++ )
1494  {
1495  auto const& thepointj2 = theelt.point( j2 );
1496  // compare barycenters
1497  bool find2=true;
1498  for (uint16_type d=0;d<nRealDim;++d)
1499  {
1500  find2 = find2 && ( std::abs( thepointj2(d)-nodePointRecv[d] )<1e-9 );
1501  }
1502  if (find2) { hasFind=true;jBis=j2; }
1503  }
1504 
1505  CHECK ( hasFind ) << "[mesh::updateEntitiesCoDimensionOneGhostCell] : invalid partitioning data, ghost point cells are not available\n";
1506  // get the good face
1507  auto point_it = this->pointIterator( theelt.point( jBis ).id() );
1508  //update the face
1509  this->points().modify( point_it, detail::updateIdInOthersPartitions( proc, idPointRecv ) );
1510  } // for ( size_type j = 0; j < element_type::numLocalVertices; j++ )
1511 
1512  /*for ( size_type j = 0; j < element_type::numLocalEdges; j++ )
1513  {
1514  }*/
1515 
1516  } // for ( int cpt=0;cpt<nbMsgToSend[proc];++cpt)
1517  } // for (int proc=0; proc<M_comm.size();++proc)
1518 
1519  //------------------------------------------------------------------------------------------------//
1520 
1521  //std::cout << "[Mesh::updateEntitiesCoDimensionOneGhostCell] finish" << std::endl;
1522 
1523 } // updateEntitiesCoDimensionOneGhostCell
1524 #endif // defined(FEELPP_ENABLE_MPI_MODE)
1525 
1526 
1527 template<typename Shape, typename T, int Tag>
1528 void
1530 {
1531  if ( nDim != nRealDim )
1532  return;
1533 #if !defined( NDEBUG )
1534  VLOG(2) << "[Mesh::check] numLocalFaces = " << this->numLocalFaces() << "\n";
1535  element_iterator iv = this->beginElementWithProcessId( this->worldComm().localRank() );
1536  element_iterator en = this->endElementWithProcessId( this->worldComm().localRank() );
1537 
1538  //boost::tie( iv, en ) = this->elementsRange();
1539  for ( ; iv != en; ++iv )
1540  {
1541  element_type const& __element = *iv;
1542 
1543  for ( size_type j = 0; j < this->numLocalFaces(); j++ )
1544  {
1545  FEELPP_ASSERT( iv->facePtr( j ) )( j )( iv->id() ).error( "invalid element face check" );
1546  VLOG(2) << "------------------------------------------------------------\n";
1547  VLOG(2) << "Element : " << iv->id() << " face lid: " << j << " face gid: " << iv->face( j ).id() << "\n";
1548 
1549  }
1550 
1551  size_type counter = 0;
1552 
1553  for ( uint16_type ms=0; ms < __element.nNeighbors(); ms++ )
1554  {
1555  if ( __element.neighbor( ms ).first != invalid_size_type_value )
1556  ++counter;
1557 
1558  }
1559 
1560  VLOG(2) << "[Mesh::check] element " << __element.id() << " number of neighbors: " << counter << "\n";
1561  FEELPP_ASSERT( counter >= 1 )( __element.id() )( __element.nNeighbors() )( counter ).warn( "invalid neighboring data" );
1562 #if 0
1563 
1564  for ( size_type j = 0; j < ( size_type )element_type::numEdges; ++j )
1565  {
1566  FEELPP_ASSERT( iv->edgePtr( j ) )( j )( iv->id() ).error( "invalid element edge check" );
1567  VLOG(2) << "------------------------------------------------------------\n";
1568  VLOG(2) << "Element : " << iv->id() << " edge lid: " << j << " edge gid: " << iv->edge( j ).id() << "\n";
1569 
1570  }
1571 
1572 #endif
1573 
1574  }
1575 
1576  // faces check
1577  typedef typename super::location_face_const_iterator location_face_const_iterator;
1578  location_face_const_iterator itf = this->beginFaceOnBoundary();
1579  location_face_const_iterator ite = this->endFaceOnBoundary();
1580  std::map<int,int> nf;
1581  for ( ; itf != ite; ++ itf )
1582  {
1583  face_type const& __face = *itf;
1584  DLOG_IF( WARNING, !__face.isConnectedTo0() && !__face.isConnectedTo1() ) << "face not connected to an element face:" << __face << "\n";
1585  if ( !__face.isConnectedTo0() && !__face.isConnectedTo1() )
1586  {
1587  auto it = nf.find( (int)__face.marker().value() );
1588  if ( it == nf.end() )
1589  nf[(int)__face.marker().value()] = 1;
1590  else
1591  it->second += 1;
1592  }
1593  DLOG_IF( WARNING, !__face.isConnectedTo0() && __face.isConnectedTo1() ) << "face connected to element 1 but not element 0. face:" << __face << "\n";
1594  DLOG_IF( WARNING, __face.isConnectedTo0() && __face.isConnectedTo1() )<< "invalid boundary face (connected to 2 elements):" << __face << "\n";
1595 
1596  DLOG_IF( WARNING, __face.isConnectedTo0() && !__face.element( 0 ).facePtr( __face.pos_first() ) ) << "invalid face in element, face: " << __face << " in element " << __face.element(0) << "\n";
1597 
1598  }
1599  auto itt = nf.begin();
1600  auto ent = nf.end();
1601  for(; itt != ent; ++itt )
1602  {
1603  LOG(INFO ) << "face with marker " << itt->first << " not attached " << itt->second << "\n";
1604  }
1605 #endif
1606 }
1607 template<typename Shape, typename T, int Tag>
1608 void
1610 {
1611  // Don't need to do anything if there is
1612  // only one processor.
1613  if ( this->worldComm().localSize() == 1 )
1614  return;
1615 
1616 #ifdef FEELPP_HAS_MPI
1617 
1618  M_neighboring_processors.clear();
1619 
1620  // Get the bounding sphere for the local processor
1621  Sphere bounding_sphere = processorBoundingSphere ( *this, this->worldComm().localRank() );
1622 
1623  // Just to be sure, increase its radius by 10%. Sure would suck to
1624  // miss a neighboring processor!
1625  bounding_sphere.setRadius( bounding_sphere.radius()*1.1 );
1626 
1627  // Collect the bounding spheres from all processors, test for intersection
1628  {
1629  std::vector<double>
1630  send ( 4, 0 ),
1631  recv ( 4*this->worldComm().localSize(), 0 );
1632 
1633  send[0] = bounding_sphere.center()( 0 );
1634  send[1] = bounding_sphere.center()( 1 );
1635  send[2] = bounding_sphere.center()( 2 );
1636  send[3] = bounding_sphere.radius();
1637 
1638  MPI_Allgather ( &send[0], send.size(), MPI_DOUBLE,
1639  &recv[0], send.size(), MPI_DOUBLE,
1640  this->worldComm().localComm() );
1641 
1642 
1643  for ( unsigned int proc=0; proc<this->worldComm().localSize(); proc++ )
1644  {
1645  const Point center ( recv[4*proc+0],
1646  recv[4*proc+1],
1647  recv[4*proc+2] );
1648 
1649  const Real radius = recv[4*proc+3];
1650 
1651  const Sphere proc_sphere ( center, radius );
1652 
1653  if ( bounding_sphere.intersects( proc_sphere ) )
1654  M_neighboring_processors.push_back( proc );
1655  }
1656 
1657  // Print out the _neighboring_processors list
1658  VLOG(2) << "Processor " << this->worldComm().localRank() << " intersects:\n";
1659 
1660  for ( unsigned int p=0; p< M_neighboring_processors.size(); p++ )
1661  VLOG(2) << " - proc " << M_neighboring_processors[p] << "\n";
1662  }
1663 
1664 #endif
1665 }
1666 
1667 template<typename Shape, typename T, int Tag>
1668 void
1669 Mesh<Shape, T, Tag>::checkLocalPermutation( mpl::bool_<true> ) const
1670 {
1671  bool mesh_well_oriented = true;
1672  std::vector<size_type> list_of_bad_elts;
1673 
1674  for ( element_const_iterator elt = this->beginElement();
1675  elt != this->endElement(); ++elt )
1676  {
1677  element_type const& __element = *elt;
1678 
1679  if ( !__element.isAnticlockwiseOriented() )
1680  {
1681  mesh_well_oriented = false;
1682  list_of_bad_elts.push_back( elt->id() );
1683  }
1684  }
1685 
1686  if ( mesh_well_oriented )
1687  VLOG(2) << "Local numbering in the elements is OK . \n";
1688 
1689  else
1690  {
1691  std::for_each( list_of_bad_elts.begin(),
1692  list_of_bad_elts.end(),
1693  []( size_type const& e )
1694  {
1695  LOG_FIRST_N(WARNING,10) << "element is not anticlockwise oriented(wrong local numbering): " << e << "\n";
1696  });
1697 
1698 
1699  }
1700 }
1701 
1702 template<typename Shape, typename T, int Tag>
1703 void
1705 {
1706  element_type __element = *this->beginElement();
1707 
1708  if ( ! ( __element.isATriangleShape() ||
1709  __element.isATetrahedraShape() ) )
1710  return;
1711 
1712  static const uint16_type otn_triangle[ 10 ] =
1713  {
1714  1, 0, 2, 3, 5, 4
1715  };
1716  static const uint16_type otn_tetra[ 10 ] =
1717  {
1718  1, 0, 2, 3, 4, 6, 5, 8, 7, 9
1719  };
1720 
1721  for ( element_iterator elt = this->beginElement();
1722  elt != this->endElement(); ++elt )
1723  {
1724  element_type __element = *elt;
1725  bool is_anticlockwise = __element.isAnticlockwiseOriented();
1726  // --verbose
1727 #if 0
1728  FEELPP_ASSERT( is_anticlockwise == true )
1729  ( is_anticlockwise )
1730  ( __element.id() )
1731  ( __element.G() ).warn( "invalid element permutation, will fix it" );
1732 #endif // 0
1733 
1734  // fix permutation
1735  if ( !is_anticlockwise )
1736  {
1737  if ( __element.isATriangleShape() )
1738  __element.exchangePoints( otn_triangle );
1739 
1740  else if ( __element.isATetrahedraShape() )
1741  __element.exchangePoints( otn_tetra );
1742 
1743  else
1744  {
1745  FEELPP_ASSERT( 0 )
1746  ( __element.id() )
1747  ( __element.G() ).error( "invalid element type" );
1748  throw std::logic_error( "invalid element type" );
1749  }
1750 
1751  is_anticlockwise = __element.isAnticlockwiseOriented();
1752  FEELPP_ASSERT( is_anticlockwise == true )
1753  ( is_anticlockwise )
1754  ( __element.id() )
1755  ( __element.G() ).error( "invalid element permutation" );
1756 
1757  this->elements().replace( elt, __element );
1758  }
1759 
1760  }
1761 }
1762 
1763 template<typename Shape, typename T, int Tag>
1764 void
1766 {
1767  encode();
1768  VLOG(1) << "sending markername\n";
1769  //this->comm().send( p, tag, M_markername.size() );
1770  VLOG(1) << "sending markername size: "<< M_markername.size() << "\n";
1771  BOOST_FOREACH(auto m, M_markername )
1772  {
1773  VLOG(1) << "sending key: "<< m.first << "\n";
1774  //this->comm().send( p, tag, m.first );
1775  VLOG(1) << "sending value\n";
1776  //this->comm().send( p, tag, m.second );
1777  }
1778 }
1779 
1780 template<typename Shape, typename T, int Tag>
1781 void
1783 {
1784  VLOG(1) << "receiving markername\n";
1785  //this->comm().recv( p, tag, M_markername );
1786  int s = 0;
1787  //this->comm().recv( p, tag, s );
1788  VLOG(1) << "receiving markername size: "<< s << "\n";
1789  for( int i = 0; i < s; ++i )
1790  {
1791  std::string k;
1792  VLOG(1) << "receiving key\n";
1793  //this->comm().recv( p, tag, k );
1794  VLOG(1) << "receiving key:"<< k << "\n";
1795  std::vector<int> v;
1796  VLOG(1) << "receiving value\n";
1797  //this->comm().recv( p, tag, v );
1798  VLOG(1) << "receiving value: "<< v[0] << ","<< v[1] <<"\n";
1799  //M_markername[k]=v;
1800  }
1801 
1802 
1803  //decode();
1804 
1805 }
1806 
1807 template<typename Shape, typename T, int Tag>
1808 typename Mesh<Shape, T, Tag>::element_iterator
1809 Mesh<Shape, T, Tag>::eraseElement( element_iterator position, bool modify )
1810 {
1811  for(int i = 0; i < element_type::numTopologicalFaces; ++i )
1812  {
1813  auto fit = this->faces().iterator_to( position->face(i) );
1814  auto modified = this->faces().modify( fit, [position]( face_type& f ) { f.disconnect(*position); } );
1815  if ( fit->isOnBoundary() && !fit->isConnectedTo0() )
1816  {
1817  std::cout << "erase boundary face...\n";
1818  this->eraseFace( fit );
1819  }
1820  if ( !fit->isOnBoundary() && fit->isConnectedTo0() && !fit->isConnectedTo1() )
1821  {
1822  std::cout << "found boundary face...\n";
1823  this->faces().modify( fit, []( face_type& f ){ f.setOnBoundary( true ); } );
1824  }
1825  }
1826  this->elements().modify( position,
1827  [] ( element_type& e )
1828  {
1829  for ( int i = 0; i < e.numPoints; ++i )
1830  e.point( i ).elements().erase( e.id() );
1831  } );
1832  auto eit = this->elements().erase( position );
1833  //this->setUpdatedForUse( false );
1834  //this->updateForUse();
1835  return eit;
1836 }
1837 template<typename Shape, typename T, int Tag>
1838 void
1840 {
1841  //std::cout<<"encode= " << this->worldComm().localSize() << std::endl;
1842 
1843  M_enc_pts.clear();
1844  for( auto pt_it = this->beginPoint(), pt_en = this->endPoint(); pt_it != pt_en; ++pt_it )
1845  {
1846  std::vector<double> pts(3,0);
1847  pts[0] = pt_it->node()[0];
1848  if ( mesh_type::nRealDim >= 2 )
1849  pts[1] = pt_it->node()[1];
1850  if ( mesh_type::nRealDim >= 3 )
1851  pts[2] = pt_it->node()[2];
1852  M_enc_pts[pt_it->id()+1] = boost::make_tuple( pt_it->isOnBoundary(), pt_it->tags(), pts );
1853 
1854  }
1855  M_enc_elts.clear();
1856 
1857 
1858  auto allmarkedfaces = boundaryfaces( *this );
1859  auto face_it = allmarkedfaces.template get<1>();
1860  auto face_end = allmarkedfaces.template get<2>();
1861 
1862  int elem_number=1;
1863 
1864  GmshOrdering<element_type> ordering;
1865 
1866  GmshOrdering<face_type> ordering_face;
1867  // save the faces
1868 
1869  for ( ; face_it != face_end; ++face_it )
1870  {
1871  std::vector<int> faces;
1872  faces.push_back( ordering_face.type() );
1873  faces.push_back( 4 + face_it->numberOfPartitions() );
1874  faces.push_back( face_it->marker().value() );
1875  faces.push_back( face_it->marker2().value() );
1876  faces.push_back( face_it->numberOfPartitions() );
1877  faces.push_back( face_it->processId() );
1878  for ( size_type i=0 ; i<face_it->numberOfNeighborPartitions(); ++i )
1879  faces.push_back( -( face_it->neighborPartitionIds()[i] ) );
1880  for ( uint16_type p=0; p<face_type::numPoints; ++p )
1881  faces.push_back( face_it->point( ordering_face.fromGmshId( p ) ).id()+1 );
1882 
1883  M_enc_faces[face_it->id()] = faces;
1884  } // faces
1885 
1886 
1887  auto eltOnProccess = elements( *this );
1888  auto elt_it = eltOnProccess.template get<1>();
1889  auto elt_en = eltOnProccess.template get<2>();
1890 
1891  for ( ; elt_it != elt_en; ++elt_it )
1892  {
1893  std::vector<int> elts;
1894  elts.push_back( ordering.type() );
1895  elts.push_back( 4 + elt_it->numberOfPartitions() );
1896  elts.push_back( elt_it->marker().value() );
1897  elts.push_back( elt_it->marker2().value() );
1898  elts.push_back( elt_it->numberOfPartitions() );
1899  elts.push_back( elt_it->processId() );
1900  for ( size_type i=0 ; i<elt_it->numberOfNeighborPartitions(); ++i )
1901  elts.push_back( -( elt_it->neighborPartitionIds()[i] ) );
1902  for ( uint16_type p=0; p<element_type::numPoints; ++p )
1903  elts.push_back( elt_it->point( ordering.fromGmshId( p ) ).id()+1 );
1904 
1905  M_enc_elts[elt_it->id()] = elts;
1906  } // elements
1907  //std::cout<<"encode= " << this->worldComm().localSize() << std::endl;
1908 }
1909 
1910 template<typename Shape, typename T, int Tag>
1911 void
1913 {
1914 #if 0
1915  std::vector<int> mapWorld(this->worldComm().size());
1916  for(int cpu=0; cpu < this->worldComm().size(); ++cpu)
1917  mapWorld[cpu] = cpu;
1918  WorldComm worldcomm(mapWorld);
1919 
1920 
1921 
1922 
1923  // std::cout<<"decode= " << this->worldComm().size() << std::endl;
1924  //std::cout<<"decode= " << worldcomm.subWorldComm().localRank() << std::endl;
1925  this->setWorldComm(worldcomm.subWorldComm());
1926 #else
1927  LOG(INFO) <<"decode= " << this->worldComm().size() << "\n" ;
1928  LOG(INFO) <<"decode= " << this->worldComm().subWorldComm().localRank() << "\n";
1929 #endif
1930  static const uint16_type npoints_per_face = ( face_type::numVertices*face_type::nbPtsPerVertex+
1931  face_type::numEdges*face_type::nbPtsPerEdge+
1932  face_type::numFaces*face_type::nbPtsPerFace );
1933 
1934  static const uint16_type npoints_per_element = element_type::numPoints;
1935  for( auto pt_it = M_enc_pts.begin(), pt_en = M_enc_pts.end();
1936  pt_it != pt_en; ++pt_it )
1937  {
1938  node_type __n( nRealDim );
1939 
1940  for ( uint16_type j = 0; j < nRealDim; ++j )
1941  __n[j] = pt_it->second.template get<2>()[j];
1942 
1943  point_type __pt( pt_it->first-1,__n, pt_it->second.template get<0>() );
1944  __pt.setOnBoundary( pt_it->second.template get<0>() );
1945  __pt.setTags( pt_it->second.template get<1>() );
1946 
1947  this->addPoint( __pt );
1948  }
1949  for( auto face_it = M_enc_faces.begin(), face_en = M_enc_faces.end();
1950  face_it != face_en; ++face_it )
1951  {
1952  face_type pf;
1953  GmshOrdering<face_type> ordering;
1954  FEELPP_ASSERT( ordering.type() == face_it->second[0] )
1955  ( ordering.type() ) ( face_it->second[0] ).error( "invalid mesh face type" );
1956  std::vector<int> tags( face_it->second[1] );
1957  for(int i = 0; i < tags.size(); ++i ) tags[i] = face_it->second[2+i];
1958  pf.setTags( tags );
1959  pf.setId( this->numFaces() );
1960  pf.setProcessIdInPartition( this->worldComm().localRank() );
1961  pf.setProcessId( this->worldComm().localRank() );
1962  //pf.setIdInPartition( this->worldComm().localRank(),pf.id() );
1963 
1964  const int shift = face_it->second[1]+1;
1965  for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
1966  {
1967  pf.setPoint( ordering.fromGmshId( jj ), this->point( face_it->second[shift+jj]-1 ) );
1968  }
1969  this->addFace( pf );
1970  }
1971 
1972  for( auto elt_it = M_enc_elts.begin(), elt_en = M_enc_elts.end();
1973  elt_it != elt_en; ++elt_it )
1974  {
1975  element_type pv;
1976  GmshOrdering<element_type> ordering;
1977  FEELPP_ASSERT( ordering.type() == elt_it->second[0] )
1978  ( ordering.type() ) ( elt_it->second[0] ).error( "invalid mesh element type" );
1979  std::vector<int> tags( elt_it->second[1] );
1980  for(int i = 0; i < tags.size(); ++i ) tags[i] = elt_it->second[2+i];
1981  pv.setTags( tags );
1982  pv.setProcessIdInPartition( this->worldComm().localRank() );
1983  pv.setProcessId( this->worldComm().localRank() );
1984  //pv.setIdInPartition( this->worldComm().localRank(),pv.id() );
1985 
1986  const int shift = elt_it->second[1]+1;
1987  for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
1988  {
1989  pv.setPoint( ordering.fromGmshId( jj ), this->point( elt_it->second[shift+jj]-1 ) );
1990  }
1991  this->addElement( pv );
1992 #if 0
1993  __idGmshToFeel=pv.id();
1994  auto theelt = mesh->elementIterator( pv.id(), pv.partitionId() );
1995  mesh->elements().modify( theelt, detail::updateIdInOthersPartitions( this->worldComm().localRank(), pv.id() ) );
1996 #endif
1997  }
1998  LOG(INFO) << "distance elts: "<< std::distance( this->beginElement(), this->endElement() ) << "\n";
1999  LOG(INFO) << "distance faces: "<< std::distance( this->beginFace(), this->endFace() ) << "\n";
2000  LOG(INFO) << "distance marker faces: "<< std::distance( this->beginFaceWithMarker(), this->endFaceWithMarker() ) << "\n";
2001  LOG(INFO) << "distance marker2 faces: "<< std::distance( this->beginFaceWithMarker2(), this->endFaceWithMarker2() ) << "\n";
2002 
2003  //this->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
2004  //this->updateForUse();
2005  //std::cout<<"decode= " << this->worldComm().localSize() << std::endl;
2006 }
2007 
2008 template<typename Shape, typename T, int Tag1, int Tag2=Tag1, int TheTag=Tag1>
2009 boost::shared_ptr<Mesh<Shape, T, TheTag> >
2010 merge( boost::shared_ptr<Mesh<Shape, T, Tag1> > m1,
2011  boost::shared_ptr<Mesh<Shape, T, Tag2> > m2 )
2012 {
2013  boost::shared_ptr<Mesh<Shape, T, TheTag> > m( new Mesh<Shape, T, TheTag> );
2014 
2015  // add the points
2016 
2017  size_type shift_p = std::distance( m1->beginPoint(), m1->endPoint() );
2018  for( auto it = m1->beginPoint(), en = m1->endPoint(); it != en; ++it )
2019  {
2020  m->addPoint( *it );
2021  }
2022  for( auto it = m2->beginPoint(), en = m2->endPoint(); it != en; ++it )
2023  {
2024  auto p = *it;
2025  // shift id
2026  p.setId( shift_p+p.id() );
2027  m->addPoint( p );
2028  }
2029  typedef typename Mesh<Shape, T, TheTag>::face_iterator face_iterator;
2030  typedef typename Mesh<Shape, T, TheTag>::element_iterator element_iterator;
2031  typedef typename Mesh<Shape, T, TheTag>::face_type face_type;
2032  typedef typename Mesh<Shape, T, TheTag>::element_type element_type;
2033 
2034  auto addface = [&]( boost::shared_ptr<Mesh<Shape, T, TheTag> > m, face_iterator it, size_type shift )
2035  {
2036  face_type pf = *it;
2037  pf.disconnect();
2038  pf.setOnBoundary( true );
2039  const uint16_type npoints_per_face = ( face_type::numVertices*face_type::nbPtsPerVertex+
2040  face_type::numEdges*face_type::nbPtsPerEdge+
2041  face_type::numFaces*face_type::nbPtsPerFace );
2042  for ( uint16_type jj = 0; jj < npoints_per_face; ++jj )
2043  {
2044  pf.setPoint( jj, m->point( shift+it->point(jj).id() ) );
2045  }
2046  m->addFace( pf );
2047  };
2048  // add the faces
2049  size_type shift_f = std::distance( m1->beginFace(), m1->endFace() );
2050  for( auto it = m1->beginFace(), en = m1->endFace(); it != en; ++it )
2051  {
2052  addface( m, it, 0 );
2053  }
2054  for( auto it = m2->beginFace(), en = m2->endFace(); it != en; ++it )
2055  {
2056  // don't forget to shift the point id in the face
2057  addface( m, it, shift_p );
2058  }
2059  // add the elements
2060  auto addelement = [&]( boost::shared_ptr<Mesh<Shape, T, TheTag> > m, element_iterator it, size_type shift_f, size_type shift_p )
2061  {
2062  element_type pf = *it;
2063 
2064  static const uint16_type npoints_per_element = element_type::numPoints;
2065 
2066 #if 0
2067  for ( uint16_type jj = 0; jj < pf.numLocalFaces; ++jj )
2068  {
2069  pf.setFace( jj, m->face( shift_f+it->face(jj).id() ) );
2070  }
2071 #endif
2072  for ( uint16_type jj = 0; jj < npoints_per_element; ++jj )
2073  {
2074  pf.setPoint( jj, m->point( shift_p+it->point(jj).id() ) );
2075  }
2076  m->addElement( pf );
2077  };
2078  for( auto it = m1->beginElement(), en = m1->endElement(); it != en; ++it )
2079  {
2080  addelement( m, it, 0, 0 );
2081  }
2082  for( auto it = m2->beginElement(), en = m2->endElement(); it != en; ++it )
2083  {
2084  // don't forget to shift the point id in the face
2085  addelement( m, it, shift_f, shift_p );
2086  }
2087  m->components().set ( MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
2088  m->updateForUse();
2089  return m;
2090 }
2091 
2092 template<typename Shape, typename T, int Tag>
2093 void
2094 Mesh<Shape, T, Tag>::Inverse::distribute( bool extrapolation )
2095 {
2096  typename self_type::element_iterator el_it;
2097  typename self_type::element_iterator el_en;
2098  boost::tie( boost::tuples::ignore, el_it, el_en ) = Feel::elements( *M_mesh );
2099 
2100  typedef typename self_type::element_type element_type;
2101  typedef typename gm_type::template Context<vm::JACOBIAN|vm::KB|vm::POINT, element_type> gmc_type;
2102  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
2103  BoundingBox<> bb;
2104 
2105  typename gm_type::reference_convex_type refelem;
2106  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( M_mesh->gm(),
2107  refelem.points() ) );
2108  boost::unordered_map<size_type,bool> npt;
2109 
2110  M_dist.clear();
2111  M_ref_coords.clear();
2112  M_cvx_pts.clear();
2113  M_pts_cvx.clear();
2114  M_pts_cvx.resize( M_mesh->numElements() );
2115 
2116  KDTree::points_type boxpts;
2117  gmc_ptrtype __c( new gmc_type( M_mesh->gm(),
2118  *el_it,
2119  __geopc ) );
2120  VLOG(2) << "[Mesh::Inverse] distribute mesh points ion kdtree\n";
2121 
2122  for ( ; el_it != el_en; ++el_it )
2123  {
2124  // get geometric transformation
2125  __c->update( *el_it );
2126  gic_type gic( M_mesh->gm(), *el_it );
2127 
2128  // create bounding box
2129  //bb.make( el_it->points() );
2130  bb.make( el_it->G() );
2131 
2132  for ( size_type k=0; k < bb.min.size(); ++k )
2133  {
2134  bb.min[k] -= 1e-10;
2135  bb.max[k] += 1e-10;
2136  }
2137 
2138  DVLOG(2) << "G = " << el_it->G() << " min = " << bb.min << ", max = " << bb.max << "\n";
2139 
2140  // check if the points
2141  this->pointsInBox( boxpts, bb.min, bb.max );
2142 
2143  DVLOG(2) << "boxpts size = " << boxpts.size() << "\n";
2144 
2145 
2146  for ( size_type i = 0; i < boxpts.size(); ++i )
2147  {
2148  size_type index = boost::get<1>( boxpts[i] );
2149 
2150  if ( ( !npt[index] ) || M_dist[index] < 0 )
2151  {
2152  // check if we are in
2153  gic.setXReal( boost::get<0>( boxpts[i] ) );
2154  bool isin;
2155  value_type dmin;
2156  boost::tie( isin, dmin ) = refelem.isIn( gic.xRef() );
2157  bool tobeadded = extrapolation || isin;
2158 
2159  DVLOG(2) << "i = " << i << " index = " << index << " isin = " << ( isin >= -1e-10 )
2160  << " xref = " << gic.xRef() << " xreal = " << boost::get<0>( boxpts[i] )
2161  << " tobeadded= " << tobeadded << " dist=" << dmin<< "\n";
2162 
2163 
2164  if ( tobeadded && npt[index] )
2165  {
2166  if ( dmin > M_dist[index] )
2167  M_pts_cvx[M_cvx_pts[index]].erase( index );
2168 
2169  else
2170  tobeadded = false;
2171  }
2172 
2173  if ( tobeadded )
2174  {
2175  M_ref_coords[index] = gic.xRef();
2176  M_dist[index] = dmin;
2177  M_cvx_pts[index]=el_it->id();
2178  M_pts_cvx[el_it->id()][index]=boost::get<2>( boxpts[i] );
2179  npt[index]=true;
2180  }
2181  }
2182  }
2183 
2184  }
2185 
2186  VLOG(2) << "[Mesh::Inverse] distribute mesh points in kdtree done\n";
2187 }
2188 
2189 template<typename Shape, typename T, int Tag>
2190 void
2191 Mesh<Shape, T, Tag>::Localization::init()
2192 {
2193  if ( !M_mesh ) return;
2194 
2195  DLOG_IF( WARNING, IsInit == false ) << "You have already initialized the tool of localization\n";
2196 
2197 
2198  //clear data
2199  M_geoGlob_Elts.clear();
2200  M_kd_tree->clear();
2201 
2202  typename self_type::element_iterator el_it;
2203  typename self_type::element_iterator el_en;
2204  boost::tie( boost::tuples::ignore, el_it, el_en ) = Feel::elements( *M_mesh );
2205 
2206  for ( ; el_it != el_en; ++el_it )
2207  {
2208  for ( int i=0; i<el_it->nPoints(); ++i )
2209  {
2210  if ( boost::get<1>( M_geoGlob_Elts[el_it->point( i ).id()] ).size()==0 )
2211  {
2212  boost::get<0>( M_geoGlob_Elts[el_it->point( i ).id()] ) = el_it->point( i ).node();
2213  M_kd_tree->addPoint( el_it->point( i ).node(),el_it->point( i ).id() );
2214  }
2215 
2216  boost::get<1>( M_geoGlob_Elts[el_it->point( i ).id()] ).push_back( el_it->id() );
2217  }
2218  }
2219 
2220  IsInit=true;
2221  IsInitBoundaryFaces=false;
2222 
2223 }
2224 
2225 template<typename Shape, typename T, int Tag>
2226 void
2227 Mesh<Shape, T, Tag>::Localization::initBoundaryFaces()
2228 {
2229  if ( !M_mesh ) return;
2230 
2231 #if !defined( NDEBUG )
2232  FEELPP_ASSERT( IsInitBoundaryFaces == false )
2233  ( IsInitBoundaryFaces ).warn( "You have already initialized the tool of localization" );
2234 #endif
2235 
2236 
2237  //clear data
2238  M_geoGlob_Elts.clear();
2239  M_kd_tree->clear()
2240 ;
2241  typename self_type::location_face_iterator face_it;
2242  typename self_type::location_face_iterator face_en;
2243  boost::tie( boost::tuples::ignore, face_it, face_en ) = Feel::boundaryfaces( M_mesh );
2244 
2245  for ( ; face_it != face_en; ++face_it )
2246  {
2247  for ( int i=0; i<face_it->nPoints(); ++i )
2248  {
2249  if ( face_it->isConnectedTo0() )
2250  {
2251  if ( boost::get<1>( M_geoGlob_Elts[face_it->point( i ).id()] ).size()==0 )
2252  {
2253  boost::get<0>( M_geoGlob_Elts[face_it->point( i ).id()] ) = face_it->point( i ).node();
2254  M_kd_tree->addPoint( face_it->point( i ).node(),face_it->point( i ).id() );
2255  }
2256  boost::get<1>( M_geoGlob_Elts[face_it->point( i ).id()] ).push_back( face_it->element( 0 ).id() );
2257  }
2258  }
2259  }
2260 
2261  IsInitBoundaryFaces=true;
2262  IsInit=false;
2263 
2264 }
2265 
2266 
2267 
2268 template<typename Shape, typename T, int Tag>
2269 boost::tuple<bool,typename Mesh<Shape, T, Tag>::node_type,double>
2270 Mesh<Shape, T, Tag>::Localization::isIn( size_type _id, const node_type & _pt ) const
2271 {
2272  bool isin=false;
2273  double dmin;
2274  node_type x_ref;
2275 
2276  //get element with the id
2277  auto const& elt = M_mesh->element( _id );
2278 
2279  if ( elt.isOnBoundary() )
2280  {
2281  // get inverse geometric transformation
2282  gmc_inverse_type gic( M_mesh->gm(), elt, this->mesh()->worldComm().subWorldCommSeq() );
2283  //apply the inverse geometric transformation for the point p
2284  gic.setXReal( _pt);
2285  x_ref=gic.xRef();
2286  // the point is in the reference element ?
2287  boost::tie( isin, dmin ) = M_refelem.isIn( gic.xRef() );
2288  }
2289  else
2290  {
2291  // get inverse geometric transformation
2292  gmc1_inverse_type gic( M_mesh->gm1(), elt, mpl::int_<1>(), this->mesh()->worldComm().subWorldCommSeq() );
2293  //apply the inverse geometric transformation for the point p
2294  gic.setXReal( _pt);
2295  x_ref=gic.xRef();
2296  // the point is in the reference element ?
2297  boost::tie( isin, dmin ) = M_refelem1.isIn( gic.xRef() );
2298  }
2299 
2300  return boost::make_tuple(isin,x_ref,dmin);
2301 }
2302 
2303 template<typename Shape, typename T, int Tag>
2304 boost::tuple<uint16_type,std::vector<bool> >
2305 Mesh<Shape, T, Tag>::Localization::isIn( std::vector<size_type> _ids, const node_type & _pt )
2306 {
2307  typedef typename self_type::gm_type::reference_convex_type ref_convex_type;
2308  typedef typename self_type::gm1_type::reference_convex_type ref_convex1_type;
2309 
2310  uint16_type nbId = _ids.size();
2311  std::vector<bool> isin( _ids.size(),false );
2312  bool isin2=false;
2313  double dmin;
2314  node_type __x_ref;
2315 
2316  uint16_type nbIsIn=0;
2317 
2318  for ( uint16_type i = 0; i< nbId ; ++i )
2319  {
2320  //get element with the id
2321  auto const& elt = M_mesh->element( _ids[i] );
2322 
2323  if ( elt.isOnBoundary() )
2324  {
2325  // get inverse geometric transformation
2326  gmc_inverse_type gic( M_mesh->gm(), elt );
2327  //apply the inverse geometric transformation for the point p
2328  gic.setXReal( _pt);
2329  __x_ref=gic.xRef();
2330  // the point is in the reference element ?
2331  boost::tie( isin2, dmin ) = M_refelem.isIn( gic.xRef() );
2332  isin[i] = isin2;
2333  }
2334  else
2335  {
2336  // get inverse geometric transformation
2337  gmc1_inverse_type gic( M_mesh->gm1(), elt, mpl::int_<1>() );
2338  //apply the inverse geometric transformation for the point p
2339  gic.setXReal( _pt);
2340  __x_ref=gic.xRef();
2341  // the point is in the reference element ?
2342  boost::tie( isin2, dmin ) = M_refelem1.isIn( gic.xRef() );
2343  isin[i] = isin2;
2344  }
2345  if (isin[i]) ++nbIsIn;
2346  }
2347 
2348  return boost::make_tuple( nbIsIn,isin );
2349 }
2350 
2351 
2352 template<typename Shape, typename T, int Tag>
2353 boost::tuple<bool, size_type, typename Mesh<Shape, T, Tag>::node_type>
2354 Mesh<Shape, T, Tag>::Localization::searchElement( const node_type & p )
2355 {
2356 
2357 #if !defined( NDEBUG )
2358  FEELPP_ASSERT( IsInit == true )
2359  ( IsInit ).warn( "You don't have initialized the tool of localization" );
2360 #endif
2361  bool isin=false;double dmin=0;
2362  node_type x_ref;
2363  size_type idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2364 
2365  std::list< std::pair<size_type, uint> > ListTri;
2366  this->searchInKdTree(p,ListTri);
2367 
2368  //research the element which contains the point p
2369  auto itLT=ListTri.begin();
2370  auto itLT_end=ListTri.end();
2371 
2372 #if !defined( NDEBUG )
2373  //if(std::distance(itLT,itLT_end)==0) std::cout<<"\nListTri vide\n";
2374  FEELPP_ASSERT( std::distance( itLT,itLT_end )>0 ).error( " problem in list localization : is empty" );
2375 #endif
2376 
2377  while ( itLT != itLT_end && !isin )
2378  {
2379  //get element with the id
2380  //elt = M_mesh->element( itLT->first );
2381 
2382  // search point in this elt
2383  boost::tie(isin,x_ref,dmin) = this->isIn(itLT->first,p);
2384  // if not inside, continue the research with an other element
2385  if (!isin) ++itLT;
2386  else idEltFound=itLT->first;
2387  }
2388 
2389 
2390  if (!isin)
2391  {
2392  if( this->doExtrapolation() )
2393  {
2394  //std::cout << "WARNING EXTRAPOLATION for the point" << p << std::endl;
2395  //std::cout << "W";
2396  auto const& eltUsedForExtrapolation = this->mesh()->element(ListTri.begin()->first);
2397  gmc_inverse_type gic( this->mesh()->gm(), eltUsedForExtrapolation, this->mesh()->worldComm().subWorldCommSeq() );
2398  //apply the inverse geometric transformation for the point p
2399  gic.setXReal( p);
2400  boost::tie(isin,idEltFound,x_ref) = boost::make_tuple(true,eltUsedForExtrapolation.id(),gic.xRef());
2401  }
2402  else
2403  {
2404  idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2405  isin = false;
2406  //x_ref=?
2407  }
2408  }
2409 
2410  return boost::make_tuple( isin, idEltFound, x_ref);
2411 
2412 }
2413 
2414 template<typename Shape, typename T, int Tag>
2415 boost::tuple<std::vector<bool>, size_type>
2416 Mesh<Shape, T, Tag>::Localization::run_analysis( const matrix_node_type & m,
2417  const size_type & eltHypothetical )
2418 {
2419 
2420 #if !defined( NDEBUG )
2421  FEELPP_ASSERT( IsInit == true )
2422  ( IsInit ).warn( "You don't have initialized the tool of localization" );
2423 #endif
2424 
2425  bool find_x;
2426  size_type cv_id=eltHypothetical;
2427  node_type x_ref;double dmin=0;
2428  std::vector<bool> hasFindPts(m.size2(),false);
2429 
2430  M_resultAnalysis.clear();
2431 
2432  bool doExtrapolationAtStart = this->doExtrapolation();
2433  auto nPtMaxNearNeighborAtStart = this->kdtree()->nPtMaxNearNeighbor();
2434 
2435 
2436  // first step : no extrapolation
2437  if ( doExtrapolationAtStart ) this->setExtrapolation( false );
2438 
2439  // first currentEltHypothetical
2440  auto currentEltHypothetical = eltHypothetical;//this->mesh()->beginElement()->id();//eltHypothetical;
2441  for ( size_type i=0; i< m.size2(); ++i )
2442  {
2443  bool testHypothetical_find = false;
2444 
2445  if ( eltHypothetical!=invalid_size_type_value )
2446  {
2447  boost::tie( testHypothetical_find,x_ref,dmin ) = this->isIn( currentEltHypothetical,ublas::column( m, i ) );
2448  }
2449  if ( testHypothetical_find )
2450  {
2451  cv_id = currentEltHypothetical;
2452  M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2453  hasFindPts[i]=true;
2454  }
2455  else // search kdtree
2456  {
2457  // kdtree call
2458  boost::tie( find_x, cv_id, x_ref ) = this->searchElement(ublas::column( m, i ));
2459  // traitement
2460  if (find_x) // if find : OK
2461  {
2462  M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2463  currentEltHypothetical = cv_id;
2464  hasFindPts[i]=true;
2465  }
2466  else// if (false) // try an other method (no efficient but maybe a solution)
2467  {
2468  // search in all element
2469  this->kdtree()->nbNearNeighbor(2*nPtMaxNearNeighborAtStart /*15*/ /*this->mesh()->numElements()*/);
2470 
2471  boost::tie( find_x, cv_id, x_ref ) = this->searchElement(ublas::column( m, i ));
2472  //revert parameter
2473  this->kdtree()->nbNearNeighbor(nPtMaxNearNeighborAtStart);
2474  // if find : OK (but strange!)
2475  if (find_x)
2476  {
2477  M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2478  currentEltHypothetical = cv_id;
2479  hasFindPts[i]=true;
2480  }
2481  else if (doExtrapolationAtStart)// && this->mesh()->worldComm().localSize()==1)
2482  {
2483  this->setExtrapolation(true);
2484 
2485  boost::tie( find_x, cv_id, x_ref ) = this->searchElement(ublas::column( m, i ));
2486  // normaly is find
2487  if (find_x)
2488  {
2489  M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2490  currentEltHypothetical = cv_id;
2491  hasFindPts[i]=true;
2492  }
2493 
2494  this->setExtrapolation(false);
2495  }
2496  else
2497  {
2498  //std::cout<<"\n Il y a un GROS Probleme de Localization\n";
2499  }
2500  }
2501  } // search kdtree
2502 
2503  } // for (size_type i=0;i< m.size2();++i)
2504 
2505  //revert parameter
2506  this->setExtrapolation(doExtrapolationAtStart);
2507 
2508  return boost::make_tuple(hasFindPts,cv_id);
2509 
2510 } //run_analysis
2511 
2512 
2513 
2514 template<typename Shape, typename T, int Tag>
2515 boost::tuple<bool, std::list<boost::tuple<size_type, typename Mesh<Shape, T, Tag>::node_type> > >
2516 Mesh<Shape, T, Tag>::Localization::searchElements( const node_type & p )
2517 {
2518 
2519 #if !defined( NDEBUG )
2520  FEELPP_ASSERT( IsInit == true )
2521  ( IsInit ).warn( "You don't have initialized the tool of localization" );
2522 #endif
2523 
2524  //this->kdtree()->nbNearNeighbor(this->mesh()->numElements());
2525 
2526  std::list< std::pair<size_type, uint> > ListTri;
2527  searchInKdTree( p,ListTri );
2528 
2529  typename self_type::element_type elt;
2530  typename self_type::gm_type::reference_convex_type refelem;
2531  typename self_type::gm1_type::reference_convex_type refelem1;
2532 
2533  bool isin=false;
2534  double dmin;
2535  node_type x_ref;
2536 
2537  //research the element which contains the point p
2538  auto itLT=ListTri.begin();
2539  auto itLT_end=ListTri.end();
2540 
2541 #if !defined( NDEBUG )
2542  //if(std::distance(itLT,itLT_end)==0) std::cout<<"\nListTri vide\n";
2543  FEELPP_ASSERT( std::distance( itLT,itLT_end )>0 ).error( " problem in list localization : is empty" );
2544 #endif
2545 
2546  std::list<boost::tuple<size_type,node_type> > newlistelts;
2547  newlistelts.clear();
2548  bool find = false;
2549  bool finishSearch=false;
2550  while ( itLT != itLT_end && !finishSearch /*&& !isin*/ )
2551  {
2552 #if 0
2553  //get element with the id
2554  elt= M_mesh->element( itLT->first );
2555 
2556  if ( elt.isOnBoundary() )
2557  {
2558  // get inverse geometric transformation
2559  typename self_type::Inverse::gic_type gic( M_mesh->gm(), elt );
2560 
2561  //apply the inverse geometric transformation for the point p
2562  gic.setXReal( p );
2563  __x_ref=gic.xRef();
2564 
2565  // the point is in the reference element ?
2566  boost::tie( isin, dmin ) = refelem.isIn( gic.xRef() );
2567  }
2568 
2569  else
2570  {
2571  // get inverse geometric transformation
2572  typename self_type::Inverse::gic1_type gic( M_mesh->gm1(), elt,mpl::int_<1>() );
2573 
2574  //apply the inverse geometric transformation for the point p
2575  gic.setXReal( p );
2576  __x_ref=gic.xRef();
2577 
2578  // the point is in the reference element ?
2579  boost::tie( isin, dmin ) = refelem1.isIn( gic.xRef() );
2580  //std::cout << "gic.xRef()" << gic.xRef() << std::endl;
2581  }
2582 #else
2583  boost::tie(isin,x_ref, dmin) = this->isIn(itLT->first,p);
2584 #endif
2585  if ( isin )
2586  {
2587  newlistelts.push_back( boost::make_tuple( itLT->first,x_ref ) );
2588  find = true;
2589  // if not on boundary -> finish for this point
2590  if (dmin>1e-7) finishSearch=true;
2591  }
2592 
2593  //if (find) std::cout << elt.G() << std::endl;
2594 
2595  //if not inside, continue the research with an other element
2596  //if (!isin) ++itLT;
2597  ++itLT;
2598  }
2599 
2600  if ( !find ) std::cout << "\n WARNING EXTRAPOLATION IN SEARCHELEMENTS!!!"<<std::endl;
2601 
2602  if ( find )
2603  return boost::make_tuple( true,newlistelts );
2604 
2605  else if ( !find && !M_doExtrapolation )
2606  return boost::make_tuple( false,newlistelts );
2607 
2608  else
2609  {
2610  //std::cout << "\n WARNING EXTRAPOLATION \n";
2611  itLT=ListTri.begin();
2612  elt= M_mesh->element( itLT->first );
2613  typename self_type::Inverse::gic_type gic( M_mesh->gm(), elt );
2614  //apply the inverse geometric transformation for the point p
2615  //gic.setXReal(boost::get<0>(*ptsNN.begin()));
2616  gic.setXReal( p );
2617  x_ref=gic.xRef();
2618  //return boost::make_tuple( true, itLT->first, __x_ref);
2619  newlistelts.push_back( boost::make_tuple( itLT->first,x_ref ) );
2620  find = true;
2621  return boost::make_tuple( true,newlistelts );
2622  }
2623 } // searchElements
2624 
2625 
2626 template<typename Shape, typename T, int Tag>
2627 void
2628 Mesh<Shape, T, Tag>::Localization::searchInKdTree( const node_type & p,
2629  std::list< std::pair<size_type, uint> > & ListTri )
2630 {
2631  //search for nearest points
2632  M_kd_tree->search( p );
2633 
2634  //get the results of research
2635  typename KDTree::points_search_type ptsNN = M_kd_tree->pointsNearNeighbor();
2636 
2637  typename KDTree::points_search_const_iterator itNN = ptsNN.begin();
2638  typename KDTree::points_search_const_iterator itNN_end = ptsNN.end();
2639 
2640 #if !defined( NDEBUG )
2641  FEELPP_ASSERT( std::distance( itNN,itNN_end )>0 ).error( "none Near Neighbor Points are find" );
2642 #endif
2643 
2644  //iterator on a l(ist index element
2645  typename std::list<size_type>::iterator itL;
2646  typename std::list<size_type>::iterator itL_end;
2647 
2648  //ListTri will contain the indices of elements (size_type)
2649  //and the number of occurence(uint)
2650  //std::list< std::pair<size_type, uint> > ListTri;
2651  std::list< std::pair<size_type, uint> >::iterator itLT;
2652  std::list< std::pair<size_type, uint> >::iterator itLT_end;
2653 
2654  //create of ListTri : sort largest to smallest occurrences
2655  //In case of equality : if the point is closer than another then it will be before
2656  // if it is in the same point then the lowest index will be before
2657  for ( ; itNN != itNN_end; ++itNN )
2658  {
2659  itL= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).begin();
2660  itL_end= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).end();
2661 
2662  for ( ; itL != itL_end; ++itL )
2663  {
2664  itLT=ListTri.begin();
2665  itLT_end=ListTri.end();
2666  bool find=false;
2667 
2668  while ( itLT != itLT_end && !find )
2669  {
2670  if ( itLT->first == *itL ) find=true;
2671 
2672  else ++itLT;
2673  }
2674 
2675  if ( find )
2676  {
2677  uint nb=itLT->second+1;
2678  size_type numEl=itLT->first;
2679  ListTri.remove( *itLT );
2680  itLT=ListTri.begin();
2681  itLT_end=ListTri.end();
2682  bool find=false;
2683 
2684  while ( itLT != itLT_end && !find )
2685  {
2686  if ( itLT->second < nb ) find=true;
2687 
2688  else ++itLT;
2689  }
2690 
2691  ListTri.insert( itLT,std::make_pair( numEl,nb ) );
2692  }
2693 
2694  else ListTri.push_back( std::make_pair( *itL,1 ) );
2695  }
2696  }
2697 
2698 }
2699 
2700 template<typename Shape, typename T, int Tag>
2701 boost::tuple<bool,typename Mesh<Shape, T, Tag>::node_type,double>
2702 Mesh<Shape, T, Tag>::Localization::isIn( size_type _id,
2703  const node_type & _pt,
2704  const matrix_node_type & setPoints,
2705  mpl::int_<1> ) const
2706 {
2707  bool isin=true; // warning : start with true
2708  double dmin=0.;
2709  node_type x_ref;
2710 
2711  //get element with the id
2712  auto const& elt= M_mesh->element( _id );
2713  auto const& eltG = elt.G();
2714 
2715  // check conformity between setPoints (given) and eltG (localize)
2716  std::vector<bool> find( setPoints.size2(),false );
2717  for ( size_type i=0; i< setPoints.size2() && isin; ++i )
2718  {
2719  auto const& thePt = ublas::column( setPoints,i );
2720 
2721  for ( size_type j=0; j<eltG.size2(); ++j )
2722  {
2723  auto ptjeltG = ublas::column( eltG,j );
2724 
2725  if ( ptjeltG.size()==1 )
2726  {
2727  if ( std::abs( thePt( 0 )-ptjeltG( 0 ) )<1e-8 )
2728  find[i]=true;
2729  }
2730  else if ( ptjeltG.size()==2 )
2731  {
2732  if ( std::abs( thePt( 0 )-ptjeltG( 0 ) )<1e-8 &&
2733  std::abs( thePt( 1 )-ptjeltG( 1 ) )<1e-8 )
2734  find[i]=true;
2735  }
2736  else if ( ptjeltG.size()==3 )
2737  {
2738  if ( std::abs( thePt( 0 )-ptjeltG( 0 ) )<1e-8 &&
2739  std::abs( thePt( 1 )-ptjeltG( 1 ) )<1e-8 &&
2740  std::abs( thePt( 2 )-ptjeltG( 2 ) )<1e-8 )
2741  find[i]=true;
2742  }
2743  }
2744  // up checking
2745  isin &= find[i];
2746  }
2747 
2748  // if find -> get ref point and check
2749  if ( isin )
2750  {
2751  bool isin2=false;
2752  boost::tie(isin2,x_ref,dmin) = this->isIn(_id,_pt);
2753  LOG_IF(WARNING, !isin2) << "Mesh::Localization::isIn<Conformal> : check fail -> maybe x_ref is not correct";
2754  }
2755 
2756  return boost::make_tuple(isin,x_ref,dmin);
2757 }
2758 
2759 
2760 
2761 template<typename Shape, typename T, int Tag>
2762 boost::tuple<bool, size_type, typename Mesh<Shape, T, Tag>::node_type>
2763 Mesh<Shape, T, Tag>::Localization::searchElement( const node_type & p,
2764  const matrix_node_type & setPoints,
2765  mpl::int_<1> )
2766 {
2767  typename self_type::element_type elt;
2768  typename self_type::gm_type::reference_convex_type refelem;
2769  typename self_type::gm1_type::reference_convex_type refelem1;
2770 
2771  bool isin=false,isin2=false;
2772  double dmin=0.;
2773  node_type x_ref;
2774  size_type idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2775 
2776  std::list< std::pair<size_type, uint> > ListTri;
2777  searchInKdTree( p,ListTri );
2778 
2779  auto itLT=ListTri.begin();
2780  auto itLT_end=ListTri.end();
2781 
2782 #if !defined( NDEBUG )
2783  //if(std::distance(itLT,itLT_end)==0) std::cout<<"\nListTri vide\n";
2784  FEELPP_ASSERT( std::distance( itLT,itLT_end )>0 ).error( " problem in list localization : is empty" );
2785 #endif
2786 
2787 
2788  //research the element which contains the point p
2789  while ( itLT != itLT_end && !isin )
2790  {
2791  boost::tie(isin,x_ref,dmin) = this->isIn(itLT->first,p,setPoints,mpl::int_<1>());
2792 
2793  //if not inside, continue the research with an other element
2794  if ( !isin ) ++itLT;
2795  else idEltFound=itLT->first;
2796  } //while ( itLT != itLT_end && !isin )
2797 
2798  if (!isin)
2799  {
2800  if( this->doExtrapolation() )
2801  {
2802  std::cout << "WARNING EXTRAPOLATION for the point" << p << std::endl;
2803  //std::cout << "W";
2804  auto const& eltUsedForExtrapolation = this->mesh()->element(ListTri.begin()->first);
2805  gmc_inverse_type gic( this->mesh()->gm(), eltUsedForExtrapolation, this->mesh()->worldComm().subWorldCommSeq() );
2806  //apply the inverse geometric transformation for the point p
2807  gic.setXReal( p);
2808  boost::tie(isin,idEltFound,x_ref) = boost::make_tuple(true,eltUsedForExtrapolation.id(),gic.xRef());
2809  }
2810  else
2811  {
2812  idEltFound = this->mesh()->beginElementWithId(this->mesh()->worldComm().localRank())->id();
2813  isin = false;
2814  //x_ref=?
2815  }
2816  }
2817 
2818  return boost::make_tuple( isin, idEltFound, x_ref);
2819 
2820 } //searchElement
2821 
2822 
2823 template<typename Shape, typename T, int Tag>
2824 boost::tuple<std::vector<bool>, size_type>
2825 Mesh<Shape, T, Tag>::Localization::run_analysis( const matrix_node_type & m,
2826  const size_type & eltHypothetical,
2827  const matrix_node_type & setPoints,
2828  mpl::int_<1> )
2829 {
2830 #if !defined( NDEBUG )
2831  FEELPP_ASSERT( IsInit == true )
2832  ( IsInit ).warn( "You don't have initialized the tool of localization" );
2833 #endif
2834 
2835  bool find_x=false;
2836  size_type cv_id=eltHypothetical;
2837  node_type x_ref;
2838  double dmin;
2839  std::vector<bool> hasFindPts(setPoints.size2(),false);
2840 
2841  M_resultAnalysis.clear();
2842  auto currentEltHypothetical = eltHypothetical;
2843  for ( size_type i=0; i< m.size2(); ++i )
2844  {
2845 
2846  bool testHypothetical_find = false;
2847  if ( eltHypothetical!=invalid_size_type_value )
2848  {
2849  boost::tie( testHypothetical_find,x_ref,dmin ) = this->isIn( currentEltHypothetical,ublas::column( m, i ),setPoints,mpl::int_<1>() );
2850  }
2851  if ( testHypothetical_find )
2852  {
2853  cv_id = currentEltHypothetical;
2854  M_resultAnalysis[cv_id].push_back( boost::make_tuple(i,x_ref) );
2855  hasFindPts[i]=true;
2856  }
2857  else
2858  {
2859  boost::tie( find_x, cv_id, x_ref ) = this->searchElement( ublas::column( m, i ),setPoints,mpl::int_<1>() );
2860 
2861  if ( find_x )
2862  {
2863  M_resultAnalysis[cv_id].push_back( boost::make_tuple( i,x_ref ) );
2864  hasFindPts[i]=true;
2865  currentEltHypothetical = cv_id;
2866  }
2867  }
2868  //else std::cout<<"\nNew Probleme Localization\n" << std::endl;
2869  }
2870 
2871  return boost::make_tuple(hasFindPts,cv_id);
2872 
2873 } // run_analysis
2874 
2875 
2876 template<typename Shape, typename T, int Tag>
2877 typename Mesh<Shape, T, Tag>::node_type
2878 Mesh<Shape, T, Tag>::Localization::barycenter() const
2879 {
2880  return this->barycenter(mpl::int_<nRealDim>());
2881 }
2882 
2883 template<typename Shape, typename T, int Tag>
2884 typename Mesh<Shape, T, Tag>::node_type
2885 Mesh<Shape, T, Tag>::Localization::barycenter(mpl::int_<1> ) const
2886 {
2887  node_type res(1);
2888  res(0)=0;
2889  for (size_type i = 0 ; i<this->kdtree()->nPoints() ; ++i)
2890  {
2891  auto const& pt = this->kdtree()->points()[i].template get<0>();
2892  res(0)+=pt(0);
2893  }
2894  res(0)/=this->kdtree()->nPoints();
2895  return res;
2896 }
2897 template<typename Shape, typename T, int Tag>
2898 typename Mesh<Shape, T, Tag>::node_type
2899 Mesh<Shape, T, Tag>::Localization::barycenter(mpl::int_<2> ) const
2900 {
2901  node_type res(2);
2902  res(0)=0;res(1)=0;
2903  for (size_type i = 0 ; i<this->kdtree()->nPoints() ; ++i)
2904  {
2905  auto const& pt = this->kdtree()->points()[i].template get<0>();
2906  res(0)+=pt(0);res(1)+=pt(1);
2907  }
2908  res(0)/=this->kdtree()->nPoints();res(1)/=this->kdtree()->nPoints();
2909  return res;
2910 }
2911 template<typename Shape, typename T, int Tag>
2912 typename Mesh<Shape, T, Tag>::node_type
2913 Mesh<Shape, T, Tag>::Localization::barycenter(mpl::int_<3> ) const
2914 {
2915  node_type res(3);
2916  res(0)=0;res(1)=0;res(2)=0;
2917  for (size_type i = 0 ; i<this->kdtree()->nPoints() ; ++i)
2918  {
2919  auto const& pt = this->kdtree()->points()[i].template get<0>();
2920  res(0)+=pt(0);res(1)+=pt(1);res(2)+=pt(2);
2921  }
2922  res(0)/=this->kdtree()->nPoints();res(1)/=this->kdtree()->nPoints();res(2)/=this->kdtree()->nPoints();
2923  return res;
2924 }
2925 
2926 
2927 
2928 
2929 
2930 #if 0
2931 #if defined( FEELPP_INSTANTIATION_MODE )
2932 
2933 # define DIMS1 BOOST_PP_TUPLE_TO_LIST(1,(1))
2934 # define RDIMS1 BOOST_PP_TUPLE_TO_LIST(3,(1,2,3))
2935 # define ORDERS1 BOOST_PP_TUPLE_TO_LIST(5,(1,2,3,4,5))
2936 
2937 # define DIMS2 BOOST_PP_TUPLE_TO_LIST(1,(2))
2938 # define RDIMS2 BOOST_PP_TUPLE_TO_LIST(2,(2,3))
2939 # define ORDERS2 BOOST_PP_TUPLE_TO_LIST(5,(1,2,3,4,5))
2940 
2941 # define DIMS3 BOOST_PP_TUPLE_TO_LIST(1,(3))
2942 # define RDIMS3 BOOST_PP_TUPLE_TO_LIST(1,(3))
2943 # define ORDERS3 BOOST_PP_TUPLE_TO_LIST(4,(1,2,3,4))
2944 
2945 # define FACTORY_SIMPLEX(LDIM,LORDER,RDIM) template class Mesh<Simplex<LDIM, LORDER, RDIM> >;
2946 # define FACTORY_HYPERCUBE(LDIM,LORDER,RDIM) template class Mesh<Hypercube<LDIM, LORDER, RDIM> >;
2947 
2948 # define FACTORY_SIMPLEX_OP(_, GDO) FACTORY_SIMPLEX GDO
2949 # define FACTORY_HYPERCUBE_OP(_, GDO) FACTORY_HYPERCUBE GDO
2950 
2951 
2952 # define FACTORY_SIMPLEX_E(LDIM,LORDER,RDIM) extern template class Mesh<Simplex<LDIM, LORDER, RDIM> >;
2953 # define FACTORY_HYPERCUBE_E(LDIM,LORDER,RDIM) extern template class Mesh<Hypercube<LDIM, LORDER, RDIM> >;
2954 
2955 # define FACTORY_SIMPLEX_OP_E(_, GDO) FACTORY_HYPERCUBE GDO
2956 # define FACTORY_HYPERCUBE_OP_E(_, GDO) FACTORY_HYPERCUBE_E GDO
2957 
2958 #if !defined( FEELPP_MESH_IMPL_NOEXTERN )
2959 
2960 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_SIMPLEX_OP_E, 3, ( DIMS1, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS1 ), RDIMS1 ) )
2961 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_HYPERCUBE_OP_E, 3, ( DIMS1, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS1 ), RDIMS1 ) )
2962 
2963 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_SIMPLEX_OP_E, 3, ( DIMS2, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS2 ), RDIMS2 ) )
2964 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_HYPERCUBE_OP_E, 3, ( DIMS2, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS2 ), RDIMS2 ) )
2965 
2966 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_SIMPLEX_OP_E, 3, ( DIMS3, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS3 ), RDIMS3 ) )
2967 BOOST_PP_LIST_FOR_EACH_PRODUCT( FACTORY_HYPERCUBE_OP_E, 3, ( DIMS3, BOOST_PP_LIST_FIRST_N( FEELPP_MESH_MAX_ORDER, ORDERS3 ), RDIMS3 ) )
2968 
2969 #endif // FEELPP_MESH_IMPL_NOEXTERN
2970 #endif // FEELPP_INSTANTIATION_MODE
2971 #endif
2972 } // namespace Feel
2973 
2974 
2975 
2976 
2977 
2978 #endif // __MESHIMPL_HPP

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