autoincme.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_AUTOINCME_H
#define MECHSYS_FEM_AUTOINCME_H

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

#include "util/string.h"
#include "fem/solver/solver.h"
#include "util/lineparser.h"
#ifdef DO_DEBUG
  #include "util/numstreams.h"
#endif

namespace FEM
{

class AutoIncME : public Solver
{
public:
    AutoIncME(Array<String> const & SolverCtes, FEM::Data & data, FEM::Output & output, Array<int> & nDiv)
        : Solver             (data, output, nDiv) ,
          _dT_ini            (0.001)  ,
          _m_min             (0.01)   ,
          _m_max             (10)     ,
          _m_sf              (0.7)    ,
          _max_sub_steps     (200)    ,
          _dtol              (1.0e-4) ,
          _use_resid_correct (false)  
    {
        for (size_t i=0; i<SolverCtes.size(); ++i)
        {
            LineParser LP(SolverCtes[i]);
            LP.ReplaceAllChars('=',' ');
            String key;                    LP>>key;
                 if (key=="dTini")         LP>>_dT_ini;
            else if (key=="mMin")          LP>>_m_min;
            else if (key=="mMax")          LP>>_m_max;
            else if (key=="mSF")           LP>>_m_sf;
            else if (key=="maxSSteps")     LP>>_max_sub_steps;
            else if (key=="DTOL")          LP>>_dtol;
            else if (key=="correctResid")  LP>>_use_resid_correct;
        }
    }
    void _do_solve_for_an_increment(int                  const   increment   ,  // In
                                    LinAlg::Vector<REAL>       & dF_ext      ,  // In/Out
                                    LinAlg::Vector<REAL>       & dU_ext      ,  // In/Out
                                    REAL                 const   norm_dF_ext ,  // In
                                    LinAlg::Matrix<REAL>       & K           ,  // Out
                                    LinAlg::Vector<REAL>       & dF_int      ,  // Out
                                    LinAlg::Vector<REAL>       & R           ,  // Out
                                    REAL                       & resid       ); // Out
private:
    REAL _dT_ini;
    REAL _m_min;
    REAL _m_max;
    REAL _m_sf; // sf: safety factor
    int  _max_sub_steps;
    REAL _dtol;
    bool _use_resid_correct;
    LinAlg::Vector<REAL> _df_ext;
    LinAlg::Vector<REAL> _du_ext;
    LinAlg::Vector<REAL> _df_int;
    LinAlg::Vector<REAL>   _df_1;
    LinAlg::Vector<REAL>   _df_2;
    LinAlg::Vector<REAL>   _du_1;
    LinAlg::Vector<REAL>   _du_2;
    LinAlg::Vector<REAL>  _du_ME;
    LinAlg::Vector<REAL>  _Err_f;
    LinAlg::Vector<REAL>  _Err_u;
}; // class AutoIncME




inline void AutoIncME::_do_solve_for_an_increment(int                  const   increment   , // In // {{{
                                                  LinAlg::Vector<REAL>       & dF_ext      , // In/Out
                                                  LinAlg::Vector<REAL>       & dU_ext      , // In/Out
                                                  REAL                 const   norm_dF_ext , // In
                                                  LinAlg::Matrix<REAL>       & K           , // Out
                                                  LinAlg::Vector<REAL>       & dF_int      , // Out
                                                  LinAlg::Vector<REAL>       & R           , // Out
                                                  REAL                       & resid       ) // Out
{
#ifdef DO_DEBUG // {{{
    std::cout << "AutoIncME: Increment # " << increment+1 << std::endl;
#endif // }}}

    if (increment==0)
    {
        // Allocate temporary vectors
        int n_tot_dof = dF_ext.Size();
        _df_ext.Resize(n_tot_dof);
        _du_ext.Resize(n_tot_dof);
        _df_int.Resize(n_tot_dof);
          _df_1.Resize(n_tot_dof);
          _df_2.Resize(n_tot_dof);
          _du_1.Resize(n_tot_dof);
          _du_2.Resize(n_tot_dof);
         _du_ME.Resize(n_tot_dof);
         _Err_f.Resize(n_tot_dof);
         _Err_u.Resize(n_tot_dof);

#ifdef DO_DEBUG // {{{
           std::cout << "Sub-step" << _8s()<< "resid" << "Accep/Reject" << _8s()<< "T" << _8s()<< "m" << _8s()<< "dT" << std::endl;
#endif // }}}
    }

    // Auxiliar "pseudo-time" variables
    REAL  T = 0.0;
    REAL dT = _dT_ini;
    REAL m;

    // Loop along sub-steps (0<k<NumSubSteps(?))
    bool DidFinish=false;
    for (int k_sub_step=0; k_sub_step<_max_sub_steps; ++k_sub_step)
    {
        // Check for convergence
        if (T>=1.0) { DidFinish=true; std::cout << std::endl; break; }
        
        // Calculate scaled increment vectors for this sub-step
        LinAlg::CopyScal(dT,dF_ext, _df_ext); // _df_ext <- dT*dF_ext
        LinAlg::CopyScal(dT,dU_ext, _du_ext); // _du_ext <- dT*dU_ext;

        // Setup temporary force vectors
        LinAlg::Copy(_df_ext, _df_1); // _df_1 <- _df_ext
        LinAlg::Copy(_df_ext, _df_2); // _df_2 <- _df_ext;

        // Setup temporary displacement vectors
        LinAlg::Copy(_du_ext, _du_1); // _du_1 <- _du_ext;
        LinAlg::Copy(_du_ext, _du_2); // _du_2 <- _du_ext;

        // Backup element state
        _backup_nodes_and_elements_state(); // _U_bkp<-Node.U; Elem->BackupState();

        // Assemble K matrix and calculate dU_ext
        _assemb_K(K);
        //_inv_K_times_X(K, false, _df_ext, _du_1); // -------------- First evaluation
        _inv_K_times_X(K, false, _df_1, _du_1); // -------------- First evaluation

        // Update nodes and elements state for a Forward-Euler evaluation of displacements
        _update_nodes_and_elements_state(_du_1, _df_int); // Node.U<-Node.U+_du_1; Elem->StressUpdate(_du_1);

        //if (_use_resid_correct) {{{
        //{
            // Calculate residual and its (Euclidian) norm
            //LinAlg::Copy(dF_ext, R);      // R <- dF_ext
            //LinAlg::Axpy(-1.0,dF_int, R); // R <- -1*dF_int+R == R-dF_int (for the first iteration: R<-dF_ext-dF_int)

            //LinAlg::Copy(dU_ext, dU_unb);
            //LinAlg::Axpy(-1.0,dU_1, dU_unb); // ???  dU_unb <- dU_ext - dU_1;

            // Calculate unbalanced displacements
            //_inv_K_times_X(K, false, R, dU_unb); // In=R,U; Out=dU  =>  dU=inv[K(U)]*R

            // Update all elements state
            //_update_elements_state(dU_unb, dF_int);
        //} }}}

        // Assemble K matrix and calculate dU_ext
        _assemb_K(K);
        //_inv_K_times_X(K, false, _df_ext, _du_2); // -------------- Second evaluation
        _inv_K_times_X(K, false, _df_2, _du_2); // -------------- Second evaluation

        // Calculate Modified-Euler displacement increment vector
        LinAlg::AddScaled(0.5,_du_1, 0.5,_du_2, _du_ME); // _du_ME <- 0.5*_du_1 + 0.5*_du_2

        // Calculate local error estimative
        //LinAlg::AddScaled(0.5,_du_2, -0.5,_du_1, _Err); // _Err <- 0.5*_du_2 - 0.5*_du_1
        //REAL norm_Err = sqrt(LinAlg::Dot(_Err,_Err));
        //REAL      err = norm_Err/_norm_essential_vector();
        LinAlg::AddScaled(0.5,_df_2, -0.5,_df_1, _Err_f); // _Err <- 0.5*_du_2 - 0.5*_du_1
        LinAlg::AddScaled(0.5,_du_2, -0.5,_du_1, _Err_u); // _Err <- 0.5*_du_2 - 0.5*_du_1
        REAL norm_Err_f = sqrt(LinAlg::Dot(_Err_f,_Err_f))/sqrt(LinAlg::Dot(_df_2,_df_2));
        REAL norm_Err_u = sqrt(LinAlg::Dot(_Err_u,_Err_u))/sqrt(LinAlg::Dot(_du_2,_du_2));
        REAL err = (norm_Err_f>norm_Err_u ? norm_Err_f : norm_Err_u);

        // Multiplier for the size of next sub-increment
        m = _m_sf*sqrt(_dtol/err);

        // Restore all elements state (useful both if sstep accepted or not)
        _restore_nodes_and_elements_state(); // Node.U<-_U_bkp; Elem->RestoreState();

        // Check if this sub-step can be accepted
        if (err<=_dtol) // sub-step accepted
        {
            // Actualize "pseudo-time"
            T = T + dT;

            // Update nodes and elements state for a Modified-Euler evaluation of displacements
            _update_nodes_and_elements_state(_du_ME, _df_int); // Node.U<-Node.U+_du_ME; Elem->StressUpdate(_du_ME);

            // Calculate residual and its (Euclidian) norm
            LinAlg::AddScaled(1.0,_df_ext, -1.0,_df_int, R); // R <- dF_ext-dF_int
            resid = sqrt(LinAlg::Dot(R,R))/norm_dF_ext;      // normalized measure of the residual

            // Check upper bound of step size
            if (m>_m_max) m=_m_max;
#ifdef DO_DEBUG // {{{
            std::cout << k_sub_step << _8s()<< resid << "Accepted";
#endif // }}}
        }
        else // sub-step rejected
        {
            if (m<_m_min) m=_m_min; // Check lower bound of step size
#ifdef DO_DEBUG // {{{
            std::cout << k_sub_step << -1 << "Rejected";
#endif // }}}
        }

        // Next sub-step
        dT = m*dT;
        if (dT>1.0-T) dT=1.0-T;
//#ifdef DO_DEBUG // {{{
        std::cout << _8s()<< T << _8s()<< m << _8s()<< dT << std::endl;
        //std::cout << _6()<< T <<  m <<  dT << std::endl;
//#endif // }}}
    }
    if (!DidFinish)
        throw new Fatal(_("AutoIncME::_do_solve did not finish with T=1. k_sub_step=<%d>"),_max_sub_steps);
} // }}}




// Allocate a new AutoIncME solver
Solver * AutoIncMEMaker(Array<String> const & SolverCtes, FEM::Data & data, FEM::Output & output, Array<int> & nDiv) // {{{
{
    return new AutoIncME(SolverCtes, data, output, nDiv);
} // }}}

// Register an AutoIncME solver into SolverFactory array map
int AutoIncMERegister() // {{{
{
    SolverFactory["AutoIncME"] = AutoIncMEMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __AutoIncME_dummy_int = AutoIncMERegister();

}; // namespace FEM

#endif // MECHSYS_FEM_AUTOINCME_H

// vim:fdm=marker

---