VectorsInHeat.h

 
#pragma once
 
#if defined(VectorsInHeat_RECURSES)
#error Recursive header files inclusion detected in VectorsInHeat.h
#else // defined(VectorsInHeat_RECURSES)
#define VectorsInHeat_RECURSES
 
#if !defined VectorsInHeat_h
#define VectorsInHeat_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/base/ConstAlias.h"
#include "DGtal/math/linalg/DirichletConditions.h"
 
namespace DGtal
{
// template class VectorInHeat
template <typename TPolygonalCalculus>
class VectorsInHeat
{
    // ----------------------- Standard services ------------------------------
public:
 
    typedef TPolygonalCalculus PolygonalCalculus;
    typedef typename PolygonalCalculus::SparseMatrix SparseMatrix;
    typedef typename PolygonalCalculus::DenseMatrix DenseMatrix;
    typedef typename PolygonalCalculus::Solver Solver;
    typedef typename PolygonalCalculus::Vector Vector;
    typedef typename PolygonalCalculus::Vertex Vertex;
    typedef typename PolygonalCalculus::LinAlg LinAlgBackend;
    typedef DirichletConditions< LinAlgBackend > Conditions;
    typedef typename Conditions::IntegerVector IntegerVector;
 
    VectorsInHeat() = delete;
 
    VectorsInHeat(ConstAlias<PolygonalCalculus> calculus): myCalculus(&calculus)
    {
        myIsInit=false;
    }
 
    ~VectorsInHeat() = default;
 
    VectorsInHeat ( const VectorsInHeat & other ) = delete;
 
    VectorsInHeat ( VectorsInHeat && other ) = delete;
 
    VectorsInHeat & operator= ( const VectorsInHeat & other ) = delete;
 
    VectorsInHeat & operator= ( VectorsInHeat && other ) = delete;
 
 
 
    // ----------------------- Interface --------------------------------------
 
    void init( double dt, double lambda = 1.0,
               bool boundary_with_mixed_solution = false )
    {
        myIsInit=true;
 
        SparseMatrix laplacian = myCalculus->globalLaplaceBeltrami( lambda );
 
        SparseMatrix connectionLaplacian = myCalculus->globalConnectionLaplace( lambda );
 
        SparseMatrix mass = myCalculus->globalLumpedMassMatrix();
        SparseMatrix mass2= myCalculus->doubledGlobalLumpedMassMatrix();
        myScalarHeatOpe   =  mass - dt*laplacian;
        myVectorHeatOpe   =  mass2 - dt*connectionLaplacian;
 
        //Prefactorizing
        myScalarHeatSolver.compute(myScalarHeatOpe);
        myVectorHeatSolver.compute(myVectorHeatOpe);
 
        //empty sources
        myVectorSource          = Vector::Zero(2*myCalculus->nbVertices());
        myScalarSource          = Vector::Zero(myCalculus->nbVertices());
        myDiracSource           = Vector::Zero(myCalculus->nbVertices());
 
        // Manage boundaries
        myManageBoundary = false;
        if ( ! boundary_with_mixed_solution ) return;
        myBoundary = IntegerVector::Zero(myCalculus->nbVertices());
        const auto surfmesh = myCalculus->getSurfaceMeshPtr();
        const auto edges    = surfmesh->computeManifoldBoundaryEdges();
        for ( auto e : edges )
          {
            const auto vtcs = surfmesh->edgeVertices( e );
            myBoundary[ vtcs.first  ] = 1;
            myBoundary[ vtcs.second ] = 1;
          }
        myManageBoundary = ! edges.empty();
        if ( ! myManageBoundary ) return;
        // Prepare solver for a problem with Dirichlet conditions.
        SparseMatrix heatOpe_d = Conditions::dirichletOperator( myScalarHeatOpe, myBoundary );
        // Prefactoring
        myHeatDirichletSolver.compute( heatOpe_d );
    }
 
    void addSource(const Vertex aV,const Vector& ev)
    {
        ASSERT_MSG(aV < myCalculus->nbVertices(), "Vertex is not in the surface mesh vertex range");
        Vector v = myCalculus->Tv(aV).transpose()*ev;
        v = v.normalized()*ev.norm();
        myVectorSource( 2*aV ) = v(0);
        myVectorSource( 2*aV+1 ) = v(1);
        myScalarSource( aV ) = v.norm();
        myDiracSource( aV ) = 1;
    }
 
    void clearSource()
    {
      myVectorSource.clear();
      myScalarSource.clear();
      myDiracSource.clear();
    }
  
    Vector vectorSource() const
    {
        FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
        return myVectorSource;
    }
 
    Vector extrinsicVectorSourceAtVertex(const Vertex aV){
        FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
        return myCalculus->toExtrinsicVector(aV,intrinsicVectorSourceAtVertex(aV));
    }
 
    Vector intrinsicVectorSourceAtVertex(const Vertex aV){
        FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
        Vector s(2);
        s(0) = myVectorSource(2*aV);
        s(1) = myVectorSource(2*aV+1);
        return s;
    }
 
 
    std::vector<Vector> compute() const
    {
        FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
        //Heat diffusion
        Vector vectorHeatDiffusion = myVectorHeatSolver.solve(myVectorSource);
        Vector scalarHeatDiffusion = myScalarHeatSolver.solve(myScalarSource);
        Vector diracHeatDiffusion = myScalarHeatSolver.solve(myDiracSource);
        auto surfmesh = myCalculus->getSurfaceMeshPtr();
 
 
        // Take care of boundaries
        if ( myManageBoundary )
        {
          Vector bValues  = Vector::Zero( myCalculus->nbVertices() );
          Vector bNormSources = Conditions::dirichletVector( myScalarHeatOpe, myScalarSource,
                                                         myBoundary, bValues );
          Vector bSol     = myHeatDirichletSolver.solve( bNormSources );
          Vector heatDiffusionDirichlet
                          = Conditions::dirichletSolution( bSol, myBoundary, bValues );
          scalarHeatDiffusion = 0.5 * ( scalarHeatDiffusion + heatDiffusionDirichlet );
        }
 
        std::vector<Vector> result(surfmesh->nbVertices());
 
        for (auto v = 0;v<surfmesh->nbVertices();v++){
            Vector Y(2);
            Y(0) = vectorHeatDiffusion(2*v);
            Y(1) = vectorHeatDiffusion(2*v+1);
            Y = Y.normalized()*(scalarHeatDiffusion(v)/diracHeatDiffusion(v));
            result[v] = myCalculus->toExtrinsicVector(v,Y);
        }
 
        return result;
    }
 
 
    bool isValid() const
    {
        return myIsInit && myCalculus->isValid();
    }
 
    // ----------------------- Private --------------------------------------
 
private:
 
    const PolygonalCalculus *myCalculus;
 
    SparseMatrix myScalarHeatOpe;
    SparseMatrix myVectorHeatOpe;
 
    Solver myScalarHeatSolver;
    Solver myVectorHeatSolver;
 
    Vector myScalarSource;
    Vector myDiracSource;
    Vector myVectorSource;
 
    bool myManageBoundary;
 
    IntegerVector myBoundary;
 
    bool myIsInit;
 
    Solver myHeatDirichletSolver;
 
}; // end of class VectorsInHeat
} // namespace DGtal
 
//                                                                           //
 
#endif // !defined VectorsInHeat_h
 
#undef VectorsInHeat_RECURSES
#endif // else defined(VectorsInHeat_RECURSES)