37 auto const* surface_material_ids = materialIDs(*sfc_mesh);
38 auto* subsurface_material_ids = materialIDs(*subsurface_mesh);
40 if (!surface_material_ids)
42 OGS_FATAL(
"Surface mesh does not contain material IDs");
44 if (surface_material_ids->empty())
46 OGS_FATAL(
"Surface mesh material IDs doesn't contain any values");
48 if (!subsurface_material_ids)
50 OGS_FATAL(
"Subsurface mesh does not contain material IDs");
52 if (subsurface_material_ids->size() % surface_material_ids->size() != 0)
54 OGS_FATAL(
"Could not determine the number of subsurface layers.");
56 int const number_of_layers =
57 subsurface_material_ids->size() / surface_material_ids->size();
60 ranges::views::repeat_n(*surface_material_ids, number_of_layers) |
62 subsurface_material_ids->begin());
65int main(
int argc,
char* argv[])
68 "Creates a layered 3D OGS mesh from an existing 2D OGS mesh and a list "
69 "of raster files representing subsurface layers. Supported raster "
70 "formats are ArcGIS ascii rasters (*.asc), Surfer Grids (*.grd), or "
71 "gridded XYZ rasters (*.xyz)."
72 "Only input meshes consisting of line and triangle elements are "
73 "currently supported as mapping of quads might result in invalid mesh "
75 "OpenGeoSys-6 software, version " +
78 "Copyright (c) 2012-2025, OpenGeoSys Community "
79 "(http://www.opengeosys.org)",
82 TCLAP::SwitchArg use_ascii_arg(
"",
"ascii_output",
83 "Write VTU output in ASCII format.");
84 cmd.add(use_ascii_arg);
86 double min_thickness(std::numeric_limits<double>::epsilon());
87 TCLAP::ValueArg<double> min_thickness_arg(
89 "The minimum thickness of a layer to be integrated at any given "
91 false, min_thickness,
"floating point number");
92 cmd.add(min_thickness_arg);
94 TCLAP::ValueArg<std::string> raster_path_arg(
96 "An ascii-file containing a list of raster files, starting from top "
98 true,
"",
"file name");
99 cmd.add(raster_path_arg);
101 TCLAP::SwitchArg keep_materials_arg(
102 "",
"keep-surface-material-ids",
103 "if the argument is present the materials defined in the surface mesh "
104 "are used to set the material information for the subsurface cells",
106 cmd.add(keep_materials_arg);
108 TCLAP::ValueArg<std::string> mesh_out_arg(
109 "o",
"output-mesh-file",
"The file name of the resulting 3D mesh.",
110 true,
"",
"file name");
111 cmd.add(mesh_out_arg);
113 TCLAP::ValueArg<std::string> mesh_arg(
"i",
"input-mesh-file",
114 "The file name of the 2D input mesh.",
115 true,
"",
"file name");
118 cmd.parse(argc, argv);
122 if (min_thickness_arg.isSet())
124 min_thickness = min_thickness_arg.getValue();
125 if (min_thickness < 0)
127 ERR(
"Minimum layer thickness must be non-negative value.");
132 INFO(
"Reading mesh '{:s}' ... ", mesh_arg.getValue());
133 std::unique_ptr<MeshLib::Mesh>
const sfc_mesh(
137 ERR(
"Error reading mesh '{:s}'.", mesh_arg.getValue());
140 if (sfc_mesh->getDimension() != 2)
142 ERR(
"Input mesh must be a 2D mesh.");
147 std::vector<std::string> raster_paths =
149 if (raster_paths.size() < 2)
151 ERR(
"At least two raster files needed to create 3D mesh.");
154 std::reverse(raster_paths.begin(), raster_paths.end());
159 if (!mapper.
createLayers(*sfc_mesh, *rasters, min_thickness))
166 ERR(
"Reading raster files.");
170 auto const result_mesh = mapper.
getMesh(
"SubsurfaceMesh");
171 if (result_mesh ==
nullptr)
173 ERR(
"Mapper returned empty result for 'SubsurfaceMesh'.");
177 if (keep_materials_arg.getValue())
183 std::string output_name(mesh_out_arg.getValue());
186 output_name.append(
".vtu");
189 INFO(
"Writing mesh '{:s}' ... ", output_name);
190 auto const data_mode =
191 use_ascii_arg.getValue() ? vtkXMLWriter::Ascii : vtkXMLWriter::Binary;