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