forwardeuler.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_FORWARDEULER_H
#define MECHSYS_FEM_FORWARDEULER_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 ForwardEuler: public Solver
{
public:
    ForwardEuler(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(_("ForwardEuler::ForwardEuler: < %s > is invalid "), key.c_str());
        }
        if (_nSS<1) throw new Fatal(_("ForwardEuler::ForwardEuler: 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; // number of substeps
}; // class ForwardEuler




inline void ForwardEuler::_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)
    {
        // Assemble K matrix and calculate dU_ext
        _assemb_K      (K);
        _inv_K_times_X (K, false, dF_ext, dU_ext); // In=R,U; Out=dU  =>  dU=inv[K(U)]*R

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

        // Calculate residual (internal)
        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 ForwardEuler solver
Solver * ForwardEulerMaker(Array<String> const & SolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output) // {{{
{
    return new ForwardEuler(SolverCtes, ID, data, output);
} // }}}

// Register an ForwardEuler solver into SolverFactory array map
int ForwardEulerRegister() // {{{
{
    SolverFactory["ForwardEuler"] = ForwardEulerMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __ForwardEuler_dummy_int = ForwardEulerRegister();

}; // namespace FEM

#endif // MECHSYS_FEM_FORWARDEULER_H

// vim:fdm=marker

---