cmodifiedeuler.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_CMODIFIEDEULER_H
#define MECHSYS_FEM_CMODIFIEDEULER_H

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

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

namespace FEM
{

class CModifiedEuler: public CSolver
{
public:
    CModifiedEuler(Array<String> const & CSolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output)
        : CSolver (ID, data, output),
          _nSS   (1)
    {
        for (size_t i=0; i<CSolverCtes.size(); ++i)
        {
            LineParser LP(CSolverCtes[i]);
            LP.ReplaceAllChars('=',' ');
            String key;     LP>>key;
            if (key=="nSS") LP>>_nSS;
            else throw new Fatal(_("CModifiedEuler::CModifiedEuler: < %s > is invalid "), key.c_str());
        }
        if (_nSS<1) throw new Fatal(_("CModifiedEuler::CModifiedEuler: 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
                                    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 _nSS;

    LinAlg::Vector<REAL> _dF_2;
    LinAlg::Vector<REAL> _dU_2;
    LinAlg::Vector<REAL> _dU_ME;
}; // class CModifiedEuler




inline void CModifiedEuler::_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
{
    // 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
    REAL h = dTime/_nSS;

    // Start
    for (int i=0; i<_nSS; ++i)
    {
        // Resize vectors
        if (increment==0)
        {
            int n_tot_dof = dF_ext.Size();
            _dF_2 .Resize(n_tot_dof);
            _dU_2 .Resize(n_tot_dof);
            _dU_ME.Resize(n_tot_dof);
        }

        // Setup temporary increment vectors
        LinAlg::Copy(dF_ext, _dF_2);
        LinAlg::Copy(dU_ext, _dU_2);

        // 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_ext
        _assemb_G      (G, hKUn, h);           // G <- C + h*alpha*K
        LinAlg::Axpy   (-1.0,hKUn, dF_ext);        // dF_ext <- dF_ext - hKUn
        _inv_K_times_X (G, false, dF_ext, dU_ext); // dU_ext <- inv(G)*(dF_ext - hKUn)
        
        // Forward-Euler: update nodes and elements state
        _update_nodes_and_elements_state(dU_ext, 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)
        
        // Calculate Modified-Euler displacement increment vector
        LinAlg::AddScaled(0.5,dU_ext, 0.5,_dU_2, _dU_ME); // _dU_ME <- 0.5*dU_ext + 0.5*_dU_2

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

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

} // }}}




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

// Register an CModifiedEuler solver into CSolverFactory array map
int CModifiedEulerRegister() // {{{
{
    CSolverFactory["CModifiedEuler"] = CModifiedEulerMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __CModifiedEuler_dummy_int = CModifiedEulerRegister();

}; // namespace FEM

#endif // MECHSYS_FEM_CMODIFIEDEULER_H

// vim:fdm=marker

---