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