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