Search All of the Math Forum:

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

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: volume integration
Replies: 11   Last Post: May 7, 2014 7:25 PM

 Messages: [ Previous | Next ]
 Bruno Luong Posts: 9,822 Registered: 7/26/08
Re: volume integration
Posted: Sep 9, 2011 2:52 PM

Sorry I clean up few pollution 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));
% Determinant
D = A(:,1).*m1 - A(:,4).*m2 + A(:,7).*m3;

JAJB = dot(JA,JB,3);
JAJB = mean(JAJB,2);
jac = JAJB.*abs(D)/6;

I = sum(jac)

% Bruno

Date Subject Author
9/6/11 Marios Karaoulis
9/6/11 Bruno Luong
9/6/11 Marios Karaoulis
9/6/11 Bruno Luong
9/9/11 Marios Karaoulis
9/9/11 Bruno Luong
9/9/11 Bruno Luong
9/9/11 Marios Karaoulis
9/10/11 Bruno Luong
12/5/12 Marios
5/7/14 Liaquat Ali
9/8/11 Marios Karaoulis