Auto-solid angle calculation by Gararkin method D = vb_solid_auto_grk(D,F,Sout) 境界頂点係数行列の計算 自己重み係数：基底関数を三角面上で積分した寄与を加える (1次近似ポテンシャル) D : 境界頂点係数行列 (重み付き立体角) Sout : 三角面の面積法線 F : 境界三角面頂点インデックス 2004-02-06 Taku Yoshioka 2004-12-26 M. Sato modified Copyright (C) 2011, ATR All Rights Reserved. License : New BSD License(see VBMEG_LICENSE.txt)
0001 function D = vb_solid_auto_grk(D,F,Sout) 0002 % Auto-solid angle calculation by Gararkin method 0003 % D = vb_solid_auto_grk(D,F,Sout) 0004 % 境界頂点係数行列の計算 0005 % 自己重み係数：基底関数を三角面上で積分した寄与を加える 0006 % (1次近似ポテンシャル) 0007 % 0008 % D : 境界頂点係数行列 (重み付き立体角) 0009 % Sout : 三角面の面積法線 0010 % F : 境界三角面頂点インデックス 0011 % 0012 % 2004-02-06 Taku Yoshioka 0013 % 2004-12-26 M. Sato modified 0014 % 0015 % Copyright (C) 2011, ATR All Rights Reserved. 0016 % License : New BSD License(see VBMEG_LICENSE.txt) 0017 0018 % 0019 % 自己重み係数：基底関数を三角面上で積分した寄与 0020 % 非対角成分を含む 0021 % A : 境界頂点係数行列 (自己重み係数) 0022 0023 % 三角面の数 Number of triangle patch 0024 NF = size(F,1); 0025 NV = size(D,1); 0026 0027 S = sqrt(sum(Sout.^2,2)); 0028 0029 for i = 1:NF 0030 ix1 = F(i,1); 0031 ix2 = F(i,2); 0032 ix3 = F(i,3); 0033 S1 = S(i)/12; 0034 S2 = S(i)/24; 0035 0036 % A(ix1,ix1) = A(ix1,ix1) + S1; 0037 % A(ix2,ix2) = A(ix2,ix2) + S1; 0038 % A(ix3,ix3) = A(ix3,ix3) + S1; 0039 % 0040 % A(ix1,ix2) = A(ix1,ix2) + S2; 0041 % A(ix2,ix1) = A(ix2,ix1) + S2; 0042 % A(ix1,ix3) = A(ix1,ix3) + S2; 0043 % A(ix3,ix1) = A(ix3,ix1) + S2; 0044 % A(ix2,ix3) = A(ix2,ix3) + S2; 0045 % A(ix3,ix2) = A(ix3,ix2) + S2; 0046 0047 D(ix1,ix1) = D(ix1,ix1) + S1; 0048 D(ix2,ix2) = D(ix2,ix2) + S1; 0049 D(ix3,ix3) = D(ix3,ix3) + S1; 0050 0051 D(ix1,ix2) = D(ix1,ix2) + S2; 0052 D(ix2,ix1) = D(ix2,ix1) + S2; 0053 D(ix1,ix3) = D(ix1,ix3) + S2; 0054 D(ix3,ix1) = D(ix3,ix1) + S2; 0055 D(ix2,ix3) = D(ix2,ix3) + S2; 0056 D(ix3,ix2) = D(ix3,ix2) + S2; 0057 end 0058 0059 % This code save the memory 0060 for i=1:(NV*NV) 0061 D(i) = D(i) + (1/NV); 0062 end 0063 0064 % Following code may alocate new memory space for calculation 0065 % and this cause memory error for large number of triangles 0066 % D = D + (1/NV); 0067