diff --git a/test/discretization/cellcentered/tpfa/CMakeLists.txt b/test/discretization/cellcentered/tpfa/CMakeLists.txt
index 56fa30fedeaf0377736cecef3c0a15073ddd6893..33625fcfb58f201a6fe632108c8513319d7f6bf4 100644
--- a/test/discretization/cellcentered/tpfa/CMakeLists.txt
+++ b/test/discretization/cellcentered/tpfa/CMakeLists.txt
@@ -1,2 +1,13 @@
-add_dumux_test(test_tpfafvgeometry test_tpfafvgeometry test_tpfafvgeometry.cc
- ${CMAKE_CURRENT_BINARY_DIR}/test_tpfafvgeometry)
+dune_add_test(NAME test_tpfafvgeometry
+ SOURCES test_tpfafvgeometry.cc
+ COMMAND ./test_tpfafvgeometry)
+
+dune_add_test(NAME test_tpfafvgeometry_nonconforming
+ SOURCES test_tpfafvgeometry_nonconforming.cc
+ COMPILE_DEFINITIONS TYPETAG=TestFVGeometryNonConforming
+ COMMAND ./test_tpfafvgeometry_nonconforming)
+
+dune_add_test(NAME test_cachedtpfafvgeometry_nonconforming
+ SOURCES test_tpfafvgeometry_nonconforming.cc
+ COMPILE_DEFINITIONS TYPETAG=TestCachedFVGeometryNonConforming
+ COMMAND ./test_cachedtpfafvgeometry_nonconforming)
diff --git a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
index ca43301d6f6050ba6ff2be7a7ad35618024d1c25..e2cc9c3870a026247bb81987645118b91d74028a 100644
--- a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
+++ b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
@@ -31,20 +31,28 @@
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/common/mcmgmapper.hh>
-#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/common/basicproperties.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
-namespace Dumux
+//! Dummy flux variables class so that we can update the connectivity map
+class MockFluxVariables
{
+public:
+ template<class Element, class FvGeometry, class Scvf>
+ static std::vector<std::size_t> computeStencil(const Element& element,
+ const FvGeometry& fvGeometry,
+ const Scvf& scvf)
+ {
+ return std::vector<std::size_t>();
+ }
+};
-namespace Properties
-{
-NEW_TYPE_TAG(TestFVGeometry, INHERITS_FROM(CCTpfaModel));
-
-SET_TYPE_PROP(TestFVGeometry, Grid, Dune::YaspGrid<2>);
-
-SET_TYPE_PROP(TestFVGeometry, Problem, Dumux::MockProblem<TypeTag>);
-}
-
+namespace Dumux {
+ namespace Properties {
+ NEW_TYPE_TAG(TestFVGeometry, INHERITS_FROM(CCTpfaModel, NumericModel));
+ SET_TYPE_PROP(TestFVGeometry, Grid, Dune::YaspGrid<2>);
+ SET_TYPE_PROP(TestFVGeometry, FluxVariables, MockFluxVariables);
+ }
}
template<class T>
@@ -74,26 +82,24 @@ int main (int argc, char *argv[]) try
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
// make a grid
GlobalPosition lower(0.0);
GlobalPosition upper(1.0);
std::array<unsigned int, dim> els{{2, 2}};
std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els);
- auto leafGridView = grid->leafGridView();
- Problem problem(leafGridView);
-
- FVGridGeometry global(leafGridView);
- global.update(problem);
+ // obtain leaf and make FVGridGeometry
+ auto leafGridView = grid->leafGridView();
+ FVGridGeometry fvGridGeometry(leafGridView);
+ fvGridGeometry.update();
// iterate over elements. For every element get fv geometry and loop over scvs and scvfaces
for (const auto& element : elements(leafGridView))
{
- auto eIdx = problem.elementMapper().index(element);
+ auto eIdx = fvGridGeometry.elementMapper().index(element);
std::cout << std::endl << "Checking fvGeometry of element " << eIdx << std::endl;
- auto fvGeometry = localView(global);
+ auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bind(element);
auto range = scvs(fvGeometry);
diff --git a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5317e7c34613bc0679eeb669f56941c64faaebf1
--- /dev/null
+++ b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc
@@ -0,0 +1,413 @@
+// -*- 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 Test for the grid finite volume element geometry for non-conforming grids.
+ * A square grid with 9 elements is created of which the central element is refined.
+ * Subsequently, all directions & connectivity of the scvfs are checked for correctness.
+ */
+#include <config.h>
+
+#include <iostream>
+#include <utility>
+
+#include <dune/common/float_cmp.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/alugrid/grid.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dumux/common/basicproperties.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+//! Dummy flux variables class so that we can update the connectivity map
+class MockFluxVariables
+{
+public:
+ template<class Element, class FvGeometry, class Scvf>
+ static std::vector<std::size_t> computeStencil(const Element& element,
+ const FvGeometry& fvGeometry,
+ const Scvf& scvf)
+ {
+ return std::vector<std::size_t>();
+ }
+};
+
+namespace Dumux {
+ namespace Properties{
+ //! Test without using global caching of the geometries
+ NEW_TYPE_TAG(TestFVGeometryNonConforming, INHERITS_FROM(CCTpfaModel, NumericModel));
+ SET_TYPE_PROP(TestFVGeometryNonConforming, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
+ SET_TYPE_PROP(TestFVGeometryNonConforming, FluxVariables, MockFluxVariables);
+
+ //! Test using global geometry caching
+ NEW_TYPE_TAG(TestCachedFVGeometryNonConforming, INHERITS_FROM(CCTpfaModel, NumericModel));
+ SET_TYPE_PROP(TestCachedFVGeometryNonConforming, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
+ SET_BOOL_PROP(TestCachedFVGeometryNonConforming, EnableFVGridGeometryCache, true);
+ SET_TYPE_PROP(TestCachedFVGeometryNonConforming, FluxVariables, MockFluxVariables);
+ }
+}
+
+//! epsilon for checking direction of scvf normals
+constexpr double eps = 1e-6;
+
+//! returns true if a given position is inside the central element (before refinement)
+bool isInCentralElement(const Dune::FieldVector<double, 2>& pos)
+{
+ return pos[0] < 2.0 && pos[0] > 1.0 && pos[1] > 1.0 && pos[1] < 2.0;
+}
+
+//! returns if the element with given center position is corner element
+bool isCornerElement(const Dune::FieldVector<double, 2>& center)
+{
+ const auto distVec = center - Dune::FieldVector<double, 2>({1.5, 1.5});
+ return (abs(distVec[0]) > 1.0 - eps && abs(distVec[1]) > 1.0 - eps);
+}
+
+//! returns if the element is middle left element
+bool isMiddleLeftElement(const Dune::FieldVector<double, 2>& center)
+{
+ const auto distVec = center - Dune::FieldVector<double, 2>({1.5, 1.5});
+ return distVec[0] < -1.0 + eps && abs(distVec[1]) < eps;
+}
+
+//! returns if the element is middle right element
+bool isMiddleRightElement(const Dune::FieldVector<double, 2>& center)
+{
+ const auto distVec = center - Dune::FieldVector<double, 2>({1.5, 1.5});
+ return distVec[0] > 1.0 - eps && abs(distVec[1]) < eps;
+}
+
+//! returns if the element is middle upper element
+bool isMiddleUpperElement(const Dune::FieldVector<double, 2>& center)
+{
+ const auto distVec = center - Dune::FieldVector<double, 2>({1.5, 1.5});
+ return distVec[1] > 1.0 - eps && abs(distVec[0]) < eps;
+}
+
+//! returns if the element is middle lower element
+bool isMiddleLowerElement(const Dune::FieldVector<double, 2>& center)
+{
+ const auto distVec = center - Dune::FieldVector<double, 2>({1.5, 1.5});
+ return distVec[1] < -1.0 + eps && abs(distVec[0]) < eps;
+}
+
+//! returns for a given element center the element type name used in this test
+std::string elementTypeName(const Dune::FieldVector<double, 2>& center)
+{
+ if (isInCentralElement(center))
+ return "Central element";
+ if (isCornerElement(center))
+ return "Corner element";
+ if (isMiddleLeftElement(center))
+ return "Middle left";
+ if (isMiddleRightElement(center))
+ return "Middle right";
+ if (isMiddleUpperElement(center))
+ return "Middle upper";
+ if (isMiddleLowerElement(center))
+ return "Middle lower";
+
+ DUNE_THROW(Dune::InvalidStateException, "Element center position could not be interpreted.");
+}
+
+int main (int argc, char *argv[]) try
+{
+ // maybe initialize mpi
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::cout << "Checking the FVGeometries, SCVs and SCV faces on a non-conforming grid" << std::endl;
+
+ //! aliases
+ using TypeTag = TTAG(TYPETAG);
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using GridView = typename Grid::LeafGridView;
+
+ constexpr int dim = GridView::dimension;
+ constexpr int dimworld = GridView::dimensionworld;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ //! make a grid
+ GlobalPosition lower(0.0);
+ GlobalPosition upper(3.0);
+ std::array<unsigned int, dim> els{{3, 3}};
+ std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els);
+
+ //! refine the central element once
+ auto leafGridView = grid->leafGridView();
+ for (const auto& element : elements(leafGridView))
+ if (isInCentralElement(element.geometry().center()))
+ grid->mark(1, element);
+ grid->preAdapt();
+ grid->adapt();
+ grid->postAdapt();
+
+ //! if the leaf now does not have 12 elements, something went wrong
+ if (leafGridView.size(0) != 12)
+ DUNE_THROW(Dune::InvalidStateException, "Refined grid does not have exactly 12 elements!");
+
+ //! instantiate and update fvGridGeometry
+ FVGridGeometry fvGridGeometry(leafGridView);
+ fvGridGeometry.update();
+
+ //! We should have constructed 12 scvfs
+ if (fvGridGeometry.numScv() != 12)
+ DUNE_THROW(Dune::InvalidStateException, "FvGridGeometry does not have exactly 12 scvs!");
+
+ //! We should have constructed 52 scvfs
+ if (fvGridGeometry.numScvf() != 52)
+ DUNE_THROW(Dune::InvalidStateException, "FvGridGeometry does not have exactly 52 scvfs!");
+
+ //! iterate over elements and check for each element the number of scvfs
+ for (const auto& element : elements(leafGridView))
+ {
+ auto fvGeometry = localView(fvGridGeometry);
+ fvGeometry.bind(element);
+
+ //! For the tpfa scheme there is always one scv per element
+ if (fvGeometry.numScv() != 1)
+ DUNE_THROW(Dune::InvalidStateException, "More than one scv found in an element!");
+
+ //! make sure the scv has the same center as the element
+ const auto elementCenter = element.geometry().center();
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto scvCenter = scv.center();
+ for (unsigned int dir = 0; dir < dim; ++dir)
+ if ( !Dune::FloatCmp::eq(elementCenter[dir], scvCenter[dir], eps) )
+ DUNE_THROW(Dune::InvalidStateException, "Element center - " << elementCenter << " - and scv center - " << scvCenter << " - do not coincide");
+ }
+
+ //! check if the right number of scvfs point in the right direction
+ std::size_t numScvfsPosX = 0;
+ std::size_t numScvfsNegX = 0;
+ std::size_t numScvfsPosY = 0;
+ std::size_t numScvfsNegY = 0;
+ std::size_t numScvfsTotal = 0;
+
+ //! Also, scheck how many neighbors on the two levels the element has
+ //! and store the neighboring centers for output
+ std::size_t numL0Neighbors = 0;
+ std::size_t numL1Neighbors = 0;
+ std::vector<GlobalPosition> neighborCenters;
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ const auto n = scvf.unitOuterNormal();
+
+ //! the normal vector should point either in x- or y-direction with length 1
+ using std::abs;
+ if (abs(n[0]) > eps && abs(n[1]) > eps)
+ DUNE_THROW(Dune::InvalidStateException, "Wrong unit outer normal vector");
+
+ //! Outer normal must be pointing outwards
+ if ( n*(scvf.center() - elementCenter) < 0 )
+ DUNE_THROW(Dune::InvalidStateException, "Normal vector does not point outwards");
+
+ //! the face should never have more the one neighbor
+ if (scvf.numOutsideScvs() > 1)
+ DUNE_THROW(Dune::InvalidStateException, "Scvf has more than one neighbor");
+
+ //! center must always be between the corners
+ const auto d1 = scvf.corner(0) - scvf.center();
+ const auto d2 = scvf.corner(1) - scvf.center();
+ if ( d1 * d2 >= 0 )
+ DUNE_THROW(Dune::InvalidStateException, "Center is not between the two corners");
+
+ //! count up faces depending on direction of n
+ if ( Dune::FloatCmp::eq(n[0], 1.0, eps) )
+ numScvfsPosX++;
+ if ( Dune::FloatCmp::eq(n[0], -1.0, eps) )
+ numScvfsNegX++;
+ if ( Dune::FloatCmp::eq(n[1], 1.0, eps) )
+ numScvfsPosY++;
+ if ( Dune::FloatCmp::eq(n[1], -1.0, eps) )
+ numScvfsNegY++;
+
+ //! keep track of total number of scvfs
+ numScvfsTotal++;
+
+ //! check levels of neighbors
+ if (!scvf.boundary())
+ {
+ const auto outsideElement = fvGridGeometry.element(scvf.outsideScvIdx());
+ const auto outsideCenter = outsideElement.geometry().center();
+
+ if (isInCentralElement(outsideCenter))
+ numL1Neighbors++;
+ else
+ numL0Neighbors++;
+
+ //! check if ipGlobal & area make sense
+ const auto distVec = scvf.ipGlobal() - elementCenter;
+ if ( isInCentralElement(outsideCenter) && !isInCentralElement(elementCenter) )
+ {
+ if ( Dune::FloatCmp::eq(distVec[0], 0.0, eps) || Dune::FloatCmp::eq(distVec[1], 0.0, eps) )
+ DUNE_THROW(Dune::InvalidStateException, "Vector from element center to ipGlobal must NOT be axis-parallel " <<
+ "when going from lower to higher level");
+ if ( !Dune::FloatCmp::eq(scvf.area(), 0.5, eps) )
+ DUNE_THROW(Dune::InvalidStateException, "Area of scvf towards level 1 element is not 0.5!");
+ }
+ else
+ {
+ if ( !Dune::FloatCmp::eq(distVec[0], 0.0, eps) && !Dune::FloatCmp::eq(distVec[1], 0.0, eps) )
+ DUNE_THROW(Dune::InvalidStateException, "Vector from element center to ipGlobal must be axis-parallel " <<
+ "when neighboring levels are identical or outside level is lower");
+
+ if ( !isInCentralElement(outsideCenter) && !isInCentralElement(elementCenter) )
+ {
+ if ( !Dune::FloatCmp::eq(scvf.area(), 1.0, eps) )
+ DUNE_THROW(Dune::InvalidStateException, "Area of scvf between level 0 element is not 1.0!");
+ }
+ else
+ {
+ if ( !Dune::FloatCmp::eq(scvf.area(), 0.5, eps) )
+ DUNE_THROW(Dune::InvalidStateException, "Area of scvf between different levels is not 0.5!");
+ }
+ }
+
+ //! store outside center for output
+ neighborCenters.push_back(outsideCenter);
+ }
+ }
+
+ //! make sure the number of found scvfs makes sense
+ const std::size_t sumScvfs = numScvfsPosX+numScvfsNegX+numScvfsPosY+numScvfsNegY;
+ if (numScvfsTotal != sumScvfs)
+ DUNE_THROW(Dune::InvalidStateException, "Number of total scvfs is not equal to number of faces in individual directions.\n"
+ << "Total number: " << numScvfsTotal << ", sum individual faces: " << sumScvfs);
+
+ //! print found combination in element to terminal
+ std::cout << elementTypeName(elementCenter) << " with element center " << elementCenter << " has " << numScvfsTotal << " faces:\n"
+ << " - " << numScvfsPosX << " in positive x-direction\n"
+ << " - " << numScvfsNegX << " in negative x-direction\n"
+ << " - " << numScvfsPosY << " in positive y-direction\n"
+ << " - " << numScvfsNegY << " in negative x-direction\n"
+ << "The neighboring element centers are:" << std::endl;
+ for (const auto& p : neighborCenters) std::cout << "\t\t - " << p << std::endl;
+
+ //! in elements within the original central element or corner elements we should have only 4 scvfs (1 in each direction)
+ using std::abs;
+ if ( isInCentralElement(elementCenter) || isCornerElement(elementCenter) )
+ {
+ if (numScvfsTotal != 4)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 4 scvfs in a central/corner element");
+ if (numScvfsPosX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in a central/corner element in positive x-direction");
+ if (numScvfsPosY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in a central/corner element in positive y-direction");
+ if (numScvfsNegX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in a central/corner element in negative x-direction");
+ if (numScvfsNegY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvf1 in a central/corner element in negative y-direction");
+
+ //! In the corners we should find two level 0 neighbor elements
+ if ( isCornerElement(elementCenter) && numL1Neighbors != 0 && numL0Neighbors != 2 )
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly two level 0 neighbors and zero level 1 neighbors in corner element");
+
+ //! In the center elements we should find two level 0 and 2 level 1 neighbor elements
+ if ( isInCentralElement(elementCenter) && numL1Neighbors != 2 && numL0Neighbors != 2 )
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly two level 0 neighbors and two level 1 neighbors in central element");
+ }
+ //! In the element on the left to the center we should have 2 scvfs in positive x-direction
+ else if ( isMiddleLeftElement(elementCenter) ) //! to the left of the center
+ {
+ if (numScvfsTotal != 5)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 5 scvfs in middle left element");
+ if (numScvfsPosX != 2)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 2 scvfs in middle left element in positive x-direction");
+ if (numScvfsPosY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle left element in positive y-direction");
+ if (numScvfsNegX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle left element in negative x-direction");
+ if (numScvfsNegY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle left element in negative y-direction");
+
+ //! we should find two level 0 and 2 level 1 neighbor elements
+ if (numL1Neighbors != 2 && numL0Neighbors != 2)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly two level 0 neighbors and two level 1 neighbors in middle left element");
+ }
+ //! In the element on the right to the center we should have 2 scvfs in negative x-direction
+ else if ( isMiddleRightElement(elementCenter) ) //! to the right of the center
+ {
+ if (numScvfsTotal != 5)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 5 scvfs in middle right element");
+ if (numScvfsPosX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle right element in positive x-direction");
+ if (numScvfsPosY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle right element in positive y-direction");
+ if (numScvfsNegX != 2)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 2 scvfs in middle right element in negative x-direction");
+ if (numScvfsNegY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle right element in negative y-direction");
+
+ //! we should find two level 0 and 2 level 1 neighbor elements
+ if ( numL1Neighbors != 2 && numL0Neighbors != 2 )
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly two level 0 neighbors and two level 1 neighbors in middle right element");
+ }
+ //! In the element above the center we should have 2 scvfs in negative y-direction
+ else if ( isMiddleUpperElement(elementCenter) ) //! above the center
+ {
+ if (numScvfsTotal != 5)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 5 scvfs in middle upper element");
+ if (numScvfsPosX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle upper element in positive x-direction");
+ if (numScvfsPosY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle upper element in positive y-direction");
+ if (numScvfsNegX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle upper element in negative x-direction");
+ if (numScvfsNegY != 2)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 2 scvfs in middle upper element in negative y-direction");
+
+ //! we should find two level 0 and 2 level 1 neighbor elements
+ if ( numL1Neighbors != 2 && numL0Neighbors != 2 )
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly two level 0 neighbors and two level 1 neighbors in middle upper element");
+ }
+ //! In the element below the center we should have 2 scvfs in positive y-direction
+ else if ( isMiddleLowerElement(elementCenter) ) //! below the center
+ {
+ if (numScvfsTotal != 5)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 5 scvfs in middle lower element");
+ if (numScvfsPosX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle lower element in positive x-direction");
+ if (numScvfsPosY != 2)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle lower element in positive y-direction");
+ if (numScvfsNegX != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle lower element in negative x-direction");
+ if (numScvfsNegY != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly 1 scvfs in middle lower element in negative y-direction");
+
+ //! we should find two level 0 and 2 level 1 neighbor elements
+ if ( numL1Neighbors != 2 && numL0Neighbors != 2 )
+ DUNE_THROW(Dune::InvalidStateException, "Did not find exactly two level 0 neighbors and two level 2 neighbors in middle lower element");
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Element center position could not be interpreted.");
+ }
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (Dune::Exception e) {
+
+ std::cout << e << std::endl;
+ return 1;
+}