Nine Point Quadrature - Level Set

RFA_LS · March 17, 2020, 16:24

Hi everyone

I found an equation of numerical integration over of control volume:

$\int _{\Omega_{ij}} g\approx \frac{h^2}{24}(16g_{ij}+g_{i-1,j-1}+g_{i-1,j} +g_{i-1,j+1}+g_{i,j-1}+g_{i,j+1} +g_{i+1,j-1}+g_{i+1,j}+g_{i+1,j+1})$

I am wondering why it is 16/24 and 1/24.
I would like to get the formula for 3D numerical integration, which should based on 3x3x3 points. What the weight should be?

The equation is from this paper (SIAM J. SCI. COMPUT. Vol. 20, No. 4, pp. 1165-1191) on page 1172 http://diyhpl.us/~bryan/papers2/frey...0algorithm.pdf

FMDenaro · March 17, 2020, 18:33

Quote:

Originally Posted by RFA_LS

Hi everyone

I found an equation of numerical integration over of control volume:

$\int _{\Omega_{ij}} g\approx \frac{h^2}{24}(16g_{ij}+g_{i-1,j-1}+g_{i-1,j} +g_{i-1,j+1}+g_{i,j-1}+g_{i,j+1} +g_{i+1,j-1}+g_{i+1,j}+g_{i+1,j+1})$

I am wondering why it is 16/24 and 1/24.
I would like to get the formula for 3D numerical integration, which should based on 3x3x3 points. What the weight should be?

The equation is from this paper (SIAM J. SCI. COMPUT. Vol. 20, No. 4, pp. 1165-1191) on page 1172 http://diyhpl.us/~bryan/papers2/frey...0algorithm.pdf

You can find the formula also in the Finite Volume chapter of the textbook Ferziger, Peric and Street

Eifoehn4 · March 17, 2020, 19:07

In one dimension it should be:

$\approx \frac{h}{6}(4g_{i}+g_{i-1}+g_{i+1})$

In two dimensions it should be:

$\approx \frac{h^2}{36}(16g_{ij}+4g_{i-1,j}+4g_{i,j-1}+4g_{i,j+1}+4g_{i+1,j} +g_{i-1,j-1}+g_{i-1,j+1}+g_{i+1,j-1}+g_{i+1,j+1})$

In three dimensions something like:

$\approx \frac{h^3}{216}\left(\begin{array}{ccccccccccccccccccc} 64 & 16 & 16 & 16 & 16 & 16 & 16 & 4 & 4 & 4 & 4 & 4 & 4 & 4 & \dots & 1 & 1 & 1 & \dots \end{array}\right)$

Regards

Eifoehn4 · March 17, 2020, 19:09

double post ...

FMDenaro · March 17, 2020, 19:33

Quote:

Originally Posted by Eifoehn4

In one dimension it should be:

$\approx \frac{h}{6}(4g_{i}+g_{i-1}+g_{i+1})$

In two dimensions it should be:

$\approx \frac{h^2}{36}(16g_{ij}+4g_{i-1,j}+4g_{i,j-1}+4g_{i,j+1}+4g_{i+1,j} +g_{i-1,j-1}+g_{i-1,j+1}+g_{i+1,j-1}+g_{i+1,j+1})$

In three dimensions something like:

$\approx \frac{h^3}{216}\left(\begin{array}{ccccccccccccccccccc} 64 & 16 & 16 & 16 & 16 & 16 & 16 & 4 & 4 & 4 & 4 & 4 & 4 & 4 & \dots & 1 & 1 & 1 & \dots \end{array}\right)$

Regards

Indeee it is just the multidimensional extension of the Simpson rule.

RFA_LS · March 17, 2020, 19:50

Thank you all!

The coefficients are different from the paper. Any reason for that?

RFA_LS · March 17, 2020, 19:53

Quote:

Originally Posted by Eifoehn4

double post ...

I saw that. But no one had answered and it was a long time ago.

Eifoehn4 · March 17, 2020, 19:55

Quote:

Originally Posted by RFA_LS

Thank you all!

The coefficients are different from the paper. Any reason for that?

There are at least three possibilities:

The paper is wrong!
You have not copied the formula correctly!
The integration rule is not based on aequidistant points including the boundaries!

I think it's point two!

RFA_LS · March 17, 2020, 20:02

Quote:

Originally Posted by Eifoehn4

There are at least three possibilities:

The paper is wrong!
You have not copied the formula correctly!
The integration rule is not based on aequidistant points including the boundaries!

I think it's point two!

It's correct.

You can check on page 1172:

http://diyhpl.us/~bryan/papers2/frey...0algorithm.pdf

Eifoehn4 · March 17, 2020, 20:20

Seems that the formula in the paper is wrong. Or they may use a generic weighting of the points.

Note that the only important fact is that the sum of all points contributing to integral must sum to unity. If the weights are different than the quadrature is defined according a different norm or point distribution.

Regards

RFA_LS · March 17, 2020, 22:32

Quote:

Originally Posted by Eifoehn4

In one dimension it should be:

$\approx \frac{h}{6}(4g_{i}+g_{i-1}+g_{i+1})$

In two dimensions it should be:

$\approx \frac{h^2}{36}(16g_{ij}+4g_{i-1,j}+4g_{i,j-1}+4g_{i,j+1}+4g_{i+1,j} +g_{i-1,j-1}+g_{i-1,j+1}+g_{i+1,j-1}+g_{i+1,j+1})$

In three dimensions something like:

$\approx \frac{h^3}{216}\left(\begin{array}{ccccccccccccccccccc} 64 & 16 & 16 & 16 & 16 & 16 & 16 & 4 & 4 & 4 & 4 & 4 & 4 & 4 & \dots & 1 & 1 & 1 & \dots \end{array}\right)$

Regards

Where did you find it? Can you suggest one book?

Eifoehn4 · March 18, 2020, 03:31

I have used my brain, starting from the one dimensional case. it's actually really simple.

You have to understand how the multi-D quadrature works on tensor-product elements.

FMDenaro · March 18, 2020, 03:51

The general rule is to write the quadratic polinomial in x,y,z by factored 1D lagrangian interpolation, then the integral is analytic. You get the formula for the 3x3x3 stencil.

praveen · March 18, 2020, 09:03

The integral in the paper is over $[x_i-h/2,x_i+h/2] \times [y_j-h/2,y_j+h/2]$ . The answers given by other people here are taking the domain to be $[x_i-h, x_i+h] \times [y_j-h,y_j+h]$ . Thats why you see a difference in the coefficients.

Eifoehn4 · March 18, 2020, 10:55

Quote:

Originally Posted by praveen

The integral in the paper is over $[x_i-h/2,x_i+h/2] \times [y_j-h/2,y_j+h/2]$ . The answers given by other people here are taking the domain to be $[x_i-h, x_i+h] \times [y_j-h,y_j+h]$ . Thats why you see a difference in the coefficients.

Dear praveen,

considering $[x_i-h/2,x_i+h/2] \times [y_j-h/2,y_j+h/2]$ with a quadratic polynomial.

This would result in one dimension

$\approx \frac{h}{24}(22g_{i}+g_{i-1}+g_{i+1})$

and in two dimensions

$\approx \frac{h^2}{576}(484g_{ij}+22g_{i-1,j}+22g_{i,j-1}+22g_{i,j+1}+22g_{i+1,j} +g_{i-1,j-1}+g_{i-1,j+1}+\dots).$

The integration rule mentioned in the paper is not well explained. Perhaps you can enlighten me. In my opinion it's something generic.

Regards

RFA_LS · March 18, 2020, 11:00

Quote:

Originally Posted by praveen

The integral in the paper is over $[x_i-h/2,x_i+h/2] \times [y_j-h/2,y_j+h/2]$ . The answers given by other people here are taking the domain to be $[x_i-h, x_i+h] \times [y_j-h,y_j+h]$ . Thats why you see a difference in the coefficients.

It doesn't explain the 24. Can you show?

Eifoehn4 · March 18, 2020, 11:04

Quote:

Originally Posted by RFA_LS

It doesn't explain the 24. Can you show?

Dear RFA_LS,

the 24 is only the result of the scaling. The sum over all DOF's must sum to unity.

RFA_LS · March 18, 2020, 11:48

I did using Simpson's rule for $x,y\pm h$ .

$\int^{1}_{-1}\int^{1}_{-1}f(x,y)dxdy=\int^{1}_{-1}\frac{\Delta x}{3}[f(x_{i-1},y)+4f(x_i,y)+f(x_{i+1},y)]dy-O(h^4) dy =$

$\int^{1}_{-1}\frac{\Delta x}{3}f(x_{i-1},y)dy+\frac{4\Delta x}{3}\int^{1}_{-1}f(x_i,y)dy+\frac{\Delta x}{3}\int^{1}_{-1}f(x_{i+1},y)dy-O(h^4) =$

$\frac{\Delta x\Delta y}{9}[f(x_{i-1},y_{i-1})+4f(x_{i-1},y_{i})+f(x_{i-1},y_{i+1})]+\frac{4\Delta x\Delta y}{9}[f(x_{i},y_{i-1})+$

$4f(x_{i},y_{i})+f(x_{i},y_{i+1})]+\frac{\Delta x\Delta y}{9}[f(x_{i+1},y_{i-1})+4f(x_{i+1},y_{i})+f(x_{i+1},y_{i+1})]$

Resulting:

$\int^{1}_{-1}\int^{1}_{-1}f(x,y)dxdy=\frac{\Delta x\Delta y}{9}[16f(x_{i},y_{i})+4f(x_{i-1},y_{i})+4f(x_{i},y_{i-1})+4f(x_{i+1},y_{i})$

$+4f(x_{i},y_{i+1})+f(x_{i-1},y_{i-1})+f(x_{i-1},y_{i+1})+f(x_{i+1},y_{i-1})+f(x_{i+1},y_{i+1})]$

I believe that changing from $x,y\pm h$ to $x,y\pm h/2$ it would change from 9 to 36 and from $\pm 1$ to $\pm 1/2$

praveen · March 18, 2020, 11:56

The formula in the paper is only exact for an affine function of the form

a + b x + c y

So it is not based on Simpson. There are many possibilities and its not obvious what the authors have done.

One suggestion: do a least squares fit of the nine values to the function

a + b (x - x_i) + c (y - y_j)

Then the integral is

a * h^2

So what is a ?

Eifoehn4 · March 18, 2020, 12:01

Yes I think we agree on that.

March 17, 2020, 16:24	Nine Point Quadrature - Level Set	#1
RFA_LS New Member Rodrigo Join Date: Mar 2020 Posts: 9 Rep Power: 6	Hi everyone I found an equation of numerical integration over of control volume: $\int _{\Omega_{ij}} g\approx \frac{h^2}{24}(16g_{ij}+g_{i-1,j-1}+g_{i-1,j} +g_{i-1,j+1}+g_{i,j-1}+g_{i,j+1} +g_{i+1,j-1}+g_{i+1,j}+g_{i+1,j+1})$ I am wondering why it is 16/24 and 1/24. I would like to get the formula for 3D numerical integration, which should based on 3x3x3 points. What the weight should be? The equation is from this paper (SIAM J. SCI. COMPUT. Vol. 20, No. 4, pp. 1165-1191) on page 1172 http://diyhpl.us/~bryan/papers2/frey...0algorithm.pdf

March 17, 2020, 19:07		#3
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 13	In one dimension it should be: $\approx \frac{h}{6}(4g_{i}+g_{i-1}+g_{i+1})$ In two dimensions it should be: $\approx \frac{h^2}{36}(16g_{ij}+4g_{i-1,j}+4g_{i,j-1}+4g_{i,j+1}+4g_{i+1,j} +g_{i-1,j-1}+g_{i-1,j+1}+g_{i+1,j-1}+g_{i+1,j+1})$ In three dimensions something like: $\approx \frac{h^3}{216}\left(\begin{array}{ccccccccccccccccccc} 64 & 16 & 16 & 16 & 16 & 16 & 16 & 4 & 4 & 4 & 4 & 4 & 4 & 4 & \dots & 1 & 1 & 1 & \dots \end{array}\right)$ Regards RFA_LS likes this.

March 18, 2020, 03:31		#12
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 13	I have used my brain, starting from the one dimensional case. it's actually really simple. You have to understand how the multi-D quadrature works on tensor-product elements. RFA_LS likes this.

March 18, 2020, 03:51		#13
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,783 Rep Power: 71	The general rule is to write the quadratic polinomial in x,y,z by factored 1D lagrangian interpolation, then the integral is analytic. You get the formula for the 3x3x3 stencil. RFA_LS likes this.

March 18, 2020, 09:03		#14
praveen Super Moderator Praveen. C Join Date: Mar 2009 Location: Bangalore Posts: 342 Blog Entries: 6 Rep Power: 18	The integral in the paper is over $[x_i-h/2,x_i+h/2] \times [y_j-h/2,y_j+h/2]$ . The answers given by other people here are taking the domain to be $[x_i-h, x_i+h] \times [y_j-h,y_j+h]$ . Thats why you see a difference in the coefficients. __________________ http://cpraveen.github.io http://twitter.com/cfdlab http://github.com/cpraveen

March 17, 2020, 19:09		#4
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 13	double post ...

March 17, 2020, 19:50		#6
RFA_LS New Member Rodrigo Join Date: Mar 2020 Posts: 9 Rep Power: 6	Thank you all! The coefficients are different from the paper. Any reason for that?

March 17, 2020, 20:20		#10
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 13	Seems that the formula in the paper is wrong. Or they may use a generic weighting of the points. Note that the only important fact is that the sum of all points contributing to integral must sum to unity. If the weights are different than the quadrature is defined according a different norm or point distribution. Regards

March 18, 2020, 11:48		#18
RFA_LS New Member Rodrigo Join Date: Mar 2020 Posts: 9 Rep Power: 6	I did using Simpson's rule for $x,y\pm h$ . $\int^{1}_{-1}\int^{1}_{-1}f(x,y)dxdy=\int^{1}_{-1}\frac{\Delta x}{3}[f(x_{i-1},y)+4f(x_i,y)+f(x_{i+1},y)]dy-O(h^4) dy =$ $\int^{1}_{-1}\frac{\Delta x}{3}f(x_{i-1},y)dy+\frac{4\Delta x}{3}\int^{1}_{-1}f(x_i,y)dy+\frac{\Delta x}{3}\int^{1}_{-1}f(x_{i+1},y)dy-O(h^4) =$ $\frac{\Delta x\Delta y}{9}[f(x_{i-1},y_{i-1})+4f(x_{i-1},y_{i})+f(x_{i-1},y_{i+1})]+\frac{4\Delta x\Delta y}{9}[f(x_{i},y_{i-1})+$ $4f(x_{i},y_{i})+f(x_{i},y_{i+1})]+\frac{\Delta x\Delta y}{9}[f(x_{i+1},y_{i-1})+4f(x_{i+1},y_{i})+f(x_{i+1},y_{i+1})]$ Resulting: $\int^{1}_{-1}\int^{1}_{-1}f(x,y)dxdy=\frac{\Delta x\Delta y}{9}[16f(x_{i},y_{i})+4f(x_{i-1},y_{i})+4f(x_{i},y_{i-1})+4f(x_{i+1},y_{i})$ $+4f(x_{i},y_{i+1})+f(x_{i-1},y_{i-1})+f(x_{i-1},y_{i+1})+f(x_{i+1},y_{i-1})+f(x_{i+1},y_{i+1})]$ I believe that changing from $x,y\pm h$ to $x,y\pm h/2$ it would change from 9 to 36 and from $\pm 1$ to $\pm 1/2$

March 18, 2020, 11:56		#19
praveen Super Moderator Praveen. C Join Date: Mar 2009 Location: Bangalore Posts: 342 Blog Entries: 6 Rep Power: 18	The formula in the paper is only exact for an affine function of the form a + b x + c y So it is not based on Simpson. There are many possibilities and its not obvious what the authors have done. One suggestion: do a least squares fit of the nine values to the function a + b (x - x_i) + c (y - y_j) Then the integral is a * h^2 So what is a ? RFA_LS likes this. __________________ http://cpraveen.github.io http://twitter.com/cfdlab http://github.com/cpraveen

March 18, 2020, 12:01		#20
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 13	Yes I think we agree on that.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[snappyHexMesh] meshing of very small patches in comparison with the overall geometry	christos	OpenFOAM Meshing & Mesh Conversion	3	December 17, 2014 16:55
[snappyHexMesh] Adding layers goes wrong with SnappyHexMesh	Elise	OpenFOAM Meshing & Mesh Conversion	1	April 22, 2013 02:32
the initialization of the coupled level set method and vof model in FLUENT 13.0	dutliang	FLUENT	3	December 13, 2011 21:28
Problem on high density ratio in Level Set method	Kai Yan	Main CFD Forum	10	December 25, 2007 06:12
[Gmsh] Gmsh and samplesurface	touf	OpenFOAM Meshing & Mesh Conversion	2	December 10, 2007 02:27