% FEAcomp. Solution of Example 6.2.
%
% J.A. Mayugo, WVU, 2006


% To run this M-file need the functions:
%       fitfunFEAcomp(xdata, ydata,ifun)
%       invlapFEAcomp(time,func)
%       PMMViscoMatrix( time, V_f, imod_m, imod_f, Elas_prop_m, Elas_prop_f, Visc_prop_m)
% available in .\Examples\Ch7ex\      

clear all; close all;format short g;

%% domain time
time = [0];   % not used in elastic solution

% Volume fraction 
V_f = 0.4;

%% Matrix Vinylester 
E_m=3.12; % GPa  
nu_m = 0.38;

% Elastic
imod_m=0; 

%%% Fiber: 
for imod_f=1:1;  % 1,1 only isotropic, 2:2 only transversaly, 1:2 both
    if imod_f==1 % isotropic
        E_f	= 241.; %GPa 24.42;    % msi 
        nu_f  = 0.2;	
    elseif imod_f==2 % transversely isotropic																					
        E_a  = 241.;  E_t = 24.8; % GPa 3.6;     % msi 
        nu_a = 0.2;  nu_t = 0.005; 
        G_a  = 44.1e6; % GPa 6.4;    
        G_t = E_t/2/(1+nu_t); % GPa, calculated
    end

    if imod_f==1       % Isotropic fibers, requires imod_f =1
        Elas_prop_f=[E_f nu_f];     
    elseif imod_f==2    % Transversaly Isotropic fibers, requires imod_f =2
        Elas_prop_f=[E_a E_t nu_a nu_t G_a G_t];     
    end
    if imod_m==1         % Maxell model, requires imod_m = 1
        Visc_prop_m=[eta_m];
    else
        Visc_prop_m=[0];
    end
    Elas_prop_m=[E_m nu_m];
    [Elas_prop, Visc_prop,model] = PMMViscoMatrix( time, V_f, imod_m, imod_f,...
                                             Elas_prop_m, Elas_prop_f, Visc_prop_m);
    %% Elastic properties
    Eo_1=Elas_prop(1); Eo_2=Elas_prop(2);
    Eo_12=Elas_prop(3); Eo_23=Elas_prop(4);
    nuo_12=Elas_prop(5); nuo_23=Elas_prop(6);

    disp('Composite: Elastic properties');
    disp('Eo_1, Eo_2, Go_12, Go_23, nuo_12, nuo_23');
    disp([Eo_1 Eo_2 Eo_12 Eo_23 nuo_12 nuo_23]);
    
end

%%% Calculate RVE dimensions
syms a2 a3 real

d_f = 7 

eq1 = V_f - pi*(d_f/2)^2/(2*a2*a3);
a3 = a2*tan(pi/3);

a2= solve(eval(eq1),a2);

a2sol = double(a2)
a3sol = double(eval(a3))

%% Solution Example_6_3

C11 =    0.1012950;
C21 =    0.5009453E-02;
C31 =    0.5008952E-02;

C12 =    0.5010452E-02;
C22 =    0.1081295E-01;
C32 =    0.5956267E-02;

C13 =    0.5008050E-02;
C23 =    0.5954785E-02;
C33 =    0.1080788E-01;

%! PARAMETER EL =    0.9830064E-01
% PARAMETER NUL =    0.2988150
% PARAMETER ET =    0.7482400E-02
% PARAMETER NUT =    0.5401690
% PARAMETER GT =    0.2429084E-02

C44 =    0.2422171E-02;
C55 =    0.2582067E-02;
C66 =    0.2583918E-02;


C_Ex63= 1000 * [C11 C12 C13 0 0 0; C21 C22 C23 0 0 0 ; C31 C32 C33 0 0 0 ;
        0 0 0 C44 0 0;0 0 0 0 C55 0 ; 0 0 0 0 0 C66  ] ;

% COMPUTE TRANSFORMATION MATRIX (Clokwise rotation of 12 respect to xy)
theta = pi/4;

a = [1  0           0;
     0  cos(theta)  sin(theta) ;
     0  -sin(theta) cos(theta) ]; 
 
     
for i1=1:3
    for j1=1:3
        if i1==j1 
            alpha=i1;
        else
            alpha=9-i1-j1;
        end
        for p1=1:3
            for q1=1:3
                if p1==q1 
                    beta=p1;
                else
                    beta=9-p1-q1;
                end
 %               te(alpha,beta)=0
                if (alpha<=3)&(beta<=3)
                    T(alpha,beta)=a(i1,p1)*a(i1,p1);
                elseif (alpha>3)&(beta<=3)
                    T(alpha,beta)=a(i1,p1)*a(j1,p1);
                elseif (alpha<=3)&(beta>3)
                    T(alpha,beta)=a(i1,q1)*a(i1,p1)+a(i1,p1)*a(i1,q1);
                elseif (alpha>3)&(beta>3)
                    T(alpha,beta)=a(i1,p1)*a(j1,q1)+a(i1,q1)*a(j1,p1);
                end
            end
        end
    end
end
%T
% COMPUTE C,S in new coordiante system

C=inv(T)*C_Ex63*inv(T)';
S=inv(C);

E1 = 1/S(1,1);
E2 = 1/S(2,2);
E3 = 1/S(3,3);
nu12=-S(2,1)*E1;
nu21=-S(1,2)*E2;
nu13=-S(3,1)*E1;
nu31=-S(1,3)*E3;
nu32=-S(2,3)*E3;
nu23=-S(3,2)*E2;

disp('Composite: Elastic properties Ex_63 (rotated)');
disp('E1, E2, E3, nu12, nu13, nu23, G23, G13, G12');
disp([E1, E2, E3, nu12, nu13, nu23 1/S(4,4) 1/S(5,5) 1/S(6,6)]); 

S=inv(C_Ex63) ;

E1 = 1/S(1,1);
E2 = 1/S(2,2);
E3 = 1/S(3,3);
nu12=-S(2,1)*E1;
nu21=-S(1,2)*E2;
nu13=-S(3,1)*E1;
nu31=-S(1,3)*E3;
nu32=-S(2,3)*E3;
nu23=-S(3,2)*E2;

disp('Composite: Elastic properties Ex_63 (no rotated)');
disp('E1, E2, E3, nu12, nu13, nu23, G23, G13, G12');
disp([E1, E2, E3, nu12, nu13, nu23 1/S(4,4) 1/S(5,5) 1/S(6,6)]); 
