Note: CREATING 2-D GEOMETRIES IS THE SAME FOR THE RECOMMENDED AND THE LEGACY WORKFLOWS. |
For basic information on 2-D geometry construction, see Create 2-D Geometry
A geometry function describes the curves that bound the geometry regions. A curve is a parametrized function (x(t),y(t)). The variable t ranges over a fixed interval. For best results, t should be proportional to arc length plus a constant.
For each region you should have at least two curves. For example,
the 'circleg'
geometry function, which ships with
the toolbox, uses four curves to describe a circle.
Curves can intersect only at the beginning or end of parameter intervals.
Toolbox functions query your geometry function by passing in 0, 1, or 2 arguments. Conditionalize your geometry function based on the number of input arguments to return the following:
Number of Input Arguments | Returned Data |
---|---|
0 (ne = pdegeom ) | ne is the number of edges in the geometry. |
1 (d = pdegeom(bs) ) | bs is a vector of edge segments. Your function
returns d as a matrix with one column for each
edge segment specified in bs . The rows of d are:
0 . |
2 ([x,y] = pdegeom(bs,s) ) | s is an array of arc lengths, and bs is
a scalar or an array the same size as s giving
edge numbers. If bs is a scalar, then it applies
to every element in s . Your function returns x and y ,
which are the x and y coordinates
of the edge segments specified in bs at parameter
value s . The x and y arrays
have the same size as s . |
This figure shows how the direction of parameter increase relates to label numbering. The arrows in the following figure show the directions of increasing parameter values. The black dots indicate curve beginning and end points. The red numbers indicate region labels. The red 0 in the center of the figure indicates that the center square is a hole.
The arrows by curves 1
and 2
show
region 1
to the left and region 0
to
the right.
The arrows by curves 3
and 4
show
region 0
to the left and region 1
to
the right.
The arrows by curves 5
and 6
show
region 0
to the left and region 1
to
the right.
The arrows by curves 7
and 8
show
region 1
to the left and region 0
to
the right.
xs = [0,1,1,0,0]; xt = [0.25,0.75,0.75,0.25,0.25]; ys = [0,0,1,1,0]; yt = [0.25,0.25,0.75,0.75,0.25]; plot(xs,ys,'b-') hold on axis equal ylim([-0.1,1.1]) plot(xt,yt,'b-') plot(xs,ys,'k.','MarkerSize',12) plot(xt,yt,'k.','MarkerSize',12) text(-0.1,0.5,'0','Color','r') text(0.5,0.5,'0','Color','r') text(1.1,0.5,'0','Color','r') text(0.15,0.15,'1','Color','r') text(0.85,0.85,'1','Color','r') text(0.5,0,'1') text(0.99,0.5,'2') text(0.5,1,'3') text(-0.01,0.5,'4') text(0.5,0.25,'5') text(0.74,0.5,'6') text(0.5,0.75,'7') text(0.24,0.5,'8') annotation('arrow',[0.6,0.7],[0.18,0.18]) annotation('arrow',[0.6,0.7],[0.86,0.86]) annotation('arrow',[0.26,0.26],[0.7,0.8]) annotation('arrow',[0.77,0.77],[0.7,0.8]) annotation('arrow',[0.53,0.63],[0.35,0.35]) annotation('arrow',[0.53,0.63],[0.69,0.69]) annotation('arrow',[0.39,0.39],[0.55,0.65]) annotation('arrow',[0.645,0.645],[0.55,0.65])
This example shows how to use a geometry function to create a circular region. Of course, you could just as easily use a circle basic shape.
You can parametrize a circle with radius 1 centered at the origin (0,0) as follows:
A geometry function needs to have at least two segments. So break up the circle into four segments:
Now that you have a parametrization, write the geometry function. Save this function file as circlefunction.m
on your MATLAB path.
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
This geometry was particularly simple to create because the parameterization did not change depending on the segment number.
Visualize the geometry, edge numbers, and domain label.
pdegplot(@circlefunction,'EdgeLabels','on','SubdomainLabels','on') axis equal
The built-in function circleg
gives a slightly different parameterization of the circle. You might find it instructive to compare the two approaches.
This example shows how to create a cardioid geometry using four distinct techniques. The techniques are ways to parametrize your geometry using arc length calculations. The cardioid satisfies the equation .
ezpolar('2*(1+cos(Phi))')
This example shows four approaches to giving a parametrization of the cardioid as a function of arc length.
Use the pdearcl
function with a polygonal approximation to the geometry. This approach is general, accurate enough, and computationally fast.
Use the integral and fzero functions to compute the arc length. This approach is more computationally costly, but can be accurate without you having to choose an arbitrary polygon.
Use a parametrization that is not proportional to arc length plus a constant. This approach is simplest, but can yield a distorted mesh that does not give the most accurate solution to your PDE problem.
Polygonal Approximation
The finite element method uses a triangular mesh to approximate the solution to a PDE numerically. So there is no loss in accuracy by taking a sufficiently fine polygonal approximation to the geometry. The pdearcl
function maps between parametrization and arc length in a form well-suited to a geometry function. Here is a geometry function for the cardioid.
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
Plot the geometry function.
pdegplot('cardioid1','EdgeLabels','on') axis equal
With 400 line segments, the geometry looks smooth.
The built-in cardg
function gives a slightly different version of this technique.
Integral for Arc Length
You can write an integral for the arc length of a curve. If the parametrization is in terms of and , then the arclength is
So for a given value
, you can find
as the root of the equation
. The fzero
function solves this type of nonlinear equation.
For the present example of a cardioid, here is the calculation.
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 cardiod 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
Plot the geometry function.
pdegplot('cardioid2','EdgeLabels','on') axis equal
The geometry looks identical to the polygonal approximation. This integral version takes much longer to calculate than the polygonal version.
Analytic Arc Length
If you are handy with integrals, or have Symbolic Math Toolbox, you can find an analytic expression for the arc length as a function of the parametrization. Then you can give the parametrization in terms of arc length. Here is an approach using 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)
So you see that, in terms of arclength s
, the parameter t
is t = 2*asin(s/8)
where s
ranges from 0 to 8, corresponding to t
ranging from 0 to
. For s
between 8 and 16, by symmetry if the cardioid, t = pi + 2*asin((16-s)/8)
. Furthermore, you can express x and y in terms of s by the following analytic calculations.
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
Now that you have analytic expressions for x and y in terms of the arclength s, you can write the geometry function. 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 2 2 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 cardiod 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
Plot the geometry function.
pdegplot('cardioid3','EdgeLabels','on') axis equal
This analytic geometry looks slightly smoother than the previous versions. However, the difference is inconsequential in terms of calculations.
Geometry Not Proportional to Arc Length
You can write a geometry function where the parameter is not proportional to arc length. This can yield a distorted 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
Plot the geometry function.
pdegplot('cardioid4','EdgeLabels','on') axis equal
The labels are not evenly spaced on the edges because the parameter is not proportional to arc length.
Examine the default mesh for each of the four methods of creating geometry.
subplot(2,2,1) [p,e,t] = initmesh(@cardioid1); pdeplot(p,e,t) title('Polygons') axis equal subplot(2,2,2) [p,e,t] = initmesh(@cardioid2); pdeplot(p,e,t) title('Integral') axis equal subplot(2,2,3) [p,e,t] = initmesh(@cardioid3); pdeplot(p,e,t) title('Analytic') axis equal subplot(2,2,4) [p,e,t] = initmesh(@cardioid4); pdeplot(p,e,t) title('Distorted') axis equal
While the "Distorted" mesh looks a bit less regular than the other meshes (it has some very narrow triangles near the cusp of the cardioid), all of the meshes appear to be usable.
This example shows how to create a geometry file for a region with subdomains and a hole. It uses the "Analytic" cardioid example from "Arc Length Calculations for a Geometry Function" and a variant of the circle function from "Geometry Function for a Circle".
The geometry consists of an outer cardioid that is divided into an upper half called subdomain 1 and a lower half called subdomain 2. Also, the lower half has a circular hole centered at (1,–1) and of radius 1/2. Here is the code of the geometry function.
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 cardiod 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
Plot the geometry, including edge labels and subdomain labels.
pdegplot(@cardg3,'EdgeLabels','on','SubdomainLabels','on') axis equal