/****************************************************************************
*
*                                 M U E S L I   v 1.8
*
*
*     Copyright 2020 IMDEA Materials Institute, Getafe, Madrid, Spain
*     Contact: muesli.materials@imdea.org
*     Author: Ignacio Romero (ignacio.romero@imdea.org)
*
*     This file is part of MUESLI.
*
*     MUESLI 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 3 of the License, or
*     (at your option) any later version.
*
*     MUESLI 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 MUESLI.  If not, see <http://www.gnu.org/licenses/>.
*
****************************************************************************/


#include <iostream>
#include <iomanip>
#include <cstdlib>
#include "muesli/muesli.h"

/* this function runs a cyclic test on an elastoplastic material
   point, and outputs results to a data for later plotting
*/
int main(int argc, char **argv)
{
    std::ofstream os("cyclictest.data");

    const double E = 210e9; // Young's modulus
    const double nu = 0.33;   // Poisson's ratio
    const double rho = 1.0;   // density
    const double Hiso = 1e9;  // isotropic hardening modulus
    const double Hkin = 2e9;  // kinematic hardening modulus
    const double Y0 = 100e6;  // yield stress
    const double alpha = 0.0; // not required
    const std::string yieldf = "mises";
    const std::string name = "an splastic material";
    
    muesli::splasticMaterial mat(name, E, nu, rho, Hiso, Hkin, Y0, alpha, yieldf);
    muesli::smallStrainMP* p = mat.createMaterialPoint();

    istensor strain;
    istensor shear;
    shear(0,1) = shear(1,0) = 0.02;

    unsigned npercycle = 20;
    unsigned ncycles   = 5;
    for (unsigned i=0; i<ncycles*npercycle; i++)
      {
    double t = sin(((double)i)/npercycle*2.0*M_PI);
    strain = shear*t;
    p->updateCurrentState(i, strain);

    istensor sigma; p->stress(sigma);
    double   tau = sigma(1,0);
    double   g   = p->plasticSlip();
    p->commitCurrentState();

    os << "\n" << std::setw(8) << std::scientific << t << " " << strain(0,1) << " " << tau << " " << g;
      }

    os.close();
    delete p;

    return 0;
}


