gsGeometry.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsFunction.h>
#include <gsCore/gsBoundary.h>


#define GISMO_BASIS_ACCESSORS \
    Basis & basis() { return static_cast<Basis&>(*this->m_basis); } \
    const Basis & basis() const { return static_cast<const Basis&>(*this->m_basis); }
    // bool isProjective() const{ return Basis::IsRational; }

namespace gismo
{

template<class T>
class gsGeometry : public gsFunction<T>
{

public:
    typedef memory::shared_ptr< gsGeometry > Ptr;

    typedef memory::unique_ptr< gsGeometry > uPtr;

    typedef T Scalar_t;

public:

    
    gsGeometry() :m_basis( NULL ), m_id(0)
    { }

    gsGeometry( const gsBasis<T> & basis, gsMatrix<Scalar_t> coefs) :
    m_basis(basis.clone().release()), m_id(0)
    {
        m_coefs.swap(coefs);
        GISMO_ASSERT( basis.size() == m_coefs.rows(), 
                      "The coefficient matrix of the geometry (rows="<<m_coefs.rows()
                      <<") does not match the number of basis functions in its basis("
                      << basis.size() <<").");
    }

    gsGeometry(const gsGeometry & o);

    gsGeometry& operator=( const gsGeometry & o);
    
    virtual ~gsGeometry();


#if EIGEN_HAS_RVALUE_REFERENCES
    gsGeometry(gsGeometry&& other) 
    : m_coefs(std::move(other.m_coefs)), m_basis(other.m_basis), 
      m_id(std::move(other.m_id))
    {
        other.m_basis = NULL;
    }
    gsGeometry & operator=(gsGeometry&& other)
    {
        m_coefs.swap(other.m_coefs); other.m_coefs.clear();
        delete m_basis;
        m_basis = other.m_basis; other.m_basis = NULL;
        m_id = std::move(other.m_id);
        return *this;
    }
#endif

public:

    // Look at gsFunction class for documentation
    void eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const;

    // Look at gsFunction class for documentation
    virtual void deriv_into(const gsMatrix<T>& u, gsMatrix<T>& result) const;


    // Look at gsFunctionSet class for documentation
    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;

    // Look at gsFunctionSet for documentation
    virtual void compute(const gsMatrix<T> & in, gsFuncData<T> & out) const;


    int orientation() const
    {
        if ( parDim() == geoDim() )
        {
            const T val = gsFunction<T>::jacobian( this->parameterCenter() ).determinant();
            return (T(0) < val) - (val < (T)(0));
        }
        return 1;
    }

    virtual void toggleOrientation();

    /*************************************************************************/


    virtual const gsBasis<T> & basis() const = 0;

    virtual       gsBasis<T> & basis()       = 0;

    short_t targetDim() const { return this->coefDim(); }

    short_t coefDim() const { return static_cast<short_t>(m_coefs.cols()); }

    short_t geoDim() const { return this->coefDim(); }

    virtual short_t domainDim() const;

    short_t parDim() const;

    short_t coDim() const;

    gsMatrix<T> support() const;

    gsMatrix<T> parameterRange() const;

    //boxSide sideOf(const gsVector<T> & u); //


    std::vector<boxSide> locateOn(const gsMatrix<T> & u, gsVector<bool> & onG2, gsMatrix<T> & preIm, bool lookForBoundary = false, real_t tol = 1.e-6) const; //

    // Whether the coefficients of this geometry are stored in projective or affine form
    //virtual bool isProjective() const = 0;


    /*************************************************************************/

    // todo: coefsSize() x (geoDim() + 1) if projective
          gsMatrix<T> & coefs()       { return this->m_coefs; }

    const gsMatrix<T> & coefs() const { return this->m_coefs; }

    typename gsMatrix<T>::RowXpr       coef(index_t i)       { return m_coefs.row(i); }

    typename gsMatrix<T>::ConstRowXpr  coef(index_t i) const { return m_coefs.row(i); }

    T & coef(index_t i, index_t j)
    {
        GISMO_ASSERT( ((i < m_coefs.rows()) && (j < m_coefs.cols()) ),
                      "Coefficient or coordinate which is out of range requested.");
        return m_coefs(i,j);
    }

    const T coef(index_t i, index_t j) const
    {
        GISMO_ASSERT( ((i < m_coefs.rows()) && (j < m_coefs.cols()) ),
                      "Coefficient or coordinate which is out of range requested.");
        return m_coefs(i,j);
    }

    void setCoefs(gsMatrix<T> cc) { this->m_coefs.swap(cc); }

    index_t coefsSize() const { return m_coefs.rows(); }
    // Warning: This can cause some clash while using periodic basis, since the ghost coefs are stored but ignored.




    /*************************************************************************/

    void linearTransform(const gsMatrix<T>& mat)
    {
        this->m_coefs = this->m_coefs * mat.transpose();
    }

    void rotate(T angle, const gsVector<T,3> & axis )
    {
        assert( geoDim() == 3 );
        gsEigen::Transform<T,3,gsEigen::Affine> 
            rot( gsEigen::AngleAxis<T> (angle,axis.normalized()) );
        // To do: Simpler way to use transforms ?
        this->m_coefs = (this->m_coefs.rowwise().homogeneous() * 
                         rot.matrix().transpose() ).leftCols(3) ;
    }

    void rotate(T angle)
    {
        GISMO_ASSERT( geoDim() == 2, "Only for 2D");
        gsEigen::Rotation2D<T> rot(angle);
        this->m_coefs *= rot.matrix().transpose();
    }

    void scale(T s, int coord = -1)
    {
        if ( coord == -1) // Uniform scaling
            this->m_coefs *= s;
        else if ( coord <geoDim() )//scale coordinate coord
            this->m_coefs.col(coord) *= s;
    }

    void scale(gsVector<T> const & v)
    {
        GISMO_ASSERT( v.rows() == this->m_coefs.cols(), "Sizes do not agree." );
        this->m_coefs.array().rowwise() *= v.array().transpose();
    }

    void translate(gsVector<T> const & v)
    {
        this->m_coefs.rowwise() += v.transpose();
    }

    typename gsMatrix<T>::RowXpr
    coefAtCorner(boxCorner const & c);

    typename gsMatrix<T>::ConstRowXpr
    coefAtCorner(boxCorner const & c) const;


    /*************************************************************************/


    virtual void uniformRefine(int numKnots = 1, int mul=1, int dir=-1); // todo: int dir = -1

    virtual void uniformCoarsen(int numKnots = 1); // todo: int dir = -1

    void refineElements( std::vector<index_t> const & boxes );

    void unrefineElements( std::vector<index_t> const & boxes );

    typename gsGeometry::uPtr coord(const index_t c) const;
    
    void embed3d()
    {
        embed(3);
    }

    void embed(index_t N, bool pad_right = true)
    {
        GISMO_ASSERT( N > 0, "Embed dimension must be positive");

        const index_t nc = N - m_coefs.cols();
        if ( nc == 0 ) return;

        if (!pad_right && nc<0)
            m_coefs.leftCols(N) = m_coefs.rightCols(N);
        m_coefs.conservativeResize(gsEigen::NoChange, N);

        if ( nc > 0 )
        {
            if (pad_right)
                m_coefs.rightCols(nc).setZero();
            else
            {
                m_coefs.rightCols(N-nc) = m_coefs.leftCols(N-nc);
                m_coefs.leftCols(nc).setZero();
            }
        }
#       ifndef NDEBUG
        else gsWarn<<"Coefficients projected (deleted)..\n";
#       endif
    }

    short_t degree(const short_t & i) const;

    virtual void insertKnot( T knot, index_t dir, index_t i = 1);

    virtual void degreeElevate(short_t const i = 1, short_t const dir = -1);

    virtual void degreeIncrease(short_t const i = 1, short_t const dir = -1);

    virtual void degreeReduce(short_t const i = 1, short_t const dir = -1);

    virtual void degreeDecrease(short_t const i = 1, short_t const dir = -1);
    
    virtual void hessian_into(const gsMatrix<T>& u, gsMatrix<T> & result,
                              index_t coord) const;
    
    void controlNet( gsMesh<T> & mesh) const;

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

    std::vector<gsGeometry *> boundary() const;

    virtual typename gsGeometry::uPtr boundary(boxSide const& s) const;

    virtual typename gsGeometry::uPtr iface(const boundaryInterface & bi,
                                            const gsGeometry & other) const;

    typename gsGeometry::uPtr component(boxComponent const& bc) const;

    GISMO_UPTR_FUNCTION_PURE(gsGeometry, clone)

    
    virtual std::ostream &print(std::ostream &os) const;

    virtual void merge( gsGeometry * other );

    virtual void toMesh(gsMesh<T> & msh, int npoints) const;

    typename gsGeometry::uPtr approximateLinearly() const;

    void evaluateMesh(gsMesh<T>& mesh) const;


    virtual std::vector<gsGeometry<T>* > uniformSplit(index_t dir =-1 ) const;


    gsGeometrySlice<T> getIsoParametricSlice(index_t dir_fixed, T par) const;

    T closestPointTo(const gsVector<T> & pt,
                        gsVector<T> & result,
                        const T accuracy = 1e-6,
                        const bool useInitialPoint = false) const;

    T directedHausdorffDistance(const gsGeometry & other,
                                const index_t nsamples = 1000,
                                const T accuracy = 1e-6) const;

    T HausdorffDistance(        const gsGeometry & other,
                                const index_t nsamples = 1000,
                                const T accuracy = 1e-6,
                                const bool directed=false) const;

    void setId(const size_t i) { m_id = i; }

    size_t id() const { return m_id; }

protected:
    void swap(gsGeometry & other)
    {
        std::swap(m_basis, other.m_basis);
        m_coefs.swap(other.m_coefs);
        std::swap(m_id, other.m_id);
    }

protected:

    //todo: coefsSize() x (geoDim() + 1) if projective
    gsMatrix<T> m_coefs;

    gsBasis<T> * m_basis;

    size_t m_id;

}; // class gsGeometry

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

} // namespace gismo

#include <gsCore/gsCurve.h>
#include <gsCore/gsSurface.h>
#include <gsCore/gsVolume.h>
#include <gsCore/gsBulk.h>

namespace gismo
{

// Generic traits for geometry with dimension known at compile time
template <short_t d, typename T>
struct gsGeoTraits
{
    typedef gsGeometry<T> GeometryBase;
};

// Traits for curve
template <typename T>
struct gsGeoTraits<1,T>
{
    typedef gsCurve<T> GeometryBase;
};

// Traits for surface
template <typename T>
struct gsGeoTraits<2,T>
{
    typedef gsSurface<T> GeometryBase;
};

// Traits for volume
template <typename T>
struct gsGeoTraits<3,T>
{
    typedef gsVolume<T> GeometryBase;
};

// Traits for bulk
template <typename T>
struct gsGeoTraits<4,T>
{
    typedef gsBulk<T> GeometryBase;
};

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsGeometry(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


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