Question

Solve 2-dim Helmholtz equation in an infinitely long perfect conductor waveguide, with dimensions a > b

a) Solve for the analytical form using separation of variable techniques

b) Plot the first 6 base functions (m, n = 0, 1, 2, with m and n not equal to zero simultaneously). Attempt some 3d plots (use MATLAB)

c) Find their frequencies

Answer 1

Dhf

God

gf_workspace('clear all');

disp('2D scalar wave equation (helmholtz) demonstration');

disp(' we present three approaches for the solution of the helmholtz problem')

disp(' - the first one is to use the new getfem "model bricks"')

disp(' - the second one is to use the old getfem "model bricks"')

disp(' - the third one is to use the "low level" approach, i.e. to assemble')

disp(' and solve the linear systems.')

disp('The result is the wave scattered by a disc, the incoming wave beeing a plane wave coming from the top');

disp(' \delta u + k^2 = 0');

disp(' u = -uinc on the interior boundary');

disp(' \partial_n u + iku = 0 on the exterior boundary');

%PK = 10; gt_order = 6; k = 7; use_hierarchical = 0; load_the_mesh=0;

PK=3; gt_order = 3; k = 1; use_hierarchical = 1; load_the_mesh=1;

if (use_hierarchical) s = 'hierarchical'; else s = 'classical'; end;

disp(sprintf('using %s P%d FEM with geometric transformations of degree %d',s,PK,gt_order));

if (load_the_mesh),

disp('the mesh is loaded from a file, gt_order ignored');

end;

if load_the_mesh == 0,

% a quadrangular mesh is generated, with a high degree geometric transformation

% number of cells for the regular mesh

Nt=10; Nr=8;

m=gfMesh('empty',2);

dtheta=2*pi*1/Nt; R=1+9*(0:Nr-1)/(Nr-1);

gt=gfGeoTrans(sprintf('GT_PRODUCT(GT_PK(1,%d),GT_PK(1,1))',gt_order));

ddtheta=dtheta/gt_order;

for i=1:Nt;

for j=1:Nr-1;

ti=(i-1)*dtheta:ddtheta:i*dtheta;

X = [R(j)*cos(ti) R(j+1)*cos(ti)];

Y = [R(j)*sin(ti) R(j+1)*sin(ti)];

m.set('add convex',gt,[X;Y]);

end;

end;

fem_u=gfFem(sprintf('FEM_QK(2,%d)',PK));

fem_d=gfFem(sprintf('FEM_QK(2,%d)',PK));

mfu=gfMeshFem(m,1);

mfd=gfMeshFem(m,1);  

mfu.set('fem',fem_u);

mfd.set('fem',fem_d);

sIM=sprintf('IM_GAUSS_PARALLELEPIPED(2,%d)',gt_order+2*PK);

mim=gfMeshIm(m, gfInteg(sIM));

else

% the mesh is loaded

m=gfMesh('import','gid','../meshes/holed_disc_with_quadratic_2D_triangles.msh');

if (use_hierarchical),

% hierarchical basis improve the condition number

% of the final linear system

fem_u=gfFem(sprintf('FEM_PK_HIERARCHICAL(2,%d)',PK));

%fem_u=gfFem('FEM_HCT_TRIANGLE');

%fem_u=gfFem('FEM_HERMITE(2)');

else,

fem_u=gfFem(sprintf('FEM_PK(2,%d)',PK));

end;

fem_d=gfFem(sprintf('FEM_PK(2,%d)',PK));

mfu=gfMeshFem(m,1);

mfd=gfMeshFem(m,1);  

set(mfu,'fem',fem_u);

set(mfd,'fem',fem_d);

mim=gfMeshIm(m,gfInteg('IM_TRIANGLE(13)'));

end;

nbdu=mfu.nbdof;

nbdd=mfd.nbdof;

% identify the inner and outer boundaries

P=m.pts; % get list of mesh points coordinates

pidobj=find(sum(P.^2) < 1*1+1e-6);

pidout=find(sum(P.^2) > 10*10-1e-2);

% build the list of faces from the list of points

fobj=get(m,'faces from pid',pidobj);

fout=get(m,'faces from pid',pidout);

set(m,'boundary',1,fobj);

set(m,'boundary',2,fout);

% expression of the incoming wave

wave_expr=sprintf('cos(%f*y+.2)+1i*sin(%f*y+.2)',k,k);

Uinc=get(mfd,'eval',{wave_expr});

%

% we present three approaches for the solution of the Helmholtz problem

% - the first one is to use the new getfem "model bricks"

% - the second one is to use the old getfem "model bricks"

% - the third one is to use the "low level" approach, i.e. to assemble

% and solve the linear systems.

if 1,

t0=cputime;

% solution using new model bricks

md=gf_model('complex');

gf_model_set(md, 'add fem variable', 'u', mfu);

gf_model_set(md, 'add initialized data', 'k', [k]);

gf_model_set(md, 'add Helmholtz brick', mim, 'u', 'k');

gf_model_set(md, 'add initialized data', 'Q', [1i*k]);

gf_model_set(md, 'add Fourier Robin brick', mim, 'u', 'Q', 2);

gf_model_set(md, 'add initialized fem data', 'DirichletData', mfd, Uinc);

gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', mfd, 1, 'DirichletData');

% gf_model_set(md, 'add Dirichlet condition with penalization', mim, 'u', 1e12, 1, 'DirichletData');

gf_model_get(md, 'solve');

U = gf_model_get(md, 'variable', 'u');

disp(sprintf('solve do

Fj