diff --git a/CMakeLists.txt b/CMakeLists.txt
index 201de083f0971f7e00d7f7ec278384ed5bce8666..57e78b71729483cb0be62688a438df0316b469fb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,6 +25,7 @@ dune_enable_all_packages()
add_subdirectory(doc)
add_subdirectory(cmake/modules)
add_subdirectory(appl)
+add_subdirectory(dumux)
# finalize the dune project, e.g. generating config.h etc.
finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..88dedfd800890b8f320070dcb07b574541976767
--- /dev/null
+++ b/dumux/CMakeLists.txt
@@ -0,0 +1,4 @@
+#install headers
+install(FILES
+poregeometries.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/)
diff --git a/dumux/poregeometries.hh b/dumux/poregeometries.hh
new file mode 100644
index 0000000000000000000000000000000000000000..114fdf67a7162ea15c1c5ddedd35ff7780efc013
--- /dev/null
+++ b/dumux/poregeometries.hh
@@ -0,0 +1,276 @@
+// -*- 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 Freeflow Porous Media Flow
+ * \ingroup Pore Geometry
+ * \brief The problem file for the evaluation of a pore geometry adjacent to a freeflow with RANS models
+ */
+
+#ifndef DUMUX_PORE_GEOMETRIES_HH
+#define DUMUX_PORE_GEOMETRIES_HH
+
+namespace Dumux::PoreGeometries {
+
+template <class PoreBlockList, class Scalar, class DimensionPair, class BoundingBox>
+void buildUniformPoreGeometry(PoreBlockList& poreBlockList,
+ Scalar porosity,
+ DimensionPair numVolumes,
+ BoundingBox porousMediaBox)
+{
+ using PoreBlock = typename std::remove_reference<decltype(poreBlockList)>::type::value_type;
+ using GlobalPosition = typename PoreBlock::GlobalPosition;
+ GlobalPosition volumeSize;
+ volumeSize[0] = (porousMediaBox.second[0] - porousMediaBox.first[0]) / numVolumes[0];
+ volumeSize[1] = (porousMediaBox.second[1] - porousMediaBox.first[1]) / numVolumes[1];
+ const Scalar volumeLength = volumeSize[1];
+ const Scalar poreLength = std::sqrt(volumeLength * volumeLength - (volumeLength * volumeLength * porosity));
+ const Scalar poreOffset = (volumeLength - poreLength) / 2.0;
+
+ for (int i = 0; i < numVolumes[0]; i++)
+ {
+ for (int j = 0; j < numVolumes[1]; j++)
+ {
+ PoreBlock poreBlock;
+ poreBlock.gridBoxMin[0] = porousMediaBox.first[0] + (volumeSize[0] * i);
+ poreBlock.gridBoxMin[1] = porousMediaBox.first[1] + (volumeSize[1] * j);
+ poreBlock.gridBoxMax[0] = porousMediaBox.first[0] + (volumeSize[0] * (i+1));
+ poreBlock.gridBoxMax[1] = porousMediaBox.first[1] + (volumeSize[1] * (j+1));
+
+ poreBlock.poreBlockMin[0] = poreBlock.gridBoxMin[0] + (poreOffset);
+ poreBlock.poreBlockMax[0] = poreBlock.poreBlockMin[0] + poreLength;
+ poreBlock.poreBlockMin[1] = poreBlock.gridBoxMin[1] + (poreOffset);
+ poreBlock.poreBlockMax[1] = poreBlock.poreBlockMin[1] + poreLength;
+
+ poreBlockList.push_back(poreBlock);
+ }
+ }
+}
+
+template <class PoreBlockList, class Scalar, class DimensionPair, class BoundingBox>
+void buildStaggeredPoreGeometry(PoreBlockList& poreBlockList,
+ Scalar porosity,
+ DimensionPair numVolumes,
+ BoundingBox porousMediaBox)
+{
+ using PoreBlock = typename std::remove_reference<decltype(poreBlockList)>::type::value_type;
+ using GlobalPosition = typename PoreBlock::GlobalPosition;
+ GlobalPosition volumeSize;
+ numVolumes[1] = numVolumes[1] + 1;
+ porousMediaBox.first[1] = porousMediaBox.first[1] - ((porousMediaBox.second[0] - porousMediaBox.first[0]) / numVolumes[0]);
+ volumeSize[0] = (porousMediaBox.second[0] - porousMediaBox.first[0]) / numVolumes[0];
+ volumeSize[1] = (porousMediaBox.second[1] - porousMediaBox.first[1]) / numVolumes[1];
+ const Scalar volumeLength = volumeSize[1];
+ const Scalar poreLength = std::sqrt(volumeLength * volumeLength - (volumeLength * volumeLength * porosity));
+ const Scalar poreOffset = (volumeLength - poreLength) / 2.0;
+
+ // Shift staggered pores vertically
+ const Scalar staggeredOffset = volumeSize[1] / 2.0;
+
+ for (int i = 0; i < numVolumes[0]; i++)
+ {
+ for (int j = 0; j < numVolumes[1]; j++)
+ {
+ PoreBlock poreBlock;
+ poreBlock.gridBoxMin[0] = porousMediaBox.first[0] + (volumeSize[0] * i);
+ poreBlock.gridBoxMin[1] = porousMediaBox.first[1] + (volumeSize[1] * j);
+ poreBlock.gridBoxMax[0] = porousMediaBox.first[0] + (volumeSize[0] * (i+1));
+ poreBlock.gridBoxMax[1] = porousMediaBox.first[1] + (volumeSize[1] * (j+1));
+
+ poreBlock.poreBlockMin[0] = poreBlock.gridBoxMin[0] + (poreOffset);
+ poreBlock.poreBlockMax[0] = poreBlock.poreBlockMin[0] + poreLength;
+ poreBlock.poreBlockMin[1] = poreBlock.gridBoxMin[1] + (poreOffset);
+ poreBlock.poreBlockMax[1] = poreBlock.poreBlockMin[1] + poreLength;
+
+ if ( i % 2 != 0)
+ {
+ poreBlock.isStaggered = true;
+ poreBlock.poreBlockMin[1] = (poreBlock.poreBlockMin[1] + staggeredOffset);
+ poreBlock.poreBlockMax[1] = (poreBlock.poreBlockMax[1] + staggeredOffset);
+ }
+ poreBlockList.push_back(poreBlock);
+ }
+ }
+}
+
+
+template <class PoreBlockList, class Scalar, class DimensionPair, class BoundingBox>
+void buildRandomConstantPorosityPoreGeometry(PoreBlockList& poreBlockList,
+ Scalar porosity,
+ DimensionPair numVolumes,
+ BoundingBox porousMediaBox)
+{
+ using PoreBlock = typename std::remove_reference<decltype(poreBlockList)>::type::value_type;
+ using GlobalPosition = typename PoreBlock::GlobalPosition;
+ GlobalPosition volumeSize;
+ volumeSize[0] = (porousMediaBox.second[0] - porousMediaBox.first[0]) / numVolumes[0];
+ volumeSize[1] = (porousMediaBox.second[1] - porousMediaBox.first[1]) / numVolumes[1];
+ const Scalar volumeLength = volumeSize[1];
+ const Scalar poreLength = std::sqrt(volumeLength * volumeLength - (volumeLength * volumeLength * porosity));
+
+ for (int i = 0; i < numVolumes[0]; i++)
+ {
+ for (int j = 0; j < numVolumes[1]; j++)
+ {
+ PoreBlock poreBlock;
+ poreBlock.gridBoxMin[0] = porousMediaBox.first[0] + (volumeSize[0] * i);
+ poreBlock.gridBoxMin[1] = porousMediaBox.first[1] + (volumeSize[1] * j);
+ poreBlock.gridBoxMax[0] = porousMediaBox.first[0] + (volumeSize[0] * (i+1));
+ poreBlock.gridBoxMax[1] = porousMediaBox.first[1] + (volumeSize[1] * (j+1));
+ poreBlock.centerPoint[0] = (poreBlock.gridBoxMax[0] + poreBlock.gridBoxMin[0]) / 2.0;
+ poreBlock.centerPoint[1] = (poreBlock.gridBoxMax[1] + poreBlock.gridBoxMin[1]) / 2.0;
+
+ double r1 = ( (double) std::rand() / (RAND_MAX));
+ poreBlock.poreBlockMin[0] = poreBlock.gridBoxMin[0] + ((poreBlock.centerPoint[0] - poreBlock.gridBoxMin[0]) * r1);
+ poreBlock.poreBlockMax[0] = poreBlock.poreBlockMin[0] + poreLength;
+ double r2 = ( (double) std::rand() / (RAND_MAX));
+ poreBlock.poreBlockMin[1] = poreBlock.gridBoxMin[1] + ((poreBlock.centerPoint[1] - poreBlock.gridBoxMin[1]) * r2);
+ poreBlock.poreBlockMax[1] = poreBlock.poreBlockMin[1] + poreLength;
+
+ poreBlockList.push_back(poreBlock);
+ }
+ }
+}
+
+template <class PoreBlockList, class Scalar, class DimensionPair, class BoundingBox>
+void buildRandomPoreGeometry(PoreBlockList& poreBlockList,
+ Scalar porosity,
+ DimensionPair numVolumes,
+ BoundingBox porousMediaBox)
+{
+ using PoreBlock = typename std::remove_reference<decltype(poreBlockList)>::type::value_type;
+ using GlobalPosition = typename PoreBlock::GlobalPosition;
+ GlobalPosition volumeSize;
+ volumeSize[0] = (porousMediaBox.second[0] - porousMediaBox.first[0]) / numVolumes[0];
+ volumeSize[1] = (porousMediaBox.second[1] - porousMediaBox.first[1]) / numVolumes[1];
+
+ for (int i = 0; i < numVolumes[0]; i++)
+ {
+ for (int j = 0; j < numVolumes[1]; j++)
+ {
+ PoreBlock poreBlock;
+ poreBlock.gridBoxMin[0] = porousMediaBox.first[0] + (volumeSize[0] * i);
+ poreBlock.gridBoxMin[1] = porousMediaBox.first[1] + (volumeSize[1] * j);
+ poreBlock.gridBoxMax[0] = porousMediaBox.first[0] + (volumeSize[0] * (i+1));
+ poreBlock.gridBoxMax[1] = porousMediaBox.first[1] + (volumeSize[1] * (j+1));
+ poreBlock.centerPoint[0] = (poreBlock.gridBoxMax[0] + poreBlock.gridBoxMin[0]) / 2.0;
+ poreBlock.centerPoint[1] = (poreBlock.gridBoxMax[1] + poreBlock.gridBoxMin[1]) / 2.0;
+
+ double r1 = ( (double) std::rand() / (RAND_MAX));
+ poreBlock.poreBlockMin[0] = poreBlock.gridBoxMin[0] + ((poreBlock.centerPoint[0] - poreBlock.gridBoxMin[0]) * r1);
+ double r2 = ( (double) std::rand() / (RAND_MAX));
+ poreBlock.poreBlockMin[1] = poreBlock.gridBoxMin[1] + ((poreBlock.centerPoint[1] - poreBlock.gridBoxMin[1]) * r2);
+
+ double r3 = ( (double) std::rand() / (RAND_MAX));
+ poreBlock.poreBlockMax[0] = poreBlock.centerPoint[0] + ((poreBlock.gridBoxMax[0] - poreBlock.centerPoint[0]) * r3);
+ double r4 = ( (double) std::rand() / (RAND_MAX));
+ poreBlock.poreBlockMax[1] = poreBlock.centerPoint[1] + ((poreBlock.gridBoxMax[1] - poreBlock.centerPoint[1]) * r4);
+
+ poreBlockList.push_back(poreBlock);
+ }
+ }
+}
+
+template <class PoreBlockList>
+void readInPoreGeometry(PoreBlockList& poreBlockList,
+ std::string gridInputFileName)
+{
+ // Open the file and make sure that the file is open
+ std::ifstream infile(gridInputFileName, std::ios::in);
+ if (!infile)
+ DUNE_THROW(Dune::IOError, "Unable to open file: " << gridInputFileName << " !\n");
+
+ using PoreBlock = typename std::remove_reference<decltype(poreBlockList)>::type::value_type;
+ using Scalar = typename PoreBlock::Scalar;
+
+ std::vector<std::vector<Scalar>> coodinates;
+ int lineIdx = 0;
+ while (infile)
+ {
+ lineIdx++;
+ std::string s;
+ // If there is no further line, end loop
+ if (!getline(infile, s))
+ break;
+
+ // Read each line to a string
+ if (s[0] != '#') // ignore lines that start with the comment indicator (#)
+ {
+ std::istringstream ss(s);
+ std::vector<Scalar> blockRecord;
+
+ while (ss)
+ {
+ std::string line;
+ if (!getline(ss, line, ',')) // delimited with a comma
+ break;
+
+ blockRecord.push_back(std::stof(line));
+ }
+
+ coodinates.push_back(blockRecord);
+ }
+ }
+
+ if (!infile.eof())
+ DUNE_THROW(Dune::IOError, "Could not read file " << gridInputFileName << "!\n");
+
+ for (int i = 0; i < coodinates.size(); i++)
+ {
+ PoreBlock poreBlock;
+ poreBlock.poreBlockMin[0] = coodinates[i][0];
+ poreBlock.poreBlockMax[0] = coodinates[i][1];
+ poreBlock.poreBlockMin[1] = coodinates[i][2];
+ poreBlock.poreBlockMax[1] = coodinates[i][3];
+ poreBlockList.push_back(poreBlock);
+ }
+ std::cout << "The Pore geometry has been sucessfully created from the provided grid file. \n";
+}
+
+
+template <class PoreBlockList>
+void writeOutPoreGeometryFileToCSV(const PoreBlockList& poreBlockList,
+ std::string gridOutputFileName)
+{
+ std::string gridListFileName = "PoreGeometry_" + gridOutputFileName + ".csv";
+ std::ofstream gridOutputFile;
+
+ gridOutputFile.open(gridListFileName, std::ios::app);
+ gridOutputFile << "#PoreMin[0]" << ","
+ << "PoreMax[0]" << ","
+ << "PoreMin[1]" << ","
+ << "PoreMax[1]" << "\n";
+ gridOutputFile.close();// Export the coordinates
+
+ for (int i = 0; i < poreBlockList.size(); i++)
+ {
+ gridOutputFile.open(gridListFileName, std::ios::app);
+ gridOutputFile << poreBlockList[i].poreBlockMin[0] << ","
+ << poreBlockList[i].poreBlockMax[0] << ","
+ << poreBlockList[i].poreBlockMin[1] << ","
+ << poreBlockList[i].poreBlockMax[1] << "\n";
+ gridOutputFile.close();
+ }
+ std::cout << "The Pore geometry description is writen to the " << gridListFileName << " file. \n";
+}
+
+} // end namespace Dumux::PoreGeometries
+
+#endif