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