129int main(
int argc,
char* argv[])
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 " +
139 "Copyright (c) 2012-2024, OpenGeoSys Community "
140 "(http://www.opengeosys.org)",
142 TCLAP::SwitchArg gml_arg(
"",
"gml",
"Write found nodes to gml file.");
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);
152 TCLAP::ValueArg<std::string> bc_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",
158 "process type as string (LIQUID_FLOW or GROUNDWATER_FLOW (default))");
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);
168 TCLAP::ValueArg<std::string> geometry_fname(
169 "i",
"input-geometry",
170 "the name of the file containing the input geometry",
true,
"",
172 cmd.add(geometry_fname);
174 TCLAP::ValueArg<std::string> mesh_arg(
175 "m",
"mesh-file",
"the name of the file containing the mesh",
true,
"",
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);
184 cmd.parse(argc, argv);
187 MPI_Init(&argc, &argv);
191 INFO(
"Reading mesh '{:s}' ... ", mesh_arg.getValue());
192 std::unique_ptr<MeshLib::Mesh> subsurface_mesh(
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(
202 subsurface_mesh.reset(
nullptr);
207 gmsh_path_arg.getValue());
213 std::vector<GeoLib::Polyline*>
const* plys(
217 ERR(
"Could not get vector of polylines out of geometry '{:s}'.",
225 auto search_length_strategy =
226 std::make_unique<MeshGeoToolsLib::SearchLength>();
227 if (search_length_arg.isSet())
229 search_length_strategy.reset(
235 *surface_mesh, std::move(search_length_strategy),
237 for (std::size_t k(0); k < plys->size(); k++)
244 std::string polyline_name(
"Polyline-" + std::to_string(k));
251 if (geo_names.empty())
253 ERR(
"Did not find mesh nodes along polylines.");
260 std::string merge_name(
"AllMeshNodesAlongPolylines");
263 merge_name = geo_names[0];
267 auto const& merged_pnts(pnt_vec->
getVector());
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)
273 pnts_with_id.emplace_back(*(merged_pnts[k]), k);
276 std::sort(pnts_with_id.begin(), pnts_with_id.end(),
278 { return p0 < p1; });
280 double const eps(std::numeric_limits<double>::epsilon());
281 std::vector<GeoLib::Point*> surface_pnts;
285 surface_pnts.push_back(
288 std::string element_name;
290 name_id_map[element_name] = 0;
292 for (std::size_t k(1); k < n_merged_pnts; ++k)
296 if (std::abs(p0[0] - p1[0]) > eps || std::abs(p0[1] - p1[1]) > eps)
298 surface_pnts.push_back(
300 std::string element_name;
302 name_id_map[element_name] = surface_pnts.size() - 1;
307 "-MeshNodesAlongPolylines");
308 geometry_sets.
addPointVec(std::move(surface_pnts), surface_name,
309 std::move(name_id_map), 1e-6);
312 std::string
const base_fname(
315 bc_type.getValue(), gml_arg.getValue());