DiscreteExteriorCalculus.h

#pragma once

#if defined(DiscreteExteriorCalculus_RECURSES)
#error Recursive header files inclusion detected in DiscreteExteriorCalculus.h
#else // defined(DiscreteExteriorCalculus_RECURSES)

#define DiscreteExteriorCalculus_RECURSES

#if !defined DiscreteExteriorCalculus_h

#define DiscreteExteriorCalculus_h

// Inclusions
#include <iostream>
#include <vector>
#include <map>
#include <list>
#include <boost/array.hpp>
#include <boost/unordered_map.hpp>
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/base/Common.h"
#include "DGtal/topology/KhalimskySpaceND.h"
#include "DGtal/dec/Duality.h"
#include "DGtal/dec/KForm.h"
#include "DGtal/dec/LinearOperator.h"
#include "DGtal/dec/VectorField.h"
#include "DGtal/base/ConstAlias.h"
#include "DGtal/topology/CanonicSCellEmbedder.h"

#include <DGtal/kernel/sets/CDigitalSet.h>
#include <DGtal/math/linalg/CDynamicMatrix.h>
#include <DGtal/math/linalg/CDynamicVector.h>
#include <DGtal/math/linalg/CLinearAlgebra.h>

namespace DGtal
{
  // forward factory declaration
  template <typename TLinearAlgebraBackend, typename TInteger>
  class DiscreteExteriorCalculusFactory;

  template <Dimension dim, typename TInteger>
  size_t
  hash_value(const KhalimskyCell<dim, TInteger>& cell);

  // template class DiscreteExteriorCalculus
  template <Dimension dimEmbedded, Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger = DGtal::int32_t>
  class DiscreteExteriorCalculus
  {
    // ----------------------- Standard services ------------------------------
  public:

    friend class DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>;

    typedef DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger> Self;

    typedef TLinearAlgebraBackend LinearAlgebraBackend;
    typedef typename LinearAlgebraBackend::DenseVector::Index Index;
    typedef typename LinearAlgebraBackend::DenseVector::Scalar Scalar;
    typedef typename LinearAlgebraBackend::DenseVector DenseVector;
    typedef typename LinearAlgebraBackend::DenseMatrix DenseMatrix;
    typedef typename LinearAlgebraBackend::SparseMatrix SparseMatrix;

    BOOST_CONCEPT_ASSERT(( concepts::CInteger<TInteger> ));

    BOOST_CONCEPT_ASSERT(( concepts::CDynamicVector<DenseVector> ));
    BOOST_CONCEPT_ASSERT(( concepts::CDynamicMatrix<DenseMatrix> ));
    BOOST_CONCEPT_ASSERT(( concepts::CDynamicMatrix<SparseMatrix> ));
    BOOST_CONCEPT_ASSERT(( concepts::CLinearAlgebra<DenseVector, SparseMatrix> ));
    BOOST_CONCEPT_ASSERT(( concepts::CLinearAlgebra<DenseVector, DenseMatrix> ));

    BOOST_STATIC_ASSERT(( boost::is_same<Dimension, Order>::value ));

    BOOST_STATIC_ASSERT(( dimAmbient >=  dimEmbedded ));

    BOOST_STATIC_CONSTANT( Dimension, dimensionEmbedded = dimEmbedded );
    BOOST_STATIC_CONSTANT( Dimension, dimensionAmbient = dimAmbient );

    typedef DGtal::KhalimskySpaceND<dimAmbient, TInteger> KSpace;
    typedef typename KSpace::Cell Cell;
    typedef typename KSpace::SCell SCell;
    typedef typename KSpace::Point Point;

    struct Property
    {
        Scalar primal_size;
        Scalar dual_size;
        Index index;
        bool flipped;
    };

    typedef boost::unordered_map<Cell, Property> Properties;

    typedef std::vector<SCell> SCells;
    typedef boost::array<SCells, dimEmbedded+1> IndexedSCells;

    typedef VectorField<Self, PRIMAL> PrimalVectorField;
    typedef VectorField<Self, DUAL> DualVectorField;

    typedef KForm<Self, 0, PRIMAL> PrimalForm0;
    typedef KForm<Self, 1, PRIMAL> PrimalForm1;
    typedef KForm<Self, 2, PRIMAL> PrimalForm2;
    typedef KForm<Self, 3, PRIMAL> PrimalForm3;
    typedef KForm<Self, 0, DUAL> DualForm0;
    typedef KForm<Self, 1, DUAL> DualForm1;
    typedef KForm<Self, 2, DUAL> DualForm2;
    typedef KForm<Self, 3, DUAL> DualForm3;

    typedef LinearOperator<Self, 0, PRIMAL, 1, PRIMAL> PrimalDerivative0;
    typedef LinearOperator<Self, 1, PRIMAL, 2, PRIMAL> PrimalDerivative1;
    typedef LinearOperator<Self, 2, PRIMAL, 3, PRIMAL> PrimalDerivative2;
    typedef LinearOperator<Self, 0, DUAL, 1, DUAL> DualDerivative0;
    typedef LinearOperator<Self, 1, DUAL, 2, DUAL> DualDerivative1;
    typedef LinearOperator<Self, 2, DUAL, 3, DUAL> DualDerivative2;

    typedef LinearOperator<Self, 1, PRIMAL, 0, PRIMAL> PrimalAntiderivative1;
    typedef LinearOperator<Self, 2, PRIMAL, 1, PRIMAL> PrimalAntiderivative2;
    typedef LinearOperator<Self, 3, PRIMAL, 2, PRIMAL> PrimalAntiderivative3;
    typedef LinearOperator<Self, 1, DUAL, 0, DUAL> DualAntiderivative1;
    typedef LinearOperator<Self, 2, DUAL, 1, DUAL> DualAntiderivative2;
    typedef LinearOperator<Self, 3, DUAL, 2, DUAL> DualAntiderivative3;

    typedef LinearOperator<Self, 0, PRIMAL, dimEmbedded-0, DUAL> PrimalHodge0;
    typedef LinearOperator<Self, 1, PRIMAL, dimEmbedded-1, DUAL> PrimalHodge1;
    typedef LinearOperator<Self, 2, PRIMAL, dimEmbedded-2, DUAL> PrimalHodge2;
    typedef LinearOperator<Self, 3, PRIMAL, dimEmbedded-3, DUAL> PrimalHodge3;
    typedef LinearOperator<Self, 0, DUAL, dimEmbedded-0, PRIMAL> DualHodge0;
    typedef LinearOperator<Self, 1, DUAL, dimEmbedded-1, PRIMAL> DualHodge1;
    typedef LinearOperator<Self, 2, DUAL, dimEmbedded-2, PRIMAL> DualHodge2;
    typedef LinearOperator<Self, 3, DUAL, dimEmbedded-3, PRIMAL> DualHodge3;

    typedef LinearOperator<Self, 0, PRIMAL, 0, PRIMAL> PrimalIdentity0;
    typedef LinearOperator<Self, 1, PRIMAL, 1, PRIMAL> PrimalIdentity1;
    typedef LinearOperator<Self, 2, PRIMAL, 2, PRIMAL> PrimalIdentity2;
    typedef LinearOperator<Self, 3, PRIMAL, 3, PRIMAL> PrimalIdentity3;
    typedef LinearOperator<Self, 0, DUAL, 0, DUAL> DualIdentity0;
    typedef LinearOperator<Self, 1, DUAL, 1, DUAL> DualIdentity1;
    typedef LinearOperator<Self, 2, DUAL, 2, DUAL> DualIdentity2;
    typedef LinearOperator<Self, 3, DUAL, 3, DUAL> DualIdentity3;

    DiscreteExteriorCalculus();

    template <typename TDomain>
    void
    initKSpace(ConstAlias<TDomain> domain);

    // ----------------------- Iterators on property map -----------------------

    typedef typename Properties::const_iterator ConstIterator;

    ConstIterator begin() const;

    ConstIterator end() const;

    typedef typename Properties::iterator Iterator;

    Iterator begin();

    Iterator end();

    // ----------------------- Interface --------------------------------------
  public:

    KSpace myKSpace;

    void selfDisplay(std::ostream& out) const;

    std::string className() const;

    bool
    insertSCell(const SCell& signed_cell);

    bool
    insertSCell(const SCell& signed_cell, const Scalar& primal_size, const Scalar& dual_size);

    bool
    eraseCell(const Cell& cell);

    void
    updateIndexes();

    const Properties&
    getProperties() const;

    template <Order order, Duality duality>
    const SCells&
    getIndexedSCells() const;

    template <Order order, Duality duality, typename TConstIterator>
    LinearOperator<Self, order, duality, order, duality>
    reorder(const TConstIterator& begin_range, const TConstIterator& end_range) const;

    template <Order order, Duality duality>
    LinearOperator<Self, order, duality, order, duality>
    identity() const;

    template <Order order, Duality duality>
    LinearOperator<Self, order, duality, order+1, duality>
    derivative() const;

    template <Order order, Duality duality>
    LinearOperator<Self, order, duality, order-1, duality>
    antiderivative() const;

    template <Duality duality>
    LinearOperator<Self, 0, duality, 0, duality>
    laplace() const;

    template <Duality duality>
    LinearOperator<Self, 0, duality, 0, duality>
    heatLaplace(const typename DenseVector::Scalar& h, const typename DenseVector::Scalar& t, const typename DenseVector::Scalar& K) const;

    template <Order order, Duality duality>
    LinearOperator<Self, order, duality, dimEmbedded-order, OppositeDuality<duality>::duality>
    hodge() const;

    template <Duality duality>
    KForm<Self, 1, duality>
    flat(const VectorField<Self, duality>& vector_field) const;

    template <Duality duality>
    LinearOperator<Self, 0, duality, 1, duality>
    flatDirectional(const Dimension& dir) const;

    template <Duality duality>
    VectorField<Self, duality>
    sharp(const KForm<Self, 1, duality>& one_form) const;

    template <Duality duality>
    LinearOperator<Self, 1, duality, 0, duality>
    sharpDirectional(const Dimension& dir) const;

    SCell
    getSCell(const Order& order, const Duality& duality, const Index& index) const;

    bool
    isCellFlipped(const Cell& cell) const;

    bool
    containsCell(const Cell& cell) const;

    Index
    getCellIndex(const Cell& cell) const;

    Index
    kFormLength(const Order& order, const Duality& duality) const;

    Order
    actualOrder(const Order& order, const Duality& duality) const;

    Scalar
    hodgeSign(const Cell& cell, const Duality& duality) const;

    Dimension
    edgeDirection(const Cell& cell, const Duality& duality) const;

    void
    resetSizes();

    bool
    isValid() const;

    // ------------------------- Private Datas --------------------------------
  private:

    Properties myCellProperties;

    IndexedSCells myIndexSignedCells;

    boost::array<boost::array<SparseMatrix, dimAmbient>, 2> myFlatOperatorMatrixes;

    boost::array<boost::array<SparseMatrix, dimAmbient>, 2> mySharpOperatorMatrixes;

    bool myCachedOperatorsNeedUpdate;

    bool myIndexesNeedUpdate;


    // ------------------------- Hidden services ------------------------------
  protected:

    // ------------------------- Internals ------------------------------------
  private:

    void
    updateCachedOperators();

    template <Duality duality>
    void
    updateFlatOperator();

    template <Duality duality>
    void
    updateSharpOperator();

  }; // end of class DiscreteExteriorCalculus


  template <Dimension dimEmbedded, Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
  std::ostream&
  operator<<(std::ostream& out, const DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>& object);

} // namespace DGtal


// Includes inline functions.
#include "DGtal/dec/DiscreteExteriorCalculus.ih"

//                                                                           //

#endif // !defined DiscreteExteriorCalculus_h

#undef DiscreteExteriorCalculus_RECURSES
#endif // else defined(DiscreteExteriorCalculus_RECURSES)