gsTensorDomainBoundaryIterator.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsCore/gsDomainIterator.h>
 
#include <gsAssembler/gsGaussRule.h>
 
#include <gsUtils/gsCombinatorics.h>
 
namespace gismo
{
 
// Class which enables iteration over all elements of a tensor product parameter domain
// Documentation in gsDomainIterator.h
 
template<class T, int D, typename uiter>
//template<class T, int D=Dynamic, typename uiter = typename gsKnotVector<T>::const_uiterator>
class gsTensorDomainBoundaryIterator : public gsDomainIterator<T>
{
public:
 
    gsTensorDomainBoundaryIterator( const std::vector< std::vector<T> > & breaks_, const boxSide & s )
    : d( breaks_.size() ),
      lower ( gsVector<T, D>::Zero(d) ),
      upper ( gsVector<T, D>::Zero(d) )
    {
        center =  gsVector<T, D>::Zero(d);
        par = s.parameter();
        dir = s.direction();
        meshBegin.resize(d);
        meshEnd.resize(d);
        curElement.resize(d);
        breaks = breaks_;
 
        for (int i=0; i < dir; ++i)
        {
            meshEnd[i]   = breaks[i].end() - 1;
            meshBegin[i] = curElement[i] = breaks[i].begin();
            if (meshEnd[i] == curElement[i])
                m_isGood = false;
        }
 
        meshEnd[dir]    = ( par ? breaks[dir].end() - 1 : breaks[dir].begin() + 1 );
        curElement[dir] =
        meshBegin[dir]  = ( par ? breaks[dir].end() - 2 : breaks[dir].begin()     );
        tindex = curElement[dir] - breaks[dir].begin();
 
        for (int i=dir+1; i < d; ++i)
        {
            meshEnd[i]   = breaks[i].end() - 1;
            meshBegin[i] = curElement[i] = breaks[i].begin();
 
            if (meshEnd[i] == curElement[i])
                m_isGood = false;
        }
 
        // Set to one quadrature point by default
        m_quadrature.setNodes( gsVector<index_t>::Ones(d) );
 
        if (m_isGood)
            update();
    }
 
    gsTensorDomainBoundaryIterator( const gsBasis<T>& b, const boxSide & s )
    : gsDomainIterator<T>(b, s),
      d( m_basis->dim() ),
      lower ( gsVector<T, D>::Zero(d) ),
      upper ( gsVector<T, D>::Zero(d) )
    {
        center =  gsVector<T, D>::Zero(d);
        par = s.parameter();
        dir = s.direction();
        meshBegin.resize(d);
        meshEnd.resize(d);
        curElement.resize(d);
        breaks.reserve(d);
 
        for (int i=0; i < dir; ++i) 
        {
            breaks.push_back( m_basis->component(i).domain()->breaks() );
            meshEnd[i]   = breaks[i].end() - 1;
            meshBegin[i] = curElement[i] = breaks[i].begin();
            //meshEnd[i]    = m_basis->component(i).domain()->uend() - 1;
            //meshBegin[i]  = 
            //curElement[i] = m_basis->component(i).knots().ubegin();
 
            if (meshEnd[i] == curElement[i])
                m_isGood = false;
        }
 
        // Fixed direction
        breaks.push_back( m_basis->component(dir).domain()->breaks() );
 
        meshEnd[dir]    = ( par ? breaks[dir].end() - 1 : breaks[dir].begin() + 1 );
        curElement[dir] =
        meshBegin[dir]  = ( par ? breaks[dir].end() - 2 : breaks[dir].begin()     );
        tindex = curElement[dir] - breaks[dir].begin();
        //meshEnd[dir]    = ( par ? m_basis->component(dir).knots().uend() - 1   : 
        //                          m_basis->component(dir).knots().ubegin() + 1 );
        //curElement[dir] =
        //meshBegin[dir]  = ( par ? 
        //                    m_basis->component(dir).knots().uend() - 2 : 
        //                    m_basis->component(dir).knots().ubegin()   );
        //tindex = curElement[dir] - m_basis->component(i).knots().ubegin();
 
        for (int i=dir+1; i < d; ++i) 
        {
            breaks.push_back( m_basis->component(i).domain()->breaks() );
            meshEnd[i]   = breaks[i].end() - 1;
            meshBegin[i] = curElement[i] = breaks[i].begin();
            //meshEnd[i]    = m_basis->component(i).knots().uend() - 1;
            //meshBegin[i]  = 
            //curElement[i] = m_basis->component(i).knots().ubegin();
 
            if (meshEnd[i] == curElement[i])
                m_isGood = false;
        }
 
        // Set to one quadrature point by default
        m_quadrature.setNodes( gsVector<index_t>::Ones(d) );
 
        if (m_isGood)
            update();
    }
 
    // ---> Documentation in gsDomainIterator.h
    // proceed to the next element; returns true if end not reached yet
    bool next()
    {
        m_isGood = m_isGood && nextLexicographic(curElement, meshBegin, meshEnd);
        if (m_isGood)
            update();
        return m_isGood;
    }
 
    // ---> Documentation in gsDomainIterator.h
    // proceed to the next element (skipping #increment elements);
    // returns true if end not reached yet
    bool next(index_t increment)
    {
        for (index_t i = 0; i < increment; i++)
            m_isGood = m_isGood && nextLexicographic(curElement, meshBegin, meshEnd);
        if (m_isGood)
            update();
        return m_isGood;
    }
 
    void reset()
    {
        curElement=meshBegin;
        m_isGood = true;
        for(int i=0; i < d; ++i)
        {
            if (i!=dir && curElement[i]==meshEnd[i])
                m_isGood=false;
        }
        if (m_isGood)
            update();
    }
 
    gsVector<unsigned, D> index() const
    {
        gsVector<unsigned, D> curr_index(d);  
        for (int i = 0; i < dir; ++i)
            curr_index[i]  = curElement[i] - meshBegin[i];
        for (int i = dir+1; i < d; ++i)
            curr_index[i]  = curElement[i] - meshBegin[i];
        curr_index[dir]  = tindex;
        return curr_index; 
    }
 
    const gsVector<T> & lowerCorner() const
    { return lower; }
 
    const gsVector<T> & upperCorner() const
    { return upper; }
 
    const T getPerpendicularCellSize() const
    {
        return *(curElement[dir]+1) - *curElement[dir];
    }
 
    size_t numElements() const
    {
        size_t result = 1;
        for (short_t i = 0; i < dir; ++i)
            result *= breaks[i].size() - 1;
        for (short_t i = dir+1; i < d; ++i)
            result *= breaks[i].size() - 1;
        
        return result;
    }
 
    void adjacent( const gsVector<bool> & orient, 
                   gsDomainIterator<T>  & other )
    {
        // 2D only for now
 
        gsTensorDomainBoundaryIterator & other_ = 
            static_cast< gsTensorDomainBoundaryIterator &>(other);
 
        int a1 = !dir;
        int a2 = !other_.dir;
 
        other_.curElement[a2] = std::lower_bound( 
            other_.breaks[a2].begin(), other_.breaks[a2].end(), 
            orient[0] ? *curElement[a1] : *(curElement[a1]+1) );
        other_.update();
    }
 
    void setBreaks(std::vector<T> newBreaks, index_t i) // i: direction
    {
        breaks[i].swap(newBreaks);
        meshEnd[i]   = breaks[i].end() - 1;
        meshBegin[i] = curElement[i] = breaks[i].begin();
        reset();
    }
 
private:
 
    void update()
    {
        for (int i = 0; i < dir ; ++i)
        {
            lower[i]  = *curElement[i];
            upper[i]  = *(curElement[i]+1);
            center[i] = (T)(0.5) * (lower[i] + upper[i]);
        }
        lower[dir]  = 
        upper[dir]  =
        center[dir] = (par ? *(curElement[dir]+1) : *curElement[dir] );
        for (int i = dir+1; i < d; ++i)
        {
            lower[i]  = *curElement[i];
            upper[i]  = *(curElement[i]+1);
            center[i] = (T)(0.5) * (lower[i] + upper[i]);
        }
 
        //gsDebug<<"lower: "<< lower.transpose() <<", upper="<<upper.transpose() <<"\n";
    }
 
// Data members
protected:
    using gsDomainIterator<T>::m_basis;
    using gsDomainIterator<T>::m_isGood;
    using gsDomainIterator<T>::center;
    using gsDomainIterator<T>::m_side;
 
private:
 
    // the dimension of the parameter space
    short_t d;
 
    // Boundary parameters
    short_t  dir;
    bool par;
    unsigned tindex;
 
    // coordinates of the grid cell boundaries
    std::vector< std::vector<T> > breaks;
 
    // Quadrature rule
    gsGaussRule<T> m_quadrature;
 
    // First mesh-line on the tensor grid
    gsVector<uiter, D> meshBegin;
 
    // Last mesh-line on the tensor grid
    gsVector<uiter, D> meshEnd;
 
    // Current element as pointers to it's supporting mesh-lines
    gsVector<uiter, D> curElement;
 
    // parameter coordinates of current grid cell
    gsVector<T> lower, upper;
 
public:
#   define Eigen gsEigen
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
#   undef Eigen
}; // class gsTensorDomainBoundaryIterator
 
 
} // namespace gismo