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

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

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

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

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

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

Функция геометрии для круга

В этом примере показано, как записать, что геометрия функционирует для создания круговой области. Параметризуйте круг с радиусом 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" "Вычислений длины дуги для примера" Функции Геометрии и варианта круговой функции от "Функции геометрии для Круга". Геометрия состоит из внешней кардиоиды, которая разделена на верхнюю половину названного субдомена 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
       % Region 1 to the left in
       % the upper half, 2 in the lower
       1   1   2   2 
       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 % Empty to the left
        2      2      2        2]; % Region 2 to the right
    % 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

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

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

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

Пример решает уравнение Пуассона с нулем граничные условия Дирихле на всех контурах. Геометрия является кардиоидой с эллиптическим отверстием, которое имеет центр в (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)