Для основной информации о создании 2D геометрии см. Три Способа Создать 2D Геометрию.
Функция геометрии описывает кривые, которые связали области геометрии. Кривая является параметрической функцией (x (t), y (t)). Переменная t передвигается за фиксированный интервал. Для лучших результатов t должен быть пропорционален длине дуги плюс константа.
Необходимо задать по крайней мере две кривые для каждой геометрической области. Например, функция геометрии 'circleg'
, которая доступна в Partial Differential Equation Toolbox™, использует четыре кривые, чтобы описать круг. Кривые могут пересечься только вначале или конец интервалов параметра.
Функции тулбокса запрашивают вашу функцию геометрии путем передачи в 0, 1, или 2 аргумента. Conditionalize ваша геометрия функционируют на основе количества входных параметров, чтобы возвратить данные, описанные в этой таблице.
Количество входных параметров | Возвращенные данные |
---|---|
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)
, можно следующим образом:
Функция геометрии должна иметь по крайней мере два сегмента. Чтобы удовлетворить это требование, разбейте круг в четыре сегмента.
Теперь, когда у вас есть параметризация, запишите функцию геометрии. Сохраните этот файл функции как 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
с многоугольным приближением к геометрии. Этот подход является общим, достаточно точным, и в вычислительном отношении быстро.
Используйте интеграл и функции 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 = 8*sin(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 - s^2/16
x = r*cos(th); x = simplify(expand(x))
x = s^4/512 - (3*s^2)/16 + 4
y = r*sin(th); y = simplify(expand(y))
y = (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 выглядит немного менее регулярной, чем другие сетки. Это имеет некоторые очень узкие треугольники около острого выступа кардиоиды. Тем не менее, все сетки, кажется, применимы.
Этот пример показывает, как создать файл геометрии для области с субдоменами и дыры. Это использует раздел "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)