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