gsFuncCoordinate.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsFunction.h>
#include <gsCore/gsFunctionSet.h>

namespace gismo
{

template<class T>
class gsFuncCoordinate : public gsFunction<T>
{
public:
    gsFuncCoordinate() : m_index(1) { }

    typedef memory::shared_ptr< gsFuncCoordinate > Ptr;

    typedef memory::unique_ptr< gsFuncCoordinate > uPtr;

    gsFuncCoordinate(const gsFunctionSet<T> & function, unsigned const i);

    gsFuncCoordinate(const gsFunctionSet<T> & function, gsVector<index_t> ind)
    : m_function(&function), m_index(give(ind)) { }

    ~gsFuncCoordinate() { } //destructor

    GISMO_CLONE_FUNCTION(gsFuncCoordinate)

public:

    short_t domainDim () const {return m_function->domainDim();}

    short_t targetDim () const {return (short_t)m_index.size();} // It is a coordinate of a vector function

    gsMatrix<T> support() const;

    void eval_into (const gsMatrix<T>& u, gsMatrix<T>& result ) const;

    void deriv_into(const gsMatrix<T>& u, gsMatrix<T>& result ) const;

    void deriv2_into(const gsMatrix<T>& u, gsMatrix<T>& result ) const;

    void setCoordinate( unsigned const & i );

    void first();

    bool next();

    bool valid();

    index_t index() const { return m_index.value(); }

    // // temporary hack
    virtual const gsFuncCoordinate & piece(const index_t k) const
    {
        auto it = m_pieces.find(k);
        if (m_pieces.end()!=it) return it->second;
        std::pair<typename std::map<index_t,gsFuncCoordinate>::iterator,bool> ins;            
#       pragma omp critical (m_pieces_touch)
        ins = m_pieces.emplace(k, gsFuncCoordinate(m_function->piece(k),m_index));
        return ins.first->second;
    }

// Data members
private:
    const gsFunctionSet<T> * m_function;
    gsVector<index_t> m_index;

    mutable std::map<index_t,gsFuncCoordinate> m_pieces;
}; // class gsFuncCoordinate


template<class T>
gsFuncCoordinate<T>::gsFuncCoordinate(const gsFunctionSet<T> & function, const unsigned i )
: m_function(&function), m_index(1)
{
    GISMO_ASSERT( i<unsigned(m_function->targetDim()),"Invalid coordinate" );
    m_index[0] = i;
}

template<class T> void
gsFuncCoordinate<T>::setCoordinate( unsigned const & i )
{
    GISMO_ASSERT( i<unsigned(m_function->targetDim()),"Invalid coordinate" );
    m_index[0] = i;
}

template<class T> void
gsFuncCoordinate<T>::first()
{
    m_index.resize(1);
    m_index[0] = 0;
}

template<class T> bool
gsFuncCoordinate<T>::next()
{
    return ( ++m_index[0]  < static_cast<index_t>(m_function->targetDim()) );
}

template<class T> bool
gsFuncCoordinate<T>::valid()
{
    return ( m_index[0]  < static_cast<index_t>(m_function->targetDim()) );
}

template<class T> gsMatrix<T>
gsFuncCoordinate<T>::support()  const
{
    return m_function->support();
}

template<class T> void
gsFuncCoordinate<T>::eval_into(const gsMatrix<T>& u, gsMatrix<T>& result )  const
{
    gsMatrix<T> tmp;
    m_function->eval_into(u, tmp);
    result = tmp(m_index, gsEigen::all);
}

template<class T> void
gsFuncCoordinate<T>::deriv_into(const gsMatrix<T>& u, gsMatrix<T>& result )  const
{
    gsMatrix<T> tmp;
    m_function->deriv_into(u, tmp);
    const index_t stride = domainDim();
    result.resize(m_index.size()*stride, u.cols() );
    for (index_t i = 0; i != m_index.size(); ++i)
        result.middleRows(i*stride,stride) = tmp.middleRows(stride*m_index[i], stride);
}

template<class T> void
gsFuncCoordinate<T>::deriv2_into(const gsMatrix<T>& u, gsMatrix<T>& result )  const
{
    gsMatrix<T> tmp;
    m_function->deriv2_into(u, tmp);
    const index_t dim = domainDim();
    const index_t stride = dim*(dim+1)/2;
    result.resize(m_index.size()*stride, u.cols() );
    for (index_t i = 0; i != m_index.size(); ++i)
        result.middleRows(i*stride,stride) = tmp.middleRows(stride*m_index[i], stride);
}


} // namespace gismo