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.
71 lines
2.8 KiB
71 lines
2.8 KiB
function [xquad,yquad,wquad] = SPECTRALPE_quads(Cellintel,pts,pp,kg)
|
|
% Calculates quadrature points for the region defined by Cellintel using
|
|
% the "SPECTRALPE" method from the paper "spectral mesh-free quadrature for
|
|
% planar regions bounded by rational parametric curves."
|
|
% Input:
|
|
% Cellintel Oriented boundaries of the region
|
|
% pts Number of intermediate quadrature points to use
|
|
% pp Number of poles to use in the intermediate quadrature rule
|
|
% (calculated as (k+3)*m, where k is the degree of
|
|
% polynomial exactness desired and m is the degree of the
|
|
% component curves
|
|
% kg Number of antiderivative quadrature points to use
|
|
% (calculated as ceil((k+1)/2), where k is the degree of
|
|
% polynomial exactness desired.
|
|
%
|
|
% Output:
|
|
% [xquad, yquad] Quadrature points locations
|
|
% wquad Quadrature weights
|
|
xquad=zeros(0,1); yquad=zeros(0,1); wquad=zeros(0,1);
|
|
for i=1:length(Cellintel)
|
|
cp=Cellintel{i}(1:3:end,:);
|
|
for j=1:3:(size(Cellintel{i},1))
|
|
rts=roots(fliplr(BernsteinToMonomial(Cellintel{i}(j+2,:))));
|
|
sgl=repmat(rts,ceil(pp/(size(Cellintel{1},2)-1)),1);
|
|
sgl((length(sgl)+1):(pts),:)=[Inf];
|
|
sglout = (2*sgl-1);
|
|
if isempty(rts)
|
|
[xg, wg ]=GaussLegendreQuad1D(pp+pts,0,1);
|
|
x=[xg';wg];
|
|
else
|
|
x = rfejer(sglout);
|
|
if any(x(2,:)<0)
|
|
[xg, wg ]=GaussLegendreQuad1D(pp+pts,0,1);
|
|
x=[xg';wg];
|
|
elseif all(abs(Cellintel{i}(j+2,:)-1)<1e-15)
|
|
[xg, wg ]=GaussLegendreQuad1D(pp+pts,0,1);
|
|
x=[xg';wg];
|
|
else
|
|
x(2,:)=x(2,:)*.5;
|
|
x(1,:)=(x(1,:)+1)/2;
|
|
end
|
|
end
|
|
[xc,yc,yp]= RatGreensFunctionQuad(x(1,:),Cellintel{i}(j:(j+2),:));
|
|
clear xq; clear wq; clear yq;
|
|
for ii=1:length(xc)
|
|
[xxq,wwq]=GaussLegendreQuad1D(kg,min(min(cp(:))),xc(ii));
|
|
xq(ii,:)=xxq;
|
|
wq(ii,:)=wwq.*yp(ii).*x(2,ii);
|
|
end
|
|
yq=repmat(yc,1,kg);
|
|
xquad=[xquad; xq(:)];
|
|
yquad=[yquad; yq(:)];
|
|
wquad=[wquad; wq(:)];
|
|
end
|
|
end
|
|
end
|
|
|
|
function [x,y,yp] = RatGreensFunctionQuad(t,Side)
|
|
% Finds x(t),y(t), and y'(t) on a Bezier Curve defined in Side
|
|
% Input:
|
|
% t Parameter values at which evaluation should occur
|
|
% Side Rational Bezier Curve
|
|
%
|
|
% Output:
|
|
% [x,y,yp] x,y, and dy/ds values for "Side" at t-values
|
|
for i=1:length(t)
|
|
xy=dCR_eval(Side,t(i));
|
|
x(i,:)=xy(:,1)';y(i,:)=xy(:,2)';
|
|
yp(i,:) = dCR_eval_dt(Side(2:3,:),t(i));
|
|
end
|
|
end
|