gsFunctionSet.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
//#include <gsCore/gsForwardDeclarations.h> // included by gsLinearAlgebra.h

// The following macros are a design pattern to solve the problem with
// unique::memory pointers as return value of virtual functions in base/derived
// classes. It is expected that a class where this macros are used is derived
// from gsFunctionSet or its derivatives. It assumes that a concrete
// implementation has the suffix "_impl". From outside that class, some can
// call that function by its name and get back a pointer inside a
// memory::unique_ptr (aka. uPtr) of the correct type. If casts are needed
// afterward, use memory::convert_ptr<toType>(from).

// NOTE: One has to add for the DOXYGEN documentation the signature by hand
// and "comment out" the _impl method in the hpp file. This can be done with
// the preprocessor definition __DOXYGEN__. Commenting out is done with @cond
// till @endcond as a doxygen comment.

// Example:

// #ifdef __DOXYGEN__
// /// @briefSome Method
// /// @param a
// /// @param b
// /// @return gsBasis<T>::uPtr
// gsBasis<T>::uPtr someMethod(int a, const int b)
// #endif
//
// GISMO_UPTR_FUNCTION_PURE(gsBasis<T>, someMethod, int, const int)

// /// @cond
// gsBasis<T>* someClass<T>::someMethod_impl(int a, const int b) { .... }
// /// @endcond

// Helper macros for counting arguments, works till highest number in PP_RSEQ_N
// Call with PP_NARG(__VA_ARGS__)
// For upgrade, add values to PP_ARG_N, PP_RSEQ_N and PP_COMMASEQ_N
#define PP_EXPAND(x) x

#define PP_ARG_N(_1,_2,_3,N,...) N
#define PP_RSEQ_N() 3,2,1,0
#define PP_NARG_(...) PP_EXPAND(PP_ARG_N(__VA_ARGS__))
#define PP_COMMASEQ_N()  1,1,0,0
#define PP_COMMA(...) ,
#define PP_HASCOMMA(...) PP_NARG_(__VA_ARGS__,PP_COMMASEQ_N())
#define PP_NARG(...) PP_NARG_HELPER1(PP_HASCOMMA(__VA_ARGS__),PP_HASCOMMA(PP_COMMA __VA_ARGS__ ()),PP_NARG_(__VA_ARGS__, PP_RSEQ_N()))
#define PP_NARG_HELPER1(a,b,N) PP_NARG_HELPER2(a, b, N)
#define PP_NARG_HELPER2(a,b,N) PP_NARG_HELPER3_ ## a ## b(N)
#define PP_NARG_HELPER3_01(N) 0
#define PP_NARG_HELPER3_00(N) 1
#define PP_NARG_HELPER3_11(N) N

// Declaration prototypes. Followed by: { ";" , "= 0;" , "{ ... }" }
#define __DECn(n, type, name, ...)  __DEC ## n(type, name, __VA_ARGS__)
#define __DEC0(type, name, void)    private: virtual type * name##_impl() const
#define __DEC1(type, name, t1)      private: virtual type * name##_impl(t1 n1) const
#define __DEC2(type, name, t1, t2)  private: virtual type * name##_impl(t1 n1, t2 n2) const

// Definition prototypes
#define __DEFn(n, type, name, ...)  __DEF ## n(type, name, __VA_ARGS__)
#define __DEF0(type, name, void)    public:  inline memory::unique_ptr< type > name() const { return memory::unique_ptr< type >(name##_impl()); }
#define __DEF1(type, name, t1)      public:  inline memory::unique_ptr< type > name(t1 n1) const { return memory::unique_ptr< type >(name##_impl(n1)); }
#define __DEF2(type, name, t1, t2)  public:  inline memory::unique_ptr< type > name(t1 n1, t2 n2) const { return memory::unique_ptr< type >(name##_impl(n1, n2)); }

// Declaration of virtual function
// 1st: return type
// 2nd: function name
// 3rd: types of parameter arguments
// Definition of real the implementation in the form "type * name_impl(...)"
// should be done in the corresponding .hpp file.
#define GISMO_UPTR_FUNCTION_DEC(type, name, ...) \
        GISMO_UPTR_FUNCTION_DEC_(PP_NARG(__VA_ARGS__), type, name, __VA_ARGS__)
#define GISMO_UPTR_FUNCTION_DEC_(n, type, name, ...) \
        __DECn(n, type, name, __VA_ARGS__); \
        __DEFn(n, type, name, __VA_ARGS__)

// Declaration and start of definition of virtual function
// 1st: return type
// 2nd: function name
// 3rd: types of parameter arguments
// must be finished with a block of { return type * your implementation }
// don't forget that your in "private:" afterward
#define GISMO_UPTR_FUNCTION_DEF(type, name, ...) \
        GISMO_UPTR_FUNCTION_DEF_(PP_NARG(__VA_ARGS__), type, name, __VA_ARGS__)
#define GISMO_UPTR_FUNCTION_DEF_(n, type, name, ...) \
        __DEFn(n, type, name, __VA_ARGS__) \
        __DECn(n, type, name, __VA_ARGS__)

// Declaration of pure virtual function
// 1st: return type
// 2nd: function name
// 3rd: types of parameter arguments
#define GISMO_UPTR_FUNCTION_PURE(type, name, ...) \
        GISMO_UPTR_FUNCTION_PURE_(PP_NARG(__VA_ARGS__), type, name, __VA_ARGS__)
#define GISMO_UPTR_FUNCTION_PURE_(n, type, name, ...) \
        __DECn(n, type, name, __VA_ARGS__) = 0; \
        __DEFn(n, type, name, __VA_ARGS__)

// Declaration and definition with GISMO_NO_IMPLEMENTATION
// 1st: return type
// 2nd: function name
// 3rd: types of parameter arguments
#define GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION(type, name, ...) \
        GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION_(PP_NARG(__VA_ARGS__), type, name, __VA_ARGS__)
#define GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION_(n, type, name, ...) \
        __DECn(n, type, name, __VA_ARGS__) { GISMO_NO_IMPLEMENTATION } \
        __DEFn(n, type, name, __VA_ARGS__)

// Declaration, definition and implementation of clone function
// 1st: return type
#define GISMO_CLONE_FUNCTION(type) \
        __DEC0(type, clone, void) { return new type(*this); } \
        __DEF0(type, clone, void)

// Declaration, definition and implementation of clone function which overrides
// 1st: return type
#define GISMO_OVERRIDE_CLONE_FUNCTION(type) \
        __DEC0(type, clone, void) override { return new type(*this); } \
        __DEF0(type, clone, void)


namespace gismo {

template <typename T>
class gsFunctionSet
{
public:

    typedef memory::shared_ptr< gsFunctionSet > Ptr;

    typedef memory::unique_ptr< gsFunctionSet > uPtr;

public:

    virtual ~gsFunctionSet();

#ifdef __DOXYGEN__
    uPtr clone();
#endif
    GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION(gsFunctionSet, clone)

    
    virtual const gsFunctionSet & piece(const index_t) const {return *this;}

    const gsFunction<T> & function(const index_t k) const;

    const gsBasis<T> & basis(const index_t k) const;

public:

    /*
     @brief Returns (a bounding box for) the support of the function(s).
    
     Returns either a zero-sized matrix or a dx2 matrix, containing
     the two diagonally extreme corners of a hypercube.

     If the returned matrix is empty, it is understood that the domain
     of definition is the whole of \f$R^{domainDim}\f$
    */
    virtual gsMatrix<T> support() const;
    virtual gsMatrix<T> supportOf(const index_t & i) const;

    virtual void active_into (const gsMatrix<T>  & u, gsMatrix<index_t> &result) const;
    
public:
    virtual void eval_into      (const gsMatrix<T>  & u, gsMatrix<T> &result) const;

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

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

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

    std::vector<gsMatrix<T> > evalAllDers(const gsMatrix<T> & u, int n, bool sameElement = false) const;

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

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

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


    gsMatrix<index_t> active(const gsMatrix<T> & u) const
    {
        gsMatrix<index_t> rvo;
        this->active_into(u, rvo);
        return rvo;
    }

    /*
       @brief Divergence (of vector-valued functions).

       This makes sense only for vector valued functions whose target dimension is a multiple of the domain dimension
       \f[\left[
       \begin{array}{ccccc}
       \text{div} f_1(p_1) & \text{div} f_1(p_2) & \ldots & \text{div} f_1(p_N)\\
       \text{div} f_2(p_1) & \text{div} f_2(p_2) & \ldots & \text{div} f_2(p_N)\\
       \vdots       & \vdots       &        & \vdots\\
       \text{div} f_S(p_1) & \text{div} f_S(p_2) & \ldots & \text{div} f_S(p_N)
       \end{array}
       \right]
       \f]
       If the target dimension is an integer multiple of the domain dimension, \f$\text{div} f(p)\f$ contains
       <em>(target dim)/(domain dim)</em> rows, each row is the divergence of a <em>(domain dim)</em>-tuple of entries of \f$f\f$.
       Equivalently, \f$f\f$ is interpreted as a matrix valued function in column maior standard and
       the divergence corresponding vector valued divergence is returned.

       If the target dimension is not an integer multiple of the domain dimension calling this function results in an error.

       @param u
       @param result
    virtual void div_into       (const gsMatrix<T>  & u, gsMatrix<T> &result) const;
    */


    /*
       @brief Curl (of vector-valued functions).

       This makes sense only for vector-valued functions whose target dimension is a multiple of the domain dimension.
       \f[\left[
       \begin{array}{ccccc}
       \text{curl} f_1(p_1) & \text{curl} f_1(p_2) & \ldots & \text{curl} f_1(p_N)\\
       \text{curl} f_2(p_1) & \text{curl} f_2(p_2) & \ldots & \text{curl} f_2(p_N)\\
       \vdots       & \vdots       &        & \vdots\\
       \text{curl} f_S(p_1) & \text{curl} f_S(p_2) & \ldots & \text{curl} f_S(p_N)
       \end{array}
       \right]
       \f]
       If the target dimension is an integer multiple of the domain dimension, \f$\text{curl} f(p)\f$ contains
       <em>(target dim)/(domain dim)</em> rows, each row is the rotor of a <em>(domain dim)</em>-tuple of entries of \f$f\f$.
       Equivalently, \f$ f\f$ is interpreted as a matrix valued function in column maior standard and
       the rotor corresponding vector valued rotor is returned.

       If the target dimension is not an integer multiple of the domain dimension calling this function results in an error.
       @param u
       @param result
    virtual void curl_into      (const gsMatrix<T> &  u, gsMatrix<T> &result) const;
    */


    /*
       @brief Laplacian (of vector-valued functions).

       \f[\left[
       \begin{array}{ccccc}
       \Delta f_1(p_1) & \Delta f_1(p_2) & \ldots & \Delta f_1(p_N)\\
       \Delta f_2(p_1) & \Delta f_2(p_2) & \ldots & \Delta f_2(p_N)\\
       \vdots       & \vdots       &        & \vdots\\
       \Delta f_S(p_1) & \Delta f_S(p_2) & \ldots & \Delta f_S(p_N)
       \end{array}
       \right]
       \f]
       If the target dimension is different from 1, then each \f$\Delta f_i(p_j)\f$ is a block of target
       dimension rows containing the component wise Laplacian.
       @param u
       @param result
    virtual void laplacian_into (const gsMatrix<T>  &  u, gsMatrix<T> &result) const;
    */



    virtual void compute(const gsMatrix<T> & in, gsFuncData<T> & out) const;

public:
    virtual short_t domainDim () const = 0;

    virtual short_t targetDim () const {return 1;}

    /*
      @brief Dimension of domain and target
      @return the pair of integers domainDim() and targetDim()
    */
    std::pair<short_t, short_t> dimensions() const {return std::make_pair(domainDim(),targetDim());}
    
    virtual index_t size() const //= 0;
    {GISMO_NO_IMPLEMENTATION}

    virtual index_t nPieces() const {return 1;}

    virtual std::ostream &print(std::ostream &os) const// = 0;
    {
        os << "gsFunctionSet\n";
        return os; 
    }

};

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

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsFunctionSet(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


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