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