#include<iostream>
using namespace std;

#include "mat.h"
using namespace ariadne;

int main()
{
   cout << "# Uniaxial strain\n";
   cout << "# Ariadne V" << ariadne::version() << endl;
   cout << "# WXMATHS V" << wxmaths::version() << endl;
   char matname[80]; cin >> matname; cout << "# material " << matname << endl;
   mat_type *a_mat;
   cin >> a_mat;
// cout << "mat*:\n" << a_mat << endl << "mat:\n" << *a_mat << endl;
   mat_type::state *st0 = a_mat->new_state(); cin >> *st0;
   cout.precision(10);
   cout << "# initial state = " << *st0 << endl;
   cout << "# p = " << a_mat->p(st0) << '\n';
   cout << "# c^2 = " << a_mat->csqv(st0) << '\n';
   cout << "# T = " << a_mat->t(st0) << '\n';

   Matrix3d gradu; cin >> gradu;
// cout << "# gradu = " << gradu << endl;

   mat_type::state *stdot = a_mat->new_state();
   a_mat->statedot_mech(st0,gradu,stdot);

   cout << "# initial stdot = " << *stdot << endl;

   double dt; cin >> dt; cout << "# dt = " << dt << endl;

   mat_type::state *st = a_mat->new_state(); a_mat->copy(st0,st);
   mat_type::state *st1 = a_mat->new_state();
   cout << "# t rho e stress T csqv cv p	st\n";
   double rho = a_mat->mass_density(st), e = a_mat->e(st),
          pres = a_mat->p(st), temp = a_mat->t(st),
          csqv = a_mat->csqv(st), cv = a_mat->cv(st);
   Matrix3d_symm_voigt tau = a_mat->stress(st,gradu);
   cout << 0.0 << '\t' << rho << ' ' << e << ' ' << tau << ' ' << temp << ' '
        << csqv << ' ' << cv << ' ' << pres << '\t' << *st << endl;
   int n; cin >> n;
   for (int i=0;i<n;i++) {
      double t = i * dt;

      // predictor
      a_mat->add_update(st,stdot,dt/2,st1);
      a_mat->evolve_internal(st1,dt/2);
      a_mat->statedot_mech(st1,gradu,stdot);

      // corrector
      a_mat->add_update(st,stdot,dt,st1);
      a_mat->evolve_internal(st1,dt);
      a_mat->copy(st1,st);

      rho = a_mat->mass_density(st); e = a_mat->e(st);
      pres = a_mat->p(st); temp = a_mat->t(st);
      csqv = a_mat->csqv(st); cv = a_mat->cv(st);
      tau = a_mat->stress(st,gradu);
      cout << t << '\t' << rho << ' ' << e << ' ' << tau << ' ' << temp << ' '
           << csqv << ' ' << cv << ' ' << pres << '\t' << *st << endl;

      a_mat->statedot_mech(st,gradu,stdot);
   }
}