#include "mat.h"

int main()
{
   mat_type *a_mat;
   cin >> a_mat;
   cout << "mat*:\n" << a_mat << endl << "mat:\n" << *a_mat << endl;
   mat_type::state *st0 = a_mat->get_state(); cin >> *st0;
   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 = \n" << gradu << endl;

   mat_type::state *stdot = a_mat->get_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->get_state(); a_mat->Copy(st0,st);
   int n; cin >> n;
   cout << "# t st p T stress_deviator\n";
   for (int i=0;i<n;i++) {
      double t = i * dt;
      a_mat->add_update(st,stdot,dt,st);
      double pres = a_mat->p(st), temp = a_mat->t(st);
      Vector sig = deviators(a_mat->stress(st,gradu));
      cout << t << '\t' << *st << '\t' << pres << ' ' << temp << '\t';
      sig.Write_elements(cout); cout << endl;
      a_mat->statedot_mech(st,gradu,stdot);
   }
}