MATLAB Examples

Nested Function for Geometry with Additional Parameters

This example shows how to include additional parameters into a function for creating a 2-D geometry.

When a 2-D geometry function requires additional parameters, you cannot use a standard anonymous function approach because geometry functions return a varying number of arguments. Instead, you can use global variables or nested functions. In most cases, the recommended approach is to use nested functions.

The example solves a Poisson's equation with zero Dirichlet boundary conditions on all boundaries. The geometry is a cardioid with an elliptic hole that has a center at (1,-1) and variable semiaxes. To set up and solve the PDE model with this geometry, use a nested function. Here, the parent function accepts the lengths of the semiaxes, rr and ss, as input parameters. The reason to nest cardioidWithEllipseGeom within cardioidWithEllipseModel is that nested functions share the workspace of their parent functions. Therefore, the cardioidWithEllipseGeom function can access the values of rr and ss that you pass to 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

When calling cardioidWithEllipseModel, ensure that the semiaxes values are small enough, so that the elliptic hole appears entirely within the outer cardioid. Otherwise, the geometry becomes invalid.

For example, call the function for the ellipse with the major semiaxis rr = 0.5 and the minor semiaxis ss = 0.25. This function call returns the following geometry and the solution.

cardioidWithEllipseModel(0.5,0.25)