Incorrect integral

[stla]

Hello,

In this blog post, using R, I calculate an integral on a polyhedron by decomposing this polyhedron into simplices.

The Julia code below is my attempt to reproduce this R example with the help of your package. But I don't find the same result. However, using another method for the calculation of the integrals, I find the same result as R.

So I believe there's a mistake or a bug in your package, or there's something bad in my code but I don't see what.

using Polyhedra
using MiniQhull # to use delaunay()
using GrundmannMoeller

# define the polyhedron Ax <= b and get its vertices
# see https://laustep.github.io/stlahblog/posts/integralPolyhedron.html
A = [
    -1  0  0;
     1  0  0;
     0 -1  0;
     1  1  0;
     0  0 -1;
     1  1  1.0
]
b = [5; 4; 5; 3; 10; 6.0]
H = hrep(A, b)
P = polyhedron(H)
pts = points(P)
vertices = collect(pts)

# decompose the polyhedron into simplices (tetrahedra)
indices = delaunay(hcat(vertices...))
_, ntetrahedra = size(indices)
tetrahedra = Vector{Vector{Vector{Float64}}}(undef, ntetrahedra)
for j in 1:ntetrahedra
    ids = vec(indices[:, j])
    tetrahedra[j] = vertices[ids]
end

# function to be integrated
function f(x)
    return x[1] + x[2]*x[3]
end
# integral is sum of integrals over the tetrahedra
scheme = grundmann_moeller(Float64, Val(3), 5)
integral = 0
for i in 1:ntetrahedra
    global integral = integral + integrate(f, scheme, tetrahedra[i])
end
integral
# -11879.7
# should be -5358.3

[stla]

Here is the calculation of the volume of the 3-dimensional unit simplex. It should be 1/6, but I get -1/6.

using GrundmannMoeller

tetrahedron = [[1.0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]]

# function to be integrated
function f(x)
    return 1
end
# integral 
scheme = grundmann_moeller(Float64, Val(3), 5)
integrate(f, scheme, tetrahedron)

[stla]

Indeed, the integrals are always the same as mine up to the sign.

If I change the order of the vertices of tetrahedron, I get 1/6 as expected. Does your package take into account the orientation of the simplices?

[stla]

I think I found the problem: in the function calc_vol_vertices, one should take the absolute value of the determinant.