Геометрическая функция описывает кривые, которые связывают геометрические области. Кривая является параметризованной функцией (x (t), y (t)). Переменная t находится в диапазоне фиксированного интервала. Для достижения наилучших результатов t должно быть пропорционально длине дуги плюс константа.
Необходимо указать не менее двух кривых для каждой геометрической области. Например, 'circleg' геометрическая функция, доступная в дифференциальном уравнении 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. x и y массивы имеют тот же размер, что и s. |
На следующем рисунке показано, как направление увеличения параметра соотносится с нумерацией меток. Стрелки на рисунке показывают направления увеличения значений параметров. Черные точки указывают начальную и конечную точки кривой. Красные цифры обозначают метки областей. Красный 0 в центре фигуры указывает на то, что центральный квадрат является отверстием.
Стрелки по кривым 1 и 2 показывают область 1 слева и область 0 справа.
Стрелки по кривым 3 и 4 показывают область 0 слева и область 1 справа.
Стрелки по кривым 5 и 6 показывают область 0 слева и область 1 справа.
Стрелки по кривым 7 и 8 показывают область 1 слева и область 0 справа.
В этом примере показано, как записать функцию геометрии для создания круговой области. Параметризовать окружность с радиусом 1 с центром в начале координат (0,0), следующим образом:
t),
t),
Геометрическая функция должна иметь не менее двух сегментов. Чтобы удовлетворить этому требованию, разбейте круг на четыре сегмента.
Теперь, когда у вас есть параметризация, напишите функцию геометрии. Сохранить этот файл функции как 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
В этом примере показано, как создать кардиоидную геометрию с помощью четырех различных методов. Методы представляют собой способы параметризации геометрии с помощью расчетов длины дуги. Кардиоид удовлетворяет уравнению cos (Start)).
ezpolar('2*(1+cos(Phi))')
Следующие четыре способа параметризовать кардиоид как функцию длины дуги:
Используйте pdearcl с многоугольным приближением к геометрии. Этот подход является общим, достаточно точным и вычислительно быстрым.
Используйте integral и fzero для вычисления длины дуги. Этот подход является более затратным с точки зрения вычислений, но может быть точным, не требуя выбора произвольного многоугольника.
Используйте аналитический расчет длины дуги. Этот подход является наилучшим, когда он применяется, но есть много случаев, когда он не применяется.
Используйте параметризацию, которая не пропорциональна длине дуги плюс константа. Этот подход является самым простым, но может дать искаженную сетку, которая не дает наиболее точного решения вашей проблемы PDE.
Многоугольное приближение
Метод конечных элементов использует треугольную сетку для аппроксимации решения к PDE численно. Можно избежать потери точности, приняв достаточно точное многоугольное приближение к геометрии. 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 функция даёт несколько другую версию этого приёма.
Интеграл для длины дуги
Можно записать интеграл для длины дуги кривой. Если параметризация в терминах ) и u), то дуги
dydu) 2du.
Для заданного значения можно найти в качестве корня уравнения s0. 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
Геометрия выглядит идентично многоугольной аппроксимации. Для вычисления этой интегральной версии требуется гораздо больше времени, чем для многоугольной версии.
Аналитическая длина дуги
Можно также найти аналитическое выражение для длины дуги как функции параметризации. Затем можно задать параметризацию в терминах длины дуги. Например, найдите аналитическое выражение для длины дуги с помощью символьного математического 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
Эта аналитическая геометрия выглядит несколько более гладкой, чем предыдущие версии. Однако разница несущественна с точки зрения расчетов.
Геометрия, не пропорциональная длине дуги
Можно также записать геометрическую функцию, где параметр не пропорционален длине дуги. Такой подход может привести к искаженной сетке.
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
Искаженная сетка выглядит немного менее регулярной, чем другие сетки. У него очень узкие треугольники рядом с бугорком кардиоида. Тем не менее, все сетки кажутся пригодными для использования.
В этом примере показано, как создать файл геометрии для области с поддоменами и отверстием. Он использует раздел «Аналитическая длина дуги» примера «Вычисления длины дуги для функции геометрии» и вариант функции окружности из «Геометрическая функция для окружности». Геометрия состоит из внешнего кардиоида, который делится на верхнюю половину, называемую поддоменом 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)
