Logo Search packages:      
Sourcecode: feel++ version File versions  Download package

doftable.hpp

Go to the documentation of this file.
/*
  This file is part of the Feel library

  Copyright (C) 2010 Université de Grenoble 1

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 3.0 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/
/**
 * \file doftable.hpp
 * \author Christophe Prud'homme
 */
#ifndef _DOFTABLE_HH
#define _DOFTABLE_HH

#include <set>
#include <map>
#include <vector>
#include <algorithm>

#include <boost/foreach.hpp>
#include <boost/multi_array.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_comparison.hpp>
#include <boost/tuple/tuple_io.hpp>
#include <boost/fusion/algorithm/iteration/accumulate.hpp>

#include <Eigen/Core>
#include<Eigen/StdVector>

#include <feel/feelcore/feel.hpp>
#include <feel/feelalg/glas.hpp>
#include <feel/feelpoly/mapped.hpp>
#include <feel/feelalg/datamap.hpp>
#include <feel/feeldiscr/dof.hpp>

namespace Feel
{
// import fusion namespace in Feel
namespace fusion = boost::fusion;

/**
 * \class DofTable
 * \ingroup SpaceTime
 * \brief Local-to-global Degree of Freedom table
 *
 * \author Christophe Prud'homme
 * \author Goncalo Pena
 */
template<typename MeshType,  typename FEType, typename PeriodicityType>
00062 class DofTable : public DataMap
{
    typedef DataMap super;
public:

    /**
     * mesh type
     */
00070     typedef MeshType mesh_type;
    typedef FEType fe_type;
    typedef boost::shared_ptr<FEType> fe_ptrtype;


    typedef typename mesh_type::pid_element_const_iterator pid_element_const_iterator;
    typedef typename mesh_type::element_const_iterator element_const_iterator;
    typedef typename mesh_type::element_type element_type;
    typedef typename mesh_type::face_type face_type;
    typedef typename mesh_type::gm_ptrtype gm_ptrtype;
    typedef typename mesh_type::gm_type gm_type;

    typedef typename fe_type::matrix_type matrix_type;
    typedef typename fe_type::value_type value_type;
    typedef typename fe_type::reference_convex_type reference_convex_type;
    typedef typename fe_type::points_type points_type;

    typedef typename reference_convex_type::super convex_type;

    typedef typename element_type::edge_permutation_type edge_permutation_type;
    typedef typename element_type::face_permutation_type face_permutation_type;

    static const uint16_type nOrder = fe_type::nOrder;
    static const uint16_type nDim = mesh_type::nDim;
    static const uint16_type Shape = mesh_type::Shape;
    static const uint16_type nComponents = fe_type::nComponents;
    static const uint16_type nComponents1 = fe_type::nComponents1;
    static const uint16_type nComponents2 = fe_type::nComponents2;


    static const bool is_continuous = fe_type::isContinuous;
    static const bool is_discontinuous_locally = fe_type::continuity_type::is_discontinuous_locally;
    static const bool is_discontinuous_totally = fe_type::continuity_type::is_discontinuous_totally;

    static const bool is_scalar = FEType::is_scalar;
    static const bool is_vectorial = FEType::is_vectorial;
    static const bool is_tensor2 = FEType::is_tensor2;
    static const bool is_modal = FEType::is_modal;
    static const bool is_product = FEType::is_product;

    static const bool is_p0_continuous = ( (nOrder == 0) && is_continuous );

    static const uint16_type nDofPerElement = mpl::if_<mpl::bool_<is_product>, mpl::int_<FEType::nLocalDof*nComponents1>, mpl::int_<FEType::nLocalDof> >::type::value;

    typedef PeriodicityType periodicity_type;
    static const bool is_periodic = periodicity_type::is_periodic;

    //typedef ContinuityType continuity_type;
    typedef typename fe_type::continuity_type continuity_type;

    typedef DofTable<MeshType, FEType, PeriodicityType> self_type;


    /**
     * A global dof is defined by its index in the global table
     * and by its sign.
     **/

00128     typedef boost::tuple<size_type, int16_type, bool> global_dof_type;


    /**
     * A global dof from face is defined by
     * -its index in the global table
     * -its sign
     * -bool : periodicity
     * -local number in the element
     **/

00139     typedef boost::tuple<size_type, int16_type, bool, int16_type> global_dof_fromface_type;


    /**
     * Type for the localToGlobal table.
     */
00145     typedef boost::multi_array<global_dof_type,2> Container;
    typedef boost::multi_array<global_dof_fromface_type,2> Container_fromface;

    typedef typename Container::array_view<1>::type indices_per_element_type;

    typedef typename node<value_type>::type node_type;

    typedef boost::array<Container::index,2> dims_type;
    typedef typename Container::index_range range_type;
    typedef boost::tuple<node_type, size_type, uint16_type > dof_point_type;
    typedef std::vector<dof_point_type> dof_points_type;
    typedef typename std::vector<dof_point_type>::iterator dof_points_iterator;
    typedef typename std::vector<dof_point_type>::const_iterator dof_points_const_iterator;

    /**
     * Tuple that holds a size_type \p elt 1 uint16_type \p l and 1
     * uint16_type ent
     * \p elt shall be an element index
     * \p l shall be the local index of the dof in the element
     * \p ent shall be the entity the dof belongs to (0: vertex, 1: edge, 2: face, 3: volume)
     */
00166     typedef boost::tuple<size_type, uint16_type, uint16_type, uint16_type> local_dof_type;

    typedef boost::tuple<uint16_type&,size_type&> ref_shift_type;

    /**
     * Type that hold the map between a global dof and the elements
     */
00173     typedef std::map<size_type, std::list<local_dof_type> >  dof_element_type;
    typedef typename dof_element_type::iterator dof_iterator;
    typedef typename dof_element_type::const_iterator dof_const_iterator;

    typedef std::list<local_dof_type>::const_iterator ldof_const_iterator;

    typedef boost::tuple<uint16_type,uint16_type,size_type> dof_type;
    typedef std::map<dof_type, size_type> dof_map_type;
    typedef std::map<dof_type, size_type>::iterator dof_map_iterator;
    typedef std::map<dof_type, size_type, size_type>::const_iterator dof_map_const_iterator;

    typedef std::map<size_type, std::set<size_type> > dof_procset_type;
    /**
     * This type is useful to construct the sign map in the modal case
     **/

00189     typedef ublas::vector<bool> face_sign_info_type;


    typedef Eigen::Matrix<int, nDofPerElement, 1> localglobal_indices_type;

    /**
     * Type for the permutations to be done in the faces
     **/

00198     typedef ublas::vector<uint16_type> permutation_vector_type;

    typedef boost::tuple<element_type const*, face_type const*> element_face_pair_type;
    typedef std::list<element_face_pair_type> periodic_element_list_type;
    typedef typename periodic_element_list_type::iterator periodic_element_list_iterator;
    typedef typename periodic_element_list_type::const_iterator periodic_element_list_const_iterator;
    typedef boost::tuple<size_type /*element id*/, uint16_type /*lid*/, size_type /*gDof*/, uint16_type /*type*/> periodic_dof_type;
    typedef std::multimap<size_type /*gid*/, periodic_dof_type> periodic_dof_map_type;

    /**
     * @brief The minimal constructor
     *
     * @param _fe reference element
     *
     */
    DofTable( fe_ptrtype const& _fe, periodicity_type const& periodicity );

    /**
     * copy constructor
     *
     * @param dof2 a dof object instance
     */
    DofTable( const DofTable & dof2 );

    /**
     * @brief  Constructor accepting a mesh as parameter
     *
     *  @param mesh a RegionMesh3D
     *  @param _fe reference element
     */
    DofTable( mesh_type& mesh, fe_ptrtype const& _fe, periodicity_type const& periodicity );

    /**
     * \return the number of dof for faces on the boundary
     */
00233     uint16_type nDofPerFaceOnBoundary() const
    {
        return _M_n_dof_per_face_on_bdy;
    }

    indices_per_element_type indices( size_type id_el ) const
    {
        return _M_el_l2g[ boost::indices[id_el][range_type()] ];
    }

    std::vector<size_type> getIndices( size_type id_el ) const
    {
        size_type nldof =
            fe_type::nDofPerVolume * element_type::numVolumes +
            fe_type::nDofPerFace * element_type::numGeometricFaces +
            fe_type::nDofPerEdge * element_type::numEdges +
            fe_type::nDofPerVertex * element_type::numVertices;
        const size_type ntdof = is_product?nComponents*nldof:nldof;
        std::vector<size_type> ind( ntdof );
        for( size_type i = 0; i < ntdof; ++i )
            ind[i] = boost::get<0>( _M_el_l2g[ id_el][ i ] );
        return ind;
    }

    /**
     * \return the dof view
     */
00260     dof_container_type const& dofView() const { return M_dof_view; }


    /**
     * \return the number of processors \c dof belongs to
     */
00266     size_type dofNProc( size_type dof ) const { return _M_dof_procset.find( dof )->second.size(); }

    /**
     * \return the set of processors \c dof belongs to
     */
00271     std::set<size_type> dofProcSet( size_type dof ) const { return _M_dof_procset.find( dof )->second; }

    /**
     * @return the coordinates of the nodal dofs associated with the
     * element \p el
     */
00277     const dof_point_type& dofPoint( size_type i ) const { return M_dof_points[i]; }

    /**
     * @return an iterator at the beginning of dof points
     */
00282     dof_points_const_iterator dofPointBegin() const { return M_dof_points.begin(); }

    /**
     * @return an iterator at the beginning of dof points
     */
00287     dof_points_iterator dofPointBegin() { return M_dof_points.begin(); }

    /**
     * @return an iterator at the end of dof points
     */
00292     dof_points_const_iterator dofPointEnd() const { return M_dof_points.end(); }

    /**
     * @return an iterator at the end of dof points
     */
00297     dof_points_iterator dofPointEnd() { return M_dof_points.end(); }


    /**
     * insted of creating the dof indices on the fly, get them from a
     * vector. The situation typically arises when we want to have dof
     * correspondance between two spaces
     *
     * \see OperatorLagrangeP1
     */
00307     void setDofIndices( std::vector<boost::tuple<size_type, uint16_type, size_type> > const& dof )
    {
        M_dof_indices.resize( dof.size() );
        std::copy( dof.begin(), dof.end(), M_dof_indices.begin() );

        if ( dof.empty() )
            return ;
        size_type nldof =
            fe_type::nDofPerVolume * element_type::numVolumes +
            fe_type::nDofPerFace * element_type::numGeometricFaces +
            fe_type::nDofPerEdge * element_type::numEdges +
            fe_type::nDofPerVertex * element_type::numVertices;
        typedef boost::tuple<size_type, uint16_type,size_type> thedof_type;
#if 1
        std::set<size_type> eltid;
        std::set<size_type> dofs;

        BOOST_FOREACH( thedof_type thedof,  dof )
            {
                eltid.insert( boost::get<0>( thedof ) );
                dofs.insert( boost::get<2>( thedof ) );
            }
#endif
        const size_type ntdof = is_product?nComponents*nldof:nldof;
        _M_el_l2g.resize( boost::extents[eltid.size()][ntdof] );


        BOOST_FOREACH( thedof_type thedof,  dof )
            {
                _M_el_l2g[ boost::get<0>( thedof ) ][ boost::get<1>( thedof ) ] = boost::get<2>( thedof );

            }
        this->_M_first_df[0] = 0;
        this->_M_last_df[0] = dofs.size()-1;
        this->_M_n_dofs = dofs.size();
    }

    /**
     * \return the dof index
     */
00347     size_type dofIndex( size_type dof ) const
    {
        return dof;
#if 0
        if ( M_dof_indices.empty() )
            return dof;
        FEEL_ASSERT( dof < M_dof_indices.size() )( dof )( M_dof_indices.size() ).warn( "invalid dof index" );
        return M_dof_indices[dof];
#endif
    }

    /**
     * \return the local to global indices
     */
00361     localglobal_indices_type const& localToGlobalIndices( size_type ElId ) { return M_locglob_indices[ElId]; }

    /**
     * \return the signs of the global dof (=1 in nodal case, +-1 in modal case)
     */
00366     localglobal_indices_type const& localToGlobalSigns( size_type ElId ) { return M_locglob_signs[ElId]; }

    /**
     * \return the specified entries of the localToGlobal table
     *
     * \param ElId the element ID
     * \param localNode the local DOF numbering (starting from 1)
     * \param c the component index, default is 0-th component
     *
     * \return the global numbering of a DOF, given an element and the local numbering
     */

00378     global_dof_type localToGlobal( const size_type ElId,
                                   const uint16_type localNode,
                                   const uint16_type c = 0 ) const
    {
        return _M_el_l2g[ ElId][ fe_type::nLocalDof * c  + localNode ];
    }

    global_dof_fromface_type faceLocalToGlobal( const size_type ElId,
                                                const uint16_type localNode,
                                                const uint16_type c = 0 ) const
    {
        const size_type nDofF = ( face_type::numVertices * fe_type::nDofPerVertex +
                                  face_type::numEdges * fe_type::nDofPerEdge +
                                  face_type::numFaces * fe_type::nDofPerFace );
        return _M_face_l2g[ ElId][ nDofF*c+localNode ];
    }

00395     struct element_access
    {
        element_access( DofTable const& __d )
            :
            _M_d( __d )
        {}
        global_dof_type operator()( size_type __id, uint16_type __loc, uint16_type c = 0 ) const
        {
            return _M_d._M_el_l2g[ __id][ fe_type::nLocalDof * c+ __loc ];
        }
        uint16_type localDofInElement(size_type __id, uint16_type __loc, uint16_type c = 0 ) const
        {
            return __loc;
        }
        DofTable const& _M_d;
    };
    friend class element_access;

00413     struct face_access
    {
        static const uint16_type nDofF = ( face_type::numVertices * fe_type::nDofPerVertex +
                                           face_type::numEdges * fe_type::nDofPerEdge +
                                           face_type::numFaces * fe_type::nDofPerFace );

        face_access( DofTable const& __d )
            :
            _M_d( __d )
        {}
        global_dof_fromface_type operator()( size_type __id, uint16_type __loc, uint16_type c = 0 ) const
        {
            return _M_d._M_face_l2g[ __id][nDofF*c+__loc];
        }

        uint16_type localDofInElement(size_type __id, uint16_type __loc, uint16_type c = 0 ) const
        {
            return boost::get<3>(_M_d._M_face_l2g[ __id][nDofF*c+__loc]);
        }

        DofTable const& _M_d;
    };
    friend class face_access;

    /**
     * @brief local to global mapping
     */
    template<typename Elem>
    typename mpl::if_<mpl::equal_to<mpl::int_<Elem::nDim>,mpl::int_<nDim> >,
                      global_dof_type,
                      global_dof_fromface_type >::type
00444     localToGlobal( Elem const& El, const uint16_type localNode, const uint16_type c = 0 ) const
    {
        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;
        //Debug( 5005 ) << "dof:(" << El.id() << ", " << localNode << ")= "
        //<< access_type(*this)( El.id(), localNode, c ) << "\n";
        return access_type(*this)( El.id(), localNode, c );
    }

    template<typename Elem>
    uint16_type
    localDofInElement(Elem const& El, const uint16_type localNode, const uint16_type c = 0 ) const
    {
        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;

        return (access_type(*this)).localDofInElement(El.id(), localNode, c );
    }

    /**
     * Number of elements in mesh
     */
00464     size_type numElements() const
    {
        return _M_n_el;
    }

    /**
     * Number of local vertices (in an element)
     */
00472     uint16_type numLocalVertices() const
    {
        return fe_type::numVertices;
    }

    /**
     * Number of local edges (in an element)
     */
00480     uint16_type numLocalEdges() const
    {
        return fe_type::numEdges;
    }

    /**
     * Number of local faces (in an element)
     */
00488     uint16_type numLocalFaces() const
    {
        return fe_type::numFaces;
    }

    /**
     * show some information about the dof instance
     */
    void showMe() const;

    void dump() const
    {
#if 0
        for ( size_type __i = 0;__i < _M_face_l2g.nrows();++__i )
            {
                for ( size_type __l = 0;__l < _M_face_l2g.ncols();++__l )
                    {
                        std::cout << "face " << __i << " local " << __l
                                  << " to global " << _M_face_l2g[ __i][ __l ] << "\n";
                    }
            }
#endif // 0
    }

    /**
     * set element to done
     */
00515     void setElementDone( size_type elt )
    {
        _M_elt_done[elt] = true;
    }
    bool isElementDone( size_type elt, int c = 0 )
    {
        const size_type nldof =
            fe_type::nDofPerVolume * element_type::numVolumes +
            fe_type::nDofPerFace * element_type::numGeometricFaces +
            fe_type::nDofPerEdge * element_type::numEdges +
            fe_type::nDofPerVertex * element_type::numVertices;
        const int ncdof = is_product?nComponents:1;
        for( int c = 0; c < ncdof; ++c )
        {
            for( uint16_type l = 0; l < nldof; ++l )
                if (boost::get<0>(_M_el_l2g[elt][c*nldof+l])==invalid_size_type_value)
                    return false;
        }
        return true;
    }
    void addDofFromElement( element_type const& __elt,
                            size_type& next_free_dof,
                            size_type processor = 0,
                            size_type shift = 0 )
    {

        size_type nldof =
            fe_type::nDofPerVolume * element_type::numVolumes +
            fe_type::nDofPerFace * element_type::numGeometricFaces +
            fe_type::nDofPerEdge * element_type::numEdges +
            fe_type::nDofPerVertex * element_type::numVertices;

        if ( !this->isElementDone( __elt.id() ) )
        {
            Debug(11111) << "adding dof from element " << __elt.id() << "\n";
            size_type gdofcount = shift;
            Debug( 5005 ) << "next_free_dof " << next_free_dof  << "\n";
            Debug( 5005 ) << "current dof " << dofIndex( next_free_dof ) << "\n";


            /*
             * Only in the continuous , we need to have the ordering [vertex,edge,face,volume]
             */
            if ( is_continuous || is_discontinuous_locally )
            {

                /* idem as above but for local element
                   numbering except that it is
                   reset to 0 after each element */
                uint16_type ldofcount = 0;

                /* pack the shifts into a tuple */
                boost::tuple<uint16_type&,size_type&> shifts = boost::make_tuple( boost::ref(ldofcount),
                                                                                  boost::ref(gdofcount) );

                /* \warning: the order of function calls is
                   crucial here we order the degrees of freedom
                   wrt the topological entities of the mesh
                   elements from lowest dimension (vertex) to
                   highest dimension (element)
                */
                addVertexDof( __elt, processor, next_free_dof, shifts  );
                addEdgeDof( __elt, processor, next_free_dof, shifts );
                addFaceDof( __elt, processor, next_free_dof, shifts );
                addVolumeDof( __elt, processor, next_free_dof, shifts );
            }
            else
            {

                size_type ie = __elt.id();

                const int ncdof = is_product?nComponents:1;
                for ( uint16_type l = 0; l < nldof; ++l )
                {
                    for( int c = 0; c < ncdof; ++c, ++next_free_dof )
                    {
                        _M_el_l2g[ ie][ fe_type::nLocalDof*c + l ] = boost::make_tuple(( dofIndex(next_free_dof)) , 1, false );
                    }
                }
            }
        }
        else
        {
            Debug(11111) << "element " << __elt.id() << "has already been taken care of\n";
        }
    }
    /**
     * \return the local to global map
     */
    std::map<size_type, std::pair<size_type, size_type> >
00605     localToGlobalMap() const
    { return _M_local2global;}

    /**
     * local(on processor) dof to global dof
     */
00611     size_type localToGlobal( size_type l ) const
    {
        return _M_local2global.find( l )->second.first;
    }

    /**
     * Add a new periodic dof to the dof map and to the list of periodic dof \p
     * periodic_dof for a given tag \p tag
     *
     * \arg __elt geometric element
     * \arg __face  face of the element that holds the periodic dof
     * \arg next_free_dof integer that will get possibilty incremented and holds the current number of dof
     * \arg periodic_dof table of periodic dof
     * \arg tag the boundary tag associated with the periodic dof
     */
00626     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 )
    {
        addVertexPeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag, mpl::bool_<fe_type::nDofPerVertex>() );
    }
    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> ) {}
    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> );

    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 )
    {
        addEdgePeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag , mpl::bool_<fe_type::nDofPerEdge>(), mpl::int_<nDim>() );
    }
    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> ) {}
    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> ) {}
    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> ) {}
    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> ) {}
    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> );
    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> );


    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 )
    {
        addFacePeriodicDof( __elt, __face, next_free_dof, periodic_dof, tag, mpl::bool_<fe_type::nDofPerFace>() );
    }
    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> ) {}
    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> );

    /**
     * Initialize the dof map table
     */
    void initDofMap( mesh_type& M );

    /**
     * build the dof map
     */
00668     void build( mesh_type* M ) { this->build( *M ); }

    /**
     * build the dof map
     */
00673     void build( boost::shared_ptr<mesh_type>& M ) { this->build( *M ); }

    /**
     * build the dof map
     */
00678     void build( mesh_type& M )
    {
        _M_elt_done.resize( M.numElements() );
        std::fill( _M_elt_done.begin(), _M_elt_done.end(), false );

        Debug( 5015 ) << "[Dof::build] initDofMap\n";
        this->initDofMap( M );

        Debug( 5015 ) << "[Dof::build] start building dof map\n";
        size_type start_next_free_dof = 0;
        Debug( 5015 ) << "[Dof::build] start_next_free_dof = " << start_next_free_dof << "\n";
        if ( is_periodic )
            {
                Debug(5015) << "[build] call buildPeriodicDofMap()\n";
                start_next_free_dof = this->buildPeriodicDofMap( M );
                Debug( 5015 ) << "[Dof::build] start_next_free_dof(after periodic) = " << start_next_free_dof << "\n";
            }

        if ( is_discontinuous_locally )
            {
                Debug(5015) << "[build] call buildLocallyDiscontinuousDofMap()\n";
                start_next_free_dof = this->buildLocallyDiscontinuousDofMap( M, start_next_free_dof );
                Debug( 5015 ) << "[Dof::build] start_next_free_dof(after local discontinuities) = " << start_next_free_dof << "\n";
            }
        Debug(5015) << "[build] call buildDofMap()\n";
        this->buildDofMap( M, start_next_free_dof );

#if !defined(NDEBUG)
        Debug(5015) << "[build] check that all elements dof were assigned()\n";
        element_const_iterator fit, fen;
        boost::tie( fit, fen ) = M.elementsRange();
        std::vector<boost::tuple<size_type,uint16_type,size_type> > em;
        for( ; fit != fen; ++fit )
            {
                const int ncdof = is_product?nComponents:1;
                for ( uint16_type c = 0; c < ncdof; ++c )
                    if (!this->isElementDone( fit->id(), c ) )
                    {
                        em.push_back( boost::make_tuple( fit->id(), c, fit->marker().value() ) );
                    }
            }
        if ( !em.empty() )
            {
                Debug(5015) << "[build] some element dof were not assigned\n";
                for( size_type i = 0; i < em.size(); ++i )
                    {
                        Debug(5015) << " - element " << boost::get<0>( em[i] ) << " c=" << boost::get<1>( em[i] )
                                    << " m=" << boost::get<2>( em[i] ) << "\n";
                    }
            }
        else
            {
                Debug(5015) << "[build] check that all elements dof were assigned: OK\n";
            }
#endif // NDEBUG
        Debug( 5015 ) << "[Dof::build] n_dof = " << this->nLocalDof() << "\n";

        Debug(5015) << "[build] call buildBoundaryDofMap()\n";
        this->buildBoundaryDofMap( M );

        Debug(5015) << "[Dof::build] done building the map\n";
    }

    /**
     * build dof map associated to the periodic dof, must be called
     * before buildDofMap
     */
    size_type buildPeriodicDofMap( mesh_type& M );

    /**
     * build dof associated to local discontinuities
     */
    size_type buildLocallyDiscontinuousDofMap( mesh_type& M, size_type start_next_free_dof );

    /**
     * @brief Build the localToGlobal table
     *
     *  \param mesh mesh
     */
    void buildDofMap( mesh_type& mesh, size_type start_next_free_dof = 0 );

    /**
     * @brief Build the localToGlobal table for the boundary
     *
     * \param mesh A mesh
     */
    void buildBoundaryDofMap( mesh_type& mesh );


    /**
     * \return the dictionnary for the global dof
     */
00770     dof_map_type const& mapGDof() const { return map_gdof; }

    /**
     * \return the dictionnary for the global dof
     */
00775     dof_map_type& mapGDof() { return map_gdof; }

    /**
     * clear the dictionnary
     */
00780     void clearMapGDof() { map_gdof.clear(); }

    /**
     * set the dictionnary for the dictionary of the global dof
     */
00785     void setMapGDof( dof_map_type const& mapdof ) { map_gdof = mapdof; }

    /**
     * The dof are ordered such that they are contiguous per element
     * and components. This way an extraction of the dof indices in
     * one element allows to extract a view of the corresponding
     * coefficient in a given basis which is then very helpful for
     * interpolation for example.
     *
     * \param ie index of the element
     * \param lc_dof index of the dof in the element
     * \param lc local index of the entity associated with the dof in the element
     * \param gDof global dof index
     * \param pDof dof index in the processor
     *
     * \return the index of the next free dof in the processor
     */
00802     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 )
    {
        bool res = true;
        const int ncdof = is_product?nComponents:1;
        for( int c = 0; c < ncdof; ++c )
        {
            gDof.template get<1>() = c;
            uint16_type lc_dof = fe_type::nLocalDof*c+l_dof;
            Feel::detail::ignore_unused_variable_warning(lc);
            dof_map_iterator itdof = map_gdof.find( gDof );
            dof_map_iterator endof = map_gdof.end();
            bool __inserted = false;
            if ( itdof == endof )
            {
                Debug( 5005 ) << "[dof] dof (" << gDof.get<0>() << "," << gDof.get<1>() << "," << gDof.get<2>() << ") not yet inserted in map\n";
                boost::tie( itdof, __inserted ) = map_gdof.insert( std::make_pair( gDof, dofIndex( pDof ) ) );

                pDof += 1;

                FEEL_ASSERT( __inserted == true )( ie )( lc_dof )
                    (gDof.get<0>())(gDof.get<1>())( gDof.get<2>() )
                    ( processor )( itdof->second ).error( "dof should have been inserted");
            }
            else
            {
                Debug( 5005 ) << "[dof] dof (" << gDof.get<0>() << ","
                              << gDof.get<1>()
                              << "," << gDof.get<2>()
                              << ") already inserted in map with dof_id = " << itdof->second << "\n";
            }

#if !defined( NDEBUG )
            Debug( 5005 ) << "global dof = " << itdof->second
                          << " local dof = " << fe_type::nLocalDof*itdof->first.get<1>() + lc_dof
                          << " element = " << ie
                          << " entity = " << itdof->first.get<0>()
                          << " component = " << itdof->first.get<1>()
                          << " index = " << itdof->first.get<2>() << "\n";
#endif
            // make sure that no already created dof is overwritten here (may be done alsewhere)
            if ( boost::get<0>( _M_el_l2g[ ie][ lc_dof ] ) == invalid_size_type_value )
            {

                FEEL_ASSERT( itdof->first == gDof ).error( "very bad logical error in insertDof" );

                FEEL_ASSERT( lc_dof >= fe_type::nLocalDof*itdof->first.get<1>() &&
                             lc_dof < fe_type::nLocalDof*(itdof->first.get<1>()+1) )
                    ( lc_dof )
                    ( fe_type::nLocalDof*itdof->first.get<1>() ).error( "invalid local dof index" );


                // add processor to the set of processors this dof belongs to
                _M_dof_procset[ itdof->second+shift ].insert( processor );

                _M_el_l2g[ ie][ lc_dof ] = boost::make_tuple(itdof->second+shift, sign, is_dof_periodic );

#if !defined(NDEBUG)
                _M_dof2elt[itdof->second+shift].push_back( boost::make_tuple( ie, lc_dof, lc, itdof->first.get<0>() ) );
#endif
#if 0
                M_dof_view.insert( Dof( itdof->second+shift,      // global index
                                        sign,                     // sign
                                        itdof->first.get<0>(),    // entity type
                                        false,                    // is on boundary ?
                                        0                         // marker
                                       ) );
#endif
#if 0
                typedef Container::index index;

                const int ntdof = is_product?nComponents:1;
                for ( index i2 = 0; i2 < ntdof*fe_type::nLocalDof; ++i2 )
                    Debug() << "dof table( " << ie << ", " << lc  << ")=" << boost::get<0>(_M_el_l2g[ ie][ i2 ]) << "\n";
#endif
            }
            res = res && ( __inserted || ( ( boost::get<0>( _M_el_l2g[ ie][ lc_dof ] ) == invalid_size_type_value ) && shift ));
        }
        return res;
    }

private:

    void addVertexDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                       ref_shift_type& shifts  )
    {
        addVertexDof( __elt, processor, next_free_dof, shifts, mpl::bool_<fe_type::nDofPerVertex>() );
    }
    void addVertexDof( element_type const& /*M*/, uint16_type /*processor*/,  size_type& /*next_free_dof*/,
                       ref_shift_type& /*shifts*/, mpl::bool_<false> )
    {}
    void addVertexDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                       ref_shift_type& shifts, mpl::bool_<true> )
    {
        uint16_type local_shift;
        size_type global_shift;
        boost::tie( local_shift, global_shift ) = shifts;


        size_type ie = __elt.id();

        uint16_type lc = local_shift;
        for ( uint16_type i = 0; i < element_type::numVertices; ++i )
            {
                for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l, ++lc )
                    {
                        //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
                        const size_type gDof = ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
                        this->insertDof( ie, lc, i, boost::make_tuple(0, 0, gDof),
                                         processor, next_free_dof, 1, false, global_shift );
                    }
            }
        // update shifts
        shifts.get<0>() = lc;

#if !defined(NDEBUG)
        Debug( 5005 ) << "[Dof::updateVolumeDof(addVertexDof] vertex proc" << processor << " next_free_dof = " << next_free_dof << "\n";
#endif
    }
    void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                     ref_shift_type& shifts )
    {
        return addEdgeDof( __elt,
                           processor,
                           next_free_dof,
                           shifts,
                           mpl::int_<fe_type::nDim>(),
                           mpl::bool_<fe_type::nDofPerEdge>() );
    }
    void addEdgeDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
                     ref_shift_type& /*shifts*/, mpl::int_<1>, mpl::bool_<false> )
    {}

    void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                     ref_shift_type& shifts, mpl::int_<1>, mpl::bool_<true> )
    {
        uint16_type local_shift;
        size_type global_shift;
        boost::tie( local_shift, global_shift ) = shifts;

        size_type ie = __elt.id();
        uint16_type lc = local_shift;
        for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
            {
                const size_type gDof = is_p0_continuous? l:ie * fe_type::nDofPerEdge + l;
                this->insertDof( ie, lc, l, boost::make_tuple(1, 0, gDof), processor, next_free_dof, 1, false, global_shift );
            }
        // update shifts
        shifts.get<0>() = lc;
#if !defined(NDEBUG)
        Debug( 5005 ) << "[Dof::addEdgeDof(1)] element proc" << processor << " next_free_dof = " << next_free_dof << "\n";
#endif
    }
    void addEdgeDof( element_type const& /*__elt*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
                     ref_shift_type& /*shifts*/, mpl::int_<2>, mpl::bool_<false> )
    {}
00965     void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                     ref_shift_type& shifts, mpl::int_<2>, mpl::bool_<true> )
    {
        uint16_type local_shift;
        size_type global_shift;
        boost::tie( local_shift, global_shift ) = shifts;

        size_type ie = __elt.id();
        uint16_type lc = local_shift;

        /** The boundary dofs are constructed in the same way if the basis is modal **/

        for ( uint16_type i = 0; i < element_type::numEdges; ++i )
            {
                for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
                    {
                        size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;
                        int32_type sign = 1;


                        if ( __elt.edgePermutation(i).value()  == edge_permutation_type::IDENTITY )
                            {
                                gDof += l ; // both nodal and modal case
                            }
                        else if ( __elt.edgePermutation(i).value()  == edge_permutation_type::REVERSE_PERMUTATION)
                            {

                                if(fe_type::is_modal)
                                    {
                                        //only half of the modes (odd polynomial order) are negative.
                                        sign = (l%2)?(-1):(1);
                                        gDof += l;
                                    }
                                else
                                    gDof += fe_type::nDofPerEdge - 1 - l ;
                            }
                        else
                            FEEL_ASSERT( 0 ).error ( "invalid edge permutation" );
                        this->insertDof( ie, lc, i, boost::make_tuple(1, 0, gDof), processor, next_free_dof, sign, false, global_shift );
                    }
            }

        // update shifts
        shifts.get<0>() = lc;
#if !defined(NDEBUG)
        Debug( 5005 ) << "[Dof::addEdgeDof] edge proc" << processor << " next_free_dof = " << next_free_dof << "\n";
#endif
    }

    void addEdgeDof( element_type const& /*__elt*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
                     ref_shift_type& /*shifts*/, mpl::int_<3>, mpl::bool_<false> )
    {}

    void addEdgeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                     ref_shift_type& shifts, mpl::int_<3>, mpl::bool_<true> )
    {
        uint16_type local_shift;
        size_type global_shift;
        boost::tie( local_shift, global_shift ) = shifts;

        size_type ie = __elt.id();
        uint16_type lc = local_shift;
        for ( uint16_type i = 0; i < element_type::numEdges; ++i )
            {
                for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l, ++lc )
                    {
                        size_type gDof = __elt.edge( i ).id() * fe_type::nDofPerEdge;

                        int32_type sign = 1;

                        if ( __elt.edgePermutation(i).value()  == edge_permutation_type::IDENTITY )
                            {
                                gDof += l ; // both nodal and modal case
                            }
                        else if ( __elt.edgePermutation(i).value()  == edge_permutation_type::REVERSE_PERMUTATION)
                            {

                                if(fe_type::is_modal)
                                    {
                                        //only half of the modes (odd polynomial order) are negative.
                                        sign = (l%2)?(-1):(1);
                                        gDof += l;
                                    }
                                else
                                    gDof += fe_type::nDofPerEdge - 1 - l ;
                            }
                        else
                            FEEL_ASSERT( 0 ).error ( "invalid edge permutation" );

                        this->insertDof( ie, lc, i, boost::make_tuple(1, 0, gDof), processor, next_free_dof, sign, false, global_shift );
                    }
            }
        // update shifts
        shifts.get<0>() = lc;
#if !defined(NDEBUG)
        Debug( 5005 ) << "[Dof::addEdgeDof] edge proc" << processor << " next_free_dof = " << next_free_dof << "\n";
#endif
    }


    void addFaceDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                     ref_shift_type& shifts )
    {
        return addFaceDof( __elt, processor, next_free_dof, shifts, mpl::int_<fe_type::nDim>(), mpl::bool_<fe_type::nDofPerFace>() );
    }
    void addFaceDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
                     ref_shift_type& /*shifts*/, mpl::int_<1>, mpl::bool_<false> )
    {}
    void addFaceDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
                     ref_shift_type& /*shifts*/, mpl::int_<2>, mpl::bool_<false> )
    {}
    void addFaceDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                     ref_shift_type& shifts, mpl::int_<2>, mpl::bool_<true> )
    {
        uint16_type local_shift;
        size_type global_shift;
        boost::tie( local_shift, global_shift ) = shifts;

        size_type ie = __elt.id();
        uint16_type lc = local_shift;
        for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
            {
                const size_type gDof = is_p0_continuous? l:ie * fe_type::nDofPerFace + l;
                this->insertDof( ie, lc, l, boost::make_tuple(2, 0, gDof), processor, next_free_dof, 1, false, global_shift );
            }
        // update shifts
        shifts.get<0>() = lc;
#if !defined(NDEBUG)
        Debug( 5005 ) << "[Dof::addFaceDof(2,true)] face proc" << processor << " next_free_dof = " << next_free_dof << "\n";
#endif
    }
    void addFaceDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
                     ref_shift_type& /*shifts*/, mpl::int_<3>, mpl::bool_<false> )
    {}
    void addFaceDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                     ref_shift_type& shifts, mpl::int_<3>, mpl::bool_<true> )
    {
        uint16_type local_shift;
        size_type global_shift;
        boost::tie( local_shift, global_shift ) = shifts;

        size_type ie = __elt.id();

        uint16_type lc = local_shift;
        for ( uint16_type i = 0; i < element_type::numFaces; ++i )
            {
                face_permutation_type permutation = __elt.facePermutation(i);
                FEEL_ASSERT( permutation != face_permutation_type(0) ).error ( "invalid face permutation" );

                // Polynomial order in each direction
                uint16_type p=1;
                uint16_type q=0;

                // MaxOrder = Order - 2
                int MaxOrder = int((3 + std::sqrt(1+8*fe_type::nDofPerFace))/2) - 2;

                for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l, ++lc )
                    {

                        // TODO: orient the dof indices such
                        // that they match properly the faces
                        // dof of the connected faces. There
                        // are a priori many permutations of
                        // the dof face indices
                        size_type gDof = __elt.face( i ).id() * fe_type::nDofPerFace;
                        int32_type sign = 1;

                        q=q+1;
                        if(q > MaxOrder)
                            {
                                q = 1;
                                p = p+1;
                                MaxOrder = MaxOrder-1;
                            }

                        if ( !fe_type::is_modal )
                            {
                                // no need of permutation is identity or only one dof on face
                                if ( permutation  == face_permutation_type(1) || fe_type::nDofPerFace == 1 )
                                    gDof += l;
                                else
                                    gDof += vector_permutation[permutation][l];
                            }
                        else
                            {
                                gDof += l;
                                if ( permutation == face_permutation_type(2) )
                                    {
                                        // Reverse sign if polynomial order in
                                        // eta_1 direction is odd

                                        if(p%2 == 0)
                                            sign = -1;

                                    }
                            }
                        this->insertDof( ie, lc, i, boost::make_tuple(2, 0, gDof), processor, next_free_dof, sign, false, global_shift );

                    }
            }
        // update shifts
        shifts.get<0>() = lc;
#if !defined(NDEBUG)
        Debug( 5005 ) << "[Dof::addFaceDof<3>] face proc" << processor << " next_free_dof = " << next_free_dof << "\n";
#endif
    }
    void addVolumeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                       ref_shift_type& shifts )
    {
        return addVolumeDof( __elt, processor, next_free_dof, shifts, mpl::bool_<fe_type::nDofPerVolume>() );
    }
    void addVolumeDof( element_type const& /*M*/, uint16_type /*processor*/, size_type& /*next_free_dof*/,
                       ref_shift_type& /*shifts*/, mpl::bool_<false> )
    {}
    void addVolumeDof( element_type const& __elt, uint16_type processor, size_type& next_free_dof,
                       ref_shift_type& shifts, mpl::bool_<true> )
    {
        BOOST_STATIC_ASSERT( element_type::numVolumes );
        uint16_type local_shift;
        size_type global_shift;
        boost::tie( local_shift, global_shift ) = shifts;

        size_type ie = __elt.id();
        uint16_type lc = local_shift;
        for ( uint16_type l = 0; l < fe_type::nDofPerVolume; ++l, ++lc )
            {
                const size_type gDof = is_p0_continuous? l:ie * fe_type::nDofPerVolume + l;
                this->insertDof( ie, lc, l, boost::make_tuple(3, 0, gDof), processor, next_free_dof, 1, false, global_shift );
            }
        // update shifts
        shifts.get<0>() = lc;
#if !defined(NDEBUG)
        Debug( 5005 ) << "[Dof::updateVolumeDof(<2>)] element proc" << processor << " next_free_dof = " << next_free_dof << "\n";
#endif
    }

    template<typename FaceIterator>
    void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc )
    {
        addVertexBoundaryDof( __face_it, lc, mpl::bool_<fe_type::nDofPerVertex>(), mpl::int_<nDim>() );
    }
    template<typename FaceIterator> void addVertexBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<1> ) {}
    template<typename FaceIterator> void addVertexBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<2> ) {}
    template<typename FaceIterator> void addVertexBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<3> ) {}
    template<typename FaceIterator>
    void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<1>  )
    {
        BOOST_STATIC_ASSERT( face_type::numVertices );

        // id of the element adjacent to the face
        // \warning NEED TO INVESTIGATE THIS
        size_type iElAd = __face_it->ad_first();
        FEEL_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );

        // local id of the face in its adjacent element
        uint16_type iFaEl = __face_it->pos_first();
        FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");

        // Loop number of Dof per vertex
        const int ncdof = is_product?nComponents:1;
        for( int c = 0; c < ncdof; ++c )
        {
            for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
            {
                auto temp= this->localToGlobal( iElAd,
                                                iFaEl * fe_type::nDofPerVertex + l,
                                                c );
                _M_face_l2g[ __face_it->id()][ lc++ ] = boost::make_tuple( boost::get<0>(temp),boost::get<1>(temp),boost::get<2>(temp),
                                                                           iFaEl * fe_type::nDofPerVertex + l );
            }
        }
    }
    template<typename FaceIterator>
    void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<2>  )
    {
        BOOST_STATIC_ASSERT( face_type::numVertices );

        // id of the element adjacent to the face
        // \warning NEED TO INVESTIGATE THIS
        size_type iElAd = __face_it->ad_first();
        FEEL_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );

        // local id of the face in its adjacent element
        uint16_type iFaEl = __face_it->pos_first();
        FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");

        size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
                            face_type::numEdges * fe_type::nDofPerEdge +
                            face_type::numFaces * fe_type::nDofPerFace );


        //_M_dof2elt[gDof].push_back( boost::make_tuple( iElAd, lc-1, 48, 0 ) );
        // loop on face vertices
        const int ncdof = is_product?nComponents:1;
        for( int c = 0; c < ncdof; ++c )
        {
            uint16_type lcc=c*ndofF;

            for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
            {
                // local vertex number (in element)
                uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );

                FEEL_ASSERT( iVeEl != invalid_uint16_type_value ).error( "invalid local dof" );

                // Loop number of Dof per vertex
                for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
                {
                    auto temp = this->localToGlobal( iElAd,
                                                     iVeEl * fe_type::nDofPerVertex + l,
                                                     c );
                    _M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>(temp),boost::get<1>(temp),boost::get<2>(temp),
                                                                                iVeEl * fe_type::nDofPerVertex + l );
                }
            }
        }
    }
    template<typename FaceIterator>
    void addVertexBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<3>  )
    {
        addVertexBoundaryDof( __face_it, lc, mpl::bool_<true>(), mpl::int_<2>() );
    }
    template<typename FaceIterator>
    void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc )
    {
        addEdgeBoundaryDof( __face_it, lc, mpl::bool_<fe_type::nDofPerEdge*face_type::numEdges>(), mpl::int_<nDim>() );
    }
    template<typename FaceIterator>
    void addEdgeBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<1> ){}
    template<typename FaceIterator>
    void addEdgeBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<2> ){}
    template<typename FaceIterator>
    void addEdgeBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false>, mpl::int_<3> ){}
    template<typename FaceIterator>
    void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<2> )
    {
        // id of the element adjacent to the face
        // \warning NEED TO INVESTIGATE THIS
        size_type iElAd = __face_it->ad_first();
        FEEL_ASSERT( iElAd != invalid_size_type_value )
            ( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );

        // local id of the face in its adjacent element
        uint16_type iFaEl = __face_it->pos_first();
        FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");
#if !defined(NDEBUG)
        Debug( 5005 ) << " local face id : " << iFaEl << "\n";
#endif
        size_type nVerticesF = face_type::numVertices * fe_type::nDofPerVertex;
        size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
                            face_type::numEdges * fe_type::nDofPerEdge +
                            face_type::numFaces * fe_type::nDofPerFace );


        const int ncdof = is_product?nComponents:1;
        for( int c = 0; c < ncdof; ++c )
        {
            uint16_type lcc=nVerticesF+c*ndofF;

            // Loop number of Dof per edge
            for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
            {
                auto temp = this->localToGlobal( iElAd,
                                                 element_type::numVertices*fe_type::nDofPerVertex +
                                                 iFaEl * fe_type::nDofPerEdge + l,
                                                 c );
                _M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>(temp),boost::get<1>(temp),boost::get<2>(temp),
                                                                            element_type::numVertices*fe_type::nDofPerVertex +
                                                                            iFaEl * fe_type::nDofPerEdge + l );
            }
        }
    }
    template<typename FaceIterator>
    void addEdgeBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true>, mpl::int_<3> )
    {
        //BOOST_STATIC_ASSERT( face_type::numEdges );

        // id of the element adjacent to the face
        // \warning NEED TO INVESTIGATE THIS
        size_type iElAd = __face_it->ad_first();
        FEEL_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );

        // local id of the face in its adjacent element
        uint16_type iFaEl = __face_it->pos_first();
        FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");
#if !defined(NDEBUG)
        Debug( 5005 ) << " local face id : " << iFaEl << "\n";
#endif
        size_type nVerticesF = face_type::numVertices * fe_type::nDofPerVertex;
        size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
                            face_type::numEdges * fe_type::nDofPerEdge +
                            face_type::numFaces * fe_type::nDofPerFace );

        const int ncdof = is_product?nComponents:1;
        for( int c = 0; c < ncdof; ++c )
        {
            uint16_type lcc=nVerticesF+c*ndofF;
            // loop on face vertices
            for ( uint16_type iEdFa = 0; iEdFa < face_type::numEdges; ++iEdFa )
            {
                // local edge number (in element)
                uint16_type iEdEl = element_type::fToE( iFaEl, iEdFa );

                FEEL_ASSERT( iEdEl != invalid_uint16_type_value ).error( "invalid local dof" );

                // Loop number of Dof per edge
                for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
                {
                    auto temp = this->localToGlobal( iElAd,
                                                     element_type::numVertices*fe_type::nDofPerVertex +
                                                     iEdEl * fe_type::nDofPerEdge + l,
                                                     c );
                    _M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>(temp),boost::get<1>(temp),boost::get<2>(temp),
                                                                                element_type::numVertices*fe_type::nDofPerVertex +
                                                                                iEdEl * fe_type::nDofPerEdge + l );
                }
            }
        }
    }


    template<typename FaceIterator>
    void addFaceBoundaryDof( FaceIterator __face_it, uint16_type& lc )
    {
        addFaceBoundaryDof( __face_it, lc, mpl::bool_<face_type::numFaces*fe_type::nDofPerFace>() );
    }
    template<typename FaceIterator>
    void addFaceBoundaryDof( FaceIterator /*__face_it*/, uint16_type& /*lc*/, mpl::bool_<false> )
    {
    }
    template<typename FaceIterator>
    void addFaceBoundaryDof( FaceIterator __face_it, uint16_type& lc, mpl::bool_<true> )
    {
        // id of the element adjacent to the face
        // \warning NEED TO INVESTIGATE THIS
        size_type iElAd = __face_it->ad_first();
        FEEL_ASSERT( iElAd != invalid_size_type_value )( __face_it->id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );

        // local id of the face in its adjacent element
        uint16_type iFaEl = __face_it->pos_first();
        FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");
#if !defined(NDEBUG)
        Debug( 5005 ) << " local face id : " << iFaEl << "\n";
#endif
        size_type nVerticesAndEdgeF = ( face_type::numVertices * fe_type::nDofPerVertex +
                                        face_type::numEdges * fe_type::nDofPerEdge);
        size_type ndofF = ( face_type::numVertices * fe_type::nDofPerVertex +
                            face_type::numEdges * fe_type::nDofPerEdge +
                            face_type::numFaces * fe_type::nDofPerFace );

        const int ncdof = is_product?nComponents:1;
        for( int c = 0; c < ncdof; ++c )
        {
            uint16_type lcc=nVerticesAndEdgeF+c*ndofF;
            // Loop on number of Dof per face
            for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l )
            {
                auto temp = this->localToGlobal( iElAd,
                                                 element_type::numVertices*fe_type::nDofPerVertex +
                                                 element_type::numEdges*fe_type::nDofPerEdge +
                                                 iFaEl * fe_type::nDofPerFace + l,
                                                 c );
                _M_face_l2g[ __face_it->id()][ lcc++ ] = boost::make_tuple( boost::get<0>(temp),boost::get<1>(temp),boost::get<2>(temp),
                                                                            element_type::numVertices*fe_type::nDofPerVertex +
                                                                            element_type::numEdges*fe_type::nDofPerEdge +
                                                                            iFaEl * fe_type::nDofPerFace + l );
            }
        }
    }

    /**
     * @brief Checks if the dofs associated with entity_id N are continuous
     *
     * \param M A mesh
     */
    template<int N>
01441     void checkDofEntity (mesh_type& M)
    {
        Feel::detail::ignore_unused_variable_warning(M);
        if ( !is_scalar )
            return;

#if !defined(NDEBUG)

        using namespace Feel;

        gm_ptrtype _M_gm_ptr = M.gm();

        fe_type _M_basis;

        //value_type tol = value_type(100.0)*type_traits<value_type>::epsilon();
        value_type tol = value_type(10.0)*type_traits<double>::epsilon();

        bool global_signs_good = 1;

        std::vector<size_type> bad_dof;

        for (uint16_type gDof = 0; gDof < nDof(); ++gDof)
            {
                uint16_type _numEntities = _M_dof2elt[gDof].size();
                uint16_type _ent = _M_dof2elt[gDof].begin()->get<3>();

                if ( _numEntities > 1 && _ent == mpl::int_<N>() )
                    {
                        bool signs_good = 1;

                        std::vector< ublas::vector<value_type> > basis_eval;

                        std::vector< points_type > real_coordinates;

                        ldof_const_iterator __ldofit = _M_dof2elt[gDof].begin();
                        ldof_const_iterator __ldofen = _M_dof2elt[gDof].end();

                        while ( __ldofit != __ldofen )
                            {
                                size_type entity_element_id = __ldofit->get<0>();
                                uint16_type entity_local_dof_id = __ldofit->get<1>();
                                uint16_type entity_local_id = __ldofit->get<2>();

                                PointSetMapped<element_type, convex_type, nOrder> test( M.element(entity_element_id) );

                                points_type Gt = test.pointsBySubEntity(N, entity_local_id);

                                real_coordinates.push_back( test.pointsBySubEntity(N, entity_local_id, 0, 1) );

                                int sign = boost::get<1>(localToGlobal( entity_element_id, entity_local_dof_id ));

                                basis_eval.push_back( value_type(sign)*ublas::row( _M_basis.evaluate(Gt), entity_local_dof_id ) );

                                ++__ldofit;
                            }

                        for (uint16_type i=1; i < _numEntities; i++)
                            {

                                FEEL_ASSERT( ublas::norm_inf( real_coordinates[i] - real_coordinates[0] ) < tol  )
                                    ( gDof )
                                    ( real_coordinates[0] )
                                    ( real_coordinates[i] ).error( "Reference points aren't being mapped to the same real one's" );

                                if ( ublas::norm_inf(basis_eval[i] - basis_eval[0]) > tol )
                                    {
                                        signs_good = 0;
                                        global_signs_good = 0;
                                    }
                            }

                        basis_eval.empty();
                        real_coordinates.empty();

                        if (signs_good == 0)
                            bad_dof.push_back(gDof);
                    }
            }

        if ( !bad_dof.empty() )
            {
                for (uint16_type i = 0; i < bad_dof.size(); ++i)
                    Warning() << bad_dof[i] << "\n";
                if (mpl::int_<N>() == 1)
                    Warning() << "Edges: ";
                else
                    Warning() << "Faces: ";
                Warning() << "Bad dof signs. \n";
            }
#endif
    }

    void checkDofContinuity( mesh_type& /*mesh*/, mpl::int_<1> ) {}

    void checkDofContinuity( mesh_type& mesh, mpl::int_<2> ) { checkDofEntity<1>(mesh); }

    void checkDofContinuity( mesh_type& mesh, mpl::int_<3> )
    {
        checkDofContinuity(mesh, mpl::int_<2>() );
        checkDofEntity<2>(mesh);
    }

    void generateFacePermutations ( mesh_type& /*mesh*/, mpl::bool_<false> ) {}

    void generateFacePermutations ( mesh_type& mesh, mpl::bool_<true> )
    {
        element_type _elt = *mesh.beginElement();
        PointSetMapped<element_type, convex_type, nOrder> pts( _elt );

        for (uint16_type i = 2; i < face_permutation_type::N_PERMUTATIONS; i++)
            vector_permutation[face_permutation_type(i)] = pts.getVectorPermutation(face_permutation_type(i));
    }
    void generateDofPoints( mesh_type& M );
    void generatePeriodicDofPoints( mesh_type& M, periodic_element_list_type const& periodic_elements, dof_points_type& periodic_dof_points );
private:

    fe_ptrtype _M_fe;

    reference_convex_type _M_convex_ref;

    size_type _M_n_el;
    std::vector<bool> _M_elt_done;

    uint16_type _M_n_dof_per_face_on_bdy;
    uint16_type _M_n_dof_per_face;

    Container _M_el_l2g;
    Container_fromface _M_face_l2g;

    dof_element_type _M_dof2elt;

    dof_map_type map_gdof;

    std::map<face_permutation_type, permutation_vector_type> vector_permutation;

    /**
     * map between local(on processor) dof and global dofs
     */
01579     std::map<size_type, std::pair<size_type, size_type> > _M_local2global;


    face_sign_info_type _M_face_sign;

    /**
     * for each dof, store the set of processor it belongs to
     */
01587     dof_procset_type _M_dof_procset;

    /**
     * coordinates of the nodal dofs
     */
01592     dof_points_type M_dof_points;

    std::vector<boost::tuple<size_type, uint16_type, size_type> > M_dof_indices;

    periodicity_type M_periodicity;

    /// a view of the dof container
01599     dof_container_type M_dof_view;

    typedef typename std::vector<localglobal_indices_type,Eigen::aligned_allocator<localglobal_indices_type> > vector_indices_type;


    vector_indices_type M_locglob_indices;
    vector_indices_type M_locglob_signs;
public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

template<typename MeshType, typename FEType, typename PeriodicityType>
const uint16_type DofTable<MeshType, FEType, PeriodicityType>::nComponents;

template<typename MeshType, typename FEType, typename PeriodicityType>
01614 DofTable<MeshType, FEType, PeriodicityType>::DofTable( mesh_type& mesh, fe_ptrtype const& _fe, periodicity_type const& periodicity )
    :
    super(),
    _M_fe( _fe ),
    _M_n_el( invalid_size_type_value ),
    _M_n_dof_per_face_on_bdy( invalid_uint16_type_value ),
    _M_n_dof_per_face( invalid_uint16_type_value ),
    _M_el_l2g(),
    _M_face_l2g(),
    map_gdof(),
    _M_face_sign(),
    M_dof_indices(),
    M_periodicity( periodicity )
{
    Debug( 5015 ) << "[dof] is_periodic = " << is_periodic << "\n";
    size_type start_next_free_dof = 0;
    if ( is_periodic )
        start_next_free_dof = buildPeriodicDofMap( mesh );
    buildDofMap( mesh, start_next_free_dof );
    buildBoundaryDofMap( mesh );
}

template<typename MeshType, typename FEType, typename PeriodicityType>
01637 DofTable<MeshType, FEType, PeriodicityType>::DofTable( fe_ptrtype const& _fe, periodicity_type const& periodicity )
    :
    super(),
    _M_fe( _fe ),
    _M_n_el( 0 ),
    _M_n_dof_per_face_on_bdy( invalid_uint16_type_value ),
    _M_n_dof_per_face( invalid_uint16_type_value ),
    _M_el_l2g(),
    _M_face_l2g(),
    map_gdof(),
    _M_face_sign(),
    M_dof_indices(),
    M_periodicity( periodicity )

{
}

template<typename MeshType, typename FEType, typename PeriodicityType>
01655 DofTable<MeshType, FEType, PeriodicityType>::DofTable( const DofTable<MeshType, FEType, PeriodicityType> & dof2 )
    :
    super( dof2 ),
    _M_fe( dof2._M_fe ),
    _M_n_el( dof2._M_n_el ),
    _M_n_dof_per_face_on_bdy( dof2._M_n_dof_per_face_on_bdy ),
    _M_n_dof_per_face( dof2._M_n_dof_per_face ),
    _M_el_l2g( dof2._M_el_l2g ),
    _M_face_l2g( dof2._M_face_l2g ),
    map_gdof( dof2.map_gdof ),
    _M_face_sign( dof2._M_face_sign),
    M_dof_indices( dof2.M_dof_indices ),
    M_periodicity( dof2.M_periodicity )
{
}

template<typename MeshType, typename FEType, typename PeriodicityType>
void
01673 DofTable<MeshType, FEType, PeriodicityType>::showMe() const
{
    Log()  << " Degree of Freedom (DofTable) Object" << "\n";
    //if ( verbose )
    {
        Log() <<  "* nDof = " << this->nLocalDof() << "\n";
        Log()  << "************************************************************" << "\n";
        Log()  << "           Local to Global DOF table" << "\n";
        Log()  << "************************************************************" << "\n";
        Log()  << "Element Id    Loc. N.    Global N.   Sign#    Element Id   Loc. N.  Global N.  Sign" << "\n";

        for ( size_type i = 0; i < _M_n_el;++i )
            {

                for ( size_type j = 0; j < nDofPerElement;++j )
                    {

                        Log()<< "elt id " << i << " : "
                                   << "(local/global/sign dof : " << j << " : "
                                   << boost::get<0>(localToGlobal( i  , j )) << " : "
                                   << boost::get<1>(localToGlobal( i  , j )) << "\n";
                    }

            }
        Log()  << "\n";

        Log()  << "************************************************************" << "\n";
        Log()  << " Boundary  Local to Global DOF table" << "\n";
        Log()  << "************************************************************" << "\n";
        typedef typename Container_fromface::const_iterator const_iterator;
        const_iterator it = _M_face_l2g.begin();
        const_iterator en = _M_face_l2g.end();
        for(size_type f = 0;it!=en;++it,++f)
            {
                std::ostringstream ostr;
                ostr  << "face id " << f << " : ";
                typedef typename Container_fromface::template subarray<1>::type::const_iterator const_iterator2;
                const_iterator2 it2 = it->begin();
                const_iterator2 en2 = it->end();
                for(size_type l = 0;it2!=en2;++it2,++l)
                    {
                        ostr << "(local/global/sign dof : " << l << " : "
                             << boost::get<0>(*it2)<< " : "
                             << boost::get<1>(*it2) << "\n";
                    }
                Log() << ostr.str() << "\n";
            }
    }

}

template<typename MeshType, typename FEType, typename PeriodicityType>
void
DofTable<MeshType, FEType, PeriodicityType>::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> )
{
    // store the element and local dof id for further
    // reference when inserting the associated global dof
    // id of the element adjacent to the face

    size_type iElAd = __face.ad_first();
    FEEL_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[periodic]invalid face/element in face" );
    Feel::detail::ignore_unused_variable_warning(iElAd);

    // local id of the face in its adjacent element
    uint16_type iFaEl = __face.pos_first();
    FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");

    // loop on face vertices
    for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
        {
            // local vertex number (in element)
            uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
            Feel::detail::ignore_unused_variable_warning(iVeEl);

            FEEL_ASSERT( iVeEl != invalid_uint16_type_value ).error( "invalid local dof" );

            // Loop number of Dof per vertex
            for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
                {
                    uint16_type lid = iVeEl * fe_type::nDofPerVertex + l;
                    //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
                    const size_type gDof = ( __elt.point( iVeEl ).id() ) * fe_type::nDofPerVertex + l;

                    Debug( 5015 ) << "add vertex periodic doc " << next_free_dof << " in element " << __elt.id() << " lid = " << lid << "\n";
                    size_type dof_id = next_free_dof;
                    // next_free_dof might be incremented if a new dof is created
                    bool inserted = this->insertDof( __elt.id(), lid, iVeEl, boost::make_tuple(0, 0, gDof), 0, next_free_dof, 1, true );
                    Debug( 5015 ) << "vertex periodic dof inserted : " << inserted << "\n";

                    // add the pair (elt, lid) to the map associated
                    // with dof_id, one dof can be shared by several
                    // elements
                    dof_id = boost::get<0>(localToGlobal( __elt.id(), lid, 0 ));
                    periodic_dof[tag].insert( std::make_pair( dof_id, boost::make_tuple( __elt.id(), lid, gDof, 0 ) ) );

                    Debug( 5015 ) << "added vertex periodic dof " <<  __elt.id() << ", " <<  lid << ", " << boost::get<0>(localToGlobal( __elt.id(), lid, 0 )) << "\n";
                }

        }

}

template<typename MeshType, typename FEType, typename PeriodicityType>
void
DofTable<MeshType, FEType, PeriodicityType>::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> )
{
#if 0
    // id of the element adjacent to the face
    // \warning NEED TO INVESTIGATE THIS
    size_type iElAd = __face.ad_first();
    FEEL_ASSERT( iElAd != invalid_size_type_value )
        ( __face.id() ).error( "[DofTable::buildBoundaryDof] invalid face/element in face" );
#endif // 0

    // local id of the face in its adjacent element
    uint16_type iFaEl = __face.pos_first();
    FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");
#if !defined(NDEBUG)
    Debug( 5005 ) << " local face id : " << iFaEl << "\n";
#endif
    // Loop number of DofTable per edge
    for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
        {
            uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
            //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
            const size_type gDof = ( __elt.edge( iFaEl ).id() ) * fe_type::nDofPerEdge + l;

            Debug( 5015 ) << "add edge periodic dof " << next_free_dof << " in element " << __elt.id() << " lid = " << lid << "\n";
            size_type dof_id = next_free_dof;
            // next_free_dof might be incremented if a new dof is created
            bool inserted = this->insertDof( __elt.id(), lid, iFaEl, boost::make_tuple(1, 0, gDof), 0, next_free_dof, 1, true );
            Debug( 5015 ) << "edge periodic dof inserted (1 or 0) : " << inserted << "\n";

            // add the pair (elt, lid) to the map associated
            // with dof_id, one dof can be shared by several
            // elements
            dof_id = boost::get<0>(localToGlobal( __elt.id(), lid, 0 ));
            periodic_dof[tag].insert( std::make_pair( dof_id, boost::make_tuple( __elt.id(), lid, gDof, 1 ) ) );

            Debug( 5015 ) << "added edge periodic dof " <<  __elt.id() << ", " <<  lid << ", " << boost::get<0>(localToGlobal( __elt.id(), lid, 0 )) << "\n";

        }
}
template<typename MeshType, typename FEType, typename PeriodicityType>
void
DofTable<MeshType, FEType, PeriodicityType>::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> )
{
    //BOOST_STATIC_ASSERT( face_type::numEdges );

    // id of the element adjacent to the face
    // \warning NEED TO INVESTIGATE THIS
    size_type iElAd = __face.ad_first();
    FEEL_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );

    // local id of the face in its adjacent element
    uint16_type iFaEl = __face.pos_first();
    FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");
#if !defined(NDEBUG)
    Debug( 5005 ) << " local face id : " << iFaEl << "\n";
#endif

    // loop on face vertices
    for ( uint16_type iEdFa = 0; iEdFa < face_type::numEdges; ++iEdFa )
        {
            // local edge number (in element)
            uint16_type iEdEl = element_type::fToE( iFaEl, iEdFa );

            FEEL_ASSERT( iEdEl != invalid_uint16_type_value ).error( "invalid local dof" );

            // Loop number of Dof per edge
            for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
                {

                    uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
                    //const size_type gDof = global_shift + ( __elt.point( i ).id() ) * fe_type::nDofPerVertex + l;
                    size_type gDof = __elt.edge( l ).id() * fe_type::nDofPerEdge;
                    if ( __elt.edgePermutation(l).value()  == edge_permutation_type::IDENTITY )
                        {
                            gDof += l ; // both nodal and modal case
                        }
                    else if ( __elt.edgePermutation(l).value()  == edge_permutation_type::REVERSE_PERMUTATION)
                        {
                            gDof += fe_type::nDofPerEdge - 1 - l ;
                        }

                    Debug( 5015 ) << "add periodic doc " << next_free_dof << " in element " << __elt.id() << " lid = " << lid << "\n";
                    size_type dof_id = next_free_dof;
                    // next_free_dof might be incremented if a new dof is created
                    bool inserted = this->insertDof( __elt.id(), lid, l, boost::make_tuple(1, 0, gDof), 0, next_free_dof, 1, true );
                    Debug( 5015 ) << "periodic dof inserted : " << inserted << "\n";

                    // add the pair (elt, lid) to the map associated
                    // with dof_id, one dof can be shared by several
                    // elements
                    dof_id = boost::get<0>(localToGlobal( __elt.id(), lid, 0 ));
                    periodic_dof[tag].insert( std::make_pair( dof_id, boost::make_tuple( __elt.id(), lid, gDof, 1 ) ) );

                    Debug( 5015 ) << "added " <<  __elt.id() << ", " <<  lid << ", " << boost::get<0>(localToGlobal( __elt.id(), lid, 0 )) << "\n";

                }
        }

}
template<typename MeshType, typename FEType, typename PeriodicityType>
void
DofTable<MeshType, FEType, PeriodicityType>::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> )
{
#if 0
    // id of the element adjacent to the face
    // \warning NEED TO INVESTIGATE THIS
    size_type iElAd = __face.ad_first();
    FEEL_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[Dof::buildBoundaryDof] invalid face/element in face" );

    // local id of the face in its adjacent element
    uint16_type iFaEl = __face.pos_first();
    FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");
#if !defined(NDEBUG)
    Debug( 5005 ) << " local face id : " << iFaEl << "\n";
#endif
    // Loop on number of Dof per face
    for ( uint16_type l = 0; l < fe_type::nDofPerFace; ++l )
        {
            auto temp = this->localToGlobal( iElAd,
                                             element_type::numVertices*fe_type::nDofPerVertex +
                                             element_type::numEdges*fe_type::nDofPerEdge +
                                             iFaEl * fe_type::nDofPerFace + l,
                                             c );
            _M_face_l2g[ __face_it->id()][ lc++ ] = boost::make_tuple( boost::get<0>(temp),boost::get<1>(temp),boost::get<2>(temp),
                                                                       element_type::numVertices*fe_type::nDofPerVertex +
                                                                       element_type::numEdges*fe_type::nDofPerEdge +
                                                                       iFaEl * fe_type::nDofPerFace + l );
        }
#endif // 0
}
template<typename MeshType, typename FEType, typename PeriodicityType>
void
01931 DofTable<MeshType, FEType, PeriodicityType>::initDofMap( mesh_type& M )
{
    _M_n_el = M.numElements();

    size_type nldof =
        fe_type::nDofPerVolume * element_type::numVolumes +
        fe_type::nDofPerFace * element_type::numGeometricFaces +
        fe_type::nDofPerEdge * element_type::numEdges +
        fe_type::nDofPerVertex * element_type::numVertices;

    Debug(5015) << "==============================\n";
    Debug(5015) << "[initDofMap]\n";
    Debug(5015) << "nldof                   = "  << int(nldof) << "\n";
    Debug(5015) << "fe_type::nLocalDof     = "  << int(fe_type::nLocalDof) << "\n";
    Debug(5015) << "fe_type::nDofPerVolume = "  << int(fe_type::nDofPerVolume) << "\n";
    Debug(5015) << "fe_type::nDofPerFace   = "  << int(fe_type::nDofPerFace) << "\n";
    Debug(5015) << "fe_type::nDofPerEdge   = "  << int(fe_type::nDofPerEdge) << "\n";
    Debug(5015) << "fe_type::nDofPerVertex = "  << int(fe_type::nDofPerVertex) << "\n";
    Debug(5015) << "element_type::numVolumes= "  << int(element_type::numVolumes) << "\n";
    Debug(5015) << "element_type::numFaces= "    << int(element_type::numFaces) << "\n";
    Debug(5015) << "element_type::numEdges= "    << int(element_type::numEdges) << "\n";
    Debug(5015) << "element_type::numVertices= " << int(element_type::numVertices) << "\n";
    Debug(5015) << "==============================\n";

    FEEL_ASSERT( nldof == fe_type::nLocalDof )
        ( nldof )
        ( fe_type::nLocalDof ).error( "Something wrong in FE specification" ) ;

    // initialize the local to global map and fill it with invalid
    // values that will allow to check whether we have a new dof or
    // not when building the table
    const size_type nV = M.numElements();
    int ntldof = is_product?nComponents*nldof:nldof;
    _M_el_l2g.resize( boost::extents[nV][ntldof] );
    M_locglob_indices.resize( nV );
    M_locglob_signs.resize( nV );
    typedef Container::index index;
    for ( index i1 = 0; i1 < index(nV); ++i1 )
        for ( index i2 = 0; i2 < index(ntldof); ++i2 )
            _M_el_l2g[i1][i2] = boost::make_tuple(invalid_size_type_value,0,false); // 0 is the invalid value for the sign !

    _M_face_sign = ublas::scalar_vector<bool>(M.numFaces(), false);
    const bool doperm = (((Shape == SHAPE_TETRA) && (nOrder > 2) ) || ((Shape == SHAPE_HEXA) && (nOrder > 1) ));
    Debug( 5005 ) << "generateFacePermutations: " << doperm << "\n";
    generateFacePermutations( M, mpl::bool_<doperm>() );
}
template<typename MeshType, typename FEType, typename PeriodicityType>
size_type
01979 DofTable<MeshType, FEType, PeriodicityType>::buildPeriodicDofMap( mesh_type& M )
{
    size_type nldof =
        fe_type::nDofPerVolume * element_type::numVolumes +
        fe_type::nDofPerFace * element_type::numGeometricFaces +
        fe_type::nDofPerEdge * element_type::numEdges +
        fe_type::nDofPerVertex * element_type::numVertices;

    FEEL_ASSERT( nldof == fe_type::nLocalDof )
        ( nldof )
        ( fe_type::nLocalDof ).error( "Something wrong in FE specification" ) ;

    typedef Container::index index;

    const size_type n_proc  = M.comm().size();

    //! list of elements which have a periodic face Tag2
    periodic_element_list_type periodic_elements;

    for (size_type processor=0; processor<n_proc; processor++)
        {
            // compute the number of dof on current processor
            element_const_iterator it_elt = M.beginElementWithProcessId( processor );
            element_const_iterator en_elt = M.endElementWithProcessId( processor );
            size_type n_elts = std::distance( it_elt, en_elt);
            Debug( 5015 ) << "[buildDofMap] n_elts =  " << n_elts << " on processor " << processor << "\n";
            //this->_M_first_df[processor] = next_free_dof;

            it_elt = M.beginElementWithProcessId( processor );

            it_elt = M.beginElementWithProcessId( processor );
            Debug( 5015 ) << "[buildDofMap] starting with elt " << it_elt->id() << "\n";
            for (;it_elt!=en_elt; ++it_elt )
                {
                    element_type const& __elt = *it_elt;
                    //Debug( 5015 ) << "next_free_dof " << next_free_dof  << "\n";
                    //Debug( 5015 ) << "current dof " << dofIndex( next_free_dof ) << "\n";

                    typename element_type::face_const_iterator it, en;
                    boost::tie( it, en ) = it_elt->faces();

                    bool found_periodic_face_in_element = false;
                    for( ;it != en; ++it )
                        {
                            if ( (*it)->marker().value() == periodicity_type::tag2 ||
                                 (*it)->marker().value() == periodicity_type::tag1 )
                                {
                                    // store the element reference for the end, the associated
                                    // dof on the periodic face is in fact already taken care of.
                                    // the "internal" dof or on not periodic face will be added
                                    periodic_elements.push_back( boost::make_tuple( boost::addressof(__elt), *it ) );
                                    found_periodic_face_in_element = true;
                                    break;
                                }
                        }
                }
        }
    Debug(5015) << "[buildPeriodicDofMap] built periodic_elements " << periodic_elements.size() << "\n";
    std::map<size_type,periodic_dof_map_type> periodic_dof;
    /*
     * Generate the periodic dof, assign a gid to the tag1 dof and set
     * the tag2 dof to invalid_size_type_value for now.
     */
    periodic_element_list_iterator it_periodic = periodic_elements.begin();
    periodic_element_list_iterator en_periodic = periodic_elements.end();
    size_type next_free_dof = 0;
    while( it_periodic != en_periodic )
        {
            element_type const& __elt = *it_periodic->template get<0>();
            face_type const& __face = *it_periodic->template get<1>();
            if ( __face.marker().value() == periodicity_type::tag1 )
                {
                    addVertexPeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
                    addEdgePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
                    addFacePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
                }

            ++it_periodic;
        }
    it_periodic = periodic_elements.begin();
    while( it_periodic != en_periodic )
        {
            element_type const& __elt = *it_periodic->template get<0>();
            face_type const& __face = *it_periodic->template get<1>();
            if ( __face.marker().value() == periodicity_type::tag2 )
                {
                    addVertexPeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
                    addEdgePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
                    addFacePeriodicDof( __elt, __face, next_free_dof, periodic_dof, __face.marker().value() );
                }

            ++it_periodic;
        }

    Debug( 5015 ) << "[periodic dof table] next_free_dof : " << next_free_dof << "\n";
    Debug( 5015 ) << "[periodic dof table] number of periodic dof : " << periodic_dof[periodicity_type::tag1].size() << "\n";

    dof_points_type periodic_dof_points( next_free_dof );
    generatePeriodicDofPoints( M, periodic_elements, periodic_dof_points );

    Debug( 5015 ) << "[periodic dof table] generated dof points\n";
    Debug( 5015 ) << "[periodic dof table] start matching the dof points\n";

    size_type max_gid = 0;
    std::pair<size_type,periodic_dof_type> dof;
    BOOST_FOREACH( dof, periodic_dof[periodicity_type::tag1] )
        {
            size_type gid = dof.first;
            max_gid = (max_gid > gid)?max_gid:gid;
        }

    size_type max_gid2 = 0;
    std::pair<size_type,periodic_dof_type> dof2;
    BOOST_FOREACH( dof2, periodic_dof[periodicity_type::tag2] )
        {
            size_type gid2 = dof2.first;
            FEEL_ASSERT( gid2 > max_gid )( gid2 )( max_gid ).error( "invalid dof index" );
            max_gid2 = (max_gid2 > gid2)?max_gid2:gid2;
        }
    FEEL_ASSERT( (max_gid+1) == (max_gid2+1-(max_gid+1)) )( max_gid )( max_gid2 ).error( "[periodic] invalid periodic setup" );

    std::vector<bool> periodic_dof_done( max_gid+1 );
    std::fill( periodic_dof_done.begin(), periodic_dof_done.end(), false );

    BOOST_FOREACH( dof, periodic_dof[periodicity_type::tag1] )
        {
            size_type gid = dof.first;
            if ( periodic_dof_done[gid] )
                continue;
            node_type x1 = periodic_dof_points[gid].template get<0>();
            bool match = false;
            typename periodic_dof_map_type::iterator it_dof2 = periodic_dof[periodicity_type::tag2].begin();
            typename periodic_dof_map_type::iterator en_dof2 = periodic_dof[periodicity_type::tag2].end();

            for( ; it_dof2 != en_dof2; ++ it_dof2 )
                {
                    size_type gid2 = it_dof2->first;
                    FEEL_ASSERT( gid2 < next_free_dof )( gid )( gid2 )( next_free_dof ).error( "[periodic] invalid dof id" );
                    node_type x2 = periodic_dof_points[gid2].template get<0>();
                    //FEEL_ASSERT( math::abs( x2[0]-M_periodicity.translation()[0]) < 1e-10 )
                    //( x1 )( x2 )( M_periodicity.translation() ).error( "[periodic] invalid periodic setup");
                }
            it_dof2 = periodic_dof[periodicity_type::tag2].begin();
            size_type corresponding_gid = invalid_size_type_value;
            for( ; it_dof2 != en_dof2; ++ it_dof2 )
                {
                    size_type gid2 = it_dof2->first;
                    FEEL_ASSERT( gid2 < next_free_dof )( gid )( gid2 )( next_free_dof ).error( "[periodic] invalid dof id" );
                    node_type x2 = periodic_dof_points[gid2].template get<0>();
                    FEEL_ASSERT( ( x1.size() == x2.size() ) &&
                                 ( x1.size() == M_periodicity.translation().size() ) )
                        ( gid )( dof.second.template get<0>() )( dof.second.template get<1>())
                        ( gid2 )( it_dof2->second.template get<0>() )( it_dof2->second.template get<1>())
                        ( x1 )( x2 )( M_periodicity.translation() ).error( "invalid point size" );
                    if ( ublas::norm_2( x1-(x2-M_periodicity.translation()) ) < 1e-10 )
                        {
                            // loop on each pair (element, lid) which
                            // has a global id gid2 and set it to gid
                            corresponding_gid = gid2;
                            match = true;
                            break;
                        }
                }
            // if we have --- actually we must have one --- a match, remove the
            // iterator from dof2 to quicken the search for the next dof1 match
            if ( match )
                {

                    it_dof2 = periodic_dof[periodicity_type::tag2].lower_bound( corresponding_gid );
                    en_dof2 = periodic_dof[periodicity_type::tag2].upper_bound( corresponding_gid );
                    Debug( 5015 ) << "distance = " << std::distance( it_dof2, en_dof2 ) << "\n";
                    while( it_dof2 != en_dof2 )
                        {

                            size_type ie = it_dof2->second.template get<0>();
                            size_type lid = it_dof2->second.template get<1>();
                            size_type gDof = it_dof2->second.template get<2>();
                            uint16_type dof2_type = it_dof2->second.template get<3>();
                            uint16_type dof1_type = dof.second.template get<3>();

                            FEEL_ASSERT( dof1_type == dof2_type )
                                ( gid )( it_dof2->first )( gDof )( lid )( ie )
                                ( dof1_type )( dof2_type ).error ( "invalid dof" );

                            Debug( 5015 ) << "link " <<  boost::get<0>( _M_el_l2g[ ie][ lid ] )  << " -> " << gid << "\n";

                            // gid is given by dof1
                            _M_el_l2g[ ie][ lid ] = boost::make_tuple(gid, 1, true );


                            // warning: must modify the data structure that allows to
                            // generate unique global dof ids
                            FEEL_ASSERT( (map_gdof[  boost::make_tuple(dof2_type, 0, gDof) ] == corresponding_gid ) ||
                                         ( map_gdof[ boost::make_tuple(dof2_type, 0, gDof) ] == gid ) )
                                ( corresponding_gid )( dof2_type )( gDof )( gid )
                                ( map_gdof[ boost::make_tuple(dof2_type, 0, gDof) ] ).error ("invalid gid" );

                            Debug( 5015 ) << "link mapgdof " <<   map_gdof[ boost::make_tuple(dof2_type, 0, gDof) ]  << " -> " << gid << "\n";
                            map_gdof[ boost::make_tuple(dof2_type, 0, gDof) ] = gid;

                            FEEL_ASSERT( map_gdof[ boost::make_tuple(dof2_type, 0, gDof) ] == gid )
                                ( corresponding_gid )( dof2_type )( gDof )( gid )
                                ( map_gdof[ boost::make_tuple(dof2_type, 0, gDof) ]) .error ("invalid gid" );

                            ++it_dof2;
                        }
                    it_dof2 = periodic_dof[periodicity_type::tag2].lower_bound( corresponding_gid );
                    periodic_dof[periodicity_type::tag2].erase( it_dof2, en_dof2 );
                    periodic_dof_done[gid] =  true;
                }
            else
                {
                    // we have a problem, no match was found, this should not happen
                    Debug( 5015 ) << "[periodic] invalid point/dof matching\n";
                    Debug( 5015 ) << "[periodic] n = " << x1 << "\n";
                }

        }
    Debug( 5015 ) << "[periodic dof table] done matching the dof points\n";
    Debug( 5015 ) << "[periodic dof table] is empty : " << periodic_dof[periodicity_type::tag2].empty() << "\n";

    // ensure that periodic_dof[periodicity_type::tag2] is empty
    if ( !periodic_dof[periodicity_type::tag2].empty() )
        {
            Debug( 5015 ) << "[periodic] periodic conditions not set properly, some periodic dof were not assigned\n";
            typename periodic_dof_map_type::iterator it_dof2 = periodic_dof[periodicity_type::tag2].begin();
            typename periodic_dof_map_type::iterator en_dof2 = periodic_dof[periodicity_type::tag2].end();
            while( it_dof2 != en_dof2 )
                {

                    size_type ie = it_dof2->second.template get<0>();
                    size_type lid = it_dof2->second.template get<1>();

                    Debug( 5015 ) << "[periodic] dof " << it_dof2->first << " not assigned, "
                                  << "x = " << periodic_dof_points[it_dof2->first].template get<0>() << " "
                                  << "elt = " << ie << ", lid= " << lid << "\n";



                    ++it_dof2;
                }
        }
    else
        {
            Debug( 5015 ) << "[periodic] periodic condition done\n";
        }

    return max_gid+1;
}

template<typename MeshType, typename FEType, typename PeriodicityType>
size_type
02231 DofTable<MeshType, FEType, PeriodicityType>::buildLocallyDiscontinuousDofMap( mesh_type& M, size_type start_next_free_dof )
{
    typedef typename continuity_type::template apply<MeshType, self_type> builder;
    return fusion::accumulate( typename continuity_type::discontinuity_markers_type(), start_next_free_dof,  builder(M, *this ) );
}
template<typename MeshType, typename FEType, typename PeriodicityType>
void
02238 DofTable<MeshType, FEType, PeriodicityType>::buildDofMap( mesh_type& M, size_type start_next_free_dof )
{
    if ( !M_dof_indices.empty() )
        {
            generateDofPoints( M );
            return;
        }

    size_type nldof =
        fe_type::nDofPerVolume * element_type::numVolumes +
        fe_type::nDofPerFace * element_type::numGeometricFaces +
        fe_type::nDofPerEdge * element_type::numEdges +
        fe_type::nDofPerVertex * element_type::numVertices;

    FEEL_ASSERT( nldof == fe_type::nLocalDof )
        ( nldof )
        ( fe_type::nLocalDof ).error( "Something wrong in FE specification" ) ;

    typedef Container::index index;

    const size_type n_proc  = M.comm().size();

    //! list of elements which have a periodic face Tag2
    std::list<std::pair<element_type const*, face_type const*> > periodic_elements;

    size_type next_free_dof = start_next_free_dof;
    for (size_type processor=0; processor<n_proc; processor++)
        {
            // compute the number of dof on current processor
            element_const_iterator it_elt = M.beginElementWithProcessId( processor );
            element_const_iterator en_elt = M.endElementWithProcessId( processor );
            size_type n_elts = std::distance( it_elt, en_elt);
            Debug( 5005 ) << "[buildDofMap] n_elts =  " << n_elts << " on processor " << processor << "\n";
            this->_M_first_df[processor] = next_free_dof;
            if ( is_periodic || is_discontinuous_locally )
                this->_M_first_df[processor] =  0;
            it_elt = M.beginElementWithProcessId( processor );

            Debug( 5005 ) << "[buildDofMap] starting with elt " << it_elt->id() << "\n";
            for (;it_elt!=en_elt; ++it_elt )
                {
                    this->addDofFromElement( *it_elt, next_free_dof, processor );
                } // elements loop

    // printing Dof table only in debug mode
#if !defined( NDEBUG )
            const int ncdof = is_product?nComponents:1;
            for( int c = 0; c < ncdof; ++c )
            {
                it_elt = M.beginElementWithProcessId( processor );
                for (;it_elt!=en_elt; ++it_elt )
                {
                    element_type const& __elt = *it_elt;
                    std::ostringstream ostr;
                    ostr << "element " << __elt.id() << " : ";
                    for ( uint16_type l = 0; l < nldof; ++l )
                        {
                            ostr << _M_el_l2g[ __elt.id()][ fe_type::nLocalDof*c + l ] << " ";
                        }
                    Debug( 5005 ) << ostr.str() << "\n";
                }
        }
#endif
    this->_M_last_df[processor] = next_free_dof-1;
        }
    this->_M_n_dofs = next_free_dof;

#if !defined(NDEBUG)
    Debug( 5005 ) << " n global dof " << nDof() << "\n";
    Debug( 5005 ) << " n local dof " << nLocalDof() << "\n";
#endif
    for (size_type processor=0; processor<M.comm().size(); processor++)
        {
            Debug( 5005 ) << "o processor " << processor << "\n";
            Debug( 5005 ) << "  - n dof on proc " << nDofOnProcessor(processor) << "\n";
            Debug( 5005 ) << "  - first dof " << firstDof(processor) << "\n";
            Debug( 5005 ) << "  - last dof " << lastDof(processor) << "\n";
        }

    //if ( is_continuous )
    //checkDofContinuity( M, mpl::int_<fe_type::nDim>() );

    for (size_type processor=0; processor<n_proc; processor++)
    {
        auto it_elt = M.beginElementWithProcessId( processor );
        auto en_elt = M.endElementWithProcessId( processor );
        for( ; it_elt != en_elt; ++it_elt )
        {
            size_type elid= it_elt->id();
            for( int i = 0; i < FEType::nLocalDof; ++i )
            {
                int nc1 = (is_product?nComponents1:1);
                for( int c1 =0; c1 < nc1; ++c1 )
                {
                    int ind = FEType::nLocalDof*c1+i;
                    boost::tie( M_locglob_indices[elid][ind],
                                M_locglob_signs[elid][ind], boost::tuples::ignore) =
                        localToGlobal( elid, i, c1 );
                }
            }
        }
    }

    generateDofPoints( M );
}

template<typename MeshType, typename FEType, typename PeriodicityType>
void
02346 DofTable<MeshType, FEType, PeriodicityType>::buildBoundaryDofMap( mesh_type& M )
{
    size_type nDofF = ( face_type::numVertices * fe_type::nDofPerVertex +
                        face_type::numEdges * fe_type::nDofPerEdge +
                        face_type::numFaces * fe_type::nDofPerFace );
    _M_n_dof_per_face_on_bdy = nDofF;
    Debug( 5005 ) << "vertex dof : " <<  face_type::numVertices * fe_type::nDofPerVertex << "\n";
    Debug( 5005 ) << "edge dof : " <<  face_type::numEdges * fe_type::nDofPerEdge << "\n";
    Debug( 5005 ) << "face dof : " << face_type::numFaces * fe_type::nDofPerFace  << "\n";
    Debug( 5005 ) << "number of Dof on an Element Face : " << nDofF << "\n";

    //
    // Face dof
    //
    typename mesh_type::face_const_iterator __face_it = M.beginFace();
    typename mesh_type::face_const_iterator __face_en = M.endFace();

    const size_type nF = M.faces().size();
    int ntldof = is_product?nComponents*nDofF:nDofF;
    _M_face_l2g.resize( boost::extents[nF][ntldof] );
    typedef Container_fromface::index index;
    global_dof_fromface_type default_dof = boost::make_tuple(invalid_size_type_value,0,false,0);
    for ( index i1 = 0; i1 < index(nF); ++i1 )
        for ( index i2 = 0; i2 < index(ntldof); ++i2 )
            // 0 is the invalid value for the sign !
            _M_face_l2g[i1][i2] = default_dof;

    Debug( 5005 ) << "[buildBoundaryDofMap] nb faces : " << nF << "\n";
    Debug( 5005 ) << "[buildBoundaryDofMap] nb dof faces : " << nDofF*nComponents << "\n";

    __face_it = M.beginFace();
    for ( size_type nf = 0; __face_it != __face_en; ++__face_it, ++nf )
    {
        FEEL_ASSERT( __face_it->isConnectedTo0() )
            ( __face_it->id() )
            ( __face_it->marker() )
            ( __face_it->isConnectedTo0() )
            ( __face_it->isConnectedTo1() ).warn( "[Dof::buildFaceDofMap] face not connected" );

        if ( !__face_it->isConnectedTo0() ) continue;

#if !defined(NDEBUG)
        if (  __face_it->isOnBoundary() )
            Debug( 5005 ) << "[buildBoundaryDofMap] boundary global face id : " << __face_it->id()
                          << " marker: " << __face_it->marker()<< "\n";
        else
            Debug( 5005 ) << "[buildBoundaryDofMap] global face id : " << __face_it->id() << "\n";
#endif
        uint16_type lcVertex = 0;
        uint16_type lcEdge = 0;
        uint16_type lcFace = 0;


        addVertexBoundaryDof( __face_it, lcVertex );
        addEdgeBoundaryDof( __face_it, lcEdge );
        addFaceBoundaryDof( __face_it, lcFace );

    }
#if !defined(NDEBUG)
    for ( index face_id = 0; face_id < index(nF); ++face_id )
        for ( index face_dof_id = 0; face_dof_id < index(ntldof); ++face_dof_id )
            FEEL_ASSERT( boost::get<0>(_M_face_l2g[face_id][face_dof_id]) != invalid_size_type_value )( face_id )( face_dof_id ).warn( "invalid dof table: initialized dof entries" );
#endif
}    // updateBoundaryDof

template<typename MeshType, typename FEType, typename PeriodicityType>
void
DofTable<MeshType, FEType, PeriodicityType>::generateDofPoints(  mesh_type& M )
{
    if ( !M_dof_points.empty() )
        return;
    if ( fe_type::is_modal )
        return;
    Debug( 5005 ) << "[Dof::generateDofPoints] generating dof coordinates\n";
    typedef typename gm_type::template Context<vm::POINT, element_type> gm_context_type;
    typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;

    typedef typename fe_type::template Context<vm::POINT, fe_type, gm_type, element_type> fecontext_type;

    gm_ptrtype gm( new gm_type );
    fe_type fe;
    //
    // Precompute some data in the reference element for
    // geometric mapping and reference finite element
    //
    typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( gm, fe.points() ) );

    //const uint16_type ndofv = fe_type::nDof;

    element_const_iterator it_elt = M.beginElementWithProcessId( M.comm().rank() );
    element_const_iterator en_elt = M.endElementWithProcessId( M.comm().rank() );

    if ( it_elt == en_elt )
        return;

    gm_context_ptrtype __c( new gm_context_type( gm, *it_elt, __geopc ) );

    std::vector<bool> dof_done( nLocalDof() );
    M_dof_points.resize( nLocalDof() );
    std::fill( dof_done.begin(), dof_done.end(), false );
    for ( size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
    {
        __c->update( *it_elt );
        for( uint16_type l =0; l < fe_type::nLocalDof; ++l )
            {
                int ncdof  = is_product?nComponents:1;
                for( uint16_type c1 = 0; c1 < ncdof; ++c1 )
                    {
                        size_type thedof = boost::get<0>(localToGlobal( it_elt->id(), l, c1 ));
                        if ( ( thedof >= firstDof() ) && ( thedof <= lastDof() ) )
                            {
                                // get only the local dof
                                //size_type thedofonproc = thedof - firstDof();
                                thedof -= firstDof();
                                FEEL_ASSERT( thedof < nLocalDof() )( thedof )( nLocalDof() )( firstDof() )( lastDof() )( it_elt->id() )( l )( c1 ).error ( "invalid local dof index");
                                if ( dof_done[ thedof ] == false )
                                    {
                                        std::set<uint16_type> lid;
                                        lid.insert( l );
                                        //M_dof_points[dof_id] = boost::make_tuple( thedof, __c->xReal( l ) );
                                        M_dof_points[thedof] = boost::make_tuple( __c->xReal( l ), firstDof()+thedof, c1 );
                                        dof_done[thedof] = true;
                                        ++dof_id;
                                    }
                            }
                    }
            }
    }
    for ( size_type dof_id = 0; dof_id < nLocalDof() ; ++dof_id )
        {
            FEEL_ASSERT( boost::get<1>(M_dof_points[dof_id]) >= firstDof() &&
                         boost::get<1>(M_dof_points[dof_id]) <= lastDof() )
                ( dof_id )( firstDof() )( lastDof() )( nLocalDof() )
                ( boost::get<1>(M_dof_points[dof_id]) )
                ( boost::get<0>(M_dof_points[dof_id]) ).error( "invalid dof point" );
            FEEL_ASSERT( dof_done[dof_id] == true )( dof_id )( nLocalDof() )( firstDof() )( lastDof() )( fe_type::nDim )( fe_type::nLocalDof ).warn( "invalid dof point" );
        }
    Debug( 5005 ) << "[Dof::generateDofPoints] generating dof coordinates done\n";
}
template<typename MeshType, typename FEType, typename PeriodicityType>
void
DofTable<MeshType, FEType, PeriodicityType>::generatePeriodicDofPoints(  mesh_type& M,
                                                                         periodic_element_list_type const& periodic_elements,
                                                                         dof_points_type& periodic_dof_points )
{
    if ( fe_type::is_modal )
        return;
    Debug( 5005 ) << "[Dof::generateDofPoints] generating dof coordinates\n";
    typedef typename gm_type::template Context<vm::POINT, element_type> gm_context_type;
    typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;

    typedef typename fe_type::template Context<vm::POINT, fe_type, gm_type, element_type> fecontext_type;

    gm_ptrtype gm( new gm_type );
    fe_type fe;
    //
    // Precompute some data in the reference element for
    // geometric mapping and reference finite element
    //
    typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( gm, fe.points() ) );

    //const uint16_type ndofv = fe_type::nDof;

    periodic_element_list_const_iterator it_elt = periodic_elements.begin();
    periodic_element_list_const_iterator en_elt = periodic_elements.end();

    if ( it_elt == en_elt )
        return;

    gm_context_ptrtype __c( new gm_context_type( gm, *it_elt->template get<0>(), __geopc ) );

    std::vector<bool> dof_done( periodic_dof_points.size() );
    std::fill( dof_done.begin(), dof_done.end(), false );
    for ( size_type dof_id = 0; it_elt!=en_elt ; ++it_elt )
    {
        __c->update( *it_elt->template get<0>() );

        face_type const& __face = *it_elt->template get<1>();

        size_type iElAd = __face.ad_first();
        FEEL_ASSERT( iElAd != invalid_size_type_value )( __face.id() ).error( "[periodic]invalid face/element in face" );
        Feel::detail::ignore_unused_variable_warning(iElAd);

        // local id of the face in its adjacent element
        uint16_type iFaEl = __face.pos_first();
        FEEL_ASSERT( iFaEl != invalid_uint16_type_value ).error ("invalid element index in face");

        // loop on face vertices
        for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
        {
            // local vertex number (in element)
            uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
            Feel::detail::ignore_unused_variable_warning(iVeEl);

            FEEL_ASSERT( iVeEl != invalid_uint16_type_value ).error( "invalid local dof" );

            // Loop number of Dof per vertex
            for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
            {
                uint16_type lid = iVeEl * fe_type::nDofPerVertex + l;

                size_type thedof = boost::get<0>(localToGlobal( it_elt->template get<0>()->id(), lid, 0 ));
                FEEL_ASSERT( thedof < dof_done.size() )
                    ( thedof )
                    ( dof_done.size() )
                    ( it_elt->template get<0>()->id() )
                    ( lid ).error ("[generatePeriodicDofPoints] invalid dof id" );
                if ( dof_done[ thedof ] == false )
                {
                    periodic_dof_points[thedof] = boost::make_tuple( __c->xReal( lid ), thedof, 0 );
                    // these tests are problem specific x=0 and x=translation
#if 0
                    if ( __face.marker().value() == periodicity_type::tag1 )
                        FEEL_ASSERT( math::abs( __c->xReal( lid )[0] ) < 1e-10 )( __c->xReal( lid ) ).warn( "[periodic] invalid p[eriodic point tag1");
                    if ( __face.marker().value() == periodicity_type::tag2 )
                        FEEL_ASSERT( math::abs( __c->xReal( lid )[0] - M_periodicity.translation()[0] ) < 1e-10 )
                            ( __c->xReal( lid ) )( M_periodicity.translation()).warn( "[periodic] invalid p[eriodic point tag1");
#endif
                    dof_done[thedof] = true;
                    ++dof_id;
                }
            }
        }
        // loop on edge vertices
        for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
        {
            uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
            size_type thedof = boost::get<0>(localToGlobal( it_elt->template get<0>()->id(), lid, 0 ));
            FEEL_ASSERT( thedof < dof_done.size() )
                ( thedof )
                ( dof_done.size() )
                ( it_elt->template get<0>()->id() )
                ( lid ).error ("[generatePeriodicDofPoints] invalid dof id" );
            if ( dof_done[ thedof ] == false )
            {
                periodic_dof_points[thedof] = boost::make_tuple( __c->xReal( lid ), thedof, 0 );
                // these tests are problem specific x=0 and x=translation
#if 0
                if ( __face.marker().value() == periodicity_type::tag1 )
                    FEEL_ASSERT( math::abs( __c->xReal( lid )[1] +1 ) < 1e-10 )( __c->xReal( lid ) ).warn( "[periodic] invalid p[eriodic point tag1");
                if ( __face.marker().value() == periodicity_type::tag2 )
                    FEEL_ASSERT( math::abs( __c->xReal( lid )[1] - (M_periodicity.translation()[1]-1) ) < 1e-10 )
                        ( __c->xReal( lid ) )( M_periodicity.translation()).warn( "[periodic] invalid p[eriodic point tag1");
#endif
                dof_done[thedof] = true;
                ++dof_id;
            }
        }
    }
    for ( size_type dof_id = 0; dof_id < periodic_dof_points.size() ; ++dof_id )
    {
        FEEL_ASSERT( boost::get<1>(periodic_dof_points[dof_id]) >= 0 &&
                     boost::get<1>(periodic_dof_points[dof_id]) < periodic_dof_points.size() )
            ( dof_id )( periodic_dof_points.size() )
            ( boost::get<1>(periodic_dof_points[dof_id]) )
            ( boost::get<0>(periodic_dof_points[dof_id]) ).error( "invalid dof point" );
        FEEL_ASSERT( dof_done[dof_id] == true )( dof_id ).error( "invalid dof point" );
    }
}

} // namespace Feel
#endif

Generated by  Doxygen 1.6.0   Back to index