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.

57 lines
2.4 KiB

1 year ago
function [Nquads,value] = recurse_quad(CPWmat,bdbx,funct,p,initlvls,lvls,tol)
%recurse_quad Performs recursive quadtree integration
% Inputs:
% CPWmat One bezier object representing one boundary of domain
% bdbx bdbx a quad which intersects the boundary of the domain
% funct Function to be integrated f = @(x,y)...
% p Order of tensor-product Gaussian quad rule to be used
% lvls Levels of adaptive quadtree subdivision to perform
%
% Outputs:
% Nquads Totol number of quadrature points used
% value Estimate of value of integral over domain using quadtree
Nquads = 0;
value = 0;
if initlvls>0
daughtbxs{1} = [bdbx(1) (bdbx(1)+bdbx(2))/2 bdbx(3) (bdbx(3) + bdbx(4))/2];
daughtbxs{2} = [(bdbx(1)+bdbx(2))/2 bdbx(2) bdbx(3) (bdbx(3) + bdbx(4))/2];
daughtbxs{3} = [(bdbx(1)+bdbx(2))/2 bdbx(2) (bdbx(3) + bdbx(4))/2 bdbx(4)];
daughtbxs{4} = [bdbx(1) (bdbx(1)+bdbx(2))/2 (bdbx(3) + bdbx(4))/2 bdbx(4)];
for i=1:4
[nqt,vt] = recurse_quad(CPWmat,daughtbxs{i},funct,p,initlvls-1,lvls,tol);
value = value + vt; Nquads = Nquads + nqt;
end
else
if lvls==0
quadpts = gauss2Dbdspts(bdbx,p);
quadpts(:,1:2)=shapIEvalpts(quadpts(:,1:2),CPWmat,
shapeIEvalpts([x y],funct,CPWmat,1,tol);
quadpts = gauss2Dbds(Ifun, bdbx, p);
Nquads = p^2;
else
daughtbxs{1} = [bdbx(1) (bdbx(1)+bdbx(2))/2 bdbx(3) (bdbx(3) + bdbx(4))/2];
daughtbxs{2} = [(bdbx(1)+bdbx(2))/2 bdbx(2) bdbx(3) (bdbx(3) + bdbx(4))/2];
daughtbxs{3} = [(bdbx(1)+bdbx(2))/2 bdbx(2) (bdbx(3) + bdbx(4))/2 bdbx(4)];
daughtbxs{4} = [bdbx(1) (bdbx(1)+bdbx(2))/2 (bdbx(3) + bdbx(4))/2 bdbx(4)];
for i=1:4
if lvls==3
bdpts = [daughtbxs{i}(1) daughtbxs{i}(3);
daughtbxs{i}(1) daughtbxs{i}(4);
daughtbxs{i}(2) daughtbxs{i}(3);
daughtbxs{i}(2) daughtbxs{i}(4);];
% scatter(bdpts(:,1),bdpts(:,2),100,'r.');
end
isContID=isContainedBx(daughtbxs{i},CPWmat,1,tol,p);
if isContID == 0
value = value + gauss2Dbds(funct,daughtbxs{i},p);
Nquads = Nquads + p^2;
elseif isContID == 1
[nqt,vt] = recurse_quad(CPWmat,daughtbxs{i},funct,p,0,lvls-1,tol);
value = value + vt; Nquads = Nquads + nqt;
end
end
end
end
end