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++ )
            {
                if ( bivec[i].type() != interaction::contact)
                {
                    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