# Compute integral for DelaunayTriangulation.jl to be able to solve PDE with FiniteVolumeMethod.jl

**URL:** <https://discourse.julialang.org/t/compute-integral-for-delaunaytriangulation-jl-to-be-able-to-solve-pde-with-finitevolumemethod-jl/124310>\
**Category:** Numerics\
**Tags:** pde, numerics, mesh, integral\
**Created:** [December 31, 2024, 1:33pm UTC](https://discourse.julialang.org/t/compute-integral-for-delaunaytriangulation-jl-to-be-able-to-solve-pde-with-finitevolumemethod-jl/124310 "2024-12-31T13:33:10Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![B4rtel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/b4rtel/32/214171_2.png) [@B4rtel](https://discourse.julialang.org/u/B4rtel)\
**Post date:** [December 31, 2024, 1:33pm UTC](https://discourse.julialang.org/t/compute-integral-for-delaunaytriangulation-jl-to-be-able-to-solve-pde-with-finitevolumemethod-jl/124310/1 "2024-12-31T13:33:10Z")

</div>

Hi :), I would like to solve a advection PDE of the form \frac{\partial}{\partial t} f(x,t) = - \nabla\_x \cdot \left(\left(\int [\phi(p\_1(x,y),\phi(p\_2(x,y)]\odot (y-x) f(y,t) \text{d} y\right) f(x,t)\right) for p\_a(x\_i,x\_j):= \beta |x\_a-y\_a| + (1-\beta)\sum\_{b=1, b\neq a}^2\alpha\_b |x\_b-y\_b|, \phi a scalar function and x,y \in R^2, with the finite volume method by using the FiniteVolumeMethod.jl package. To be able to do that I need to compute the integral. Is there a way to do that when I discretise the space with the DelaunyTriangulation.jl package as

```julia
tri = triangulate_rectangle(a, b, c, d, nx, ny, single_boundary=true)
mesh = FVMGeometry(tri)

```

(which is the format I need it to be discretised with to then use FiniteVolumeMethod.jl to solve the PDE).  
Here you find the code I have so far (which is based on the example of FiniteVolumeMethod.jl) where I marked what is missing.

```julia
using FiniteVolumeMethod, DelaunayTriangulation, CairoMakie, OrdinaryDiffEq

N=64
T=2.0
xmin = -1.0
xmax = 1.0

# choose parameters
const beta = 0.5
const alpha_fixed = [0.5, 0.5]

# make the mesh
a, b, c, d = xmin, xmax, xmin, xmax
nx, ny = N, N
tri = triangulate_rectangle(a, b, c, d, nx, ny, single_boundary=true)
mesh = FVMGeometry(tri)
# set boundary conditions
bc = (x, y, t, u, p) -> zero(u)
BCs = BoundaryConditions(mesh, bc, Dirichlet)
# initialise f
f = (x, y) -> y ≤ 0.0 ? 50.0 : 0.0
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]

# Define the function φ
ϕ_choice = "ifelse_r_less_3_10"
ϕ(r) = ifelse(r < 3.0/10, 1.0, 0.0) # this function returns 1 if r is less than 3/10 and 0 otherwise

# Define the p_a norm
function p_a(x_i, x_j, alpha_i, beta, a)
    return beta * abs(x_i[a] - x_j[a]) + (1 - beta) * sum(alpha_i[b] * abs(x_i[b] - x_j[b]) for b in 1:2 if b != a)
end

function integr(f, x_fix, y_fix, tri, integ, ϕ_vec, beta, alpha_fixed)
  THIS PART IS MISSING:
  It shoud return the numerically computerd integral over all points (x,y) in the mesh of the function
    ϕ(p_a([x_fix, y_fix],[x,y], alpha_fixed, beta, a)) *([x,y]-[x_fix, y_fix])* f(x, y) dx dy
    where p_a and ϕ are defined above
end

flux_function = (x, y, t, α, β, γ, p) -> begin
    u = α * x + β * y + γ
    ν = integr(u, x, y, p.tri, p.integ, p.ϕ_vec, beta, alpha_fixed)
    qx = p.ν * u
    qy = 0
    return (qx, qy)
end
# initialise integ and \phi_vec (that will be used in the integration) here for faster runtime
integ = zeros(2)
ϕ_vec = zeros(2)

# Define the flux parameters
flux_parameters = (ν=0.05, tri=tri,integ, ϕ_vec)

prob = FVMProblem(mesh, BCs;
    initial_condition,
    flux_function,
    flux_parameters,
    final_time=T)

# solve the system of ODEs with a condition ensuring non-negativity
sol = solve(prob, Tsit5(), isoutofdomain = (f,p,t) -> any(x->x<0, f), saveat=0.001)

# plit the solution as an animation
u = Observable(sol.u[1])
fig, ax, sc = tricontourf(tri, u, levels=0:5:50, colormap=:matter)
tightlimits!(ax)
record(fig, "animation", eachindex(sol)) do i
    u[] = sol.u[i]
end

```

I would appreciate your help with this a lot! 🙂

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [December 31, 2024, 2:42pm UTC](https://discourse.julialang.org/t/compute-integral-for-delaunaytriangulation-jl-to-be-able-to-solve-pde-with-finitevolumemethod-jl/124310/2 "2024-12-31T14:42:51Z")

</div>

Can you just split the integral into a sum over each triangle and then loop over each triangle and then sum accordingly? i.e. use `each_solid_triangle` and approximate the integral over each

---

<div class="post-metadata">

**Author:** ![B4rtel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/b4rtel/32/214171_2.png) [@B4rtel](https://discourse.julialang.org/u/B4rtel)\
**Post date:** [December 31, 2024, 4:47pm UTC](https://discourse.julialang.org/t/compute-integral-for-delaunaytriangulation-jl-to-be-able-to-solve-pde-with-finitevolumemethod-jl/124310/3 "2024-12-31T16:47:22Z")

</div>

Thx a lot for your answer! `each_solid_triangle` is definitely useful for my problem. 🙂 My idea was now to just evaluate the integral at the centroid of each triangle. I’m struggling to get the actual values though from the triangles I’m looping over since their points are related to the number of approximation points and not to the actual point the in the plane where they correspond to (so one entry of each\_solid\_triangle is smth like `(4,20,29) and not a smth like `([0.5,0.2],[-0.8,0.1],[0.3,-0.9])`). Do you know though how to work around that (so far I could only obtain all points of the triangluation via `get\_points(tri)` but not those of just one triangle)? Or maybe even how to efficiently integrate over the triangles?

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [December 31, 2024, 4:49pm UTC](https://discourse.julialang.org/t/compute-integral-for-delaunaytriangulation-jl-to-be-able-to-solve-pde-with-finitevolumemethod-jl/124310/4 "2024-12-31T16:49:30Z")

</div>

Usually you do something like

```julia
for T in each_solid_triangle(tri)
    p, q, r = get_point(tri, triangle_vertices(T)...)
end

```

then `p`, `q`, and `r` are the coordinates of the vertices of `T` in counter-clockwise order. Note that the centroid of each triangle can also be given by `DelaunayTriangulation.triangle_centroid(p, q, r)` once you have `T` (I thought there was a `triangle_centroid(::Triangulation, T)` method, but apparently not).

Here’s an example where I use this to approximate \int\_{[0, 1]^2} \cos(x - y)\,\mathrm dA:

```julia
julia> using DelaunayTriangulation

julia> f = (x, y) -> cos(x - y);

julia> tri = triangulate_rectangle(0, 1, 0, 1, 25, 25);

julia> int = 0.0
0.0

julia> for T in each_solid_triangle(tri)
           p, q, r = get_point(tri, triangle_vertices(T)...)
           c = DelaunayTriangulation.triangle_centroid(p, q, r)
           cx, cy = getxy(c)
           int += DelaunayTriangulation.triangle_area(p, q, r) * f(cx, cy)
       end

julia> int
0.9195284141899974

julia> 2 - 2cos(1)
0.9193953882637205

```

---

<div class="post-metadata">

**Author:** ![B4rtel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/b4rtel/32/214171_2.png) [@B4rtel](https://discourse.julialang.org/u/B4rtel)\
**Post date:** [January 2, 2025, 3:20pm UTC](https://discourse.julialang.org/t/compute-integral-for-delaunaytriangulation-jl-to-be-able-to-solve-pde-with-finitevolumemethod-jl/124310/5 "2025-01-02T15:20:13Z")

</div>

Thx a lot for the explanation! 🙂