gsFieldCreator.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsField.h>

namespace gismo
{


template<class T>
class gsAbsError : public gsFunction<T> 
{
public:
    typedef memory::shared_ptr< gsAbsError > Ptr;

    typedef memory::unique_ptr< gsAbsError > uPtr;

    gsAbsError( gsFunction<T> const & f1, 
                gsGeometry<T> const & geo, 
                gsFunction<T> const & f2,
                bool _f2param = false)
    : m_geo(geo), m_f1(f1), m_f2(f2), f2param(_f2param)
    { 
            
    }

    GISMO_CLONE_FUNCTION(gsAbsError)

    void eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
    {
        gsMatrix<T> tmp;
        m_geo.eval_into(u, tmp);

        if(!f2param)
            result.noalias() = ( m_f1.eval(u) - m_f2.eval(tmp) ).cwiseAbs();
        else
            result.noalias() = ( m_f1.eval(u) - m_f2.eval(u) ).cwiseAbs(); // f2 parametric function

    }

    short_t domainDim() const { return m_f1.domainDim(); }
    short_t targetDim() const { return m_f1.targetDim(); }

    std::ostream &print(std::ostream &os) const
    { os << "Absolute error.\n"; return os; };
private:
    const gsGeometry<T> & m_geo;

    const gsFunction<T> & m_f1;
    const gsFunction<T> & m_f2;

    bool f2param;

private:
//    gsAbsError() { }
};


template<class T>
class gsGradientField : public gsFunction<T> 
{
public:
    typedef memory::shared_ptr< gsGradientField > Ptr;

    typedef memory::unique_ptr< gsGradientField > uPtr;

    gsGradientField( gsGeometry<T> const & geo, gsFunction<T> const & f)
    : m_geo(geo), m_f(f)
    { }

    GISMO_CLONE_FUNCTION(gsGradientField)

    void eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
    {
        gsMatrix<T> tmp;
        m_geo.eval_into(u, tmp);
        m_f.deriv_into(tmp, result);
    }

    short_t domainDim() const { return m_geo.parDim(); }
    short_t targetDim() const { return m_geo.geoDim(); }

    std::ostream &print(std::ostream &os) const
    { os << "Gradient field of "; return os; };
private:
    const gsGeometry<T> & m_geo;
    const gsFunction<T> & m_f  ;
};

template<class T>
class gsJacDetField : public gsFunction<T> 
{
public:
    typedef memory::shared_ptr< gsJacDetField > Ptr;

    typedef memory::unique_ptr< gsJacDetField > uPtr;

    gsJacDetField( gsGeometry<T> const & geo)
    : m_geo(geo), m_dim(m_geo.parDim())
    { 
        GISMO_ENSURE(m_dim == m_geo.geoDim(), "Not extended to surface case yet." );
    }

    GISMO_CLONE_FUNCTION(gsJacDetField)

    void eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
    {
        result.resize(1, u.cols() );
        gsMatrix<T> tmp;
        for (index_t k=0; k!= u.cols(); ++k)
            result(0,k) = m_geo.deriv(u.col(k)).reshape(m_dim,m_dim).determinant();
    }

    short_t domainDim() const { return m_geo.parDim(); }
    short_t targetDim() const { return 1; }

    std::ostream &print(std::ostream &os) const
    { os << "Jacobian determinant field of "; return os; };
private:
    const gsGeometry<T> & m_geo;
    const int m_dim;
};

template<class T>
class gsNormalField : public gsFunction<T> 
{
public:
    typedef memory::shared_ptr< gsNormalField > Ptr;

    typedef memory::unique_ptr< gsNormalField > uPtr;

    gsNormalField( gsGeometry<T> const & geo) : m_geo(geo)
    { 
            
    }

    GISMO_CLONE_FUNCTION(gsNormalField)

    void eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
    {
        const short_t ParDim = m_geo.parDim();
        result.resize(ParDim+1, u.cols()) ;

        gsMatrix<T> Jk;

        for( index_t j=0; j < u.cols(); ++j )
        {
            Jk = m_geo.jacobian( u.col(j) );
            T alt_sgn(1.0);
            gsMatrix<T> mm(ParDim,ParDim);
            for (int i = 0; i <= ParDim; ++i) // for all components of the normal vector
            {
                Jk.rowMinor(i, mm);
                result(i,j) = alt_sgn * mm.determinant();
                alt_sgn = -alt_sgn;
            }
        }
    }

    short_t domainDim() const { return m_geo.parDim(); }
    short_t targetDim() const { return m_geo.geoDim();}

    std::ostream &print(std::ostream &os) const
    { os << "NormalField"; return os; };
private:
    const gsGeometry<T> & m_geo;

private:
//    gsNormalField() { }
};

template<class T>
class gsParamField : public gsFunction<T> 
{
public:
    typedef memory::shared_ptr< gsParamField > Ptr;

    typedef memory::unique_ptr< gsParamField > uPtr;

    gsParamField(gsGeometry<T> const & geo)
    : m_geo(geo)
    { 
            
    }

    GISMO_CLONE_FUNCTION(gsParamField)

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

    short_t domainDim() const { return m_geo.domainDim(); }
    short_t targetDim() const { return m_geo.domainDim(); }

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

private:
    const gsGeometry<T> & m_geo;
};


template<class T>
class gsPatchIdField : public gsFunction<T>
{
public:
    typedef memory::shared_ptr< gsPatchIdField > Ptr;

    typedef memory::unique_ptr< gsPatchIdField > uPtr;

    explicit gsPatchIdField(gsGeometry<T> const & geo_)
    : geo(geo_), m_supp(geo.support())
    { }

    GISMO_CLONE_FUNCTION(gsPatchIdField)

    void eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
    { 
        result.setZero(2, u.cols() );
        result.row(0).setConstant(geo.id());

        const int d = geo.parDim();
        for (boxSide c=boxSide::getFirst(d); c<boxSide::getEnd(d); ++c)
        {
            const index_t dir = c.direction();
            const T par = m_supp(dir, c.parameter() );

            for (index_t v = 0; v != u.cols(); ++v) // for all columns of u
            {
                if ( math::abs( u(dir,v) - par ) < 1e-2 )
                    result(1,v) = static_cast<T>(c);
            }
        }
    }

    short_t domainDim() const { return geo.domainDim(); }
    short_t targetDim() const { return 2; }

    std::ostream &print(std::ostream &os) const
    { os << "Boundary side indicator field.\n"; return os; };

private:
    const gsGeometry<T> & geo;
    gsMatrix<T>           m_supp;
};



template<class T>
struct gsFieldCreator
{

    static gsField<T> absError(gsField<T> const & field, gsFunction<T> const & f, bool fparam = false)
    {
        const gsMultiPatch<T> & mp = field.patches();
    
        gsPiecewiseFunction<T> * nFields = new gsPiecewiseFunction<T>(mp.nPatches());
        
        for (size_t k=0; k< mp.nPatches(); ++k)
            nFields->addPiecePointer( new gsAbsError<T>(field.function(k), mp.patch(k), f, fparam) );
        
        return gsField<T>(mp, typename gsPiecewiseFunction<T>::Ptr(nFields), true );
    }
    
    
    static gsField<T> gradient(gsMultiPatch<T> const & mp, gsFunction<T> const & f)
    {
        gsPiecewiseFunction<T> * nFields = new gsPiecewiseFunction<T>(mp.nPatches());
        
        
        for (size_t k=0; k< mp.nPatches(); ++k)
            nFields->addPiecePointer( new gsGradientField<T>(mp.patch(k), f) );
        
        return gsField<T>(mp, typename gsPiecewiseFunction<T>::Ptr(nFields), true );
    }

    static gsField<T> normal(gsMultiPatch<T> const & mp)
    {
        gsPiecewiseFunction<T> * nFields = new gsPiecewiseFunction<T>(mp.nPatches());
        
        
        for (size_t k=0; k< mp.nPatches(); ++k)
            nFields->addPiecePointer( new gsNormalField<T>(mp.patch(k)) );
        
        return gsField<T>(mp, typename gsPiecewiseFunction<T>::Ptr(nFields), true );
    }


    static gsField<T> jacDet(gsMultiPatch<T> const & mp)
    {
        gsPiecewiseFunction<T> * nFields = new gsPiecewiseFunction<T>(mp.nPatches());
        
        
        for (size_t k=0; k< mp.nPatches(); ++k)
            nFields->addPiecePointer( new gsJacDetField<T>(mp.patch(k)) );

        return gsField<T>(mp, typename gsPiecewiseFunction<T>::Ptr(nFields), true );
    }

    static gsField<T> parameters(gsMultiPatch<T> const & mp)
    {
        gsPiecewiseFunction<T> * nFields = new gsPiecewiseFunction<T>(mp.nPatches());
        
        
        for (size_t k=0; k< mp.nPatches(); ++k)
            nFields->addPiecePointer( new gsParamField<T>(mp.patch(k)) );

        return gsField<T>(mp, typename gsPiecewiseFunction<T>::Ptr(nFields), true );
    }

    static gsField<T> patchIds(gsMultiPatch<T> const & mp)
    {
        gsPiecewiseFunction<T> * nFields = new gsPiecewiseFunction<T>(mp.nPatches());
        for (size_t k=0; k< mp.nPatches(); ++k)
            nFields->addPiecePointer( new gsPatchIdField<T>(mp.patch(k)) );
        return gsField<T>(mp, typename gsPiecewiseFunction<T>::Ptr(nFields), true );
    }

}; // struct gsFieldCreator



} // namespace gismo