% FEAcomp. Solution of Example 8.5
% J.A. Mayugo & E.J. Barbero (c) 2006, 2007, 2013

clc; clear all; close all;format short g; szefont=20;

f0=figure;
subplot(1,2,1);
ezplot('sqrt(0.6104816314e0 * Y2s + 0.4047803847e1 * Y2s ^ 2 + 0.3849035026e0 * Y6s ^ 2) - 1'...
    ,[0,.6,0,1.75]);
xlabel('Y_2','Fontsize',szefont);ylabel('Y_6','Fontsize',szefont);
title('');set(gca,'Fontsize',szefont);
hold on;
plot([0,.6],[0,0],'--','Linewidth',1.4);
plot([0,0],[-2,2],'--','Linewidth',1.4);
hold off;
subplot(1,2,2);
ezplot('sqrt(0.6723365985e-4 * sigma2 ^ 2 + 0.4047803847e1 * (0.1101321586e-3 * sigma2 ^ 2 ) ^ 2 + 0.1375436417e-7 * sigma6 ^ 4)-1' ...
    ,[0,80,-100,100]);
xlabel('\sigma_2 [MPa]','Fontsize',szefont);ylabel('\sigma_6 [MPa]','Fontsize',szefont);
title('');set(gca,'Fontsize',szefont);
hold on;
plot([0,80],[0,0],'--','Linewidth',2.4,'color','k');
plot([0,0],[-100,100],'--','Linewidth',2.4,'color','k');
hold off;
saveas(f0,['figEx85a'],'eps');

% FEA Graph solution sigma2 TEST (PRVAR.LIS)
% time, SY, EPELY
FEsig2 =    [ 
  0.20000E-01     3.63200      0.400000E-03;
  0.40000E-01     7.26400      0.800000E-03;
  0.60000E-01     10.8960      0.120000E-02;
  0.80000E-01     14.5280      0.160000E-02;
  0.10000         18.1600      0.200000E-02;
  0.12000         21.7920      0.240000E-02;
  0.14000         25.4240      0.280000E-02;
  0.16000         29.0560      0.320000E-02;
  0.18000         32.6880      0.360000E-02;
  0.20000         36.3200      0.400000E-02;
  0.22000         39.9520      0.440000E-02;
  0.24000         43.5840      0.480000E-02;
  0.26000         47.2160      0.520000E-02;
  0.28000         50.8480      0.560000E-02;
  0.30000         54.4800      0.600000E-02;
  0.32000         58.1120      0.640000E-02;
  0.34000         61.7440      0.680000E-02;
  0.36000         55.3705      0.720000E-02;
  0.38000         48.3341      0.760000E-02;
  0.40000         42.4207      0.800000E-02;
  0.42000         37.4106      0.840000E-02;
  0.44000         33.1410      0.880000E-02;
  0.46000         29.4826      0.920000E-02;
  0.48000         26.3320      0.960000E-02;
  0.50000         23.6056      0.100000E-01;
  0.52000         21.2357      0.104000E-01;
  0.54000         19.1668      0.108000E-01;
  0.56000         17.3534      0.112000E-01;
  0.58000         15.7659      0.116000E-01;
  0.60000         14.3563      0.120000E-01;
  0.62000         13.1071      0.124000E-01;
  0.64000         11.9966      0.128000E-01;
  0.66000         11.0064      0.132000E-01;
  0.68000         10.1207      0.136000E-01;
  0.70000         9.32641      0.140000E-01;
  0.72000         8.61205      0.144000E-01;
  0.74000         7.96797      0.148000E-01;
  0.76000         7.38583      0.152000E-01;
  0.78000         6.85843      0.156000E-01;
  0.80000         6.37955      0.160000E-01;
  0.82000         5.94380      0.164000E-01;
  0.84000         5.54647      0.168000E-01;
  0.86000         5.18345      0.172000E-01;
  0.88000         4.85115      0.176000E-01;
  0.90000         4.54869      0.180000E-01;
  0.92000         4.26862      0.184000E-01;
  0.94000         4.01085      0.188000E-01;
  0.96000         3.77327      0.192000E-01;
  0.98000         3.55395      0.196000E-01;
   1.0000         3.35119      0.200000E-01];

% FEA Graph solution SHEAR TEST (PRVAR.LIS)
% time, SXY, EPELXY
FEshear =    [
  0.20000E-01     4.08290      0.771814E-03;
  0.40000E-01     8.16579      0.154363E-02;
  0.60000E-01     12.2487      0.231544E-02;
  0.80000E-01     16.3316      0.308726E-02;
  0.10000         20.4145      0.385907E-02;
  0.12000         24.4974      0.463088E-02;
  0.14000         28.5803      0.540270E-02;
  0.16000         32.6632      0.617451E-02;
  0.18000         36.7461      0.694633E-02;
  0.20000         40.8290      0.771814E-02;
  0.22000         44.9119      0.848996E-02;
  0.24000         48.9948      0.926177E-02;
  0.26000         53.0777      0.100336E-01;
  0.28000         57.1606      0.108054E-01;
  0.30000         61.2435      0.115772E-01;
  0.32000         65.3263      0.123490E-01;
  0.34000         69.4092      0.131208E-01;
  0.36000         73.4921      0.138927E-01;
  0.38000         77.5750      0.146645E-01;
  0.40000         81.6579      0.154363E-01;
  0.42000         85.7408      0.162081E-01;
  0.44000         89.8237      0.169799E-01;
  0.46000         93.9066      0.177517E-01;
  0.48000         97.9895      0.185235E-01;
  0.50000         102.072      0.192954E-01;
  0.52000         106.155      0.200672E-01;
  0.54000         110.238      0.208390E-01;
  0.56000         114.321      0.216108E-01;
  0.58000         118.404      0.223826E-01;
  0.60000         122.487      0.231544E-01;
  0.62000         126.570      0.239262E-01;
  0.64000         130.530      0.247114E-01;
  0.66000         122.904      0.255634E-01;
  0.68000         115.903      0.264112E-01;
  0.70000         109.450      0.272546E-01;
  0.72000         103.494      0.280948E-01;
  0.74000         97.9756      0.289320E-01;
  0.76000         92.8470      0.297667E-01;
  0.78000         88.0811      0.305989E-01;
  0.80000         83.6397      0.314290E-01;
  0.82000         79.4894      0.322572E-01;
  0.84000         75.6135      0.330835E-01;
  0.86000         71.9848      0.339081E-01;
  0.88000         68.5796      0.347312E-01;
  0.90000         65.3869      0.355530E-01;
  0.92000         62.3870      0.363734E-01;
  0.94000         59.5620      0.371926E-01;
  0.96000         56.9048      0.380107E-01;
  0.98000         54.4003      0.388278E-01;
   1.0000         52.0376      0.396439E-01];

[m,n]=size (FEshear');

time  = FEshear(:,1)';
epsilon = (FEshear(:,3)*100/2)';
sigma   = FEshear(:,2)';

epsilon2 = (FEsig2(:,3)*100)';
sigma2   = FEsig2(:,2)';

f1=figure;
hold on;
plot([0 epsilon],[0 sigma],'-','LineWidth',2.4,'color','k');
plot([0 epsilon2],[0 sigma2],'--','LineWidth',2.4,'color','k'); 
%plot(1.5,sigma_cr,'o','LineWidth',1.5,'MarkerSize',10);
legend('Transverse stress','In-plane shear stress','Location','NorthWest');legend('boxoff');
xlabel('Strain [%]','Fontsize',szefont);ylabel('Stress [MPa]','Fontsize',szefont);
%linetxt1 = ['    \leftarrow \sigma_{cr}=3060 MPa']
%linetxt2 = ['\angle (1-D_{cr})E_f']
%linetxt3 = ['    where D_{cr}=10.63%']
%text(1.5,sigma_cr,linetxt1,'HorizontalAlignment','left','Fontsize',szefont)
%text(0.8,1500,linetxt2,'HorizontalAlignment','left','Fontsize',szefont)
%text(0.8,1200,linetxt3,'HorizontalAlignment','left','Fontsize',szefont)
set(gca,'Fontsize',szefont);
axis([0 2 0 140]);
hold off;
saveas(f1,['figEx85b'],'pdf');

