OGS
CreateBoundaryConditionsAlongPolylines.cpp File Reference
#include <tclap/CmdLine.h>
#include <fstream>
#include <map>
#include <string>
#include <vector>
#include "Applications/FileIO/readGeometryFromFile.h"
#include "Applications/FileIO/writeGeometryToFile.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/MPI.h"
#include "BaseLib/TCLAPArguments.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 36 of file CreateBoundaryConditionsAlongPolylines.cpp.

40{
41 // copy data
42 std::vector<GeoLib::Point*> pnts;
44 std::size_t cnt(0);
45 for (std::size_t id : node_ids)
46 {
47 pnts.push_back(new GeoLib::Point(*(nodes[id]), cnt));
48 pnt_names[geo_name + "-PNT-" + std::to_string(cnt)] = cnt;
49 cnt++;
50 }
51
52 // create data structures for geometry
53 geometry_sets.addPointVec(std::move(pnts), geo_name, std::move(pnt_names));
54}
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 128 of file CreateBoundaryConditionsAlongPolylines.cpp.

129{
130 TCLAP::CmdLine cmd(
131 "Creates boundary conditions for mesh nodes along polylines."
132 "The documentation is available at "
133 "https://docs.opengeosys.org/docs/tools/model-preparation/"
134 "create-boundary-conditions-along-a-polyline.\n\n"
135 "OpenGeoSys-6 software, version " +
137 ".\n"
138 "Copyright (c) 2012-2025, OpenGeoSys Community "
139 "(http://www.opengeosys.org)",
141 TCLAP::SwitchArg gml_arg("", "gml", "Write found nodes to gml file.");
142 cmd.add(gml_arg);
143
144 TCLAP::ValueArg<std::string> output_base_fname(
145 "o", "output-base-file-name",
146 "Output. The base name of the file the output (geometry (gli) and "
147 "boundary "
148 "condition (bc)) will be written to",
149 true, "", "BASE_FILENAME_OUTPUT");
150 cmd.add(output_base_fname);
151
152 std::vector<std::string> allowed_types_vector{"LIQUID_FLOW",
153 "GROUNDWATER_FLOW"};
154 TCLAP::ValuesConstraint<std::string> allowed_types(allowed_types_vector);
155 TCLAP::ValueArg<std::string> bc_type(
156 "t", "type",
157 "the process type the boundary condition will be written for currently "
158 "LIQUID_FLOW (primary variable PRESSURE1) and GROUNDWATER_FLOW "
159 "(primary variable HEAD, default) are supported, ",
160 true, "", &allowed_types);
161 cmd.add(bc_type);
162
163 TCLAP::ValueArg<double> search_length_arg(
164 "s", "search-length", "The size of the search length ", false,
165 std::numeric_limits<double>::epsilon(), "SEARCH_LENGTH");
166 cmd.add(search_length_arg);
167
168 TCLAP::ValueArg<std::string> geometry_fname(
169 "i", "input-geometry",
170 "Input (.gml | .gli). The name of the input file containing the "
171 "geometry",
172 true, "", "INPUT_FILE");
173 cmd.add(geometry_fname);
174
175 TCLAP::ValueArg<std::string> mesh_arg(
176 "m", "mesh-file",
177 "Input (.vtu). The name of the input file containing the mesh", true,
178 "", "INPUT_FILE");
179 cmd.add(mesh_arg);
180
181 TCLAP::ValueArg<std::string> gmsh_path_arg(
182 "g", "gmsh-path", "Input (.msh). The path to the gmsh binary", false,
183 "", "INPUT_FILE");
184 cmd.add(gmsh_path_arg);
185
186 auto log_level_arg = BaseLib::makeLogLevelArg();
187 cmd.add(log_level_arg);
188 cmd.parse(argc, argv);
189
190 BaseLib::MPI::Setup mpi_setup(argc, argv);
191 BaseLib::initOGSLogger(log_level_arg.getValue());
192
193 // *** read mesh
194 INFO("Reading mesh '{:s}' ... ", mesh_arg.getValue());
195 std::unique_ptr<MeshLib::Mesh> subsurface_mesh(
196 MeshLib::IO::readMeshFromFile(mesh_arg.getValue()));
197 INFO("done.");
198 INFO("Extracting top surface of mesh '{:s}' ... ", mesh_arg.getValue());
199 Eigen::Vector3d const dir({0, 0, -1});
200 double const angle(90);
201 std::unique_ptr<MeshLib::Mesh> surface_mesh(
203 dir, angle));
204 INFO("done.");
205 subsurface_mesh.reset(nullptr);
206
207 // *** read geometry
208 GeoLib::GEOObjects geometries;
209 FileIO::readGeometryFromFile(geometry_fname.getValue(), geometries,
210 gmsh_path_arg.getValue());
211
212 auto const geo_name = geometries.getGeometryNames()[0];
213
214 // *** check if the data is usable
215 // *** get vector of polylines
216 std::vector<GeoLib::Polyline*> const* plys(
217 geometries.getPolylineVec(geo_name));
218 if (!plys)
219 {
220 ERR("Could not get vector of polylines out of geometry '{:s}'.",
221 geo_name);
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 return EXIT_FAILURE;
255 }
256
257 std::string merge_name("AllMeshNodesAlongPolylines");
258 if (geometry_sets.mergeGeometries(geo_names, merge_name) == 2)
259 {
260 merge_name = geo_names[0];
261 }
262
263 GeoLib::PointVec const* pnt_vec(geometry_sets.getPointVecObj(merge_name));
264 auto const& merged_pnts(pnt_vec->getVector());
265
266 std::vector<GeoLib::Point> pnts_with_id;
267 const std::size_t n_merged_pnts(merged_pnts.size());
268 for (std::size_t k(0); k < n_merged_pnts; ++k)
269 {
270 pnts_with_id.emplace_back(*(merged_pnts[k]), k);
271 }
272
273 std::sort(pnts_with_id.begin(), pnts_with_id.end(),
274 [](GeoLib::Point const& p0, GeoLib::Point const& p1)
275 { return p0 < p1; });
276
277 double const eps(std::numeric_limits<double>::epsilon());
278 std::vector<GeoLib::Point*> surface_pnts;
279 GeoLib::PointVec::NameIdMap name_id_map;
280
281 // insert first point
282 surface_pnts.push_back(
283 new GeoLib::Point(pnts_with_id[0], surface_pnts.size()));
284 {
285 std::string element_name;
286 pnt_vec->getNameOfElementByID(0, element_name);
287 name_id_map[element_name] = 0;
288 }
289 for (std::size_t k(1); k < n_merged_pnts; ++k)
290 {
291 const GeoLib::Point& p0(pnts_with_id[k - 1]);
292 const GeoLib::Point& p1(pnts_with_id[k]);
293 if (std::abs(p0[0] - p1[0]) > eps || std::abs(p0[1] - p1[1]) > eps)
294 {
295 surface_pnts.push_back(
296 new GeoLib::Point(pnts_with_id[k], surface_pnts.size()));
297 std::string element_name;
298 pnt_vec->getNameOfElementByID(k, element_name);
299 name_id_map[element_name] = surface_pnts.size() - 1;
300 }
301 }
302
303 std::string surface_name(BaseLib::dropFileExtension(mesh_arg.getValue()) +
304 "-MeshNodesAlongPolylines");
305 geometry_sets.addPointVec(std::move(surface_pnts), surface_name,
306 std::move(name_id_map), 1e-6);
307
308 // write the BCs and the merged geometry set to file
309 std::string const base_fname(
310 BaseLib::dropFileExtension(output_base_fname.getValue()));
311 writeBCsAndGeometry(geometry_sets, surface_name, base_fname,
312 bc_type.getValue(), gml_arg.getValue());
313 return EXIT_SUCCESS;
314}
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:36
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
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="")
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:64
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:227

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(), BaseLib::initOGSLogger(), BaseLib::makeLogLevelArg(), 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 86 of file CreateBoundaryConditionsAlongPolylines.cpp.

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

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

Referenced by writeBCsAndGeometry().

◆ writeLiquidFlowPointBC()

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

Definition at line 70 of file CreateBoundaryConditionsAlongPolylines.cpp.

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

Referenced by writeBCsAndGeometry().