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 <mpi.h>
#include "BaseLib/FileTools.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 36 of file MeshMapping.cpp.

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

Referenced by main().

◆ main()

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

Definition at line 55 of file MeshMapping.cpp.

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