equipment.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_LABTESTSIM_EQUIPMENT_H
#define MECHSYS_LABTESTSIM_EQUIPMENT_H

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

#include <cmath>
#include <fstream>

#include "models/equilibmodel.h"
#include "tensors/functions.h"
#include "tensors/tensors.h"
#include "util/numstreams.h"
#include "util/exception.h"
#include "util/string.h"
#include "util/array.h"
#include "util/util.h" // for Sort

namespace LabTestSim
{

// Auxiliar structure
struct State
{
    REAL Sa;
    REAL Sr;
    REAL St;
    REAL Pp; // Pore-pressure
}; // struct State

template<typename T_Model>
class Equipment
{
public:

    // Constructor
    Equipment (T_Model * pModel, String const & strEquipment, bool StopOnFailure, bool SaveIntStates, bool Silent=false);

    // Destructor
    virtual ~Equipment() {}

    // Methods
    void                     Run       (int nDiv, Array<State> const & States);
    void                     WriteFile (String const & Filename, bool PrintPrincVals=false, bool PrintHead=true);
    int                      nStates   ()            const { return _SIG.size();   }
    Tensors::Tensor2 const & Sig       (int i_state) const { return _SIG[i_state]; }
    Tensors::Tensor2 const & Eps       (int i_state) const { return _EPS[i_state]; }

protected:
    // Data
    T_Model * _model;
    String    _equipment;
    bool      _stop_on_failure;
    bool      _save_int_states; // Save intermediary states?
    bool      _silent;          // Do not show warning messages

    // Methods
    virtual int  _exx_idx() const =0;
    virtual bool _do_update_one_increment(State const & DeltaState) =0;

private:
    // Data
    Array<Tensors::Tensor2> _SIG; // Stress tensor's        size == num states
    Array<Tensors::Tensor2> _EPS; // Strain tensor's        size == num states
    Array<Array<REAL>     > _ISV; // Internal State Values  size == num states
    Array<REAL            > _PPR; // Pore-pressures

}; // class Equipment




template<typename T_Model>
inline Equipment<T_Model>::Equipment(T_Model * pModel, String const & strEquipment, bool StopOnFailure, bool SaveIntStates, bool Silent) // {{{
    : _model           (pModel)       ,
      _equipment       (strEquipment) ,
      _stop_on_failure (StopOnFailure),
      _save_int_states (SaveIntStates),
      _silent          (Silent)
{
} // }}}

template<typename T_Model>
inline void Equipment<T_Model>::Run(int nDiv, Array<State> const & States ) // {{{
{
    // Number of states
    size_t n_states = States.size();

    // Re-allocate arrays
    if (_save_int_states)
    {
        _SIG.resize(1); // will use push_back
        _EPS.resize(1);
        _ISV.resize(1);
#ifdef USE_COUPLEDMODEL
        _PPR.resize(1); // pore-pressures
#endif
    }
    else
    {
        _SIG.resize(n_states);
        _EPS.resize(n_states);
        _ISV.resize(n_states);
#ifdef USE_COUPLEDMODEL
        _PPR.resize(n_states);
#endif
    }

    // Save initial state
    Array<REAL> tmp_isv; _model->InternalStateValues(tmp_isv);
    _SIG[0] = _model->Sig();
    _EPS[0] = _model->Eps();
    _ISV[0] = tmp_isv;
#ifdef USE_COUPLEDMODEL
    _PPR[0] = _model->Pp();
#endif

    // Loop along final states
    for (size_t i_state=1; i_state<n_states; ++i_state)
    {
        // Delta
        Tensors::Tensor2 const & sig = _model->Sig();
#ifdef USE_COUPLEDMODEL
        REAL ppr = _model->Pp();
#endif
        State delta_state;
        delta_state.Sa  = (States[i_state].Sa  - sig(0))/nDiv;
        delta_state.Sr  = (States[i_state].Sr  - sig(1))/nDiv;
        delta_state.St  = (States[i_state].St  - sig(2))/nDiv;
#ifdef USE_COUPLEDMODEL
        delta_state.Pp  = (States[i_state].Pp  - ppr   )/nDiv;
#endif

        // Loop along divisions
        for (int j=0; j<nDiv; ++j)
        {
            bool ok = _do_update_one_increment(delta_state);
            if (_save_int_states)
            {
                _model->InternalStateValues(tmp_isv);
                _SIG.push_back(_model->Sig());
                _EPS.push_back(_model->Eps());
                _ISV.push_back(tmp_isv);
#ifdef USE_COUPLEDMODEL
                _PPR.push_back(_model->Pp());
#endif
            }
            if ((!ok) && _stop_on_failure) return;
            if (!ok) break;
        }

        // Save state
        if (!_save_int_states)
        {
            _model->InternalStateValues(tmp_isv);
            _SIG[i_state] = _model->Sig();
            _EPS[i_state] = _model->Eps();
            _ISV[i_state] = tmp_isv;
#ifdef USE_COUPLEDMODEL
            _PPR[i_state] = _model->Pp();
#endif
        }
    }
} // }}}

template<typename T_Model>
inline void Equipment<T_Model>::WriteFile(String const & Filename, bool PrintPrincVals, bool PrintHead) // {{{
{

    using Util::_8s;

    std::ofstream of(Filename.GetSTL().c_str(), std::ios::out);
    if (!of.fail())
    {
        if (PrintHead)
        {
            // Heading
            of << _8s()<< "p"  << _8s()<< "q"  << _8s()<< "sin3th" << _8s()<< "Ev" << _8s()<< "Ed";
            of << _8s()<< "Sa" << _8s()<< "Sr" << _8s()<< "St"     << _8s()<< "Sat";
            of << _8s()<< "Ea" << _8s()<< "Er" << _8s()<< "Et"     << _8s()<< "Eat";

            // Internal state values names
            Array<String> names;  _model->InternalStateNames(names);
            for (size_t i=0; i<names.size(); ++i)
                of << _8s()<< names[i];

            // Principal values
            if (PrintPrincVals)
            {
                of << _8s()<< "S1" << _8s()<< "S2" << _8s()<< "S3";
                of << _8s()<< "E1" << _8s()<< "E2" << _8s()<< "E3";
            }

#ifdef USE_COUPLEDMODEL
            // Pore-pressure
            of << _8s()<< "Pp";
#endif

            // Next line
            of << std::endl;
        }

        REAL sq2 = sqrt(2.0);

        for (size_t i=0; i<_EPS.size(); ++i)
        {
            // Invariants
            REAL p,q,sin3th; Tensors::Tensor2 S; Tensors::Stress_p_q_S_sin3th(_SIG[i], p,q,S,sin3th);
            REAL Ev,Ed;                          Tensors::Strain_Ev_Ed       (_EPS[i], Ev,Ed);

            // Output values
            of << _8s()<< p                << _8s()<< q                << _8s()<< sin3th           << _8s()<< Ev*100.0 << _8s()<< Ed*100.0;
            of << _8s()<< _SIG[i](0)       << _8s()<< _SIG[i](1)       << _8s()<< _SIG[i](2)       << _8s()<< _SIG[i](5)/sq2;
            of << _8s()<< _EPS[i](0)*100.0 << _8s()<< _EPS[i](1)*100.0 << _8s()<< _EPS[i](2)*100.0 << _8s()<< _EPS[i](5)*100.0/sq2;

            // Output internal values
            for (size_t j=0; j<_ISV[i].size(); ++j)
                of << _8s()<< _ISV[i][j];

            // Principal values
            if (PrintPrincVals)
            {
                REAL S[3]; Tensors::Eigenvals(_SIG[i], S);  Util::Sort(S,3);
                REAL E[3]; Tensors::Eigenvals(_EPS[i], E);  Util::Sort(E,3);
                of << _8s()<< S[2] << _8s()<< S[1] << _8s()<< S[0];
                of << _8s()<< E[2] << _8s()<< E[1] << _8s()<< E[0];
            }

#ifdef USE_COUPLEDMODEL
            // Pore-pressure
            of << _8s()<< _PPR[i];
#endif

            // Next line
            of << std::endl;
        }

        of.close();
    }
    else
        throw new Fatal(_("Equipment::WriteFile: Could not open file < %s > for output"), Filename.c_str());

} // }}}

}; // namespace LabTestSim

#endif // MECHSYS_LABTESTSIM_EQUIPMENT_H

// vim:fdm=marker

---