MRMesh/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,  // true Laplacian mode when initial mesh shape is remembered and copied in apply
        No    // ignore initial mesh shape in the region and just position vertices smoothly in the region
    };
 
    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:
    // updates solver_ only
    void updateSolver_();
 
    // updates rhs_ only
    void updateRhs_();
 
    template <typename I, typename G, typename S>
    void prepareRhs_( I && iniRhs, G && g, S && s );
 
    const MeshTopology & topology_;
    VertCoords & points_;
 
    // all initially free vertices and the first layer of vertices around them
    VertBitSet region_;
 
    // currently free vertices
    VertBitSet freeVerts_;
 
    // fixed vertices where no smoothness is required
    VertBitSet fixedSharpVertices_;
 
    // fixed vertices from the first layer around free vertices
    VertBitSet firstLayerFixedVerts_;
 
    // for all vertices in the region
    struct Equation
    {
        Vector3d rhs;           // equation right hand side
        double centerCoeff = 0; // coefficient on matrix diagonal
        int firstElem = 0;      // index in nonZeroElements_
    };
    std::vector<Equation> equations_;
 
    struct Element
    {
        double coeff = 0;
        VertId neiVert;
    };
    std::vector<Element> nonZeroElements_;
 
    // map from vertex index to matrix row/col
    Vector< int, VertId > regionVert2id_;
    Vector< int, VertId > freeVert2id_;
 
    using SparseMatrix = Eigen::SparseMatrix<double,Eigen::RowMajor>;
    SparseMatrix M_;
 
    // if true then we do not need to recompute solver_ in the apply
    bool solverValid_ = false;
    using SparseMatrixColMajor = Eigen::SparseMatrix<double,Eigen::ColMajor>;
 
    // interface needed to hide implementation headers
    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_;
 
    // if true then we do not need to recompute rhs_ in the apply
    bool rhsValid_ = false;
    Eigen::VectorXd rhs_[3];
};
 
} //namespace MR