OGS
IntegrateBoreholesIntoMesh.cpp File Reference

Detailed Description

Definition in file IntegrateBoreholesIntoMesh.cpp.

#include <iterator>
#include <limits>
#include <memory>
#include <string>
#include <vector>
#include <tclap/CmdLine.h>
#include <mpi.h>
#include "GeoLib/GEOObjects.h"
#include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
#include "InfoLib/GitInfo.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Elements/Line.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/Node.h"
#include "MeshLib/Utils/DuplicateMeshComponents.h"
Include dependency graph for IntegrateBoreholesIntoMesh.cpp:

Go to the source code of this file.

Functions

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)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ getNodes()

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 at line 34 of file IntegrateBoreholesIntoMesh.cpp.

40{
41 std::vector<std::size_t> pnt_nodes;
42 for (auto node : nodes)
43 {
44 double const eps = std::numeric_limits<double>::epsilon();
45 if (std::abs((*node)[0] - pnt[0]) < eps &&
46 std::abs((*node)[1] - pnt[1]) < eps)
47 {
48 auto const& elems = mesh.getElementsConnectedToNode(*node);
49 for (auto e : elems)
50 {
51 if (mat_ids[e->getID()] >= mat_limits.first &&
52 mat_ids[e->getID()] <= mat_limits.second &&
53 (*node)[2] >= elevation_limits.first &&
54 (*node)[2] <= elevation_limits.second)
55 {
56 pnt_nodes.push_back(node->getID());
57 break;
58 }
59 }
60 }
61 }
62 if (pnt_nodes.size() < 2)
63 {
64 WARN("No nodes found for point {:d}...", pnt.getID());
65 }
66 else
67 {
68 // sort found nodes from top to bottom (required for BHE simulations)
69 std::sort(pnt_nodes.begin(), pnt_nodes.end(),
70 [nodes](std::size_t a, std::size_t b)
71 { return (*nodes[a])[2] > (*nodes[b])[2]; });
72 }
73 return pnt_nodes;
74}
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40

References MeshLib::Mesh::getElementsConnectedToNode(), MathLib::Point3dWithID::getID(), and WARN().

Referenced by fixVtkInconsistencies(), main(), reverseNodeOrder(), and MeshLib::setMeshSpaceDimension().

◆ main()

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

Definition at line 76 of file IntegrateBoreholesIntoMesh.cpp.

77{
78 TCLAP::CmdLine cmd(
79 "Integrates line elements representing boreholes into a pre-existing "
80 "3D mesh. Corresponding nodes matching the (x,y)-coordinates given in "
81 "the gml-file are found in the mesh and connected from top to bottom "
82 "via line elements. Each borehole (i.e. all points at a given "
83 "(x,y)-location but at different depths) is assigned a unique material "
84 "ID. Vertical limits of boreholes can be specified via Material IDs "
85 "and/or elevation. Point not matching any mesh nodes or located outside"
86 " the mesh are ignored.\n\n"
87 "OpenGeoSys-6 software, version " +
89 ".\n"
90 "Copyright (c) 2012-2024, OpenGeoSys Community "
91 "(http://www.opengeosys.org)",
93
94 double const dmax = std::numeric_limits<double>::max();
95 TCLAP::ValueArg<double> max_elevation_arg(
96 "", "max-elevation", "Maximum elevation for an integrated borehole",
97 false, 0, "a number");
98 cmd.add(max_elevation_arg);
99 TCLAP::ValueArg<double> min_elevation_arg(
100 "", "min-elevation", "Minimum elevation for an integrated borehole",
101 false, 0, "a number");
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, "a number");
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, "a number");
110 cmd.add(min_id_arg);
111 TCLAP::ValueArg<std::string> geo_arg("g", "geo",
112 "Name of the geometry file (*.gml)",
113 true, "", "geometry file name");
114 cmd.add(geo_arg);
115 TCLAP::ValueArg<std::string> output_arg("o", "output",
116 "Name of the output mesh (*.vtu)",
117 true, "", "output file name");
118 cmd.add(output_arg);
119 TCLAP::ValueArg<std::string> input_arg("i", "input",
120 "Name of the input mesh (*.vtu)",
121 true, "", "input file name");
122 cmd.add(input_arg);
123 cmd.parse(argc, argv);
124
125#ifdef USE_PETSC
126 MPI_Init(&argc, &argv);
127#endif
128
129 std::pair<int, int> mat_limits(0, std::numeric_limits<int>::max());
130 std::pair<double, double> elevation_limits(
131 std::numeric_limits<double>::lowest(), dmax);
132
133 if (min_id_arg.isSet() != max_id_arg.isSet())
134 {
135 ERR("If minimum MaterialID is set, maximum ID must be set, too (and "
136 "vice versa).");
137#ifdef USE_PETSC
138 MPI_Finalize();
139#endif
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#ifdef USE_PETSC
155 MPI_Finalize();
156#endif
157 return EXIT_FAILURE;
158 }
159 if (min_elevation_arg.isSet() != max_elevation_arg.isSet())
160 {
161 ERR("If minimum elevation is set, maximum elevation must be set, too "
162 "(and vice versa).");
163#ifdef USE_PETSC
164 MPI_Finalize();
165#endif
166 return EXIT_FAILURE;
167 }
168 if (min_elevation_arg.isSet() && max_elevation_arg.isSet())
169 {
170 elevation_limits = std::make_pair(min_elevation_arg.getValue(),
171 max_elevation_arg.getValue());
172 }
173 if (elevation_limits.first > elevation_limits.second)
174 {
175 std::swap(elevation_limits.first, elevation_limits.second);
176 }
177
178 std::string const& mesh_name = input_arg.getValue();
179 std::string const& output_name = output_arg.getValue();
180 std::string const& geo_name = geo_arg.getValue();
181
184 if (!xml_io.readFile(geo_name))
185 {
186 ERR("Failed to read geometry file `{:s}'.", geo_name);
187#ifdef USE_PETSC
188 MPI_Finalize();
189#endif
190 return EXIT_FAILURE;
191 }
192 std::vector<GeoLib::Point*> const& points =
193 *geo.getPointVec(geo.getGeometryNames()[0]);
194
195 std::unique_ptr<MeshLib::Mesh> const mesh(
197 if (mesh == nullptr)
198 {
199 ERR("Failed to read input mesh file `{:s}'.", mesh_name);
200#ifdef USE_PETSC
201 MPI_Finalize();
202#endif
203 return EXIT_FAILURE;
204 }
205 if (mesh->getDimension() != 3)
206 {
207 ERR("Method can only be applied to 3D meshes.");
208#ifdef USE_PETSC
209 MPI_Finalize();
210#endif
211 return EXIT_FAILURE;
212 }
213
214 auto const& nodes = mesh->getNodes();
215 auto const& mat_ids = MeshLib::materialIDs(*mesh);
216 if (mat_ids == nullptr)
217 {
218 ERR("Mesh is required to have MaterialIDs");
219#ifdef USE_PETSC
220 MPI_Finalize();
221#endif
222 return EXIT_FAILURE;
223 }
224
225 auto const& elems = mesh->getElements();
227 auto new_mat_ids = props.createNewPropertyVector<int>(
228 "MaterialIDs", MeshLib::MeshItemType::Cell);
229 std::copy(mat_ids->cbegin(), mat_ids->cend(),
230 std::back_inserter(*new_mat_ids));
231 int const max_id = *std::max_element(mat_ids->begin(), mat_ids->end());
232 std::vector<MeshLib::Node*> new_nodes = MeshLib::copyNodeVector(nodes);
233 std::size_t const n_points = points.size();
234 std::vector<MeshLib::Element*> new_elems =
235 MeshLib::copyElementVector(elems, new_nodes);
236
237 for (std::size_t i = 0; i < n_points; ++i)
238 {
239 std::vector<std::size_t> const& line_nodes = getNodes(
240 *points[i], nodes, *mat_ids, mat_limits, elevation_limits, *mesh);
241 std::size_t const n_line_nodes = line_nodes.size();
242 if (n_line_nodes < 2)
243 {
244 continue;
245 }
246 for (std::size_t j = 0; j < n_line_nodes - 1; ++j)
247 {
248 new_elems.push_back(new MeshLib::Line(
249 {new_nodes[line_nodes[j]], new_nodes[line_nodes[j + 1]]},
250 elems.size()));
251 new_mat_ids->push_back(max_id + i + 1);
252 }
253 }
254
255 MeshLib::Mesh const result("result", new_nodes, new_elems,
256 true /* compute_element_neighbors */, props);
257 MeshLib::IO::VtuInterface vtu(&result);
258 vtu.writeToFile(output_name);
259#ifdef USE_PETSC
260 MPI_Finalize();
261#endif
262 return EXIT_SUCCESS;
263}
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)
Definition Logging.h:45
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
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....
Definition Properties.h:36
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)
Definition Mesh.cpp:268
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)

References MeshLib::Cell, MeshLib::copyElementVector(), MeshLib::copyNodeVector(), MeshLib::Properties::createNewPropertyVector(), ERR(), GeoLib::GEOObjects::getGeometryNames(), getNodes(), GeoLib::GEOObjects::getPointVec(), MeshLib::materialIDs(), GitInfoLib::GitInfo::ogs_version, GeoLib::IO::BoostXmlGmlInterface::readFile(), MeshLib::IO::readMeshFromFile(), and MeshLib::IO::VtuInterface::writeToFile().