17 #include <tclap/CmdLine.h>
31 GeoLib::Point const& pnt, std::vector<MeshLib::Node*>
const& nodes,
33 std::pair<int, int>
const& mat_limits,
34 std::pair<double, double>
const& elevation_limits,
37 std::vector<std::size_t> pnt_nodes;
38 for (
auto node : nodes)
40 double const eps = std::numeric_limits<double>::epsilon();
41 if (std::abs((*node)[0] - pnt[0]) < eps &&
42 std::abs((*node)[1] - pnt[1]) < eps)
47 if (mat_ids[e->getID()] >= mat_limits.first &&
48 mat_ids[e->getID()] <= mat_limits.second &&
49 (*node)[2] >= elevation_limits.first &&
50 (*node)[2] <= elevation_limits.second)
52 pnt_nodes.push_back(node->getID());
58 if (pnt_nodes.size() < 2)
60 WARN(
"No nodes found for point {:d}...", pnt.
getID());
65 std::sort(pnt_nodes.begin(), pnt_nodes.end(),
66 [nodes](std::size_t a, std::size_t b)
67 { return (*nodes[a])[2] > (*nodes[b])[2]; });
72 int main(
int argc,
char* argv[])
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 " +
86 "Copyright (c) 2012-2021, OpenGeoSys Community "
87 "(http://www.opengeosys.org)",
90 double const dmax = std::numeric_limits<double>::max();
91 TCLAP::ValueArg<double> max_elevation_arg(
92 "",
"max-elevation",
"Maximum elevation for an integrated borehole",
93 false, 0,
"a number");
94 cmd.add(max_elevation_arg);
95 TCLAP::ValueArg<double> min_elevation_arg(
96 "",
"min-elevation",
"Minimum elevation for an integrated borehole",
97 false, 0,
"a number");
98 cmd.add(min_elevation_arg);
99 TCLAP::ValueArg<int> max_id_arg(
100 "",
"max-id",
"Maximum MaterialID for an integrated borehole",
false,
103 TCLAP::ValueArg<int> min_id_arg(
104 "",
"min-id",
"Minimum MaterialID for an integrated borehole",
false,
107 TCLAP::ValueArg<std::string> geo_arg(
"g",
"geo",
108 "Name of the geometry file (*.gml)",
109 true,
"",
"geometry file name");
111 TCLAP::ValueArg<std::string> output_arg(
"o",
"output",
112 "Name of the output mesh (*.vtu)",
113 true,
"",
"output file name");
115 TCLAP::ValueArg<std::string> input_arg(
"i",
"input",
116 "Name of the input mesh (*.vtu)",
117 true,
"",
"input file name");
119 cmd.parse(argc, argv);
121 std::pair<int, int> mat_limits(0, std::numeric_limits<int>::max());
122 std::pair<double, double> elevation_limits(
123 std::numeric_limits<double>::lowest(), dmax);
125 if (min_id_arg.isSet() != max_id_arg.isSet())
127 ERR(
"If minimum MaterialID is set, maximum ID must be set, too (and "
131 if (min_id_arg.isSet() && max_id_arg.isSet())
134 std::make_pair(min_id_arg.getValue(), max_id_arg.getValue());
136 if (mat_limits.first > mat_limits.second)
138 std::swap(mat_limits.first, mat_limits.second);
140 if (min_id_arg.isSet() && (mat_limits.first < 0 || mat_limits.second < 0))
142 ERR(
"Specified MaterialIDs must have non-negative values.");
145 if (min_elevation_arg.isSet() != max_elevation_arg.isSet())
147 ERR(
"If minimum elevation is set, maximum elevation must be set, too "
148 "(and vice versa).");
151 if (min_elevation_arg.isSet() && max_elevation_arg.isSet())
153 elevation_limits = std::make_pair(min_elevation_arg.getValue(),
154 max_elevation_arg.getValue());
156 if (elevation_limits.first > elevation_limits.second)
158 std::swap(elevation_limits.first, elevation_limits.second);
161 std::string
const& mesh_name = input_arg.getValue();
162 std::string
const& output_name = output_arg.getValue();
163 std::string
const& geo_name = geo_arg.getValue();
169 ERR(
"Failed to read geometry file `{:s}'.", geo_name);
172 std::vector<GeoLib::Point*>
const& points =
175 std::unique_ptr<MeshLib::Mesh>
const mesh(
179 ERR(
"Failed to read input mesh file `{:s}'.", mesh_name);
182 if (mesh->getDimension() != 3)
184 ERR(
"Method can only be applied to 3D meshes.");
188 auto const& nodes = mesh->getNodes();
190 if (mat_ids ==
nullptr)
192 ERR(
"Mesh is required to have MaterialIDs");
196 auto const& elems = mesh->getElements();
200 std::copy(mat_ids->cbegin(), mat_ids->cend(),
201 std::back_inserter(*new_mat_ids));
202 int const max_id = *std::max_element(mat_ids->begin(), mat_ids->end());
204 std::size_t
const n_points = points.size();
205 std::vector<MeshLib::Element*> new_elems =
208 for (std::size_t i = 0; i < n_points; ++i)
210 std::vector<std::size_t>
const& line_nodes =
getNodes(
211 *points[i], nodes, *mat_ids, mat_limits, elevation_limits, *mesh);
212 std::size_t
const n_line_nodes = line_nodes.size();
213 if (n_line_nodes < 2)
217 for (std::size_t j = 0; j < n_line_nodes - 1; ++j)
220 {new_nodes[line_nodes[j]], new_nodes[line_nodes[j + 1]]},
222 new_mat_ids->push_back(max_id + i + 1);
226 MeshLib::Mesh const result(
"result", new_nodes, new_elems, props);
Definition of the BoostXmlGmlInterface class.
Definition of Duplicate functions.
Definition of the Element class.
Definition of the GEOObjects class.
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(char const *fmt, Args const &... args)
void WARN(char const *fmt, Args const &... args)
Definition of the Mesh class.
Definition of the Node class.
Implementation of the VtuInterface class.
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
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
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
PropertyVector< T > * createNewPropertyVector(std::string const &name, MeshItemType mesh_item_type, std::size_t n_components=1)
GITINFOLIB_EXPORT const std::string ogs_version
void copy(PETScVector const &x, PETScVector &y)
MeshLib::Mesh * readMeshFromFile(const std::string &file_name)
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)
Definition of readMeshFromFile function.