diff --git a/dumux/common/geometry/makegeometry.hh b/dumux/common/geometry/makegeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e5ed691cb46c0c8a186875f3cdaa9a223c93072a
--- /dev/null
+++ b/dumux/common/geometry/makegeometry.hh
@@ -0,0 +1,167 @@
+/*****************************************************************************
+ * 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 3 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Create Dune geometries from user-specified points
+ */
+#ifndef DUMUX_MAKE_GEOMETRY_HH
+#define DUMUX_MAKE_GEOMETRY_HH
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace Dumux {
+
+//! Checks if four points lie within the same plane.
+template
+bool pointsAreCoplanar(const std::vector>& points, Scalar eps = 1e-20)
+{
+ assert(points.size() == 4);
+ // (see "Real-Time Collision Detection" by Christer Ericson)
+ Dune::FieldMatrix M;
+ for(int i = 0; i < 3; ++i )
+ M[i] = {points[0][i], points[1][i], points[2][i], points[3][i]};
+ M[3] = {1.0, 1.0, 1.0, 1.0};
+
+ using std::abs;
+ return abs(M.determinant()) < eps;
+}
+
+/*!
+ * \brief Returns a vector of points following the dune ordering.
+ * Convenience method that creates a temporary object in case no array of orientations is desired.
+ *
+ * \param points The user-specified vector of points (potentially in wrong order).
+ */
+template
+std::vector> getReorderedPoints(const std::vector>& points)
+{
+ std::array tmp;
+ return getReorderedPoints(points, tmp);
+}
+
+/*!
+ * \brief Returns a vector of points following the dune ordering.
+ *
+ * \param points The user-specified vector of points (potentially in wrong order).
+ * \param orientations An array of orientations that can be useful for further processing.
+ */
+template
+std::vector> getReorderedPoints(const std::vector>& points,
+ std::array& orientations)
+{
+ if(points.size() == 4)
+ {
+ auto& p0 = points[0];
+ auto& p1 = points[1];
+ auto& p2 = points[2];
+ auto& p3 = points[3];
+
+ // check if the points define a proper quadrilateral
+ const auto normal = crossProduct((p1 - p0), (p2 - p0));
+
+ orientations = { getOrientation(p0, p3, p2, normal),
+ getOrientation(p0, p3, p1, normal),
+ getOrientation(p2, p1, p0, normal),
+ getOrientation(p2, p1, p3, normal) };
+
+
+ // check if the points follow the dune ordering (see http://www.dcs.gla.ac.uk/~pat/52233/slides/Geometry1x1.pdf)
+ const bool diagonalsIntersect = (orientations[0] != orientations[1]) && (orientations[2] != orientations[3]);
+
+ // the points conform with the dune ordering
+ if(diagonalsIntersect)
+ return points;
+
+ // the points do not conform with the dune ordering, re-order
+ using GlobalPosition = Dune::FieldVector;
+ if(!diagonalsIntersect && orientations[0] == 1)
+ return std::vector{p1, p0, p2, p3};
+ else if(!diagonalsIntersect && orientations[0] == -1)
+ return std::vector{p3, p1, p0, p2};
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Could not reorder points");
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented, "Reorder for " << points.size() << " points.");
+}
+
+/*!
+ * \brief Creates a dune quadrilateral geometry given 4 corner points.
+ *
+ * \tparam Scalar The Scalar type.
+ * \tparam enableSanityCheck Turn on/off sanity check and reordering of points
+ * \param points The user-specified vector of points (potentially in wrong order).
+ */
+template
+auto makeDuneQuadrilaterial(const std::vector>& points)
+{
+ assert(points.size() == 4 && "A quadrilateral needs 4 corner points!");
+
+ using GlobalPosition = Dune::FieldVector;
+ static constexpr auto coordDim = GlobalPosition::dimension;
+ static constexpr auto dim = coordDim-1;
+ using GeometryType = Dune::MultiLinearGeometry;
+
+ // if no sanity check if desired, use the given points directly to construct the geometry
+ if(!enableSanityCheck)
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ return GeometryType(Dune::GeometryTypes::quadrilateral, points);
+#else
+ {
+ static Dune::GeometryType gt(Dune::GeometryType::cube, dim);
+ return GeometryType(gt, points);
+ }
+#endif
+
+ // otherwise, perform a number of checks and corrections
+ if(!pointsAreCoplanar(points))
+ DUNE_THROW(Dune::InvalidStateException, "Points do not lie within a plane");
+
+ // make sure points conform with dune ordering
+ std::array orientations;
+ auto corners = getReorderedPoints(points, orientations);
+
+ if(std::any_of(orientations.begin(), orientations.end(), [](auto i){ return i == 0; }))
+ DUNE_THROW(Dune::InvalidStateException, "More than two points lie on the same line.");
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ const auto quadrilateral = GeometryType(Dune::GeometryTypes::quadrilateral, corners);
+#else
+ static Dune::GeometryType gt(Dune::GeometryType::cube, dim);
+ const auto quadrilateral = GeometryType(gt, points);
+#endif
+
+ const auto eps = 1e-20;
+ if(quadrilateral.volume() < eps)
+ DUNE_THROW(Dune::InvalidStateException, "Something went wrong, geometry has area of zero");
+
+ return quadrilateral;
+}
+
+
+
+} // end namespace Dumux
+
+#endif // DUMUX_MAKE_GEOMETRY_HH
diff --git a/dumux/freeflow/navierstokes/staggered/CMakeLists.txt b/dumux/freeflow/navierstokes/staggered/CMakeLists.txt
index 4e0e1f7ca324ec44abdd3a925848802768348a6e..d12ae0320aa19d72907b30cb8945ef2968c7f4d4 100644
--- a/dumux/freeflow/navierstokes/staggered/CMakeLists.txt
+++ b/dumux/freeflow/navierstokes/staggered/CMakeLists.txt
@@ -1,7 +1,7 @@
#install headers
install(FILES
-fluxoverplane.hh
+fluxoversurface.hh
fluxvariables.hh
localresidual.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/navierstokes/staggered)
diff --git a/dumux/freeflow/navierstokes/staggered/fluxoverplane.hh b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
similarity index 66%
rename from dumux/freeflow/navierstokes/staggered/fluxoverplane.hh
rename to dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
index cb1420e263a59d045d949bdfe93fc7aebe3ef82a..0f2c4f1448094d6748625b8d8d9bfb1fca4ad1da 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxoverplane.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
@@ -19,31 +19,35 @@
/*!
* \file
* \ingroup NavierStokesModel
- * \copydoc Dumux::FluxOverPlane
+ * \copydoc Dumux::FluxOverSurface
*/
-#ifndef DUMUX_FLUX_OVER_PLANE_STAGGERED_HH
-#define DUMUX_FLUX_OVER_PLANE_STAGGERED_HH
+#ifndef DUMUX_FLUX_OVER_SURFACE_STAGGERED_HH
+#define DUMUX_FLUX_OVER_SURFACE_STAGGERED_HH
#include
+#include
+#include
+#include
#include
#include
+#include
#include
-#include
+#include
#include
#include
#include
-#include
+#include
namespace Dumux {
/*!
* \ingroup NavierStokesModel
- * \brief Class used to calculate fluxes over planes. This only works for the staggered grid discretization.
+ * \brief Class used to calculate fluxes over surfaces. This only works for the staggered grid discretization.
*/
template
-class FluxOverPlane
+class FluxOverSurface
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -64,54 +68,53 @@ class FluxOverPlane
dimWorld = GridView::dimensionworld
};
- using GlobalPosition = Dune::FieldVector;
-
- static constexpr auto planeDim = dim - 1;
- using PlaneGeometryType = Dune::AffineGeometry< Scalar, planeDim, dim >;
+ using GlobalPosition = Dune::FieldVector;
+ static constexpr auto surfaceDim = dimWorld - 1;
+ using SurfaceGeometryType = Dune::MultiLinearGeometry< Scalar, surfaceDim, dimWorld >;
/*!
- * \brief Auxiliary class that holds the plane-specific data.
+ * \brief Auxiliary class that holds the surface-specific data.
*/
template
- class PlaneData
+ class SurfaceData
{
- using Geo = Dune::AffineGeometry< Scalar, mydim, coordDim >;
+ using Geo = SurfaceGeometryType;
public:
- PlaneData() {}
+ SurfaceData() {}
- PlaneData(Geo&& geo)
+ SurfaceData(Geo&& geo)
{
values_.resize(1);
geometries_.push_back(std::move(geo));
}
- PlaneData(std::vector&& geo)
+ SurfaceData(std::vector&& geo)
{
values_.resize(geo.size());
geometries_ = std::forward(geo);
}
- void addValue(int planeIdx, const CellCenterPrimaryVariables& value)
+ void addValue(int surfaceIdx, const CellCenterPrimaryVariables& value)
{
- values_[planeIdx] += value;
+ values_[surfaceIdx] += value;
}
- void addSubPlane(Geo&& geo)
+ void addSubSurface(Geo&& geo)
{
values_.emplace_back(0.0);
geometries_.push_back(std::move(geo));
}
- auto& subPlanes() const
+ auto& subSurfaces() const
{
return geometries_;
}
- CellCenterPrimaryVariables& value(int planeIdx) const
+ CellCenterPrimaryVariables& value(int surfaceIdx) const
{
- return values_[planeIdx];
+ return values_[surfaceIdx];
}
auto& values() const
@@ -119,11 +122,11 @@ class FluxOverPlane
return values_;
}
- void printPlaneBoundaries(int planeIdx) const
+ void printSurfaceBoundaries(int surfaceIdx) const
{
- const auto& geometry = geometries_[planeIdx];
+ const auto& geometry = geometries_[surfaceIdx];
for(int i = 0; i < geometry.corners(); ++i)
- std::cout << geometry.corner(i) << " ";
+ std::cout << "(" << geometry.corner(i) << ") ";
}
void resetValues()
@@ -140,7 +143,7 @@ class FluxOverPlane
public:
- using PlaneList = std::vector;
+ using SurfaceList = std::vector;
/*!
* \brief The constructor
@@ -149,101 +152,87 @@ public:
* \param sol The solution vector
*/
template
- FluxOverPlane(const Assembler& assembler,
- const SolutionVector& sol)
- : problem_(assembler.problem())
- , gridVariables_(assembler.gridVariables())
- , sol_(sol)
+ FluxOverSurface(const Assembler& assembler,
+ const SolutionVector& sol)
+ : problem_(assembler.problem()),
+ gridVariables_(assembler.gridVariables()),
+ sol_(sol)
{
- verbose_ = getParamFromGroup(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "FluxOverPlane.Verbose", false);
+ verbose_ = getParamFromGroup(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "FluxOverSurface.Verbose", false);
}
/*!
- * \brief Add a collection of sub planes under a given name
+ * \brief Add a collection of sub surfaces under a given name
*
- * \param name The name of the plane
- * \param planes The list of sub planes
+ * \param name The name of the surface
+ * \param surfaces The list of sub surfaces
*/
- void addPlane(const std::string& name, PlaneList&& planes )
+ void addSurface(const std::string& name, SurfaceList&& surfaces)
{
- planes_[name] = PlaneData(std::forward(planes));
+ surfaces_[name] = SurfaceData(std::forward(surfaces));
}
/*!
- * \brief Add a plane under a given name, specifing the plane's corner points.
- * This is a specialization for 2D, therefore the plane is actually a line.
+ * \brief Add a surface under a given name, specifying the surface's corner points.
+ * This is a specialization for 2D, therefore the surface is actually a line.
*
- * \param name The name of the plane
+ * \param name The name of the surface
* \param p0 The first corner
* \param p1 The second corner
*/
- void addPlane(const std::string& name, const GlobalPosition& p0, const GlobalPosition& p1)
+ void addSurface(const std::string& name, const GlobalPosition& p0, const GlobalPosition& p1)
{
- planes_[name].addSubPlane(makePlane(p0, p1));
+ surfaces_[name].addSubSurface(makeSurface(std::vector{p0, p1}));
}
/*!
- * \brief Add a plane under a given name, specifing the plane's corner points.
+ * \brief Add a surface under a given name, specifying the surface's corner points.
* This is a specialization for 3D.
*
- * \param name The name of the plane
+ * \param name The name of the surface
* \param p0 The first corner
* \param p1 The second corner
* \param p2 The third corner
* \param p3 The fourth corner
*/
- void addPlane(const std::string& name,
- const GlobalPosition& p0,
- const GlobalPosition& p1,
- const GlobalPosition& p2,
- const GlobalPosition& p3)
+ void addSurface(const std::string& name,
+ const GlobalPosition& p0,
+ const GlobalPosition& p1,
+ const GlobalPosition& p2,
+ const GlobalPosition& p3)
{
- planes_[name].addSubPlane(makePlane(p0, p1, p2, p3));
+ surfaces_[name].addSubSurface(makeSurface(std::vector{p0, p1, p2, p3}));
}
/*!
- * \brief Creates a geometrical plane object.
- * This is a specialization for 2D, therefore the plane is actually a line.
+ * \brief Creates a geometrical surface object for (2D).
*
- * \param p0 The first corner
- * \param p1 The second corner
+ * \param corners The vector storing the surface's corners
*/
- static PlaneGeometryType makePlane(const GlobalPosition& p0, const GlobalPosition& p1)
+ static SurfaceGeometryType makeSurface(const std::vector>& corners)
{
- const std::vector< Dune::FieldVector< Scalar, dim > > corners = {p0, p1};
#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
- return PlaneGeometryType(Dune::GeometryTypes::line, corners);
+ return SurfaceGeometryType(Dune::GeometryTypes::line, corners);
#else
- static Dune::GeometryType gt(Dune::GeometryType::simplex, dim-1);
- return PlaneGeometryType(gt, corners);
+ {
+ static Dune::GeometryType gt2d(Dune::GeometryType::simplex, dim-1);
+ return SurfaceGeometryType(gt2d, corners);
+ }
#endif
}
/*!
- * \brief Creates a geometrical plane object.
- * This is a specialization for 3D.
+ * \brief Creates a geometrical surface object for (3D).
*
- * \param p0 The first corner
- * \param p1 The second corner
- * \param p2 The third corner
- * \param p3 The fourth corner
+ * \param corners The vector storing the surface's corners
*/
- static PlaneGeometryType makePlane(const GlobalPosition& p0,
- const GlobalPosition& p1,
- const GlobalPosition& p2,
- const GlobalPosition& p3)
+ static SurfaceGeometryType makeSurface(const std::vector>& corners)
{
- const std::vector< Dune::FieldVector< Scalar, dim > > corners = {p0, p1, p2, p3};
-#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
- return PlaneGeometryType(Dune::GeometryTypes::quadrilateral, corners);
-#else
- static Dune::GeometryType gt(Dune::GeometryType::cube, dim-1);
- return PlaneGeometryType(gt, corners);
-#endif
+ return makeDuneQuadrilaterial(corners);
}
/*!
- * \brief Calculate the mass or mole fluxes over all planes
+ * \brief Calculate the mass or mole fluxes over all surfaces
*/
void calculateMassOrMoleFluxes()
{
@@ -262,8 +251,7 @@ public:
}
/*!
- * \brief Calculate the volume fluxes over all planes.
- * This method simply averages the densities between two adjacent cells.
+ * \brief Calculate the volume fluxes over all surfaces.
*/
void calculateVolumeFluxes()
{
@@ -272,7 +260,7 @@ public:
}
/*!
- * \brief Calculate the fluxes over all planes for a given flux type.
+ * \brief Calculate the fluxes over all surfaces for a given flux type.
*
* \param fluxType The flux type. This can be a lambda of the following form:
* [](const auto& problem,
@@ -287,9 +275,9 @@ public:
template
void calculateFluxes(const FluxType& fluxType)
{
- // make sure to reset all the values of the planes, in case this method has been called already before
- for(auto&& plane : planes_)
- plane.second.resetValues();
+ // make sure to reset all the values of the surfaces, in case this method has been called already before
+ for(auto&& surface : surfaces_)
+ surface.second.resetValues();
// make sure not to iterate over the same dofs twice
std::vector dofVisited(problem_.fvGridGeometry().numFaceDofs(), false);
@@ -317,26 +305,26 @@ public:
dofVisited[dofIdx] = true;
- // iterate through all planes and check if the flux at the given position
- // should be accounted for in the respective plane
- for(auto&& plane : planes_)
+ // iterate through all surfaces and check if the flux at the given position
+ // should be accounted for in the respective surface
+ for(auto&& surface : surfaces_)
{
- const auto& subPlanes = plane.second.subPlanes();
+ const auto& subSurfaces = surface.second.subSurfaces();
- for(int planeIdx = 0; planeIdx < subPlanes.size(); ++planeIdx)
+ for(int surfaceIdx = 0; surfaceIdx < subSurfaces.size(); ++surfaceIdx)
{
- if(intersectsPointGeometry(scvf.center(), subPlanes[planeIdx]))
+ if(intersectsPointGeometry(scvf.center(), subSurfaces[surfaceIdx]))
{
const auto result = fluxType(problem_, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
- plane.second.addValue(planeIdx, result);
+ surface.second.addValue(surfaceIdx, result);
if(verbose_)
{
- std::cout << "Flux at face " << scvf.center() << " (" << plane.first << "): " << result;
- std::cout << " (directionIndex: " << scvf.directionIndex() << "; plane boundaries: ";
- plane.second.printPlaneBoundaries(planeIdx);
- std::cout << ", planeIdx " << planeIdx << ")" << std::endl;
+ std::cout << "Flux at face " << scvf.center() << " (" << surface.first << "): " << result;
+ std::cout << " (directionIndex: " << scvf.directionIndex() << "; surface boundaries: ";
+ surface.second.printSurfaceBoundaries(surfaceIdx);
+ std::cout << ", surfaceIdx " << surfaceIdx << ")" << std::endl;
}
}
}
@@ -346,30 +334,30 @@ public:
}
/*!
- * \brief Return the fluxes of the individual sub planes of a given name.
+ * \brief Return the fluxes of the individual sub surface of a given name.
*
- * \param name The name of the plane
+ * \param name The name of the surface
*/
auto& values(const std::string& name) const
{
- return planes_.at(name).values();
+ return surfaces_.at(name).values();
}
/*!
- * \brief Return the cumulative net fluxes of a plane of a given name.
+ * \brief Return the cumulative net fluxes of a surface of a given name.
*
- * \param name The name of the plane
+ * \param name The name of the surface
*/
auto netFlux(const std::string& name) const
{
- const auto& planeResults = values(name);
- return std::accumulate(planeResults.begin(), planeResults.end(), CellCenterPrimaryVariables(0.0));
+ const auto& surfaceResults = values(name);
+ return std::accumulate(surfaceResults.begin(), surfaceResults.end(), CellCenterPrimaryVariables(0.0));
}
private:
/*!
- * \brief Calculate the volume fluxes over all planes for compositional models.
+ * \brief Calculate the volume fluxes over all surfaces for compositional models.
* This method simply averages the densities between two adjacent cells.
*/
void calculateVolumeFluxesImpl_(std::true_type)
@@ -421,7 +409,7 @@ private:
}
/*!
- * \brief Calculate the volume fluxes over all planes for non-compositional models.
+ * \brief Calculate the volume fluxes over all surfaces for non-compositional models.
* This method simply averages the densities between two adjacent cells.
*/
void calculateVolumeFluxesImpl_(std::false_type)
@@ -451,7 +439,7 @@ private:
calculateFluxes(fluxType);
}
- std::map > planes_;
+ std::map > surfaces_;
const Problem& problem_;
const GridVariables& gridVariables_;
const SolutionVector& sol_;
diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index c127993f4d07a3e12686521d6259c384c7dcb76e..b20b36ca512d99a7494b45cdde85454227704f9e 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -1,12 +1,14 @@
dune_add_test(SOURCES test_1d3d_intersection.cc)
dune_add_test(SOURCES test_2d3d_intersection.cc)
dune_add_test(SOURCES test_graham_convex_hull.cc)
+dune_add_test(SOURCES test_makegeometry.cc)
#install sources
install(FILES
test_1d3d_intersection.cc
test_2d3d_intersection.cc
test_graham_convex_hull.cc
+test_makegeometry.cc
writetriangulation.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/common/geometry)
diff --git a/test/common/geometry/test_makegeometry.cc b/test/common/geometry/test_makegeometry.cc
new file mode 100644
index 0000000000000000000000000000000000000000..78821515ad733c23a2f1869aba998c1f20a20493
--- /dev/null
+++ b/test/common/geometry/test_makegeometry.cc
@@ -0,0 +1,213 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Test for the makegeometry method
+ */
+ #include
+
+ #include
+ #include
+ #include
+
+#include
+#include
+#include
+
+//! test if the ordering of the input points has an effect on the resulting geometry
+template
+void permutatePointsAndTest(const std::vector& cornerPoints,
+ const std::vector& pointsWithinGeometry,
+ const std::vector& pointsOutsideGeometry,
+ const bool verbose = false)
+{
+ std::array s = {0,1,2,3};
+
+ const auto area = Dumux::makeDuneQuadrilaterial(cornerPoints).volume();
+
+ do {
+ std::vector tmp = {cornerPoints[s[0]], cornerPoints[s[1]], cornerPoints[s[2]], cornerPoints[s[3]]};
+
+ if(verbose)
+ {
+ std::cout << "input corner points: " << std::endl;
+ for(auto&& i : tmp)
+ std::cout << "(" << i << ") ";
+ std::cout << std::endl;
+ }
+
+ const auto quad = Dumux::makeDuneQuadrilaterial(tmp);
+
+ if(verbose)
+ {
+ std::cout << "resulting corners \n";
+ for(int i = 0; i < quad.corners(); ++i)
+ std::cout << "(" << quad.corner(i) << ") ";
+ std::cout << std::endl;
+
+ std::cout << "actual area: " << area << ", area after permuation: " << quad.volume() << std::endl;
+ }
+
+ if(!Dune::FloatCmp::eq(quad.volume(), area))
+ DUNE_THROW(Dune::InvalidStateException, "Area of quadrilateral after permuation of input points is wrong");
+
+ if(verbose)
+ std::cout << "checking point intersection: \n";
+
+ for(const auto& p: pointsWithinGeometry)
+ {
+ if(Dumux::intersectsPointGeometry(p, quad))
+ {
+ if(verbose)
+ std::cout << "point " << p << " lies within the quadrilateral" << std::endl;
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Check for point inside geometry failed. Point " << p << " does not lie within the geometry!");
+ }
+
+ for(const auto& p : pointsOutsideGeometry)
+ {
+ if(!Dumux::intersectsPointGeometry(p, quad))
+ {
+ if(verbose)
+ std::cout << "point " << p << " lies outside of the quadrilateral" << std::endl;
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Check for point outside geometry failed. Point " << p << " does lie within the geometry!");
+ }
+
+ } while(std::next_permutation(s.begin(), s.end()));
+}
+
+template
+void checkAxisAlignedGeometry(std::vector& cornerPoints,
+ std::vector& pointsWithinGeometry,
+ std::vector& pointsOutsideGeometry,
+ const int normalDirection)
+{
+ static const char dim[]= "xyz";
+ std::cout << "testing for quadrilateral with normal in " << dim[normalDirection] << " direction" << std::endl;
+
+ // check if points are coplanar
+ if(!Dumux::pointsAreCoplanar(cornerPoints))
+ DUNE_THROW(Dune::InvalidStateException, "False positive, points are actually coplanar!");
+
+ cornerPoints[0][normalDirection] += 1e-9;
+ if(Dumux::pointsAreCoplanar(cornerPoints))
+ DUNE_THROW(Dune::InvalidStateException, "Points are not coplanar!");
+
+ cornerPoints[0][normalDirection] -= 1e-9;
+
+ permutatePointsAndTest(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry);
+
+ const auto origCornerPoints = cornerPoints;
+
+ std::vector pointsWithinGeometryForRandomTest = { {0.6, 0.6, 0.6},
+ {0.4, 0.4, 0.4} };
+ for(auto& p : pointsWithinGeometryForRandomTest)
+ p[normalDirection] = 0.0;
+
+ auto pointsOutsideGeometryForRandomTest = pointsOutsideGeometry;
+ pointsOutsideGeometryForRandomTest[0] = {1.5, 1.5, 1.5};
+ pointsOutsideGeometryForRandomTest[0][normalDirection] = 0.0;
+
+ for(int i = 0; i< 10; i++)
+ {
+ // uniform random number generator
+ std::random_device rd;
+ std::mt19937 generator(rd());
+ std::uniform_real_distribution<> uniformdist(-0.3, 0.3);
+
+ for(auto&& p : cornerPoints)
+ {
+ for(int x = 0; x < p.size(); ++x)
+ {
+ if(x != normalDirection)
+ p[x] += uniformdist(generator);
+ }
+ }
+
+ permutatePointsAndTest(cornerPoints, pointsWithinGeometryForRandomTest, pointsOutsideGeometryForRandomTest);
+
+ cornerPoints = origCornerPoints;
+ }
+
+
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+ using GlobalPosition = Dune::FieldVector;
+
+ GlobalPosition p0 = {0,0,0};
+ GlobalPosition p1 = {1,0,0};
+ GlobalPosition p2 = {0,1,0};
+ GlobalPosition p3 = {1,1,0};
+
+ std::vector cornerPoints = {p0, p1, p2, p3};
+
+ std::vector pointsWithinGeometry = { GlobalPosition{0.5, 0.5, 0.0},
+ GlobalPosition{cornerPoints[0][0] + 1e-3, cornerPoints[0][1] + 1e-3, 0.0},
+ GlobalPosition{cornerPoints[3][0] - 1e-3, cornerPoints[3][1] - 1e-3, 0.0} };
+
+ std::vector pointsOutsideGeometry = { GlobalPosition{cornerPoints[0][0] - 1e-3, cornerPoints[0][1] - 1e-3, 0.0},
+ GlobalPosition{0.5, 0.5, 1e-3} };
+
+ // do the checks for a quadrilateral parallel to the x and y axis
+ checkAxisAlignedGeometry(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry, 2);
+
+ // rotate the quadrilateral to make it parallel to the other axes and test again
+ for(int i = 1; i >=0; --i)
+ {
+ for(auto& p : cornerPoints)
+ std::rotate(p.begin(), p.begin() + 1, p.end());
+ for(auto& p : pointsWithinGeometry)
+ std::rotate(p.begin(), p.begin() + 1, p.end());
+ for(auto& p : pointsOutsideGeometry)
+ std::rotate(p.begin(), p.begin() + 1, p.end());
+
+ checkAxisAlignedGeometry(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry, i);
+ }
+
+ std::cout << "testing for non axis-aligned quadrilateral" << std::endl;
+
+ cornerPoints[0][0] += 0.5;
+ cornerPoints[1][0] -= 0.5;
+ cornerPoints[2][0] += 0.5;
+ cornerPoints[3][0] -= 0.5;
+
+ GlobalPosition pointToCheck5 = {0.0, 0.5, 0.5};
+
+ pointsWithinGeometry = {pointToCheck5};
+
+ pointsOutsideGeometry[1] = pointsOutsideGeometry[0];
+
+ permutatePointsAndTest(cornerPoints, pointsWithinGeometry, pointsOutsideGeometry);
+
+ return 0;
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (const Dune::Exception& e) {
+ std::cout << e << std::endl;
+ return 1;
+}
diff --git a/test/freeflow/navierstokes/test_channel.cc b/test/freeflow/navierstokes/test_channel.cc
index 1ac6b724a67d80c8e9ebf14c080486b9efe0308a..505abff099579731be8faf0be1c58f59b609a3f2 100644
--- a/test/freeflow/navierstokes/test_channel.cc
+++ b/test/freeflow/navierstokes/test_channel.cc
@@ -50,7 +50,7 @@
#include
-#include
+#include
/*!
* \brief Provides an interface for customizing error messages associated with
@@ -153,7 +153,7 @@ int main(int argc, char** argv) try
auto gridVariables = std::make_shared(problem, fvGridGeometry);
gridVariables->init(x, xOld);
- // intialize the vtk output module
+ // initialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
StaggeredVtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
@@ -171,8 +171,8 @@ int main(int argc, char** argv) try
using NewtonSolver = Dumux::NewtonSolver;
NewtonSolver nonLinearSolver(assembler, linearSolver);
- // set up two planes over which fluxes are calculated
- FluxOverPlane flux(*assembler, x);
+ // set up two surfaces over which fluxes are calculated
+ FluxOverSurface flux(*assembler, x);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using GlobalPosition = Dune::FieldVector;
@@ -181,23 +181,23 @@ int main(int argc, char** argv) try
const Scalar yMin = fvGridGeometry->bBoxMin()[1];
const Scalar yMax = fvGridGeometry->bBoxMax()[1];
- // The first plane shall be placed at the middle of the channel.
+ // The first surface shall be placed at the middle of the channel.
// If we have an odd number of cells in x-direction, there would not be any cell faces
- // at the postion of the plane (which is required for the flux calculation).
+ // at the position of the surface (which is required for the flux calculation).
// In this case, we add half a cell-width to the x-position in order to make sure that
- // the cell faces lie on the plane. This assumes a regular cartesian grid.
+ // the cell faces lie on the surface. This assumes a regular cartesian grid.
const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin);
const int numCellsX = getParam>("Grid.Cells")[0];
const Scalar offsetX = (numCellsX % 2 == 0) ? 0.0 : 0.5*((xMax - xMin) / numCellsX);
const auto p0middle = GlobalPosition{planePosMiddleX + offsetX, yMin};
const auto p1middle = GlobalPosition{planePosMiddleX + offsetX, yMax};
- flux.addPlane("middle", p0middle, p1middle);
+ flux.addSurface("middle", p0middle, p1middle);
- // The second plane is placed at the outlet of the channel.
+ // The second surface is placed at the outlet of the channel.
const auto p0outlet = GlobalPosition{xMax, yMin};
const auto p1outlet = GlobalPosition{xMax, yMax};
- flux.addPlane("outlet", p0outlet, p1outlet);
+ flux.addSurface("outlet", p0outlet, p1outlet);
// time loop
timeLoop->start(); do
diff --git a/test/freeflow/rans/test_pipe_laufer.cc b/test/freeflow/rans/test_pipe_laufer.cc
index f01a3c253ac3a549ed21c235224ec78c8c21effd..cffaff6fbce264b84a3b3eb2c301d1b0a721a16b 100644
--- a/test/freeflow/rans/test_pipe_laufer.cc
+++ b/test/freeflow/rans/test_pipe_laufer.cc
@@ -54,8 +54,6 @@
#include
#include
-#include
-
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.