gsMappedBasis.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsBoxTopology.h>
#include <gsCore/gsMultiPatch.h>

#include <gsMSplines/gsMappedSingleBasis.h>
#include <gsMSplines/gsWeightMapper.h>
#include <gsCore/gsDomainIterator.h>
#include <gsCore/gsPiecewiseFunction.h>
#include <gsCore/gsBasisFun.h>

namespace gismo
{

// forward declarations of the mapper classes
template<short_t d,class T> class gsMappedSingleBasis;

template<short_t d,class T>
class gsMappedBasis  : public gsFunctionSet<T>
{
    friend class gsMappedSingleBasis<d,T>;

private:
    typedef std::vector<std::pair<int,int> >::iterator step_iter;
    typedef gsBasis<T> BasisType;
    typedef typename std::vector<BasisType *>::const_iterator ConstBasisIter;
    typedef typename std::vector<BasisType *>::iterator BasisIter;
    typedef typename std::vector<gsMatrix<T> *>::const_iterator ConstMatrixPtrIter;
    typedef typename gsSparseMatrix<T,0,index_t>::InnerIterator InIterMat;
    typedef typename std::vector<index_t> IndexContainer;
    typedef typename std::vector<index_t>::const_iterator ConstIndexIter;
    typedef typename std::vector<T> WeightContainer;
    typedef typename std::vector<T>::const_iterator ConstWeightIter;
    typedef gsEigen::PermutationMatrix<Dynamic,Dynamic,index_t> gsPermutationMatrix;

    typedef memory::unique_ptr< gsDomainIterator<T> > domainIter;

public:
    typedef memory::shared_ptr< gsMappedBasis > Ptr;

    typedef memory::unique_ptr< gsMappedBasis > uPtr;

    gsMappedBasis() : m_mapper(nullptr)
    { }

    gsMappedBasis(gsMultiBasis<T> const & mb, const gsSparseMatrix<T> & m)
    : m_mapper(nullptr)
    {
        init(mb,m);
    }

    gsMappedBasis( const gsMappedBasis& other );

    virtual ~gsMappedBasis();

    void init(gsMultiBasis<T> const & mb, const gsSparseMatrix<T> & m)
    {
        GISMO_ASSERT(mb.domainDim()==d, "Error in dimensions");
        m_topol  = mb.topology();

        delete m_mapper;
        m_mapper = new gsWeightMapper<T>(m);

        freeAll(m_bases);
        m_bases.reserve(mb.nBases());
        cloneAll(mb.patchBases(),m_bases);
        m_sb.clear();
        m_sb.reserve(m_bases.size());
        for ( size_t q=0; q!=m_bases.size(); ++q )
            m_sb.push_back( gsMappedSingleBasis<d,T>(this,q) );

        m_mapper->optimize(gsWeightMapper<T>::optSourceToTarget);
    }

    index_t nPieces() const {return m_topol.nBoxes();}

public:
    // getters for the private fields

    BasisType const & getBase(index_t i) const
    { return *m_bases[i]; }

    BasisType & getBase(index_t i)
    { return *m_bases[i]; }

    const std::vector<BasisType*> getBases() const;

    std::vector<BasisType*> getBasesCopy() const
    {
        std::vector<BasisType*> bases;
        cloneAll(m_bases,bases);
        return bases;
    }

    gsWeightMapper<T> const & getMapper() const
    { return *m_mapper; }

    gsWeightMapper<T> * getMapPointer() const
    { return m_mapper; }

    gsBoxTopology const & getTopol() const
    { return m_topol; }

public:
    // general functions for interacting with this class

    GISMO_CLONE_FUNCTION(gsMappedBasis)

    // Prints the object to the stream
    //std::ostream & print(std::ostream & os) const { }

    short_t domainDim() const
    {
        GISMO_ASSERT(m_bases.size()>0,"there should be at least one basis provided.");
        return m_bases[0]->domainDim();
    }

    index_t size(const index_t index) const;
    index_t size() const { return m_mapper->getNrOfTargets(); }

    index_t globalSize() const { return m_mapper->getNrOfTargets(); }

    index_t localSize() const { return m_mapper->getNrOfSources(); }

    index_t localSize(const index_t index) const
    { return m_bases[index]->size(); }

    short_t maxDegree() const;

    short_t degree(const index_t patch,const short_t i) const
    { return m_bases[patch]->degree(i); }

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

private:
    void addLocalIndicesOfPatchSide(const patchSide& ps,index_t offset,std::vector<index_t>& locals) const;

public:
    void reorderDofs(const gsPermutationMatrix& permMatrix)
    {
        (*m_mapper)*=permMatrix;
        m_mapper->optimize(gsWeightMapper<T>::optSourceToTarget);
    }

    void reorderDofs_withCoef(const gsPermutationMatrix& permMatrix,gsMatrix<T>& coefs)
    {
        reorderDofs(permMatrix);
        coefs*=permMatrix;
    }

    void boundary(std::vector<index_t> & indices,index_t offset = 0) const;

    void innerBoundaries(std::vector<index_t> & indices,index_t offset = 0) const;

    gsMappedSingleBasis<d,T> & getMappedSingleBasis(const index_t i) const
    {
        // if (m_sb.empty())
        // {
        //     m_sb.reserve(m_bases.size());
        //     for (size_t q = 0; q!=m_bases.size(); ++q)
        //         m_sb.push_back( gsMappedSingleBasis<d,T>(this,q) );
        // }
        return m_sb[i];
    }

    const gsMappedSingleBasis<d,T> & piece(const index_t k) const { return getMappedSingleBasis(k); }
    //const gsFunctionSet & piece(const index_t k) const { return m_sb[k]; }

    domainIter makeDomainIterator(const index_t patch,const boxSide & s) const
    { return m_bases[patch]->makeDomainIterator(s); }

    gsGeometry<T>* exportPatch(const index_t i,gsMatrix<T> const & localCoef) const;

    gsMultiPatch<T> exportToPatches(gsMatrix<T> const & localCoef) const;

private:
    // Avoid warnings for hidden overloads w.r.t gsFunctionSet
    void active_into(const gsMatrix<T> &,gsMatrix<index_t>&) const
    { GISMO_NO_IMPLEMENTATION; }
    void eval_into(const gsMatrix<T> &,gsMatrix<T>&) const
    { GISMO_NO_IMPLEMENTATION; }
    void deriv_into(const gsMatrix<T> &,gsMatrix<T>&) const
    { GISMO_NO_IMPLEMENTATION; }
    void deriv2_into(const gsMatrix<T> &,gsMatrix<T>&) const
    { GISMO_NO_IMPLEMENTATION; }
    void evalAllDers_into(const gsMatrix<T> &, int,
                          std::vector<gsMatrix<T> >&, bool) const
    { GISMO_NO_IMPLEMENTATION; }

public:
    // functions for evaluating and derivatives

    void active_into(const index_t patch, const gsMatrix<T> & u,
                     gsMatrix<index_t>& result) const; //global BF active on patch at point

    virtual void numActive_into(const index_t patch,const gsMatrix<T> & u, gsVector<index_t>& result) const
    {
        GISMO_UNUSED(patch); GISMO_UNUSED(u); GISMO_UNUSED(result);
        GISMO_NO_IMPLEMENTATION
    }


    void eval_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T>& result ) const;
    void deriv_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T>& result ) const;
    void deriv2_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T>& result ) const;

    gsPiecewiseFunction<T> basisFunction(index_t global_BF)
    {
        const size_t np = nPatches();
        gsPiecewiseFunction<T> bf(np);
        for( size_t i = 0; i!=np; ++i)
            bf.addPiece( getMappedSingleBasis(i).function(global_BF) );
        return bf;
    }

    void evalSingle_into(const index_t patch, const index_t global_BF, const gsMatrix<T> & u, gsMatrix<T>& result ) const;
    void derivSingle_into(const index_t patch, const index_t global_BF, const gsMatrix<T> & u, gsMatrix<T>& result ) const;
    void deriv2Single_into(const index_t patch, const index_t global_BF, const gsMatrix<T> & u, gsMatrix<T>& result ) const;

    gsMatrix<T> evalSingle(const index_t patch, const index_t global_BF, const gsMatrix<T> & u) const
    {
        gsMatrix<T> result;
        this->evalSingle_into(patch, global_BF, u, result);
        return result;
    }

    gsMatrix<T> derivSingle(const index_t patch, const index_t global_BF, const gsMatrix<T> & u) const
    {
        gsMatrix<T> result;
        this->derivSingle_into(patch, global_BF, u, result);
        return result;
    }

    gsMatrix<T> deriv2Single(const index_t patch, const index_t global_BF, const gsMatrix<T> & u) const
    {
        gsMatrix<T> result;
        this->deriv2Single_into(patch, global_BF, u, result);
        return result;
    }

    void evalAllDers_into(const index_t patch, const gsMatrix<T> & u,
                          const index_t n, std::vector<gsMatrix<T> >& result,
                          bool sameElement = false) const;

    void evalAllDersSingle_into(const index_t patch,const index_t global_BF, const gsMatrix<T> & u,const index_t n,gsMatrix<T> & result ) const;


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

public:
    // functions for converting global to local coefs and back

    void global_coef_to_local_coef(gsMatrix<T> const & globalCoefs,gsMatrix<T> & localCoefs) const
    { m_mapper->mapToSourceCoefs(globalCoefs,localCoefs); }

    void local_coef_to_global_coef(gsMatrix<T> const & localCoefs,gsMatrix<T> & globalCoefs) const
    {  m_mapper->mapToTargetCoefs(localCoefs,globalCoefs); }

    index_t getPatch(index_t const localIndex)
    {
        return _getPatch(localIndex);
    }
    index_t getPatchIndex(index_t const localIndex)
    {
        return _getPatchIndex(localIndex);
    }

protected:
    // functions for working with Indexes

    index_t _getPatch(index_t const localIndex) const;

    index_t _getPatchIndex(const index_t localIndex) const;

    index_t _getLocalIndex(index_t const patch,index_t const patchIndex) const
    { return _getFirstLocalIndex(patch)+patchIndex; }

    index_t _getFirstLocalIndex(index_t const patch) const;

    index_t _getLastLocalIndex(index_t const patch) const
    { return _getFirstLocalIndex(patch)+m_bases[patch]->size()-1; }

    // Data members
protected:
    gsBoxTopology m_topol;

    std::vector<BasisType *> m_bases;

    gsWeightMapper<T> * m_mapper;
    // gsSparseMatrix<T> r:C, c:B

    mutable std::vector<gsMappedSingleBasis<d,T> > m_sb;

    // Make gsMultiBasis a member instead of m_bases and m_topol?
};

#ifdef GISMO_WITH_PYBIND11

  // void pybind11_init_gsMappedBasis1(pybind11::module &m);
  void pybind11_init_gsMappedBasis2(pybind11::module &m);
  // void pybind11_init_gsMappedBasis3(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

}

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