%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% NAME OF FILE: FEM2DL_Cyl.m
%
% PURPOSE: This Matlab code solves a two-dimensional scattering problem
% using the finite element method. The scatterer is a perfectly
% conducting circular cylinder of radius robj; all the dimensions are given
% in terms of the free-space wavelength. A TM-to-z incident plane wave is
% scattered from the circular cylinder and propagates undisturbed toward
% infinity. To simulate this undisturbed propagation of the scattered field
% toward infinity, a first-order absorbing boundary condition (ABC) was
% imposed at a distance rho away from the center of the cylinder. The farther
% the ABC boundary is placed, the more accurate the ABC is. The free space
% between the circular cylinder and the ABC boundary is subdivided into
% triangles governed by linear interpolation functions. Dirichlet boundary
% conditions are applied on the surface of the cylinder; i.e., the tangential
% electric field, in this case Ez, is set to zero for all the nodes that
% coincide with the surface of the circular cylinder.
%
% The finite element method is applied to the homogeneous scalar wave
% equation, otherwise known as the homogeneous Helmholtz equation. The
% primary unknown quantity is the total electric field in the z-direction
% which is given by the incident field plus the scattered field. The
% direction of the incident field is set toward the positive x-axis (phi_i=0)
% whereas the total field is evaluated at a distance half-way between the
% scatterer and the ABC boundary for all observation angles between 0 and 360
% degrees. The numerical solution is compared with the exact analytical
% solution.
%
% The user is allowed to set the following input parameters:
%
% rhoj = radius of the scatterer (circular cylinder) in wavelengths
% rho = radius of the ABC boundary in wavelengths
% h = discritization size in wavelengths
% E0 = amplitude of the incident electric field
%
% IMPORTANT: Depending on the number of nodes in the domain and the clock
% speed of your computer, the finite element code may take several minutes
% to execute. Try not to exceed 5,000 nodes otherwise you have to wait a
% significant amount of time to get results!
%
% Written by Anastasis Polycarpou (Last updated: 8/12/2005)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear % Clear all variables
%
% Define the Input parameters
% ===========================
robj=0.5; % Set the radius of the circular PEC cylinder in wavelengths
rho=1.5; % Set the radius of the outer circular ABC boundary in wavelengths
h=0.04; % Set the discretization size in wavelengths
E0=1; % Set the amplitude of the Incident field
% ===========================
%
% Determine the number of annular rings
% =====================================
Nseg=ceil((rho-robj)/(sqrt(3)*h/2));
dr=(rho-robj)/Nseg;
%
% Create the nodes on each annular ring
% =====================================
Nnods=0;
for J=1:Nseg+1
r=robj+(J-1)*dr;
Npoints=ceil(2*pi*r/h);
dphi=2*pi/Npoints;
for JJ=1:Npoints
Nnods=Nnods+1;
phi=(JJ-1)*dphi;
x(Nnods)=r*cos(phi);
y(Nnods)=r*sin(phi);
end
end
%
% Triangulate using Delaunay triangulation
% ========================================
tri=delaunay(x,y);
%
% Eliminate the triangles within the perfectly conducting circular cylinder
% =========================================================================
Nelms=0;
for J=1:length(tri)
Innods=0;
for JJ=1:3
r0=sqrt(x(tri(J,JJ))^2+y(tri(J,JJ))^2);
if abs(robj-r0) < 0.0001*robj
Innods=Innods+1;
end
end
if Innods <= 2
Nelms=Nelms+1;
elmnod(Nelms,:)=tri(J,:);
end
end
%
% Print the number of nodes and the number of elements
% ====================================================
fprintf('The number of nodes is: %6i\n',Nnods)
fprintf('The number of elements is: %6i\n',Nelms)
%
% Set all triangles in a counter-clockwise sense of numbering
% ===========================================================
for e=1:Nelms
x21=x(elmnod(e,2))-x(elmnod(e,1));
x31=x(elmnod(e,3))-x(elmnod(e,1));
y21=y(elmnod(e,2))-y(elmnod(e,1));
y31=y(elmnod(e,3))-y(elmnod(e,1));
Ae=0.5*(x21*y31-x31*y21);
if Ae < 0
temp=elmnod(e,2);
elmnod(e,2)=elmnod(e,3);
elmnod(e,3)=temp;
end
end
%
% Plot the mesh
% =============
figure(1)
triplot(elmnod,x,y);
xlabel('x (wavelengths)');
ylabel('y (wavelengths)');
axis([-rho rho -rho rho]);
axis square;
%
% Set the Dirichlet & Absorbing BCs on the cylinder and outer boundary,
% respectively
% =====================================================================
NEBC=0;
NABC=0;
for JJ=1:Nnods
r0=sqrt(x(JJ)^2+y(JJ)^2);
if abs(robj-r0) < 0.0001*robj
NEBC=NEBC+1;
EBC(NEBC)=JJ;
elseif abs(rho-r0) < 0.0001*rho
NABC=NABC+1;
ABC(NABC)=JJ;
end
end
for e=1:Nelms
for I=1:3
Nd=elmnod(e,I);
flag(e,I)=0;
for J=1:NABC
if Nd == ABC(J)
flag(e,I)=1;
end
end
end
end
%
% Print the number of Dirichlet BCs and the number of ABCs
% ========================================================
fprintf('The number of nodes on a Dirichlet boundary is: %6i\n',NEBC)
fprintf('The number of nodes on an ABC boundary is: %6i\n',NABC)
%
% Assign the constants for each element
% =====================================
mur0=1;
epsr0=1;
k0=2*pi;
gamma=k0*j+1/(2*rho);
for e=1:Nelms
alphax(e)=1/mur0;
alphay(e)=1/mur0;
beta(e)=k0^2*epsr0;
g(e)=0;
end
%
% Initialization of the global K matrix and right-hand side vector
% ================================================================
K=sparse(Nnods,Nnods);
b=zeros(Nnods,1);
%
% Form the element matrices and assemble to the global matrix
% ===========================================================
for e=1:Nelms
x21=x(elmnod(e,2))-x(elmnod(e,1));
x31=x(elmnod(e,3))-x(elmnod(e,1));
x32=x(elmnod(e,3))-x(elmnod(e,2));
x13=x(elmnod(e,1))-x(elmnod(e,3));
y12=y(elmnod(e,1))-y(elmnod(e,2));
y21=y(elmnod(e,2))-y(elmnod(e,1));
y31=y(elmnod(e,3))-y(elmnod(e,1));
y23=y(elmnod(e,2))-y(elmnod(e,3));
Ae=0.5*(x21*y31-x31*y21);
% Evaluation of the element K matrix
% ==================================
Me(1,1)=-(alphax(e)*y23^2+alphay(e)*x32^2)/(4*Ae);
Me(1,2)=-(alphax(e)*y23*y31+alphay(e)*x32*x13)/(4*Ae);
Me(2,1)=Me(1,2);
Me(1,3)=-(alphax(e)*y23*y12+alphay(e)*x32*x21)/(4*Ae);
Me(3,1)=Me(1,3);
Me(2,2)=-(alphax(e)*y31^2+alphay(e)*x13^2)/(4*Ae);
Me(2,3)=-(alphax(e)*y31*y12+alphay(e)*x13*x21)/(4*Ae);
Me(3,2)=Me(2,3);
Me(3,3)=-(alphax(e)*y12^2+alphay(e)*x21^2)/(4*Ae);
% Evaluation of the element T matrix
% ==================================
Te(1,1)=beta(e)*Ae/6;
Te(1,2)=beta(e)*Ae/12;
Te(2,1)=beta(e)*Ae/12;
Te(1,3)=beta(e)*Ae/12;
Te(3,1)=beta(e)*Ae/12;
Te(2,2)=beta(e)*Ae/6;
Te(2,3)=beta(e)*Ae/12;
Te(3,2)=beta(e)*Ae/12;
Te(3,3)=beta(e)*Ae/6;
% Sum the element matrices Me and Te
% ==================================
for I=1:3
for J=1:3
Ke(I,J)=Me(I,J)+Te(I,J);
end
end
% Evaluation of element vector ge
% ===============================
ge(1)=g(e)*Ae/3;
ge(2)=g(e)*Ae/3;
ge(3)=g(e)*Ae/3;
% Evaluation of the element vector pe & update of the element K-matrix
% ====================================================================
pe=zeros(3,1);
if flag(e,1) == 1 && flag(e,2) == 1
x1=x(elmnod(e,1));
y1=y(elmnod(e,1));
x2=x(elmnod(e,2));
