OGS
ipDataToPointCloud.cpp
Go to the documentation of this file.
1
10#include <tclap/CmdLine.h>
11
12#include <unordered_map>
13
14#include "BaseLib/Logging.h"
15#include "BaseLib/MPI.h"
17#include "InfoLib/GitInfo.h"
20#include "MeshLib/Mesh.h"
24
25template <typename MeshElement>
26std::vector<std::array<double, 3>> interpolateElementCoords(
27 MeshLib::Element const& e, NumLib::ShapeMatrixCache const& sm_cache)
28{
29 constexpr int dim = MeshElement::dimension;
30 using ShapeFunction =
33
34 auto const& Ns = sm_cache.NsHigherOrder<MeshElement>();
35 std::vector<std::array<double, 3>> coords;
36 coords.reserve(Ns.size());
37
38 for (auto const& N : Ns)
39 {
40 coords.push_back(
42 }
43
44 return coords;
45}
46
47// A template metafunction "returning" the cell type of a mesh element
48template <typename Element>
50
51template <typename ElementRule>
52struct CellTypeOfTemplateElement<MeshLib::TemplateElement<ElementRule>>
53{
54 static constexpr MeshLib::CellType value = ElementRule::cell_type;
55};
56
57std::unordered_map<MeshLib::CellType, std::vector<std::array<double, 3>> (*)(
58 MeshLib::Element const&,
61{
62 std::unordered_map<MeshLib::CellType,
63 std::vector<std::array<double, 3>> (*)(
64 MeshLib::Element const&,
66 map_cell_type_to_element_coords_interpolator;
67
68 map_cell_type_to_element_coords_interpolator.reserve(
69 boost::mp11::mp_size<NumLib::AllElementTraitsLagrange>::value);
70
71 boost::mp11::mp_for_each<NumLib::AllElementTraitsLagrange>(
72 [&map_cell_type_to_element_coords_interpolator]<typename ETL>(ETL)
73 {
74 using MeshElement = typename ETL::Element;
75 auto constexpr cell_type =
77
78 auto const [it, newly_inserted] =
79 map_cell_type_to_element_coords_interpolator.emplace(
81
82 if (!newly_inserted)
83 {
85 "Coordinate interpolator for cell type {} is already "
86 "present",
87 MeshLib::CellType2String(cell_type));
88 }
89 });
90
91 return map_cell_type_to_element_coords_interpolator;
92}
93
94std::vector<MeshLib::Node*> computePointCloudNodes(
95 MeshLib::Mesh const& mesh, unsigned const integration_order)
96{
97 NumLib::ShapeMatrixCache sm_cache{integration_order};
98 auto const map_cell_type_to_element_coords_interpolator =
100
101 std::vector<MeshLib::Node*> nodes;
102
103 for (auto const* element : mesh.getElements())
104 {
105 auto const cell_type = element->getCellType();
106
107 auto const it =
108 map_cell_type_to_element_coords_interpolator.find(cell_type);
109 if (it == map_cell_type_to_element_coords_interpolator.end())
110 {
111 OGS_FATAL("Unsupported cell type {} for element #{}",
112 MeshLib::CellType2String(cell_type), element->getID());
113 }
114
115 auto const& element_coords_interpolator = it->second;
116 auto const coords = element_coords_interpolator(*element, sm_cache);
117
118 for (auto& cs : coords)
119 {
120 nodes.push_back(new MeshLib::Node(cs, nodes.size()));
121 }
122 }
123
124 return nodes;
125}
126
128{
129 auto const& properties = mesh.getProperties();
130 auto const ip_meta_data = MeshLib::getIntegrationPointMetaData(properties);
131
132 std::optional<unsigned> integration_order;
133
134 for (auto const& [name, property] : properties)
135 {
136 if (property->getMeshItemType() !=
138 {
139 continue;
140 }
141
142 if (name == "IntegrationPointMetaData" || name == "OGS_VERSION")
143 {
144 continue;
145 }
146
147 auto const order =
150
151 if (!integration_order)
152 {
153 integration_order = order;
154 }
155 else if (integration_order != order)
156 {
157 OGS_FATAL("Integration orders differ: {} != {}", *integration_order,
158 order);
159 }
160 }
161
162 if (!integration_order)
163 {
164 OGS_FATAL("Integration order could not be determined.");
165 }
166
167 return *integration_order;
168}
169
171 MeshLib::Properties& props_out,
172 std::size_t const num_ips)
173{
174 for (auto const& [prop_name, prop_in] : props_in)
175 {
176 auto const mesh_item_type = prop_in->getMeshItemType();
177
178 if (mesh_item_type != MeshLib::MeshItemType::IntegrationPoint)
179 {
180 continue;
181 }
182
183 auto const num_comp = prop_in->getNumberOfGlobalComponents();
184 auto const* prop_in_double =
185 dynamic_cast<MeshLib::PropertyVector<double> const*>(prop_in);
186
187 if (prop_in_double == nullptr)
188 {
189 INFO(
190 "Ignoring non-double-valued property '{}' with {} components "
191 "on {}",
192 prop_name, num_comp, toString(mesh_item_type));
193 continue;
194 }
195
196 DBUG("Converting property '{}' with {} components on {}", prop_name,
197 num_comp, toString(mesh_item_type));
198
199 if (props_out.existsPropertyVector<double>(prop_name))
200 {
201 OGS_FATAL(
202 "A property vector with name '{}' already exists. Not "
203 "adding it again",
204 prop_name);
205 }
206
207 if (auto const num_ips_actual = prop_in_double->getNumberOfTuples();
208 num_ips_actual != num_ips)
209 {
210 WARN(
211 "Property vector '{}' has {} tuples, which differs from "
212 "the number of integration points ({}). Skipping this "
213 "property.",
214 prop_name, num_ips_actual, num_ips);
215 }
216
217 auto* prop_out = props_out.createNewPropertyVector<double>(
218 prop_name, MeshLib::MeshItemType::Node, num_comp);
219
220 prop_out->assign(*prop_in_double);
221 }
222}
223
224int main(int argc, char** argv)
225{
226 // TODO future additions to this tool might include:
227 // -C --copy-cell-data
228 // -N --interpolate-node-data
229 // -I --add-cell-ids
230 // -O --integration-order
231 // --natural-coords add natural coordinates of integration points, or better
232 // not?
233 // --no-data
234 // --linear-shape-functions={pressure, ...} for the interpolation of nodal
235 // data
236 // PVD support
237 // handle other data than only double data
238
239 TCLAP::CmdLine cmd(
240 "Creates a point cloud mesh for the integration point data existing on "
241 "a given input mesh.\n\n"
242 "OpenGeoSys-6 software, version " +
244 ".\n"
245 "Copyright (c) 2012-2025, OpenGeoSys Community "
246 "(http://www.opengeosys.org)",
248 TCLAP::ValueArg<std::string> arg_out_file(
249 "o", "output-file", "Output (.vtu). The output mesh (point cloud)",
250 true, "", "OUTPUT_FILE");
251 cmd.add(arg_out_file);
252 TCLAP::ValueArg<std::string> arg_in_file("i", "input-file",
253 "Input (.vtu). The input mesh",
254 true, "", "INPUT_FILE");
255 cmd.add(arg_in_file);
256
257 auto log_level_arg = BaseLib::makeLogLevelArg();
258 cmd.add(log_level_arg);
259 cmd.parse(argc, argv);
260
261 BaseLib::MPI::Setup mpi_setup(argc, argv);
262 BaseLib::initOGSLogger(log_level_arg.getValue());
263
264 std::unique_ptr<MeshLib::Mesh const> mesh_in(
265 MeshLib::IO::readMeshFromFile(arg_in_file.getValue()));
266
267 auto const integration_order = determineIntegrationOrder(*mesh_in);
268 auto nodes = computePointCloudNodes(*mesh_in, integration_order);
269
270 MeshLib::Mesh point_cloud("point_cloud", std::move(nodes), {});
271
272 copyDoubleValuedFieldDataToPointCloud(mesh_in->getProperties(),
273 point_cloud.getProperties(),
274 point_cloud.getNumberOfNodes());
275
276 MeshLib::IO::writeMeshToFile(point_cloud, arg_out_file.getValue());
277
278 return EXIT_SUCCESS;
279}
#define OGS_FATAL(...)
Definition Error.h:26
Git information.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:42
Definition of the Mesh class.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
Properties & getProperties()
Definition Mesh.h:136
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition Properties.h:33
bool existsPropertyVector(std::string_view name) const
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
boost::mp11::mp_at< ShapeFunctionsHigherOrder, boost::mp11::mp_find< MeshElements, MeshElement > > ShapeFunctionHigherOrder
auto const & NsHigherOrder() const
unsigned determineIntegrationOrder(MeshLib::Mesh const &mesh)
int main(int argc, char **argv)
std::vector< MeshLib::Node * > computePointCloudNodes(MeshLib::Mesh const &mesh, unsigned const integration_order)
std::unordered_map< MeshLib::CellType, std::vector< std::array< double, 3 > >(*)(MeshLib::Element const &, NumLib::ShapeMatrixCache const &)> createElementCoordInterpolatorsForAllElementTypes()
void copyDoubleValuedFieldDataToPointCloud(MeshLib::Properties const &props_in, MeshLib::Properties &props_out, std::size_t const num_ips)
std::vector< std::array< double, 3 > > interpolateElementCoords(MeshLib::Element const &e, NumLib::ShapeMatrixCache const &sm_cache)
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:64
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
int writeMeshToFile(const MeshLib::Mesh &mesh, std::filesystem::path const &file_path, std::set< std::string > variable_output_names)
CellType
Types of mesh elements supported by OpenGeoSys.
Definition MeshEnums.h:64
std::optional< IntegrationPointMetaData > getIntegrationPointMetaData(MeshLib::Properties const &properties)
std::string CellType2String(const CellType t)
Given a MeshElemType this returns the appropriate string.
IntegrationPointMetaDataSingleField getIntegrationPointMetaDataSingleField(std::optional< IntegrationPointMetaData > const &ip_meta_data, std::string const &field_name)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition of readMeshFromFile function.