|
|
Re: volume integration
Posted:
Sep 9, 2011 2:37 PM
|
|
You should try to vectorize the code rather than using PAR-FOR on a slow code.
xyz = unx(:,1:3);
JA = unx(:,4:6);
JB = unx(:,7:9);
T = DelaunayTri(xyz);
T = T.Triangulation;
JA = reshape(JA(T,:),[],4,3);
JB = reshape(JB(T,:),[],4,3);
Tetra = reshape(xyz(T,:),[],4,3);
Tetra = permute(Tetra, [2 3 1]);
A = cat(1,Tetra(2,:,:)-Tetra(1,:,:), ...
Tetra(3,:,:)-Tetra(2,:,:), ...
Tetra(4,:,:)-Tetra(3,:,:));
A = reshape(A, 9, []).';
% This code compute determinant hacked from here:
% http://www.mathworks.com/matlabcentral/fileexchange/27762-small-size-linear-solver/content/SmallSizeSolver/inv3.m
% minors
m1 = (A(:,5).*A(:,9) - A(:,8).*A(:,6));
m2 = (A(:,2).*A(:,9) - A(:,8).*A(:,3));
m3 = (A(:,2).*A(:,6) - A(:,5).*A(:,3));
X = [+m1 ...
-m2 ...
+m3 ...
-(A(:,4).*A(:,9)-A(:,7).*A(:,6)) ...
+(A(:,1).*A(:,9)-A(:,7).*A(:,3)) ...
-(A(:,1).*A(:,6)-A(:,4).*A(:,3)) ...
+(A(:,4).*A(:,8)-A(:,7).*A(:,5)) ...
-(A(:,1).*A(:,8)-A(:,7).*A(:,2)) ...
+(A(:,1).*A(:,5)-A(:,4).*A(:,2))];
% Determinant
D = A(:,1).*m1 - A(:,4).*m2 + A(:,7).*m3;
D = abs(D);
JAJB = dot(JA,JB,3);
JAJB = mean(JAJB,2);
jac = JAJB.*abs(D)/6;
I = sum(jac)
% Bruno
|
|
