-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathquadpts.m
150 lines (148 loc) · 7.75 KB
/
quadpts.m
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
function [lambda,weight] = quadpts(order)
%% QUADPTS quadrature points in 2-D.
%
% [lambda,weight] = quadpts(order) return quadrature points with given
% order (up to 9) in the barycentric coordinates.
%
% The output lambda is a matrix of size nQ by 3, where nQ is the number of
% quadrature points. lambda(i,:) is three barycentric coordinate of the
% i-th quadrature point and lambda(:,j) is the j-th barycentric coordinate
% of all quadrature points. The x-y coordinate of the p-th quadrature point
% can be computed as
%
% pxy = lambda(p,1)*node(elem(:,1),:) ...
% + lambda(p,2)*node(elem(:,2),:) ...
% + lambda(p,3)*node(elem(:,3),:);
%
% The weight of p-th quadrature point is given by weight(p). See
% verifyquadpts for the usage of qudrature rules to compute integrals over
% triangles.
%
% References:
%
% * David Dunavant. High degree efficient symmetrical Gaussian
% quadrature rules for the triangle. International journal for numerical
% methods in engineering. 21(6):1129--1148, 1985.
% * John Burkardt. DUNAVANT Quadrature Rules for the Triangle.
% http://people.sc.fsu.edu/~burkardt/m_src/dunavant/dunavant.html
%
% See also quadpts1, quadpts3, verifyquadpts
%
% Order 6 - 9 is added by Huayi Wei, modify by Jie Zhou
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.
if order>9
order = 9;
end
switch order
case 1 % Order 1, nQuad 1
lambda = [1/3, 1/3, 1/3];
weight = 1;
case 2 % Order 2, nQuad 3
lambda = [2/3, 1/6, 1/6; ...
1/6, 2/3, 1/6; ...
1/6, 1/6, 2/3];
weight = [1/3, 1/3, 1/3];
case 3 % Order 3, nQuad 4
lambda = [1/3, 1/3, 1/3; ...
0.6, 0.2, 0.2; ...
0.2, 0.6, 0.2; ...
0.2, 0.2, 0.6];
weight = [-27/48, 25/48, 25/48, 25/48];
case 4 % Order 4, nQuad 6
lambda = [0.108103018168070, 0.445948490915965, 0.445948490915965; ...
0.445948490915965, 0.108103018168070, 0.445948490915965; ...
0.445948490915965, 0.445948490915965, 0.108103018168070; ...
0.816847572980459, 0.091576213509771, 0.091576213509771; ...
0.091576213509771, 0.816847572980459, 0.091576213509771; ...
0.091576213509771, 0.091576213509771, 0.816847572980459];
weight = [0.223381589678011, 0.223381589678011, 0.223381589678011, ...
0.109951743655322, 0.109951743655322, 0.109951743655322];
case 5 % Order 5, nQuad 7
alpha1 = 0.059715871789770; beta1 = 0.470142064105115;
alpha2 = 0.797426985353087; beta2 = 0.101286507323456;
lambda = [ 1/3, 1/3, 1/3; ...
alpha1, beta1, beta1; ...
beta1, alpha1, beta1; ...
beta1, beta1, alpha1; ...
alpha2, beta2, beta2; ...
beta2, alpha2, beta2; ...
beta2, beta2, alpha2];
weight = [0.225, 0.132394152788506, 0.132394152788506, 0.132394152788506, ...
0.125939180544827, 0.125939180544827, 0.125939180544827];
case 6
A =[0.249286745170910 0.249286745170910 0.116786275726379
0.249286745170910 0.501426509658179 0.116786275726379
0.501426509658179 0.249286745170910 0.116786275726379
0.063089014491502 0.063089014491502 0.050844906370207
0.063089014491502 0.873821971016996 0.050844906370207
0.873821971016996 0.063089014491502 0.050844906370207
0.310352451033784 0.636502499121399 0.082851075618374
0.636502499121399 0.053145049844817 0.082851075618374
0.053145049844817 0.310352451033784 0.082851075618374
0.636502499121399 0.310352451033784 0.082851075618374
0.310352451033784 0.053145049844817 0.082851075618374
0.053145049844817 0.636502499121399 0.082851075618374];
lambda = [A(:,[1,2]), 1 - sum(A(:,[1,2]),2)];
weight = A(:,3)';
case 7
A =[0.333333333333333 0.333333333333333 -0.149570044467682
0.260345966079040 0.260345966079040 0.175615257433208
0.260345966079040 0.479308067841920 0.175615257433208
0.479308067841920 0.260345966079040 0.175615257433208
0.065130102902216 0.065130102902216 0.053347235608838
0.065130102902216 0.869739794195568 0.053347235608838
0.869739794195568 0.065130102902216 0.053347235608838
0.312865496004874 0.638444188569810 0.077113760890257
0.638444188569810 0.048690315425316 0.077113760890257
0.048690315425316 0.312865496004874 0.077113760890257
0.638444188569810 0.312865496004874 0.077113760890257
0.312865496004874 0.048690315425316 0.077113760890257
0.048690315425316 0.638444188569810 0.077113760890257];
lambda = [A(:,[1,2]), 1 - sum(A(:,[1,2]),2)];
weight = A(:,3)';
case 8
A =[0.333333333333333 0.333333333333333 0.144315607677787
0.081414823414554 0.459292588292723 0.095091634267285
0.459292588292723 0.081414823414554 0.095091634267285
0.459292588292723 0.459292588292723 0.095091634267285
0.658861384496480 0.170569307751760 0.103217370534718
0.170569307751760 0.658861384496480 0.103217370534718
0.170569307751760 0.170569307751760 0.103217370534718
0.898905543365938 0.050547228317031 0.032458497623198
0.050547228317031 0.898905543365938 0.032458497623198
0.050547228317031 0.050547228317031 0.032458497623198
0.008394777409958 0.263112829634638 0.027230314174435
0.008394777409958 0.728492392955404 0.027230314174435
0.263112829634638 0.008394777409958 0.027230314174435
0.728492392955404 0.008394777409958 0.027230314174435
0.263112829634638 0.728492392955404 0.027230314174435
0.728492392955404 0.263112829634638 0.027230314174435];
lambda = [A(:,[1,2]), 1 - sum(A(:,[1,2]),2)];
weight = A(:,3)';
case 9
A =[0.333333333333333 0.333333333333333 0.097135796282799
0.020634961602525 0.489682519198738 0.031334700227139
0.489682519198738 0.020634961602525 0.031334700227139
0.489682519198738 0.489682519198738 0.031334700227139
0.125820817014127 0.437089591492937 0.07782754100474
0.437089591492937 0.125820817014127 0.07782754100474
0.437089591492937 0.437089591492937 0.07782754100474
0.623592928761935 0.188203535619033 0.079647738927210
0.188203535619033 0.623592928761935 0.079647738927210
0.188203535619033 0.188203535619033 0.079647738927210
0.910540973211095 0.044729513394453 0.025577675658698
0.044729513394453 0.910540973211095 0.025577675658698
0.044729513394453 0.044729513394453 0.025577675658698
0.036838412054736 0.221962989160766 0.043283539377289
0.036838412054736 0.741198598784498 0.043283539377289
0.221962989160766 0.036838412054736 0.043283539377289
0.741198598784498 0.036838412054736 0.043283539377289
0.221962989160766 0.741198598784498 0.043283539377289
0.741198598784498 0.221962989160766 0.043283539377289];
lambda = [A(:,[1,2]), 1 - sum(A(:,[1,2]),2)];
weight = A(:,3)';
end
%% Verification
% The order of the quadrature rule is verified by the function
% verifyquadpts. See <matlab:ifem('verifyquadpts') verifyquadpts>.