diff --git a/mrst_api.c b/mrst_api.c
index b51580f9..a188a6c4 100644
--- a/mrst_api.c
+++ b/mrst_api.c
@@ -1,10 +1,7 @@
 /* "API" */
-#include <stdlib.h>
 #include <mex.h>
 #include "mrst_api.h"
 
-#define Assert mxAssert
-
 /* ------------------------------------------------------------------ */
 static mxArray*
 getField(const mxArray *a, char *field, char *subfield)
@@ -25,8 +22,8 @@ static int *
 extractIntMatrix(const mxArray *a)
 /* ------------------------------------------------------------------ */
 {
-   int n = mxGetNumberOfElements(a);
-   int *q         = malloc(n * sizeof *q);
+   int n  = mxGetNumberOfElements(a);
+   int *q = mxMalloc(n * sizeof *q);
    if (q != NULL)
    {
       if (mxIsInt32(a))
@@ -35,8 +32,6 @@ extractIntMatrix(const mxArray *a)
          int i,j;
          for (i=0; i<n; ++i)
          {
-            mxAssert (p[i] <= INT_MAX,
-                      "Matrix entry exceeds INT_MAX");
             q[i] = p[i]-1;
          }      
       }
@@ -46,7 +41,7 @@ extractIntMatrix(const mxArray *a)
          int i,j;
          for (i=0; i<n; ++i)
          {
-            mxAssert (p[i] <= INT_MAX,
+            mxAssert ((1 <= p[i]) && (p[i] <= INT_MAX),
                       "Matrix entry exceeds INT_MAX");
             q[i] = p[i]-1;
          }      
@@ -61,48 +56,80 @@ static int *
 extractIntMatrixTranspose(const mxArray *a)
 /* ------------------------------------------------------------------ */
 {
-   int *q = extractIntMatrix(a);
+   int *p = extractIntMatrix(a);
    int  M = mxGetM(a);
    int  N = mxGetN(a);
-
-   int i,j;
-   for(i=0; i<M; ++i)
+   
+   int *q = mxMalloc(M * N * sizeof *q);
+   if (q != NULL)
    {
-      for(j=0; j<i; ++j)
+      int i,j;
+      for(i=0; i<M; ++i)
       {
-         int tmp  = q[i+M*j];
-         q[i+M*j] = q[j+N*i];
-         q[j+N*i] = tmp;
+         for(j=0; j<N; ++j)
+         {
+            q[i*N+j] = p[i+M*j];
+         }
+      }
+   }
+   mxFree(p);
+   return q;
+}
+
+/* ------------------------------------------------------------------ */
+static double *
+extractDoubleMatrixTranspose(const mxArray *a)
+/* ------------------------------------------------------------------ */
+{
+   int  M    = mxGetM(a);
+   int  N    = mxGetN(a);
+   double *q = mxMalloc(M * N * sizeof *q);
+   if (q != NULL)
+   {
+      double *p = mxGetPr(a);
+      int i,j;
+      for(i=0; i<M; ++i)
+      {
+         for(j=0; j<N; ++j)
+         {
+            q[i*N+j] = p[i+M*j];
+         }
       }
    }
    return q;
 }
 
 /* ------------------------------------------------------------------ */
-int  getNumberOfDimensions(const mxArray *G)
+int  
+getNumberOfDimensions(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {   
    return mxGetN(getField(G, "nodes", "coords"));
 }
 
 /* ------------------------------------------------------------------ */
-void getLocal2GlobalCellMap(const mxArray *G)
+void 
+getLocal2GlobalCellMap(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-   Assert(0, "Not implemented!");
+   mxAssert(0, "Not implemented!");
 }
 
 /* ------------------------------------------------------------------ */
-int getNumberOfNodes(const mxArray *G)
+int 
+getNumberOfNodes(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
    return mxGetM(getField(G, "nodes", "coords"));
 }
 
 /* ------------------------------------------------------------------ */
-double *getNodeCoordinates(const mxArray *G)
+double *
+getNodeCoordinates(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
+   return extractDoubleMatrixTranspose(getField(G, "nodes", "coords"));
+#if 0
    mxArray *p1, *p2;
 
    p1 = mxGetField(G , 0, "nodes" );
@@ -111,7 +138,7 @@ double *getNodeCoordinates(const mxArray *G)
    const int n = getNumberOfNodes(G);
    const int d = getNumberOfDimensions(G);
 
-   double *v = malloc(n * d * sizeof *v);
+   double *v = mxMalloc(n * d * sizeof *v);
    if (v != NULL)
    {
 
@@ -126,31 +153,36 @@ double *getNodeCoordinates(const mxArray *G)
       }
    }
    return v;
+#endif
 }
 
 /* ------------------------------------------------------------------ */
-int  getNumberOfFaces(const mxArray *G)
+int  
+getNumberOfFaces(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
    return mxGetNumberOfElements(getField(G, "faces", "nodePos"))-1;
 }
 
 /* ------------------------------------------------------------------ */
-int *getFaceNodePos(const mxArray *G)
+int *
+getFaceNodePos(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
    return extractIntMatrix(getField(G, "faces", "nodePos"));
 }
 
 /* ------------------------------------------------------------------ */
-int getNumberOfFaceNodes(const mxArray *G)
+int 
+getNumberOfFaceNodes(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
      return mxGetNumberOfElements(getField(G, "faces", "nodes"));
 }
 
 /* ------------------------------------------------------------------ */
-int *getFaceNodes(const mxArray *G)
+int *
+getFaceNodes(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
      return extractIntMatrix(getField(G, "faces", "nodes"));
@@ -158,74 +190,184 @@ int *getFaceNodes(const mxArray *G)
 
 
 /* ------------------------------------------------------------------ */
-int *getFaceCellNeighbors(const mxArray *G)
+int *
+getFaceCellNeighbors(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
      return extractIntMatrixTranspose(getField(G, "faces", "neighbors"));
 }
 
 /* ------------------------------------------------------------------ */
-void getFaceAreas(const mxArray *G, double **v)
+double * 
+getFaceAreas(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     *v = mxGetPr(getField(G, "faces", "areas"));
+   return mxGetPr(getField(G, "faces", "areas"));
 }
 
 /* ------------------------------------------------------------------ */
-void getFaceNormals(const mxArray *G, double **v)
+double *
+getFaceNormals(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     *v = mxGetPr(getField(G, "faces", "normals"));
+   return extractDoubleMatrixTranspose(getField(G, "faces", "normals"));
 }
 
 /* ------------------------------------------------------------------ */
-void getFaceCentroids(const mxArray *G, double **v)
+double *
+getFaceCentroids(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     *v = mxGetPr(getField(G, "faces", "centroids"));
+   return extractDoubleMatrixTranspose(getField(G, "faces", "centroids"));
 }
 
 /* ------------------------------------------------------------------ */
-int  getNumberOfCells(const mxArray *G)
+int  
+getNumberOfCells(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     return mxGetNumberOfElements(getField(G, "cells", "facePos"))-1;   
+   return mxGetNumberOfElements(getField(G, "cells", "facePos"))-1;   
 }
 
 /* ------------------------------------------------------------------ */
 int *getCellFacePos(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     return extractIntMatrix(getField(G, "cells", "facePos"));
+   return extractIntMatrix(getField(G, "cells", "facePos"));
 }
 
 /* ------------------------------------------------------------------ */
 int getNumberOfCellFaces(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     return mxGetNumberOfElements(getField(G, "cells", "faces"))-1;   
+   return mxGetNumberOfElements(getField(G, "cells", "faces"))-1;   
 }
 
 /* ------------------------------------------------------------------ */
 int *getCellFaces(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     return extractIntMatrix(getField(G, "cells", "faces"));
+   return extractIntMatrix(getField(G, "cells", "faces"));
 }
 
 /* ------------------------------------------------------------------ */
-void getCellVolumes(const mxArray *G, double **v)
+double *
+getCellVolumes(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     *v = mxGetPr(getField(G, "cells", "volumes"));
+   return mxGetPr(getField(G, "cells", "volumes"));
 }
 
 
 
 /* ------------------------------------------------------------------ */
-void getCellCentroids(const mxArray *G, double **v)
+double *
+getCellCentroids(const mxArray *G)
 /* ------------------------------------------------------------------ */
 {
-     *v = mxGetPr(getField(G, "cells", "centroids"));
+   return extractDoubleMatrixTranspose(getField(G, "cells", "centroids"));
 }
 
+
+
+/*
+ *
+ *
+ *
+ */
+
+
+/* ------------------------------------------------------------------ */
+static int
+allocate_perm_data(int ncells, int d, double **K)
+/* ------------------------------------------------------------------ */
+{
+    int ret;
+    double *perm;
+
+    perm = mxMalloc(ncells * d * d * sizeof *perm);
+
+    if (perm != NULL) {
+        *K = perm;
+        ret = 1;
+    } else {
+        ret = 0;
+    }
+
+    return ret;
+}
+
+
+/* ------------------------------------------------------------------ */
+double *
+getPermeability(const mxArray *perm, int d)
+/* ------------------------------------------------------------------ */
+{
+    int ncells, ncomp, alloc_ok;
+
+    int c, i, off;
+
+    double *k, *tensor;
+
+    ncells = mxGetM(perm);
+    ncomp  = mxGetN(perm);
+
+    alloc_ok = allocate_perm_data(ncells, d, &k);
+
+    if (alloc_ok) {
+        for (i = 0; i < ncells * d * d; i++) {
+            k[i] = 0.0;
+        }
+
+        tensor = mxGetPr(perm);
+
+        if (ncomp == 1) {
+            /* Isotropic (scalar) tensor */
+            for (c = 0; c < ncells; c++) {
+                off = c * d * d;
+                for (i = 0; i < d; i++) {
+                    k[i*(d + 1) + off] = tensor[c];
+                }
+            }
+        } else if (ncomp == d) {
+            /* Diagonal tensor */
+            for (c = 0; c < ncells; c++) {
+                off = c * d * d;
+                for (i = 0; i < d; i++) {
+                    k[i*(d + 1) + off] = tensor[c + i*ncells];
+                }
+            }
+        } else if (d == 2) {
+            /* Full 2D tensor */
+           mxAssert (ncomp == 3, "");
+
+            for (c = 0; c < ncells; c++) {
+                off = c * d * d;
+                k[0 + off] = tensor[c + 0*ncells];
+                k[1 + off] = tensor[c + 1*ncells];
+                k[2 + off] = tensor[c + 1*ncells];
+                k[3 + off] = tensor[c + 2*ncells];
+            }
+        } else {
+            /* Full 3D tensor */
+            mxAssert ((d == 3) && (ncomp == 6), "");
+
+            for (c = 0; c < ncells; c++) {
+                off = c * d * d;
+
+                k[0 + off] = tensor[c + 0*ncells];
+                k[1 + off] = tensor[c + 1*ncells];
+                k[2 + off] = tensor[c + 2*ncells];
+
+                k[3 + off] = tensor[c + 1*ncells];
+                k[4 + off] = tensor[c + 3*ncells];
+                k[5 + off] = tensor[c + 4*ncells];
+
+                k[6 + off] = tensor[c + 2*ncells];
+                k[7 + off] = tensor[c + 4*ncells];
+                k[8 + off] = tensor[c + 5*ncells];
+            }
+        }
+    }
+    return k;
+}
diff --git a/mrst_api.h b/mrst_api.h
index 8f36183b..2351e72c 100644
--- a/mrst_api.h
+++ b/mrst_api.h
@@ -20,9 +20,9 @@ int    *getFaceNodes          (const mxArray *G);    /* copy */
 int    *getFaceCellNeighbors  (const mxArray *G);    /* copy */
 
 /* Face geometry */
-void    getFaceAreas          (const mxArray *G, double **v);
-void    getFaceNormals        (const mxArray *G, double **v);
-void    getFaceCentroids      (const mxArray *G, double **v);
+double *getFaceAreas          (const mxArray *G);
+double *getFaceNormals        (const mxArray *G);    /* copy */
+double *getFaceCentroids      (const mxArray *G);    /* copy */
 
 /* Cell topology */
 int     getNumberOfCells      (const mxArray *G);
@@ -31,7 +31,12 @@ int    *getCellFacePos        (const mxArray *G);    /* copy */
 int    *getCellFaces          (const mxArray *G);    /* copy */
 
 /* Cell geometry */
-void    getCellVolumes        (const mxArray *G, double **v);
-void    getCellCentroids      (const mxArray *G, double **v);
+double *getCellVolumes        (const mxArray *G);
+double *getCellCentroids      (const mxArray *G);    /* copy */
+
+
+/* Rock properties */
+double *
+getPermeability(const mxArray *perm, int d);    /* copy */
 
 #endif /* MRST_API_H_INCLUDED */