Документация

Необходимый синтаксис

Функция геометрии описывает кривые, которые связывают области геометрии. Кривая является параметризованной функцией (x (t), y (t)). Переменная t находится в областях значений через фиксированный интервал. Для наилучших результатов t должны быть пропорциональны длине дуги плюс константа.

Необходимо задать по крайней мере две кривые для каждой геометрической области. Для примера, 'circleg' функция геометрии, которая доступна в Partial Differential Equation Toolbox™, использует четыре кривые для описания окружности. Кривые могут пересекаться только в начале или конце интервалов параметров.

Функции тулбокса запрашивают геометрическую функцию путем передачи аргументов 0, 1 или 2. Обусловьте свою геометрическую функцию на основе количества входных параметров, чтобы вернуть данные, описанные в этой таблице.

Отношение между параметризацией и метками области

Геометрическая функция для круга

В этом примере показано, как записать геометрическую функцию для создания круговой области. Параметризируйте окружность с радиусом 1, центрированным в источник (0,0), следующим образом:

Геометрическая функция должна иметь как минимум два сегмента. Чтобы удовлетворить этому требованию, разделите круг на четыре сегмента.

Теперь, когда у вас есть параметризация, запишите геометрическую функцию. Сохраните этот файл функции как circlefunction.m на пути MATLAB ®. Эта геометрия проста в создании, поскольку параметризация не меняется в зависимости от номера сегмента.

function [x,y] = circlefunction(bs,s)
% Create a unit circle centered at (0,0) using four segments.
switch nargin
    case 0
        x = 4; % four edge segments
        return
    case 1
        A = [0,pi/2,pi,3*pi/2; % start parameter values
             pi/2,pi,3*pi/2,2*pi; % end parameter values
             1,1,1,1; % region label to left
             0,0,0,0]; % region label to right
        x = A(:,bs); % return requested columns
        return
    case 2
        x = cos(s);
        y = sin(s);
end

Постройте график геометрии с указанием номеров ребер и меток граней.

pdegplot(@circlefunction,'EdgeLabels','on','FaceLabels','on')
axis equal

Вычисление длины дуги для геометрической функции

Этот пример показывает, как создать кардиоидную геометрию с помощью четырех различных методов. Методы являются способами параметризации геометрии с помощью вычислений длины дуги. Кардиоид удовлетворяет уравнению $$r=2(1+\cos \left(\Phi \right))$$.

ezpolar('2*(1+cos(Phi))')

Ниже приведены четыре способа параметризации кардиоида как функции от длины дуги:

function [x,y] = cardioid1(bs,s) 
% CARDIOID1 Geometry file defining the geometry of a cardioid. 

if nargin == 0  
  x = 4; % four segments in boundary
  return 
end

if nargin == 1
  dl = [0    pi/2   pi       3*pi/2
        pi/2   pi     3*pi/2   2*pi
        1      1      1        1
        0      0      0        0];
  x = dl(:,bs);   
  return 
end 

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

nth = 400; % fine polygon, 100 segments per quadrant
th = linspace(0,2*pi,nth); % parametrization
r = 2*(1 + cos(th));
xt = r.*cos(th); % Points for interpolation of arc lengths
yt = r.*sin(th);
% Compute parameters corresponding to the arc length values in s
th = pdearcl(th,[xt;yt],s,0,2*pi); % th contains the parameters
% Now compute x and y for the parameters th
r = 2*(1 + cos(th));
x(:) = r.*cos(th);
y(:) = r.*sin(th);
end

pdegplot('cardioid1','EdgeLabels','on')
axis equal

function [x,y] = cardioid2(bs,s) 
% CARDIOID2 Geometry file defining the geometry of a cardioid. 

if nargin == 0  
  x = 4; % four segments in boundary
  return 
end

if nargin == 1
  dl = [0    pi/2   pi       3*pi/2
        pi/2   pi     3*pi/2   2*pi
        1      1      1        1
        0      0      0        0];
  x = dl(:,bs);   
  return 
end 

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

cbs = find(bs < 3); % upper half of cardioid
fun = @(ss)integral(@(t)sqrt(4*(1 + cos(t)).^2 + 4*sin(t).^2),0,ss);
sscale  = fun(pi);
for ii = cbs(:)' % ensure a row vector
    myfun = @(rr)fun(rr)-s(ii)*sscale/pi;
    theta = fzero(myfun,[0,pi]);
    r = 2*(1 + cos(theta));
    x(ii) = r*cos(theta);
    y(ii) = r*sin(theta);
end
cbs = find(bs >= 3); % lower half of cardioid
s(cbs) = 2*pi - s(cbs);
for ii = cbs(:)'
    theta = fzero(@(rr)fun(rr)-s(ii)*sscale/pi,[0,pi]);
    r = 2*(1 + cos(theta));
    x(ii) = r*cos(theta);
    y(ii) = -r*sin(theta);
end
end

pdegplot('cardioid2','EdgeLabels','on')
axis equal

syms t real
r = 2*(1+cos(t));
x = r*cos(t);
y = r*sin(t);
arcl = simplify(sqrt(diff(x)^2+diff(y)^2));
s = int(arcl,t,0,t,'IgnoreAnalyticConstraints',true)

s = 
$$8\hspace{0.17em}\sin \left(\frac{t}{2}\right)$$

syms s real
th = 2*asin(s/8);
r = 2*(1 + cos(th));
r = expand(r)

r = 
$$4-\frac{s^2}{16}$$

x = r*cos(th);
x = simplify(expand(x))

x = 
$$\frac{s^4}{512}-\frac{3\hspace{0.17em}{s}^2}{16}+4$$

y = r*sin(th);
y = simplify(expand(y))

y = 
$$\frac{s\hspace{0.17em}{\left(64-s^2\right)}^{3/2}}{512}$$

function [x,y] = cardioid3(bs,s) 
% CARDIOID3 Geometry file defining the geometry of a cardioid. 

if nargin == 0  
  x = 4; % four segments in boundary
  return 
end

if nargin == 1
dl = [0   4   8  12
      4   8  12  16
      1   1   1   1
      0   0   0   0];
  x = dl(:,bs);   
  return 
end 

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

cbs = find(bs < 3); % upper half of cardioid
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); % lower half
s(cbs) = 16 - s(cbs); % take the reflection
x(cbs) = s(cbs).^4/512 - 3*s(cbs).^2/16 + 4;
y(cbs) = -s(cbs).*(64 - s(cbs).^2).^(3/2)/512; % negate y
end

pdegplot('cardioid3','EdgeLabels','on')
axis equal

function [x,y] = cardioid4(bs,s) 
% CARDIOID4 Geometry file defining the geometry of a cardioid. 

if nargin == 0  
  x = 4; % four segments in boundary
  return 
end

if nargin == 1
  dl = [0    pi/2   pi       3*pi/2
        pi/2   pi     3*pi/2   2*pi
        1      1      1        1
        0      0      0        0];
  x = dl(:,bs);   
  return 
end 

r = 2*(1 + cos(s)); % s is not proportional to arc length
x = r.*cos(s);
y = r.*sin(s);
end

pdegplot('cardioid4','EdgeLabels','on')
axis equal

subplot(2,2,1)
model = createpde;
geometryFromEdges(model,@cardioid1);
generateMesh(model);
pdeplot(model)
title('Polygons')
axis equal

subplot(2,2,2)
model = createpde;
geometryFromEdges(model,@cardioid2);
generateMesh(model);
pdeplot(model)
title('Integral')
axis equal

subplot(2,2,3)
model = createpde;
geometryFromEdges(model,@cardioid3);
generateMesh(model);
pdeplot(model)
title('Analytic')
axis equal

subplot(2,2,4)
model = createpde;
geometryFromEdges(model,@cardioid4);
generateMesh(model);
pdeplot(model)
title('Distorted')
axis equal

Пример функции геометрии с поддоменами и отверстием

В этом примере показано, как создать файл геометрии для области с поддоменами и отверстием. Он использует раздел «Analytic Arc Length» примера «Arc Length Calculations for a Geometry Function» и вариант функции круга из «Geometry Function for a Circle». Геометрия состоит из внешнего кардиоида, который разделяется на верхнюю половину, называемую субдоменом 1 и нижнюю половину, называемую субдоменом 2. Кроме того, нижняя половина имеет круглое отверстие с центром (1, -1) и радиусом 1/2. Ниже приведен код геометрической функции.

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 cardioid
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

Постройте график геометрии, включая метки ребер и метки поддомена.

pdegplot(@cardg3,'EdgeLabels','on','FaceLabels','on')
axis equal

Вложенная функция для геометрии с дополнительными параметрами

В этом примере показано, как включить дополнительные параметры в функцию для создания 2-D геометрии.

Когда функция 2-D геометрии требует дополнительных параметров, вы не можете использовать стандартный подход анонимной функции, потому что функции геометрии возвращают различное количество аргументов. Вместо этого можно использовать глобальные переменные или вложенные функции. В большинстве случаев рекомендуемый подход заключается в использовании вложенных функций.

Пример решает уравнение Пуассона с нулем граничных условий Дирихле на всех контурах. Геометрия является кардиоидом с эллиптическим отверстием, которое имеет центр в (1, -1) и переменные полуоси. Чтобы настроить и решить модель PDE с этой геометрией, используйте вложенную функцию. Здесь родительская функция принимает длины полуосей, rr и ss, как входные параметры. Причина гнездиться cardioidWithEllipseGeom в пределах cardioidWithEllipseModel это то, что вложенные функции совместно используют рабочую область своих родительских функций. Поэтому cardioidWithEllipseGeom функция может получить доступ к значениям rr и ss что вы передаете в cardioidWithEllipseModel.

function cardioidWithEllipseModel(rr,ss) 

if (rr > 0) & (ss > 0)
  model = createpde(); 
  geometryFromEdges(model,@cardioidWithEllipseGeom); 
  pdegplot(model,'EdgeLabels','on','FaceLabels','on') 
  axis equal 

  applyBoundaryCondition(model,'dirichlet','Edge',1:8,'u',0); 
  specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1);

  generateMesh(model); 
  u = solvepde(model); 
  figure 
  pdeplot(model,'XYData',u.NodalSolution) 
  axis equal

else
    display('Semiaxes values must be positive numbers.')
end 

function [x,y] = cardioidWithEllipseGeom(bs,s) 

if nargin == 0 
 x = 8; % eight segments in boundary
 return 
end 

if nargin == 1 
    % Cardioid 
    dlc = [ 0   4   8  12 
            4   8  12  16 
            1   1   1   1 
            0   0   0   0]; 
    % Ellipse
    dle = [0      pi/2   pi       3*pi/2 
           pi/2   pi     3*pi/2   2*pi 
           0      0      0        0 
           1      1      1        1];
    % Combine the edge matrices 
    dl = [dlc,dle]; 
    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 cardioid 
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 > 4); % Inner ellipse center (1,-1) axes rr and ss 
x(cbs) = 1 + rr*cos(s(cbs)); 
y(cbs) = -1 + ss*sin(s(cbs)); 

end  

end 

При вызове cardioidWithEllipseModel, убедитесь, что значения полуосей достаточно малы, так что эллиптическое отверстие появляется полностью внутри внешнего кардиоида. В противном случае геометрия становится недопустимой.

Для примера вызовите функцию для эллипса с большой полуосью rr = 0.5 и малой полуоси ss = 0.25. Этот вызов функции возвращает следующую геометрию и решение.

cardioidWithEllipseModel(0.5,0.25)