Algoim is a collection of high-order accurate numerical methods and C++ algorithms for working with implicitly defined geometry and level set methods. Motivated by multi-phase multi-physics applications, particularly those with evolving dynamic interfaces, these algorithms target core, fundamental techniques in level set methods. They have been designed with a view to standard finite difference implementations as well as more advanced finite element and discontinuous Galerkin implementations, multi-threading and massively parallel MPI computation. The collection includes:



Many of the numerical algorithms implemented in Algoim are templated on the spatial dimension N, allowing one to develop numerical schemes for any number of spatial dimensions. To assist with this functionality, Algoim makes use of the open source C++ library blitz++ for high-performance fixed-length vector arithmetic, through template<typename T, int N> blitz::TinyVector<T,N>. Therefore, blitz++ should first be downloaded, configured for your compiler, and installed such that it can be found in the appropriate include directories, i.e., so that #include <blitz/array.h> resolves correctly.

In addition, a C++ compiler supporting standard C++14 or higher is required.


Clone via

git clone

or visit the Algoim GitHub page.


Algoim is a header-only C++ library. Except for small example/demonstration applications, all of the files are C++ .hpp header files. As such, it requires minimal installation effort: simply download and configure so that the appropriate header driver can be found by your compiler when you include it in your C++ program, e.g., #include "algoim/src/algoim_quad.hpp".

High-Order Quadrature Algorithms for Implicitly Defined Domains

An implicitly defined domain is either a volumetric region or codimension-one surface whose shape is characterised implicitly by an isosurface of a continuous scalar function. A variety of applications involving implicitly defined geometry require the evaluation of integrals over such domains, including level set methods for propagating interfaces in computational physics, embedded boundary methods for solving partial differential equations on curved domains, and in treating jump conditions and singular source terms in weak formulations. A specific example is that of implicit mesh discontinuous Galerkin methods, which have been developed to facilitate high-order accurate modelling of interfacial fluid dynamics.

In practice, to calculate integrals over implicitly defined domains, a quadrature scheme must be computed. In R. I. Saye, High-Order Quadrature Methods for Implicitly Defined Surfaces and Volumes in Hyperrectangles, SIAM Journal on Scientific Computing, 37(2), A993-A1019 (2015), a general purpose, high-order accurate quadrature algorithm has been developed. Based on the idea of converting the implicitly defined geometry into the graph of an implicitly defined height function, the approach requires only one-dimensional root finding and simple one-dimensional Gaussian quadrature schemes. These algorithms produce quadrature schemes with strictly positive weights and inherits the high-order accuracy of Gaussian quadrature, e.g., with 4 points (per dimension), 8th order accuracy can be achieved. Examples of generated quadrature schemes are shown in the figure.

Some examples illustrating the quadrature schemes constructed by the algorithm (based on a one-dimensional Gaussian quadrature scheme of order 4) for surface integrals (left column) and volume integrals (middle/right columns). The weights are coloured according to a scale that is normalised for each particular case: pale indicates a low-valued weight and dark blue indicates a high-valued weight.


If you make use of these quadrature algorithms in your research, in any documentation or publication please cite - in addition to Algoim - the original paper describing these algorithms:


The driver header file for the quadrature algorithms is located at algoim/src/algoim_quad.hpp. There is only one driver routine, quadGen which is templated on the level set function object phi and the dimension N:

template<int N, typename F>
Algoim::QuadratureRule<N> Algoim::quadGen(const F& phi, const Algoim::BoundingBox<Real,N>& xrange,
                                          int dim, int side, int qo);

The provided parameters are as follows:

The output of quadGen is an Algoim::QuadratureRule<N> object. This object is essentially a std::vector listing the set of quadrature points (as blitz::TinyVector<Real,N>) and their corresponding weights. A QuadratureRule object is also a functor – its associated templated member function operator()(const F& f) can be applied to any user-specified integrand function to compute the result of applying the quadrature rule to a functional.

Requirements on phi. As mentioned above, the phi functor must be able to accept arguments of type blitz::TinyVector<Algoim::Interval<N>,N>. Here, Algoim::Interval<N> is a special type whose purpose is to calculate a first-order Taylor series with bounded remainder, and shares concepts in common with automatic differentiation (see the paper cited above for more details). Interval<N> implements operator+, operator*, operator/, etc., and can be used in a variety of ways, most commonly in evaluating polynomial expressions. Common unary operators are also implemented, e.g., sin(Interval<N>) and exp(Interval<N>). However, one cannot straightforwardly apply Interval<N> arithmetic to max or min statements, if conditions, and other non-standard or non-smooth functions. This relates to the fact that it is very difficult and subtle to compute quadrature schemes for shapes which have corners, holes, or cusps, etc. In using Algoim quadrature schemes for the first time, it is recommended they be applied to smooth level set functions made out of polynomial expressions and common smooth functions like sin, exp, etc.


All of these examples are implemented in algoim/examples/examples_quad.cpp in a single main() entry point; compile them with a C++14 compatible compiler of your choice.

The quadrature algorithms of Algoim are here demonstrated with a level set function describing an ellipse (in N = 2 dimensions) or ellipsoid (in N = 3 dimensions). First, we define a function object implementing the level set function:

template<int N>
struct Ellipsoid
    template<typename T>
    T operator() (const blitz::TinyVector<T,N>& x) const
        if (N == 2)
            return x(0)*x(0) + 4.0*x(1)*x(1) - 1.0;
            return x(0)*x(0) + 4.0*x(1)*x(1) + 9.0*x(2)*x(2) - 1.0;

    template<typename T>
    blitz::TinyVector<T,N> grad(const blitz::TinyVector<T,N>& x) const
        if (N == 2)
            return blitz::TinyVector<T,N>(2.0*x(0), 8.0*x(1));
            return blitz::TinyVector<T,N>(2.0*x(0), 8.0*x(1), 18.0*x(2));

To compute the area of the ellipse in 2D using a scheme with qo = 4, apply quadGen to a bounding box encapsulating the extent of the ellipse, and then apply the resulting quadrature rule to the integrand f(x) = 1:

Ellipsoid<2> phi;
auto q = Algoim::quadGen<2>(phi, Algoim::BoundingBox<double,2>(-1.1, 1.1), -1, -1, 4);
double area = q([](const auto& x) { return 1.0; });
// area ≈ 1.5708

To compute the volume of the ellipsoid in 3D, it is as simple as changing N:

Ellipsoid<3> phi;
auto q = Algoim::quadGen<3>(phi, Algoim::BoundingBox<double,3>(-1.1, 1.1), -1, -1, 4);
double volume = q([](const auto& x) { return 1.0; });
// volume ≈ 0.69813

Each of the two examples above relied on the quadrature scheme’s ability to automatically subdivide the given bounding box until the interface geometry can be represented as the graph of a well-defined height function. The subdivision routine terminates as soon as an internal criterion is met (discussed more in the paper cited above). An alternative method revealing asymptotic high-order accuracy is to subdivide the ellipsoid via a Cartesian grid and apply quadGen to individual grid cells, as follows:

int n = 16;
double dx = 2.2 / n;
Ellipsoid<2> phi;
double area = 0.0;
for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j)
    blitz::TinyVector<double,2> xmin = {-1.1 + i*dx, -1.1 + j*dx};
    blitz::TinyVector<double,2> xmax = {-1.1 + i*dx + dx, -1.1 + j*dx + dx};
    area += Algoim::quadGen<2>(phi, Algoim::BoundingBox<double,2>(xmin, xmax), -1, -1, 4).sumWeights();
// area ≈ 1.570796327

To compute the measure of an implicitly defined codimension-one surface, change the dim parameter to equal N:

Ellipsoid<3> phi;
auto q = Algoim::quadGen<3>(phi, Algoim::BoundingBox<double,3>(-1.1, 1.1), 3, -1, 4);
double surface_area = q.sumWeights();
// surface_area ≈ 4.4007

To visualise a quadrature scheme produced by Algoim::quadGen, a tool is provided in algoim/src/algoim_quad.hpp which outputs a scheme in the format of an XML VTK .vtp file, which one can open with, e.g., ParaView. The routine takes as input a user-defined stream (e.g., a std::ofstream) and writes XML to the stream, interpreting the quadrature scheme as a scattered set of points with associated weights:

#include <fstream>
Algoim::QuadratureRule<N> q = ...;
std::ofstream f("scheme.vtp");
Algoim::outputQuadratureRuleAsVtpXML(q, f);
// Open scheme.vtp in ParaView

More examples of using Algoim’s quadrature code will be added in the future or as requested.

Advanced: Double-Double and Quadruple-Double Precision

Algoim is by default configured to use double-precision arithmetic (Algoim::Real is a typedef for double). For very high precision applications of the associated high-order accurate algorithms, Algoim can interface to the open-source QD library which implements software-driven double-double and quadruple-double arithmetic. To enable this feature:

  1. Download and install QD so that your compiler can resolve the necessary includes, e.g., #include <qd/dd_real.h>.
  2. When compiling code using Algoim, define a global-scope preprocessor directive as follows:
    • To use double-double precision, define ALGOIM_DDREAL. For example, for GCC or ICC, use the compiler argument -DALGOIM_DDREAL.
    • To use quadruple-double precision, define ALGOIM_QDREAL.

When either of the above preprocessor directives are enabled, Algoim::Real ceases to be a typedef for double and is replaced by a typedef for dd_real or qd_real. This also affects other internal workings of Algoim, for example, the precomputed Gaussian quadrature schemes in algoim/src/algoim_gaussquad.hpp are replaced with QD-compatible constants having over 64 digits of accuracy.


Algoim was developed by Robert Saye in the Mathematics Group at the Lawrence Berkeley National Laboratory as part of his research in numerical methods for high-order accurate interface dynamics with applications in multi-phase multi-physics. For more information and applications of these algorithms, visit his research page.

This work was supported in part by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program through the FASTMath Institute under Contract No. DE-AC02-05CH11231 at Lawrence Berkeley National Laboratory.