Code covered by the BSD License  

Highlights from
Accelerating Finite Element Analysis (FEA) in MATLAB

image thumbnail

Accelerating Finite Element Analysis (FEA) in MATLAB

by

 

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;
        [flag loadSplit] = is_point_in_triangle(P,A,B,C);
        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;
        tempF(idx) = tempF(idx) + extForces(ip).compX*loadSplit(ii);
        tempF(idx+1) = tempF(idx+1) + extForces(ip).compY*loadSplit(ii);
    end
    F = fcn(F,tempF);
end
perfStats.LoadMatrixAssemblyTimeInSec = toc(LATtic);
disp(sprintf('Completed Load Matrix (F) Assembly. Time taken:%d sec',...
    perfStats.LoadMatrixAssemblyTimeInSec))

% 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

% Calculate element load matrix
function Felem = get_element_loading(A,B,C,th)
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
function [flag loadSplit] = is_point_in_triangle(P,A,B,C)
flag = false;
loadSplit = [0 0 0];

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;
        loadSplit = [1 0 0];
        return;
    elseif (abs(1-av) < tol)
        flag = true;
        loadSplit = [0 1 0];
        return;
    elseif (v < 1)
        flag = true;
        loadSplit = [0.5 0.5 0];
        return;
    end
elseif (abs(1-au) < tol)
    if (av < tol)
        flag = true;
        loadSplit = [0 0 1];
        return;
    end
elseif (u<1)
    if (av < tol)
        flag = true;
        loadSplit = [0.5 0 0.5];
        return;
    end
end

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

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

function nodeIdx = nearest_node(p,nodes)
P = [p(1)*ones(size(nodes,1),1) p(2)*ones(size(nodes,1),1)];
D = P-nodes;
[minDist nodeIdx] = min(sum(D.*D,2));
end

Contact us