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
mesh.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: 2005-07-05
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2006-2012 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __mesh_H
31 #define __mesh_H 1
32 
33 #include <boost/unordered_map.hpp>
34 
35 #include <boost/foreach.hpp>
36 #include <boost/signals2/signal.hpp>
37 #include <boost/serialization/vector.hpp>
38 #include <boost/serialization/array.hpp>
39 #include <boost/serialization/base_object.hpp>
40 #include <boost/archive/text_oarchive.hpp>
41 #include <boost/archive/text_iarchive.hpp>
42 #include <boost/archive/binary_oarchive.hpp>
43 #include <boost/archive/binary_iarchive.hpp>
44 
45 #include <feel/feelcore/context.hpp>
46 //#include <feel/feelcore/worldcomm.hpp>
47 
48 #include <feel/feelmesh/mesh0d.hpp>
49 #include <feel/feelmesh/mesh1d.hpp>
50 #include <feel/feelmesh/mesh2d.hpp>
51 #include <feel/feelmesh/mesh3d.hpp>
52 #include <feel/feelmesh/meshutil.hpp>
54 #include <feel/feelmesh/enums.hpp>
55 
56 #include <feel/feelpoly/geomap.hpp>
59 
60 
61 
62 
63 #include <boost/preprocessor/comparison/less.hpp>
64 #include <boost/preprocessor/logical/and.hpp>
65 #include <boost/preprocessor/control/if.hpp>
66 #include <boost/preprocessor/list/at.hpp>
67 #include <boost/preprocessor/list/cat.hpp>
68 #include <boost/preprocessor/list/for_each_product.hpp>
69 #include <boost/preprocessor/logical/or.hpp>
70 #include <boost/preprocessor/tuple/to_list.hpp>
71 #include <boost/preprocessor/tuple/eat.hpp>
72 #include <boost/preprocessor/facilities/empty.hpp>
73 #include <boost/preprocessor/punctuation/comma.hpp>
74 #include <boost/preprocessor/facilities/identity.hpp>
75 
76 #include <boost/enable_shared_from_this.hpp>
77 
78 namespace Feel
79 {
80 const size_type EXTRACTION_KEEP_POINTS_IDS = ( 1<<0 );
81 const size_type EXTRACTION_KEEP_EDGES_IDS = ( 1<<1 );
82 const size_type EXTRACTION_KEEP_FACES_IDS = ( 1<<2 );
83 const size_type EXTRACTION_KEEP_VOLUMES_IDS = ( 1<<3 );
84 const size_type EXTRACTION_KEEP_ALL_IDS = ( EXTRACTION_KEEP_POINTS_IDS |
85  EXTRACTION_KEEP_EDGES_IDS |
86  EXTRACTION_KEEP_FACES_IDS |
87  EXTRACTION_KEEP_VOLUMES_IDS );
88 const size_type EXTRACTION_KEEP_MESH_RELATION = ( 1<<4 );
89 }
91 
92 namespace Feel
93 {
94 
95 struct MeshMarkerName
96 {
97  std::string name;
98  std::vector<int> ids;
99 
100 };
101 
102 std::vector<MeshMarkerName> markerMap( int Dim );
103 
104 
105 
109 template<typename Mesh> class Partitioner;
110 
117 // template <typename GeoShape, typename T > class Mesh;
118 
119 template <typename GeoShape, typename T = double, int Tag = 0>
120 class Mesh
121  :
122 public mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<0> >,
123  mpl::identity<Mesh0D<GeoShape > >,
124  typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<1> >,
125  mpl::identity<Mesh1D<GeoShape > >,
126  typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<2> >,
127  mpl::identity<Mesh2D<GeoShape> >,
128  mpl::identity<Mesh3D<GeoShape> > >::type>::type>::type::type,
129 public boost::addable<Mesh<GeoShape,T,Tag> >,
130 public boost::enable_shared_from_this< Mesh<GeoShape,T,Tag> >
131 {
132  typedef typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<0> >,
133  mpl::identity<Mesh0D<GeoShape> >,
134  typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<1> >,
135  mpl::identity<Mesh1D<GeoShape> >,
136  typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<2> >,
137  mpl::identity<Mesh2D<GeoShape> >,
138  mpl::identity<Mesh3D<GeoShape> > >::type>::type>::type::type super;
139 public:
140 
141 
145 
146  static const uint16_type nDim = GeoShape::nDim;
147  static const uint16_type nRealDim = GeoShape::nRealDim;
148  static const uint16_type Shape = GeoShape::Shape;
149  static const uint16_type nOrder = GeoShape::nOrder;
150  static const uint16_type tag = Tag;
151 
153 
156  typedef T value_type;
157  typedef GeoShape shape_type;
158  typedef typename super::return_type return_type;
159  typedef typename node<double>::type node_type;
160 
161  typedef typename super::super_elements super_elements;
162  typedef typename super::elements_type elements_type;
163  typedef typename super::element_type element_type;
164  typedef typename super::element_iterator element_iterator;
165  typedef typename super::element_const_iterator element_const_iterator;
166 
167  typedef typename super::super_faces super_faces;
168  typedef typename super::faces_type faces_type;
169  typedef typename super::face_type face_type;
170  typedef typename super::face_iterator face_iterator;
171  typedef typename super::face_const_iterator face_const_iterator;
172 
173  typedef typename super::edge_type edge_type;
174 
175  typedef typename super::points_type points_type;
176  typedef typename super::point_type point_type;
177  typedef typename super::point_iterator point_iterator;
178  typedef typename super::point_const_iterator point_const_iterator;
179 
180  typedef typename element_type::gm_type gm_type;
181  typedef boost::shared_ptr<gm_type> gm_ptrtype;
182 
183  typedef typename element_type::gm1_type gm1_type;
184  typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
185 
186  template<size_type ContextID>
187  struct gmc
188  {
189  typedef typename gm_type::template Context<ContextID, element_type> type;
190  typedef boost::shared_ptr<type> ptrtype;
191  };
192 
193  typedef Mesh<shape_type, T, Tag> self_type;
194  typedef self_type mesh_type;
195  typedef boost::shared_ptr<self_type> self_ptrtype;
196  typedef self_ptrtype mesh_ptrtype;
197 
198  typedef typename element_type::template reference_convex<T>::type reference_convex_type;
199 
200  //typedef Partitioner<self_type> partitioner_type;
201  //typedef boost::shared_ptr<partitioner_type> partitioner_ptrtype;
202 
203  typedef typename super::face_processor_type face_processor_type;
204  typedef typename super::face_processor_type element_edge_type;
205 
206  typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
207  mpl::identity< Mesh< Simplex< GeoShape::nDim,1,GeoShape::nRealDim>, value_type, Tag > >,
208  mpl::identity< Mesh< Hypercube<GeoShape::nDim,1,GeoShape::nRealDim>,value_type, Tag > > >::type::type P1_mesh_type;
209 
210  typedef boost::shared_ptr<P1_mesh_type> P1_mesh_ptrtype;
211 
212  template<int TheTag>
213  struct trace_mesh
214  {
215  typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
216  mpl::identity< Mesh< Simplex< GeoShape::nDim-1,nOrder,GeoShape::nRealDim>, value_type, TheTag > >,
217  mpl::identity< Mesh< Hypercube<GeoShape::nDim-1,nOrder,GeoShape::nRealDim>,value_type, TheTag > > >::type::type type;
218  typedef boost::shared_ptr<type> ptrtype;
219  typedef boost::shared_ptr<const type> const_ptrtype;
220  };
221  typedef typename trace_mesh<Tag>::type trace_mesh_type;
222  typedef typename trace_mesh<Tag>::ptrtype trace_mesh_ptrtype;
223 
224 
225  template<int TheTag>
226  struct trace_trace_mesh
227  {
228  static const uint16_type nDim = (GeoShape::nDim==1)?GeoShape::nDim-1:GeoShape::nDim-2;
229  typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
230  mpl::identity< Mesh< Simplex<nDim,nOrder,GeoShape::nRealDim>, value_type, TheTag > >,
231  mpl::identity< Mesh< Hypercube<nDim,nOrder,GeoShape::nRealDim>,value_type, TheTag > > >::type::type type;
232 
233  typedef boost::shared_ptr<type> ptrtype;
234  typedef boost::shared_ptr<const type> const_ptrtype;
235  };
236  typedef typename trace_trace_mesh<Tag>::type trace_trace_mesh_type;
237  typedef typename trace_trace_mesh<Tag>::ptrtype trace_trace_mesh_ptrtype;
238 
239 
241 
245  Mesh( WorldComm const& worldComm = Environment::worldComm() );
246 
247  ~Mesh()
248  {
249  M_gm.reset();
250  M_gm1.reset();
251  M_tool_localization.reset();
252  }
257  static mesh_ptrtype New()
258  {
259  return mesh_ptrtype(new mesh_type);
260  }
261 
262  self_type& operator+=( self_type const& m );
263 
267 
271 #if 0
272  size_type numElements() const
273  {
274 
275  return std::distance( this->beginElementWithProcessId( this->worldComm().rank() ),
276  this->endElementWithProcessId( this->worldComm().rank() ) );
277  }
278 #endif
279 
280 
282  {
283  //int ne = numElements();
284  int ne = std::distance( this->beginElementWithProcessId( this->worldComm().rank() ),
285  this->endElementWithProcessId( this->worldComm().rank() ) );
286  int gne;
287  mpi::all_reduce( this->worldComm(), ne, gne, [] ( int x, int y )
288  {
289  return x + y;
290  } );
291  return gne;
292  }
296  uint16_type dimension() const
297  {
298  return nDim;
299  }
300 
304  gm_ptrtype const& gm() const
305  {
306  return M_gm;
307  }
308 
312  gm1_ptrtype const& gm1() const
313  {
314  return M_gm1;
315  }
316 
320  gm_ptrtype& gm()
321  {
322  return M_gm;
323  }
324 
328  gm1_ptrtype& gm1()
329  {
330  return M_gm1;
331  }
332 
337  reference_convex_type referenceConvex() const
338  {
339  return reference_convex_type();
340  }
341 
345  face_processor_type const& localFaceId( element_type const& e,
346  size_type const n ) const
347  {
348  return M_e2f.find(std::make_pair(e.id(),n))->second;
349  }
350 
354  face_processor_type const& localFaceId( size_type const e,
355  size_type const n ) const
356  {
357  return M_e2f.find(std::make_pair(e,n))->second;
358  }
359 #if 0
360 
363  WorldComm const& worldComm() const
364  {
365  return M_worldComm;
366  }
367 
368  void setWorldComm( WorldComm const& _worldComm )
369  {
370  M_worldComm = _worldComm;
371  }
372 
373  mpi::communicator const& comm() const
374  {
375  return M_worldComm.localComm();
376  }
377 #endif
378 
379 
383 
387  void setPartitioner( std::string partitioner )
388  {
389  //M_part = partitioner_ptrtype( partitioner_type::New( partitioner ) );
390  }
391 
395  face_processor_type& localFaceId( element_type const& e,
396  size_type const n )
397  {
398  return M_e2f[std::make_pair(e.id(),n)];
399  }
400 
404  face_processor_type& localFaceId( size_type const e,
405  size_type const n )
406  {
407  return M_e2f[std::make_pair(e,n)];
408  }
409 
413  size_type markerName( std::string const& marker ) const
414  {
415  auto mit = M_markername.find( marker );
416  if ( mit != M_markername.end() )
417  return mit->second[0];
419  }
423  size_type markerDim( std::string const& marker ) const
424  {
425  auto mit = M_markername.find( marker );
426  if ( mit != M_markername.end() )
427  return mit->second[1];
429  }
430 
434  std::map<std::string, std::vector<size_type> > markerNames() const
435  {
436  return M_markername;
437  }
438 
442  struct Localization;
443  boost::shared_ptr<typename self_type::Localization> tool_localization()
444  {
445  return M_tool_localization;
446  }
447 
451  value_type measure() const
452  {
453  return M_meas;
454  }
455 
459  value_type measureBoundary() const
460  {
461  return M_measbdy;
462  }
463 
465 
469 
473  void addMarkerName( std::pair<std::string, std::vector<size_type> > const& marker )
474  {
475  M_markername.insert( marker );
476  }
477 
481  void addMarkerName( std::string __name, int __id ,int __topoDim )
482  {
483  std::vector<size_type> data(2);
484  data[0]=__id;
485  data[1]=__topoDim;
486  M_markername[__name]=data;
487  }
488 
489  flag_type markerId( boost::any const& marker );
490 
501  element_iterator eraseElement( element_iterator position, bool modify=true );
502 
503 
504  template<typename RangeT, int TheTag>
505  typename trace_mesh<TheTag>::ptrtype
506  trace( RangeT const& range, mpl::int_<TheTag> ) const;
507 
520  template<int TheTag=Tag>
521  typename trace_mesh<TheTag>::ptrtype
522  trace() const
523  {
524  return trace(boundaryfaces(this->shared_from_this()),mpl::int_<TheTag>());
525  }
526 
527  template<typename RangeT>
528  typename trace_mesh<Tag>::ptrtype
529  trace( RangeT const& range ) const;
530 
531 
532  template<int TheTag=Tag>
533  typename trace_mesh<TheTag>::ptrtype
534  wireBasket() const
535  {
536  return wireBasket(boundaryfaces(this->shared_from_this()),mpl::int_<TheTag>());
537  }
538 
539  template<typename RangeT, int TheTag>
540  typename trace_trace_mesh<TheTag>::ptrtype
541  wireBasket( RangeT const& range, mpl::int_<TheTag> ) const;
542 
543  template<typename RangeT>
544  typename trace_trace_mesh<Tag>::ptrtype
545  wireBasket( RangeT const& range ) const;
546 
547 
548  template<typename Iterator>
549  void createSubmesh( self_type& mesh,
550  Iterator const& begin_elt,
551  Iterator const& end_elt,
552  size_type extraction_policies = EXTRACTION_KEEP_ALL_IDS ) const;
553 
561  void createSubmeshByProcessId( self_type& mesh, uint16_type pid ) const
562  {
563  this->createSubmesh( mesh,
564  this->beginElementWithProcessId( pid ),
565  this->endElementWithProcessId( pid ) );
566  }
567 
571  P1_mesh_ptrtype createP1mesh() const;
572 
578  template<typename IteratorRange>
579  void updateMarker2WithRange( IteratorRange const& range,flag_type flag )
580  {
581  const size_type iDim = boost::tuples::template element<0, IteratorRange>::type::value;
582  this->updateMarker2WithRange( range,flag,mpl::int_<iDim>() );
583  }
584 
588  template<typename IteratorRange>
589  void updateMarker2WithRange( IteratorRange const& range,flag_type flag, mpl::int_<MESH_ELEMENTS> )
590  {
591  const size_type iDim = boost::tuples::template element<0, IteratorRange>::type::value;
592  this->updateMarker2WithRangeElements( range,flag );
593  }
594 
598  template<typename IteratorRange>
599  void updateMarker2WithRange( IteratorRange const& range,flag_type flag, mpl::int_<MESH_FACES> )
600  {
601  this->updateMarker2WithRangeFaces( range,flag );
602  }
603 
609  template<typename IteratorRange>
610  void updateMarker3WithRange( IteratorRange const& range,flag_type flag )
611  {
612  const size_type iDim = boost::tuples::template element<0, IteratorRange>::type::value;
613  this->updateMarker3WithRange( range,flag,mpl::int_<iDim>() );
614  }
615 
619  template<typename IteratorRange>
620  void updateMarker3WithRange( IteratorRange const& range,flag_type flag, mpl::int_<MESH_ELEMENTS> )
621  {
622  this->updateMarker3WithRangeElements( range,flag );
623  }
624 
628  template<typename IteratorRange>
629  void updateMarker3WithRange( IteratorRange const& range,flag_type flag, mpl::int_<MESH_FACES> )
630  {
631  this->updateMarker3WithRangeFaces( range,flag );
632  }
633 
637  void partition ( const uint16_type n_parts = 1 );
638 
647  void renumber()
648  {
649  renumber( mpl::bool_<( nDim > 1 )>() );
650  }
651  void renumber( std::vector<size_type> const& node_map, mpl::int_<1> );
652  void renumber( std::vector<size_type> const& node_map, mpl::int_<2> );
653  void renumber( std::vector<size_type> const& node_map, mpl::int_<3> );
654 
662  void localrenumber();
663 
667  void checkAndFixPermutation();
668 
672  void send( int p, int tag );
673 
677  void recv( int p, int tag );
678 
684  void encode();
685 
691  void decode();
692 
693  BOOST_PARAMETER_MEMBER_FUNCTION( ( void ),
694  save,
695  tag,
696  ( required
697  ( name,(std::string) )
698  ( path,* ) )
699  ( optional
700  ( type,( std::string ),std::string( "binary" ) )
701  ( suffix,( std::string ),std::string( "" ) )
702  ( sep,( std::string ),std::string( "" ) )
703  ) )
704  {
705  Feel::detail::ignore_unused_variable_warning( args );
706 
707  if ( !fs::exists( fs::path( path ) ) )
708  {
709  fs::create_directories( fs::path( path ) );
710  }
711 
712  std::ostringstream os1;
713  os1 << name << sep << suffix << "-" << this->worldComm().globalSize() << "." << this->worldComm().globalRank() << ".fdb";
714  fs::path p = fs::path( path ) / os1.str();
715  fs::ofstream ofs( p );
716 
717  if ( type == "binary" )
718  {
719  boost::archive::binary_oarchive oa( ofs );
720  oa << *this;
721  }
722 
723  else if ( type == "text" )
724  {
725  boost::archive::text_oarchive oa( ofs );
726  oa << *this;
727  }
728 
729  else if ( type == "xml" )
730  {
731  //boost::archive::xml_oarchive oa(ofs);
732  //oa << *this;
733  }
734  }
735  BOOST_PARAMETER_MEMBER_FUNCTION(
736  ( bool ),
737  load,
738  tag,
739  ( required
740  ( name,(std::string) )
741  ( path,* ) )
742  ( optional
743  ( update,(size_type), MESH_CHECK|MESH_UPDATE_EDGES|MESH_UPDATE_FACES )
744  ( type,( std::string ),std::string( "binary" ) )
745  ( suffix,( std::string ),std::string( "" ) )
746  ( sep,( std::string ),std::string( "" ) )
747  )
748  )
749  {
750  Feel::detail::ignore_unused_variable_warning( args );
751  std::ostringstream os1;
752  os1 << name << sep << suffix << "-" << this->worldComm().globalSize() << "." << this->worldComm().globalRank() << ".fdb";
753  fs::path p = fs::path( path ) / os1.str();
754  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
755  std::cout << "try loading " << p.native() << "\n";
756  if ( !fs::exists( p ) )
757  {
758  LOG(INFO) << "[mesh::load] failed loading " << p.native() << "\n";
759  std::ostringstream os2;
760  os2 << name << sep << suffix << "-" << this->worldComm().globalSize() << "." << this->worldComm().globalRank();
761  p = fs::path( path ) / os2.str();
762  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
763  std::cout << " now try loading " << p.native() << "\n";
764 
765  if ( !fs::exists( p ) )
766  {
767  LOG(INFO) << "[mesh::load] failed loading " << p.native() << "\n";
768  return false;
769  }
770  }
771 
772  if ( !fs::is_regular_file( p ) )
773  {
774  LOG(INFO) << "[mesh::load] failed loading " << p.native() << "\n";
775  return false;
776  }
777 
778  fs::ifstream ifs( p );
779 
780  if ( type == "binary" )
781  {
782  boost::archive::binary_iarchive ia( ifs );
783  ia >> *this;
784  }
785 
786  else if ( type == "text" )
787  {
788  boost::archive::text_iarchive ia( ifs );
789  ia >> *this;
790  }
791 
792  else if ( type == "xml" )
793  {
794  //boost::archive::xml_iarchive ia(ifs);
795  //ia >> *this;
796  }
797  if ( update )
798  {
799  this->components().reset();
800  this->components().set( update );
801  this->updateForUse();
802  }
803 
804  else
805  {
806  this->components().reset();
807  }
808  return true;
809  }
810 
811 
814 
815  //private:
816 
826 
832  void checkLocalPermutation( mpl::bool_<false> ) const {}
833  void checkLocalPermutation( mpl::bool_<true> ) const;
834 
835 
843  void updateForUse();
844 
845  void meshModified()
846  {
847  for( auto fit = this->beginFace(), fen = this->endFace(); fit != fen; ++fit )
848  {
849  if ( fit->isOnBoundary() && !fit->isConnectedTo0() )
850  {
851  std::cout << "erase boundary face...\n";
852  this->eraseFace( fit );
853  }
854  if ( !fit->isOnBoundary() && fit->isConnectedTo0() && !fit->isConnectedTo1() )
855  {
856  std::cout << "found boundary face...\n";
857  this->faces().modify( fit, []( face_type& f ){ f.setOnBoundary( true ); } );
858  }
859  }
860  this->setUpdatedForUse( false );
861  this->updateForUse();
862  }
863 private:
864 
865  void propagateMarkers( mpl::int_<1> ) {}
866  void propagateMarkers( mpl::int_<2> );
867  void propagateMarkers( mpl::int_<3> );
868 
869  friend class boost::serialization::access;
870  template<class Archive>
871  void serialize( Archive & ar, const unsigned int version )
872  {
873  if ( Archive::is_saving::value )
874  {
875  DVLOG(2) << "Serializing mesh(saving) ...\n";
876  DVLOG(2) << "encoding...\n";
877  encode();
878  DVLOG(2) << "loading markers...\n";
879  ar & M_markername;
880  DVLOG(2) << "loading pts...\n";
881  ar & M_enc_pts;
882  DVLOG(2) << "loading faces...\n";
883  ar & M_enc_faces;
884  DVLOG(2) << "loading elts...\n";
885  ar & M_enc_elts;
886  }
887  if ( Archive::is_loading::value )
888  {
889  DVLOG(2) << "Serializing mesh(loading) ...\n";
890  DVLOG(2) << "loading markers...\n";
891  ar & M_markername;
892  DVLOG(2) << "loading pts...\n";
893  ar & M_enc_pts;
894  DVLOG(2) << "loading faces...\n";
895  ar & M_enc_faces;
896  DVLOG(2) << "loading elts...\n";
897  ar & M_enc_elts;
898  decode();
899 
900  }
901 
902  }
903 public:
904  struct Inverse
905  :
906  public mpl::if_<mpl::bool_<GeoShape::is_simplex>,
907  mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Simplex> >,
908  mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Hypercube> > >::type::type
909  {
910  typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
911  mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Simplex> >,
912  mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Hypercube> > >::type::type super;
913  typedef typename super::gic_type gic_type;
914 
915  typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
916  mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Simplex> >,
917  mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Hypercube> > >::type::type super1;
918  typedef typename super1::gic_type gic1_type;
919 
920  Inverse( boost::shared_ptr<self_type> const& m )
921  :
922  super(),
923  M_mesh ( m )
924  {}
925 
926  ~Inverse()
927  {}
928 
929  size_type nPointsInConvex( size_type i ) const
930  {
931  return M_pts_cvx[i].size();
932  }
933  void pointsInConvex( size_type i, std::vector<boost::tuple<size_type, uint16_type > > &itab ) const
934  {
935  itab.resize( M_pts_cvx[i].size() );
936  size_type j = 0;
937 
938  for ( map_iterator it = M_pts_cvx[i].begin(); it != M_pts_cvx[i].end(); ++it )
939  itab[j++] = boost::make_tuple( it->first, it->second );
940  }
941 
942  const boost::unordered_map<size_type,node_type> &referenceCoords( void )
943  {
944  return M_ref_coords;
945  }
946 
955  void distribute( bool extrapolation = false );
956 
957 
958 
959  private:
960  boost::shared_ptr<self_type> M_mesh;
961  std::vector<std::map<size_type,uint16_type > > M_pts_cvx;
962  typedef typename std::map<size_type, uint16_type >::const_iterator map_iterator;
963  //typedef typename node<value_type>::type node_type;
964  boost::unordered_map<size_type,node_type> M_ref_coords;
965  boost::unordered_map<size_type,double> M_dist;
966  boost::unordered_map<size_type,size_type> M_cvx_pts;
967 
968  }; // Inverse
969 
970  struct Localization
971  {
972 
973  typedef Localization localization_type;
974  typedef boost::shared_ptr<localization_type> localization_ptrtype;
975 
976  typedef typename matrix_node<typename node_type::value_type>::type matrix_node_type;
977 
978  typedef KDTree kdtree_type;
979  typedef typename boost::shared_ptr<KDTree> kdtree_ptrtype;
980 
981  // a node x => a list of id elt which contain the node x
982  typedef boost::tuple<node_type, std::list<size_type> > node_elem_type;
983 
984  typedef typename std::list<boost::tuple<size_type,node_type> > container_output_type;
985  typedef typename std::list<boost::tuple<size_type,node_type> >::iterator container_output_iterator_type;
986  //map between element id and list of node described in the reference elt
987  //typedef std::map<size_type, std::list<boost::tuple<size_type,node_type> > > container_search_type
988  typedef std::map<size_type, container_output_type > container_search_type;
989  typedef typename container_search_type::const_iterator container_search_const_iterator_type;
990  typedef typename container_search_type::iterator container_search_iterator_type;
991 
992  // geomap inverse
993  typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
994  mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Simplex> >,
995  mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Hypercube> > >::type::type gm_inverse_type;
996  typedef typename gm_inverse_type::gic_type gmc_inverse_type;
997 
998  typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
999  mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Simplex> >,
1000  mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Hypercube> > >::type::type gm1_inverse_type;
1001  typedef typename gm1_inverse_type::gic_type gmc1_inverse_type;
1002 
1003  // reference convex
1004  typedef typename self_type::gm_type::reference_convex_type ref_convex_type;
1005  typedef typename self_type::gm1_type::reference_convex_type ref_convex1_type;
1006 
1010  Localization() :
1011  M_mesh (),
1012  M_kd_tree( new kdtree_type() ),
1013  IsInit( false ),
1014  IsInitBoundaryFaces( false ),
1015  M_doExtrapolation( true )
1016  {
1017  DVLOG(2) << "[Mesh::Localization] create Localization tool\n";
1018  M_kd_tree->nbNearNeighbor( 15 );
1019  M_resultAnalysis.clear();
1020  DVLOG(2) << "[Mesh::Localization] create Localization tool done\n";
1021  }
1022 
1023  Localization( boost::shared_ptr<self_type> m, bool init_b = true ) :
1024  M_mesh ( m ),
1025  IsInit( init_b ),
1026  IsInitBoundaryFaces( false ),
1027  M_doExtrapolation( true )
1028  {
1029  if ( IsInit )
1030  init();
1031 
1032  M_kd_tree->nbNearNeighbor( 15 );
1033  M_resultAnalysis.clear();
1034  }
1035 
1036  Localization( Localization const & L ) :
1037  M_mesh( L.M_mesh ),
1038  M_kd_tree( new kdtree_type( *( L.M_kd_tree ) ) ),
1039  M_geoGlob_Elts( L.M_geoGlob_Elts ),
1040  IsInit( L.IsInit ),
1041  IsInitBoundaryFaces( L.IsInitBoundaryFaces ),
1042  M_resultAnalysis( L.M_resultAnalysis ),
1043  M_doExtrapolation( L.M_doExtrapolation )
1044  {}
1045 
1046  /*--------------------------------------------------------------
1047  * Define the mesh whith or not init
1048  */
1049  void
1050  setMesh( boost::shared_ptr<self_type> m,bool b=true )
1051  {
1052  M_mesh=m;
1053 
1054  if ( b )
1055  init();
1056 
1057  else IsInit=b;
1058 
1059  M_resultAnalysis.clear();
1060  }
1061 
1062  /*--------------------------------------------------------------
1063  * Define if necessary to use extrapolation
1064  */
1065  void
1066  setExtrapolation( bool b )
1067  {
1068  M_doExtrapolation = b;
1069  }
1070 
1071  /*--------------------------------------------------------------
1072  * Run the init function if necessary
1073  */
1074  void updateForUse()
1075  {
1076  if ( IsInit==false )
1077  init();
1078  }
1079 
1080  /*--------------------------------------------------------------
1081  * Run the init function if necessary
1082  */
1083  void updateForUseBoundaryFaces()
1084  {
1085  if ( IsInitBoundaryFaces==false )
1086  initBoundaryFaces();
1087  }
1088 
1089  /*--------------------------------------------------------------
1090  * Access
1091  */
1092  bool isInit() const
1093  {
1094  return IsInit;
1095  }
1096 
1097  bool isInitBoundaryFaces() const
1098  {
1099  return IsInitBoundaryFaces;
1100  }
1101 
1102  bool doExtrapolation() const
1103  {
1104  return M_doExtrapolation;
1105  }
1106 
1107  boost::shared_ptr<self_type> mesh()
1108  {
1109  return M_mesh;
1110  }
1111  boost::shared_ptr<self_type> const& mesh() const
1112  {
1113  return M_mesh;
1114  }
1115 
1116  kdtree_ptrtype kdtree()
1117  {
1118  return M_kd_tree;
1119  }
1120 
1121  kdtree_ptrtype const& kdtree() const
1122  {
1123  return M_kd_tree;
1124  }
1125 
1126  node_type barycenter() const;
1127  node_type barycenter(mpl::int_<1> ) const;
1128  node_type barycenter(mpl::int_<2> ) const;
1129  node_type barycenter(mpl::int_<3> ) const;
1130 
1131  container_search_type const & result_analysis() const { return M_resultAnalysis;}
1132 
1133  container_search_iterator_type result_analysis_begin()
1134  {
1135  return M_resultAnalysis.begin();
1136  }
1137  container_search_iterator_type result_analysis_end()
1138  {
1139  return M_resultAnalysis.end();
1140  }
1141 
1142  /*---------------------------------------------------------------
1143  * True if the node p is in mesh->element(id)
1144  */
1145  boost::tuple<bool,node_type,double> isIn( size_type _id, const node_type & _pt ) const;
1146  boost::tuple<bool,node_type,double> isIn( size_type _id, const node_type & _pt, const matrix_node_type & setPoints, mpl::int_<1> ) const;
1147  boost::tuple<uint16_type,std::vector<bool> > isIn( std::vector<size_type> _ids, const node_type & _pt );
1148 
1149  /*---------------------------------------------------------------
1150  * Research only one element wich contains the node p
1151  */
1152  boost::tuple<bool, size_type,node_type> searchElement(const node_type & p);
1153  /*---------------------------------------------------------------
1154  * Research only one element wich contains the node p
1155  */
1156  boost::tuple<bool, size_type,node_type> searchElement(const node_type & p,
1157  const matrix_node_type & setPoints,
1158  mpl::int_<0> )
1159  {
1160  return searchElement( p );
1161  }
1162 
1163  /*---------------------------------------------------------------
1164  * Research only one element wich contains the node p and which this elt have as geometric point contain setPoints
1165  */
1166  boost::tuple<bool, size_type,node_type> searchElement( const node_type & p,
1167  const matrix_node_type & setPoints,
1168  mpl::int_<1> );
1169 
1170  /*---------------------------------------------------------------
1171  * Research all elements wich contains the node p
1172  */
1173  boost::tuple<bool, std::list<boost::tuple<size_type,node_type> > > searchElements( const node_type & p );
1174 
1175  /*---------------------------------------------------------------
1176  * Research the element wich contains the node p, forall p in the
1177  * matrix_node_type m. The result is save by this object
1178  */
1179  boost::tuple<std::vector<bool>, size_type> run_analysis(const matrix_node_type & m,
1180  const size_type & eltHypothetical);
1181 
1182  /*---------------------------------------------------------------
1183  * Research the element wich contains the node p, forall p in the
1184  * matrix_node_type m. The result is save by this object
1185  */
1186  boost::tuple<std::vector<bool>, size_type> run_analysis(const matrix_node_type & m,
1187  const size_type & eltHypothetical,
1188  const matrix_node_type & setPoints,
1189  mpl::int_<0> )
1190  {
1191  return run_analysis( m,eltHypothetical );
1192  }
1193 
1194  /*---------------------------------------------------------------
1195  * Research the element wich contains the node p, forall p in the
1196  * matrix_node_type m. The result is save by this object
1197  */
1198  boost::tuple<std::vector<bool>, size_type> run_analysis(const matrix_node_type & m,
1199  const size_type & eltHypothetical,
1200  const matrix_node_type & setPoints,
1201  mpl::int_<1> );
1202 
1203  /*---------------------------------------------------------------
1204  * Reset all data
1205  */
1206  void reset()
1207  {
1208  IsInit=false;
1209  init();
1210  }
1211 
1212  void resetBoundaryFaces()
1213  {
1214  IsInitBoundaryFaces=false;
1215  initBoundaryFaces();
1216  }
1217 
1218  private :
1219 
1220  /*---------------------------------------------------------------
1221  *initializes the kd tree and the map between node and list elements(all elements)
1222  */
1223  void init();
1224 
1225  /*---------------------------------------------------------------
1226  *initializes the kd tree and the map between node and list elements(only on boundary)
1227  */
1228  void initBoundaryFaces();
1229 
1230  /*---------------------------------------------------------------
1231  *search near elt in kd tree and get a sorted list
1232  */
1233  void searchInKdTree( const node_type & p,
1234  std::list< std::pair<size_type, uint> > & listTri );
1235 
1236  private:
1237 
1238  boost::shared_ptr<self_type> M_mesh;
1239  //KDTree M_kd_tree;
1240  kdtree_ptrtype M_kd_tree;
1241  //map between node and list elements
1242  std::map<size_type, node_elem_type > M_geoGlob_Elts;
1243  bool IsInit,IsInitBoundaryFaces;
1244  container_search_type M_resultAnalysis;
1245  bool M_doExtrapolation;
1246 
1247  ref_convex_type M_refelem;
1248  ref_convex1_type M_refelem1;
1249  };
1250 
1251 
1255 
1256 
1260  boost::signals2::signal<void ( MESH_CHANGES )> meshChanged;
1261 
1262  template<typename Observer>
1263  void addObserver( Observer& obs )
1264  {
1265  meshChanged.connect( obs );
1266  }
1267 
1268  void removeFacesFromBoundary( std::initializer_list<uint16_type> markers );
1269 
1270 
1271 
1273 
1274 protected:
1279  void updateEntitiesCoDimensionOne(mpl::bool_<true>);
1280  void updateEntitiesCoDimensionOne(mpl::bool_<false>);
1281 
1286 
1290  void check() const;
1291 
1292 
1293 private:
1294 
1298  void renumber( mpl::bool_<false> ) {}
1299 
1303  void renumber( mpl::bool_<true> );
1304 
1308  void modifyEdgesOnBoundary( face_iterator& face, mpl::bool_<true> );
1309 
1313  void modifyEdgesOnBoundary( face_iterator& face, mpl::bool_<false> );
1314 
1318  bool modifyElementOnBoundaryFromEdge( element_iterator& e, mpl::bool_<false> );
1319 
1323  bool modifyElementOnBoundaryFromEdge( element_iterator& e, mpl::bool_<true> );
1324 
1328  void updateOnBoundary();
1329 
1330 private:
1331 
1333  //WorldComm M_worldComm;
1334 
1335  gm_ptrtype M_gm;
1336  gm1_ptrtype M_gm1;
1337 
1339  value_type M_meas;
1340 
1342  value_type M_measbdy;
1343 
1348  std::vector<uint16_type> M_neighboring_processors;
1349 
1350  //partitioner_ptrtype M_part;
1351 
1355  boost::unordered_map<std::pair<int,int>,face_processor_type> M_e2f;
1356 
1360  boost::multi_array<element_edge_type,2> M_e2e;
1361 
1367  std::map<std::string, std::vector<size_type> > M_markername;
1368 
1372  std::map<uint64_type, boost::tuple<bool, std::vector<int>, std::vector<double> > > M_enc_pts;
1373 
1377  std::map<uint64_type, std::vector<int> > M_enc_faces;
1378 
1382  std::map<uint64_type, std::vector<int> > M_enc_elts;
1383 
1387  boost::shared_ptr<Localization> M_tool_localization;
1388 
1389 };
1390 
1391 template<typename Shape, typename T, int Tag>
1392 const uint16_type Mesh<Shape, T, Tag>::nDim;
1393 template<typename Shape, typename T, int Tag>
1394 const uint16_type Mesh<Shape, T, Tag>::nOrder;
1395 
1396 template<typename Shape, typename T, int Tag>
1397 template<typename RangeT>
1398 typename Mesh<Shape, T, Tag>::template trace_mesh<Tag>::ptrtype
1399 Mesh<Shape, T, Tag>::trace( RangeT const& range ) const
1400 {
1401  DVLOG(2) << "[trace] extracting " << range.template get<0>() << " nb elements :"
1402  << std::distance(range.template get<1>(),range.template get<2>()) << "\n";
1403  return Feel::createSubmesh<const mesh_type,RangeT,Tag>( this->shared_from_this(), range );
1404 
1405 }
1406 
1407 template<typename Shape, typename T, int Tag>
1408 template<typename RangeT>
1409 typename Mesh<Shape, T, Tag>::template trace_trace_mesh<Tag>::ptrtype
1410 Mesh<Shape, T, Tag>::wireBasket( RangeT const& range ) const
1411 {
1412  DVLOG(2) << "[trace] extracting " << range.template get<0>() << " nb elements :"
1413  << std::distance(range.template get<1>(),range.template get<2>()) << "\n";
1414  return Feel::createSubmesh<const mesh_type,RangeT,Tag>( this->shared_from_this(), range );
1415 
1416 }
1417 
1418 
1419 template<int TheTag> struct Tag : public mpl::int_<TheTag> {};
1420 
1421 template<typename Shape, typename T, int Tag>
1422 template<typename RangeT,int TheTag>
1423 typename Mesh<Shape, T, Tag>::template trace_mesh<TheTag>::ptrtype
1424 Mesh<Shape, T, Tag>::trace( RangeT const& range, mpl::int_<TheTag> ) const
1425 {
1426  DVLOG(2) << "[trace] extracting " << range.template get<0>() << " nb elements :"
1427  << std::distance(range.template get<1>(),range.template get<2>()) << "\n";
1428  return Feel::createSubmesh<const mesh_type,RangeT,TheTag>( this->shared_from_this(), range );
1429 
1430 }
1431 
1432 template<typename Shape, typename T, int Tag>
1433 template<typename RangeT,int TheTag>
1434 typename Mesh<Shape, T, Tag>::template trace_trace_mesh<TheTag>::ptrtype
1435 Mesh<Shape, T, Tag>::wireBasket( RangeT const& range, mpl::int_<TheTag> ) const
1436 {
1437  DVLOG(2) << "[trace] extracting " << range.template get<0>() << " nb elements :"
1438  << std::distance(range.template get<1>(),range.template get<2>()) << "\n";
1439  return Feel::createSubmesh<const mesh_type,RangeT,TheTag>( this->shared_from_this(), range );
1440 
1441 }
1442 
1443 template<typename Shape, typename T, int Tag>
1444 template<typename Iterator>
1445 void
1446 Mesh<Shape, T, Tag>::createSubmesh( self_type& new_mesh,
1447  Iterator const& begin_elt,
1448  Iterator const& end_elt,
1449  size_type extraction_policies ) const
1450 {
1451 #if 0
1452  new_mesh = Feel::createSubmesh( this->shared_from_this(), boost::make_tuple(mpl::int_<MESH_ELEMENTS>(),begin_elt, end_elt), extraction_policies );
1453 #else
1454  Context policies( extraction_policies );
1455 
1456  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] start\n";
1457  // make sure it is all empty
1458  new_mesh.clear();
1459  // inherit all the markers
1460  new_mesh.M_markername = this->markerNames();
1461  BOOST_FOREACH( auto marker, new_mesh.M_markername )
1462  {
1463  LOG(INFO) << "marker name " << marker.first
1464  << " id: " << marker.second[0]
1465  << " geoe: " << marker.second[1] << "\n";
1466 
1467  }
1468  // How the nodes on this mesh will be renumbered to nodes
1469  // on the new_mesh.
1470  std::vector<size_type> new_node_numbers ( this->numPoints() );
1471  std::vector<size_type> new_vertex ( this->numPoints() );
1472 
1473  std::fill ( new_node_numbers.begin(),
1474  new_node_numbers.end(),
1476 
1477  std::fill ( new_vertex.begin(),
1478  new_vertex.end(),
1479  0 );
1480 
1481 
1482 
1483  // the number of nodes on the new mesh, will be incremented
1484  unsigned int n_new_nodes = 0;
1485  unsigned int n_new_elem = 0;
1486  size_type n_new_faces = 0;
1487 
1488  for ( Iterator it = begin_elt; it != end_elt; ++it )
1489  {
1490 
1491  element_type const& old_elem = *it;
1492 
1493  // copy element so that we can modify it
1494  element_type new_elem = old_elem;
1495 
1496  // Loop over the nodes on this element.
1497  for ( unsigned int n=0; n < old_elem.nPoints(); n++ )
1498  {
1499  FEELPP_ASSERT ( old_elem.point( n ).id() < new_node_numbers.size() ).error( "invalid point id()" );
1500 
1501  if ( new_node_numbers[old_elem.point( n ).id()] == invalid_size_type_value )
1502  {
1503  new_node_numbers[old_elem.point( n ).id()] = n_new_nodes;
1504 
1505  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] insert point " << old_elem.point( n ) << "\n";
1506 
1507  point_type pt( old_elem.point( n ) );
1508  pt.setId( n_new_nodes );
1509 
1510  // Add this node to the new mesh
1511  new_mesh.addPoint ( pt );
1512 
1513  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] number of points " << new_mesh.numPoints() << "\n";
1514 
1515  // Increment the new node counter
1516  n_new_nodes++;
1517 
1518  if ( n < element_type::numVertices )
1519  {
1520  FEELPP_ASSERT( new_vertex[old_elem.point( n ).id()] == 0 ).error( "already seen this point?" );
1521  new_vertex[old_elem.point( n ).id()]=1;
1522  }
1523  }
1524 
1525  // Define this element's connectivity on the new mesh
1526  FEELPP_ASSERT ( new_node_numbers[old_elem.point( n ).id()] < new_mesh.numPoints() ).error( "invalid connectivity" );
1527 
1528  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] adding point old(" << old_elem.point( n ).id()
1529  << ") as point new(" << new_node_numbers[old_elem.point( n ).id()]
1530  << ") in element " << new_elem.id() << "\n";
1531 
1532  new_elem.setPoint( n, new_mesh.point( new_node_numbers[old_elem.point( n ).id()] ) );
1533 
1534  }
1535 
1536  // set id of element
1537  new_elem.setId ( n_new_elem );
1538 
1539  // increment the new element counter
1540  n_new_elem++;
1541 
1542 
1543  // Add an equivalent element type to the new_mesh
1544  new_mesh.addElement( new_elem );
1545 
1546 
1547  // Maybe add faces for this element
1548  for ( unsigned int s=0; s<old_elem.numTopologicalFaces; s++ )
1549  {
1550  if ( !old_elem.facePtr( s ) ) continue;
1551 
1552 #if 0
1553  std::cout << "local face id: " << s
1554  << " global face id: " << old_elem.face( s ).id() << "\n";
1555 #endif
1556  // only add face on the boundary: they have some data
1557  // (boundary ids) which cannot be retrieved otherwise
1558  //if ( old_elem.neighbor(s) == invalid_size_type_value )
1559  size_type global_face_id = old_elem.face( s ).id();
1560 
1561  if ( this->hasFace( global_face_id ) )
1562  {
1563 
1564  //std::cout << "found face " << global_face_id << "\n";
1565  // get the corresponding face
1566  face_type const& old_face = old_elem.face( s );
1567  face_type new_face = old_face;
1568 
1569  // disconnect from elements of old mesh,
1570  // the connection will be redone in
1571  // \c updateForUse()
1572  new_face.disconnect();
1573 
1574  //std::cout << "disconnect face\n";
1575  // update points info
1576  for ( uint16_type p = 0; p < new_face.nPoints(); ++p )
1577  {
1578 #if 0
1579  std::cout << "add new point " << new_face.point( p ).id() << " to face \n";
1580  std::cout << "add old point " << old_face.point( p ).id() << " to face \n";
1581  std::cout << "new point id " << new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).id()] << "\n";
1582 #endif
1583  new_face.setPoint( p, new_mesh.point( new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).id()] ) );
1584 
1585  }
1586 
1587  new_face.setId( n_new_faces++ );
1588 #if 0
1589  std::cout << "face id" << new_face.id()
1590  << " marker1 : " << new_face.marker()
1591  << " old marker1 : " << old_face.marker()
1592  << "\n";
1593 #endif
1594  // add it to the list of faces
1595  new_mesh.addFace( new_face );
1596  }
1597  }
1598  }
1599 
1600  new_mesh.setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0 ) );
1601 
1602  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] update face/edge info if necessary\n";
1603  // Prepare the new_mesh for use
1604  new_mesh.components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
1605  new_mesh.updateForUse();
1606 
1607  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] stop\n";
1608 #endif
1609 }
1610 
1611 
1612 template<typename Shape, typename T, int Tag>
1613 typename Mesh<Shape, T, Tag>::P1_mesh_ptrtype
1615 {
1616 
1617  P1_mesh_ptrtype new_mesh( new P1_mesh_type( this->worldComm() ) );
1618 
1619  // How the nodes on this mesh will be renumbered to nodes
1620  // on the new_mesh.
1621  boost::unordered_map<size_type,size_type> new_node_numbers;// ( this->numPoints() );
1622  boost::unordered_map<size_type,int> new_vertex;// ( this->numPoints() );
1623  std::vector<size_type> new_element_numbers ( this->numElements() );
1624 
1625  // the number of nodes on the new mesh, will be incremented
1626  unsigned int n_new_nodes = 0;
1627  unsigned int n_new_elem = 0;
1628  size_type n_new_faces = 0;
1629 
1630  // inherit the table of markersName
1631  BOOST_FOREACH( auto itMark, this->markerNames() )
1632  {
1633  new_mesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
1634  }
1635 
1636 
1637 
1638  const int proc_id = this->worldComm().localRank();
1639  const int nProc = this->worldComm().localSize();
1640  std::vector< std::list<boost::tuple<size_type,size_type> > > memory_ghostid( nProc );
1641  std::vector< std::vector<size_type> > memory_id( nProc );
1642  std::vector< std::vector<size_type> > vecToSend( nProc );
1643  std::vector< std::vector<size_type> > vecToRecv( nProc );
1644 
1645 
1646  //auto it = this->beginElementWithProcessId( this->worldComm().localRank() );
1647  //auto en = this->endElementWithProcessId( this->worldComm().localRank() );
1648  auto it = this->beginElement();
1649  auto en = this->endElement();
1650 
1651  for ( ; it != en; ++it )
1652  {
1653  element_type const& old_elem = *it;
1654 
1655  // create a new element
1656  typename P1_mesh_type::element_type new_elem;
1657  // set id of element
1658  new_elem.setId ( n_new_elem );
1659  new_element_numbers[old_elem.id()]= n_new_elem;
1660  // set element markers
1661  new_elem.setMarker( old_elem.marker().value() );
1662  new_elem.setMarker2( old_elem.marker2().value() );
1663  new_elem.setMarker3( old_elem.marker3().value() );
1664  // partitioning update
1665  new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
1666  new_elem.setNumberOfPartitions(old_elem.numberOfPartitions());
1667  new_elem.setProcessId(old_elem.processId());
1668  //new_elem.setIdInPartition( old_elem.pidInPartition(), n_new_elem );
1669  new_elem.setNeighborPartitionIds(old_elem.neighborPartitionIds());
1670 
1671  // Loop over the P1 nodes on this element.
1672  for ( uint16_type n=0; n < element_type::numVertices; n++ )
1673  {
1674  auto const& old_point = old_elem.point( n );
1675 
1676  //if ( !new_node_numbers[old_point.id()] )
1677  if ( new_node_numbers.find( old_point.id() ) == new_node_numbers.end() )
1678  {
1679  new_node_numbers[old_point.id()] = n_new_nodes;
1680  DVLOG(2) << "[Mesh<Shape,T>::createP1mesh] insert point " << old_point << "\n";
1681  typename P1_mesh_type::point_type pt( old_point );
1682  pt.setId( n_new_nodes );
1683  pt.setProcessId(old_point.processId());
1684  // Add this node to the new mesh
1685  new_mesh->addPoint( pt );
1686  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] number of points " << new_mesh->numPoints() << "\n";
1687  // Increment the new node counter
1688  n_new_nodes++;
1689  FEELPP_ASSERT( !new_vertex[old_point.id()] ).error( "already seen this point?" );
1690  new_vertex[old_point.id()]=1;
1691  }
1692  // Define this element's connectivity on the new mesh
1693  //FEELPP_ASSERT ( new_node_numbers[old_elem.point( n ).id()] < new_mesh->numPoints() ).error( "invalid connectivity" );
1694  DVLOG(2) << "[Mesh<Shape,T>::createP1mesh] adding point old(" << old_point.id()
1695  << ") as point new(" << new_node_numbers[old_point.id()]
1696  << ") in element " << new_elem.id() << "\n";
1697  // add point in element
1698  new_elem.setPoint( n, new_mesh->point( new_node_numbers[old_point.id()] ) );
1699  } //for ( uint16_type n=0; n < element_type::numVertices; n++ )
1700 
1701  // Add an equivalent element type to the new_mesh
1702  new_mesh->addElement( new_elem );
1703  // increment the new element counter
1704  n_new_elem++;
1705 
1706  // Maybe add faces for this element
1707  for ( unsigned int s=0; s<old_elem.numTopologicalFaces; s++ )
1708  {
1709  if ( !old_elem.facePtr( s ) ) continue;
1710  // only add face on the boundary: they have some data
1711  // (boundary ids) which cannot be retrieved otherwise
1712  const size_type global_face_id = old_elem.face( s ).id();
1713  if ( this->hasFace( global_face_id ) )
1714  {
1715  // get the corresponding face
1716  face_type const& old_face = old_elem.face( s );
1717  typename P1_mesh_type::face_type new_face;
1718  // disconnect from elements of old mesh,
1719  // the connection will be redone in updateForUse()
1720  new_face.disconnect();
1721  // is on boundary
1722  new_face.setOnBoundary( old_face.isOnBoundary() );
1723  // set id of face
1724  new_face.setId( n_new_faces );
1725  // set face markers
1726  new_face.setMarker( old_face.marker().value() );
1727  new_face.setMarker2( old_face.marker2().value() );
1728  new_face.setMarker2( old_face.marker3().value() );
1729  // partitioning update
1730  new_face.setProcessIdInPartition( old_face.pidInPartition() );
1731  new_face.setNumberOfPartitions(old_face.numberOfPartitions());
1732  new_face.setProcessId(old_face.processId());
1733  //new_face.setIdInPartition( old_face.pidInPartition(), n_new_faces );
1734  new_face.setNeighborPartitionIds(old_face.neighborPartitionIds());
1735  // update P1 points info
1736  for ( uint16_type p = 0; p < face_type::numVertices; ++p )
1737  {
1738  new_face.setPoint( p, new_mesh->point( new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).id()] ) );
1739  }
1740  // add it to the list of faces
1741  new_mesh->addFace( new_face );
1742  // increment the new face counter
1743  ++n_new_faces;
1744  } // if ( this->hasFace( global_face_id ) )
1745  } // for ( unsigned int s=0; s<old_elem.numTopologicalFaces; s++ )
1746 
1747 
1748  if ( it->isGhostCell() )
1749  {
1750  DVLOG(2) << "element " << it->id() << " is a ghost cell\n";
1751  for (auto it_pid=it->idInOthersPartitions().begin(),en_pid=it->idInOthersPartitions().end() ; it_pid!=en_pid ; ++it_pid)
1752  {
1753  DVLOG(2) << " " << it_pid->first << "-" << it_pid->second << "-"<<it->pidInPartition()<<"-"<<new_mesh->worldComm().localRank();
1754  const int procToSend=it_pid->first;
1755  if (procToSend!=it->pidInPartition())
1756  {
1757  memory_ghostid[procToSend].push_back(boost::make_tuple(new_elem.id()/*it->id()*/,it_pid->second));
1758  }
1759  }
1760  }
1761  } // end for it
1762 
1763 
1764  if ( nProc > 1 )
1765  {
1766  // init container to send
1767  for ( int proc=0; proc<nProc; ++proc )
1768  {
1769  vecToSend[proc].resize(memory_ghostid[proc].size());
1770  memory_id[proc].resize(memory_ghostid[proc].size());
1771  auto it_mem = memory_ghostid[proc].begin();
1772  auto en_mem = memory_ghostid[proc].end();
1773  for ( int k = 0 ; it_mem!=en_mem ; ++k,++it_mem )
1774  {
1775  vecToSend[proc][k] = it_mem->get<1>();
1776  memory_id[proc][k] = it_mem->get<0>();
1777  }
1778  }
1779 
1780  //------------------------------------------------------
1781  this->worldComm().localComm().barrier();
1782  //------------------------------------------------------
1783 
1784  std::vector<size_type> nDataSize_vec(nProc);
1785  for ( int proc=0; proc<nProc; ++proc )
1786  {
1787  const auto nDataSize = vecToSend[proc].size();
1788  mpi::all_gather( this->worldComm().localComm(),
1789  nDataSize,
1790  nDataSize_vec );
1791 
1792  for ( int proc2=0; proc2<nProc; ++proc2 )
1793  {
1794  if ( nDataSize_vec[proc2] > 0 && proc!=proc2 )
1795  {
1796  if (proc_id==proc2) // send
1797  {
1798  this->worldComm().localComm().send( proc, 0, vecToSend[proc] );
1799  }
1800  else if (proc_id==proc) // recv
1801  {
1802  this->worldComm().localComm().recv( proc2, 0, vecToRecv[proc2] );
1803  }
1804  }
1805  }
1806  }
1807 
1808  //------------------------------------------------------
1809  this->worldComm().localComm().barrier();
1810  //------------------------------------------------------
1811 
1812  for ( int proc=0 ; proc<nProc; ++proc )
1813  {
1814  vecToSend[proc].resize(vecToRecv[proc].size());
1815  for (int k=0;k<vecToRecv[proc].size();++k)
1816  {
1817  const size_type idRequest = vecToRecv[proc][k];
1818  vecToSend[proc][k]= new_element_numbers[idRequest];
1819  }
1820  }
1821 
1822  //------------------------------------------------------
1823  this->worldComm().localComm().barrier();
1824  //------------------------------------------------------
1825 
1826  for ( int proc=0; proc<nProc; ++proc )
1827  {
1828  const auto nDataSize = vecToSend[proc].size();
1829  mpi::all_gather( this->worldComm().localComm(),
1830  nDataSize,
1831  nDataSize_vec );
1832 
1833  for ( int proc2=0; proc2<nProc; ++proc2 )
1834  {
1835  if ( nDataSize_vec[proc2] > 0 && proc!=proc2 )
1836  {
1837  if (proc_id==proc2) // send
1838  {
1839  this->worldComm().localComm().send( proc, 0, vecToSend[proc] );
1840  }
1841  else if (proc_id==proc) // recv
1842  {
1843  this->worldComm().localComm().recv( proc2, 0, vecToRecv[proc2] );
1844  }
1845  }
1846  }
1847  }
1848 
1849 
1850  //------------------------------------------------------
1851  this->worldComm().localComm().barrier();
1852  //------------------------------------------------------
1853 
1854  for ( int proc=0 ; proc<nProc; ++proc )
1855  {
1856  //std::cout << " vecToRecv[proc].size() " << vecToRecv[proc].size() << std::endl;
1857  for (int k=0;k<vecToRecv[proc].size();++k)
1858  {
1859  auto elttt = new_mesh->elementIterator( memory_id[proc][k], /*new_mesh->worldComm().localRank()*/ proc );
1860  new_mesh->elements().modify( elttt, detail::updateIdInOthersPartitions( proc, vecToRecv[proc][k] ) );
1861  }
1862  }
1863 
1864 
1865  } // if ( nProc > 1 )
1866 
1867 
1868 
1869  new_mesh->setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0,
1870  //[]( size_type lhs, std::pair<size_type,int> const& rhs )
1871  []( size_type lhs, std::pair<size_type,int> const& rhs )
1872  {
1873  return lhs+rhs.second;
1874  } ) );
1875 
1876 
1877  DVLOG(2) << "[Mesh<Shape,T>::createP1mesh] update face/edge info if necessary\n";
1878  // Prepare the new_mesh for use
1879  new_mesh->components().set ( /*MESH_RENUMBER|*/MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
1880  // run intensive job
1881  new_mesh->updateForUse();
1882 
1883  return new_mesh;
1884 }
1885 namespace detail
1886 {
1887 template<typename T>
1888 struct MeshPoints
1889 {
1890  template<typename MeshType, typename IteratorType>
1891  MeshPoints( MeshType*, IteratorType it, IteratorType en, const bool outer = false, const bool renumber = false );
1892 
1893  std::vector<int> ids;
1894  std::map<int,int> new2old;
1895  std::map<int,int> old2new;
1896  std::map<int,int> nodemap;
1897  std::vector<T> coords;
1898 };
1899 template<typename T>
1900 template<typename MeshType, typename IteratorType>
1901 MeshPoints<T>::MeshPoints( MeshType* mesh, IteratorType it, IteratorType en, const bool outer, const bool renumber )
1902 {
1903  std::set<int> nodeset;
1904  size_type p = 0;
1905  for( ; it != en; ++it )
1906  {
1907  for ( size_type j = 0; j < it->numLocalVertices; j++ )
1908  {
1909  int pid = it->point( j ).id();
1910  auto ins = nodeset.insert( pid );
1911  if ( ins.second )
1912  {
1913  if ( renumber )
1914  ids.push_back( p+1 );
1915  else
1916  ids.push_back( pid );
1917  old2new[pid]=ids[p];
1918  new2old[ids[p]]=pid;
1919  nodemap[pid] = p;
1920  ++p;
1921  }
1922  }
1923  }
1924  CHECK( p == ids.size() ) << "Invalid number of points " << ids.size() << "!=" << p;
1925  int nv = ids.size();
1926 
1927  coords.resize( 3*nv, 0 );
1928 
1929  auto pit = ids.begin();
1930  auto pen = ids.end();
1931  //for( auto i = 0; i < nv; ++i )
1932  for( int i = 0; pit != pen; ++pit, ++i )
1933  {
1934  CHECK( *pit > 0 ) << "invalid id " << *pit;
1935  //LOG(INFO) << "p " << i << "/" << nv << " =" << *pit;
1936  //int pid = (renumber)?nodemap[*pit]+1:*pit;
1937  int pid = *pit;
1938 
1939  auto const& p = mesh->point( new2old[*pit] );
1940  if ( outer )
1941  coords[i] = ( T ) p.node()[0];
1942  else
1943  coords[3*i] = ( T ) p.node()[0];
1944 
1945  if ( MeshType::nRealDim >= 2 )
1946  {
1947  if ( outer )
1948  coords[nv+i] = ( T ) p.node()[1];
1949  else
1950  coords[3*i+1] = ( T ) p.node()[1];
1951  }
1952  if ( MeshType::nRealDim >= 3 )
1953  {
1954  if ( outer )
1955  coords[2*nv+i] = T( p.node()[2] );
1956  else
1957  coords[3*i+2] = T( p.node()[2] );
1958  }
1959  }
1960 }
1961 }
1962 
1963 } // Feel
1964 
1965 
1966 //#if !defined(FEELPP_INSTANTIATION_MODE)
1967 # include <feel/feeldiscr/meshimpl.hpp>
1968 //#endif //
1969 
1970 #endif /* __mesh_H */
reference_convex_type referenceConvex() const
Definition: mesh.hpp:337
void encode()
Definition: meshimpl.hpp:1839
void updateMarker2WithRange(IteratorRange const &range, flag_type flag, mpl::int_< MESH_ELEMENTS >)
Definition: mesh.hpp:589
simplex of dimension Dim
Definition: simplex.hpp:73
element_iterator eraseElement(element_iterator position, bool modify=true)
Definition: meshimpl.hpp:1809
void findNeighboringProcessors()
Definition: meshimpl.hpp:1609
void send(int p, int tag)
Definition: meshimpl.hpp:1765
void check() const
Definition: meshimpl.hpp:1529
void addMarkerName(std::pair< std::string, std::vector< size_type > > const &marker)
Definition: mesh.hpp:473
Definition: mesh.hpp:109
size_type numGlobalElements() const
Definition: mesh.hpp:281
void createSubmeshByProcessId(self_type &mesh, uint16_type pid) const
Definition: mesh.hpp:561
void localrenumber()
Definition: meshimpl.hpp:709
void updateMarker3WithRange(IteratorRange const &range, flag_type flag, mpl::int_< MESH_FACES >)
Definition: mesh.hpp:629
trace_mesh< TheTag >::ptrtype trace() const
Definition: mesh.hpp:522
face_processor_type & localFaceId(element_type const &e, size_type const n)
Definition: mesh.hpp:395
void updateMarker3WithRange(IteratorRange const &range, flag_type flag, mpl::int_< MESH_ELEMENTS >)
Definition: mesh.hpp:620
value_type measureBoundary() const
Definition: mesh.hpp:459
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
void addMarkerName(std::string __name, int __id, int __topoDim)
Definition: mesh.hpp:481
void recv(int p, int tag)
Definition: meshimpl.hpp:1782
face_processor_type const & localFaceId(size_type const e, size_type const n) const
Definition: mesh.hpp:354
Definition: glas.hpp:125
void updateEntitiesCoDimensionOneGhostCell()
void updateForUse()
Definition: meshimpl.hpp:96
unifying mesh class
Definition: createsubmesh.hpp:41
void updateEntitiesCoDimensionOne()
Definition: meshimpl.hpp:753
static mesh_ptrtype New()
Definition: mesh.hpp:257
void updateMarker2WithRange(IteratorRange const &range, flag_type flag)
Definition: mesh.hpp:579
void setPartitioner(std::string partitioner)
Definition: mesh.hpp:387
void updateMarker2WithRange(IteratorRange const &range, flag_type flag, mpl::int_< MESH_FACES >)
Definition: mesh.hpp:599
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::pid_face_const_iterator, typename MeshTraits< MeshType >::pid_face_const_iterator > faces(MeshType const &mesh)
Definition: filters.hpp:933
void checkLocalPermutation(mpl::bool_< false >) const
Definition: mesh.hpp:832
Mesh(WorldComm const &worldComm=Environment::worldComm())
Definition: meshimpl.hpp:45
std::map< std::string, std::vector< size_type > > markerNames() const
Definition: mesh.hpp:434
size_type markerName(std::string const &marker) const
Definition: mesh.hpp:413
gm1_ptrtype const & gm1() const
Definition: mesh.hpp:312
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
uint16_type dimension() const
Definition: mesh.hpp:296
void updateMarker3WithRange(IteratorRange const &range, flag_type flag)
Definition: mesh.hpp:610
void renumber()
Definition: mesh.hpp:647
boost::signals2::signal< void(MESH_CHANGES)> meshChanged
Definition: mesh.hpp:1260
gm_ptrtype & gm()
Definition: mesh.hpp:320
gm_ptrtype const & gm() const
Definition: mesh.hpp:304
face_processor_type & localFaceId(size_type const e, size_type const n)
Definition: mesh.hpp:404
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::location_face_const_iterator, typename MeshTraits< MeshType >::location_face_const_iterator > boundaryfaces(MeshType const &mesh)
Definition: filters.hpp:1158
P1_mesh_ptrtype createP1mesh() const
Definition: mesh.hpp:1614
face_processor_type const & localFaceId(element_type const &e, size_type const n) const
Definition: mesh.hpp:345
gm1_ptrtype & gm1()
Definition: mesh.hpp:328
Context class.
Definition: feelcore/context.hpp:63
size_type markerDim(std::string const &marker) const
Definition: mesh.hpp:423
value_type measure() const
Definition: mesh.hpp:451
#define FEELPP_DEFINE_VISITABLE()
Definition: visitor.hpp:319
void partition(const uint16_type n_parts=1)
Definition: meshimpl.hpp:59
void checkAndFixPermutation()
Definition: meshimpl.hpp:1704
void decode()
Definition: meshimpl.hpp:1912

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