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.

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
% E = -grad(potential)
[ex,ey] = pdegrad(nodes,triangles,potential);
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

Contact us