diff --git a/Makefile.am b/Makefile.am
index d10c0bad..a1f38153 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -74,6 +74,7 @@ opm/core/pressure/msmfem/ifsh_ms.c \
opm/core/pressure/msmfem/dfs.c \
opm/core/pressure/msmfem/coarse_sys.c \
opm/core/pressure/ifsh.c \
+opm/core/pressure/IncompTpfa.cpp \
opm/core/pressure/mimetic/mimetic.c \
opm/core/pressure/mimetic/hybsys_global.c \
opm/core/pressure/mimetic/hybsys.c \
@@ -157,6 +158,7 @@ opm/core/GridAdapter.hpp \
opm/core/GridManager.hpp \
opm/core/pressure/fsh.h \
opm/core/pressure/HybridPressureSolver.hpp \
+opm/core/pressure/IncompTpfa.hpp \
opm/core/pressure/TPFAPressureSolver.hpp \
opm/core/pressure/tpfa/compr_quant.h \
opm/core/pressure/tpfa/compr_source.h \
diff --git a/opm/core/pressure/IncompTpfa.cpp b/opm/core/pressure/IncompTpfa.cpp
new file mode 100644
index 00000000..10a2b7ab
--- /dev/null
+++ b/opm/core/pressure/IncompTpfa.cpp
@@ -0,0 +1,112 @@
+/*
+ Copyright 2012 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 .
+*/
+
+#include
+#include
+#include
+#include
+#include
+
+
+namespace Opm
+{
+
+
+
+ /// Construct solver.
+ /// \param[in] g A 2d or 3d grid.
+ /// \param[in] permeability Array of permeability tensors, the array
+ /// should have size N*D^2, if D == g.dimensions
+ /// and N == g.number_of_cells.
+ /// \param[in] gravity Gravity vector. If nonzero, the array should
+ /// have D elements.
+ IncompTpfa::IncompTpfa(const UnstructuredGrid& g,
+ const double* permeability,
+ const double* gravity)
+ : grid_(g),
+ htrans_(g.cell_facepos[ g.number_of_cells ]),
+ trans_ (g.number_of_faces),
+ gpress_(g.cell_facepos[ g.number_of_cells ], 0.0),
+ gpress_omegaweighted_(g.cell_facepos[ g.number_of_cells ], 0.0)
+ {
+ UnstructuredGrid* gg = const_cast(&grid_);
+ tpfa_htrans_compute(gg, permeability, &htrans_[0]);
+ if (gravity) {
+ mim_ip_compute_gpress(gg->number_of_cells, gg->dimensions, gravity,
+ gg->cell_facepos, gg->cell_faces,
+ gg->face_centroids, gg->cell_centroids,
+ &gpress_[0]);
+ }
+ h_ = ifs_tpfa_construct(gg);
+ }
+
+
+
+
+ /// Destructor.
+ IncompTpfa::~IncompTpfa()
+ {
+ ifs_tpfa_destroy(h_);
+ }
+
+
+
+ /// Assemble and solve pressure system.
+ /// \param[in] totmob Must contain N total mobility values (one per cell).
+ /// totmob = \sum_{p} kr_p/mu_p.
+ /// \param[in] omega Must contain N mobility-weighted density values (one per cell).
+ /// omega = \frac{\sum_{p} mob_p rho_p}{\sum_p rho_p}.
+ /// \param[in] src Must contain N souce rates (one per cell).
+ /// Positive values represent total inflow rates,
+ /// negative values represent total outflow rates.
+ /// \param[out] pressure Will contain N cell-pressure values.
+ /// \param[out] faceflux Will contain F signed face flux values.
+ void IncompTpfa::solve(const std::vector& totmob,
+ const std::vector& omega,
+ const std::vector& src,
+ std::vector& pressure,
+ std::vector& faceflux)
+ {
+ UnstructuredGrid* gg = const_cast(&grid_);
+ tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
+
+ if (!omega.empty()) {
+ mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
+ &omega[0],
+ &gpress_[0], &gpress_omegaweighted_[0]);
+ }
+
+ ifs_tpfa_assemble(gg, &trans_[0], &src[0], &gpress_omegaweighted_[0], h_);
+
+ using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
+ CSRMatrixUmfpackSolver linsolve;
+ linsolve.solve(h_->A, h_->b, h_->x);
+
+ pressure.resize(grid_.number_of_cells);
+ faceflux.resize(grid_.number_of_faces);
+
+ ifs_tpfa_press_flux(gg, &trans_[0], h_,
+ &pressure[0],
+ &faceflux[0]);
+ }
+
+
+
+
+} // namespace Opm
diff --git a/opm/core/pressure/IncompTpfa.hpp b/opm/core/pressure/IncompTpfa.hpp
new file mode 100644
index 00000000..ce12b7d0
--- /dev/null
+++ b/opm/core/pressure/IncompTpfa.hpp
@@ -0,0 +1,83 @@
+/*
+ Copyright 2012 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 .
+*/
+
+#ifndef OPM_INCOMPTPFA_HEADER_INCLUDED
+#define OPM_INCOMPTPFA_HEADER_INCLUDED
+
+
+#include
+
+struct UnstructuredGrid;
+struct ifs_tpfa_data;
+
+
+namespace Opm
+{
+
+ /// Encapsulating a tpfa pressure solver for the incompressible case.
+ /// Supports gravity and simple sources as driving forces.
+ /// Below we use the shortcuts D for the number of dimensions, N
+ /// for the number of cells and F for the number of faces.
+ /// Note: we intend to add wells in the future.
+ class IncompTpfa
+ {
+ public:
+ /// Construct solver.
+ /// \param[in] g A 2d or 3d grid.
+ /// \param[in] permeability Array of permeability tensors, the array
+ /// should have size N*D^2, if D == g.dimensions
+ /// and N == g.number_of_cells.
+ /// \param[in] gravity Gravity vector. If nonzero, the array should
+ /// have D elements.
+ IncompTpfa(const UnstructuredGrid& g,
+ const double* permeability,
+ const double* gravity);
+
+ /// Destructor.
+ ~IncompTpfa();
+
+ /// Assemble and solve pressure system.
+ /// \param[in] totmob Must contain N total mobility values (one per cell).
+ /// totmob = \sum_{p} kr_p/mu_p.
+ /// \param[in] omega Must contain N mobility-weighted density values (one per cell).
+ /// omega = \frac{\sum_{p} mob_p rho_p}{\sum_p rho_p}.
+ /// \param[in] src Must contain N souce rates (one per cell).
+ /// Positive values represent total inflow rates,
+ /// negative values represent total outflow rates.
+ /// \param[out] pressure Will contain N cell-pressure values.
+ /// \param[out] faceflux Will contain F signed face flux values.
+ void solve(const std::vector& totmob,
+ const std::vector& omega,
+ const std::vector& src,
+ std::vector& pressure,
+ std::vector& faceflux);
+
+ private:
+ const UnstructuredGrid& grid_;
+ ::std::vector htrans_;
+ ::std::vector trans_ ;
+ ::std::vector gpress_;
+ ::std::vector gpress_omegaweighted_;
+
+ struct ifs_tpfa_data* h_;
+ };
+
+} // namespace Opm
+
+#endif // OPM_INCOMPTPFA_HEADER_INCLUDED