particlegrid.h

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2013   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


#ifndef particlegrid_h
#define particlegrid_h

#include "floatarray.h"
#include "intarray.h"
#include "mathfem.h"

#include <list>

namespace oofem
{

template< class Point >
class ParticleGridIterator;

template <class Point>
class OOFEM_EXPORT ParticleGrid
{
public:
    typedef ParticleGridIterator<Point> iterator;

protected:
    struct RefinedParticlePoint {
        Point *point;
        ParticleGrid<Point> *subgrid;
    };

    int n;
    IntArray res;
    int total;
    FloatArray bb0,bb1;
    FloatArray dx;
    RefinedParticlePoint **data;
    IntArray res_prod;

    int giveIndex ( const IntArray &pos ) const;
    void givePosition ( IntArray &pos, int index ) const;

    void init ( const IntArray &res, const FloatArray &bb0, const FloatArray &bb1 );

public:
    ParticleGrid ( const IntArray &res, const FloatArray &bb0, const FloatArray &bb1 );
    ParticleGrid ( const ParticleGrid *old );
    ~ParticleGrid();

    int giveDimensions() {
        return n;
    }
    int getTotal() const {
        return this->total;
    }
    int getNumberOfPoints() const;
    bool isEmpty() const;
    iterator begin();
    iterator beginAt ( const FloatArray &x0, const FloatArray &x1 );

    int getResolution ( int i ) {
        return this->res ( i );
    }
    double getGridStep ( int i ) const {
        return this->dx ( i );
    }

    void getGridCoord ( FloatArray &answer, const IntArray &pos ) const;
    void getGridCoord ( FloatArray &answer, int ind ) const;

    bool getPoint ( Point *&answer, const IntArray &pos ) const;
    bool getPoint ( Point *&answer, int ind ) const;

    bool getSubGrid ( ParticleGrid<Point> *&subgrid, const IntArray &pos ) const;
    bool getSubGrid ( ParticleGrid<Point> *&subgrid, int ind ) const;

    void createSubGrid ( ParticleGrid<Point> *&subgrid, const IntArray &res, const IntArray &pos );

    void setPoint ( Point *point, const IntArray &pos );
    void setPoint ( Point *point, int ind );

    void clearPosition ( const IntArray &pos );
    void clearPosition ( int ind );

    void getPointsWithin ( std::list<Point*> &answer, const FloatArray &x0, const FloatArray &x1 );
    void appendAllPoints ( std::list<Point*> &answer );

    void getBoundingBox ( const FloatArray &x0, const FloatArray &x1, IntArray &ind0, IntArray &ind1 ) const;

    void getPosition ( const FloatArray &x, IntArray &pos ) const;

    std :: string errorInfo(const char *func) const { return std :: string("ParticleGrid :: ") + func; }

    friend class ParticleGridIterator<Point>;
};

template <class Point>
ParticleGrid<Point> :: ParticleGrid ( const IntArray &res, const FloatArray &bb0, const FloatArray &bb1 )
{
    this->init ( res, bb0, bb1 );
}

template <class Point>
ParticleGrid<Point> :: ParticleGrid ( const ParticleGrid *old )
{
    this->init ( old->res, old->bb0, old->bb1 );
}

template <class Point>
void ParticleGrid<Point> :: init ( const IntArray &res, const FloatArray &bb0, const FloatArray &bb1 )
{
    this->res = res;
    this->bb0 = bb0;
    this->bb1 = bb1;

    this->n = this->res.giveSize();
    this->dx.resize ( this->n );
    this->total = 1;
    this->res_prod.resize ( this->n );
    this->res_prod ( 0 ) = 1;

    for ( int i = 0; i < this->n; ++i ) {
        this->dx ( i ) = ( this->bb1 ( i ) - this->bb0 ( i ) ) / ( this->res ( i )-1 );

        this->total *= this->res ( i );

        this->res_prod ( i ) = 1;
        for ( int j = 0; j < i; ++j ) {
            this->res_prod ( i ) *= this->res ( j );
        }
    }

    this->data = new RefinedParticlePoint*[this->total];
    for ( int i = 0; i < this->total; ++i ) {
        this->data[i] = NULL;
    }
}

template <class Point>
ParticleGrid<Point> :: ~ParticleGrid()
{
    // Clear all points
    for ( int i = 0; i < this->total; ++i ) {
        this->clearPosition ( i );
    }
    delete[] this->data;
}

template <class Point>
int ParticleGrid<Point> :: getNumberOfPoints() const
{
    RefinedParticlePoint *p;
    int c = 0;
    for ( int i = 0; i < this->total; ++i ) {
        p = this->data[i];
        if ( p ) {
            if ( p->point ) {
                c += 1;
            } else if ( p->subgrid ) {
                c += p->subgrid->getNumberOfPoints();
            }
        }
    }
    return c;
}

template <class Point>
bool ParticleGrid<Point> :: isEmpty() const
{
    RefinedParticlePoint *p;
    for ( int i = 0; i < this->total; ++i ) {
        p = this->data[i];
        if ( p ) {
            if ( p->point ) {
                return false;
            } else if ( p->subgrid ) {
                if ( !p->subgrid->isEmpty() ) {
                    return false;
                }
            }
        }
    }
    return true;
}

template <class Point>
int ParticleGrid<Point> :: giveIndex ( const IntArray &pos ) const
{
    int index = 0;
    for ( int i = 0; i < this->n; ++i ) {
        index += this->res_prod ( i ) *pos ( i );
    }
    return index;
}

template <class Point>
void ParticleGrid<Point> :: givePosition ( IntArray &pos, int index ) const
{
    pos.resize ( this->n );
    for ( int i = 0; i < this->n; ++i ) {
        pos ( i ) = ( index/res_prod ( i ) ) % res ( i );
    }
}

template <class Point>
void ParticleGrid<Point> :: getGridCoord ( FloatArray &answer, const IntArray &pos ) const
{
    answer.resize ( this->n );
    for ( int i = 0; i < this->n; ++i ) {
        answer ( i ) = pos ( i ) *this->dx ( i ) + this->bb0 ( i );
    }
}

template <class Point>
void ParticleGrid<Point> :: getGridCoord ( FloatArray &answer, int ind ) const
{
    IntArray pos;
    this->givePosition ( pos, ind );
    this->getGridCoord ( answer, pos );
}

template <class Point>
void ParticleGrid<Point> :: getBoundingBox ( const FloatArray &x0, const FloatArray &x1, IntArray &ind0, IntArray &ind1 ) const
{
    ind0.resize ( this->n );
    ind1.resize ( this->n );
    for ( int i = 0; i < this->n; ++i ) {
        ind0 ( i ) = max ( ( int ) floor ( ( x0 ( i ) - this->bb0 ( i ) ) /this->dx ( i ) ), 0 );
        ind1 ( i ) = min ( ( int ) ceil ( ( x1 ( i ) - this->bb0 ( i ) ) /this->dx ( i ) ), this->res ( i ) );
    }
}

template <class Point>
void ParticleGrid<Point> :: getPosition ( const FloatArray &x, IntArray &pos ) const
{
    pos.resize ( this->n );
    for ( int i = 0; i < this->n; ++i ) {
        //pos(i) = (int)nearest((x(i) - this->bb0(i))/this->dx(i));
        pos ( i ) = ( int ) floor ( ( x ( i ) - this->bb0 ( i ) ) /this->dx ( i ) + 0.5 );
        pos ( i ) = min ( max ( pos ( i ),0 ),this->res ( i )-1 );
    }
}

template <class Point>
typename ParticleGrid<Point>::iterator ParticleGrid<Point> :: begin()
{
    iterator it ( this );
    return it;
}

template <class Point>
typename ParticleGrid<Point>::iterator ParticleGrid<Point> :: beginAt ( const FloatArray &x0, const FloatArray &x1 )
{
    IntArray ind0, ind1;
    this->getBoundingBox ( x0, x1, ind0, ind1 );
    iterator it ( this, ind0, ind1 );
    return it;
}

template <class Point>
void ParticleGrid<Point> :: getPointsWithin ( std::list<Point*> &answer, const FloatArray &x0, const FloatArray &x1 )
{
    answer.clear();

    // First find a suitable bounding box of positions (note: this isn't perfect in regards to sub grids)
    IntArray p0, p1;
    this->getBoundingBox ( x0, x1, p0, p1 );

    // Get all grid points in that region
    IntArray p ( this->n );
    if ( this->n == 2 ) {
        for ( int x = p0 ( 0 ); x < p1 ( 0 ); x++ ) {
            p ( 0 ) = x;
            for ( int y = p0 ( 1 ); y < p1 ( 1 ); y++ ) {
                p ( 1 ) = y;
                RefinedParticlePoint *rp = this->data[this->giveIndex ( p )];
                if ( rp ) {
                    if ( rp->subgrid ) {
                        rp->subgrid->appendAllPoints ( answer );
                    } else {
                        answer.push_front ( rp->point );
                    }
                }
            }
        }
    } else if ( this->n == 3 ) {
        OOFEM_ERROR ( "3D grids not implemented yet" );
    }
}

template <class Point>
void ParticleGrid<Point> :: appendAllPoints ( std::list<Point*> &answer )
{
    for ( int i = 0; i < this->total; ++i ) {
        RefinedParticlePoint *rp = this->data[i];
        if ( rp ) {
            if ( rp->subgrid ) {
                rp->subgrid->appendAllPoints ( answer );
            } else if ( rp->point ) {
                answer.push_front ( rp->point );
            } else {
                OOFEM_ERROR ( "Refined particle does not contain subgrid or point." );
            }
        }
    }
}

template <class Point>
bool ParticleGrid<Point> :: getPoint ( Point *&answer, const IntArray &pos ) const
{
    return this->getPoint ( answer, this->giveIndex ( pos ) );
}

template <class Point>
void ParticleGrid<Point> :: setPoint ( Point *point, const IntArray &pos )
{
    this->setPoint ( point, this->giveIndex ( pos ) );
}

template <class Point>
bool ParticleGrid<Point> :: getPoint ( Point *&answer, int ind ) const
{
    RefinedParticlePoint *rp = this->data[ind];
    if ( rp == NULL ) { // Not active
        answer = NULL;
        return false;
    } else {
        if ( rp->point ) { // Active and has a point
            answer = rp->point;
            return true;
        } else { // Active and has a sub grid
            answer = NULL;
            return false;
        }
    }
}

template <class Point>
void ParticleGrid<Point> :: setPoint ( Point *point, int ind )
{
    this->clearPosition ( ind );
    this->data[ind] = new RefinedParticlePoint;
    this->data[ind]->point = point;
    this->data[ind]->subgrid = NULL;
}

template <class Point>
void ParticleGrid<Point> :: clearPosition ( const IntArray &pos )
{
    return this->clearPosition ( this->giveIndex ( pos ) );
}

template <class Point>
void ParticleGrid<Point> :: clearPosition ( int ind )
{
    RefinedParticlePoint *rp = this->data[ind];
    if ( rp ) {
        if ( rp->point ) {
            delete rp->point;
        }
        if ( rp->subgrid ) {
            delete rp->subgrid;
        }
        delete rp;
    }
    this->data[ind] = NULL;
}

template <class Point>
bool ParticleGrid<Point> :: getSubGrid ( ParticleGrid<Point> *&subgrid, int ind ) const
{
    RefinedParticlePoint *rp = this->data[ind];
    if ( rp == NULL ) {
        subgrid = NULL;
        return false;
    } else {
        if ( rp->subgrid ) {
            subgrid = rp->subgrid;
            return true;
        } else {
            subgrid = NULL;
            return false;
        }
    }
}

template <class Point>
bool ParticleGrid<Point> :: getSubGrid ( ParticleGrid<Point> *&subgrid, const IntArray &pos ) const
{
    return this->getSubGrid ( subgrid, this->giveIndex ( pos ) );
}

template <class Point>
void ParticleGrid<Point> :: createSubGrid ( ParticleGrid<Point> *&subgrid, const IntArray &res, const IntArray &pos )
{
    // The bounding boxes for the sub domains.
    FloatArray grid_point;
    this->getGridCoord ( grid_point, pos );
    FloatArray sub_bb0, sub_bb1;
    sub_bb0 = grid_point;
    sub_bb0.add ( -0.5, dx );
    sub_bb1 = grid_point;
    sub_bb1.add ( 0.5, dx );
    subgrid = new ParticleGrid ( res, sub_bb0, sub_bb1 );

    int ind = this->giveIndex ( pos );
    this->clearPosition ( ind );
    this->data[ind] = new RefinedParticlePoint;
    this->data[ind]->point = NULL;
    this->data[ind]->subgrid = subgrid;
}

template <class Point>
class ParticleGridIterator
{
protected:
    int index, endind;
    ParticleGridIterator<Point> *sub_it;
    ParticleGrid<Point> *grid;

    IntArray pos;
    const IntArray ind0, ind1;
    bool limited;

public:
    ParticleGridIterator() :
        index ( 0 ), endind ( 0 ), sub_it ( NULL ), grid ( NULL ), ind0(), ind1(), limited ( false ) {}
    ParticleGridIterator ( ParticleGrid<Point> *g ) :
        index ( 0 ), endind ( g->total ), sub_it ( NULL ), grid ( g ), ind0(), ind1(), limited ( false ) {
        if ( this->grid->data[this->index] && this->grid->data[this->index]->subgrid ) {
            this->sub_it = new ParticleGridIterator<Point> ( this->grid->data[this->index]->subgrid );
        }
    }
    ParticleGridIterator ( ParticleGrid<Point> *g, IntArray ind0, IntArray ind1 ) :
        index ( 0 ), sub_it ( NULL ), grid ( g ), ind0 ( ind0 ), ind1 ( ind1 ), limited ( true ) {
        IntArray end = ind1;
        end.add ( -1 );
        this->endind = this->grid->giveIndex ( end ) + 1;
        this->pos = ind0;
        this->index = this->grid->giveIndex ( ind0 );
        if ( this->grid->data[this->index] && this->grid->data[this->index]->subgrid ) {
            this->sub_it = new ParticleGridIterator<Point> ( this->grid->data[this->index]->subgrid );
        }
    }
    ParticleGridIterator ( const ParticleGridIterator<Point> &x ) :
        index ( x.index ), sub_it ( x.sub_it ), grid ( x.grid ), ind0 ( x.ind0 ), ind1 ( x.ind1 ), limited ( x.limited )  { }
    ~ParticleGridIterator() {
        if ( this->sub_it ) {
            delete sub_it;
        }
    }

    bool end() const {
        return this->index >= this->endind;
    }

    Point *getPoint() {
        if ( this->end() ) {
            return NULL;
        }
        if ( this->subGridActive() ) {
            return this->sub_it->getPoint();
        }
        if ( this->grid->data[this->index] != NULL ) {
            return this->grid->data[this->index]->point;
        }
        return NULL;
    }

    int getIndex() {
        return this->index;
    }

    void setPoint ( Point *p ) {
        if ( this->end() ) {
            OOFEM_ERROR ( "Can't set element in outside grid" );
        }
        if ( this->subGridActive() ) {
            this->sub_it->setPoint ( p );
        }
        this->grid->setPoint ( p, this->index );
    }

    bool subGridActive() const {
        return this->sub_it != NULL;
    }

    void getGridPoint ( FloatArray &x ) {
        if ( this->subGridActive() ) {
            this->sub_it->getGridPoint ( x );
        } else {
            this->grid->getGridCoord ( x, this->index );
        }
    }

    void operator++() {
        if ( this->subGridActive() ) {
            this->sub_it++;
            if ( !this->sub_it->end() ) {
                return;
            } else {
                delete this->sub_it;
                this->sub_it = NULL;
            }
        }
        this->index++;
        if ( this->limited ) {
            this->pos ( 0 ) ++;
            for ( int i = 0; i < this->pos.giveSize() - 1; ++i ) {
                if ( this->pos ( i ) >= this->ind1 ( i ) ) {
                    this->pos ( i ) = this->ind0 ( i );
                    this->pos ( i+1 ) ++;
                } else {
                    break;
                }
            }
            index = this->grid->giveIndex ( pos );
        }
        if ( !this->end() && this->grid->data[this->index] && this->grid->data[this->index]->subgrid ) {
            this->sub_it = new ParticleGridIterator<Point> ( this->grid->data[this->index]->subgrid );
        }
    }

    std :: string errorInfo(const char *func) const { return std :: string("ParticleGridIterator :: ") + func; }

    friend class ParticleGrid<Point>;
};

} // end namespace oofem
#endif // particlegrid_h