书接上回《高斯-勒朗德积分公式》
需求：在给定空间三角形\(\Delta ABC\)中，\(A(x_1,y_1,z_1),B(x_2,y_2,z_2),C(x_3,y_3,z_3)\)，已知函数\(f(x,y,z)\)，求利用数值方法求解积分：\(\iint_{\Delta ABC}f(x,y,z)\text dS\)。
解决方法：参考triangle_lyness_rule给出的积分方法，具体细节也不是太懂，但是思路上与高斯积分类似，计算平面上的积分点与系数权重进行积分
triangle_lyness_rule中给出了不同阶数下，在标准三角形中的系数点位置与权重系数。例如下图中展示\(Rule=10\)时的积分点位置与权重系数。下表中显示了不同\(Rule\)下的积分精度\(Precision\)积分点数目\(order\)以及积分点是否包含三角形中心\(center\)。
| Rule | Order | Precision | Center |
|---|---|---|---|
| 0 | 1 | 1 | YES |
| 1 | 3 | 2 | NO |
| 2 | 4 | 2 | YES |
| 3 | 4 | 3 | YES |
| 4 | 7 | 3 | YES |
| 5 | 6 | 4 | NO |
| 6 | 10 | 4 | YES |
| 7 | 9 | 4 | NO |
| 8 | 7 | 5 | YES |
| 9 | 10 | 5 | YES |
| 10 | 12 | 6 | NO |
| 11 | 16 | 6 | YES |
| 12 | 13 | 6 | YES |
| 13 | 13 | 7 | YES |
| 14 | 16 | 7 | YES |
| 15 | 16 | 8 | YES |
| 16 | 21 | 8 | NO |
| 17 | 16 | 8 | YES |
| 18 | 19 | 9 | YES |
| 19 | 22 | 9 | YES |
| 20 | 27 | 11 | NO |
| 21 | 28 | 11 | YES |
采用与之前文章中《高斯-勒朗德积分公式》形函数方式计算坐标转换关系。得到三节点形函数为，剩下步骤与《高斯-勒朗德积分公式》中类似。
在KY师兄指点下，以上步骤可以进一步简化。原因在于三角形坐标变换的形函数简单，可以直接进行坐标运算，Jacobi系数直接等于三角形面积，具体见代码。
以下为测试积分函数，其中\(LYNESS_RULE.txt\)存储的数据太长了，就放到Gitee:链接待更新仓库了。
%% 测试三角形积分clc;clear;global TriCoeff% 导入积分系数TriCoeff=loadLynessFromTxT("LYNESS_RULE.txt");P1=[0,0,0];P2=[2,0,0];P3=[0,3,0];% 积分函数func=@(x,y,z) (x^6+y^3+1);count=1;for rule=0:1:21 [P_W] = getTrianglePoints([P1;P2;P3],rule); [N,~]=size(P_W); res=0; for i=1:1:N res=res+func(P_W(i,1),P_W(i,2),P_W(i,3))*P_W(i,4); end resA(count,1)=res; resA(count,2)=rule; resA(count,3)=N; count=count+1;end%% matlab 自带积分函数pfun = @(x,y) (x.^6+y.^3+1);xmin = 0;xmax = 2;ymin = 0;ymax = @(x) -3/2*x+3;r = integral2(pfun,xmin,xmax,ymin,ymax);%% plotfigure(22)plot(resA(:,2),resA(:,1),'r-o');hold on;plot(resA(:,2),r*ones(22,1),'b-');grid on;xticks([0:2:22]);xlim([0,22]);legend("TRIANGLE LYNESS RULE 积分","Matlab integral2积分");text(10,14,"积分函数：(x^6+y^3+1)")text(10,12,"积分区域：(0,0,0),(2,0,0),(0,3,0)");xlabel("Lyness Rule");ylabel("积分数值");积分结果对比
function [P_W] = getTrianglePoints(Triangle,Rule)% getTrianglePoints 三角形面元积分% https://people.sc.fsu.edu/~jburkardt/cpp_src/triangle_lyness_rule/triangle_lyness_rule.html% 输入：% Triangle(3,3)：三角形面元三个点% Rule：triangle_lyness_rule% 输出：% P_W(:,4)：P_W(:,1:3)积分点、P_W(:,4)权重系数%% 任意空间三角形 =》平面直角三角形 坐标转换% 形函数N1=@(s,t) -s-t+1;N2=@(s,t) s;N3=@(s,t) t;N1_s=@(s,t) -1;N2_s=@(s,t) 1;N3_s=@(s,t) 0;N1_t=@(s,t) -1;N2_t=@(s,t) 0;N3_t=@(s,t) 1;P1=Triangle(1,:);P2=Triangle(2,:);P3=Triangle(3,:);global TriCoeff;data=TriCoeff{Rule+1,1};[order,~]=size(data);P_W=zeros(order,4);for i=1:1:order P_W(i,1:3)=Loc2Glo(data(i,1:2)); P_W(i,4)=data(i,3)*Jacobi(data(i,1:2));end function Pglobal=Loc2Glo(loc) %loc(1,2) Pglobal=N1(loc(1),loc(2))*P1+... N2(loc(1),loc(2))*P2+... N3(loc(1),loc(2))*P3; end function J=Jacobi(Loc) s=N1_s(Loc(1),Loc(2))*P1+... N2_s(Loc(1),Loc(2))*P2+... N3_s(Loc(1),Loc(2))*P3; t=N1_t(Loc(1),Loc(2))*P1+... N2_t(Loc(1),Loc(2))*P2+... N3_t(Loc(1),Loc(2))*P3; %三角形，这里多除了一个2 J=norm(cross(s,t))/2; endendfunction [P_W] = getTrianglePointsSimplified(Triangle,Rule) global TriCoeff; points = TriCoeff{Rule+1}; weights = points(:,3)'; points(:,3) = 1-points(:,1)-points(:,2); P_W = zeros(size(points,1),4); P_W(:,1:3)=points*Triangle; area = 0.5*norm(cross(Triangle(1,:)-Triangle(2,:),Triangle(1,:)-Triangle(3,:))); P_W(:,4)=weights*area;endfunction TriCoeff = loadLynessFromTxT(filename)%LOADLYNESSFROMTXT 加载系数TriCoeff=cell(22,1);fp=fopen(filename,'r');data=textscan(fp,"%f,%f,%f");fclose(fp);ALL=[data{1,1},data{1,2},data{1,3}];[N,~]=size(ALL);i=1;while i<N rule=ALL(i,1); order=ALL(i,2); TriCoeff{rule+1,1}=ALL(i+1:i+order,:); i=i+order+1;endend