From 8cd8367c7e1eafce3095e0fd56a58394cf3b4d40 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Mon, 3 Dec 2012 15:30:22 +0000
Subject: [PATCH] implicit branch: unify Box- and CC-Model to ImplicitModel.
Delete ccmodel.hh. Deprecate BoxModel in favor of ImplicitModel. Adapt
includes and names.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9755 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/common/boxmodel.hh | 57 +-
dumux/freeflow/stokes/stokeslocalresidual.hh | 2 +-
dumux/freeflow/stokes/stokesmodel.hh | 2 +-
dumux/implicit/1p/1pmodel.hh | 2 +-
dumux/implicit/1p2c/1p2clocalresidual.hh | 2 +-
dumux/implicit/1p2c/1p2cmodel.hh | 2 +-
dumux/implicit/2p/2plocalresidual.hh | 2 +-
dumux/implicit/2p2c/2p2clocalresidual.hh | 2 +-
dumux/implicit/2p2c/2p2cvolumevariables.hh | 2 +-
dumux/implicit/3p3c/3p3clocalresidual.hh | 2 +-
dumux/implicit/3p3c/3p3cvolumevariables.hh | 2 +-
dumux/implicit/box/boxmodel.hh | 926 ------------------
dumux/implicit/box/boxproperties.hh | 2 +
dumux/implicit/box/boxpropertydefaults.hh | 4 -
dumux/implicit/cellcentered/ccmodel.hh | 755 --------------
.../cellcentered/ccpropertydefaults.hh | 4 -
dumux/implicit/common/implicitassembler.hh | 6 +-
dumux/implicit/common/implicitmodel.hh | 127 ++-
dumux/implicit/common/implicitproblem.hh | 1 +
.../common/implicitpropertydefaults.hh | 4 +
dumux/implicit/mpnc/mpnclocalresidual.hh | 2 +-
dumux/implicit/mpnc/mpncmodel.hh | 2 +-
dumux/implicit/richards/richardsmodel.hh | 2 +-
dumux/nonlinear/newtoncontroller.hh | 8 +-
test/implicit/2p/lensproblem.hh | 2 +-
25 files changed, 156 insertions(+), 1766 deletions(-)
delete mode 100644 dumux/implicit/box/boxmodel.hh
delete mode 100644 dumux/implicit/cellcentered/ccmodel.hh
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index c0b766d422..5beaf2c43f 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -1,3 +1,56 @@
-#warning This file is deprecated. Include dumux/implicit/box/boxmodel.hh instead.
+// -*- 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 Base class for models using box discretization
+ */
+#ifndef DUMUX_BOX_MODEL_HH
+#define DUMUX_BOX_MODEL_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/implicit/common/implicitmodel.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup BoxModel
+ *
+ * \brief The base class for the vertex centered finite volume
+ * discretization scheme.
+ */
+template<class TypeTag>
+class BoxModel : public ImplicitModel<TypeTag>
+{
+ typedef ImplicitModel<TypeTag> ParentType;
+
+ // copying a model is not a good idea
+ BoxModel(const BoxModel &);
+
+public:
+ /*!
+ * \brief The constructor.
+ */
+ DUNE_DEPRECATED_MSG("Use ImplicitModel from dumux/implicit/common/implicitmodel.hh instead.")
+ BoxModel() : ParentType()
+ {}
+};
+}
+
+#endif
diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh
index 93ba061821..4f5dd41ac0 100644
--- a/dumux/freeflow/stokes/stokeslocalresidual.hh
+++ b/dumux/freeflow/stokes/stokeslocalresidual.hh
@@ -26,7 +26,7 @@
#ifndef DUMUX_STOKES_LOCAL_RESIDUAL_BASE_HH
#define DUMUX_STOKES_LOCAL_RESIDUAL_BASE_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
//#include <dumux/implicit/box/boxcouplinglocalresidual.hh>
#include "stokesproperties.hh"
diff --git a/dumux/freeflow/stokes/stokesmodel.hh b/dumux/freeflow/stokes/stokesmodel.hh
index 79545f0a87..586fb1486f 100644
--- a/dumux/freeflow/stokes/stokesmodel.hh
+++ b/dumux/freeflow/stokes/stokesmodel.hh
@@ -24,7 +24,7 @@
* \brief Base class for all models which use the Stokes box model.
*/
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include "stokeslocalresidual.hh"
#include "stokesnewtoncontroller.hh"
diff --git a/dumux/implicit/1p/1pmodel.hh b/dumux/implicit/1p/1pmodel.hh
index 14f31be286..40afc03848 100644
--- a/dumux/implicit/1p/1pmodel.hh
+++ b/dumux/implicit/1p/1pmodel.hh
@@ -27,7 +27,7 @@
#ifndef DUMUX_1P_MODEL_HH
#define DUMUX_1P_MODEL_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include "1plocalresidual.hh"
diff --git a/dumux/implicit/1p2c/1p2clocalresidual.hh b/dumux/implicit/1p2c/1p2clocalresidual.hh
index 7412f97448..9dbcde20fe 100644
--- a/dumux/implicit/1p2c/1p2clocalresidual.hh
+++ b/dumux/implicit/1p2c/1p2clocalresidual.hh
@@ -27,7 +27,7 @@
#define DUMUX_ONEP_TWOC_LOCAL_RESIDUAL_HH
#define VELOCITY_OUTPUT 1 //1 turns velocity output on, 0 turns it off
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/implicit/1p2c/1p2cproperties.hh>
#include <dumux/implicit/1p2c/1p2cvolumevariables.hh>
diff --git a/dumux/implicit/1p2c/1p2cmodel.hh b/dumux/implicit/1p2c/1p2cmodel.hh
index 31f295e3f2..3b70302f85 100644
--- a/dumux/implicit/1p2c/1p2cmodel.hh
+++ b/dumux/implicit/1p2c/1p2cmodel.hh
@@ -30,7 +30,7 @@
#include "1p2cproperties.hh"
#include "1p2clocalresidual.hh"
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
namespace Dumux
{
diff --git a/dumux/implicit/2p/2plocalresidual.hh b/dumux/implicit/2p/2plocalresidual.hh
index 83b1ceac8b..c72e5046c8 100644
--- a/dumux/implicit/2p/2plocalresidual.hh
+++ b/dumux/implicit/2p/2plocalresidual.hh
@@ -24,7 +24,7 @@
#ifndef DUMUX_TWOP_LOCAL_RESIDUAL_BASE_HH
#define DUMUX_TWOP_LOCAL_RESIDUAL_BASE_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include "2pproperties.hh"
diff --git a/dumux/implicit/2p2c/2p2clocalresidual.hh b/dumux/implicit/2p2c/2p2clocalresidual.hh
index 9f45f06a3c..d36494883f 100644
--- a/dumux/implicit/2p2c/2p2clocalresidual.hh
+++ b/dumux/implicit/2p2c/2p2clocalresidual.hh
@@ -26,7 +26,7 @@
#ifndef DUMUX_2P2C_LOCAL_RESIDUAL_BASE_HH
#define DUMUX_2P2C_LOCAL_RESIDUAL_BASE_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>
#include "2p2cproperties.hh"
diff --git a/dumux/implicit/2p2c/2p2cvolumevariables.hh b/dumux/implicit/2p2c/2p2cvolumevariables.hh
index 0546e6258f..ac9719c58b 100644
--- a/dumux/implicit/2p2c/2p2cvolumevariables.hh
+++ b/dumux/implicit/2p2c/2p2cvolumevariables.hh
@@ -25,7 +25,7 @@
#ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH
#define DUMUX_2P2C_VOLUME_VARIABLES_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>
#include <dune/common/collectivecommunication.hh>
diff --git a/dumux/implicit/3p3c/3p3clocalresidual.hh b/dumux/implicit/3p3c/3p3clocalresidual.hh
index aa39b2dacb..03b0e2474f 100644
--- a/dumux/implicit/3p3c/3p3clocalresidual.hh
+++ b/dumux/implicit/3p3c/3p3clocalresidual.hh
@@ -25,7 +25,7 @@
#ifndef DUMUX_3P3C_LOCAL_RESIDUAL_HH
#define DUMUX_3P3C_LOCAL_RESIDUAL_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>
#include "3p3cproperties.hh"
diff --git a/dumux/implicit/3p3c/3p3cvolumevariables.hh b/dumux/implicit/3p3c/3p3cvolumevariables.hh
index f4f4ef3570..0049fe5290 100644
--- a/dumux/implicit/3p3c/3p3cvolumevariables.hh
+++ b/dumux/implicit/3p3c/3p3cvolumevariables.hh
@@ -25,7 +25,7 @@
#ifndef DUMUX_3P3C_VOLUME_VARIABLES_HH
#define DUMUX_3P3C_VOLUME_VARIABLES_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>
#include <dune/common/collectivecommunication.hh>
diff --git a/dumux/implicit/box/boxmodel.hh b/dumux/implicit/box/boxmodel.hh
deleted file mode 100644
index a73f4741e5..0000000000
--- a/dumux/implicit/box/boxmodel.hh
+++ /dev/null
@@ -1,926 +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 Base class for models using box discretization
- */
-#ifndef DUMUX_BOX_MODEL_HH
-#define DUMUX_BOX_MODEL_HH
-
-#include "boxproperties.hh"
-#include "boxpropertydefaults.hh"
-
-#include "boxelementvolumevariables.hh"
-#include "boxlocalresidual.hh"
-
-#include <dumux/parallel/vertexhandles.hh>
-
-#include <dune/grid/common/geometry.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup BoxModel
- *
- * \brief The base class for the vertex centered finite volume
- * discretization scheme.
- */
-template<class TypeTag>
-class BoxModel
-{
- typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, NumEq),
- dim = GridView::dimension
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) LocalJacobian;
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
-
- typedef typename GridView::ctype CoordScalar;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef typename GridView::template Codim<dim>::Entity Vertex;
- typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
- typedef typename GridView::IntersectionIterator IntersectionIterator;
-
- typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
-
- // copying a model is not a good idea
- BoxModel(const BoxModel &);
-
-public:
- /*!
- * \brief The constructor.
- */
- BoxModel()
- {
- enableHints_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableHints);
- }
-
- ~BoxModel()
- { delete jacAsm_; }
-
- /*!
- * \brief Apply the initial conditions to the model.
- *
- * \param problem The object representing the problem which needs to
- * be simulated.
- */
- void init(Problem &problem)
- {
- problemPtr_ = &problem;
-
- updateBoundaryIndices_();
-
- int nDofs = asImp_().numDofs();
- uCur_.resize(nDofs);
- uPrev_.resize(nDofs);
- boxVolume_.resize(nDofs);
-
- localJacobian_.init(problem_());
- jacAsm_ = new JacobianAssembler();
- jacAsm_->init(problem_());
-
- asImp_().applyInitialSolution_();
-
- // resize the hint vectors
- if (enableHints_) {
- int nVerts = gridView_().size(dim);
- curHints_.resize(nVerts);
- prevHints_.resize(nVerts);
- hintsUsable_.resize(nVerts);
- std::fill(hintsUsable_.begin(),
- hintsUsable_.end(),
- false);
- }
-
- // also set the solution of the "previous" time step to the
- // initial solution.
- uPrev_ = uCur_;
- }
-
- void setHints(const Element &element,
- ElementVolumeVariables &prevVolVars,
- ElementVolumeVariables &curVolVars) const
- {
- if (!enableHints_)
- return;
-
- int n = element.template count<dim>();
- prevVolVars.resize(n);
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i) {
- int globalIdx = vertexMapper().map(element, i, dim);
-
- if (!hintsUsable_[globalIdx]) {
- curVolVars[i].setHint(NULL);
- prevVolVars[i].setHint(NULL);
- }
- else {
- curVolVars[i].setHint(&curHints_[globalIdx]);
- prevVolVars[i].setHint(&prevHints_[globalIdx]);
- }
- }
- }
-
- void setHints(const Element &element,
- ElementVolumeVariables &curVolVars) const
- {
- if (!enableHints_)
- return;
-
- int n = element.template count<dim>();
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i) {
- int globalIdx = vertexMapper().map(element, i, dim);
-
- if (!hintsUsable_[globalIdx])
- curVolVars[i].setHint(NULL);
- else
- curVolVars[i].setHint(&curHints_[globalIdx]);
- }
- }
-
- void updatePrevHints()
- {
- if (!enableHints_)
- return;
-
- prevHints_ = curHints_;
- }
-
- void updateCurHints(const Element &element,
- const ElementVolumeVariables &elemVolVars) const
- {
- if (!enableHints_)
- return;
-
- for (unsigned int i = 0; i < elemVolVars.size(); ++i) {
- int globalIdx = vertexMapper().map(element, i, dim);
- curHints_[globalIdx] = elemVolVars[i];
- if (!hintsUsable_[globalIdx])
- prevHints_[globalIdx] = elemVolVars[i];
- hintsUsable_[globalIdx] = true;
- }
- }
-
-
- /*!
- * \brief Compute the global residual for an arbitrary solution
- * vector.
- *
- * \param residual Stores the result
- * \param u The solution for which the residual ought to be calculated
- */
- Scalar globalResidual(SolutionVector &residual,
- const SolutionVector &u)
- {
- SolutionVector tmp(curSol());
- curSol() = u;
- Scalar res = globalResidual(residual);
- curSol() = tmp;
- return res;
- }
-
- /*!
- * \brief Compute the global residual for the current solution
- * vector.
- *
- * \param residual Stores the result
- */
- Scalar globalResidual(SolutionVector &residual)
- {
- residual = 0;
-
- ElementIterator elemIt = gridView_().template begin<0>();
- const ElementIterator elemEndIt = gridView_().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- localResidual().eval(*elemIt);
-
- for (int i = 0; i < elemIt->template count<dim>(); ++i) {
- int globalI = vertexMapper().map(*elemIt, i, dim);
- residual[globalI] += localResidual().residual(i);
- }
- }
-
- // calculate the square norm of the residual
- Scalar result2 = residual.two_norm2();
- if (gridView_().comm().size() > 1)
- result2 = gridView_().comm().sum(result2);
-
- // add up the residuals on the process borders
- if (gridView_().comm().size() > 1) {
- VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
- sumHandle(residual, vertexMapper());
- gridView_().communicate(sumHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
- }
-
- return std::sqrt(result2);
- }
-
- /*!
- * \brief Compute the integral over the domain of the storage
- * terms of all conservation quantities.
- *
- * \param storage Stores the result
- */
- void globalStorage(PrimaryVariables &storage)
- {
- storage = 0;
-
- ElementIterator elemIt = gridView_().template begin<0>();
- const ElementIterator elemEndIt = gridView_().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- localResidual().evalStorage(*elemIt);
-
- for (int i = 0; i < elemIt->template count<dim>(); ++i)
- storage += localResidual().storageTerm()[i];
- }
-
- if (gridView_().comm().size() > 1)
- storage = gridView_().comm().sum(storage);
- }
-
- /*!
- * \brief Returns the volume \f$\mathrm{[m^3]}\f$ of a given control volume.
- *
- * \param globalIdx The global index of the control volume's
- * associated vertex
- */
- Scalar boxVolume(const int globalIdx) const
- { return boxVolume_[globalIdx][0]; }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- const SolutionVector &curSol() const
- { return uCur_; }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- SolutionVector &curSol()
- { return uCur_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- const SolutionVector &prevSol() const
- { return uPrev_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- SolutionVector &prevSol()
- { return uPrev_; }
-
- /*!
- * \brief Returns the operator assembler for the global jacobian of
- * the problem.
- */
- JacobianAssembler &jacobianAssembler()
- { return *jacAsm_; }
-
- /*!
- * \copydoc jacobianAssembler()
- */
- const JacobianAssembler &jacobianAssembler() const
- { return *jacAsm_; }
-
- /*!
- * \brief Returns the local jacobian which calculates the local
- * stiffness matrix for an arbitrary element.
- *
- * The local stiffness matrices of the element are used by
- * the jacobian assembler to produce a global linerization of the
- * problem.
- */
- LocalJacobian &localJacobian()
- { return localJacobian_; }
- /*!
- * \copydoc localJacobian()
- */
- const LocalJacobian &localJacobian() const
- { return localJacobian_; }
-
- /*!
- * \brief Returns the local residual function.
- */
- LocalResidual &localResidual()
- { return localJacobian().localResidual(); }
- /*!
- * \copydoc localResidual()
- */
- const LocalResidual &localResidual() const
- { return localJacobian().localResidual(); }
-
- /*!
- * \brief Returns the relative error between two vectors of
- * primary variables.
- *
- * \param vertexIdx The global index of the control volume's
- * associated vertex
- * \param priVars1 The first vector of primary variables
- * \param priVars2 The second vector of primary variables
- *
- * \todo The vertexIdx argument is pretty hacky. it is required by
- * models with pseudo primary variables (i.e. the primary
- * variable switching models). the clean solution would be
- * to access the pseudo primary variables from the primary
- * variables.
- */
- Scalar relativeErrorVertex(const int vertexIdx,
- const PrimaryVariables &priVars1,
- const PrimaryVariables &priVars2)
- {
- Scalar result = 0.0;
- for (int j = 0; j < numEq; ++j) {
- Scalar eqErr = std::abs(priVars1[j] - priVars2[j]);
- eqErr /= std::max<Scalar>(1.0, std::abs(priVars1[j] + priVars2[j])/2);
-
- result = std::max(result, eqErr);
- }
- return result;
- }
-
- /*!
- * \brief Try to progress the model to the next timestep.
- *
- * \param solver The non-linear solver
- * \param controller The controller which specifies the behaviour
- * of the non-linear solver
- */
- bool update(NewtonMethod &solver,
- NewtonController &controller)
- {
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i)
- Valgrind::CheckDefined(curSol()[i]);
-#endif // HAVE_VALGRIND
-
- asImp_().updateBegin();
-
- bool converged = solver.execute(controller);
- if (converged) {
- asImp_().updateSuccessful();
- }
- else
- asImp_().updateFailed();
-
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i) {
- Valgrind::CheckDefined(curSol()[i]);
- }
-#endif // HAVE_VALGRIND
-
- return converged;
- }
-
-
- /*!
- * \brief Called by the update() method before it tries to
- * apply the newton method. This is primary a hook
- * which the actual model can overload.
- */
- void updateBegin()
- { }
-
-
- /*!
- * \brief Called by the update() method if it was
- * successful. This is primary a hook which the actual
- * model can overload.
- */
- void updateSuccessful()
- { }
-
- /*!
- * \brief Called by the update() method if it was
- * unsuccessful. This is primary a hook which the actual
- * model can overload.
- */
- void updateFailed()
- {
- // Reset the current solution to the one of the
- // previous time step so that we can start the next
- // update at a physically meaningful solution.
- uCur_ = uPrev_;
- curHints_ = prevHints_;
-
- jacAsm_->reassembleAll();
- }
-
- /*!
- * \brief Called by the problem if a time integration was
- * successful, post processing of the solution is done and
- * the result has been written to disk.
- *
- * This should prepare the model for the next time integration.
- */
- void advanceTimeLevel()
- {
- // make the current solution the previous one.
- uPrev_ = uCur_;
- prevHints_ = curHints_;
-
- updatePrevHints();
- }
-
- /*!
- * \brief Serializes the current state of the model.
- *
- * \tparam Restarter The type of the serializer class
- *
- * \param res The serializer object
- */
- template <class Restarter>
- void serialize(Restarter &res)
- { res.template serializeEntities<dim>(asImp_(), this->gridView_()); }
-
- /*!
- * \brief Deserializes the state of the model.
- *
- * \tparam Restarter The type of the serializer class
- *
- * \param res The serializer object
- */
- template <class Restarter>
- void deserialize(Restarter &res)
- {
- res.template deserializeEntities<dim>(asImp_(), this->gridView_());
- prevSol() = curSol();
- }
-
- /*!
- * \brief Write the current solution for a vertex to a restart
- * file.
- *
- * \param outstream The stream into which the vertex data should
- * be serialized to
- * \param vertex The DUNE Codim<dim> entity which's data should be
- * serialized
- */
- void serializeEntity(std::ostream &outstream,
- const Vertex &vertex)
- {
- int vertIdx = dofMapper().map(vertex);
-
- // write phase state
- if (!outstream.good()) {
- DUNE_THROW(Dune::IOError,
- "Could not serialize vertex "
- << vertIdx);
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- outstream << curSol()[vertIdx][eqIdx] << " ";
- }
- }
-
- /*!
- * \brief Reads the current solution variables for a vertex from a
- * restart file.
- *
- * \param instream The stream from which the vertex data should
- * be deserialized from
- * \param vertex The DUNE Codim<dim> entity which's data should be
- * deserialized
- */
- void deserializeEntity(std::istream &instream,
- const Vertex &vertex)
- {
- int vertIdx = dofMapper().map(vertex);
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (!instream.good())
- DUNE_THROW(Dune::IOError,
- "Could not deserialize vertex "
- << vertIdx);
- instream >> curSol()[vertIdx][eqIdx];
- }
- }
-
- /*!
- * \brief Returns the number of global degrees of freedoms (DOFs)
- */
- size_t numDofs() const
- { return gridView_().size(dim); }
-
- /*!
- * \brief Mapper for the entities where degrees of freedoms are
- * defined to indices.
- *
- * This usually means a mapper for vertices.
- */
- const DofMapper &dofMapper() const
- { return problem_().vertexMapper(); }
-
- /*!
- * \brief Mapper for vertices to indices.
- */
- const VertexMapper &vertexMapper() const
- { return problem_().vertexMapper(); }
-
- /*!
- * \brief Mapper for elements to indices.
- */
- const ElementMapper &elementMapper() const
- { return problem_().elementMapper(); }
-
- /*!
- * \brief Resets the Jacobian matrix assembler, so that the
- * boundary types can be altered.
- */
- void resetJacobianAssembler ()
- {
- delete jacAsm_;
- jacAsm_ = new JacobianAssembler;
- jacAsm_->init(problem_());
- }
-
- /*!
- * \brief Update the weights of all primary variables within an
- * element given the complete set of volume variables
- *
- * \param element The DUNE codim 0 entity
- * \param volVars All volume variables for the element
- */
- void updatePVWeights(const Element &element,
- const ElementVolumeVariables &volVars) const
- { }
-
- /*!
- * \brief Add the vector fields for analysing the convergence of
- * the newton method to the a VTK multi writer.
- *
- * \tparam MultiWriter The type of the VTK multi writer
- *
- * \param writer The VTK multi writer object on which the fields should be added.
- * \param u The solution function
- * \param deltaU The delta of the solution function before and after the Newton update
- */
- template <class MultiWriter>
- void addConvergenceVtkFields(MultiWriter &writer,
- const SolutionVector &u,
- const SolutionVector &deltaU)
- {
- typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
-
- SolutionVector residual(u);
- asImp_().globalResidual(residual, u);
-
- // create the required scalar fields
- unsigned numVertices = this->gridView_().size(dim);
-
- // global defect of the two auxiliary equations
- ScalarField* def[numEq];
- ScalarField* delta[numEq];
- ScalarField* x[numEq];
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- x[eqIdx] = writer.allocateManagedBuffer(numVertices);
- delta[eqIdx] = writer.allocateManagedBuffer(numVertices);
- def[eqIdx] = writer.allocateManagedBuffer(numVertices);
- }
-
- VertexIterator vIt = this->gridView_().template begin<dim>();
- VertexIterator vEndIt = this->gridView_().template end<dim>();
- for (; vIt != vEndIt; ++ vIt)
- {
- int globalIdx = vertexMapper().map(*vIt);
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- (*x[eqIdx])[globalIdx] = u[globalIdx][eqIdx];
- (*delta[eqIdx])[globalIdx] = - deltaU[globalIdx][eqIdx];
- (*def[eqIdx])[globalIdx] = residual[globalIdx][eqIdx];
- }
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- std::ostringstream oss;
- oss.str(""); oss << "x_" << eqIdx;
- writer.attachVertexData(*x[eqIdx], oss.str());
- oss.str(""); oss << "delta_" << eqIdx;
- writer.attachVertexData(*delta[eqIdx], oss.str());
- oss.str(""); oss << "defect_" << eqIdx;
- writer.attachVertexData(*def[eqIdx], oss.str());
- }
-
- asImp_().addOutputVtkFields(u, writer);
- }
-
- /*!
- * \brief Add the quantities of a time step which ought to be written to disk.
- *
- * This should be overwritten by the actual model if any secondary
- * variables should be written out. Read: This should _always_ be
- * overwritten by well behaved models!
- *
- * \tparam MultiWriter The type of the VTK multi writer
- *
- * \param sol The global vector of primary variable values.
- * \param writer The VTK multi writer where the fields should be added.
- */
- template <class MultiWriter>
- void addOutputVtkFields(const SolutionVector &sol,
- MultiWriter &writer)
- {
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
-
- // create the required scalar fields
- unsigned numVertices = this->gridView_().size(dim);
-
- // global defect of the two auxiliary equations
- ScalarField* x[numEq];
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- x[eqIdx] = writer.allocateManagedBuffer(numVertices);
- }
-
- VertexIterator vIt = this->gridView_().template begin<dim>();
- VertexIterator vEndIt = this->gridView_().template end<dim>();
- for (; vIt != vEndIt; ++ vIt)
- {
- int globalIdx = vertexMapper().map(*vIt);
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- (*x[eqIdx])[globalIdx] = sol[globalIdx][eqIdx];
- }
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- std::ostringstream oss;
- oss << "primaryVar_" << eqIdx;
- writer.attachVertexData(*x[eqIdx], oss.str());
- }
- }
-
- /*!
- * \brief Reference to the grid view of the spatial domain.
- */
- const GridView &gridView() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns true if the vertex with 'globalVertIdx' is
- * located on the grid's boundary.
- *
- * \param globalVertIdx The global index of the control volume's
- * associated vertex
- */
- bool onBoundary(const int globalVertIdx) const
- { return boundaryIndices_[globalVertIdx]; }
-
- /*!
- * \brief Returns true if a vertex is located on the grid's
- * boundary.
- *
- * \param element A DUNE Codim<0> entity which contains the control
- * volume's associated vertex.
- * \param vIdx The local vertex index inside element
- */
- bool onBoundary(const Element &element, const int vIdx) const
- { return onBoundary(vertexMapper().map(element, vIdx, dim)); }
-
- /*!
- * \brief Fill the fluid state according to the primary variables.
- *
- * Taking the information from the primary variables,
- * the fluid state is filled with every information that is
- * necessary to evaluate the model's local residual.
- *
- * \param priVars The primary variables of the model.
- * \param problem The problem at hand.
- * \param element The current element.
- * \param fvGeometry The finite volume element geometry.
- * \param scvIdx The index of the subcontrol volume.
- * \param fluidState The fluid state to fill.
- */
- template <class FluidState>
- static void completeFluidState(const PrimaryVariables& priVars,
- const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const int scvIdx,
- FluidState& fluidState)
- {
- VolumeVariables::completeFluidState(priVars, problem, element,
- fvGeometry, scvIdx, fluidState);
- }
-protected:
- /*!
- * \brief A reference to the problem on which the model is applied.
- */
- Problem &problem_()
- { return *problemPtr_; }
- /*!
- * \copydoc problem_()
- */
- const Problem &problem_() const
- { return *problemPtr_; }
-
- /*!
- * \brief Reference to the grid view of the spatial domain.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Reference to the local residal object
- */
- LocalResidual &localResidual_()
- { return localJacobian_.localResidual(); }
-
- /*!
- * \brief Applies the initial solution for all vertices of the grid.
- */
- void applyInitialSolution_()
- {
- // first set the whole domain to zero
- uCur_ = Scalar(0.0);
- boxVolume_ = Scalar(0.0);
-
- FVElementGeometry fvGeometry;
-
- // iterate through leaf grid and evaluate initial
- // condition at the center of each sub control volume
- //
- // TODO: the initial condition needs to be unique for
- // each vertex. we should think about the API...
- ElementIterator eIt = gridView_().template begin<0>();
- const ElementIterator &eEndIt = gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt) {
- // deal with the current element
- fvGeometry.update(gridView_(), *eIt);
-
- // loop over all element vertices, i.e. sub control volumes
- for (int scvIdx = 0; scvIdx < fvGeometry.numVertices; scvIdx++)
- {
- // map the local vertex index to the global one
- int globalIdx = vertexMapper().map(*eIt,
- scvIdx,
- dim);
-
- // let the problem do the dirty work of nailing down
- // the initial solution.
- PrimaryVariables initPriVars;
- Valgrind::SetUndefined(initPriVars);
- problem_().initial(initPriVars,
- *eIt,
- fvGeometry,
- scvIdx);
- Valgrind::CheckDefined(initPriVars);
-
- // add up the initial values of all sub-control
- // volumes. If the initial values disagree for
- // different sub control volumes, the initial value
- // will be the arithmetic mean.
- initPriVars *= fvGeometry.subContVol[scvIdx].volume;
- boxVolume_[globalIdx] += fvGeometry.subContVol[scvIdx].volume;
- uCur_[globalIdx] += initPriVars;
- Valgrind::CheckDefined(uCur_[globalIdx]);
- }
- }
-
- // add up the primary variables and the volumes of the boxes
- // which cross process borders
- if (gridView_().comm().size() > 1) {
- VertexHandleSum<Dune::FieldVector<Scalar, 1>,
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> >,
- VertexMapper> sumVolumeHandle(boxVolume_, vertexMapper());
- gridView_().communicate(sumVolumeHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
-
- VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
- sumPVHandle(uCur_, vertexMapper());
- gridView_().communicate(sumPVHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
- }
-
- // divide all primary variables by the volume of their boxes
- int n = gridView_().size(dim);
- for (int i = 0; i < n; ++i) {
- uCur_[i] /= boxVolume(i);
- }
- }
-
- /*!
- * \brief Find all indices of boundary vertices.
- *
- * For this we need to loop over all intersections (which is slow
- * in general). If the DUNE grid interface would provide a
- * onBoundary() method for entities this could be done in a much
- * nicer way (actually this would not be necessary)
- */
- void updateBoundaryIndices_()
- {
- boundaryIndices_.resize(numDofs());
- std::fill(boundaryIndices_.begin(), boundaryIndices_.end(), false);
-
- ElementIterator eIt = gridView_().template begin<0>();
- ElementIterator eEndIt = gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt) {
- Dune::GeometryType geoType = eIt->geometry().type();
- const ReferenceElement &refElement = ReferenceElements::general(geoType);
-
- IntersectionIterator isIt = gridView_().ibegin(*eIt);
- IntersectionIterator isEndIt = gridView_().iend(*eIt);
- for (; isIt != isEndIt; ++isIt) {
- if (isIt->boundary()) {
- // add all vertices on the intersection to the set of
- // boundary vertices
- int faceIdx = isIt->indexInInside();
- int numFaceVerts = refElement.size(faceIdx, 1, dim);
- for (int faceVertIdx = 0;
- faceVertIdx < numFaceVerts;
- ++faceVertIdx)
- {
- int elemVertIdx = refElement.subEntity(faceIdx,
- 1,
- faceVertIdx,
- dim);
- int globalVertIdx = vertexMapper().map(*eIt, elemVertIdx, dim);
- boundaryIndices_[globalVertIdx] = true;
- }
- }
- }
- }
- }
-
- // the hint cache for the previous and the current volume
- // variables
- mutable std::vector<bool> hintsUsable_;
- mutable std::vector<VolumeVariables> curHints_;
- mutable std::vector<VolumeVariables> prevHints_;
-
- // the problem we want to solve. defines the constitutive
- // relations, matxerial laws, etc.
- Problem *problemPtr_;
-
- // calculates the local jacobian matrix for a given element
- LocalJacobian localJacobian_;
- // Linearizes the problem at the current time step using the
- // local jacobian
- JacobianAssembler *jacAsm_;
-
- // the set of all indices of vertices on the boundary
- std::vector<bool> boundaryIndices_;
-
- // cur is the current iterative solution, prev the converged
- // solution of the previous time step
- SolutionVector uCur_;
- SolutionVector uPrev_;
-
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
-
-private:
- /*!
- * \brief Returns whether messages should be printed
- */
- bool verbose_() const
- { return gridView_().comm().rank() == 0; }
-
- Implementation &asImp_()
- { return *static_cast<Implementation*>(this); }
- const Implementation &asImp_() const
- { return *static_cast<const Implementation*>(this); }
-
- bool enableHints_;
-};
-}
-
-#endif
diff --git a/dumux/implicit/box/boxproperties.hh b/dumux/implicit/box/boxproperties.hh
index 3bd1812dd8..75b1fc3733 100644
--- a/dumux/implicit/box/boxproperties.hh
+++ b/dumux/implicit/box/boxproperties.hh
@@ -59,4 +59,6 @@ NEW_PROP_TAG(ImplicitUseTwoPointFlux); //! indicates whether two-point flux shou
// \}
+#include "boxpropertydefaults.hh"
+
#endif
diff --git a/dumux/implicit/box/boxpropertydefaults.hh b/dumux/implicit/box/boxpropertydefaults.hh
index 6e592e0b45..1487bce6d4 100644
--- a/dumux/implicit/box/boxpropertydefaults.hh
+++ b/dumux/implicit/box/boxpropertydefaults.hh
@@ -29,7 +29,6 @@
#include <dumux/implicit/common/implicitpropertydefaults.hh>
#include "boxassembler.hh"
-#include "boxmodel.hh"
#include "boxfvelementgeometry.hh"
#include "boxelementboundarytypes.hh"
#include "boxlocalresidual.hh"
@@ -60,9 +59,6 @@ SET_TYPE_PROP(BoxModel, DofMapper, typename GET_PROP_TYPE(TypeTag, VertexMapper)
//! Set the BaseLocalResidual to BoxLocalResidual
SET_TYPE_PROP(BoxModel, BaseLocalResidual, Dumux::BoxLocalResidual<TypeTag>);
-//! Set the BaseModel to BoxModel
-SET_TYPE_PROP(BoxModel, BaseModel, Dumux::BoxModel<TypeTag>);
-
//! An array of secondary variable containers
SET_TYPE_PROP(BoxModel, ElementVolumeVariables, Dumux::BoxElementVolumeVariables<TypeTag>);
diff --git a/dumux/implicit/cellcentered/ccmodel.hh b/dumux/implicit/cellcentered/ccmodel.hh
deleted file mode 100644
index 00e30d5407..0000000000
--- a/dumux/implicit/cellcentered/ccmodel.hh
+++ /dev/null
@@ -1,755 +0,0 @@
-// -*- 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-2010 by Andreas Lauser *
- * Copyright (C) 2008-2010 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 Base class for models using box discretization
- */
-#ifndef DUMUX_CC_MODEL_HH
-#define DUMUX_CC_MODEL_HH
-
-#include "ccproperties.hh"
-#include "ccpropertydefaults.hh"
-
-#include "ccelementvolumevariables.hh"
-#include "cclocalresidual.hh"
-
-#include <dumux/implicit/common/implicitlocaljacobian.hh>
-#include <dumux/parallel/vertexhandles.hh>
-
-#include <dune/grid/common/geometry.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup CCModel
- *
- * \brief The base class for the vertex centered finite volume
- * discretization scheme.
- */
-template<class TypeTag>
-class CCModel
-{
- typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, NumEq),
- dim = GridView::dimension
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) LocalJacobian;
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
-
- typedef typename GridView::ctype CoordScalar;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef typename GridView::template Codim<dim>::Entity Vertex;
- typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
- typedef typename GridView::IntersectionIterator IntersectionIterator;
-
- typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
-
- // copying a model is not a good idea
- CCModel(const CCModel &);
-
-public:
- /*!
- * \brief The constructor.
- */
- CCModel()
- {}
-
- ~CCModel()
- { delete jacAsm_; }
-
- /*!
- * \brief Apply the initial conditions to the model.
- *
- * \param prob The object representing the problem which needs to
- * be simulated.
- */
- void init(Problem &prob)
- {
- problemPtr_ = &prob;
-
- updateBoundaryIndices_();
-
- int nDofs = asImp_().numDofs();
- uCur_.resize(nDofs);
- uPrev_.resize(nDofs);
-
- localJacobian_.init(problem_());
- jacAsm_ = new JacobianAssembler();
- jacAsm_->init(problem_());
-
- asImp_().applyInitialSolution_();
-
- // also set the solution of the "previous" time step to the
- // initial solution.
- uPrev_ = uCur_;
- }
-
- // functions for interface consistence
- void setHints(const Element &elem,
- ElementVolumeVariables &prevVolVars,
- ElementVolumeVariables &curVolVars) const
- {}
-
- void setHints(const Element &elem,
- ElementVolumeVariables &curVolVars) const
- {}
-
- void updatePrevHints()
- {}
-
- void updateCurHints(const Element &elem,
- const ElementVolumeVariables &ev) const
- {}
-
- /*!
- * \brief Compute the global residual for an arbitrary solution
- * vector.
- *
- * \param dest Stores the result
- * \param u The solution for which the residual ought to be calculated
- */
- Scalar globalResidual(SolutionVector &dest,
- const SolutionVector &u)
- {
- SolutionVector tmp(curSol());
- curSol() = u;
- Scalar res = globalResidual(dest);
- curSol() = tmp;
- return res;
- }
-
- /*!
- * \brief Compute the global residual for the current solution
- * vector.
- *
- * \param dest Stores the result
- */
- Scalar globalResidual(SolutionVector &dest)
- {
- dest = 0;
-
- ElementIterator elemIt = gridView_().template begin<0>();
- const ElementIterator elemEndIt = gridView_().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- localResidual().eval(*elemIt);
-
- int globalI = elementMapper().map(*elemIt);
- dest[globalI] = localResidual().residual(0);
- }
-
- // calculate the square norm of the residual
- Scalar result2 = dest.two_norm2();
- if (gridView_().comm().size() > 1)
- result2 = gridView_().comm().sum(result2);
-
- return std::sqrt(result2);
- }
-
- /*!
- * \brief Compute the integral over the domain of the storage
- * terms of all conservation quantities.
- *
- * \param dest Stores the result
- */
- void globalStorage(PrimaryVariables &dest)
- {
- dest = 0;
-
- ElementIterator elemIt = gridView_().template begin<0>();
- const ElementIterator elemEndIt = gridView_().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- localResidual().evalStorage(*elemIt);
-
- dest += localResidual().storageTerm()[0];
- };
-
- if (gridView_().comm().size() > 1)
- dest = gridView_().comm().sum(dest);
- }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- const SolutionVector &curSol() const
- { return uCur_; }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- SolutionVector &curSol()
- { return uCur_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- const SolutionVector &prevSol() const
- { return uPrev_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- SolutionVector &prevSol()
- { return uPrev_; }
-
- /*!
- * \brief Returns the operator assembler for the global jacobian of
- * the problem.
- */
- JacobianAssembler &jacobianAssembler()
- { return *jacAsm_; }
-
- /*!
- * \copydoc jacobianAssembler()
- */
- const JacobianAssembler &jacobianAssembler() const
- { return *jacAsm_; }
-
- /*!
- * \brief Returns the local jacobian which calculates the local
- * stiffness matrix for an arbitrary element.
- *
- * The local stiffness matrices of the element are used by
- * the jacobian assembler to produce a global linerization of the
- * problem.
- */
- LocalJacobian &localJacobian()
- { return localJacobian_; }
- /*!
- * \copydoc localJacobian()
- */
- const LocalJacobian &localJacobian() const
- { return localJacobian_; }
-
- /*!
- * \brief Returns the local residual function.
- */
- LocalResidual &localResidual()
- { return localJacobian().localResidual(); }
- /*!
- * \copydoc localResidual()
- */
- const LocalResidual &localResidual() const
- { return localJacobian().localResidual(); }
-
- /*!
- * \brief Returns the relative weight of a primary variable for
- * calculating relative errors.
- *
- * \param dofIdx The global index of the control volume
- * \param pvIdx The index of the primary variable
- */
- Scalar primaryVarWeight(int dofIdx, int pvIdx) const
- {
- return 1.0/std::max(std::abs(this->prevSol()[dofIdx][pvIdx]), 1.0);
- }
-
- /*!
- * \brief Returns the relative error between two vectors of
- * primary variables.
- *
- * \param vertexIdx The global index of the control volume's
- * associated vertex
- * \param pv1 The first vector of primary variables
- * \param pv2 The second vector of primary variables
- *
- * \todo The vertexIdx argument is pretty hacky. it is required by
- * models with pseudo primary variables (i.e. the primary
- * variable switching models). the clean solution would be
- * to access the pseudo primary variables from the primary
- * variables.
- */
- Scalar relativeErrorVertex(int vertexIdx,
- const PrimaryVariables &pv1,
- const PrimaryVariables &pv2)
- {
- Scalar result = 0.0;
- for (int j = 0; j < numEq; ++j) {
- //Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
- //Scalar eqErr = std::abs(pv1[j] - pv2[j])*weight;
- Scalar eqErr = std::abs(pv1[j] - pv2[j]);
- eqErr /= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
-
- result = std::max(result, eqErr);
- }
- return result;
- }
-
- /*!
- * \brief Try to progress the model to the next timestep.
- *
- * \param solver The non-linear solver
- * \param controller The controller which specifies the behaviour
- * of the non-linear solver
- */
- bool update(NewtonMethod &solver,
- NewtonController &controller)
- {
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i)
- Valgrind::CheckDefined(curSol()[i]);
-#endif // HAVE_VALGRIND
-
- asImp_().updateBegin();
-
- bool converged = solver.execute(controller);
- if (converged) {
- asImp_().updateSuccessful();
- }
- else
- asImp_().updateFailed();
-
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i) {
- Valgrind::CheckDefined(curSol()[i]);
- }
-#endif // HAVE_VALGRIND
-
- return converged;
- }
-
-
- /*!
- * \brief Called by the update() method before it tries to
- * apply the newton method. This is primary a hook
- * which the actual model can overload.
- */
- void updateBegin()
- { }
-
-
- /*!
- * \brief Called by the update() method if it was
- * successful. This is primary a hook which the actual
- * model can overload.
- */
- void updateSuccessful()
- { };
-
- /*!
- * \brief Called by the update() method if it was
- * unsuccessful. This is primary a hook which the actual
- * model can overload.
- */
- void updateFailed()
- {
- // Reset the current solution to the one of the
- // previous time step so that we can start the next
- // update at a physically meaningful solution.
- uCur_ = uPrev_;
-
- jacAsm_->reassembleAll();
- };
-
- /*!
- * \brief Called by the problem if a time integration was
- * successful, post processing of the solution is done and
- * the result has been written to disk.
- *
- * This should prepare the model for the next time integration.
- */
- void advanceTimeLevel()
- {
- // make the current solution the previous one.
- uPrev_ = uCur_;
-
- updatePrevHints();
- }
-
- /*!
- * \brief Serializes the current state of the model.
- *
- * \tparam Restarter The type of the serializer class
- *
- * \param res The serializer object
- */
- template <class Restarter>
- void serialize(Restarter &res)
- { res.template serializeEntities<0>(asImp_(), this->gridView_()); }
-
- /*!
- * \brief Deserializes the state of the model.
- *
- * \tparam Restarter The type of the serializer class
- *
- * \param res The serializer object
- */
- template <class Restarter>
- void deserialize(Restarter &res)
- {
- res.template deserializeEntities<0>(asImp_(), this->gridView_());
- prevSol() = curSol();
- }
-
- /*!
- * \brief Write the current solution for a vertex to a restart
- * file.
- *
- * \param outstream The stream into which the vertex data should
- * be serialized to
- * \param vert The DUNE Codim<dim> entity which's data should be
- * serialized
- */
- void serializeEntity(std::ostream &outstream,
- const Element &element)
- {
- int elemIdx = dofMapper().map(element);
-
- // write phase state
- if (!outstream.good()) {
- DUNE_THROW(Dune::IOError,
- "Could not serialize element "
- << elemIdx);
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- outstream << curSol()[elemIdx][eqIdx] << " ";
- }
- }
-
- /*!
- * \brief Reads the current solution variables for a vertex from a
- * restart file.
- *
- * \param instream The stream from which the vertex data should
- * be deserialized from
- * \param vert The DUNE Codim<dim> entity which's data should be
- * deserialized
- */
- void deserializeEntity(std::istream &instream,
- const Element &element)
- {
- int elemIdx = dofMapper().map(element);
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (!instream.good())
- DUNE_THROW(Dune::IOError,
- "Could not deserialize element "
- << elemIdx);
- instream >> curSol()[elemIdx][eqIdx];
- }
- }
-
- /*!
- * \brief Returns the number of global degrees of freedoms (DOFs)
- */
- size_t numDofs() const
- { return gridView_().size(0); }
-
- /*!
- * \brief Mapper for the entities where degrees of freedoms are
- * defined to indices.
- *
- * Here, this means a mapper for elements.
- */
- const DofMapper &dofMapper() const
- { return problem_().elementMapper(); };
-
- /*!
- * \brief Mapper for vertices to indices.
- */
- const VertexMapper &vertexMapper() const
- { return problem_().vertexMapper(); };
-
- /*!
- * \brief Mapper for elements to indices.
- */
- const ElementMapper &elementMapper() const
- { return problem_().elementMapper(); };
-
- /*!
- * \brief Resets the Jacobian matrix assembler, so that the
- * boundary types can be altered.
- */
- void resetJacobianAssembler ()
- {
- delete jacAsm_;
- jacAsm_ = new JacobianAssembler;
- jacAsm_->init(problem_());
- }
-
- /*!
- * \brief Update the weights of all primary variables within an
- * element given the complete set of volume variables
- *
- * \param element The DUNE codim 0 entity
- * \param volVars All volume variables for the element
- */
- void updatePVWeights(const Element &element,
- const ElementVolumeVariables &volVars) const
- { };
-
- /*!
- * \brief Add the vector fields for analysing the convergence of
- * the newton method to the a VTK multi writer.
- *
- * \tparam MultiWriter The type of the VTK multi writer
- *
- * \param writer The VTK multi writer object on which the fields should be added.
- * \param u The solution function
- * \param deltaU The delta of the solution function before and after the Newton update
- */
- template <class MultiWriter>
- void addConvergenceVtkFields(MultiWriter &writer,
- const SolutionVector &u,
- const SolutionVector &deltaU)
- {}
-
- /*!
- * \brief Add the quantities of a time step which ought to be written to disk.
- *
- * This should be overwritten by the actual model if any secondary
- * variables should be written out. Read: This should _always_ be
- * overwritten by well behaved models!
- *
- * \tparam MultiWriter The type of the VTK multi writer
- *
- * \param sol The global vector of primary variable values.
- * \param writer The VTK multi writer where the fields should be added.
- */
- template <class MultiWriter>
- void addOutputVtkFields(const SolutionVector &sol,
- MultiWriter &writer)
- {
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
-
- // create the required scalar fields
- unsigned numElements = this->gridView_().size(0);
-
- // global defect of the two auxiliary equations
- ScalarField* x[numEq];
- for (int i = 0; i < numEq; ++i) {
- x[i] = writer.allocateManagedBuffer(numElements);
- }
-
- ElementIterator eIt = this->gridView_().template begin<0>();
- ElementIterator eEndIt = this->gridView_().template end<0>();
- for (; eIt != eEndIt; ++ eIt)
- {
- int globalIdx = elementMapper().map(*eIt);
- for (int i = 0; i < numEq; ++i) {
- (*x[i])[globalIdx] = sol[globalIdx][i];
- }
- }
-
- for (int i = 0; i < numEq; ++i) {
- std::ostringstream oss;
- oss << "primaryVar_" << i;
- writer.attachCellData(*x[i], oss.str());
- }
- }
-
- /*!
- * \brief Reference to the grid view of the spatial domain.
- */
- const GridView &gridView() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns true if the control volume touches
- * the grid's boundary.
- *
- * \param globalIdx The global index of the control volume
- */
- bool onBoundary(int globalIdx) const
- { return boundaryIndices_[globalIdx]; }
-
- /*!
- * \brief Returns true if the control volume touches
- * the grid's boundary.
- *
- * \param elem A DUNE Codim<0> entity coinciding with the control
- * volume.
- */
- bool onBoundary(const Element &elem) const
- { return onBoundary(elementMapper().map(elem)); }
-
- /*!
- * \brief Fill the fluid state according to the primary variables.
- *
- * Taking the information from the primary variables,
- * the fluid state is filled with every information that is
- * necessary to evaluate the model's local residual.
- *
- * \param primaryVariables The primary variables of the model.
- * \param problem The problem at hand.
- * \param element The current element.
- * \param elementGeometry The finite volume element geometry.
- * \param scvIdx The index of the subcontrol volume.
- * \param fluidState The fluid state to fill.
- */
- template <class FluidState>
- static void completeFluidState(const PrimaryVariables& primaryVariables,
- const Problem& problem,
- const Element& element,
- const FVElementGeometry& elementGeometry,
- FluidState& fluidState)
- {
- VolumeVariables::completeFluidState(primaryVariables, problem, element,
- elementGeometry, fluidState);
- }
-protected:
- /*!
- * \brief A reference to the problem on which the model is applied.
- */
- Problem &problem_()
- { return *problemPtr_; }
- /*!
- * \copydoc problem_()
- */
- const Problem &problem_() const
- { return *problemPtr_; }
-
- /*!
- * \brief Reference to the grid view of the spatial domain.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Reference to the local residal object
- */
- LocalResidual &localResidual_()
- { return localJacobian_.localResidual(); }
-
- /*!
- * \brief Applies the initial solution for all vertices of the grid.
- */
- void applyInitialSolution_()
- {
- // first set the whole domain to zero
- uCur_ = Scalar(0.0);
-
- FVElementGeometry fvElemGeom;
-
- // iterate through leaf grid and evaluate initial
- // condition at the center of each sub control volume
- //
- // TODO: the initial condition needs to be unique for
- // each vertex. we should think about the API...
- ElementIterator it = gridView_().template begin<0>();
- const ElementIterator &eendit = gridView_().template end<0>();
- for (; it != eendit; ++it) {
- // deal with the current element
- fvElemGeom.update(gridView_(), *it);
-
- // map the local vertex index to the global one
- int globalIdx = elementMapper().map(*it);
-
- // let the problem do the dirty work of nailing down
- // the initial solution.
- PrimaryVariables initVal;
- Valgrind::SetUndefined(initVal);
- problem_().initial(initVal,
- *it,
- fvElemGeom,
- 0);
- Valgrind::CheckDefined(initVal);
-
- uCur_[globalIdx] = initVal;
- Valgrind::CheckDefined(uCur_[globalIdx]);
- }
- }
-
- /*!
- * \brief Find all indices of boundary vertices.
- *
- * For this we need to loop over all intersections (which is slow
- * in general). If the DUNE grid interface would provide a
- * onBoundary() method for entities this could be done in a much
- * nicer way (actually this would not be necessary)
- */
- void updateBoundaryIndices_()
- {
- boundaryIndices_.resize(numDofs());
- std::fill(boundaryIndices_.begin(), boundaryIndices_.end(), false);
-
- ElementIterator eIt = gridView_().template begin<0>();
- ElementIterator eEndIt = gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt) {
- IntersectionIterator isIt = gridView_().ibegin(*eIt);
- IntersectionIterator isEndIt = gridView_().iend(*eIt);
- for (; isIt != isEndIt; ++isIt) {
- if (!isIt->boundary())
- continue;
-
- int globalIdx = elementMapper().map(*eIt);
- boundaryIndices_[globalIdx] = true;
- }
- }
- }
-
- // the problem we want to solve. defines the constitutive
- // relations, matxerial laws, etc.
- Problem *problemPtr_;
-
- // calculates the local jacobian matrix for a given element
- LocalJacobian localJacobian_;
- // Linearizes the problem at the current time step using the
- // local jacobian
- JacobianAssembler *jacAsm_;
-
- // the set of all indices of vertices on the boundary
- std::vector<bool> boundaryIndices_;
-
- // cur is the current iterative solution, prev the converged
- // solution of the previous time step
- SolutionVector uCur_;
- SolutionVector uPrev_;
-
-private:
- /*!
- * \brief Returns whether messages should be printed
- */
- bool verbose_() const
- { return gridView_().comm().rank() == 0; };
-
- Implementation &asImp_()
- { return *static_cast<Implementation*>(this); }
- const Implementation &asImp_() const
- { return *static_cast<const Implementation*>(this); }
-
- bool enableHints_;
-};
-}
-
-#endif
diff --git a/dumux/implicit/cellcentered/ccpropertydefaults.hh b/dumux/implicit/cellcentered/ccpropertydefaults.hh
index 57904da7ed..3f63566d08 100644
--- a/dumux/implicit/cellcentered/ccpropertydefaults.hh
+++ b/dumux/implicit/cellcentered/ccpropertydefaults.hh
@@ -33,7 +33,6 @@
#include <dumux/implicit/common/implicitpropertydefaults.hh>
#include "ccassembler.hh"
-#include "ccmodel.hh"
#include "ccfvelementgeometry.hh"
#include "ccelementboundarytypes.hh"
#include "cclocalresidual.hh"
@@ -59,9 +58,6 @@ SET_TYPE_PROP(CCModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper)
//! Set the BaseLocalResidual to CCLocalResidual
SET_TYPE_PROP(CCModel, BaseLocalResidual, Dumux::CCLocalResidual<TypeTag>);
-//! Set the BaseModel to CCModel
-SET_TYPE_PROP(CCModel, BaseModel, Dumux::CCModel<TypeTag>);
-
//! An array of secondary variable containers
SET_TYPE_PROP(CCModel, ElementVolumeVariables, Dumux::CCElementVolumeVariables<TypeTag>);
diff --git a/dumux/implicit/common/implicitassembler.hh b/dumux/implicit/common/implicitassembler.hh
index 5ca3d76c91..9c731fe616 100644
--- a/dumux/implicit/common/implicitassembler.hh
+++ b/dumux/implicit/common/implicitassembler.hh
@@ -290,9 +290,9 @@ public:
// we need to add the distance the solution was moved for
// this vertex
- Scalar dist = model_().relativeErrorVertex(i,
- currentPriVars,
- nextPriVars);
+ Scalar dist = model_().relativeErrorDof(i,
+ currentPriVars,
+ nextPriVars);
delta_[i] += std::abs(dist);
}
diff --git a/dumux/implicit/common/implicitmodel.hh b/dumux/implicit/common/implicitmodel.hh
index 087dbf778a..842e9e910e 100644
--- a/dumux/implicit/common/implicitmodel.hh
+++ b/dumux/implicit/common/implicitmodel.hh
@@ -24,17 +24,13 @@
#ifndef DUMUX_IMPLICIT_MODEL_HH
#define DUMUX_IMPLICIT_MODEL_HH
-#include "implicitproperties.hh"
-#include "implicitpropertydefaults.hh"
-
-#include "implicitelementvolumevariables.hh"
-#include "implicitlocaljacobian.hh"
-#include "implicitlocalresidual.hh"
+#include <dune/grid/common/geometry.hh>
+#include <dune/istl/bvector.hh>
+#include "implicitproperties.hh"
+#include <dumux/common/valgrind.hh>
#include <dumux/parallel/vertexhandles.hh>
-#include <dune/grid/common/geometry.hh>
-
namespace Dumux
{
@@ -246,7 +242,7 @@ public:
else
{
int globalI = elementMapper().map(*elemIt);
- dest[globalI] = localResidual().residual(0);
+ residual[globalI] = localResidual().residual(0);
}
}
@@ -311,7 +307,7 @@ public:
}
else
{
- DUNE_THROW(Dune::InvalidStateError,
+ DUNE_THROW(Dune::InvalidStateException,
"requested box volume for cell-centered model");
}
}
@@ -545,7 +541,7 @@ public:
const Entity &entity)
{
int dofIdx = dofMapper().map(entity);
-
+
// write phase state
if (!outstream.good()) {
DUNE_THROW(Dune::IOError,
@@ -573,6 +569,7 @@ public:
const Entity &entity)
{
int dofIdx = dofMapper().map(entity);
+
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!instream.good())
DUNE_THROW(Dune::IOError,
@@ -590,7 +587,7 @@ public:
if (isBox)
return gridView_().size(dim);
else
- return gridView_().size(dim);
+ return gridView_().size(0);
}
/*!
@@ -601,12 +598,15 @@ public:
* for vertices, if the cell centered method is used,
* this means a mapper for elements.
*/
- const DofMapper &dofMapper() const
+ template <class T = TypeTag>
+ const typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox), VertexMapper>::type &dofMapper() const
{
- if (isBox)
- return problem_().vertexMapper();
- else
- return problem_().elementMapper();
+ return problem_().vertexMapper();
+ }
+ template <class T = TypeTag>
+ const typename std::enable_if<!GET_PROP_VALUE(T, ImplicitIsBox), ElementMapper>::type &dofMapper() const
+ {
+ return problem_().elementMapper();
}
/*!
@@ -780,7 +780,7 @@ public:
if (isBox)
return onBoundary(vertexMapper().map(element, vIdx, dim));
else
- DUNE_THROW(Dune::InvalidStateError,
+ DUNE_THROW(Dune::InvalidStateException,
"requested for cell-centered model");
}
@@ -797,7 +797,7 @@ public:
if (!isBox)
return onBoundary(elementMapper().map(elem));
else
- DUNE_THROW(Dune::InvalidStateError,
+ DUNE_THROW(Dune::InvalidStateException,
"requested for box model");
}
@@ -873,12 +873,20 @@ protected:
fvGeometry.update(gridView_(), *eIt);
// loop over all element vertices, i.e. sub control volumes
- for (int scvIdx = 0; scvIdx < fvGeometry.numVertices; scvIdx++)
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; scvIdx++)
{
- // map the local vertex index to the global one
- int globalIdx = vertexMapper().map(*eIt,
- scvIdx,
- dim);
+ int globalIdx;
+
+ if (isBox)
+ {
+ // map the local vertex index to the global one
+ globalIdx = vertexMapper().map(*eIt, scvIdx, dim);
+ }
+ else
+ {
+ // get the global index of the element
+ globalIdx = elementMapper().map(*eIt);
+ }
// let the problem do the dirty work of nailing down
// the initial solution.
@@ -890,12 +898,16 @@ protected:
scvIdx);
Valgrind::CheckDefined(initPriVars);
- // add up the initial values of all sub-control
- // volumes. If the initial values disagree for
- // different sub control volumes, the initial value
- // will be the arithmetic mean.
- initPriVars *= fvGeometry.subContVol[scvIdx].volume;
- boxVolume_[globalIdx] += fvGeometry.subContVol[scvIdx].volume;
+ if (isBox)
+ {
+ // add up the initial values of all sub-control
+ // volumes. If the initial values disagree for
+ // different sub control volumes, the initial value
+ // will be the arithmetic mean.
+ initPriVars *= fvGeometry.subContVol[scvIdx].volume;
+ boxVolume_[globalIdx] += fvGeometry.subContVol[scvIdx].volume;
+ }
+
uCur_[globalIdx] += initPriVars;
Valgrind::CheckDefined(uCur_[globalIdx]);
}
@@ -903,7 +915,7 @@ protected:
// add up the primary variables and the volumes of the boxes
// which cross process borders
- if (gridView_().comm().size() > 1) {
+ if (isBox && gridView_().comm().size() > 1) {
VertexHandleSum<Dune::FieldVector<Scalar, 1>,
Dune::BlockVector<Dune::FieldVector<Scalar, 1> >,
VertexMapper> sumVolumeHandle(boxVolume_, vertexMapper());
@@ -918,20 +930,17 @@ protected:
Dune::ForwardCommunication);
}
- // divide all primary variables by the volume of their boxes
- int n = gridView_().size(dim);
- for (int i = 0; i < n; ++i) {
- uCur_[i] /= boxVolume(i);
+ if (isBox)
+ {
+ // divide all primary variables by the volume of their boxes
+ for (int i = 0; i < uCur_.size(); ++i) {
+ uCur_[i] /= boxVolume(i);
+ }
}
}
/*!
- * \brief Find all indices of boundary vertices.
- *
- * For this we need to loop over all intersections (which is slow
- * in general). If the DUNE grid interface would provide a
- * onBoundary() method for entities this could be done in a much
- * nicer way (actually this would not be necessary)
+ * \brief Find all indices of boundary vertices (box) / elements (cell centered).
*/
void updateBoundaryIndices_()
{
@@ -948,20 +957,28 @@ protected:
IntersectionIterator isEndIt = gridView_().iend(*eIt);
for (; isIt != isEndIt; ++isIt) {
if (isIt->boundary()) {
- // add all vertices on the intersection to the set of
- // boundary vertices
- int faceIdx = isIt->indexInInside();
- int numFaceVerts = refElement.size(faceIdx, 1, dim);
- for (int faceVertIdx = 0;
- faceVertIdx < numFaceVerts;
- ++faceVertIdx)
+ if (isBox)
{
- int elemVertIdx = refElement.subEntity(faceIdx,
- 1,
- faceVertIdx,
- dim);
- int globalVertIdx = vertexMapper().map(*eIt, elemVertIdx, dim);
- boundaryIndices_[globalVertIdx] = true;
+ // add all vertices on the intersection to the set of
+ // boundary vertices
+ int faceIdx = isIt->indexInInside();
+ int numFaceVerts = refElement.size(faceIdx, 1, dim);
+ for (int faceVertIdx = 0;
+ faceVertIdx < numFaceVerts;
+ ++faceVertIdx)
+ {
+ int elemVertIdx = refElement.subEntity(faceIdx,
+ 1,
+ faceVertIdx,
+ dim);
+ int globalVertIdx = vertexMapper().map(*eIt, elemVertIdx, dim);
+ boundaryIndices_[globalVertIdx] = true;
+ }
+ }
+ else
+ {
+ int globalIdx = elementMapper().map(*eIt);
+ boundaryIndices_[globalIdx] = true;
}
}
}
@@ -1010,4 +1027,6 @@ private:
};
}
+#include "implicitpropertydefaults.hh"
+
#endif
diff --git a/dumux/implicit/common/implicitproblem.hh b/dumux/implicit/common/implicitproblem.hh
index 492ffd3d34..aa499baa09 100644
--- a/dumux/implicit/common/implicitproblem.hh
+++ b/dumux/implicit/common/implicitproblem.hh
@@ -24,6 +24,7 @@
#define DUMUX_IMPLICIT_PROBLEM_HH
#include "implicitproperties.hh"
+#include "implicitmodel.hh"
#include <dumux/io/vtkmultiwriter.hh>
#include <dumux/io/restart.hh>
diff --git a/dumux/implicit/common/implicitpropertydefaults.hh b/dumux/implicit/common/implicitpropertydefaults.hh
index 032ab4faa5..5076a3c27e 100644
--- a/dumux/implicit/common/implicitpropertydefaults.hh
+++ b/dumux/implicit/common/implicitpropertydefaults.hh
@@ -33,6 +33,7 @@
#include <dumux/common/timemanager.hh>
#include "implicitproperties.hh"
+#include "implicitmodel.hh"
#include "implicitlocaljacobian.hh"
#include "implicitvolumevariables.hh"
@@ -73,6 +74,9 @@ SET_TYPE_PROP(ImplicitBase,
Dune::MultipleCodimMultipleGeomTypeMapper<typename GET_PROP_TYPE(TypeTag, GridView),
Dune::MCMGElementLayout>);
+//! Set the BaseModel to ImplicitModel
+SET_TYPE_PROP(ImplicitBase, BaseModel, Dumux::ImplicitModel<TypeTag>);
+
//! The volume variable class, to be overloaded by the model
SET_TYPE_PROP(ImplicitBase, VolumeVariables, Dumux::ImplicitVolumeVariables<TypeTag>);
diff --git a/dumux/implicit/mpnc/mpnclocalresidual.hh b/dumux/implicit/mpnc/mpnclocalresidual.hh
index c2623eb26c..5513b44529 100644
--- a/dumux/implicit/mpnc/mpnclocalresidual.hh
+++ b/dumux/implicit/mpnc/mpnclocalresidual.hh
@@ -24,7 +24,7 @@
#include "energy/mpnclocalresidualenergy.hh"
#include "mass/mpnclocalresidualmass.hh"
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>
namespace Dumux
diff --git a/dumux/implicit/mpnc/mpncmodel.hh b/dumux/implicit/mpnc/mpncmodel.hh
index 5b31ce8871..e8ba3bb925 100644
--- a/dumux/implicit/mpnc/mpncmodel.hh
+++ b/dumux/implicit/mpnc/mpncmodel.hh
@@ -22,7 +22,7 @@
#include "mpncproperties.hh"
#include "mpncvtkwriter.hh"
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
namespace Dumux
{
diff --git a/dumux/implicit/richards/richardsmodel.hh b/dumux/implicit/richards/richardsmodel.hh
index c66248d2b4..ff71be4f6f 100644
--- a/dumux/implicit/richards/richardsmodel.hh
+++ b/dumux/implicit/richards/richardsmodel.hh
@@ -25,7 +25,7 @@
#ifndef DUMUX_RICHARDS_MODEL_HH
#define DUMUX_RICHARDS_MODEL_HH
-#include <dumux/implicit/box/boxmodel.hh>
+#include <dumux/boxmodels/common/boxmodel.hh>
#include "richardslocalresidual.hh"
#include "richardsproblem.hh"
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index b6d3ca41b4..987bbb3faa 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -346,10 +346,10 @@ public:
PrimaryVariables uNewI = uLastIter[i];
uNewI -= deltaU[i];
- Scalar vertError = model_().relativeErrorVertex(i,
- uLastIter[i],
- uNewI);
- error_ = std::max(error_, vertError);
+ Scalar dofError = model_().relativeErrorDof(i,
+ uLastIter[i],
+ uNewI);
+ error_ = std::max(error_, dofError);
}
diff --git a/test/implicit/2p/lensproblem.hh b/test/implicit/2p/lensproblem.hh
index b9c19b5205..0f95a12ab5 100644
--- a/test/implicit/2p/lensproblem.hh
+++ b/test/implicit/2p/lensproblem.hh
@@ -255,7 +255,7 @@ public:
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
- values = 0;
+ values = 0.0;
}
// \}
--
GitLab