|
|
Re: volume integration
Posted:
Sep 9, 2011 12:04 PM
|
|
On Sep 6, 4:25 pm, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Marios Karaoulis <marios.karaou...@gmail.com> wrote in message <5210aa44-a37a-435e-b109-177774624...@h14g2000yqi.googlegroups.com>...
> > By DelaunayTri, you create triangles. Should't I need delaunayn ?
>
> In 3d DelaunayTri returns tetrahedra.
>
> Bruno
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
|
|
