solver.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_SOLVER_H
#define MECHSYS_FEM_SOLVER_H

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

#include "fem/solver/intsolverdata.h"
#include "util/string.h"
#include "fem/data.h"
#include "fem/output.h"
#include "fem/node.h"
#include "fem/ele/element.h"
#include "linalg/vector.h"
#include "linalg/matrix.h"
#include "linalg/lawrap.h"

namespace FEM
{


class Solver
{
public:
    // Constructor
    Solver(FEM::InputData const & ID, FEM::Data & data, FEM::Output & output);
    // Destructor
    virtual ~Solver() {}
    // Methods
    void Solve(int iStage);
protected:
    // Data
    FEM::InputData            const & _idat;
    FEM::Data                       & _data;    // Nodes and elements data
    FEM::Output                     & _output;  // Nodes and elements data
    Array<FEM::Node::DOFVarsStruct*>  _a_udofs; // Array of ptrs to unknown DOFs.    size = total number of unknown dofs of a current stage
    Array<FEM::Node::DOFVarsStruct*>  _a_pdofs; // Array of ptrs to prescribed DOFs. size = total number of prescribed dofs of a current stage
    LinAlg::Vector<REAL>              _DF_ext;  // External increment of natural dofs.    size = total number of DOFs
    LinAlg::Vector<REAL>              _DU_ext;  // External increment of essential dofs.  size = total number of DOFs

    // Local
    void _realloc_and_setup_dof_vectors();
    void _assemb_natural_stage_inc();         // Allocate and assemble global (stage) increment of external force vector for the current stage
    void _assemb_essential_stage_inc();       // Allocate and assemble global (stage) increment displacements vector for the current stage
    void _assemb_K(LinAlg::Matrix<REAL> & K); // Assemble matrix K for the actual state of elements/nodes inside _data

    void _inv_K_times_X(LinAlg::Matrix<REAL> & K           ,  // In/Out: (if DoPreserveK==true => just In)  Out => unknown state (values)
                        bool                   DoPreserveK ,  // In:
                        LinAlg::Vector<REAL> & X           ,  // In/Out:
                        LinAlg::Vector<REAL> & Y           ); // In/Out:   Y <- inv(K)*X   (K*Y=X)
    int _n_tot_dof() const { return _a_udofs.size() + _a_pdofs.size(); }

    void _backup_nodes_and_elements_state();
    void _update_nodes_and_elements_state(LinAlg::Vector<REAL> const & dU, LinAlg::Vector<REAL> & dF_int);
    void _restore_nodes_and_elements_state();

    REAL _norm_essential_vector();
    REAL _norm_natural_vector();

    // To be overriden by derived methods
    virtual 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      ) =0; // Out
private:
    LinAlg::Vector<REAL> _U_bkp; // A place to store a temporary state of displacements (backup)
    LinAlg::Vector<REAL> _F_bkp; // A place to store a temporary state of internal forces (backup)
    // Local functions
    void _do_scatter(LinAlg::Matrix<REAL> const & K,   // {{{
                     LinAlg::Vector<REAL> const & X,
                     LinAlg::Vector<REAL> const & Y,
                     LinAlg::Matrix<REAL>       & K11, // Unknown_size x Unknown_size
                     LinAlg::Matrix<REAL>       & K12, // Unknown_size x Prescribed_size
                     LinAlg::Matrix<REAL>       & K21, // Prescribed_size x Unknown_size
                     LinAlg::Matrix<REAL>       & K22, // Prescribed_size x Prescribed_size
                     LinAlg::Vector<REAL>       & Y2,  // Prescribed_size
                     LinAlg::Vector<REAL>       & X1); // Unknown_size }}}

    void _do_gather(LinAlg::Vector<REAL> const & Y1, // Unknown_size {{{
                    LinAlg::Vector<REAL> const & Y2, // Prescribed_size
                    LinAlg::Vector<REAL> const & X1, // Unknown_size
                    LinAlg::Vector<REAL> const & X2, // Prescribed_size
                    LinAlg::Vector<REAL>       & Y ,
                    LinAlg::Vector<REAL>       & X); // }}}
}; // class Solver




inline Solver::Solver(FEM::InputData const & ID, FEM::Data & data, FEM::Output & output) // {{{
    : _idat   (ID),
      _data   (data),
      _output (output)
{
} // }}}

inline void Solver::Solve(int iStage) // {{{
{
    // - Realloc DOF vectors: _a_pdofs and _a_pdofs
    // - Give an equation ID number to all DOFs inside all nodes
    _realloc_and_setup_dof_vectors(); // EqID. OBS.: _n_tot_dof() will work only after this call

    // Allocate auxiliar vectors for backup displacements and internal forces state
    _U_bkp.Resize(_n_tot_dof());
    _F_bkp.Resize(_n_tot_dof());

    // Allocate and assemble global (stage) increment of external force vector for the current stage
    _assemb_natural_stage_inc(); // assemb DF_ext

    // Allocate and assemble global (stage) increment displacements vector for the current stage
    _assemb_essential_stage_inc(); // assemb DU_ext

    // ------------------------------------------------------------------------------------ do solve ---

    // Total number of DOF == number of equations
    int n_tot_dof = _n_tot_dof();

    // Number of divisions for each increment
    int num_div = _idat.nDIV[iStage];

    // Scale increments according to a step controller
    REAL lam=1.0/num_div;        // Step controller
    LinAlg::Scal(lam, _DF_ext);  // dF_ext <- lam*_DF_ext
    LinAlg::Scal(lam, _DU_ext);  // dF_ext <- lam*_DF_ext

    // Allocate K stifness (dense) matrix
    LinAlg::Matrix<REAL> K(n_tot_dof, n_tot_dof);

    // Allocate a vector of residual forces = dF_int - dF_ext
    LinAlg::Vector<REAL> dF_int   (n_tot_dof);
    LinAlg::Vector<REAL> Rinternal(n_tot_dof);

    // Output initial state
    if (iStage==0) _output.OutputIncrement(iStage, -1);

    // Loop over increments (0<increment<NumDiv)
    for (int increment=0; increment<num_div; ++increment)
    {
        // Internal solver data
        IntSolverData isd;

        // Solve for increment
        _do_solve_for_an_increment(increment, _DF_ext, _DU_ext, K, dF_int, Rinternal, isd);

        // Output results for increment
        _output.OutputIncrement(iStage, increment, isd);
    }

    // Output end of Stage
    _output.OutputStage(iStage);
    // ----------------------------------------------------------------------------------- end solve ---
} // }}}

inline void Solver::_realloc_and_setup_dof_vectors() // {{{
{
    // ID of equation
    int eq_id=0;

    // Deallocate vectors
    _a_udofs.resize(0);
    _a_pdofs.resize(0);

    // Loop along all nodes
    for (size_t i=0; i<_data.Nodes.size(); ++i)
    {
        // Declare a reference to DOFVars array inside Node [i]
        Array<FEM::Node::DOFVarsStruct> & dofs = _data.Nodes[i].DOFVars();

        // Loop along DOFs inside Node [i]
        for (size_t j=0; j<dofs.size(); ++j)
        {
            // Fill DOF vectors
            if (dofs[j].IsEssenPresc) _a_pdofs.push_back(&dofs[j]); // prescribed
            else                      _a_udofs.push_back(&dofs[j]); // unknown
            
            // Assing an equation ID to DOF [j]
            dofs[j].EqID = eq_id;
            eq_id++;
        }
    }
} // }}}

inline void Solver::_assemb_natural_stage_inc() // Allocate and assemble global (stage) increment of external force vector for the current stage {{{
{
    // Total number of DOF == number of equations
    int n_tot_dof = _n_tot_dof();

    // Resize vector DF_ext
    _DF_ext.Resize(n_tot_dof);

    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
        _DF_ext(_a_udofs[iu]->EqID) = _a_udofs[iu]->NaturalBry;

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
        _DF_ext(_a_pdofs[ip]->EqID) = _a_pdofs[ip]->NaturalBry;
} // }}}

inline void Solver::_assemb_essential_stage_inc() // Allocate and assemble global (stage) increment displacements vector for the current stage {{{
{
    // Total number of DOF == number of equations
    int n_tot_dof = _n_tot_dof();

    // Resize vector DU_ext
    _DU_ext.Resize(n_tot_dof);

    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
        _DU_ext(_a_udofs[iu]->EqID) = _a_udofs[iu]->EssentialBry;

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
        _DU_ext(_a_pdofs[ip]->EqID) = _a_pdofs[ip]->EssentialBry;
} // }}}

inline void Solver::_assemb_K(LinAlg::Matrix<REAL> & K) // Assemble matrix K for the actual state of elements/nodes inside _data {{{
{
    // Clear K (global) stifness matrix
    K.SetValues(0.0);

    // Loop along elements
    for (size_t i_ele=0; i_ele<_data.Elements.size(); ++i_ele)
    {
        // Declare a pointer to Element [i_ele]
        FEM::EquilibElem const * const E = dynamic_cast<FEM::EquilibElem *>(_data.Elements[i_ele]);
        if (!E) throw new Fatal(_("FEM::Solver: This solver only works for equilibrium type elements"));
    
        // Allocate an array of position of DOF components inside Ke of Element [i_ele] {{{
        /* Ex.:
         *               a_map_Ke_to_K_pos[0] = EqID   =>   K[EqID,EqID] = Ke[0,0]
         *               a_map_Ke_to_K_pos[1] = EqID   =>   K[EqID,EqID] = Ke[1,1]
         *                                      ................
         *     a_map_Ke_to_K_pos[Ke.Rows()-1] = EqID   =>   K[EqID,EqID] = Ke[Ke.Rows()-1, Ke.Rows()-1]
         *
         * Ex.: 
         *       Node.NGL=2, Ele.NNODE=2 => Ke.Rows() = 4
         *
         *             Ke                K(global)_pos
         *   _                     _   
         *  |  Ke11 Ke12 Ke13 Ke14  | => Ele.Node[0].DOF[0].EqID
         *  |  Ke21 Ke22 Ke23 Ke24  | => Ele.Node[0].DOF[1].EqID
         *  |  Ke31 Ke32 Ke33 Ke34  | => Ele.Node[1].DOF[0].EqID
         *  |_ Ke41 Ke42 Ke43 Ke44 _| => Ele.Node[1].DOF[1].EqID
         }}} */
            
        if (E->IsActive())
        {
            // Calculate the stiffness matrix of Element [i_ele]
            LinAlg::Matrix<REAL> Ke;
            Array<size_t> a_map_Ke_to_K_pos; // size=Ke.Rows()=Ke.Cols()
            E->Stiffness(Ke, a_map_Ke_to_K_pos);


            // Assemble local (Ke) DOFs into global (K) matrix
            // Loop along Ke (local) matrix Ke.rows and Ke.cols
            for (int i_row_Ke=0; i_row_Ke<Ke.Rows(); ++i_row_Ke)
            for (int j_col_Ke=0; j_col_Ke<Ke.Rows(); ++j_col_Ke)
                K(a_map_Ke_to_K_pos[i_row_Ke], a_map_Ke_to_K_pos[j_col_Ke]) += Ke(i_row_Ke, j_col_Ke);
        }
    }
} // }}}

inline void Solver::_inv_K_times_X(LinAlg::Matrix<REAL> & K           , // In/Out: (if DoPreserveK==true => just In)  Out => unknown state (values) {{{
                                   bool                   DoPreserveK , // In:
                                   LinAlg::Vector<REAL> & X           , // In/Out:
                                   LinAlg::Vector<REAL> & Y           ) // In/Out:
{
    // Modify K for essential boundary conditions
    //if (_do_split)
    //{
        /*   _               _  
         *  |  [K11] | [K12]  | / {Y1}=? \   / {X1}   \
         *  |  - - - - - - -  | | - - -  | = | - - -  |
         *  |_ [K21] | [K22] _| \ {Y2}   /   \ {X2}=? /
         *
         *  {Y2} and {X1} = prescribed
         *
         *  K11*Y1 + K12*Y2 = X1
         *  K21*Y1 + K22*Y2 = X2
         *
         *  Y1 = inv(K11)*(X1 - K12*Y2)
         *  X2 = K21*Y1 + K22*Y2
         */

        // Allocate temporary matrices and vectors
        LinAlg::Matrix<REAL> K11,K12;
        LinAlg::Matrix<REAL> K21,K22;
        LinAlg::Vector<REAL> Y2;
        LinAlg::Vector<REAL> X1;

        // Scatter global K, X and Y values into smaller groups
        _do_scatter(K,X,Y, K11,K12,K21,K22, Y2, X1);

        // Solve for {Y1} and {X2}
        LinAlg::Vector<REAL> TMP(X1);                // TMP <- X1
        LinAlg::Gemv(-1.0,K12,Y2,1.0,TMP);           // TMP <- -1.0*K12*Y2 + 1.0*X1

        // Call a linear solver for: TMP <- inv(K11)*TMP (K11 is lost)
        if (_idat.linSOLVER==InputData::lsLAPACK)
        {
#if defined (USE_LAPACK) || defined (USE_GOTOBLAS)
            LinAlg::Gesv(K11, TMP); // TMP <- inv(K11)*TMP (K11 is lost)
#else
            throw new Fatal(_("Solver::_inv_K_times_X: LAPACK (or GOTOBLAS) is not available (USE_LAPACK=false, or USE_GOTOBLAS=false)"));
#endif
        }
        else if (_idat.linSOLVER==InputData::lsUMFPACK)
        {
            throw new Fatal(_("Solver::_inv_K_times_X: UMFPACK is not working yet"));
        }
        else if (_idat.linSOLVER==InputData::lsSUPERLU)
        {
            throw new Fatal(_("Solver::_inv_K_times_X: SUPERLU is not working yet"));
        }
        else throw new Fatal(_("Solver::_inv_K_times_X: linSOLVER is invalid"));

        LinAlg::Vector<REAL> & Y1 = TMP;             //  Y1 <- TMP (reference)
        LinAlg::Vector<REAL>   X2(Y2.Size());        // allocate X2
        LinAlg::Gemvpmv(1.0,K21,Y1, 1.0,K22,Y2, X2); //  X2 <- 1*K21*Y1 + 1*K22*Y2

        // Gather Y1, Y2, X1 and X2 into Y and X
        _do_gather(Y1,Y2, X1,X2, Y,X);
    //}
    //else if (_do_penalty)
    //{
    // Copy X to Res; neccessary for Gesv
    //   LinAlg::Copy(X, Y); // Y <- X
    //   K(ipresc,jpresc) <- Vpres*BIGNUM
    //   Gesv(K,Y) Y <- inv(K)*Y
    //}
} // }}}

inline void Solver::_do_scatter(LinAlg::Matrix<REAL> const & K,   // {{{
                                LinAlg::Vector<REAL> const & X,
                                LinAlg::Vector<REAL> const & Y,
                                LinAlg::Matrix<REAL>       & K11, // Unknown_size x Unknown_size
                                LinAlg::Matrix<REAL>       & K12, // Unknown_size x Prescribed_size
                                LinAlg::Matrix<REAL>       & K21, // Prescribed_size x Unknown_size
                                LinAlg::Matrix<REAL>       & K22, // Prescribed_size x Prescribed_size
                                LinAlg::Vector<REAL>       & Y2,  // Prescribed_size
                                LinAlg::Vector<REAL>       & X1)  // Unknown_size
{
    // Unknown_size = us
    // Prescribed_size = ps
    int us = _a_udofs.size(); // unknown
    int ps = _a_pdofs.size(); // prescribed

    // Allocate matrices/vectors
    K11.Resize(us,us);
    K12.Resize(us,ps);
    K21.Resize(ps,us);
    K22.Resize(ps,ps);
    Y2.Resize(ps);
    X1.Resize(us);
    // Loop along Prescribed_size
    for (int ip=0; ip<ps; ++ip)

    // Loop along Unknown_size
    for (int iu=0; iu<us; ++iu)
    {
        // K11: Loop along Unknown_size
        for (int ju=0; ju<us; ++ju)
            K11(iu,ju) = K(_a_udofs[iu]->EqID, _a_udofs[ju]->EqID);

        // K12: Loop along Prescribed_size
        for (int jp=0; jp<ps; ++jp)
            K12(iu,jp) = K(_a_udofs[iu]->EqID, _a_pdofs[jp]->EqID);

        // X1
        X1(iu) = X(_a_udofs[iu]->EqID);
    }

    // Loop along Prescribed_size
    for (int ip=0; ip<ps; ++ip)
    {
        // K21: Loop along Unknown_size
        for (int ju=0; ju<us; ++ju)
            K21(ip,ju) = K(_a_pdofs[ip]->EqID, _a_udofs[ju]->EqID);

        // K22: Loop along Prescribed_size
        for (int jp=0; jp<ps; ++jp)
            K22(ip,jp) = K(_a_pdofs[ip]->EqID, _a_pdofs[jp]->EqID);

        // Y2
        Y2(ip) = Y(_a_pdofs[ip]->EqID);
    }
} // }}}

inline void Solver::_do_gather(LinAlg::Vector<REAL> const & Y1, // Unknown_size {{{
                               LinAlg::Vector<REAL> const & Y2, // Prescribed_size
                               LinAlg::Vector<REAL> const & X1, // Unknown_size
                               LinAlg::Vector<REAL> const & X2, // Prescribed_size
                               LinAlg::Vector<REAL>       & Y ,
                               LinAlg::Vector<REAL>       & X) 
{
    // Unknown_size = us
    // Prescribed_size = ps
    int us = _a_udofs.size(); // unknown
    int ps = _a_pdofs.size(); // prescribed

    // Loop along Unknown_size
    for (int iu=0; iu<us; ++iu)
    {
        // Y1
        Y(_a_udofs[iu]->EqID) = Y1(iu);

        // X1
        X(_a_udofs[iu]->EqID) = X1(iu);
    }

    // Loop along Prescribed_size
    for (int ip=0; ip<ps; ++ip)
    {
        // Y2
        Y(_a_pdofs[ip]->EqID) = Y2(ip);

        // X2
        X(_a_pdofs[ip]->EqID) = X2(ip);
    }
} // }}}

inline void Solver::_backup_nodes_and_elements_state() // {{{
{
    // ### Backup all nodes ############################################################################
    
    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
    {
         _U_bkp(_a_udofs[iu]->EqID) = _a_udofs[iu]->EssentialVal;
         _F_bkp(_a_udofs[iu]->EqID) = _a_udofs[iu]->NaturalVal;
    }

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
    {   
         _U_bkp(_a_pdofs[ip]->EqID) = _a_pdofs[ip]->EssentialVal;
         _F_bkp(_a_pdofs[ip]->EqID) = _a_pdofs[ip]->NaturalVal;
    }

    // ### Backup all elements #########################################################################
    
    // Loop along elements
    for (size_t i=0; i<_data.Elements.size(); ++i)
        _data.Elements[i]->BackupState();

    // #################################################################################################
} // }}}

inline void Solver::_update_nodes_and_elements_state(LinAlg::Vector<REAL> const & dU, LinAlg::Vector<REAL> & dF_int) // {{{
{
    // ### Update all nodes #############################################################################
    
    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
        _a_udofs[iu]->EssentialVal += dU(_a_udofs[iu]->EqID);

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
        _a_pdofs[ip]->EssentialVal += dU(_a_pdofs[ip]->EqID);

    // ### Update all elements ##########################################################################
    
    // --- Setup dU values for all nodes and Zeroes _IncNatural values ----------------------------------
    
    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
    {
        _a_udofs[iu]->_IncEssenVal = dU(_a_udofs[iu]->EqID);
        _a_udofs[iu]->_IncNaturVal = 0.0; // Need to be set zero because will be output
    }

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
    {
        _a_pdofs[ip]->_IncEssenVal = dU(_a_pdofs[ip]->EqID);
        _a_pdofs[ip]->_IncNaturVal = 0.0; // Need to be set zero because will be output
    }

    // --- Update all elements --------------------------------------------------------------------------
    
    // Loop along elements
    for (size_t i=0; i<_data.Elements.size(); ++i)
        _data.Elements[i]->UpdateState(); // Read nodal values _IncEssential and write to _IncNatural

    // --- Read dF values (just calculated) from all nodes ----------------------------------------------
    
    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
        dF_int(_a_udofs[iu]->EqID) = _a_udofs[iu]->_IncNaturVal;

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
        dF_int(_a_pdofs[ip]->EqID) = _a_pdofs[ip]->_IncNaturVal;

    // ##################################################################################################
} // }}}

inline void Solver::_restore_nodes_and_elements_state() // {{{
{
    // ### Restore all nodes ###########################################################################
    
    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
    {
        _a_udofs[iu]->EssentialVal = _U_bkp(_a_udofs[iu]->EqID);
          _a_udofs[iu]->NaturalVal = _F_bkp(_a_udofs[iu]->EqID);
    }

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
    {
        _a_pdofs[ip]->EssentialVal = _U_bkp(_a_pdofs[ip]->EqID);
          _a_pdofs[ip]->NaturalVal = _F_bkp(_a_pdofs[ip]->EqID);
    }

    // ### Restore all elements ########################################################################
    
    // Loop along elements
    for (size_t i=0; i<_data.Elements.size(); ++i)
        _data.Elements[i]->RestoreState();

    // #################################################################################################
} // }}}

inline REAL Solver::_norm_essential_vector() // {{{
{
    REAL tmp=0.0;
    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
        tmp += pow(_a_udofs[iu]->EssentialVal,2.0);

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
        tmp += pow(_a_pdofs[ip]->EssentialVal,2.0);

    return sqrt(tmp);
} // }}}

inline REAL Solver::_norm_natural_vector() // {{{
{
    REAL tmp=0.0;
    // Loop along Unknown_size
    for (size_t iu=0; iu<_a_udofs.size(); ++iu)
        tmp += pow(_a_udofs[iu]->NaturalVal,2.0);

    // Loop along Prescribed_size
    for (size_t ip=0; ip<_a_pdofs.size(); ++ip)
        tmp += pow(_a_pdofs[ip]->NaturalVal,2.0);

    return sqrt(tmp);
} // }}}



// Define a pointer to a function that makes (allocate) a new solver
typedef Solver * (*SolverMakerPtr)(Array<String> const & SolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output);

// Typdef of the array map that contains all the pointers to the functions that makes solvers
typedef std::map<String, SolverMakerPtr, std::less<String> > SolverFactory_t;

// Instantiate the array map that contains all the pointers to the functions that makes solvers
SolverFactory_t SolverFactory;

// Allocate a new solver according to a string giving the name of the solver
Solver * AllocSolver(String const & SolverName, Array<String> const & SolverCtes, FEM::InputData const & ID, FEM::Data & data, FEM::Output & output) // {{{
{
    SolverMakerPtr ptr=NULL;
    ptr = SolverFactory[SolverName];
    if (ptr==NULL)
        throw new Fatal(_("FEM::AllocSolver: There is no < %s > solver implemented in this library"), SolverName.c_str());
    return (*ptr)(SolverCtes, ID, data, output);
} // }}}

}; // namespace FEM

#endif // MECHSYS_FEM_SOLVER_H

// vim:fdm=marker

---