Параметрическая функция для 2D создания геометрии

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

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

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

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

Количество входных параметров Возвращенные данные

0 (ne = pdegeom) ne количество ребер в геометрии.

Количество входных параметров	Возвращенные данные
0 (`ne = pdegeom`)	`ne` количество ребер в геометрии.
1 (`d = pdegeom(bs)`)	`bs` вектор из сегментов ребра. Ваша функция возвращает `d` когда матрица с одним столбцом для каждого сегмента ребра задана в `bs`. Строки `d` : Запустите значение параметров Закончите значение параметров Оставленная метка области, где “оставлено” относительно направления от запуска в конец значение параметров Правильная метка области Метка области совпадает с номером субдомена. Меткой области внешнего вида геометрии является `0`.
2 (`[x,y] = pdegeom(bs,s)`)	`s` массив длин дуги и `bs` скаляр или массив одного размера с `s` это дает числа ребра. Если `bs` скаляр, затем он применяется к каждому элементу в `s`. Ваша функция возвращает `x` и `y`, которые являются `x` и `y` координаты сегментов ребра заданы в `bs` в значении параметров `s`. `x` и `y` массивы имеют тот же размер как `s`.

1 (d = pdegeom(bs))

bs вектор из сегментов ребра. Ваша функция возвращает d когда матрица с одним столбцом для каждого сегмента ребра задана в bs. Строки d :

Запустите значение параметров
Закончите значение параметров
Оставленная метка области, где “оставлено” относительно направления от запуска в конец значение параметров
Правильная метка области

Метка области совпадает с номером субдомена. Меткой области внешнего вида геометрии является 0.

2 ([x,y] = pdegeom(bs,s)) s массив длин дуги и bs скаляр или массив одного размера с s это дает числа ребра. Если bs скаляр, затем он применяется к каждому элементу в s. Ваша функция возвращает x и y, которые являются x и y координаты сегментов ребра заданы в bs в значении параметров s. x и y массивы имеют тот же размер как s.

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

Следующий рисунок показывает, как направление увеличения параметра имеет отношение к нумерации метки. Стрелки в фигуре показывают направления увеличения значений параметров. Черные точки указывают на начало кривой и конечные точки. Красные числа указывают на метки области. Красный 0 в центре фигуры указывает, что центральный квадрат является отверстием.

Стрелки кривыми 1 и 2 показывают область 1 налево и область 0 направо.
Стрелки кривыми 3 и 4 показывают область 0 налево и область 1 направо.
Стрелки кривыми 5 и 6 показывают область 0 налево и область 1 направо.
Стрелки кривыми 7 и 8 показывают область 1 налево и область 0 направо.

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

Скрипт Open Live Script

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

$x = \cos (t),$

$y = \sin (t),$

$0 \leq t \leq 2 π .$

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

$0 \leq t \leq π / 2$
$π / 2 \leq t \leq π$
$π \leq t \leq 3 π / 2$
$3 π / 2 \leq t \leq 2 π$

Теперь, когда у вас есть параметризация, запишите функцию геометрии. Сохраните этот файл функции как 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

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

Этот пример использует:

Скрипт Open Live Script

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

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

Следующее является этими четырьмя способами параметризовать кардиоиду как функцию длины дуги:

Используйте pdearcl функция многоугольным приближением к геометрии. Этот подход является общим, достаточно точным, и в вычислительном отношении быстро.
Используйте integral и fzero функции, чтобы вычислить длину дуги. Этот подход является более в вычислительном отношении дорогостоящим, но может быть точным, не требуя, чтобы вы выбрали произвольный многоугольник.
Используйте аналитическое вычисление длины дуги. Этот подход является лучшим, когда он применяется, но существует много случаев, где он не применяется.
Используйте параметризацию, которая не пропорциональна длине дуги плюс константа. Этот подход является самым простым, но может дать к искаженной mesh, которая не дает самое точное решение вашей проблемы УЧП.

Многоугольное приближение

Метод конечных элементов использует треугольную mesh, чтобы аппроксимировать решение УЧП численно. Можно избежать потери в точности путем взятия достаточно прекрасного многоугольного приближения к геометрии. pdearcl функционируйте карты между параметризацией и длиной дуги в форме, хорошо подходящей для функции геометрии. Запишите следующую функцию геометрии для кардиоиды.

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

С 400 линейными сегментами выглядит гладко геометрия.

Встроенный cardg функция дает немного отличающуюся версию этого метода.

Интеграл для длины дуги

Можно записать интеграл для длины дуги кривой. Если параметризация в терминах $x (u)$ и $y (u)$ , затем длина дуги $s (t)$

$s (t) = \int_{0}^{t} \sqrt{{(\frac{d x}{d u})}^{2} + {(\frac{d y}{d u})}^{2}} d u .$

Для данного значения $s 0$ , можно найти $t$ как корень уравнения $s (t) = s 0$ . fzero функция решает этот тип нелинейного уравнения.

Запишите следующую функцию геометрии для кардиоидного примера.

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

Геометрия выглядит идентичной многоугольному приближению. Эта интегральная версия берет намного дольше, чтобы вычислить, чем многоугольная версия.

Аналитическая длина дуги

Также можно найти аналитическое выражение для длины дуги как функция параметризации. Затем можно дать параметризацию в терминах длины дуги. Например, найдите аналитическое выражение для длины дуги при помощи Symbolic Math Toolbox™.

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 \sin (\frac{t}{2})$

В терминах длины дуги s, параметр t t = 2*asin(s/8), где s диапазоны от 0 до 8, соответствуя t в пределах от 0 к $π$ . Для s между 8 и 16, симметрией кардиоиды, t = pi + 2*asin((16-s)/8). Кроме того, можно описать x и y в терминах s этими аналитическими вычислениями.

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 s^{2}}{16} + 4$

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

y = 
 $\frac{s {(64 - s^{2})}^{3 / 2}}{512}$

Теперь, когда у вас есть аналитические выражения для x и y в терминах длины дуги s, запишите функцию геометрии.

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

Эта аналитическая геометрия выглядит немного более сглаженной, чем предыдущие версии. Однако различие несущественно в терминах вычислений.

Геометрия, не пропорциональная длине дуги

Также можно записать функцию геометрии, где параметр не пропорционален длине дуги. Этот подход может дать к искаженной mesh.

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

Метки равномерно не расположены с интервалами на ребрах, потому что параметр не пропорционален длине дуги.

Исследуйте mesh по умолчанию на каждый из четырех методов создания геометрии.

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

Искаженная mesh выглядит немного менее регулярной, чем другие сетки. Это имеет некоторые очень узкие треугольники около острого выступа кардиоиды. Тем не менее, все сетки, кажется, применимы.

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

Скрипт Open Live Script

В этом примере показано, как создать файл геометрии для области с субдоменами и отверстия. Это использует раздел "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
          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

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

Скрипт Open Live Script

В этом примере показано, как включать дополнительные параметры в функцию для создания 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)

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