gsMultiBasis.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsGeometry.h>
#include <gsCore/gsBasis.h>
#include <gsCore/gsDofMapper.h>
#include <gsCore/gsBoxTopology.h>
#include <gsPde/gsBoundaryConditions.h>
#include <gsAssembler/gsAssemblerOptions.h>


namespace gismo
{

template<class T>
class gsMultiBasis : public gsFunctionSet<T>
{
public:
    typedef gsFunctionSet<T> Base;

    typedef memory::shared_ptr<gsMultiBasis> Ptr;
    typedef memory::unique_ptr<gsMultiBasis> uPtr;

    typedef std::vector<gsBasis<T> *> BasisContainer;

public:

    typedef typename BasisContainer::iterator iterator;
    typedef typename BasisContainer::const_iterator const_iterator;

    typedef          gsBasis<T> & reference;
    typedef          const gsBasis<T> & const_reference;

public:

    gsMultiBasis() { }

    explicit gsMultiBasis( const gsMultiPatch<T> & mpatch,
                           bool numeratorOnly = false);

    gsMultiBasis(BasisContainer& bases, const gsBoxTopology & topology)
        : m_topology( topology )
    {
        m_bases.swap(bases);// consumes the pointers
    }

    gsMultiBasis( const gsBasis<T> & geo, bool numeratorOnly = false);

    gsMultiBasis( BasisContainer & bases,
                  const std::vector<patchSide>& boundary,
                  const std::vector<boundaryInterface>& interfaces )
        : m_topology( bases[0]->dim(), bases.size(), boundary, interfaces )
    {
        m_bases.swap(bases);// consumes the pointers
    }

    ~gsMultiBasis();

    gsMultiBasis( const gsMultiBasis& other );

#if EIGEN_HAS_RVALUE_REFERENCES
    gsMultiBasis(gsMultiBasis&& other) : m_bases(give(other.m_bases)), m_topology(give(other.m_topology)) {}

    gsMultiBasis& operator= ( const gsMultiBasis& other );

    gsMultiBasis& operator= ( gsMultiBasis&& other )
    {
        freeAll(m_bases);
        m_bases = give(other.m_bases);
        m_topology = give(other.m_topology);
        return *this;
    }
#else
    gsMultiBasis& operator= ( gsMultiBasis other )
    {
        this->swap( other );
        return *this;
    }
#endif

    GISMO_CLONE_FUNCTION(gsMultiBasis)

public:

    const_iterator begin() const
    {
        return m_bases.begin();
    }

    const_iterator end() const
    {
        return m_bases.end();
    }

    iterator begin() {
        return m_bases.begin();
    }

    iterator end() {
        return m_bases.end();
    }

    void clear()
    {
        m_topology.clearAll();
        m_bases   .clear();
    }

    const gsBoxTopology & topology() const { return m_topology; }

public:

    const_reference at(size_t i) const
    {return *m_bases.at(i);}

    //reference at(size_t i)
    //{return *m_bases.at(i);}

    const_reference operator[](size_t i) const
    {return *m_bases[i];}

    reference operator[](size_t i)
    {return *m_bases[i];}

    // reference front()
    // {return *m_bases.front();}

    const_reference front() const
    {return *m_bases.front();}

    // reference back()
    // {return *m_bases.back();}

    const_reference back() const
    {return *m_bases.back();}

    void pop_back()
    {m_bases.pop_back();}

public:

    short_t domainDim () const {return m_bases.empty() ? 0 : m_bases.front()->domainDim();}

    short_t targetDim () const {return m_bases.empty() ? 0 : m_bases.front()->targetDim();}

    void swap(gsMultiBasis& other)
    {
        m_topology.swap( other.m_topology );

        m_bases.swap( other.m_bases );
    }

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

    short_t dim() const { return m_bases[0]->dim();}

    short_t degree(size_t i = 0, short_t comp = 0) const
    {
        GISMO_ASSERT( i < m_bases.size(),
                      "Invalid patch index "<<i<<" requested from gsMultiBasis" );
        return m_bases[i]->degree(comp);
    }

    short_t maxDegree(short_t k) const;

    short_t minDegree(short_t k) const;

    short_t maxCwiseDegree() const;

    short_t minCwiseDegree() const;

    int size(size_t i) const
    {
        GISMO_ASSERT( i < m_bases.size(),
                      "Invalid patch index "<<i<<" requested from gsMultiBasis" );
        return m_bases[i]->size();
    }

    index_t size() const
    {
        //GISMO_ERROR("call gsMultiBasis::nBases() instead.");
        return totalSize();
    }

    size_t totalSize() const
    {
        size_t sum = 0;
        for (size_t k = 0; k < m_bases.size(); ++k)
            sum += m_bases[k]->size();
        return sum;
    }

    size_t totalElements() const
    {
        size_t sum = 0;
        for (size_t k = 0; k < m_bases.size(); ++k)
            sum += m_bases[k]->numElements();
        return sum;
    }

    size_t nBases() const          { return m_bases.size(); }

    const gsBasis<T> & basis(const  size_t i ) const
    {
        GISMO_ASSERT( i < m_bases.size(),
                      "Invalid patch index"<<i<<" requested from gsMultiBasis" );
        return *m_bases[i];
    }

    const gsBasis<T> & piece(const index_t i) const
    {
        GISMO_ASSERT( static_cast<size_t>(i) < m_bases.size(),
                      "Invalid patch index"<<i<<" requested from gsMultiBasis" );
        return *m_bases[i];
    }

    index_t nPieces() const { return static_cast<index_t>(m_bases.size()); }

    gsBasis<T> & basis(const size_t i )
    {
        GISMO_ASSERT( i < m_bases.size(),
                      "Invalid patch index"<<i<<" requested from gsMultiBasis" );
        return *m_bases[i];
    }

    void addBasis( gsBasis<T> * g );

    void addBasis(typename gsBasis<T>::uPtr g);

    int findBasisIndex( gsBasis<T>* g ) const;

    void addInterface( gsBasis<T>* g1, boxSide s1,
                       gsBasis<T>* g2, boxSide s2 );

    void addPatchBoundary( gsBasis<T>* g, boxSide s )
    {
        const int p =findBasisIndex( g );
        m_topology.addBoundary( patchSide( p, s ) );
    }

    void uniformRefine(int numKnots = 1, int mul = 1, int dir = -1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
        {
            m_bases[k]->uniformRefine(numKnots,mul,dir);
        }
    }

    static void combineTransferMatrices(
        const std::vector< gsSparseMatrix<T, RowMajor> >& localTransferMatrices,
        const gsDofMapper& coarseMapper,
        const gsDofMapper& fineMapper,
        gsSparseMatrix<T, RowMajor>& transferMatrix
    );

    void uniformRefine_withTransfer(
        gsSparseMatrix<T, RowMajor>& transfer,
        const gsBoundaryConditions<T>& boundaryConditions,
        const gsOptionList& assemblerOptions,
        int numKnots = 1,
        int mul = 1,
        index_t unk = 0
        );

    void uniformRefineComponent(int comp, int numKnots = 1, int mul = 1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
        {
            m_bases[k]->component(comp).uniformRefine(numKnots,mul);
        }
    }

    // @brief Refine the boxes defined by "boxes"
    void refine(int k, gsMatrix<T> const & boxes)
    {
        m_bases[k]->refine(boxes);
    }

    // @brief Refine the boxes defined by "boxes"
    void unrefine(int k, gsMatrix<T> const & boxes)
    {
        m_bases[k]->unrefine(boxes);
    }

    void refineElements(int k, std::vector<index_t> const & boxes)
    {
        m_bases[k]->refineElements(boxes);
    }

    void unrefineElements(int k, std::vector<index_t> const & boxes)
    {
        m_bases[k]->unrefineElements(boxes);
    }

    void refine(size_t k, gsMatrix<T> const & boxes, int refExt)
    {
        GISMO_ASSERT( k < m_bases.size(),
                      "Invalid patch index "<<k<<" requested from gsMultiBasis" );
        m_bases[k]->refine( boxes, refExt);
    }
    void unrefine(size_t k, gsMatrix<T> const & boxes, int refExt)
    {
        GISMO_ASSERT( k < m_bases.size(),
                      "Invalid patch index "<<k<<" requested from gsMultiBasis" );
        m_bases[k]->unrefine( boxes, refExt);
    }

    void uniformCoarsen(int numKnots = 1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
        {
            m_bases[k]->uniformCoarsen(numKnots);
        }
    }

    void uniformCoarsen_withTransfer(
        gsSparseMatrix<T, RowMajor>& transfer,
        const gsBoundaryConditions<T>& boundaryConditions,
        const gsOptionList& assemblerOptions,
        int numKnots = 1,
        index_t unk = 0
        );

    typename gsBasis<T>::uPtr componentBasis(patchComponent p) const
    { return m_bases[p.patch()]->componentBasis(p); }

    typename gsBasis<T>::uPtr componentBasis_withIndices(
        patchComponent pc,
        const gsDofMapper& dm,
        gsMatrix<index_t>& indices,
        bool no_lower = true
    ) const;

    std::vector<typename gsBasis<T>::uPtr> componentBasis_withIndices(
        const std::vector<patchComponent>& pc,
        const gsDofMapper& dm,
        gsMatrix<index_t>& indices,
        bool no_lower = true
    ) const;

    void repairInterfaces( const std::vector< boundaryInterface > & bivec )
    {
        size_t kmax = 2*bivec.size();
        size_t k = 0;
        bool sthChanged = false;
        bool change = false;
        do
        {
            sthChanged = false;
            for( size_t i = 0; i < bivec.size(); i++ )
            {
                change = repairInterface( bivec[i] );
                sthChanged = sthChanged || change;
            }
            k++; // just to be sure this cannot go on infinitely
        }
        while( sthChanged && k <= kmax );
    }

    bool repairInterface2d( const boundaryInterface & bi );

    bool repairInterface( const boundaryInterface & bi );

    template<short_t d>
    bool repairInterfaceFindElements( const boundaryInterface & bi,
                                      std::vector<index_t> & refEltsFirst,
                                      std::vector<index_t> & refEltsSecond );

    void degreeElevate(short_t const i = 1, short_t const dir = -1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
            m_bases[k]->degreeElevate(i,dir);
    }

    void degreeIncrease(short_t const i = 1, short_t const dir = -1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
            m_bases[k]->degreeIncrease(i,dir);
    }

    void degreeDecrease(short_t const i = 1, short_t const dir = -1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
            m_bases[k]->degreeDecrease(i,dir);
    }

    void degreeReduce(short_t const i = 1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
            m_bases[k]->degreeReduce(i);
    }

    void setDegree(short_t const& i)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
            m_bases[k]->setDegree(i);
    }

    void reduceContinuity(int const i = 1)
    {
        for (size_t k = 0; k < m_bases.size(); ++k)
            m_bases[k]->reduceContinuity(i);
    }

    BasisContainer & patchBases()
    {
        return m_bases;
    }

    const BasisContainer & patchBases() const
    {
        return m_bases;
    }

    void setTopology(const gsBoxTopology & tpl)
    {
        m_topology = tpl;
    }


    void getMapper(bool conforming,
                   const gsBoundaryConditions<T> & bc,
                   int unk,
                   gsDofMapper & mapper,
                   bool finalize = true) const;

    void getMapper(bool conforming,
                   const gsBoundaryConditions<T> & bc,
                   gsDofMapper & mapper,
                   bool finalize = true) const
    { getMapper(conforming, bc, 0, mapper, finalize); }

    void getMapper(iFace::strategy is,
                   const gsBoundaryConditions<T> & bc,
                   gsDofMapper & mapper,
                   int unk,
                   bool finalize = true) const
    { getMapper(is==iFace::glue, bc, unk, mapper, finalize); }

    void getMapper(dirichlet::strategy ds,
                   iFace::strategy is,
                   const gsBoundaryConditions<T> & bc,
                   gsDofMapper & mapper,
                   int unk,
                   bool finalize = true) const
    {
        if ( ds == dirichlet::elimination )
            getMapper(is==iFace::glue, bc, unk, mapper, finalize);
        else
            getMapper(is==iFace::glue,        mapper, finalize);
    }

    gsDofMapper getMapper(dirichlet::strategy ds,
                          iFace::strategy is,
                          const gsBoundaryConditions<T> & bc,
                          int unk,
                          bool finalize = true) const
    {
        gsDofMapper mapper;
        if ( ds == dirichlet::elimination )
            getMapper(is==iFace::glue, bc, unk, mapper, finalize);
        else
            getMapper(is==iFace::glue,        mapper, finalize);
        return mapper;
    }


    // to remove
    void getMapper(bool conforming, gsDofMapper & mapper, bool finalize = true) const;

    void getMapper(iFace::strategy is, gsDofMapper & mapper, bool finalize = true) const
    { getMapper(is==iFace::glue, mapper, finalize); }


    //private: // to do

    void matchInterface(const boundaryInterface & bi,
                        gsDofMapper & mapper) const;

    void tileParameters();

private:

    BasisContainer m_bases;

    gsBoxTopology m_topology;

}; // class gsMultiBasis


template<class T>
std::ostream& operator<<( std::ostream& os, const gsMultiBasis<T>& b )
{
    return b.print( os );
}

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsMultiBasis(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


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