cautome.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_CAUTOME_H
#define MECHSYS_FEM_CAUTOME_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/csolver.h"

namespace FEM
{

class CAutoME: public CSolver
{
public:
    CAutoME(Array<String> const & CSolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output)
        : CSolver (ID, data, output),
          _maxSS (200)   ,
          _STOL  (1.0e-2),
          _dTini (0.001) ,
          _mMin  (0.01)  ,
          _mMax  (10)    ,
          _mCoef (0.7)   ,
          _ZTOL  (1e-5)  ,
          _check_convergence (true)
    {
        for (size_t i=0; i<CSolverCtes.size(); ++i)
        {
            LineParser LP(CSolverCtes[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=="checkConvergence") LP>>_check_convergence;
            else throw new Fatal(_("CAutoME::CAutoME: < %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
                                    REAL                 const   dTime    ,  // In
                                    LinAlg::Matrix<REAL>       & G        ,  // (no needed outside)
                                    LinAlg::Vector<REAL>       & hKUn     ,  // (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 _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;
}; // class CAutoME




inline void CAutoME::_do_solve_for_an_increment(int                  const   increment,  // In // {{{
                                                LinAlg::Vector<REAL> const & DF_ext   ,  // In (may be changed)
                                                LinAlg::Vector<REAL> const & DU_ext   ,  // In (may be changed)
                                                REAL                 const   dTime    ,  // In
                                                LinAlg::Matrix<REAL>       & G        ,  // (no needed outside)
                                                LinAlg::Vector<REAL>       & hKUn     ,  // (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);
    }

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

        // Sub-divide timestep
        REAL h = dT*dTime;

        // 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 G matrix and calculate _dU_1
        _assemb_G      (G, hKUn, h);             // G <- C + h*alpha*K
        LinAlg::Axpy   (-1.0,hKUn, _dF_1);       // _dF_1 <- _dF_1 - hKUn
        _inv_K_times_X (G, false, _dF_1, _dU_1); // _dU_1 <- inv(G)*(dF_ext - hKUn)
    
        // Forward-Euler: update nodes and elements state
        _update_nodes_and_elements_state(_dU_1, dF_int, h);

        // Modified-Euler: Assemble G matrix and calculate dU_2
        _assemb_G      (G, hKUn, h);             // G <- C + h*alpha*K
        LinAlg::Axpy   (-1.0,hKUn, _dF_2);       // _dF_2 <- _dF_2 - hKUn
        _inv_K_times_X (G, false, _dF_2, _dU_2); // _dU_2 <- inv(G)*(dF_ext - hKun)
    
        // 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(_("CAutoME::_do_solve_for_an_increment: k=%d: normF=%e cannot be equal to ZTOL (%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, h);

            // 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(_("CAutoME::_do_solve_for_an_increment: did not converge for %d substeps"), _maxSS);

} // }}}




// Allocate a new CAutoME solver
CSolver * CAutoMEMaker(Array<String> const & CSolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output) // {{{
{
    return new CAutoME(CSolverCtes, ID, data, output);
} // }}}

// Register an CAutoME solver into CSolverFactory array map
int CAutoMERegister() // {{{
{
    CSolverFactory["CAutoME"] = CAutoMEMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __CAutoME_dummy_int = CAutoMERegister();

}; // namespace FEM

#endif // MECHSYS_FEM_CAUTOME_H

// vim:fdm=marker

---