gsNewtonIterator.h

Go to the documentation of this file.

#pragma once

#include <gsAssembler/gsAssembler.h>


namespace gismo
{

template <class T> 
class gsNewtonIterator
{
public:

    gsNewtonIterator(  gsAssembler<T> & assembler,
                       const gsMultiPatch<T> & initialSolution)
    : m_assembler(assembler),
      m_curSolution(initialSolution),
      m_numIterations(0),
      m_maxIterations(100),
      m_tolerance(1e-12),
      m_converged(false)
    { 

    }

    gsNewtonIterator(gsAssembler<T> & assembler)
    : m_assembler(assembler),
      m_numIterations(0),
      m_maxIterations(100),
      m_tolerance(1e-12),
      m_converged(false)
    { 

    }


public:

    void solve();

    void firstIteration();

    void nextIteration();

public:

    const gsMultiPatch<T> & solution() const { return m_curSolution; }

    bool converged() const {return m_converged;}

    index_t numIterations() const { return m_numIterations; }

    T tolerance() const {return m_tolerance;}

    T residue()   const {return m_residue;}

    void setMaxIterations(index_t nIter) {m_maxIterations = nIter;}

    void setTolerance(T tol) {m_tolerance = tol;}

protected:

    virtual void solveLinearProblem(gsMatrix<T> &updateVector);

    virtual void solveLinearProblem(const gsMultiPatch<T> & currentSol, gsMatrix<T> &updateVector);

    virtual T getResidue() {return m_assembler.rhs().norm();}
protected:

    gsAssembler<T> & m_assembler;

    gsMultiPatch<T>     m_curSolution;

    gsMatrix<T>         m_updateVector;

    //gsSparseSolver<>::LU  m_solver;
    //typename gsSparseSolver<T>::BiCGSTABDiagonal m_solver;
    //typename gsSparseSolver<>::CGDiagonal m_solver;
    gsSparseSolver<>::LU  m_solver;

protected:

    index_t m_numIterations;

    index_t m_maxIterations;

    T       m_tolerance;

protected:

    bool m_converged;

    T m_residue;

        T m_updnorm;

};


} // namespace gismo


namespace gismo
{

template <class T>
void gsNewtonIterator<T>::solveLinearProblem(gsMatrix<T>& updateVector)
{
    // Construct the linear system
    m_assembler.assemble();

    // gsDebugVar( m_assembler.matrix().toDense() );
    // gsDebugVar( m_assembler.rhs().transpose() );

    // Compute the newton update
    m_solver.compute( m_assembler.matrix() );
    updateVector = m_solver.solve( m_assembler.rhs() );
    
    // gsDebugVar(updateVector);
}

template <class T>
void gsNewtonIterator<T>::solveLinearProblem(const gsMultiPatch<T> & currentSol, gsMatrix<T>& updateVector)
{
    // Construct linear system for next iteration
    m_assembler.assemble(currentSol);

    // gsDebugVar( m_assembler.matrix().toDense() );
    // gsDebugVar( m_assembler.rhs().transpose() );
    
    // Compute the newton update
    m_solver.compute( m_assembler.matrix() );
    updateVector = m_solver.solve( m_assembler.rhs() );

    // gsDebugVar(updateVector);
}


template <class T> 
void gsNewtonIterator<T>::solve()
{
    firstIteration();

    const T initResidue = m_residue;
        const T initUpdate = m_updnorm;

    // ----- Iterations start -----
    for (m_numIterations = 1; m_numIterations < m_maxIterations; ++m_numIterations)
    {
        gsDebug << "Newton iteration " << m_numIterations 
            << " residue " << math::abs(m_residue)
            << ".\n";
        nextIteration();
        
        // termination criteria
        if ( math::abs(m_updnorm / initUpdate)  < m_tolerance ||
             math::abs(m_residue / initResidue) < m_tolerance )
        {
            m_converged = true;
            break;
        }
    }
}


template <class T> 
void gsNewtonIterator<T>::firstIteration()
{
    // ----- First iteration -----
    m_converged = false;

    // Solve 
    solveLinearProblem(m_updateVector);

    // Construct initial solution
    m_assembler.constructSolution(m_updateVector, m_curSolution);

    // Homogenize Dirichlet dofs (values are now copied in m_curSolution)
    m_assembler.homogenizeFixedDofs(-1);

    // Compute initial residue
    m_residue = getResidue();
    m_updnorm = m_updateVector   .norm();

        gsDebug<<"Iteration: "<< 0
               <<", residue: "<< m_residue
               <<", update norm: "<< m_updnorm
               <<"\n";
}

template <class T> 
void gsNewtonIterator<T>::nextIteration()
{
    // Solve the linaer system of the current iteration
    solveLinearProblem(m_curSolution, m_updateVector);

    // Update the deformed solution
    m_assembler.updateSolution(m_updateVector, m_curSolution);
    
    // Compute residue
    m_residue = getResidue();
    m_updnorm = m_updateVector.norm();
    
    gsDebug<<"Iteration: "<< m_numIterations
           <<", residue: "<< m_residue
           <<", update norm: "<< m_updnorm
           <<"\n";
}



} // namespace gismo