OGS
CreateBoundaryConditionsAlongPolylines.cpp File Reference
#include <tclap/CmdLine.h>
#include <mpi.h>
#include <fstream>
#include <map>
#include <string>
#include <vector>
#include "Applications/FileIO/readGeometryFromFile.h"
#include "Applications/FileIO/writeGeometryToFile.h"
#include "BaseLib/FileTools.h"
#include "GeoLib/GEOObjects.h"
#include "GeoLib/Point.h"
#include "InfoLib/GitInfo.h"
#include "MeshGeoToolsLib/MeshNodeSearcher.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/IO/writeMeshToFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshSurfaceExtraction.h"
Include dependency graph for CreateBoundaryConditionsAlongPolylines.cpp:

Go to the source code of this file.

Functions

void convertMeshNodesToGeometry (std::vector< MeshLib::Node * > const &nodes, std::vector< std::size_t > const &node_ids, std::string &geo_name, GeoLib::GEOObjects &geometry_sets)
 
void writeGroundwaterFlowPointBC (std::ostream &bc_out, std::string const &pnt_name, double head_value)
 
void writeLiquidFlowPointBC (std::ostream &bc_out, std::string const &pnt_name)
 
void writeBCsAndGeometry (GeoLib::GEOObjects &geometry_sets, std::string const &geo_name, std::string const &out_fname, std::string const &bc_type, bool const write_gml)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ convertMeshNodesToGeometry()

void convertMeshNodesToGeometry ( std::vector< MeshLib::Node * > const & nodes,
std::vector< std::size_t > const & node_ids,
std::string & geo_name,
GeoLib::GEOObjects & geometry_sets )

Definition at line 37 of file CreateBoundaryConditionsAlongPolylines.cpp.

41{
42 // copy data
43 std::vector<GeoLib::Point*> pnts;
45 std::size_t cnt(0);
46 for (std::size_t id : node_ids)
47 {
48 pnts.push_back(new GeoLib::Point(*(nodes[id]), cnt));
49 pnt_names[geo_name + "-PNT-" + std::to_string(cnt)] = cnt;
50 cnt++;
51 }
52
53 // create data structures for geometry
54 geometry_sets.addPointVec(std::move(pnts), geo_name, std::move(pnt_names));
55}
void addPointVec(std::vector< Point * > &&points, std::string &name, PointVec::NameIdMap &&pnt_id_name_map, double const eps=std::sqrt(std::numeric_limits< double >::epsilon()))
std::map< std::string, std::size_t > NameIdMap
Definition TemplateVec.h:41

References GeoLib::GEOObjects::addPointVec().

Referenced by main().

◆ main()

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

Definition at line 129 of file CreateBoundaryConditionsAlongPolylines.cpp.

130{
131 TCLAP::CmdLine cmd(
132 "Creates boundary conditions for mesh nodes along polylines."
133 "The documentation is available at "
134 "https://docs.opengeosys.org/docs/tools/model-preparation/"
135 "create-boundary-conditions-along-a-polyline.\n\n"
136 "OpenGeoSys-6 software, version " +
138 ".\n"
139 "Copyright (c) 2012-2024, OpenGeoSys Community "
140 "(http://www.opengeosys.org)",
142 TCLAP::SwitchArg gml_arg("", "gml", "Write found nodes to gml file.");
143 cmd.add(gml_arg);
144
145 TCLAP::ValueArg<std::string> output_base_fname(
146 "o", "output-base-file-name",
147 "the base name of the file the output (geometry (gli) and boundary "
148 "condition (bc)) will be written to",
149 true, "", "file name");
150 cmd.add(output_base_fname);
151
152 TCLAP::ValueArg<std::string> bc_type(
153 "t", "type",
154 "the process type the boundary condition will be written for currently "
155 "LIQUID_FLOW (primary variable PRESSURE1) and GROUNDWATER_FLOW "
156 "(primary variable HEAD, default) are supported",
157 true, "",
158 "process type as string (LIQUID_FLOW or GROUNDWATER_FLOW (default))");
159 cmd.add(bc_type);
160
161 TCLAP::ValueArg<double> search_length_arg(
162 "s", "search-length",
163 "The size of the search length. The default value is "
164 "std::numeric_limits<double>::epsilon()",
165 false, std::numeric_limits<double>::epsilon(), "floating point number");
166 cmd.add(search_length_arg);
167
168 TCLAP::ValueArg<std::string> geometry_fname(
169 "i", "input-geometry",
170 "the name of the file containing the input geometry", true, "",
171 "file name");
172 cmd.add(geometry_fname);
173
174 TCLAP::ValueArg<std::string> mesh_arg(
175 "m", "mesh-file", "the name of the file containing the mesh", true, "",
176 "file name");
177 cmd.add(mesh_arg);
178
179 TCLAP::ValueArg<std::string> gmsh_path_arg("g", "gmsh-path",
180 "the path to the gmsh binary",
181 false, "", "path as string");
182 cmd.add(gmsh_path_arg);
183
184 cmd.parse(argc, argv);
185
186#ifdef USE_PETSC
187 MPI_Init(&argc, &argv);
188#endif
189
190 // *** read mesh
191 INFO("Reading mesh '{:s}' ... ", mesh_arg.getValue());
192 std::unique_ptr<MeshLib::Mesh> subsurface_mesh(
193 MeshLib::IO::readMeshFromFile(mesh_arg.getValue()));
194 INFO("done.");
195 INFO("Extracting top surface of mesh '{:s}' ... ", mesh_arg.getValue());
196 Eigen::Vector3d const dir({0, 0, -1});
197 double const angle(90);
198 std::unique_ptr<MeshLib::Mesh> surface_mesh(
200 dir, angle));
201 INFO("done.");
202 subsurface_mesh.reset(nullptr);
203
204 // *** read geometry
205 GeoLib::GEOObjects geometries;
206 FileIO::readGeometryFromFile(geometry_fname.getValue(), geometries,
207 gmsh_path_arg.getValue());
208
209 auto const geo_name = geometries.getGeometryNames()[0];
210
211 // *** check if the data is usable
212 // *** get vector of polylines
213 std::vector<GeoLib::Polyline*> const* plys(
214 geometries.getPolylineVec(geo_name));
215 if (!plys)
216 {
217 ERR("Could not get vector of polylines out of geometry '{:s}'.",
218 geo_name);
219#ifdef USE_PETSC
220 MPI_Finalize();
221#endif
222 return EXIT_FAILURE;
223 }
224
225 auto search_length_strategy =
226 std::make_unique<MeshGeoToolsLib::SearchLength>();
227 if (search_length_arg.isSet())
228 {
229 search_length_strategy.reset(
230 new MeshGeoToolsLib::SearchLength(search_length_arg.getValue()));
231 }
232
233 GeoLib::GEOObjects geometry_sets;
235 *surface_mesh, std::move(search_length_strategy),
237 for (std::size_t k(0); k < plys->size(); k++)
238 {
239 auto const& ids = mesh_searcher.getMeshNodeIDs(*((*plys)[k]));
240 if (ids.empty())
241 {
242 continue;
243 }
244 std::string polyline_name("Polyline-" + std::to_string(k));
245 convertMeshNodesToGeometry(surface_mesh->getNodes(), ids, polyline_name,
246 geometry_sets);
247 }
248
249 // merge all together
250 auto const geo_names = geometry_sets.getGeometryNames();
251 if (geo_names.empty())
252 {
253 ERR("Did not find mesh nodes along polylines.");
254#ifdef USE_PETSC
255 MPI_Finalize();
256#endif
257 return EXIT_FAILURE;
258 }
259
260 std::string merge_name("AllMeshNodesAlongPolylines");
261 if (geometry_sets.mergeGeometries(geo_names, merge_name) == 2)
262 {
263 merge_name = geo_names[0];
264 }
265
266 GeoLib::PointVec const* pnt_vec(geometry_sets.getPointVecObj(merge_name));
267 auto const& merged_pnts(pnt_vec->getVector());
268
269 std::vector<GeoLib::Point> pnts_with_id;
270 const std::size_t n_merged_pnts(merged_pnts.size());
271 for (std::size_t k(0); k < n_merged_pnts; ++k)
272 {
273 pnts_with_id.emplace_back(*(merged_pnts[k]), k);
274 }
275
276 std::sort(pnts_with_id.begin(), pnts_with_id.end(),
277 [](GeoLib::Point const& p0, GeoLib::Point const& p1)
278 { return p0 < p1; });
279
280 double const eps(std::numeric_limits<double>::epsilon());
281 std::vector<GeoLib::Point*> surface_pnts;
282 GeoLib::PointVec::NameIdMap name_id_map;
283
284 // insert first point
285 surface_pnts.push_back(
286 new GeoLib::Point(pnts_with_id[0], surface_pnts.size()));
287 {
288 std::string element_name;
289 pnt_vec->getNameOfElementByID(0, element_name);
290 name_id_map[element_name] = 0;
291 }
292 for (std::size_t k(1); k < n_merged_pnts; ++k)
293 {
294 const GeoLib::Point& p0(pnts_with_id[k - 1]);
295 const GeoLib::Point& p1(pnts_with_id[k]);
296 if (std::abs(p0[0] - p1[0]) > eps || std::abs(p0[1] - p1[1]) > eps)
297 {
298 surface_pnts.push_back(
299 new GeoLib::Point(pnts_with_id[k], surface_pnts.size()));
300 std::string element_name;
301 pnt_vec->getNameOfElementByID(k, element_name);
302 name_id_map[element_name] = surface_pnts.size() - 1;
303 }
304 }
305
306 std::string surface_name(BaseLib::dropFileExtension(mesh_arg.getValue()) +
307 "-MeshNodesAlongPolylines");
308 geometry_sets.addPointVec(std::move(surface_pnts), surface_name,
309 std::move(name_id_map), 1e-6);
310
311 // write the BCs and the merged geometry set to file
312 std::string const base_fname(
313 BaseLib::dropFileExtension(output_base_fname.getValue()));
314 writeBCsAndGeometry(geometry_sets, surface_name, base_fname,
315 bc_type.getValue(), gml_arg.getValue());
316#ifdef USE_PETSC
317 MPI_Finalize();
318#endif
319 return EXIT_SUCCESS;
320}
void writeBCsAndGeometry(GeoLib::GEOObjects &geometry_sets, std::string const &geo_name, std::string const &out_fname, std::string const &bc_type, bool const write_gml)
void convertMeshNodesToGeometry(std::vector< MeshLib::Node * > const &nodes, std::vector< std::size_t > const &node_ids, std::string &geo_name, GeoLib::GEOObjects &geometry_sets)
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
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 PointVec * getPointVecObj(const std::string &name) const
int mergeGeometries(std::vector< std::string > const &geo_names, std::string &merged_geo_name)
const std::vector< Polyline * > * getPolylineVec(const std::string &name) const
This class manages pointers to Points in a std::vector along with a name. It also handles the deletio...
Definition PointVec.h:36
static MeshLib::Mesh * getMeshSurface(const MeshLib::Mesh &subsfc_mesh, Eigen::Vector3d const &dir, double angle, std::string_view subsfc_node_id_prop_name="", std::string_view subsfc_element_id_prop_name="", std::string_view face_id_prop_name="")
std::string dropFileExtension(std::string const &filename)
void readGeometryFromFile(std::string const &fname, GeoLib::GEOObjects &geo_objs, std::string const &gmsh_path)
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:225

References GeoLib::GEOObjects::addPointVec(), convertMeshNodesToGeometry(), BaseLib::dropFileExtension(), ERR(), GeoLib::GEOObjects::getGeometryNames(), MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeIDs(), MeshToolsLib::MeshSurfaceExtraction::getMeshSurface(), GeoLib::TemplateVec< T >::getNameOfElementByID(), GeoLib::GEOObjects::getPointVecObj(), GeoLib::GEOObjects::getPolylineVec(), GeoLib::TemplateVec< T >::getVector(), INFO(), GeoLib::GEOObjects::mergeGeometries(), GitInfoLib::GitInfo::ogs_version, FileIO::readGeometryFromFile(), MeshLib::IO::readMeshFromFile(), writeBCsAndGeometry(), and MeshGeoToolsLib::Yes.

◆ writeBCsAndGeometry()

void writeBCsAndGeometry ( GeoLib::GEOObjects & geometry_sets,
std::string const & geo_name,
std::string const & out_fname,
std::string const & bc_type,
bool const write_gml )

Definition at line 87 of file CreateBoundaryConditionsAlongPolylines.cpp.

91{
92 if (write_gml)
93 {
94 INFO("write points to '{:s}.gml'.", geo_name);
95 FileIO::writeGeometryToFile(geo_name, geometry_sets,
96 out_fname + ".gml");
97 }
98 FileIO::writeGeometryToFile(geo_name, geometry_sets, out_fname + ".gli");
99
100 bool liquid_flow(false);
101 if (bc_type == "LIQUID_FLOW")
102 {
103 liquid_flow = true;
104 }
105
106 GeoLib::PointVec const* pnt_vec_objs(
107 geometry_sets.getPointVecObj(geo_name));
108 auto const& pnts(pnt_vec_objs->getVector());
109 std::ofstream bc_out(out_fname + ".bc");
110 for (std::size_t k(0); k < pnts.size(); k++)
111 {
112 std::string const& pnt_name(pnt_vec_objs->getItemNameByID(k));
113 if (!pnt_name.empty())
114 {
115 if (liquid_flow)
116 {
117 writeLiquidFlowPointBC(bc_out, pnt_name);
118 }
119 else
120 {
121 writeGroundwaterFlowPointBC(bc_out, pnt_name, (*pnts[k])[2]);
122 }
123 }
124 }
125 bc_out << "#STOP\n";
126 bc_out.close();
127}
void writeGroundwaterFlowPointBC(std::ostream &bc_out, std::string const &pnt_name, double head_value)
void writeLiquidFlowPointBC(std::ostream &bc_out, std::string const &pnt_name)
void writeGeometryToFile(std::string const &geo_name, GeoLib::GEOObjects &geo_objs, std::string const &fname)

References GeoLib::PointVec::getItemNameByID(), GeoLib::GEOObjects::getPointVecObj(), GeoLib::TemplateVec< T >::getVector(), INFO(), FileIO::writeGeometryToFile(), writeGroundwaterFlowPointBC(), and writeLiquidFlowPointBC().

Referenced by main().

◆ writeGroundwaterFlowPointBC()

void writeGroundwaterFlowPointBC ( std::ostream & bc_out,
std::string const & pnt_name,
double head_value )

Definition at line 57 of file CreateBoundaryConditionsAlongPolylines.cpp.

59{
60 bc_out << "#BOUNDARY_CONDITION\n";
61 bc_out << " $PCS_TYPE\n";
62 bc_out << " GROUNDWATER_FLOW\n";
63 bc_out << " $PRIMARY_VARIABLE\n";
64 bc_out << " HEAD\n";
65 bc_out << " $GEO_TYPE\n";
66 bc_out << " POINT " << pnt_name << "\n";
67 bc_out << " $DIS_TYPE\n";
68 bc_out << " CONSTANT " << head_value << "\n";
69}

Referenced by writeBCsAndGeometry().

◆ writeLiquidFlowPointBC()

void writeLiquidFlowPointBC ( std::ostream & bc_out,
std::string const & pnt_name )

Definition at line 71 of file CreateBoundaryConditionsAlongPolylines.cpp.

72{
73 bc_out << "#BOUNDARY_CONDITION\n";
74 bc_out << " $PCS_TYPE\n";
75 bc_out << " LIQUID_FLOW\n";
76 bc_out << " $PRIMARY_VARIABLE\n";
77 bc_out << " PRESSURE1\n";
78 bc_out << " $GEO_TYPE\n";
79 bc_out << " POINT " << pnt_name << "\n";
80 bc_out << " $DIS_TYPE\n";
81 bc_out << " CONSTANT 0.0\n";
82}

Referenced by writeBCsAndGeometry().