% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 494;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/09/19 16:09:48 $';
fem.version = vrsn;

% Geometry
g1=rect2(0.24,0.24,'base','center','pos',[-0.13,0.07]);
g2=rect2('5.3/1000','5.3/1000','base','center','pos',{'0','0'},'rot','0');
g3=ellip2(5.5E-4,3.0E-4,'base','corner','pos',[-8.0E-4,0]);
g4=ellip2('0.8/1000','0.8/1000','base','center','pos',{'0','0'},'rot','0');
g5=geomcomp({g2,g4},'ns',{'R1','E1'},'sf','R1-E1','edge','none');

% Analyzed geometry
clear s
s.objs={g5};
s.name={'CO1'};
s.tags={'g5'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Constants
fem.const = {'R','5.3/1000', ...
  'a','0.8/1000'};
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $)

%%% Loop for the wave vector k
n=20   %Number of k values
for ii=1:n
    ii
    kx(ii)=(ii-1)/(n-1)*50
    ky=0
    % Constants
fem.const = {
    %'d','5.3/1000', ...
%  'r','0.8/1000', ...
  'kx' ,num2str(kx(ii)), ......
  'ky',num2str(ky)}
%  'k2','0', ...
%  'a1x','d', ...
%  'a1y','0', ...
%  'a2x','0', ...
%  'a2y','d'};
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $)

% Geometry
g5=move(g5,[5.3/1000/2,5.3/1000/2]);

% Analyzed geometry
clear s
s.objs={g5};
s.name={'CO1'};
s.tags={'g5'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $)

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');
% 
% % Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');
% 
% % Refine mesh
% fem.mesh=meshrefine(fem, ...
%                     'mcase',0, ...
%                     'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'AcoPressure';
appl.module = 'ACO';
appl.assignsuffix = '_acpr';
clear prop
prop.analysis='eigen';
appl.prop = prop;
clear bnd
bnd.type = 'SH';
bnd.symmetryx = {{}};
bnd.symmetryy = {{}};
bnd.fartype = {{}};
bnd.symtypey = {{}};
bnd.symtypex = {{}};
bnd.ind = [1,1,1,1,1,1,1,1];
appl.bnd = bnd;
clear equ
equ.rho = 1000;
equ.cs = 1500;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
%fem.expr = {'kx','k*(k1*b1x+k2*b2x)', ...
 % 'ky','k*(k1*b1y+k2*b2y)', ...
 % 'b1x','2*pi*a2y/(a1x*a2y-a1y*a2x)', ...
 % 'b1y','-2*pi*a2x/(a1x*a2y-a1y*a2x)', ...
 % 'b2x','-2*pi*a1y/(a1x*a2y-a1y*a2x)', ...
 % 'b2y','2*pi*a1x/(a1x*a2y-a1y*a2x)'};

% Descriptions
clear descr
descr.expr= {'kx','1st k-component for periodic condition','ky','2nd k-component for periodic condition','b1y','y-component of 1st reciprocal lattice vector','b2x','x-component of 2nd reciprocal lattice vector','b2y','y-component of 2nd reciprocal lattice vector','b1x','x-component of 1st reciprocal lattice vector'};
fem.descr = descr;

% Coupling variable elements
clear elemcpl
% Extrusion coupling variables
clear elem
elem.elem = 'elcplextr';
elem.g = {'1'};
src = cell(1,1);
clear bnd
bnd.expr = {{'p',{},{}},{{},'p',{}}};
bnd.map = {{'1','1','1'},{'1','1','1'}};
bnd.ind = {{'1'},{'2'},{'3','4','5','6','7','8'}};
src{1} = {{},bnd,{}};
elem.src = src;
geomdim = cell(1,1);
clear bnd
bnd.map = {{{},{},'2'},{{},'3',{}}};
bnd.ind = {{'1','2','5','6','7','8'},{'3'},{'4'}};
geomdim{1} = {{},bnd,{}};
elem.geomdim = geomdim;
elem.var = {'pconstr1','pconstr2'};
map = cell(1,3);
clear submap
submap.type = 'unit';
map{1} = submap;
clear submap
submap.type = 'linear';
submap.sg = '1';
submap.sv = {'7','8'};
submap.dg = '1';
submap.dv = {'1','2'};
map{2} = submap;
clear submap
submap.type = 'linear';
submap.sg = '1';
submap.sv = {'2','8'};
submap.dg = '1';
submap.dv = {'1','7'};
map{3} = submap;
elem.map = map;
elemcpl{1} = elem;
% Point constraint variables (used for periodic conditions)
clear elem
elem.elem = 'elpconstr';
elem.g = {'1'};
clear bnd
bnd.constr = {{'0','pconstr2-(p*exp(i*ky*y))'},{'pconstr1-(p*exp(i*kx*x))', ...
  '0'}};
bnd.cpoints = {{'2','2'},{'2','2'}};
bnd.ind = {{'3'},{'4'}};
elem.geomdim = {{{},bnd,{}}};
elemcpl{2} = elem;
fem.elemcpl = elemcpl;

% Descriptions
descr = fem.descr;
descr.const= {'a2y','y component of 2st lattice vector','d','periodicidad','a1y','y component of 1st lattice vector','a1x','x component of 1st lattice vector','r','radio','k','Fraction of wave vector magnitude','k1','1st component of wave direction vector','a2x','x component of 2nd lattice vector','k2','2nd component of wave direction vector'};
fem.descr = descr;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femeig(fem, ...
               'solcomp',{'p'}, ...
               'outcomp',{'p'}, ...
               'blocksize','auto', ...
               'neigs',1, ...
               'shift',0-2*pi*140000*i, ...
               'eigref','-i*2*pi*140000');
               
% Save current fem structure for restart purposes
fem0=fem;

%%% Save eigenfrequencies
freky(1:1,ii)=real(fem.sol.lambda(1:1)*i/2/pi)

% Plot solution
% postplot(fem, ...
%          'tridata',{'abs(p)','cont','internal','unit','Pa'}, ...
%          'trimap','jet(1024)', ...
%          'solnum',1, ...
%          'title','freq_acpr(1)=-5.695122e-6+21.190758i    Surface: abs(p) [Pa]');%, ...
% %%%         'axis',[-0.0017892932446615213,0.007044900157276357,-1.4675477294442702E-4,0.005824094536405545]);
     
end
%%% Save eigenfrequencies
d=str2num(cell2mat(fem.const(2)));

%%% Plot dispersion relation
plot(max(freky,[],1))