Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: volume integration
Replies: 10   Last Post: Dec 5, 2012 5:57 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Bruno Luong

Posts: 8,572
Registered: 7/26/08
Re: volume integration
Posted: Sep 9, 2011 2:37 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.