OGS
MeshSurfaceExtraction.cpp
Go to the documentation of this file.
1 
15 #include "MeshSurfaceExtraction.h"
16 
17 #include <boost/math/constants/constants.hpp>
18 #include <memory>
19 
20 #include "BaseLib/Logging.h"
21 #include "MeshLib/Elements/Line.h"
22 #include "MeshLib/Elements/Point.h"
23 #include "MeshLib/Elements/Quad.h"
24 #include "MeshLib/Elements/Tri.h"
25 #include "MeshLib/Mesh.h"
29 
30 namespace MeshLib
31 {
32 template <typename T>
34  std::vector<std::size_t> const& id_map,
35  MeshLib::Mesh& sfc_mesh)
36 {
37  auto const number_of_components = property.getNumberOfGlobalComponents();
38 
39  auto sfc_prop = getOrCreateMeshProperty<T>(
40  sfc_mesh, property.getPropertyName(), property.getMeshItemType(),
41  number_of_components);
42  sfc_prop->clear();
43  sfc_prop->reserve(id_map.size());
44 
45  for (auto bulk_id : id_map)
46  {
47  std::copy_n(&property.getComponent(bulk_id, 0 /*component_id*/),
48  number_of_components, back_inserter(*sfc_prop));
49  }
50 }
51 
53  MeshLib::Properties const& properties,
54  std::vector<std::size_t> const& node_ids_map,
55  std::vector<std::size_t> const& element_ids_map)
56 {
57  std::size_t const n_elems(sfc_mesh.getNumberOfElements());
58  std::size_t const n_nodes(sfc_mesh.getNumberOfNodes());
59  if (node_ids_map.size() != n_nodes)
60  {
61  ERR("createSfcMeshProperties() - Incorrect number of node IDs ({:d}) "
62  "compared to actual number of surface nodes ({:d}).",
63  node_ids_map.size(), n_nodes);
64  return false;
65  }
66 
67  if (element_ids_map.size() != n_elems)
68  {
69  ERR("createSfcMeshProperties() - Incorrect number of element IDs "
70  "({:d}) compared to actual number of surface elements ({:d}).",
71  element_ids_map.size(), n_elems);
72  return false;
73  }
74  std::map<MeshLib::MeshItemType, std::vector<std::size_t> const*> const
75  id_maps = {{MeshItemType::Cell, &element_ids_map},
76  {MeshItemType::Node, &node_ids_map}};
77 
78  std::size_t vectors_copied(0);
79  std::size_t vectors_skipped(0);
80  for (auto [name, property] : properties)
81  {
82  if (property->getMeshItemType() != MeshItemType::Cell &&
83  property->getMeshItemType() != MeshItemType::Node)
84  {
85  WARN(
86  "Skipping property vector '{:s}' - not defined on cells or "
87  "nodes.",
88  name);
89  vectors_skipped++;
90  continue;
91  }
92 
93  auto const& id_map = *id_maps.at(property->getMeshItemType());
94  if (auto p = dynamic_cast<PropertyVector<double>*>(property))
95  {
96  processPropertyVector(*p, id_map, sfc_mesh);
97  vectors_copied++;
98  }
99  else if (auto p = dynamic_cast<PropertyVector<float>*>(property))
100  {
101  processPropertyVector(*p, id_map, sfc_mesh);
102  vectors_copied++;
103  }
104  else if (auto p = dynamic_cast<PropertyVector<int>*>(property))
105  {
106  processPropertyVector(*p, id_map, sfc_mesh);
107  vectors_copied++;
108  }
109  else if (auto p = dynamic_cast<PropertyVector<unsigned>*>(property))
110  {
111  processPropertyVector(*p, id_map, sfc_mesh);
112  vectors_copied++;
113  }
114  else if (auto p = dynamic_cast<PropertyVector<long>*>(property))
115  {
116  processPropertyVector(*p, id_map, sfc_mesh);
117  vectors_copied++;
118  }
119  else if (auto p = dynamic_cast<PropertyVector<long long>*>(property))
120  {
121  processPropertyVector(*p, id_map, sfc_mesh);
122  vectors_copied++;
123  }
124  else if (auto p =
125  dynamic_cast<PropertyVector<unsigned long>*>(property))
126  {
127  processPropertyVector(*p, id_map, sfc_mesh);
128  vectors_copied++;
129  }
130  else if (auto p = dynamic_cast<PropertyVector<unsigned long long>*>(
131  property))
132  {
133  processPropertyVector(*p, id_map, sfc_mesh);
134  vectors_copied++;
135  }
136  else if (auto p = dynamic_cast<PropertyVector<std::size_t>*>(property))
137  {
138  processPropertyVector(*p, id_map, sfc_mesh);
139  vectors_copied++;
140  }
141  else if (auto p = dynamic_cast<PropertyVector<char>*>(property))
142  {
143  processPropertyVector(*p, id_map, sfc_mesh);
144  vectors_copied++;
145  }
146  else
147  {
148  WARN(
149  "Skipping property vector '{:s}' - no matching data type "
150  "'{:s}' found.",
151  name, typeid(*property).name());
152  vectors_skipped++;
153  }
154  }
155  INFO("{:d} property vectors copied, {:d} vectors skipped.", vectors_copied,
156  vectors_skipped);
157  return true;
158 }
159 
160 std::tuple<std::vector<MeshLib::Node*>, std::vector<std::size_t>>
161 createNodesAndIDMapFromElements(std::vector<MeshLib::Element*> const& elements,
162  std::size_t const n_all_nodes)
163 {
164  std::vector<MeshLib::Node const*> tmp_nodes(n_all_nodes, nullptr);
165  for (auto const* elem : elements)
166  {
167  auto const n_nodes = elem->getNumberOfNodes();
168  for (unsigned j = 0; j < n_nodes; ++j)
169  {
170  const MeshLib::Node* node(elem->getNode(j));
171  tmp_nodes[node->getID()] = node;
172  }
173  }
174 
175  std::vector<MeshLib::Node*> nodes;
176  std::vector<std::size_t> node_id_map(n_all_nodes);
177  for (unsigned i = 0; i < n_all_nodes; ++i)
178  {
179  if (tmp_nodes[i])
180  {
181  node_id_map[i] = nodes.size();
182  nodes.push_back(new MeshLib::Node(*tmp_nodes[i]));
183  }
184  }
185  return {nodes, node_id_map};
186 }
187 
189  const MeshLib::Mesh& mesh)
190 {
191  std::vector<double> node_area_vec;
192  if (mesh.getDimension() != 2)
193  {
194  ERR("Error in MeshSurfaceExtraction::getSurfaceAreaForNodes() - Given "
195  "mesh is no surface mesh (dimension != 2).");
196  return node_area_vec;
197  }
198 
199  double total_area(0);
200 
201  // for each node, a vector containing all the element idget every element
202  const std::vector<MeshLib::Node*>& nodes = mesh.getNodes();
203  const std::size_t nNodes(mesh.getNumberOfNodes());
204  for (std::size_t n = 0; n < nNodes; ++n)
205  {
206  double node_area(0);
207 
208  auto const conn_elems = mesh.getElementsConnectedToNode(*nodes[n]);
209  const std::size_t nConnElems(conn_elems.size());
210 
211  for (std::size_t i = 0; i < nConnElems; ++i)
212  {
213  const MeshLib::Element* elem(conn_elems[i]);
214  const unsigned nElemParts =
215  (elem->getGeomType() == MeshElemType::TRIANGLE) ? 3 : 4;
216  const double area = conn_elems[i]->getContent() / nElemParts;
217  node_area += area;
218  total_area += area;
219  }
220 
221  node_area_vec.push_back(node_area);
222  }
223 
224  INFO("Total surface Area: {:f}", total_area);
225 
226  return node_area_vec;
227 }
228 
230  const MeshLib::Mesh& subsfc_mesh, Eigen::Vector3d const& dir, double angle,
231  std::string const& subsfc_node_id_prop_name,
232  std::string const& subsfc_element_id_prop_name,
233  std::string const& face_id_prop_name)
234 {
235  // allow slightly greater angles than 90 degrees for numerical errors
236  if (angle < 0 || angle > 91)
237  {
238  ERR("Supported angle between 0 and 90 degrees only.");
239  return nullptr;
240  }
241 
242  INFO("Extracting mesh surface...");
243  std::vector<MeshLib::Element*> sfc_elements;
244  std::vector<std::size_t> element_ids_map;
245  std::vector<std::size_t> face_ids_map;
246  get2DSurfaceElements(subsfc_mesh.getElements(), sfc_elements,
247  element_ids_map, face_ids_map, dir, angle,
248  subsfc_mesh.getDimension());
249 
250  if (sfc_elements.empty())
251  {
252  return nullptr;
253  }
254 
255  auto [sfc_nodes, node_id_map] = createNodesAndIDMapFromElements(
256  sfc_elements, subsfc_mesh.getNumberOfNodes());
257 
258  // create new elements vector with newly created nodes (and delete
259  // temp-elements)
260  auto new_elements =
261  MeshLib::copyElementVector(sfc_elements, sfc_nodes, &node_id_map);
262  std::for_each(sfc_elements.begin(), sfc_elements.end(),
263  [](MeshLib::Element* e) { delete e; });
264 
265  std::vector<std::size_t> id_map;
266  id_map.reserve(sfc_nodes.size());
267  std::transform(begin(sfc_nodes), end(sfc_nodes), std::back_inserter(id_map),
268  [](MeshLib::Node* const n) { return n->getID(); });
269 
270  MeshLib::Mesh* result(
271  new Mesh(subsfc_mesh.getName() + "-Surface", sfc_nodes, new_elements));
272 
273  addBulkIDPropertiesToMesh(*result, subsfc_node_id_prop_name, id_map,
274  subsfc_element_id_prop_name, element_ids_map,
275  face_id_prop_name, face_ids_map);
276 
277  if (!createSfcMeshProperties(*result, subsfc_mesh.getProperties(), id_map,
278  element_ids_map))
279  {
280  ERR("Transferring subsurface properties failed.");
281  }
282  return result;
283 }
284 
286  const std::vector<MeshLib::Element*>& all_elements,
287  std::vector<MeshLib::Element*>& sfc_elements,
288  std::vector<std::size_t>& element_to_bulk_element_id_map,
289  std::vector<std::size_t>& element_to_bulk_face_id_map,
290  Eigen::Vector3d const& dir, double angle, unsigned mesh_dimension)
291 {
292  if (mesh_dimension < 2 || mesh_dimension > 3)
293  {
294  ERR("Cannot handle meshes of dimension {:i}", mesh_dimension);
295  }
296 
297  bool const complete_surface = (dir.dot(dir) == 0);
298 
299  double const pi(boost::math::constants::pi<double>());
300  double const cos_theta(std::cos(angle * pi / 180.0));
301  Eigen::Vector3d const norm_dir(dir.normalized());
302 
303  for (auto const* elem : all_elements)
304  {
305  const unsigned element_dimension(elem->getDimension());
306  if (element_dimension < mesh_dimension)
307  {
308  continue;
309  }
310 
311  if (element_dimension == 2)
312  {
313  if (!complete_surface)
314  {
315  if (FaceRule::getSurfaceNormal(*elem).normalized().dot(
316  norm_dir) > cos_theta)
317  {
318  continue;
319  }
320  }
321  sfc_elements.push_back(elem->clone());
322  element_to_bulk_element_id_map.push_back(elem->getID());
323  element_to_bulk_face_id_map.push_back(0);
324  }
325  else
326  {
327  if (!elem->isBoundaryElement())
328  {
329  continue;
330  }
331  const unsigned nFaces(elem->getNumberOfFaces());
332  for (unsigned j = 0; j < nFaces; ++j)
333  {
334  if (elem->getNeighbor(j) != nullptr)
335  {
336  continue;
337  }
338 
339  auto const face =
340  std::unique_ptr<MeshLib::Element const>{elem->getFace(j)};
341  if (!complete_surface)
342  {
343  if (FaceRule::getSurfaceNormal(*face).normalized().dot(
344  norm_dir) < cos_theta)
345  {
346  continue;
347  }
348  }
349  if (face->getGeomType() == MeshElemType::TRIANGLE)
350  {
351  sfc_elements.push_back(new MeshLib::Tri(
352  *static_cast<const MeshLib::Tri*>(face.get())));
353  }
354  else
355  {
356  sfc_elements.push_back(new MeshLib::Quad(
357  *static_cast<const MeshLib::Quad*>(face.get())));
358  }
359  element_to_bulk_element_id_map.push_back(elem->getID());
360  element_to_bulk_face_id_map.push_back(j);
361  }
362  }
363  }
364 }
365 
366 std::vector<MeshLib::Node*> MeshSurfaceExtraction::getSurfaceNodes(
367  const MeshLib::Mesh& mesh, Eigen::Vector3d const& dir, double angle)
368 {
369  INFO("Extracting surface nodes...");
370  std::vector<MeshLib::Element*> sfc_elements;
371  std::vector<std::size_t> element_to_bulk_element_id_map;
372  std::vector<std::size_t> element_to_bulk_face_id_map;
373 
375  mesh.getElements(), sfc_elements, element_to_bulk_element_id_map,
376  element_to_bulk_face_id_map, dir, angle, mesh.getDimension());
377 
378  std::vector<MeshLib::Node*> surface_nodes;
379  std::tie(surface_nodes, std::ignore) =
381 
382  for (auto e : sfc_elements)
383  {
384  delete e;
385  }
386 
387  return surface_nodes;
388 }
389 
391  MeshLib::Element const& surface_element,
392  std::vector<MeshLib::Element*>& surface_elements,
393  std::vector<std::size_t>& element_to_bulk_element_id_map,
394  std::vector<std::size_t>& element_to_bulk_face_id_map)
395 {
396  const unsigned n_faces(surface_element.getNumberOfBoundaries());
397  for (unsigned j = 0; j < n_faces; ++j)
398  {
399  if (surface_element.getNeighbor(j) != nullptr)
400  {
401  continue;
402  }
403 
404  surface_elements.push_back(
405  const_cast<MeshLib::Element*>(surface_element.getBoundary(j)));
406  element_to_bulk_face_id_map.push_back(j);
407  element_to_bulk_element_id_map.push_back(surface_element.getID());
408  }
409 }
410 
411 std::tuple<std::vector<MeshLib::Element*>, std::vector<std::size_t>,
412  std::vector<std::size_t>>
414 {
415  std::vector<std::size_t> element_to_bulk_element_id_map;
416  std::vector<std::size_t> element_to_bulk_face_id_map;
417  std::vector<MeshLib::Element*> surface_elements;
418 
419  auto const& bulk_elements = bulk_mesh.getElements();
420  auto const mesh_dimension = bulk_mesh.getDimension();
421 
422  for (auto const* elem : bulk_elements)
423  {
424  const unsigned element_dimension(elem->getDimension());
425  if (element_dimension < mesh_dimension)
426  {
427  continue;
428  }
429 
430  if (!elem->isBoundaryElement())
431  {
432  continue;
433  }
434  createSurfaceElementsFromElement(*elem, surface_elements,
435  element_to_bulk_element_id_map,
436  element_to_bulk_face_id_map);
437  }
438  return {surface_elements, element_to_bulk_element_id_map,
439  element_to_bulk_face_id_map};
440 }
441 
442 namespace BoundaryExtraction
443 {
444 std::unique_ptr<MeshLib::Mesh> getBoundaryElementsAsMesh(
445  MeshLib::Mesh const& bulk_mesh,
446  std::string const& subsfc_node_id_prop_name,
447  std::string const& subsfc_element_id_prop_name,
448  std::string const& face_id_prop_name)
449 {
450  auto const mesh_dimension = bulk_mesh.getDimension();
451  if (mesh_dimension < 2 || mesh_dimension > 3)
452  {
453  ERR("Cannot handle meshes of dimension {:i}", mesh_dimension);
454  }
455 
456  // create boundary elements based on the subsurface nodes
457  auto [boundary_elements, element_to_bulk_element_id_map,
458  element_to_bulk_face_id_map] = createBoundaryElements(bulk_mesh);
459 
460  // create new nodes needed for the new surface elements
461  auto [boundary_nodes, node_id_map] = createNodesAndIDMapFromElements(
462  boundary_elements, bulk_mesh.getNumberOfNodes());
463 
464  // create new elements using newly created nodes and delete temp-elements
465  auto new_elements = MeshLib::copyElementVector(
466  boundary_elements, boundary_nodes, &node_id_map);
467  for (auto* e : boundary_elements)
468  {
469  delete e;
470  }
471 
472  std::vector<std::size_t> nodes_to_bulk_nodes_id_map;
473  nodes_to_bulk_nodes_id_map.reserve(boundary_nodes.size());
474  std::transform(begin(boundary_nodes), end(boundary_nodes),
475  std::back_inserter(nodes_to_bulk_nodes_id_map),
476  [](MeshLib::Node* const n) { return n->getID(); });
477 
478  std::unique_ptr<MeshLib::Mesh> boundary_mesh(new Mesh(
479  bulk_mesh.getName() + "-boundary", boundary_nodes, new_elements));
480 
482  *boundary_mesh, subsfc_node_id_prop_name, nodes_to_bulk_nodes_id_map,
483  subsfc_element_id_prop_name, element_to_bulk_element_id_map,
484  face_id_prop_name, element_to_bulk_face_id_map);
485 
487  if (!createSfcMeshProperties(*boundary_mesh, bulk_mesh.getProperties(),
488  nodes_to_bulk_nodes_id_map,
489  element_to_bulk_element_id_map))
490  {
491  ERR("Transferring subsurface properties failed.");
492  }
493  return boundary_mesh;
494 }
495 
496 } // namespace BoundaryExtraction
497 
499  MeshLib::Mesh& surface_mesh,
500  std::string const& node_to_bulk_node_id_map_name,
501  std::vector<std::size_t> const& node_to_bulk_node_id_map,
502  std::string const& element_to_bulk_element_id_map_name,
503  std::vector<std::size_t> const& element_to_bulk_element_id_map,
504  std::string const& element_to_bulk_face_id_map_name,
505  std::vector<std::size_t> const& element_to_bulk_face_id_map)
506 {
507  // transmit the original node ids of the bulk mesh as a property
508  if (!node_to_bulk_node_id_map_name.empty())
509  {
510  MeshLib::addPropertyToMesh(surface_mesh, node_to_bulk_node_id_map_name,
512  node_to_bulk_node_id_map);
513  }
514 
515  // transmit the original bulk element ids as a property
516  if (!element_to_bulk_element_id_map_name.empty())
517  {
519  surface_mesh, element_to_bulk_element_id_map_name,
520  MeshLib::MeshItemType::Cell, 1, element_to_bulk_element_id_map);
521  }
522 
523  // transmit the face id of the original bulk element as a property
524  if (!element_to_bulk_face_id_map_name.empty())
525  {
527  surface_mesh, element_to_bulk_face_id_map_name,
528  MeshLib::MeshItemType::Cell, 1, element_to_bulk_face_id_map);
529  }
530 }
531 
532 } // end namespace MeshLib
Definition of Duplicate functions.
Definition of the Line class.
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Definition of the MeshSurfaceExtraction class.
Definition of the Mesh class.
Definition of the Quad class.
Definition of the Tri class.
std::size_t getID() const
Definition: Point3dWithID.h:62
virtual MeshElemType getGeomType() const =0
virtual const Element * getBoundary(unsigned i) const =0
virtual const Element * getNeighbor(unsigned i) const =0
Get the specified neighbor.
virtual unsigned getNumberOfBoundaries() const =0
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
static Eigen::Vector3d getSurfaceNormal(Element const &e)
Returns the surface normal of a 2D element.
Definition: FaceRule.cpp:40
static void get2DSurfaceElements(const std::vector< MeshLib::Element * > &all_elements, std::vector< MeshLib::Element * > &sfc_elements, std::vector< std::size_t > &element_to_bulk_element_id_map, std::vector< std::size_t > &element_to_bulk_face_id_map, Eigen::Vector3d const &dir, double angle, unsigned mesh_dimension)
Functionality needed for getSurfaceNodes() and getMeshSurface()
static MeshLib::Mesh * getMeshSurface(const MeshLib::Mesh &subsfc_mesh, Eigen::Vector3d const &dir, double angle, std::string const &subsfc_node_id_prop_name="", std::string const &subsfc_element_id_prop_name="", std::string const &face_id_prop_name="")
static std::vector< MeshLib::Node * > getSurfaceNodes(const MeshLib::Mesh &mesh, Eigen::Vector3d const &dir, double angle)
Returns the surface nodes of a mesh.
static std::vector< double > getSurfaceAreaForNodes(const MeshLib::Mesh &mesh)
Returns a vector of the areas assigned to each node on a surface mesh.
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:95
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:71
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
const std::string getName() const
Get name of the mesh.
Definition: Mesh.h:92
Properties & getProperties()
Definition: Mesh.h:123
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:89
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:232
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:86
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition: Properties.h:36
std::string const & getPropertyName() const
MeshItemType getMeshItemType() const
PROP_VAL_TYPE & getComponent(std::size_t tuple_index, int component)
Returns the value for the given component stored in the given tuple.
std::unique_ptr< MeshLib::Mesh > getBoundaryElementsAsMesh(MeshLib::Mesh const &bulk_mesh, std::string const &subsfc_node_id_prop_name, std::string const &subsfc_element_id_prop_name, std::string const &face_id_prop_name)
void createSurfaceElementsFromElement(MeshLib::Element const &surface_element, std::vector< MeshLib::Element * > &surface_elements, std::vector< std::size_t > &element_to_bulk_element_id_map, std::vector< std::size_t > &element_to_bulk_face_id_map)
bool createSfcMeshProperties(MeshLib::Mesh &sfc_mesh, MeshLib::Properties const &properties, std::vector< std::size_t > const &node_ids_map, std::vector< std::size_t > const &element_ids_map)
void processPropertyVector(MeshLib::PropertyVector< T > const &property, std::vector< std::size_t > const &id_map, MeshLib::Mesh &sfc_mesh)
std::tuple< std::vector< MeshLib::Node * >, std::vector< std::size_t > > createNodesAndIDMapFromElements(std::vector< MeshLib::Element * > const &elements, std::size_t const n_all_nodes)
void addBulkIDPropertiesToMesh(MeshLib::Mesh &surface_mesh, std::string const &node_to_bulk_node_id_map_name, std::vector< std::size_t > const &node_to_bulk_node_id_map, std::string const &element_to_bulk_element_id_map_name, std::vector< std::size_t > const &element_to_bulk_element_id_map, std::string const &element_to_bulk_face_id_map_name, std::vector< std::size_t > const &element_to_bulk_face_id_map)
void addPropertyToMesh(Mesh &mesh, std::string const &name, MeshItemType item_type, std::size_t number_of_components, std::vector< T > const &values)
Definition: Mesh.h:193
std::vector< Element * > copyElementVector(std::vector< Element * > const &elements, std::vector< Node * > const &new_nodes, std::vector< std::size_t > const *const node_id_map)
std::tuple< std::vector< MeshLib::Element * >, std::vector< std::size_t >, std::vector< std::size_t > > createBoundaryElements(MeshLib::Mesh const &bulk_mesh)