cellbase.h

Go to the documentation of this file.

// Copyright (C) 1997-2010 by Pawel Pilarczyk.
//
// This file is part of the Homology Library.  This library 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 library 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 software; see the file "license.txt".  If not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started in January 2002. Last revision: August 4, 2010.


#ifndef _CHOMP_CUBES_CELLBASE_H_
#define _CHOMP_CUBES_CELLBASE_H_

#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/cubes/pointset.h"
#include "chomp/homology/chains.h"
#include "chomp/struct/bitfield.h"
#include "chomp/struct/integer.h"
#include "chomp/struct/hashsets.h"
#include "chomp/homology/gcomplex.h"
#include "chomp/cubes/pointbas.h"
#include "chomp/cubes/cubebase.h"
#include "chomp/cubes/cellmain.h"

#include <iostream>
#include <fstream>
#include <cstdlib>

namespace chomp {
namespace homology {


// --------------------------------------------------
// ----------- CubicalCell with PointBase -----------
// --------------------------------------------------

template <class coordtype>
class tCellBase
{
public:
        typedef coordtype CoordType;

        static const int MaxDim = tCubeBase<coordtype>::MaxDim;

        static const int MaxSpaceDim = tCubeBase<coordtype>::MaxDim;

        typedef typename tCubeBase<coordtype>::PointBase PointBase;

        tCellBase ();

        tCellBase (const coordtype *c1, const coordtype *c2,
                int spcdim, int celldim = -1);

        tCellBase (const tCubeBase<coordtype> &q1,
                const tCubeBase<coordtype> &q2);

        tCellBase (const tCubeBase<coordtype> &q, int facedim);

        explicit tCellBase (const tCubeBase<coordtype> &q);

        tCellBase (const tCellBase<coordtype> &q, int offset, int ncoords);

        tCellBase (const tCellBase<coordtype> &c);

        tCellBase<coordtype> &operator =
                (const tCellBase<coordtype> &c);

        int dim () const;

        int spacedim () const;

        coordtype *leftcoord (coordtype *c) const;

        coordtype *rightcoord (coordtype *c) const;

        int_t hashkey1 () const;

        int_t hashkey2 () const;

        static const char *name ();

        static const char *pluralname ();

        enum OutputBitValues
        {
                BitProduct = 0x01,      // unset => two vertices
                BitSpace = 0x02
        };

        static int OutputBits;

        // --- friends: ---

        friend inline int operator == (const tCellBase<coordtype> &c1,
                const tCellBase<coordtype> &c2)
        {
                return (c1. n1 == c2. n1) && (c1. n2 == c2. n2);
        } /* operator == */

private:
        int_t num1 () const;

        int_t num2 () const;

        int_t n1;

        int_t n2;

        void initialize (const coordtype *c1, const coordtype *c2,
                int spcdim, int celldim = -1);

}; /* class tCellBase */

// --------------------------------------------------

template<class coordtype>
int tCellBase<coordtype>::OutputBits = 0;

// --------------------------------------------------

template <class coordtype>
inline int_t tCellBase<coordtype>::num1 () const
{
        return (n1 & NumMask);
} /* tCellBase::num1 */

template <class coordtype>
inline int_t tCellBase<coordtype>::num2 () const
{
        return (n2 & NumMask);
} /* tCellBase::num2 */

template <class coordtype>
inline void tCellBase<coordtype>::initialize (const coordtype *c1,
        const coordtype *c2, int spcdim, int celldim)
{
        // test the validity of the dimension
        if (spcdim <= 0)
                throw "Non-positive dimension of the space.";
        else if (spcdim >= MaxDim)
                throw "Too high dimension of a cell.";

        // get the number of the left vertex
        n1 = PointBase::number (c1, spcdim);
        if (n1 < 0)
                throw "Negative number of the left vertex.";

        // get the number of the right vertex
        n2 = PointBase::number (c2, spcdim);
        if (n2 < 0)
                throw "Negative number of the right vertex.";

        // calculate the dimension of the cubical cell if it is unknown
        if (celldim < 0)
        {
                celldim = 0;
                c1 = PointBase::coord (n1, spcdim);
                c2 = PointBase::coord (n2, spcdim);
                for (int i = 0; i < spcdim; ++ i)
                {
                        if (c2 [i] != c1 [i])
                                ++ celldim;
                }
        }

        // add the space dimension and the cell dimension to the numbers
        n1 |= (static_cast<int_t> (spcdim) << NumBits);
        n2 |= (static_cast<int_t> (celldim) << NumBits);
        return;
} /* tCellBase::initialize */

// --------------------------------------------------

template <class coordtype>
inline tCellBase<coordtype>::tCellBase ()
: n1 (0), n2 (0)
{
        return;
} /* tCellBase::tCellBase */

template <class coordtype>
inline tCellBase<coordtype>::tCellBase
        (const coordtype *c1, const coordtype *c2, int spcdim, int celldim)
{
        initialize (c1, c2, spcdim, celldim);
        return;
} /* tCellBase::tCellBase */

template <class coordtype>
inline tCellBase<coordtype>::tCellBase
        (const tCubeBase<coordtype> &q1, const tCubeBase<coordtype> &q2)
{
        // get the dimension of the space and check for consistency
        int spcdim = q1. dim ();
        if (spcdim != q2. dim ())
                throw "Trying to intersect cubes of different dimension.";

        // get the coordinates of minimal vertices of both cubes
        coordtype c1 [MaxDim];
        q1. coord (c1);
        coordtype c2 [MaxDim];
        q2. coord (c2);

        // prepare tables for the new coordinates of the cubical cell
        coordtype left [MaxDim];
        coordtype right [MaxDim];

        // compute the new coordinates of the cubical cell
        int celldim = CommonCell (left, right, c1, c2, spcdim,
                PointBase::getwrapping (spcdim));
        
        // create the cell as computed
        initialize (left, right, spcdim, celldim);

        return;
} /* tCellBase::tCellBase */

template <class coordtype>
inline tCellBase<coordtype>::tCellBase (const tCubeBase<coordtype> &q,
        int facedim)
: n1 (q. n)
{
        // check the dimensions
        int spcdim = static_cast<int> (n1 >> NumBits);
        if (facedim < 0)
                throw "Negative dimension of a face requested.";
        if (facedim > spcdim)
                throw "Too high dimension of a face requested.";

        // create the opposite corner of the cell
        const coordtype *c1 = PointBase::coord (num1 (), spcdim);
        coordtype c2 [MaxDim];
        for (int i = 0; i < spcdim; ++ i)
                c2 [i] = c1 [i] + ((i < facedim) ? 1 : 0);

        // determine the number of the opposite corner and add cell dimension
        n2 = PointBase::number (c2, spcdim) |
                (static_cast<int_t> (facedim) << NumBits);
        return;
} /* tCellBase::tCellBase */

template <class coordtype>
inline tCellBase<coordtype>::tCellBase (const tCubeBase<coordtype> &q)
: n1 (q. n)
{
        // create the opposite corner of the cell
        int spcdim = static_cast<int> (n1 >> NumBits);
        const coordtype *c1 = PointBase::coord (num1 (), spcdim);
        coordtype c2 [MaxDim];
        for (int i = 0; i < spcdim; ++ i)
                c2 [i] = c1 [i] + 1;

        // determine the number of the opposite corner and add cell dimension
        n2 = PointBase::number (c2, spcdim) |
                (static_cast<int_t> (spcdim) << NumBits);
        return;
} /* tCellBase::tCellBase */

template <class coordtype>
inline tCellBase<coordtype>::tCellBase (const tCellBase<coordtype> &q,
        int offset, int ncoords)
{
        int spcdim = static_cast<int> (q. n1 >> NumBits);
        if ((offset < 0) || (ncoords <= 0) || (offset + ncoords > spcdim))
                throw "Wrong cell projection requested.";
        const coordtype *c1 = PointBase::coord (q. n1 & NumMask, spcdim);
        const coordtype *c2 = PointBase::coord (q. n2 & NumMask, spcdim);
        initialize (c1 + offset, c2 + offset, ncoords);
        return;
} /* tCellBase::tCellBase */

template <class coordtype>
inline tCellBase<coordtype>::tCellBase (const tCellBase<coordtype> &q)
: n1 (q. n1), n2 (q. n2)
{
        return;
} /* tCellBase::tCellBase */

template <class coordtype>
inline tCellBase<coordtype> &tCellBase<coordtype>::operator =
        (const tCellBase<coordtype> &q)
{
        n1 = q. n1;
        n2 = q. n2;
        return *this;
} /* tCellBase::operator = */

template <class coordtype>
inline int tCellBase<coordtype>::dim () const
{
        return static_cast<int> (n2 >> NumBits);
} /* tCellBase::dim */

template <class coordtype>
inline int tCellBase<coordtype>::spacedim () const
{
        return (n1 >> NumBits);
} /* tCellBase::spacedim */

template <class coordtype>
inline coordtype *tCellBase<coordtype>::leftcoord (coordtype *c) const
{
        const coordtype *leftc = PointBase::coord (num1 (), spacedim ());
        if (!c)
                throw "Null pointer to save left coordinates.";
        copycoord (c, leftc, spacedim ());
        return c;
} /* tCellBase::leftcoord */

template <class coordtype>
inline coordtype *tCellBase<coordtype>::rightcoord (coordtype *c) const
{
        const coordtype *rightc = PointBase::coord (num2 (), spacedim ());
        if (!c)
                throw "Null pointer to save right coordinates.";
        copycoord (c, rightc, spacedim ());
        return c;
} /* tCellBase::rightcoord */

template <class coordtype>
inline int_t tCellBase<coordtype>::hashkey1 () const
{
        return ((n1 ^ 0x55555555u) << 17) ^ ((n1 ^ 0xAAAAAAAAu) << 7) ^
                ((n2 ^ 0xAAAAAAAAu) >> 7);
} /* tCellBase::hashkey1 */

template <class coordtype>
inline int_t tCellBase<coordtype>::hashkey2 () const
{
        return ((n2 ^ 0xAAAAAAAAu) << 21) ^ ((n2 ^ 0x55555555u) << 10) ^
                ((n1 ^ 0x55555555u) >> 9);
} /* tCellBase::hashkey2 */

template <class coordtype>
const char *tCellBase<coordtype>::name ()
{
        return "cubical cell";
} /* tCellBase::name */

template <class coordtype>
const char *tCellBase<coordtype>::pluralname ()
{
        return "cubical cells";
} /* tCellBase::pluralname */

// --------------------------------------------------

template <class coordtype>
inline int operator != (const tCellBase<coordtype> &c1,
        const tCellBase<coordtype> &c2)
{
        return !(c1 == c2);
} /* operator != */

// --------------------------------------------------

template <class coordtype>
inline tCellBase<coordtype> operator *
        (const tCellBase<coordtype> &c1,
        const tCellBase<coordtype> &c2)
{
        // get the underlying space dimensions for both cells
        int d1 = c1. spacedim (), d2 = c2. spacedim ();
        if (d1 + d2 >= tCellBase<coordtype>::MaxDim)
                throw "Too high dimension of a Cartesian product of cells.";

        // prepare arrays for the coordinates of the cell to create
        coordtype left [tCellBase<coordtype>::MaxDim];
        coordtype right [tCellBase<coordtype>::MaxDim];

        // extract the coordinates of the first cell
        c1. leftcoord (left);
        c1. rightcoord (right);

        // extract the coordinates of the second cell
        c2. leftcoord (left + d1);
        c2. rightcoord (right + d1);

        // create the Cartesian product of the cells
        return tCellBase<coordtype> (left, right, d1 + d2,
                c1. dim () + c2. dim ());
} /* operator * */

template <class coordtype>
inline std::ostream &operator << (std::ostream &out,
        const tCellBase<coordtype> &c)
{
        return WriteCubicalCell (out, c);
} /* operator << */

template <class coordtype>
inline std::istream &operator >> (std::istream &in,
        tCellBase<coordtype> &c)
{
        return ReadCubicalCell (in, c);
} /* operator >> */

// --------------------------------------------------

template <class coordtype>
inline tCellBase<coordtype> boundarycell (const tCellBase<coordtype> &q,
        int i, bool onlyexisting)
{
        return CubicalBoundaryCell (q, i, onlyexisting);
} /* boundarycell */

template <class coordtype>
inline tCellBase<coordtype> boundarycell (const tCellBase<coordtype> &q,
        int i)
{
        return CubicalBoundaryCell (q, i);
} /* boundarycell */
        
template <class coordtype>
inline int boundarylength (const tCellBase<coordtype> &q)
{
        return CubicalBoundaryLength (q);
} /* boundarylength */

template <class coordtype>
inline int boundarycoef (const tCellBase<coordtype> &q, int i)
{
        return CubicalBoundaryCoef (q, i);
} /* boundarycoef */


} // namespace homology
} // namespace chomp

#endif // _CHOMP_CUBES_CELLBASE_H_