function Calc_ee_se_vvf_eta
  u = dlmread('VoidedCube_u.dat');
  f = dlmread('VoidedCube_f.dat');
  solid_vol = dlmread('VoidedCube_solid_vol.dat');

  vvf = Calc_VVF(u, solid_vol);
  ee = Calc_ee(u);
  [se, eta] = Calc_se(u, f);

  d = [ee, se, vvf, eta];
  dlmwrite('VoidedCube_ee_se_vvf_eta.dat', d, 'delimiter', ' ', ...
    'precision', '%.18e');
end

function vvf = Calc_VVF(u, vol_sol)
  % calculates void volume fraction

  % displaced coordinates of master corner node
  y = ones(size(u)) + u;

  % total volume of deformed cube
  vol_tot = y(:, 1).*y(:, 2).*y(:, 3);

  % volume of void
  vol_void = vol_tot - vol_sol;

  % void volume fraction
  vvf = vol_void./vol_tot;
end

function ee = Calc_ee(u)
  % calculates effective strain

  % get logarithmic strains E1, E2, E3
  % Ei = log(Li/Li0)
  %   Li = i^th deformed length
  %   Li0 = i^th reference length

  % deformed length
  Li = ones(size(u)) + u;

  % reference length = 1.0 (cube is 1.0^3)
  % therefore Li/Li0 == Li

  % logarithmic strains
  Ei = log(Li);

  E1 = Ei(:, 1);
  E2 = Ei(:, 2);
  E3 = Ei(:, 3);

  ee = sqrt(2.0)/3.0 * sqrt((E1 - E2).^2 + (E2 - E3).^2 + (E3 - E1).^2);
end

function [se, eta] = Calc_se(u, f)
  % calculates effective stress

  % true stress s_{ii} = load f_i / deformed area_i

  % loads
  f1 = f(:, 1);
  f2 = f(:, 2);
  f3 = f(:, 3);

  % get deformed coords
  y1 = 1.0 + u(:, 1);
  y2 = 1.0 + u(:, 2);
  y3 = 1.0 + u(:, 3);

  % s11
  % deformed area = y2*y3
  a1 = y2.*y3;
  s11 = f1./a1;

  % s22
  % deformed area = y1*y3
  a2 = y1.*y3;
  s22 = f2./a2;

  % s33
  % deformed area = y1*y2
  a3 = y1.*y2;
  s33 = f3./a3;

  % pressure stress
  p = -(s11 + s22 + s33)/3.0;

  % effective stress
  se = sqrt((s11 - s22).^2 + (s22 - s33).^2 + (s33 - s11).^2)/sqrt(2.0);

  % triaxiality
  eta = -p./se;
end
