newtonraphson.h

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#ifndef MECHSYS_FEM_NEWTON_H
#define MECHSYS_FEM_NEWTON_H

#ifdef HAVE_CONFIG_H
  #include "config.h"
#else
  #ifndef REAL
    #define REAL double
  #endif
#endif

// MechSys
#include "fem/solver/solver.h"

namespace FEM
{

class NewtonRaphson : public Solver
{
public:
    NewtonRaphson(Array<String> const & SolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output)
        : Solver              (ID, data, output),
          _itol               (1.0e-10) ,
          _max_it             (20)      ,
          _use_modfied_newton (false)   ,
          _nSS                (1)
    {
        for (size_t i=0; i<SolverCtes.size(); ++i)
        {
            LineParser LP(SolverCtes[i]);
            LP.ReplaceAllChars('=',' ');
            String key;               LP>>key;
                 if (key=="ITOL")     LP>>_itol;
            else if (key=="maxIt")    LP>>_max_it;
            else if (key=="modified") LP>>_use_modfied_newton;
            else if (key=="nSS")      LP>>_nSS;
            else throw new Fatal(_("AutoME::AutoME: < %s > is invalid "), key.c_str());
        }
        if (_nSS<1) throw new Fatal(_("NewtonRaphson::NewtonRaphson: The number of substeps (nSS=%d) must be greater than or equal to 1"), _nSS);
    }
    void _do_solve_for_an_increment(int                  const   increment,  // In
                                    LinAlg::Vector<REAL> const & DF_ext   ,  // In
                                    LinAlg::Vector<REAL> const & DU_ext   ,  // In
                                    LinAlg::Matrix<REAL>       & K        ,  // (no needed outside)
                                    LinAlg::Vector<REAL>       & dF_int   ,  // (no needed outside)
                                    LinAlg::Vector<REAL>       & Rinternal,  // (no needed outside)
                                    IntSolverData              & ISD      ); // Out
private:
    REAL _itol;
    int  _max_it;
    bool _use_modfied_newton;
    int  _nSS;
    LinAlg::Vector<REAL> _dU;
}; // class NewtonRaphson




inline void NewtonRaphson::_do_solve_for_an_increment(int                  const   increment, // In // {{{
                                                      LinAlg::Vector<REAL> const & DF_ext   , // In
                                                      LinAlg::Vector<REAL> const & DU_ext   , // In
                                                      LinAlg::Matrix<REAL>       & K        , // (no needed outside)
                                                      LinAlg::Vector<REAL>       & dF_int   , // (no needed outside)
                                                      LinAlg::Vector<REAL>       & Rinternal, // (no needed outside)
                                                      IntSolverData              & ISD      ) // Out

{
    // Allocate and fill dF_ext and dU_ext
    int n_tot_dof = DF_ext.Size();
    LinAlg::Vector<REAL> dF_ext; dF_ext.Resize(n_tot_dof);
    LinAlg::Vector<REAL> dU_ext; dU_ext.Resize(n_tot_dof);
    LinAlg::CopyScal(1.0/_nSS,DF_ext, dF_ext); // dF_ext <- DF_ext/_nss
    LinAlg::CopyScal(1.0/_nSS,DU_ext, dU_ext); // dU_ext <- DU_ext/_nss

    // Start
    for (int i=0; i<_nSS; ++i)
    {
        // Check for the first increment
        if (increment==0)
        {
            // Allocate a vector of increment of displacements
            _dU.Resize(dF_ext.Size());
        }

        // Calculate the first residual (=dF_ext)
        LinAlg::Copy(dF_ext, Rinternal); // R <- dF_ext
        REAL resid = 1.0+_itol;          // resid <- norm_R/den

        // Calculate the first (external) disp. increment
        LinAlg::Copy(dU_ext, _dU);  // _dU <- dU_ext

        // Check for Modified-Newton-Raphson scheme
        if (_use_modfied_newton)
            _assemb_K(K);

        // Loop over iterations
        bool converged=false;
        for (int k=0; k<_max_it; ++k)
        {
            // Check for convergence
            if (resid<=_itol) { converged=true; break;}
            
            // Calculate dU_ext using a Forward Euler step
            if (_use_modfied_newton)
                _inv_K_times_X(K, true, Rinternal, _dU); // In=R,U; Out=_dU  =>  _dU=inv[K(U)]*R
            else
            {
                _assemb_K(K);
                _inv_K_times_X(K, false, Rinternal, _dU); // In=R,U; Out=_dU  =>  _dU=inv[K(U)]*R
            }
            
            // Update Nodes and Elements state
            _update_nodes_and_elements_state(_dU, dF_int);

            // Calculate residual (internal)
            LinAlg::Axpy(-1.0,dF_int, Rinternal);  //  R <- R - dF_int (for the first iteration: R<-dF_ext-dF_int)
            _dU.SetValues(0.0);                    // _dU <- 0.0 TODO: check this: after the first iteration there will be no more increment of disp.

            // Calc resid
            REAL denom = 0.0;                      // Normalizer
            for (int i=0; i<Rinternal.Size(); i++) denom += pow((dF_ext(i)+dF_int(i))/2.0, 2.0);
            resid = LinAlg::Norm(Rinternal)/sqrt(denom);
            ISD.BE_resid(resid); 

        }
        if (!converged)
            throw new Fatal(_("NewtonRaphson::_do_solve did not converge for %d iterations"),_max_it);
    }

} // }}}




// Allocate a new NewtonRaphson solver
Solver * NewtonRaphsonMaker(Array<String> const & SolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output) // {{{
{
    return new NewtonRaphson(SolverCtes, ID, data, output);
} // }}}

// Register an NewtonRaphson solver into SolverFactory array map
int NewtonRaphsonRegister() // {{{
{
    SolverFactory["NewtonRaphson"] = NewtonRaphsonMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __NewtonRaphson_dummy_int = NewtonRaphsonRegister();

}; // namespace FEM

#endif // MECHSYS_FEM_NEWTON_H

// vim:fdm=marker

---