From 288062511e57fa3c6ffa72ddc1fa0686a8e32bcb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Fri, 29 Oct 2010 14:44:10 +0200
Subject: [PATCH] Implement tools for computing accumulation term.

---
 trans_tpfa.c | 36 ++++++++++++++++++++++++++++++++++++
 trans_tpfa.h |  8 ++++++++
 2 files changed, 44 insertions(+)

diff --git a/trans_tpfa.c b/trans_tpfa.c
index 6e75b4fb..e3665b78 100644
--- a/trans_tpfa.c
+++ b/trans_tpfa.c
@@ -231,3 +231,39 @@ tpfa_compr_htran(grid_t       *G         ,
     free(v);
     free(luAc);
 }
+
+
+/* ---------------------------------------------------------------------- */
+void
+tpfa_compr_accum(grid_t       *G         ,
+                 int           np        ,
+                 size_t        max_ngconn,
+                 const double *Ac        ,
+                 const double *xf        ,
+                 double       *src)
+/* ---------------------------------------------------------------------- */
+{
+    int         c, i;
+    size_t      ngconn_tot;
+    double     *luAc, *v, *qf;
+    MAT_SIZE_T *ipiv;
+
+    ngconn_tot = G->cell_facepos[ G->number_of_cells ];
+
+    qf   = malloc(ngconn_tot      * sizeof *qf);
+    luAc = malloc(np * np         * sizeof *luAc);
+    v    = malloc(np * max_ngconn * sizeof *v);
+    ipiv = malloc(np              * sizeof *ipiv);
+
+    if ((qf != NULL) && (luAc != NULL) && (v != NULL) && (ipiv != NULL)) {
+        compr_htran_core(G, np, Ac, xf, qf, luAc, v, ipiv);
+
+        for (c = i = 0; c < G->number_of_cells; c++) {
+            for (; i < G->cell_facepos[c + 1]; i++) {
+                src[c] -= qf[i];
+            }
+        }
+    }
+
+    free(ipiv);   free(v);   free(luAc);   free(qf);
+}
diff --git a/trans_tpfa.h b/trans_tpfa.h
index 1be204c2..59723da1 100644
--- a/trans_tpfa.h
+++ b/trans_tpfa.h
@@ -52,6 +52,14 @@ tpfa_compr_htran(grid_t       *G         ,
                  const double *xf        ,
                  double       *htrans);
 
+void
+tpfa_compr_accum(grid_t       *G         ,
+                 int           np        ,
+                 size_t        max_ngconn,
+                 const double *Ac        ,
+                 const double *xf        ,
+                 double       *src);
+
 #ifdef __cplusplus
 }
 #endif