%%this function gets coordinates of a CST (constant strain triangle) element
%and populates the N, B, and D matrices for plane stress or plane strain
%BY ABIGAIL OWENS 11/1/17


disp('enter the coordinates of nodes 1, 2, and 3 of the triangle, up to four decimal places')
disp('**********************************************************')
disp('counterclockwise, start with point 1.')

x1=input('enter X1 ');
y1=input('enter Y1 ');
fprintf('point 1 is (%0.4f , %0.4f)\n\n',x1,y1)
%tells user what point they've entered

disp('point 2:')
x2=input('enter X2 ');
y2=input('enter Y2 ');
fprintf('point 2 is (%0.4f , %0.4f)\n\n',x2,y2)
%tells user point 2

disp('point 3:')
x3=input('enter X3 ');
y3=input('enter Y3 ');
fprintf('point 3 is (%0.4f , %0.4f)\n\n',x3,y3)
%tells user point 3

%disp('what is the effective area of the triangle?')
%A=input('area: ');

disp('modulus of the element? make sure your units are consistent with the coordinates you entered')
E=input('modulus: ');

disp('poissons ratio?')
nu=input('ratio: ');

disp('element thickness?')
t=input('thickness: ');

    
Dtype = input('Will problem be plane stress [1] or plane strain [2]? ');
while (Dtype ~= 1 && Dtype ~= 2)  % Ensure the user inputs a correct number
    Dtype = input('Invalid Value! Try again.\n Plane Stress [1] or Plane Strain [2]? '); 
end

%point 1 N calcs
alpha1=(x2*y3)-(x3*y2);
beta1=(y2-y3);
gamma1=-(x2-x3);

%point 2 N calcs
alpha2=(x3*y1)-(x1*y3);
beta2=(y3-y1);
gamma2=-(x3-x1);

%point 3 N calcs
alpha3=(x1*y2)-(x2*y1);
beta3=(y1-y2);
gamma3=-(x1-x2);

%calculating area 
A=abs((x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2);

%making N matrix
N1=0.5*A*(alpha1+beta1+gamma1);
N2=0.5*A*(alpha2+beta2+gamma2);
N3=0.5*A*(alpha3+beta3+gamma3);
N=[N1 0 N2 0 N3 0;0 N1 0 N2 0 N3];

%making B matrix
B=(1/(2*A))*[beta1 0 beta2 0 beta3 0;0 gamma1 0 gamma2 0 gamma3; gamma1 beta1 gamma2 beta2 gamma3 beta3];

  if Dtype==1 %plane stress
    D= (E/(1-(nu)^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2]
   end
  if Dtype==2 %plane strain
        D= (E/((1-nu)*(1-2*nu)))*[(1-nu) nu 0;nu (1-nu) 0;0 0 (1-2*nu)/2]
   end 
    
  %making entire stiffness matrix
  transpB= transpose (B); %doing transpose of B matrix in different step just because
  k=t*A*transpB*D*B
  
  %stress matrices 
  d=input('Do you have a displacement matrix? enter yes [1] or no  [2] ');
  while (d ~= 1 && d ~= 2)  % Ensure the user inputs a valid entry
    d = input('Invalid Value! Try again.\n Displacement matrix known [1] or no [2]? '); 
end

if d==1
    disp('Enter displacement values:')%gets displacement values from user
    u1=input('U1: ');
    v1=input('V1: ');
    u2=input('U2: ');
    v2=input('V2: ');
    u3=input('U3: ');
    v3=input('V3: ');
    
    %put values into displacement matrix
    displacement=[u1;v1;u2;v2;u3;v3];
    
    %calculate stress and tau xy
    stress=D*B*displacement;
    sigmax=stress(1,1); %assigning variable names to components of stress matrix
    sigmay=stress(2,1);
    tauxy=stress(3,1);
    
    %tell stress to user
    fprintf('stress in the x direction is %0.4f \n',sigmax)
    fprintf('stress in the y direction is %0.4f \n',sigmay)
    fprintf('stress in the xy direction is %0.4f \n',tauxy)
    
    %calculate principal stresses and principle angle
    princ1=((sigmax+sigmay)/2)+sqrt(((sigmax-sigmay)/2)^2+(tauxy)^2);
    princ2=((sigmax+sigmay)/2)-sqrt(((sigmax-sigmay)/2)^2+(tauxy)^2);
    theta=(atan((2*tauxy)/(sigmax-sigmay)))/2;
    
    %tell to user
    fprintf('principal stress 1: %0.4f \n',princ1)
    fprintf('principal stress 2: %0.4f \n',princ2)
    fprintf('principal angle: %0.4f \n', theta)
    
    %von Mises stresses for plane stress (sigma z = 0)
    if Dtype==1 
        sigmaz=0;
        sigmavm=sqrt((1/(sqrt(2)))*((sigmax-sigmay)^2+(sigmax-sigmaz)^2+(sigmay-sigmaz)^2+6*(tauxy^2)));
        fprintf('von Mises stress: %0.4f \n',sigmavm)
    end 
    %von Mises stress for plane strain (sigma z = something)
    if Dtype==2
        sigmaz=nu*(sigmax+sigmay);
        sigmavm=(1/(sqrt(2)))*sqrt(((sigmax-sigmay)^2+(sigmax-sigmaz)^2+(sigmay-sigmaz)^2+6*(tauxy^2)));
        fprintf('von Mises stress: %0.4f \n',sigmavm)
    end
end

%if no displacement, end program
if d==2
    disp('sounds good, see you next time')
end 
