% FEM program of elasticity in 2D for quadrilateral linear elements
% Makes use of: 
% natural_deriv_2D4N, jacobian_2D4N, Bmatrix_2D4N, element_stiffness_2D4N

clear;

% 0.- PREPROCESSING _____________________________________________________
% 0.1  Define nodes (coordinate list)
% This data do not enter  in the stiffness matrix in case of springs 
nodeCoordinates=(1/sqrt(2))*[0 -1;1 0;0 1;-1 0];

% 0.2  Define connectivities (element definition)
elementNodes=[1 2 3 4];

% Size of model
[numnodes nDof]=size(nodeCoordinates); %nē de nodos/nēde dof por nodo
[numelements nNodeElem]=size(elementNodes); % nē elementos/nē nodos por elemento
nDofsGlobal=nDof*numnodes; % nē nodos globales=(nē dof por nodo)x(número nodos)

% 0.3 Define material assignment for each element 
elementMat(1)=1;    %elementMat(2) para otro material 'elementMat=[1;2];'

% 0.4 Define propertie of each material
propMatlist.E=200E9;
propMatlist.nu=0.2;

% 0.5 Boundary conditions 
% 0.5.1. Define  dof with prescribed displacement 0 or force
prescribedDispNodes=[4 4 2 2];    % prescribed Displacement Nodes
prescribedDispDof=[1 2 1 2];      % prescribed Displacement DOF (nDof)
prescribedForcNodes=[];
prescribedForcDof=[];

% 0.5.2 Forces
u=0.005*sqrt(2);
forces=[];
displacements =[0 0 u 0];

% 1.- SOLVER _____________________________________________________________
% Elementary stiffness matrix
%1.1 Form stiffness matrix
stiffglob=zeros(nDofsGlobal);
displacement=zeros(nDofsGlobal,1);
force=zeros(nDofsGlobal,1);
ndofs_global=zeros(nDof*nNodeElem,1);
force(2*(prescribedForcNodes-1)+prescribedForcDof)=forces % applied forces
displacement(2*(prescribedDispNodes-1)+prescribedDispDof)=displacements % applied displacement

% for ielem=1:numelements %Loop in elements
% 	nodos_global=elementNodes(ielem,:); %Obtain dof numbers for element ielem
% 	j=0;
% 	for idof=nodos_global
% 		j=j+1;
% 		ndofs_global(j)=(2*idof)-1;
% 		j=j+1;
% 		ndofs_global(j)=(2*idof);
% 	end
% 	xnod(1,1)=nodeCoordinates(elementNodes(ielem,1),1);
% 	xnod(2,1)=nodeCoordinates(elementNodes(ielem,1),2);
% 	xnod(3,1)=nodeCoordinates(elementNodes(ielem,2),1);
% 	xnod(4,1)=nodeCoordinates(elementNodes(ielem,2),2);
% 	xnod(5,1)=nodeCoordinates(elementNodes(ielem,3),1);
% 	xnod(6,1)=nodeCoordinates(elementNodes(ielem,3),2);
% 	xnod(7,1)=nodeCoordinates(elementNodes(ielem,4),1);
% 	xnod(8,1)=nodeCoordinates(elementNodes(ielem,4),2);
% 	propMat=propMatlist(elementMat(ielem))
% 	stiffelem=element_stiffness_2D4N(propMat,xnod,2);
% 	stiffglob(ndofs_global, ndofs_global)=stiffglob(ndofs_global,...
%     ndofs_global)+stiffelem; % Add element stiffness matrix to global one (assembly)
% end
xnod=(1/sqrt(2))*[0;-1;1;0;0;1;-1;0];
stiffglob=element_stiffness_2D4N(propMatlist,xnod,2)

Fut=-stiffglob*displacement;
prescribedDofs=[2*(prescribedDispNodes-1)+prescribedDispDof]

%1.2 Eliminate equations in which applied u=0
activeDofs=setdiff([1:nDofsGlobal],[prescribedDofs])
% unactiveDofs=setdiff([1:nDofsGlobal],activeDofs) % =prescribedDofs

%1.3 Solve linear system
U=stiffglob(activeDofs,activeDofs)\(force(activeDofs)+Fut(activeDofs));

displacement(activeDofs)=U;
displacement(prescribedDofs)=displacements
force=stiffglob*displacement