www.gusucode.com > pde 案例源码 matlab代码程序 > pde/cardg3.m

    function [x,y] = cardg3(bs,s) 
% CARDG3 Geometry File defining the geometry of a cardioid with two
% subregions and a hole.

if nargin == 0  
  x = 9; % 9 segments
  return 
end

if nargin == 1
    % Outer cardioid
    dl = [0   4   8  12
          4   8  12  16
          1   1   2   2 % Region 1 to the left in the upper half, 2 in the lower
          0   0   0   0];
    % Dividing line between top and bottom
    dl2 = [0
        4
        1 % Region 1 to the left
        2]; % Region 2 to the right
    % Inner circular hole
    dl3 = [0      pi/2   pi       3*pi/2
           pi/2   pi     3*pi/2   2*pi
           0      0      0        0 % To the left is empty
           2      2      2        2]; % To the right is region 2
    % Combine the three edge matrices
    dl = [dl,dl2,dl3];

    x = dl(:,bs);   
    return 
end 

x = zeros(size(s)); 
y = zeros(size(s)); 
if numel(bs) == 1 % Does bs need scalar expansion?
    bs = bs*ones(size(s)); % Expand bs
end


cbs = find(bs < 3); % Upper half of cardiod
x(cbs) = s(cbs).^4/512 - 3*s(cbs).^2/16 + 4;
y(cbs) = s(cbs).*(64 - s(cbs).^2).^(3/2)/512;
cbs = find(bs >= 3 & bs <= 4); % Lower half of cardioid
s(cbs) = 16 - s(cbs);
x(cbs) = s(cbs).^4/512 - 3*s(cbs).^2/16 + 4;
y(cbs) = -s(cbs).*(64 - s(cbs).^2).^(3/2)/512;
cbs = find(bs == 5); % Index of straight line
x(cbs) = s(cbs);
y(cbs) = zeros(size(cbs));
cbs = find(bs > 5); % Inner circle radius 0.25 center (1,-1)
x(cbs) = 1 + 0.25*cos(s(cbs));
y(cbs) = -1 + 0.25*sin(s(cbs));
end