chaincon/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-2011 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: March 6, 2011.


#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/combchain.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 at the moment; the general formula
/// will be implemented soon.
template <class CoordT>
inline tCombTensor<tCubCell<CoordT> > AWdiagonal (const tCubCell<CoordT> &c)
{
        // abort the program if the dimension is different than two
        if (c. dim () != 2)
                throw "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 indices [2];
        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];
                        indices [index ++] = i;
                }
                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<tSimplex<int_t> > s;
        int_t numbers [3];
        numbers [0] = num [0];
        numbers [1] = num [1];
        numbers [2] = num [3];
        std::sort (numbers, numbers + 3);
        s. add (tSimplex<int_t> (2, numbers));
        numbers [0] = num [0];
        numbers [1] = num [2];
        numbers [2] = num [3];
        std::sort (numbers, numbers + 3);
        s. add (tSimplex<int_t> (2, numbers));

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

        // translate the simplices into cubical cells
        // note: extract only the 1-dimensional cells
        tCombTensor<tCubCell<CoordT> > tensor;
        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
                tSimplex<int_t> 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<tCubCell<CoordT> > cc [2];
                for (int side = 0; side < 2; ++ side)
                {
                        if ((vert [side] [0] == 0) && (vert [side] [1] == 3))
                        {
                                cc [side]. add (tCubCell<CoordT> (dim,
                                        coord [0], coord [1]));
                                cc [side]. add (tCubCell<CoordT> (dim,
                                        coord [1], coord [3]));
                        }
                        else
                        {
                                cc [side]. add (tCubCell<CoordT> (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)
                        {
                                tensor. add (cc [0] [i], cc [1] [j]);
                        }
                }
        }

        return tensor;
} /* AWdiagonal */


#endif // _CHAINCON_AWDIAGCUBS_H_