From 06cf7da25c2145cfa4aa7adcc062a10f81c52735 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 10 Apr 2018 17:04:47 +0200
Subject: [PATCH] [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 <numeric>
+#include <functional>
+#include <type_traits>
+#include <dune/common/exceptions.hh>
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
#include <dune/geometry/type.hh>
-#include <dune/geometry/affinegeometry.hh>
+#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
-#include <dumux/common/geometry/intersectspointgeometry.hh>
+#include <dumux/common/geometry/makegeometry.hh>
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 TypeTag>
-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<Scalar, dim>;
-
- static constexpr auto planeDim = dim - 1;
- using PlaneGeometryType = Dune::AffineGeometry< Scalar, planeDim, dim >;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ 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<int mydim, int coordDim, class Scalar=double>
- 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>&& geo)
+ SurfaceData(std::vector<Geo>&& geo)
{
values_.resize(geo.size());
geometries_ = std::forward<decltype(geo)>(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<PlaneGeometryType>;
+ using SurfaceList = std::vector<SurfaceGeometryType>;
/*!
* \brief The constructor
@@ -149,101 +152,87 @@ public:
* \param sol The solution vector
*/
template<class Assembler>
- 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<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "FluxOverPlane.Verbose", false);
+ verbose_ = getParamFromGroup<bool>(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<planeDim,dim>(std::forward<decltype(planes)>(planes));
+ surfaces_[name] = SurfaceData<surfaceDim,dim>(std::forward<decltype(surfaces)>(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<GlobalPosition>{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<GlobalPosition>{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<Dune::FieldVector<Scalar, 2>>& 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<Dune::FieldVector<Scalar, 3>>& 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<class FluxType>
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<bool> 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<std::string, PlaneData<planeDim ,dim> > planes_;
+ std::map<std::string, SurfaceData<surfaceDim ,dim> > 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 <dumux/io/staggeredvtkoutputmodule.hh>
-#include <dumux/freeflow/navierstokes/staggered/fluxoverplane.hh>
+#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
/*!
* \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<GridVariables>(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<TypeTag, GET_PROP_VALUE(TypeTag, PhaseIdx)> 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<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
- // set up two planes over which fluxes are calculated
- FluxOverPlane<TypeTag> flux(*assembler, x);
+ // set up two surfaces over which fluxes are calculated
+ FluxOverSurface<TypeTag> flux(*assembler, x);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
@@ -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<std::vector<int>>("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 <dumux/io/gnuplotinterface.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
-#include <dumux/freeflow/navierstokes/staggered/fluxoverplane.hh>
-
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
--
GitLab