# FiniteVolumeMethod.jl v1.0: Solving conservation PDEs in 2D, now with systems and internal conditions

**URL:** <https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715>\
**Category:** Package Announcements\
**Tags:** package, numerics, differentialequation\
**Created:** [September 10, 2023, 7:59am UTC](https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715 "2023-09-10T07:59:50Z")\
**Posts on this page:** 6\
**Page:** 1

<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:** [September 10, 2023, 7:59am UTC](https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715/1 "2023-09-10T07:59:50Z")

</div>

I have just released v1.0 of my [FiniteVolumeMethod.jl](https://github.com/DanielVandH/FiniteVolumeMethod.jl) package, which is a package for solving PDEs of the form

\frac{\partial u(\boldsymbol x, t)}{\partial t} + \boldsymbol\nabla\boldsymbol\cdot\boldsymbol q(\boldsymbol x, t, u) = S(\boldsymbol x, t, u),

with additional support provided for systems of PDEs of the above form, steady state problems where \partial u/\partial t = 0. Previous versions of FiniteVolumeMethod.jl only support scalar PDEs instead of systems. The package uses my other package [DelaunayTriangulation.jl](https://github.com/DanielVandH/DelaunayTriangulation.jl) for defining the mesh. Additionally, support is now added for inhomogeneous Neumann boundary conditions and for internal conditions. All the new changes are listed in the [release](https://github.com/DanielVandH/FiniteVolumeMethod.jl/releases/tag/v1.0.0). I’ve also rewrote the [mathematical details section](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/math/) to be (hopefully) a lot clearer.

There are 15 examples in documentation which cover a range of examples, covering e.g. domains with holes, mixed boundary conditions, interpolation with [NaturalNeighbours.jl](https://github.com/DanielVandH/NaturalNeighbours.jl). Here are two examples showing the support for internal conditions and for systems of PDEs. For an explanation, see the corresponding tutorials in the [docs](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/).

# Laplace’s Equation with Internal Conditions

As a first example, I’ll solve Laplace’s equation with an internal constraint so that the solution is zero on some line inside the domain. In particular:

\begin{equation\*} \begin{aligned} \boldsymbol\nabla^2 u &= 0 & \quad (x, y) \in [0, 1]^2, \\ u(0, y) &= 100y, & 0 \leq y \leq 1, \\ u(1, y) &= 100y, & 0 \leq y \leq 1, \\ u(x, 0) &= 0, & 0 \leq x \leq 1, \\ u(x, 1) &= 100, & 0 \leq x\leq 1, \\ u(1/2, y) &= 0, & 0 \leq y \leq 2/5. \end{aligned} \end{equation\*}

This last condition u(1/2, y)=0 for 0 \leq y \leq 2/5 is the interior constraint. Here is how we solve it. The main difference from previous versions is that `InternalConditions` is needed to enforce constraints at nodes in the interior.

```julia
# The mesh. Need to add points on the line to enforce the interior constraints.
using DelaunayTriangulation, FiniteVolumeMethod
tri = triangulate_rectangle(0, 1, 0, 1, 50, 50, single_boundary=false)
new_points = LinRange(0, 2 / 5, 250)
for y in new_points
    add_point!(tri, 1 / 2, y)
end
refine!(tri, max_area=1e-4)
mesh = FVMGeometry(tri)

# The boundary conditions 
bc_bot = (x, y, t, u, p) -> zero(u)
bc_right = (x, y, t, u, p) -> oftype(u, 100y) # helpful to have each bc return the same type
bc_top = (x, y, t, u, p) -> oftype(u, 100)
bc_left = (x, y, t, u, p) -> oftype(u, 100y)
bcs = (bc_bot, bc_right, bc_top, bc_left)
types = (Dirichlet, Dirichlet, Dirichlet, Dirichlet)
BCs = BoundaryConditions(mesh, bcs, types)

# The internal conditions 
function find_all_points_on_line(tri)
    vertices = Int[]
    for i in each_solid_vertex(tri)
        x, y = get_point(tri, i)
        if x == 1 / 2 && 0 ≤ y ≤ 2 / 5
            push!(vertices, i)
        end
    end
    return vertices
end
vertices = find_all_points_on_line(tri)
ICs = InternalConditions((x, y, t, u, p) -> zero(u),
    dirichlet_nodes=Dict(vertices .=> 1))

# The problem 
initial_condition = zeros(num_points(tri)) # initial guess for the steady state
for i in each_solid_vertex(tri)
    x, y = get_point(tri, i)
    initial_condition[i] = ifelse(x == 1 / 2 && 0 ≤ y ≤ 2 / 5, 0, 100y)
end
diffusion_function = (x, y, t, u, p) -> one(u) # ∇²u = ∇⋅[D∇u], D = 1
final_time = Inf
prob = FVMProblem(mesh, BCs, ICs;
    diffusion_function,
    initial_condition,
    final_time)
steady_prob = SteadyFVMProblem(prob)

# The solution 
using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve=KLUFactorization())))
using CairoMakie 
fig, ax, sc = tricontourf(tri, sol.u, levels=LinRange(0, 100, 28))
tightlimits!(ax)
fig

```

 ![Untitled](https://global.discourse-cdn.com/julialang/original/3X/a/8/a8b517acc93d5c03e9f2859f01b7e3744e6ad401.png)

# Keller-Segel Model of Chemotaxis

Now I’ll solve the following Keller-Segel model of chemotaxis (this example takes quite a long time to run, so I wouldn’t recommend you run it yourself - just has a nice animation). For systems, you need to provide an `FVMProblem` for each equation.

\begin{equation\*} \begin{aligned} \frac{\partial u}{\partial t} &= \boldsymbol\nabla^2 u - \boldsymbol\nabla\boldsymbol\cdot\left(\frac{cu}{1+u^2}\boldsymbol\nabla v\right) + u(1-u), \\ \frac{\partial v}{\partial t} &= D\boldsymbol\nabla^2 v + u - av, \\ \end{aligned} \end{equation\*}

where (x, y) \in [0, 100]^2 and we have homogeneous Neumann boundary conditions. Note that I’ve had to implement this problem by rewriting in the conservation form.

```math
# Defining the problem
using FiniteVolumeMethod, DelaunayTriangulation
tri = triangulate_rectangle(0, 100, 0, 100, 250, 250, single_boundary=true)
mesh = FVMGeometry(tri)
bc_u = (x, y, t, (u, v), p) -> zero(u)
bc_v = (x, y, t, (u, v), p) -> zero(v)
BCs_u = BoundaryConditions(mesh, bc_u, Neumann)
BCs_v = BoundaryConditions(mesh, bc_v, Neumann)
q_u = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    u = αu * x + βu * y + γu
    ∇u = (αu, βu)
    ∇v = (αv, βv)
    χu = p.c * u / (1 + u^2)
    _q = χu .* ∇v .- ∇u
    return _q
end
q_v = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    ∇v = (αv, βv)
    _q = -p.D .* ∇v
    return _q
end
S_u = (x, y, t, (u, v), p) -> begin
    return u * (1 - u)
end
S_v = (x, y, t, (u, v), p) -> begin
    return u - p.a * v
end
q_u_parameters = (c=4.0,)
q_v_parameters = (D=1.0,)
S_v_parameters = (a=0.1,)
u_initial_condition = 0.01rand(num_points(tri))
v_initial_condition = zeros(num_points(tri))
final_time = 1000.0
u_prob = FVMProblem(mesh, BCs_u;
    flux_function=q_u, flux_parameters=q_u_parameters,
    source_function=S_u,
    initial_condition=u_initial_condition, final_time=final_time)
v_prob = FVMProblem(mesh, BCs_v;
    flux_function=q_v, flux_parameters=q_v_parameters,
    source_function=S_v, source_parameters=S_v_parameters,
    initial_condition=v_initial_condition, final_time=final_time)
prob = FVMSystem(u_prob, v_prob)

# Solving and visualising
using OrdinaryDiffEq, LinearSolve, CairoMakie
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=1.0) # very slow
fig = Figure(fontsize=44)
x = LinRange(0, 100, 250)
y = LinRange(0, 100, 250)
i = Observable(1)
axu = Axis(fig[1, 1], width=600, height=600,
    title=map(i -> L"u(x,~ y,~ %$(sol.t[i]))", i), xlabel=L"x", ylabel=L"y")
axv = Axis(fig[1, 2], width=600, height=600,
    title=map(i -> L"v(x,~ y,~ %$(sol.t[i]))", i), xlabel=L"x", ylabel=L"y")
u = map(i -> reshape(sol.u[i][1, :], 250, 250), i)
v = map(i -> reshape(sol.u[i][2, :], 250, 250), i)
heatmap!(axu, x, y, u, colorrange=(0.0, 2.5), colormap=:turbo)
heatmap!(axv, x, y, v, colorrange=(0.0, 10.0), colormap=:turbo)
resize_to_layout!(fig)
record(fig, "keller_segel_chemotaxis.mp4", eachindex(sol);
    framerate=60) do _i
    i[] = _i
end

```

The complete list of examples that I give in the docs are below.

- [Diffusion Equation on a Square Plate](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Diffusion-Equation-on-a-Square-Plate)
- [Diffusion Equation in a Wedge with Mixed Boundary Conditions](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Diffusion-Equation-in-a-Wedge-with-Mixed-Boundary-Conditions)
- [Reaction Diffusion Equation with a Time-dependent Dirichlet Boundary Condition on a Disk](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Reaction-Diffusion-Equation-with-a-Time-dependent-Dirichlet-Boundary-Condition-on-a-Disk)
- [Porous-Medium Equation](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Porous-Medium-Equation)
- [Porous-Fisher Equation and Travelling Waves](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Porous-Fisher-Equation-and-Travelling-Waves)
- [Piecewise Linear and Natural Neighbour Interpolation for an Advection-Diffusion Equation](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Piecewise-Linear-and-Natural-Neighbour-Interpolation-for-an-Advection-Diffusion-Equation)
- [Helmholtz Equation with Inhomogeneous Boundary Conditions](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Helmholtz-Equation-with-Inhomogeneous-Boundary-Conditions)
- [Laplace’s Equation with Internal Dirichlet Conditions](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Laplace's-Equation-with-Internal-Dirichlet-Conditions)
- [Equilibrium Temperature Distribution with Mixed Boundary Conditions and using EnsembleProblems](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Equilibrium-Temperature-Distribution-with-Mixed-Boundary-Conditions-and-using-EnsembleProblems)
- [A Reaction-Diffusion Brusselator System of PDEs](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#A-Reaction-Diffusion-Brusselator-System-of-PDEs)
- [Gray-Scott Model: Turing Patterns from a Coupled Reaction-Diffusion System](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Gray-Scott-Model:-Turing-Patterns-from-a-Coupled-Reaction-Diffusion-System)
- [Diffusion Equation on an Annulus](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Diffusion-Equation-on-an-Annulus)
- [Mean Exit Time](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Mean-Exit-Time)
- [Solving Mazes with Laplace’s Equation](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Solving-Mazes-with-Laplace's-Equation)
- [Keller-Segel Model of Chemotaxis](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/tutorials/overview/#Keller-Segel-Model-of-Chemotaxis)

---

<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:** [September 10, 2023, 8:02am UTC](https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715/2 "2023-09-10T08:02:04Z")

</div>

Also just another fun example: [Solving mazes with Laplace’s equation](https://danielvandh.github.io/FiniteVolumeMethod.jl/stable/tutorials/solving_mazes_with_laplaces_equation/):

[https://global.discourse-cdn.com/julialang/original/3X/a/7/a7844f8fec6bc5632febf71d93139a5d5495b3ba.mp4](https://global.discourse-cdn.com/julialang/original/3X/a/7/a7844f8fec6bc5632febf71d93139a5d5495b3ba.mp4)

---

<div class="post-metadata">

**Author:** ![TheLateKronos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thelatekronos/32/12824_2.png) [@TheLateKronos](https://discourse.julialang.org/u/TheLateKronos)\
**Post date:** [September 10, 2023, 1:57pm UTC](https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715/3 "2023-09-10T13:57:27Z")

</div>

Great work 👏

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [September 13, 2023, 12:47pm UTC](https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715/4 "2023-09-13T12:47:48Z")

</div>

cool

---

<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:** [September 16, 2023, 4:50pm UTC](https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715/5 "2023-09-16T16:50:20Z")

</div>

I’ve just updated the package so that I also provide some [specialised solvers for specific problems](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/wyos/overview/), namely:

1. `DiffusionEquation`s: \partial\_tu = \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u].
2. `MeanExitTimeProblem`s: \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla T(\boldsymbol x)] = -1.
3. `LinearReactionDiffusionEquation`s: \partial\_tu = \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] + f(\boldsymbol x)u.
4. `PoissonsEquation`: \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] = f(\boldsymbol x).
5. `LaplacesEquation`: \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] = 0.

This list will hopefully include semilinear systems in the future.

I also derive, in the linked section, the details for all of these methods, so that it’s also clear if anyone ever wants to make their own specialised solver for their problem. These are significantly faster than if you use the more general `FVMProblem` approach, since the linearity in the diffusion term can be easily exploited. I also fixed a design issue that was accidentally causing FVMSystem problems to take 50x longer than they should… so now the Keller-Segel model problem above can be solved relatively quickly (if you also use `CVODE_BDF(linear_solver=:GMRES)` from Sundials.jl.

Here’s an example of one of these template problems (see [here](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/wyos/poissons_equation/)), solving Poisson’s equation describing a problem with two parallel capacitor plates with a space-varying dielectric function and charge density so that \boldsymbol\nabla\boldsymbol\cdot[\epsilon(\boldsymbol x)\boldsymbol\nabla V(\boldsymbol x)] = -\rho(\boldsymbol x)/\epsilon\_0. Here is the setup of the problem. It’s a bit long since \epsilon and \rho have complicated definitions and we need the indices of the two plates for the internal Dirichlet conditions; none of this is new in 1.1.

> **Setting up the problem**
>
> ```julia
> using DelaunayTriangulation, FiniteVolumeMethod, CairoMakie, LinearAlgebra, SpecialFunctions
> ## The mesh
> a, b, c, d = 0.0, 10.0, 0.0, 10.0
> e, f = (2.0, 3.0), (8.0, 3.0)
> g, h = (2.0, 7.0), (8.0, 7.0)
> points = [(a, c), (b, c), (b, d), (a, d), e, f, g, h]
> boundary_nodes = [[1, 2], [2, 3], [3, 4], [4, 1]]
> edges = Set(((5, 6), (7, 8))) # the two capacitor plates
> tri = triangulate(points; boundary_nodes, edges, delete_ghosts=false)
> refine!(tri; max_area=1e-4get_total_area(tri))
> mesh = FVMGeometry(tri)
> ## The boundary and internal conditions
> zero_f = (x, y, t, u, p) -> zero(x)
> bc_f = (zero_f, zero_f, zero_f, zero_f)
> bc_types = (Dirichlet, Neumann, Dirichlet, Neumann) # bot, right, top, left
> BCs = BoundaryConditions(mesh, bc_f, bc_types)
> function find_vertices_on_plates(tri)
> lower_plate = Int[]
> upper_plate = Int[]
> for i in each_solid_vertex(tri)
> x, y = get_point(tri, i)
> in_range = 2 ≤ x ≤ 8
> if in_range
> y == 3 && push!(lower_plate, i)
> y == 7 && push!(upper_plate, i)
> end
> end
> return lower_plate, upper_plate
> end
> lower_plate, upper_plate = find_vertices_on_plates(tri)
> dirichlet_nodes = Dict(
> (lower_plate .=> 1)...,
> (upper_plate .=> 2)...
> )
> internal_f = ((x, y, t, u, p) -> -one(x), (x, y, t, u, p) -> one(x))
> ICs = InternalConditions(internal_f, dirichlet_nodes=dirichlet_nodes)
> ## The dielectric function and charge density
> function dist_to_line(p, a, b)
> ℓ² = norm(a .- b)^2
> t = max(0.0, min(1.0, dot(p .- a, b .- a) / ℓ²))
> pv = a .+ t .* (b .- a)
> return norm(p .- pv)
> end
> function dist_to_plates(x, y)
> p = (x, y)
> a1, b1 = (2.0, 3.0), (8.0, 3.0)
> a2, b2 = (2.0, 7.0), (8.0, 7.0)
> d1 = dist_to_line(p, a1, b1)
> d2 = dist_to_line(p, a2, b2)
> return min(d1, d2)
> end
> function dielectric_function(x, y, p)
> r = dist_to_plates(x, y)
> return 1 + (p.ϵ₀ - 1) * (1 + erf(r / p.Δ)) / 2
> end
> function charge_density(x, y, p)
> r = dist_to_plates(x, y)
> return p.Q / (2π) * exp(-r^2 / 2)
> end
> function plate_source_function(x, y, p)
> ρ = charge_density(x, y, p)
> return -ρ / p.ϵ₀
> end
> diffusion_parameters = (ϵ₀=8.8541878128e-12, Δ=4.0)
> source_parameters = (ϵ₀=8.8541878128e-12, Q=1e-6)
> 
> ```

Here is the new part. Rather than using `FVMProblem`, you can use `PoissonsEquation`:

```julia
diffusion_parameters = (ϵ₀=8.8541878128e-12, Δ=4.0)
source_parameters = (ϵ₀=8.8541878128e-12, Q=1e-6)
prob = PoissonsEquation(mesh, BCs, ICs;
    diffusion_function=dielectric_function,
    diffusion_parameters=diffusion_parameters,
    source_function=plate_source_function,
    source_parameters=source_parameters)
sol = solve(prob, KLUFactorization()) # steady state problems are a linear solve

```

This takes about 9.68 ms, while the corresponding FVMProblem takes 352 ms - a great speedup! In the [docs](https://danielvandh.github.io/FiniteVolumeMethod.jl/dev/wyos/poissons_equation/), I also show how the voltage and electric fields can be visualised, using [NaturalNeighbours.jl](https://github.com/DanielVandH/NaturalNeighbours.jl) to compute the gradients, giving the figure below.

 ![poissons_equation-47](https://global.discourse-cdn.com/julialang/original/3X/8/a/8a1c7c85cc41c24072ae472b86e78b2ae3425973.png)

---

<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:** [May 11, 2024, 10:56am UTC](https://discourse.julialang.org/t/finitevolumemethod-jl-v1-0-solving-conservation-pdes-in-2d-now-with-systems-and-internal-conditions/103715/6 "2024-05-11T10:56:15Z")

</div>

Now that DelaunayTriangulation v1.0 has been released, support for curve-bounded domains now exists. The documentation of FiniteVolumeMethod.jl has now been updated, where applicable, to use these new features. This allows for us to specify, for example, a circle as a function rather than manually discretising it. For example, instead of specifying the circle as

```julia
using DelaunayTriangulation, CairoMakie
R = 2.0
θ = LinRange(0, 2π, 250)
x = @. R * cos(θ)
y = @. R * sin(θ)
x[end] = x[1]
y[end] = y[1]
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
tri = triangulate(points; boundary_nodes)
refine!(tri; max_area=1e-3 * π * R^2)

```

you can do

```julia
circle = CircularArc((0.0, R), (0.0, R), (0.0, 0.0))
points = NTuple{2,Float64}[] # still need to provide some points, even if empty
tri2 = triangulate(points; boundary_nodes = [circle])
refine!(tri2; max_area=1e-3 * π * R^2)

```

which can significantly improve the quality of the mesh and, thus, the solution of your PDE - just see how the quality of the mesh around the boundary changes here:

```julia
fig = Figure()
ax=Axis(fig[1,1],width=400,height=400,title="Manual discretisation")
triplot!(ax, tri, color=:black)
ax=Axis(fig[1,2],width=400,height=400,title="Automatic discretisation")
triplot!(ax, tri2, color=:black)
resize_to_layout!(fig)
fig

```

 ![Yf1T9uh](https://global.discourse-cdn.com/julialang/original/3X/3/5/3570e2c62c5b51c0c1e2550c64361bbc1bbc4e36.jpeg)

Please see the docs for more examples, and the docs of DelaunayTriangulation.jl for more examples of specifying curve-bounded domains.