diff --git a/opm/core/pressure/mimetic/mimetic.h b/opm/core/pressure/mimetic/mimetic.h
index 03d396f95..50f4ddc3a 100644
--- a/opm/core/pressure/mimetic/mimetic.h
+++ b/opm/core/pressure/mimetic/mimetic.h
@@ -20,78 +20,250 @@
#ifndef OPM_MIMETIC_HEADER_INCLUDED
#define OPM_MIMETIC_HEADER_INCLUDED
+/**
+ * \file
+ * Routines to assist mimetic discretisations of the flow equation.
+ */
+
#ifdef __cplusplus
extern "C" {
#endif
-void mim_ip_span_nullspace(int nf, int nconn, int d,
- double *C,
- double *A,
- double *X,
- double *work, int lwork);
-
-void mim_ip_linpress_exact(int nf, int nconn, int d,
- double vol, double *K,
- double *N,
- double *Binv,
- double *work, int lwork);
-
-void mim_ip_simple(int nf, int nconn, int d,
- double v, double *K, double *C,
- double *A, double *N,
- double *Binv,
- double *work, int lwork);
-
-
-/** Compute the mimetic inner products given a grid and cellwise
- * permeability tensors.
+/**
+ * Form linear operator to span the null space of the normal vectors
+ * of a grid cell.
*
- * @param ncells Number of cells in grid.
- * @param d Number of space dimensions.
- * @param max_ncf Maximum number of faces per cell.
- * @param pconn Start indices in conn for each cell, plus end
- * marker. The size of pconn is (ncells + 1), and for a
- * cell i, [conn[pconn[i]], conn[pconn[i+1]]) is a
- * half-open interval containing the indices of faces
- * adjacent to i.
- * @param conn Cell to face mapping. Size shall be equal to the sum of
- * ncf. See pconn for explanation.
- * @param fneighbour Face to cell mapping. Its size shall be equal to
- * the number of faces times 2. For each face, the
- * two entries are either a cell number or -1
- * (signifying the outer boundary). The face normal
- * points out of the first cell and into the second.
- * @param fcentroid Face centroids. Size shall be equal to the number
- * of faces times d.
- * @param fnormal Face normale. Size shall be equal to the number
- * of faces times d.
- * @param farea Face areas.
- * @param ccentroid Cell centroids. Size shall be ncells*d.
- * @param cvol Cell volumes.
- * @param perm Permeability. Size shall be ncells*d*d, storing a
- * d-by-d positive definite tensor per cell.
- * @param[out] Binv This is where the inner product will be
- * stored. Its size shall be equal to \f$\sum_i
- * n_i^2\f$.
+ * Specifically,
+ * \f[
+ * \begin{aligned}
+ * X &= \operatorname{diag}(A) (I - QQ^\mathsf{T})
+ * \operatorname{diag}(A), \\
+ * Q &= \operatorname{orth}(\operatorname{diag}(A) C)
+ * \end{aligned}
+ * \f]
+ * in which \f$\operatorname{orth}(M)\f$ denotes an orthonormal
+ * basis for the colum space (range) of the matrix \f$M\f$,
+ * represented as a matrix.
+ *
+ * @param[in] nf Number of faces connected to single grid cell.
+ * @param[in] nconn Total number of grid cell connections.
+ * Typically equal to @c nf.
+ * @param[in] d Number of physical dimensions.
+ * Assumed less than four.
+ * @param[in,out] C Centroid vectors. Specifically,
+ * \f$c_{ij} = \Bar{x}_{ij} - \Bar{x}_{cj}\f$.
+ * Array of size \f$\mathit{nf}\times d\f$
+ * in column major (Fortran) order.
+ * Contents destroyed on output.
+ * @param[in] A Interface areas.
+ * @param[out] X Null space linear operator. Array of size
+ * \f$\mathit{nconn}\times\mathit{nconn}\f$
+ * in column major (Fortran) order. On output,
+ * the upper left \f$\mathit{nf}\times\mathit{nf}\f$
+ * sub-matrix contains the required null space
+ * linear operator.
+ * @param[out] work Scratch array of size at least @c nconn.
+ * @param[in] lwork Actual size of scratch array.
*/
-void mim_ip_simple_all(int ncells, int d, int max_ncf,
- int *pconn, int *conn,
- int *fneighbour, double *fcentroid, double *fnormal,
- double *farea, double *ccentroid, double *cvol,
- double *perm, double *Binv);
+void
+mim_ip_span_nullspace(int nf, int nconn, int d,
+ double *C,
+ double *A,
+ double *X,
+ double *work, int lwork);
+/**
+ * Form (inverse) mimetic inner product that reproduces linear
+ * pressure drops (constant velocity) on general polyhedral cells.
+ *
+ * Specifically
+ * \f[
+ * B^{-1} = \frac{1}{v} \big(NKN^\mathsf{T} + \frac{6t}{d}\,X\big)
+ * \f]
+ * in which \f$t = \operatorname{tr}(K)\f$ is the trace of \f$K\f$
+ * and \f$X\f$ is the result of function mim_ip_span_nullspace().
+ *
+ * @param[in] nf Number of faces connected to single grid cell.
+ * @param[in] nconn Total number of grid cell connections.
+ * Typically equal to @c nf.
+ * @param[in] d Number of physical dimensions.
+ * Assumed less than four.
+ * @param[in] vol Cell volume.
+ * @param[in] K Permeability. A \f$d\times d\f$ matrix in
+ * column major (Fortran) order.
+ * @param[in] N Normal vectors. An \f$\mathit{nf}\times d\f$
+ * matrix in column major (Fortran) order.
+ * @param[in,out] Binv Inverse inner product result. An
+ * \f$\mathit{nconn}\times\mathit{nconn}\f$
+ * matrix in column major format. On input,
+ * the result of mim_ip_span_nullspace(). On
+ * output, the upper left
+ * \f$\mathit{nf}\times\mathit{nf}\f$ sub-matrix
+ * will be overwritten with \f$B^{-1}\f$.
+ * @param[in,out] work Scratch array of size at least nf * d.
+ * @param[in] lwork Actual size of scratch array.
+ */
+void
+mim_ip_linpress_exact(int nf, int nconn, int d,
+ double vol, double *K,
+ double *N,
+ double *Binv,
+ double *work, int lwork);
+
+
+/**
+ * Convenience wrapper around the function pair mim_ip_span_nullspace()
+ * and mim_ip_linpress_exact().
+ *
+ * @param[in] nf Number of faces connected to single grid cell.
+ * @param[in] nconn Total number of grid cell connections.
+ * Typically equal to @c nf.
+ * @param[in] d Number of physical dimensions.
+ * Assumed less than four.
+ * @param[in] v Cell volume.
+ * @param[in] K Permeability. A \f$d\times d\f$ matrix in
+ * column major (Fortran) order.
+ * @param[in,out] C Centroid vectors. Specifically,
+ * \f$c_{ij} = \Bar{x}_{ij} - \Bar{x}_{cj}\f$.
+ * Array of size \f$\mathit{nf}\times d\f$
+ * in column major (Fortran) order.
+ * Contents destroyed on output.
+ * @param[in] A Interface areas.
+ * @param[in] N Outward normal vectors.
+ * An \f$\mathit{nf}\times d\f$ matrix in
+ * column major (Fortran) order.
+ * @param[out] Binv Inverse inner product result. An
+ * \f$\mathit{nconn}\times\mathit{nconn}\f$
+ * matrix in column major format. On
+ * output, the upper left
+ * \f$\mathit{nf}\times\mathit{nf}\f$ sub-matrix
+ * will be overwritten with \f$B^{-1}\f$
+ * defined by function mim_ip_linpress_exact().
+ * @param[in,out] work Scratch array of size at least nf * d.
+ * @param[in] lwork Actual size of scratch array.
+ */
+void
+mim_ip_simple(int nf, int nconn, int d,
+ double v, double *K, double *C,
+ double *A, double *N,
+ double *Binv,
+ double *work, int lwork);
+
+
+/**
+ * Compute the mimetic inner products given a grid and cell-wise
+ * permeability tensors.
+ *
+ * This function applies mim_ip_simple() to all specified cells.
+ *
+ * @param[in] ncells Number of cells.
+ * @param[in] d Number of physical dimensions.
+ * @param[in] max_ncf Maximum number of connections (faces)
+ * of any individual cell.
+ * @param[in] pconn Start pointers of cell-to-face topology
+ * mapping.
+ * @param[in] conn Actual cell-to-face topology mapping.
+ * @param[in] fneighbour Face-to-cell mapping.
+ * @param[in] fcentroid Face centroids.
+ * @param[in] fnormal Face normals.
+ * @param[in] farea Face areas.
+ * @param[in] ccentroid Cell centroids.
+ * @param[in] cvol Cell volumes.
+ * @param[in] perm Cell permeability.
+ * @param[out] Binv Inverse inner product result. Must point
+ * to an array of size at least
+ * \f$\sum_c n_c^2\f$ when \f$n_c\f$ denotes
+ * the number of connections (faces) of
+ * cell \f$c\f$.
+ */
+void
+mim_ip_simple_all(int ncells, int d, int max_ncf,
+ int *pconn, int *conn,
+ int *fneighbour, double *fcentroid, double *fnormal,
+ double *farea, double *ccentroid, double *cvol,
+ double *perm, double *Binv);
+
+/**
+ * Compute local, static gravity pressure contributions to Darcy
+ * flow equation discretised using a mimetic finite-difference method.
+ *
+ * The pressure contribution of local face \f$i\f$ in cell \f$c\f$ is
+ * \f[
+ * \mathit{gpress}_{\mathit{pconn}_c + i} =
+ * \vec{g}\cdot (\Bar{x}_{\mathit{conn}_{\mathit{pconn}_c + i}}
+ * - \Bar{x}_c)
+ * \f]
+ *
+ * @param[in] nc Number of cells.
+ * @param[in] d Number of physcial dimensions.
+ * @param[in] grav Gravity vector. Array of size @c d.
+ * @param[in] pconn Start pointers of cell-to-face topology
+ * mapping.
+ * @param[in] conn Actual cell-to-face topology mapping.
+ * @param[in] fcentroid Face centroids.
+ * @param[in] ccentroid Cell centroids.
+ * @param[out] gpress Gravity pressure result. Array of size
+ * at least pconn[nc].
+ */
void
mim_ip_compute_gpress(int nc, int d, const double *grav,
const int *pconn, const int *conn,
const double *fcentroid, const double *ccentroid,
double *gpress);
-/* inv(B) <- \lambda_t(s)*inv(B)_0 */
+
+/**
+ * Incorporate effects of multiple phases in mimetic discretisation of
+ * flow equations.
+ *
+ * Specifically, update the (inverse) inner products \f$B^{-1}\f$
+ * previously computed using function mim_ip_linpress_exact() according
+ * to the rule
+ * \f[
+ * \Tilde{B}_c^{-1} = \frac{1}{\lambda_{T,c}} B_c^{-1},
+ * \quad i=0,\dots,\mathit{nc}-1
+ * \f]
+ * in which \f$B_c^{-1}\f$ denotes the result of mim_ip_linpress_exact()
+ * for cell \f$c\f$ and \f$\lambda_{T,c}\f$ denotes the total mobility
+ * of cell \f$c\f$.
+ *
+ * @param[in] nc Number of cells.
+ * @param[in] pconn Start pointers of cell-to-face topology
+ * mapping.
+ * @param[in] totmob Total mobility for all cells. Array of size @c nc.
+ * @param[in] Binv0 Inverse inner product results for all cells.
+ * @param[out] Binv Inverse inner product results incorporating
+ * effects of multiple fluid phases.
+ */
void
mim_ip_mobility_update(int nc, const int *pconn, const double *totmob,
const double *Binv0, double *Binv);
-/* G <- \sum_i \rho_i f_i(s) * G_0 */
+
+/**
+ * Incorporate effects of multiple fluid phases into existing, local,
+ * static mimetic discretisations of gravity pressure.
+ *
+ * Specifically, update the result of mim_ip_compute_gpress()
+ * according to the rule
+ * \f[
+ * \Tilde{G}_{\mathit{pconn}_c + i} = \omega_c\cdot
+ * G_{\mathit{pconn}_c + i}, \quad i=\mathit{pconn}_c, \dots,
+ * \mathit{pconn}_{c+1}-1, \quad c=0,\dots,\mathit{nc}-1
+ * \f]
+ * in which \f$\omega_c = (\sum_\alpha \lambda_{\alpha,c}
+ * \rho_\alpha)/\lambda_{T,c}\f$ and \f$\Tilde{G}\f$ denotes the result
+ * of function mim_ip_compute_gpress().
+ *
+ * @param[in] nc Number of cells.
+ * @param[in] pconn Start pointers of cell-to-face topology
+ * mapping.
+ * @param[in] omega Sum of phase densities weighted by
+ * fractional flow.
+ * @param[in] gpress0 Result of mim_ip_compute_gpress().
+ * @param[out] gpress Gravity pressure incorporating effects
+ * of multiple fluid phases.
+ */
void
mim_ip_density_update(int nc, const int *pconn, const double *omega,
const double *gpress0, double *gpress);