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 12

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 12

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 12

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 12

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 12

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 12

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 12

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 12

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 12

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 12

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 12

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.

Subject: volume integration

From: Liaquat Ali

Date: 7 May, 2014 23:25:10

Message: 12 of 12

"Marios Karaoulis" <marios_hellas_23@yahoo.gr> wrote in message <k9ojg4$pk6$1@newscl01ah.mathworks.com>...
> 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.

Hi Marios,

I have the same problem i get two vector fields D and F in 3D spcae(x,y,z) now firstly i need to normalized the D then calculate the volume integral of the D and F.

I have the data from COMSOL in this format
x y z dx dy dz fx fy fz

Kindly tell me how to do that integral

Tags for 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