gsDofMapper.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsForwardDeclarations.h>
#include <gsCore/gsBoundary.h>
#include <gsCore/gsExport.h>

namespace gismo
{

#define MAPPER_PATCH_DOF(a,b,c) m_dofs[c][m_offset[b]+a]

class GISMO_EXPORT gsDofMapper
{
public:

    gsDofMapper();

    template<class T>
    gsDofMapper(
        const gsMultiBasis<T>         &bases,
        const gsBoundaryConditions<T> &dirichlet,
        int unk = 0
        ) : m_shift(0), m_bshift(0)
    {
      init(bases, dirichlet, unk); //obsolete, one component
    }

    template<class T>
    gsDofMapper(const gsMultiBasis<T> & bases, index_t nComp = 1)
    {
        init(bases, nComp);
    }

    template<class T>
    gsDofMapper(
        std::vector<const gsMultiBasis<T> *> const & bases
        )
    {
        init(bases);
    }


    template<class T>
    gsDofMapper(const gsBasis<T> & basis, index_t nComp = 1)
    {
      initSingle(basis, nComp);
    }

    gsDofMapper(const gsVector<index_t> &patchDofSizes, index_t nComp = 1)
    {
        initPatchDofs(patchDofSizes, nComp);
    }

    template <typename T>
    void init(const gsMultiBasis<T> & bases, index_t nComp = 1);

    template <typename T>
    void init( std::vector<const gsMultiBasis<T> *> const & bases);

    template<class T>
    void init(const gsMultiBasis<T>         &basis,
              const gsBoundaryConditions<T> &dirichlet, int unk = 0);

    void swap(gsDofMapper & other)
    {
        m_dofs  .swap(other.m_dofs);
        m_offset.swap(other.m_offset);

        std::swap(m_shift      , other.m_shift);
        std::swap(m_bshift     , other.m_bshift);
        std::swap(m_numFreeDofs, other.m_numFreeDofs);
        std::swap(m_numElimDofs, other.m_numElimDofs);
        std::swap(m_numCpldDofs, other.m_numCpldDofs);
        std::swap(m_curElimId  , other.m_curElimId);
        std::swap(m_tagged     , other.m_tagged);
    }

private:

    template <typename T>
    void initSingle( const gsBasis<T> & basis, index_t nComp = 1);

    void initPatchDofs(const gsVector<index_t> & patchDofSizes,
                       index_t nComp = 1);

public:

    gsVector<index_t> asVector(index_t comp = 0) const;

    gsVector<index_t> inverseAsVector(index_t comp = 0) const;

    void setMatchingInterfaces(const gsBoxTopology & mp);

    void colapseDofs(index_t k, const gsMatrix<unsigned> & b, index_t comp = 0);

    void matchDof( index_t u, index_t i, index_t v, index_t j, index_t comp = 0);

    void matchDofs(index_t u, const gsMatrix<index_t> & b1,
                   index_t v, const gsMatrix<index_t> & b2,
                           index_t comp = 0);

    void markCoupled(index_t i, index_t k, index_t comp = 0);

    void markTagged(index_t i, index_t k, index_t comp = 0);

    void markCoupledAsTagged();

    // to do: put k at the end
    void markBoundary(index_t k, const gsMatrix<index_t> & boundaryDofs, index_t comp = 0);

    void eliminateDof(index_t i, index_t k, index_t comp = 0);

    void finalize();

    bool isFinalized() const { return m_curElimId>=0; }

    bool isPermutation() const { return static_cast<size_t>(size())==mapSize(); }

    std::ostream& print( std::ostream& os = gsInfo ) const;

    void setIdentity(index_t nPatches, size_t nDofs, size_t nComp = 1);

    void setShift(index_t shift);

    void addShift(index_t shift);

    void permuteFreeDofs(const gsVector<index_t>& permutation, index_t comp = 0);

    index_t firstIndex(index_t comp = 0) const
    { return m_numFreeDofs[comp] + m_numElimDofs[comp] + m_shift; }

    index_t lastIndex() const { return m_shift + freeSize(); }

    void setBoundaryShift(index_t shift);

    void localToGlobal(const gsMatrix<index_t>& locals,
                       index_t patchIndex,
                       gsMatrix<index_t>& globals, 
                               index_t comp = 0) const;

    void localToGlobal2(const gsMatrix<index_t>& locals,
                        index_t patchIndex,
                        gsMatrix<index_t>& globals,
                        index_t & numFree, 
                                index_t comp = 0) const;

    inline index_t freeIndex(index_t i, index_t k = 0, index_t c = 0) const
    {
        GISMO_ASSERT(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return MAPPER_PATCH_DOF(i,k,c);
    }

 index_t componentOf(index_t gl) const
 {
     GISMO_ASSERT(m_curElimId>=0,"finalize() was not called on gsDofMapper");
     //index_t c = 1; // could do some binary search
     //in place
     //while (gl >= m_numFreeDofs[c]+m_numElimDofs[c]) { ++c; }
     //elim
     return (gl<m_numFreeDofs.back() ?
       std::distance(m_numFreeDofs.begin(), std::upper_bound(m_numFreeDofs.begin(), m_numFreeDofs.end(), gl))
             : std::distance(m_numElimDofs.begin(),std::upper_bound(m_numElimDofs.begin(), m_numElimDofs.end(), gl-m_numFreeDofs.back())) ) - 1;
     //while (gl >= m_numFreeDofs[c] + m_shift) { ++c; } return c-1;
 }

    inline index_t index(index_t i, index_t k = 0, index_t c = 0) const
    {
        GISMO_ASSERT(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return MAPPER_PATCH_DOF(i,k,c)+m_shift;
    }

    inline index_t bindex(index_t i, index_t k = 0, index_t c = 0) const
    {
        GISMO_ASSERT(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return MAPPER_PATCH_DOF(i,k,c) - m_numFreeDofs.back()
            //- m_numElimDofs[c]
            + m_bshift;
    }

    bool allFree() const
    { return m_numFreeDofs.back()+m_numElimDofs.back()==m_curElimId; }

    inline index_t cindex(index_t i, index_t k = 0, index_t c = 0) const
    {
        GISMO_ASSERT(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return MAPPER_PATCH_DOF(i,k,c) - m_numFreeDofs[c+1]
          + m_numCpldDofs[c+1];
    }

    inline index_t tindex(index_t i, index_t k = 0, index_t c = 0) const
    {
        GISMO_ASSERT(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return std::distance(m_tagged.begin(),std::lower_bound(m_tagged.begin(),m_tagged.end(),MAPPER_PATCH_DOF(i,k,c)));
    }

    inline index_t global_to_bindex(index_t gl) const
    {
        GISMO_ASSERT( is_boundary_index( gl ),
                      "global_to_bindex(): dof "<<gl<<" is not on the boundary");
        return gl - m_numFreeDofs.back() - m_shift + m_bshift;
    //gl -= m_numFreeDofs.back() + m_shift;
        //const index_t c = std::distance(m_numElimDofs.begin(),
    //    std::upper_bound(m_numElimDofs.begin(), m_numElimDofs.end(),gl)) -1;
        //return gl + m_bshift; // - m_numElimDofs[c]
    }

    inline bool is_free_index(index_t gl) const
    {
      return gl < m_curElimId + m_shift;
    }

    inline bool is_free( index_t i, index_t k = 0, index_t c = 0) const
    { return is_free_index( index(i, k, c) ); }

    inline bool is_boundary_index( index_t gl ) const
    { return gl >= m_numFreeDofs.back() + m_shift; }

    inline bool is_boundary(index_t i, index_t k = 0, index_t c = 0) const
    {return is_boundary_index( index(i, k, c) );}

    inline bool is_coupled( index_t i, index_t k = 0, index_t c = 0) const
    { return  is_coupled_index( index(i, k, c) ); }

    inline bool is_coupled_index(index_t gl) const
    {
      const index_t gc = componentOf(gl);
      const index_t vv = m_numFreeDofs[gc+1] + m_shift;
      return  (gl < vv && // is a free dof and
              (gl + m_numCpldDofs[gc+1] + 1 > vv) );  // is not standard dof
    }

    inline bool is_tagged(index_t i, index_t k = 0, index_t c = 0) const
    { return  is_tagged_index( index(i, k, c) ); }

    inline bool is_tagged_index(index_t gl) const
    {
        return std::binary_search(m_tagged.begin(),m_tagged.end(),gl);
    }

    inline index_t numComponents() const
    { return static_cast<index_t>(m_dofs.size()); }

    inline index_t size() const
    {
        GISMO_ENSURE(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return freeSize() + boundarySize();
    }

    inline index_t size(index_t comp) const
    {
        GISMO_ENSURE(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return m_numFreeDofs[comp+1]-m_numFreeDofs[comp]
          + m_numElimDofs[comp+1]-m_numElimDofs[comp];
    }

    inline index_t freeSize() const
    {
      return m_curElimId;
    }

    inline index_t freeSize(index_t comp) const
    {
      return m_numFreeDofs[comp+1]-m_numFreeDofs[comp] +
        (allFree() ? m_numElimDofs[comp+1]-m_numElimDofs[comp] : 0 );
    }

    const std::vector<index_t> & getTagged() const { return m_tagged; }

    index_t coupledSize() const;

    index_t taggedSize() const;

    inline index_t boundarySize() const
    {
      GISMO_ENSURE(m_curElimId>=0, "finalize() was not called on gsDofMapper");
        return m_numElimDofs.back();
    }

    index_t boundarySizeWithDuplicates() const;

    size_t offset(int k) const {return m_offset[k];}

    size_t numPatches() const {return m_offset.size();}

    size_t mapSize() const
    { return (m_dofs.empty()?0:m_dofs.size() * m_dofs.front().size()); }

    size_t componentsSize() const {return m_dofs.size();}

    size_t patchSize(const index_t k, const index_t c = 0) const
    {
        const size_t k1(k+1);
        GISMO_ASSERT(k1<=numPatches(), "Invalid patch index "<< k <<" >= "<< numPatches() );
        if ( 1==m_offset.size() ) return  m_dofs[c].size();
        else if ( k1==m_offset.size() ) return (m_dofs[c].size() - m_offset.back());
        else return (m_offset[k1]-m_offset[k]);
    }

    size_t totalSize(const index_t c = 0) const
    {
        return m_dofs[c].size();
    }

    void preImage(index_t gl, std::vector<std::pair<index_t,index_t> > & result) const;

    std::pair<index_t,index_t> anyPreImage(index_t gl) const;

    std::vector<std::pair<index_t,index_t> > anyPreImages(index_t comp = 0) const;
    
    std::map<index_t,index_t> inverseOnPatch(const index_t k) const;

    bool indexOnPatch(const index_t gl, const index_t k) const;

    index_t mapIndex(index_t n) const
    {
        return m_dofs[n/m_dofs.front().size()]
            [n%m_dofs.front().size()] + m_shift;
    }

     gsVector<index_t> findBoundary(const index_t k) const;

     gsVector<index_t> findFree(const index_t k) const;

     gsVector<index_t> findCoupled(const index_t k, const index_t j = -1) const;

     gsVector<index_t> findFreeUncoupled(const index_t k) const;

     gsVector<index_t> findTagged(const index_t k) const;

private:

    template<class Predicate, class Iterator>
      static gsVector<index_t> find_impl(Iterator istart, Iterator iend, Predicate pred);

    void finalizeComp(const index_t comp);

    // replace all references to oldIdx by newIdx
    inline void replaceDofGlobally(index_t oldIdx, index_t newIdx);
    inline void replaceDofGlobally(index_t oldIdx, index_t newIdx, index_t comp);

    void mergeDofsGlobally(index_t dof1, index_t dof2);
    void mergeDofsGlobally(index_t dof1, index_t dof2, index_t comp);

// Data members
private:

    // m_dofs/m_patchDofs stores for each patch the mapping from local to global dofs.
    //
    // During setup, the entries have a different meaning:
    //   0        -- regular free dof
    //   negative -- an eliminated dof
    //   positive -- a coupling dof
    // For nonzero entries, the value is an id which identifies the eliminated/coupling
    // group of the dof. Dofs with the same id will get the same dof index in the
    // final numbering stage in finalize().

    // Representation of each component as a single vector plus
    // offsets for patch-local indices
    std::vector<std::vector<index_t> >  m_dofs;

    std::vector<size_t> m_offset;

    index_t m_shift;

    index_t m_bshift;

    std::vector<index_t> m_numFreeDofs;
    std::vector<index_t> m_numElimDofs;
    std::vector<index_t> m_numCpldDofs;

    // used during setup: running id for current eliminated dof
    // After finalize() is called m_curElimId takes positive value.
    index_t m_curElimId;

    std::vector<index_t> m_tagged;

}; // class gsDofMapper

inline std::ostream& operator<<( std::ostream& os, const gsDofMapper& b )
{
    return b.print( os );
}

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsDofMapper(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsDofMapper.hpp)
#endif