29 #include <unordered_map>
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>
46 #include<Eigen/StdVector>
48 #include <feel/feelcore/feel.hpp>
52 #include <feel/feelalg/datamap.hpp>
58 namespace fusion = boost::fusion;
59 namespace bimaps = boost::bimaps;
68 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
78 typedef FEType fe_type;
79 typedef boost::shared_ptr<FEType> fe_ptrtype;
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;
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;
94 typedef typename reference_convex_type::super convex_type;
96 typedef typename element_type::edge_permutation_type edge_permutation_type;
97 typedef typename element_type::face_permutation_type face_permutation_type;
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;
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;
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;
118 static const bool is_p0_continuous = ( ( nOrder == 0 ) && is_continuous );
120 static const uint16_type nDofPerElement = mpl::if_<mpl::bool_<is_product>, mpl::int_<FEType::nLocalDof*nComponents1>, mpl::int_<FEType::nLocalDof> >::type::value;
122 typedef PeriodicityType periodicity_type;
123 static const bool is_periodic = periodicity_type::is_periodic;
126 typedef typename fe_type::continuity_type continuity_type;
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;
160 typedef typename std::map<int,global_dof_type> indices_per_element_type;
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;
177 typedef boost::tuple<size_type, uint16_type, uint16_type, uint16_type>
local_dof_type;
178 typedef LocalDofSet local_dof_set_type;
180 typedef boost::tuple<uint16_type&,size_type&> ref_shift_type;
187 typedef boost::bimap<size_type,boost::bimaps::multiset_of<size_type> > dof_marker_type;
188 typedef typename dof_marker_type::value_type dof2marker;
190 typedef typename dof_element_type::iterator dof_iterator;
191 typedef typename dof_element_type::const_iterator dof_const_iterator;
193 typedef std::list<local_dof_type>::const_iterator ldof_const_iterator;
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;
200 typedef std::map<size_type, std::set<size_type> > dof_procset_type;
208 typedef Eigen::Matrix<int, nDofPerElement, 1> localglobal_indices_type;
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 , uint16_type , uint16_type ,
size_type , uint16_type > periodic_dof_type;
221 typedef std::multimap<
size_type , periodic_dof_type> periodic_dof_map_type;
223 DofTable( WorldComm
const& _worldComm )
234 DofTable( fe_ptrtype
const& _fe, periodicity_type
const& periodicity, WorldComm
const& _worldComm );
249 DofTable(
mesh_type& mesh, fe_ptrtype
const& _fe, periodicity_type
const& periodicity, WorldComm
const& _worldComm );
255 M_dof_points.clear();
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);
268 constexpr
size_type nLocalDofOnFace(
bool per_component =
false )
const
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 );
277 local_dof_set_type
const&
280 return M_local_dof_set.update( eid );
283 mesh_type* mesh()
const {
return M_mesh; }
290 return M_n_dof_per_face_on_bdy;
293 indices_per_element_type indices(
size_type id_el )
const
295 indices_per_element_type ind;
296 BOOST_FOREACH( LocalDof
const& ldof, localDofSet( id_el ) )
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;
305 constexpr
size_type getIndicesSize()
const
309 std::vector<size_type> getIndices(
size_type id_el )
const
311 std::vector<size_type> ind(
nLocalDof() );
312 getIndicesSet( id_el, ind );
317 std::vector<size_type> getIndices(
size_type id_el, mpl::size_t<MESH_ELEMENTS> )
const
319 return getIndices( id_el );
322 void getIndicesSet(
size_type id_el, std::vector<size_type>& ind )
const
324 BOOST_FOREACH( LocalDof
const& ldof, this->localDofSet( id_el ) )
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();
332 std::vector<size_type> getIndices(
size_type id_el, mpl::size_t<MESH_FACES> )
const
334 std::vector<size_type> ind( nLocalDofOnFace() );
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 );
347 void getIndicesSetOnGlobalCluster(
size_type id_el, std::vector<size_type>& ind )
const
349 BOOST_FOREACH( LocalDof
const& ldof, this->localDofSet( id_el ) )
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() ];
358 std::vector<size_type> getIndicesOnGlobalCluster(
size_type id_el )
const
361 std::vector<size_type> ind( s );
362 getIndicesSetOnGlobalCluster( id_el, ind );
372 return M_dof_procset.find( dof )->second.size();
380 return M_dof_procset.find( dof )->second;
389 return M_dof_points[i];
397 return M_dof_points.begin();
405 return M_dof_points.begin();
413 return M_dof_points.end();
421 return M_dof_points.end();
424 periodic_element_list_const_iterator beginPeriodicElements()
const {
return periodic_elements.begin(); }
425 periodic_element_list_const_iterator endPeriodicElements()
const {
return periodic_elements.end(); }
436 M_dof_indices.resize( dof.size() );
437 std::copy( dof.begin(), dof.end(), M_dof_indices.begin() );
443 std::set<size_type> eltid;
444 std::set<size_type> dofs;
446 BOOST_FOREACH(
Dof thedof, dof )
448 eltid.insert( boost::get<0>( thedof ) );
449 dofs.insert( boost::get<2>( thedof ) );
453 BOOST_FOREACH(
Dof thedof, dof )
455 M_el_l2g.insert( dof_relation( LocalDof( boost::get<0>( thedof ), boost::get<1>( thedof ) ),
456 Dof( boost::get<2>( thedof ), 0,
false ) ) );
459 int processor = this->
worldComm().localRank();
462 this->
M_last_df[processor] = dofs.size()-1;
479 if ( M_dof_indices.empty() )
482 FEELPP_ASSERT( dof < M_dof_indices.size() )( dof )( M_dof_indices.size() ).warn(
"invalid dof index" );
483 return M_dof_indices[dof];
492 return M_locglob_indices[ElId];
500 return M_locglobOnCluster_indices[ElId];
508 return M_locglob_signs[ElId];
521 const uint16_type
id )
const
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();
538 auto it = M_el_l2g.right.find(
Dof( dof ) );
539 DCHECK( it != M_el_l2g.right.end() ) <<
"Invalid global dof entry ( " << dof <<
")";
543 uint16_type localDofId( uint16_type
const lid, uint16_type
const c = 0 )
const
545 return fe_type::nLocalDof * c +
lid;
548 const uint16_type localNode,
549 const uint16_type c = 0 )
const
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 <<
")";
558 const uint16_type localNode,
559 const uint16_type c = 0 )
const
561 Dof resloc = M_el_l2g.left.find( LocalDof( ElId, fe_type::nLocalDof * c + localNode ) )->second;
562 resloc.setIndex( this->mapGlobalProcessToGlobalCluster()[resloc.index()] );
567 const uint16_type localNode,
568 const uint16_type c = 0 )
const
570 const size_type nDofF = nLocalDofOnFace(
true );
571 return M_face_l2g.find( ElId )->second.find( nDofF*c+localNode )->second;
574 struct element_access
576 element_access( DofTable
const& __d )
582 return M_d.M_el_l2g.left.find( LocalDof(__id, M_d.nLocalDof(
true) * c+ __loc ) )->second;
584 uint16_type localDofInElement(
size_type __id, uint16_type __loc, uint16_type c = 0 )
const
590 friend struct element_access;
595 face_access( DofTable
const& __d )
601 return M_d.M_face_l2g.find( __id)->second.find(M_d.nLocalDofOnFace(
true )*c+__loc)->second;
604 uint16_type localDofInElement(
size_type __id, uint16_type __loc, uint16_type c = 0 )
const
606 return boost::get<3>( M_d.M_face_l2g.find( __id )->second.find( M_d.nLocalDofOnFace(
true )*c+__loc)->second );
611 friend struct face_access;
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
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;
625 return access_type( *
this )( El.id(), localNode, c );
628 template<
typename Elem>
630 localDofInElement( Elem
const& El,
const uint16_type localNode,
const uint16_type c = 0 )
const
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;
634 return ( access_type( *
this ) ).localDofInElement( El.id(), localNode, c );
650 return fe_type::numVertices;
658 return fe_type::numEdges;
666 return fe_type::numFaces;
678 for (
size_type __i = 0; __i < M_face_l2g.nrows(); ++__i )
680 for (
size_type __l = 0; __l < M_face_l2g.ncols(); ++__l )
682 std::cout <<
"face " << __i <<
" local " << __l
683 <<
" to global " << M_face_l2g[ __i][ __l ] <<
"\n";
695 M_elt_done[elt] =
true;
697 bool isElementDone(
size_type elt,
int c = 0 )
700 BOOST_FOREACH( LocalDof
const& local_dof, this->localDofSet( elt ) )
702 auto it = M_el_l2g.left.find( local_dof );
703 if ( it == M_el_l2g.left.end() )
708 void addDofFromElement( element_type
const& __elt,
716 if ( !this->isElementDone( __elt.id() ) )
718 DVLOG(3) <<
"adding dof from element " << __elt.id() <<
"\n";
720 DVLOG(3) <<
"next_free_dof " << next_free_dof <<
"\n";
721 DVLOG(3) <<
"current dof " <<
dofIndex( next_free_dof ) <<
"\n";
727 if ( is_continuous || is_discontinuous_locally )
733 uint16_type ldofcount = 0;
736 boost::tuple<uint16_type&,size_type&> shifts = boost::make_tuple( boost::ref( ldofcount ),
737 boost::ref( gdofcount ) );
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 );
756 const int ncdof = is_product?nComponents:1;
758 for ( uint16_type l = 0; l < nldof; ++l )
760 for (
int c = 0; c < ncdof; ++c, ++next_free_dof )
762 M_el_l2g.insert( dof_relation( LocalDof( ie, fe_type::nLocalDof*c + l ),
763 Dof( (
dofIndex( next_free_dof ) ) , 1,
false ) ) );
771 DVLOG(3) <<
"element " << __elt.id() <<
"has already been taken care of\n";
777 std::map<size_type, std::pair<size_type, size_type> >
780 return M_local2global;
788 return M_local2global.find( l )->second.first;
802 face_type
const& __face,
804 std::map<size_type,periodic_dof_map_type>& periodic_dof,
807 addVertexPeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag, mpl::bool_<(fe_type::nDofPerVertex>0)>() );
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> );
812 void addEdgePeriodicDof( element_type
const& __elt,
813 face_type
const& __face,
815 std::map<size_type,periodic_dof_map_type>& periodic_dof,
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>() );
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> );
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 )
831 static const bool cond = fe_type::nDofPerFace > 0;
832 addFacePeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag, mpl::bool_<cond>() );
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> );
853 void build( boost::shared_ptr<mesh_type>& M )
863 M_mesh = boost::addressof( M );
865 M_elt_done.resize( M.numElements() );
866 std::fill( M_elt_done.begin(), M_elt_done.end(), false );
868 VLOG(2) <<
"[Dof::build] initDofMap\n";
871 VLOG(2) <<
"[Dof::build] start building dof map\n";
873 VLOG(2) <<
"[Dof::build] start_next_free_dof = " << start_next_free_dof <<
"\n";
877 VLOG(2) <<
"[build] call buildPeriodicDofMap()\n";
879 VLOG(2) <<
"[Dof::build] start_next_free_dof(after periodic) = " << start_next_free_dof <<
"\n";
882 if ( is_discontinuous_locally )
884 VLOG(2) <<
"[build] call buildLocallyDiscontinuousDofMap()\n";
886 VLOG(2) <<
"[Dof::build] start_next_free_dof(after local discontinuities) = " << start_next_free_dof <<
"\n";
889 VLOG(2) <<
"[build] call buildDofMap()\n";
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;
899 for ( ; fit != fen; ++fit )
901 const int ncdof = is_product?nComponents:1;
903 for ( uint16_type c = 0; c < ncdof; ++c )
904 if ( !this->isElementDone( fit->id(), c ) )
906 em.push_back( boost::make_tuple( fit->id(), c, fit->marker().value() ) );
912 VLOG(2) <<
"[build] some element dof were not assigned\n";
914 for (
size_type i = 0; i < em.size(); ++i )
916 VLOG(3) <<
" - element " << boost::get<0>( em[i] ) <<
" c=" << boost::get<1>( em[i] )
917 <<
" m=" << boost::get<2>( em[i] ) <<
"\n";
923 VLOG(2) <<
"[build] check that all elements dof were assigned: OK\n";
929 VLOG(2) <<
"[build] call buildBoundaryDofMap()\n";
938 bool isP0continuous =
true;
939 for ( ; it_nlocDof!=en_nlocDof && isP0continuous ; ++it_nlocDof )
940 isP0continuous = isP0continuous && (*it_nlocDof==1);
942 if ( !isP0continuous )
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";
949 std::cerr <<
"ERROR : FEELPP_ENABLE_MPI_MODE is OFF" << std::endl;
955 for (
int proc=0; proc<this->
worldComm().globalSize(); ++proc )
989 #if defined(FEELPP_ENABLE_MPI_MODE)
1002 VLOG(2) <<
"[Dof::build] done building the map\n";
1030 #if defined(FEELPP_ENABLE_MPI_MODE)
1034 void buildGhostDofMap(
mesh_type& mesh );
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();
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 );
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 );
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 );
1096 typename dof_marker_type::right_range_type
1097 markerToDof( boost::any
const& marker )
1099 using namespace boost::bimaps;
1100 int id = M_mesh->markerId( marker );
1101 return M_dof_marker.right.range(
id <= _key, _key<
id+1 );
1104 typename dof_marker_type::right_range_type
1105 markerToDofLessThan( boost::any
const& marker )
1107 using namespace boost::bimaps;
1108 int id = M_mesh->markerId( marker );
1109 return M_dof_marker.right.range( unbounded, _key<
id );
1111 typename dof_marker_type::right_range_type
1112 markerToDofGreaterThan( boost::any
const& marker )
1114 using namespace boost::bimaps;
1115 int id = M_mesh->markerId( marker );
1116 return M_dof_marker.right.range(
id<_key, unbounded );
1119 void printDofMarker(std::string
const& filename )
1121 std::ofstream ofs( filename.c_str() );
1122 BOOST_FOREACH(
auto dof, M_dof_marker.left )
1124 ofs << dof.first <<
" " << dof.second <<
"\n";
1146 uint16_type processor,
1148 int32_type sign = 1,
1149 bool is_dof_periodic =
false,
1151 Marker1
const& marker = Marker1( 0 ) )
1154 const int ncdof = is_product?nComponents:1;
1156 for (
int c = 0; c < ncdof; ++c )
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;
1165 if ( itdof == endof )
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 ) ) );
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" );
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";
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";
1193 auto eit = M_el_l2g.left.find( LocalDof(ie,lc_dof) );
1195 if ( eit == M_el_l2g.left.end() )
1197 DCHECK( itdof->first == gDof ) <<
"very bad logical error in insertDof";
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>();
1206 M_dof_procset[ itdof->second+shift ].insert( processor );
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() ) );
1215 #if !defined(NDEBUG)
1216 M_dof2elt[itdof->second+shift].push_back( boost::make_tuple( ie, lc_dof, lc, itdof->first.template get<0>() ) );
1219 M_dof_view.insert(
Dof( itdof->second+shift,
1221 itdof->first.template get<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";
1234 size_type _dof = M_el_l2g.left.find(LocalDof(ie,lc_dof))->second.index();
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";
1244 res = res && ( __inserted || ( ( M_el_l2g.left.find(LocalDof(ie,lc_dof)) != M_el_l2g.left.end() ) && shift ) );
1255 M_dof_points.clear();
1256 this->generateDofPoints(M);
1262 std::pair<std::map<size_type,size_type>,std::map<size_type,size_type> >
1266 void addVertexDof( element_type
const& __elt, uint16_type processor,
size_type& next_free_dof,
1267 ref_shift_type& shifts )
1269 addVertexDof( __elt, processor, next_free_dof, shifts, mpl::bool_<(fe_type::nDofPerVertex>0)>() );
1271 void addVertexDof( element_type
const& , uint16_type ,
size_type& ,
1272 ref_shift_type& , mpl::bool_<false> )
1274 void addVertexDof( element_type
const& __elt, uint16_type processor,
size_type& next_free_dof,
1275 ref_shift_type& shifts, mpl::bool_<true> )
1277 uint16_type local_shift;
1279 boost::tie( local_shift, global_shift ) = shifts;
1284 uint16_type lc = local_shift;
1286 for ( uint16_type i = 0; i < element_type::numVertices; ++i )
1288 for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l, ++lc )
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() );
1298 shifts.template get<0>() = lc;
1300 #if !defined(NDEBUG)
1301 DVLOG(4) <<
"[Dof::updateVolumeDof(addVertexDof] vertex proc" << processor <<
" next_free_dof = " << next_free_dof <<
"\n";
1304 void addEdgeDof( element_type
const& __elt, uint16_type processor,
size_type& next_free_dof,
1305 ref_shift_type& shifts )
1307 static const bool cond = fe_type::nDofPerEdge > 0;
1308 return addEdgeDof( __elt,
1312 mpl::int_<fe_type::nDim>(),
1313 mpl::bool_<cond>() );
1315 void addEdgeDof( element_type
const& , uint16_type ,
size_type& ,
1316 ref_shift_type& , mpl::int_<1>, mpl::bool_<false> )
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> )
1322 uint16_type local_shift;
1324 boost::tie( local_shift, global_shift ) = shifts;
1327 uint16_type lc = local_shift;
1329 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
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() );
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";
1341 void addEdgeDof( element_type
const& , uint16_type ,
size_type& ,
1342 ref_shift_type& , mpl::int_<2>, mpl::bool_<false> )
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> )
1347 uint16_type local_shift;
1349 boost::tie( local_shift, global_shift ) = shifts;
1352 uint16_type lc = local_shift;
1356 for ( uint16_type i = 0; i < element_type::numEdges; ++i )
1358 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
1360 size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;
1361 int32_type sign = 1;
1364 if ( __elt.edgePermutation( i ).value() == edge_permutation_type::IDENTITY )
1369 else if ( __elt.edgePermutation( i ).value() == edge_permutation_type::REVERSE_PERMUTATION )
1372 if ( fe_type::is_modal )
1375 sign = ( l%2 )?( -1 ):( 1 );
1380 gDof += fe_type::nDofPerEdge - 1 - l ;
1384 FEELPP_ASSERT( 0 ).error (
"invalid edge permutation" );
1386 this->
insertDof( ie, lc, i, boost::make_tuple( 1, 0, gDof ), processor, next_free_dof, sign,
false, global_shift, __elt.edge( i ).marker() );
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";
1397 void addEdgeDof( element_type
const& , uint16_type ,
size_type& ,
1398 ref_shift_type& , mpl::int_<3>, mpl::bool_<false> )
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> )
1404 uint16_type local_shift;
1406 boost::tie( local_shift, global_shift ) = shifts;
1409 uint16_type lc = local_shift;
1411 for ( uint16_type i = 0; i < element_type::numEdges; ++i )
1413 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
1415 size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;
1417 int32_type sign = 1;
1419 if ( __elt.edgePermutation( i ).value() == edge_permutation_type::IDENTITY )
1424 else if ( __elt.edgePermutation( i ).value() == edge_permutation_type::REVERSE_PERMUTATION )
1427 if ( fe_type::is_modal )
1430 sign = ( l%2 )?( -1 ):( 1 );
1435 gDof += fe_type::nDofPerEdge - 1 - l ;
1439 FEELPP_ASSERT( 0 ).error (
"invalid edge permutation" );
1441 this->
insertDof( ie, lc, i, boost::make_tuple( 1, 0, gDof ), processor, next_free_dof, sign,
false, global_shift, __elt.edge( i ).marker() );
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";
1453 void addFaceDof( element_type
const& __elt, uint16_type processor,
size_type& next_free_dof,
1454 ref_shift_type& shifts )
1456 return addFaceDof( __elt, processor, next_free_dof, shifts, mpl::int_<fe_type::nDim>(), mpl::bool_<(fe_type::nDofPerFace > 0)>() );
1458 void addFaceDof( element_type
const& , uint16_type ,
size_type& ,
1459 ref_shift_type& , mpl::int_<1>, mpl::bool_<false> )
1461 void addFaceDof( element_type
const& , uint16_type ,
size_type& ,
1462 ref_shift_type& , mpl::int_<2>, mpl::bool_<false> )
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> )
1467 uint16_type local_shift;
1469 boost::tie( local_shift, global_shift ) = shifts;
1472 uint16_type lc = local_shift;
1474 for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
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() );
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";
1486 void addFaceDof( element_type
const& , uint16_type ,
size_type& ,
1487 ref_shift_type& , mpl::int_<3>, mpl::bool_<false> )
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> )
1492 uint16_type local_shift;
1494 boost::tie( local_shift, global_shift ) = shifts;
1498 uint16_type lc = local_shift;
1500 for ( uint16_type i = 0; i < element_type::numFaces; ++i )
1502 face_permutation_type permutation = __elt.facePermutation( i );
1503 FEELPP_ASSERT( permutation != face_permutation_type( 0 ) ).error (
"invalid face permutation" );
1510 int MaxOrder = int( ( 3 + std::sqrt( 1+8*fe_type::nDofPerFace ) )/2 ) - 2;
1512 for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
1520 size_type gDof = __elt.face( i ).id() * fe_type::nDofPerFace;
1521 int32_type sign = 1;
1529 MaxOrder = MaxOrder-1;
1532 if ( !fe_type::is_modal )
1535 if ( permutation == face_permutation_type( 1 ) || fe_type::nDofPerFace == 1 )
1539 gDof += vector_permutation[permutation][l];
1546 if ( permutation == face_permutation_type( 2 ) )
1557 this->
insertDof( ie, lc, i, boost::make_tuple( 2, 0, gDof ), processor, next_free_dof, sign,
false, global_shift,__elt.face( i ).marker() );
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";
1568 void addVolumeDof( element_type
const& __elt, uint16_type processor,
size_type& next_free_dof,
1569 ref_shift_type& shifts )
1571 return addVolumeDof( __elt, processor, next_free_dof, shifts, mpl::bool_<(fe_type::nDofPerVolume>0)>() );
1573 void addVolumeDof( element_type
const& , uint16_type ,
size_type& ,
1574 ref_shift_type& , mpl::bool_<false> )
1576 void addVolumeDof( element_type
const& __elt, uint16_type processor,
size_type& next_free_dof,
1577 ref_shift_type& shifts, mpl::bool_<true> )
1579 BOOST_STATIC_ASSERT( element_type::numVolumes );
1580 uint16_type local_shift;
1582 boost::tie( local_shift, global_shift ) = shifts;
1585 uint16_type lc = local_shift;
1587 for ( uint16_type l = 0; l < fe_type::nDofPerVolume; ++l, ++lc )
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() );
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";
1600 template<
typename FaceIterator>
1601 void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc )
1603 addVertexBoundaryDof( __face_it, lc, mpl::bool_<(fe_type::nDofPerVertex>0)>(), mpl::int_<nDim>() );
1605 template<
typename FaceIterator>
void addVertexBoundaryDof( FaceIterator , uint16_type& , mpl::bool_<false>, mpl::int_<1> ) {}
1606 template<
typename FaceIterator>
void addVertexBoundaryDof( FaceIterator , uint16_type& , mpl::bool_<false>, mpl::int_<2> ) {}
1607 template<
typename FaceIterator>
void addVertexBoundaryDof( FaceIterator , uint16_type& , 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> )
1611 BOOST_STATIC_ASSERT( face_type::numVertices );
1612 #if !defined(FEELPP_ENABLE_MPI_MODE)
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" );
1619 uint16_type iFaEl = __face_it->pos_first();
1625 if ( __face_it->processId() == __face_it->proc_first() )
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" );
1631 iFaEl = __face_it->pos_first();
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" );
1641 iFaEl = __face_it->pos_second();
1647 const int ncdof = is_product?nComponents:1;
1649 for (
int c = 0; c < ncdof; ++c )
1651 for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
1653 auto temp= this->localToGlobal( iElAd,
1654 iFaEl * fe_type::nDofPerVertex + l,
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 );
1661 template<
typename FaceIterator>
1662 void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<2> )
1664 BOOST_STATIC_ASSERT( face_type::numVertices );
1665 #if !defined(FEELPP_ENABLE_MPI_MODE)
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" );
1672 uint16_type iFaEl = __face_it->pos_first();
1678 if ( __face_it->processId() == __face_it->proc_first() )
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" );
1684 iFaEl = __face_it->pos_first();
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" );
1694 iFaEl = __face_it->pos_second();
1700 size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
1701 face_type::numEdges * fe_type::nDofPerEdge +
1702 face_type::numFaces * fe_type::nDofPerFace );
1707 const int ncdof = is_product?nComponents:1;
1709 for (
int c = 0; c < ncdof; ++c )
1711 uint16_type lcc=c*ndofF;
1713 for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
1716 uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
1721 for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
1723 auto temp = this->localToGlobal( iElAd,
1724 iVeEl * fe_type::nDofPerVertex + l,
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 );
1732 template<
typename FaceIterator>
1733 void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<3> )
1735 addVertexBoundaryDof( __face_it, lc, mpl::bool_<true>(), mpl::int_<2>() );
1737 template<
typename FaceIterator>
1738 void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc )
1740 static const bool cond = fe_type::nDofPerEdge*face_type::numEdges > 0;
1741 addEdgeBoundaryDof( __face_it, lc, mpl::bool_<cond>(), mpl::int_<nDim>() );
1743 template<
typename FaceIterator>
1744 void addEdgeBoundaryDof( FaceIterator , uint16_type& , mpl::bool_<false>, mpl::int_<1> ) {}
1745 template<
typename FaceIterator>
1746 void addEdgeBoundaryDof( FaceIterator , uint16_type& , mpl::bool_<false>, mpl::int_<2> ) {}
1747 template<
typename FaceIterator>
1748 void addEdgeBoundaryDof( FaceIterator , uint16_type& , 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> )
1752 #if !defined(FEELPP_ENABLE_MPI_MODE)
1755 size_type iElAd = __face_it->ad_first();
1757 ( __face_it->id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
1760 uint16_type iFaEl = __face_it->pos_first();
1765 if ( __face_it->processId() == __face_it->proc_first() )
1767 iElAd = __face_it->ad_first();
1769 ( __face_it->id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
1772 iFaEl = __face_it->pos_first();
1777 iElAd = __face_it->ad_second();
1779 ( __face_it->id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
1782 iFaEl = __face_it->pos_second();
1789 #if !defined(NDEBUG)
1790 DVLOG(4) <<
" local face id : " << iFaEl <<
"\n";
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 );
1798 const int ncdof = is_product?nComponents:1;
1800 for (
int c = 0; c < ncdof; ++c )
1802 uint16_type lcc=nVerticesF+c*ndofF;
1805 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
1807 auto temp = this->localToGlobal( iElAd,
1808 element_type::numVertices*fe_type::nDofPerVertex +
1809 iFaEl * fe_type::nDofPerEdge + l,
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 );
1817 template<
typename FaceIterator>
1818 void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<3> )
1821 #if !defined(FEELPP_ENABLE_MPI_MODE)
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" );
1828 uint16_type iFaEl = __face_it->pos_first();
1834 if ( __face_it->processId() == __face_it->proc_first() )
1836 iElAd = __face_it->ad_first();
1838 ( __face_it->id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
1841 iFaEl = __face_it->pos_first();
1846 iElAd = __face_it->ad_second();
1848 ( __face_it->id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
1851 iFaEl = __face_it->pos_second();
1856 #if !defined(NDEBUG)
1857 DVLOG(4) <<
" local face id : " << iFaEl <<
"\n";
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 );
1864 const int ncdof = is_product?nComponents:1;
1866 for (
int c = 0; c < ncdof; ++c )
1868 uint16_type lcc=nVerticesF+c*ndofF;
1871 for ( uint16_type iEdFa = 0; iEdFa < face_type::numEdges; ++iEdFa )
1874 uint16_type iEdEl = element_type::fToE( iFaEl, iEdFa );
1879 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
1881 auto temp = this->localToGlobal( iElAd,
1882 element_type::numVertices*fe_type::nDofPerVertex +
1883 iEdEl * fe_type::nDofPerEdge + l,
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 );
1894 template<
typename FaceIterator>
1895 void addFaceBoundaryDof( FaceIterator __face_it, uint16_type& lc )
1897 addFaceBoundaryDof( __face_it, lc, mpl::bool_<(face_type::numFaces*fe_type::nDofPerFace > 0)>() );
1899 template<
typename FaceIterator>
1900 void addFaceBoundaryDof( FaceIterator , uint16_type& , mpl::bool_<false> )
1903 template<
typename FaceIterator>
1904 void addFaceBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true> )
1906 #if !defined(FEELPP_ENABLE_MPI_MODE)
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" );
1913 uint16_type iFaEl = __face_it->pos_first();
1919 if ( __face_it->processId() == __face_it->proc_first() )
1921 iElAd = __face_it->ad_first();
1923 ( __face_it->id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
1926 iFaEl = __face_it->pos_first();
1931 iElAd = __face_it->ad_second();
1933 ( __face_it->id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
1936 iFaEl = __face_it->pos_second();
1941 #if !defined(NDEBUG)
1942 DVLOG(4) <<
" local face id : " << iFaEl <<
"\n";
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 );
1950 const int ncdof = is_product?nComponents:1;
1952 for (
int c = 0; c < ncdof; ++c )
1954 uint16_type lcc=nVerticesAndEdgeF+c*ndofF;
1957 for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l )
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,
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 );
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> );
1989 Feel::detail::ignore_unused_variable_warning( M );
1994 #if !defined(NDEBUG)
1996 using namespace Feel;
1998 gm_ptrtype M_gm_ptr = M.gm();
2003 value_type tol = value_type( 10.0 )*type_traits<double>::epsilon();
2005 bool global_signs_good = 1;
2007 std::vector<size_type> bad_dof;
2009 for ( uint16_type gDof = 0; gDof <
nDof(); ++gDof )
2011 uint16_type _numEntities = M_dof2elt[gDof].size();
2012 uint16_type _ent = M_dof2elt[gDof].begin()->template get<3>();
2014 if ( _numEntities > 1 && _ent == mpl::int_<N>() )
2016 bool signs_good = 1;
2018 std::vector< ublas::vector<value_type> > basis_eval;
2020 std::vector< points_type > real_coordinates;
2022 ldof_const_iterator __ldofit = M_dof2elt[gDof].begin();
2023 ldof_const_iterator __ldofen = M_dof2elt[gDof].end();
2025 while ( __ldofit != __ldofen )
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>();
2031 PointSetMapped<element_type, convex_type, nOrder> test( M.element( entity_element_id ) );
2033 points_type Gt = test.pointsBySubEntity( N, entity_local_id );
2035 real_coordinates.push_back( test.pointsBySubEntity( N, entity_local_id, 0, 1 ) );
2037 int sign = boost::get<1>( localToGlobal( entity_element_id, entity_local_dof_id ) );
2039 basis_eval.push_back( value_type( sign )*ublas::row( M_basis.evaluate( Gt ), entity_local_dof_id ) );
2044 for ( uint16_type i=1; i < _numEntities; i++ )
2047 FEELPP_ASSERT( ublas::norm_inf( real_coordinates[i] - real_coordinates[0] ) < tol )
2049 ( real_coordinates[0] )
2050 ( real_coordinates[i] ).error(
"Reference points aren't being mapped to the same real one's" );
2052 if ( ublas::norm_inf( basis_eval[i] - basis_eval[0] ) > tol )
2055 global_signs_good = 0;
2060 real_coordinates.empty();
2062 if ( signs_good == 0 )
2063 bad_dof.push_back( gDof );
2067 if ( !bad_dof.empty() )
2069 for ( uint16_type i = 0; i < bad_dof.size(); ++i )
2070 Warning() << bad_dof[i] <<
"\n";
2072 if ( mpl::int_<N>() == 1 )
2073 Warning() <<
"Edges: ";
2076 Warning() <<
"Faces: ";
2078 Warning() <<
"Bad dof signs. \n";
2084 void checkDofContinuity(
mesh_type& , mpl::int_<1> ) {}
2086 void checkDofContinuity(
mesh_type& mesh, mpl::int_<2> )
2088 checkDofEntity<1>( mesh );
2091 void checkDofContinuity(
mesh_type& mesh, mpl::int_<3> )
2093 checkDofContinuity( mesh, mpl::int_<2>() );
2094 checkDofEntity<2>( mesh );
2097 void generateFacePermutations (
mesh_type& , mpl::bool_<false> ) {}
2099 void generateFacePermutations (
mesh_type& mesh, mpl::bool_<true> )
2101 element_type _elt = *mesh.beginElement();
2102 PointSetMapped<element_type, convex_type, nOrder> pts( _elt );
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 ) );
2108 void generatePeriodicDofPoints(
mesh_type& M, periodic_element_list_type
const& periodic_elements, dof_points_type& periodic_dof_points );
2117 reference_convex_type M_convex_ref;
2120 std::vector<bool> M_elt_done;
2122 uint16_type M_n_dof_per_face_on_bdy;
2123 uint16_type M_n_dof_per_face;
2126 Container_fromface M_face_l2g;
2128 mutable local_dof_set_type M_local_dof_set;
2130 dof_marker_type M_dof_marker;
2132 dof_map_type map_gdof;
2134 std::map<face_permutation_type, permutation_vector_type> vector_permutation;
2139 std::map<size_type, std::pair<size_type, size_type> > M_local2global;
2147 dof_procset_type M_dof_procset;
2152 dof_points_type M_dof_points;
2154 std::vector<Dof> M_dof_indices;
2156 periodicity_type M_periodicity;
2158 periodic_element_list_type periodic_elements;
2163 typedef typename std::vector<localglobal_indices_type,Eigen::aligned_allocator<localglobal_indices_type> > vector_indices_type;
2166 vector_indices_type M_locglob_indices;
2167 vector_indices_type M_locglob_signs;
2170 vector_indices_type M_locglobOnCluster_indices;
2171 vector_indices_type M_locglobOnCluster_signs;
2174 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
2177 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2180 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2183 super( _worldComm ),
2194 M_periodicity( periodicity )
2196 VLOG(2) <<
"[dof] is_periodic = " << is_periodic <<
"\n";
2206 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2209 super( _worldComm ),
2216 M_local_dof_set( 0, nLocalDof() ),
2220 M_periodicity( periodicity )
2224 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
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 )
2242 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2246 LOG(INFO) <<
" Degree of Freedom (DofTable) Object" <<
"\n";
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";
2255 for (
size_type i = 0; i < M_n_el; ++i )
2258 for (
size_type j = 0; j < nDofPerElement; ++j )
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";
2271 LOG(INFO) <<
"************************************************************" <<
"\n";
2272 LOG(INFO) <<
" Boundary Local to Global DOF table" <<
"\n";
2273 LOG(INFO) <<
"************************************************************" <<
"\n";
2275 auto it = M_face_l2g.begin();
2276 auto en = M_face_l2g.end();
2278 for (
size_type f = 0; it!=en; ++it,++f )
2280 std::ostringstream ostr;
2281 ostr <<
"face id " << it->first <<
" : ";
2283 auto it2 = it->second.begin();
2284 auto en2 = it->second.end();
2286 for (
size_type l = 0; it2!=en2; ++it2,++l )
2288 ostr <<
"(local/global/sign dof : " << it2->first <<
" : "
2289 << boost::get<0>( it2->second )<<
" : "
2290 << boost::get<1>( it2->second ) <<
"\n";
2293 LOG(INFO) << ostr.str() <<
"\n";
2299 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2302 face_type
const& __face,
2304 std::map<size_type,periodic_dof_map_type>& periodic_dof,
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 );
2317 uint16_type iFaEl = __face.pos_first();
2321 for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
2324 uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
2325 Feel::detail::ignore_unused_variable_warning( iVeEl );
2330 for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
2332 uint16_type lid = iVeEl * fe_type::nDofPerVertex + l;
2334 const size_type gDof = ( __elt.point( iVeEl ).id() ) * fe_type::nDofPerVertex + l;
2336 VLOG(2) <<
"add vertex periodic doc " << next_free_dof <<
" in element " << __elt.id() <<
" lid = " << lid <<
"\n";
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";
2342 const int ncdof = is_product?nComponents:1;
2343 for (
int c1 = 0; c1 < ncdof; ++c1 )
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 ) ) );
2351 VLOG(2) <<
"added vertex periodic dof " << __elt.id() <<
", " << lid <<
", " << boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) ) <<
"\n";
2359 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2362 face_type
const& __face,
2364 std::map<size_type,periodic_dof_map_type>& periodic_dof,
2374 ( __face.id() ).error(
"[DofTable::buildBoundaryDof] invalid face/element in face" );
2378 uint16_type iFaEl = __face.pos_first();
2380 #if !defined(NDEBUG)
2381 DVLOG(4) <<
" local face id : " << iFaEl <<
"\n";
2385 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
2387 uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
2389 const size_type gDof = ( __elt.edge( iFaEl ).id() ) * fe_type::nDofPerEdge + l;
2391 DVLOG(4) <<
"add edge periodic dof " << next_free_dof <<
" in element " << __elt.id() <<
" lid = " << lid <<
"\n";
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";
2397 const int ncdof = is_product?nComponents:1;
2398 for (
int c1 = 0; c1 < ncdof; ++c1 )
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 ) ) );
2406 DVLOG(4) <<
"added edge periodic dof " << __elt.id() <<
", " << lid <<
", " << boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) ) <<
"\n";
2411 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2414 face_type
const& __face,
2416 std::map<size_type,periodic_dof_map_type>& periodic_dof,
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" );
2430 uint16_type iFaEl = __face.pos_first();
2432 #if !defined(NDEBUG)
2433 DVLOG(4) <<
" local face id : " << iFaEl <<
"\n";
2437 for ( uint16_type iEdFa = 0; iEdFa < face_type::numEdges; ++iEdFa )
2440 uint16_type iEdEl = element_type::fToE( iFaEl, iEdFa );
2441 Feel::detail::ignore_unused_variable_warning( iEdEl );
2445 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
2448 uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
2450 size_type gDof = __elt.edge( l ).id() * fe_type::nDofPerEdge;
2452 if ( __elt.edgePermutation( l ).value() == edge_permutation_type::IDENTITY )
2457 else if ( __elt.edgePermutation( l ).value() == edge_permutation_type::REVERSE_PERMUTATION )
2459 gDof += fe_type::nDofPerEdge - 1 - l ;
2462 DVLOG(4) <<
"add periodic doc " << next_free_dof <<
" in element " << __elt.id() <<
" lid = " << lid <<
"\n";
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";
2468 const int ncdof = is_product?nComponents:1;
2469 for (
int c1 = 0; c1 < ncdof; ++c1 )
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 ) ) );
2477 DVLOG(4) <<
"added " << __elt.id() <<
", " << lid <<
", " << boost::get<0>( localToGlobal( __elt.id(), lid, c1 ) ) <<
"\n";
2483 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2486 face_type
const& __face,
2488 std::map<size_type,periodic_dof_map_type>& periodic_dof,
2496 FEELPP_ASSERT( iElAd !=
invalid_size_type_value )( __face.id() ).error(
"[Dof::buildBoundaryDof] invalid face/element in face" );
2499 uint16_type iFaEl = __face.pos_first();
2501 #if !defined(NDEBUG)
2502 DVLOG(2) <<
" local face id : " << iFaEl <<
"\n";
2506 for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l )
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,
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 );
2521 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2525 M_n_el = M.numElements();
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;
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";
2547 FEELPP_ASSERT( nldof == fe_type::nLocalDof )
2549 ( fe_type::nLocalDof ).error(
"Something wrong in FE specification" ) ;
2555 int ntldof = is_product?nComponents*nldof:nldof;
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 );
2563 M_face_sign = ublas::scalar_vector<bool>( M.numFaces(), false );
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>() );
2569 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
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;
2579 FEELPP_ASSERT( nldof == fe_type::nLocalDof )
2581 ( fe_type::nLocalDof ).error(
"Something wrong in FE specification" ) ;
2583 const size_type n_proc = M.worldComm().localSize();
2586 for (
size_type processor=0; processor<n_proc; 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";
2595 it_elt = M.beginElementWithProcessId( processor );
2597 it_elt = M.beginElementWithProcessId( processor );
2598 VLOG(2) <<
"[buildDofMap] starting with elt " << it_elt->id() <<
"\n";
2600 for ( ; it_elt!=en_elt; ++it_elt )
2602 element_type
const& __elt = *it_elt;
2606 typename element_type::face_const_iterator it, en;
2607 boost::tie( it, en ) = it_elt->faces();
2610 for ( ; it != en; ++it )
2612 if ( ( *it )->marker().value() == M_periodicity.tag2() ||
2613 ( *it )->marker().value() == M_periodicity.tag1() )
2618 periodic_elements.push_back( boost::make_tuple( boost::addressof( __elt ), *it ) );
2626 VLOG(2) <<
"[buildPeriodicDofMap] built periodic_elements " << periodic_elements.size() <<
"\n";
2627 std::map<size_type,periodic_dof_map_type> periodic_dof;
2632 periodic_element_list_iterator it_periodic = periodic_elements.begin();
2633 periodic_element_list_iterator en_periodic = periodic_elements.end();
2636 while ( it_periodic != en_periodic )
2638 element_type
const& __elt = *it_periodic->template get<0>();
2639 face_type
const& __face = *it_periodic->template get<1>();
2641 if ( __face.marker().value() == M_periodicity.tag1() )
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() );
2651 it_periodic = periodic_elements.begin();
2653 while ( it_periodic != en_periodic )
2655 element_type
const& __elt = *it_periodic->template get<0>();
2656 face_type
const& __face = *it_periodic->template get<1>();
2658 if ( __face.marker().value() == M_periodicity.tag2() )
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() );
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";
2671 dof_points_type periodic_dof_points( next_free_dof );
2672 generatePeriodicDofPoints( M, periodic_elements, periodic_dof_points );
2674 VLOG(2) <<
"[periodic dof table] generated dof points\n";
2675 VLOG(2) <<
"[periodic dof table] start matching the dof points\n";
2678 std::pair<size_type,periodic_dof_type> dof;
2679 BOOST_FOREACH( dof, periodic_dof[M_periodicity.tag1()] )
2682 max_gid = ( max_gid > gid )?max_gid:gid;
2686 std::pair<size_type,periodic_dof_type> dof2;
2687 BOOST_FOREACH( dof2, periodic_dof[M_periodicity.tag2()] )
2690 FEELPP_ASSERT( gid2 > max_gid )( gid2 )( max_gid ).error(
"invalid dof index" );
2691 max_gid2 = ( max_gid2 > gid2 )?max_gid2:gid2;
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";
2700 std::vector<bool> periodic_dof_done( max_gid+1 );
2701 std::fill( periodic_dof_done.begin(), periodic_dof_done.end(), false );
2703 BOOST_FOREACH( dof, periodic_dof[M_periodicity.tag1()] )
2707 if ( periodic_dof_done[gid] )
2710 node_type x1 = periodic_dof_points[gid].template get<0>();
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();
2715 for ( ; it_dof2 != en_dof2; ++ it_dof2 )
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>();
2724 it_dof2 = periodic_dof[M_periodicity.tag2()].begin();
2727 for ( ; it_dof2 != en_dof2; ++ it_dof2 )
2731 if ( it_dof2->second.template get<2>() != dof.second.template get<2>() )
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" );
2742 if ( ublas::norm_2( x1-( x2-M_periodicity.translation() ) ) < 1e-10 )
2746 corresponding_gid = gid2;
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>();
2762 VLOG(2) <<
"matching dof id " << gid <<
" with dof id=" << corresponding_gid <<
"\n";
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";
2768 while ( it_dof2 != en_dof2 )
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>();
2779 FEELPP_ASSERT( dof1_type == dof2_type )
2780 ( gid )( it_dof2->first )( gDof )( lid )( c2) ( ie )
2781 ( dof1_type )( dof2_type ).error (
"invalid dof" );
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";
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 ) ) );
2791 CHECK( successful_modify ) <<
"modify periodic dof table failed: element id "
2792 << ie <<
" local dof id " << lid <<
" component " << c2;
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
2801 <<
", map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ]= "
2802 << map_gdof[ boost::make_tuple( dof2_type, c2, gDof ) ] <<
"\n";
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;
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" );
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;
2822 VLOG(2) <<
"[periodic] invalid point/dof matching\n";
2823 VLOG(2) <<
"[periodic] n = " << x1 <<
"\n";
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";
2831 if ( !periodic_dof[M_periodicity.tag2()].empty() )
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();
2837 while ( it_dof2 != en_dof2 )
2840 size_type ie = it_dof2->second.template get<0>();
2841 size_type lid = it_dof2->second.template get<1>();
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";
2855 VLOG(2) <<
"[periodic] periodic condition done\n";
2861 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
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 ) );
2868 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
2872 if ( !M_dof_indices.empty() )
2874 generateDofPoints( M );
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;
2884 FEELPP_ASSERT( nldof == fe_type::nLocalDof )
2886 ( fe_type::nLocalDof ).error(
"Something wrong in FE specification" ) ;
2888 #if !defined(FEELPP_ENABLE_MPI_MODE) // sequential if (this->worldComm().size()==1)
2890 const size_type n_proc = M.worldComm().localSize();
2893 std::list<std::pair<element_type const*, face_type const*> > periodic_elements;
2895 size_type next_free_dof = start_next_free_dof;
2897 for (
size_type processor=0; processor<n_proc; 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;
2906 if ( is_periodic || is_discontinuous_locally )
2907 this->M_first_df[processor] = 0;
2909 it_elt = M.beginElementWithProcessId( processor );
2911 DVLOG(2) <<
"[buildDofMap] starting with elt " << it_elt->id() <<
"\n";
2913 for ( ; it_elt!=en_elt; ++it_elt )
2915 this->addDofFromElement( *it_elt, next_free_dof, processor );
2918 #if !defined( NDEBUG )
2919 const int ncdof = is_product?nComponents:1;
2921 for (
int c = 0; c < ncdof; ++c )
2923 it_elt = M.beginElementWithProcessId( processor );
2925 for ( ; it_elt!=en_elt; ++it_elt )
2927 element_type
const& __elt = *it_elt;
2928 std::ostringstream ostr;
2929 ostr <<
"element " << __elt.id() <<
" : ";
2931 for ( uint16_type l = 0; l < nldof; ++l )
2933 ostr << M_el_l2g.left.find( LocalDof( __elt.id(), fe_type::nLocalDof*c + l ) )->second.index() <<
" ";
2936 DVLOG(2) << ostr.str() <<
"\n";
2941 this->M_last_df[processor] = next_free_dof-1;
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];
2950 this->M_n_dofs = next_free_dof;
2952 #if !defined(NDEBUG)
2953 DVLOG(2) <<
" n global dof " << nDof() <<
"\n";
2954 DVLOG(2) <<
" n local dof " << nLocalDof() <<
"\n";
2957 for (
size_type processor=0; processor<M.worldComm().localSize(); processor++ )
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";
2968 for (
size_type processor=0; processor<n_proc; processor++ )
2970 auto it_elt = M.beginElementWithProcessId( processor );
2971 auto en_elt = M.endElementWithProcessId( processor );
2973 for ( ; it_elt != en_elt; ++it_elt )
2977 for (
int i = 0; i < FEType::nLocalDof; ++i )
2979 int nc1 = ( is_product?nComponents1:1 );
2981 for (
int c1 =0; c1 < nc1; ++c1 )
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 );
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 );
3001 size_type theFirstDf = start_next_free_dof;
3003 if ( is_periodic || is_discontinuous_locally )
3006 mpi::all_gather( M.worldComm().localComm(),
3013 size_type next_free_dof = start_next_free_dof;
3015 for ( ; it_elt!=en_elt; ++it_elt )
3017 this->addDofFromElement( *it_elt, next_free_dof, M.worldComm().localRank() );
3021 for (
auto mit = M_dof_marker.right.begin(), men = M_dof_marker.right.end() ; mit != men ; ++mit )
3023 LOG(INFO) <<
"marker " << mit->first <<
" dof id " << mit->second;
3028 LOG(INFO) <<
"local to global view";
3029 for(
auto it = M_el_l2g.left.begin(), en = M_el_l2g.left.end();
3032 LOG(INFO) <<
"local dof (" << it->first.elementId()<<
","<< it->first.localDof() <<
") --> global dof " << it->second.index();
3034 LOG(INFO) <<
"global to local view";
3035 for(
auto it = M_el_l2g.right.begin(), en = M_el_l2g.right.end();
3038 LOG(INFO) <<
"global dof " << it->first.index() <<
" --> local dof (" << it->second.elementId()<<
","<< it->second.localDof() <<
")";
3041 const size_type thelastDof = ( !hasNoElt )?next_free_dof-1:0;
3043 mpi::all_gather( M.worldComm().localComm(),
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()];
3053 mpi::all_gather( M.worldComm().localComm(),
3055 this->M_n_localWithGhost_df );
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()
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;
3069 it_elt = M.beginElementWithProcessId( M.worldComm().localRank() );
3071 for ( ; it_elt != en_elt; ++it_elt )
3075 for (
int i = 0; i < FEType::nLocalDof; ++i )
3077 int nc1 = ( is_product?nComponents1:1 );
3079 for (
int c1 =0; c1 < nc1; ++c1 )
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();
3090 generateDofPoints( M );
3094 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
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";
3105 if ( nDofF == 0 )
return;
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;
3114 typename mesh_type::face_const_iterator __face_it = M.beginFace();
3115 typename mesh_type::face_const_iterator __face_en = M.endFace();
3119 int ntldof = nLocalDofOnFace();
3121 DVLOG(2) <<
"[buildBoundaryDofMap] nb faces : " << nF <<
"\n";
3122 DVLOG(2) <<
"[buildBoundaryDofMap] nb dof faces : " << nDofF*nComponents <<
"\n";
3124 for (
size_type nf = 0; __face_it != __face_en; ++__face_it, ++nf )
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();
3132 if ( !__face_it->isConnectedTo0() )
continue;
3134 #if !defined(NDEBUG)
3136 if ( __face_it->isOnBoundary() )
3137 DVLOG(4) <<
"[buildBoundaryDofMap] boundary global face id : " << __face_it->id()
3138 <<
" marker: " << __face_it->marker()<<
"\n";
3141 DVLOG(4) <<
"[buildBoundaryDofMap] global face id : " << __face_it->id() <<
"\n";
3144 uint16_type lcVertex = 0;
3145 uint16_type lcEdge = 0;
3146 uint16_type lcFace = 0;
3149 addVertexBoundaryDof( __face_it, lcVertex );
3150 addEdgeBoundaryDof( __face_it, lcEdge );
3151 addFaceBoundaryDof( __face_it, lcFace );
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" );
3165 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3169 if ( !M_dof_points.empty() )
3172 if ( fe_type::is_modal )
3175 DVLOG(2) <<
"[Dof::generateDofPoints] generating dof coordinates\n";
3177 typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;
3181 gm_ptrtype gm(
new gm_type );
3187 typename gm_type::precompute_ptrtype __geopc(
new typename gm_type::precompute_type( gm, fe.points() ) );
3191 element_const_iterator it_elt = M.beginElementWithProcessId( M.worldComm().localRank() );
3192 element_const_iterator en_elt = M.endElementWithProcessId( M.worldComm().localRank() );
3194 if ( it_elt == en_elt )
3197 gm_context_ptrtype __c(
new gm_context_type( gm, *it_elt, __geopc ) );
3199 std::vector<bool> dof_done( nLocalDofWithGhost() );
3200 M_dof_points.resize( nLocalDofWithGhost() );
3201 std::fill( dof_done.begin(), dof_done.end(), false );
3203 for (
size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
3205 __c->update( *it_elt );
3207 for ( uint16_type l =0; l < fe_type::nLocalDof; ++l )
3209 int ncdof = is_product?nComponents:1;
3211 for ( uint16_type c1 = 0; c1 < ncdof; ++c1 )
3213 size_type thedof = boost::get<0>( localToGlobal( it_elt->id(), l, c1 ) );
3215 if ( ( thedof >= firstDof() ) && ( thedof <= lastDof() ) )
3219 thedof -= firstDof();
3220 DCHECK( thedof < nLocalDofWithGhost() )
3221 <<
"invalid local dof index "
3222 << thedof <<
", " << nLocalDofWithGhost() <<
"," << firstDof() <<
","
3223 << lastDof() <<
"," << it_elt->id() <<
"," << l <<
"," << c1;
3225 if ( dof_done[ thedof ] ==
false )
3227 std::set<uint16_type> lid;
3230 M_dof_points[thedof] = boost::make_tuple( __c->xReal( l ), firstDof()+thedof, c1 );
3231 dof_done[thedof] =
true;
3239 for (
size_type dof_id = 0; dof_id < nLocalDofWithGhost() ; ++dof_id )
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;
3253 DVLOG(2) <<
"[Dof::generateDofPoints] generating dof coordinates done\n";
3255 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3258 periodic_element_list_type
const& periodic_elements,
3259 dof_points_type& periodic_dof_points )
3261 if ( fe_type::is_modal )
3264 DVLOG(2) <<
"[Dof::generateDofPoints] generating dof coordinates\n";
3266 typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;
3270 gm_ptrtype gm(
new gm_type );
3276 typename gm_type::precompute_ptrtype __geopc(
new typename gm_type::precompute_type( gm, fe.points() ) );
3280 periodic_element_list_const_iterator it_elt = periodic_elements.begin();
3281 periodic_element_list_const_iterator en_elt = periodic_elements.end();
3283 if ( it_elt == en_elt )
3286 gm_context_ptrtype __c(
new gm_context_type( gm, *it_elt->template get<0>(), __geopc ) );
3288 std::vector<bool> dof_done( periodic_dof_points.size() );
3289 std::fill( dof_done.begin(), dof_done.end(), false );
3291 for (
size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
3293 __c->update( *it_elt->template get<0>() );
3295 face_type
const& __face = *it_elt->template get<1>();
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 );
3302 uint16_type iFaEl = __face.pos_first();
3305 int ncdof = is_product?nComponents:1;
3307 for ( uint16_type c1 = 0; c1 < ncdof; ++c1 )
3310 for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
3313 uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
3314 Feel::detail::ignore_unused_variable_warning( iVeEl );
3319 for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
3321 uint16_type lid = iVeEl * fe_type::nDofPerVertex + l;
3323 size_type thedof = boost::get<0>( localToGlobal( it_elt->template get<0>()->id(), lid, c1 ) );
3324 FEELPP_ASSERT( thedof < dof_done.size() )
3328 ( it_elt->template get<0>()->
id() )
3329 ( lid ).error (
"[generatePeriodicDofPoints] invalid dof id" );
3331 if ( dof_done[ thedof ] ==
false )
3333 periodic_dof_points[thedof] = boost::make_tuple( __c->xReal( lid ), thedof, c1 );
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" );
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" );
3345 dof_done[thedof] =
true;
3350 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
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() )
3358 ( it_elt->template get<0>()->
id() )
3359 ( lid ).error (
"[generatePeriodicDofPoints] invalid dof id" );
3361 if ( dof_done[ thedof ] ==
false )
3363 periodic_dof_points[thedof] = boost::make_tuple( __c->xReal( lid ), thedof, c1 );
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" );
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" );
3375 dof_done[thedof] =
true;
3382 for (
size_type dof_id = 0; dof_id < periodic_dof_points.size() ; ++dof_id )
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" );
3394 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3398 addSubstructuringDofVertex( M, next_free_dof );
3399 addSubstructuringDofEdge( M, next_free_dof, mpl::int_<nDim>() );
3400 addSubstructuringDofFace( M, next_free_dof, mpl::int_<nDim>() );
3403 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
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";
3411 std::cout <<
"n wb: " << std::distance( M.beginEdgeWithMarker( M.markerName(
"WireBasket") ), M.endEdgeWithMarker( M.markerName(
"WireBasket") ) ) <<
"\n";
3415 for(
auto pit = M.beginPointWithMarker( M.markerName(
"CrossPoints") ),
3416 pen = M.endPointWithMarker( M.markerName(
"CrossPoints") );
3420 auto __elt = M.element( *pit->elements().begin() );
3423 for ( uint16_type i = 0; i < element_type::numVertices; ++i )
3425 for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l, ++lc )
3427 if (__elt.point( i ).id()==pit->id() )
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";
3439 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3446 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3453 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3460 for(
auto pit = M.beginEdgeWithMarker( M.markerName(
"WireBasket") ),
3461 pen = M.endEdgeWithMarker( M.markerName(
"WireBasket") );
3464 auto __elt = M.element( *pit->elements().begin() );
3465 std::cout <<
"Adding wirebasket edge " << pit->id() <<
" using element " << __elt.id() <<
"\n";
3469 for ( uint16_type i = 0; i < element_type::numEdges; ++i )
3471 for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
3473 if (__elt.edge( i ).id()==pit->id() )
3475 size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;
3476 int32_type sign = 1;
3478 if ( __elt.edgePermutation( i ).value() == edge_permutation_type::IDENTITY )
3482 else if ( __elt.edgePermutation( i ).value() == edge_permutation_type::REVERSE_PERMUTATION )
3485 if ( fe_type::is_modal )
3488 sign = ( l%2 )?( -1 ):( 1 );
3493 gDof += fe_type::nDofPerEdge - 1 - l ;
3496 FEELPP_ASSERT( 0 ).error (
"invalid edge permutation" );
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";
3506 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3513 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3520 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3526 std::vector<std::string>
faces = assign::list_of(
"TOP")(
"BOTTOM")(
"NORTH")(
"EAST")(
"WEST")(
"SOUTH");
3527 BOOST_FOREACH(
auto face, faces )
3531 for(
auto pit = faces.template get<1>(), pen = faces.template get<2>(); pit!=pen; ++pit )
3533 auto __elt = M.element( *pit->elements().begin() );
3534 std::cout <<
"Adding face " << pit->id() <<
" with marker " << face <<
" using element " << __elt.id() <<
"\n";
3538 for ( uint16_type i = 0; i < element_type::numFaces; ++i )
3540 face_permutation_type permutation = __elt.facePermutation( i );
3541 FEELPP_ASSERT( permutation != face_permutation_type( 0 ) ).error (
"invalid face permutation" );
3548 int MaxOrder = int( ( 3 + std::sqrt( 1+8*fe_type::nDofPerFace ) )/2 ) - 2;
3550 for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
3552 if (__elt.face( i ).id()==pit->id() )
3559 size_type gDof = __elt.face( i ).id() * fe_type::nDofPerFace;
3560 int32_type sign = 1;
3568 MaxOrder = MaxOrder-1;
3571 if ( !fe_type::is_modal )
3574 if ( permutation == face_permutation_type( 1 ) || fe_type::nDofPerFace == 1 )
3578 gDof += vector_permutation[permutation][l];
3585 if ( permutation == face_permutation_type( 2 ) )
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";
3607 template<
typename MeshType,
typename FEType,
typename PeriodicityType>
3608 std::pair<std::map<size_type,size_type>,std::map<size_type,size_type> >
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() );
3615 if ( it_elt == en_elt )
3616 return std::make_pair(doftopid,pidtodof);
3618 for (
size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
3620 for ( uint16_type i = 0; i < M_mesh->numLocalVertices(); ++i )
3622 int ncdof = is_product?nComponents:1;
3623 for ( uint16_type c1 = 0; c1 < ncdof; ++c1 )
3625 const size_type gDof = ( it_elt->point( i ).id() );
3626 size_type thedof = boost::get<0>( localToGlobal( it_elt->id(), i, c1 ) );
3628 pidtodof[ncdof*gDof+c1] = thedof;
3629 doftopid[thedof] = ncdof*gDof+c1;
3634 if ( !fname.empty() )
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 )
3645 ofs << p.first <<
" " << p.second <<
"\n";
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 )
3653 ofs2 << p.first <<
" " << p.second <<
"\n";
3657 return std::make_pair(doftopid,pidtodof);
3663 #if defined(FEELPP_ENABLE_MPI_MODE)
3667 #endif //_DOFTABLE_HH