cubproduct.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A function that converts a cubical cell into a simplicial set
/// that corresponds to the cartesian product of intervals.
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2009-2016 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 <http://www.gnu.org/licenses/>.

// Started on March 24, 2009. Last revision: August 21, 2013.


#ifndef _CHAINCON_CUBPRODUCT_H_
#define _CHAINCON_CUBPRODUCT_H_


// include relevant local header files
#include "chaincon/prodcell.h"
//#include "chaincon/simplex.h"
#include "chaincon/simplset.h"

// include selected header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/struct/hashsets.h"

// include some standard C++ header files
#include <istream>
#include <ostream>
#include <algorithm>
#include <vector>


// --------------------------------------------------
// ---------------- cube to product -----------------
// --------------------------------------------------

/// Transforms an elementary cube into a simplicial set that represents
/// the Cartesian product of the simplicial sets corresponding to the edges
/// of the cube.
template <class CubCellT, class SetT>
inline void cube2product (const CubCellT &q, SetT &cells)
{
        using chomp::homology::sbug;

        typedef typename SetT::value_type ProdT;
        typedef typename ProdT::CellType SimCellT;
        typedef typename CubCellT::CoordType CoordT;

        // note: SimCellT == tSimplex<CoordT,*,*>

        // create the simplicial sets corresponding to the edges of the cube
        int spaceDim (q. spaceDim ());
        std::vector<tSimplSet<SimCellT> > edgeComplexes (spaceDim);
        for (int i = 0; i < spaceDim; ++ i)
        {
                CoordT coords [2];
                coords [0] = q. left (i);
                coords [1] = q. right (i);
                SimCellT cell (1, coords);
                edgeComplexes [i]. add (cell);
                edgeComplexes [i]. addFaces (1);
                edgeComplexes [i]. addDegenerate (0, spaceDim);
        }

        // create product cells and normalize the result
        std::vector<int_t> counters (spaceDim);
        std::vector<int_t> maxCounters (spaceDim);
        for (int i = 0; i < spaceDim; ++ i)
                maxCounters [i] = edgeComplexes [i]. size ();
//      for (int i = 0; i < spaceDim; ++ i)
//      {
//              sbug << "DEBUG: edgeComplexes [" << i << "]:\n";
//              sbug << edgeComplexes [i];
//      }
        for (int dim = 0; dim <= spaceDim; ++ dim)
        {
        //      sbug << "DEBUG: --- dimension " << dim << " ---\n";
                // prepare the initial counters
                for (int cur = 0; cur < spaceDim; ++ cur)
                {
                        counters [cur] = 0;
                        while ((counters [cur] < maxCounters [cur]) &&
                                (edgeComplexes [cur] [counters [cur]].
                                dim () != dim))
                        {
                                ++ (counters [cur]);
                        }
                        if (counters [cur] >= maxCounters [cur])
                                throw "Missing cells in cube2product init.";
                }
                while (1)
                {
                        // move to the nearest combination of cells of this dim
                        int cur = 0;
                        while ((counters [cur] < maxCounters [cur]) &&
                                (edgeComplexes [cur] [counters [cur]].
                                dim () != dim))
                        {
                                ++ (counters [cur]);
                        }
                        while (counters [cur] >= maxCounters [cur])
                        {
                                counters [cur] = 0;
                                while ((counters [cur] < maxCounters [cur]) &&
                                        (edgeComplexes [cur] [counters [cur]].
                                        dim () != dim))
                                {
                                        ++ (counters [cur]);
                                }
                                if (counters [cur] >= maxCounters [cur])
                                        throw "Missing cells in cube2product.";
                                ++ cur;
                                if (cur >= spaceDim)
                                        break;
                                ++ (counters [cur]);
                                while ((counters [cur] < maxCounters [cur]) &&
                                        (edgeComplexes [cur] [counters [cur]].
                                        dim () != dim))
                                {
                                        ++ (counters [cur]);
                                }
                        }
                        if (cur >= spaceDim)
                                break;

                        // create the product cell
                //      sbug << "DEBUG: counters [0] = " << counters [0] << ".\n";
                //      sbug << "DEBUG: cell [0] = " << edgeComplexes [0] [counters [0]] << ".\n";
                        ProdT prodCell (edgeComplexes [0] [counters [0]]);
                        for (int i = 1; i < spaceDim; ++ i)
                        {
                        //      sbug << "DEBUG: counters [" << i << "] = " << counters [i] << ".\n";
                        //      sbug << "DEBUG: cell [" << i << "] = " << edgeComplexes [i] [counters [i]] << ".\n";
                                ProdT prodNew (prodCell, ProdT
                                        (edgeComplexes [i] [counters [i]]));
                                prodCell. swap (prodNew);
                        }
                        cells. add (prodCell);

                        // take the next combination of cells
                        ++ (counters [0]);
                }
        }

        return;
} /* cube2product */


#endif // _CHAINCON_CUBPRODUCT_H_