From ced66b0c8040bb7c6ae27e03b33b400a07531152 Mon Sep 17 00:00:00 2001
From: Markus Wolff <markus.wolff@twt-gmbh.de>
Date: Fri, 6 Sep 2013 15:00:24 +0000
Subject: [PATCH] replaced mimetic finite difference implementation of the 2p
pressure equation by the newer implementations from dumux-devel
- includes gravity terms but no capillary pressure terms!
- additional implementation for adaptive grids
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11406 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../{croperator.hh => croperator2p.hh} | 66 +-
.../diffusion/mimetic/croperator2padaptive.hh | 564 +++++++++++++
.../2p/diffusion/mimetic/mimetic2p.hh | 747 +++++++++++++++++
.../2p/diffusion/mimetic/mimetic2padaptive.hh | 775 ++++++++++++++++++
.../diffusion/mimetic/mimeticgroundwater.hh | 436 ----------
.../2p/diffusion/mimetic/mimeticoperator.hh | 126 ---
.../2p/diffusion/mimetic/mimeticoperator2p.hh | 300 +++++++
.../mimetic/mimeticoperator2padaptive.hh | 308 +++++++
.../2p/diffusion/mimetic/mimeticpressure2p.hh | 481 ++++++-----
.../mimetic/mimeticpressure2padaptive.hh | 481 +++++++++++
.../mimetic/mimeticpressureproperties2p.hh | 23 +-
.../mimeticpressureproperties2padaptive.hh | 79 ++
12 files changed, 3589 insertions(+), 797 deletions(-)
rename dumux/decoupled/2p/diffusion/mimetic/{croperator.hh => croperator2p.hh} (85%)
create mode 100644 dumux/decoupled/2p/diffusion/mimetic/croperator2padaptive.hh
create mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimetic2p.hh
create mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimetic2padaptive.hh
delete mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
delete mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimeticoperator.hh
create mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2p.hh
create mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2padaptive.hh
create mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
create mode 100644 dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2padaptive.hh
diff --git a/dumux/decoupled/2p/diffusion/mimetic/croperator.hh b/dumux/decoupled/2p/diffusion/mimetic/croperator2p.hh
similarity index 85%
rename from dumux/decoupled/2p/diffusion/mimetic/croperator.hh
rename to dumux/decoupled/2p/diffusion/mimetic/croperator2p.hh
index 8325fd1567..6b07041189 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/croperator.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/croperator2p.hh
@@ -16,8 +16,8 @@
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
-#ifndef DUMUX_CROPERATOR_HH
-#define DUMUX_CROPERATOR_HH
+#ifndef DUMUX_CROPERATOR2P_HH
+#define DUMUX_CROPERATOR2P_HH
#include<iostream>
#include<vector>
@@ -25,6 +25,7 @@
#include<map>
#include<cassert>
+#include<dune/common/timer.hh>
#include<dune/common/fvector.hh>
#include<dune/common/fmatrix.hh>
#include<dune/common/exceptions.hh>
@@ -46,7 +47,7 @@
namespace Dumux
{
/*!
- * \ingroup MimeticPressure2p
+ * \ingroup Mimetic2p
*/
/**
* @brief defines a class for Crozieux-Raviart piecewise linear finite element functions
@@ -68,8 +69,8 @@ namespace Dumux
* \tparam Scalar The field type used in the elements of the global stiffness matrix
* \tparam GridView The grid view of the simulation grid
*/
-template<class Scalar, class GridView>
-class CROperatorAssembler
+template<class TypeTag>
+class CROperatorAssemblerTwoP
{
template<int dim>
struct FaceLayout
@@ -79,13 +80,11 @@ class CROperatorAssembler
return gt.dim() == dim-1;
}
};
-
- typedef typename GridView::Grid Grid;
- enum {dim=Grid::dimension};
- typedef typename Grid::template Codim<0>::Entity Entity;
- typedef typename GridView::template Codim<0>::Iterator Iterator;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ enum {dim=GridView::dimension};
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IndexSet IS;
- typedef typename Grid::template Codim<0>::EntityPointer EEntityPointer;
typedef Dune::FieldMatrix<Scalar,1,1> BlockType;
typedef Dune::BCRSMatrix<BlockType> MatrixType;
typedef typename MatrixType::block_type MBlockType;
@@ -95,9 +94,10 @@ class CROperatorAssembler
typedef Dune::BlockVector< Dune::FieldVector<double,1> > SatType;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,FaceLayout> FaceMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum
{
- pressEqIdx = 0
+ pressEqIdx = Indices::pressEqIdx,
};
//! a function to approximately compute the number of nonzeros
@@ -109,13 +109,12 @@ class CROperatorAssembler
public:
typedef MatrixType RepresentationType;
- CROperatorAssembler (const GridView& gridview)
+ CROperatorAssemblerTwoP (const GridView& gridview)
: gridView_(gridview), faceMapper_(gridView_), size_(faceMapper_.size()),
A_(size_, size_, nnz(), RepresentationType::random)
{}
- //!Initialize the global matrix of the system of equations to solve
- void initializeMatrix()
+ void initialize()
{
faceMapper_.update();
@@ -132,8 +131,8 @@ public:
std::vector<bool> visited(size_, false);
// LOOP 1 : Compute row sizes
- Iterator eEndIt = gridView_.template end<0>();
- for (Iterator eIt = gridView_.template begin<0>(); eIt != eEndIt; ++eIt)
+ ElementIterator eEndIt = gridView_.template end<0>();
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt != eEndIt; ++eIt)
{
int numFaces = eIt->template count<1>();
@@ -160,7 +159,7 @@ public:
visited.assign(size_, false);
// LOOP 2 : insert the nonzeros
- for (Iterator eIt = gridView_.template begin<0>(); eIt!=eEndIt; ++eIt)
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt!=eEndIt; ++eIt)
{
int numFaces = eIt->template count<1>();
@@ -197,6 +196,16 @@ public:
return A_;
}
+ const FaceMapper& faceMapper()
+ {
+ return faceMapper_;
+ }
+
+ const FaceMapper& faceMapper() const
+ {
+ return faceMapper_;
+ }
+
/*! @brief Assemble global stiffness matrix
This method takes an object that can compute local stiffness matrices and
@@ -223,7 +232,7 @@ public:
// check size
if (u.N()!=A_.M() || f.N()!=A_.N())
- DUNE_THROW(Dune::MathError,"CROperatorAssembler::assemble(): size mismatch");
+ DUNE_THROW(Dune::MathError,"CROperatorAssemblerTwoP::assemble(): size mismatch");
// clear global stiffness matrix and right hand side
A_ = 0;
f = 0;
@@ -234,11 +243,11 @@ public:
essential[i][0] = BoundaryConditions::neumann;
// local to global id mapping (do not ask vertex mapper repeatedly
- int local2Global[2*GridView::dimension];
+ Dune::FieldVector<int, 2*dim> local2Global(0);
// run over all leaf elements
- Iterator eEndIt = gridView_.template end<0>();
- for (Iterator eIt = gridView_.template begin<0>(); eIt!=eEndIt; ++eIt)
+ ElementIterator eEndIt = gridView_.template end<0>();
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt!=eEndIt; ++eIt)
{
unsigned int numFaces = eIt->template count<1>();
@@ -274,6 +283,19 @@ public:
f[local2Global[i]][0] += loc.rhs(i)[0];
}
}
+ // run over all leaf elements
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt!=eEndIt; ++eIt)
+ {
+ unsigned int numFaces = eIt->template count<1>();
+
+ // get local to global id map
+ for (int k = 0; k < numFaces; k++)
+ {
+ int alpha = faceMapper_.map(*eIt, k, 1);
+ local2Global[k] = alpha;
+ }
+ loc.completeRHS(*eIt, local2Global, f);
+ }
// put in essential boundary conditions
rowiterator endi=A_.end();
diff --git a/dumux/decoupled/2p/diffusion/mimetic/croperator2padaptive.hh b/dumux/decoupled/2p/diffusion/mimetic/croperator2padaptive.hh
new file mode 100644
index 0000000000..f8cb84ee07
--- /dev/null
+++ b/dumux/decoupled/2p/diffusion/mimetic/croperator2padaptive.hh
@@ -0,0 +1,564 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_CROPERATOR2PADAPTIVE_HH
+#define DUMUX_CROPERATOR2PADAPTIVE_HH
+
+#include<iostream>
+#include<vector>
+#include<set>
+#include<map>
+#include<cassert>
+#include<stdio.h>
+#include<stdlib.h>
+
+#include<dune/common/timer.hh>
+#include<dune/common/fvector.hh>
+#include<dune/common/fmatrix.hh>
+#include<dune/common/exceptions.hh>
+#include<dune/geometry/type.hh>
+#include<dune/grid/common/grid.hh>
+#include<dune/grid/common/mcmgmapper.hh>
+
+#include<dune/istl/bvector.hh>
+#include<dune/istl/operators.hh>
+#include<dune/istl/bcrsmatrix.hh>
+#include<dumux/common/boundaryconditions.hh>
+#include<dumux/decoupled/2p/diffusion/mimetic/localstiffness.hh>
+
+/**
+ * @file
+ * @brief defines a class for piecewise linear finite element functions
+ */
+
+namespace Dumux
+{
+
+//forward declaration
+template<class GridView>
+class IntersectionMapper;
+
+/*!
+ * \ingroup Mimetic2p
+ */
+/**
+ * @brief defines a class for Crozieux-Raviart piecewise linear finite element functions
+ *
+ */
+
+/*! @brief A class for mapping a CR function to a CR function
+
+ This class sets up a compressed row storage matrix with connectivity for CR elements.
+
+ This class does not fill any entries into the matrix.
+
+ The template parameter TypeTag describes what kind of Assembler we are. There two choices:
+ <dt>LevelTag</dt> We assemble on a grid level.
+ <dt>LeafTag</dt> We assemble on the leaf entities of the grid
+*/
+/*! @brief Extends CROperatorBase by a generic methods to assemble global stiffness matrix from local stiffness matrices
+ *
+ *
+ * The template parameter TypeTag describes what kind of Assembler we are. There two choices:
+ * <dt>LevelTag</dt> We assemble on a grid level.
+ * <dt>LeafTag</dt> We assemble on the leaf entities of the grid
+ */
+template<class TypeTag>
+class CROperatorAssemblerTwoPAdaptive
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ enum {dim=GridView::dimension};
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::IndexSet IS;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef Dune::FieldMatrix<Scalar,1,1> BlockType;
+ typedef Dune::BCRSMatrix<BlockType> MatrixType;
+ typedef typename MatrixType::block_type MBlockType;
+ typedef typename MatrixType::RowIterator rowiterator;
+ typedef typename MatrixType::ColIterator coliterator;
+ typedef Dune::array<BoundaryConditions::Flags,1> BCBlockType; // componentwise boundary conditions
+ typedef Dune::BlockVector< Dune::FieldVector<double,1> > SatType;
+ typedef Dumux::IntersectionMapper<GridView> IntersectionMapper;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum
+ {
+ pressEqIdx = Indices::pressEqIdx,
+ };
+
+ // return number of rows/columns
+ int size () const
+ {
+ return intersectionMapper_.size();
+ }
+
+public:
+ typedef MatrixType RepresentationType;
+
+ CROperatorAssemblerTwoPAdaptive (const GridView& gridview)
+ : gridView_(gridview), is_(gridView_.indexSet()), intersectionMapper_(gridView_)
+ {
+ A_.setBuildMode(MatrixType::random);
+ }
+
+ void initialize()
+ {
+ adapt();
+ }
+
+ void adapt()
+ {
+ intersectionMapper_.update();
+ A_.setSize(size(), size());
+ updateMatrix();
+ }
+
+ //! return const reference to operator matrix
+ const RepresentationType& operator* () const
+ {
+ return A_;
+ }
+
+ //! return reference to operator matrix
+ RepresentationType& operator* ()
+ {
+ return A_;
+ }
+
+ const IntersectionMapper& intersectionMapper()
+ {
+ return intersectionMapper_;
+ }
+
+ const IntersectionMapper& intersectionMapper() const
+ {
+ return intersectionMapper_;
+ }
+
+ /*! @brief Assemble global stiffness matrix
+
+ This method takes an object that can compute local stiffness matrices and
+ assembles the global linear system Au=f.
+
+ @param[in] loc the local assembler providing element stiffness and boundary conditions for all elements
+ @param[in,out] u solution, contains initial values on input, Dirichlet values are set. The
+ type of boundary condition for a node is inferred from the values returned
+ by the local assembler. A node is of Neumann type if all elements referring
+ to that node report a Neumann boundary condition, it is set to Dirichlet
+ if a least one element reports a process or Dirichlet boundary condition. The difference
+ between process and Dirichlet is that process always denotes a homogeneous Dirichlet
+ value.
+ @param[in] f right hand side is filled by this method
+
+ Note that the rows corresponding to nodes at the Dirichlet boundary are filled
+ with trivial equations of the form \f[1\cdot u_i = f_i \f] where \f$u_i\f$ and \f$f_i\f$ are both set to the
+ Dirichlet value at the \f$i\f$th node.
+
+ */
+ template <class LocalStiffness, class Vector>
+ void assemble (LocalStiffness& loc, Vector& u, Vector& f);
+
+ void updateMatrix();
+
+protected:
+ const GridView& gridView_;
+ const IS& is_;
+ IntersectionMapper intersectionMapper_;
+ RepresentationType A_;
+};
+
+template<class TypeTag>
+void CROperatorAssemblerTwoPAdaptive<TypeTag>::updateMatrix()
+{
+ // set size of all rows to zero
+ for (unsigned int i = 0; i < size(); i++)
+ A_.setrowsize(i, 0);
+
+ // build needs a flag for all entities of all codims
+ std::vector<bool> visited(size(), false);
+
+ int numElem = gridView_.size(0);
+
+ for (int elemIdx = 0; elemIdx < numElem; elemIdx++)
+ {
+ int numFaces = intersectionMapper_.size(elemIdx);
+ for (int faceIdx = 0; faceIdx < numFaces; faceIdx++)
+ {
+ int faceIdxGlobal = intersectionMapper_.map(elemIdx, faceIdx);
+ if (!visited[faceIdxGlobal])
+ {
+ A_.incrementrowsize(faceIdxGlobal);
+ visited[faceIdxGlobal] = true;
+ }
+ for (int k = 0; k < numFaces-1; k++) {
+ A_.incrementrowsize(faceIdxGlobal);
+ }
+
+ }
+
+ }
+
+ A_.endrowsizes();
+
+ // clear the flags for the next round, actually that is not necessary because addindex takes care of this
+ for (int i = 0; i < size(); i++)
+ visited[i] = false;
+
+ for (int elemIdx = 0; elemIdx < numElem; elemIdx++)
+ {
+ int numFaces = intersectionMapper_.size(elemIdx);
+ for (int faceIdx = 0; faceIdx < numFaces; faceIdx++)
+ {
+ int faceIdxGlobalI = intersectionMapper_.map(elemIdx, faceIdx);
+ if (!visited[faceIdxGlobalI])
+ {
+ A_.addindex(faceIdxGlobalI,faceIdxGlobalI);
+ visited[faceIdxGlobalI] = true;
+ }
+ for (int k = 0; k < numFaces; k++)
+ if (k != faceIdx) {
+ int faceIdxGlobalJ = intersectionMapper_.map(elemIdx, k);
+ A_.addindex(faceIdxGlobalI, faceIdxGlobalJ);
+ }
+
+ }
+
+ }
+
+ A_.endindices();
+}
+
+template<class TypeTag>
+template <class LocalStiffness, class Vector>
+void CROperatorAssemblerTwoPAdaptive<TypeTag>::assemble(LocalStiffness& loc, Vector& u, Vector& f)
+{
+
+ // check size
+ if (u.N()!=A_.M() || f.N()!=A_.N())
+ DUNE_THROW(Dune::MathError,"CROperatorAssemblerTwoPAdaptive::assemble(): size mismatch");
+ // clear global stiffness matrix and right hand side
+ A_ = 0;
+ f = 0;
+
+ // allocate flag vector to hold flags for essential boundary conditions
+ std::vector<BCBlockType> essential(intersectionMapper_.size());
+ for (typename std::vector<BCBlockType>::size_type i=0; i<essential.size(); i++)
+ essential[i][0] = BoundaryConditions::neumann;
+
+ // local to global id mapping (do not ask vertex mapper repeatedly
+ std::vector<int> local2Global(2*dim);
+
+ // run over all leaf elements
+ ElementIterator eendit = gridView_.template end<0>();
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt!=eendit; ++eIt)
+ {
+ // build local stiffness matrix for CR elements
+ // inludes rhs and boundary condition information
+ loc.assemble(*eIt, 1); // assemble local stiffness matrix
+
+ int globalIdx = intersectionMapper_.map(*eIt);
+
+ unsigned int numFaces = intersectionMapper_.size(globalIdx);
+ local2Global.resize(numFaces);
+
+ for (int i = 0; i < numFaces; i++)
+ {
+ int idx = intersectionMapper_.map(*eIt, i);
+ local2Global[i] = idx;
+ }
+
+ // accumulate local matrix into global matrix for non-hanging nodes
+ for (int i=0; i<numFaces; i++) // loop over rows, i.e. test functions
+ {
+ // accumulate matrix
+ for (int j=0; j<numFaces; j++)
+ {
+ // the standard entry
+ A_[local2Global[i]][local2Global[j]] += loc.mat(i,j);
+ }
+
+ // essential boundary condition and rhs
+ if (loc.bc(i).isDirichlet(pressEqIdx))
+ {
+ essential[local2Global[i]][0] = BoundaryConditions::dirichlet;
+ f[local2Global[i]][0] = loc.rhs(i)[0];
+ }
+ else
+ f[local2Global[i]][0] += loc.rhs(i)[0];
+ }
+
+ }
+ // run over all leaf elements
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt!=eendit; ++eIt)
+ {
+ int globalIdx = intersectionMapper_.map(*eIt);
+
+ unsigned int numFaces = intersectionMapper_.size(globalIdx);
+ local2Global.resize(numFaces);
+
+ for (int i = 0; i < numFaces; i++)
+ {
+ int idx = intersectionMapper_.map(*eIt, i);
+ local2Global[i] = idx;
+ }
+
+ loc.completeRHS(*eIt, local2Global, f);
+ }
+
+ // put in essential boundary conditions
+ rowiterator endi=A_.end();
+ for (rowiterator i=A_.begin(); i!=endi; ++i)
+ {
+ // muck up extra rows
+ if ((int) i.index() >= (int) size())
+ {
+ coliterator endj=(*i).end();
+ for (coliterator j=(*i).begin(); j!=endj; ++j)
+ {
+ (*j) = 0;
+ if (j.index()==i.index())
+ (*j)[0][0] = 1;
+ }
+ f[i.index()] = 0;
+ continue;
+ }
+
+ // insert dirichlet ans processor boundary conditions
+ if (essential[i.index()][0]!=BoundaryConditions::neumann)
+ {
+ coliterator endj=(*i).end();
+ for (coliterator j=(*i).begin(); j!=endj; ++j)
+ if (j.index()==i.index())
+ {
+ (*j)[0][0] = 1;
+ }
+ else
+ {
+ (*j)[0][0] = 0;
+ }
+ u[i.index()][0] = f[i.index()][0];
+ }
+ }
+}
+
+template<class GridView>
+class IntersectionMapper
+{
+ // mapper: one data element in every entity
+ template<int dim>
+ struct ElementLayout
+ {
+ bool contains (Dune::GeometryType gt)
+ {
+ return gt.dim() == dim;
+ }
+ };
+
+ typedef typename GridView::Grid Grid;
+ enum {dim=Grid::dimension};
+ typedef typename Grid::template Codim<0>::Entity Element;
+ typedef typename Grid::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::Intersection Intersection;
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,ElementLayout> ElementMapper;
+
+
+
+public:
+ IntersectionMapper(const GridView& gridview)
+ : gridView_(gridview), elementMapper_(gridView_), size_(gridView_.size(1)), intersectionMapGlobal_(gridView_.size(0)), intersectionMapLocal_(gridView_.size(0))
+ {
+ ElementIterator eIt = gridView_.template begin<0>();
+
+ int faceIdx = 0;
+
+ IntersectionIterator isItEnd = gridView_.iend(*eIt);
+ for (IntersectionIterator isIt = gridView_.ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ int idxInInside = isIt->indexInInside();
+
+ standardLocalIdxMap_[idxInInside] = faceIdx;
+
+ faceIdx++;
+ }
+ }
+
+ const ElementMapper& elementMapper() const
+ {
+ return elementMapper_;
+ }
+
+ int map(const Element& element) const
+ {
+ return elementMapper_.map(element);
+ }
+
+ int map(int elemIdx, int faceIdx)
+ {
+ return intersectionMapGlobal_[elemIdx][faceIdx];
+ }
+
+ int map(int elemIdx, int faceIdx) const
+ {
+ return (intersectionMapGlobal_[elemIdx].find(faceIdx))->second;//use find() for const function!
+ }
+
+ int map(const Element& element, int faceIdx)
+ {
+ return intersectionMapGlobal_[map(element)][faceIdx];
+ }
+
+ int map(const Element& element, int faceIdx) const
+ {
+ return intersectionMapGlobal_[map(element)].find(faceIdx)->second;//use find() for const function!
+ }
+
+ int maplocal(int elemIdx, int faceIdx)
+ {
+ return intersectionMapLocal_[elemIdx][faceIdx];
+ }
+
+ int maplocal(int elemIdx, int faceIdx) const
+ {
+ return (intersectionMapLocal_[elemIdx].find(faceIdx))->second;//use find() for const function!
+ }
+
+
+ int maplocal(const Element& element, int faceIdx)
+ {
+ return intersectionMapLocal_[map(element)][faceIdx];
+ }
+
+ int maplocal(const Element& element, int faceIdx) const
+ {
+ return (intersectionMapLocal_[map(element)].find(faceIdx))->second;//use find() for const function!
+ }
+
+ // return number intersections
+ unsigned int size () const
+ {
+ return size_;
+ }
+
+ unsigned int size (int elemIdx) const
+ {
+ return intersectionMapLocal_[elemIdx].size();
+ }
+
+ unsigned int size (const Element& element) const
+ {
+ return intersectionMapLocal_[map(element)].size();
+ }
+
+ void update()
+ {
+ elementMapper_.update();
+
+ intersectionMapGlobal_.clear();
+ intersectionMapGlobal_.resize(elementMapper_.size());
+ intersectionMapLocal_.clear();
+ intersectionMapLocal_.resize(elementMapper_.size());
+
+ ElementIterator eItEnd = gridView_.template end<0>();
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = map(*eIt);
+
+ int faceIdx = 0;
+ // run through all intersections with neighbors
+ IntersectionIterator isItEnd = gridView_.iend(*eIt);
+ for (IntersectionIterator isIt = gridView_.ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ int indexInInside = isIt->indexInInside();
+ intersectionMapLocal_[globalIdx][faceIdx] = indexInInside;
+
+ faceIdx++;
+ }
+ }
+
+ int globalIntersectionIdx = 0;
+ for (ElementIterator eIt = gridView_.template begin<0>(); eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = map(*eIt);
+
+ int faceIdx = 0;
+ // run through all intersections with neighbors
+ IntersectionIterator isItEnd = gridView_.iend(*eIt);
+ for (IntersectionIterator isIt = gridView_.ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ if (isIt->neighbor())
+ {
+ ElementPointer neighbor = isIt->outside();
+ int globalIdxNeighbor = map(*neighbor);
+
+ if (eIt->level() > neighbor->level() || (eIt->level() == neighbor->level() && globalIdx < globalIdxNeighbor))
+ {
+
+ int faceIdxNeighbor = 0;
+ if (size(globalIdxNeighbor) == 2 * dim)
+ {
+ faceIdxNeighbor = standardLocalIdxMap_[isIt->indexInOutside()];
+ }
+ else
+ {
+ IntersectionIterator isItNEnd = gridView_.iend(*neighbor);
+ for (IntersectionIterator isItN = gridView_.ibegin(*neighbor); isItN != isItNEnd; ++isItN)
+ {
+ if (isItN->neighbor())
+ {
+ if (isItN->outside() == eIt)
+ {
+ break;
+ }
+ }
+ faceIdxNeighbor++;
+ }
+ }
+
+ intersectionMapGlobal_[globalIdx][faceIdx] = globalIntersectionIdx;
+ intersectionMapGlobal_[globalIdxNeighbor][faceIdxNeighbor] = globalIntersectionIdx;
+ globalIntersectionIdx ++;
+ }
+ }
+ else
+ {
+ intersectionMapGlobal_[globalIdx][faceIdx] = globalIntersectionIdx;
+ globalIntersectionIdx ++;
+ }
+ faceIdx++;
+ }
+ }
+ size_ = globalIntersectionIdx;
+ }
+
+protected:
+ const GridView& gridView_;
+ ElementMapper elementMapper_;
+ unsigned int size_;
+ std::vector<std::map<int, int> > intersectionMapGlobal_;
+ std::vector<std::map<int, int> > intersectionMapLocal_;
+ std::map<int, int> standardLocalIdxMap_;
+};
+
+}
+
+#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimetic2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimetic2p.hh
new file mode 100644
index 0000000000..c0e936867e
--- /dev/null
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimetic2p.hh
@@ -0,0 +1,747 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Markus Wolff *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Local stiffness matrix for the diffusion equation discretized by mimetic FD
+ */
+#ifndef DUMUX_MIMETIC2P_HH
+#define DUMUX_MIMETIC2P_HH
+
+#include<map>
+#include<iostream>
+#include<iomanip>
+#include<fstream>
+#include<vector>
+#include<sstream>
+
+#include<dune/common/exceptions.hh>
+#include<dune/grid/common/grid.hh>
+#include<dune/geometry/type.hh>
+#include<dune/geometry/quadraturerules.hh>
+
+#include<dumux/common/boundaryconditions.hh>
+
+#include <dune/common/dynvector.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mimetic2p
+ */
+/**
+ * @brief compute local stiffness matrix for conforming finite elements for the full 2-phase pressure equation
+ *
+ */
+
+//! A class for computing local stiffness matrices
+/*! A class for computing local stiffness matrix for the full 2-phase pressure equation
+ */
+template<class TypeTag>
+class MimeticTwoPLocalStiffness: public LocalStiffness<TypeTag, 1>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ // grid types
+ enum
+ {
+ dim = GridView::dimension
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+ pressEqIdx = Indices::pressEqIdx,
+ satEqIdx = Indices::satEqIdx,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation), //!< gives kind of pressure used (\f$ 0 = p_w\f$, \f$ 1 = p_n\f$, \f$ 2 = p_{global}\f$)
+ saturationType = GET_PROP_VALUE(TypeTag, SaturationFormulation) //!< gives kind of saturation used (\f$ 0 = S_w\f$, \f$ 1 = S_n\f$)
+ };
+
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+
+public:
+ // define the number of components of your system, this is used outside
+ // to allocate the correct size of (dense) blocks with a FieldMatrix
+ enum
+ {
+ m = 1
+ };
+ enum
+ {
+ size = 2 * dim
+ };
+
+ //! Constructor
+ MimeticTwoPLocalStiffness(Problem& problem, bool levelBoundaryAsDirichlet,
+ const GridView& gridView, bool procBoundaryAsDirichlet = true) :
+ problem_(problem), gridView_(gridView), maxError_(0), timeStep_(1)
+ {
+ ErrorTermFactor_ = GET_PARAM(TypeTag, Scalar, ImpetErrorTermFactor);
+ ErrorTermLowerBound_ = GET_PARAM(TypeTag, Scalar, ImpetErrorTermLowerBound);
+ ErrorTermUpperBound_ = GET_PARAM(TypeTag, Scalar, ImpetErrorTermUpperBound);
+
+ density_[wPhaseIdx] = 0.0;
+ density_[nPhaseIdx] = 0.0;
+ viscosity_[wPhaseIdx] =0.0;
+ viscosity_[nPhaseIdx] = 0.0;
+ }
+
+ void initialize()
+ {
+ ElementIterator element = problem_.gridView().template begin<0>();
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setPressure(nPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setTemperature(problem_.temperature(*element));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
+ density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
+ int size = gridView_.size(0);
+ rhs_.resize(size , 0.);
+ W_.resize(size);
+ }
+
+ void reset()
+ {
+ rhs_ = 0;
+ maxError_ = 0;
+ timeStep_ = 1;
+ }
+
+ void setErrorInfo(Scalar maxErr, Scalar dt)
+ {
+ maxError_ = maxErr;
+ timeStep_ = dt;
+ }
+
+ //! assemble local stiffness matrix for given element and order
+ /*! On exit the following things have been done:
+ - The stiffness matrix for the given entity and polynomial degree has been assembled and is
+ accessible with the mat() method.
+ - The boundary conditions have been evaluated and are accessible with the bc() method
+ - The right hand side has been assembled. It contains either the value of the essential boundary
+ condition or the assembled source term and neumann boundary condition. It is accessible via the rhs() method.
+ @param[in] element a codim 0 entity reference
+ @param[in] k order of CR basis (only k = 1 is implemented)
+ */
+ void assemble(const Element& element, int k = 1)
+ {
+ unsigned int numFaces = element.template count<1>();
+ this->setcurrentsize(numFaces);
+
+ // clear assemble data
+ for (int i = 0; i < numFaces; i++)
+ {
+ this->b[i] = 0;
+ this->bctype[i].reset();
+ for (int j = 0; j < numFaces; j++)
+ this->A[i][j] = 0;
+ }
+
+ assembleV(element, k);
+ assembleBC(element, k);
+ }
+
+ // TODO/FIXME: this is only valid for linear problems where
+ // the local stiffness matrix is independend of the current
+ // solution. We need to implement this properly, but this
+ // should at least make the thing compile...
+ typedef Dune::FieldVector<Scalar, m> VBlockType;
+ void assemble(const Element &cell, const Dune::BlockVector<VBlockType>& localSolution, int orderOfShapeFns = 1)
+ {
+ assemble(cell, orderOfShapeFns);
+ }
+
+ //! assemble only boundary conditions for given element
+ /*! On exit the following things have been done:
+ - The boundary conditions have been evaluated and are accessible with the bc() method
+ - The right hand side contains either the value of the essential boundary
+ condition or the assembled neumann boundary condition. It is accessible via the rhs() method.
+ @param[in] element a codim 0 entity reference
+ @param[in] k order of CR basis
+ */
+ void assembleBoundaryCondition(const Element& element, int k = 1)
+ {
+ unsigned int numFaces = element.template count<1>();
+ this->setcurrentsize(numFaces);
+
+ // clear assemble data
+ for (int i = 0; i < numFaces; i++)
+ {
+ this->b[i] = 0;
+ this->bctype[i].reset();
+ }
+
+ assembleBC(element, k);
+ }
+
+ template<class Vector>
+ void completeRHS(const Element& element, Dune::FieldVector<int, 2*dim>& local2Global, Vector& f)
+ {
+ int globalIdx = problem_.variables().index(element);
+ unsigned int numFaces = element.template count<1>();
+
+ Dune::FieldVector<Scalar, 2 * dim> F(0.);
+ Scalar dInv = 0.;
+ computeReconstructionMatrices(element, W_[globalIdx], F, dInv);
+
+ for (int i = 0; i < numFaces; i++)
+ {
+ if (!this->bc(i).isDirichlet(pressEqIdx))
+ f[local2Global[i]][0] += (dInv * F[i] * rhs_[globalIdx]);
+ }
+ }
+
+ Scalar constructPressure(const Element& element, Dune::FieldVector<Scalar,2*dim>& pressTrace)
+ {
+ int globalIdx = problem_.variables().index(element);
+ unsigned int numFaces = element.template count<1>();
+ Scalar volume = element.geometry().volume();
+
+ PrimaryVariables sourceVec(0.0);
+ problem_.source(sourceVec, element);
+ Scalar qmean = volume * (sourceVec[wPhaseIdx]/density_[wPhaseIdx] + sourceVec[nPhaseIdx]/density_[nPhaseIdx]);
+ qmean += rhs_[globalIdx];
+
+ Dune::FieldVector<Scalar, 2 * dim> F(0.);
+ Scalar dInv = 0.;
+ computeReconstructionMatrices(element, W_[globalIdx], F, dInv);
+
+ return (dInv*(qmean + (F*pressTrace)));
+ }
+
+ void constructVelocity(const Element& element, Dune::FieldVector<Scalar,2*dim>& vel, Dune::FieldVector<Scalar,2*dim>& pressTrace, Scalar press)
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ Dune::FieldVector<Scalar, 2 * dim> faceVol(0);
+ IntersectionIterator isItEnd = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ faceVol[isIt->indexInInside()] = isIt->geometry().volume();
+ }
+
+ vel = 0;
+ for (int i = 0; i < 2*dim; i++)
+ for (int j = 0; j < 2*dim; j++)
+ vel[i] += W_[globalIdx][i][j]*faceVol[j]*(press - pressTrace[j]);
+ }
+
+ void constructVelocity(const Element& element, int faceIdx, Scalar& vel, Dune::FieldVector<Scalar,2*dim>& pressTrace, Scalar press)
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ Dune::FieldVector<Scalar, 2 * dim> faceVol(0);
+ IntersectionIterator isItEnd = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ faceVol[isIt->indexInInside()] = isIt->geometry().volume();
+ }
+
+ vel = 0;
+ for (int j = 0; j < 2*dim; j++)
+ vel += W_[globalIdx][faceIdx][j]*faceVol[j]*(press - pressTrace[j]);
+ }
+
+ void computeReconstructionMatrices(const Element& element, const Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim>& W, Dune::FieldVector<Scalar, 2 * dim>& F, Scalar& dInv)
+ {
+ Dune::FieldVector<Scalar, 2 * dim> c(0);
+ Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> Pi(0);
+
+ IntersectionIterator isItEnd = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ // local number of facet
+ int idx = isIt->indexInInside();
+
+ Scalar faceVol = isIt->geometry().volume();
+
+ // Corresponding to the element under consideration,
+ // calculate the part of the matrix C coupling velocities and element pressures.
+ // This is just a row vector of size numFaces.
+ // scale with volume
+ c[idx] = faceVol;
+ // Set up the element part of the matrix \Pi coupling velocities
+ // and pressure-traces. This is a diagonal matrix with entries given by faceVol.
+ Pi[idx][idx] = faceVol;
+ }
+
+ // Calculate the element part of the matrix D^{-1} = (c W c^T)^{-1} which is just a scalar value.
+ Dune::FieldVector<Scalar, 2 * dim> Wc(0);
+ W.umv(c, Wc);
+ dInv = 1.0 / (c * Wc);
+
+ // Calculate the element part of the matrix F = Pi W c^T which is a column vector.
+ F = 0;
+ Pi.umv(Wc, F);
+ }
+
+ void assembleElementMatrices(const Element& element, Dune::FieldVector<Scalar, 2 * dim>& faceVol,
+ Dune::FieldVector<Scalar, 2 * dim>& F,
+ Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim>& Pi,
+ Scalar& dInv,
+ Scalar& qmean);
+
+
+ const Problem& problem() const
+ {
+ return problem_;
+ }
+
+private:
+ void assembleV(const Element& element, int k);
+
+ void assembleBC(const Element& element, int k);
+
+ Scalar evaluateErrorTerm(CellData& cellData)
+ {
+ //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
+ // error reduction routine: volumetric error is damped and inserted to right hand side
+ Scalar sat = 0;
+ switch (saturationType)
+ {
+ case Sw:
+ sat = cellData.saturation(wPhaseIdx);
+ break;
+ case Sn:
+ sat = cellData.saturation(nPhaseIdx);
+ break;
+ }
+
+ Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
+ if (sat < 0.0)
+ {
+ error = sat;
+ }
+ error /= timeStep_;
+
+ Scalar errorAbs = std::abs(error);
+
+ if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
+ && (!problem_.timeManager().willBeFinished()))
+ {
+ return ErrorTermFactor_ * error;
+ }
+ return 0.0;
+ }
+
+private:
+ Problem& problem_;
+ GridView gridView_;
+ Dune::DynamicVector<Scalar> rhs_;
+ std::vector<Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> > W_;
+ Scalar maxError_;
+ Scalar timeStep_;
+ Scalar ErrorTermFactor_; //!< Handling of error term: relaxation factor
+ Scalar ErrorTermLowerBound_; //!< Handling of error term: lower bound for error dampening
+ Scalar ErrorTermUpperBound_; //!< Handling of error term: upper bound for error dampening
+
+ Scalar density_[numPhases];
+ Scalar viscosity_[numPhases];
+};
+
+template<class TypeTag>
+void MimeticTwoPLocalStiffness<TypeTag>::assembleV(const Element& element, int k = 1)
+{
+ unsigned int numFaces = element.template count<1>();
+ this->setcurrentsize(numFaces);
+
+ int globalIdx = problem_.variables().index(element);
+
+ // The notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
+ // The matrix W developed here corresponds to one element-associated
+ // block of the matrix B^{-1} there.
+ Dune::FieldVector<Scalar, 2 * dim> faceVol(0);
+ Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> Pi(0.);
+ Dune::FieldVector<Scalar, 2 * dim> F(0.);
+ Scalar dInv = 0.;
+ Scalar qmean = 0.;
+ this->assembleElementMatrices(element, faceVol, F, Pi, dInv, qmean);
+
+ // Calculate the element part of the matrix Pi W Pi^T.
+ Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> PiWPiT(W_[globalIdx]);
+ PiWPiT.rightmultiply(Pi);
+ PiWPiT.leftmultiply(Pi);
+
+ // Calculate the element part of the matrix F D^{-1} F^T.
+ Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> FDinvFT(0);
+ for (int i = 0; i < numFaces; i++)
+ for (int j = 0; j < numFaces; j++)
+ FDinvFT[i][j] = dInv * F[i] * F[j];
+
+ // Calculate the element part of the matrix S = Pi W Pi^T - F D^{-1} F^T.
+ for (int i = 0; i < numFaces; i++)
+ for (int j = 0; j < numFaces; j++)
+ this->A[i][j] = PiWPiT[i][j] - FDinvFT[i][j];
+
+ // Calculate the source term F D^{-1} f
+ // NOT WORKING AT THE MOMENT
+ Scalar factor = dInv * qmean;
+ for (int i = 0; i < numFaces; i++)
+ this->b[i] = F[i] * factor;
+
+ // std::cout << "faceVol = " << faceVol << std::endl << "W = " << std::endl << W << std::endl
+ // << "c = " << c << std::endl << "Pi = " << std::endl << Pi << std::endl
+ // << "dinv = " << dinv << std::endl << "F = " << F << std::endl;
+ // std::cout << "dinvF = " << dinvF << ", q = " << qmean
+ // << ", b = " << this->b[0] << ", " << this->b[1] << ", " << this->b[2] << ", " << this->b[3] << std::endl;
+}
+
+template<class TypeTag>
+void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element& element,
+ Dune::FieldVector<Scalar, 2 * dim>& faceVol,
+ Dune::FieldVector<Scalar, 2 * dim>& F,
+ Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim>& Pi,
+ Scalar& dInv,
+ Scalar& qmean)
+{
+ unsigned int numFaces = element.template count<1>();
+ this->setcurrentsize(numFaces);
+
+ // get global coordinate of cell center
+ Dune::FieldVector<Scalar, dim> centerGlobal = element.geometry().center();
+
+ int globalIdx = problem_.variables().index(element);
+
+ CellData& cellData = problem_.variables().cellData(globalIdx);
+
+ // cell volume
+ Scalar volume = element.geometry().volume();
+
+ // build the matrices R and ~N
+ Dune::FieldMatrix<Scalar, 2 * dim, dim> R(0), N(0);
+
+ // std::cout << "element " << elemId << ": center " << centerGlobal << std::endl;;
+
+ //collect information needed for calculation of fluxes due to capillary-potential (pc + gravity!)
+ Scalar gravPotFace[2*dim];
+
+ Scalar gravPot = (problem_.bBoxMax() - centerGlobal) * problem_.gravity() * (density_[nPhaseIdx] - density_[wPhaseIdx]);
+
+ int i = -1;
+ IntersectionIterator isItEnd = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ // local number of facet
+ i = isIt->indexInInside();
+
+ Dune::FieldVector<Scalar, dim> faceGlobal = isIt->geometry().center();
+ faceVol[i] = isIt->geometry().volume();
+
+ // get normal vector
+ const Dune::FieldVector<Scalar, dim>& unitOuterNormal = isIt->centerUnitOuterNormal();
+
+ N[i] = unitOuterNormal;
+
+ for (int k = 0; k < dim; k++)
+ // move origin to the center of gravity
+ R[i][k] = faceVol[i] * (faceGlobal[k] - centerGlobal[k]);
+
+ gravPotFace[i] = (problem_.bBoxMax() - faceGlobal) * problem_.gravity() * (density_[nPhaseIdx] - density_[wPhaseIdx]);
+ }
+
+ // proceed along the lines of Algorithm 1 from
+ // Brezzi/Lipnikov/Simonicini M3AS 2005
+ // (1) orthonormalize columns of the matrix R
+ Scalar norm = R[0][0] * R[0][0];
+ for (int i = 1; i < numFaces; i++)
+ norm += R[i][0] * R[i][0];
+ norm = sqrt(norm);
+ for (int i = 0; i < numFaces; i++)
+ R[i][0] /= norm;
+ Scalar weight = R[0][1] * R[0][0];
+ for (int i = 1; i < numFaces; i++)
+ weight += R[i][1] * R[i][0];
+ for (int i = 0; i < numFaces; i++)
+ R[i][1] -= weight * R[i][0];
+ norm = R[0][1] * R[0][1];
+ for (int i = 1; i < numFaces; i++)
+ norm += R[i][1] * R[i][1];
+ norm = sqrt(norm);
+ for (int i = 0; i < numFaces; i++)
+ R[i][1] /= norm;
+ if (dim == 3)
+ {
+ Scalar weight1 = R[0][2] * R[0][0];
+ Scalar weight2 = R[0][2] * R[0][1];
+ for (int i = 1; i < numFaces; i++)
+ {
+ weight1 += R[i][2] * R[i][0];
+ weight2 += R[i][2] * R[i][1];
+ }
+ for (int i = 0; i < numFaces; i++)
+ R[i][2] -= weight1 * R[i][0] + weight2 * R[i][1];
+ norm = R[0][2] * R[0][2];
+ for (int i = 1; i < numFaces; i++)
+ norm += R[i][2] * R[i][2];
+ norm = sqrt(norm);
+ for (int i = 0; i < numFaces; i++)
+ R[i][2] /= norm;
+ }
+ // std::cout << "~R =\dim" << R;
+
+ // (2) Build the matrix ~D
+ Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> D(0);
+ for (int s = 0; s < numFaces; s++)
+ {
+ Dune::FieldVector<Scalar, 2 * dim> es(0);
+ es[s] = 1;
+ for (int k = 0; k < numFaces; k++)
+ {
+ D[k][s] = es[k];
+ for (int i = 0; i < dim; i++)
+ {
+ D[k][s] -= R[s][i] * R[k][i];
+ }
+ }
+ }
+
+// std::cout<<"D = "<<D<<"\n";
+
+ // eval diffusion tensor, ASSUMING to be constant over each cell
+ Dune::FieldMatrix<Scalar, dim, dim> mobility(problem_.spatialParams().intrinsicPermeability(element));
+ mobility *= (cellData.mobility(wPhaseIdx) + cellData.mobility(nPhaseIdx));
+
+ Scalar traceMob = mobility[0][0];
+ for (int i = 1; i < dim; i++)
+ traceMob += mobility[i][i];
+ D *= 2.0 * traceMob / volume;
+ // std::cout << "u~D =\dim" << D;
+
+ // (3) Build the matrix W = Minv
+ Dune::FieldMatrix<Scalar, 2 * dim, dim> NK(N);
+ NK.rightmultiply(mobility);
+
+ for (int i = 0; i < numFaces; i++)
+ {
+ for (int j = 0; j < numFaces; j++)
+ {
+ W_[globalIdx][i][j] = NK[i][0] * N[j][0];
+ for (int k = 1; k < dim; k++)
+ W_[globalIdx][i][j] += NK[i][k] * N[j][k];
+ }
+ }
+
+ W_[globalIdx] /= volume;
+ W_[globalIdx] += D;
+
+ // std::cout << "W = \dim" << W;
+
+ // Now the notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
+ // The matrix W developed so far corresponds to one element-associated
+ // block of the matrix B^{-1} there.
+
+ // Corresponding to the element under consideration,
+ // calculate the part of the matrix C coupling velocities and element pressures.
+ // This is just a row vector of size numFaces.
+ // scale with volume
+ Dune::FieldVector<Scalar, 2 * dim> c(0);
+ for (int i = 0; i < numFaces; i++)
+ c[i] = faceVol[i];
+
+ // Set up the element part of the matrix \Pi coupling velocities
+ // and pressure-traces. This is a diagonal matrix with entries given by faceVol.
+ for (int i = 0; i < numFaces; i++)
+ Pi[i][i] = faceVol[i];
+
+ // Calculate the element part of the matrix D^{-1} = (c W c^T)^{-1} which is just a scalar value.
+ Dune::FieldVector<Scalar, 2 * dim> Wc(0);
+ W_[globalIdx].umv(c, Wc);
+ dInv = 1.0 / (c * Wc);
+
+ // Calculate the element part of the matrix F = Pi W c^T which is a column vector.
+ F = 0;
+ Pi.umv(Wc, F);
+ // std::cout << "Pi = \dim" << Pi << "c = " << c << ", F = " << F << std::endl;
+
+ //accumulate fluxes due to capillary potential (pc + gravity!)
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ int idx = isIt->indexInInside();
+
+ Scalar fracFlow = 0;
+
+ Scalar flux = 0;
+ for (int j = 0; j < 2 * dim; j++)
+ flux += W_[globalIdx][idx][j] * faceVol[j] * (gravPot - gravPotFace[j]);
+
+ //it is enough to evaluate the capillary/gravity flux for neighbors -> not needed anyway at the boundaries!
+ if (isIt->neighbor())
+ {
+ ElementPointer neighbor = isIt->outside();
+ int globalIdxNeighbor = problem_.variables().index(*neighbor);
+ if (flux >= 0.)
+ {
+ switch (pressureType)
+ {
+ case pw:
+ {
+ fracFlow = cellData.fracFlowFunc(nPhaseIdx);
+ break;
+ }
+ case pn:
+ {
+ fracFlow = -cellData.fracFlowFunc(wPhaseIdx);
+ break;
+ }
+ }
+
+
+ rhs_[globalIdx] -= (faceVol[idx] * fracFlow * flux);
+ rhs_[globalIdxNeighbor] += (faceVol[idx] * fracFlow * flux);
+ }
+ }
+ else if (isIt->boundary())
+ {
+ BoundaryTypes bctype;
+ problem_.boundaryTypes(bctype, *isIt);
+
+ if (bctype.isDirichlet(pressEqIdx))
+ {
+ if (flux > 0. || !bctype.isDirichlet(satEqIdx))
+ {
+ switch (pressureType)
+ {
+ case pw:
+ {
+ fracFlow = cellData.fracFlowFunc(nPhaseIdx);
+ break;
+ }
+ case pn:
+ {
+ fracFlow = -cellData.fracFlowFunc(wPhaseIdx);
+ break;
+ }
+ }
+ }
+ else if (flux < 0. && bctype.isDirichlet(satEqIdx))
+ {
+ PrimaryVariables boundValues(0.0);
+ problem_.dirichlet(boundValues, *isIt);
+
+ Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
+ boundValues[saturationIdx]);
+ Scalar krn = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element),
+ boundValues[saturationIdx]);
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ fracFlow = krn / viscosity_[nPhaseIdx]
+ / (krw / viscosity_[wPhaseIdx] + krn / viscosity_[nPhaseIdx]);
+ break;
+ }
+ case pn:
+ {
+ fracFlow = -krw / viscosity_[wPhaseIdx]
+ / (krw / viscosity_[wPhaseIdx] + krn / viscosity_[nPhaseIdx]);
+ break;
+ }
+ }
+ }
+
+ rhs_[globalIdx] -= (faceVol[idx] * fracFlow * flux);
+ }
+ }
+ }
+
+ PrimaryVariables sourceVec(0.0);
+ problem_.source(sourceVec, element);
+ qmean = volume * (sourceVec[wPhaseIdx]/density_[wPhaseIdx] + sourceVec[nPhaseIdx]/density_[nPhaseIdx]);
+
+ qmean += evaluateErrorTerm(cellData) * volume;
+}
+
+template<class TypeTag>
+void MimeticTwoPLocalStiffness<TypeTag>::assembleBC(const Element& element, int k = 1)
+{
+ // evaluate boundary conditions via intersection iterator
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ IntersectionIterator endIsIt = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != endIsIt; ++isIt)
+ {
+ int faceIndex = isIt->indexInInside();
+ if (isIt->boundary())
+ {
+ problem_.boundaryTypes(this->bctype[faceIndex], *isIt);
+ PrimaryVariables boundValues(0.0);
+
+ if (this->bctype[faceIndex].isNeumann(pressEqIdx))
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ problem_.neumann(boundValues, *isIt);
+ Scalar J = (boundValues[wPhaseIdx]/density_[wPhaseIdx] + boundValues[nPhaseIdx]/density_[nPhaseIdx]);
+ this->b[faceIndex] -= J * isIt->geometry().volume();
+ }
+ else if (this->bctype[faceIndex].isDirichlet(pressEqIdx))
+ {
+ problem_.dirichlet(boundValues, *isIt);
+ if (pressureType == pw)
+ this->b[faceIndex] = boundValues[pressureIdx] + (problem_.bBoxMax() - isIt->geometry().center()) * problem_.gravity() * density_[wPhaseIdx];
+ else if (pressureType == pn)
+ this->b[faceIndex] = boundValues[pressureIdx] + (problem_.bBoxMax() - isIt->geometry().center()) * problem_.gravity() * density_[nPhaseIdx];
+ else
+ this->b[faceIndex] = boundValues[pressureIdx];
+ }
+ }
+ }
+
+}
+
+/** @} */
+}
+#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimetic2padaptive.hh b/dumux/decoupled/2p/diffusion/mimetic/mimetic2padaptive.hh
new file mode 100644
index 0000000000..90cb1be901
--- /dev/null
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimetic2padaptive.hh
@@ -0,0 +1,775 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Markus Wolff *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Local stiffness matrix for the diffusion equation discretized by mimetic FD
+ */
+#ifndef DUMUX_MIMETIC2PADAPTIVE_HH
+#define DUMUX_MIMETIC2PADAPTIVE_HH
+
+#include<map>
+#include<iostream>
+#include<iomanip>
+#include<fstream>
+#include<vector>
+#include<sstream>
+
+#include<dune/common/exceptions.hh>
+#include<dune/grid/common/grid.hh>
+#include<dune/geometry/type.hh>
+#include<dune/geometry/quadraturerules.hh>
+
+#include<dumux/common/boundaryconditions.hh>
+
+#include <dune/common/dynvector.hh>
+#include <dune/common/dynmatrix.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mimetic2p
+ */
+/**
+ * @brief compute local stiffness matrix for conforming finite elements for the full 2-phase pressure equation
+ *
+ */
+
+//! A class for computing local stiffness matrices
+/*! A class for computing local stiffness matrix for the full 2-phase pressure equation
+ */
+template<class TypeTag>
+class MimeticTwoPLocalStiffnessAdaptive: public LocalStiffness<TypeTag, 1>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ // grid types
+ enum
+ {
+ dim = GridView::dimension
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+ pressEqIdx = Indices::pressEqIdx,
+ satEqIdx = Indices::satEqIdx,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation), //!< gives kind of pressure used (\f$ 0 = p_w\f$, \f$ 1 = p_n\f$, \f$ 2 = p_{global}\f$)
+ saturationType = GET_PROP_VALUE(TypeTag, SaturationFormulation) //!< gives kind of saturation used (\f$ 0 = S_w\f$, \f$ 1 = S_n\f$)
+ };
+
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+
+ typedef Dumux::IntersectionMapper<GridView> IntersectionMapper;
+
+public:
+ // define the number of components of your system, this is used outside
+ // to allocate the correct size of (dense) blocks with a FieldMatrix
+ enum
+ {
+ m = 1
+ };
+ enum
+ {
+ size = 2 * dim
+ };
+
+ //! Constructor
+ MimeticTwoPLocalStiffnessAdaptive(Problem& problem, bool levelBoundaryAsDirichlet,
+ const GridView& gridView, const IntersectionMapper& intersectionMapper,bool procBoundaryAsDirichlet = true) :
+ problem_(problem), gridView_(gridView), intersectionMapper_(intersectionMapper), maxError_(0), timeStep_(1)
+ {
+ int size = gridView_.size(0);
+ rhs_.resize(size , 0.);
+ W_.resize(size);
+ ErrorTermFactor_ = GET_PARAM(TypeTag, Scalar, ImpetErrorTermFactor);
+ ErrorTermLowerBound_ = GET_PARAM(TypeTag, Scalar, ImpetErrorTermLowerBound);
+ ErrorTermUpperBound_ = GET_PARAM(TypeTag, Scalar, ImpetErrorTermUpperBound);
+
+ density_[wPhaseIdx] = 0.0;
+ density_[nPhaseIdx] = 0.0;
+ viscosity_[wPhaseIdx] =0.0;
+ viscosity_[nPhaseIdx] = 0.0;
+ }
+
+ void initialize()
+ {
+ ElementIterator element = problem_.gridView().template begin<0>();
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setPressure(nPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setTemperature(problem_.temperature(*element));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
+ density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
+ adapt();
+ }
+
+ void adapt()
+ {
+ int size = gridView_.size(0);
+ rhs_.resize(size , 0.);
+ W_.resize(size);
+ }
+
+ void reset()
+ {
+ rhs_ = 0;
+ maxError_ = 0;
+ timeStep_ = 1;
+ }
+
+ void setErrorInfo(Scalar maxErr, Scalar dt)
+ {
+ maxError_ = maxErr;
+ timeStep_ = dt;
+ }
+
+ int numberOfFaces(int globalIdx)
+ {
+ return W_[globalIdx].N();
+ }
+
+ //! assemble local stiffness matrix for given element and order
+ /*! On exit the following things have been done:
+ - The stiffness matrix for the given entity and polynomial degree has been assembled and is
+ accessible with the mat() method.
+ - The boundary conditions have been evaluated and are accessible with the bc() method
+ - The right hand side has been assembled. It contains either the value of the essential boundary
+ condition or the assembled source term and neumann boundary condition. It is accessible via the rhs() method.
+ @param[in] element a codim 0 entity reference
+ @param[in] k order of CR basis (only k = 1 is implemented)
+ */
+ void assemble(const Element& element, int k = 1)
+ {
+ int numFaces = intersectionMapper_.size(element);
+
+ this->setcurrentsize(numFaces);
+
+ // clear assemble data
+ for (int i = 0; i < numFaces; i++)
+ {
+ this->b[i] = 0;
+ this->bctype[i].reset();
+ for (int j = 0; j < numFaces; j++)
+ this->A[i][j] = 0;
+ }
+
+ assembleV(element, numFaces, k);
+ assembleBC(element, k);
+ }
+
+ // TODO/FIXME: this is only valid for linear problems where
+ // the local stiffness matrix is independend of the current
+ // solution. We need to implement this properly, but this
+ // should at least make the thing compile...
+ typedef Dune::FieldVector<Scalar, m> VBlockType;
+ void assemble(const Element &cell, const Dune::BlockVector<VBlockType>& localSolution, int orderOfShapeFns = 1)
+ {
+ assemble(cell, orderOfShapeFns);
+ }
+
+ //! assemble only boundary conditions for given element
+ /*! On exit the following things have been done:
+ - The boundary conditions have been evaluated and are accessible with the bc() method
+ - The right hand side contains either the value of the essential boundary
+ condition or the assembled neumann boundary condition. It is accessible via the rhs() method.
+ @param[in] element a codim 0 entity reference
+ @param[in] k order of CR basis
+ */
+ void assembleBoundaryCondition(const Element& element, int k = 1)
+ {
+ int numFaces = intersectionMapper_.size(element);
+
+ // clear assemble data
+ for (int i = 0; i < numFaces; i++)
+ {
+ this->b[i] = 0;
+ this->bctype[i].reset();
+ }
+
+ assembleBC(element, k);
+ }
+
+ template<class Vector>
+ void completeRHS(const Element& element, std::vector<int>& local2Global, Vector& f)
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ int numFaces = numberOfFaces(globalIdx);
+
+ Dune::DynamicVector<Scalar> F(numFaces);
+ Scalar dInv = 0.;
+ computeReconstructionMatrices(element, W_[globalIdx], F, dInv);
+
+ for (int i = 0; i < numFaces; i++)
+ {
+// if (!this->bctype[i].isSet())
+ f[local2Global[i]][0] += (dInv * F[i] * rhs_[globalIdx]);
+ }
+ }
+
+ Scalar constructPressure(const Element& element, Dune::DynamicVector<Scalar>& pressTrace)
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ int numFaces = numberOfFaces(globalIdx);
+
+ Scalar volume = element.geometry().volume();
+
+ PrimaryVariables sourceVec(0.0);
+ problem_.source(sourceVec, element);
+ Scalar qmean = volume * (sourceVec[wPhaseIdx]/density_[wPhaseIdx] + sourceVec[nPhaseIdx]/density_[nPhaseIdx]);
+ qmean += rhs_[globalIdx];
+
+ Dune::DynamicVector<Scalar> F(numFaces, 0.);
+ Scalar dInv = 0.;
+ computeReconstructionMatrices(element, W_[globalIdx], F, dInv);
+
+ return (dInv*(qmean + (F*pressTrace)));
+ }
+
+ void constructVelocity(const Element& element, Dune::DynamicVector<Scalar>& vel, Dune::DynamicVector<Scalar>& pressTrace, Scalar press)
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ int numFaces = numberOfFaces(globalIdx);
+
+ Dune::DynamicVector<Scalar> faceVol(numFaces);
+
+ Dune::FieldVector<Scalar, 2*dim> faceVolumeReal(0.0);
+
+ int idx = 0;
+ IntersectionIterator isItEnd = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ faceVol[idx] = isIt->geometry().volume();
+ int indexInInside = isIt->indexInInside();
+ faceVolumeReal[indexInInside] += faceVol[idx];
+
+ idx++;
+ }
+
+ vel = 0;
+ for (int i = 0; i < numFaces; i++)
+ for (int j = 0; j < numFaces; j++)
+ vel[i] += W_[globalIdx][i][j]*faceVol[j]*(press - pressTrace[j]);
+
+ for (int i = 0; i < numFaces; i++)
+ {
+ vel[i] *= faceVol[i]/faceVolumeReal[intersectionMapper_.maplocal(globalIdx, i)];
+ }
+
+ }
+
+ void computeReconstructionMatrices(const Element& element, const Dune::DynamicMatrix<Scalar>& W, Dune::DynamicVector<Scalar>& F, Scalar& dInv)
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ int numFaces = numberOfFaces(globalIdx);
+
+ Dune::DynamicVector<Scalar> c(numFaces);
+ Dune::DynamicMatrix<Scalar> Pi(numFaces, numFaces);
+
+ int idx = 0;
+ IntersectionIterator isItEnd = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ Scalar faceVol = isIt->geometry().volume();
+
+ // Corresponding to the element under consideration,
+ // calculate the part of the matrix C coupling velocities and element pressures.
+ // This is just a row vector of size numFaces.
+ // scale with volume
+ c[idx] = faceVol;
+ // Set up the element part of the matrix \Pi coupling velocities
+ // and pressure-traces. This is a diagonal matrix with entries given by faceVol.
+ Pi[idx][idx] = faceVol;
+ idx++;
+ }
+
+ // Calculate the element part of the matrix D^{-1} = (c W c^T)^{-1} which is just a scalar value.
+ Dune::DynamicVector<Scalar> Wc(numFaces);
+ W.umv(c, Wc);
+ dInv = 1.0 / (c * Wc);
+
+ // Calculate the element part of the matrix F = Pi W c^T which is a column vector.
+ F = 0;
+ Pi.umv(Wc, F);
+ }
+
+ void assembleElementMatrices(const Element& element,
+ Dune::DynamicVector<Scalar>& faceVol,
+ Dune::DynamicVector<Scalar>& F,
+ Dune::DynamicMatrix<Scalar>& Pi,
+ Scalar& dInv,
+ Scalar& qmean);
+
+
+ const Problem& problem() const
+ {
+ return problem_;
+ }
+
+private:
+ void assembleV(const Element& element, int numFaces, int k);
+
+ void assembleBC(const Element& element, int k);
+
+ Scalar evaluateErrorTerm(CellData& cellData)
+ {
+ //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
+ // error reduction routine: volumetric error is damped and inserted to right hand side
+ Scalar sat = 0;
+ switch (saturationType)
+ {
+ case Sw:
+ sat = cellData.saturation(wPhaseIdx);
+ break;
+ case Sn:
+ sat = cellData.saturation(nPhaseIdx);
+ break;
+ }
+
+ Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
+ if (sat < 0.0)
+ {
+ error = sat;
+ }
+ error /= timeStep_;
+
+ Scalar errorAbs = std::abs(error);
+
+ if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
+ && (!problem_.timeManager().willBeFinished()))
+ {
+ return ErrorTermFactor_ * error;
+ }
+ return 0.0;
+ }
+
+private:
+ Problem& problem_;
+ const GridView& gridView_;
+ const IntersectionMapper& intersectionMapper_;
+ Dune::DynamicVector<Scalar> rhs_;
+ std::vector<Dune::DynamicMatrix<Scalar> > W_;
+ Scalar maxError_;
+ Scalar timeStep_;
+ Scalar ErrorTermFactor_; //!< Handling of error term: relaxation factor
+ Scalar ErrorTermLowerBound_; //!< Handling of error term: lower bound for error dampening
+ Scalar ErrorTermUpperBound_; //!< Handling of error term: upper bound for error dampening
+
+ Scalar density_[numPhases];
+ Scalar viscosity_[numPhases];
+};
+
+template<class TypeTag>
+void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleV(const Element& element, int numFaces, int k = 1)
+{
+ int globalIdx = problem_.variables().index(element);
+
+ // The notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
+ // The matrix W developed here corresponds to one element-associated
+ // block of the matrix B^{-1} there.
+ W_[globalIdx].resize(numFaces, numFaces);
+ Dune::DynamicVector<Scalar> faceVol(numFaces);
+ Dune::DynamicMatrix<Scalar> Pi(numFaces, numFaces);
+ Dune::DynamicVector<Scalar> F(numFaces);
+ Scalar dInv = 0.;
+ Scalar qmean = 0.;
+ this->assembleElementMatrices(element, faceVol, F, Pi, dInv, qmean);
+
+ // Calculate the element part of the matrix Pi W Pi^T.
+ Dune::DynamicMatrix<Scalar> PiWPiT(W_[globalIdx]);
+ PiWPiT.rightmultiply(Pi);
+ PiWPiT.leftmultiply(Pi);
+
+ // Calculate the element part of the matrix F D^{-1} F^T.
+ Dune::DynamicMatrix<Scalar> FDinvFT(numFaces, numFaces);
+ for (int i = 0; i < numFaces; i++)
+ for (int j = 0; j < numFaces; j++)
+ FDinvFT[i][j] = dInv * F[i] * F[j];
+
+ // Calculate the element part of the matrix S = Pi W Pi^T - F D^{-1} F^T.
+ for (int i = 0; i < numFaces; i++)
+ for (int j = 0; j < numFaces; j++)
+ this->A[i][j] = PiWPiT[i][j] - FDinvFT[i][j];
+
+ // Calculate the source term F D^{-1} f
+ // NOT WORKING AT THE MOMENT
+ Scalar factor = dInv * qmean;
+ for (int i = 0; i < numFaces; i++)
+ this->b[i] = F[i] * factor;
+
+ // std::cout << "faceVol = " << faceVol << std::endl << "W = " << std::endl << W << std::endl
+ // << "c = " << c << std::endl << "Pi = " << std::endl << Pi << std::endl
+ // << "dinv = " << dinv << std::endl << "F = " << F << std::endl;
+ // std::cout << "dinvF = " << dinvF << ", q = " << qmean
+ // << ", b = " << this->b[0] << ", " << this->b[1] << ", " << this->b[2] << ", " << this->b[3] << std::endl;
+}
+
+template<class TypeTag>
+void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const Element& element,
+ Dune::DynamicVector<Scalar>& faceVol,
+ Dune::DynamicVector<Scalar>& F,
+ Dune::DynamicMatrix<Scalar>& Pi,
+ Scalar& dInv,
+ Scalar& qmean)
+{
+ // get global coordinate of cell center
+ const Dune::FieldVector<Scalar, dim>& centerGlobal = element.geometry().center();
+
+ int globalIdx = problem_.variables().index(element);
+
+ CellData& cellData = problem_.variables().cellData(globalIdx);
+
+ int numFaces = intersectionMapper_.size(globalIdx);
+
+ // cell volume
+ Scalar volume = element.geometry().volume();
+
+ // build the matrices R and ~N
+ Dune::DynamicMatrix<Scalar> R(numFaces, dim, 0.);
+ Dune::DynamicMatrix<Scalar> N(numFaces, dim, 0.);
+
+ // std::cout << "element " << elemId << ": center " << centerGlobal << std::endl;;
+
+ //collect information needed for calculation of fluxes due to capillary-potential (pc + gravity!)
+ std::vector<Scalar> pcPotFace(numFaces);
+ Scalar pcPot = (problem_.bBoxMax() - element.geometry().center()) * problem_.gravity() * (density_[nPhaseIdx] - density_[wPhaseIdx]);
+
+ int idx = 0;
+ IntersectionIterator isItEnd = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ // local number of facet
+
+ Dune::FieldVector<Scalar, dim> faceGlobal(isIt->geometry().center());
+ faceVol[idx] = isIt->geometry().volume();
+
+ // get normal vector
+ const Dune::FieldVector<Scalar, dim>& unitOuterNormal = isIt->centerUnitOuterNormal();
+
+
+ for (int k = 0; k < dim; k++)
+ {
+ // move origin to the center of gravity
+ R[idx][k] = faceVol[idx] * (faceGlobal[k] - centerGlobal[k]);
+ N[idx][k] = unitOuterNormal[k];
+ }
+
+ pcPotFace[idx] = (problem_.bBoxMax() - faceGlobal) * problem_.gravity() * (density_[nPhaseIdx] - density_[wPhaseIdx]);
+
+ idx++;
+ }
+
+ // proceed along the lines of Algorithm 1 from
+ // Brezzi/Lipnikov/Simonicini M3AS 2005
+ // (1) orthonormalize columns of the matrix R (Gram-Schmidt orthonormalization)
+ Scalar norm = R[0][0] * R[0][0];
+ for (int i = 1; i < numFaces; i++)
+ norm += R[i][0] * R[i][0];
+ norm = sqrt(norm);
+ for (int i = 0; i < numFaces; i++)
+ R[i][0] /= norm;
+ Scalar weight = R[0][1] * R[0][0];
+ for (int i = 1; i < numFaces; i++)
+ weight += R[i][1] * R[i][0];
+ for (int i = 0; i < numFaces; i++)
+ R[i][1] -= weight * R[i][0];
+ norm = R[0][1] * R[0][1];
+ for (int i = 1; i < numFaces; i++)
+ norm += R[i][1] * R[i][1];
+ norm = sqrt(norm);
+ for (int i = 0; i < numFaces; i++)
+ R[i][1] /= norm;
+ if (dim == 3)
+ {
+ Scalar weight1 = R[0][2] * R[0][0];
+ Scalar weight2 = R[0][2] * R[0][1];
+ for (int i = 1; i < numFaces; i++)
+ {
+ weight1 += R[i][2] * R[i][0];
+ weight2 += R[i][2] * R[i][1];
+ }
+ for (int i = 0; i < numFaces; i++)
+ R[i][2] -= weight1 * R[i][0] + weight2 * R[i][1];
+ norm = R[0][2] * R[0][2];
+ for (int i = 1; i < numFaces; i++)
+ norm += R[i][2] * R[i][2];
+ norm = sqrt(norm);
+ for (int i = 0; i < numFaces; i++)
+ R[i][2] /= norm;
+ }
+ // std::cout << "~R =\dim" << R;
+
+ // (2) Build the matrix ~D
+ Dune::DynamicMatrix<Scalar> D(numFaces, numFaces, 0.);
+ for (int s = 0; s < numFaces; s++)
+ {
+ Dune::DynamicVector<Scalar> es(numFaces, 0.);
+ es[s] = 1;
+ for (int k = 0; k < numFaces; k++)
+ {
+ D[k][s] = es[k];
+ for (int i = 0; i < dim; i++)
+ {
+ D[k][s] -= R[s][i] * R[k][i];
+ }
+ }
+ }
+// std::cout<<"D = "<<D<<"\n";
+
+ // eval diffusion tensor, ASSUMING to be constant over each cell
+ Dune::FieldMatrix<Scalar, dim, dim> mobility(problem_.spatialParams().intrinsicPermeability(element));
+ mobility *= (cellData.mobility(wPhaseIdx) + cellData.mobility(nPhaseIdx));
+
+ Scalar traceMob = mobility[0][0];
+ for (int i = 1; i < dim; i++)
+ traceMob += mobility[i][i];
+ D *= 2.0 * traceMob / volume;
+ // std::cout << "u~D =\dim" << D;
+
+ // (3) Build the matrix W = Minv
+ Dune::DynamicMatrix<Scalar> NK(N);
+ NK.rightmultiply(mobility);
+
+// std::cout<<"NK = "<<NK<<"\n";
+
+ for (int i = 0; i < numFaces; i++)
+ {
+ for (int j = 0; j < numFaces; j++)
+ {
+ W_[globalIdx][i][j] = NK[i][0] * N[j][0];
+ for (int k = 1; k < dim; k++)
+ W_[globalIdx][i][j] += NK[i][k] * N[j][k];
+ }
+ }
+
+ W_[globalIdx] /= volume;
+ W_[globalIdx] += D;
+
+ // std::cout << "W = \dim" << W;
+
+ // Now the notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
+ // The matrix W developed so far corresponds to one element-associated
+ // block of the matrix B^{-1} there.
+
+ // Corresponding to the element under consideration,
+ // calculate the part of the matrix C coupling velocities and element pressures.
+ // This is just a row vector of size numFaces.
+ // scale with volume
+ Dune::DynamicVector<Scalar> c(numFaces);
+ for (int i = 0; i < numFaces; i++)
+ c[i] = faceVol[i];
+
+ // Set up the element part of the matrix \Pi coupling velocities
+ // and pressure-traces. This is a diagonal matrix with entries given by faceVol.
+ for (int i = 0; i < numFaces; i++)
+ Pi[i][i] = faceVol[i];
+
+ // Calculate the element part of the matrix D^{-1} = (c W c^T)^{-1} which is just a scalar value.
+ Dune::DynamicVector<Scalar> Wc(numFaces);
+ W_[globalIdx].umv(c, Wc);
+ dInv = 1.0 / (c * Wc);
+
+ // Calculate the element part of the matrix F = Pi W c^T which is a column vector.
+ F = 0;
+ Pi.umv(Wc, F);
+ // std::cout << "Pi = \dim" << Pi << "c = " << c << ", F = " << F << std::endl;
+
+ //accumulate fluxes due to capillary potential (pc + gravity!)
+ idx = 0;
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isItEnd; ++isIt)
+ {
+ Scalar fracFlow = 0;
+
+ Scalar flux = 0;
+ for (int j = 0; j < numFaces; j++)
+ flux += W_[globalIdx][idx][j] * faceVol[j] * (pcPot - pcPotFace[j]);
+
+ if (isIt->neighbor())
+ {
+
+ ElementPointer neighbor = isIt->outside();
+ int globalIdxNeighbor = problem_.variables().index(*neighbor);
+ if (flux > 0.)
+ {
+ switch (pressureType)
+ {
+ case pw:
+ {
+ fracFlow = cellData.fracFlowFunc(nPhaseIdx);
+ break;
+ }
+ case pn:
+ {
+ fracFlow = -cellData.fracFlowFunc(wPhaseIdx);
+ break;
+ }
+ }
+
+ rhs_[globalIdx] -= (faceVol[idx] * fracFlow * flux);
+ rhs_[globalIdxNeighbor] += (faceVol[idx] * fracFlow * flux);
+ }
+ }
+ else if (isIt->boundary())
+ {
+ BoundaryTypes bctype;
+ problem_.boundaryTypes(bctype, *isIt);
+
+ if (bctype.isDirichlet(pressEqIdx))
+ {
+ if (flux > 0. || !bctype.isDirichlet(satEqIdx))
+ {
+ switch (pressureType)
+ {
+ case pw:
+ {
+ fracFlow = cellData.fracFlowFunc(nPhaseIdx);
+ break;
+ }
+ case pn:
+ {
+ fracFlow = -cellData.fracFlowFunc(wPhaseIdx);
+ break;
+ }
+ }
+ }
+ else if (flux < 0. && bctype.isDirichlet(satEqIdx))
+ {
+ PrimaryVariables boundValues(0.0);
+ problem_.dirichlet(boundValues, *isIt);
+
+ Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
+ boundValues[saturationIdx]);
+ Scalar krn = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element),
+ boundValues[saturationIdx]);
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ fracFlow = krn / viscosity_[nPhaseIdx]
+ / (krw / viscosity_[wPhaseIdx] + krn / viscosity_[nPhaseIdx]);
+ break;
+ }
+ case pn:
+ {
+ fracFlow = -krw / viscosity_[wPhaseIdx]
+ / (krw / viscosity_[wPhaseIdx] + krn / viscosity_[nPhaseIdx]);
+ break;
+ }
+ }
+ }
+
+ rhs_[globalIdx] -= (faceVol[idx] * fracFlow * flux);
+ }
+ }
+ idx++;
+ }
+
+ PrimaryVariables sourceVec(0.0);
+ problem_.source(sourceVec, element);
+ qmean = volume * (sourceVec[wPhaseIdx]/density_[wPhaseIdx] + sourceVec[nPhaseIdx]/density_[nPhaseIdx]);
+
+ qmean += evaluateErrorTerm(cellData) * volume;
+}
+
+template<class TypeTag>
+void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleBC(const Element& element, int k = 1)
+{
+ // evaluate boundary conditions via intersection iterator
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ unsigned int faceIndex = 0;
+ IntersectionIterator endIsIt = gridView_.iend(element);
+ for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != endIsIt; ++isIt)
+ {
+ if (isIt->neighbor())
+ {
+
+ }
+ else if (isIt->boundary())
+ {
+ problem_.boundaryTypes(this->bctype[faceIndex], *isIt);
+ PrimaryVariables boundValues(0.0);
+
+ if (this->bctype[faceIndex].isNeumann(pressEqIdx))
+ {
+ int globalIdx = problem_.variables().index(element);
+
+ problem_.neumann(boundValues, *isIt);
+ Scalar J = (boundValues[wPhaseIdx]/density_[wPhaseIdx] + boundValues[nPhaseIdx]/density_[nPhaseIdx]);
+ this->b[faceIndex] -= J * isIt->geometry().volume();
+ }
+ else if (this->bctype[faceIndex].isDirichlet(pressEqIdx))
+ {
+ problem_.dirichlet(boundValues, *isIt);
+ if (pressureType == pw)
+ this->b[faceIndex] = boundValues[pressureIdx] + (problem_.bBoxMax() - isIt->geometry().center()) * problem_.gravity() * density_[wPhaseIdx];
+ else if (pressureType == pn)
+ this->b[faceIndex] = boundValues[pressureIdx] + (problem_.bBoxMax() - isIt->geometry().center()) * problem_.gravity() * density_[nPhaseIdx];
+ else
+ this->b[faceIndex] = boundValues[pressureIdx];
+ }
+ }
+ faceIndex++;
+ }
+}
+
+/** @} */
+}
+#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
deleted file mode 100644
index 286f1f83a9..0000000000
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
+++ /dev/null
@@ -1,436 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program 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 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Local stiffness matrix for the diffusion equation discretized by mimetic FD
- */
-#ifndef DUMUX_MIMETICGROUNDWATER_HH
-#define DUMUX_MIMETICGROUNDWATER_HH
-
-#include<map>
-#include<iostream>
-#include<iomanip>
-#include<fstream>
-#include<vector>
-#include<sstream>
-
-#include<dune/geometry/type.hh>
-#include<dune/geometry/quadraturerules.hh>
-#include<dune/grid/common/grid.hh>
-
-#include<dumux/common/boundaryconditions.hh>
-#include "localstiffness.hh"
-#include <dumux/decoupled/2p/diffusion/diffusionproperties2p.hh>
-namespace Dumux
-{
-/*!
- * \ingroup MimeticPressure2p
- */
-/**
- * @brief compute local stiffness matrix for conforming finite elements for diffusion equation
- *
- */
-
-//! A class for computing local stiffness matrices
-/*! A class for computing local stiffness matrix for the diffusion equation
- * \f{eqnarray*}{
- * \text{div} \boldsymbol v_{total} &=& q \\
- * \boldsymbol v_{total} &=& - \lambda \boldsymbol K \textbf{grad} p
- * \f}
- * where \f$ p = p_D \f$ on \f$ \Gamma_{Dirichlet} \f$, and \f$ \boldsymbol v_{total} \cdot \boldsymbol n = q_N \f$ on \f$ \Gamma_{Neumann} \f$.
- *
- * \tparam TypeTag The problem TypeTag
- */
-template<class TypeTag>
-class MimeticGroundwaterEquationLocalStiffness: public LocalStiffness<TypeTag, 1>
-{
- typedef typename GET_PROP_TYPE(TypeTag, GridView)GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-
- // grid types
- enum
- {
- dim = GridView::dimension
- };
- enum
- {
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
- pressureIdx = Indices::pressureIdx,
- pressEqIdx = Indices::pressureEqIdx,
- numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
- };
- typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
-
- typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
- typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
-
- typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
-
-public:
- // define the number of components of your system, this is used outside
- // to allocate the correct size of (dense) blocks with a FieldMatrix
- enum
- { m=1};
- enum
- { size=2*dim};
-
- //! Constructor
- MimeticGroundwaterEquationLocalStiffness (Problem& problem,
- bool levelBoundaryAsDirichlet, const GridView& gridView,
- bool procBoundaryAsDirichlet=true)
- : problem_(problem), gridView_(gridView)
- {}
-
- //! For initialization
- void initialize()
- {
- ElementIterator element = problem_.gridView().template begin<0> ();
- FluidState fluidState;
- fluidState.setPressure(wPhaseIdx, problem_.referencePressure(*element));
- fluidState.setPressure(nPhaseIdx, problem_.referencePressure(*element));
- fluidState.setTemperature(problem_.temperature(*element));
- fluidState.setSaturation(wPhaseIdx, 1.);
- fluidState.setSaturation(nPhaseIdx, 0.);
- density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
- density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
- }
-
- //! assemble local stiffness matrix for given element and order
- /*! On exit the following things have been done:
- - The stiffness matrix for the given entity and polynomial degree has been assembled and is
- accessible with the mat() method.
- - The boundary conditions have been evaluated and are accessible with the bc() method
- - The right hand side has been assembled. It contains either the value of the essential boundary
- condition or the assembled source term and neumann boundary condition. It is accessible via the rhs() method.
- @param[in] element a codim 0 entity reference
- @param[in] k order of CR basis (only k = 1 is implemented)
- */
- void assemble (const Element& element, int k=1)
- {
- unsigned int numFaces = element.template count<1>();
- this->setcurrentsize(numFaces);
-
- // clear assemble data
- for (unsigned int i=0; i<numFaces; i++)
- {
- this->b[i] = 0;
- this->bctype[i].reset();
- for (unsigned int j=0; j<numFaces; j++)
- this->A[i][j] = 0;
- }
-
- assembleV(element,k);
- assembleBC(element,k);
- }
-
- // TODO/FIXME: this is only valid for linear problems where
- // the local stiffness matrix is independend of the current
- // solution. We need to implement this properly, but this
- // should at least make the thing compile...
- typedef Dune::FieldVector<Scalar, m> VBlockType;
- void assemble(const Element &cell, const Dune::BlockVector<VBlockType>& localSolution, int orderOfShapeFns = 1)
- {
- assemble(cell, orderOfShapeFns);
- }
-
- //! assemble only boundary conditions for given element
- /*! On exit the following things have been done:
- - The boundary conditions have been evaluated and are accessible with the bc() method
- - The right hand side contains either the value of the essential boundary
- condition or the assembled neumann boundary condition. It is accessible via the rhs() method.
- @param[in] element a codim 0 entity reference
- @param[in] k order of CR basis
- */
- void assembleBoundaryCondition (const Element& element, int k=1)
- {
- unsigned int numFaces = element.template count<1>();
- this->setcurrentsize(numFaces);
-
- // clear assemble data
- for (unsigned int i=0; i<numFaces; i++)
- {
- this->b[i] = 0;
- this->bctype[i].reset();
- }
-
- assembleBC(element,k);
- }
-
- void assembleElementMatrices(const Element& element, Dune::FieldVector<Scalar,2*dim>& faceVol,
- Dune::FieldMatrix<Scalar,2*dim,2*dim>& W, Dune::FieldVector<Scalar,2*dim>& c,
- Dune::FieldMatrix<Scalar,2*dim,2*dim>& Pi, Scalar& dinv, Dune::FieldVector<Scalar,2*dim>& F, Scalar& qmean)
- {
- unsigned int numFaces = element.template count<1>();
- this->setcurrentsize(numFaces);
-
- // get global coordinate of cell center
- const Dune::FieldVector<Scalar,dim>& centerGlobal = element.geometry().center();
-
- int globalIdx = problem_.variables().index(element);
-
- CellData& cellData = problem_.variables().cellData(globalIdx);
-
- // eval diffusion tensor, ASSUMING to be constant over each cell
- Dune::FieldMatrix<Scalar,dim,dim> K(0);
- K = problem_.spatialParams().intrinsicPermeability(element);
-
- K *= (cellData.mobility(wPhaseIdx) + cellData.mobility(nPhaseIdx));
-
- // cell volume
- Scalar volume = element.geometry().volume();
-
- // build the matrices R and ~N
- Dune::FieldMatrix<Scalar,2*dim,dim> R(0), N(0);
-
- // std::cout << "element " << elemId << ": center " << centerGlobal << std::endl;;
-
- typedef typename GridView::IntersectionIterator IntersectionIterator;
-
- int i = -1;
- IntersectionIterator isEndIt = gridView_.iend(element);
- for (IntersectionIterator isIt = gridView_.ibegin(element); isIt!=isEndIt; ++isIt)
- {
- // local number of facet
- i = isIt->indexInInside();
-
- const Dune::FieldVector<Scalar,dim>& faceGlobal = isIt->geometry().center();
- faceVol[i] = isIt->geometry().volume();
-
- // get normal vector
- const Dune::FieldVector<Scalar,dim>& unitOuterNormal = isIt->centerUnitOuterNormal();
-
- N[i] = unitOuterNormal;
-
- for (int k = 0; k < dim; k++)
- // move origin to the center of gravity
- R[i][k] = faceVol[i]*(faceGlobal[k] - centerGlobal[k]);
- }
-
- // std::cout << "N =\dim" << N;
- // std::cout << "R =\dim" << R;
-
- // proceed along the lines of Algorithm 1 from
- // Brezzi/Lipnikov/Simonicini M3AS 2005
- // (1) orthonormalize columns of the matrix R
- Scalar norm = R[0][0]*R[0][0];
- for (unsigned int i = 1; i < numFaces; i++)
- norm += R[i][0]*R[i][0];
- norm = sqrt(norm);
- for (unsigned int i = 0; i < numFaces; i++)
- R[i][0] /= norm;
- Scalar weight = R[0][1]*R[0][0];
- for (unsigned int i = 1; i < numFaces; i++)
- weight += R[i][1]*R[i][0];
- for (unsigned int i = 0; i < numFaces; i++)
- R[i][1] -= weight*R[i][0];
- norm = R[0][1]*R[0][1];
- for (unsigned int i = 1; i < numFaces; i++)
- norm += R[i][1]*R[i][1];
- norm = sqrt(norm);
- for (unsigned int i = 0; i < numFaces; i++)
- R[i][1] /= norm;
- if (dim == 3)
- {
- Scalar weight1 = R[0][2]*R[0][0];
- Scalar weight2 = R[0][2]*R[0][1];
- for (unsigned int i = 1; i < numFaces; i++)
- {
- weight1 += R[i][2]*R[i][0];
- weight2 += R[i][2]*R[i][1];
- }
- for (unsigned int i = 0; i < numFaces; i++)
- R[i][1] -= weight1*R[i][0] + weight2*R[i][1];
- norm = R[0][2]*R[0][2];
- for (unsigned int i = 1; i < numFaces; i++)
- norm += R[i][2]*R[i][2];
- norm = sqrt(norm);
- for (unsigned int i = 0; i < numFaces; i++)
- R[i][2] /= norm;
- }
- // std::cout << "~R =\dim" << R;
-
- // (2) Build the matrix ~D
- Dune::FieldMatrix<Scalar,2*dim,2*dim> D(0);
- for (unsigned int s = 0; s < numFaces; s++)
- {
- Dune::FieldVector<Scalar,2*dim> es(0);
- es[s] = 1;
- for (unsigned int k = 0; k < numFaces; k++)
- {
- D[k][s] = es[k];
- for (int i = 0; i < dim; i++)
- {
- D[k][s] -= R[s][i]*R[k][i];
- }
- }
- }
-
- Scalar traceK = K[0][0];
- for (int i = 1; i < dim; i++)
- traceK += K[i][i];
- D *= 2.0*traceK/volume;
- // std::cout << "u~D =\dim" << D;
-
- // (3) Build the matrix W = Minv
- Dune::FieldMatrix<Scalar,2*dim,dim> NK(N);
- NK.rightmultiply (K);
- for (unsigned int i = 0; i < numFaces; i++)
- {
- for (unsigned int j = 0; j < numFaces; j++)
- {
- W[i][j] = NK[i][0]*N[j][0];
- for (int k = 1; k < dim; k++)
- W[i][j] += NK[i][k]*N[j][k];
- }
- }
-
- W /= volume;
- W += D;
- // std::cout << "W = \dim" << W;
- // std::cout << D[2][2] << ", " << D[2][0] << ", " << D[2][3] << ", " << D[2][1] << std::endl;
- // std::cout << D[0][2] << ", " << D[0][0] << ", " << D[0][3] << ", " << D[0][1] << std::endl;
- // std::cout << D[3][2] << ", " << D[3][0] << ", " << D[3][3] << ", " << D[3][1] << std::endl;
- // std::cout << D[1][2] << ", " << D[1][0] << ", " << D[1][3] << ", " << D[1][1] << std::endl;
-
- // Now the notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
- // The matrix W developed so far corresponds to one element-associated
- // block of the matrix B^{-1} there.
-
- // Corresponding to the element under consideration,
- // calculate the part of the matrix C coupling velocities and element pressures.
- // This is just a row vector of size numFaces.
- // scale with volume
- for (unsigned int i = 0; i < numFaces; i++)
- c[i] = faceVol[i];
-
- // Set up the element part of the matrix \Pi coupling velocities
- // and pressure-traces. This is a diagonal matrix with entries given by faceVol.
- for (unsigned int i = 0; i < numFaces; i++)
- Pi[i][i] = faceVol[i];
-
- // Calculate the element part of the matrix D^{-1} = (c W c^T)^{-1} which is just a scalar value.
- Dune::FieldVector<Scalar,2*dim> Wc(0);
- W.umv(c, Wc);
- dinv = 1.0/(c*Wc);
-
- // Calculate the element part of the matrix F = Pi W c^T which is a column vector.
- F = 0;
- Pi.umv(Wc, F);
- // std::cout << "Pi = \dim" << Pi << "c = " << c << ", F = " << F << std::endl;
-
- PrimaryVariables sourceVec(0.0);
- problem_.source(sourceVec, element);
- qmean = volume*(sourceVec[wPhaseIdx]/density_[wPhaseIdx] + sourceVec[nPhaseIdx]/density_[nPhaseIdx]);
- }
-
-private:
- void assembleV (const Element& element, int k=1)
- {
- unsigned int numFaces = element.template count<1>();
- this->setcurrentsize(numFaces);
-
- // The notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
- // The matrix W developed here corresponds to one element-associated
- // block of the matrix B^{-1} there.
- Dune::FieldVector<Scalar,2*dim> faceVol(0);
- Dune::FieldMatrix<Scalar,2*dim,2*dim> W(0);
- Dune::FieldVector<Scalar,2*dim> c(0);
- Dune::FieldMatrix<Scalar,2*dim,2*dim> Pi(0);
- Dune::FieldVector<Scalar,2*dim> F(0);
- Scalar dinv;
- Scalar qmean;
- this->assembleElementMatrices(element, faceVol, W, c, Pi, dinv, F, qmean);
-
- // Calculate the element part of the matrix Pi W Pi^T.
- Dune::FieldMatrix<Scalar,2*dim,2*dim> PiWPiT(W);
- PiWPiT.rightmultiply(Pi);
- PiWPiT.leftmultiply(Pi);
-
- // Calculate the element part of the matrix F D^{-1} F^T.
- Dune::FieldMatrix<Scalar,2*dim,2*dim> FDinvFT(0);
- for (unsigned int i = 0; i < numFaces; i++)
- for (unsigned int j = 0; j < numFaces; j++)
- FDinvFT[i][j] = dinv*F[i]*F[j];
-
- // Calculate the element part of the matrix S = Pi W Pi^T - F D^{-1} F^T.
- for (unsigned int i = 0; i < numFaces; i++)
- for (unsigned int j = 0; j < numFaces; j++)
- this->A[i][j] = PiWPiT[i][j] - FDinvFT[i][j];
-
- // Calculate the source term F D^{-1} f
- // NOT WORKING AT THE MOMENT
- Scalar factor = dinv*qmean;
- for (unsigned int i = 0; i < numFaces; i++)
- this->b[i] = F[i]*factor;
-
- // std::cout << "faceVol = " << faceVol << std::endl << "W = " << std::endl << W << std::endl
- // << "c = " << c << std::endl << "Pi = " << std::endl << Pi << std::endl
- // << "dinv = " << dinv << std::endl << "F = " << F << std::endl;
- // std::cout << "dinvF = " << dinvF << ", q = " << qmean
- // << ", b = " << this->b[0] << ", " << this->b[1] << ", " << this->b[2] << ", " << this->b[3] << std::endl;
- }
-
- void assembleBC (const Element& element, int k=1)
- {
- // evaluate boundary conditions via intersection iterator
- typedef typename GridView::IntersectionIterator IntersectionIterator;
-
- IntersectionIterator endIsIt = gridView_.iend(element);
- for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != endIsIt; ++isIt)
- {
- int faceIndex = isIt->indexInInside();
- if (isIt->boundary())
- {
- problem_.boundaryTypes(this->bctype[faceIndex], *isIt);
- PrimaryVariables boundValues(0.0);
-
- if (this->bctype[faceIndex].isNeumann(pressEqIdx))
- {
- problem_.neumann(boundValues, *isIt);
- Scalar J = (boundValues[wPhaseIdx]/density_[wPhaseIdx] + boundValues[nPhaseIdx]/density_[nPhaseIdx]);
- this->b[faceIndex] -= J * isIt->geometry().volume();
- }
- else if (this->bctype[faceIndex].isDirichlet(pressEqIdx))
- {
- problem_.dirichlet(boundValues, *isIt);
- this->b[faceIndex] = boundValues[pressureIdx];
- }
- }
- }
- }
-
-private:
- Problem& problem_;
- const GridView gridView_;
-
- Scalar density_[numPhases];
-};
-
-/** @} */
-}
-#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator.hh
deleted file mode 100644
index 7646d1a815..0000000000
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator.hh
+++ /dev/null
@@ -1,126 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program 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 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief An assembler for the Jacobian matrix based on mimetic FD.
- */
-#ifndef DUMUX_MIMETICOPERATOR_HH
-#define DUMUX_MIMETICOPERATOR_HH
-
-#include<iostream>
-#include<map>
-
-#include<dune/common/timer.hh>
-#include<dune/common/fvector.hh>
-#include<dune/common/fmatrix.hh>
-#include<dune/geometry/type.hh>
-#include<dune/grid/common/grid.hh>
-#include<dune/grid/common/mcmgmapper.hh>
-#include<dune/istl/bvector.hh>
-#include<dune/istl/operators.hh>
-#include<dune/istl/bcrsmatrix.hh>
-
-#include<dumux/common/boundaryconditions.hh>
-#include"croperator.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup MimeticPressure2p
- * @brief Levelwise assembler
-
- This class serves as a base class for local assemblers. It provides
- space and access to the local stiffness matrix. The actual assembling is done
- in a derived class via the virtual assemble method.
-
- \tparam Scalar The field type used in the elements of the stiffness matrix
- \tparam GridView The grid view of the simulation grid
-*/
-template<class Scalar, class GridView>
-class MimeticOperatorAssembler : public CROperatorAssembler<Scalar, GridView>
-{
- template<int dim>
- struct ElementLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim;
- }
- };
-
- enum {dim=GridView::dimension};
- typedef Dune::BlockVector< Dune::FieldVector<Scalar,2*dim> > VType;
- typedef Dune::BlockVector< Dune::FieldVector<Scalar,1> > PType;
- typedef typename GridView::template Codim<0>::Iterator Iterator;
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,ElementLayout> ElementMapper;
-
-public:
-
- MimeticOperatorAssembler (const GridView& gridView)
- : CROperatorAssembler<Scalar, GridView>(gridView), elementMapper(gridView)
- {}
-
- template<class LocalStiffness, class Vector>
- void calculatePressure (LocalStiffness& loc, Vector& u,
- VType& velocity, PType& pressure)
- {
- // run over all level elements
- Iterator eEndIt = this->gridView_.template end<0>();
- for (Iterator eIt = this->gridView_.template begin<0>(); eIt!=eEndIt; ++eIt)
- {
- unsigned int numFaces = eIt->template count<1>();
-
- int elemId = elementMapper.map(*eIt);
-
- // get local to global id map and pressure traces
- Dune::FieldVector<Scalar,2*dim> pressTrace(0);
- for (unsigned int k = 0; k < numFaces; k++)
- {
- pressTrace[k] = u[this->faceMapper_.map(*eIt, k, 1)];
- }
-
- // The notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
- // The matrix W developed here corresponds to one element-associated
- // block of the matrix B^{-1} there.
- Dune::FieldVector<Scalar,2*dim> faceVol(0);
- Dune::FieldMatrix<Scalar,2*dim,2*dim> W(0);
- Dune::FieldVector<Scalar,2*dim> c(0);
- Dune::FieldMatrix<Scalar,2*dim,2*dim> Pi(0);
- Dune::FieldVector<Scalar,2*dim> F(0);
- Scalar dinv = 0;
- Scalar qmean = 0;
- loc.assembleElementMatrices(*eIt, faceVol, W, c, Pi, dinv, F, qmean);
-
- pressure[elemId] = dinv*(qmean + (F*pressTrace));
-
- Dune::FieldVector<Scalar,2*dim> v(0);
- for (int i = 0; i < 2*dim; i++)
- for (int j = 0; j < 2*dim; j++)
- v[i] += W[i][j]*faceVol[j]*(pressure[elemId] - pressTrace[j]);
-
- velocity[elemId] = v;
- }
- }
-
-private:
- ElementMapper elementMapper;
-};
-}
-#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2p.hh
new file mode 100644
index 0000000000..2a176b74a0
--- /dev/null
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2p.hh
@@ -0,0 +1,300 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Markus Wolff *
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute eIt and/or modify *
+ * eIt under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that eIt 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief An assembler for the Jacobian matrix based on mimetic FD.
+ */
+#ifndef DUMUX_MIMETICOPERATOR2P_HH
+#define DUMUX_MIMETICOPERATOR2P_HH
+
+#include"croperator2p.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mimetic2P
+ * @brief Levelwise assembler
+
+ This class serves as a base class for local assemblers. It provides
+ space and access to the local stiffness matrix. The actual assembling is done
+ in a derived class via the virtual assemble method.
+
+ The template parameters are:
+
+ - Scalar The field type used in the elements of the stiffness matrix
+ */
+template<class TypeTag>
+class MimeticOperatorAssemblerTwoP: public CROperatorAssemblerTwoP<TypeTag>
+{
+ typedef CROperatorAssemblerTwoP<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld,
+ };
+ typedef typename GET_PROP_TYPE(TypeTag, LocalStiffness) LocalStiffness;
+
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNw,
+ pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation),
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ saturationIdx = Indices::saturationIdx,
+ satEqIdx = Indices::satEqIdx,
+ pressEqIdx = Indices::pressEqIdx
+ };
+
+ typedef Dune::FieldVector<Scalar, dimWorld> FieldVector;
+
+public:
+
+ MimeticOperatorAssemblerTwoP(const GridView& gridView) :
+ ParentType(gridView)
+ {
+ }
+
+ template<class Vector>
+ void calculatePressure(LocalStiffness& loc, Vector& u, Problem& problem)
+ {
+ Dune::FieldVector<Scalar, 2 * dim> velocityW(0);
+ Dune::FieldVector<Scalar, 2 * dim> velocityNw(0);
+ Dune::FieldVector<Scalar, 2 * dim> pressTrace(0);
+ Dune::FieldVector<Scalar, 2 * dim> gravPotTrace(0);
+
+ // run over all level elements
+ ElementIterator eIt = this->gridView_.template begin<0>();
+ ElementIterator eItEnd = this->gridView_.template end<0>();
+
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem.referencePressure(*eIt));
+ fluidState.setPressure(nPhaseIdx, problem.referencePressure(*eIt));
+ fluidState.setTemperature(problem.temperature(*eIt));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ Scalar densityDiff = FluidSystem::density(fluidState, nPhaseIdx) - FluidSystem::density(fluidState, wPhaseIdx);
+ Scalar viscosityW = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ Scalar viscosityNw = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
+ //reset velocity
+ for (int i = 0; i < problem.gridView().size(0); i++)
+ {
+ problem.variables().cellData(i).fluxData().resetVelocity();
+ }
+
+ for (; eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = problem.variables().index(*eIt);
+
+ CellData& cellData = problem.variables().cellData(globalIdx);
+ FieldVector globalPos = eIt->geometry().center();
+
+ // get local to global id map and pressure traces
+ IntersectionIterator isIt = problem.gridView().template ibegin(*eIt);
+ const IntersectionIterator &isItEnd = problem.gridView().template iend(*eIt);
+ for (; isIt != isItEnd; ++isIt)
+ {
+ int indexInInside = isIt->indexInInside();
+
+ int globalIdxFace = this->faceMapper_.map(*eIt, indexInInside, 1);
+
+ pressTrace[indexInInside] = u[globalIdxFace];
+
+ gravPotTrace[indexInInside] = (problem.bBoxMax() - isIt->geometry().center()) * problem.gravity() * densityDiff;
+ }
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ Scalar potW = loc.constructPressure(*eIt, pressTrace);
+ Scalar gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * densityDiff;
+ Scalar potNw = potW + gravPot;
+
+ cellData.setPotential(wPhaseIdx, potW);
+ cellData.setPotential(nPhaseIdx, potNw);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, wPhaseIdx);
+
+ cellData.setPressure(wPhaseIdx, potW - gravPot);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, nPhaseIdx);
+
+ cellData.setPressure(nPhaseIdx, potNw - gravPot);
+
+ break;
+ }
+ case pn:
+ {
+ Scalar potNw = loc.constructPressure(*eIt, pressTrace);
+ Scalar gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * densityDiff;
+ Scalar potW = potNw - gravPot;
+
+ cellData.setPotential(nPhaseIdx, potNw);
+ cellData.setPotential(wPhaseIdx, potW);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, wPhaseIdx);
+
+ cellData.setPressure(wPhaseIdx, potW - gravPot);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, nPhaseIdx);
+
+ cellData.setPressure(nPhaseIdx, potNw - gravPot);
+
+ break;
+ }
+ }
+
+ //velocity reconstruction: !!! The velocity which is not reconstructed from the primary pressure variable can be slightly wrong and not conservative!!!!
+ // -> Should not be used for transport!!
+ switch (pressureType)
+ {
+ case pw:
+ {
+ loc.constructVelocity(*eIt, velocityW, pressTrace, cellData.potential(wPhaseIdx));
+ pressTrace += gravPotTrace;
+ loc.constructVelocity(*eIt, velocityNw, pressTrace, cellData.potential(nPhaseIdx));
+
+ break;
+ }
+ case pn:
+ {
+ loc.constructVelocity(*eIt, velocityW, pressTrace, cellData.potential(nPhaseIdx));
+ pressTrace -= gravPotTrace;
+ loc.constructVelocity(*eIt, velocityNw, pressTrace, cellData.potential(wPhaseIdx));
+
+ break;
+ }
+ }
+
+ isIt = problem.gridView().template ibegin(*eIt);
+ for (; isIt != isItEnd; ++isIt)
+ {
+ int idxInInside = isIt->indexInInside();
+
+ cellData.fluxData().setUpwindPotential(wPhaseIdx, idxInInside, velocityW[idxInInside]);
+ cellData.fluxData().setUpwindPotential(nPhaseIdx, idxInInside, velocityNw[idxInInside]);
+
+ Scalar mobilityW = 0;
+ Scalar mobilityNw = 0;
+
+ if (isIt->neighbor())
+ {
+ int neighborIdx = problem.variables().index(*(isIt->outside()));
+
+ CellData& cellDataNeighbor = problem.variables().cellData(neighborIdx);
+
+ mobilityW =
+ (velocityW[idxInInside] >= 0.) ? cellData.mobility(wPhaseIdx) :
+ cellDataNeighbor.mobility(wPhaseIdx);
+ mobilityNw =
+ (velocityNw[idxInInside] >= 0.) ? cellData.mobility(nPhaseIdx) :
+ cellDataNeighbor.mobility(nPhaseIdx);
+
+ int idxInOutside = isIt->indexInOutside();
+
+ if (velocityW[idxInInside] >= 0.)
+ {
+ FieldVector velocity(isIt->centerUnitOuterNormal());
+ velocity *= mobilityW/(mobilityW+mobilityNw) * velocityW[idxInInside];
+ cellData.fluxData().addVelocity(wPhaseIdx, idxInInside, velocity);
+ cellDataNeighbor.fluxData().addVelocity(wPhaseIdx, isIt->indexInOutside(), velocity);
+ }
+ if (velocityNw[idxInInside] >= 0.)
+ {
+ FieldVector velocity(isIt->centerUnitOuterNormal());
+ velocity *= mobilityNw/(mobilityW+mobilityNw) * velocityNw[idxInInside];
+ cellData.fluxData().addVelocity(nPhaseIdx, idxInInside, velocity);
+ cellDataNeighbor.fluxData().addVelocity(nPhaseIdx, isIt->indexInOutside(), velocity);
+ }
+
+ cellData.fluxData().setVelocityMarker(idxInInside);
+ }
+ else
+ {
+ BoundaryTypes bctype;
+ problem.boundaryTypes(bctype, *isIt);
+ if (bctype.isDirichlet(satEqIdx))
+ {
+ PrimaryVariables boundValues(0.0);
+ problem.dirichlet(boundValues, *isIt);
+
+ if (velocityW[idxInInside] >= 0.)
+ {
+ mobilityW = cellData.mobility(wPhaseIdx);
+ }
+ else
+ {
+ mobilityW = MaterialLaw::krw(problem.spatialParams().materialLawParams(*eIt),
+ boundValues[saturationIdx]) / viscosityW;
+ }
+
+ if (velocityNw[idxInInside] >= 0.)
+ {
+ mobilityNw = cellData.mobility(nPhaseIdx);
+ }
+ else
+ {
+ mobilityNw = MaterialLaw::krn(problem.spatialParams().materialLawParams(*eIt),
+ boundValues[saturationIdx]) / viscosityNw;
+ }
+ }
+ else
+ {
+ mobilityW = cellData.mobility(wPhaseIdx);
+ mobilityNw = cellData.mobility(nPhaseIdx);
+ }
+
+ FieldVector velocity(isIt->centerUnitOuterNormal());
+ velocity *= mobilityW / (mobilityW + mobilityNw) * velocityW[idxInInside];
+ cellData.fluxData().setVelocity(wPhaseIdx, idxInInside, velocity);
+
+ velocity = 0;
+ velocity = isIt->centerUnitOuterNormal();
+ velocity *= mobilityNw / (mobilityW + mobilityNw) * velocityNw[idxInInside];
+ cellData.fluxData().setVelocity(nPhaseIdx, idxInInside, velocity);
+ cellData.fluxData().setVelocityMarker(idxInInside);
+ }
+ }
+ }
+ }
+};
+}
+#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2padaptive.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2padaptive.hh
new file mode 100644
index 0000000000..8819201010
--- /dev/null
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2padaptive.hh
@@ -0,0 +1,308 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Markus Wolff *
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief An assembler for the Jacobian matrix based on mimetic FD.
+ */
+#ifndef DUMUX_MIMETICOPERATOR2PADAPTIVE_HH
+#define DUMUX_MIMETICOPERATOR2PADAPTIVE_HH
+
+#include"croperator2padaptive.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mimetic2P
+ * @brief Levelwise assembler
+
+ This class serves as a base class for local assemblers. It provides
+ space and access to the local stiffness matrix. The actual assembling is done
+ in a derived class via the virtual assemble method.
+
+ The template parameters are:
+
+ - Scalar The field type used in the elements of the stiffness matrix
+*/
+template<class TypeTag>
+class MimeticOperatorAssemblerTwoPAdaptive : public CROperatorAssemblerTwoPAdaptive<TypeTag>
+{
+ typedef CROperatorAssemblerTwoPAdaptive<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ enum
+ {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld
+ };
+ typedef typename GET_PROP_TYPE(TypeTag, LocalStiffness) LocalStiffness;
+
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNw,
+ pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation),
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ saturationIdx = Indices::saturationIdx,
+ satEqIdx = Indices::satEqIdx
+ };
+
+ typedef Dune::FieldVector<Scalar, dimWorld> FieldVector;
+public:
+
+ MimeticOperatorAssemblerTwoPAdaptive (const GridView& gridView)
+ : ParentType(gridView)
+ {}
+
+ template<class Vector>
+ void calculatePressure(LocalStiffness& loc, Vector& u, Problem& problem)
+ {
+ Dune::DynamicVector<Scalar> velocityW(2*dim);
+ Dune::DynamicVector<Scalar> velocityNw(2*dim);
+ Dune::DynamicVector<Scalar> pressTraceW(2*dim);
+ Dune::DynamicVector<Scalar> pressTraceNw(2*dim);
+
+ // run over all level elements
+ ElementIterator eIt = this->gridView_.template begin<0>();
+ ElementIterator eItEnd = this->gridView_.template end<0>();
+
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem.referencePressure(*eIt));
+ fluidState.setPressure(nPhaseIdx, problem.referencePressure(*eIt));
+ fluidState.setTemperature(problem.temperature(*eIt));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ Scalar densityDiff = FluidSystem::density(fluidState, nPhaseIdx) - FluidSystem::density(fluidState, wPhaseIdx);
+ Scalar viscosityW = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ Scalar viscosityNw = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
+ //reset velocity
+ for (int i = 0; i < problem.gridView().size(0); i++)
+ {
+ problem.variables().cellData(i).fluxData().resetVelocity();
+ }
+
+ for (; eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = problem.variables().index(*eIt);
+
+ unsigned int numFaces = this->intersectionMapper_.size(globalIdx);
+
+ // get local to global id map and pressure traces
+ velocityW.resize(numFaces);
+ velocityNw.resize(numFaces);
+ pressTraceW.resize(numFaces);
+ pressTraceNw.resize(numFaces);
+
+ CellData& cellData = problem.variables().cellData(globalIdx);
+ FieldVector globalPos = eIt->geometry().center();
+
+ int intersectionIdx = -1;
+ // get local to global id map and pressure traces
+ IntersectionIterator isIt = problem.gridView().template ibegin(*eIt);
+ const IntersectionIterator &isItEnd = problem.gridView().template iend(*eIt);
+ for (; isIt != isItEnd; ++isIt)
+ {
+ ++intersectionIdx;
+
+ int globalIdxFace = this->intersectionMapper_.map(*eIt, intersectionIdx);
+
+ Scalar pcPotFace = (problem.bBoxMax() - isIt->geometry().center()) * problem.gravity() * densityDiff;
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressTraceW[intersectionIdx] = u[globalIdxFace];
+ pressTraceNw[intersectionIdx] = u[globalIdxFace] + pcPotFace;
+ break;
+ }
+ case pn:
+ {
+ pressTraceNw[intersectionIdx] = u[globalIdxFace];
+ pressTraceW[intersectionIdx] = u[globalIdxFace] - pcPotFace;
+
+ break;
+ }
+ }
+ }
+
+ switch (pressureType)
+ {
+ case pw:
+ {
+ Scalar potW = loc.constructPressure(*eIt, pressTraceW);
+ Scalar gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * densityDiff;
+ Scalar potNw = potW + gravPot;
+
+ cellData.setPotential(wPhaseIdx, potW);
+ cellData.setPotential(nPhaseIdx, potNw);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, wPhaseIdx);
+
+ cellData.setPressure(wPhaseIdx, potW - gravPot);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, nPhaseIdx);
+
+ cellData.setPressure(nPhaseIdx, potNw - gravPot);
+
+ break;
+ }
+ case pn:
+ {
+ Scalar potNw = loc.constructPressure(*eIt, pressTraceNw);
+ Scalar gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * densityDiff;
+ Scalar potW = potNw - gravPot;
+
+ cellData.setPotential(nPhaseIdx, potNw);
+ cellData.setPotential(wPhaseIdx, potW);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, wPhaseIdx);
+
+ cellData.setPressure(wPhaseIdx, potW - gravPot);
+
+ gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, nPhaseIdx);
+
+ cellData.setPressure(nPhaseIdx, potNw - gravPot);
+
+ break;
+ }
+ }
+
+ //velocity reconstruction: !!! The velocity which is not reconstructed from the primary pressure variable can be slightly wrong and not conservative!!!!
+ // -> Should not be used for transport!!
+ loc.constructVelocity(*eIt, velocityW, pressTraceW, cellData.potential(wPhaseIdx));
+ loc.constructVelocity(*eIt, velocityNw, pressTraceNw, cellData.potential(nPhaseIdx));
+
+ intersectionIdx = -1;
+ isIt = problem.gridView().template ibegin(*eIt);
+ for (; isIt != isItEnd; ++isIt)
+ {
+ ++intersectionIdx;
+ int idxInInside = isIt->indexInInside();
+
+ cellData.fluxData().addUpwindPotential(wPhaseIdx, idxInInside, velocityW[intersectionIdx]);
+ cellData.fluxData().addUpwindPotential(nPhaseIdx, idxInInside, velocityNw[intersectionIdx]);
+
+ Scalar mobilityW = 0;
+ Scalar mobilityNw = 0;
+
+ if (isIt->neighbor())
+ {
+ int neighborIdx = problem.variables().index(*(isIt->outside()));
+
+ CellData& cellDataNeighbor = problem.variables().cellData(neighborIdx);
+
+ mobilityW =
+ (velocityW[intersectionIdx] >= 0.) ? cellData.mobility(wPhaseIdx) :
+ cellDataNeighbor.mobility(wPhaseIdx);
+ mobilityNw =
+ (velocityNw[intersectionIdx] >= 0.) ? cellData.mobility(nPhaseIdx) :
+ cellDataNeighbor.mobility(nPhaseIdx);
+
+ int idxInOutside = isIt->indexInOutside();
+
+ if (velocityW[intersectionIdx] >= 0.)
+ {
+ FieldVector velocity(isIt->centerUnitOuterNormal());
+ velocity *= mobilityW/(mobilityW+mobilityNw) * velocityW[intersectionIdx];
+ cellData.fluxData().addVelocity(wPhaseIdx, idxInInside, velocity);
+ cellDataNeighbor.fluxData().addVelocity(wPhaseIdx, isIt->indexInOutside(), velocity);
+ }
+ if (velocityNw[intersectionIdx] >= 0.)
+ {
+ FieldVector velocity(isIt->centerUnitOuterNormal());
+ velocity *= mobilityNw/(mobilityW+mobilityNw) * velocityNw[intersectionIdx];
+ cellData.fluxData().addVelocity(nPhaseIdx, idxInInside, velocity);
+ cellDataNeighbor.fluxData().addVelocity(nPhaseIdx, isIt->indexInOutside(), velocity);
+ }
+
+ cellData.fluxData().setVelocityMarker(idxInInside);
+ }
+ else
+ {
+ BoundaryTypes bctype;
+ problem.boundaryTypes(bctype, *isIt);
+ if (bctype.isDirichlet(satEqIdx))
+ {
+ PrimaryVariables boundValues(0.0);
+ problem.dirichlet(boundValues, *isIt);
+
+ if (velocityW[intersectionIdx] >= 0.)
+ {
+ mobilityW = cellData.mobility(wPhaseIdx);
+ }
+ else
+ {
+ mobilityW = MaterialLaw::krw(problem.spatialParams().materialLawParams(*eIt),
+ boundValues[saturationIdx]) / viscosityW;
+ }
+
+ if (velocityNw[intersectionIdx] >= 0.)
+ {
+ mobilityNw = cellData.mobility(nPhaseIdx);
+ }
+ else
+ {
+ mobilityNw = MaterialLaw::krn(problem.spatialParams().materialLawParams(*eIt),
+ boundValues[saturationIdx]) / viscosityNw;
+ }
+ }
+ else
+ {
+ mobilityW = cellData.mobility(wPhaseIdx);
+ mobilityNw = cellData.mobility(nPhaseIdx);
+ }
+
+ FieldVector velocity(isIt->centerUnitOuterNormal());
+ velocity *= mobilityW/(mobilityW+mobilityNw) * velocityW[intersectionIdx];
+ cellData.fluxData().addVelocity(wPhaseIdx, idxInInside, velocity);
+
+
+ velocity = isIt->centerUnitOuterNormal();
+ velocity *= mobilityNw/(mobilityW+mobilityNw) * velocityNw[intersectionIdx];
+ cellData.fluxData().addVelocity(nPhaseIdx, idxInInside, velocity);
+ cellData.fluxData().setVelocityMarker(idxInInside);
+ }
+ }
+ }
+ }
+};
+}
+#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
index bcd2c99920..3e519d7520 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
@@ -1,7 +1,9 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
- * See the file COPYING for full copying permissions. *
+ * Copyright (C) 2011 by Markus Wolff *
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
@@ -10,7 +12,7 @@
* *
* This program 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 *
+ * 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 *
@@ -24,48 +26,40 @@
#ifndef DUMUX_MIMETICPRESSURE2P_HH
#define DUMUX_MIMETICPRESSURE2P_HH
-// dune environent:
-#include <dune/istl/bvector.hh>
-#include <dune/istl/operators.hh>
-#include <dune/istl/solvers.hh>
-#include <dune/istl/preconditioners.hh>
-
// dumux environment
-#include <dumux/decoupled/2p/diffusion/diffusionproperties2p.hh>
-#include <dumux/decoupled/common/mimetic/mimeticproperties.hh>
-#include <dumux/decoupled/2p/diffusion/mimetic/mimeticoperator.hh>
+#include "dumux/decoupled/common/mimetic/mimeticproperties.hh"
+#include "dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2p.hh"
+#include "dumux/decoupled/2p/diffusion/mimetic/mimetic2p.hh"
namespace Dumux
{
-//! \ingroup MimeticPressure2p
-/*! \brief Mimetic finite differences discretization of a two-phase pressure equation of the sequential IMPES model.
- *
- * This class provides a mimetic finite differences implementation for solving equations of the form
- * \f[
- * \text{div}\, \boldsymbol{v}_{total} = q.
- * \f]
- * The total velocity \f$\boldsymbol{v}_{total}\f$ is defined using a global pressure approach. This leads to
- * \f[
- * - \text{div}\, \left(\lambda \boldsymbol K \textbf{grad}\, p_{global}\right) = q.
- * \f]
- * Here, \f$ p_{global} \f$ is the global pressure of a classical fractional flow formulation
- * (see e.g. P. Binning and M. A. Celia, “Practical implementation of the fractional flow approach to multi-phase flow simulation,†Advances in water resources, vol. 22, no. 5, pp. 461-478, 1999.),
- * \f$ \boldsymbol K \f$ the absolute permeability, \f$ \lambda = \lambda_w + \lambda_n \f$ the total mobility depending on the
- * saturation (\f$ \lambda_\alpha = k_{r_\alpha} / \mu_\alpha \f$) and \f$ q \f$ the source term. Gravity is neglected in this implementation.
+
+/*! \ingroup Mimetic2p
*
- * \f$ p = p_D \f$ on \f$ \Gamma_{Dirichlet} \f$, and \f$ \boldsymbol v_{total} \cdot \boldsymbol n = q_N \f$ on \f$ \Gamma_{Neumann} \f$.
+ * \brief mimetic method for the pressure equation
*
- * Remark: gravity is neglected!
+ * Provides a mimetic implementation for the evaluation
+ * of equations of the form
+ * \f[\text{div}\, \boldsymbol{v}_{total} = q.\f]
+ * The definition of the total velocity \f$\boldsymbol{v}_total\f$ depends on the kind of pressure chosen. This could be a wetting (w) phase pressure leading to
+ * \f[ - \text{div}\, \left[\lambda \boldsymbol{K} \left(\text{grad}\, p_w + f_n \text{grad}\, p_c + \sum f_\alpha \rho_\alpha g \text{grad}\, z\right)\right] = q, \f]
+ * a non-wetting (n) phase pressure yielding
+ * \f[ - \text{div}\, \left[\lambda \boldsymbol{K} \left(\text{grad}\, p_n - f_w \text{grad}\, p_c + \sum f_\alpha \rho_\alpha g \text{grad}\, z\right)\right] = q, \f]
+ * or a global pressure leading to
+ * \f[ - \text{div}\, \left[\lambda \boldsymbol{K} \left(\text{grad}\, p_{global} + \sum f_\alpha \rho_\alpha g \text{grad}\, z\right)\right] = q.\f]
+ * Here, \f$p\f$ denotes a pressure, \f$\boldsymbol{K}\f$ the absolute permeability, \f$\lambda\f$ the total mobility, possibly depending on the
+ * saturation,\f$f\f$ the fractional flow function of a phase, \f$\rho\f$ a phase density, \f$g\f$ the gravity constant and \f$q\f$ the source term.
+ * For all cases, \f$p = p_D\f$ on \f$\Gamma_{Neumann}\f$, and \f$\boldsymbol{v}_{total} = q_N\f$
+ * on \f$\Gamma_{Dirichlet}\f$.
*
- * \tparam TypeTag The Type Tag
+ *\tparam TypeTag The Type Tag
*/
template<class TypeTag> class MimeticPressure2P
{
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Variables) Variables;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
typedef typename SpatialParams::MaterialLaw MaterialLaw;
@@ -75,16 +69,21 @@ template<class TypeTag> class MimeticPressure2P
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
- typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
-
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
};
enum
{
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
pGlobal = Indices::pressureGlobal,
- sw = Indices::saturationW
+ sw = Indices::saturationW,
+ sn = Indices::saturationNW,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation), //!< gives kind of pressure used (\f$ 0 = p_w\f$, \f$ 1 = p_n\f$, \f$ 2 = p_{global}\f$)
+ saturationType = GET_PROP_VALUE(TypeTag, SaturationFormulation), //!< gives kind of saturation used (\f$ 0 = S_w\f$, \f$ 1 = S_n\f$)
};
enum
{
@@ -95,20 +94,25 @@ template<class TypeTag> class MimeticPressure2P
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Grid Grid;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
+
typedef typename GET_PROP_TYPE(TypeTag, LocalStiffness) LocalStiffness;
typedef Dune::BlockVector< Dune::FieldVector<Scalar, 1> > TraceType;
- typedef Dune::BlockVector< Dune::FieldVector<Scalar, 2*dim> > NormalVelType;
- typedef MimeticOperatorAssembler<Scalar,GridView> OperatorAssembler;
+ typedef MimeticOperatorAssemblerTwoP<TypeTag> OperatorAssembler;
- ///@cond false
- typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolution;
- ///@endcond
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) Matrix;
- typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) RHSVector;
- typedef typename GET_PROP_TYPE(TypeTag, PressureSolutionVector) PressureSolution;
+ typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) Vector;
+
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElementContainer;
+ typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
//initializes the matrix to store the system of equations
void initializeMatrix();
@@ -116,9 +120,33 @@ template<class TypeTag> class MimeticPressure2P
//function which assembles the system of equations to be solved
void assemble(bool first)
{
- LocalStiffness lstiff(problem_, false, problem_.gridView());
- lstiff.initialize();
- A_.assemble(lstiff, pressTrace_, f_);
+ Scalar timeStep = problem_.timeManager().timeStepSize();
+ Scalar maxError = 0.0;
+ int size = problem_.gridView().size(0);
+ for (int i = 0; i < size; i++)
+ {
+ Scalar sat = 0;
+ switch (saturationType)
+ {
+ case sw:
+ sat = problem_.variables().cellData(i).saturation(wPhaseIdx);
+ break;
+ case sn:
+ sat = problem_.variables().cellData(i).saturation(nPhaseIdx);
+ break;
+ }
+ if (sat > 1.0)
+ {
+ maxError = std::max(maxError, (sat - 1.0) / timeStep);
+ }
+ if (sat < 0.0)
+ {
+ maxError = std::max(maxError, (-sat) / timeStep);
+ }
+ }
+
+ lstiff_.setErrorInfo(maxError, timeStep);
+ A_.assemble(lstiff_, pressTrace_, f_);
return;
}
@@ -127,36 +155,15 @@ template<class TypeTag> class MimeticPressure2P
void postprocess()
{
- LocalStiffness lstiff(problem_, false, problem_.gridView());
- lstiff.initialize();
- A_.calculatePressure(lstiff, pressTrace_, normalVelocity_, pressure_);
- return;
- }
-
-protected:
- //! \cond \private
- Problem& problem()
- {
- return problem_;
- }
+ A_.calculatePressure(lstiff_, pressTrace_, problem_);
- const Problem& problem() const
- {
- return problem_;
+ return;
}
- //! \endcond
public:
//constitutive functions are initialized and stored in the variables object
void updateMaterialLaws();
-
- /*! \brief Initializes the pressure model
- *
- * Initializes pressure and velocity field.
- *
- * \param solveTwice indicates if more than one iteration is allowed to get an initial pressure solution
- */
void initialize(bool solveTwice = true)
{
ElementIterator element = problem_.gridView().template begin<0> ();
@@ -172,140 +179,238 @@ public:
viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
updateMaterialLaws();
- A_.initializeMatrix();
- f_.resize(problem_.gridView().size(1));//resize to make sure the final grid size (after the problem was completely built) is used!
- pressure_.resize(problem_.gridView().size(0));
+ A_.initialize();
pressTrace_.resize(problem_.gridView().size(1));
- normalVelocity_.resize(problem_.gridView().size(0)),
- pressure_ = 0;
- pressTrace_ = 0;
+ f_.resize(problem_.gridView().size(1));
+ lstiff_.initialize();
+ lstiff_.reset();
+
+ pressTrace_ = 0.0;
f_ = 0;
- normalVelocity_ = 0.0;
+
assemble(true);
solve();
postprocess();
- storePressureSolution();
- storeVelocity();
-
return;
}
- /*! \brief Pressure and velocity update */
+ void updateVelocity()
+ {
+ updateMaterialLaws();
+ postprocess();
+ }
+
void update()
{
+ lstiff_.reset();
+
assemble(false);
solve();
postprocess();
- storePressureSolution();
- storeVelocity();
-
return;
}
- /*! \brief Globally stores the pressure solution*/
- void storePressureSolution()
+ //! \brief Write data files
+ /* \param name file name */
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
{
int size = problem_.gridView().size(0);
- for (int i = 0; i < size; i++)
+ ScalarSolutionType *potential = writer.allocateManagedBuffer(size);
+ ScalarSolutionType *pressure = 0;
+ ScalarSolutionType *pressureSecond = 0;
+ ScalarSolutionType *potentialSecond = 0;
+ Dune::BlockVector < DimVector > *velocityWetting = 0;
+ Dune::BlockVector < DimVector > *velocityNonwetting = 0;
+
+ if (vtkOutputLevel_ > 0)
{
- CellData& cellData = problem_.variables().cellData(i);
- storePressureSolution(i, cellData);
- cellData.fluxData().resetVelocity();
+ pressure = writer.allocateManagedBuffer(size);
+ pressureSecond = writer.allocateManagedBuffer(size);
+ potentialSecond = writer.allocateManagedBuffer(size);
+ velocityWetting = writer.template allocateManagedBuffer<Scalar, dim>(size);
+ velocityNonwetting = writer.template allocateManagedBuffer<Scalar, dim>(size);
}
- }
- /*! \brief Stores the pressure solution of a cell
- *
- * \param globalIdx Global cell index
- * \param cellData A CellData object
- */
- void storePressureSolution(int globalIdx, CellData& cellData)
- {
- cellData.setGlobalPressure(pressure_[globalIdx]);
- }
-
- // Globally stores the velocity solution at cell-cell interfaces
- void storeVelocity();
+ ElementIterator eItBegin = problem_.gridView().template begin<0>();
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = eItBegin; eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = problem_.variables().index(*eIt);
+ CellData& cellData = problem_.variables().cellData(globalIdx);
+
+ if (pressureType == pw)
+ {
+ (*potential)[globalIdx] = cellData.potential(wPhaseIdx);
+ }
+
+ if (pressureType == pn)
+ {
+ (*potential)[globalIdx] = cellData.potential(nPhaseIdx);
+ }
+
+ if (vtkOutputLevel_ > 0)
+ {
+
+ if (pressureType == pw)
+ {
+ (*pressure)[globalIdx] = cellData.pressure(wPhaseIdx);
+ (*potentialSecond)[globalIdx] = cellData.potential(nPhaseIdx);
+ (*pressureSecond)[globalIdx] = cellData.pressure(nPhaseIdx);
+ }
+
+ if (pressureType == pn)
+ {
+ (*pressure)[globalIdx] = cellData.pressure(nPhaseIdx);
+ (*potentialSecond)[globalIdx] = cellData.potential(wPhaseIdx);
+ (*pressureSecond)[globalIdx] = cellData.pressure(wPhaseIdx);
+ }
+
+ Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
+ Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
+
+ // run through all intersections with neighbors and boundary
+ IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ int isIndex = isIt->indexInInside();
+
+ fluxW[isIndex] += isIt->geometry().volume()
+ * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
+ fluxNw[isIndex] += isIt->geometry().volume()
+ * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
+ }
+
+ DimVector refVelocity(0);
+ refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
+ refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+ refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
+
+ const DimVector& localPos = ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
+
+ // get the transposed Jacobian of the element mapping
+ const DimMatrix& jacobianT = eIt->geometry().jacobianTransposed(localPos);
+
+ // calculate the element velocity by the Piola transformation
+ DimVector elementVelocity(0);
+ jacobianT.umtv(refVelocity, elementVelocity);
+ elementVelocity /= eIt->geometry().integrationElement(localPos);
+
+ (*velocityWetting)[globalIdx] = elementVelocity;
+
+ refVelocity = 0;
+ refVelocity[0] = 0.5 * (fluxNw[1] - fluxNw[0]);
+ refVelocity[1] = 0.5 * (fluxNw[3] - fluxNw[2]);
+ refVelocity[2] = 0.5 * (fluxNw[5] - fluxNw[4]);
+
+ // calculate the element velocity by the Piola transformation
+ elementVelocity = 0;
+ jacobianT.umtv(refVelocity, elementVelocity);
+ elementVelocity /= eIt->geometry().integrationElement(localPos);
+
+ (*velocityNonwetting)[globalIdx] = elementVelocity;
+ }
+ }
+
+ if (pressureType == pw)
+ {
+ writer.attachCellData(*potential, "wetting potential");
+ }
+
+ if (pressureType == pn)
+ {
+ writer.attachCellData(*potential, "nonwetting potential");
+ }
+
+ if (vtkOutputLevel_ > 0)
+ {
+ if (pressureType == pw)
+ {
+ writer.attachCellData(*pressure, "wetting pressure");
+ writer.attachCellData(*pressureSecond, "nonwetting pressure");
+ writer.attachCellData(*potentialSecond, "nonwetting potential");
+ }
+
+ if (pressureType == pn)
+ {
+ writer.attachCellData(*pressure, "nonwetting pressure");
+ writer.attachCellData(*pressureSecond, "wetting pressure");
+ writer.attachCellData(*potentialSecond, "wetting potential");
+ }
+
+ writer.attachCellData(*velocityWetting, "wetting-velocity", dim);
+ writer.attachCellData(*velocityNonwetting, "non-wetting-velocity", dim);
+ }
+ }
- /*! \brief Function for serialization of the pressure field.
- *
- * Function needed for restart option. Writes the pressure of a grid element to a restart file.
- *
- * \param outstream Stream into the restart file.
- * \param element Grid element
- */
+ /*! \name general methods for serialization, output */
+ //@{
+ // serialization methods
+ //! Function needed for restart option.
void serializeEntity(std::ostream &outstream, const Element &element)
{
- int globalIdx = problem_.variables().index(element);
- outstream << pressure_[globalIdx][0];
+ int numFaces = element.template count<1>();
+ for (int i=0; i < numFaces; i++)
+ {
+ int globalIdx = A_.faceMapper().map(element, i, 1);
+ outstream << pressTrace_[globalIdx][0];
+ }
}
- /*! \brief Function for deserialization of the pressure field.
- *
- * Function needed for restart option. Reads the pressure of a grid element from a restart file.
- *
- * \param instream Stream from the restart file.
- * \param element Grid element
- */
void deserializeEntity(std::istream &instream, const Element &element)
{
- int globalIdx = problem_.variables().index(element);
- instream >> pressure_[globalIdx][0];
- }
-
- /*! \brief Adds pressure output to the output file
- *
- * Adds the global pressure to the output.
- *
- * \tparam MultiWriter Class defining the output writer
- * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
- *
- */
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
- {
- ScalarSolution *pressure = writer.allocateManagedBuffer (problem_.gridView().size(0));
-
- *pressure = pressure_;
-
- writer.attachCellData(*pressure, "global pressure");
-
- return;
+ int numFaces = element.template count<1>();
+ for (int i=0; i < numFaces; i++)
+ {
+ int globalIdx = A_.faceMapper().map(element, i, 1);
+ instream >> pressTrace_[globalIdx][0];
+ }
}
+ //@}
//! Constructs a MimeticPressure2P object
/**
- * \param problem A problem class object
+ * \param problem The Dumux problem
*/
MimeticPressure2P(Problem& problem) :
problem_(problem),
- A_(problem.gridView())
+ A_(problem.gridView()), lstiff_(problem_, false, problem_.gridView())
{
- if (pressureType != pGlobal)
+ if (pressureType != pw && pressureType != pn)
{
DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
}
- if (saturationType != sw)
+ if (saturationType != sw && saturationType != sn)
{
DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
}
+ if (GET_PROP_VALUE(TypeTag, EnableCompressibility))
+ {
+ DUNE_THROW(Dune::NotImplemented, "Compressibility not supported!");
+ }
+
+ density_[wPhaseIdx] = 0.0;
+ density_[nPhaseIdx] = 0.0;
+ viscosity_[wPhaseIdx] = 0.0;
+ viscosity_[nPhaseIdx] = 0.0;
+
+ vtkOutputLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, Vtk, OutputLevel);
}
private:
Problem& problem_;
- ScalarSolution pressure_;
- PressureSolution pressTrace_; //!< vector of pressure traces
- NormalVelType normalVelocity_;
- RHSVector f_;
+ TraceType pressTrace_; //!< vector of pressure traces
+ TraceType f_;
OperatorAssembler A_;
+ LocalStiffness lstiff_;
Scalar density_[numPhases];
Scalar viscosity_[numPhases];
- static const int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
- static const int saturationType = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$ 0 = S_w \f$, \f$ 1 = S_n \f$)
+
+ int vtkOutputLevel_;
};
//solves the system of equations to get the spatial distribution of the pressure
@@ -319,74 +424,44 @@ void MimeticPressure2P<TypeTag>::solve()
if (verboseLevelSolver)
std::cout << "MimeticPressure2P: solve for pressure" << std::endl;
+// printmatrix(std::cout, *A_, "global stiffness matrix", "row", 11, 3);
+// printvector(std::cout, f_, "right hand side", "row", 10, 1, 3);
+
Solver solver(problem_);
solver.solve(*A_, pressTrace_, f_);
+
return;
}
-/*! \brief Updates constitutive relations and stores them in the variable class
- *
- * Stores mobility, fractional flow function and capillary pressure for all grid cells.
- */
+//constitutive functions are updated once if new saturations are calculated and stored in the variables object
template<class TypeTag>
void MimeticPressure2P<TypeTag>::updateMaterialLaws()
{
- // iterate through leaf grid an evaluate c0 at cell center
- ElementIterator eItEnd = problem_.gridView().template end<0>();
- for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
- {
- int globalIdx = problem_.variables().index(*eIt);
-
- CellData& cellData = problem_.variables().cellData(globalIdx);
-
- //determine phase saturations from primary saturation variable
+ // iterate through leaf grid an evaluate c0 at cell center
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = problem_.variables().index(*eIt);
- Scalar satW = cellData.saturation(wPhaseIdx);
- Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(*eIt), satW);
+ CellData& cellData = problem_.variables().cellData(globalIdx);
- cellData.setCapillaryPressure(pc);
+ Scalar satW = cellData.saturation(wPhaseIdx);
- // initialize mobilities
- Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx];
- Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx];
+ // initialize mobilities
+ Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW)
+ / viscosity_[wPhaseIdx];
+ Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW)
+ / viscosity_[nPhaseIdx];
- // initialize mobilities
- cellData.setMobility(wPhaseIdx, mobilityW);
- cellData.setMobility(nPhaseIdx, mobilityNW);
+ // initialize mobilities
+ cellData.setMobility(wPhaseIdx, mobilityW);
+ cellData.setMobility(nPhaseIdx, mobilityNW);
- //initialize fractional flow functions
- cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNW));
- cellData.setFracFlowFunc(nPhaseIdx, mobilityNW / (mobilityW + mobilityNW));
- }
- return;
-}
-
-/*! \brief Globally stores the velocity solution at cell-cell interfaces*/
-template<class TypeTag>
-void MimeticPressure2P<TypeTag>::storeVelocity()
-{
- // iterate through leaf grid an evaluate c0 at cell center
- const ElementIterator &eItEnd = this->problem_.gridView().template end<0>();
- for (ElementIterator eIt = this->problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
- {
- int globalIdx = problem_.variables().index(*eIt);
-
- CellData& cellData = problem_.variables().cellData(globalIdx);
-
- IntersectionIterator isIt = problem_.gridView().ibegin(*eIt);
- const IntersectionIterator &isItEnd = problem_.gridView().iend(*eIt);
- for (; isIt != isItEnd; ++isIt)
- {
- int idxInInside = isIt->indexInInside();
- Dune::FieldVector<Scalar, dimWorld> velocity(isIt->centerUnitOuterNormal());
- velocity*= normalVelocity_[globalIdx][idxInInside];
- cellData.fluxData().setVelocity(wPhaseIdx, idxInInside, velocity);
- cellData.fluxData().setUpwindPotential(wPhaseIdx, idxInInside, normalVelocity_[globalIdx][idxInInside]);
- cellData.fluxData().setUpwindPotential(nPhaseIdx, idxInInside, normalVelocity_[globalIdx][idxInInside]);
+ //initialize fractional flow functions
+ cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNW));
+ cellData.setFracFlowFunc(nPhaseIdx, mobilityNW / (mobilityW + mobilityNW));
}
- }
-// printvector(std::cout, problem_.variables().velocity(), "velocity", "row", 4, 1, 3);
- return;
+ return;
}
}
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
new file mode 100644
index 0000000000..0f1f462629
--- /dev/null
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
@@ -0,0 +1,481 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Markus Wolff *
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Model for the pressure equation discretized by mimetic FD.
+ */
+#ifndef DUMUX_MIMETICPRESSURE2PADAPTIVE_HH
+#define DUMUX_MIMETICPRESSURE2PADAPTIVE_HH
+
+// dumux environment
+#include "dumux/decoupled/common/mimetic/mimeticproperties.hh"
+#include "dumux/decoupled/2p/diffusion/mimetic/mimeticoperator2padaptive.hh"
+#include "dumux/decoupled/2p/diffusion/mimetic/mimetic2padaptive.hh"
+
+namespace Dumux
+{
+
+/*! \ingroup Mimetic2p
+ *
+ * \brief mimetic method for the pressure equation
+ *
+ * Provides a mimetic implementation for the evaluation
+ * of equations of the form
+ * \f[\text{div}\, \boldsymbol{v}_{total} = q.\f]
+ * The definition of the total velocity \f$\boldsymbol{v}_total\f$ depends on the kind of pressure chosen. This could be a wetting (w) phase pressure leading to
+ * \f[ - \text{div}\, \left[\lambda \boldsymbol{K} \left(\text{grad}\, p_w + f_n \text{grad}\, p_c + \sum f_\alpha \rho_\alpha g \text{grad}\, z\right)\right] = q, \f]
+ * a non-wetting (n) phase pressure yielding
+ * \f[ - \text{div}\, \left[\lambda \boldsymbol{K} \left(\text{grad}\, p_n - f_w \text{grad}\, p_c + \sum f_\alpha \rho_\alpha g \text{grad}\, z\right)\right] = q, \f]
+ * or a global pressure leading to
+ * \f[ - \text{div}\, \left[\lambda \boldsymbol{K} \left(\text{grad}\, p_{global} + \sum f_\alpha \rho_\alpha g \text{grad}\, z\right)\right] = q.\f]
+ * Here, \f$p\f$ denotes a pressure, \f$\boldsymbol{K}\f$ the absolute permeability, \f$\lambda\f$ the total mobility, possibly depending on the
+ * saturation,\f$f\f$ the fractional flow function of a phase, \f$\rho\f$ a phase density, \f$g\f$ the gravity constant and \f$q\f$ the source term.
+ * For all cases, \f$p = p_D\f$ on \f$\Gamma_{Neumann}\f$, and \f$\boldsymbol{v}_{total} = q_N\f$
+ * on \f$\Gamma_{Dirichlet}\f$.
+ *
+ *\tparam TypeTag The Type Tag
+ */
+template<class TypeTag> class MimeticPressure2PAdaptive
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename SpatialParams::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation), //!< gives kind of pressure used (\f$ 0 = p_w\f$, \f$ 1 = p_n\f$, \f$ 2 = p_{global}\f$)
+ saturationType = GET_PROP_VALUE(TypeTag, SaturationFormulation), //!< gives kind of saturation used (\f$ 0 = S_w\f$, \f$ 1 = S_n\f$)
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, LocalStiffness) LocalStiffness;
+ typedef Dune::BlockVector< Dune::FieldVector<Scalar, 1> > TraceType;
+ typedef MimeticOperatorAssemblerTwoPAdaptive<TypeTag> OperatorAssembler;
+
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) Vector;
+
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElementContainer;
+ typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
+
+ //initializes the matrix to store the system of equations
+ void initializeMatrix();
+
+ //function which assembles the system of equations to be solved
+ void assemble(bool first)
+ {
+ Scalar timeStep = problem_.timeManager().timeStepSize();
+ Scalar maxError = 0.0;
+ int size = problem_.gridView().size(0);
+ for (int i = 0; i < size; i++)
+ {
+ Scalar sat = 0;
+ switch (saturationType)
+ {
+ case Sw:
+ sat = problem_.variables().cellData(i).saturation(wPhaseIdx);
+ break;
+ case Sn:
+ sat = problem_.variables().cellData(i).saturation(nPhaseIdx);
+ break;
+ }
+ if (sat > 1.0)
+ {
+ maxError = std::max(maxError, (sat - 1.0) / timeStep);
+ }
+ if (sat < 0.0)
+ {
+ maxError = std::max(maxError, (-sat) / timeStep);
+ }
+ }
+
+ lstiff_.setErrorInfo(maxError, timeStep);
+ A_.assemble(lstiff_, pressTrace_, f_);
+ return;
+ }
+
+ //solves the system of equations to get the spatial distribution of the pressure
+ void solve();
+
+ void postprocess()
+ {
+ A_.calculatePressure(lstiff_, pressTrace_, problem_);
+
+ return;
+ }
+
+public:
+ //constitutive functions are initialized and stored in the variables object
+ void updateMaterialLaws();
+
+ void initialize(bool solveTwice = true)
+ {
+ ElementIterator element = problem_.gridView().template begin<0> ();
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setPressure(nPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setTemperature(problem_.temperature(*element));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
+ density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
+ updateMaterialLaws();
+ adapt();
+ lstiff_.initialize();
+ lstiff_.reset();
+
+ assemble(true);
+ solve();
+ postprocess();
+
+ return;
+ }
+
+ void adapt()
+ {
+ A_.adapt();
+ pressTrace_.resize(A_.intersectionMapper().size());
+ f_.resize(A_.intersectionMapper().size());
+ pressTrace_ = 0.0;
+ f_ = 0;
+ lstiff_.adapt();
+ }
+
+ void updateVelocity()
+ {
+ updateMaterialLaws();
+ postprocess();
+ }
+
+ void update()
+ {
+ if (problem_.gridAdapt().wasAdapted())
+ {
+ adapt();
+ }
+
+ lstiff_.reset();
+ assemble(false);
+
+ solve();
+
+ postprocess();
+
+ return;
+ }
+
+
+ //! \brief Write data files
+ /* \param name file name */
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ int size = problem_.gridView().size(0);
+ ScalarSolutionType *potential = writer.allocateManagedBuffer(size);
+ ScalarSolutionType *pressure = 0;
+ ScalarSolutionType *pressureSecond = 0;
+ ScalarSolutionType *potentialSecond = 0;
+ Dune::BlockVector < DimVector > *velocityWetting = 0;
+ Dune::BlockVector < DimVector > *velocityNonwetting = 0;
+
+ if (vtkOutputLevel_ > 0)
+ {
+ pressure = writer.allocateManagedBuffer(size);
+ pressureSecond = writer.allocateManagedBuffer(size);
+ potentialSecond = writer.allocateManagedBuffer(size);
+ velocityWetting = writer.template allocateManagedBuffer<Scalar, dim>(size);
+ velocityNonwetting = writer.template allocateManagedBuffer<Scalar, dim>(size);
+ }
+
+
+ ElementIterator eItBegin = problem_.gridView().template begin<0>();
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = eItBegin; eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = problem_.variables().index(*eIt);
+ CellData& cellData = problem_.variables().cellData(globalIdx);
+
+ if (pressureType == pw)
+ {
+ (*potential)[globalIdx] = cellData.potential(wPhaseIdx);
+ }
+
+ if (pressureType == pn)
+ {
+ (*potential)[globalIdx] = cellData.potential(nPhaseIdx);
+ }
+
+ if (vtkOutputLevel_ > 0)
+ {
+
+ if (pressureType == pw)
+ {
+ (*pressure)[globalIdx] = cellData.pressure(wPhaseIdx);
+ (*potentialSecond)[globalIdx] = cellData.potential(nPhaseIdx);
+ (*pressureSecond)[globalIdx] = cellData.pressure(nPhaseIdx);
+ }
+
+ if (pressureType == pn)
+ {
+ (*pressure)[globalIdx] = cellData.pressure(nPhaseIdx);
+ (*potentialSecond)[globalIdx] = cellData.potential(wPhaseIdx);
+ (*pressureSecond)[globalIdx] = cellData.pressure(wPhaseIdx);
+ }
+
+ Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
+ Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
+
+ // run through all intersections with neighbors and boundary
+ IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ int isIndex = isIt->indexInInside();
+
+ fluxW[isIndex] += isIt->geometry().volume()
+ * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
+ fluxNw[isIndex] += isIt->geometry().volume()
+ * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
+ }
+
+ DimVector refVelocity(0);
+ refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
+ refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+ refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
+
+ const DimVector& localPos = ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
+
+ // get the transposed Jacobian of the element mapping
+ const DimMatrix& jacobianT = eIt->geometry().jacobianTransposed(localPos);
+
+ // calculate the element velocity by the Piola transformation
+ DimVector elementVelocity(0);
+ jacobianT.umtv(refVelocity, elementVelocity);
+ elementVelocity /= eIt->geometry().integrationElement(localPos);
+
+ (*velocityWetting)[globalIdx] = elementVelocity;
+
+ refVelocity = 0;
+ refVelocity[0] = 0.5 * (fluxNw[1] - fluxNw[0]);
+ refVelocity[1] = 0.5 * (fluxNw[3] - fluxNw[2]);
+ refVelocity[2] = 0.5 * (fluxNw[5] - fluxNw[4]);
+
+ // calculate the element velocity by the Piola transformation
+ elementVelocity = 0;
+ jacobianT.umtv(refVelocity, elementVelocity);
+ elementVelocity /= eIt->geometry().integrationElement(localPos);
+
+ (*velocityNonwetting)[globalIdx] = elementVelocity;
+ }
+ }
+
+ if (pressureType == pw)
+ {
+ writer.attachCellData(*potential, "wetting potential");
+ }
+
+ if (pressureType == pn)
+ {
+ writer.attachCellData(*potential, "nonwetting potential");
+ }
+
+ if (vtkOutputLevel_ > 0)
+ {
+ if (pressureType == pw)
+ {
+ writer.attachCellData(*pressure, "wetting pressure");
+ writer.attachCellData(*pressureSecond, "nonwetting pressure");
+ writer.attachCellData(*potentialSecond, "nonwetting potential");
+ }
+
+ if (pressureType == pn)
+ {
+ writer.attachCellData(*pressure, "nonwetting pressure");
+ writer.attachCellData(*pressureSecond, "wetting pressure");
+ writer.attachCellData(*potentialSecond, "wetting potential");
+ }
+
+ writer.attachCellData(*velocityWetting, "wetting-velocity", dim);
+ writer.attachCellData(*velocityNonwetting, "non-wetting-velocity", dim);
+ }
+ }
+
+ /*! \name general methods for serialization, output */
+ //@{
+ // serialization methods
+ //! Function needed for restart option.
+ void serializeEntity(std::ostream &outstream, const Element &element)
+ {
+ int numFaces = element.template count<1>();
+ for (int i=0; i < numFaces; i++)
+ {
+ int globalIdx = A_.intersectionMapper().map(element, i);
+ outstream << pressTrace_[globalIdx][0];
+ }
+ }
+
+ void deserializeEntity(std::istream &instream, const Element &element)
+ {
+ int numFaces = element.template count<1>();
+ for (int i=0; i < numFaces; i++)
+ {
+ int globalIdx = A_.intersectionMapper().map(element, i);
+ instream >> pressTrace_[globalIdx][0];
+ }
+ }
+ //@}
+
+ //! Constructs a MimeticPressure2PAdaptive object
+ /**
+ * \param problem The Dumux problem
+ */
+ MimeticPressure2PAdaptive(Problem& problem) :
+ problem_(problem),
+ A_(problem.gridView()), lstiff_(problem_, false, problem_.gridView(), A_.intersectionMapper())
+ {
+ if (pressureType != pw && pressureType != pn)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
+ }
+ if (saturationType != Sw && saturationType != Sn)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
+ }
+ if (GET_PROP_VALUE(TypeTag, EnableCompressibility))
+ {
+ DUNE_THROW(Dune::NotImplemented, "Compressibility not supported!");
+ }
+
+ density_[wPhaseIdx] = 0.0;
+ density_[nPhaseIdx] = 0.0;
+ viscosity_[wPhaseIdx] = 0.0;
+ viscosity_[nPhaseIdx] = 0.0;
+
+ vtkOutputLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, Vtk, OutputLevel);
+ }
+
+private:
+ Problem& problem_;
+ TraceType pressTrace_; //!< vector of pressure traces
+ TraceType f_;
+ OperatorAssembler A_;
+ LocalStiffness lstiff_;
+
+ Scalar density_[numPhases];
+ Scalar viscosity_[numPhases];
+
+ int vtkOutputLevel_;
+};
+
+//solves the system of equations to get the spatial distribution of the pressure
+template<class TypeTag>
+void MimeticPressure2PAdaptive<TypeTag>::solve()
+{
+ typedef typename GET_PROP_TYPE(TypeTag, LinearSolver) Solver;
+
+ int verboseLevelSolver = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
+
+ if (verboseLevelSolver)
+ std::cout << "MimeticPressure2PAdaptive: solve for pressure" << std::endl;
+
+ Solver solver(problem_);
+ solver.solve(*A_, pressTrace_, f_);
+
+// printmatrix(std::cout, *A_, "global stiffness matrix", "row", 11, 3);
+// printvector(std::cout, f_, "right hand side", "row", 200, 1, 3);
+// printvector(std::cout, pressTrace_, "pressure", "row", 200, 1, 3);
+ return;
+}
+
+//constitutive functions are updated once if new saturations are calculated and stored in the variables object
+template<class TypeTag>
+void MimeticPressure2PAdaptive<TypeTag>::updateMaterialLaws()
+{
+ // iterate through leaf grid an evaluate c0 at cell center
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = problem_.variables().index(*eIt);
+
+ CellData& cellData = problem_.variables().cellData(globalIdx);
+
+ Scalar satW = cellData.saturation(wPhaseIdx);
+
+ // initialize mobilities
+ Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW)
+ / viscosity_[wPhaseIdx];
+ Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW)
+ / viscosity_[nPhaseIdx];
+
+ // initialize mobilities
+ cellData.setMobility(wPhaseIdx, mobilityW);
+ cellData.setMobility(nPhaseIdx, mobilityNW);
+
+ //initialize fractional flow functions
+ cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNW));
+ cellData.setFracFlowFunc(nPhaseIdx, mobilityNW / (mobilityW + mobilityNW));
+ }
+ return;
+}
+
+
+
+}
+#endif
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2p.hh
index 0be450181a..39b197ff63 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2p.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2p.hh
@@ -1,7 +1,12 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
- * See the file COPYING for full copying permissions. *
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
@@ -10,20 +15,20 @@
* *
* This program 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 *
+ * 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
- * \ingroup MimeticPressure2p
- * \ingroup IMPETProperties
+ * \ingroup IMPES
+ * \ingroup Properties
*/
/*!
* \file
*
- * \brief Defines the properties required for a two-phase mimetic finite differences model.
+ * \brief Defines the properties required for (immiscible) twophase sequential models.
*/
#ifndef DUMUX_MIMETICPROPERTIES2P_DECOUPLED_HH
@@ -49,7 +54,7 @@ namespace Properties
// Type tags
//////////////////////////////////////////////////////////////////
-//! The type tag for two-phase problems using a mimetic finite differences method.
+//! The type tag for the two-phase problems
NEW_TYPE_TAG(MimeticPressureTwoP, INHERITS_FROM(PressureTwoP, Mimetic))
;
@@ -60,16 +65,14 @@ NEW_TYPE_TAG(MimeticPressureTwoP, INHERITS_FROM(PressureTwoP, Mimetic))
}
#include <dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh>
-#include <dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh>
+#include <dumux/decoupled/2p/diffusion/mimetic/mimetic2p.hh>
namespace Dumux
{
namespace Properties
{
-//! Set mimetic finite differences implementation of the two-phase pressure equation as default pressure model
SET_TYPE_PROP(MimeticPressureTwoP, PressureModel, MimeticPressure2P<TypeTag>);
-//! Set the local stiffness implementation for the two-phase model
-SET_TYPE_PROP(MimeticPressureTwoP, LocalStiffness, MimeticGroundwaterEquationLocalStiffness<TypeTag>);
+SET_TYPE_PROP(MimeticPressureTwoP, LocalStiffness, MimeticTwoPLocalStiffness<TypeTag>);
}
}
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2padaptive.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2padaptive.hh
new file mode 100644
index 0000000000..6962b76c0d
--- /dev/null
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2padaptive.hh
@@ -0,0 +1,79 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup IMPES
+ * \ingroup Properties
+ */
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for (immiscible) twophase sequential models.
+ */
+
+#ifndef DUMUX_MIMETICPROPERTIES2PADAPTIVE_DECOUPLED_HH
+#define DUMUX_MIMETICPROPERTIES2PADAPTIVE_DECOUPLED_HH
+
+//Dumux-includes
+#include <dumux/decoupled/2p/diffusion/diffusionproperties2p.hh>
+#include <dumux/decoupled/common/mimetic/mimeticproperties.hh>
+namespace Dumux
+{
+
+////////////////////////////////
+// forward declarations
+////////////////////////////////
+
+
+////////////////////////////////
+// properties
+////////////////////////////////
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the two-phase problems
+NEW_TYPE_TAG(MimeticPressureTwoPAdaptive, INHERITS_FROM(PressureTwoP, Mimetic))
+;
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+}
+}
+
+#include <dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh>
+#include <dumux/decoupled/2p/diffusion/mimetic/mimetic2padaptive.hh>
+
+namespace Dumux
+{
+namespace Properties
+{
+SET_TYPE_PROP(MimeticPressureTwoPAdaptive, PressureModel, MimeticPressure2PAdaptive<TypeTag>);
+SET_TYPE_PROP(MimeticPressureTwoPAdaptive, LocalStiffness, MimeticTwoPLocalStiffnessAdaptive<TypeTag>);
+}
+}
+
+#endif
--
GitLab