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