## Code

### High-order algorithms for computing closest points on implicitly-defined surfaces

A powerful technique for representing curves in two dimensions (N = 2) and surfaces in three dimension (N = 3) is to define them implicitly as a level set of a continuous, N-dimensional scalar function. This technique of embedding the surface geometry in a higher-dimensional function (which is often called the "level set function") leads to a wide array of mathematical and computational advantages, as exemplified by the level set method for moving interface problems, and the Voronoi implicit interface method for computing multiphase interface evolution.

Numerical methods making use of this idea often require accurate approximations of minimum distances to and closest points on implicitly defined surfaces. I have developed high-order accurate algorithms for this purpose and these are described in the paper:

• R. I. Saye, High-order methods for computing distances to implicitly defined surfaces, Communications in Applied Mathematics and Computational Science, 9(1), 107–141 (2014). Link pdf

Provided on this page is C++ code implementing the algorithms developed in the paper. The code is freely available (subject to the license) and can be used to implement, for example, high-order accurate reinitialisation/redistancing algorithms in level set methods. The code mainly applies to the case that the level set function is defined on a rectangular Cartesian grid. However, as discussed in the paper, it is possible to extend the algorithms to the case of unstructured grids – one could use the code as a starting point. For example, the implementations of the k-d tree and Newton's method could be used as-is without modification in such an adaptation.

#### C++ code

Most of the code is in C++ header files (`.hpp`). Except for two demonstration programs, there is only one `.cpp` file (`poly.cpp`) that needs to be compiled. This file contains the precomputed pseudoinverses of the interpolation matrices and the various stencils for performing least-squares polynomial interpolation on a rectangular Cartesian grid.

The package contains:

• `kdtree.hpp` – Implements a k-d tree data structure for a given collection of points in RN. This particular k-d tree has been optimised for clouds of points that originate from smooth codimension-one surfaces using coordinate transformations that result in "tight" bounding boxes for nodes in the tree.
• `hocp.hpp` – Main driver routines for the high-order closest point algorithm applied to rectangular Cartesian grids.
• `newtoncp.hpp` – Newton's method for the constrained minimum-distance optimisation problem applied to finding closest points on the zero level set of a function (typically a polynomial).
• `poly.hpp` and `poly.cpp` – Implements the least-squares polynomial interpolation algorithms for 10 different classes of polynomial in 2D and 3D.
• `hocp_util.hpp` – Helper methods including `MultiLoop` which can be used to write N-dimensional nested for loops, where N is a template parameter, for iterating over grid points in a N-dimensional array.
• `testSimple.cpp` – A very simple demonstration of how to reinitialise a level set function in 2D and 3D using the high-order closest point algorithm.
• `testConvergence.cpp` – Measures convergence rates of the algorithm on several test interface geometries (circle, sphere, square, cube, ellipse, ellipsoid, rounded rectangle and rounded cylinder), where the test problem is specified as a command-line parameter.

#### Instructions and code dependencies

Much of the code is templated on the dimension N. To assist with this functionality, the code makes use of blitz++, an open-source implementation of N-dimensional arrays and fixed-length vectors in C++ with convenient expression template techniques.

To compile and test the code:

3. Compile:
1. If you are using blitz-0.9, uncomment the line "`#include <blitz/tinyvec-et.h>`" in `hocp_util.hpp`. If you are using blitz-0.10, you do not need to do this.
2. Modify the `Makefile` to correctly point to the location of the `blitz` directory in the include command.
3. `make`, or compile according to your favourite methods.
4. Run:
1. The program "`testSimple`" demonstrates a very simple problem of reinitialising a level set function in 2D and 3D. Type in a dimension (2 or 3) and the number of grid points in each axis to have the program reinitialise a spherical interface level set function and report the maximum-norm error of the computed signed distance function.
2. The program "`testConvergence`" implements a series of convergence tests on different interface geometries and can be used to replicate the results in the above mentioned paper. Use it by specifying the test problem via command-line parameters – for example, `./testConvergence 2 sphere` or `./testConvergence 3 ellipsoid`.