P03_非凹四边形面积计算、截断式四棱柱体积计算

非凹四边形面积计算，可分成两个三角形，三角形计算面积可以用两向量叉乘的模一半来表示（规则四边形两边的叉乘的模就是其面积）。

二截断式四棱柱可以分为两个截断式三菱柱计算，每个截断式三菱柱可以分为一个三菱柱+两个三菱锥（或者一个三菱柱+一个四凌锥）

计算案例区域是由一系列节点网格组成的，文本信息由节点信息（ N(x y z) ) 和面元素信息( E(1 2 3 4) )组成，计算蓝色网格的面积，及其对应的体积，求z的平均值。

---

%S_mesh_quad是每个四边形的面积，S_sum_quad是总网格的面积
%V_mesh_quad是每个类四棱柱的体积，V_sum_quad是总体积

clc,clear

%写入节点信息
N = importdata('Note.dat');
E = importdata('Element.dat');
%写入完毕

S_mesh_quad = zeros(size(E, 1), 1);
V_mesh_quad = zeros(size(E, 1), 1);

for i_sv = 1:size(E, 1)
    N_quad = N(E(i_sv, :), :);
    %拆分为两个三角形tri1(241)和tri2（243）来计算XOY平面的四边形网格面积，tecplot四边形云图是以2、4节点连接的
    N_quadxy=N_quad;
    N_quadxy(:,3)=0;
    L24_tri = N_quadxy(4, :) - N_quadxy(2, :);
    L21_tri = N_quadxy(1, :) - N_quadxy(2, :);
    L23_tri = N_quadxy(3, :) - N_quadxy(2, :);
    S_tri1 = 0.5 * norm(cross(L24_tri, L21_tri));
    S_tri2 = 0.5 * norm(cross(L24_tri, L23_tri));
    S_quad=S_tri1+S_tri2;

    %计算第一个三角形的投影体积，拆分为一个三菱柱column和两个三菱锥cone1、cone2
    N_tri1=N_quad;
    N_tri1(3,:)=[];     %第一个三角形去掉第三点数据
    N_tri1=sortrows(N_tri1,3);  %按z值对点数据排列，利于计算体积，最小z值确定三菱柱体积，2、3、4z值确定两个三菱锥的底面

    tri1_cone2=[N_tri1;N_tri1(3,:)];
    tri1_cone2(4,3)=N_tri1(1,3);       %顶部三菱锥点位信息，第四点要把z值设置到最小值

    tri1_cone2_L12=tri1_cone2(2,:)-tri1_cone2(1,:);
    tri1_cone2_L13=tri1_cone2(3,:)-tri1_cone2(1,:);
    tri1_cone2_L14=tri1_cone2(4,:)-tri1_cone2(1,:);  %最底部三菱锥的三条边，便可求此三菱锥面积

    V_tri1_cone1=S_tri1*(N_tri1(2,3)-N_tri1(1,3))/3;
    V_tri1_cone2=abs(dot(tri1_cone2_L14,cross(tri1_cone2_L12,tri1_cone2_L13))/6);  %底面积的法向量点积另外一条向量再/2/3
    V_tri1=S_tri1*N_tri1(1,3)+V_tri1_cone1+V_tri1_cone2;  %第一部分类三菱柱总体积为V_column+V_tri1_cone1+V_tri1_cone2

    N_tri2=N_quad;
    N_tri2(1,:)=[];     %第二个三角形去掉第一点数据
    N_tri2=sortrows(N_tri2,3);  %按z值对点数据排列，利于计算体积

    tri2_cone2=[N_tri2;N_tri2(3,:)];
    tri2_cone2(4,3)=N_tri2(1,3);       %顶部三菱锥点位信息，第四点要把z值降到最小值

    tri2_cone2_L12=tri2_cone2(2,:)-tri2_cone2(1,:);
    tri2_cone2_L13=tri2_cone2(3,:)-tri2_cone2(1,:);
    tri2_cone2_L14=tri2_cone2(4,:)-tri2_cone2(1,:);  %最底部三菱锥的三条边，便可求此三菱锥面积

    V_tri2_cone1=S_tri2*(N_tri2(2,3)-N_tri2(1,3))/3;
    V_tri2_cone2=abs(dot(tri2_cone2_L14,cross(tri2_cone2_L12,tri2_cone2_L13))/6);  %底面积的法向量点积另外一条向量再/2/3
    V_tri2=S_tri2*N_tri2(1,3)+V_tri2_cone1+V_tri2_cone2;  %第一部分类三菱柱总体积为V_column+V_tri2_cone1+V_tri2_cone2
    V_quad=V_tri1+V_tri2;

    %写入单个网格的面积及投影体积
    S_mesh_quad(i_sv,1)=S_quad;
    V_mesh_quad(i_sv,1)=V_quad;

end

%计算总面积、总体积
S_sum_quad=sum(S_mesh_quad(:,1));
V_sum_quad=sum(V_mesh_quad(:,1));
A_average=V_sum_quad/S_sum_quad;
%计算完毕

%写入总面积、总体积