diff --git a/dumux/boxmodels/common/boxelementboundarytypes.hh b/dumux/boxmodels/common/boxelementboundarytypes.hh
index 44998ebbf5042548f657495b91aa0182ccf8e8e6..c3e229974bff92ab3b2b5cd8f681e3f02fa95f67 100644
--- a/dumux/boxmodels/common/boxelementboundarytypes.hh
+++ b/dumux/boxmodels/common/boxelementboundarytypes.hh
@@ -29,8 +29,6 @@
#include "boxproperties.hh"
-#include <dune/grid/common/geometry.hh>
-
#include <dumux/common/valgrind.hh>
namespace Dumux
@@ -55,11 +53,6 @@ class BoxElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTa
enum { dim = GridView::dimension };
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
- typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef typename GridView::ctype CoordScalar;
-
- typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
public:
/*!
@@ -97,50 +90,14 @@ public:
int numVerts = element.template count<dim>();
this->resize(numVerts);
- int nBoundary = 0;
- std::vector<bool> onBoundary(numVerts, false);
-
- // loop over all intersections of the element and mark all
- // vertices in these intersections
- Dune::GeometryType geoType = element.geometry().type();
- const ReferenceElement &refElem = ReferenceElements::general(geoType);
- IntersectionIterator isIt = problem.gridView().ibegin(element);
- IntersectionIterator isEndIt = problem.gridView().iend(element);
- for (; isIt != isEndIt; ++isIt) {
- if (!isIt->boundary())
- continue; // intersection is not on grid boundary
-
- // mark all vertices on the intersection
- 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);
- if (!onBoundary[elemVertIdx]) {
- ++ nBoundary;
- onBoundary[elemVertIdx] = true;
- }
- }
- }
-
hasDirichlet_ = false;
hasNeumann_ = false;
hasOutflow_ = false;
- if (nBoundary == 0) {
- for (int i = 0; i < numVerts; ++i)
- (*this)[i].reset();
- return;
- }
for (int i = 0; i < numVerts; ++i) {
(*this)[i].reset();
- if (!onBoundary[i])
+ if (!problem.model().onBoundary(element, i))
continue;
const VertexPointer vptr = element.template subEntity<dim>(i);
@@ -150,7 +107,7 @@ public:
hasNeumann_ = hasNeumann_ || (*this)[i].hasNeumann();
hasOutflow_ = hasOutflow_ || (*this)[i].hasOutflow();
}
- };
+ }
/*!
* \brief Returns whether the element has a vertex which contains
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index 2463f7e1d202f8ea236d5696c153fa4cbcf1035f..73621ebf8569464701d27340b7cc491127ce5df2 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -37,6 +37,7 @@
#include <dumux/parallel/vertexhandles.hh>
+#include <dune/grid/common/geometry.hh>
namespace Dumux
{
@@ -74,10 +75,15 @@ class BoxModel
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 &);
@@ -104,6 +110,8 @@ public:
{
problemPtr_ = &prob;
+ updateBoundaryIndices_();
+
int nDofs = asImp_().numDofs();
uCur_.resize(nDofs);
uPrev_.resize(nDofs);
@@ -705,6 +713,27 @@ public:
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(int globalVertIdx) const
+ { return boundaryIndices_[globalVertIdx]; }
+
+ /*!
+ * \brief Returns true if a vertex is located on the grid's
+ * boundary.
+ *
+ * \param elem A DUNE Codim<0> entity which contains the control
+ * volume's associated vertex.
+ * \param vIdx The local vertex index inside elem
+ */
+ bool onBoundary(const Element &elem, int vIdx) const
+ { return onBoundary(vertexMapper().map(elem, vIdx, dim)); }
+
template <class FluidState>
static void completeFluidState(const PrimaryVariables& primaryVariables,
const Problem& problem,
@@ -815,23 +844,55 @@ protected:
}
}
+ /*!
+ * \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 &refElem = ReferenceElements::general(geoType);
+
+ IntersectionIterator isIt = gridView_().ibegin(*eIt);
+ IntersectionIterator isEndIt = gridView_().iend(*eIt);
+ for (; isIt != isEndIt; ++isIt) {
+ if (!isIt->boundary())
+ continue;
+ // add all vertices on the intersection to the set of
+ // boundary vertices
+ 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 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_;
- /*!
- * \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); }
-
// the problem we want to solve. defines the constitutive
// relations, matxerial laws, etc.
Problem *problemPtr_;
@@ -842,6 +903,9 @@ protected:
// 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_;
@@ -850,6 +914,17 @@ protected:
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_;
};
}