OGS
MeshMapping.cpp File Reference

Detailed Description

2012/03/07 KR Initial implementation

Definition in file MeshMapping.cpp.

#include <tclap/CmdLine.h>
#include <fstream>
#include <memory>
#include <string>
#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/MPI.h"
#include "BaseLib/TCLAPArguments.h"
#include "GeoLib/AABB.h"
#include "GeoLib/IO/AsciiRasterInterface.h"
#include "GeoLib/Raster.h"
#include "InfoLib/GitInfo.h"
#include "MathLib/MathTools.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/IO/writeMeshToFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshSearch/MeshElementGrid.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshEditing/ProjectPointOnMesh.h"
#include "MeshToolsLib/MeshGenerators/MeshLayerMapper.h"
Include dependency graph for MeshMapping.cpp:

Go to the source code of this file.

Functions

double getClosestPointElevation (MeshLib::Node const &p, std::vector< MeshLib::Node * > const &nodes, double const &max_dist)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ getClosestPointElevation()

double getClosestPointElevation ( MeshLib::Node const & p,
std::vector< MeshLib::Node * > const & nodes,
double const & max_dist )

Definition at line 35 of file MeshMapping.cpp.

38{
39 double sqr_shortest_dist(max_dist);
40 double elevation(p[2]);
41 for (MeshLib::Node* node : nodes)
42 {
43 double sqr_dist = (p[0] - (*node)[0]) * (p[0] - (*node)[0]) +
44 (p[1] - (*node)[1]) * (p[1] - (*node)[1]);
45 if (sqr_dist < sqr_shortest_dist)
46 {
47 sqr_shortest_dist = sqr_dist;
48 elevation = (*node)[2];
49 }
50 }
51 return elevation;
52}

Referenced by main().

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 54 of file MeshMapping.cpp.

55{
56 TCLAP::CmdLine cmd(
57 "Changes the elevation of 2D mesh nodes based on either raster data or "
58 "another 2D mesh. In addition, a low pass filter can be applied to "
59 "node elevation based connected nodes.\n\n"
60 "OpenGeoSys-6 software, version " +
62 ".\n"
63 "Copyright (c) 2012-2025, OpenGeoSys Community "
64 "(http://www.opengeosys.org)",
66 auto log_level_arg = BaseLib::makeLogLevelArg();
67 cmd.add(log_level_arg);
68 TCLAP::SwitchArg lowpass_arg(
69 "", "lowpass",
70 "Applies a lowpass filter to elevation over connected nodes.", false);
71 cmd.add(lowpass_arg);
72 TCLAP::ValueArg<double> map_static_arg(
73 "s", "static",
74 "Static elevation to map the input file to. This can be combined with "
75 "mapping based on rasters or other meshes to deal with locations where "
76 "no corresponding data exists.",
77 false, 0, "MAP_STATIC");
78 cmd.add(map_static_arg);
79 TCLAP::ValueArg<double> max_dist_arg(
80 "d", "distance",
81 "Maximum distance to search for mesh nodes if there is no "
82 "corresponding data for input mesh nodes on the mesh it should be "
83 "mapped on. (min = 0)",
84 false, 1, "MAX_DISTANCE");
85 cmd.add(max_dist_arg);
86 TCLAP::ValueArg<std::string> map_mesh_arg(
87 "m", "mesh", "Input (.vtu). 2D mesh file to map the input file on",
88 false, "", "INPUT_FILE");
89 cmd.add(map_mesh_arg);
90 TCLAP::ValueArg<std::string> map_raster_arg(
91 "r", "raster",
92 "Input (.asc | .grd | .xyz) Raster file to map the input file on",
93 false, "", "INPUT_FILE");
94 cmd.add(map_raster_arg);
95 TCLAP::ValueArg<std::string> output_arg(
96 "o", "output", "Output (.vtu) mesh file", true, "", "OUTPUT_FILE");
97 cmd.add(output_arg);
98 TCLAP::ValueArg<std::string> input_arg(
99 "i", "input", "Input (.vtu | .msh) mesh file", true, "", "INPUT_FILE");
100 cmd.add(input_arg);
101 cmd.parse(argc, argv);
102
103 BaseLib::MPI::Setup mpi_setup(argc, argv);
104 BaseLib::initOGSLogger(log_level_arg.getValue());
105
106 std::unique_ptr<MeshLib::Mesh> mesh(
107 MeshLib::IO::readMeshFromFile(input_arg.getValue()));
108 if (mesh == nullptr)
109 {
110 ERR("Error reading mesh file.");
111 return EXIT_FAILURE;
112 }
113
114 if (!(map_static_arg.isSet() || map_raster_arg.isSet() ||
115 map_mesh_arg.isSet()))
116 {
117 ERR("Nothing to do. Please choose mapping based on a raster or mesh "
118 "file, or to a static value.");
119 return EXIT_FAILURE;
120 }
121
122 if (map_raster_arg.isSet() && map_mesh_arg.isSet())
123 {
124 ERR("Please select mapping based on *either* a mesh or a raster file.");
125 return EXIT_FAILURE;
126 }
127
128 // Maps the elevation of mesh nodes to static value
129 if (map_static_arg.isSet())
130 {
132 *mesh, map_static_arg.getValue());
133 }
134
135 // Maps the elevation of mesh nodes according to raster
136 if (map_raster_arg.isSet())
137 {
138 std::string const raster_path = map_raster_arg.getValue();
139 std::ifstream file_stream(raster_path, std::ifstream::in);
140 if (!file_stream.good())
141 {
142 ERR("Opening raster file {} failed.", raster_path);
143 return EXIT_FAILURE;
144 }
145 file_stream.close();
146
147 std::unique_ptr<GeoLib::Raster> const raster(
149 MeshToolsLib::MeshLayerMapper::layerMapping(*mesh, *raster, 0, true);
150 }
151
152 // Maps the elevation of mesh nodes according to a ground truth mesh
153 if (map_mesh_arg.isSet())
154 {
155 std::unique_ptr<MeshLib::Mesh> ground_truth(
156 MeshLib::IO::readMeshFromFile(map_mesh_arg.getValue()));
157 if (ground_truth == nullptr)
158 {
159 ERR("Error reading mesh file.");
160 return EXIT_FAILURE;
161 }
162
163 std::vector<MeshLib::Node*> const& nodes = mesh->getNodes();
164 MeshLib::MeshElementGrid const grid(*ground_truth);
165 auto const edgeLengths = minMaxEdgeLength(mesh->getElements());
166 double const max_edge = edgeLengths.second;
167 double const max_dist(pow(max_dist_arg.getValue(), 2));
168
169 for (MeshLib::Node* node : nodes)
170 {
171 MathLib::Point3d min_vol{{(*node)[0] - max_edge,
172 (*node)[1] - max_edge,
173 -std::numeric_limits<double>::max()}};
174 MathLib::Point3d max_vol{{(*node)[0] + max_edge,
175 (*node)[1] + max_edge,
176 std::numeric_limits<double>::max()}};
177 std::vector<const MeshLib::Element*> const& elems =
178 grid.getElementsInVolume(min_vol, max_vol);
179 auto const* element =
181 *node);
182 if (element != nullptr)
183 {
185 *element, *node);
186 }
187 else
188 {
189 (*node)[2] = getClosestPointElevation(
190 *node, ground_truth->getNodes(), max_dist);
191 }
192 }
193 }
194
195 // a simple lowpass filter for the elevation of mesh nodes using the
196 // elevation of each node weighted by 2 and the elevation of each connected
197 // node weighted by 1
198 if (lowpass_arg.isSet())
199 {
200 INFO("lowpass");
201 const std::size_t nNodes(mesh->getNumberOfNodes());
202 std::vector<MeshLib::Node*> nodes(mesh->getNodes());
203
204 std::vector<double> elevation(nNodes);
205 for (std::size_t i = 0; i < nNodes; i++)
206 {
207 elevation[i] = (*nodes[i])[2];
208 }
209
210 auto const& connections =
212 for (std::size_t i = 0; i < nNodes; i++)
213 {
214 auto const& conn_nodes(connections[nodes[i]->getID()]);
215 const unsigned nConnNodes(conn_nodes.size());
216 elevation[i] = (2 * (*nodes[i])[2]);
217 for (std::size_t j = 0; j < nConnNodes; ++j)
218 {
219 elevation[i] += (*conn_nodes[j])[2];
220 }
221 elevation[i] /= (nConnNodes + 2);
222 }
223
224 for (std::size_t i = 0; i < nNodes; i++)
225 {
226 (*nodes[i])[2] = elevation[i];
227 }
228 }
229
230 if (MeshLib::IO::writeMeshToFile(*mesh, output_arg.getValue()) != 0)
231 {
232 return EXIT_FAILURE;
233 }
234
235 INFO("Result successfully written.");
236 return EXIT_SUCCESS;
237}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
double getClosestPointElevation(MeshLib::Node const &p, std::vector< MeshLib::Node * > const &nodes, double const &max_dist)
static GeoLib::Raster * readRaster(std::string const &fname)
static bool layerMapping(MeshLib::Mesh const &mesh, const GeoLib::Raster &raster, double nodata_replacement=0.0, bool const ignore_nodata=false)
static bool mapToStaticValue(MeshLib::Mesh const &mesh, double value)
Maps the elevation of all mesh nodes to the specified static value.
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:64
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
int writeMeshToFile(const MeshLib::Mesh &mesh, std::filesystem::path const &file_path, std::set< std::string > variable_output_names)
std::vector< std::vector< Node * > > calculateNodesConnectedByElements(Mesh const &mesh)
Definition Mesh.cpp:309
std::pair< double, double > minMaxEdgeLength(std::vector< Element * > const &elements)
Returns the minimum and maximum edge length for given elements.
Definition Mesh.cpp:190
double getElevation(MeshLib::Element const &element, MathLib::Point3d const &node)
MeshLib::Element const * getProjectedElement(std::vector< const MeshLib::Element * > const &elements, MathLib::Point3d const &node)

References MeshLib::calculateNodesConnectedByElements(), ERR(), getClosestPointElevation(), MeshLib::MeshElementGrid::getElementsInVolume(), MeshToolsLib::ProjectPointOnMesh::getElevation(), MeshToolsLib::ProjectPointOnMesh::getProjectedElement(), INFO(), BaseLib::initOGSLogger(), MeshToolsLib::MeshLayerMapper::layerMapping(), BaseLib::makeLogLevelArg(), MeshToolsLib::MeshLayerMapper::mapToStaticValue(), GitInfoLib::GitInfo::ogs_version, MeshLib::IO::readMeshFromFile(), FileIO::AsciiRasterInterface::readRaster(), and MeshLib::IO::writeMeshToFile().