%
% COMSOL / Matlab script
% Calculation of the acoustic radiation impedance of an unflanged pipe 
% with the FEM (Helmholtz equation)
%
% Antoine Lefebvre, McGill University
%
% Copyright (c) 2010 Antoine Lefebvre
%
flclear fem

rho = 1.25;
c = 343;
a = 0.01;
L = 0.03;
N = 5;

ka = [0.01 colon(0.1,0.1,1.0)];
k = ka/a;

block = block3(2*N*a,2*N*a,2*N*a, 'pos',[0,0,-N*a]);
sphere = sphere3(N*a)*block;
pipe = cylinder3(a,L,'pos', [0,0,-L])*block;
radiation = sphere-pipe;

% Analyzed geometry
clear s
s.objs={radiation,pipe};
s.name={'rad','pipe'};
s.tags={'rad','pipe'};

fem.draw=struct('s',s);
fem.geom=geomgroup(fem,'imprint','off','paircand','none');


edge=0.001;
fem.mesh=meshinit(fem, 'hmax',0.005, 'hmaxfact',0.8, 'hcurve',0.5, ...
                       'hgrad',1.3, 'hcutoff',0.02, 'hnarrow',0.6, ...
                       'hmaxedg', [13,edge,23,edge], 'point',[], ...
                       'edge',[], 'face',[], 'subdomain',[2]);

fem.mesh=meshcopy(fem, 'source',11, 'target',5, 'mcase',0);


fem.mesh=meshinit(fem, 'hmax',0.025, 'hmaxfact',0.8, 'hcurve',0.5, ...
                       'hgrad',1.3, 'hcutoff',0.02, 'hnarrow',0.6, ...
                       'hmaxedg', [13,edge,23,edge], 'point',[], ...
                       'edge',[], 'face',[], 'subdomain',[1], ...
                       'meshstart',fem.mesh);


% Constants
fem.const = {'v0','1 [m/s]', 'a','0.01 [m]', 'S','pi*a^2', 'c','343 [m/s]', ...
             'rho','1.25 [kg/m^3]'};

clear appl
appl.mode.class = 'AcoPressure';
appl.module = 'ACO';
appl.sshape = 3;
appl.assignsuffix = '_acpr';
clear prop
prop.elemdefault='Lag3';
appl.prop=prop;
clear bnd
bnd.p0 = {0,0,1,0};
bnd.type = {'SH','RAD','RAD','NA'};
bnd.wavetype = {'PL','SPH','PL','PL'};
bnd.ind = [4,4,2,1,1,2,1,4,4,3,1,1];
clear pair
pair.type = 'cont';
pair.pair = 'open end';
bnd.pair = pair;
appl.bnd = bnd;
appl.var = {'freq','freq'};
fem.appl{1} = appl;
fem.border = 1;

% Boundary pairs
clear pair
pair{1}.type= 'identity';
pair{1}.name= 'open end';
pair{1}.src.dl = [11];
pair{1}.src.operator = 'src2dst_ip1';
pair{1}.dst.dl = [5];
pair{1}.dst.operator = 'dst2src_ip1';
bnd.pair = pair;

fem.bnd = bnd;
fem.expr = {'freq','ka*c/(2*pi*a)','a0','i*omega_acpr*v0'};
fem=multiphysics(fem);
fem.xmesh=meshextend(fem);
fem.sol=femstatic(fem, 'solcomp',{'p'}, 'outcomp',{'p'}, 'blocksize','auto', ...
		       'pname','ka', 'plist',ka, 'oldcomp',{});

pin = postint(fem,'p/S','Dl',[10],'Edim',[2],'Solnum',[1:length(ka)]);
ZcUin = -postint(fem,'rho*c*nv_acpr/S','Dl',[10],'Edim',[2],'Solnum',[1:length(ka)]);

Zin = pin./ZcUin;
Zr = (Zin.*cos(k*L)-i*sin(k*L))./(-i*Zin.*sin(k*L)+cos(k*L));
flsave('unflanged_pipe_quarter.mph', fem);
save('unflanged_pipe_quarter.mat', 'Zr','a','ka');