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

    %% Geometry Function Example with Subdomains and a Hole
% This example shows how to create a geometry file for a region with
% subdomains and a hole. It uses the “Analytic” cardioid
% example from "Arc Length Calculations for a Geometry Function" and a
% variant of the circle function from "Geometry Function for a Circle".
%
% The geometry consists of an outer cardioid that is divided into an upper
% half called subdomain 1 and a lower half called subdomain 2. Also, the
% lower half has a circular hole centered at (1,–1) and of radius
% 1/2. Here is the code of the geometry function.
%
%  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
%%
% Plot the geometry, including edge labels and subdomain labels.
pdegplot(@cardg3,'EdgeLabels','on','SubdomainLabels','on')
axis equal