Функция геометрии описывает кривые, которые связали области геометрии. Кривая является параметрической функцией (x (t), y (t)). Переменная t передвигается за фиксированный интервал. Для лучших результатов t должен быть пропорционален длине дуги плюс константа.
Необходимо задать по крайней мере две кривые для каждой геометрической области. Например, 'circleg' функция геометрии, которая доступна в Partial Differential Equation Toolbox™, использует четыре кривые, чтобы описать круг. Кривые могут пересечься только вначале или конец интервалов параметра.
Функции тулбокса запрашивают вашу функцию геометрии путем передачи в 0, 1, или 2 аргумента. Conditionalize ваша геометрия функционируют на основе количества входных параметров, чтобы возвратить данные, описанные в этой таблице.
| Количество входных параметров | Возвращенные данные |
|---|---|
0ne = pdegeom) | ne количество ребер в геометрии. |
1d = pdegeom(bs)) |
Метка области совпадает с номером субдомена. Меткой области внешнего вида геометрии является |
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 направо.
В этом примере показано, как записать, что геометрия функционирует для создания круговой области. Параметризуйте круг с радиусом 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
В этом примере показано, как создать кардиоидную геометрию с помощью четырех отличных методов. Методы являются способами параметризовать вашу геометрию с помощью вычислений длины дуги. Кардиоида удовлетворяет уравнению .
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 функция дает немного отличающуюся версию этого метода.
Интеграл для длины дуги
Можно записать интеграл для длины дуги кривой. Если параметризация в терминах и , затем длина дуги
Для данного значения , можно найти как корень уравнения . 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 =
В терминах длины дуги 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 =
x = r*cos(th); x = simplify(expand(x))
x =
y = r*sin(th); y = simplify(expand(y))
y =
Теперь, когда у вас есть аналитические выражения для 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 выглядит немного менее регулярной, чем другие сетки. Это имеет некоторые очень узкие треугольники около острого выступа кардиоиды. Тем не менее, все сетки, кажется, применимы.
В этом примере показано, как создать файл геометрии для области с субдоменами и отверстия. Это использует раздел "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
В этом примере показано, как включать дополнительные параметры в функцию для создания 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)
