|
Re: volume integration
Posted:
Sep 8, 2011 5:52 PM
|
|
Hi again,
Let's say I have the unx matrix where
1st column is the x-coordinate 2nd column is the y-coordinate 3rd column is the z=coordinate
4th column is the x component of J1 5th column is the y component of J1 6th column is the z component of J1
7th column is the x component of J2 8th column is the y component of J2 9th column is the z component of J2
I create the tetrahedrals DT=DelaunayTri(unx(:,1),unx(:,2),unx(:,3)); num_tet=length(DT(:,1));
% find volume
% find volume of each tetrahedra voule_tet=zeros(num_tet,1); parfor i=1:num_tet
n1=DT(i,1); n2=DT(i,2); n3=DT(i,3); n4=DT(i,4);
P1=[unx(n1,1) unx(n1,2) unx(n1,3)]; P2=[unx(n2,1) unx(n2,2) unx(n2,3)]; P3=[unx(n3,1) unx(n3,2) unx(n3,3)]; P4=[unx(n4,1) unx(n4,2) unx(n4,3)];
V1=P2-P1; V2=P3-P2; V3=P4-P3;
volume_tet(i) = abs(det([V1(:) V2(:) V3(:)]))/6;
end
% find integration
jac=zeros(1,num_tet);
parfor k=1:num_tet % take appropriate solution
n1=DT(k,1); n2=DT(k,2); n3=DT(k,3); n4=DT(k,4);
%take from A (Jx,Jy,Jz) for each 4 points JA1=[unx(n1,4:6)]; JA2=[unx(n2,4:6)]; JA3=[unx(n3,4:6)]; JA4=[unx(n4,4:6)];
JB1=[unx(n1,7:9)]; JB2=[unx(n2,7:9)]; JB3=[unx(n3,7:9)]; JB4=[unx(n4,7:9)];
jac(1,k)=0.25*volume(k)*(dot(JA1,JB1)+dot(JA2,JB2)+dot(JA3,JB3)+dot(JA4,JB4));
end
Does this look right? Do you have any suggestions how to speed this up?
Thanks in adnvance
|
|