MRLaplacian.h

Go to the documentation of this file.

#pragma once
 
#include "MRBitSet.h"
#include "MRVector.h"
#include "MRVector3.h"
#include "MREnums.h"
#include <MRPch/MREigenSparseCore.h>
 
namespace MR
{
 
 
class Laplacian
{
public:
    enum class RememberShape
    {
        Yes,  
        No    
    };
 
    MRMESH_API explicit Laplacian( Mesh & mesh );
    Laplacian( const MeshTopology & topology, VertCoords & points ) : topology_( topology ), points_( points ) { }
 
    MRMESH_API void init( const VertBitSet & freeVerts, EdgeWeights weights, VertexMass vmass = VertexMass::Unit,
        RememberShape rem = Laplacian::RememberShape::Yes );
 
    MRMESH_API void fixVertex( VertId v, bool smooth = true );
 
    MRMESH_API void fixVertex( VertId v, const Vector3f & fixedPos, bool smooth = true );
 
    MRMESH_API void updateSolver();
 
    MRMESH_API void apply();
 
    MRMESH_API void applyToScalar( VertScalars & scalarField );
 
    [[nodiscard]] const VertBitSet & region() const { return region_; }
 
    [[nodiscard]] const VertBitSet & freeVerts() const { return freeVerts_; }
 
    [[nodiscard]] const VertBitSet & firstLayerFixedVerts() const { assert( solverValid_ ); return firstLayerFixedVerts_; }
 
private:
    void updateSolver_();
 
    void updateRhs_();
 
    template <typename I, typename G, typename S>
    void prepareRhs_( I && iniRhs, G && g, S && s );
 
    const MeshTopology & topology_;
    VertCoords & points_;
 
    VertBitSet region_;
 
    VertBitSet freeVerts_;
 
    VertBitSet fixedSharpVertices_;
 
    VertBitSet firstLayerFixedVerts_;
 
    struct Equation
    {
        Vector3d rhs;           
        double centerCoeff = 0; 
        int firstElem = 0;      
    };
    std::vector<Equation> equations_;
 
    struct Element
    {
        double coeff = 0;
        VertId neiVert;
    };
    std::vector<Element> nonZeroElements_;
 
    Vector< int, VertId > regionVert2id_;
    Vector< int, VertId > freeVert2id_;
 
    using SparseMatrix = Eigen::SparseMatrix<double,Eigen::RowMajor>;
    SparseMatrix M_;
 
    bool solverValid_ = false;
    using SparseMatrixColMajor = Eigen::SparseMatrix<double,Eigen::ColMajor>;
 
    class Solver
    {
    public:
        virtual ~Solver() = default;
        virtual void compute( const SparseMatrixColMajor& A ) = 0;
        virtual Eigen::VectorXd solve( const Eigen::VectorXd& rhs ) = 0;
    };
    std::unique_ptr<Solver> solver_;
 
    bool rhsValid_ = false;
    Eigen::VectorXd rhs_[3];
};
 
}