Функция геометрии описывает кривые, которые связывают области геометрии. Кривая является параметризованной функцией (x (t), y (t)). Переменная t находится в областях значений через фиксированный интервал. Для наилучших результатов t должны быть пропорциональны длине дуги плюс константа.
Необходимо задать по крайней мере две кривые для каждой геометрической области. Для примера, 'circleg' функция геометрии, которая доступна в Partial Differential Equation Toolbox™, использует четыре кривые для описания окружности. Кривые могут пересекаться только в начале или конце интервалов параметров.
Функции тулбокса запрашивают геометрическую функцию путем передачи аргументов 0, 1 или 2. Обусловьте свою геометрическую функцию на основе количества входных параметров, чтобы вернуть данные, описанные в этой таблице.
| Количество Входных параметров | Возвращенные данные |
|---|---|
0 (ne = pdegeom) | ne - количество ребер в геометрии. |
1 (d = pdegeom(bs)) |
Метка области совпадает с числом поддомена. Метка области внешней геометрии |
2 ([x,y] = pdegeom(bs,s)) | s является массивом длин дуг и bs является скаляром или массивом того же размера, что и s который задает ребро числа. Если bs является скаляром, затем применяется к каждому элементу в s. Ваша функция возвращается x и y, которые являются x и y координаты сегментов ребра, указанные в bs при значении параметров s. The 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, чтобы аппроксимировать решение к УЧП численно. Можно избежать потерь точности, приняв достаточно тонкое полигональное приближение к геометрии. The 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 функция дает несколько другую версию этого метода.
Интеграл для длины дуги
Можно записать интеграл для длины дуги кривой. Если параметризация в терминах и , затем длина дуги является
Для заданного значения , вы можете найти как корень уравнения . The 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» примера «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)
