Commit b0db5d0e authored by Andreas Lauser's avatar Andreas Lauser
Browse files

big cleanup

timemanager API, restart, boxmodel are affected

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@3955 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 7d86dff8
......@@ -95,8 +95,7 @@ AC_CONFIG_FILES([dumux.pc
dumux/boxmodels/2p2c/Makefile
dumux/boxmodels/2p2cni/Makefile
dumux/boxmodels/2pni/Makefile
dumux/boxmodels/boxscheme/Makefile
dumux/boxmodels/pdelab/Makefile
dumux/boxmodels/common/Makefile
dumux/boxmodels/richards/Makefile
dumux/common/Makefile
dumux/decoupled/Makefile
......
SUBDIRS = 1p 1p2c 2p 2p2c 2p2cni 2pni boxscheme pdelab richards
SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni 2pni richards
boxmodelsdir = $(includedir)/dumux/boxmodels
......
boxschemedir = $(includedir)/dumux/boxmodels/boxscheme
boxscheme_HEADERS = *.hh
include $(top_srcdir)/am/global-rules
This diff is collapsed.
// $Id: fvelementgeometry.hh 3774 2010-06-24 06:22:44Z bernd $
/*****************************************************************************
* Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
* 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, as long as this copyright notice *
* is included in its original form. *
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
#ifndef DUMUX_FVELEMENTGEOMETRY_HH
#define DUMUX_FVELEMENTGEOMETRY_HH
#if HAVE_DUNE_PDELAB
#include "fvelementgeometry-pdelab.hh"
#else
#include "fvelementgeometry-disc.hh"
#endif
#endif
// $Id$
/*****************************************************************************
* Copyright (C) 2010 by Andreas Lauser *
* 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, as long as this copyright notice *
* is included in its original form. *
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
#ifndef DUMUX_BOX_ELEMENT_BOUNDARY_TYPES_HH
#define DUMUX_BOX_ELEMENT_BOUNDARY_TYPES_HH
#include "boxproperties.hh"
#include <vector>
#include <dumux/common/boundarytypes.hh>
namespace Dumux
{
/*!
* \ingroup BoxModel
*
* \brief This class stores an array of BoundaryTypes objects
*/
template<class TypeTag>
class BoxElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) >
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef std::vector<BoundaryTypes> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
typedef typename RefElemProp::Container ReferenceElements;
typedef typename RefElemProp::ReferenceElement ReferenceElement;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
enum { dim = GridView::dimension };
public:
/*!
* \brief The constructor.
*/
BoxElementBoundaryTypes()
{ }
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvElemGeom)
{
Dune::GeometryType geoType = element.geometry().type();
const ReferenceElement &refElem = ReferenceElements::general(geoType);
int numVerts = element.template count<dim>();
this->resize(numVerts);
for (int i = 0; i < numVerts; ++i)
(*this)[i].reset();
// evaluate boundary conditions
IntersectionIterator isIt = problem.gridView().template ibegin(element);
const IntersectionIterator &endIt = problem.gridView().template iend(element);
for (; isIt != endIt; ++isIt) {
// Ignore non- boundary faces.
if (!isIt->boundary())
continue;
// Set the boundary type for all vertices of the face
int faceIdx = isIt->indexInInside();
int numFaceVerts = refElem.size(faceIdx, 1, dim);
for (int faceVertIdx = 0;
faceVertIdx < numFaceVerts;
faceVertIdx++)
{
int elemVertIdx = refElem.subEntity(faceIdx,
1,
faceVertIdx,
dim);
int boundaryFaceIdx =
fvElemGeom.boundaryFaceIndex(faceIdx,
faceVertIdx);
// set the boundary types
problem.boundaryTypes((*this)[elemVertIdx],
element,
fvElemGeom,
*isIt,
elemVertIdx,
boundaryFaceIdx);
(*this)[elemVertIdx].checkWellPosed();
Valgrind::CheckDefined((*this)[elemVertIdx]);
}
}
};
};
} // namespace Dumux
#endif
// $Id: boundarytypespdelab.hh 3736 2010-06-15 09:52:10Z lauser $
// $Id$
/*****************************************************************************
* Copyright (C) 2009-2010 by Bernd Flemisch, Andreas Lauser *
* Copyright (C) 2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
......@@ -13,37 +13,75 @@
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
#if HAVE_DUNE_PDELAB && !defined DUMUX_BOUNDARYTYPESPDELAB_HH
#define DUMUX_BOUNDARYTYPESPDELAB_HH
#ifndef DUMUX_BOX_ELEMENT_SECONDARY_VARS_HH
#define DUMUX_BOX_ELEMENT_SECONDARY_VARS_HH
#include<dune/common/fvector.hh>
#include<dune/pdelab/common/function.hh>
#include "boxproperties.hh"
template <class TypeTag>
class BoundaryIndexHelperPDELab
#include <vector>
#include <dumux/common/boundarytypes.hh>
namespace Dumux
{
/*!
* \ingroup BoxModel
*
* \brief This class stores an array of SecondaryVars objects
*/
template<class TypeTag>
class BoxElementSecondaryVars : public std::vector<typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) >
{
enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) SecondaryVars;
typedef std::vector<SecondaryVars> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVarVector)) PrimaryVarVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
enum { dim = GridView::dimension };
public:
template <int i>
struct Child {
struct Type {
enum { isLeaf = true };
enum { eqIdx = i };
};
};
/*!
* \brief The constructor.
*/
BoxElementSecondaryVars()
{ }
template <int i>
const typename Child<i>::Type &getChild() const
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvElemGeom,
bool oldSol)
{
static typename Child<i>::Type dummy;
return dummy;
const SolutionVector &globalSol =
oldSol?
problem.model().prevSol():
problem.model().curSol();
const VertexMapper &vertexMapper = problem.vertexMapper();
// we assert that the i-th shape function is
// associated to the i-th vert of the element.
int n = element.template count<dim>();
this->resize(n);
for (int i = 0; i < n; i++) {
const PrimaryVarVector &solI
= globalSol[vertexMapper.map(element, i, dim)];
(*this)[i].update(solI,
problem,
element,
fvElemGeom,
i,
oldSol);
}
};
enum { isLeaf = false };
enum { CHILDREN = numEq };
BoundaryIndexHelperPDELab()
{}
};
} // namespace Dumux
#endif
// $Id: fvelementgeometry-pdelab.hh 3783 2010-06-24 11:33:53Z bernd $
// $Id$
/*****************************************************************************
* Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
* Institute of Hydraulic Engineering *
......@@ -13,13 +13,8 @@
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
#ifndef DUMUX_FVELEMENTGEOMETRY_PDELAB_HH
#define DUMUX_FVELEMENTGEOMETRY_PDELAB_HH
#ifdef DUMUX_FVELEMENTGEOMETRY_DISC_HH
#error "you can not use dune-disc and dune-pdelab at the same time!"
#endif
#ifndef DUMUX_BOX_FV_ELEMENTGEOMETRY_HH
#define DUMUX_BOX_FV_ELEMENTGEOMETRY_HH
#include <dune/grid/common/intersectioniterator.hh>
#include <dune/grid/common/capabilities.hh>
......@@ -37,7 +32,7 @@ NEW_PROP_TAG(Scalar);
/////////////////////
// HACK: Functions to initialize the subcontrol volume data
// structures of FVElementGeomerty.
// structures of BoxFVElementGeomerty.
//
// this is a work around to be able to to use specialization for
// doing this. it is required since it is not possible to
......@@ -46,24 +41,24 @@ NEW_PROP_TAG(Scalar);
/** \todo Please doc me! */
template <typename FVElementGeometry, int dim>
class _FVElemGeomHelper
template <typename BoxFVElementGeometry, int dim>
class _BoxFVElemGeomHelper
{
public:
static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
DUNE_THROW(Dune::NotImplemented, "_FVElemGeomHelper::fillSubContVolData dim = " << dim);
DUNE_THROW(Dune::NotImplemented, "_BoxFVElemGeomHelper::fillSubContVolData dim = " << dim);
};
};
/** \todo Please doc me! */
template <typename FVElementGeometry>
class _FVElemGeomHelper<FVElementGeometry, 1>
template <typename BoxFVElementGeometry>
class _BoxFVElemGeomHelper<BoxFVElementGeometry, 1>
{
public:
enum { dim = 1 };
static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
// 1D
eg.subContVol[0].volume = 0.5*eg.elementVolume;
......@@ -73,13 +68,13 @@ public:
/** \todo Please doc me! */
template <typename FVElementGeometry>
class _FVElemGeomHelper<FVElementGeometry, 2>
template <typename BoxFVElementGeometry>
class _BoxFVElemGeomHelper<BoxFVElementGeometry, 2>
{
public:
enum { dim = 2 };
static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
switch (numVertices) {
case 3: // 2D, triangle
......@@ -88,7 +83,7 @@ public:
eg.subContVol[2].volume = eg.elementVolume/3;
break;
case 4: // 2D, quadrilinear
if (!Dune::Capabilities::IsUnstructured<typename FVElementGeometry::Grid>::v) {
if (!Dune::Capabilities::IsUnstructured<typename BoxFVElementGeometry::Grid>::v) {
eg.subContVol[0].volume = eg.elementVolume/4;
eg.subContVol[1].volume = eg.elementVolume/4;
eg.subContVol[2].volume = eg.elementVolume/4;
......@@ -102,20 +97,20 @@ public:
}
break;
default:
DUNE_THROW(Dune::NotImplemented, "_FVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
DUNE_THROW(Dune::NotImplemented, "_BoxFVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
}
}
};
/** \todo Please doc me! */
template <typename FVElementGeometry>
class _FVElemGeomHelper<FVElementGeometry, 3>
template <typename BoxFVElementGeometry>
class _BoxFVElemGeomHelper<BoxFVElementGeometry, 3>
{
public:
enum { dim = 3 };
static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
switch (numVertices) {
case 4: // 3D, tetrahedron
......@@ -278,7 +273,7 @@ public:
eg.edgeCoord[9]);
break;
default:
DUNE_THROW(Dune::NotImplemented, "_FVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
DUNE_THROW(Dune::NotImplemented, "_BoxFVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
}
}
};
......@@ -288,18 +283,18 @@ public:
/** \todo Please doc me! */
template<class TypeTag>
class FVElementGeometry
class BoxFVElementGeometry
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) LocalFEMSpace;
typedef FVElementGeometry<TypeTag> ThisType;
typedef BoxFVElementGeometry<TypeTag> ThisType;
/** \todo Please doc me! */
friend class _FVElemGeomHelper<ThisType, Grid::dimension>;
friend class _BoxFVElemGeomHelper<ThisType, Grid::dimension>;
typedef _FVElemGeomHelper<ThisType, Grid::dimension> FVElemGeomHelper;
typedef _BoxFVElemGeomHelper<ThisType, Grid::dimension> BoxFVElemGeomHelper;
enum{dim = Grid::dimension};
enum{maxNC = (dim < 3 ? 4 : 8)};
......@@ -443,7 +438,7 @@ class FVElementGeometry
rightFace = edgeToFaceHex[1][k];
break;
default:
DUNE_THROW(Dune::NotImplemented, "FVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
DUNE_THROW(Dune::NotImplemented, "BoxFVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
break;
}
}
......@@ -525,7 +520,7 @@ class FVElementGeometry
rightEdge = faceAndVertexToRightEdgeHex[face][vert];
break;
default:
DUNE_THROW(Dune::NotImplemented, "FVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
DUNE_THROW(Dune::NotImplemented, "BoxFVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
break;
}
}
......@@ -578,13 +573,12 @@ public:
int numFaces; //!< number of faces (0 in < 3D)
const LocalFEMSpace feMap_;
const GridView& gridView_;
FVElementGeometry(const GridView& gridView)
: feMap_(), gridView_(gridView)
BoxFVElementGeometry()
: feMap_()
{}
void update(const Element& e)
void update(const GridView& gridView, const Element& e)
{
const Geometry& geometry = e.geometry();
Dune::GeometryType gt = geometry.type();
......@@ -623,9 +617,9 @@ public:
// fill sub control volume data use specialization for this
// \todo maybe it would be a good idea to migrate everything
// which is dependend of the grid's dimension to
// _FVElemGeomHelper in order to benefit from more aggressive
// _BoxFVElemGeomHelper in order to benefit from more aggressive
// compiler optimizations...
FVElemGeomHelper::fillSubContVolData(*this, numVertices);
BoxFVElemGeomHelper::fillSubContVolData(*this, numVertices);
// fill sub control volume face data:
for (int k = 0; k < numEdges; k++) { // begin loop over edges / sub control volume faces
......@@ -700,8 +694,8 @@ public:
} // end loop over edges / sub control volume faces
// fill boundary face data:
IntersectionIterator endit = gridView_.template iend(e);
for (IntersectionIterator it = gridView_.template ibegin(e); it != endit; ++it)
IntersectionIterator endit = gridView.template iend(e);
for (IntersectionIterator it = gridView.template ibegin(e); it != endit; ++it)
if (it->boundary())
{
int face = it->indexInInside();
......@@ -711,7 +705,7 @@ public:
int vertInElement = referenceElement.subEntity(face, 1, vertInFace, dim);
int bfIdx = boundaryFaceIndex(face, vertInFace);
subContVol[vertInElement].inner = false;
switch (dim) {
switch ((short) dim) {
case 1:
boundaryFace[bfIdx].ipLocal = referenceElement.position(vertInElement, dim);
boundaryFace[bfIdx].area = 1.0;
......@@ -735,7 +729,7 @@ public:
edgeCoord[rightEdge], faceCoord[face], edgeCoord[leftEdge]);
break;
default:
DUNE_THROW(Dune::NotImplemented, "FVElementGeometry for dim = " << dim);
DUNE_THROW(Dune::NotImplemented, "BoxFVElementGeometry for dim = " << dim);
}
boundaryFace[bfIdx].ipGlobal = geometry.global(boundaryFace[bfIdx].ipLocal);
......@@ -830,6 +824,5 @@ public:
}
#endif
// $Id$
/*****************************************************************************
* Copyright (C) 2007 by Peter Bastian *
* Institute of Parallel and Distributed System *
* Department Simulation of Large Systems *
* University of Stuttgart, Germany *
* *
* Copyright (C) 2008-2009 by Andreas Lauser *
* Copyright (C) 2007-2009 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, as long as this copyright notice *
* is included in its original form. *
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
/*!
* \file
* \brief Caculates the jacobian of models based on the box scheme element-wise.
*/
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH
#include <dune/common/exceptions.hh>
#include <dumux/common/valgrind.hh>
#include <dune/grid/common/genericreferenceelements.hh>
#include <boost/format.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
#include "boxelementsecondaryvars.hh"
#include "boxfvelementgeometry.hh"
#include "boxlocalresidual.hh"
#include "boxproperties.hh"
namespace Dumux
{
/*!
* \ingroup BoxModel
* \brief Element-wise caculation of the jacobian matrix for models
* based on the box scheme .
*
* \todo Please doc me more!
*/
template<class TypeTag>
class BoxLocalJacobian
{
private:
typedef BoxLocalJacobian<TypeTag> ThisType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) LocalResidual;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GridView::Grid::ctype CoordScalar;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename Element::EntityPointer ElementPointer;
typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
typedef typename RefElemProp::Container ReferenceElements;
typedef typename RefElemProp::ReferenceElement ReferenceElement;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename Element::Geometry Geometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVarVector)) PrimaryVarVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) SecondaryVars;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSecondaryVars)) ElementSecondaryVars;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
public: