generateStructuredMesh.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include <tclap/CmdLine.h>
 
#include <algorithm>
#include <memory>
#include <string>
#include <vector>
 
#include "BaseLib/Error.h"
#include "BaseLib/Logging.h"
#include "BaseLib/MPI.h"
#include "BaseLib/Subdivision.h"
#include "BaseLib/TCLAPArguments.h"
#include "InfoLib/GitInfo.h"
#include "MathLib/Point3d.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/IO/writeMeshToFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshEnums.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshGenerators/MeshGenerator.h"
 
namespace
{
unsigned getDimension(MeshLib::MeshElemType eleType)
{
    switch (eleType)
    {
        case MeshLib::MeshElemType::LINE:
            return 1;
        case MeshLib::MeshElemType::QUAD:
        case MeshLib::MeshElemType::TRIANGLE:
            return 2;
        case MeshLib::MeshElemType::HEXAHEDRON:
        case MeshLib::MeshElemType::PRISM:
        case MeshLib::MeshElemType::PYRAMID:
        case MeshLib::MeshElemType::TETRAHEDRON:
            return 3;
        case MeshLib::MeshElemType::POINT:
        case MeshLib::MeshElemType::INVALID:
            return 0;
    }
    return 0;
}
 
}  // end namespace
 
int main(int argc, char* argv[])
{
    TCLAP::CmdLine cmd(
        "Structured mesh generator.\n"
        "The documentation is available at "
        "https://docs.opengeosys.org/docs/tools/meshing/"
        "structured-mesh-generation.\n\n"
        "OpenGeoSys-6 software, version " +
            GitInfoLib::GitInfo::ogs_version +
            ".\n"
            "Copyright (c) 2012-2026, OpenGeoSys Community "
            "(http://www.opengeosys.org)",
        ' ', GitInfoLib::GitInfo::ogs_version);
 
    std::vector<std::string> allowed_ele_types;
    allowed_ele_types.emplace_back("line");
    allowed_ele_types.emplace_back("tri");
    allowed_ele_types.emplace_back("quad");
    allowed_ele_types.emplace_back("hex");
    allowed_ele_types.emplace_back("prism");
    allowed_ele_types.emplace_back("tet");
    allowed_ele_types.emplace_back("pyramid");
    TCLAP::ValuesConstraint<std::string> allowedVals(allowed_ele_types);
    TCLAP::ValueArg<std::string> eleTypeArg("e", "element-type",
                                            "element type to be created", true,
                                            "line", &allowedVals);
    cmd.add(eleTypeArg);
    TCLAP::ValueArg<std::string> mesh_out(
        "o", "mesh-output-file",
        "Output (.vtu). The name of the file the mesh will be written to", true,
        "", "OUTPUT_FILE");
    cmd.add(mesh_out);
    TCLAP::ValueArg<double> lengthXArg("", "lx",
                                       "length of a domain in x direction, "
                                       "(min = 0)",
                                       false, 10.0, "LENGTH_X");
    cmd.add(lengthXArg);
    TCLAP::ValueArg<double> lengthYArg("", "ly",
                                       "length of a domain in y direction, "
                                       "(min = 0)",
                                       false, 10.0, "LENGTH_Y");
    cmd.add(lengthYArg);
    TCLAP::ValueArg<double> lengthZArg("", "lz",
                                       "length of a domain in z direction, "
                                       "(min = 0)",
                                       false, 10.0, "LENGTH_Z");
    cmd.add(lengthZArg);
    TCLAP::ValueArg<unsigned> nsubdivXArg(
        "", "nx",
        "the number of subdivision in x direction, "
        "(min = 0)",
        false, 10, "SUBDIVISIONS_X");
    cmd.add(nsubdivXArg);
    TCLAP::ValueArg<unsigned> nsubdivYArg(
        "", "ny",
        "the number of subdivision in y direction, "
        "(min = 0)",
        false, 10, "SUBDIVISIONS_Y");
    cmd.add(nsubdivYArg);
    TCLAP::ValueArg<unsigned> nsubdivZArg(
        "", "nz",
        "the number of subdivision in z direction, "
        "(min = 0)",
        false, 10, "SUBDIVISIONS_Z");
    cmd.add(nsubdivZArg);
    // in case of gradual refinement
    TCLAP::ValueArg<double> d0XArg("", "dx0",
                                   "initial cell length in x direction, "
                                   "(min = 0)",
                                   false, 1, "INITIAL_X");
    cmd.add(d0XArg);
    TCLAP::ValueArg<double> d0YArg("", "dy0",
                                   "initial cell length in y direction, "
                                   "(min = 0)",
                                   false, 1, "INITIAL_Y");
    cmd.add(d0YArg);
    TCLAP::ValueArg<double> d0ZArg("", "dz0",
                                   "initial cell length in z direction, "
                                   "(min = 0)",
                                   false, 1, "INITIAL_Z");
    cmd.add(d0ZArg);
    TCLAP::ValueArg<double> dmaxXArg("", "dx-max",
                                     "maximum cell length in x direction, "
                                     "(min = 0)",
                                     false, std::numeric_limits<double>::max(),
                                     "MAX_X");
    cmd.add(dmaxXArg);
    TCLAP::ValueArg<double> dmaxYArg("", "dy-max",
                                     "maximum cell length in y direction, "
                                     "(min = 0)",
                                     false, std::numeric_limits<double>::max(),
                                     "MAX_Y");
    cmd.add(dmaxYArg);
    TCLAP::ValueArg<double> dmaxZArg("", "dz-max",
                                     "maximum cell length in z direction, "
                                     "(min = 0)",
                                     false, std::numeric_limits<double>::max(),
                                     "MAX_Z");
    cmd.add(dmaxZArg);
    TCLAP::ValueArg<double> multiXArg("", "mx", "multiplier in x direction",
                                      false, 1, "MULTIPLIER_X");
    cmd.add(multiXArg);
    TCLAP::ValueArg<double> multiYArg("", "my", "multiplier in y direction",
                                      false, 1, "MULTIPLIER_Y");
    cmd.add(multiYArg);
    TCLAP::ValueArg<double> multiZArg("", "mz", "multiplier in z direction",
                                      false, 1, "MULTIPLIER_Z");
    cmd.add(multiZArg);
    TCLAP::ValueArg<double> originXArg(
        "", "ox", "mesh origin (lower left corner) in x direction", false, 0,
        "ORIGIN_X");
    cmd.add(originXArg);
    TCLAP::ValueArg<double> originYArg(
        "", "oy", "mesh origin (lower left corner) in y direction", false, 0,
        "ORIGIN_Y");
    cmd.add(originYArg);
    TCLAP::ValueArg<double> originZArg(
        "", "oz", "mesh origin (lower left corner) in z direction", false, 0,
        "ORIGIN_Z");
    cmd.add(originZArg);
 
    // parse arguments
    auto log_level_arg = BaseLib::makeLogLevelArg();
    cmd.add(log_level_arg);
    cmd.parse(argc, argv);
    BaseLib::MPI::Setup mpi_setup(argc, argv);
    BaseLib::initOGSLogger(log_level_arg.getValue());
 
    const std::string eleTypeName(eleTypeArg.getValue());
    const MeshLib::MeshElemType eleType =
        MeshLib::String2MeshElemType(eleTypeName);
    const unsigned dim = getDimension(eleType);
 
    bool dim_used[3] = {false};
    for (unsigned i = 0; i < dim; i++)
    {
        dim_used[i] = true;
    }
 
    std::vector<TCLAP::ValueArg<double>*> vec_lengthArg = {
        &lengthXArg, &lengthYArg, &lengthZArg};
    std::vector<TCLAP::ValueArg<unsigned>*> vec_ndivArg = {
        &nsubdivXArg, &nsubdivYArg, &nsubdivZArg};
    std::vector<TCLAP::ValueArg<double>*> vec_d0Arg = {&d0XArg, &d0YArg,
                                                       &d0ZArg};
    std::vector<TCLAP::ValueArg<double>*> vec_dMaxArg = {&dmaxXArg, &dmaxYArg,
                                                         &dmaxZArg};
    std::vector<TCLAP::ValueArg<double>*> vec_multiArg = {
        &multiXArg, &multiYArg, &multiZArg};
    MathLib::Point3d const origin(
        {originXArg.getValue(), originYArg.getValue(), originZArg.getValue()});
 
    const bool isLengthSet =
        std::any_of(vec_lengthArg.begin(), vec_lengthArg.end(),
                    [&](TCLAP::ValueArg<double>* arg) { return arg->isSet(); });
    if (!isLengthSet)
    {
        ERR("Missing input: Length information is not provided at all.");
        return EXIT_FAILURE;
    }
    for (unsigned i = 0; i < 3; i++)
    {
        if (dim_used[i] && !vec_lengthArg[i]->isSet())
        {
            ERR("Missing input: Length for dimension [{:d}] is required but "
                "missing.",
                i);
            return EXIT_FAILURE;
        }
    }
 
    std::vector<double> length(dim);
    std::vector<unsigned> n_subdivision(dim);
    std::vector<double> vec_dx(dim);
    for (unsigned i = 0; i < dim; i++)
    {
        length[i] = vec_lengthArg[i]->getValue();
        n_subdivision[i] = vec_ndivArg[i]->getValue();
        vec_dx[i] = length[i] / n_subdivision[i];
    }
 
    std::vector<std::unique_ptr<BaseLib::ISubdivision>> vec_div;
    vec_div.reserve(dim);
    for (unsigned i = 0; i < dim; i++)
    {
        if (vec_multiArg[i]->isSet())
        {
            if (vec_ndivArg[i]->isSet())
            {
                // number of partitions in direction is specified
                if (vec_d0Arg[i]->isSet())
                {
                    OGS_FATAL(
                        "Specifying all of --m?, --d?0 and --n? for coordinate "
                        "'?' is not supported.");
                }
                vec_div.emplace_back(new BaseLib::GradualSubdivisionFixedNum(
                    length[i], vec_ndivArg[i]->getValue(),
                    vec_multiArg[i]->getValue()));
            }
            else
            {
                vec_div.emplace_back(new BaseLib::GradualSubdivision(
                    length[i], vec_d0Arg[i]->getValue(),
                    vec_dMaxArg[i]->getValue(), vec_multiArg[i]->getValue()));
            }
        }
        else
        {
            vec_div.emplace_back(
                new BaseLib::UniformSubdivision(length[i], n_subdivision[i]));
        }
    }
 
    // generate a mesh
    std::unique_ptr<MeshLib::Mesh> mesh;
    switch (eleType)
    {
        case MeshLib::MeshElemType::LINE:
            mesh.reset(MeshToolsLib::MeshGenerator::generateLineMesh(
                *vec_div[0], origin));
            break;
        case MeshLib::MeshElemType::TRIANGLE:
            mesh.reset(MeshToolsLib::MeshGenerator::generateRegularTriMesh(
                *vec_div[0], *vec_div[1], origin));
            break;
        case MeshLib::MeshElemType::QUAD:
            mesh.reset(MeshToolsLib::MeshGenerator::generateRegularQuadMesh(
                *vec_div[0], *vec_div[1], origin));
            break;
        case MeshLib::MeshElemType::HEXAHEDRON:
            mesh.reset(MeshToolsLib::MeshGenerator::generateRegularHexMesh(
                *vec_div[0], *vec_div[1], *vec_div[2], origin));
            break;
        case MeshLib::MeshElemType::PRISM:
            mesh.reset(MeshToolsLib::MeshGenerator::generateRegularPrismMesh(
                length[0], length[1], length[2], n_subdivision[0],
                n_subdivision[1], n_subdivision[2], origin));
            break;
        case MeshLib::MeshElemType::TETRAHEDRON:
            mesh.reset(MeshToolsLib::MeshGenerator::generateRegularTetMesh(
                *vec_div[0], *vec_div[1], *vec_div[2], origin));
            break;
        case MeshLib::MeshElemType::PYRAMID:
            mesh.reset(MeshToolsLib::MeshGenerator::generateRegularPyramidMesh(
                *vec_div[0], *vec_div[1], *vec_div[2], origin));
            break;
        default:
            ERR("Given element type is not supported.");
            break;
    }
 
    if (mesh)
    {
        INFO("Mesh created: {:d} nodes, {:d} elements.",
             mesh->getNumberOfNodes(), mesh->getNumberOfElements());
 
        // write into a file
        MeshLib::IO::writeMeshToFile(
            *(mesh.get()), std::filesystem::path(mesh_out.getValue()));
    }
 
    return EXIT_SUCCESS;
}