Code covered by the BSD License

Accelerating Finite Element Analysis (FEA) in MATLAB

Vaishali Shrivastava (view profile)

Accelerate computationally intensive part of FEA by using Parallel computing.

mems_cantilever_spmode2(xCoords,yCoords,extForces,boundaryCondn, ...
function [perfStats output fea_points fea_tri] = mems_cantilever_spmode2(xCoords,yCoords,extForces,boundaryCondn, ...
elemLength,refineFactor,parmode)
% 0. Assign default arguments if necessary
switch nargin
case 0
errordlg('Missing X and Y coordinates of the geometry, external forces and boundary conditions');
return;
case 1
errordlg('Missing Y coordinates of the geometry, external forces and boundary conditions');
return;
case 2
errordlg('Missing external forces and boundary conditions');
return;
case 3
errordlg('Missing boundary conditions');
return;
case 4
elemLength = 0.1;
refineFactor = 0;
parmode = false;
case 5
refineFactor = 0;
parmode = false;
case 6
parmode = false;
otherwise
end
TETtic = tic;
perfStats.ParalllelMode = parmode;
% 1. Define the domain (geometry, external forces, boundary conditions)
%% 1.1 Geometry Parameters
th = 1;

%% 1.2 Edge nodes (Boundary nodes)
%%% notched beam
bdryNodes = [xCoords yCoords];

%% 1.3 Edges defining boundaries
bdryEdges = [(1:size(bdryNodes,1))' [2:size(bdryNodes,1) 1]'];

%% 1.4 Define external forces on the system
nF = length(extForces); % Number of external point forces applied

% 2. Create the mesh (discretize domain)
%% 2.1 Mesh the domain
hdata.hmax = elemLength;
[nodes tri meshStats] = mesh2d(bdryNodes,bdryEdges,hdata);
for ii = 1:refineFactor
[nodes tri] = refine(nodes,tri);
end
drawnow
close all;
NElems = size(tri,1);
NDOF = 2*size(nodes,1);
disp(sprintf('Completed Mesh Generation. NELems = %i, NDOF = %i',NElems,NDOF))
perfStats.NumberOfElems = NElems;
perfStats.DegreesOfFreedom = NDOF;

% 3. Calculate element level matrices
% 4. Assemble element matrices to global matrices
%% 4.1 Initialize global stiffness and load matrices
K = sparse(NDOF,NDOF);
F = zeros(NDOF,1);
fcn = @plus;

ElemVertices = zeros(NElems,6);
for ie = 1:NElems
lnodes= tri(ie,:);
ElemVertices(ie,:) = [nodes(lnodes(1),:) nodes(lnodes(2),:) nodes(lnodes(3),:)];
end

%% 4.2 Loop through all elements and fill in the global stiffness matrix
if parmode
disp(sprintf('Running in parallel mode'))

% dummy loop to not include parallel setup in tic-toc for KAT
parfor i=1:32
i*2;
end
ATtic = tic;
KATtic = tic;
tempP1 = tic;
parfor ie = 1:NElems
%%% Extract the element vertices
%%% get element stiffness matrix
kelem = get_element_stiffness(ElemVertices(ie,:),th, 170e9, 0.34);
%%% get the global stiffness matrix mapping indices
[rowMap colMap] = get_node_mapping_matrix(tri(ie,:));
K = fcn(K,sparse(rowMap,colMap,kelem,NDOF,NDOF));
end
perfStats.StiffnessMatrixAssemblyTimeInSec = toc(tempP1);
else
disp(sprintf('Running in serial mode'))
ATtic = tic;
KATtic = tic;
for ie = 1:NElems
%%% Extract the element vertices
%%% get element stiffness matrix
kelem = get_element_stiffness(ElemVertices(ie,:),th, 170e9, 0.34);
%%% get the global stiffness matrix mapping indices
[rowMap colMap] = get_node_mapping_matrix(tri(ie,:));
K = fcn(K,sparse(rowMap,colMap,kelem,NDOF,NDOF));
end
perfStats.StiffnessMatrixAssemblyTimeInSec = toc(KATtic);
end
disp(sprintf('Completed Stiffness Matrix (K) Assembly. Time taken:%d sec',...
perfStats.StiffnessMatrixAssemblyTimeInSec))

%% 4.3 Loop through all point forces and fill in the global loading matrix
LATtic = tic;
parfor ip = 1:nF
%%% Identify the element on/within where this point force is being
%%% applied
flag = false;
ie = 1;
while(~flag && ie<=NElems)
%%% Extract the element vertices
A = ElemVertices(ie,1:2);
B = ElemVertices(ie,3:4);
C = ElemVertices(ie,5:6);
P = extForces(ip).location;
ie = ie+1;
end

%%% Error out if the flag is still false after parsing all elements
if (~flag)
error(sprintf(['The point force %i is not applied within' ...
'the domain boundary. Check its location: %i %i'],ip,...
extForces(ip).location(1), extForces(ip).location(2)));
end

ie = ie-1;
%%% fill the load matrix with forces
tempF = zeros(NDOF,1);
for ii = 1:3
idx = 2*tri(ie,ii)-1;
end
F = fcn(F,tempF);
end
disp(sprintf('Completed Load Matrix (F) Assembly. Time taken:%d sec',...

% 5. Apply boundary conditions
%% loop through all boundary condition points and apply the BCs on nearest
%% nodes available
BCTtic = tic;
CONST = full(max(max(K)))*1e6;
for ibc = 1:numel(boundaryCondn)
%% Find the closest node to apply the boundaryCondn
nodeIdx = nearest_node(boundaryCondn(ibc).location,nodes);

%% Modify K and F to include the boundaryCondn
if (~isempty(boundaryCondn(ibc).X))
K(2*nodeIdx-1,2*nodeIdx-1) = K(2*nodeIdx-1,2*nodeIdx-1) + CONST;
F(2*nodeIdx-1) = F(2*nodeIdx-1) + boundaryCondn(ibc).X*CONST;
end
if (~isempty(boundaryCondn(ibc).Y))
K(2*nodeIdx,2*nodeIdx) = K(2*nodeIdx,2*nodeIdx) + CONST;
F(2*nodeIdx) = F(2*nodeIdx) + boundaryCondn(ibc).Y*CONST;
end
end
perfStats.BoundaryConditionsTimeInSec = toc(BCTtic);
disp(sprintf('Applied Boundary Conditions successfully. Time taken:%d sec',...
perfStats.BoundaryConditionsTimeInSec))

tt = toc(ATtic);
perfStats.AssemblyTimeInSec = tt;
disp(sprintf('Completed assembly operations. Time taken:%d sec',...
perfStats.AssemblyTimeInSec))

% 6. Solve the system of equations
perfStats.K_SparsePercent = 100*nnz(K)/(NDOF*NDOF);
kinfo = whos('K');
perfStats.Ksize_bytes = kinfo.bytes;
X = K\F;
disp(sprintf('Solving KX=F complete'))
perfStats.TotalExecutionTimeInSec = toc(TETtic);
disp(sprintf('Total execution time:%d sec',...
perfStats.TotalExecutionTimeInSec))

perfStats.KAT_TET_Ratio = perfStats.AssemblyTimeInSec/perfStats.TotalExecutionTimeInSec;
newNodes = [X(1:2:end)+nodes(:,1) X(2:2:end)+nodes(:,2)];

% 7. Post-processing data
output = X;
fea_points = nodes;
fea_tri = tri;
end

% Calculate element area
function area = get_element_area(A,B,C)
% area of the triangle
area = .5*det([1 A(1) A(2);...
1 B(1) B(2);...
1 C(1) C(2)]);
end

% Calculate element stiffness
function Kelem = get_element_stiffness(elemVer,th, YM, PR)
% YM = 1e6; % Youngs Modulus
% PR = 0.3; % Poisson Ratio
% Stress/Strain matrix for plane-stress condition
% stress = D x strain
D = YM/(1-PR^2).*[1 PR 0;PR 1 0;0 0 (1-PR)/2];

A = elemVer(1:2);
B = elemVer(3:4);
C = elemVer(5:6);

area = get_element_area(A,B,C);

% strain/displacement matrix
% strain = B x disp
B = 1/(2*area).*[B(2)-C(2) 0 C(2)-A(2) 0 A(2)-B(2) 0;...
0 B(1)-C(1) 0 C(1)-A(1) 0 A(1)-B(1);...
B(1)-C(1) B(2)-C(2) C(1)-A(1) C(2)-A(2) A(1)-B(1) A(2)-B(2)];

% element stiffness matrix
Kelem = th*area.*B'*D*B;
end

% Calculate element mass
function Melem = get_element_mass(A,B,C,th)
density = 1e3; % Density of the material
area = get_element_area(A,B,C);
Melem = (density*area*th)/12.*[2 0 1 0 1 0;...
0 2 0 1 0 1;...
1 0 2 0 1 0;...
0 1 0 2 0 1;...
1 0 1 0 2 0;...
0 1 0 1 0 2];
end

Felem = zeros(6,1);
end

% Get local node to global node mapping matrix
function [rowMap colMap] = get_node_mapping_matrix(tri)
rowMap = zeros(1,36);
colMap = zeros(1,36);
kk = 1;
for ii = 1:6
rowIdx = 2*tri(ceil(ii/2)) - mod(ii,2);
for jj = 1:6
colIdx = 2*tri(ceil(jj/2)) - mod(jj,2);
rowMap(kk) = rowIdx;
colMap(kk) = colIdx;
kk = kk+1;
end
end
end
function map = get_node_mapping_matrix_dummy(tri)
map = cell(6);
for ii = 1:6
rowIdx = 2*tri(ceil(ii/2)) - mod(ii,2);
for jj = 1:6
colIdx = 2*tri(ceil(jj/2)) - mod(jj,2);
map{ii,jj} = [rowIdx colIdx];
end
end
end

% Check if a point is inside a triangle
flag = false;

v0 = C-A;
v1 = B-A;
v2 = P-A;

num = dot(v1,v1)*dot(v2,v0) - dot(v1,v0)*dot(v2,v1);
den = dot(v0,v0)*dot(v1,v1) - dot(v0,v1)*dot(v1,v0);
u = num/den;
num = dot(v0,v0)*dot(v2,v1) - dot(v0,v1)*dot(v2,v0);
v = num/den;

au = abs(u);
av = abs(v);
tol = 1e-8;

if (au < tol)
if (av < tol)
flag = true;
return;
elseif (abs(1-av) < tol)
flag = true;
return;
elseif (v < 1)
flag = true;
return;
end
elseif (abs(1-au) < tol)
if (av < tol)
flag = true;
return;
end
elseif (u<1)
if (av < tol)
flag = true;
return;
end
end

if (u>0 && v>0 && abs(1-(u+v))< tol)
flag = true;
return;
end

% Condition for point to be strictly inside the triangle
if ((u>0) && (v>0) && (u+v<1))
flag = true;