73{
74 TCLAP::CmdLine cmd(
75 "Integrates line elements representing boreholes into a pre-existing "
76 "3D mesh. Corresponding nodes matching the (x,y)-coordinates given in "
77 "the gml-file are found in the mesh and connected from top to bottom "
78 "via line elements. Each borehole (i.e. all points at a given "
79 "(x,y)-location but at different depths) is assigned a unique material "
80 "ID. Vertical limits of boreholes can be specified via Material IDs "
81 "and/or elevation. Point not matching any mesh nodes or located outside"
82 " the mesh are ignored.\n\n"
83 "OpenGeoSys-6 software, version " +
85 ".\n"
86 "Copyright (c) 2012-2025, OpenGeoSys Community "
87 "(http://www.opengeosys.org)",
89
90 double const dmax = std::numeric_limits<double>::max();
91 TCLAP::ValueArg<double> max_elevation_arg(
92 "", "max-elevation",
93 "Maximum elevation for an integrated borehole, "
94 "(min = 0)",
95 false, 0, "MAX_ELEVATION");
96 cmd.add(max_elevation_arg);
97 TCLAP::ValueArg<double> min_elevation_arg(
98 "", "min-elevation",
99 "Minimum elevation for an integrated borehole, "
100 "(min = 0)",
101 false, 0, "MIN_ELEVATION");
102 cmd.add(min_elevation_arg);
103 TCLAP::ValueArg<int> max_id_arg(
104 "", "max-id", "Maximum MaterialID for an integrated borehole", false,
105 -1, "MAX_ID");
106 cmd.add(max_id_arg);
107 TCLAP::ValueArg<int> min_id_arg(
108 "", "min-id", "Minimum MaterialID for an integrated borehole", false,
109 -1, "MIN_ID");
110 cmd.add(min_id_arg);
111 TCLAP::ValueArg<std::string> geo_arg(
112 "g", "geo", "Input (.gml). Name of the geometry file", true, "",
113 "INPUT_FILE");
114 cmd.add(geo_arg);
115 TCLAP::ValueArg<std::string> output_arg(
116 "o", "output", "Output (.vtu). Name of the mesh file", true, "",
117 "OUTPUT_FILE");
118 cmd.add(output_arg);
119 TCLAP::ValueArg<std::string> input_arg(
120 "i", "input", "Input (.vtu). Name of the mesh file", true, "",
121 "INPUT_FILE");
122 cmd.add(input_arg);
123 cmd.parse(argc, argv);
124
126
127 std::pair<int, int> mat_limits(0, std::numeric_limits<int>::max());
128 std::pair<double, double> elevation_limits(
129 std::numeric_limits<double>::lowest(), dmax);
130
131 if (min_id_arg.isSet() != max_id_arg.isSet())
132 {
133 ERR(
"If minimum MaterialID is set, maximum ID must be set, too (and "
134 "vice versa).");
135 return EXIT_FAILURE;
136 }
137 if (min_id_arg.isSet() && max_id_arg.isSet())
138 {
139 mat_limits =
140 std::make_pair(min_id_arg.getValue(), max_id_arg.getValue());
141 }
142 if (mat_limits.first > mat_limits.second)
143 {
144 std::swap(mat_limits.first, mat_limits.second);
145 }
146 if (min_id_arg.isSet() && (mat_limits.first < 0 || mat_limits.second < 0))
147 {
148 ERR(
"Specified MaterialIDs must have non-negative values.");
149 return EXIT_FAILURE;
150 }
151 if (min_elevation_arg.isSet() != max_elevation_arg.isSet())
152 {
153 ERR(
"If minimum elevation is set, maximum elevation must be set, too "
154 "(and vice versa).");
155 return EXIT_FAILURE;
156 }
157 if (min_elevation_arg.isSet() && max_elevation_arg.isSet())
158 {
159 elevation_limits = std::make_pair(min_elevation_arg.getValue(),
160 max_elevation_arg.getValue());
161 }
162 if (elevation_limits.first > elevation_limits.second)
163 {
164 std::swap(elevation_limits.first, elevation_limits.second);
165 }
166
167 std::string const& mesh_name = input_arg.getValue();
168 std::string const& output_name = output_arg.getValue();
169 std::string const& geo_name = geo_arg.getValue();
170
173 if (!xml_io.readFile(geo_name))
174 {
175 ERR(
"Failed to read geometry file `{:s}'.", geo_name);
176 return EXIT_FAILURE;
177 }
178 std::vector<GeoLib::Point*> const& points =
180
181 std::unique_ptr<MeshLib::Mesh> const mesh(
183 if (mesh == nullptr)
184 {
185 ERR(
"Failed to read input mesh file `{:s}'.", mesh_name);
186 return EXIT_FAILURE;
187 }
188 if (mesh->getDimension() != 3)
189 {
190 ERR(
"Method can only be applied to 3D meshes.");
191 return EXIT_FAILURE;
192 }
193
194 auto const& nodes = mesh->getNodes();
196 if (mat_ids == nullptr)
197 {
198 ERR(
"Mesh is required to have MaterialIDs");
199 return EXIT_FAILURE;
200 }
201
202 auto const& elems = mesh->getElements();
206 std::copy(mat_ids->cbegin(), mat_ids->cend(),
207 std::back_inserter(*new_mat_ids));
208 int const max_id = *std::max_element(mat_ids->begin(), mat_ids->end());
210 std::size_t const n_points = points.size();
211 std::vector<MeshLib::Element*> new_elems =
213
214 for (std::size_t i = 0; i < n_points; ++i)
215 {
216 std::vector<std::size_t>
const& line_nodes =
getNodes(
217 *points[i], nodes, *mat_ids, mat_limits, elevation_limits, *mesh);
218 std::size_t const n_line_nodes = line_nodes.size();
219 if (n_line_nodes < 2)
220 {
221 continue;
222 }
223 for (std::size_t j = 0; j < n_line_nodes - 1; ++j)
224 {
226 {new_nodes[line_nodes[j]], new_nodes[line_nodes[j + 1]]},
227 elems.size()));
228 new_mat_ids->push_back(max_id + i + 1);
229 }
230 }
231
233 true , props);
235 vtu.writeToFile(output_name);
236 return EXIT_SUCCESS;
237}
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)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Container class for geometric objects.
std::vector< std::string > getGeometryNames() const
Returns the names of all geometry vectors.
const std::vector< Point * > * getPointVec(const std::string &name) const
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
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)
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)