tbarx.cpp

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

#include <iostream>

#include "util/array.h"
#include "util/exception.h"
#include "models/barcelonax.h"
#include "tensors/tensors.h"
#include "tensors/functions.h"
#include "fem/solver/intsolverdata.h"

using Tensors::Tensor2;
using Tensors::Tensor4;
using std::cout;
using std::endl;

int main(int argc, char **argv) try
{
    REAL OCR    = 1.0;
    REAL OSI    = 10.0;
    REAL Pp_ini = -10.0;
    REAL DSigX  = 10;
    REAL DSigY  = 10;
    REAL DSigZ  = 10;
    REAL DPp    = -20;
    int  FEndiv = 1;
    if (argc>=2) OCR    = atof(argv[1]);
    if (argc>=3) OSI    = atof(argv[2]);
    if (argc>=4) Pp_ini = atof(argv[3]);
    if (argc>=5) DSigX  = atof(argv[4]);
    if (argc>=6) DSigY  = atof(argv[5]);
    if (argc>=7) DSigZ  = atof(argv[6]);
    if (argc>=8) DPp    = atof(argv[7]);
    if (argc>=9) FEndiv = atoi(argv[8]);

    cout << "OCR    = " << OCR    << endl;
    cout << "OSI    = " << OSI    << endl;
    cout << "Pp_ini = " << Pp_ini << endl;
    cout << "DSigX  = " << DSigX  << endl;
    cout << "DSigY  = " << DSigY  << endl;
    cout << "DSigZ  = " << DSigZ  << endl;
    cout << "DPp    = " << DPp    << endl;
    cout << "FEndiv = " << FEndiv << endl;
    cout << endl << endl;

    // Parameters
    Array<REAL> prms; prms.resize(18);
    prms[0]  = 0.25  ; // lam0
    prms[1]  = 0.05  ; // kap
    prms[2]  = 20.0  ; // phics (deg)
    prms[3]  = 10000 ; // G (kPa)
    prms[4]  = 0.9   ; // r
    prms[5]  = 0.1   ; // bet
    prms[6]  = 0.5   ; // k
    prms[7]  = 0.002 ; // lams
    prms[8]  = 0.001 ; // kaps
    prms[9]  = 101.3 ; // patm (kPa)
    prms[10] = 1     ; // pref
    prms[11] = 10    ; // B
    // Hydraulic
    prms[12] = 0.1;    // a
    prms[13] = 0.5;    // b
    prms[14] = 1.0;    // c
    prms[15] = 2.0;    // m
    prms[16] = 10e-8;  // ks
    prms[17] = 10;     // gamaW (kN/m³)

    // Initial Data
    Array<REAL> inidata; inidata.resize(7); //{ SigX SigY SigZ vini OCR OSI Pp }
    inidata[0] = 24;     // =z0 (kPa)
    inidata[1] = 24;     // =z0 (kPa)
    inidata[2] = 24;     // =z0 (kPa)
    inidata[3] = 1.8;    // v_ini
    inidata[4] = OCR;    // OCR
    inidata[5] = OSI;    // OSI (kPa)
    inidata[6] = Pp_ini; // Pp (kPa);  ==> Suc = 0-(-10) = +10 (kPa)

    // Barcelona Model
    BarcelonaX barx(prms, inidata);

    // Set Integration Schemes Constants Structure
    IntegSchemesCtes isc;
    isc.Type     (IntegSchemesCtes::FE);
    isc.FE_ndiv  (FEndiv);
    isc.ME_maxSS (20);
    isc.ME_STOL  (1e-1);
    isc.ME_dTini (0.1);
    isc.EdMax    (20);
    CoupledModel::SetIntegSchemesCtes(isc);

    // TgStiffness
    Tensor4 Dep;
    Tensor2 dep;
    barx.TgStiffness(Dep,dep);
    cout << "Dep =\n" << Dep << endl;
    cout << "dep ="   << dep << endl;
    cout << endl << endl;

    // Actualize
    Tensor2 DSig; DSig=DSigX,DSigY,DSigZ, 0,0,0;
    Tensor2 DEps;
    REAL    DnSr;
    cout << ".... Actualize ...\n";
    barx.BackupState();
    barx.Actualize(DSig,DPp, DEps,DnSr);
    barx.RestoreState();

    // StressUpdate
    Tensor2 dsig;
    REAL    dnsr;
    cout << ".... StressUpdate ...\n";
    barx.StressUpdate(DEps,DPp, dsig,dnsr);

    // Output
    cout << "DSig     = " << DSig*1     << endl;
    cout << "DEps (%) = " << DEps*100.0 << endl;
    cout << "dsig     = " << dsig*1     << endl;
    cout << "DnSr     = " << DnSr       << endl;
    cout << "dnsr     = " << dnsr       << endl;

    return 0;
}
catch (Exception * e) // {{{
{
    e->Cout();
    if (e->IsFatal()) exit(1);
    delete e;
}
catch (char const * s)
{
    cout << " Fatal:" << s << "\n";
    exit(1);
} // }}}

// vim:fdm=marker

---