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
doftable.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  This file is part of the Feel library
3 
4  Copyright (C) 2010 Université de Grenoble 1
5 
6  This library is free software; you can redistribute it and/or
7  modify it under the terms of the GNU Lesser General Public
8  License as published by the Free Software Foundation; either
9  version 3.0 of the License, or (at your option) any later version.
10 
11  This library is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public
17  License along with this library; if not, write to the Free Software
18  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
24 #ifndef _DOFTABLE_HH
25 #define _DOFTABLE_HH
26 
27 #include <set>
28 #include <map>
29 #include <unordered_map>
30 #include <vector>
31 #include <algorithm>
32 
33 #include <boost/foreach.hpp>
34 #include <boost/multi_array.hpp>
35 #include <boost/tuple/tuple.hpp>
36 #include <boost/tuple/tuple_comparison.hpp>
37 #include <boost/tuple/tuple_io.hpp>
38 #include <boost/fusion/algorithm/iteration/accumulate.hpp>
39 #include <boost/bimap.hpp>
40 #include <boost/bimap/support/lambda.hpp>
41 #include <boost/bimap/unordered_set_of.hpp>
42 #include <boost/bimap/multiset_of.hpp>
43 #include <boost/bimap/set_of.hpp>
44 
45 #include <Eigen/Core>
46 #include<Eigen/StdVector>
47 
48 #include <feel/feelcore/feel.hpp>
50 #include <feel/feelalg/glas.hpp>
51 #include <feel/feelpoly/mapped.hpp>
52 #include <feel/feelalg/datamap.hpp>
53 #include <feel/feeldiscr/dof.hpp>
54 
55 namespace Feel
56 {
57 // import fusion namespace in Feel
58 namespace fusion = boost::fusion;
59 namespace bimaps = boost::bimaps;
68 template<typename MeshType, typename FEType, typename PeriodicityType>
69 class DofTable : public DataMap
70 {
71  typedef DataMap super;
72 public:
73 
77  typedef MeshType mesh_type;
78  typedef FEType fe_type;
79  typedef boost::shared_ptr<FEType> fe_ptrtype;
80 
81 
82  typedef typename mesh_type::pid_element_const_iterator pid_element_const_iterator;
83  typedef typename mesh_type::element_const_iterator element_const_iterator;
84  typedef typename mesh_type::element_type element_type;
85  typedef typename mesh_type::face_type face_type;
86  typedef typename mesh_type::gm_ptrtype gm_ptrtype;
87  typedef typename mesh_type::gm_type gm_type;
88 
89  typedef typename fe_type::matrix_type matrix_type;
90  typedef typename fe_type::value_type value_type;
91  typedef typename fe_type::reference_convex_type reference_convex_type;
92  typedef typename fe_type::points_type points_type;
93 
94  typedef typename reference_convex_type::super convex_type;
95 
96  typedef typename element_type::edge_permutation_type edge_permutation_type;
97  typedef typename element_type::face_permutation_type face_permutation_type;
98 
99  static const uint16_type nOrder = fe_type::nOrder;
100  static const uint16_type nDim = mesh_type::nDim;
101  static const uint16_type nRealDim = mesh_type::nRealDim;
102  static const uint16_type Shape = mesh_type::Shape;
103  static const uint16_type nComponents = fe_type::nComponents;
104  static const uint16_type nComponents1 = fe_type::nComponents1;
105  static const uint16_type nComponents2 = fe_type::nComponents2;
106 
107 
108  static const bool is_continuous = fe_type::isContinuous;
109  static const bool is_discontinuous_locally = fe_type::continuity_type::is_discontinuous_locally;
110  static const bool is_discontinuous_totally = fe_type::continuity_type::is_discontinuous_totally;
111 
112  static const bool is_scalar = FEType::is_scalar;
113  static const bool is_vectorial = FEType::is_vectorial;
114  static const bool is_tensor2 = FEType::is_tensor2;
115  static const bool is_modal = FEType::is_modal;
116  static const bool is_product = FEType::is_product;
117 
118  static const bool is_p0_continuous = ( ( nOrder == 0 ) && is_continuous );
119 
120  static const uint16_type nDofPerElement = mpl::if_<mpl::bool_<is_product>, mpl::int_<FEType::nLocalDof*nComponents1>, mpl::int_<FEType::nLocalDof> >::type::value;
121 
122  typedef PeriodicityType periodicity_type;
123  static const bool is_periodic = periodicity_type::is_periodic;
124 
125  //typedef ContinuityType continuity_type;
126  typedef typename fe_type::continuity_type continuity_type;
127 
129 
130 
136  //typedef boost::tuple<size_type, int16_type, bool> global_dof_type;
138 
139 
148  typedef boost::tuple<size_type, int16_type, bool, int16_type> global_dof_fromface_type;
149 
150 
154  //typedef std::unordered_map<int,std::map<int,global_dof_type> > Container;
155  //typedef typename std::map<int,global_dof_type>::iterator local_map_iterator;
156  typedef boost::bimap<bimaps::unordered_set_of<LocalDof>, bimaps::multiset_of<Dof> > dof_table;
157  typedef dof_table::value_type dof_relation;
158  typedef std::unordered_map<int,std::map<int,global_dof_fromface_type> > Container_fromface;
159 
160  typedef typename std::map<int,global_dof_type> indices_per_element_type;
161 
162  typedef typename node<value_type>::type node_type;
163 
164 
165  typedef boost::tuple<node_type, size_type, uint16_type > dof_point_type;
166  typedef std::vector<dof_point_type> dof_points_type;
167  typedef typename std::vector<dof_point_type>::iterator dof_points_iterator;
168  typedef typename std::vector<dof_point_type>::const_iterator dof_points_const_iterator;
169 
177  typedef boost::tuple<size_type, uint16_type, uint16_type, uint16_type> local_dof_type;
178  typedef LocalDofSet local_dof_set_type;
179 
180  typedef boost::tuple<uint16_type&,size_type&> ref_shift_type;
181 
185  typedef std::map<size_type, std::list<local_dof_type> > dof_element_type;
186 
187  typedef boost::bimap<size_type,boost::bimaps::multiset_of<size_type> > dof_marker_type;
188  typedef typename dof_marker_type::value_type dof2marker;
189 
190  typedef typename dof_element_type::iterator dof_iterator;
191  typedef typename dof_element_type::const_iterator dof_const_iterator;
192 
193  typedef std::list<local_dof_type>::const_iterator ldof_const_iterator;
194 
195  typedef boost::tuple<uint16_type,uint16_type,size_type> dof_type;
196  typedef std::map<dof_type, size_type> dof_map_type;
197  typedef std::map<dof_type, size_type>::iterator dof_map_iterator;
198  typedef std::map<dof_type, size_type, size_type>::const_iterator dof_map_const_iterator;
199 
200  typedef std::map<size_type, std::set<size_type> > dof_procset_type;
205  typedef ublas::vector<bool> face_sign_info_type;
206 
207 
208  typedef Eigen::Matrix<int, nDofPerElement, 1> localglobal_indices_type;
209 
214  typedef ublas::vector<uint16_type> permutation_vector_type;
215 
216  typedef boost::tuple<element_type const*, face_type const*> element_face_pair_type;
217  typedef std::list<element_face_pair_type> periodic_element_list_type;
218  typedef typename periodic_element_list_type::iterator periodic_element_list_iterator;
219  typedef typename periodic_element_list_type::const_iterator periodic_element_list_const_iterator;
220  typedef boost::tuple<size_type /*element id*/, uint16_type /*lid*/, uint16_type /*c*/, size_type /*gDof*/, uint16_type /*type*/> periodic_dof_type;
221  typedef std::multimap<size_type /*gid*/, periodic_dof_type> periodic_dof_map_type;
222 
223  DofTable( WorldComm const& _worldComm )
224  :
225  super( _worldComm )
226  {}
227 
234  DofTable( fe_ptrtype const& _fe, periodicity_type const& periodicity, WorldComm const& _worldComm );
235 
241  DofTable( const DofTable & dof2 );
242 
249  DofTable( mesh_type& mesh, fe_ptrtype const& _fe, periodicity_type const& periodicity, WorldComm const& _worldComm );
250 
251  ~DofTable()
252  {
253  M_el_l2g.clear();
254  M_face_l2g.clear();
255  M_dof_points.clear();
256  }
257  constexpr size_type nLocalDof( bool per_component = false ) const
258  {
259  return (is_product&&!per_component)?(nComponents*(fe_type::nDofPerVolume * element_type::numVolumes +
260  fe_type::nDofPerFace * element_type::numGeometricFaces +
261  fe_type::nDofPerEdge * element_type::numEdges +
262  fe_type::nDofPerVertex * element_type::numVertices)):
263  (fe_type::nDofPerVolume * element_type::numVolumes +
264  fe_type::nDofPerFace * element_type::numGeometricFaces +
265  fe_type::nDofPerEdge * element_type::numEdges +
266  fe_type::nDofPerVertex * element_type::numVertices);
267  }
268  constexpr size_type nLocalDofOnFace( bool per_component = false ) const
269  {
270  return (is_product&&!per_component)?(nComponents*( face_type::numVertices * fe_type::nDofPerVertex +
271  face_type::numEdges * fe_type::nDofPerEdge +
272  face_type::numFaces * fe_type::nDofPerFace )):
273  ( face_type::numVertices * fe_type::nDofPerVertex +
274  face_type::numEdges * fe_type::nDofPerEdge +
275  face_type::numFaces * fe_type::nDofPerFace );
276  }
277  local_dof_set_type const&
278  localDofSet( size_type eid ) const
279  {
280  return M_local_dof_set.update( eid );
281  }
282  mesh_type* mesh() { return M_mesh; }
283  mesh_type* mesh() const { return M_mesh; }
284 
288  uint16_type nDofPerFaceOnBoundary() const
289  {
290  return M_n_dof_per_face_on_bdy;
291  }
292 
293  indices_per_element_type indices( size_type id_el ) const
294  {
295  indices_per_element_type ind;
296  BOOST_FOREACH( LocalDof const& ldof, localDofSet( id_el ) )
297  {
298  auto it = M_el_l2g.left.find( ldof );
299  DCHECK( it != M_el_l2g.left.end() ) << "Invalid element id " << id_el;
300  ind[ldof.localDof()] = *it;
301  }
302  return ind;
303  }
304 
305  constexpr size_type getIndicesSize() const
306  {
307  return nLocalDof();
308  }
309  std::vector<size_type> getIndices( size_type id_el ) const
310  {
311  std::vector<size_type> ind( nLocalDof() );
312  getIndicesSet( id_el, ind );
313 
314  return ind;
315  }
316 
317  std::vector<size_type> getIndices( size_type id_el, mpl::size_t<MESH_ELEMENTS> ) const
318  {
319  return getIndices( id_el );
320  }
321 
322  void getIndicesSet( size_type id_el, std::vector<size_type>& ind ) const
323  {
324  BOOST_FOREACH( LocalDof const& ldof, this->localDofSet( id_el ) )
325  {
326  auto it = M_el_l2g.left.find( ldof );
327  DCHECK(it != M_el_l2g.left.end() ) << "Invalid element id " << id_el;
328  ind[ldof.localDof()] = it->second.index();
329  }
330  }
331 
332  std::vector<size_type> getIndices( size_type id_el, mpl::size_t<MESH_FACES> ) const
333  {
334  std::vector<size_type> ind( nLocalDofOnFace() );
335 
336  auto eit = M_face_l2g.find( id_el );
337  DCHECK( eit != M_face_l2g.end() ) << "Invalid face id " << id_el;
338  auto it = eit->second.begin();
339  auto en = eit->second.end();
340  for( ; it != en; ++ it )
341  ind[it->first] = boost::get<0>( it->second );
342 
343  return ind;
344  }
345 
346 
347  void getIndicesSetOnGlobalCluster( size_type id_el, std::vector<size_type>& ind ) const
348  {
349  BOOST_FOREACH( LocalDof const& ldof, this->localDofSet( id_el ) )
350  {
351  auto it = M_el_l2g.left.find( ldof );
352  DCHECK( it != M_el_l2g.left.end() ) << "Invalid element id " << id_el;
353  ind[ldof.localDof()] =this->mapGlobalProcessToGlobalCluster()[ it->second.index() ];
354 
355  }
356  }
357 
358  std::vector<size_type> getIndicesOnGlobalCluster( size_type id_el ) const
359  {
360  const size_type s = getIndicesSize();
361  std::vector<size_type> ind( s );
362  getIndicesSetOnGlobalCluster( id_el, ind );
363  return ind;
364  }
365 
366 
371  {
372  return M_dof_procset.find( dof )->second.size();
373  }
374 
378  std::set<size_type> dofProcSet( size_type dof ) const
379  {
380  return M_dof_procset.find( dof )->second;
381  }
382 
387  const dof_point_type& dofPoint( size_type i ) const
388  {
389  return M_dof_points[i];
390  }
391 
395  dof_points_const_iterator dofPointBegin() const
396  {
397  return M_dof_points.begin();
398  }
399 
403  dof_points_iterator dofPointBegin()
404  {
405  return M_dof_points.begin();
406  }
407 
411  dof_points_const_iterator dofPointEnd() const
412  {
413  return M_dof_points.end();
414  }
415 
419  dof_points_iterator dofPointEnd()
420  {
421  return M_dof_points.end();
422  }
423 
424  periodic_element_list_const_iterator beginPeriodicElements() const { return periodic_elements.begin(); }
425  periodic_element_list_const_iterator endPeriodicElements() const { return periodic_elements.end(); }
426 
434  void setDofIndices( std::vector<Dof> const& dof )
435  {
436  M_dof_indices.resize( dof.size() );
437  std::copy( dof.begin(), dof.end(), M_dof_indices.begin() );
438 
439  if ( dof.empty() )
440  return ;
441 
442 #if 1
443  std::set<size_type> eltid;
444  std::set<size_type> dofs;
445 
446  BOOST_FOREACH( Dof thedof, dof )
447  {
448  eltid.insert( boost::get<0>( thedof ) );
449  dofs.insert( boost::get<2>( thedof ) );
450  }
451 #endif
452 
453  BOOST_FOREACH( Dof thedof, dof )
454  {
455  M_el_l2g.insert( dof_relation( LocalDof( boost::get<0>( thedof ), boost::get<1>( thedof ) ),
456  Dof( boost::get<2>( thedof ), 0, false ) ) );
457 
458  }
459  int processor = this->worldComm().localRank();
460 
461  this->M_first_df[processor] = 0;
462  this->M_last_df[processor] = dofs.size()-1;
463  this->M_n_dofs = dofs.size();
464 
465  this->M_n_localWithGhost_df[processor] = this->M_last_df[processor] - this->M_first_df[processor] + 1;
466  this->M_n_localWithoutGhost_df[processor]=this->M_n_localWithGhost_df[processor];
467  this->M_first_df_globalcluster[processor]=this->M_first_df[processor];
468  this->M_last_df_globalcluster[processor]=this->M_last_df[processor];
469  }
470 
475  {
476  return dof;
477 #if 0
478 
479  if ( M_dof_indices.empty() )
480  return dof;
481 
482  FEELPP_ASSERT( dof < M_dof_indices.size() )( dof )( M_dof_indices.size() ).warn( "invalid dof index" );
483  return M_dof_indices[dof];
484 #endif
485  }
486 
490  localglobal_indices_type const& localToGlobalIndices( size_type ElId )
491  {
492  return M_locglob_indices[ElId];
493  }
494 
498  localglobal_indices_type const& localToGlobalIndicesOnCluster( size_type ElId )
499  {
500  return M_locglobOnCluster_indices[ElId];
501  }
502 
506  localglobal_indices_type const& localToGlobalSigns( size_type ElId )
507  {
508  return M_locglob_signs[ElId];
509  }
510 
521  const uint16_type id ) const
522  {
523  auto it = M_el_l2g.left.find( LocalDof(ElId, id ) );
524  DCHECK( it != M_el_l2g.left.end() ) << "Invalid dof entry ( " << ElId << ", " << id << ")";
525  DCHECK( it->second.index() < nDof() ) << "Invalid Dof Entry: " << it->second.index() << " > " << this->nDof();
526  return it->second.index();
527  }
528 
536  LocalDof const& globalToLocal( size_type dof ) const
537  {
538  auto it = M_el_l2g.right.find( Dof( dof ) );
539  DCHECK( it != M_el_l2g.right.end() ) << "Invalid global dof entry ( " << dof << ")";
540  return it->second;
541  }
542 
543  uint16_type localDofId( uint16_type const lid, uint16_type const c = 0 ) const
544  {
545  return fe_type::nLocalDof * c + lid;
546  }
547  global_dof_type const& localToGlobal( const size_type ElId,
548  const uint16_type localNode,
549  const uint16_type c = 0 ) const
550  {
551  auto it = M_el_l2g.left.find( LocalDof(ElId,fe_type::nLocalDof * c + localNode ) );
552  DCHECK( it != M_el_l2g.left.end() ) << "Invalid dof entry ( " << ElId << ", " << fe_type::nLocalDof * c + localNode << ")";
553  //DCHECK( it->second.index() < nDof() && nDof() > 0 ) << "Invalid Dof Entry: " << it->second.index() << " > " << this->nDof();
554  return it->second;
555  }
556 
557  global_dof_type localToGlobalOnCluster( const size_type ElId,
558  const uint16_type localNode,
559  const uint16_type c = 0 ) const
560  {
561  Dof resloc = M_el_l2g.left.find( LocalDof( ElId, fe_type::nLocalDof * c + localNode ) )->second;
562  resloc.setIndex( this->mapGlobalProcessToGlobalCluster()[resloc.index()] );
563  return resloc;
564  }
565 
566  global_dof_fromface_type const& faceLocalToGlobal( const size_type ElId,
567  const uint16_type localNode,
568  const uint16_type c = 0 ) const
569  {
570  const size_type nDofF = nLocalDofOnFace( true );
571  return M_face_l2g.find( ElId )->second.find( nDofF*c+localNode )->second;
572  }
573 
574  struct element_access
575  {
576  element_access( DofTable const& __d )
577  :
578  M_d( __d )
579  {}
580  global_dof_type const& operator()( size_type __id, uint16_type __loc, uint16_type c = 0 ) const
581  {
582  return M_d.M_el_l2g.left.find( LocalDof(__id, M_d.nLocalDof(true) * c+ __loc ) )->second;
583  }
584  uint16_type localDofInElement( size_type __id, uint16_type __loc, uint16_type c = 0 ) const
585  {
586  return __loc;
587  }
588  DofTable const& M_d;
589  };
590  friend struct element_access;
591 
592  struct face_access
593  {
594 
595  face_access( DofTable const& __d )
596  :
597  M_d( __d )
598  {}
599  global_dof_fromface_type operator()( size_type __id, uint16_type __loc, uint16_type c = 0 ) const
600  {
601  return M_d.M_face_l2g.find( __id)->second.find(M_d.nLocalDofOnFace( true )*c+__loc)->second;
602  }
603 
604  uint16_type localDofInElement( size_type __id, uint16_type __loc, uint16_type c = 0 ) const
605  {
606  return boost::get<3>( M_d.M_face_l2g.find( __id )->second.find( M_d.nLocalDofOnFace( true )*c+__loc)->second );
607  }
608 
609  DofTable const& M_d;
610  };
611  friend struct face_access;
612 
616  template<typename Elem>
617  typename mpl::if_<mpl::equal_to<mpl::int_<Elem::nDim>,mpl::int_<nDim> >,
620  localToGlobal( Elem const& El, const uint16_type localNode, const uint16_type c = 0 ) const
621  {
622  typedef typename mpl::if_<mpl::equal_to<mpl::int_<Elem::nDim>,mpl::int_<nDim> >,mpl::identity<element_access>,mpl::identity<face_access> >::type::type access_type;
623  //DVLOG(2) << "dof:(" << El.id() << ", " << localNode << ")= "
624  //<< access_type(*this)( El.id(), localNode, c ) << "\n";
625  return access_type( *this )( El.id(), localNode, c );
626  }
627 
628  template<typename Elem>
629  uint16_type
630  localDofInElement( Elem const& El, const uint16_type localNode, const uint16_type c = 0 ) const
631  {
632  typedef typename mpl::if_<mpl::equal_to<mpl::int_<Elem::nDim>,mpl::int_<nDim> >,mpl::identity<element_access>,mpl::identity<face_access> >::type::type access_type;
633 
634  return ( access_type( *this ) ).localDofInElement( El.id(), localNode, c );
635  }
636 
641  {
642  return M_n_el;
643  }
644 
648  uint16_type numLocalVertices() const
649  {
650  return fe_type::numVertices;
651  }
652 
656  uint16_type numLocalEdges() const
657  {
658  return fe_type::numEdges;
659  }
660 
664  uint16_type numLocalFaces() const
665  {
666  return fe_type::numFaces;
667  }
668 
672  void showMe() const;
673 
674  void dump() const
675  {
676 #if 0
677 
678  for ( size_type __i = 0; __i < M_face_l2g.nrows(); ++__i )
679  {
680  for ( size_type __l = 0; __l < M_face_l2g.ncols(); ++__l )
681  {
682  std::cout << "face " << __i << " local " << __l
683  << " to global " << M_face_l2g[ __i][ __l ] << "\n";
684  }
685  }
686 
687 #endif // 0
688  }
689 
694  {
695  M_elt_done[elt] = true;
696  }
697  bool isElementDone( size_type elt, int c = 0 )
698  {
699  bool done = true;
700  BOOST_FOREACH( LocalDof const& local_dof, this->localDofSet( elt ) )
701  {
702  auto it = M_el_l2g.left.find( local_dof );
703  if ( it == M_el_l2g.left.end() )
704  return false;
705  }
706  return done;
707  }
708  void addDofFromElement( element_type const& __elt,
709  size_type& next_free_dof,
710  size_type processor = 0,
711  size_type shift = 0 )
712  {
713 
714  size_type nldof = this->nLocalDof( true );
715 
716  if ( !this->isElementDone( __elt.id() ) )
717  {
718  DVLOG(3) << "adding dof from element " << __elt.id() << "\n";
719  size_type gdofcount = shift;
720  DVLOG(3) << "next_free_dof " << next_free_dof << "\n";
721  DVLOG(3) << "current dof " << dofIndex( next_free_dof ) << "\n";
722 
723 
724  /*
725  * Only in the continuous , we need to have the ordering [vertex,edge,face,volume]
726  */
727  if ( is_continuous || is_discontinuous_locally )
728  {
729 
730  /* idem as above but for local element
731  numbering except that it is
732  reset to 0 after each element */
733  uint16_type ldofcount = 0;
734 
735  /* pack the shifts into a tuple */
736  boost::tuple<uint16_type&,size_type&> shifts = boost::make_tuple( boost::ref( ldofcount ),
737  boost::ref( gdofcount ) );
738 
739  /* \warning: the order of function calls is
740  crucial here we order the degrees of freedom
741  wrt the topological entities of the mesh
742  elements from lowest dimension (vertex) to
743  highest dimension (element)
744  */
745  addVertexDof( __elt, processor, next_free_dof, shifts );
746  addEdgeDof( __elt, processor, next_free_dof, shifts );
747  addFaceDof( __elt, processor, next_free_dof, shifts );
748  addVolumeDof( __elt, processor, next_free_dof, shifts );
749  }
750 
751  else
752  {
753 
754  size_type ie = __elt.id();
755 
756  const int ncdof = is_product?nComponents:1;
757 
758  for ( uint16_type l = 0; l < nldof; ++l )
759  {
760  for ( int c = 0; c < ncdof; ++c, ++next_free_dof )
761  {
762  M_el_l2g.insert( dof_relation( LocalDof( ie, fe_type::nLocalDof*c + l ),
763  Dof( ( dofIndex( next_free_dof ) ) , 1, false ) ) );
764  }
765  }
766  }
767  }
768 
769  else
770  {
771  DVLOG(3) << "element " << __elt.id() << "has already been taken care of\n";
772  }
773  }
777  std::map<size_type, std::pair<size_type, size_type> >
779  {
780  return M_local2global;
781  }
782 
787  {
788  return M_local2global.find( l )->second.first;
789  }
790 
801  void addVertexPeriodicDof( element_type const& __elt,
802  face_type const& __face,
803  size_type& next_free_dof,
804  std::map<size_type,periodic_dof_map_type>& periodic_dof,
805  size_type tag )
806  {
807  addVertexPeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag, mpl::bool_<(fe_type::nDofPerVertex>0)>() );
808  }
809  void addVertexPeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof,size_type tag, mpl::bool_<false> ) {}
810  void addVertexPeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof,size_type tag, mpl::bool_<true> );
811 
812  void addEdgePeriodicDof( element_type const& __elt,
813  face_type const& __face,
814  size_type& next_free_dof,
815  std::map<size_type,periodic_dof_map_type>& periodic_dof,
816  size_type tag )
817  {
818  static const bool cond = fe_type::nDofPerEdge > 0;
819  addEdgePeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag , mpl::bool_<cond>(), mpl::int_<nDim>() );
820  }
821  void addEdgePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof, size_type tag, mpl::bool_<false>, mpl::int_<1> ) {}
822  void addEdgePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof, size_type tag, mpl::bool_<false>, mpl::int_<2> ) {}
823  void addEdgePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof, size_type tag, mpl::bool_<false>, mpl::int_<3> ) {}
824  void addEdgePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof, size_type tag, mpl::bool_<true>, mpl::int_<1> ) {}
825  void addEdgePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof, size_type tag, mpl::bool_<true>, mpl::int_<2> );
826  void addEdgePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof, size_type tag, mpl::bool_<true>, mpl::int_<3> );
827 
828 
829  void addFacePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof,size_type tag )
830  {
831  static const bool cond = fe_type::nDofPerFace > 0;
832  addFacePeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag, mpl::bool_<cond>() );
833  }
834  void addFacePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof,size_type tag, mpl::bool_<false> ) {}
835  void addFacePeriodicDof( element_type const& __elt,face_type const& __face,size_type& next_free_dof,std::map<size_type,periodic_dof_map_type>& periodic_dof,size_type tag, mpl::bool_<true> );
836 
840  void initDofMap( mesh_type& M );
841 
845  void build( mesh_type* M )
846  {
847  this->build( *M );
848  }
849 
853  void build( boost::shared_ptr<mesh_type>& M )
854  {
855  this->build( *M );
856  }
857 
861  void build( mesh_type& M )
862  {
863  M_mesh = boost::addressof( M );
864 
865  M_elt_done.resize( M.numElements() );
866  std::fill( M_elt_done.begin(), M_elt_done.end(), false );
867 
868  VLOG(2) << "[Dof::build] initDofMap\n";
869  this->initDofMap( M );
870 
871  VLOG(2) << "[Dof::build] start building dof map\n";
872  size_type start_next_free_dof = 0;
873  VLOG(2) << "[Dof::build] start_next_free_dof = " << start_next_free_dof << "\n";
874 
875  if ( is_periodic )
876  {
877  VLOG(2) << "[build] call buildPeriodicDofMap()\n";
878  start_next_free_dof = this->buildPeriodicDofMap( M );
879  VLOG(2) << "[Dof::build] start_next_free_dof(after periodic) = " << start_next_free_dof << "\n";
880  }
881 
882  if ( is_discontinuous_locally )
883  {
884  VLOG(2) << "[build] call buildLocallyDiscontinuousDofMap()\n";
885  start_next_free_dof = this->buildLocallyDiscontinuousDofMap( M, start_next_free_dof );
886  VLOG(2) << "[Dof::build] start_next_free_dof(after local discontinuities) = " << start_next_free_dof << "\n";
887  }
888 
889  VLOG(2) << "[build] call buildDofMap()\n";
890  this->buildDofMap( M, start_next_free_dof );
891  //std::cout << "[build] callFINISH buildDofMap() with god rank " << this->worldComm().godRank() <<"\n";
892 
893 #if !defined(NDEBUG)
894  VLOG(2) << "[build] check that all elements dof were assigned()\n";
895  element_const_iterator fit, fen;
896  boost::tie( fit, fen ) = M.elementsRange();
897  std::vector<boost::tuple<size_type,uint16_type,size_type> > em;
898 
899  for ( ; fit != fen; ++fit )
900  {
901  const int ncdof = is_product?nComponents:1;
902 
903  for ( uint16_type c = 0; c < ncdof; ++c )
904  if ( !this->isElementDone( fit->id(), c ) )
905  {
906  em.push_back( boost::make_tuple( fit->id(), c, fit->marker().value() ) );
907  }
908  }
909 
910  if ( !em.empty() )
911  {
912  VLOG(2) << "[build] some element dof were not assigned\n";
913 
914  for ( size_type i = 0; i < em.size(); ++i )
915  {
916  VLOG(3) << " - element " << boost::get<0>( em[i] ) << " c=" << boost::get<1>( em[i] )
917  << " m=" << boost::get<2>( em[i] ) << "\n";
918  }
919  }
920 
921  else
922  {
923  VLOG(2) << "[build] check that all elements dof were assigned: OK\n";
924  }
925 
926 #endif // NDEBUG
927  VLOG(2) << "[Dof::build] n_dof = " << this->nLocalDofWithGhost() << "\n";
928 
929  VLOG(2) << "[build] call buildBoundaryDofMap()\n";
930  this->buildBoundaryDofMap( M );
931 
932 
933  // multi process
934  if ( this->worldComm().localSize()>1 )
935  {
936  auto it_nlocDof = this->M_n_localWithGhost_df.begin();
937  auto const en_nlocDof = this->M_n_localWithGhost_df.end();
938  bool isP0continuous = true;
939  for ( ; it_nlocDof!=en_nlocDof && isP0continuous ; ++it_nlocDof )
940  isP0continuous = isP0continuous && (*it_nlocDof==1);
941 
942  if ( !isP0continuous )//this->M_n_dofs>1 )
943  {
944 #if defined(FEELPP_ENABLE_MPI_MODE)
945  VLOG(2) << "[build] call buildGhostDofMap () with god rank " << this->worldComm().godRank() << "\n";
946  this->buildGhostDofMap( M );
947  VLOG(2) << "[build] callFINISH buildGhostDofMap () with god rank " << this->worldComm().godRank() << "\n";
948 #else
949  std::cerr << "ERROR : FEELPP_ENABLE_MPI_MODE is OFF" << std::endl;
950  //throw std::logic_error( "ERROR : FEELPP_ENABLE_MPI_MODE is OFF" );
951 #endif
952  }
953  else
954  {
955  for ( int proc=0; proc<this->worldComm().globalSize(); ++proc )
956  {
957  if (proc==0)
958  {
959  this->M_n_localWithoutGhost_df[proc] = 1;
960  this->M_first_df_globalcluster[proc] = 0;
961  this->M_last_df_globalcluster[proc] = 0;
962  }
963  else
964  {
965  this->M_n_localWithoutGhost_df[proc] = 0;
966  this->M_first_df_globalcluster[proc] = 25;// 0;
967  this->M_last_df_globalcluster[proc] = 25; //0;
968  }
969  this->M_n_localWithGhost_df[proc] = 1;
970  }
971 
972  if (this->worldComm().globalRank()==0)
973  {
974  this->M_mapGlobalClusterToGlobalProcess.resize( 1 );
976  }
977  else
978  {
979  this->M_mapGlobalClusterToGlobalProcess.resize( 0 );
980  }
981 
982  this->M_mapGlobalProcessToGlobalCluster.resize( 1 );
984  }
985  }
986 
987  else
988  {
989 #if defined(FEELPP_ENABLE_MPI_MODE)
990  // in sequential : identity map
991  const size_type s = this->M_n_localWithGhost_df[this->comm().rank()];
992  this->M_mapGlobalProcessToGlobalCluster.resize( s );
993  this->M_mapGlobalClusterToGlobalProcess.resize( s );
994 
995  for ( size_type i=0; i<s ; ++i ) this->M_mapGlobalProcessToGlobalCluster[i]=i;
996 
997  for ( size_type i=0; i<s ; ++i ) this->M_mapGlobalClusterToGlobalProcess[i]=i;
998 
999 #endif
1000  }
1001 
1002  VLOG(2) << "[Dof::build] done building the map\n";
1003  }
1004 
1010 
1014  size_type buildLocallyDiscontinuousDofMap( mesh_type& M, size_type start_next_free_dof );
1015 
1021  void buildDofMap( mesh_type& mesh, size_type start_next_free_dof = 0 );
1022 
1028  void buildBoundaryDofMap( mesh_type& mesh );
1029 
1030 #if defined(FEELPP_ENABLE_MPI_MODE)
1031 
1034  void buildGhostDofMap( mesh_type& mesh );
1035 
1039  void buildGhostInterProcessDofMap( mesh_type& mesh,
1040  std::map<size_type,boost::tuple<size_type,size_type> > & mapInterProcessDof );
1041  void buildGlobalProcessToGlobalClusterDofMapContinuous( mesh_type& mesh );
1042  void buildGlobalProcessToGlobalClusterDofMapContinuousActifDof( mesh_type& mesh,
1043  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > & listToSend );
1044  void buildGlobalProcessToGlobalClusterDofMapContinuousGhostDof( mesh_type& mesh,
1045  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > const& listToSend );
1046  void buildGlobalProcessToGlobalClusterDofMapDiscontinuous();
1047 
1048  void buildGhostInterProcessDofMapInit( mesh_type& mesh,
1049  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > & listToSend );
1050  boost::tuple<bool, std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > >
1051  buildGhostInterProcessDofMapRecursive( mesh_type& mesh,
1052  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > const& listToSend,
1053  std::map<size_type,boost::tuple<size_type,size_type> > & mapInterProcessDof,
1054  std::vector< std::set<size_type > > & memoryFace );
1055 
1056  void buildDofNotPresent( std::map<size_type,boost::tuple<size_type,size_type> > const & mapInterProcessDof,
1057  std::map<size_type,boost::tuple<size_type,size_type> > & setInterProcessDofNotPresent );
1058 
1059  void buildGlobalProcessToGlobalClusterDofMap( mesh_type& mesh,
1060  std::map<size_type,boost::tuple<size_type,size_type> > const& setInterProcessDofNotPresent );
1061  void updateGhostGlobalDof( std::map<size_type,boost::tuple<size_type,size_type> > const& setInterProcessDofNotPresent );
1062 
1063 #endif
1064 
1067  dof_map_type const& mapGDof() const
1068  {
1069  return map_gdof;
1070  }
1071 
1075  dof_map_type& mapGDof()
1076  {
1077  return map_gdof;
1078  }
1079 
1084  {
1085  map_gdof.clear();
1086  }
1087 
1091  void setMapGDof( dof_map_type const& mapdof )
1092  {
1093  map_gdof = mapdof;
1094  }
1095 
1096  typename dof_marker_type::right_range_type
1097  markerToDof( boost::any const& marker )
1098  {
1099  using namespace boost::bimaps;
1100  int id = M_mesh->markerId( marker );
1101  return M_dof_marker.right.range( id <= _key, _key<id+1 );
1102  }
1103 
1104  typename dof_marker_type::right_range_type
1105  markerToDofLessThan( boost::any const& marker )
1106  {
1107  using namespace boost::bimaps;
1108  int id = M_mesh->markerId( marker );
1109  return M_dof_marker.right.range( unbounded, _key<id );
1110  }
1111  typename dof_marker_type::right_range_type
1112  markerToDofGreaterThan( boost::any const& marker )
1113  {
1114  using namespace boost::bimaps;
1115  int id = M_mesh->markerId( marker );
1116  return M_dof_marker.right.range( id<_key, unbounded );
1117  }
1118 
1119  void printDofMarker(std::string const& filename )
1120  {
1121  std::ofstream ofs( filename.c_str() );
1122  BOOST_FOREACH( auto dof, M_dof_marker.left )
1123  {
1124  ofs << dof.first << " " << dof.second << "\n";
1125  }
1126  }
1143  uint16_type l_dof,
1144  uint16_type lc,
1145  dof_type gDof,
1146  uint16_type processor,
1147  size_type& pDof,
1148  int32_type sign = 1,
1149  bool is_dof_periodic = false,
1150  size_type shift = 0,
1151  Marker1 const& marker = Marker1( 0 ) )
1152  {
1153  bool res = true;
1154  const int ncdof = is_product?nComponents:1;
1155 
1156  for ( int c = 0; c < ncdof; ++c )
1157  {
1158  gDof.template get<1>() = c;
1159  uint16_type lc_dof = fe_type::nLocalDof*c+l_dof;
1160  Feel::detail::ignore_unused_variable_warning( lc );
1161  dof_map_iterator itdof = map_gdof.find( gDof );
1162  dof_map_iterator endof = map_gdof.end();
1163  bool __inserted = false;
1164 
1165  if ( itdof == endof )
1166  {
1167  DVLOG(4) << "[dof] dof (" << gDof.template get<0>() << "," << gDof.template get<1>() << "," << gDof.template get<2>() << ") not yet inserted in map\n";
1168  boost::tie( itdof, __inserted ) = map_gdof.insert( std::make_pair( gDof, dofIndex( pDof ) ) );
1169 
1170  pDof += 1;
1171 
1172  FEELPP_ASSERT( __inserted == true )( ie )( lc_dof )
1173  ( gDof.template get<0>() )( gDof.template get<1>() )( gDof.template get<2>() )
1174  ( processor )( itdof->second ).error( "dof should have been inserted" );
1175  }
1176 
1177  else
1178  {
1179  DVLOG(4) << "[dof] dof (" << gDof.template get<0>() << ","
1180  << gDof.template get<1>()
1181  << "," << gDof.template get<2>()
1182  << ") already inserted in map with dof_id = " << itdof->second << "\n";
1183  }
1184 
1185 #if !defined( NDEBUG )
1186  DVLOG(4) << "global dof = " << itdof->second
1187  << " local dof = " << fe_type::nLocalDof*itdof->first.template get<1>() + lc_dof
1188  << " element = " << ie
1189  << " entity = " << itdof->first.template get<0>()
1190  << " component = " << itdof->first.template get<1>()
1191  << " index = " << itdof->first.template get<2>() << "\n";
1192 #endif
1193  auto eit = M_el_l2g.left.find( LocalDof(ie,lc_dof) );
1194  // make sure that no already created dof is overwritten here (may be done alsewhere)
1195  if ( eit == M_el_l2g.left.end() )
1196  {
1197  DCHECK( itdof->first == gDof ) << "very bad logical error in insertDof";
1198 
1199  DCHECK( lc_dof >= fe_type::nLocalDof*itdof->first.template get<1>() &&
1200  lc_dof < fe_type::nLocalDof*( itdof->first.template get<1>()+1 ) )
1201  << "invalid local dof index"
1202  << lc_dof << ", " << fe_type::nLocalDof*itdof->first.template get<1>();
1203 
1204 
1205  // add processor to the set of processors this dof belongs to
1206  M_dof_procset[ itdof->second+shift ].insert( processor );
1207 
1208  auto res = M_el_l2g.insert( dof_relation( LocalDof( ie, lc_dof ),
1209  Dof( itdof->second+shift, sign, is_dof_periodic, 0, 0, marker.value() ) ) );
1210  DCHECK( res.second ) << "global dof " << itdof->second+shift << " not inserted in local dof (" <<
1211  ie << "," << lc_dof << ")";
1212  M_dof_marker.insert( dof2marker( itdof->second+shift, marker.value() ) );
1213 
1214 
1215 #if !defined(NDEBUG)
1216  M_dof2elt[itdof->second+shift].push_back( boost::make_tuple( ie, lc_dof, lc, itdof->first.template get<0>() ) );
1217 #endif
1218 #if 0
1219  M_dof_view.insert( Dof( itdof->second+shift, // global index
1220  sign, // sign
1221  itdof->first.template get<0>(), // entity type
1222  false, // is on boundary ?
1223  0 // marker
1224  ) );
1225 #endif
1226 #if 0
1227  for ( index i2 = 0; i2 < nLocalDof(); ++i2 )
1228  VLOG(1) << "dof table( " << ie << ", " << lc << ")=" << M_el_l2g.left.find(LocalDof(ie,i2))->second.index() << "\n";
1229 
1230 #endif
1231  }
1232  else
1233  {
1234  size_type _dof = M_el_l2g.left.find(LocalDof(ie,lc_dof))->second.index();
1235 #if 0
1236  CHECK( M_dof_marker[_dof] == marker.value() ) << "Invalid dof marker, element id: " << ie
1237  << ", local dof id: " << lc_dof
1238  << ", global dof id: "<< _dof
1239  << ", dof marker: " << M_dof_marker[_dof]
1240  << ", marker: " << marker.value() << "\n";
1241 #endif
1242  }
1243 
1244  res = res && ( __inserted || ( ( M_el_l2g.left.find(LocalDof(ie,lc_dof)) != M_el_l2g.left.end() ) && shift ) );
1245  }
1246 
1247  return res;
1248  }
1249 
1254  {
1255  M_dof_points.clear();
1256  this->generateDofPoints(M);
1257  }
1258 
1262  std::pair<std::map<size_type,size_type>,std::map<size_type,size_type> >
1263  pointIdToDofRelation(std::string fname="") const;
1264 private:
1265 
1266  void addVertexDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1267  ref_shift_type& shifts )
1268  {
1269  addVertexDof( __elt, processor, next_free_dof, shifts, mpl::bool_<(fe_type::nDofPerVertex>0)>() );
1270  }
1271  void addVertexDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1272  ref_shift_type& /*shifts*/, mpl::bool_<false> )
1273  {}
1274  void addVertexDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1275  ref_shift_type& shifts, mpl::bool_<true> )
1276  {
1277  uint16_type local_shift;
1278  size_type global_shift;
1279  boost::tie( local_shift, global_shift ) = shifts;
1280 
1281 
1282  size_type ie = __elt.id();
1283 
1284  uint16_type lc = local_shift;
1285 
1286  for ( uint16_type i = 0; i < element_type::numVertices; ++i )
1287  {
1288  for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l, ++lc )
1289  {
1290  //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
1291  const size_type gDof = ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
1292  this->insertDof( ie, lc, i, boost::make_tuple( 0, 0, gDof ),
1293  processor, next_free_dof, 1, false, global_shift, __elt.point( i ).marker() );
1294  }
1295  }
1296 
1297  // update shifts
1298  shifts.template get<0>() = lc;
1299 
1300 #if !defined(NDEBUG)
1301  DVLOG(4) << "[Dof::updateVolumeDof(addVertexDof] vertex proc" << processor << " next_free_dof = " << next_free_dof << "\n";
1302 #endif
1303  }
1304  void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1305  ref_shift_type& shifts )
1306  {
1307  static const bool cond = fe_type::nDofPerEdge > 0;
1308  return addEdgeDof( __elt,
1309  processor,
1310  next_free_dof,
1311  shifts,
1312  mpl::int_<fe_type::nDim>(),
1313  mpl::bool_<cond>() );
1314  }
1315  void addEdgeDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1316  ref_shift_type& /*shifts*/, mpl::int_<1>, mpl::bool_<false> )
1317  {}
1318 
1319  void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1320  ref_shift_type& shifts, mpl::int_<1>, mpl::bool_<true> )
1321  {
1322  uint16_type local_shift;
1323  size_type global_shift;
1324  boost::tie( local_shift, global_shift ) = shifts;
1325 
1326  size_type ie = __elt.id();
1327  uint16_type lc = local_shift;
1328 
1329  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
1330  {
1331  const size_type gDof = is_p0_continuous? l:ie * fe_type::nDofPerEdge + l;
1332  this->insertDof( ie, lc, l, boost::make_tuple( 1, 0, gDof ), processor, next_free_dof, 1, false, global_shift, __elt.marker() );
1333  }
1334 
1335  // update shifts
1336  shifts.template get<0>() = lc;
1337 #if !defined(NDEBUG)
1338  DVLOG(4) << "[Dof::addEdgeDof(1)] element proc" << processor << " next_free_dof = " << next_free_dof << "\n";
1339 #endif
1340  }
1341  void addEdgeDof( element_type const& /*__elt*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1342  ref_shift_type& /*shifts*/, mpl::int_<2>, mpl::bool_<false> )
1343  {}
1344  void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1345  ref_shift_type& shifts, mpl::int_<2>, mpl::bool_<true> )
1346  {
1347  uint16_type local_shift;
1348  size_type global_shift;
1349  boost::tie( local_shift, global_shift ) = shifts;
1350 
1351  size_type ie = __elt.id();
1352  uint16_type lc = local_shift;
1353 
1356  for ( uint16_type i = 0; i < element_type::numEdges; ++i )
1357  {
1358  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
1359  {
1360  size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;
1361  int32_type sign = 1;
1362 
1363 
1364  if ( __elt.edgePermutation( i ).value() == edge_permutation_type::IDENTITY )
1365  {
1366  gDof += l ; // both nodal and modal case
1367  }
1368 
1369  else if ( __elt.edgePermutation( i ).value() == edge_permutation_type::REVERSE_PERMUTATION )
1370  {
1371 
1372  if ( fe_type::is_modal )
1373  {
1374  //only half of the modes (odd polynomial order) are negative.
1375  sign = ( l%2 )?( -1 ):( 1 );
1376  gDof += l;
1377  }
1378 
1379  else
1380  gDof += fe_type::nDofPerEdge - 1 - l ;
1381  }
1382 
1383  else
1384  FEELPP_ASSERT( 0 ).error ( "invalid edge permutation" );
1385 
1386  this->insertDof( ie, lc, i, boost::make_tuple( 1, 0, gDof ), processor, next_free_dof, sign, false, global_shift, __elt.edge( i ).marker() );
1387  }
1388  }
1389 
1390  // update shifts
1391  shifts.template get<0>() = lc;
1392 #if !defined(NDEBUG)
1393  DVLOG(4) << "[Dof::addEdgeDof] edge proc" << processor << " next_free_dof = " << next_free_dof << "\n";
1394 #endif
1395  }
1396 
1397  void addEdgeDof( element_type const& /*__elt*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1398  ref_shift_type& /*shifts*/, mpl::int_<3>, mpl::bool_<false> )
1399  {}
1400 
1401  void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1402  ref_shift_type& shifts, mpl::int_<3>, mpl::bool_<true> )
1403  {
1404  uint16_type local_shift;
1405  size_type global_shift;
1406  boost::tie( local_shift, global_shift ) = shifts;
1407 
1408  size_type ie = __elt.id();
1409  uint16_type lc = local_shift;
1410 
1411  for ( uint16_type i = 0; i < element_type::numEdges; ++i )
1412  {
1413  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
1414  {
1415  size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;
1416 
1417  int32_type sign = 1;
1418 
1419  if ( __elt.edgePermutation( i ).value() == edge_permutation_type::IDENTITY )
1420  {
1421  gDof += l ; // both nodal and modal case
1422  }
1423 
1424  else if ( __elt.edgePermutation( i ).value() == edge_permutation_type::REVERSE_PERMUTATION )
1425  {
1426 
1427  if ( fe_type::is_modal )
1428  {
1429  //only half of the modes (odd polynomial order) are negative.
1430  sign = ( l%2 )?( -1 ):( 1 );
1431  gDof += l;
1432  }
1433 
1434  else
1435  gDof += fe_type::nDofPerEdge - 1 - l ;
1436  }
1437 
1438  else
1439  FEELPP_ASSERT( 0 ).error ( "invalid edge permutation" );
1440 
1441  this->insertDof( ie, lc, i, boost::make_tuple( 1, 0, gDof ), processor, next_free_dof, sign, false, global_shift, __elt.edge( i ).marker() );
1442  }
1443  }
1444 
1445  // update shifts
1446  shifts.template get<0>() = lc;
1447 #if !defined(NDEBUG)
1448  DVLOG(4) << "[Dof::addEdgeDof] edge proc" << processor << " next_free_dof = " << next_free_dof << "\n";
1449 #endif
1450  }
1451 
1452 
1453  void addFaceDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1454  ref_shift_type& shifts )
1455  {
1456  return addFaceDof( __elt, processor, next_free_dof, shifts, mpl::int_<fe_type::nDim>(), mpl::bool_<(fe_type::nDofPerFace > 0)>() );
1457  }
1458  void addFaceDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1459  ref_shift_type& /*shifts*/, mpl::int_<1>, mpl::bool_<false> )
1460  {}
1461  void addFaceDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1462  ref_shift_type& /*shifts*/, mpl::int_<2>, mpl::bool_<false> )
1463  {}
1464  void addFaceDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1465  ref_shift_type& shifts, mpl::int_<2>, mpl::bool_<true> )
1466  {
1467  uint16_type local_shift;
1468  size_type global_shift;
1469  boost::tie( local_shift, global_shift ) = shifts;
1470 
1471  size_type ie = __elt.id();
1472  uint16_type lc = local_shift;
1473 
1474  for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
1475  {
1476  const size_type gDof = is_p0_continuous? l:ie * fe_type::nDofPerFace + l;
1477  this->insertDof( ie, lc, l, boost::make_tuple( 2, 0, gDof ), processor, next_free_dof, 1, false, global_shift, __elt.marker() );
1478  }
1479 
1480  // update shifts
1481  shifts.template get<0>() = lc;
1482 #if !defined(NDEBUG)
1483  DVLOG(4) << "[Dof::addFaceDof(2,true)] face proc" << processor << " next_free_dof = " << next_free_dof << "\n";
1484 #endif
1485  }
1486  void addFaceDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1487  ref_shift_type& /*shifts*/, mpl::int_<3>, mpl::bool_<false> )
1488  {}
1489  void addFaceDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1490  ref_shift_type& shifts, mpl::int_<3>, mpl::bool_<true> )
1491  {
1492  uint16_type local_shift;
1493  size_type global_shift;
1494  boost::tie( local_shift, global_shift ) = shifts;
1495 
1496  size_type ie = __elt.id();
1497 
1498  uint16_type lc = local_shift;
1499 
1500  for ( uint16_type i = 0; i < element_type::numFaces; ++i )
1501  {
1502  face_permutation_type permutation = __elt.facePermutation( i );
1503  FEELPP_ASSERT( permutation != face_permutation_type( 0 ) ).error ( "invalid face permutation" );
1504 
1505  // Polynomial order in each direction
1506  uint16_type p=1;
1507  uint16_type q=0;
1508 
1509  // MaxOrder = Order - 2
1510  int MaxOrder = int( ( 3 + std::sqrt( 1+8*fe_type::nDofPerFace ) )/2 ) - 2;
1511 
1512  for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
1513  {
1514 
1515  // TODO: orient the dof indices such
1516  // that they match properly the faces
1517  // dof of the connected faces. There
1518  // are a priori many permutations of
1519  // the dof face indices
1520  size_type gDof = __elt.face( i ).id() * fe_type::nDofPerFace;
1521  int32_type sign = 1;
1522 
1523  q=q+1;
1524 
1525  if ( q > MaxOrder )
1526  {
1527  q = 1;
1528  p = p+1;
1529  MaxOrder = MaxOrder-1;
1530  }
1531 
1532  if ( !fe_type::is_modal )
1533  {
1534  // no need of permutation is identity or only one dof on face
1535  if ( permutation == face_permutation_type( 1 ) || fe_type::nDofPerFace == 1 )
1536  gDof += l;
1537 
1538  else
1539  gDof += vector_permutation[permutation][l];
1540  }
1541 
1542  else
1543  {
1544  gDof += l;
1545 
1546  if ( permutation == face_permutation_type( 2 ) )
1547  {
1548  // Reverse sign if polynomial order in
1549  // eta_1 direction is odd
1550 
1551  if ( p%2 == 0 )
1552  sign = -1;
1553 
1554  }
1555  }
1556 
1557  this->insertDof( ie, lc, i, boost::make_tuple( 2, 0, gDof ), processor, next_free_dof, sign, false, global_shift,__elt.face( i ).marker() );
1558 
1559  }
1560  }
1561 
1562  // update shifts
1563  shifts.template get<0>() = lc;
1564 #if !defined(NDEBUG)
1565  DVLOG(4) << "[Dof::addFaceDof<3>] face proc" << processor << " next_free_dof = " << next_free_dof << "\n";
1566 #endif
1567  }
1568  void addVolumeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1569  ref_shift_type& shifts )
1570  {
1571  return addVolumeDof( __elt, processor, next_free_dof, shifts, mpl::bool_<(fe_type::nDofPerVolume>0)>() );
1572  }
1573  void addVolumeDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
1574  ref_shift_type& /*shifts*/, mpl::bool_<false> )
1575  {}
1576  void addVolumeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
1577  ref_shift_type& shifts, mpl::bool_<true> )
1578  {
1579  BOOST_STATIC_ASSERT( element_type::numVolumes );
1580  uint16_type local_shift;
1581  size_type global_shift;
1582  boost::tie( local_shift, global_shift ) = shifts;
1583 
1584  size_type ie = __elt.id();
1585  uint16_type lc = local_shift;
1586 
1587  for ( uint16_type l = 0; l < fe_type::nDofPerVolume; ++l, ++lc )
1588  {
1589  const size_type gDof = is_p0_continuous? l:ie * fe_type::nDofPerVolume + l;
1590  this->insertDof( ie, lc, l, boost::make_tuple( 3, 0, gDof ), processor, next_free_dof, 1, false, global_shift, __elt.marker() );
1591  }
1592 
1593  // update shifts
1594  shifts.template get<0>() = lc;
1595 #if !defined(NDEBUG)
1596  DVLOG(4) << "[Dof::updateVolumeDof(<2>)] element proc" << processor << " next_free_dof = " << next_free_dof << "\n";
1597 #endif
1598  }
1599 
1600  template<typename FaceIterator>
1601  void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc )
1602  {
1603  addVertexBoundaryDof( __face_it, lc, mpl::bool_<(fe_type::nDofPerVertex>0)>(), mpl::int_<nDim>() );
1604  }
1605  template<typename FaceIterator> void addVertexBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<1> ) {}
1606  template<typename FaceIterator> void addVertexBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<2> ) {}
1607  template<typename FaceIterator> void addVertexBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<3> ) {}
1608  template<typename FaceIterator>
1609  void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<1> )
1610  {
1611  BOOST_STATIC_ASSERT( face_type::numVertices );
1612 #if !defined(FEELPP_ENABLE_MPI_MODE)
1613  // id of the element adjacent to the face
1614  // \warning NEED TO INVESTIGATE THIS
1615  size_type iElAd = __face_it->ad_first();
1616  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1617 
1618  // local id of the face in its adjacent element
1619  uint16_type iFaEl = __face_it->pos_first();
1620  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1621 #else // MPI
1622  uint16_type iFaEl;
1623  size_type iElAd;
1624 
1625  if ( __face_it->processId() == __face_it->proc_first() )
1626  {
1627  iElAd = __face_it->ad_first();
1628  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1629 
1630  // local id of the face in its adjacent element
1631  iFaEl = __face_it->pos_first();
1632  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1633  }
1634 
1635  else
1636  {
1637  iElAd = __face_it->ad_second();
1638  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1639 
1640  // local id of the face in its adjacent element
1641  iFaEl = __face_it->pos_second();
1642  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1643  }
1644 
1645 #endif
1646  // Loop number of Dof per vertex
1647  const int ncdof = is_product?nComponents:1;
1648 
1649  for ( int c = 0; c < ncdof; ++c )
1650  {
1651  for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
1652  {
1653  auto temp= this->localToGlobal( iElAd,
1654  iFaEl * fe_type::nDofPerVertex + l,
1655  c );
1656  M_face_l2g[ __face_it->id()][ lc++ ] = boost::make_tuple( boost::get<0>( temp ),boost::get<1>( temp ),boost::get<2>( temp ),
1657  iFaEl * fe_type::nDofPerVertex + l );
1658  }
1659  }
1660  }
1661  template<typename FaceIterator>
1662  void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<2> )
1663  {
1664  BOOST_STATIC_ASSERT( face_type::numVertices );
1665 #if !defined(FEELPP_ENABLE_MPI_MODE)
1666  // id of the element adjacent to the face
1667  // \warning NEED TO INVESTIGATE THIS
1668  size_type iElAd = __face_it->ad_first();
1669  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1670 
1671  // local id of the face in its adjacent element
1672  uint16_type iFaEl = __face_it->pos_first();
1673  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1674 #else // MPI
1675  uint16_type iFaEl;
1676  size_type iElAd;
1677 
1678  if ( __face_it->processId() == __face_it->proc_first() )
1679  {
1680  iElAd = __face_it->ad_first();
1681  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1682 
1683  // local id of the face in its adjacent element
1684  iFaEl = __face_it->pos_first();
1685  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1686  }
1687 
1688  else
1689  {
1690  iElAd = __face_it->ad_second();
1691  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1692 
1693  // local id of the face in its adjacent element
1694  iFaEl = __face_it->pos_second();
1695  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1696  }
1697 
1698 #endif
1699 
1700  size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
1701  face_type::numEdges * fe_type::nDofPerEdge +
1702  face_type::numFaces * fe_type::nDofPerFace );
1703 
1704 
1705  //M_dof2elt[gDof].push_back( boost::make_tuple( iElAd, lc-1, 48, 0 ) );
1706  // loop on face vertices
1707  const int ncdof = is_product?nComponents:1;
1708 
1709  for ( int c = 0; c < ncdof; ++c )
1710  {
1711  uint16_type lcc=c*ndofF;
1712 
1713  for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
1714  {
1715  // local vertex number (in element)
1716  uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
1717 
1718  FEELPP_ASSERT( iVeEl != invalid_uint16_type_value ).error( "invalid local dof" );
1719 
1720  // Loop number of Dof per vertex
1721  for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
1722  {
1723  auto temp = this->localToGlobal( iElAd,
1724  iVeEl * fe_type::nDofPerVertex + l,
1725  c );
1726  M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>( temp ),boost::get<1>( temp ),boost::get<2>( temp ),
1727  iVeEl * fe_type::nDofPerVertex + l );
1728  }
1729  }
1730  }
1731  }
1732  template<typename FaceIterator>
1733  void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<3> )
1734  {
1735  addVertexBoundaryDof( __face_it, lc, mpl::bool_<true>(), mpl::int_<2>() );
1736  }
1737  template<typename FaceIterator>
1738  void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc )
1739  {
1740  static const bool cond = fe_type::nDofPerEdge*face_type::numEdges > 0;
1741  addEdgeBoundaryDof( __face_it, lc, mpl::bool_<cond>(), mpl::int_<nDim>() );
1742  }
1743  template<typename FaceIterator>
1744  void addEdgeBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<1> ) {}
1745  template<typename FaceIterator>
1746  void addEdgeBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<2> ) {}
1747  template<typename FaceIterator>
1748  void addEdgeBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<3> ) {}
1749  template<typename FaceIterator>
1750  void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<2> )
1751  {
1752 #if !defined(FEELPP_ENABLE_MPI_MODE)
1753  // id of the element adjacent to the face
1754  // \warning NEED TO INVESTIGATE THIS
1755  size_type iElAd = __face_it->ad_first();
1756  FEELPP_ASSERT( iElAd != invalid_size_type_value )
1757  ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1758 
1759  // local id of the face in its adjacent element
1760  uint16_type iFaEl = __face_it->pos_first();
1761 #else // MPI
1762  uint16_type iFaEl;
1763  size_type iElAd;
1764 
1765  if ( __face_it->processId() == __face_it->proc_first() )
1766  {
1767  iElAd = __face_it->ad_first();
1768  FEELPP_ASSERT( iElAd != invalid_size_type_value )
1769  ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1770 
1771  // local id of the face in its adjacent element
1772  iFaEl = __face_it->pos_first();
1773  }
1774 
1775  else
1776  {
1777  iElAd = __face_it->ad_second();
1778  FEELPP_ASSERT( iElAd != invalid_size_type_value )
1779  ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1780 
1781  // local id of the face in its adjacent element
1782  iFaEl = __face_it->pos_second();
1783  }
1784 
1785 #endif
1786 
1787 
1788  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1789 #if !defined(NDEBUG)
1790  DVLOG(4) << " local face id : " << iFaEl << "\n";
1791 #endif
1792  size_type nVerticesF = face_type::numVertices * fe_type::nDofPerVertex;
1793  size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
1794  face_type::numEdges * fe_type::nDofPerEdge +
1795  face_type::numFaces * fe_type::nDofPerFace );
1796 
1797 
1798  const int ncdof = is_product?nComponents:1;
1799 
1800  for ( int c = 0; c < ncdof; ++c )
1801  {
1802  uint16_type lcc=nVerticesF+c*ndofF;
1803 
1804  // Loop number of Dof per edge
1805  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
1806  {
1807  auto temp = this->localToGlobal( iElAd,
1808  element_type::numVertices*fe_type::nDofPerVertex +
1809  iFaEl * fe_type::nDofPerEdge + l,
1810  c );
1811  M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>( temp ),boost::get<1>( temp ),boost::get<2>( temp ),
1812  element_type::numVertices*fe_type::nDofPerVertex +
1813  iFaEl * fe_type::nDofPerEdge + l );
1814  }
1815  }
1816  }
1817  template<typename FaceIterator>
1818  void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<3> )
1819  {
1820  //BOOST_STATIC_ASSERT( face_type::numEdges );
1821 #if !defined(FEELPP_ENABLE_MPI_MODE)
1822  // id of the element adjacent to the face
1823  // \warning NEED TO INVESTIGATE THIS
1824  size_type iElAd = __face_it->ad_first();
1825  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1826 
1827  // local id of the face in its adjacent element
1828  uint16_type iFaEl = __face_it->pos_first();
1829  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1830 #else // MPI
1831  uint16_type iFaEl;
1832  size_type iElAd;
1833 
1834  if ( __face_it->processId() == __face_it->proc_first() )
1835  {
1836  iElAd = __face_it->ad_first();
1837  FEELPP_ASSERT( iElAd != invalid_size_type_value )
1838  ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1839 
1840  // local id of the face in its adjacent element
1841  iFaEl = __face_it->pos_first();
1842  }
1843 
1844  else
1845  {
1846  iElAd = __face_it->ad_second();
1847  FEELPP_ASSERT( iElAd != invalid_size_type_value )
1848  ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1849 
1850  // local id of the face in its adjacent element
1851  iFaEl = __face_it->pos_second();
1852  }
1853 
1854 #endif
1855 
1856 #if !defined(NDEBUG)
1857  DVLOG(4) << " local face id : " << iFaEl << "\n";
1858 #endif
1859  size_type nVerticesF = face_type::numVertices * fe_type::nDofPerVertex;
1860  size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
1861  face_type::numEdges * fe_type::nDofPerEdge +
1862  face_type::numFaces * fe_type::nDofPerFace );
1863 
1864  const int ncdof = is_product?nComponents:1;
1865 
1866  for ( int c = 0; c < ncdof; ++c )
1867  {
1868  uint16_type lcc=nVerticesF+c*ndofF;
1869 
1870  // loop on face vertices
1871  for ( uint16_type iEdFa = 0; iEdFa < face_type::numEdges; ++iEdFa )
1872  {
1873  // local edge number (in element)
1874  uint16_type iEdEl = element_type::fToE( iFaEl, iEdFa );
1875 
1876  FEELPP_ASSERT( iEdEl != invalid_uint16_type_value ).error( "invalid local dof" );
1877 
1878  // Loop number of Dof per edge
1879  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
1880  {
1881  auto temp = this->localToGlobal( iElAd,
1882  element_type::numVertices*fe_type::nDofPerVertex +
1883  iEdEl * fe_type::nDofPerEdge + l,
1884  c );
1885  M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>( temp ),boost::get<1>( temp ),boost::get<2>( temp ),
1886  element_type::numVertices*fe_type::nDofPerVertex +
1887  iEdEl * fe_type::nDofPerEdge + l );
1888  }
1889  }
1890  }
1891  }
1892 
1893 
1894  template<typename FaceIterator>
1895  void addFaceBoundaryDof( FaceIterator __face_it, uint16_type& lc )
1896  {
1897  addFaceBoundaryDof( __face_it, lc, mpl::bool_<(face_type::numFaces*fe_type::nDofPerFace > 0)>() );
1898  }
1899  template<typename FaceIterator>
1900  void addFaceBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false> )
1901  {
1902  }
1903  template<typename FaceIterator>
1904  void addFaceBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true> )
1905  {
1906 #if !defined(FEELPP_ENABLE_MPI_MODE)
1907  // id of the element adjacent to the face
1908  // \warning NEED TO INVESTIGATE THIS
1909  size_type iElAd = __face_it->ad_first();
1910  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1911 
1912  // local id of the face in its adjacent element
1913  uint16_type iFaEl = __face_it->pos_first();
1914  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
1915 #else // MPI
1916  uint16_type iFaEl;
1917  size_type iElAd;
1918 
1919  if ( __face_it->processId() == __face_it->proc_first() )
1920  {
1921  iElAd = __face_it->ad_first();
1922  FEELPP_ASSERT( iElAd != invalid_size_type_value )
1923  ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1924 
1925  // local id of the face in its adjacent element
1926  iFaEl = __face_it->pos_first();
1927  }
1928 
1929  else
1930  {
1931  iElAd = __face_it->ad_second();
1932  FEELPP_ASSERT( iElAd != invalid_size_type_value )
1933  ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
1934 
1935  // local id of the face in its adjacent element
1936  iFaEl = __face_it->pos_second();
1937  }
1938 
1939 #endif
1940 
1941 #if !defined(NDEBUG)
1942  DVLOG(4) << " local face id : " << iFaEl << "\n";
1943 #endif
1944  size_type nVerticesAndEdgeF = ( face_type::numVertices * fe_type::nDofPerVertex +
1945  face_type::numEdges * fe_type::nDofPerEdge );
1946  size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
1947  face_type::numEdges * fe_type::nDofPerEdge +
1948  face_type::numFaces * fe_type::nDofPerFace );
1949 
1950  const int ncdof = is_product?nComponents:1;
1951 
1952  for ( int c = 0; c < ncdof; ++c )
1953  {
1954  uint16_type lcc=nVerticesAndEdgeF+c*ndofF;
1955 
1956  // Loop on number of Dof per face
1957  for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l )
1958  {
1959  auto temp = this->localToGlobal( iElAd,
1960  element_type::numVertices*fe_type::nDofPerVertex +
1961  element_type::numEdges*fe_type::nDofPerEdge +
1962  iFaEl * fe_type::nDofPerFace + l,
1963  c );
1964  M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>( temp ),boost::get<1>( temp ),boost::get<2>( temp ),
1965  element_type::numVertices*fe_type::nDofPerVertex +
1966  element_type::numEdges*fe_type::nDofPerEdge +
1967  iFaEl * fe_type::nDofPerFace + l );
1968  }
1969  }
1970  }
1971 
1972  void addSubstructuringDofMap( mesh_type const& M, size_type next_free_dof );
1973  void addSubstructuringDofVertex( mesh_type const& M, size_type next_free_dof );
1974  void addSubstructuringDofEdge( mesh_type const& M, size_type next_free_dof, mpl::int_<1> );
1975  void addSubstructuringDofEdge( mesh_type const& M, size_type next_free_dof, mpl::int_<2> );
1976  void addSubstructuringDofEdge( mesh_type const& M, size_type next_free_dof, mpl::int_<3> );
1977  void addSubstructuringDofFace( mesh_type const& M, size_type next_free_dof, mpl::int_<1> );
1978  void addSubstructuringDofFace( mesh_type const& M, size_type next_free_dof, mpl::int_<2> );
1979  void addSubstructuringDofFace( mesh_type const& M, size_type next_free_dof, mpl::int_<3> );
1980 
1986  template<int N>
1987  void checkDofEntity ( mesh_type& M )
1988  {
1989  Feel::detail::ignore_unused_variable_warning( M );
1990 
1991  if ( !is_scalar )
1992  return;
1993 
1994 #if !defined(NDEBUG)
1995 
1996  using namespace Feel;
1997 
1998  gm_ptrtype M_gm_ptr = M.gm();
1999 
2000  fe_type M_basis;
2001 
2002  //value_type tol = value_type(100.0)*type_traits<value_type>::epsilon();
2003  value_type tol = value_type( 10.0 )*type_traits<double>::epsilon();
2004 
2005  bool global_signs_good = 1;
2006 
2007  std::vector<size_type> bad_dof;
2008 
2009  for ( uint16_type gDof = 0; gDof < nDof(); ++gDof )
2010  {
2011  uint16_type _numEntities = M_dof2elt[gDof].size();
2012  uint16_type _ent = M_dof2elt[gDof].begin()->template get<3>();
2013 
2014  if ( _numEntities > 1 && _ent == mpl::int_<N>() )
2015  {
2016  bool signs_good = 1;
2017 
2018  std::vector< ublas::vector<value_type> > basis_eval;
2019 
2020  std::vector< points_type > real_coordinates;
2021 
2022  ldof_const_iterator __ldofit = M_dof2elt[gDof].begin();
2023  ldof_const_iterator __ldofen = M_dof2elt[gDof].end();
2024 
2025  while ( __ldofit != __ldofen )
2026  {
2027  size_type entity_element_id = __ldofit->template get<0>();
2028  uint16_type entity_local_dof_id = __ldofit->template get<1>();
2029  uint16_type entity_local_id = __ldofit->template get<2>();
2030 
2031  PointSetMapped<element_type, convex_type, nOrder> test( M.element( entity_element_id ) );
2032 
2033  points_type Gt = test.pointsBySubEntity( N, entity_local_id );
2034 
2035  real_coordinates.push_back( test.pointsBySubEntity( N, entity_local_id, 0, 1 ) );
2036 
2037  int sign = boost::get<1>( localToGlobal( entity_element_id, entity_local_dof_id ) );
2038 
2039  basis_eval.push_back( value_type( sign )*ublas::row( M_basis.evaluate( Gt ), entity_local_dof_id ) );
2040 
2041  ++__ldofit;
2042  }
2043 
2044  for ( uint16_type i=1; i < _numEntities; i++ )
2045  {
2046 
2047  FEELPP_ASSERT( ublas::norm_inf( real_coordinates[i] - real_coordinates[0] ) < tol )
2048  ( gDof )
2049  ( real_coordinates[0] )
2050  ( real_coordinates[i] ).error( "Reference points aren't being mapped to the same real one's" );
2051 
2052  if ( ublas::norm_inf( basis_eval[i] - basis_eval[0] ) > tol )
2053  {
2054  signs_good = 0;
2055  global_signs_good = 0;
2056  }
2057  }
2058 
2059  basis_eval.empty();
2060  real_coordinates.empty();
2061 
2062  if ( signs_good == 0 )
2063  bad_dof.push_back( gDof );
2064  }
2065  }
2066 
2067  if ( !bad_dof.empty() )
2068  {
2069  for ( uint16_type i = 0; i < bad_dof.size(); ++i )
2070  Warning() << bad_dof[i] << "\n";
2071 
2072  if ( mpl::int_<N>() == 1 )
2073  Warning() << "Edges: ";
2074 
2075  else
2076  Warning() << "Faces: ";
2077 
2078  Warning() << "Bad dof signs. \n";
2079  }
2080 
2081 #endif
2082  }
2083 
2084  void checkDofContinuity( mesh_type& /*mesh*/, mpl::int_<1> ) {}
2085 
2086  void checkDofContinuity( mesh_type& mesh, mpl::int_<2> )
2087  {
2088  checkDofEntity<1>( mesh );
2089  }
2090 
2091  void checkDofContinuity( mesh_type& mesh, mpl::int_<3> )
2092  {
2093  checkDofContinuity( mesh, mpl::int_<2>() );
2094  checkDofEntity<2>( mesh );
2095  }
2096 
2097  void generateFacePermutations ( mesh_type& /*mesh*/, mpl::bool_<false> ) {}
2098 
2099  void generateFacePermutations ( mesh_type& mesh, mpl::bool_<true> )
2100  {
2101  element_type _elt = *mesh.beginElement();
2102  PointSetMapped<element_type, convex_type, nOrder> pts( _elt );
2103 
2104  for ( uint16_type i = 2; i < face_permutation_type::N_PERMUTATIONS; i++ )
2105  vector_permutation[face_permutation_type( i )] = pts.getVectorPermutation( face_permutation_type( i ) );
2106  }
2107  void generateDofPoints( mesh_type& M );
2108  void generatePeriodicDofPoints( mesh_type& M, periodic_element_list_type const& periodic_elements, dof_points_type& periodic_dof_points );
2109 
2110 
2111 
2112 private:
2113 
2114  mesh_type* M_mesh;
2115  fe_ptrtype M_fe;
2116 
2117  reference_convex_type M_convex_ref;
2118 
2119  size_type M_n_el;
2120  std::vector<bool> M_elt_done;
2121 
2122  uint16_type M_n_dof_per_face_on_bdy;
2123  uint16_type M_n_dof_per_face;
2124 
2125  dof_table M_el_l2g;
2126  Container_fromface M_face_l2g;
2127 
2128  mutable local_dof_set_type M_local_dof_set;
2129  dof_element_type M_dof2elt;
2130  dof_marker_type M_dof_marker;
2131 
2132  dof_map_type map_gdof;
2133 
2134  std::map<face_permutation_type, permutation_vector_type> vector_permutation;
2135 
2139  std::map<size_type, std::pair<size_type, size_type> > M_local2global;
2140 
2141 
2142  face_sign_info_type M_face_sign;
2143 
2147  dof_procset_type M_dof_procset;
2148 
2152  dof_points_type M_dof_points;
2153 
2154  std::vector<Dof> M_dof_indices;
2155 
2156  periodicity_type M_periodicity;
2158  periodic_element_list_type periodic_elements;
2159 
2161  //dof_container_type M_dof_view;
2162 
2163  typedef typename std::vector<localglobal_indices_type,Eigen::aligned_allocator<localglobal_indices_type> > vector_indices_type;
2164 
2165 
2166  vector_indices_type M_locglob_indices;
2167  vector_indices_type M_locglob_signs;
2168 
2169  // multi process
2170  vector_indices_type M_locglobOnCluster_indices;
2171  vector_indices_type M_locglobOnCluster_signs;
2172 
2173 public:
2174  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
2175 };
2176 
2177 template<typename MeshType, typename FEType, typename PeriodicityType>
2179 
2180 template<typename MeshType, typename FEType, typename PeriodicityType>
2181 DofTable<MeshType, FEType, PeriodicityType>::DofTable( mesh_type& mesh, fe_ptrtype const& _fe, periodicity_type const& periodicity, WorldComm const& _worldComm )
2182  :
2183  super( _worldComm ),
2184  M_fe( _fe ),
2185  M_n_el( invalid_size_type_value ),
2186  M_n_dof_per_face_on_bdy( invalid_uint16_type_value ),
2187  M_n_dof_per_face( invalid_uint16_type_value ),
2188  M_el_l2g(),
2189  M_face_l2g(),
2190  M_local_dof_set( 0, nLocalDof() ),
2191  map_gdof(),
2192  M_face_sign(),
2193  M_dof_indices(),
2194  M_periodicity( periodicity )
2195 {
2196  VLOG(2) << "[dof] is_periodic = " << is_periodic << "\n";
2197  size_type start_next_free_dof = 0;
2198 
2199  if ( is_periodic )
2200  start_next_free_dof = buildPeriodicDofMap( mesh );
2201 
2202  buildDofMap( mesh, start_next_free_dof );
2203  buildBoundaryDofMap( mesh );
2204 }
2205 
2206 template<typename MeshType, typename FEType, typename PeriodicityType>
2207 DofTable<MeshType, FEType, PeriodicityType>::DofTable( fe_ptrtype const& _fe, periodicity_type const& periodicity, WorldComm const& _worldComm )
2208  :
2209  super( _worldComm ),
2210  M_fe( _fe ),
2211  M_n_el( 0 ),
2212  M_n_dof_per_face_on_bdy( invalid_uint16_type_value ),
2213  M_n_dof_per_face( invalid_uint16_type_value ),
2214  M_el_l2g(),
2215  M_face_l2g(),
2216  M_local_dof_set( 0, nLocalDof() ),
2217  map_gdof(),
2218  M_face_sign(),
2219  M_dof_indices(),
2220  M_periodicity( periodicity )
2221 {
2222 }
2223 
2224 template<typename MeshType, typename FEType, typename PeriodicityType>
2226  :
2227  super( dof2 ),
2228  M_fe( dof2.M_fe ),
2229  M_n_el( dof2.M_n_el ),
2230  M_n_dof_per_face_on_bdy( dof2.M_n_dof_per_face_on_bdy ),
2231  M_n_dof_per_face( dof2.M_n_dof_per_face ),
2232  M_el_l2g( dof2.M_el_l2g ),
2233  M_face_l2g( dof2.M_face_l2g ),
2234  M_local_dof_set( dof2.M_local_dof_set),
2235  map_gdof( dof2.map_gdof ),
2236  M_face_sign( dof2.M_face_sign ),
2237  M_dof_indices( dof2.M_dof_indices ),
2238  M_periodicity( dof2.M_periodicity )
2239 {
2240 }
2241 
2242 template<typename MeshType, typename FEType, typename PeriodicityType>
2243 void
2245 {
2246  LOG(INFO) << " Degree of Freedom (DofTable) Object" << "\n";
2247  //if ( verbose )
2248  {
2249  LOG(INFO) << "* nDof = " << this->nLocalDof() << "\n";
2250  LOG(INFO) << "************************************************************" << "\n";
2251  LOG(INFO) << " Local to Global DOF table" << "\n";
2252  LOG(INFO) << "************************************************************" << "\n";
2253  LOG(INFO) << "Element Id Loc. N. Global N. Sign# Element Id Loc. N. Global N. Sign" << "\n";
2254 
2255  for ( size_type i = 0; i < M_n_el; ++i )
2256  {
2257 
2258  for ( size_type j = 0; j < nDofPerElement; ++j )
2259  {
2260 
2261  LOG(INFO)<< "elt id " << i << " : "
2262  << "(local/global/sign dof : " << j << " : "
2263  << boost::get<0>( localToGlobal( i , j ) ) << " : "
2264  << boost::get<1>( localToGlobal( i , j ) ) << "\n";
2265  }
2266 
2267  }
2268 
2269  LOG(INFO) << "\n";
2270 
2271  LOG(INFO) << "************************************************************" << "\n";
2272  LOG(INFO) << " Boundary Local to Global DOF table" << "\n";
2273  LOG(INFO) << "************************************************************" << "\n";
2274 
2275  auto it = M_face_l2g.begin();
2276  auto en = M_face_l2g.end();
2277 
2278  for ( size_type f = 0; it!=en; ++it,++f )
2279  {
2280  std::ostringstream ostr;
2281  ostr << "face id " << it->first << " : ";
2282 
2283  auto it2 = it->second.begin();
2284  auto en2 = it->second.end();
2285 
2286  for ( size_type l = 0; it2!=en2; ++it2,++l )
2287  {
2288  ostr << "(local/global/sign dof : " << it2->first << " : "
2289  << boost::get<0>( it2->second )<< " : "
2290  << boost::get<1>( it2->second ) << "\n";
2291  }
2292 
2293  LOG(INFO) << ostr.str() << "\n";
2294  }
2295  }
2296 
2297 }
2298 
2299 template<typename MeshType, typename FEType, typename PeriodicityType>
2300 void
2302  face_type const& __face,
2303  size_type& next_free_dof,
2304  std::map<size_type,periodic_dof_map_type>& periodic_dof,
2305  size_type tag,
2306  mpl::bool_<true> )
2307 {
2308  // store the element and local dof id for further
2309  // reference when inserting the associated global dof
2310  // id of the element adjacent to the face
2311 
2312  size_type iElAd = __face.ad_first();
2313  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[periodic]invalid face/element in face" );
2314  Feel::detail::ignore_unused_variable_warning( iElAd );
2315 
2316  // local id of the face in its adjacent element
2317  uint16_type iFaEl = __face.pos_first();
2318  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
2319 
2320  // loop on face vertices
2321  for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
2322  {
2323  // local vertex number (in element)
2324  uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
2325  Feel::detail::ignore_unused_variable_warning( iVeEl );
2326 
2327  FEELPP_ASSERT( iVeEl != invalid_uint16_type_value ).error( "invalid local dof" );
2328 
2329  // Loop number of Dof per vertex
2330  for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
2331  {
2332  uint16_type lid = iVeEl * fe_type::nDofPerVertex + l;
2333  //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
2334  const size_type gDof = ( __elt.point( iVeEl ).id() ) * fe_type::nDofPerVertex + l;
2335 
2336  VLOG(2) << "add vertex periodic doc " << next_free_dof << " in element " << __elt.id() << " lid = " << lid << "\n";
2337  size_type dof_id = next_free_dof;
2338  // next_free_dof might be incremented if a new dof is created
2339  bool inserted = this->insertDof( __elt.id(), lid, iVeEl, boost::make_tuple( 0, 0, gDof ), 0, next_free_dof, 1, true, 0, __elt.point( iVeEl ).marker() );
2340  VLOG(2) << "vertex periodic dof inserted : " << inserted << "\n";
2341 
2342  const int ncdof = is_product?nComponents:1;
2343  for ( int c1 = 0; c1 < ncdof; ++c1 )
2344  {
2345  // add the pair (elt, lid) to the map associated
2346  // with dof_id, one dof can be shared by several
2347  // elements
2348  dof_id = localToGlobal( __elt.id(), lid, c1 ).index();
2349  periodic_dof[tag].insert( std::make_pair( dof_id, boost::make_tuple( __elt.id(), lid, c1, gDof, 0 ) ) );
2350 
2351  VLOG(2) << "added vertex periodic dof " << __elt.id() << ", " << lid << ", " << boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) ) << "\n";
2352  }
2353  }
2354 
2355  }
2356 
2357 }
2358 
2359 template<typename MeshType, typename FEType, typename PeriodicityType>
2360 void
2362  face_type const& __face,
2363  size_type& next_free_dof,
2364  std::map<size_type,periodic_dof_map_type>& periodic_dof,
2365  size_type tag,
2366  mpl::bool_<true>,
2367  mpl::int_<2> )
2368 {
2369 #if 0
2370  // id of the element adjacent to the face
2371  // \warning NEED TO INVESTIGATE THIS
2372  size_type iElAd = __face.ad_first();
2373  FEELPP_ASSERT( iElAd != invalid_size_type_value )
2374  ( __face.id() ).error( "[DofTable::buildBoundaryDof] invalid face/element in face" );
2375 #endif // 0
2376 
2377  // local id of the face in its adjacent element
2378  uint16_type iFaEl = __face.pos_first();
2379  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
2380 #if !defined(NDEBUG)
2381  DVLOG(4) << " local face id : " << iFaEl << "\n";
2382 #endif
2383 
2384  // Loop number of DofTable per edge
2385  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
2386  {
2387  uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
2388  //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
2389  const size_type gDof = ( __elt.edge( iFaEl ).id() ) * fe_type::nDofPerEdge + l;
2390 
2391  DVLOG(4) << "add edge periodic dof " << next_free_dof << " in element " << __elt.id() << " lid = " << lid << "\n";
2392  size_type dof_id = next_free_dof;
2393  // next_free_dof might be incremented if a new dof is created
2394  bool inserted = this->insertDof( __elt.id(), lid, iFaEl, boost::make_tuple( 1, 0, gDof ), 0, next_free_dof, 1, true, 0, __elt.edge( iFaEl ).marker() );
2395  DVLOG(4) << "edge periodic dof inserted (1 or 0) : " << inserted << "\n";
2396 
2397  const int ncdof = is_product?nComponents:1;
2398  for ( int c1 = 0; c1 < ncdof; ++c1 )
2399  {
2400  // add the pair (elt, lid) to the map associated
2401  // with dof_id, one dof can be shared by several
2402  // elements
2403  dof_id = boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) );
2404  periodic_dof[tag].insert( std::make_pair( dof_id, boost::make_tuple( __elt.id(), lid, c1, gDof, 1 ) ) );
2405 
2406  DVLOG(4) << "added edge periodic dof " << __elt.id() << ", " << lid << ", " << boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) ) << "\n";
2407  }
2408 
2409  }
2410 }
2411 template<typename MeshType, typename FEType, typename PeriodicityType>
2412 void
2414  face_type const& __face,
2415  size_type& next_free_dof,
2416  std::map<size_type,periodic_dof_map_type>& periodic_dof,
2417  size_type tag,
2418  mpl::bool_<true>,
2419  mpl::int_<3> )
2420 {
2421  //BOOST_STATIC_ASSERT( face_type::numEdges );
2422 
2423  // id of the element adjacent to the face
2424  // \warning NEED TO INVESTIGATE THIS
2425  size_type iElAd = __face.ad_first();
2426  Feel::detail::ignore_unused_variable_warning( iElAd );
2427  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
2428 
2429  // local id of the face in its adjacent element
2430  uint16_type iFaEl = __face.pos_first();
2431  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
2432 #if !defined(NDEBUG)
2433  DVLOG(4) << " local face id : " << iFaEl << "\n";
2434 #endif
2435 
2436  // loop on face vertices
2437  for ( uint16_type iEdFa = 0; iEdFa < face_type::numEdges; ++iEdFa )
2438  {
2439  // local edge number (in element)
2440  uint16_type iEdEl = element_type::fToE( iFaEl, iEdFa );
2441  Feel::detail::ignore_unused_variable_warning( iEdEl );
2442  FEELPP_ASSERT( iEdEl != invalid_uint16_type_value ).error( "invalid local dof" );
2443 
2444  // Loop number of Dof per edge
2445  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
2446  {
2447 
2448  uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
2449  //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
2450  size_type gDof = __elt.edge( l ).id() * fe_type::nDofPerEdge;
2451 
2452  if ( __elt.edgePermutation( l ).value() == edge_permutation_type::IDENTITY )
2453  {
2454  gDof += l ; // both nodal and modal case
2455  }
2456 
2457  else if ( __elt.edgePermutation( l ).value() == edge_permutation_type::REVERSE_PERMUTATION )
2458  {
2459  gDof += fe_type::nDofPerEdge - 1 - l ;
2460  }
2461 
2462  DVLOG(4) << "add periodic doc " << next_free_dof << " in element " << __elt.id() << " lid = " << lid << "\n";
2463  size_type dof_id = next_free_dof;
2464  // next_free_dof might be incremented if a new dof is created
2465  bool inserted = this->insertDof( __elt.id(), lid, l, boost::make_tuple( 1, 0, gDof ), 0, next_free_dof, 1, true, 0, __elt.edge( l ).marker() );
2466  DVLOG(4) << "periodic dof inserted : " << inserted << "\n";
2467 
2468  const int ncdof = is_product?nComponents:1;
2469  for ( int c1 = 0; c1 < ncdof; ++c1 )
2470  {
2471  // add the pair (elt, lid) to the map associated
2472  // with dof_id, one dof can be shared by several
2473  // elements
2474  dof_id = boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) );
2475  periodic_dof[tag].insert( std::make_pair( dof_id, boost::make_tuple( __elt.id(), lid, c1, gDof, 1 ) ) );
2476 
2477  DVLOG(4) << "added " << __elt.id() << ", " << lid << ", " << boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) ) << "\n";
2478  }
2479  }
2480  }
2481 
2482 }
2483 template<typename MeshType, typename FEType, typename PeriodicityType>
2484 void
2486  face_type const& __face,
2487  size_type& next_free_dof,
2488  std::map<size_type,periodic_dof_map_type>& periodic_dof,
2489  size_type tag,
2490  mpl::bool_<true> )
2491 {
2492 #if 0
2493  // id of the element adjacent to the face
2494  // \warning NEED TO INVESTIGATE THIS
2495  size_type iElAd = __face.ad_first();
2496  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );
2497 
2498  // local id of the face in its adjacent element
2499  uint16_type iFaEl = __face.pos_first();
2500  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
2501 #if !defined(NDEBUG)
2502  DVLOG(2) << " local face id : " << iFaEl << "\n";
2503 #endif
2504 
2505  // Loop on number of Dof per face
2506  for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l )
2507  {
2508  auto temp = this->localToGlobal( iElAd,
2509  element_type::numVertices*fe_type::nDofPerVertex +
2510  element_type::numEdges*fe_type::nDofPerEdge +
2511  iFaEl * fe_type::nDofPerFace + l,
2512  c );
2513  M_face_l2g[ __face_it->id()][ lc++ ] = boost::make_tuple( boost::get<0>( temp ),boost::get<1>( temp ),boost::get<2>( temp ),
2514  element_type::numVertices*fe_type::nDofPerVertex +
2515  element_type::numEdges*fe_type::nDofPerEdge +
2516  iFaEl * fe_type::nDofPerFace + l );
2517  }
2518 
2519 #endif // 0
2520 }
2521 template<typename MeshType, typename FEType, typename PeriodicityType>
2522 void
2524 {
2525  M_n_el = M.numElements();
2526 
2527  size_type nldof =
2528  fe_type::nDofPerVolume * element_type::numVolumes +
2529  fe_type::nDofPerFace * element_type::numGeometricFaces +
2530  fe_type::nDofPerEdge * element_type::numEdges +
2531  fe_type::nDofPerVertex * element_type::numVertices;
2532 
2533  VLOG(2) << "==============================\n";
2534  VLOG(2) << "[initDofMap]\n";
2535  VLOG(2) << "nldof = " << int( nldof ) << "\n";
2536  VLOG(2) << "fe_type::nLocalDof = " << int( fe_type::nLocalDof ) << "\n";
2537  VLOG(2) << "fe_type::nDofPerVolume = " << int( fe_type::nDofPerVolume ) << "\n";
2538  VLOG(2) << "fe_type::nDofPerFace = " << int( fe_type::nDofPerFace ) << "\n";
2539  VLOG(2) << "fe_type::nDofPerEdge = " << int( fe_type::nDofPerEdge ) << "\n";
2540  VLOG(2) << "fe_type::nDofPerVertex = " << int( fe_type::nDofPerVertex ) << "\n";
2541  VLOG(2) << "element_type::numVolumes= " << int( element_type::numVolumes ) << "\n";
2542  VLOG(2) << "element_type::numFaces= " << int( element_type::numFaces ) << "\n";
2543  VLOG(2) << "element_type::numEdges= " << int( element_type::numEdges ) << "\n";
2544  VLOG(2) << "element_type::numVertices= " << int( element_type::numVertices ) << "\n";
2545  VLOG(2) << "==============================\n";
2546 
2547  FEELPP_ASSERT( nldof == fe_type::nLocalDof )
2548  ( nldof )
2549  ( fe_type::nLocalDof ).error( "Something wrong in FE specification" ) ;
2550 
2551  // initialize the local to global map and fill it with invalid
2552  // values that will allow to check whether we have a new dof or
2553  // not when building the table
2554  const size_type nV = M.numElements();
2555  int ntldof = is_product?nComponents*nldof:nldof;//this->getIndicesSize();
2556  M_locglob_indices.resize( nV );
2557  M_locglob_signs.resize( nV );
2558 #if defined(FEELPP_ENABLE_MPI_MODE)
2559  M_locglobOnCluster_indices.resize( nV );
2560  M_locglobOnCluster_signs.resize( nV );
2561 #endif
2562 
2563  M_face_sign = ublas::scalar_vector<bool>( M.numFaces(), false );
2564 
2565  const bool doperm = ( ( ( Shape == SHAPE_TETRA ) && ( nOrder > 2 ) ) || ( ( Shape == SHAPE_HEXA ) && ( nOrder > 1 ) ) );
2566  DVLOG(2) << "generateFacePermutations: " << doperm << "\n";
2567  generateFacePermutations( M, mpl::bool_<doperm>() );
2568 }
2569 template<typename MeshType, typename FEType, typename PeriodicityType>
2570 size_type
2572 {
2573  size_type nldof =
2574  fe_type::nDofPerVolume * element_type::numVolumes +
2575  fe_type::nDofPerFace * element_type::numGeometricFaces +
2576  fe_type::nDofPerEdge * element_type::numEdges +
2577  fe_type::nDofPerVertex * element_type::numVertices;
2578 
2579  FEELPP_ASSERT( nldof == fe_type::nLocalDof )
2580  ( nldof )
2581  ( fe_type::nLocalDof ).error( "Something wrong in FE specification" ) ;
2582 
2583  const size_type n_proc = M.worldComm().localSize();
2584 
2585 
2586  for ( size_type processor=0; processor<n_proc; processor++ )
2587  {
2588  // compute the number of dof on current processor
2589  element_const_iterator it_elt = M.beginElementWithProcessId( processor );
2590  element_const_iterator en_elt = M.endElementWithProcessId( processor );
2591  size_type n_elts = std::distance( it_elt, en_elt );
2592  VLOG(2) << "[buildDofMap] n_elts = " << n_elts << " on processor " << processor << "\n";
2593  //this->M_first_df[processor] = next_free_dof;
2594 
2595  it_elt = M.beginElementWithProcessId( processor );
2596 
2597  it_elt = M.beginElementWithProcessId( processor );
2598  VLOG(2) << "[buildDofMap] starting with elt " << it_elt->id() << "\n";
2599 
2600  for ( ; it_elt!=en_elt; ++it_elt )
2601  {
2602  element_type const& __elt = *it_elt;
2603  //VLOG(2) << "next_free_dof " << next_free_dof << "\n";
2604  //VLOG(2) << "current dof " << dofIndex( next_free_dof ) << "\n";
2605 
2606  typename element_type::face_const_iterator it, en;
2607  boost::tie( it, en ) = it_elt->faces();
2608 
2609  //bool found_periodic_face_in_element = false;
2610  for ( ; it != en; ++it )
2611  {
2612  if ( ( *it )->marker().value() == M_periodicity.tag2() ||
2613  ( *it )->marker().value() == M_periodicity.tag1() )
2614  {
2615  // store the element reference for the end, the associated
2616  // dof on the periodic face is in fact already taken care of.
2617  // the "internal" dof or on not periodic face will be added
2618  periodic_elements.push_back( boost::make_tuple( boost::addressof( __elt ), *it ) );
2619  //found_periodic_face_in_element = true;
2620  break;
2621  }
2622  }
2623  }
2624  }
2625 
2626  VLOG(2) << "[buildPeriodicDofMap] built periodic_elements " << periodic_elements.size() << "\n";
2627  std::map<size_type,periodic_dof_map_type> periodic_dof;
2628  /*
2629  * Generate the periodic dof, assign a gid to the tag1 dof and set
2630  * the tag2 dof to invalid_size_type_value for now.
2631  */
2632  periodic_element_list_iterator it_periodic = periodic_elements.begin();
2633  periodic_element_list_iterator en_periodic = periodic_elements.end();
2634  size_type next_free_dof = 0;
2635 
2636  while ( it_periodic != en_periodic )
2637  {
2638  element_type const& __elt = *it_periodic->template get<0>();
2639  face_type const& __face = *it_periodic->template get<1>();
2640 
2641  if ( __face.marker().value() == M_periodicity.tag1() )
2642  {
2643  addVertexPeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
2644  addEdgePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
2645  addFacePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
2646  }
2647 
2648  ++it_periodic;
2649  }
2650 
2651  it_periodic = periodic_elements.begin();
2652 
2653  while ( it_periodic != en_periodic )
2654  {
2655  element_type const& __elt = *it_periodic->template get<0>();
2656  face_type const& __face = *it_periodic->template get<1>();
2657 
2658  if ( __face.marker().value() == M_periodicity.tag2() )
2659  {
2660  addVertexPeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
2661  addEdgePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
2662  addFacePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
2663  }
2664 
2665  ++it_periodic;
2666  }
2667 
2668  VLOG(2) << "[periodic dof table] next_free_dof : " << next_free_dof << "\n";
2669  VLOG(2) << "[periodic dof table] number of periodic dof : " << periodic_dof[M_periodicity.tag1()].size() << "\n";
2670 
2671  dof_points_type periodic_dof_points( next_free_dof );
2672  generatePeriodicDofPoints( M, periodic_elements, periodic_dof_points );
2673 
2674  VLOG(2) << "[periodic dof table] generated dof points\n";
2675  VLOG(2) << "[periodic dof table] start matching the dof points\n";
2676 
2677  size_type max_gid = 0;
2678  std::pair<size_type,periodic_dof_type> dof;
2679  BOOST_FOREACH( dof, periodic_dof[M_periodicity.tag1()] )
2680  {
2681  size_type gid = dof.first;
2682  max_gid = ( max_gid > gid )?max_gid:gid;
2683  }
2684 
2685  size_type max_gid2 = 0;
2686  std::pair<size_type,periodic_dof_type> dof2;
2687  BOOST_FOREACH( dof2, periodic_dof[M_periodicity.tag2()] )
2688  {
2689  size_type gid2 = dof2.first;
2690  FEELPP_ASSERT( gid2 > max_gid )( gid2 )( max_gid ).error( "invalid dof index" );
2691  max_gid2 = ( max_gid2 > gid2 )?max_gid2:gid2;
2692  }
2693  CHECK( ( max_gid+1 ) == ( max_gid2+1-( max_gid+1 ) ) )
2694  << "[periodic] invalid periodic setup"
2695  << " max_gid+1 = " << max_gid+1
2696  << ", ( max_gid2+1-( max_gid+1 ) =" << ( max_gid2+1-( max_gid+1 ) )
2697  << ", max gid = " << max_gid
2698  << ", max_gid2 = " << max_gid2 << "\n";
2699 
2700  std::vector<bool> periodic_dof_done( max_gid+1 );
2701  std::fill( periodic_dof_done.begin(), periodic_dof_done.end(), false );
2702 
2703  BOOST_FOREACH( dof, periodic_dof[M_periodicity.tag1()] )
2704  {
2705  size_type gid = dof.first;
2706 
2707  if ( periodic_dof_done[gid] )
2708  continue;
2709 
2710  node_type x1 = periodic_dof_points[gid].template get<0>();
2711  bool match = false;
2712  typename periodic_dof_map_type::iterator it_dof2 = periodic_dof[M_periodicity.tag2()].begin();
2713  typename periodic_dof_map_type::iterator en_dof2 = periodic_dof[M_periodicity.tag2()].end();
2714 #if 0
2715  for ( ; it_dof2 != en_dof2; ++ it_dof2 )
2716  {
2717  size_type gid2 = it_dof2->first;
2718  FEELPP_ASSERT( gid2 < next_free_dof )( gid )( gid2 )( next_free_dof ).error( "[periodic] invalid dof id" );
2719  node_type x2 = periodic_dof_points[gid2].template get<0>();
2720  //FEELPP_ASSERT( math::abs( x2[0]-M_periodicity.translation()[0]) < 1e-10 )
2721  //( x1 )( x2 )( M_periodicity.translation() ).error( "[periodic] invalid periodic setup");
2722  }
2723 #endif
2724  it_dof2 = periodic_dof[M_periodicity.tag2()].begin();
2725  size_type corresponding_gid = invalid_size_type_value;
2726 
2727  for ( ; it_dof2 != en_dof2; ++ it_dof2 )
2728  {
2729  // make sure that we iterate over dof belonging to the same function
2730  // component (e.g. in vectorial)
2731  if ( it_dof2->second.template get<2>() != dof.second.template get<2>() )
2732  continue;
2733  size_type gid2 = it_dof2->first;
2734  FEELPP_ASSERT( gid2 < next_free_dof )( gid )( gid2 )( next_free_dof ).error( "[periodic] invalid dof id" );
2735  node_type x2 = periodic_dof_points[gid2].template get<0>();
2736  FEELPP_ASSERT( ( x1.size() == x2.size() ) &&
2737  ( x1.size() == M_periodicity.translation().size() ) )
2738  ( gid )( dof.second.template get<0>() )( dof.second.template get<1>() )
2739  ( gid2 )( it_dof2->second.template get<0>() )( it_dof2->second.template get<1>() )
2740  ( x1 )( x2 )( M_periodicity.translation() ).error( "invalid point size" );
2741 
2742  if ( ublas::norm_2( x1-( x2-M_periodicity.translation() ) ) < 1e-10 )
2743  {
2744  // loop on each pair (element, lid) which
2745  // has a global id gid2 and set it to gid
2746  corresponding_gid = gid2;
2747  match = true;
2748  break;
2749  }
2750  }
2751 
2752  // if we have --- actually we must have one --- a match, remove the
2753  // iterator from dof2 to quicken the search for the next dof1 match
2754  if ( match )
2755  {
2756  size_type ie1 = dof.second.template get<0>();
2757  size_type lid1 = dof.second.template get<1>();
2758  size_type c1 = dof.second.template get<2>();
2759  size_type gDof1 = dof.second.template get<3>();
2760  uint16_type dof1_type = dof.second.template get<4>();
2761 
2762  VLOG(2) << "matching dof id " << gid << " with dof id=" << corresponding_gid << "\n";
2763 
2764  it_dof2 = periodic_dof[M_periodicity.tag2()].lower_bound( corresponding_gid );
2765  en_dof2 = periodic_dof[M_periodicity.tag2()].upper_bound( corresponding_gid );
2766  VLOG(2) << "distance = " << std::distance( it_dof2, en_dof2 ) << "\n";
2767 
2768  while ( it_dof2 != en_dof2 )
2769  {
2770 
2771  size_type ie = it_dof2->second.template get<0>();
2772  size_type lid = it_dof2->second.template get<1>();
2773  size_type c2 = it_dof2->second.template get<2>();
2774  CHECK( c1 == c2 ) << "[periodic] invalid dof component, c1 = " << c1 << ", c2 = " << c2 << "\n";
2775  size_type gDof = it_dof2->second.template get<3>();
2776  uint16_type dof2_type = it_dof2->second.template get<4>();
2777  uint16_type dof1_type = dof.second.template get<4>();
2778 
2779  FEELPP_ASSERT( dof1_type == dof2_type )
2780  ( gid )( it_dof2->first )( gDof )( lid )( c2) ( ie )
2781  ( dof1_type )( dof2_type ).error ( "invalid dof" );
2782 
2783  VLOG(2) << "link " << M_el_l2g.left.find( LocalDof( ie, localDofId(lid,c2) ) )->second.index() << " -> " << gid << "\n"
2784  << "element id1: " << ie1 << ", lid1: " << lid1 << ", c1: " << c1 << ", gDof1: " << gDof1 << ", type1: " << dof1_type << "\n"
2785  << "element id2: " << ie << ", lid2: " << lid << ", c2: " << c2 << ", gDof2: " << gDof << ", type: " << dof2_type << "\n";
2786 
2787  // gid is given by dof1
2788  auto it = M_el_l2g.left.find(LocalDof( ie, localDofId(lid,c2) ));
2789  bool successful_modify = M_el_l2g.left.modify_data( it, bimaps::_data = Dof( boost::make_tuple( gid, 1, true ) ) );
2790 
2791  CHECK( successful_modify ) << "modify periodic dof table failed: element id "
2792  << ie << " local dof id " << lid << " component " << c2;
2793  // warning: must modify the data structure that allows to
2794  // generate unique global dof ids
2795  CHECK( ( map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] == corresponding_gid ) ||
2796  ( map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] == gid ) )
2797  << "[periodic] invalid matching periodic gid, "
2798  << "corresponding_gid = " << corresponding_gid << ", dof2_type = " << dof2_type
2799  << ", gDof = " << gDof << ", gid=" << gid
2800  << ", c2 = " << c2
2801  << ", map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ]= "
2802  << map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] << "\n";
2803 
2804  VLOG(2) << "link mapgdof " << map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] << " -> " << gid << "\n";
2805  map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] = gid;
2806 
2807  FEELPP_ASSERT( map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] == gid )
2808  ( corresponding_gid )( dof2_type )( gDof )( gid )
2809  ( map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] ) .error ( "invalid gid" );
2810 
2811  ++it_dof2;
2812  }
2813 
2814  it_dof2 = periodic_dof[M_periodicity.tag2()].lower_bound( corresponding_gid );
2815  periodic_dof[M_periodicity.tag2()].erase( it_dof2, en_dof2 );
2816  periodic_dof_done[gid] = true;
2817  }
2818 
2819  else
2820  {
2821  // we have a problem, no match was found, this should not happen
2822  VLOG(2) << "[periodic] invalid point/dof matching\n";
2823  VLOG(2) << "[periodic] n = " << x1 << "\n";
2824  }
2825 
2826  }
2827  VLOG(2) << "[periodic dof table] done matching the dof points\n";
2828  VLOG(2) << "[periodic dof table] is empty : " << periodic_dof[M_periodicity.tag2()].empty() << "\n";
2829 
2830  // ensure that periodic_dof[M_periodicity.tag2()] is empty
2831  if ( !periodic_dof[M_periodicity.tag2()].empty() )
2832  {
2833  VLOG(2) << "[periodic] periodic conditions not set properly, some periodic dof were not assigned\n";
2834  typename periodic_dof_map_type::iterator it_dof2 = periodic_dof[M_periodicity.tag2()].begin();
2835  typename periodic_dof_map_type::iterator en_dof2 = periodic_dof[M_periodicity.tag2()].end();
2836 
2837  while ( it_dof2 != en_dof2 )
2838  {
2839 
2840  size_type ie = it_dof2->second.template get<0>();
2841  size_type lid = it_dof2->second.template get<1>();
2842 
2843  VLOG(2) << "[periodic] dof " << it_dof2->first << " not assigned, "
2844  << "x = " << periodic_dof_points[it_dof2->first].template get<0>() << " "
2845  << "elt = " << ie << ", lid= " << lid << "\n";
2846 
2847 
2848 
2849  ++it_dof2;
2850  }
2851  }
2852 
2853  else
2854  {
2855  VLOG(2) << "[periodic] periodic condition done\n";
2856  }
2857 
2858  return max_gid+1;
2859 }
2860 
2861 template<typename MeshType, typename FEType, typename PeriodicityType>
2862 size_type
2864 {
2865  typedef typename continuity_type::template apply<MeshType, self_type> builder;
2866  return fusion::accumulate( typename continuity_type::discontinuity_markers_type(), start_next_free_dof, builder( M, *this ) );
2867 }
2868 template<typename MeshType, typename FEType, typename PeriodicityType>
2869 void
2871 {
2872  if ( !M_dof_indices.empty() )
2873  {
2874  generateDofPoints( M );
2875  return;
2876  }
2877 
2878  size_type nldof =
2879  fe_type::nDofPerVolume * element_type::numVolumes +
2880  fe_type::nDofPerFace * element_type::numGeometricFaces +
2881  fe_type::nDofPerEdge * element_type::numEdges +
2882  fe_type::nDofPerVertex * element_type::numVertices;
2883 
2884  FEELPP_ASSERT( nldof == fe_type::nLocalDof )
2885  ( nldof )
2886  ( fe_type::nLocalDof ).error( "Something wrong in FE specification" ) ;
2887 
2888 #if !defined(FEELPP_ENABLE_MPI_MODE) // sequential if (this->worldComm().size()==1)
2889 
2890  const size_type n_proc = M.worldComm().localSize();
2891 
2893  std::list<std::pair<element_type const*, face_type const*> > periodic_elements;
2894 
2895  size_type next_free_dof = start_next_free_dof;
2896 
2897  for ( size_type processor=0; processor<n_proc; processor++ )
2898  {
2899  // compute the number of dof on current processor
2900  element_const_iterator it_elt = M.beginElementWithProcessId( processor );
2901  element_const_iterator en_elt = M.endElementWithProcessId( processor );
2902  size_type n_elts = std::distance( it_elt, en_elt );
2903  DVLOG(2) << "[buildDofMap] n_elts = " << n_elts << " on processor " << processor << "\n";
2904  this->M_first_df[processor] = next_free_dof;
2905 
2906  if ( is_periodic || is_discontinuous_locally )
2907  this->M_first_df[processor] = 0;
2908 
2909  it_elt = M.beginElementWithProcessId( processor );
2910 
2911  DVLOG(2) << "[buildDofMap] starting with elt " << it_elt->id() << "\n";
2912 
2913  for ( ; it_elt!=en_elt; ++it_elt )
2914  {
2915  this->addDofFromElement( *it_elt, next_free_dof, processor );
2916  } // elements loop
2917  // printing Dof table only in debug mode
2918 #if !defined( NDEBUG )
2919  const int ncdof = is_product?nComponents:1;
2920 
2921  for ( int c = 0; c < ncdof; ++c )
2922  {
2923  it_elt = M.beginElementWithProcessId( processor );
2924 
2925  for ( ; it_elt!=en_elt; ++it_elt )
2926  {
2927  element_type const& __elt = *it_elt;
2928  std::ostringstream ostr;
2929  ostr << "element " << __elt.id() << " : ";
2930 
2931  for ( uint16_type l = 0; l < nldof; ++l )
2932  {
2933  ostr << M_el_l2g.left.find( LocalDof( __elt.id(), fe_type::nLocalDof*c + l ) )->second.index() << " ";
2934  }
2935 
2936  DVLOG(2) << ostr.str() << "\n";
2937  }
2938  }
2939 
2940 #endif
2941  this->M_last_df[processor] = next_free_dof-1;
2942 
2943  // update info
2944  this->M_n_localWithGhost_df[processor] = this->M_last_df[processor] - this->M_first_df[processor] + 1;
2945  this->M_n_localWithoutGhost_df[processor]=this->M_n_localWithGhost_df[processor];
2946  this->M_first_df_globalcluster[processor]=this->M_first_df[processor];
2947  this->M_last_df_globalcluster[processor]=this->M_last_df[processor];
2948  }
2949 
2950  this->M_n_dofs = next_free_dof;
2951 
2952 #if !defined(NDEBUG)
2953  DVLOG(2) << " n global dof " << nDof() << "\n";
2954  DVLOG(2) << " n local dof " << nLocalDof() << "\n";
2955 #endif
2956 
2957  for ( size_type processor=0; processor<M.worldComm().localSize(); processor++ )
2958  {
2959  DVLOG(2) << "o processor " << processor << "\n";
2960  DVLOG(2) << " - n dof on proc " << nDofOnProcessor( processor ) << "\n";
2961  DVLOG(2) << " - first dof " << firstDof( processor ) << "\n";
2962  DVLOG(2) << " - last dof " << lastDof( processor ) << "\n";
2963  }
2964 
2965  //if ( is_continuous )
2966  //checkDofContinuity( M, mpl::int_<fe_type::nDim>() );
2967 
2968  for ( size_type processor=0; processor<n_proc; processor++ )
2969  {
2970  auto it_elt = M.beginElementWithProcessId( processor );
2971  auto en_elt = M.endElementWithProcessId( processor );
2972 
2973  for ( ; it_elt != en_elt; ++it_elt )
2974  {
2975  size_type elid= it_elt->id();
2976 
2977  for ( int i = 0; i < FEType::nLocalDof; ++i )
2978  {
2979  int nc1 = ( is_product?nComponents1:1 );
2980 
2981  for ( int c1 =0; c1 < nc1; ++c1 )
2982  {
2983  int ind = FEType::nLocalDof*c1+i;
2984  boost::tie( M_locglob_indices[elid][ind],
2985  M_locglob_signs[elid][ind], boost::tuples::ignore ) =
2986  localToGlobal( elid, i, c1 );
2987  }
2988  }
2989  }
2990  }
2991 
2992 #else // MPI_MODE
2993  // compute the number of dof on current processor
2994  auto it_elt = M.beginElementWithProcessId( M.worldComm().localRank() );
2995  auto en_elt = M.endElementWithProcessId( M.worldComm().localRank() );
2996  bool hasNoElt = ( it_elt == en_elt );
2997 
2998  //size_type n_elts = std::distance( it_elt, en_elt);
2999  //DVLOG(2) << "[buildDofMap] n_elts = " << n_elts << " on processor " << processor << "\n";
3000 
3001  size_type theFirstDf = start_next_free_dof;
3002 
3003  if ( is_periodic || is_discontinuous_locally )
3004  theFirstDf = 0;
3005 
3006  mpi::all_gather( M.worldComm().localComm(),
3007  theFirstDf,//start_next_free_dof,
3008  this->M_first_df );
3009 
3010  //if ( is_periodic || is_discontinuous_locally )
3011  // this->M_first_df[processor] = 0;
3012 
3013  size_type next_free_dof = start_next_free_dof;
3014 
3015  for ( ; it_elt!=en_elt; ++it_elt )
3016  {
3017  this->addDofFromElement( *it_elt, next_free_dof, M.worldComm().localRank() );
3018  } // elements loop
3019 
3020 #if 0
3021  for ( auto mit = M_dof_marker.right.begin(), men = M_dof_marker.right.end() ; mit != men ; ++mit )
3022  {
3023  LOG(INFO) << "marker " << mit->first << " dof id " << mit->second;
3024  }
3025 #endif
3026 
3027 #if 0
3028  LOG(INFO) << "local to global view";
3029  for( auto it = M_el_l2g.left.begin(), en = M_el_l2g.left.end();
3030  it != en; ++it )
3031  {
3032  LOG(INFO) << "local dof (" << it->first.elementId()<< ","<< it->first.localDof() << ") --> global dof " << it->second.index();
3033  }
3034  LOG(INFO) << "global to local view";
3035  for( auto it = M_el_l2g.right.begin(), en = M_el_l2g.right.end();
3036  it != en; ++it )
3037  {
3038  LOG(INFO) << "global dof " << it->first.index() << " --> local dof (" << it->second.elementId()<< ","<< it->second.localDof() << ")";
3039  }
3040 #endif
3041  const size_type thelastDof = ( !hasNoElt )?next_free_dof-1:0;
3042 
3043  mpi::all_gather( M.worldComm().localComm(),
3044  thelastDof,//next_free_dof-1,
3045  this->M_last_df );
3046 
3047  // access to M_n_localWithGhost_df for each process
3048  size_type mynDofWithGhost = ( !hasNoElt )?
3049  this->M_last_df[M.worldComm().localRank()] - this->M_first_df[M.worldComm().localRank()] + 1 :
3050  this->M_first_df[M.worldComm().localRank()];
3051 
3052  // update info with all_gather
3053  mpi::all_gather( M.worldComm().localComm(),
3054  mynDofWithGhost,
3055  this->M_n_localWithGhost_df );
3056 #if 0
3057  std::cout << "\n build Dof Map --2---with god rank " << this->worldComm().godRank()
3058  << " local rank DofT " << this->worldComm().localRank()
3059  << " local rank mesh " << M.worldComm().localRank()
3060  << std::endl;
3061 #endif
3062 
3063  // only true in sequential, redefine in buildDofGhostMap
3064  this->M_n_localWithoutGhost_df[this->worldComm().localRank()]=this->M_n_localWithGhost_df[this->worldComm().localRank()];
3065  this->M_first_df_globalcluster[this->worldComm().localRank()]=this->M_first_df[this->worldComm().localRank()];
3066  this->M_last_df_globalcluster[this->worldComm().localRank()]=this->M_last_df[this->worldComm().localRank()];
3067  this->M_n_dofs = next_free_dof;
3068 
3069  it_elt = M.beginElementWithProcessId( M.worldComm().localRank() );
3070 
3071  for ( ; it_elt != en_elt; ++it_elt )
3072  {
3073  size_type elid= it_elt->id();
3074 
3075  for ( int i = 0; i < FEType::nLocalDof; ++i )
3076  {
3077  int nc1 = ( is_product?nComponents1:1 );
3078 
3079  for ( int c1 =0; c1 < nc1; ++c1 )
3080  {
3081  int ind = FEType::nLocalDof*c1+i;
3082  auto const& dof = localToGlobal( elid, i, c1 );
3083  M_locglob_indices[elid][ind] = dof.index();
3084  M_locglob_signs[elid][ind] = dof.sign();
3085  }
3086  }
3087  }
3088 
3089 #endif // MPI_MODE
3090  generateDofPoints( M );
3091 
3092 }
3093 
3094 template<typename MeshType, typename FEType, typename PeriodicityType>
3095 void
3097 {
3098  size_type nDofF = nLocalDofOnFace(true);
3099  M_n_dof_per_face_on_bdy = nDofF;
3100  DVLOG(2) << "vertex dof : " << face_type::numVertices * fe_type::nDofPerVertex << "\n";
3101  DVLOG(2) << "edge dof : " << face_type::numEdges * fe_type::nDofPerEdge << "\n";
3102  DVLOG(2) << "face dof : " << face_type::numFaces * fe_type::nDofPerFace << "\n";
3103  DVLOG(2) << "number of Dof on an Element Face : " << nDofF << "\n";
3104 
3105  if ( nDofF == 0 ) return;
3106 
3107  //
3108  // Face dof
3109  //
3110 #if defined(FEELPP_ENABLE_MPI_MODE)
3111  auto __face_it = M.facesWithProcessId( M.worldComm().localRank() ).first;
3112  auto __face_en = M.facesWithProcessId( M.worldComm().localRank() ).second;
3113 #else
3114  typename mesh_type::face_const_iterator __face_it = M.beginFace();
3115  typename mesh_type::face_const_iterator __face_en = M.endFace();
3116 #endif
3117 
3118  const size_type nF = M.faces().size();
3119  int ntldof = nLocalDofOnFace();
3120 
3121  DVLOG(2) << "[buildBoundaryDofMap] nb faces : " << nF << "\n";
3122  DVLOG(2) << "[buildBoundaryDofMap] nb dof faces : " << nDofF*nComponents << "\n";
3123 
3124  for ( size_type nf = 0; __face_it != __face_en; ++__face_it, ++nf )
3125  {
3126  LOG_IF(WARNING, !__face_it->isConnectedTo0() )
3127  << "face " << __face_it->id() << " not connected"
3128  << " marker : " << __face_it->marker()
3129  << " connectedTo0 : " << __face_it->isConnectedTo0()
3130  << " connectedTo1 : " << __face_it->isConnectedTo1();
3131 
3132  if ( !__face_it->isConnectedTo0() ) continue;
3133 
3134 #if !defined(NDEBUG)
3135 
3136  if ( __face_it->isOnBoundary() )
3137  DVLOG(4) << "[buildBoundaryDofMap] boundary global face id : " << __face_it->id()
3138  << " marker: " << __face_it->marker()<< "\n";
3139 
3140  else
3141  DVLOG(4) << "[buildBoundaryDofMap] global face id : " << __face_it->id() << "\n";
3142 
3143 #endif
3144  uint16_type lcVertex = 0;
3145  uint16_type lcEdge = 0;
3146  uint16_type lcFace = 0;
3147 
3148 
3149  addVertexBoundaryDof( __face_it, lcVertex );
3150  addEdgeBoundaryDof( __face_it, lcEdge );
3151  addFaceBoundaryDof( __face_it, lcFace );
3152 
3153  }
3154 
3155 #if !defined(NDEBUG)
3156  __face_it = M.facesWithProcessId( M.worldComm().localRank() ).first;
3157  __face_en = M.facesWithProcessId( M.worldComm().localRank() ).second;
3158  for ( ; __face_it != __face_en; ++__face_it )
3159  for ( int face_dof_id = 0; face_dof_id < int( ntldof ); ++face_dof_id )
3160  FEELPP_ASSERT( boost::get<0>( M_face_l2g[__face_it->id()][face_dof_id] ) != invalid_size_type_value )( __face_it->id() )( face_dof_id ).warn( "invalid dof table: initialized dof entries" );
3161 
3162 #endif
3163 } // updateBoundaryDof
3164 
3165 template<typename MeshType, typename FEType, typename PeriodicityType>
3166 void
3168 {
3169  if ( !M_dof_points.empty() )
3170  return;
3171 
3172  if ( fe_type::is_modal )
3173  return;
3174 
3175  DVLOG(2) << "[Dof::generateDofPoints] generating dof coordinates\n";
3176  typedef typename gm_type::template Context<vm::POINT, element_type> gm_context_type;
3177  typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;
3178 
3179  typedef typename fe_type::template Context<vm::POINT, fe_type, gm_type, element_type> fecontext_type;
3180 
3181  gm_ptrtype gm( new gm_type );
3182  fe_type fe;
3183  //
3184  // Precompute some data in the reference element for
3185  // geometric mapping and reference finite element
3186  //
3187  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( gm, fe.points() ) );
3188 
3189  //const uint16_type ndofv = fe_type::nDof;
3190 
3191  element_const_iterator it_elt = M.beginElementWithProcessId( M.worldComm().localRank() );
3192  element_const_iterator en_elt = M.endElementWithProcessId( M.worldComm().localRank() );
3193 
3194  if ( it_elt == en_elt )
3195  return;
3196 
3197  gm_context_ptrtype __c( new gm_context_type( gm, *it_elt, __geopc ) );
3198 
3199  std::vector<bool> dof_done( nLocalDofWithGhost() );
3200  M_dof_points.resize( nLocalDofWithGhost() );
3201  std::fill( dof_done.begin(), dof_done.end(), false );
3202 
3203  for ( size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
3204  {
3205  __c->update( *it_elt );
3206 
3207  for ( uint16_type l =0; l < fe_type::nLocalDof; ++l )
3208  {
3209  int ncdof = is_product?nComponents:1;
3210 
3211  for ( uint16_type c1 = 0; c1 < ncdof; ++c1 )
3212  {
3213  size_type thedof = boost::get<0>( localToGlobal( it_elt->id(), l, c1 ) );
3214 
3215  if ( ( thedof >= firstDof() ) && ( thedof <= lastDof() ) )
3216  {
3217  // get only the local dof
3218  //size_type thedofonproc = thedof - firstDof();
3219  thedof -= firstDof();
3220  DCHECK( thedof < nLocalDofWithGhost() )
3221  << "invalid local dof index "
3222  << thedof << ", " << nLocalDofWithGhost() << "," << firstDof() << ","
3223  << lastDof() << "," << it_elt->id() << "," << l << "," << c1;
3224 
3225  if ( dof_done[ thedof ] == false )
3226  {
3227  std::set<uint16_type> lid;
3228  lid.insert( l );
3229  //M_dof_points[dof_id] = boost::make_tuple( thedof, __c->xReal( l ) );
3230  M_dof_points[thedof] = boost::make_tuple( __c->xReal( l ), firstDof()+thedof, c1 );
3231  dof_done[thedof] = true;
3232  ++dof_id;
3233  }
3234  }
3235  }
3236  }
3237  }
3238 
3239  for ( size_type dof_id = 0; dof_id < nLocalDofWithGhost() ; ++dof_id )
3240  {
3241  CHECK( boost::get<1>( M_dof_points[dof_id] ) >= firstDof() &&
3242  boost::get<1>( M_dof_points[dof_id] ) <= lastDof() )
3243  << "invalid dof point "
3244  << dof_id << ", " << firstDof() << ", " << lastDof() << ", " << nLocalDofWithGhost()
3245  << ", " << boost::get<1>( M_dof_points[dof_id] )
3246  << ", " << boost::get<0>( M_dof_points[dof_id] ) ;
3247  CHECK( dof_done[dof_id] == true )
3248  << "invalid dof point"
3249  << dof_id << ", " << nLocalDofWithGhost() << ", " << firstDof() << ", "
3250  << lastDof() << ", " << fe_type::nDim << ", " << fe_type::nLocalDof;
3251  }
3252 
3253  DVLOG(2) << "[Dof::generateDofPoints] generating dof coordinates done\n";
3254 }
3255 template<typename MeshType, typename FEType, typename PeriodicityType>
3256 void
3258  periodic_element_list_type const& periodic_elements,
3259  dof_points_type& periodic_dof_points )
3260 {
3261  if ( fe_type::is_modal )
3262  return;
3263 
3264  DVLOG(2) << "[Dof::generateDofPoints] generating dof coordinates\n";
3265  typedef typename gm_type::template Context<vm::POINT, element_type> gm_context_type;
3266  typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;
3267 
3268  typedef typename fe_type::template Context<vm::POINT, fe_type, gm_type, element_type> fecontext_type;
3269 
3270  gm_ptrtype gm( new gm_type );
3271  fe_type fe;
3272  //
3273  // Precompute some data in the reference element for
3274  // geometric mapping and reference finite element
3275  //
3276  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( gm, fe.points() ) );
3277 
3278  //const uint16_type ndofv = fe_type::nDof;
3279 
3280  periodic_element_list_const_iterator it_elt = periodic_elements.begin();
3281  periodic_element_list_const_iterator en_elt = periodic_elements.end();
3282 
3283  if ( it_elt == en_elt )
3284  return;
3285 
3286  gm_context_ptrtype __c( new gm_context_type( gm, *it_elt->template get<0>(), __geopc ) );
3287 
3288  std::vector<bool> dof_done( periodic_dof_points.size() );
3289  std::fill( dof_done.begin(), dof_done.end(), false );
3290 
3291  for ( size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
3292  {
3293  __c->update( *it_elt->template get<0>() );
3294 
3295  face_type const& __face = *it_elt->template get<1>();
3296 
3297  size_type iElAd = __face.ad_first();
3298  FEELPP_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[periodic]invalid face/element in face" );
3299  Feel::detail::ignore_unused_variable_warning( iElAd );
3300 
3301  // local id of the face in its adjacent element
3302  uint16_type iFaEl = __face.pos_first();
3303  FEELPP_ASSERT( iFaEl != invalid_uint16_type_value ).error ( "invalid element index in face" );
3304 
3305  int ncdof = is_product?nComponents:1;
3306 
3307  for ( uint16_type c1 = 0; c1 < ncdof; ++c1 )
3308  {
3309  // loop on face vertices
3310  for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
3311  {
3312  // local vertex number (in element)
3313  uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
3314  Feel::detail::ignore_unused_variable_warning( iVeEl );
3315 
3316  FEELPP_ASSERT( iVeEl != invalid_uint16_type_value ).error( "invalid local dof" );
3317 
3318  // Loop number of Dof per vertex
3319  for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
3320  {
3321  uint16_type lid = iVeEl * fe_type::nDofPerVertex + l;
3322 
3323  size_type thedof = boost::get<0>( localToGlobal( it_elt->template get<0>()->id(), lid, c1 ) );
3324  FEELPP_ASSERT( thedof < dof_done.size() )
3325  ( thedof )
3326  ( dof_done.size() )
3327  ( c1 )
3328  ( it_elt->template get<0>()->id() )
3329  ( lid ).error ( "[generatePeriodicDofPoints] invalid dof id" );
3330 
3331  if ( dof_done[ thedof ] == false )
3332  {
3333  periodic_dof_points[thedof] = boost::make_tuple( __c->xReal( lid ), thedof, c1 );
3334  // these tests are problem specific x=0 and x=translation
3335 #if 0
3336 
3337  if ( __face.marker().value() == M_periodicity.tag1() )
3338  FEELPP_ASSERT( math::abs( __c->xReal( lid )[0] ) < 1e-10 )( __c->xReal( lid ) ).warn( "[periodic] invalid p[eriodic point tag1" );
3339 
3340  if ( __face.marker().value() == M_periodicity.tag2() )
3341  FEELPP_ASSERT( math::abs( __c->xReal( lid )[0] - M_periodicity.translation()[0] ) < 1e-10 )
3342  ( __c->xReal( lid ) )( M_periodicity.translation() ).warn( "[periodic] invalid p[eriodic point tag1" );
3343 
3344 #endif
3345  dof_done[thedof] = true;
3346  ++dof_id;
3347  }
3348  }
3349  // loop on edge
3350  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
3351  {
3352  uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
3353  size_type thedof = boost::get<0>( localToGlobal( it_elt->template get<0>()->id(), lid, c1 ) );
3354  FEELPP_ASSERT( thedof < dof_done.size() )
3355  ( thedof )
3356  ( dof_done.size() )
3357  ( c1 )
3358  ( it_elt->template get<0>()->id() )
3359  ( lid ).error ( "[generatePeriodicDofPoints] invalid dof id" );
3360 
3361  if ( dof_done[ thedof ] == false )
3362  {
3363  periodic_dof_points[thedof] = boost::make_tuple( __c->xReal( lid ), thedof, c1 );
3364  // these tests are problem specific x=0 and x=translation
3365 #if 0
3366 
3367  if ( __face.marker().value() == M_periodicity.tag1() )
3368  FEELPP_ASSERT( math::abs( __c->xReal( lid )[1] +1 ) < 1e-10 )( __c->xReal( lid ) ).warn( "[periodic] invalid p[eriodic point tag1" );
3369 
3370  if ( __face.marker().value() == M_periodicity.tag2() )
3371  FEELPP_ASSERT( math::abs( __c->xReal( lid )[1] - ( M_periodicity.translation()[1]-1 ) ) < 1e-10 )
3372  ( __c->xReal( lid ) )( M_periodicity.translation() ).warn( "[periodic] invalid p[eriodic point tag1" );
3373 
3374 #endif
3375  dof_done[thedof] = true;
3376  ++dof_id;
3377  }
3378  }
3379  }
3380  }
3381  }
3382  for ( size_type dof_id = 0; dof_id < periodic_dof_points.size() ; ++dof_id )
3383  {
3384  FEELPP_ASSERT( boost::get<1>( periodic_dof_points[dof_id] ) >= 0 &&
3385  boost::get<1>( periodic_dof_points[dof_id] ) < periodic_dof_points.size() )
3386  ( dof_id )( periodic_dof_points.size() )
3387  ( boost::get<1>( periodic_dof_points[dof_id] ) )
3388  ( boost::get<0>( periodic_dof_points[dof_id] ) ).error( "invalid dof point" );
3389  FEELPP_ASSERT( dof_done[dof_id] == true )( dof_id ).error( "invalid dof point" );
3390  }
3391 }
3392 
3393 
3394 template<typename MeshType, typename FEType, typename PeriodicityType>
3395 void
3397 {
3398  addSubstructuringDofVertex( M, next_free_dof );
3399  addSubstructuringDofEdge( M, next_free_dof, mpl::int_<nDim>() );
3400  addSubstructuringDofFace( M, next_free_dof, mpl::int_<nDim>() );
3401 }
3402 
3403 template<typename MeshType, typename FEType, typename PeriodicityType>
3404 void
3406  size_type next_free_dof )
3407 {
3408  std::cout << "found CrossPoints and WireBasket\n";
3409  std::cout << "n cp: " << std::distance( M.beginPointWithMarker( M.markerName("CrossPoints") ), M.endPointWithMarker( M.markerName("CrossPoints") ) ) << "\n";
3410 #if 0
3411  std::cout << "n wb: " << std::distance( M.beginEdgeWithMarker( M.markerName("WireBasket") ), M.endEdgeWithMarker( M.markerName("WireBasket") ) ) << "\n";
3412 #endif
3413  // go through all the crosspoints and add them to the dof table
3414 
3415  for( auto pit = M.beginPointWithMarker( M.markerName("CrossPoints") ),
3416  pen = M.endPointWithMarker( M.markerName("CrossPoints") );
3417  pit!=pen; ++pit )
3418  {
3419  // get one element
3420  auto __elt = M.element( *pit->elements().begin() );
3421  size_type ie = __elt.id();
3422  int lc = 0;
3423  for ( uint16_type i = 0; i < element_type::numVertices; ++i )
3424  {
3425  for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l, ++lc )
3426  {
3427  if (__elt.point( i ).id()==pit->id() )
3428  {
3429  const size_type gDof = ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
3430  this->insertDof( ie, lc, i, boost::make_tuple( 0, 0, gDof ),
3431  M.worldComm().localRank(), next_free_dof, 1, false, 0 );
3432  std::cout << "Adding crosspoint " << pit->id() << " with dof " << next_free_dof << "\n";
3433  }
3434  }
3435  }
3436  }
3437 }
3438 
3439 template<typename MeshType, typename FEType, typename PeriodicityType>
3440 void
3442  size_type next_free_dof,
3443  mpl::int_<1> )
3444 {}
3445 
3446 template<typename MeshType, typename FEType, typename PeriodicityType>
3447 void
3449  size_type next_free_dof,
3450  mpl::int_<2> )
3451 {}
3452 
3453 template<typename MeshType, typename FEType, typename PeriodicityType>
3454 void
3456  size_type next_free_dof,
3457  mpl::int_<3> )
3458 {
3459  // go through all Wirebasket edges
3460  for( auto pit = M.beginEdgeWithMarker( M.markerName("WireBasket") ),
3461  pen = M.endEdgeWithMarker( M.markerName("WireBasket") );
3462  pit!=pen; ++pit )
3463  {
3464  auto __elt = M.element( *pit->elements().begin() );
3465  std::cout << "Adding wirebasket edge " << pit->id() << " using element " << __elt.id() << "\n";
3466  size_type ie = __elt.id();
3467  uint16_type lc = 0;
3468 
3469  for ( uint16_type i = 0; i < element_type::numEdges; ++i )
3470  {
3471  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
3472  {
3473  if (__elt.edge( i ).id()==pit->id() )
3474  {
3475  size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;
3476  int32_type sign = 1;
3477 
3478  if ( __elt.edgePermutation( i ).value() == edge_permutation_type::IDENTITY )
3479  {
3480  gDof += l ; // both nodal and modal case
3481  }
3482  else if ( __elt.edgePermutation( i ).value() == edge_permutation_type::REVERSE_PERMUTATION )
3483  {
3484 
3485  if ( fe_type::is_modal )
3486  {
3487  //only half of the modes (odd polynomial order) are negative.
3488  sign = ( l%2 )?( -1 ):( 1 );
3489  gDof += l;
3490  }
3491 
3492  else
3493  gDof += fe_type::nDofPerEdge - 1 - l ;
3494  }
3495  else
3496  FEELPP_ASSERT( 0 ).error ( "invalid edge permutation" );
3497 
3498  this->insertDof( ie, lc, i, boost::make_tuple( 1, 0, gDof ), M.worldComm().localRank(), next_free_dof, sign, false, 0 );
3499  std::cout << "Adding wirebasket edge " << pit->id() << " with dof " << next_free_dof << "\n";
3500  }
3501  }
3502  }
3503 
3504  }
3505 }
3506 template<typename MeshType, typename FEType, typename PeriodicityType>
3507 void
3509  size_type next_free_dof,
3510  mpl::int_<1> )
3511 {}
3512 
3513 template<typename MeshType, typename FEType, typename PeriodicityType>
3514 void
3516  size_type next_free_dof,
3517  mpl::int_<2> )
3518 {}
3519 
3520 template<typename MeshType, typename FEType, typename PeriodicityType>
3521 void
3523  size_type next_free_dof,
3524  mpl::int_<3> )
3525 {
3526  std::vector<std::string> faces = assign::list_of("TOP")("BOTTOM")("NORTH")("EAST")("WEST")("SOUTH");
3527  BOOST_FOREACH( auto face, faces )
3528  {
3529  auto faces = markedfaces( &M, face );
3530 
3531  for( auto pit = faces.template get<1>(), pen = faces.template get<2>(); pit!=pen; ++pit )
3532  {
3533  auto __elt = M.element( *pit->elements().begin() );
3534  std::cout << "Adding face " << pit->id() << " with marker " << face << " using element " << __elt.id() << "\n";
3535  size_type ie = __elt.id();
3536  uint16_type lc = 0;
3537 
3538  for ( uint16_type i = 0; i < element_type::numFaces; ++i )
3539  {
3540  face_permutation_type permutation = __elt.facePermutation( i );
3541  FEELPP_ASSERT( permutation != face_permutation_type( 0 ) ).error ( "invalid face permutation" );
3542 
3543  // Polynomial order in each direction
3544  uint16_type p=1;
3545  uint16_type q=0;
3546 
3547  // MaxOrder = Order - 2
3548  int MaxOrder = int( ( 3 + std::sqrt( 1+8*fe_type::nDofPerFace ) )/2 ) - 2;
3549 
3550  for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
3551  {
3552  if (__elt.face( i ).id()==pit->id() )
3553  {
3554  // TODO: orient the dof indices such
3555  // that they match properly the faces
3556  // dof of the connected faces. There
3557  // are a priori many permutations of
3558  // the dof face indices
3559  size_type gDof = __elt.face( i ).id() * fe_type::nDofPerFace;
3560  int32_type sign = 1;
3561 
3562  q=q+1;
3563 
3564  if ( q > MaxOrder )
3565  {
3566  q = 1;
3567  p = p+1;
3568  MaxOrder = MaxOrder-1;
3569  }
3570 
3571  if ( !fe_type::is_modal )
3572  {
3573  // no need of permutation is identity or only one dof on face
3574  if ( permutation == face_permutation_type( 1 ) || fe_type::nDofPerFace == 1 )
3575  gDof += l;
3576 
3577  else
3578  gDof += vector_permutation[permutation][l];
3579  }
3580 
3581  else
3582  {
3583  gDof += l;
3584 
3585  if ( permutation == face_permutation_type( 2 ) )
3586  {
3587  // Reverse sign if polynomial order in
3588  // eta_1 direction is odd
3589 
3590  if ( p%2 == 0 )
3591  sign = -1;
3592 
3593  }
3594  }
3595 
3596  this->insertDof( ie, lc, i, boost::make_tuple( 2, 0, gDof ), M.worldComm().localRank(), next_free_dof, sign, false, 0 );
3597  std::cout << "Adding face " << pit->id() << " with dof " << next_free_dof << "\n";
3598  }
3599  }
3600  }
3601 
3602  }
3603  }
3604 
3605 }
3606 
3607 template<typename MeshType, typename FEType, typename PeriodicityType>
3608 std::pair<std::map<size_type,size_type>,std::map<size_type,size_type> >
3610 {
3611  std::map<size_type,size_type> pidtodof,doftopid;
3612  element_const_iterator it_elt = M_mesh->beginElementWithProcessId( M_mesh->worldComm().localRank() );
3613  element_const_iterator en_elt = M_mesh->endElementWithProcessId( M_mesh->worldComm().localRank() );
3614 
3615  if ( it_elt == en_elt )
3616  return std::make_pair(doftopid,pidtodof);
3617 
3618  for ( size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
3619  {
3620  for ( uint16_type i = 0; i < M_mesh->numLocalVertices(); ++i )
3621  {
3622  int ncdof = is_product?nComponents:1;
3623  for ( uint16_type c1 = 0; c1 < ncdof; ++c1 )
3624  {
3625  const size_type gDof = ( it_elt->point( i ).id() );
3626  size_type thedof = boost::get<0>( localToGlobal( it_elt->id(), i, c1 ) );
3627  //pidtodof[ncdof*it_elt->point(l).id()+c1] = thedof;
3628  pidtodof[ncdof*gDof+c1] = thedof;
3629  doftopid[thedof] = ncdof*gDof+c1;
3630 
3631  }
3632  }
3633  }
3634  if ( !fname.empty() )
3635  {
3636  std::ostringstream os1,os2;
3637  os1 << fs::path( fname ).stem().string() << "_pidtodof" << fs::path( fname ).extension().string();
3638  os2 << fs::path( fname ).stem().string() << "_doftopid" << fs::path( fname ).extension().string();
3639  std::ofstream ofs( os1.str().c_str() );
3640  auto it = pidtodof.begin();
3641  auto en = pidtodof.end();
3642  std::for_each( it, en,
3643  [&ofs]( std::pair<size_type, size_type> const& p )
3644  {
3645  ofs << p.first << " " << p.second << "\n";
3646  });
3647  std::ofstream ofs2( os2.str().c_str() );
3648  it = doftopid.begin();
3649  en = doftopid.end();
3650  std::for_each( it, en,
3651  [&ofs2]( std::pair<size_type, size_type> const& p )
3652  {
3653  ofs2 << p.first << " " << p.second << "\n";
3654  });
3655 
3656  }
3657  return std::make_pair(doftopid,pidtodof);
3658 }
3659 } // namespace Feel
3660 
3661 
3662 
3663 #if defined(FEELPP_ENABLE_MPI_MODE)
3665 #endif
3666 
3667 #endif //_DOFTABLE_HH
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::marker_face_const_iterator, typename MeshTraits< MeshType >::marker_face_const_iterator > markedfaces(MeshType const &mesh)
Definition: filters.hpp:969
mpl::if_< mpl::equal_to< mpl::int_< Elem::nDim >, mpl::int_< nDim > >, global_dof_type, global_dof_fromface_type >::type localToGlobal(Elem const &El, const uint16_type localNode, const uint16_type c=0) const
local to global mapping
Definition: doftable.hpp:620
size_type dofNProc(size_type dof) const
Definition: doftable.hpp:370
void setDofIndices(std::vector< Dof > const &dof)
Definition: doftable.hpp:434
bool insertDof(size_type ie, uint16_type l_dof, uint16_type lc, dof_type gDof, uint16_type processor, size_type &pDof, int32_type sign=1, bool is_dof_periodic=false, size_type shift=0, Marker1 const &marker=Marker1(0))
Definition: doftable.hpp:1142
void setElementDone(size_type elt)
Definition: doftable.hpp:693
Dof global_dof_type
Definition: doftable.hpp:137
size_type localToGlobalId(const size_type ElId, const uint16_type id) const
Definition: doftable.hpp:520
const uint16_type invalid_uint16_type_value
Definition: feelcore/feel.hpp:338
void build(mesh_type *M)
Definition: doftable.hpp:845
LocalDof const & globalToLocal(size_type dof) const
Definition: doftable.hpp:536
localglobal_indices_type const & localToGlobalIndicesOnCluster(size_type ElId)
Definition: doftable.hpp:498
boost::tuple< size_type, uint16_type, uint16_type, uint16_type > local_dof_type
Definition: doftable.hpp:177
std::vector< size_type > M_mapGlobalProcessToGlobalCluster
Definition: datamap.hpp:474
uint16_type numLocalEdges() const
Definition: doftable.hpp:656
size_type nDof() const
Definition: datamap.hpp:99
MeshType mesh_type
Definition: doftable.hpp:77
void clearMapGDof()
Definition: doftable.hpp:1083
ublas::vector< uint16_type > permutation_vector_type
Definition: doftable.hpp:214
size_type numElements() const
Definition: doftable.hpp:640
void addVertexPeriodicDof(element_type const &__elt, face_type const &__face, size_type &next_free_dof, std::map< size_type, periodic_dof_map_type > &periodic_dof, size_type tag)
Definition: doftable.hpp:801
void build(mesh_type &M)
Definition: doftable.hpp:861
void initDofMap(mesh_type &M)
Definition: doftable.hpp:2523
void rebuildDofPoints(mesh_type &M)
Definition: doftable.hpp:1253
dof_map_type & mapGDof()
Definition: doftable.hpp:1075
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
std::vector< size_type > M_n_localWithoutGhost_df
Definition: datamap.hpp:440
std::vector< size_type > M_mapGlobalClusterToGlobalProcess
Definition: datamap.hpp:479
std::vector< size_type > M_first_df
Definition: datamap.hpp:450
Definition: glas.hpp:125
localglobal_indices_type const & localToGlobalIndices(size_type ElId)
Definition: doftable.hpp:490
dof_points_const_iterator dofPointEnd() const
Definition: doftable.hpp:411
class that represents a degree of freedom
Definition: dof.hpp:48
ublas::vector< bool > face_sign_info_type
Definition: doftable.hpp:205
void build(boost::shared_ptr< mesh_type > &M)
Definition: doftable.hpp:853
void buildBoundaryDofMap(mesh_type &mesh)
Build the localToGlobal table for the boundary.
Definition: doftable.hpp:3096
size_type lid(size_type GID) const
Returns local ID of global ID, return invalid_size_type_value if not found on this processor...
Definition: datamap.hpp:263
size_type nLocalDof() const
Definition: datamap.hpp:107
uint16_type nDofPerFaceOnBoundary() const
Definition: doftable.hpp:288
std::map< size_type, std::list< local_dof_type > > dof_element_type
Definition: doftable.hpp:185
boost::tuple< size_type, int16_type, bool, int16_type > global_dof_fromface_type
Definition: doftable.hpp:148
data layout in a multi-processor environnement
Definition: datamap.hpp:46
void buildDofMap(mesh_type &mesh, size_type start_next_free_dof=0)
Build the localToGlobal table.
Definition: doftable.hpp:2870
boost::bimap< bimaps::unordered_set_of< LocalDof >, bimaps::multiset_of< Dof > > dof_table
Definition: doftable.hpp:156
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
localglobal_indices_type const & localToGlobalSigns(size_type ElId)
Definition: doftable.hpp:506
uint16_type numLocalVertices() const
Definition: doftable.hpp:648
std::set< size_type > dofProcSet(size_type dof) const
Definition: doftable.hpp:378
const dof_point_type & dofPoint(size_type i) const
Definition: doftable.hpp:387
std::map< size_type, std::pair< size_type, size_type > > localToGlobalMap() const
Definition: doftable.hpp:778
WorldComm const & worldComm() const
Definition: datamap.hpp:396
size_type dofIndex(size_type dof) const
Definition: doftable.hpp:474
Local-to-global Degree of Freedom table.
Definition: doftable.hpp:69
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
dof_points_iterator dofPointEnd()
Definition: doftable.hpp:419
dof_map_type const & mapGDof() const
Definition: doftable.hpp:1067
size_type localToGlobal(size_type l) const
Definition: doftable.hpp:786
std::vector< size_type > M_first_df_globalcluster
Definition: datamap.hpp:460
uint16_type numLocalFaces() const
Definition: doftable.hpp:664
dof_points_iterator dofPointBegin()
Definition: doftable.hpp:403
std::vector< size_type > M_n_localWithGhost_df
Definition: datamap.hpp:445
std::vector< size_type > M_last_df_globalcluster
Definition: datamap.hpp:465
std::vector< size_type > M_last_df
Definition: datamap.hpp:455
void setMapGDof(dof_map_type const &mapdof)
Definition: doftable.hpp:1091
Context class.
Definition: feelcore/context.hpp:63
dof_points_const_iterator dofPointBegin() const
Definition: doftable.hpp:395
void showMe() const
Definition: doftable.hpp:2244
size_type buildLocallyDiscontinuousDofMap(mesh_type &M, size_type start_next_free_dof)
Definition: doftable.hpp:2863
size_type M_n_dofs
Definition: datamap.hpp:435
std::pair< std::map< size_type, size_type >, std::map< size_type, size_type > > pointIdToDofRelation(std::string fname="") const
Definition: doftable.hpp:3609
size_type buildPeriodicDofMap(mesh_type &M)
Definition: doftable.hpp:2571
size_type nLocalDofWithGhost() const
Definition: datamap.hpp:131

Generated on Sun Dec 22 2013 13:10:58 for Feel++ by doxygen 1.8.5