You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
432 lines
17 KiB
432 lines
17 KiB
function [avgDAT,DAT] = convergenceAnalysis(filename,orientations,functs,poldeg,timings)
|
|
% Performs convergence analysis of the integral of funct over the
|
|
% interior of the region defined by the set of oriented NURBS cuves
|
|
% given by CP, W, and Xi using six strategies:
|
|
% 1. Quadtree integration
|
|
% Uses my own adaptive quadtree implementation
|
|
% 2. Linear meshing and order 3 Gaussian quadrature under h-refinement
|
|
% using a modified version of mesh2d (xmesh provided in TRIGA)
|
|
% 3. Exact rational meshing and order 3 Gaussian quadrature under
|
|
% p-refinement using xmesh2d from TRIGA
|
|
% 4. Polynomial Spline approximation of the boundary followed by
|
|
% Green's theorem-based quadrature using splinegauss
|
|
% 5. SPECTRAL quadrature using our method
|
|
% 6. SPECTRALPE quadrature using our method
|
|
%
|
|
% Input:
|
|
% filename A string containing the filename where the boundary
|
|
% NURBS curves are stored, with the extension _NURBS
|
|
% removed. (e.g. two_rotors will load the file
|
|
% two_rotors_NURBS). These NURBS curves bound the
|
|
% integraiton region.
|
|
% orientations A vectors of 1's and -1's representing the
|
|
% orientations of each of the boundary loops in the
|
|
% region in filename
|
|
% functs A list of functions for which integrals will be
|
|
% computed.
|
|
% poldegs If the functions are polynomial, this contains an
|
|
% integer for each functs which describes its degree.
|
|
% If the functions are nonpolynomial, this contains a
|
|
% non-integer representing what degree of polynomial
|
|
% exactness should be guaranteed for the SPECTRALPE
|
|
% method.
|
|
% timings 0 if timings are not desired. 1 if timings are
|
|
% Note that timings==1 takes much longer.
|
|
% Output:
|
|
% avgDAT Contains four lists of data for each
|
|
% method:
|
|
% 1. Number of quad points
|
|
% 2. Average quadrature error for all functions in
|
|
% functs for the associated # of quad points
|
|
% 3. Pre-processing time for the associated # of quad
|
|
% points (if timings==1)
|
|
% 4. Average evaluation time for all functions in
|
|
% functs for the associate # of quadrature points
|
|
% (if timings ==1)
|
|
% DAT Formatted the same as avgDAT, except contains one cell
|
|
% for each function (i.e., is not an average over the
|
|
% functs).
|
|
|
|
% Turn off irritating warnings
|
|
warning('off','all');
|
|
|
|
% This hard-codes how many iterations of each strategy should be
|
|
% employed.
|
|
maxmeshQuad=10;
|
|
if strcmp(filename,'treble_clef') || strcmp(filename,'guitar')
|
|
maxquadmeshlin=5;
|
|
maxGGquadssq=35;
|
|
numQuadquad=3;
|
|
elseif strcmp(filename,'two_rotors')
|
|
maxquadmeshlin=5;
|
|
maxGGquadssq=55;
|
|
numQuadquad=4;
|
|
else
|
|
maxquadmeshlin=6;
|
|
maxGGquadssq=25;
|
|
numQuadquad=5;
|
|
end
|
|
maxGGquads=ceil(sqrt(maxGGquadssq));
|
|
if mod(poldeg,1)==0
|
|
maxSGquads=1;
|
|
maxVianquad=10;
|
|
maxVianquadLin=10;
|
|
else
|
|
maxSGquads=maxGGquads;
|
|
maxVianquad=9;
|
|
maxVianquadLin=9;
|
|
end
|
|
|
|
%% On off switch
|
|
% If only some of the comparison methods are desired, the others can be
|
|
% turned off using this onoff list (0= off, 1= on). The order is the
|
|
% same as in the description of the function.
|
|
onoff=[0 0 0 0 1 1 0];
|
|
|
|
for jj=1:1
|
|
if ~onoff(1)
|
|
maxmeshQuad=0;
|
|
end
|
|
if ~onoff(2)
|
|
maxquadmeshlin=0;
|
|
end
|
|
if ~onoff(3)
|
|
numQuadquad=0;
|
|
end
|
|
if ~onoff(4)
|
|
maxVianquad=0;
|
|
end
|
|
if ~onoff(5)
|
|
maxGGquads=0;
|
|
end
|
|
if ~onoff(6)
|
|
maxSGquads=0;
|
|
end
|
|
if ~onoff(7)
|
|
maxVianquadLin=0;
|
|
end
|
|
end
|
|
|
|
%% Initialize Output
|
|
avgDAT=cell(1,7);
|
|
avgDAT{1}=zeros(4,maxmeshQuad);
|
|
avgDAT{2}=zeros(4,maxquadmeshlin);
|
|
avgDAT{3}=zeros(4,numQuadquad);
|
|
avgDAT{4}=zeros(4,maxVianquad);
|
|
avgDAT{5}=zeros(4,maxGGquads);
|
|
avgDAT{6}=zeros(4,maxSGquads);
|
|
avgDAT{7}=zeros(4,maxVianquadLin);
|
|
for funi=1:length(functs)
|
|
DAT{funi}=avgDAT;
|
|
end
|
|
|
|
% Load the NURBS curves.
|
|
load([filename,'_NURBS'],'CP','W','Xi');
|
|
|
|
% Convert the NURBS format to my format
|
|
CPW_my_format=extract_my_format_from_NURBS(CP,W,Xi);
|
|
|
|
% Save the NURBS format as a spline file for TRIGA
|
|
toSplineFile(filename,CP,W,Xi);
|
|
[xtrue,ytrue,wtrue] = SPECTRAL_quads(CPW_my_format,55,55);
|
|
|
|
|
|
for funi=1:length(functs)
|
|
% Creates new function which can track number of evaluations
|
|
funct=functs{funi};
|
|
|
|
% Integrate using szego-gauss
|
|
fprintf('Performing our Fejer-Gauss method \n');
|
|
p=length(Xi{1})-length(W{1})-1;
|
|
if mod(poldeg,1)==0
|
|
timingstrue=zeros(1,1);
|
|
timingsAtrue=zeros(1,1);
|
|
f = @() SPECTRALPE_quads(CPW_my_format,0,p*(poldeg+3),ceil((poldeg+1)/2));
|
|
[xq,yq,wq]=f();
|
|
nQuadsSzego=length(xq);
|
|
g= @() applyRule(xq,yq,wq,funct);
|
|
truev= g();
|
|
if timings
|
|
timingsAtrue=timeit(g);
|
|
timingstrue=timeit(f);
|
|
end
|
|
else
|
|
truev=zeros(1,maxSGquads);
|
|
nQuadsSzego=zeros(1,maxSGquads);
|
|
timingstrue=zeros(1,maxSGquads);
|
|
timingsAtrue=zeros(1,maxSGquads);
|
|
for ii=1:maxSGquads
|
|
f = @() SPECTRALPE_quads(CPW_my_format,ii^2,p*(floor(poldeg)+3),ii^2);
|
|
[xq,yq,wq]=f();
|
|
g= @() applyRule(xq,yq,wq,funct);
|
|
truev(ii)= g();
|
|
if timings
|
|
timingstrue(ii)=timeit(f);
|
|
timingsAtrue(ii)=timeit(g);
|
|
end
|
|
nQuadsSzego(ii)=length(xq);
|
|
end
|
|
end
|
|
|
|
|
|
% Integrate using green-gauss
|
|
nQuadsGG=zeros(1,maxGGquads); valueGG=zeros(1,maxGGquads); timingsGG=zeros(1,maxGGquads);
|
|
timingsAGG=zeros(1,maxGGquads);
|
|
for ii=1:maxGGquads
|
|
fprintf('Performing our Gauss-Green method with %d orders of accuracy \n',ii^2);
|
|
tic
|
|
if mod(poldeg,1)==0
|
|
f = @() SPECTRAL_quads(CPW_my_format,ii^2,ceil((poldeg+1)/2));
|
|
% valueGG(ii)=-SPECTRALquads(CPW_my_format,field,ii^2,ceil((poldeg+1)/2));
|
|
else
|
|
f = @() SPECTRAL_quads(CPW_my_format,ii^2,ii^2);
|
|
% valueGG(ii)=-SPECTRALquads(CPW_my_format,field,ii^2,ii^2);
|
|
end
|
|
[xq,yq,wq]=f();
|
|
nQuadsGG(ii)=length(xq);
|
|
g= @() applyRule(xq,yq,wq,funct);
|
|
valueGG(ii)= g();
|
|
if timings
|
|
timingsAGG(ii)=timeit(g);
|
|
timingsGG(ii)=timeit(f);
|
|
end
|
|
end
|
|
|
|
fprintf("Computing true value using SPECTRAL quadrature with 55 quad points \n");
|
|
trueI(funi) = applyRule(xtrue,ytrue,wtrue,funct);
|
|
DAT{funi}([5,6])={...
|
|
[nQuadsGG; valueGG; timingsGG; timingsAGG],...
|
|
[nQuadsSzego; truev; timingstrue; timingsAtrue]};
|
|
end
|
|
|
|
|
|
% Integrate using quadtree decomposition using 2nd order
|
|
% Gaussian quadrature
|
|
warning('off','all');
|
|
nQuadsquad = zeros(length(functs),numQuadquad); valuequad = zeros(length(functs),numQuadquad);
|
|
timingsQuadquad = zeros(1,numQuadquad); timingsAQuadquad=zeros(length(functs),numQuadquad);
|
|
% initlvls=0;
|
|
for intlvls=1:numQuadquad
|
|
tol=1e-8;
|
|
fprintf('Performing quadtree quadrature with %d maximum levels \n',intlvls);
|
|
f=@() quadtreeIntegratePts(CPW_my_format,intlvls,3,tol,orientations);
|
|
[xq,yq,wq]=f();
|
|
nQuadsquad(:,intlvls)=length(xq);
|
|
if timings
|
|
timingsQuadquad(intlvls)=timeit(f);
|
|
end
|
|
for funi=1:length(functs)
|
|
funct=functs{funi};
|
|
g=@() applyRule(xq,yq,wq,funct);
|
|
valuequad(funi,intlvls)=g();
|
|
if timings
|
|
timingsAQuadquad(funi,intlvls)=timeit(g);
|
|
end
|
|
end
|
|
end
|
|
for funi=1:length(functs)
|
|
DAT{funi}(3)= {[nQuadsquad(funi,:); valuequad(funi,:); timingsQuadquad; timingsAQuadquad(funi,:);]};
|
|
end
|
|
|
|
% Integrate using Sommariva/Vianello strategy (approximate boundary by
|
|
% polynomial spline, then effectively use gauss-green)
|
|
% integrateVian(CPWmat,numSplinepts,
|
|
nQuadsVian=zeros(length(functs),maxVianquad); valueVian=zeros(length(functs),maxVianquad);
|
|
timingsVian=zeros(1,maxVianquad); timingsAVian=zeros(length(functs),maxVianquad);
|
|
for ll=3:maxVianquad
|
|
kk=2^(ll-1);
|
|
fprintf('Performing cubic spline approximation quadrature with %d sample points per curve \n',kk);
|
|
if mod(poldeg,1)==0
|
|
N_2009=10;
|
|
else
|
|
N_2009=22;
|
|
end
|
|
f=@() performVianQuadPts(CPW_my_format,kk,N_2009);
|
|
[xq,yq,wq]=f();
|
|
for funi=1:length(functs)
|
|
funct=functs{funi};
|
|
g=@()applyRule(xq,yq,wq,funct);
|
|
valueVian(funi,ll)= g();
|
|
if timings
|
|
timingsAVian(funi,ll)=timeit(g);
|
|
end
|
|
end
|
|
nQuadsVian(:,ll)=length(xq);
|
|
if timings
|
|
timingsVian(ll)=timeit(f);
|
|
end
|
|
end
|
|
for funi=1:length(functs)
|
|
DAT{funi}(4)= {[nQuadsVian(funi,:); valueVian(funi,:); timingsVian; timingsAVian(funi,:);]};
|
|
end
|
|
|
|
% Integrate using Sommariva/Vianello strategy (approximate boundary by
|
|
% polynomial spline, then effectively use gauss-green)
|
|
% integrateVian(CPWmat,numSplinepts,
|
|
nQuadsVianLin=zeros(length(functs),maxVianquadLin); valueVianLin=zeros(length(functs),maxVianquadLin);
|
|
timingsVianLin=zeros(1,maxVianquadLin); timingsAVianLin=zeros(length(functs),maxVianquadLin);
|
|
for ll=3:maxVianquadLin
|
|
kk=2^(ll-1);
|
|
fprintf('Performing linear approximation quadrature with %d sample points per curve \n',kk);
|
|
if mod(poldeg,1)==0
|
|
N_2009=5;
|
|
else
|
|
N_2009=15;
|
|
end
|
|
f=@() performVianQuadPtsLin(CPW_my_format,kk,N_2009);
|
|
[xq,yq,wq]=f();
|
|
for funi=1:length(functs)
|
|
funct=functs{funi};
|
|
g=@()applyRule(xq,yq,wq,funct);
|
|
valueVianLin(funi,ll)= g();
|
|
if timings
|
|
timingsAVianLin(funi,ll)=timeit(g);
|
|
end
|
|
end
|
|
nQuadsVianLin(:,ll)=length(xq);
|
|
if timings
|
|
timingsVianLin(ll)=timeit(f);
|
|
end
|
|
end
|
|
for funi=1:length(functs)
|
|
DAT{funi}(7)= {[nQuadsVianLin(funi,:); valueVianLin(funi,:); timingsVianLin; timingsAVianLin(funi,:);]};
|
|
end
|
|
% Exact and linear mesh the geometry, saved in filename.neu and
|
|
% filenamelin.neu
|
|
% Integrate over the exact and linear meshes using pth order Gaussian
|
|
% quadrature
|
|
nQuadsmesh= zeros(length(functs),maxmeshQuad); valuemesh= zeros(length(functs),maxmeshQuad);
|
|
timingsQuadmesh=zeros(1,maxmeshQuad); timingsAQuadmesh=zeros(length(functs),maxmeshQuad);
|
|
options.thresh=2;
|
|
if maxmeshQuad>1
|
|
fprintf('Attempting to exactly mesh the geometry \n');
|
|
tic
|
|
xmesh2d(filename,options);
|
|
meshtime=toc;
|
|
for i=1:maxmeshQuad
|
|
tic
|
|
fprintf('Performing exact rational meshing quadrature with %d Gauss points per triangle \n',i^2);
|
|
% p=length(Xi{1})-length(W{1});
|
|
for funi=1:length(functs)
|
|
funct=functs{funi};
|
|
f=@() integrate2dPts(filename,i);
|
|
[xq,yq,wq]=f();
|
|
nQuadsmesh(funi,i)=length(xq);
|
|
g=@() applyRule(xq,yq,wq,funct);
|
|
valuemesh(funi,i)=g();
|
|
if timings
|
|
timingsAQuadmesh(funi,i)=timeit(g);
|
|
end
|
|
end
|
|
% refineMeshRefLvl(filename);
|
|
if timings
|
|
timingsQuadmesh(i)=timeit(f)+meshtime;
|
|
end
|
|
end
|
|
end
|
|
for funi=1:length(functs)
|
|
DAT{funi}(1)={[nQuadsmesh(funi,:); valuemesh(funi,:);timingsQuadmesh; timingsAQuadmesh(funi,:);]};
|
|
end
|
|
|
|
|
|
|
|
nQuadsmeshlin= zeros(length(functs),maxquadmeshlin); valuemeshlin= zeros(length(functs),maxquadmeshlin);
|
|
timingsQuadmeshlin=zeros(1,maxquadmeshlin);
|
|
timingsAQuadmeshlin=zeros(length(functs),maxquadmeshlin);
|
|
for i=1:maxquadmeshlin
|
|
fprintf('Performing linear meshing quadrature with %d refinements \n',i-1);
|
|
options.thresh=1 + 1*(10^(1-i));
|
|
options.hmax=.3*(.5)^(1-i);
|
|
tic
|
|
[pts,tri]=xmesh2dlin(filename,options);
|
|
meshtime=toc;
|
|
for funi=1:length(functs)
|
|
funct=functs{funi};
|
|
f=@() integrate2dlinPts(pts,tri,p);
|
|
[xq,yq,wq] = f();
|
|
nQuadsmeshlin(funi,i)= length(xq);
|
|
g=@() applyRule(xq,yq,wq,funct);
|
|
valuemeshlin(funi,i)=g();
|
|
if timings
|
|
timingsAQuadmeshlin(funi,i)=timeit(g);
|
|
end
|
|
end
|
|
if timings
|
|
timingsQuadmeshlin(i)=timeit(f)+meshtime;
|
|
end
|
|
end
|
|
for funi=1:length(functs)
|
|
DAT{funi}(2)={[nQuadsmeshlin(funi,:); valuemeshlin(funi,:); timingsQuadmeshlin; timingsAQuadmeshlin]};
|
|
end
|
|
|
|
|
|
for lpl=1:7
|
|
if onoff(lpl)
|
|
for funi=1:length(functs)
|
|
avgDAT{lpl}(2:4,:)=avgDAT{lpl}(2:4,:)+(abs(DAT{funi}{lpl}(2:4,:))-[abs(trueI(funi));0; 0])./(length(functs)*[abs(trueI(funi));1; 1]);
|
|
end
|
|
avgDAT{lpl}(1,:)=DAT{funi}{lpl}(1,:);
|
|
end
|
|
end
|
|
|
|
end
|
|
|
|
function [xq,yq,wq]=performVianQuadPts(CPW_my_format,kk,N_2009)
|
|
xq=zeros(0,1); yq=zeros(0,1); wq=zeros(0,1);
|
|
rotation=0; P=[0 0]; Q=[0 0]; cumulative=1;
|
|
cubature_type=4; SPLtypestring='complete';
|
|
numEvalspcurve=kk;evalPts=linspace(0,1-(1/numEvalspcurve),numEvalspcurve);
|
|
% spline_order_vett{1}=[2*ones(4,1) (1:numEvalspcurve:(numEvalspcurve*4))'];
|
|
for ii=1:length(CPW_my_format)
|
|
% fst_spline_order=(bkpts{ii}-[1 bkpts{ii}(1:end-1)]);
|
|
% fst_spline_order(fst_spline_order>0)=3; fst_spline_order=fst_spline_order+1;
|
|
% spline_order_vett{ii}=[fst_spline_order' (kk*bkpts{ii}-kk+1)'];
|
|
spline_order_vett{ii}=[repmat(4,size(CPW_my_format{ii},1)/3,1) ...
|
|
(kk:kk:(kk*size(CPW_my_format{ii},1)/3))'+1];
|
|
clear vertices;
|
|
for jj=1:3:(size(CPW_my_format{ii},1))
|
|
verts=dCR_eval(CPW_my_format{ii}(jj:(jj+2),:),evalPts);
|
|
vertices((((jj-1)/3)*numEvalspcurve+1):((jj-1)/3+1)*numEvalspcurve,:)=verts;
|
|
end
|
|
% spline_order_vett{ii}=[4 length(vertices)+1];
|
|
vertices(end+1,:)=vertices(1,:);
|
|
f=@() splinegauss_2009b(N_2009,...
|
|
vertices,rotation,P,Q,spline_order_vett{ii},...
|
|
cumulative,SPLtypestring,cubature_type);
|
|
[x_nodes_2009,y_nodes_2009,weights_2009]=f();
|
|
xq=[xq; x_nodes_2009]; yq=[yq; y_nodes_2009];
|
|
wq=[wq; weights_2009];
|
|
end
|
|
end
|
|
|
|
function [xq,yq,wq]=performVianQuadPtsLin(CPW_my_format,kk,N_2009)
|
|
xq=zeros(0,1); yq=zeros(0,1); wq=zeros(0,1);
|
|
rotation=0; P=[0 0]; Q=[0 0]; cumulative=1;
|
|
cubature_type=4; SPLtypestring='complete';
|
|
numEvalspcurve=kk;evalPts=linspace(0,1-(1/numEvalspcurve),numEvalspcurve);
|
|
% spline_order_vett{1}=[2*ones(4,1) (1:numEvalspcurve:(numEvalspcurve*4))'];
|
|
for ii=1:length(CPW_my_format)
|
|
% fst_spline_order=(bkpts{ii}-[1 bkpts{ii}(1:end-1)]);
|
|
% fst_spline_order(fst_spline_order>0)=3; fst_spline_order=fst_spline_order+1;
|
|
% spline_order_vett{ii}=[fst_spline_order' (kk*bkpts{ii}-kk+1)'];
|
|
spline_order_vett{ii}=[repmat(2,size(CPW_my_format{ii},1)/3,1) ...
|
|
(kk:kk:(kk*size(CPW_my_format{ii},1)/3))'+1];
|
|
clear vertices;
|
|
for jj=1:3:(size(CPW_my_format{ii},1))
|
|
verts=dCR_eval(CPW_my_format{ii}(jj:(jj+2),:),evalPts);
|
|
vertices((((jj-1)/3)*numEvalspcurve+1):((jj-1)/3+1)*numEvalspcurve,:)=verts;
|
|
end
|
|
% spline_order_vett{ii}=[4 length(vertices)+1];
|
|
vertices(end+1,:)=vertices(1,:);
|
|
f=@() splinegauss_2009b(N_2009,...
|
|
vertices,rotation,P,Q,spline_order_vett{ii},...
|
|
cumulative,SPLtypestring,cubature_type);
|
|
[x_nodes_2009,y_nodes_2009,weights_2009]=f();
|
|
xq=[xq; x_nodes_2009]; yq=[yq; y_nodes_2009];
|
|
wq=[wq; weights_2009];
|
|
end
|
|
end
|
|
|
|
function I = applyRule(xq,yq,wq,funct)
|
|
I=wq'*funct(xq,yq);
|
|
end
|