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