diff --git a/.all-contributorsrc b/.all-contributorsrc
deleted file mode 100644
index d4e9442..0000000
--- a/.all-contributorsrc
+++ /dev/null
@@ -1,5 +0,0 @@
-{
-  "projectName": "LevelSetMethods",
-  "projectOwner": "maltezfaria",
-  "files": ["README.md", "docs/src/index.md"]
-}
diff --git a/Project.toml b/Project.toml
index 2276e2a..583f646 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,5 +20,7 @@ MakieExt = "Makie"
 [compat]
 julia = "1.9"
 LinearAlgebra = "1"
+Makie = "0.21"
 MarchingCubes = "0.1.10"
 StaticArrays = "1"
+MMG_jll = "5"
diff --git a/README.md b/README.md
index 10206dc..20c8d7a 100644
--- a/README.md
+++ b/README.md
@@ -7,28 +7,11 @@
 [![Lint workflow Status](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Lint.yml/badge.svg?branch=main)](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Lint.yml?query=branch%3Amain)
 [![Docs workflow Status](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Docs.yml/badge.svg?branch=main)](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Docs.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/maltezfaria/LevelSetMethods.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/maltezfaria/LevelSetMethods.jl)
-[![DOI](https://zenodo.org/badge/DOI/FIXME)](https://doi.org/FIXME)
-[![Contributor Covenant](https://img.shields.io/badge/Contributor%20Covenant-2.1-4baaaa.svg)](CODE_OF_CONDUCT.md)
-[![All Contributors](https://img.shields.io/github/all-contributors/maltezfaria/LevelSetMethods.jl?labelColor=5e1ec7&color=c0ffee&style=flat-square)](#contributors)
 [![BestieTemplate](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/JuliaBesties/BestieTemplate.jl/main/docs/src/assets/badge.json)](https://github.com/JuliaBesties/BestieTemplate.jl)
 
-## How to Cite
+This package provides tools to solve level-set equations on cartesian meshes. It closely
+follows the book [Level Set Methods and Dynamic Implicit
+Surfaces](https://link.springer.com/book/10.1007/b98879) by Stanley Osher and Ronald Fedkiw.
 
-If you use LevelSetMethods.jl in your work, please cite using the reference given in [CITATION.cff](https://github.com/maltezfaria/LevelSetMethods.jl/blob/main/CITATION.cff).
-
-## Contributing
-
-If you want to make contributions of any kind, please first that a look into our [contributing guide directly on GitHub](docs/src/90-contributing.md) or the [contributing page on the website](https://maltezfaria.github.io/LevelSetMethods.jl/dev/90-contributing/)
-
----
-
-### Contributors
-
-<!-- ALL-CONTRIBUTORS-LIST:START - Do not remove or modify this section -->
-<!-- prettier-ignore-start -->
-<!-- markdownlint-disable -->
-
-<!-- markdownlint-restore -->
-<!-- prettier-ignore-end -->
-
-<!-- ALL-CONTRIBUTORS-LIST:END -->
+For more details, see the [stable](https://maltezfaria.github.io/LevelSetMethods.jl/stable)
+or [development](https://maltezfaria.github.io/LevelSetMethods.jl/dev) documentations.
diff --git a/docs/make.jl b/docs/make.jl
index 596fc45..79ca5cf 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -42,8 +42,8 @@ makedocs(;
         "terms.md",
         "time-integrators.md",
         "boundary-conditions.md",
-        hide("Extensions" => "extensions.md", ["extension-makie.md", "extension-mmg.md"]),
-        hide("Examples" => "examples.md", ["example-zalesak.md", "example-shape-optim.md"]),
+        "Extensions" => ["extension-makie.md", "extension-mmg.md"],
+        "Examples" => ["example-zalesak.md", "example-shape-optim.md"],
         numbered_pages,
     ),
     pagesonly = true, # ignore .md files not in the pages list
diff --git a/docs/src/90-contributing.md b/docs/src/90-contributing.md
deleted file mode 100644
index 204001c..0000000
--- a/docs/src/90-contributing.md
+++ /dev/null
@@ -1,25 +0,0 @@
-# [Contributing guidelines](@id contributing)
-
-First of all, thanks for the interest!
-
-We welcome all kinds of contribution, including, but not limited to code, documentation, examples, configuration, issue creating, etc.
-
-Be polite and respectful, and follow the code of conduct.
-
-## Bug reports and discussions
-
-If you think you found a bug, feel free to open an [issue](https://github.com/maltezfaria/LevelSetMethods.jl/issues).
-Focused suggestions and requests can also be opened as issues.
-Before opening a pull request, start an issue or a discussion on the topic, please.
-
-## Working on an issue
-
-If you found an issue that interests you, comment on that issue what your plans are.
-If the solution to the issue is clear, you can immediately create a pull request (see below).
-Otherwise, say what your proposed solution is and wait for a discussion around it.
-
-!!! tip
-    Feel free to ping us after a few days if there are no responses.
-
-If your solution involves code (or something that requires running the package locally), check the [developer documentation](91-developer.md).
-Otherwise, you can use the GitHub interface directly to create your pull request.
diff --git a/docs/src/91-developer.md b/docs/src/91-developer.md
deleted file mode 100644
index d9f30a4..0000000
--- a/docs/src/91-developer.md
+++ /dev/null
@@ -1,163 +0,0 @@
-# [Developer documentation](@id dev_docs)
-
-!!! note "Contributing guidelines"
-    If you haven't, please read the [Contributing guidelines](90-contributing.md) first.
-
-If you want to make contributions to this package that involves code, then this guide is for you.
-
-## First time clone
-
-!!! tip "If you have writing rights"
-    If you have writing rights, you don't have to fork. Instead, simply clone and skip ahead. Whenever **upstream** is mentioned, use **origin** instead.
-
-If this is the first time you work with this repository, follow the instructions below to clone the repository.
-
-1. Fork this repo
-2. Clone your repo (this will create a `git remote` called `origin`)
-3. Add this repo as a remote:
-
-   ```bash
-   git remote add upstream https://github.com/maltezfaria/LevelSetMethods.jl
-   ```
-
-This will ensure that you have two remotes in your git: `origin` and `upstream`.
-You will create branches and push to `origin`, and you will fetch and update your local `main` branch from `upstream`.
-
-## Linting and formatting
-
-Install a plugin on your editor to use [EditorConfig](https://editorconfig.org).
-This will ensure that your editor is configured with important formatting settings.
-
-We use [https://pre-commit.com](https://pre-commit.com) to run the linters and formatters.
-In particular, the Julia code is formatted using [JuliaFormatter.jl](https://github.com/domluna/JuliaFormatter.jl), so please install it globally first:
-
-```julia-repl
-julia> # Press ]
-pkg> activate
-pkg> add JuliaFormatter
-```
-
-To install `pre-commit`, we recommend using [pipx](https://pipx.pypa.io) as follows:
-
-```bash
-# Install pipx following the link
-pipx install pre-commit
-```
-
-With `pre-commit` installed, activate it as a pre-commit hook:
-
-```bash
-pre-commit install
-```
-
-To run the linting and formatting manually, enter the command below:
-
-```bash
-pre-commit run -a
-```
-
-**Now, you can only commit if all the pre-commit tests pass**.
-
-## Testing
-
-As with most Julia packages, you can just open Julia in the repository folder, activate the environment, and run `test`:
-
-```julia-repl
-julia> # press ]
-pkg> activate .
-pkg> test
-```
-
-## Working on a new issue
-
-We try to keep a linear history in this repo, so it is important to keep your branches up-to-date.
-
-1. Fetch from the remote and fast-forward your local main
-
-   ```bash
-   git fetch upstream
-   git switch main
-   git merge --ff-only upstream/main
-   ```
-
-2. Branch from `main` to address the issue (see below for naming)
-
-   ```bash
-   git switch -c 42-add-answer-universe
-   ```
-
-3. Push the new local branch to your personal remote repository
-
-   ```bash
-   git push -u origin 42-add-answer-universe
-   ```
-
-4. Create a pull request to merge your remote branch into the org main.
-
-### Branch naming
-
-- If there is an associated issue, add the issue number.
-- If there is no associated issue, **and the changes are small**, add a prefix such as "typo", "hotfix", "small-refactor", according to the type of update.
-- If the changes are not small and there is no associated issue, then create the issue first, so we can properly discuss the changes.
-- Use dash separated imperative wording related to the issue (e.g., `14-add-tests`, `15-fix-model`, `16-remove-obsolete-files`).
-
-### Commit message
-
-- Use imperative or present tense, for instance: *Add feature* or *Fix bug*.
-- Have informative titles.
-- When necessary, add a body with details.
-- If there are breaking changes, add the information to the commit message.
-
-### Before creating a pull request
-
-!!! tip "Atomic git commits"
-    Try to create "atomic git commits" (recommended reading: [The Utopic Git History](https://blog.esciencecenter.nl/the-utopic-git-history-d44b81c09593)).
-
-- Make sure the tests pass.
-- Make sure the pre-commit tests pass.
-- Fetch any `main` updates from upstream and rebase your branch, if necessary:
-
-  ```bash
-  git fetch upstream
-  git rebase upstream/main BRANCH_NAME
-  ```
-
-- Then you can open a pull request and work with the reviewer to address any issues.
-
-## Building and viewing the documentation locally
-
-Following the latest suggestions, we recommend using `LiveServer` to build the documentation.
-Here is how you do it:
-
-1. Run `julia --project=docs` to open Julia in the environment of the docs.
-1. If this is the first time building the docs
-   1. Press `]` to enter `pkg` mode
-   1. Run `pkg> dev .` to use the development version of your package
-   1. Press backspace to leave `pkg` mode
-1. Run `julia> using LiveServer`
-1. Run `julia> servedocs()`
-
-## Making a new release
-
-To create a new release, you can follow these simple steps:
-
-- Create a branch `release-x.y.z`
-- Update `version` in `Project.toml`
-- Update the `CHANGELOG.md`:
-  - Rename the section "Unreleased" to "[x.y.z] - yyyy-mm-dd" (i.e., version under brackets, dash, and date in ISO format)
-  - Add a new section on top of it named "Unreleased"
-  - Add a new link in the bottom for version "x.y.z"
-  - Change the "[unreleased]" link to use the latest version - end of line, `vx.y.z ... HEAD`.
-- Create a commit "Release vx.y.z", push, create a PR, wait for it to pass, merge the PR.
-- Go back to main screen and click on the latest commit (link: <https://github.com/maltezfaria/LevelSetMethods.jl/commit/main>)
-- At the bottom, write `@JuliaRegistrator register`
-
-After that, you only need to wait and verify:
-
-- Wait for the bot to comment (should take < 1m) with a link to a RP to the registry
-- Follow the link and wait for a comment on the auto-merge
-- The comment should said all is well and auto-merge should occur shortly
-- After the merge happens, TagBot will trigger and create a new GitHub tag. Check on <https://github.com/maltezfaria/LevelSetMethods.jl/releases>
-- After the release is create, a "docs" GitHub action will start for the tag.
-- After it passes, a deploy action will run.
-- After that runs, the [stable docs](https://maltezfaria.github.io/LevelSetMethods.jl/stable) should be updated. Check them and look for the version number.
diff --git a/docs/src/example-shape-optim.md b/docs/src/example-shape-optim.md
index fce8d4c..1885e25 100644
--- a/docs/src/example-shape-optim.md
+++ b/docs/src/example-shape-optim.md
@@ -4,14 +4,19 @@
 
 We consider in this example the *isoperimetric inequality* which states that among all closed surfaces enclosing a fixed area with volume $V_0 > 0$, the sphere is the one with minimal perimeter.
 We show here how to demonstrate this result through numerical optimization.
-To do this, we first define the problem mathematically (1):
+
+!!! warning
+    This example is purely illustrative. The optimization method used here has not been extensively tested.
+    Coupling the [LevelSetMethods](https://github.com/maltezfaria/LevelSetMethods.jl) toolbox to any simulation package makes it possible to solve PDE-constrained optimization problems (see for instance [allaire2004structural](@cite)).
+
+To do this, we first define the problem mathematically:
 
 ```math
     \begin{array}{rl}
-        \displaystyle\min_{\Omega \subset \mathbb{R}^N} & P(\Omega)
+        \displaystyle\min_{\Omega \subset \mathbb{R}^d} & P(\Omega)
         \\
         \text{u.c.} & V(\Omega) = V_0
-    \end{array},
+    \end{array},\qquad\text{(1)}
 ```
 
 where $P(\Omega), V(\Omega)$ are the perimeter and volume of $\Omega$ defined by
@@ -23,16 +28,17 @@ where $P(\Omega), V(\Omega)$ are the perimeter and volume of $\Omega$ defined by
     .
 ```
 
-The optimization problem (1) can be solved using the augmented Lagrangian approach by minimizing iteratively the following functional (2):
+The optimization problem $\text{(1)}$ can be solved using the augmented Lagrangian approach by minimizing iteratively the following functional:
 
 ```math
     f(\Omega) = P(\Omega) + \lambda (V(\Omega) - V_0) + \frac{\mu}{2} (V(\Omega) - V_0)^2
+    \qquad\text{(2)}
 ```
 
 where $\mu$ is a parameter updated during the course of the optimization.
-To minimize (2), we will use a gradient-based algorithm.
+To minimize $\text{(2)}$, we use a gradient-based algorithm.
 For this, we need to define what a *small variation* of $\Omega$ is.
-As such, we define for any shape $\Omega \subset \mathbb{R}^N$ its deformed configuration $\Omega_{\boldsymbol{\theta}}$ by a small vector field $\boldsymbol{\theta} \in W^{1,\infty}(\mathbb{R}^N, \mathbb{R}^N)$ as:
+As such, for any shape $\Omega \subset \mathbb{R}^d$ we define (following Hadamard method) its deformation $\Omega_{\boldsymbol{\theta}}$ by a small vector field $\boldsymbol{\theta} \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ as:
 
 ```math
     \Omega_{\boldsymbol{\theta}}
@@ -60,7 +66,7 @@ In other words, using $\boldsymbol{\theta} = - (\kappa + (\lambda + \mu (V(\Omeg
 
 ## Numerical solution using the level-set method
 
-If $\Omega$ is given by the level-set function $\phi_0 : \R^N \to \R$ then one associated with $\Omega_{\tau\boldsymbol{\theta}}$ is given by $\phi(\cdot, \tau)$ solution of
+If $\Omega$ is given by the level-set function $\phi_0 : \R^d \to \R$ then one associated with $\Omega_{\tau\boldsymbol{\theta}}$ is given by $\phi(\cdot, \tau)$ solution of
 
 ```math
     \partial_t \phi - \kappa |\nabla \phi| - (\lambda + \mu (V(\Omega) - V_0)) |\nabla \phi| = 0
@@ -127,3 +133,6 @@ end
 ```
 
 ![Optimization](optimization.gif)
+
+Different values and updates for the $\lambda$ and $\mu$ coefficients can be used to control the extent to which the optimization focuses on minimizing the objective or satisfying the constraint.
+Taking a smaller time step can also limit the oscillations observed at the end of optimization.
diff --git a/docs/src/example-zalesak.md b/docs/src/example-zalesak.md
index 1f60b31..c20b17d 100644
--- a/docs/src/example-zalesak.md
+++ b/docs/src/example-zalesak.md
@@ -1 +1,119 @@
 # [Zalesak disk](@id zalesak)
+
+This example demonstrates the setup and simulation of a Zalesak disk, a standard benchmark
+in computational fluid dynamics, using the `LevelSetMethods.jl`.
+
+## Setting up the Grid and Disk
+
+```@setup zalesak_disk_example
+using LevelSetMethods
+using GLMakie
+LevelSetMethods.set_makie_theme!()
+```
+
+First, we create a Cartesian grid to represent the computational domain. We define a
+circular disk with a rectangular cut to form the Zalesak disk.
+
+```@example zalesak_disk_example
+# Define the computational grid
+grid = CartesianGrid((-1.5, -1.5), (1.5, 1.5), (100, 100))
+# Define the center and radius of the circular disk
+center = (-0.75, 0)
+radius = 0.5
+# Define the height and width of the rectangular notch
+h = 1.0
+w = 0.2
+# Create the circular disk and rectangular notch
+disk = LevelSetMethods.circle(grid; center, radius)
+rec = LevelSetMethods.rectangle(grid; center = center .- (0, radius), width = (w, h))
+# Use set difference to carve out the notch in the disk
+ϕ = setdiff(disk, rec)
+plot(ϕ)
+current_figure() # hide
+```
+
+## Setting Up the Level Set Equation
+
+We use the level set method to evolve the disk over time with a velocity field that
+simulates rotation:
+
+```@example zalesak_disk_example
+# Define the advection equation with a rotational velocity field
+eq = LevelSetEquation(;
+    levelset = ϕ,
+    terms = AdvectionTerm((x, t) -> (-x[2], x[1])),
+    bc = NeumannBC(),
+)
+```
+
+## Evolving the Zalesak Disk
+
+To evolve the Zalesak disk, we integrate the level set equation over time. We will visualize
+the evolution using `GLMakie`:
+
+```@example zalesak_disk_example
+obs = Observable(eq)
+fig = Figure()
+ax = Axis(fig[1, 1])
+plot!(ax, obs)
+framerate = 30
+t0 = current_time(eq)
+tf = 2*π
+timestamps = range(t0, tf; step = 1 / framerate)
+record(fig, joinpath(@__DIR__, "zalesak2d.gif"), timestamps) do t_
+    integrate!(eq, t_)
+    return obs[] = eq
+end
+```
+
+![Zalesak 2D](zalesak2d.gif)
+
+We see some small smearing of the disk due to numerical diffusion; this is a common issue,
+and the situation would be much worse with a low-order upwind scheme.
+
+## Three-dimensional Zalesak Disk
+
+The same example can be run in 3D, but the solution takes longer to compute and visualize.
+
+!!! warning "Performance Warning"
+    The 3D example below may take a minute or two to run.
+
+```@setup zalesak_disk_3d
+using LevelSetMethods
+using GLMakie
+LevelSetMethods.set_makie_theme!()
+```
+
+```julia
+# 3D Zalesak's sphere example
+grid = CartesianGrid((-1, -1, -1), (1, 1, 1), (50, 50, 50))
+center = (-1 / 3, 0, 0)
+radius = 0.5
+disk = LevelSetMethods.sphere(grid; center, radius)
+rec = LevelSetMethods.rectangle(
+    grid;
+    center = center .+ (0, radius, 0),
+    width = (1 / 3, 1.0, 2),
+)
+ϕ = setdiff(disk, rec)
+eq = LevelSetEquation(;
+    levelset = ϕ,
+    terms = AdvectionTerm((x, t) -> π * SVector(x[2], -x[1], 0)),
+    bc = NeumannBC(),
+)
+LevelSetMethods.set_makie_theme!()
+eq.t = 0
+obs = Observable(eq)
+fig = Figure()
+ax = Axis3(fig[1, 1])
+plot!(ax, obs)
+framerate = 30
+tf = 2
+timestamps = range(0, tf, tf * framerate)
+record(fig, joinpath(@__DIR__, "zalesak3d.gif"), timestamps) do t_
+    integrate!(eq, t_)
+    return obs[] = eq
+end
+```
+
+![Zalesak 3D](zalesak3d.gif)
diff --git a/docs/src/examples.md b/docs/src/examples.md
deleted file mode 100644
index 839372a..0000000
--- a/docs/src/examples.md
+++ /dev/null
@@ -1,6 +0,0 @@
-# [Examples](@id examples)
-
-The following examples demonstrate some applications of the package.
-
-- [Zalesak disk](@ref zalesak)
-- [Shape optimization](@ref example-shape-optim)
diff --git a/docs/src/extension-mmg.md b/docs/src/extension-mmg.md
index d66d821..6733319 100644
--- a/docs/src/extension-mmg.md
+++ b/docs/src/extension-mmg.md
@@ -29,6 +29,7 @@ GLMakie.closeall() # hide
 using Gmsh
 try
   gmsh.initialize()
+  gmsh.option.setNumber("General.Verbosity", 0)
   gmsh.open(volume_mesh2d)
   gmsh.fltk.initialize()
   gmsh.write(joinpath(@__DIR__, "volume2d.png"))
@@ -74,6 +75,7 @@ GLMakie.closeall() # hide
 using Gmsh
 try
   gmsh.initialize()
+  gmsh.option.setNumber("General.Verbosity", 0)
   gmsh.open(surf_mesh3d)
   gmsh.fltk.initialize()
   gmsh.write(joinpath(@__DIR__, "surface3d.png"))
diff --git a/docs/src/extensions.md b/docs/src/extensions.md
deleted file mode 100644
index f0ca6e7..0000000
--- a/docs/src/extensions.md
+++ /dev/null
@@ -1,6 +0,0 @@
-# [Extensions](@id extensions)
-
-The following extensions are present in the LevelSetMethods.jl package.
-
-- [Makie extension](@ref extension-makie)
-- [Mmg extension](@ref extension-mmg)
diff --git a/docs/src/index.md b/docs/src/index.md
index aca5004..1058a66 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -4,19 +4,19 @@ CurrentModule = LevelSetMethods
 
 # LevelSetMethods
 
-Documentation for [LevelSetMethods](https://github.com/maltezfaria/LevelSetMethods.jl).
-
 ## Installation
 
-LevelSetMethods.jl is not yet registered in the Julia package registry. To install it, run
-the following command on a Julia REPL:
+Run the following command on a Julia REPL:
 
 ```julia
-using Pkg; Pkg.add("https://github.com/maltezfaria/LevelSetMethods.jl")
+using Pkg; Pkg.add("LevelSetMethods")
 ```
 
 This will install the latest tagged version of the package and its dependencies.
 
+For **visualization**, you may also want to install one of the
+[Makie](https://docs.makie.org) backends (we suggest `GLMakie` for 3D plots and animations).
+
 ## Overview
 
 This package defines a [`LevelSetEquation`](@ref) type that can be used to solve partial
@@ -28,20 +28,20 @@ differential equations of the form
 
 where
 
-- ``\phi : \mathbb{R}^d \to \mathbb{R}`` is the level set function
-- ``\boldsymbol{u} :\mathbb{R}^d \to \mathbb{R}^d`` is a given (external) velocity field
-- ``v : \mathbb{R}^d \to \mathbb{R}`` is a normal speed
-- ``b : \mathbb{R}^d \to \mathbb{R}`` is a function that multiplies the curvature ``\kappa =
+- ``\phi : \mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}`` is the level set function
+- ``\boldsymbol{u} :\mathbb{R}^d \times \mathbb{R}^+  \to \mathbb{R}^d`` is a given (external) velocity field
+- ``v : \mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}`` is a normal speed
+- ``b : \mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}`` is a function that multiplies the curvature ``\kappa =
   \nabla \cdot (\nabla \phi / |\nabla \phi|)``
 
 Here is how it looks in practice to create a simple `LevelSetEquation`:
 
 ```@example ls-intro
-using LevelSetMethods, StaticArrays
+using LevelSetMethods
 grid = CartesianGrid((-1, -1), (1, 1), (100, 100))
-# ϕ    = LevelSet(x -> sqrt(2*x[1]^2 + x[2]^2) - 1/2, grid)
-ϕ    = LevelSetMethods.dumbbell(grid)
-𝐮    = MeshField(x -> SVector(-x[2], x[1]), grid)
+# ϕ    = LevelSet(x -> sqrt(2*x[1]^2 + x[2]^2) - 1/2, grid) # a disk
+ϕ    = LevelSetMethods.dumbbell(grid) # a predefined shape
+𝐮    = (x,t) -> (-x[2], x[1])
 eq   = LevelSetEquation(;
   terms = (AdvectionTerm(𝐮),),
   levelset = ϕ,
@@ -55,7 +55,7 @@ from [Makie](https://docs.makie.org):
 ```@example ls-intro
 using GLMakie # loads the MakieExt from LevelSetMethods
 LevelSetMethods.set_makie_theme!() # optional theme customization
-plot(eq)
+plot(ϕ)
 ```
 
 To step it in time, we can use the [`integrate!`](@ref) function:
@@ -67,7 +67,7 @@ integrate!(eq, 1)
 This will advance the solution up to `t = 1`, modifying `ϕ` in the process:
 
 ```@example ls-intro
-plot(eq)
+plot(ϕ)
 ```
 
 Creating an animation can be achieved by calling `integrate!` in a loop and saving the
@@ -125,7 +125,7 @@ To learn more about the package, you should also check out the following section
 - The section on [boundary conditions](@ref boundary-conditions) for a description of the
   available boundary conditions and how to use them
 
-Finally, the [examples](@ref examples) section contains a list of examples that demonstrate some
+Finally, the examples section contains a list of examples that demonstrate some
 hopefully cool applications.
 
 ## Bibliography
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
index af389e0..4f2f374 100644
--- a/docs/src/refs.bib
+++ b/docs/src/refs.bib
@@ -28,3 +28,14 @@ @article{mitchell2007toolbox
   pages   = {6},
   year    = {2007}
 }
+
+@article{allaire2004structural,
+  title={Structural optimization using sensitivity analysis and a level-set method},
+  author={Allaire, Gr{\'e}goire and Jouve, Fran{\c{c}}ois and Toader, Anca-Maria},
+  journal={Journal of computational physics},
+  volume={194},
+  number={1},
+  pages={363--393},
+  year={2004},
+  publisher={Elsevier}
+}
diff --git a/docs/src/terms.md b/docs/src/terms.md
index 8995c7c..d64958d 100644
--- a/docs/src/terms.md
+++ b/docs/src/terms.md
@@ -153,7 +153,7 @@ is proportional to the mean curvature:
   b \kappa |\nabla \phi|
 ```
 
-where ``\kappa = \nabla \cdot (\nabla / |\nabla|)`` is the mean curvature. Note that the
+where ``\kappa = \nabla \cdot (\nabla \phi / |\nabla \phi|)`` is the mean curvature. Note that the
 coefficient ``b`` should be negative; a positive value of ``b`` would yield an ill-posed
 evolution problem (akin to a negative diffusion coefficient).
 
@@ -216,9 +216,9 @@ sdf = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) # signed distance functio
 LevelSetMethods.set_makie_theme!()
 fig = Figure(; size = (800, 400))
 ax = Axis(fig[1,1], title = "Signed distance function")
-contour!(ax, sdf; levels = [0, 0.5], labels = true, labelsize = 14)
+contour!(ax, sdf; levels = [0.25, 0, 0.5], labels = true, labelsize = 14)
 ax = Axis(fig[1,2], title = "ϕ at t = 0")
-contour!(ax, ϕ, levels = [0, 0.5], labels = true, labelsize = 14)
+contour!(ax, ϕ, levels = [0.25, 0, 0.5], labels = true, labelsize = 14)
 fig
 ```
 
@@ -230,7 +230,9 @@ fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.25, 0.5, 0.75])
     integrate!(eq, t)
     ax = Axis(fig[1,n], title = "t = $t")
-    contour!(ax, LevelSetMethods.current_state(eq); levels = [0, 0.5], labels = true, labelsize = 14)
+    contour!(ax, LevelSetMethods.current_state(eq); levels = [0.25, 0, 0.5], labels = true, labelsize = 14)
 end
 fig
 ```
+
+Note that `ϕ` converges to the signed distance function `sdf` shown in the first figure.
diff --git a/docs/src/zalesak3d.gif b/docs/src/zalesak3d.gif
new file mode 100644
index 0000000..7a3d757
Binary files /dev/null and b/docs/src/zalesak3d.gif differ
diff --git a/examples/advection_curvature_example.jl b/examples/advection_curvature_example.jl
deleted file mode 100644
index ee1facc..0000000
--- a/examples/advection_curvature_example.jl
+++ /dev/null
@@ -1,36 +0,0 @@
-using Test
-using LevelSetMethods
-using LinearAlgebra
-using Plots
-
-nx, ny = 100, 100
-x      = LinRange(-1, 1, nx)
-y      = LinRange(-1, 1, ny)
-hx, hy = step(x), step(y)
-grid   = CartesianGrid(x, y)
-ϕ      = LevelSet(grid) do (x, y)
-    return 0.5 - (4 * x)^2 - y^2
-end
-𝐮      = MeshField(grid) do (x, y)
-    return SVector(1, 0)
-end
-b      = MeshField(x -> -0.5, grid)
-term1  = AdvectionTerm(; velocity = 𝐮)
-term2  = CurvatureTerm(b)
-terms  = (term1, term2)
-bc     = PeriodicBC(1)
-buffer = zero(ϕ)
-dt     = 0.5 * (min(hx, hy))^2 # stiff
-t      = 0
-pgap   = 10
-
-integrator = ForwardEuler(0.5)
-
-anim = @animate for n in 0:200
-    fill!(values(buffer), 0)
-    _, t = LevelSetMethods.evolve!(buffer, integrator, ϕ, terms, bc, t)
-    if n % pgap == 0
-        plot(ϕ; title = "t=$t")
-    end
-end
-gif(anim, "test.gif"; fps = 15)
diff --git a/examples/advection_example.jl b/examples/advection_example.jl
deleted file mode 100644
index 588fabd..0000000
--- a/examples/advection_example.jl
+++ /dev/null
@@ -1,40 +0,0 @@
-using LevelSetMethods
-using LinearAlgebra
-using GLMakie
-
-nx, ny = 100, 100
-x = range(-1, 1, nx)
-y = range(-1, 1, ny)
-hx, hy = step(x), step(y)
-grid = CartesianGrid(x, y)
-bc = PeriodicBC(3)
-ϕ = LevelSet(grid, bc) do (x, y)
-    return 0.5^2 - x^2 - y^2
-end
-𝐮 = MeshField(grid) do (x, y)
-    return SVector(1, 0)
-end
-term1 = AdvectionTerm(; velocity = 𝐮)
-terms = (term1,)
-b = zero(ϕ)
-integrator = ForwardEuler(0.5)
-eq = LevelSetEquation(; terms, integrator, state = ϕ, t = 0, buffer = b)
-
-##
-# Load a default theme
-theme = LevelSetMethods.makie_theme()
-
-anim = with_theme(theme) do
-    eq.t = 0
-    obs = Observable(eq)
-    fig = Figure()
-    ax = Axis(fig[1, 1])
-    plot!(ax, obs)
-    framerate = 30
-    tf = 2
-    timestamps = range(0, tf, tf * framerate)
-    record(fig, "advection.gif", timestamps) do t_
-        integrate!(eq, t_)
-        return obs[] = eq
-    end
-end
diff --git a/examples/curvature_example.jl b/examples/curvature_example.jl
deleted file mode 100644
index 6754347..0000000
--- a/examples/curvature_example.jl
+++ /dev/null
@@ -1,27 +0,0 @@
-using Test
-using LevelSetMethods
-using LinearAlgebra
-using Plots
-
-nx, ny = 100, 100
-x = LinRange(-2, 2, nx)
-y = LinRange(-2, 2, ny)
-hx, hy = step(x), step(y)
-grid = CartesianGrid(x, y)
-bc = PeriodicBC(1)
-ϕ = LevelSet(grid, bc) do (x, y)
-    return 0.5 - (4x)^2 - y^2
-end
-b = MeshField(x -> -0.1, grid)
-term = CurvatureTerm(b)
-terms = (term,)
-b = zero(ϕ)
-integrator = ForwardEuler(0.5)
-eq = LevelSetEquation(; terms, integrator, state = ϕ, t = 0, buffer = b)
-
-anim = @animate for n in 0:50
-    tf = 0.01 * n
-    integrate!(eq, tf)
-    plot(eq; linecolor = :black)
-end
-gif(anim, "test.gif")
diff --git a/examples/mixed_terms_example.jl b/examples/mixed_terms_example.jl
deleted file mode 100644
index b94336d..0000000
--- a/examples/mixed_terms_example.jl
+++ /dev/null
@@ -1,34 +0,0 @@
-using Test
-using LevelSetMethods
-using LinearAlgebra
-using GLMakie
-
-grid = CartesianGrid((-1, -1), (1, 1), (100, 100))
-ϕ    = LevelSet(grid) do (x, y)
-    return 1.0
-end
-add_circle!(ϕ, SVector(0.5, 0.0), 0.25)
-add_circle!(ϕ, SVector(-0.5, 0.0), 0.25)
-add_rectangle!(ϕ, SVector(0.0, 0.0), SVector(1.0, 0.1))
-plot(ϕ; levels = [0])
-
-# create terms
-𝐮 = MeshField(x -> SVector(-x[2], x[1]), grid)
-eq = LevelSetEquation(; terms = AdvectionTerm(𝐮), levelset = ϕ, t = 0, bc = PeriodicBC())
-
-theme = LevelSetMethods.makie_theme()
-
-anim = with_theme(theme) do
-    eq.t = 0
-    obs = Observable(eq)
-    fig = Figure()
-    ax = Axis(fig[1, 1])
-    plot!(ax, obs)
-    framerate = 30
-    tf = 2π
-    timestamps = range(0, tf; step = 1 / framerate)
-    record(fig, joinpath(@__DIR__, "ls_intro.gif"), timestamps) do t_
-        integrate!(eq, t_)
-        return obs[] = eq
-    end
-end
diff --git a/examples/normal_advection_example.jl b/examples/normal_advection_example.jl
deleted file mode 100644
index 878dd13..0000000
--- a/examples/normal_advection_example.jl
+++ /dev/null
@@ -1,34 +0,0 @@
-using Test
-using LevelSetMethods
-using LinearAlgebra
-using Plots
-
-nx, ny = 100, 100
-x      = LinRange(-1, 1, nx)
-y      = LinRange(-1, 1, ny)
-hx, hy = step(x), step(y)
-grid   = CartesianGrid(x, y)
-bc     = PeriodicBC(2)
-ϕ      = LevelSet(grid, bc) do (x, y)
-    return 1.0
-end
-add_circle!(ϕ, SVector(0.5, 0.0), 0.25)
-add_circle!(ϕ, SVector(-0.5, 0.0), 0.25)
-add_rectangle!(ϕ, SVector(0.0, 0.0), SVector(1.0, 0.1))
-plot(ϕ)
-v = MeshField(grid) do (x, y)
-    return -1.0
-end
-term1 = NormalMotionTerm(v)
-terms = (term1,)
-b = zero(ϕ)
-integrator = ForwardEuler(0.5)
-eq = LevelSetEquation(; terms, integrator, state = ϕ, t = 0, buffer = b)
-
-dt = 0.01
-anim = @animate for n in 0:20
-    tf = dt * n
-    integrate!(eq, tf)
-    plot(eq; linecolor = :black)
-end
-gif(anim, "test.gif")
diff --git a/examples/plot_example.jl b/examples/plot_example.jl
deleted file mode 100644
index d7227fb..0000000
--- a/examples/plot_example.jl
+++ /dev/null
@@ -1,49 +0,0 @@
-using LevelSetMethods
-using LinearAlgebra
-using GLMakie
-
-## 2D
-nx, ny = 100, 100
-x = range(-1, 1, nx)
-y = range(-1, 1, ny)
-grid = CartesianGrid(x, y)
-bc = PeriodicBC(1)
-ϕ = LevelSet(grid, bc) do (x, y)
-    return 0.5^2 - x^2 - y^2
-end
-
-plot(ϕ; levels = [0])
-
-## 3D
-nx, ny, nz = 100, 100, 100
-x = range(-1, 1, nx)
-y = range(-1, 1, ny)
-z = range(-1, 1, nz)
-grid = CartesianGrid(x, y, z)
-bc = PeriodicBC(1)
-ϕ = LevelSet(grid, bc) do (x, y, z)
-    return 0.5^2 - x^2 - y^2 - z^2
-end
-
-plot(ϕ; algorithm = :iso, isovalue = 0, alpha = 0.5)
-
-## Following a solution over time
-obs = Observable(ϕ)
-
-fig = Figure()
-ax1 = Axis3(fig[1, 1])
-ax2 = Axis3(fig[1, 2])
-ax3 = Axis3(fig[2, 1:2])
-
-plot!(ax1, obs; algorithm = :iso, isovalue = 0, alpha = 0.5)
-plot!(ax2, obs; algorithm = :iso, isovalue = 0, alpha = 0.3)
-plot!(ax3, obs; algorithm = :iso, isovalue = 0, alpha = 1.0)
-fig
-
-for t in 0:0.01:1
-    sleep(0.1)
-    ϕ = LevelSet(grid, bc) do (x, y, z)
-        return 0.5^2 - (x)^2 - (y)^2 - (z - t)^2
-    end
-    obs[] = ϕ
-end
diff --git a/examples/reinitialization_example.jl b/examples/reinitialization_example.jl
deleted file mode 100644
index e8674d9..0000000
--- a/examples/reinitialization_example.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-using Test
-using LevelSetMethods
-using LinearAlgebra
-using Plots
-
-nx, ny = 100, 100
-x      = LinRange(-1, 1, nx)
-y      = LinRange(-1, 1, ny)
-hx, hy = step(x), step(y)
-grid   = CartesianGrid(x, y)
-bc     = PeriodicBC(2)
-ϕ      = LevelSet(grid, bc) do (x, y)
-    return 1.0
-end
-add_circle!(ϕ, SVector(0.0, 0.0), 0.9)
-remove_rectangle!(ϕ, SVector(0.0, 0.5), SVector(0.5, 1.0))
-# modify the level set function such that it is no longer
-# a signed distance function
-modifier = MeshField(grid) do (x, y)
-    return cos(x * 4) * sin(y * 4)
-end
-@. ϕ.vals = min(ϕ.vals, 0.5)^3.0 * (abs(modifier.vals) + 0.01)
-
-term1 = ReinitializationTerm()
-terms = (term1,)
-b = zero(ϕ)
-integrator = ForwardEuler(0.5)
-eq = LevelSetEquation(; terms, integrator, state = ϕ, t = 0, buffer = b)
-
-# comparison between the level set before and after reinitialization
-# heatmap(eq,title="Before")
-# plot!(eq,color=:black,colorbar=true)
-integrate!(eq, 1.0)
-heatmap(eq; title = "After")
-plot!(eq; color = :black, colorbar = true)
diff --git a/examples/rotation_example.jl b/examples/rotation_example.jl
deleted file mode 100644
index dfee446..0000000
--- a/examples/rotation_example.jl
+++ /dev/null
@@ -1,39 +0,0 @@
-using Test
-using LevelSetMethods
-using LinearAlgebra
-using Plots
-
-t1 = time()
-
-nx, ny = 201, 201
-x = LinRange(-1, 1, nx)
-y = LinRange(-1, 1, ny)
-hx, hy = step(x), step(y)
-grid = CartesianGrid(x, y)
-bc = PeriodicBC(1)
-ϕ = LevelSet(grid, bc) do (x, y)
-    return sqrt((x + 0.4)^2 + y^2) - 0.35
-end
-𝐮 = MeshField(grid) do (x, y)
-    return 2π * SVector(-y, x)
-end
-term1 = AdvectionTerm(; velocity = 𝐮)
-terms = (term1,)
-b = zero(ϕ)
-integrator = ForwardEuler(0.5)
-
-eq = LevelSetEquation(; terms, integrator, state = ϕ, t = 0, buffer = b)
-fig = plot(eq)
-@time integrate!(eq, 1)
-plot!(fig, eq)
-
-t2 = time()
-
-@info "Total time : $(t2 - t1)"
-
-# anim = @animate for n ∈ 0:100
-#     tf = 0.02*n
-#     integrate!(eq,tf)
-#     plot(eq)
-# end
-# gif(anim, "test.gif")
diff --git a/examples/time_integrators_example.jl b/examples/time_integrators_example.jl
deleted file mode 100644
index b059c93..0000000
--- a/examples/time_integrators_example.jl
+++ /dev/null
@@ -1,57 +0,0 @@
-using Test
-using LevelSetMethods
-using LinearAlgebra
-using Plots
-
-nx, ny = 100, 100
-x      = LinRange(-1, 1, nx)
-y      = LinRange(-1, 1, ny)
-hx, hy = step(x), step(y)
-grid   = CartesianGrid(x, y)
-bc     = PeriodicBC(3)
-ϕ      = LevelSet(grid, bc) do (x, y)
-    return 1.0
-end
-add_circle!(ϕ, SVector(0.5, 0.0), 0.25)
-add_circle!(ϕ, SVector(-0.5, 0.0), 0.25)
-add_rectangle!(ϕ, SVector(0.0, 0.0), SVector(1.0, 0.1))
-v = MeshField(grid) do (x, y)
-    return -0.1
-end
-𝐮 = MeshField(grid) do (x, y)
-    return SVector(-y, x)
-end
-b = MeshField(grid) do (x, y)
-    return -min(hx, hy)
-end
-term1 = NormalMotionTerm(v)
-term2 = AdvectionTerm(𝐮)
-term3 = CurvatureTerm(b)
-terms = (term1, term2, term3)
-b = (zero(ϕ), zero(ϕ))
-integrator = ForwardEuler(0.5)
-eq = LevelSetEquation(; terms, integrator, state = ϕ, t = 0, buffer = b[1])
-integrator = RK2(0.5)
-eq2 =
-    LevelSetEquation(; terms, integrator, state = deepcopy(ϕ), t = 0, buffer = deepcopy(b))
-integrator = LevelSetMethods.RKLM2(0.5)
-eq3 = LevelSetEquation(;
-    terms,
-    integrator,
-    state = deepcopy(ϕ),
-    t = 0,
-    buffer = deepcopy(b[1]),
-)
-
-dt = 0.01
-anim = @animate for n in 0:50
-    tf = dt * n
-    integrate!(eq, tf)
-    integrate!(eq2, tf)
-    integrate!(eq3, tf)
-    fig = plot(eq; linecolor = :black, linestyle = :dash)
-    plot!(fig, eq2; linecolor = :blue, linestyle = :dot)
-    plot!(fig, eq3; linecolor = :red, linestyle = :dashdot)
-    # pplot(eq,eq2)
-end
-gif(anim, "test.gif")
diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl
index d8379fb..4940a8c 100644
--- a/ext/MakieExt.jl
+++ b/ext/MakieExt.jl
@@ -74,6 +74,7 @@ function Makie.plot!(p::LevelSetPlot)
     ϕ = @lift LSM.current_state($eq)
     N = @lift LSM.dimension($ϕ)
     if to_value(N) == 2
+        contourf!(p, ϕ; levels = [0], extendlow = (:lightgray, 0.5))
         contour!(p, ϕ; levels = [0], linewidth = 2, color = :black)
     elseif to_value(N) == 3
         volume!(p, ϕ; algorithm = :iso, isovalue = 0, alpha = 0.5)
@@ -94,7 +95,6 @@ function Makie.preferred_axis_type(p::LevelSetPlot)
     if to_value(dim) == 2
         return Axis
     elseif to_value(dim) == 3
-        @info "here"
         return LScene
     else
         throw(ArgumentError("Plot of $dim dimensional level-set not supported."))
diff --git a/src/levelset.jl b/src/levelset.jl
index b877689..9bd6f50 100644
--- a/src/levelset.jl
+++ b/src/levelset.jl
@@ -289,7 +289,6 @@ function zalesak_disk(grid; center = (0, 0), radius = 0.5, width = 0.25, height
         throw(ArgumentError("zalesak disk shape is only available in two dimensions"))
     disk = circle(grid; center = center, radius = radius)
     rec = rectangle(grid; center = center .- (0, radius), width = (width, height))
-    # return union(disk, rec)
     return setdiff(disk, rec)
 end
 
diff --git a/test/test-mmgext.jl b/test/test-mmgext.jl
deleted file mode 100644
index bf33211..0000000
--- a/test/test-mmgext.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-using LevelSetMethods
-using MMG_jll
-
-# if we got here is that the MMG library was loaded
-@test true == true