Code covered by the BSD License

# 3D Solid FEM

3D FEM analysis using linear tetrahedral elements. Also a modal analysis is shown.

ExternalForces0D1D2D.m
```function f=ExternalForces0D1D2D(f,PermutationList,External0DElementForce,External1DElementForce,External2DElementForce,x,y,z,DegreesOfFreedom)
% Assembly of external loads. This function stores pointloads (0D), integrates line tractions (1D), surface tractions (2D), and
% assembles (stores) them. Element volumeloads are integrated and assembled together with the
% element stiffness matrices in the routine called 'ElementMatrix'.
for i=1:numel(External0DElementForce)
node1=External0DElementForce{i}(1);
xdof1=DegreesOfFreedom*node1-2; % xdof of node of 0D element============
f(PermutationList(xdof1))  =f(PermutationList(xdof1))  +External0DElementForce{i}(2);
f(PermutationList(xdof1+1))=f(PermutationList(xdof1+1))+External0DElementForce{i}(3);
f(PermutationList(xdof1+2))=f(PermutationList(xdof1+2))+External0DElementForce{i}(4);
end
psi_Gausspoints=0;
weights=2;
for i=1:numel(External1DElementForce)
% Initialize the external 1D element loadvector=========================
f_element=zeros(2*DegreesOfFreedom,1);
% First  node of 1D element=============================================
node1=External1DElementForce{i}(1);
% Second node of 1D element=============================================
node2=External1DElementForce{i}(2);
%=======================================================================
% Determinant of Jacobian for 1D lineair element========================
LengthOfElement=norm([x(node2)-x(node1),y(node2)-y(node1),z(node2)-z(node1)]);
detJ=1/2*LengthOfElement;
for j=1:numel(weights)
psi=psi_Gausspoints(j);
N=1/2*[1-psi    0     0      1+psi 0     0;
0        1-psi 0      0     1+psi 0;
0        0     1-psi  0     0     1+psi];
f_element=f_element+(N')*...
External1DElementForce{i}(3:5)*weights(j)*detJ;
end
% f_element yields 1/2*[f1D_x f1D_y f1D_z]*length_element===============
xdof1=DegreesOfFreedom*node1-2;     % xdof of left node of 1D element
xdof2=DegreesOfFreedom*node2-2;     % xdof of right node of 1D element
f(PermutationList(xdof1))  =f(PermutationList(xdof1))  +f_element(1);
f(PermutationList(xdof1+1))=f(PermutationList(xdof1+1))+f_element(2);
f(PermutationList(xdof1+2))=f(PermutationList(xdof1+2))+f_element(3);
f(PermutationList(xdof2))  =f(PermutationList(xdof2))  +f_element(4);
f(PermutationList(xdof2+1))=f(PermutationList(xdof2+1))+f_element(5);
f(PermutationList(xdof2+2))=f(PermutationList(xdof2+2))+f_element(6);
end
% Page 200, Eugenio Oate==================================================
xi_Gausspoints =1/3;
eta_Gausspoints=1/3;
weights=1/2;
for i=1:numel(External2DElementForce)
% Initialize the external 1D element loadvector=========================
f_element=zeros(3*DegreesOfFreedom,1);
% First  node of 2D element=============================================
node1=External2DElementForce{i}(1);
% Second node of 2D element=============================================
node2=External2DElementForce{i}(2);
% Third node of 2D element==============================================
node3=External2DElementForce{i}(3);
% Determinant of Jacobian for 2D lineair element========================
% Area of element=======================================================
A=1/2*norm(cross([x(node3)-x(node1),y(node3)-y(node1),z(node3)-z(node1)],...
[x(node3)-x(node2),y(node3)-y(node2),z(node3)-z(node2)]));
% Determinant (twice the area)==========================================
detJ=2*A;
for j=1:numel(weights)
xi=xi_Gausspoints(j);
eta=eta_Gausspoints(j);
% Shape functions for 2D triangular element==========================
N=[1-xi-eta 0        0        xi 0   0  eta 0   0;
0        1-xi-eta 0        0  xi  0  0   eta 0;
0        0        1-xi-eta 0  0   xi 0   0   eta];
f_element=f_element+(N')*...
External2DElementForce{i}(4:6)*weights(j)*detJ;
end
% f_element yields 1/3*[f2D_x f2D_y f2D_z]*area_element=================
xdof1=DegreesOfFreedom*node1-2;     % xdof of node 1 of 2D element
xdof2=DegreesOfFreedom*node2-2;     % xdof of node 2 of 2D element
xdof3=DegreesOfFreedom*node3-2;     % xdof of node 3 of 2D element
f(PermutationList(xdof1))  =f(PermutationList(xdof1))  +f_element(1);
f(PermutationList(xdof1+1))=f(PermutationList(xdof1+1))+f_element(2);
f(PermutationList(xdof1+2))=f(PermutationList(xdof1+2))+f_element(3);
f(PermutationList(xdof2))  =f(PermutationList(xdof2))  +f_element(4);
f(PermutationList(xdof2+1))=f(PermutationList(xdof2+1))+f_element(5);
f(PermutationList(xdof2+2))=f(PermutationList(xdof2+2))+f_element(6);
f(PermutationList(xdof3))  =f(PermutationList(xdof3))  +f_element(7);
f(PermutationList(xdof3+1))=f(PermutationList(xdof3+1))+f_element(8);
f(PermutationList(xdof3+2))=f(PermutationList(xdof3+2))+f_element(9);
end

```