1

我正在使用 Delaunay 三角剖分将多边形拆分为三角形。我使用大代码处理 FEM，我的“检查点”之一是对称性（如果数据是对称的，则输出也必须是对称的）。但是，由于我无法控制 Delaunay 三角剖分，它使我失去了对称性。

我写了一个小代码来说明我的问题：我们考虑两个不相交的三角形和一个与它们相交的大矩形。我们想用矩形对这些三角形的减法进行三角测量：

clear all
close all
warning off % the warning is about duplicate points, not important here

figure
hold on

p =[.3 .3
.4 .3
.3 .4
.7 .6
.6 .7
.7 .7]; % coordinates of the points for the triangles

px = 1/3;
py = 1/3;
lx = 1/3;
ly = 1/3; % size and position of the rectangle

% rearrange the polygon with clockwise-ordered vertices
[x1,y1]=poly2cw([px; px+lx; px+lx; px],[py; py; py+ly; py+ly]); % rectangle

patch(x1,y1, 1, 'EdgeColor', 'k');

for i=1:2

    pc = p(3*i-2:3*i,:); % current triangle
    % rearrange the polygon with clockwise-ordered vertices
    [x0,y0]=poly2cw(pc(:,1),pc(:,2)); % triangle

    [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection
    [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction

    DT = delaunayTriangulation(x3,y3);

    triplot(DT,'Marker','o')

end
XL = xlim; xlim(XL+[-1 +1]*diff(XL)/10);
YL = ylim; ylim(YL+[-1 +1]*diff(YL)/10);
axis equal;
box on;

delaunay三角剖分

如您所见，Delaunay 三角剖分在两个三角形中的行为不同，因此失去了对称性。

有没有一种简单的方法来恢复对称性？

我使用 Matlab R2013a。

4

2 回答 2

0

看来这不是一个错误。实际上，您的结果来自您的数据。

我玩了一点你的代码

clear all
close all
warning off % the warning is about duplicate points, not important here

figure
hold on

p =[.3 .3
.4 .3
.3 .4
.7 .6
.6 .7
.7 .7]; % coordinates of the points for the triangles

px = 1/3;
py = 1/3;
lx = 1/3;
ly = 1/3; % size and position of the rectangle

% rearrange the polygon with clockwise-ordered vertices
[x1,y1]=poly2cw([px; px+lx; px+lx; px],[py; py; py+ly; py+ly]); % rectangle

patch(x1,y1, 1, 'EdgeColor', 'k');

for i=1:2

    pc = p(3*i-2:3*i,:); % current triangle
    % rearrange the polygon with clockwise-ordered vertices
    [x0,y0]=poly2cw(pc(:,1),pc(:,2)); % triangle

    [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection
    [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction

    % This is for R2013a
    %DT = delaunayTriangulation(x3,y3);
    %triplot(DT,'Marker','o');

    % This is for R2011b
    %DT = DelaunayTri(x3,y3);
    %triplot(DT,'Marker','o');

    % This is plain delaunay version
    DT = delaunay(x3,y3);
    triplot(DT,x3,y3,'Marker','o')
    
    % we break here to analyze the first triangulation
    break

end
XL = xlim; xlim(XL+[-1 +1]*diff(XL)/10);
YL = ylim; ylim(YL+[-1 +1]*diff(YL)/10);
axis equal;
box on;

% % % % % % % % % % % % % % % % % %
% Checking the triangulation
% % % % % % % % % % % % % % % % % %

% Wrong triangulation for i=2 is hard-coded
DT2     = [
    2 1 6
    6 5 2
    5 3 2
    5 4 3
    2 3 1 ];

figure;
hold all;
triplot(DT2,x3,y3,'Marker','o', 'Color','r', 'LineWidth',1)
axis equal;
axis tight;
box on;
XL = xlim; xlim(XL+[-1 +1]*diff(XL)*0.5);
YL = ylim; ylim(YL+[-1 +1]*diff(YL)*0.5);

% circumcircle: http://www.mathworks.com/matlabcentral/fileexchange/17300
ca = linspace(0,2*pi);
cx = cos(ca);
cy = sin(ca);

hl = [];
for k=1:size(DT2,1)
    tx  = x3(DT(k,:));
    ty  = y3(DT(k,:));
    [r,cn]=circumcircle([tx,ty]',0);
    if ~isempty(hl)
        %delete(hl);
    end
    fprintf('Circle %d: Center at (%.23f,%.23f); R=%.23f\n',k,cn,r);
    text( cn(1),cn(2), sprintf('c%d',k) );
    hl = plot( cx*r+cn(1), r*cy+cn(2), 'm-' );
    drawnow;
    %pause(3); %if you prefere to go slowly 
end

这是我看到的输出：

圆圈1：中心在（0.28333333333333333000000,0.35000000000000003000000）；R=0.05270462766947306400000 圆2：中心在（0.34999999999999998000000,0.34999999999999998000000）；R=0.02357022603955168100000 圆3：中心在（0.28333333333333338000000,0.34999999999999992000000）；R=0.05270462766947289800000 圆 4：中心在 (0.35000000000000003000000,0.28333333333333355000000); R=0.05270462766947290500000 圆 5：中心在 (0.35000000000000003000000,0.28333333333333333000000); R=0.05270462766947312000000

和图：

在此处输入图像描述

所以圈 1 和 3 以及圈 4 和 5 几乎相同。因此，您和我的结果之间的差异甚至可能来自舍入误差，因为相应的四个点在浮点数学精度内位于同一个圆上。你必须重新设计你的观点以获得不依赖于这些东西的可靠结果。

玩得开心;o)

于 2013-04-26T14:57:49.133 回答
0

好像您正在使用 MatLab R2013，因为在我的 R2011b 中没有delaunayTriangulation功能。为了能够运行您的代码，我稍微更改了它：

clear all
close all
warning off % the warning is about duplicate points, not important here

figure
hold on

p =[.3 .3
    .4 .3
    .3 .4
    .7 .6
    .6 .7
    .7 .7]; % coordinates of the points for the triangles

px = 1/3;
py = 1/3;
lx = 1/3;
ly = 1/3; % size and position of the rectangle

% rearrange the polygon with clockwise-ordered vertices
[x1,y1]=poly2cw([px; px+lx; px+lx; px],[py; py; py+ly; py+ly]); % rectangle

patch(x1,y1, 1, 'EdgeColor', 'k');

for i=1:2

    pc = p(3*i-2:3*i,:); % current triangle
    % rearrange the polygon with clockwise-ordered vertices
    [x0,y0]=poly2cw(pc(:,1),pc(:,2)); % triangle

    [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection
    [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction

    %DT = delaunayTriangulation(x3,y3);
    %triplot(DT)

    % This is triangulation of subtraction
    DT = delaunay(x3,y3);
    triplot(DT,x3,y3, 'Marker','.', 'Color','r')

    % This is triangulation of intersection
    DT = delaunay(x2,y2);
    triplot(DT,x2,y2, 'Marker','o', 'Color','b', 'LineWidth',1)

end
axis equal;
axis tight;
box on;

XL = xlim; xlim(XL+[-1 +1]*diff(XL)/10);
YL = ylim; ylim(YL+[-1 +1]*diff(YL)/10);

text(0.5,0.55,'triangulation of subtraction',  'HorizontalAlignment','center', 'VerticalAlignment','bottom', 'Color','r');
text(0.5,0.45,'triangulation of intersection', 'HorizontalAlignment','center', 'VerticalAlignment','top',    'Color','b');

这是我看到的结果

在此处输入图像描述

和你的一样吗？你能否在你的问题中添加一张图片来描述你得到的结果有什么问题？

于 2013-04-25T06:25:19.830 回答