diff --git a/mex_ip_simple.c b/mex_ip_simple.c
index eb67f256..a5e291b2 100644
--- a/mex_ip_simple.c
+++ b/mex_ip_simple.c
@@ -1,4 +1,6 @@
+#include <string.h>
 #include <mex.h>
+
 #include "mrst_api.h"
 #include "mimetic.h"
 
@@ -76,23 +78,30 @@ verify_structural_consistency(int nlhs, int nrhs, const mxArray *prhs[])
 {
     int rock_ok;
     int grid_ok;
+    int conn_ok;
 
-    mxAssert((nlhs == 1) && (nrhs == 2),
-             "Must be called with exactly two inputs and one output.");
+    mxAssert((nlhs == 1) && (nrhs == 4),
+             "Must be called with exactly four inputs and one output.");
 
     mxAssert(mxIsStruct(prhs[0]) && mxIsStruct(prhs[1]),
-             "Inputs must be STRUCTs.");
+             "First two inputs must be STRUCTs.");
 
     mxAssert((mxGetNumberOfElements(prhs[0]) == 1) &&
              (mxGetNumberOfElements(prhs[1]) == 1),
-             "Inputs must be 1-by-1 STRUCT arrays.");
+             "First two inputs must be 1-by-1 STRUCT arrays.");
 
     grid_ok = verify_grid_structure(prhs[0]);
 
     /* rock must contain a field 'perm' */
     rock_ok = mxGetFieldNumber(prhs[1], "perm") >= 0;
 
-    return grid_ok && rock_ok;
+    /* connection structure */
+    conn_ok = (mxIsDouble(prhs[2]) || mxIsInt32(prhs[2])) &&
+              (mxIsDouble(prhs[3]) || mxIsInt32(prhs[3])) &&
+              (mxGetNumberOfElements(prhs[2]) <
+               mxGetNumberOfElements(prhs[3]));
+    
+    return grid_ok && rock_ok && conn_ok;
 }
 
 
@@ -122,11 +131,72 @@ deallocate_face_data(int *fneighbour, double *fnormal, double *fcentroid)
 }
 
 
+/* ------------------------------------------------------------------ */
+static void
+copy_int_vector(const mxArray *M_v, int *v)
+/* ------------------------------------------------------------------ */
+{
+    size_t n, i;
+
+    int    *pi;
+    double *pd;
+
+    n = mxGetNumberOfElements(M_v);
+
+    if (mxIsDouble(M_v)) {
+        pd = mxGetPr(M_v);
+
+        for (i = 0; i < n; i++) {
+            mxAssert ((INT_MIN <= pd[i]) && (pd[i] <= INT_MAX),
+                      "Data element outside range for INT32");
+
+            v[i] = pd[i];
+        }
+    } else {
+        pi = mxGetData(M_v);
+        memcpy(v, pi, n * sizeof *v);
+    }
+}
+
+
+/* ------------------------------------------------------------------ */
+static void
+extract_connection_data(const mxArray *M_nconn, const mxArray *M_conn,
+                        const int nc, int **nconn, int **conn,
+                        int *sum_nconn, int *sum_nconn2)
+/* ------------------------------------------------------------------ */
+{
+    size_t i, n;
+
+    n = mxGetNumberOfElements(M_conn);
+
+    *nconn = mxMalloc(nc * sizeof *nconn);
+    *conn  = mxMalloc(n  * sizeof *conn );
+
+    copy_int_vector(M_nconn, *nconn);
+    copy_int_vector(M_conn,  *conn );
+
+    /* Adjust connections for 1-based indexing. */
+    for (i = 0; i < n; i++) {
+        (*conn)[i] -= 1;
+    }
+
+    /* Count connections */
+    n = nc;
+    *sum_nconn = *sum_nconn2 = 0;
+    for (i = 0; i < n; i++) {
+        *sum_nconn  += (*nconn)[i];
+        *sum_nconn2 += (*nconn)[i] * (*nconn)[i];
+    }
+}
+
+
 /* ------------------------------------------------------------------ */
 static int
-extract_cell_data(const mxArray *G, int d,
-                  int *max_ncf, int *sum_ncf, int *sum_ncf2,
-                  int **ncfaces, int **cfaces,
+extract_cell_data(const mxArray *G,
+                  const mxArray *M_nconn, const mxArray *M_conn,
+                  int d, int *max_ncf, int *sum_nconn, int *sum_nconn2,
+                  int **ncfaces, int **nconn, int **conn,
                   double **ccentroids, double **cvolumes)
 /* ------------------------------------------------------------------ */
 {
@@ -134,7 +204,10 @@ extract_cell_data(const mxArray *G, int d,
    int ncells = getNumberOfCells(G);
    int i;
 
-   *max_ncf = *sum_ncf = *sum_ncf2 = 0;
+   extract_connection_data(M_nconn, M_conn, ncells,
+                           nconn, conn, sum_nconn, sum_nconn2);
+   
+   *max_ncf = 0;
 
    for(i=0; i<ncells; ++i)
    {
@@ -142,11 +215,8 @@ extract_cell_data(const mxArray *G, int d,
       (*ncfaces)[i] = n;
       
       *max_ncf   = MAX(*max_ncf, n);
-      *sum_ncf  += n; 
-      *sum_ncf2 += n*n;       
    }
 
-   *cfaces     = getCellFaces(G);
    *ccentroids = getCellCentroids(G);
    *cvolumes   = getCellVolumes(G);
    
@@ -156,11 +226,12 @@ extract_cell_data(const mxArray *G, int d,
 
 /* ------------------------------------------------------------------ */
 static void
-deallocate_cell_data(int *ncfaces, int *cfaces, double *ccentroids)
+deallocate_cell_data(int *ncfaces, int *nconn, int *conn, double *ccentroids)
 /* ------------------------------------------------------------------ */
 {
     mxFree(ccentroids);
-    mxFree(cfaces);
+    mxFree(conn);
+    mxFree(nconn);
     mxFree(ncfaces);
 }
 
@@ -175,7 +246,7 @@ deallocate_perm_data(double *K)
 
 
 /*
- * BI = mex_ip_simple(G, rock)
+ * BI = mex_ip_simple(G, rock, nconn, conn)
  */
 
 /* ------------------------------------------------------------------ */
@@ -189,7 +260,7 @@ mexFunction(int nlhs,       mxArray *plhs[],
     int *fneighbour;
     double *farea, *fnormal, *fcentroid;
 
-    int *ncfaces, *cfaces;
+    int *ncfaces, *nconn, *conn;
     double *ccentroids, *cvolumes;
 
     double *perm;
@@ -203,25 +274,26 @@ mexFunction(int nlhs,       mxArray *plhs[],
         nfaces = extract_face_data(prhs[0], d,
                                    &fneighbour, &farea, &fnormal, &fcentroid);
 
-        ncells = extract_cell_data(prhs[0], d,
+        ncells = extract_cell_data(prhs[0], prhs[2], prhs[3], d,
                                    &max_ncf, &sum_ncf, &sum_ncf2, &ncfaces,
-                                   &cfaces, &ccentroids, &cvolumes);
+                                   &nconn, &conn, &ccentroids, &cvolumes);
 
         perm = getPermeability(mxGetField(prhs[1], 0, "perm"), d);
 
         if ((ncells > 0) && (nfaces > 0)) {
             plhs[0] = mxCreateDoubleMatrix(sum_ncf2, 1, mxREAL);
-            Binv    = mxGetPr(plhs[0]);
+            Binv    = mxGetPr(plhs[0]); /* mxCreateDoubleMatrix == ZEROS */
 
-            mim_ip_simple_all(ncells, d, max_ncf, ncfaces, cfaces, fneighbour,
-                              fcentroid, fnormal, farea, ccentroids, cvolumes,
-                              perm, Binv);
+            /* mim_ip_simple_all() requires zeroed target array... */
+            mim_ip_simple_all(ncells, d, max_ncf, ncfaces, nconn, conn,
+                              fneighbour, fcentroid, fnormal, farea,
+                              ccentroids, cvolumes, perm, Binv);
         } else {
             plhs[0] = mxCreateDoubleScalar(mxGetNaN());
         }
 
         deallocate_perm_data(perm);
-        deallocate_cell_data(ncfaces, cfaces, ccentroids);
+        deallocate_cell_data(ncfaces, nconn, conn, ccentroids);
         deallocate_face_data(fneighbour, fnormal, fcentroid);
     } else {
         plhs[0] = mxCreateDoubleScalar(mxGetNaN());
diff --git a/mex_ip_simple.m b/mex_ip_simple.m
index 728e57ca..b33e0579 100644
--- a/mex_ip_simple.m
+++ b/mex_ip_simple.m
@@ -2,16 +2,24 @@ function varargout = mex_ip_simple(varargin)
 %Compute 'ip_simple' inner product values using compiled C code.
 %
 % SYNOPSIS:
-%   BI = mex_ip_simple(G, rock)
+%   BI = mex_ip_simple(G, rock, nconn, conn)
 %
 % PARAMETERS:
-%   G    - Grid data structure.
+%   G     - Grid data structure.
 %
-%   rock - Rock data structure.  Must contain a valid field 'perm'.
+%   rock  - Rock data structure.  Must contain a valid field 'perm'.
+%
+%   nconn - Number of connections per cell.  Often coincides with
+%           DIFF(G.cells.facePos), but may be larger if any cells are
+%           perforated by one or more wells.
+%
+%   conn  - Connection data per cell.  Often coincides with
+%           G.cells.faces(:,1) but will contain additional data if a cell is
+%           perforated by one or more wells.
 %
 % RETURNS:
-%   BI   - A SUM(DIFF(G.cells.facePos) .^ 2)-by-1 array of inner product
-%          values, ordered by the cells of the input grid.
+%   BI   - A SUM(nconn .^ 2)-by-1 array of inner product values, ordered by
+%          the cells of the input grid.
 %
 % NOTE:
 %   As the return value 'BI' is but a simple data array value, it must be
@@ -30,12 +38,14 @@ function varargout = mex_ip_simple(varargin)
 %   rock.perm = convertFrom(rock.perm(G.cells.indexMap, :), ...
 %                           milli*darcy);
 %
+%   nconn = diff(G.cells.facePos);
+%   conn  = G.cells.faces(:,1);
+%
 %   t0 = tic;
-%   BI = mex_ip_simple(G, rock);
+%   BI = mex_ip_simple(G, rock, nconn, conn);
 %   toc(t0)
 %
-%   [i, j] = blockDiagIndex(diff(G.cells.facePos), ...
-%                           diff(G.cells.facePos));
+%   [i, j] = blockDiagIndex(nconn, nconn);
 %
 %   S = struct('BI', sparse(i, j, BI), 'type', 'hybrid', 'ip', 'ip_simple')
 %
diff --git a/mimetic.c b/mimetic.c
index b4107180..195ac7c3 100644
--- a/mimetic.c
+++ b/mimetic.c
@@ -1,9 +1,9 @@
+#include <assert.h>
 #include <stddef.h>
 #include <stdlib.h>
 
 #if defined(COMPILING_FOR_MATLAB) && COMPILING_FOR_MATLAB
-#include "mex.h"
-#include "matrix.h"
+#include <mex.h>
 #define MAT_SIZE_T mwSignedIndex
 #endif
 
@@ -16,8 +16,9 @@
 
 /* ------------------------------------------------------------------ */
 void
-mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
-                  int *fneighbour,  double *fcentroid, double *fnormal,
+mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
+                  int *nconn, int *cf,
+                  int *fneighbour, double *fcentroid, double *fnormal,
                   double *farea, double *ccentroid, double *cvol,
                   double *perm, double *Binv)
 /* ------------------------------------------------------------------ */
@@ -56,11 +57,11 @@ mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
                 }
             }
 
-            mim_ip_simple(nf, d, cvol[c], A, &perm[c * d * d],
-                          C, N, &Binv[fpos2], work, lwork);
+            mim_ip_simple(nf, nconn[c], d, cvol[c], &perm[c * d * d],
+                          C, A, N, &Binv[fpos2], work, lwork);
 
-            fpos  += nf;
-            fpos2 += nf * nf;
+            fpos  += nconn[c];
+            fpos2 += nconn[c] * nconn[c];
         }
     }
 
@@ -70,32 +71,41 @@ mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
 
 /* ------------------------------------------------------------------ */
 void
-mim_ip_simple(int nf, int d, double vol,
-              double *A, double *K, double *C, double *N,
-              double *Binv, double *work, int lwork)
+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)
 /* ------------------------------------------------------------------ */
 {
-    mim_ip_span_nullspace(nf, d, C, A, Binv, work, lwork);
-    mim_ip_linpress_exact(nf, d, vol, N, K, Binv, C);
+    mim_ip_span_nullspace(nf, nconn, d,    C, A, Binv, work, lwork);
+    mim_ip_linpress_exact(nf, nconn, d, v, K, N, Binv, work, lwork);
 }
 
 
 /* ------------------------------------------------------------------ */
 void
-mim_ip_span_nullspace(int nf, int d, double *C, double *A,
-                      double *X, double *work, int nwork)
+mim_ip_span_nullspace(int nf, int nconn, int d,
+                      double *C,
+                      double *A,
+                      double *X,
+                      double *work, int nwork)
 /* ------------------------------------------------------------------ */
 {
-    MAT_SIZE_T m, n, k, ld, info, lwork;
+    MAT_SIZE_T m, n, k, ldC, ldX, info, lwork;
 
     int    i, j;
     double a1, a2;
 
     double tau[3] = { 0.0 };  /* No more than 3 spatial dimensions */
 
-    /* Step 1) X <- I */
-    for (i = 0; i < nf*nf; i++) { X[i           ] = 0.0; }
-    for (i = 0; i < nf   ; i++) { X[i * (nf + 1)] = 1.0; }
+    /* Step 1) X(1:nf, 1:nf) <- I_{nf} */
+    for (j = 0; j < nf; j++) {
+        for (i = 0; i < nf; i++) {
+            X[i + j*nconn] = 0.0;
+        }
+        X[j * (nconn + 1)] = 1.0;
+    }
 
     /* Step 2) C <- orth(A * C) */
     for (j = 0; j < d; j++) {
@@ -104,24 +114,25 @@ mim_ip_span_nullspace(int nf, int d, double *C, double *A,
         }
     }
 
-    m = nf;  n = d;  ld = nf;  k = d;  lwork = nwork;
-    dgeqrf_(&m, &n,     C, &ld, tau, work, &lwork, &info);
-    dorgqr_(&m, &n, &k, C, &ld, tau, work, &lwork, &info);
+    m = nf;  n = d;  ldC = nf;  k = d;  lwork = nwork;
+    dgeqrf_(&m, &n,     C, &ldC, tau, work, &lwork, &info);
+    dorgqr_(&m, &n, &k, C, &ldC, tau, work, &lwork, &info);
 
     /* Step 3) X <- A * (X - C*C') * A */
+    ldX = nconn;
     a1 = -1.0;  a2 = 1.0;
     dsyrk_("Upper Triangular", "No Transpose",
-           &m, &n, &a1, C, &ld, &a2, X, &ld);
+           &m, &n, &a1, C, &ldC, &a2, X, &ldX);
     for (j = 0; j < nf; j++) {
         for (i = 0; i <= j; i++) {
-            X[i + j*nf] *= A[i] * A[j];
+            X[i + j*nconn] *= A[i] * A[j];
         }
     }
 
     /* Account for DSYRK only assigning upper triangular part. */
     for (j = 0; j < nf; j++) {
         for (i = j + 1; i < nf; i++) {
-            X[i + j*nf] = X[j + i*nf];
+            X[i + j*nconn] = X[j + i*nconn];
         }
     }
 }
@@ -129,14 +140,19 @@ mim_ip_span_nullspace(int nf, int d, double *C, double *A,
 
 /* ------------------------------------------------------------------ */
 void
-mim_ip_linpress_exact(int nf, int d, double vol, double *N, double *K,
-                      double *Binv, double *T)
+mim_ip_linpress_exact(int nf, int nconn, int d,
+                      double vol, double *K,
+                      double *N,
+                      double *Binv,
+                      double *work, int lwork)
 /* ------------------------------------------------------------------ */
 {
     MAT_SIZE_T m, n, k, ld1, ld2, ldBinv;
     int        i;
     double     a1, a2, t;
 
+    assert (lwork >= d * nf);
+
     t = 0.0;
     for (i = 0; i < d; i++) {
         t += K[i + i*d];
@@ -147,12 +163,12 @@ mim_ip_linpress_exact(int nf, int d, double vol, double *N, double *K,
     ld1 = nf ;  ld2 = d  ;
     a1  = 1.0;  a2  = 0.0;
     dgemm_("No Transpose", "No Transpose", &m, &n, &k,
-           &a1, N, &ld1, K, &ld2, &a2, T, &ld1);
+           &a1, N, &ld1, K, &ld2, &a2, work, &ld1);
 
     /*  Step 5) Binv <- (N*K*N' + t*X) / vol */
     a1 = 1.0     /      vol ;
     a2 = 6.0 * t / (d * vol);
-    ldBinv = nf;
+    ldBinv = nconn;
     dgemm_("No Transpose", "Transpose", &m, &m, &n,
-           &a1, T, &ld1, N, &ld1, &a2, Binv, &ldBinv);
+           &a1, work, &ld1, N, &ld1, &a2, Binv, &ldBinv);
 }
diff --git a/mimetic.h b/mimetic.h
index 7f0623ef..95a1a842 100644
--- a/mimetic.h
+++ b/mimetic.h
@@ -2,17 +2,26 @@
 #define MIMETIC_H_INCLUDED
 
 
-void mim_ip_span_nullspace(int nf, int d, double *C, double *A,
-                           double *X, double *work, int lwork);
+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 d, double vol, double *N,
-                           double *K, double *Binv, double *T);
+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 d, double v, double *A, double *K,
-                   double *C, double *N, double *Binv,
+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);
 
-void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf, int *cf,
+void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
+                       int *nconn, int *conn,
                        int *fneighbour, double *fcentroid, double *fnormal,
                        double *farea, double *ccentroid, double *cvol,
                        double *perm, double *Binv);