modifiedeuler.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_MODIFIEDEULER_H
#define MECHSYS_FEM_MODIFIEDEULER_H

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

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

namespace FEM
{

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

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



inline void ModifiedEuler::_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::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

    // 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 K matrix and calculate dU_ext
        _assemb_K      (K);
        _inv_K_times_X (K, false, dF_ext, dU_ext); // dU_ext <- inv(K)*dF_ext = dU_FE


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

        // Modified-Euler: Assemble K matrix and calculate dU_2
        _assemb_K      (K);
        _inv_K_times_X (K, false, _dF_2, _dU_2); // _dU_ori <- inv(K)*_dF_ori = dU_intermediate

        // 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_ori

        // 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);
    }

} // }}}




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

// Register an ModifiedEuler solver into SolverFactory array map
int ModifiedEulerRegister() // {{{
{
    SolverFactory["ModifiedEuler"] = ModifiedEulerMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __ModifiedEuler_dummy_int = ModifiedEulerRegister();

}; // namespace FEM

#endif // MECHSYS_FEM_MODIFIEDEULER_H

// vim:fdm=marker

---