-
Notifications
You must be signed in to change notification settings - Fork 2
/
meshUniformRefineTriangular.asv
135 lines (120 loc) · 5.42 KB
/
meshUniformRefineTriangular.asv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
function [meshData] = meshUniformRefineTriangular(meshData)
% Refine a triangular mesh uniformly.
%
% Remark: After unifirm refine, edges may not start from points with a
% lower index. However, this program guarantees that the edges and boundary
% edges have the same direction. The reason that we do not want to make
% all edge direction from smaller point indices to larger is that, this
% breaks the coarse-fine mesh relationship which is important in multigrid
% methods.
if ~strcmp(meshData.elementType,'triangular')
error('meshUniformRefineTriangular only works for triangular meshes!');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Refine
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 0: Set some quantities
% ---------------------------
np = meshData.nP;
ne = meshData.nE;
nt = meshData.nT;
nbe = meshData.nBE;
ne2 = 2*ne;
nt3 = nt*3;
% Step 1: Add new vertices
% ---------------------------
meshData.P(:, (np+1):(np+ne)) = ...
( meshData.P(:,meshData.E(1,:)) + meshData.P(:,meshData.E(2,:)) ) *0.5;
% Step 2: Add new edges
% ---------------------------
meshData.E(:, (ne+1):(ne2+nt3)) = 0;
% set two sub-edges of all original edges.
% The edge order is:
% first, ne edges with vertices [startP, centerP]
% second, ne edges with vertices [endP, centerP]
meshData.E(:, (ne+1):(ne2)) = [meshData.E(2,1:ne); meshData.nP+(1:ne)];
meshData.E(2, 1:ne) = np+(1:ne);
% Set three new edges inside each original triangle.
% Pay attention to the order and orientation.
meshData.E(:, (ne2+1):3:end) = np + abs([meshData.T2E(1,:); meshData.T2E(3,:)]);
meshData.E(:, (ne2+2):3:end) = np + abs([meshData.T2E(2,:); meshData.T2E(1,:)]);
meshData.E(:, (ne2+3):3:end) = np + abs([meshData.T2E(3,:); meshData.T2E(2,:)]);
% Step 3: Add new triangles
% ---------------------------
meshData.T(:, (nt+1):(4*nt)) = 0;
% Add three sub-triangles who have one original vertex.
% Pay attention to the orientation of these triangles.
meshData.T(:, (nt+1):3:end) = ...
[meshData.T(1,1:nt); np+abs(meshData.T2E(1,:)); ...
np+abs(meshData.T2E(3,:))];
meshData.T(:, (nt+2):3:end) = ...
[meshData.T(2,1:nt); np+abs(meshData.T2E(2,:)); ...
np+abs(meshData.T2E(1,:))];
meshData.T(:, (nt+3):3:end) = ...
[meshData.T(3,1:nt); np+abs(meshData.T2E(3,:)); ...
np+abs(meshData.T2E(2,:))];
% The new sub-triangle at the center inherits
% the original triangle's index.
meshData.T(:, 1:nt) = np+abs(meshData.T2E);
% Set subdomains
meshData.TSubdomN(1,(nt+1):(4*nt)) = 0;
meshData.TSubdomN((nt+1):3:end) = meshData.TSubdomN(1:nt);
meshData.TSubdomN((nt+2):3:end) = meshData.TSubdomN(1:nt);
meshData.TSubdomN((nt+3):3:end) = meshData.TSubdomN(1:nt);
% Step 4: Add new boundary edges
% ---------------------------
meshData.BE(:, (nbe+1):(2*nbe)) = ...
[meshData.BE(2,1:nbe); np+meshData.BE2E];
meshData.BE(2, 1:nbe) = np+meshData.BE2E(:);
meshData.BESegN((nbe+1):(2*nbe)) = meshData.BESegN(1:nbe);
meshData.BEOrient((nbe+1):(2*nbe)) = -meshData.BEOrient(1:nbe);
% Step 5: Set T2E
% ---------------------------
% Each column of temp will store all six sub-edges
% on the boundary of a triangle, in counterclockwise direction,
% starting from the one [vertex1, venterOfOriginalEdge1].
% Each sub-edge has the orientation that points to the center
% of the original edges.
temp = zeros(6, nt);
temp([1,3,5],:) = abs(meshData.T2E);
temp([2,4,6],:) = meshData.nE+abs(meshData.T2E);
ind = ( meshData.T2E(1,:)<0 );
temp([1,2], ind) = temp([2,1], ind);
ind = ( meshData.T2E(2,:)<0 );
temp([3,4], ind) = temp([4,3], ind);
ind = ( meshData.T2E(3,:)<0 );
temp([5,6], ind) = temp([6,5], ind);
meshData.T2E(:, (nt+1):(4*nt)) = 0;
meshData.T2E(:, (nt+1):3:end) = [temp(1,:); ne2+(1:3:nt3); -temp(6,:)];
meshData.T2E(:, (nt+2):3:end) = [temp(3,:); ne2+(2:3:nt3); -temp(2,:)];
meshData.T2E(:, (nt+3):3:end) = [temp(5,:); ne2+(3:3:nt3); -temp(4,:)];
meshData.T2E(:, 1:nt) = -(ne2 + [2:3:nt3; 3:3:nt3; 1:3:nt3]);
% Step 6: Set E2T and ELInd
% ---------------------------
meshData.E2T(:, 1:(ne2+nt3)) = 0;
meshData.ELInd(:, 1:(ne2+nt3)) = 0;
% For sub-edges of original edges.
% We will use the fact that all subedges are pointed in to the center
% of the original edge.
meshData.E2T(1,meshData.T2E(1,(nt+1):end)) = (nt+1):(4*nt);
meshData.E2T(2,-meshData.T2E(3,(nt+1):end)) = (nt+1):(4*nt);
meshData.ELInd(1,meshData.T2E(1,(nt+1):end)) = 1;
meshData.ELInd(2,-meshData.T2E(3,(nt+1):end)) = 3;
% For three new internal edges.
meshData.E2T(:, (ne2+1):end) = [(nt+1):(4*nt); reshape([1:nt;1:nt;1:nt],1,nt3)];
meshData.ELInd(1,(ne2+1):end) = 2;
meshData.ELInd(2,(ne2+1):3:end) = 3;
meshData.ELInd(2,(ne2+2):3:end) = 1;
meshData.ELInd(2,(ne2+3):3:end) = 2;
% Step 7: Set E2BE and BE2E
% ---------------------------
indE = find(meshData.E2BE>0);
meshData.E2BE((ne+1):(ne2+nt3)) = 0;
meshData.E2BE(ne+indE) = nbe + meshData.E2BE(indE);
meshData.BE2E((nbe+1):(2*nbe)) = ne + meshData.BE2E(1:nbe);
% Step 8: Set numbers
% ---------------------------
meshData.nP = meshData.nP+meshData.nE;
meshData.nE = 2*meshData.nE+3*meshData.nT;
meshData.nT = 4*meshData.nT;
meshData.nBE = 2*meshData.nBE;