From 23974f94bab1ef1be8ee71e9eea7d280aec43df9 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 14 Apr 2020 16:24:15 +0200
Subject: [PATCH] [fvproblem] Move applyInitialSolution to a separate header
---
dumux/assembly/CMakeLists.txt | 1 +
dumux/assembly/initialsolution.hh | 111 ++++++++++++++++++++++++++++++
dumux/common/fvproblem.hh | 59 +---------------
3 files changed, 115 insertions(+), 56 deletions(-)
create mode 100644 dumux/assembly/initialsolution.hh
diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt
index 65a7539064..2543a6631f 100644
--- a/dumux/assembly/CMakeLists.txt
+++ b/dumux/assembly/CMakeLists.txt
@@ -8,6 +8,7 @@ entitycolor.hh
fvassembler.hh
fvlocalassemblerbase.hh
fvlocalresidual.hh
+initialsolution.hh
jacobianpattern.hh
numericepsilon.hh
partialreassembler.hh
diff --git a/dumux/assembly/initialsolution.hh b/dumux/assembly/initialsolution.hh
new file mode 100644
index 0000000000..dc8bf7ebf3
--- /dev/null
+++ b/dumux/assembly/initialsolution.hh
@@ -0,0 +1,111 @@
+// -*- 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 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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Function to create initial solution vectors
+ */
+#ifndef DUMUX_ASSEMBLY_INITIAL_SOLUTION_HH
+#define DUMUX_ASSEMBLY_INITIAL_SOLUTION_HH
+
+#include <vector>
+#include <type_traits>
+
+#include <dumux/discretization/method.hh>
+
+namespace Dumux {
+
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Set a solution vector to the initial solution provided by the problem
+ */
+template<class SolutionVector, class Problem>
+void assembleInitialSolution(SolutionVector& sol, const Problem& problem)
+{
+ const auto& gg = problem.gridGeometry();
+ using GridGeometry = std::decay_t<decltype(gg)>;
+
+ // box method
+ if constexpr (GridGeometry::discMethod == DiscretizationMethod::box)
+ {
+ constexpr int dim = GridGeometry::GridView::dimension;
+ const auto numDofs = gg.vertexMapper().size();
+ const auto numVert = gg.gridView().size(dim);
+ sol.resize(numDofs);
+
+ // if there are more dofs than vertices (enriched nodal dofs), we have to
+ // call initial for all dofs at the nodes, coming from all neighboring elements.
+ if (numDofs != numVert)
+ {
+ std::vector<bool> dofVisited(numDofs, false);
+ for (const auto& element : elements(gg.gridView()))
+ {
+ for (int i = 0; i < element.subEntities(dim); ++i)
+ {
+ const auto dofIdxGlobal = gg.vertexMapper().subIndex(element, i, dim);
+ // forward to implementation if value at dof is not set yet
+ if (!dofVisited[dofIdxGlobal])
+ {
+ sol[dofIdxGlobal] = problem.initial(element.template subEntity<dim>(i));
+ dofVisited[dofIdxGlobal] = true;
+ }
+ }
+ }
+ }
+
+ // otherwise we directly loop over the vertices
+ else
+ {
+ for (const auto& vertex : vertices(gg.gridView()))
+ sol[gg.vertexMapper().index(vertex)] = problem.initial(vertex);
+ }
+ }
+
+ // staggered methods
+ else if constexpr (GridGeometry::discMethod == DiscretizationMethod::staggered)
+ {
+ problem.applyInitialSolution(sol);
+ }
+
+ // default: cell-centered methods
+ else
+ {
+ sol.resize(gg.numDofs());
+ for (const auto& element : elements(gg.gridView()))
+ sol[gg.elementMapper().index(element)] = problem.initial(element);
+ }
+}
+
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Create a solution vector filled with the initial solution provided by the problem
+ */
+template<class SolutionVector, class Problem>
+SolutionVector makeInitialSolution(const Problem& problem)
+{
+ SolutionVector sol;
+ assembleInitialSolution(sol, problem);
+ return sol;
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index 47e6f4313e..49dea16866 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -34,6 +34,8 @@
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
+#include <dumux/assembly/initialsolution.hh>
+
namespace Dumux {
/*!
@@ -510,8 +512,7 @@ public:
*/
void applyInitialSolution(SolutionVector& sol) const
{
- // set the initial values by forwarding to a specialized method
- applyInitialSolutionImpl_(sol, std::integral_constant<bool, isBox>());
+ assembleInitialSolution(sol, asImp_());
}
/*!
@@ -596,60 +597,6 @@ protected:
{ return *static_cast<const Implementation *>(this); }
private:
- /*!
- * \brief Applies the initial solution for the box method
- */
- void applyInitialSolutionImpl_(SolutionVector& sol, /*isBox=*/std::true_type) const
- {
- const auto numDofs = gridGeometry_->vertexMapper().size();
- const auto numVert = gridGeometry_->gridView().size(dim);
- sol.resize(numDofs);
-
- // if there are more dofs than vertices (enriched nodal dofs), we have to
- // call initial for all dofs at the nodes, coming from all neighboring elements.
- if (numDofs != numVert)
- {
- std::vector<bool> dofVisited(numDofs, false);
- for (const auto& element : elements(gridGeometry_->gridView()))
- {
- for (int i = 0; i < element.subEntities(dim); ++i)
- {
- const auto dofIdxGlobal = gridGeometry_->vertexMapper().subIndex(element, i, dim);
-
- // forward to implementation if value at dof is not set yet
- if (!dofVisited[dofIdxGlobal])
- {
- sol[dofIdxGlobal] = asImp_().initial(element.template subEntity<dim>(i));
- dofVisited[dofIdxGlobal] = true;
- }
- }
- }
- }
-
- // otherwise we directly loop over the vertices
- else
- {
- for (const auto& vertex : vertices(gridGeometry_->gridView()))
- {
- const auto dofIdxGlobal = gridGeometry_->vertexMapper().index(vertex);
- sol[dofIdxGlobal] = asImp_().initial(vertex);
- }
- }
- }
-
- /*!
- * \brief Applies the initial solution for cell-centered methods
- */
- void applyInitialSolutionImpl_(SolutionVector& sol, /*isBox=*/std::false_type) const
- {
- sol.resize(gridGeometry_->numDofs());
- for (const auto& element : elements(gridGeometry_->gridView()))
- {
- const auto dofIdxGlobal = gridGeometry_->elementMapper().index(element);
- sol[dofIdxGlobal] = asImp_().initial(element);
- }
- }
-
//! The finite volume grid geometry
std::shared_ptr<const GridGeometry> gridGeometry_;
--
GitLab