gsPointGrid.hpp

Go to the documentation of this file.

#pragma once


namespace gismo 
{

/* **************** Utility functions **************** */


template<typename T> //cwiseSampleCount
gsVector<unsigned> uniformSampleCount (const gsVector<T>& lower, 
                                       const gsVector<T>& upper, 
                                       int numPoints)
{
    const index_t d = lower.rows();
    assert( d == upper.rows() );

    // TO do : phys. volume criterion
    gsVector<T> span = upper - lower;
    const T volume = span.prod();
    const T h = math::pow(volume / (T)(numPoints), (T)(1) / (T)(d));

    gsVector<unsigned> np(d);

    for (index_t i = 0; i < d; ++i)
    {
        np[i] = cast<T,unsigned>(math::ceil( span[i] / h ) );
        GISMO_ASSERT( np[i] > 0, "Something went wrong, number of points is zero..");
    }

    return np;
}

template<typename T>
void uniformIntervals(const gsVector<T>& lower, 
                      const gsVector<T>& upper, 
                      std::vector< std::vector<T> >& intervals, 
                      int numIntervals)
{
    const int d = lower.rows();
    assert( d == upper.rows() );

    gsVector<unsigned> np = uniformSampleCount( lower, upper, numIntervals );

    // resize to dimension d without copying old contents, if any
    intervals.clear();
    intervals.resize(d);

    for (int i = 0; i < d; ++i)
    {
        int numInt = np[i] - 1;
        if (numInt <= 1)
            numInt = 1;

        const T h = (T)(1) / (T)(numInt);

        intervals[i].resize(numInt + 1);
        for (int j = 0; j <= numInt; ++j)
            intervals[i][j] = (T)(j) * h;
    }
}


/* **************** Uniform grids described by limits/corners **************** */


template<class T>
gsMatrix<T> gsPointGrid( gsVector<T> const & a, 
                         gsVector<T> const & b, 
                         gsVector<unsigned> const & np )
{
    gsMatrix<T> res(a.size(), np.prod() );
    gsGridIterator<T,CUBE> pt(a, b, np.cast<index_t>());
    for(index_t c = 0; pt; ++pt, ++c)
        res.col(c) = *pt;
    return res;
}

template<typename T>
gsMatrix<T> uniformPointGrid(const gsVector<T>& lower, 
                             const gsVector<T>& upper, 
                             int numPoints)
{
    const gsVector<unsigned> cwisePoints = uniformSampleCount(lower, upper, numPoints);
    return gsPointGrid(lower, upper, cwisePoints); // note: structure lost
}

};// namespace gismo