Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
volume integration

Subject: volume integration

From: Marios Karaoulis

Date: 6 Sep, 2011 18:33:38

Message: 1 of 11

Hi all,

I have two vector field J1 and J2, in 3D space (x,y,z) and calculated
at some discrete points (extracted in a txt file from comsol.

I need to volume integrate the
triple_integration of (J1 dot J2), which is expressed as

triple_integration ( J1(x)*J2(x) + J1(y)*J2(y) + J1(z)*J2(z) ).


I have found that http://www-users.math.umd.edu/~jmr/241/tripleint.html
, but in my case, I have no analytically expression of the funstion.
I guess I could interpolate using TriScatinterp, but is any other way
to do that?

Subject: volume integration

From: Bruno Luong

Date: 6 Sep, 2011 19:04:27

Message: 2 of 11

Marios Karaoulis <marios.karaoulis@gmail.com> wrote in message <8e295512-97e5-44e1-a95f-5e59af177653@o10g2000vby.googlegroups.com>...
> Hi all,
>
> I have two vector field J1 and J2, in 3D space (x,y,z) and calculated
> at some discrete points (extracted in a txt file from comsol.
>
> I need to volume integrate the
> triple_integration of (J1 dot J2), which is expressed as
>
> triple_integration ( J1(x)*J2(x) + J1(y)*J2(y) + J1(z)*J2(z) ).
>
>
> I have found that http://www-users.math.umd.edu/~jmr/241/tripleint.html
> , but in my case, I have no analytically expression of the funstion.
> I guess I could interpolate using TriScatinterp, but is any other way
> to do that?

You could you partition the volume by Delaunay tetrahedron (help DelaunayTri). Over each tetrahedron Tj, approximate the integral by:

Ij = |Tj|/4 sum(i) dot(J1,J2)(xi)

where {xi}={x1,x2,x3,x4} are four corners of Tj. and |Tj| is the volume

|Tj| = 1/6*abs(det [x2-x1 x3-x2 x4-x2 ]).

Then sum Ij for all j to get the integral on the whole volume.

Bruno

Subject: volume integration

From: Marios Karaoulis

Date: 6 Sep, 2011 21:51:26

Message: 3 of 11

On Sep 6, 1:04 pm, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Marios Karaoulis <marios.karaou...@gmail.com> wrote in message <8e295512-97e5-44e1-a95f-5e59af177...@o10g2000vby.googlegroups.com>...
> > Hi all,
>
> > I have two vector field J1 and J2, in 3D space (x,y,z) and calculated
> > at some discrete points (extracted in a txt file from comsol.
>
> > I need to volume integrate the
> > triple_integration of (J1 dot J2), which is expressed as
>
> > triple_integration ( J1(x)*J2(x) + J1(y)*J2(y) + J1(z)*J2(z) ).
>
> > I have found thathttp://www-users.math.umd.edu/~jmr/241/tripleint.html
> > , but in my case, I have no analytically expression of the funstion.
> > I guess I could interpolate using TriScatinterp, but is any other way
> > to do that?
>
> You could you partition the volume by Delaunay tetrahedron (help DelaunayTri). Over each tetrahedron Tj, approximate the integral by:
>
> Ij = |Tj|/4 sum(i) dot(J1,J2)(xi)
>
> where {xi}={x1,x2,x3,x4} are four corners of Tj. and |Tj| is the volume
>
> |Tj| = 1/6*abs(det [x2-x1 x3-x2 x4-x2 ]).
>
> Then sum Ij for all j to get the integral on the whole volume.
>
> Bruno



Thanks for the reply. I am not sure I get this, but I will give it a
shot.
By DelaunayTri, you create triangles. Should't I need delaunayn ?

Best

Subject: volume integration

From: Bruno Luong

Date: 6 Sep, 2011 22:25:28

Message: 4 of 11

Marios Karaoulis <marios.karaoulis@gmail.com> wrote in message <5210aa44-a37a-435e-b109-177774624018@h14g2000yqi.googlegroups.com>...

> By DelaunayTri, you create triangles. Should't I need delaunayn ?

In 3d DelaunayTri returns tetrahedra.

Bruno

Subject: volume integration

From: Marios Karaoulis

Date: 8 Sep, 2011 21:52:39

Message: 5 of 11

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

Subject: volume integration

From: Marios Karaoulis

Date: 9 Sep, 2011 16:04:24

Message: 6 of 11

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

Subject: volume integration

From: Bruno Luong

Date: 9 Sep, 2011 18:37:14

Message: 7 of 11

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

Subject: volume integration

From: Bruno Luong

Date: 9 Sep, 2011 18:52:29

Message: 8 of 11

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

Subject: volume integration

From: Marios Karaoulis

Date: 9 Sep, 2011 19:16:28

Message: 9 of 11

Thanks so much Bruno. You blew my mind.

One last question,

For every node, actually I need to evaluate the integral 4 times.
Using the same notation,
I have
JA JB JM JN

The integration is
triple integral ( (JAJM -JAJN -JBM +JBN))

So, if I change the last part to something like this

JAJM = dot(JA,JM,3);
JAJM = mean(JAJM,2);

JAJN = dot(JA,JN,3);
JAJN = mean(JAJN,2);


JBJM = dot(JB,JM,3);
JBJM = mean(JBJM,2);


JBJN = dot(JB,JN,3);
JBJN = mean(JBJN,2);


J_ALL=(JAJM-JAJN-JBJM+JBJN)

jac2 = J_ALL.*abs(D)/6;


it should be alright.

But is possible to write this in the following way

J_ALL=dot(JA-JB,JM-JN,3)
J_ALL=mean(J_ALL,2)

thanks in advance.

Subject: volume integration

From: Bruno Luong

Date: 10 Sep, 2011 08:33:11

Message: 10 of 11

Marios Karaoulis <marios.karaoulis@gmail.com> wrote in message <2d7de568-c9f7-4dbf-9eb1-a5c4571a8897@m18g2000vbe.googlegroups.com>...
>
> But is possible to write this in the following way
>
> J_ALL=dot(JA-JB,JM-JN,3)
> J_ALL=mean(J_ALL,2)

This second way looks shorter and better in all point of view IMO.

Bruno

Subject: volume integration

From: Marios Karaoulis

Date: 5 Dec, 2012 22:57:09

Message: 11 of 11

Hi Bruno,

I need your assistance again.
I have the same problem again but in 2D.
I have created the triangles, and need to area integrate over each triangle.
Could you please helo me vectrotizing this?

Thanks in advance.

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us