From 5f2893d1320b6e98b56b16173326b4ed5d5732d2 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 16 Apr 2018 09:45:37 +0200 Subject: [PATCH 1/2] [geometry] Add method to robustly create dune quadrilateral from 4 given points --- dumux/common/geometry/makegeometry.hh | 167 +++++++++++++++++ test/common/geometry/CMakeLists.txt | 2 + test/common/geometry/test_makegeometry.cc | 213 ++++++++++++++++++++++ 3 files changed, 382 insertions(+) create mode 100644 dumux/common/geometry/makegeometry.hh create mode 100644 test/common/geometry/test_makegeometry.cc diff --git a/dumux/common/geometry/makegeometry.hh b/dumux/common/geometry/makegeometry.hh new file mode 100644 index 0000000000..e5ed691cb4 --- /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/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt index c127993f4d..b20b36ca51 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 0000000000..78821515ad --- /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; +} -- GitLab From 06cf7da25c2145cfa4aa7adcc062a10f81c52735 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Tue, 10 Apr 2018 17:04:47 +0200 Subject: [PATCH 2/2] [staggered][navierstokes] rename FluxOverPlane to FluxOverSurface and improve --- .../navierstokes/staggered/CMakeLists.txt | 2 +- .../{fluxoverplane.hh => fluxoversurface.hh} | 198 ++++++++---------- test/freeflow/navierstokes/test_channel.cc | 20 +- test/freeflow/rans/test_pipe_laufer.cc | 2 - 4 files changed, 104 insertions(+), 118 deletions(-) rename dumux/freeflow/navierstokes/staggered/{fluxoverplane.hh => fluxoversurface.hh} (66%) diff --git a/dumux/freeflow/navierstokes/staggered/CMakeLists.txt b/dumux/freeflow/navierstokes/staggered/CMakeLists.txt index 4e0e1f7ca3..d12ae0320a 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 cb1420e263..0f2c4f1448 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/freeflow/navierstokes/test_channel.cc b/test/freeflow/navierstokes/test_channel.cc index 1ac6b724a6..505abff099 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 f01a3c253a..cffaff6fbc 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. -- GitLab