cellbase.h

Go to the documentation of this file.

 
 
 
// Copyright (C) 1997-2020 by Pawel Pilarczyk.
//
// This file is part of my research software package. This 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 3 of the License, or (at your option) any later version.
//
// This software 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,
// please, see <https://www.gnu.org/licenses/>.
 
// Started in January 2002. Last revision: October 24, 2013.
 
 
#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 hashkey1 (const tCellBase<coordtype> &c)
{
        return c. hashkey1 ();
} /* 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>
inline int_t hashkey2 (const tCellBase<coordtype> &c)
{
        return c. hashkey2 ();
} /* 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_