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.

computeMForcesfromV(nodes, edges, triangles, potential, geomtype, xyvalues)
```function perfStats = computeMForcesfromV(nodes, edges, triangles, potential, geomtype, xyvalues)
%% Function to compute forces from electrostatic analysis output from PDE
%% toolbox and do the mechanical analysis to get the deformation

% Compute the magnitude of electric field from potential difference
absE_tri_midpoints = sqrt(ex.^2+ey.^2);

% Interpolate from triangle data midpoints to node data.
absE = pdeprtni(nodes,triangles,absE_tri_midpoints);

% Compute electrostatic pressure
% permittivity of vacuum (8.854*10^-12 farad/meter)
epsilon0 = 8.854e-12;
electricFluxDensity = epsilon0*absE;
ePressure = CalculateElectrostaticPressure(electricFluxDensity);

% Interpolate the data to compute the electrostatic pressure on a uniform
% grid
xmin=min(nodes(1,triangles)); xmax=max(nodes(1,triangles));
ymin=min(nodes(2,triangles)); ymax=max(nodes(2,triangles));

nt=size(triangles,2);
%nxy=ceil(sqrt(nt/2))+1;
nxy=2*ceil(sqrt(nt/2))+1;
x=linspace(xmin,xmax,nxy);
y=linspace(ymin,ymax,nxy);

[xgrid, ygrid] = meshgrid(x,y);

% Electrostatic Pressure on a uniform grid
ePressure_grid = tri2grid(nodes,triangles,ePressure,x,y);
% Potential on a uniform grid
potential_grid = tri2grid(nodes,triangles,potential,x,y);
% Absolute value of Electric field on a uniform grid
absE_grid = tri2grid(nodes,triangles,absE,x,y);
% Electric flux denisty at the grid points on a uniform grid
electricFluxDensity_grid = tri2grid(nodes,triangles,electricFluxDensity,x,y);

zdata = zeros(size(ePressure_grid));

% Plot the data on the uniform grid
figure
subplot(2,2,1);
surf(x,y,zdata,potential_grid);
curAxis = axis;
axis(curAxis(1:4));
title('Potential Difference')

subplot(2,2,2)
surf(x,y,zdata,absE_grid);
axis(curAxis(1:4));
title('Electric Field')

subplot(2,2,3)
surf(x,y,zdata,ePressure_grid);
axis(curAxis(1:4));
title('Electrostatic Pressure')

% Iteration 1
% Calculate loads and boundary conditions for mechanical analysis
[xCoords, yCoords, pForce, BC, nEdges, nodeIdEdges] = CalculateInitialLoadsBCs(geomtype, xyvalues, ...
x, y, xgrid, ygrid, ...
ePressure_grid);

% Perform mechanical analysis
[perfStats output feaPoints feaTriangles] = mems_cantilever_spmode2(xCoords, yCoords, pForce, BC);

% Plot the results
xDisp = output(1:2:end);
yDisp = output(2:2:end);

% Iteration 2
% Calculate subsequent iteration loads and boundary conditions for
% mechanical analysis
[xCoords, yCoords, pForce, BC] = CalculateLoadsBCs(geomtype, feaPoints, feaTriangles, ...
yDisp, electricFluxDensity_grid, xyvalues, ...
nEdges, nodeIdEdges,  xgrid, ygrid, 2);

% Perform 2nd iteration of mechanical analysis
[perfStats output feaPoints2 feaTriangles] = mems_cantilever_spmode2(xCoords, yCoords, pForce, BC);
yDisp2 = output(2:2:end);

% Iteration 3
% Calculate subsequent iteration loads and boundary conditions for
% mechanical analysis
[xCoords, yCoords, pForce, BC] = CalculateLoadsBCs(geomtype, feaPoints2, feaTriangles, ...
yDisp2, electricFluxDensity_grid, xyvalues, ...
nEdges, nodeIdEdges,  xgrid, ygrid, 3);

% Save the electrostatic forces, x y coordinate and boundary conditions for
% the next iteration to a mat file
save('inputdata.mat', 'xCoords', 'yCoords', 'pForce', 'BC');

% Perform 2nd iteration of mechanical analysis
[perfStats output feaPoints3 feaTriangles] = mems_cantilever_spmode2(xCoords', yCoords', pForce, BC, 0.01);
yDisp3 = output(2:2:end);

xfea = feaPoints(:,1);
yfea = feaPoints(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);
xbd1 = xfea(bdId);
ybd1 = yfea(bdId);
yDisp_Bd1 = yDisp(bdId);

xfea = feaPoints2(:,1);
yfea = feaPoints2(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);
xbd2 = xfea(bdId);
ybd2 = yfea(bdId);
yDisp_Bd2 = yDisp2(bdId);

xfea = feaPoints3(:,1);
yfea = feaPoints3(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);
xbd3 = xfea(bdId);
ybd3 = yfea(bdId);
yDisp_Bd3 = yDisp3(bdId);

figure
plot(xyvalues(1:4,1), xyvalues(1:4,2));
hold on
plot(xbd1,ybd1+yDisp_Bd1, 'Color', 'r');
plot(xbd2,ybd2+yDisp_Bd2, 'Color', 'g');
plot(xbd3,ybd3+yDisp_Bd3, 'Color', 'm');
legend({
'Undeformed'
'Iteration 1'
'Iteration 2'
'Iteration 3'
},'Location','NorthWest');

hold off

figure
plot([max(abs(yDisp))*1e6 max(abs(yDisp2))*1e6 max(abs(yDisp3))*1e6]);

end

%%
function ePressure = CalculateElectrostaticPressure(electricFluxDensity)
% Function to compute the electrostatic pressure
% Electrostatic pressure is given by
% P = absD0^2/2*epsilon
% where: D0 = epsilon*E
%        epsilon = relative coefficient of dielectricity * epsilon0
% assuming, relative coefficient of dielectricity = 1

% permittivity of vacuum (8.854*10^-12 farad/meter)
epsilon0 = 8.854e-12;
ePressure = electricFluxDensity.^2/(2*epsilon0);

end

%%
function [xCoords, yCoords, pForce, BC, nEdges, nodeIdEdges] = CalculateInitialLoadsBCs(geomtype, xyvalues, x, y, xgrid, ygrid, ePressure_grid)

% Find the grid points on the geometry
if strcmpi(geomtype, 'rectangle')
%
% (x4,y4) ------------------------ (x3,y3)
%        |                        |
%        |                        |
%        |                        |
% (x1,y1) ------------------------ (x2,y2)
%
x1 = xyvalues(1,1);  y1 = xyvalues(1,2);
x2 = xyvalues(2,1);  y2 = xyvalues(2,2);
x3 = xyvalues(3,1);  y3 = xyvalues(3,2);
x4 = xyvalues(4,1);  y4 = xyvalues(4,2);
dx = x(2)-x(1);
dy = y(2)-y(1);

% Find points closest to the edge connecting (x1,y1) and (x2,y2)
xid = find(xgrid >= (x1-0.5*dx) & xgrid <= (x2+0.5*dx));
yid = find(ygrid >= (y1-dy) & ygrid <= y1);
e1id = intersect(xid, yid);

% Find points closest to the edge connecting (x2,y2) and (x3,y3)
xid = find(xgrid >= x2 & xgrid <= (x2+dx));
yid = find(ygrid >= (y2-0.5*dy) & ygrid <= (y3+0.5*dy));
e2id = intersect(xid, yid);

% Find points closest to the edge connecting (x3,y3) and (x4,y4)
xid = find(xgrid >= (x1-0.5*dx) & xgrid <= (x2+0.5*dx));
yid = find(ygrid >= y3 & ygrid <= (y3+dy));
e3id = intersect(xid, yid);

% Plot the Electrostic pressure at the boundary of the rectangle
subplot(2,2,4)

pForce = [];
% Force on the edge connecting (x1,y1) to (x2,y2)
% unit vector normal to this edge
mag = sqrt((x1-x2)^2+(y2-y1)^2);
f1_dir = [y2-y1 x1-x2]./mag;
% x,y coordinates of distributed load
f1_xCoord = xgrid(e1id(1))+0.5*dx:dx:xgrid(e1id(end));
f1_yCoord = y1*ones(size(f1_xCoord));
for idx=1:length(e1id)-1
f1_vals(idx) = 0.5*(ePressure_grid(e1id(idx))+ePressure_grid(e1id(idx+1)))*dx;  %*1e-6; % Assuming unit micron thickness of the beam
end
for idx=1:length(f1_xCoord)
pForce(idx).location = [f1_xCoord(idx) f1_yCoord(idx)];
pForce(idx).compX = f1_vals(idx)*f1_dir(1);
pForce(idx).compY = f1_vals(idx)*f1_dir(2);
end

% Force on the edge connecting (x2,y2) to (x3,y3)
% unit vector normal to this edge
mag = sqrt((x3-x2)^2+(y3-y2)^2);
f2_dir = [y3-y2 x2-x3]./mag;
% x,y coordinates of distributed load
f2_yCoord = ygrid(e2id(1))+0.5*dy:dy:ygrid(e2id(end));
f2_xCoord = x2*ones(size(f2_yCoord));
for idx=1:length(e2id)-1
f2_vals(idx) = 0.5*(ePressure_grid(e2id(idx))+ePressure_grid(e2id(idx+1)))*dy; %*1e-6; % Assuming unit micron thickness of the beam
end

% Force on the edge connecting (x3,y3) to (x4,y4)
% unit vector normal to this edge
mag = sqrt((x4-x3)^2+(y4-y3)^2);
f3_dir = [y4-y3 x3-x4]./mag;
% x,y coordinates of distributed load
f3_xCoord = xgrid(e3id(1))+0.5*dx:dx:xgrid(e3id(end));
f3_yCoord = y3*ones(size(f3_xCoord));
for idx=1:length(e3id)-1
f3_vals(idx) = 0.5*(ePressure_grid(e3id(idx))+ePressure_grid(e3id(idx+1)))*dx;  %*1e-6; % Assuming unit micron thickness of the beam
end
lp = length(pForce);
for idx=1:length(f3_xCoord)
pForce(lp+idx).location = [f3_xCoord(idx) f3_yCoord(idx)];
pForce(lp+idx).compX = f3_vals(idx)*f3_dir(1);
pForce(lp+idx).compY = f3_vals(idx)*f3_dir(2);
end

% Plot the Electrostatic forces on the boundary
plot3(f1_xCoord,f1_yCoord,f1_vals, 'Color', 'g', 'Marker', '*')
hold on
plot3(f2_xCoord,f2_yCoord,f2_vals, 'Color', 'g', 'Marker', '*')
plot3(f3_xCoord,f3_yCoord,f3_vals, 'Color', 'g', 'Marker', '*')
hold off

% Boundary conditions
BC(1).location = [x1 y1];
BC(1).X = 0;
BC(1).Y = 0;

BC(2).location = [x4 y4];
BC(2).X = 0;
BC(2).Y = 0;

% XY coordinates of the geometry
xCoords = xyvalues(1:4,1);
yCoords = xyvalues(1:4,2);

% Number of free edges
nEdges = 3;
% Node IDs on free edges
nodeIdEdges(1).id = e1id;
nodeIdEdges(2).id = e2id;
nodeIdEdges(3).id = e3id;
end
end

%%
function [xCoords, yCoords, pForce, BC] = CalculateLoadsBCs(geomtype, feaPoints, feaTriangles, yDisp, electricFluxDensity_grid, xyvalues, nEdges, nodeIdEdges, xgrid, ygrid, itr)

xCoords = [];
yCoords = [];
pForce = [];
BC = [];
if strcmpi(geomtype, 'rectangle')
e1id = nodeIdEdges(1).id;
e3id = nodeIdEdges(3).id;

% Find the points of the FEA mesh that lie on the boundary
xfea = feaPoints(:,1);
yfea = feaPoints(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);

% Coordinates of FEA points on the boundary
xbd = xfea(bdId);
ybd = yfea(bdId);
yDisp_Bd = yDisp(bdId);

% Find the FEA nodes close to rectangular grid
xgrid1 = xgrid(e1id);
xgrid3 = xgrid(e3id);

dx = 0.5*(abs(xgrid1(2)-xgrid1(1)));

yDisp_grid = [];
xbd_grid = [];
ybd_grid = [];
for idx=1:length(xgrid1)
nx = find(xbd >= (xgrid1(idx)-dx) & xbd <= (xgrid1(idx)+dx));
ny = find(ybd >= (xyvalues(1,2)-0.01*dx) & ybd <= (xyvalues(1,2)+0.01*dx));
n1 = intersect(nx, ny);
n11 = nearestNeighbor(dt, [xgrid1(idx), xyvalues(1,2)]);
mindist = sqrt((xgrid1(idx)-dt.X(n11,1))^2+(xyvalues(1,2)-dt.X(n11,2))^2);
xNearest = dt.X(n11,1);
yNearest = dt.X(n11,2);
yDispNearest = yDisp(n11);
for jdx=1:length(n1)
dist = sqrt((xgrid1(idx)-xbd(n1(jdx)))^2+(xyvalues(1,2)-ybd(n1(jdx)))^2);
if dist < mindist
mindist = dist;
xNearest = xbd(n1(jdx));
yNearest = ybd(n1(jdx));
yDispNearest = yDisp_Bd(n1(jdx));
end
end
yDisp_grid(idx,1) = yDispNearest; %sum(yDisp_Bd(n1))/length(n1);
xbd_grid(end+1) = xNearest;
ybd_grid(end+1) = yNearest;
end

for idx=1:length(xgrid3)
nx = find(dt.X(:,1) >= (xgrid3(idx)-dx) & dt.X(:,1) <= (xgrid3(idx)+dx));
ny = find(dt.X(:,2) >= (xyvalues(3,2)-0.01*dx) & dt.X(:,2) <= (xyvalues(3,2)+0.01*dx));
n3 = intersect(nx, ny);
n33 = nearestNeighbor(dt, [xgrid3(idx), xyvalues(3,2)]);
mindist = sqrt((xgrid3(idx)-dt.X(n33,1))^2+(xyvalues(3,2)-dt.X(n33,2))^2);
xNearest = dt.X(n33,1);
yNearest = dt.X(n33,2);
yDispNearest = yDisp(n33);
for jdx=1:length(n3)
dist = sqrt((xgrid3(idx)-dt.X(n3(jdx),1))^2+(xyvalues(3,2)-dt.X(n3(jdx),2))^2);
if dist < mindist
mindist = dist;
xNearest = dt.X(n3(jdx),1);
yNearest = dt.X(n3(jdx),2);
yDispNearest = yDisp(n3(jdx));
end
end
yDisp_grid(idx,2) = yDispNearest;  %sum(yDisp_Bd(n3))/length(n3);
xbd_grid(end+1) = xNearest;
ybd_grid(end+1) = yNearest;
end

% Assemble D0 (electric flux density) at the boundary
D0 = [electricFluxDensity_grid(e1id) electricFluxDensity_grid(e3id)];

% New electrostatic pressure at the boundary
% Ddef = D0*(Gap/Gap-y displacement)
Gap = xyvalues(5,1)*ones(size(yDisp_grid));
Ddef = D0.*(Gap./(Gap-yDisp_grid));
ePressure_Def = CalculateElectrostaticPressure(Ddef);

% XY coordinates of the geometry
if length(xbd) < length(xgrid1)
xCoords = [xbd_grid(1:length(xgrid1)) xbd_grid(end:-1:length(xgrid1)+1)];
ytempCoords = ybd_grid + [yDisp_grid(:,1)' yDisp_grid(:,2)'];
yCoords = [ytempCoords(1:length(xgrid1)) ytempCoords(end:-1:length(xgrid1)+1)];
if (xCoords(end) ~= xCoords(1) || yCoords(end) ~= yCoords(1))
xCoords(end+1) = xCoords(1);
yCoords(end+1) = yCoords(1);
end
else
xCoords = xbd;
yCoords = ybd + yDisp_Bd;
end
figure
plot(xCoords, yCoords);
hold on

% Compute forces
% Force on the edge connecting (x1,y1) to (x2,y2)
% x coordinate of the distributed load
f1_xCoord = xgrid1(1)+dx:2*dx:xgrid1(end);
ytemp = xyvalues(1,2)*ones(1,length(e1id))+yDisp_grid(:,1)';
for idx=1:length(e1id)-1
f1_vals(idx) = 0.5*(ePressure_Def(idx,1)+ePressure_Def(idx,1))*2*dx;  %*1e-6; % Assuming unit micron thickness of the beam

nx = find(xCoords >= (f1_xCoord(idx)-2*dx) & xCoords <= (f1_xCoord(idx)+2*dx));
ny = find(yCoords >= (ytemp(idx)-0.02*dx) & yCoords <= (ytemp(idx)+0.02*dx));
n1 = intersect(nx, ny);

if length(n1)>1
xdist = abs(xCoords(n1)-f1_xCoord(idx)*ones(size(xCoords(n1))));
[xdist ix] = sort(xdist);
f1_yCoord(idx) = max([yCoords(n1(ix(2))) yCoords(n1(ix(1)))]);
else
f1_yCoord(idx) = yCoords(n1);
end

% unit vector normal to this edge
mag = sqrt((xgrid1(idx)-xgrid1(idx+1))^2+(ytemp(idx+1)-ytemp(idx))^2);
f1_dir = [ytemp(idx+1)-ytemp(idx) xgrid1(idx)-xgrid1(idx+1)]./mag;

pForce(idx).location = [f1_xCoord(idx) f1_yCoord(idx)];
pForce(idx).compX = f1_vals(idx)*f1_dir(1);
pForce(idx).compY = f1_vals(idx)*f1_dir(2);
end

% Force on the edge connecting (x3,y3) to (x4,y4)
% x coordinate of the distributed load
lp = length(pForce);
f3_xCoord = xgrid3(1)+dx:2*dx:xgrid3(end);
ytemp = xyvalues(3,2)*ones(1,length(e3id))+yDisp_grid(:,2)';
for idx=1:length(e3id)-1
f3_vals(idx) = 0.5*(ePressure_Def(idx,2)+ePressure_Def(idx,2))*2*dx;  %*1e-6; % Assuming unit micron thickness of the beam

nx = find(xCoords >= (f3_xCoord(idx)-2*dx) & xCoords <= (f3_xCoord(idx)+2*dx));
ny = find(yCoords >= (ytemp(idx)-0.05*dx) & yCoords <= (ytemp(idx)+0.05*dx));
n1 = intersect(nx, ny);

if length(n1)>1
xdist = abs(xCoords(n1)-f3_xCoord(idx)*ones(size(xCoords(n1))));
[xdist ix] = sort(xdist);
f3_yCoord(idx) = min([yCoords(n1(ix(2))) yCoords(n1(ix(1)))]);
else
f3_yCoord(idx) = yCoords(n1);
end

% unit vector normal to this edge
mag = sqrt((xgrid3(idx)-xgrid3(idx+1))^2+(ytemp(idx+1)-ytemp(idx))^2);
f3_dir = [ytemp(idx+1)-ytemp(idx) xgrid3(idx)-xgrid3(idx+1)]./mag;

pForce(lp+idx).location = [f3_xCoord(idx) f3_yCoord(idx)];
pForce(lp+idx).compX = f3_vals(idx)*f3_dir(1);
pForce(lp+idx).compY = f3_vals(idx)*f3_dir(2);
end

plot(f1_xCoord,f1_yCoord, '.', 'markersize', 10);
plot(f3_xCoord,f3_yCoord, '.', 'markersize', 10);
title(['Deformed shape - Iteration ' num2str(itr)]);
hold off

% Boundary conditions
BC(1).location = [xyvalues(1,1) xyvalues(1,2)];
BC(1).X = 0;
BC(1).Y = 0;

BC(2).location = [xyvalues(4,1) xyvalues(4,2)];
BC(2).X = 0;
BC(2).Y = 0;

end
end
```