gsField.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsFunctionSet.h>
#include <gsCore/gsGeometry.h>
#include <gsCore/gsMultiPatch.h>
#include <gsCore/gsMultiBasis.h>
#include <gsUtils/gsPointGrid.h>

namespace gismo
{
template<class T>
class gsField
{

public:
    typedef memory::shared_ptr< gsField >  Ptr;// todo: remove
    typedef memory::unique_ptr< gsField >  uPtr;// todo: remove

    gsField(): m_patches(NULL) { }

    gsField( const gsFunctionSet<T> & mp,
             typename gsFunctionSet<T>::Ptr fs,
             const bool isparam)
    : m_patches(&mp), m_fields(fs), m_parametric(isparam)
    { }

    gsField( const gsGeometry<T> & sp, const gsFunctionSet<T> & pf, const bool isparam = false)
    : m_patches(&sp), m_fields(memory::make_shared_not_owned(&pf)), m_parametric(isparam)
    { }

    gsField( const gsGeometry<T> & sp, const gsGeometry<T> & pf)
    : m_patches(&sp), m_fields(memory::make_shared_not_owned(&pf)), m_parametric(true)
    { }

    gsField( const gsMultiPatch<T> & mp, const gsFunctionSet<T> & f, const bool isparam = false)
    : m_patches(&mp), m_fields(memory::make_shared_not_owned(&f)), m_parametric(isparam)
    { }

    gsField( const gsMultiPatch<T> & mp, const gsMultiPatch<T> & f)
    : m_patches(&mp), m_fields(memory::make_shared_not_owned(&f)), m_parametric(true)
    { }

public:

// TO DO:

// EVAL_physical_position
// Need to solve x = m_geometry(u) for u

    // Return a point in  the physical domain at parameter value u
    gsMatrix<T> point(const gsMatrix<T>& u, int i = 0) const
    {
        return m_patches->piece(i).eval(u);
    }

    // Return the value of the Field at parameter value u
    // TO DO: rename to evalParam()
    gsMatrix<T> value(const gsMatrix<T>& u, int i = 0)  const
    {
        return m_parametric
            ? m_fields->piece(i).eval(u)
            : m_fields->piece(i).eval( point(u, i) );
    }

    // Return the value of the Field at physical value u
    // TO DO: rename to evalPhys()
    gsMatrix<T> pvalue(const gsMatrix<T>& u, int i)  const
    {
        GISMO_ASSERT(!m_parametric, "Cannot compute physical value");
        return ( m_fields->piece(i).eval(u) );
    }

    T distanceL2(gsField<T> const & field, int numEvals= 1000) const;

    T distanceL2(gsFunctionSet<T> const & func,
                 bool isFunc_param = false,
                 int numEvals=1000) const;

    T distanceL2(gsFunctionSet<T> const & func,
                 gsMultiBasis<T> const & B,
                 bool isFunc_param = false,
                 int numEvals=1000) const;

    T distanceH1(gsFunctionSet<T> const & func,
                 bool isFunc_param = false,
                 int = 1000) const;

    T distanceH1(gsFunctionSet<T> const & func,
                 gsMultiBasis<T> const & B,
                 bool isFunc_param = false,
                 int = 1000) const;

    T distanceH2(gsFunctionSet<T> const & func,
                 bool isFunc_param = false) const;

    T distanceDG(gsFunctionSet<T> const & func,
                 bool isFunc_param = false,
                 int = 1000) const;

    std::ostream &print(std::ostream &os) const
    {
        os << ( m_parametric ? "Parameterized f" : "F")
           << "unction field.\n Defined on " << m_patches;
        return os;
    }

    short_t parDim() const { return m_patches->domainDim(); }

    short_t geoDim() const { return m_patches->targetDim(); }

    short_t dim() const { return m_fields->targetDim(); }

    GISMO_DEPRECATED index_t nPatches()  const { return m_patches->nPieces(); }

    int nPieces()  const { return m_patches->nPieces(); }

    const gsGeometry<T> & geometry() const
    {
        GISMO_ASSERT(dynamic_cast<const gsGeometry<T>*>(m_patches),
                     "No geometry in field. The domain is"<< *m_patches);
        return *static_cast<const gsGeometry<T>*>(m_patches);
    }

    const gsMultiPatch<T> & patches() const
    {
        GISMO_ASSERT(dynamic_cast<const gsMultiPatch<T>*>(m_patches),
                     "No patches in field. The field domain is "<< *m_patches);
        return *static_cast<const gsMultiPatch<T>*>(m_patches);
    }

    const gsFunctionSet<T> & fields() const { return *m_fields; }

    const gsGeometry<T> & patch(int i=0) const
    {
        GISMO_ASSERT( i<nPieces(),
                      "gsField: Invalid patch index.");
        GISMO_ASSERT(dynamic_cast<const gsGeometry<T>*>(&m_patches->piece(i)),
                     "No geometry in field. The domain is"<< m_patches->piece(i));
        return static_cast<const gsGeometry<T>&>(m_patches->piece(i));
    }

    const gsFunction<T> & function(int i=0) const
    {
        GISMO_ASSERT(dynamic_cast<const gsFunction<T>*>(&m_fields->piece(i)),
                     "No function in field. The domain is"<< m_patches->piece(i));
        return static_cast<const gsFunction<T>&>(m_fields->piece(i));
    }

    const gsGeometry<T> & igaFunction(int i=0) const
    {
        GISMO_ASSERT(i<m_fields->nPieces(),
                     "gsField: Invalid patch index.");
        GISMO_ASSERT(m_parametric,
                     "Cannot get an IGA function from non-parametric field.");
        GISMO_ASSERT(dynamic_cast<const gsGeometry<T>*>(&m_fields->piece(i)),
                     "Cannot return an igaFunction from a function of type: "<< m_fields->piece(i) );
        return static_cast<const gsGeometry<T> &>(m_fields->piece(i));
    }

    bool isParametric() const { return m_parametric; }

    bool isParametrized() const
    { return m_parametric && dynamic_cast<const gsGeometry<T>*>(&m_fields->piece(0));}

    const gsMatrix<T> & coefficientVector(int i=0) const
    {
        return igaFunction(i).coefs();
    }

// Data members
private:

    const gsFunctionSet<T> * m_patches;

    // If there are many patches, one field per patch

    typename gsFunctionSet<T>::Ptr m_fields;

    bool m_parametric;// True iff this is an Isogeometric field, living on parameter domain

}; // class gsField


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


} // namespace gismo

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