OGS
Mesh2MeshPropertyInterpolation.cpp
Go to the documentation of this file.
1
16
17#include <fstream>
18#include <vector>
19
20#include "BaseLib/Logging.h"
21#include "GeoLib/AABB.h"
22#include "GeoLib/Grid.h"
24#include "MeshLib/Mesh.h"
25#include "MeshLib/MeshEnums.h"
26#include "MeshLib/Node.h"
27
28namespace MeshToolsLib
29{
31 MeshLib::Mesh const& src_mesh, std::string const& property_name)
32 : _src_mesh(src_mesh), _property_name(property_name)
33{
34}
35
37 MeshLib::Mesh& dest_mesh) const
38{
39 if (_src_mesh.getDimension() != dest_mesh.getDimension())
40 {
41 ERR("MeshLib::Mesh2MeshPropertyInterpolation::setPropertiesForMesh() "
42 "dimension of source (dim = {:d}) and destination (dim = {:d}) "
43 "mesh does not match.",
44 _src_mesh.getDimension(), dest_mesh.getDimension());
45 return false;
46 }
47
48 if (_src_mesh.getDimension() != 2)
49 {
50 WARN(
51 "MeshLib::Mesh2MeshPropertyInterpolation::setPropertiesForMesh() "
52 "implemented only for 2D case at the moment.");
53 return false;
54 }
55
56 MeshLib::PropertyVector<double>* dest_properties;
57 if (dest_mesh.getProperties().existsPropertyVector<double>(_property_name))
58 {
59 dest_properties =
61 assert(dest_properties->size() != dest_mesh.getNumberOfElements());
62 }
63 else
64 {
65 INFO("Create new PropertyVector '{:s}' of type double.",
67 dest_properties =
68 dest_mesh.getProperties().createNewPropertyVector<double>(
70 dest_mesh.getNumberOfElements(), 1);
71 if (!dest_properties)
72 {
73 WARN(
74 "Could not get or create a PropertyVector of type double using "
75 "the given name '{:s}'.",
77 return false;
78 }
79 }
80
81 interpolatePropertiesForMesh(dest_mesh, *dest_properties);
82
83 return true;
84}
85
87 MeshLib::Mesh const& dest_mesh,
88 MeshLib::PropertyVector<double>& dest_properties) const
89{
90 std::vector<double> interpolated_src_node_properties(
93 interpolated_src_node_properties);
94
95 // idea: looping over the destination elements and calculate properties
96 // from interpolated_src_node_properties to accelerate the (source) point
97 // search construct a grid
98 std::vector<MeshLib::Node*> const& src_nodes(_src_mesh.getNodes());
99 GeoLib::Grid<MeshLib::Node> src_grid(src_nodes.begin(), src_nodes.end(),
100 64);
101
102 auto const& dest_elements(dest_mesh.getElements());
103 const std::size_t n_dest_elements(dest_elements.size());
104 for (std::size_t k(0); k < n_dest_elements; k++)
105 {
106 MeshLib::Element& dest_element(*dest_elements[k]);
107 if (dest_element.getGeomType() == MeshLib::MeshElemType::LINE)
108 {
109 continue;
110 }
111
112 // compute axis aligned bounding box around the current element
113 const GeoLib::AABB elem_aabb(
114 dest_element.getNodes(),
115 dest_element.getNodes() + dest_element.getNumberOfBaseNodes());
116
117 // request "interesting" nodes from grid
118 std::vector<std::vector<MeshLib::Node*> const*> const nodes =
120 elem_aabb.getMinPoint(), elem_aabb.getMaxPoint());
121
122 std::size_t cnt(0);
123 double average_value(0.0);
124
125 for (auto const* nodes_vec : nodes)
126 {
127 for (auto const* node : *nodes_vec)
128 {
129 if (elem_aabb.containsPointXY(*node) &&
130 MeshLib::isPointInElementXY(*node, dest_element))
131 {
132 average_value +=
133 interpolated_src_node_properties[node->getID()];
134 cnt++;
135 }
136 }
137 }
138
139 if (cnt == 0)
140 {
141 OGS_FATAL(
142 "Mesh2MeshInterpolation: Could not find values in source mesh "
143 "for the element {:d}.",
144 k);
145 }
146 dest_properties[k] = average_value / cnt;
147 }
148}
149
152 std::vector<double>& interpolated_properties) const
153{
154 // fetch the source of property values
156 {
157 WARN("Did not find PropertyVector<double> '{:s}'.", _property_name);
158 return;
159 }
160 auto const* elem_props =
162
163 std::vector<MeshLib::Node*> const& src_nodes(_src_mesh.getNodes());
164 const std::size_t n_src_nodes(src_nodes.size());
165 for (std::size_t k(0); k < n_src_nodes; k++)
166 {
167 const std::size_t n_con_elems(
168 _src_mesh.getElementsConnectedToNode(*src_nodes[k]).size());
169 interpolated_properties[k] = (*elem_props)
170 [_src_mesh.getElementsConnectedToNode(*src_nodes[k])[0]->getID()];
171 for (std::size_t j(1); j < n_con_elems; j++)
172 {
173 interpolated_properties[k] +=
174 (*elem_props)[_src_mesh
175 .getElementsConnectedToNode(*src_nodes[k])[j]
176 ->getID()];
177 }
178 interpolated_properties[k] /= n_con_elems;
179 }
180}
181
182} // namespace MeshToolsLib
Definition of the AABB class.
Definition of the Element class.
#define OGS_FATAL(...)
Definition Error.h:26
Definition of the Grid class.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
Implementation of the Mesh2MeshPropertyInterpolation class.
Definition of mesh-related Enumerations.
Definition of the Mesh class.
Definition of the Node class.
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
Eigen::Vector3d const & getMaxPoint() const
Definition AABB.h:187
Eigen::Vector3d const & getMinPoint() const
Definition AABB.h:180
bool containsPointXY(T const &pnt) const
Definition AABB.h:161
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCuboid(Eigen::Vector3d const &min_pnt, Eigen::Vector3d const &max_pnt) const
Definition Grid.h:286
virtual MeshElemType getGeomType() const =0
virtual unsigned getNumberOfBaseNodes() const =0
virtual Node *const * getNodes() const =0
Get array of element nodes.
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:90
Properties & getProperties()
Definition Mesh.h:136
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:102
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:257
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
bool existsPropertyVector(std::string_view name) const
Definition Properties.h:93
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
Definition Properties.h:17
PropertyVector< T > const * getPropertyVector(std::string_view name) const
constexpr std::size_t size() const
void interpolatePropertiesForMesh(MeshLib::Mesh const &dest_mesh, MeshLib::PropertyVector< double > &dest_properties) const
void interpolateElementPropertiesToNodeProperties(std::vector< double > &interpolated_properties) const
Mesh2MeshPropertyInterpolation(MeshLib::Mesh const &src_mesh, std::string const &property_name)
bool isPointInElementXY(MathLib::Point3d const &p, Element const &e)
Definition Element.cpp:189