%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%?NAME?OF?FILE:?FEM2DL_Box.m
%
%?PURPOSE:?This?Matlab?code?solves?a?two-dimensional?boundary-value?problem
%?using?the?finite?element?method.?The?boundary-value?problem?at?hand
%?corresponds?to?the?Laplace‘s?equation?as?applied?to?a?rectangular?domain?
%?characterized?by?a?simple?medium?（isotropic?linear?and?homogeneous）
%?with?dielectric?constant?epsr.?The?finite?element?method?is?based?on?sub-
%?dividing?the?two-dimensional?domain?into?Ne?linear?triangular?elements.?
%?Dirichlet?boundary?conditions?are?applied?to?all?four?metallic?walls:?
%
%?V=0?on?the?left?sidewall?
%?V=0?on?the?right?sidewall
%?V=0?on?the?bottom?sidewall
%?V=V0?on?the?top?wall?which?is?separated?by?tiny?gaps?from?the?two?sidewalls
%
%?The?dimensions?of?the?domain?are?WxH?where?W=Width?and?H=Hight.?
%?The?user?has?the?freedom?to?modify?any?of?the?defined?input?variables?
%?including?domain?width?and?height?（W?&?L）?top?wall?voltage?（V0）?and?number?
%?of?bricks?along?the?x-?and?y-axes.?Note?that?each?brick?is?then?subdivided?
%?into?two?triangles.
%
%?All?quantities?are?expressed?in?the?SI?system?of?units.
%
%?The?Primary?Unknown?Quantity?is?the?Electric?Potential.
%?The?Secondary?Unknown?Quantity?is?the?Electric?Field.
%
%?Written?by?Anastasis?Polycarpou?（Last?updated:?8/12/2005）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear?%?Clear?all?variables

%?Define?Rectangular?Geometry
%?===========================
W=1;
H=1;

%?Define?the?number?of?bricks?in?the?x?and?y?directions
%?=====================================================
XBRICKS=20;
YBRICKS=20;

%?Define?the?Voltage?（Electric?potential）?on?the?top?wall
%?=======================================================
V0=1;

%?Generate?the?triangular?mesh
%?============================
XLLC=0;??%?x-coordinate?of?the?Lower?Left?Corner
YLLC=0;??%?y-coordinate?of?the?Lower?Left?Corner
XURC=W;??%?x-coordinate?of?the?Upper?Right?Corner
YURC=H;??%?y-coordinate?of?the?Upper?Right?Corner
TBRICKS=XBRICKS*YBRICKS;???%?Total?number?of?bricks
Dx=（XURC-XLLC）/XBRICKS;??%?Edge?length?along?x-direction
Dy=（YURC-YLLC）/YBRICKS;??%?Edge?length?along?y-direction

TNNDS=（XBRICKS+1）*（YBRICKS+1）;????%?Total?number?of?nodes
TNELMS=2*TBRICKS;????%?Total?number?of?triangular?elements
%
%?Print?the?number?of?nodes?and?the?number?of?elements
%?====================================================
fprintf（‘The?number?of?nodes?is:?%6i\n‘TNNDS）
fprintf（‘The?number?of?elements?is:?%6i\n‘TNELMS）

for?I=1:TNNDS
????X（I）=mod（I-1XBRICKS+1）*Dx;
????Y（I）=floor（（I-1）/（XBRICKS+1））*Dy;
end

%?Connectivity?Information
%?========================
for?I=1:TBRICKS
????ELMNOD（2*I-11）=I+floor（（I-1）/XBRICKS）;
????ELMNOD（2*I-12）=ELMNOD（2*I-11）+1+（XBRICKS+1）;
????ELMNOD（2*I-13）=ELMNOD（2*I-11）+1+（XBRICKS+1）-1;
????
????ELMNOD（2*I1）=I+floor（（I-1）/XBRICKS）;
????ELMNOD（2*