gsFunction.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsFunctionSet.h>

namespace gismo
{

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

    typedef memory::shared_ptr< gsFunction > Ptr;

    typedef memory::unique_ptr< gsFunction > uPtr;
    
    using Base::support;
    using Base::domainDim;
    using Base::targetDim;

    GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION(gsFunction, clone)

    virtual const gsFunction & piece(const index_t k) const
    {
        GISMO_ENSURE(0==k, "Single function of type "<< typeid(*this).name() <<" is defined on single subdomain, received: "<<k<<". Is piece(k) implemented?" );
        return *this; 
    }

    gsFuncCoordinate<T> coord(const index_t c) const;

    void active_into (const gsMatrix<T>  & u, gsMatrix<index_t> &result) const
    { result.setConstant(1,u.cols(),0); }
    
    /*
      Member functions with non-virtual implementations
      (override the _into versions in derived classes).
    */

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

    virtual void eval_component_into(const gsMatrix<T>& u, 
                                     const index_t comp, 
                                     gsMatrix<T>& result) const;

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

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

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

    gsMatrix<T> jacobian(const gsMatrix<T>& u) const;

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

    virtual void hessian_into(const gsMatrix<T>& u, gsMatrix<T>& result,
                              index_t coord = 0) const;

    virtual gsMatrix<T> hessian(const gsMatrix<T>& u, index_t coord = 0) const
    {
        gsMatrix<T> res;
        hessian_into(u,res,coord);
        return res;
    }

    virtual gsMatrix<T> laplacian( const gsMatrix<T>& u ) const;
  

    virtual T distanceL2(gsFunction<T> const &) const
    { GISMO_NO_IMPLEMENTATION }

    virtual void invertPoints(const gsMatrix<T> & points, gsMatrix<T> & result,
                              const T accuracy = 1e-6,
                              const bool useInitialPoint = false) const;

    virtual void invertPointGrid(gsGridIterator<T,0> & git,
                                 gsMatrix<T> & result, const T accuracy = 1e-6,
                                 const bool useInitialPoint = false) const;

    int newtonRaphson(const gsVector<T> & value,
                      gsVector<T> & arg,
                      bool withSupport = true, 
                      const T accuracy = 1e-6,
                      int max_loop = 100,
                      T damping_factor = 1) const;

    gsMatrix<T> argMin(const T accuracy = 1e-6,//index_t coord = 0
                       int max_loop = 100,
                       gsMatrix<T> init = gsMatrix<T>(),
                       T damping_factor = 1) const;

    void recoverPoints(gsMatrix<T> & xyz, gsMatrix<T> & uv, index_t k,
                           const T accuracy = 1e-6) const;

    void recoverPointGrid(gsGridIterator<T,0> & git,
                          gsMatrix<T> & xyz, gsMatrix<T> & uv,
                          index_t k, const T accuracy = 1e-6) const;

    virtual gsMatrix<T> parameterCenter() const
    { 
        // default impl. assumes convex support
        gsMatrix<T> S = this->support();
        return ( S.col(0) + S.col(1) ) * (T)(0.5);
    }

    gsMatrix<T> parameterCenter( const boxCorner& bc ) const;

    gsMatrix<T> parameterCenter( const boxSide& bs ) const;

    
    virtual std::ostream &print(std::ostream &os) const
    {
        os << "gsFunction.\n"; return os; 
    }

    virtual void computeMap(gsMapData<T> & InOut) const;

    index_t size() const { return 1;}

private:

    template<int mode, int _Dim=-1>
    int newtonRaphson_impl(
        const gsVector<T> & value,
        gsVector<T> & arg, bool withSupport = true,
        const T accuracy = 1e-6, int max_loop = 100,
        T damping_factor = 1, T scale = 1.0) const;

    gsVector<T> _argMinOnGrid(index_t numpts = 20) const;

    gsVector<T> _argMinNormOnGrid(index_t numpts = 20) const;
    
}; // class gsFunction


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

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsFunction(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11


} // namespace gismo

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