%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finite Element MATLAB script for
% deflection and stress analysis of beams
% Written by: G.K. Ananthasuresh
% for MEAM 550 class in Fall 1998 to solve
% coupled electrostatic-elastostatic problem
% using the relaxation method.
% Modified in Fall 2001, Spring 2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf						% Clear graphics window
clear all				% Clear all variables
clf						% Clear command window
hold off				% No hold on the graphics window
%
% This script needs the following scripts to run
%	matcut.m	--> trims a matrix
%	veccut.m	--> trims a vector
%	femtruss.m	--> FEA for trusses
% It also needs the following input files	
%	node.dat	--> 	nodal data
%	elem.dat	-->	element data
%	disp.dat	-->	displacement boundary conditon data
%
% Enter data
length_in = 300; % in micron units
width = 15; % width of the beam in micron units
thickness = 1.0; % thickness of the beam in micron units
Y_in = 155E3; % Young's modulus of 155 GPa
% Initial gap in microns
g0 = 1.2;
% Pull-in voltage with the 1-D formula
Area = length_in*width;
k = 32*Y_in*width*thickness^3/(length_in)^3;
Vpi = sqrt( 8*k*g0^3/27/8.854e-30/Area/(1+0.65*g0/width))*1e-12
%VpiBeam = sqrt( 0.3475*k*g0^3/8.854e-30/Area*(1-0.0501*g0/thickness))*1e-12
alpha = 11.2782;
VpiBeam = sqrt( alpha*Y_in*thickness^3*g0^3/8.854e-30/length_in^4)*1e-12
%
% a = Y_in*thickness^3*g0^3
% b = 8.854e-30*length_in^4*g0/width*(11.15e12)^2
% c = 8.854e-30*length_in^4*(11.15e12)^2
a1 = 155e3*1.0^3*1.0^3;
b1 = -8.854e-30*200^4*1.0/15.0*(11.15e12)^2;
c1 = 8.854e-30*200^4*(11.15e12)^2;
a2 = 155e3*1.4^3*1.0^3;
b2 = -8.854e-30*320^4*1.0/25.0*(7.25e12)^2;
c2 = 8.854e-30*320^4*(7.25e12)^2;
alpha = (b2*c1-b1*c2)/(a1*b2-a2*b1);
beta = (c1-a1*alpha)/b1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ INPUT from files 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Included global variables from file globals.m
%
% Read in nodal and element connectivity data from
% files: nodes.dat and elem.dat
%load node.dat
nodenos = 65; % Number of nodes (do not change this!)
node(:,1) = [1:1:nodenos]';
node(:,2) = [0:length_in/(nodenos-1):length_in]';
node(:,3) = zeros(nodenos,1);
%load elem.dat
elem(:,1) = [1:1:nodenos-1]';
elem(:,2) = [1:1:nodenos-1]';
elem(:,3) = [2:1:nodenos]';
elem(:,4) = Y_in*ones(nodenos-1,1);
elem(:,5) = width*ones(nodenos-1,1);
elem(:,6) = thickness*ones(nodenos-1,1);
%
% Read in displacement boundary condition data from the file disp.dat
load dispbc.dat
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRE-PROCESSING 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify the number of nodes, X and Y Coordinates of the nodes
NNODE = size(node,1);
nx = node(:,2); 
ny = node(:,3);
%
% Identify the number of elements and form an element connectivity array,
% the cross-section and Young's modulus arrays.
NELEM = size(elem,1);
ncon = elem(:,[2 3]);
E = elem(:,4);
width = elem(:,5);
t = elem(:,6);
A = width.*t; %cross-section area of the beam
Inertia= width.*t.^3/12; %Area moment of inertia of the beam
%
% Displacement boundary conditions
Nfix = size(dispbc,1);
j = 0;
for i = 1:Nfix,
    j = j + 1;
    dispID(j) = (dispbc(i,2)-1)*3+dispbc(i,3);
end
[dispID sortIndex] = sort(dispID);
%
% Compute the lengths of the elements
for ie=1:NELEM,
    eye = ncon(ie,1);
    jay = ncon(ie,2);
    L(ie) = sqrt ( (nx(jay) - nx(eye))^2 + (ny(jay) - ny(eye))^2 );
end
%
% Specify input voltage here
V = input('Enter voltage:');

% Arrange force information into a force vector, F
F = zeros(3*NNODE,1);			% Initialization

U = zeros(3*NNODE,1);
U_prev = U;
% Initialization
max_iter = 50;
delta_U = 100.0;    % Initialized to a large value
% Force per unit gap per unit gap (in micro N)
fV = 8.854e-12 * V^2/2*1E6;    % Force per unit area per
% unit gap. (Note conversion factor of 1E6 for "micron" units).
%
%Relaxation iterations begin here.
for relax_iteration = 1:max_iter, 
    F = zeros(3*NNODE,1);	% Initialization to zero.
    if abs(delta_U) < 1e-4, 
        break
    end
    % 
    % Computation of electrostatic forces on nodes
    %
    for i = 5:1:NELEM-4,
        gi = g0 - U( 3*(i-1) + 2);
        gj = g0 - U( 3*i + 2 );
        gave = (gi+gj)/2;
        eforce = fV * width(i) * L(i) / gave^2;
        eforce = eforce * (1 + 0.65*gave/width(i));
        F( 3*(i-1) + 2 ) = F( 3*(i-1) + 2 ) + eforce/2;
        F( 3*i + 2 ) = F( 3*i + 2 ) + eforce/2;
    end
    %
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % SOLUTION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Call fembeam.m to solve for the following.
    %	Deflections: U
    %	Reaction forces at the constrained nodes: R
    %
    [U,R] = fembeam(A,L,E,nx,ny,ncon,NELEM,NNODE,F,dispID,Inertia);
    %
    
    delta_U = max( abs(U(2:3:3*NNODE)) ) - max( abs(U_prev(2:3:3*NNODE)) );
    U_prev = U;
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % POST-PROCESSING
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Plotting
    %clf
    for ip=1:NELEM,
        pt1 = ncon(ip,1); pt2 = ncon(ip,2);
        dx1 = U(3*(pt1-1)+1);
        dy1 = U(3*(pt1-1)+2);
        dx2 = U(3*(pt2-1)+1);
        dy2 = U(3*(pt2-1)+2);
        % 
        plot([nx(pt1) nx(pt2)], [ny(pt1) ny(pt2)],'-y');
        hold on
        plot([nx(pt1)+dx1 nx(pt2)+dx2], [ny(pt1)-dy1 ny(pt2)-dy2],'-b');
        hold on  
    end
    relax_iteration
    %delta_U
    %max( abs(U(2:3:3*NNODE)) )
    %pause
end
% Relaxation iterations end here.
echo on
% Converged in iterations...
relax_iteration - 1
echo off
% Plot the final solution in a different color. 
for ip=1:NELEM,
    pt1 = ncon(ip,1); pt2 = ncon(ip,2);
    dx1 = U(3*(pt1-1)+1);
    dy1 = U(3*(pt1-1)+2);
    dx2 = U(3*(pt2-1)+1);
    dy2 = U(3*(pt2-1)+2);
    % 
    plot([nx(pt1)+dx1 nx(pt2)+dx2], [ny(pt1)-dy1 ny(pt2)-dy2],'-r');
    hold on  
end
xlabel('Length along the beam (um)');
ylabel('downward Deflection (um)');
%axis('equal');
%legend('Undeformed','Deformed');
grid
hndl = plot([0 320],[-1 -1],'-r');
set(hndl,'LineWidth',3);
title('Ground electrode is shown in thick red line at the bottom');






