Add preliminary TPFA transmissibility calculator.

Makefile.am

@@ -26,6 +26,7 @@ ifsh_ms.h                                       \
 mimetic.h                                       \
 partition.h                                     \
 sparse_sys.h                                    \
+trans_tpfa.h                                    \
 well.h                                          \
 GridAdapter.hpp                                 \
 HybridPressureSolver.hpp
@@ -47,6 +48,7 @@ ifsh_ms.c                                       \
 mimetic.c                                       \
 partition.c                                     \
 sparse_sys.c                                    \
+trans_tpfa.c                                    \
 well.c
 
 

trans_tpfa.c

@@ -0,0 +1,82 @@
+#include <assert.h>
+#include <math.h>
+
+#include "blas_lapack.h"
+#include "trans_tpfa.h"
+
+/* ---------------------------------------------------------------------- */
+/* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */
+/* ---------------------------------------------------------------------- */
+void
+tpfa_htrans_compute(grid_t *G, const double *perm, double *htrans)
+/* ---------------------------------------------------------------------- */
+{
+    int    c, d, f, i, j;
+    double s, dist, denom;
+
+    double Kn[3];
+    double *cc, *fc, *n;
+    double *K;
+
+    MAT_SIZE_T nrows, ncols, ldA, incx, incy;
+    double a1, a2;
+
+    d = G->dimensions;
+
+    nrows = ncols    = ldA = d;
+    incx  = incy     = 1      ;
+    a1    = 1.0;  a2 = 0.0    ;
+
+    for (c = i = 0; c < G->number_of_cells; c++) {
+        K  = perm + (c * d * d);
+        cc = G->cell_centroids + (c * d);
+
+        for (; i < G->cell_facepos[c + 1]; i++) {
+            f = G->cell_faces[i];
+            s = 2.0*(G->face_cells[2*f + 0] == c) - 1.0;
+
+            n  = G->face_normals   + (f * d);
+            fc = G->face_centroids + (f * d);
+
+            dgemv_("No Transpose", &nrows, &ncols,
+                   &a1, K, &ldA, n, &incx, &a2, &Kn[0], &incy);
+
+            htrans[i] = denom = 0.0;
+            for (j = 0; j < d; j++) {
+                dist = fc[j] - cc[j];
+
+                htrans[i] += s * dist * Kn[j];
+                denom     +=     dist * dist;
+            }
+
+            assert (denom > 0);
+            htrans[i] /= denom;
+            htrans[i]  = fabs(htrans[i]);
+        }
+    }
+}
+
+
+/* ---------------------------------------------------------------------- */
+void
+tpfa_trans_compute(grid_t *G, const double *htrans, double *trans)
+/* ---------------------------------------------------------------------- */
+{
+    int c, i, f;
+
+    for (f = 0; f < G->number_of_faces; f++) {
+        trans[f] = 0.0;
+    }
+
+    for (c = i = 0; c < G->number_of_cells; c++) {
+        for (; i < G->cell_facepos[c + 1]; i++) {
+            f = G->cell_faces[i];
+
+            trans[f] += 1.0 / htrans[i];
+        }
+    }
+
+    for (f = 0; f < G->number_of_faces; f++) {
+        trans[f] = 1.0 / trans[f];
+    }
+}

trans_tpfa.h

@@ -0,0 +1,31 @@
+/*
+  Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_TRANS_TPFA_HEADER_INCLUDED
+#define OPM_TRANS_TPFA_HEADER_INCLUDED
+
+#include "grid.h"
+
+void
+tpfa_htrans_compute(grid_t *G, const double *perm, double *htrans);
+
+void
+tpfa_trans_compute(grid_t *G, const double *htrans, double *trans);
+
+#endif  /* OPM_TRANS_TPFA_HEADER_INCLUDED */