OGS
IntegrateBoreholesIntoMesh.cpp
Go to the documentation of this file.
1
10#include <tclap/CmdLine.h>
11
12#include <iterator>
13#include <limits>
14#include <memory>
15#include <string>
16#include <vector>
17
18#include "BaseLib/Logging.h"
19#include "BaseLib/MPI.h"
21#include "GeoLib/GEOObjects.h"
23#include "InfoLib/GitInfo.h"
28#include "MeshLib/Mesh.h"
29#include "MeshLib/Node.h"
31
32std::vector<std::size_t> getNodes(
33 GeoLib::Point const& pnt, std::vector<MeshLib::Node*> const& nodes,
34 MeshLib::PropertyVector<int> const& mat_ids,
35 std::pair<int, int> const& mat_limits,
36 std::pair<double, double> const& elevation_limits,
37 MeshLib::Mesh const& mesh)
38{
39 std::vector<std::size_t> pnt_nodes;
40 for (auto node : nodes)
41 {
42 double const eps = std::numeric_limits<double>::epsilon();
43 if (std::abs((*node)[0] - pnt[0]) < eps &&
44 std::abs((*node)[1] - pnt[1]) < eps)
45 {
46 auto const& elems = mesh.getElementsConnectedToNode(*node);
47 for (auto e : elems)
48 {
49 if (mat_ids[e->getID()] >= mat_limits.first &&
50 mat_ids[e->getID()] <= mat_limits.second &&
51 (*node)[2] >= elevation_limits.first &&
52 (*node)[2] <= elevation_limits.second)
53 {
54 pnt_nodes.push_back(node->getID());
55 break;
56 }
57 }
58 }
59 }
60 if (pnt_nodes.size() < 2)
61 {
62 WARN("No nodes found for point {:d}...", pnt.getID());
63 }
64 else
65 {
66 // sort found nodes from top to bottom (required for BHE simulations)
67 std::sort(pnt_nodes.begin(), pnt_nodes.end(),
68 [nodes](std::size_t a, std::size_t b)
69 { return (*nodes[a])[2] > (*nodes[b])[2]; });
70 }
71 return pnt_nodes;
72}
73
74int main(int argc, char* argv[])
75{
76 TCLAP::CmdLine cmd(
77 "Integrates line elements representing boreholes into a pre-existing "
78 "3D mesh. Corresponding nodes matching the (x,y)-coordinates given in "
79 "the gml-file are found in the mesh and connected from top to bottom "
80 "via line elements. Each borehole (i.e. all points at a given "
81 "(x,y)-location but at different depths) is assigned a unique material "
82 "ID. Vertical limits of boreholes can be specified via Material IDs "
83 "and/or elevation. Point not matching any mesh nodes or located outside"
84 " the mesh are ignored.\n\n"
85 "OpenGeoSys-6 software, version " +
87 ".\n"
88 "Copyright (c) 2012-2025, OpenGeoSys Community "
89 "(http://www.opengeosys.org)",
91
92 double const dmax = std::numeric_limits<double>::max();
93 TCLAP::ValueArg<double> max_elevation_arg(
94 "", "max-elevation",
95 "Maximum elevation for an integrated borehole, "
96 "(min = 0)",
97 false, 0, "MAX_ELEVATION");
98 cmd.add(max_elevation_arg);
99 TCLAP::ValueArg<double> min_elevation_arg(
100 "", "min-elevation",
101 "Minimum elevation for an integrated borehole, "
102 "(min = 0)",
103 false, 0, "MIN_ELEVATION");
104 cmd.add(min_elevation_arg);
105 TCLAP::ValueArg<int> max_id_arg(
106 "", "max-id", "Maximum MaterialID for an integrated borehole", false,
107 -1, "MAX_ID");
108 cmd.add(max_id_arg);
109 TCLAP::ValueArg<int> min_id_arg(
110 "", "min-id", "Minimum MaterialID for an integrated borehole", false,
111 -1, "MIN_ID");
112 cmd.add(min_id_arg);
113 TCLAP::ValueArg<std::string> geo_arg(
114 "g", "geo", "Input (.gml). Name of the geometry file", true, "",
115 "INPUT_FILE");
116 cmd.add(geo_arg);
117 TCLAP::ValueArg<std::string> output_arg(
118 "o", "output", "Output (.vtu). Name of the mesh file", true, "",
119 "OUTPUT_FILE");
120 cmd.add(output_arg);
121 TCLAP::ValueArg<std::string> input_arg(
122 "i", "input", "Input (.vtu). Name of the mesh file", true, "",
123 "INPUT_FILE");
124 cmd.add(input_arg);
125 auto log_level_arg = BaseLib::makeLogLevelArg();
126 cmd.add(log_level_arg);
127 cmd.parse(argc, argv);
128
129 BaseLib::MPI::Setup mpi_setup(argc, argv);
130 BaseLib::initOGSLogger(log_level_arg.getValue());
131
132 std::pair<int, int> mat_limits(0, std::numeric_limits<int>::max());
133 std::pair<double, double> elevation_limits(
134 std::numeric_limits<double>::lowest(), dmax);
135
136 if (min_id_arg.isSet() != max_id_arg.isSet())
137 {
138 ERR("If minimum MaterialID is set, maximum ID must be set, too (and "
139 "vice versa).");
140 return EXIT_FAILURE;
141 }
142 if (min_id_arg.isSet() && max_id_arg.isSet())
143 {
144 mat_limits =
145 std::make_pair(min_id_arg.getValue(), max_id_arg.getValue());
146 }
147 if (mat_limits.first > mat_limits.second)
148 {
149 std::swap(mat_limits.first, mat_limits.second);
150 }
151 if (min_id_arg.isSet() && (mat_limits.first < 0 || mat_limits.second < 0))
152 {
153 ERR("Specified MaterialIDs must have non-negative values.");
154 return EXIT_FAILURE;
155 }
156 if (min_elevation_arg.isSet() != max_elevation_arg.isSet())
157 {
158 ERR("If minimum elevation is set, maximum elevation must be set, too "
159 "(and vice versa).");
160 return EXIT_FAILURE;
161 }
162 if (min_elevation_arg.isSet() && max_elevation_arg.isSet())
163 {
164 elevation_limits = std::make_pair(min_elevation_arg.getValue(),
165 max_elevation_arg.getValue());
166 }
167 if (elevation_limits.first > elevation_limits.second)
168 {
169 std::swap(elevation_limits.first, elevation_limits.second);
170 }
171
172 std::string const& mesh_name = input_arg.getValue();
173 std::string const& output_name = output_arg.getValue();
174 std::string const& geo_name = geo_arg.getValue();
175
178 if (!xml_io.readFile(geo_name))
179 {
180 ERR("Failed to read geometry file `{:s}'.", geo_name);
181 return EXIT_FAILURE;
182 }
183 std::vector<GeoLib::Point*> const& points =
184 *geo.getPointVec(geo.getGeometryNames()[0]);
185
186 std::unique_ptr<MeshLib::Mesh> const mesh(
188 if (mesh == nullptr)
189 {
190 ERR("Failed to read input mesh file `{:s}'.", mesh_name);
191 return EXIT_FAILURE;
192 }
193 if (mesh->getDimension() != 3)
194 {
195 ERR("Method can only be applied to 3D meshes.");
196 return EXIT_FAILURE;
197 }
198
199 auto const& nodes = mesh->getNodes();
200 auto const& mat_ids = MeshLib::materialIDs(*mesh);
201 if (mat_ids == nullptr)
202 {
203 ERR("Mesh is required to have MaterialIDs");
204 return EXIT_FAILURE;
205 }
206
207 auto const& elems = mesh->getElements();
209 auto new_mat_ids = props.createNewPropertyVector<int>(
210 "MaterialIDs", MeshLib::MeshItemType::Cell);
211 std::copy(mat_ids->cbegin(), mat_ids->cend(),
212 std::back_inserter(*new_mat_ids));
213 int const max_id = *std::max_element(mat_ids->begin(), mat_ids->end());
214 std::vector<MeshLib::Node*> new_nodes = MeshLib::copyNodeVector(nodes);
215 std::size_t const n_points = points.size();
216 std::vector<MeshLib::Element*> new_elems =
217 MeshLib::copyElementVector(elems, new_nodes);
218
219 for (std::size_t i = 0; i < n_points; ++i)
220 {
221 std::vector<std::size_t> const& line_nodes = getNodes(
222 *points[i], nodes, *mat_ids, mat_limits, elevation_limits, *mesh);
223 std::size_t const n_line_nodes = line_nodes.size();
224 if (n_line_nodes < 2)
225 {
226 continue;
227 }
228 for (std::size_t j = 0; j < n_line_nodes - 1; ++j)
229 {
230 new_elems.push_back(new MeshLib::Line(
231 {new_nodes[line_nodes[j]], new_nodes[line_nodes[j + 1]]},
232 elems.size()));
233 new_mat_ids->push_back(max_id + i + 1);
234 }
235 }
236
237 MeshLib::Mesh const result("result", new_nodes, new_elems,
238 true /* compute_element_neighbors */, props);
239 MeshLib::IO::VtuInterface vtu(&result);
240 vtu.writeToFile(output_name);
241 return EXIT_SUCCESS;
242}
Definition of the BoostXmlGmlInterface class.
Definition of Duplicate functions.
Definition of the Element class.
Definition of the GEOObjects class.
Git information.
int main(int argc, char *argv[])
std::vector< std::size_t > getNodes(GeoLib::Point const &pnt, std::vector< MeshLib::Node * > const &nodes, MeshLib::PropertyVector< int > const &mat_ids, std::pair< int, int > const &mat_limits, std::pair< double, double > const &elevation_limits, MeshLib::Mesh const &mesh)
Definition of the Line class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:42
Definition of the Mesh class.
Definition of the Node class.
Implementation of the VtuInterface class.
Container class for geometric objects.
Definition GEOObjects.h:57
std::vector< std::string > getGeometryNames() const
Returns the names of all geometry vectors.
const std::vector< Point * > * getPointVec(const std::string &name) const
bool readFile(const std::string &fname) override
Reads an xml-file containing OGS geometry.
std::size_t getID() const
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
bool writeToFile(std::filesystem::path const &file_path)
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:257
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition Properties.h:33
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
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)
std::vector< Node * > copyNodeVector(const std::vector< Node * > &nodes)
Creates a deep copy of a Node vector.
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition Mesh.cpp:269
std::vector< Element * > copyElementVector(std::vector< Element * > const &elements, std::vector< Node * > const &new_nodes, std::vector< std::size_t > const *const node_id_map)
Definition of readMeshFromFile function.