% =====================================================================
% ----    Example of Finite Element Method on plain stress state   ----
% Solves an elasticity problem for a rectangular domain
% using a regular mesh of triangular linear elements
% The boundary conditions are dirichlet homogeneous and newmann.
%
% working on GNU Octave v. 3.6.4
% @ Copyleft, all wrongs reserved
%
% Date      Description                           Author
% Sep-2009  Creation                              J.M. Perez Zerpa
% Sep-2013  test on octave 3.6.4, comments added. J.M. Perez Zerpa
% Sep-2014  Inclusion vtk output functions        J.   Perlas
% ---------------------------------------------------------------------

function fem_rectangle_example

close all

% =========================================================
% ====================   Parameters   =====================

% ------------- Material properties -----------------------
E   = 200e6  ;  % Young modulus
nu  = 0.25   ;  % Poisson modulus

% -------------------- Mesh properties -------------------
L   = [1 4.0 ]' ;  % L(i) dimension in ith direction
t   = 0.2   ;  % thickness

% --- mesh ---
nx  = 10      ;  % number of divisions in direction x
ny  = 20      ;  % number of divisions in direction y
nel = [ nx ny ]' ; 

% ------------- Mesh Construction  ----------------
nnos    = nel + 1          ;  % nnos(i): number of noos in the i-th direction
nnostot = nnos(1)*nnos(2)  ;  % total nodes number

% ---------------- Boundary Conditions --------------------
q   = 1e6  ;  % distributed load on the up_edge
%~ LoadNodes = ( ( (nnos(2)-1) * nnos(1)+1 ) : ( nnos(2)*nnos(1)) ) ; 
startloadnode = (nnos(2)-1) * nnos(1)+1 ;
endinloadnode = nnos(1) * nnos(2)     ;
LoadNodes = ( startloadnode : endinloadnode )' ; 

% --- nodes cordinates matrix ---
lins1   = linspace( 0 , L(1) , nnos(1) )';
lins2   = linspace( 0 , L(2) , nnos(2) )';

% =========================================================


% =========================================================
% =================     Preprocessing     =================

% --- Constitutive Matrix ---
% matrix for plain strees state
C = E / (1-nu^2) * [ 1   nu  0           ; ...
                     nu  1   0           ; ...
                     0   0   (1-nu )/2 ] ;


nodes = [] ;  % nodes first column x-coord  second column y-coord

for i = 1:nnos(2) % discretiza primero por y
  nodes( ( (nnos(1) * (i-1) + 1) : nnos(1)*i ) ,:) = [ lins1 lins2(i)*ones(nnos(1),1) ] ;
end


% --- Conectivity Matrix ---
Mcon = []; %conectivity matrix
for j = 1:nel(2)
  for i = 1:nel(1)
    intri1 = 2*(i-1)+1+2*(j-1)*nx ;
    intri2 = intri1 + 1           ;
    Mcon( intri1 , : ) = [ j*nnos(1)+i    (j-1)*nnos(1)+i  j*nnos(1)+i+1     ] ;
    Mcon( intri2 , : ) = [ j*nnos(1)+i+1  (j-1)*nnos(1)+i  (j-1)*nnos(1)+i+1 ] ; 
  end
end
%~ Mcon
%~ nodes
% ------------------------


% --- Degrees of freedom fix and free ---
FixNodes  = ( 1:nnos(1) )'          ;
FixDofs   = nodes2dofs ( FixNodes )      ; % first add all dofs of the nodes 
FixDofs (3:2:(length(FixDofs)-2)) = [] ; % then remove the free dofs of the nodes
% free degrees of fredom
LibDofs = 1:2*nnostot ;
LibDofs ( FixDofs) = [] ;



% --- Build Elements coordinates vectors
neltotal = 2*nel(1)*nel(2); % number of elements

for i = 1:neltotal
    Xel(:,i) = (nodes( Mcon(i,:) , 1 )) ; % on j column has de x of the j-th element 
    Yel(:,i) = (nodes( Mcon(i,:) , 2 )) ; % on j column has de y of the j-th element 
end
% ===============================================================



% ===============================================================
% ==================      Processing       ======================

% ------------------   Stifness Matrix  ---------------------
KG = sparse( 2*nnostot,2*nnostot );
for i = 1:neltotal
  [Ke] = KT3(Xel(:,i),Yel(:,i),C , t );
  aux  = nodes2dofs( Mcon(i,:) )';
  KG( aux,aux ) = KG( aux,aux ) + Ke;
end
% -----------------------------------------------------------

%-----------------   Independent Term  ---------------------

areaelemsup     = L(1) / nel(1) * t ;
f               = zeros( 2*nnostot , 1 ) ;
f(2*LoadNodes(1)) = 0.5 * areaelemsup ;
for i=2:(length(LoadNodes)-1)
  %
  % fuerza horizontal 
  f(2*LoadNodes(i)-1) = 1 * areaelemsup ;            
  %
  % o
  %
  %fuerza vertical
  %~ f(2*LoadNodes(i)) = 1 * areaelemsup ;            
  %
end
f(2*LoadNodes(length(LoadNodes)))= 0.5 * areaelemsup;

F = q * f ;
%-----------------------------------------------------------

%-----------------   System Resolution  --------------------
KG ( : , FixDofs ) = [] ;
KG ( FixDofs , : ) = [] ;
K                  = KG ;

F  ( FixDofs     ) = [] ;

u = K\F ;
U = zeros( 2 * nnostot , 1 ) ; 

U( LibDofs ) = u ;

Ux = U( 1 : 2 : 2*nnostot-1 ) ;
Uy = U( 2 : 2 : 2*nnostot   ) ;

%-----------------------------------------------------------


epsilons=[];
sigmas=[];

for i=1:neltotal
	A=det( [ ones(1,3) ; Xel(:,i)' ; Yel(:,i)' ] )/2 ; % Area del elemento
	B= 1/(2*A)* [Yel(:,i)(2)-Yel(:,i)(3) 0          			 Yel(:,i)(3)-Yel(:,i)(1) 0          			 Yel(:,i)(1)-Yel(:,i)(2) 0          ;...
				 0          			 Xel(:,i)(3)-Xel(:,i)(2) 0          			 Xel(:,i)(1)-Xel(:,i)(3) 0          			 Xel(:,i)(2)-Xel(:,i)(1)  ;...
				 Xel(:,i)(3)-Xel(:,i)(2) Yel(:,i)(2)-Yel(:,i)(3) Xel(:,i)(1)-Xel(:,i)(3) Yel(:,i)(3)-Yel(:,i)(1) Xel(:,i)(2)-Xel(:,i)(1) Yel(:,i)(1)-Yel(:,i)(2) ] ;

  % determinacion de los grados de libertad del elemento
	glibelem = nodes2dofs ( Mcon(i,:) ) ;
	epsilons = [epsilons ; ( B * U(glibelem) )' ]; % calculo de las deformaciones de cada elemento
	sigmas = [sigmas ; ( C* epsilons(end,:)' )'] ; % calculo de las nuevas tensiones de cada elemento
end




% =================================================================
% ======================    Ploting     ===========================

%----------------  Undeformed Mesh  ------------------------
i=1; 
plot([Xel(:,i);Xel(1,i)],[Yel(:,i);Yel(1,i)],'r-o','LineWidth',1.4),
hold on,
for i=2:neltotal
    plot([Xel(:,i);Xel(1,i)],[Yel(:,i);Yel(1,i)],'r-o','LineWidth',1.4),
end
%-----------------------------------------------------------

% repeat nodes matrix construction with new desplacements
nodesdef(:,1)= nodes(:,1) + Ux ; 
nodesdef(:,2)= nodes(:,2) + Uy ;

% Build Elements coordinates vectors
for i=1:neltotal
    Xel(:,i)=(nodesdef( Mcon(i,:) , 1 )) ; % on j column has de x of the j-th element 
    Yel(:,i)=(nodesdef( Mcon(i,:) , 2 )) ; % on j column has de y of the j-th element 
end

%----------------  Deformed Mesh  ------------------------
for i=1:neltotal
    plot([Xel(:,i);Xel(1,i)],[Yel(:,i);Yel(1,i)],'b-o','LineWidth',1.5),%y(i))
end

axis equal, xlabel('x'), ylabel('y'),  
title('Displacements')
%===================================================================


disp('ATENCION: cree la carpeta output y si ya existe vaciela por favor!!')
intIteraciones=input('Ingrese el numero de cuadros para la animacion: ->');

for j=0:intIteraciones
	%vtk_2d (j,intIteraciones,"output","MEF_CST",nnostot,nodes,Ux,Uy,neltotal,Mcon',epsilons,sigmas)
	vtk_3d (j,intIteraciones,"output3d","MEF_CST",t,nnostot,nodes,Ux,Uy,neltotal,Mcon',epsilons,sigmas)
end
disp('listo. Abra paraview para ver los vtks')



%------------------------------------------------------------------
function v= nodes2dofs (u)
for i=1:length(u)
    v(2*i-1:2*i)=[ 2*u(i)-1 2*u(i)];        
end


%-------------------------------------------------------------
% calculates element matrices for clockwise numeration
function [K] = KT3 ( x , y , C , t )

A=det( [ ones(1,3) ; x' ; y' ] )/2; % element area
    
B= 1 / (2*A) * [ y(2)-y(3)  0          y(3)-y(1)  0          y(1)-y(2)  0          ;
                 0          x(3)-x(2)  0          x(1)-x(3)  0          x(2)-x(1)  ;
                 x(3)-x(2)  y(2)-y(3)  x(1)-x(3)  y(3)-y(1)  x(2)-x(1)  y(1)-y(2)  ];
K= B' * C * B * A * t ;
%--------------------------------------------------------


function vtk_2d (intIteracion,intIteraciones,strOutput,strTitulo,intNodos,arrNodos,arrDespX,arrDespY,intElementos,arrConect,arrDeform,arrTensiones)
	format short
	
	intCifras=ceil(log10(intIteraciones+1));	
	iterationid=int2str(intIteracion);
	for i=1:(intCifras-1)
		if intIteracion/(10^i) < 1
			iterationid = strcat(int2str(0),iterationid);
		endif
	endfor
	
	[err,msg] = mkdir ( "output" );
	strArchivo = ["./output/" strOutput "_" iterationid ".vtk"]
	fid=fopen(strArchivo,'w');

	fdisp(fid,"# vtk DataFile Version 2.0")
	fdisp(fid,strTitulo)
	fdisp(fid,"ASCII")
	fdisp(fid,"DATASET UNSTRUCTURED_GRID")
	fdisp(fid,"")

	fdisp(fid,["POINTS " int2str(intNodos) " float"])
	for i=1:intNodos
		fdisp(fid,["   " mat2str(arrNodos(i,1)+intIteracion/intIteraciones*arrDespX(i),5) "   " mat2str(arrNodos(i,2)+intIteracion/intIteraciones*arrDespY(i),5) "   0"])
	endfor

	arrConect = arrConect .- 1;
	fdisp(fid,"")
	fdisp(fid,["CELLS " int2str(intElementos) " " int2str(4*intElementos)])  
	for i=1:intElementos
		fdisp(fid,["   3   " int2str(arrConect(1,i)) "   " int2str(arrConect(2,i)) "   " int2str(arrConect(3,i))])
	endfor

	fdisp(fid,"")
	fdisp(fid,["CELL_TYPES " int2str(intElementos)])  
	for i=1:intElementos
		fdisp(fid,"   5")
	endfor

	fdisp(fid,"")
	fdisp(fid,["POINT_DATA " int2str(intNodos)])
	fdisp(fid,"VECTORS desplazamientos float")
	for i=1:intNodos
		fdisp(fid,["   " mat2str(intIteracion/intIteraciones*arrDespX(i),5) "   " mat2str(intIteracion/intIteraciones*arrDespY(i),5) "   0"])
	endfor

	fdisp(fid,"")
	fdisp(fid,["CELL_DATA " int2str(intElementos)])
	fdisp(fid,"VECTORS deformaciones float")
	for i=1:intElementos
		fdisp(fid,["   " mat2str(intIteracion/intIteraciones*arrDeform(i,1),5) "   " mat2str(intIteracion/intIteraciones*arrDeform(i,2),5) "   " mat2str(intIteracion/intIteraciones*arrDeform(i,3),5)])
	endfor

	fdisp(fid,"")
	fdisp(fid,"VECTORS sigmas float")
	for i=1:intElementos
		fdisp(fid,["   " mat2str(intIteracion/intIteraciones*arrTensiones(i,1),5) "   " mat2str(intIteracion/intIteraciones*arrTensiones(i,2),5) "   " mat2str(intIteracion/intIteraciones*arrTensiones(i,3),5)])
	endfor
	fclose(fid);



function vtk_3d (intIteracion,intIteraciones,strOutput,strTitulo,espesor,intNodos,arrNodos,arrDespX,arrDespY,intElementos,arrConect,arrDeform,arrTensiones)
	format short
	
	intCifras=ceil(log10(intIteraciones+1));
	iterationid=int2str(intIteracion);
	for i=1:(intCifras-1)
		if intIteracion/(10^i) < 1
			iterationid = strcat(int2str(0),iterationid);
		endif
	endfor
	
	[err,msg] = mkdir ( "output" );
	strArchivo = ["./output/" strOutput "_" iterationid ".vtk"] ;
	fid=fopen(strArchivo,'w');

	fdisp(fid,"# vtk DataFile Version 2.0")
	fdisp(fid,strTitulo)
	fdisp(fid,"ASCII")
	fdisp(fid,"DATASET UNSTRUCTURED_GRID")
	fdisp(fid,"")

	fdisp(fid,["POINTS " int2str(2*intNodos) " float"])
	for i=1:intNodos
		fdisp(fid,["   " mat2str(arrNodos(i,1)+intIteracion/intIteraciones*arrDespX(i),5) "   " mat2str(arrNodos(i,2)+intIteracion/intIteraciones*arrDespY(i),5) "   0"])
	endfor
	for i=1:intNodos
		fdisp(fid,["   " mat2str(arrNodos(i,1)+intIteracion/intIteraciones*arrDespX(i),5) "   " mat2str(arrNodos(i,2)+intIteracion/intIteraciones*arrDespY(i),5) "   " mat2str(espesor)])
	endfor
	
	arrConect = arrConect .- 1;
	fdisp(fid,"")
	fdisp(fid,["CELLS " int2str(intElementos) " " int2str(7*intElementos)])  
	for i=1:intElementos
		fdisp(fid,["   6   " int2str(arrConect(1,i)) "   " int2str(arrConect(2,i)) "   " int2str(arrConect(3,i)) "   " int2str(intNodos+arrConect(1,i)) "   " int2str(intNodos+arrConect(2,i)) "   " int2str(intNodos+arrConect(3,i))])
	endfor

	fdisp(fid,"")
	fdisp(fid,["CELL_TYPES " int2str(intElementos)])  
	for i=1:intElementos
		fdisp(fid,"   13")
	endfor

	fdisp(fid,"")
	fdisp(fid,["POINT_DATA " int2str(2*intNodos)])
	fdisp(fid,"VECTORS desplazamientos float")
	for i=1:intNodos
		fdisp(fid,["   " mat2str(intIteracion/intIteraciones*arrDespX(i),5) "   " mat2str(intIteracion/intIteraciones*arrDespY(i),5) "   0"])
	endfor
	for i=1:intNodos
		fdisp(fid,["   " mat2str(intIteracion/intIteraciones*arrDespX(i),5) "   " mat2str(intIteracion/intIteraciones*arrDespY(i),5) "   0"])
	endfor

	fdisp(fid,"")
	fdisp(fid,["CELL_DATA " int2str(intElementos)])
	fdisp(fid,"VECTORS deformaciones float")
	for i=1:intElementos
		fdisp(fid,["   " mat2str(intIteracion/intIteraciones*arrDeform(i,1),5) "   " mat2str(intIteracion/intIteraciones*arrDeform(i,2),5) "   " mat2str(intIteracion/intIteraciones*arrDeform(i,3),5)])
	endfor

	fdisp(fid,"")
	fdisp(fid,"VECTORS sigmas float")
	for i=1:intElementos
		fdisp(fid,["   " mat2str(intIteracion/intIteraciones*arrTensiones(i,1),5) "   " mat2str(intIteracion/intIteraciones*arrTensiones(i,2),5) "   " mat2str(intIteracion/intIteraciones*arrTensiones(i,3),5)])
	endfor
	fclose(fid);



