From 9021fdbb331386aec55c64963131db166df9bf3a Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Thu, 2 Feb 2012 12:38:56 +0000
Subject: [PATCH] revert r6411 with some modifications
the reason why the boundary indices where a global array is because it
is possible that an element's vertex is on the boundary, but the
element does not have any boundary intersections. only checking for
boundary intersections does not work for dirichlet conditions on such
grids. (for simplex grids that case is actually quite common.)
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7605 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../common/boxelementboundarytypes.hh | 47 +--------
dumux/boxmodels/common/boxmodel.hh | 97 ++++++++++++++++---
2 files changed, 88 insertions(+), 56 deletions(-)
diff --git a/dumux/boxmodels/common/boxelementboundarytypes.hh b/dumux/boxmodels/common/boxelementboundarytypes.hh
index 44998ebbf5..c3e229974b 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 2463f7e1d2..73621ebf85 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_;
};
}
--
GitLab