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, short_t const 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 recoverPoints(gsMatrix<T> & xyz,
                       gsMatrix<T> & uv, index_t k,
                       const T accuracy) 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