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)