用参数函数创建几何模型

有的研究对象或区域的形状可以用参数函数进行描述，比如心形线、椭圆形等。此时将参数函数按照工具箱的要求写成MATLAB函数，创建PDE模型时可以直接用函数名指定几何模型。

实现参数函数的MATLAB函数的结构如下面代码所示。该函数分无输入参数、有一个输入参数和有两个输入参数等3种情况进行描述。

不指定输入参数时只需指定边线段的条数，调用该MATLAB函数的函数会自动对边线分段并编号，得到下面函数的第1个输入参数。然后调用有一个输入参数的该函数进行处理。

有一个输入参数时用dl矩阵描述每条边线的起始和终止角度，以及该边线与相邻面的拓扑关系。调用该MATLAB函数的函数会自动处理得到下面函数的第2个输入参数。然后调用有两个输入参数的该函数进行处理。

有两个输入参数时，可以直接计算出每条边线细分成很多小段后各小段端点的坐标，并作为函数的返回值返回。

function [x,y] = myregionfunction(bs,s)
% 输入参数bs指定边线段的编号，s为弧长矩阵
% 输出参数x和y分别保存边线段分成很多小段后各端点的横坐标和纵坐标
% 按照下面格式自定义参数函数
switch nargin
    case 0  % 如果没有输入参数
        x = 4; % 指定边线段的条数
        return
    case 1  % 如果有一个输入参数
        dl = [];  % 描述每条边线段拓扑关系的矩阵
        x = dl(:,bs);  % 返回需要处理的列
        return
    case 2  % 如果有两个输入参数
        % 将每条边线段细分成很多小段，计算各顶点的坐标
        x = …;
        y = …;
end

上面代码中，dl矩阵具有类似下面的构造形式。这里描述的参数曲线被分成4段，dl矩阵的每列数据描述其中的某1段。每列有4个数据，前两个数据指定当前分段的起点和终点对应的值(如角度)，后两个数字描述当前分段与相邻面的拓扑关系。在当前分段的前进方向上，第3个数据为其左侧面的编号，第4个数据为其右侧面的编号。

dl = [0 pi/2 pi 3*pi/2

pi/2 pi 3*pi/2 2*pi

1 1 1 1

0 0 0 0];

如图4-5中所示，将矩形区域的边分成4段，编号分别为1, 2, 3和4，其中第1和2段的前进方向为逆时针方向，第3和4段的前进方向为顺时针方向。矩形区域内面的编号为1，区域外面的编号为0。所以，对于第1和2段边线，在其前进方向上，左侧面的编号为1，右侧面的编号为0；对于第3和4段边线，在其前进方向上，左侧面的编号为0，右侧面的编号为1。

图4-5 边线段与面之间的拓扑关系

【例10】下面用椭圆方程定义PDE模型的几何模型，用myellipse函数根据椭圆方程获取椭圆上若干顶点的坐标。椭圆是曲线，在椭圆上相隔很小的距离连续取点，然后将它们连起来，可以得到一个逼近曲线的封闭多边形。

function [x,y] = myellipse(bs,s)
% 参数函数定义的椭圆区域.
if nargin == 0
  x = 4;   % 4条边线段
  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
% 如果输入参数为2个，将每条边线段细分成若干段
% 获取各断点的横坐标和纵坐标
x = 3.*cos(s);
y = 2.*sin(s);
end

下面用pdegplot函数调用上面定义的myellipse函数，绘制椭圆方程定义的几何模型。

>> pdegplot('myellipse','EdgeLabels','on','FaceLabels','on')

生成图4-6。

图4-6 用参数函数生成的几何模型

下面创建一个PDE模型，使用geometryFromEdges函数指定PDE模型的几何模型为myellipse函数定义的椭圆形区域。

>> model=createpde;
>> gp=geometryFromEdges(model,@myellipse)

使用pdegplot函数绘制几何模型的图形。

>> pdegplot(model,'EdgeLabels','on','FaceLabels','on')

生成图4-6。