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

Specific Domains Numerics
, , ,
1

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

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.

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! :slight_smile:

2

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

3

Thx a lot for your answer! each_solid_triangle is definitely useful for my problem. :slight_smile: 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?

4

Usually you do something like

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> 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
3 Likes
5

Thx a lot for the explanation! :slight_smile: