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