autome.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_AUTOME_H
#define MECHSYS_FEM_AUTOME_H

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

// MechSys
#include "util/string.h"
#include "util/lineparser.h"
#include "fem/solver/solver.h"

namespace FEM
{

class AutoME : public Solver
{
public:
    AutoME(Array<String> const & SolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output)
        : Solver (ID, data, output),
          _maxSS (200)   ,
          _STOL  (1.0e-2),
          _dTini (0.001) ,
          _mMin  (0.01)  ,
          _mMax  (10)    ,
          _mCoef (0.7)   ,
          _ZTOL  (1e-3)  ,
          _unbal_correction  (false),
          _check_convergence (true)
    {
        for (size_t i=0; i<SolverCtes.size(); ++i)
        {
            LineParser LP(SolverCtes[i]);
            LP.ReplaceAllChars('=',' ');
            String key;            LP>>key;
                 if (key=="maxSS") LP>>_maxSS;
            else if (key=="STOL")  LP>>_STOL;
            else if (key=="dTini") LP>>_dTini;
            else if (key=="mMin")  LP>>_mMin;
            else if (key=="mMax")  LP>>_mMax;
            else if (key=="mCoef") LP>>_mCoef;
            else if (key=="ZTOL")  LP>>_ZTOL;
            else if (key=="unbcor") LP>>_unbal_correction;
            else if (key=="checkConvergence") LP>>_check_convergence;
            else throw new Fatal(_("AutoME::AutoME: < %s > is invalid "), key.c_str());
        }
    }
    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:
    int  _maxSS;
    REAL _STOL;
    REAL _dTini;
    REAL _mMin;
    REAL _mMax;
    REAL _mCoef;
    REAL _ZTOL;
    bool _unbal_correction;
    bool _check_convergence;

    LinAlg::Vector<REAL> _dF_1;
    LinAlg::Vector<REAL> _dU_1;
    LinAlg::Vector<REAL> _dF_2;
    LinAlg::Vector<REAL> _dU_2;
    LinAlg::Vector<REAL> _dU_ME;
    LinAlg::Vector<REAL> _Err_U;
    LinAlg::Vector<REAL> _Err_F;
    LinAlg::Vector<REAL> _dU_unb;
}; // class AutoME




inline void AutoME::_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::Copy(DF_ext, dF_ext);
    LinAlg::Copy(DU_ext, dU_ext);

    // Allocate auxiliar vectors
    if (increment==0)
    {
        _dF_1 .Resize(n_tot_dof);
        _dU_1 .Resize(n_tot_dof);
        _dF_2 .Resize(n_tot_dof);
        _dU_2 .Resize(n_tot_dof);
        _dU_ME.Resize(n_tot_dof);
        _Err_U.Resize(n_tot_dof);
        _Err_F.Resize(n_tot_dof);
        //if (_unbal_correction)
            //_dU_unb.Resize(n_tot_dof);
    }
    
    // Initialize unbalanced forces array
    //if (_unbal_correction)
        //_dU_unb.SetValues(0.0);

    // Start substeps
    REAL  T = 0.0;
    REAL dT = _dTini;
    for (int k=0; k<_maxSS; ++k)
    {
        if (T>=1.0) return;

        // Calculate scaled increment vectors for this sub-step
        LinAlg::CopyScal(dT,dF_ext, _dF_1); // _dF_1 <- dT*dF_ext
        LinAlg::CopyScal(dT,dU_ext, _dU_1); // _dU_1 <- dT*dU_ext;
        LinAlg::CopyScal(dT,dF_ext, _dF_2); // _dF_2 <- dT*dF_ext
        LinAlg::CopyScal(dT,dU_ext, _dU_2); // _dU_2 <- dT*dU_ext;

        // Backup element state (needed when updating disp. state for a ME increment)
        _backup_nodes_and_elements_state();

        // Forward-Euler: Assemble K matrix and calculate _dU_1
        _assemb_K(K);
        _inv_K_times_X(K, false, _dF_1, _dU_1);

        // Add unbalanced displacements
        //if (_unbal_correction)
            //LinAlg::Axpy(1.0,_dU_unb, _dU_1); // _dU_1 <- _dU_1 + _dU_unb

        // Forward-Euler: update nodes and elements state
        _update_nodes_and_elements_state(_dU_1, dF_int);

        // Modified-Euler: Assemble K matrix and calculate dU_2
        _assemb_K(K);
        _inv_K_times_X(K, false, _dF_2, _dU_2);

        // Save the norm of essential and natural vectors
        REAL normU = _norm_essential_vector(); // (%)
        REAL normF = _norm_natural_vector();

        // Calculate local error
        if (normF<=_ZTOL) throw new Message(_("AutoME::_do_solve_for_an_increment: k=%d: normF=%e cannot be equal to zero (%e)"),k,normF,_ZTOL);
        LinAlg::AddScaled(0.5,_dU_2, -0.5,_dU_1, _Err_U); // Error on U (%)
        LinAlg::AddScaled(0.5,_dF_2, -0.5,_dF_1, _Err_F); // Error on F
        REAL Rerr_U = LinAlg::Norm(_Err_U)/(1.0+normU);   // (R)elative error on U
        REAL Rerr_F = LinAlg::Norm(_Err_F)/normF;         // (R)elative error on F
        REAL   Rerr = (Rerr_U>Rerr_F ? Rerr_U : Rerr_F);  // (R)elative error
        REAL      m = _mCoef*sqrt(_STOL/Rerr);

        // Restore nodes and element for initial disp. state given at the start of the increment
        _restore_nodes_and_elements_state();

        if (Rerr<=_STOL)
        {
            // Calculate Modified-Euler displacement increment vector
            LinAlg::AddScaled(0.5,_dU_1, 0.5,_dU_2, _dU_ME);

            // Update nodes and elements state for a Modified-Euler evaluation of displacements
            _update_nodes_and_elements_state(_dU_ME, dF_int);

            /*
            // Calculate residual (internal)
            LinAlg::AddScaled (1.0,_dF_2, -1.0,dF_int, Rinternal); // Rinternal <- dF_ext - dF_int
            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);
            ISD.ME_resid(LinAlg::Norm(Rinternal)/sqrt(denom));
            */

            // Calculate unbalanced displacements
            //if (_unbal_correction)
                //_inv_K_times_X(K, false, Rinternal, _dU_unb);

            // Next pseudo time
            T = T + dT;
            if (m>_mMax) m=_mMax;
        }
        else
            if (m<_mMin) m=_mMin;

        // Next substep size
        dT = m*dT;
        if (dT>1.0-T) dT=1.0-T;

        // IntSolverData
        ISD.ME_Rerr(Rerr).ME_T(T).ME_dT(dT).ME_m(m);
    }
    if (_check_convergence)
        throw new Fatal(_("AutoME::_do_solve_for_an_increment: did not converge for %d substeps"), _maxSS);

} // }}}




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

// Register an AutoME solver into SolverFactory array map
int AutoMERegister() // {{{
{
    SolverFactory["AutoME"] = AutoMEMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __AutoME_dummy_int = AutoMERegister();

}; // namespace FEM

#endif // MECHSYS_FEM_AUTOME_H

// vim:fdm=marker

---