cforwardeuler.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_CFORWARDEULER_H
#define MECHSYS_FEM_CFORWARDEULER_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 CForwardEuler: public CSolver
{
public:
    CForwardEuler(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(_("CForwardEuler::CForwardEuler: < %s > is invalid "), key.c_str());
        }
        if (_nSS<1) throw new Fatal(_("CForwardEuler::CForwardEuler: 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;
}; // class CForwardEuler




inline void CForwardEuler::_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)
    {
        // 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)

        // Update nodes and elements state
        _update_nodes_and_elements_state(dU_ext, dF_int, h);

        // Calculate residual (internal)
        LinAlg::Axpy      (+1.0,hKUn, dF_ext);                  // dF_ext    <- dF_ext + hKUn
        LinAlg::AddScaled (1.0,dF_ext, -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.FE_resid(LinAlg::Norm(Rinternal)/sqrt(denom));
    }

} // }}}




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

// Register an CForwardEuler solver into CSolverFactory array map
int CForwardEulerRegister() // {{{
{
    CSolverFactory["CForwardEuler"] = CForwardEulerMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __CForwardEuler_dummy_int = CForwardEulerRegister();

}; // namespace FEM

#endif // MECHSYS_FEM_CFORWARDEULER_H

// vim:fdm=marker

---