[PoreGeometry] move geometry header to base

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)

dumux/CMakeLists.txt

+#install headers
+install(FILES
+poregeometries.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/)

dumux/poregeometries.hh

+// -*- 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