awdiagcubs.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A cubical version of the Alexander-Whitney diagonal,
/// based on simplicial subdivision of cubical cells.
///
/////////////////////////////////////////////////////////////////////////////

// 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: April 11, 2013.


#ifndef _CHAINCON_AWDIAGCUBS_H_
#define _CHAINCON_AWDIAGCUBS_H_


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

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

// include relevant local header files
#include "chaincon/cubcell.h"
#include "chaincon/combtensor.h"
#include "chaincon/tensor.h"
#include "chaincon/combchain.h"
#include "chaincon/chain.h"
#include "chaincon/simplex.h"
#include "chaincon/awdiagsim.h"
#include "chaincon/awdiag.h"


// --------------------------------------------------
// ------------------ AW diagonal -------------------
// --------------------------------------------------

/// Computes the Alexander-Whitney diagonal of a cubical cell
/// using a simplicial subdivision of the cell.
/// Please, note that this formula has only been implemented
/// for 2-dimensional cells. Please, use the other function
/// based on Serre's paper for the general formula.
/// This is a combinatorial version of the function (coefficients in Z_p).
template <class CoordT, class WrapT, class EmptyT>
inline void AWdiagonal (const tCubCell<CoordT,WrapT,EmptyT> &c,
        tCombTensor<tCubCell<CoordT,WrapT,EmptyT>,
        tCubCell<CoordT,WrapT,EmptyT> > &t)
{
        // prepare the data types to be used in this procedure
        typedef tCubCell<CoordT,WrapT,EmptyT> CubeType;
        typedef tSimplex<int_t,EmptyT> SimplexType;

        // abort the program if the dimension is different than two
        if (c. dim () != 2)
        {
                throw "Trying to use a version of A-W diagonal "
                        "implemented for 2-dim cubes only.";
        }

        // extract the coordinates of the four vertices of the cubical cell
        int dim = c. spaceDim ();
        CoordT coord [4] [chomp::homology::MaxBasDim];
        int index = 0;
        for (int i = 0; i < dim; ++ i)
        {
                coord [0] [i] = c. left (i);
                coord [3] [i] = c. right (i);
                if (coord [0] [i] != coord [3] [i])
                {
                        coord [1] [i] = index ? coord [3] [i] :
                                coord [0] [i];
                        coord [2] [i] = index ? coord [0] [i] :
                                coord [3] [i];
                        ++ index;
                }
                else
                {
                        coord [1] [i] = coord [0 ][i];
                        coord [2] [i] = coord [0] [i];
                }
        }

        // retrieve the numbers of the coordinates
        int_t num [4];
        for (int i = 0; i < 4; ++ i)
        {
                num [i] = chomp::homology::tPointBase<CoordT>::number
                        (coord [i], dim);
        }

        // subdivide the cubical cell into simplices
        tCombChain<SimplexType> s;
        int_t numbers [3];
        numbers [0] = num [0];
        numbers [1] = num [1];
        numbers [2] = num [3];
        std::sort (numbers, numbers + 3);
        s. add (SimplexType (2, numbers));
        numbers [0] = num [0];
        numbers [1] = num [2];
        numbers [2] = num [3];
        std::sort (numbers, numbers + 3);
        s. add (SimplexType (2, numbers));

        // compute the AW diagonal of the sum of these two simplices
        tCombTensor<SimplexType,SimplexType> simplTensor;
        AWdiagonal (s, simplTensor);

        // translate the simplices into cubical cells
        // note: extract only the 1-dimensional cells
        int_t size = simplTensor. size ();
        for (int_t n = 0; n < size; ++ n)
        {
                // extract the pair of simplices that form the n-th term
                // of the tensor product
                SimplexType sim [2];
                for (int side = 0; side < 2; ++ side)
                {
                        sim [side] = side ? simplTensor. right (n) :
                                simplTensor. left (n);
                }
                if ((sim [0]. dim () != 1) || (sim [1]. dim () != 1))
                        continue;

                // find the indices of the cubical vertices of the simplices
                int_t vert [2] [2];
                for (int side = 0; side < 2; ++ side)
                {
                        for (int v = 0; v < 2; ++ v)
                        {
                                int_t n = sim [side] [v];
                                int index = 0;
                                while (num [index] != n)
                                        ++ index;
                                vert [side] [v] = index;
                        }
                        if (vert [side] [0] > vert [side] [1])
                                std::swap (vert [side] [0], vert [side] [1]);
                }

                // convert this tensor product of simplices into cells
                tCombChain<CubeType> cc [2];
                for (int side = 0; side < 2; ++ side)
                {
                        if ((vert [side] [0] == 0) && (vert [side] [1] == 3))
                        {
                                cc [side]. add (CubeType (dim,
                                        coord [0], coord [1]));
                                cc [side]. add (CubeType (dim,
                                        coord [1], coord [3]));
                        }
                        else
                        {
                                cc [side]. add (CubeType (dim,
                                        coord [vert [side] [0]],
                                        coord [vert [side] [1]]));
                        }
                }

                // add the tensor products of the cells
                for (int_t i = 0; i < cc [0]. size (); ++ i)
                {
                        for (int_t j = 0; j < cc [1]. size (); ++ j)
                        {
                                t. add (cc [0]. getCell (i),
                                        cc [1]. getCell (j));
                        }
                }
        }

        return;
} /* AWdiagonal */

/// Computes the Alexander-Whitney diagonal of a cubical cell
/// using a simplicial subdivision of the cell.
/// A general version (for any coefficients).
template <class CoordT, class WrapT, class EmptyT, class CoefT>
inline void AWdiagonal (const tCubCell<CoordT,WrapT,EmptyT> &/*c*/,
        tTensor<tCubCell<CoordT,WrapT,EmptyT>,
        tCubCell<CoordT,WrapT,EmptyT>,CoefT> &/*t*/)
{
        throw "Trying to use a version of A-W diagonal "
                "that was not implemented.";
        return;
} /* AWdiagonal */


#endif // _CHAINCON_AWDIAGCUBS_H_