ComputeNodeAreasFromSurfaceMesh.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 <fstream>
#include <memory>
#include <numeric>
#include <string>
#include <vector>
 
#include "BaseLib/Error.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/MPI.h"
#include "BaseLib/TCLAPArguments.h"
#include "InfoLib/GitInfo.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshSurfaceExtraction.h"
 
static void writeToFile(
    std::string const& id_area_fname, std::string const& csv_fname,
    std::vector<std::pair<std::size_t, double>> const& ids_and_areas,
    std::vector<MeshLib::Node*> const& mesh_nodes)
{
    std::ofstream ids_and_area_out(id_area_fname);
    if (!ids_and_area_out)
    {
        OGS_FATAL("Unable to open the file '{:s}' - aborting.", id_area_fname);
    }
    std::ofstream csv_out(csv_fname);
    if (!csv_out)
    {
        OGS_FATAL("Unable to open the file '{:s}' - aborting.", csv_fname);
    }
 
    ids_and_area_out.precision(std::numeric_limits<double>::max_digits10);
    csv_out.precision(std::numeric_limits<double>::max_digits10);
 
    csv_out << "ID x y z area node_id\n";  // CSV header
    for (std::size_t k(0); k < ids_and_areas.size(); k++)
    {
        ids_and_area_out << ids_and_areas[k].first << " "
                         << ids_and_areas[k].second << "\n";
        csv_out << k << " " << *(mesh_nodes[k]) << ids_and_areas[k].second
                << " " << ids_and_areas[k].first << "\n";
    }
    ids_and_area_out << "\n";
    csv_out << "\n";
}
 
int main(int argc, char* argv[])
{
    TCLAP::CmdLine cmd(
        "The tool computes the area per node of the surface mesh and writes "
        "the information as txt and csv data.\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);
    TCLAP::ValueArg<std::string> mesh_in(
        "i", "mesh-input-file",
        "Input (.vtu). The name of the file containing the input mesh", true,
        "", "INPUT_FILE");
    cmd.add(mesh_in);
 
    std::vector<std::string> allowed_types_vector{
        "bulk_node_ids", "bulk_element_ids", "bulk_edge_ids", "bulk_face_ids"};
    TCLAP::ValuesConstraint<std::string> allowed_types(allowed_types_vector);
    TCLAP::ValueArg<std::string> id_prop_name(
        "", "id-prop-name",
        "the name of the property containing the id information", false,
        "bulk_node_ids", &allowed_types);
    cmd.add(id_prop_name);
    TCLAP::ValueArg<std::string> out_base_fname(
        "p", "output-base-name",
        "Output (.csv | .txt). The path and base file name the output will be "
        "written to",
        false, "", "BASE_FILENAME_OUTPUT");
    cmd.add(out_base_fname);
 
    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());
 
    std::unique_ptr<MeshLib::Mesh> surface_mesh(
        MeshLib::IO::readMeshFromFile(mesh_in.getValue()));
    INFO("Mesh read: {:d} nodes, {:d} elements.",
         surface_mesh->getNumberOfNodes(), surface_mesh->getNumberOfElements());
    // ToDo check if mesh is read correct and if the mesh is a surface mesh
 
    MeshLib::PropertyVector<std::size_t>* orig_node_ids(nullptr);
    // check if a node property containing the subsurface ids is available
    // if the node property is not available generate it
    if (!surface_mesh->getProperties().existsPropertyVector<std::size_t>(
            id_prop_name.getValue()))
    {
        orig_node_ids =
            surface_mesh->getProperties().createNewPropertyVector<std::size_t>(
                id_prop_name.getValue(), MeshLib::MeshItemType::Node,
                surface_mesh->getNumberOfNodes(), 1);
        if (!orig_node_ids)
        {
            ERR("Fatal error: could not create property.");
            return EXIT_FAILURE;
        }
        std::iota(orig_node_ids->begin(), orig_node_ids->end(), 0);
    }
    else
    {
        orig_node_ids =
            surface_mesh->getProperties().getPropertyVector<std::size_t>(
                id_prop_name.getValue());
    }
 
    std::vector<double> areas(
        MeshToolsLib::MeshSurfaceExtraction::getSurfaceAreaForNodes(
            *surface_mesh));
 
    // pack area and node id together
    std::vector<std::pair<std::size_t, double>> ids_and_areas;
    std::transform(orig_node_ids->cbegin(), orig_node_ids->cend(),
                   areas.cbegin(), std::back_inserter(ids_and_areas),
                   std::make_pair<std::size_t const&, double const&>);
 
    // generate file names for output
    std::string path(out_base_fname.getValue());
    if (path.empty())
    {
        path = BaseLib::dropFileExtension(mesh_in.getValue());
    }
    std::string const id_and_area_fname(path + ".txt");
    std::string const csv_fname(path + ".csv");
 
    writeToFile(id_and_area_fname, csv_fname, ids_and_areas,
                surface_mesh->getNodes());
 
    return EXIT_SUCCESS;
}