OGS
MeshRevision.cpp
Go to the documentation of this file.
1 
15 #include "MeshRevision.h"
16 
17 #include <numeric>
18 
19 #include "BaseLib/Algorithm.h"
20 #include "BaseLib/Logging.h"
22 #include "GeoLib/Grid.h"
25 #include "MeshLib/Mesh.h"
26 
27 namespace MeshLib
28 {
29 const std::array<unsigned, 8> MeshRevision::_hex_diametral_nodes = {
30  {6, 7, 4, 5, 2, 3, 0, 1}};
31 
33 
35 {
36  std::vector<std::size_t> id_map(this->collapseNodeIndices(eps));
37  std::size_t nNodes(id_map.size());
38  unsigned count(0);
39  for (std::size_t i = 0; i < nNodes; ++i)
40  {
41  if (i != id_map[i])
42  {
43  count++;
44  }
45  }
46  return count;
47 }
48 
49 MeshLib::Mesh* MeshRevision::simplifyMesh(const std::string& new_mesh_name,
50  double eps,
51  unsigned min_elem_dim) const
52 {
53  if (this->_mesh.getNumberOfElements() == 0)
54  {
55  return nullptr;
56  }
57 
58  std::vector<MeshLib::Element*> const& elements(this->_mesh.getElements());
59  auto const node_ids = collapseNodeIndices(eps);
60  std::vector<MeshLib::Node*> new_nodes =
61  this->constructNewNodesArray(node_ids);
62  std::vector<MeshLib::Element*> new_elements;
63  std::vector<std::size_t> element_ids;
64 
65  for (std::size_t k(0); k < elements.size(); ++k)
66  {
67  MeshLib::Element const* const elem(elements[k]);
68  unsigned n_unique_nodes(this->getNumberOfUniqueNodes(elem));
69  if (n_unique_nodes == elem->getNumberOfBaseNodes() &&
70  elem->getDimension() >= min_elem_dim)
71  {
72  ElementErrorCode e(elem->validate());
74  {
75  std::size_t const n_new_elements(
76  subdivideElement(elem, new_nodes, new_elements));
77  if (n_new_elements == 0)
78  {
79  ERR("Element {:d} has unknown element type.", k);
81  BaseLib::cleanupVectorElements(new_nodes, new_elements);
82  return nullptr;
83  }
84  element_ids.insert(element_ids.end(), n_new_elements, k);
85  }
86  else
87  {
88  new_elements.push_back(MeshLib::copyElement(elem, new_nodes));
89  element_ids.push_back(k);
90  }
91  }
92  else if (n_unique_nodes < elem->getNumberOfBaseNodes() &&
93  n_unique_nodes > 1)
94  {
95  std::size_t const n_new_elements(reduceElement(
96  elem, n_unique_nodes, new_nodes, new_elements, min_elem_dim));
97  element_ids.insert(element_ids.end(), n_new_elements, k);
98  }
99  else
100  {
101  ERR("Something is wrong, more unique nodes than actual nodes");
102  }
103  }
104 
105  auto const& props = _mesh.getProperties();
106  MeshLib::Properties const new_properties =
107  copyProperties(props, node_ids, element_ids);
108 
110  if (!new_elements.empty())
111  {
112  return new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements,
113  new_properties);
114  }
115 
116  BaseLib::cleanupVectorElements(new_nodes, new_elements);
117  return nullptr;
118 }
119 
120 std::vector<std::size_t> MeshRevision::collapseNodeIndices(double eps) const
121 {
122  const std::vector<MeshLib::Node*>& nodes(_mesh.getNodes());
123  const std::size_t nNodes(_mesh.getNumberOfNodes());
124  std::vector<std::size_t> id_map(nNodes);
125  const double half_eps(eps / 2.0);
126  const double sqr_eps(eps * eps);
127  std::iota(id_map.begin(), id_map.end(), 0);
128 
129  GeoLib::Grid<MeshLib::Node> const grid(nodes.begin(), nodes.end(), 64);
130 
131  for (std::size_t k = 0; k < nNodes; ++k)
132  {
133  MeshLib::Node const* const node(nodes[k]);
134  if (node->getID() != k)
135  {
136  continue;
137  }
138  std::vector<std::vector<MeshLib::Node*> const*> node_vectors(
139  grid.getPntVecsOfGridCellsIntersectingCube(*node, half_eps));
140 
141  const std::size_t nVectors(node_vectors.size());
142  for (std::size_t i = 0; i < nVectors; ++i)
143  {
144  const std::vector<MeshLib::Node*>& cell_vector(*node_vectors[i]);
145  const std::size_t nGridCellNodes(cell_vector.size());
146  for (std::size_t j = 0; j < nGridCellNodes; ++j)
147  {
148  MeshLib::Node const* const test_node(cell_vector[j]);
149  // are node indices already identical (i.e. nodes will be
150  // collapsed)
151  if (id_map[node->getID()] == id_map[test_node->getID()])
152  {
153  continue;
154  }
155 
156  // if test_node has already been collapsed to another node x,
157  // ignore it (if the current node would need to be collapsed
158  // with x it would already have happened when x was tested)
159  if (test_node->getID() != id_map[test_node->getID()])
160  {
161  continue;
162  }
163 
164  // calc distance
165  if (MathLib::sqrDist(*node, *test_node) < sqr_eps)
166  {
167  id_map[test_node->getID()] = node->getID();
168  }
169  }
170  }
171  }
172  return id_map;
173 }
174 
175 std::vector<MeshLib::Node*> MeshRevision::constructNewNodesArray(
176  const std::vector<std::size_t>& id_map) const
177 {
178  const std::vector<MeshLib::Node*>& nodes(_mesh.getNodes());
179  const std::size_t nNodes(nodes.size());
180  std::vector<MeshLib::Node*> new_nodes;
181  new_nodes.reserve(nNodes);
182  for (std::size_t k = 0; k < nNodes; ++k)
183  {
184  // all nodes that have not been collapsed with other nodes are copied
185  // into new array
186  if (nodes[k]->getID() == id_map[k])
187  {
188  std::size_t const id(new_nodes.size());
189  new_nodes.push_back(new MeshLib::Node(
190  (*nodes[k])[0], (*nodes[k])[1], (*nodes[k])[2], id));
191  nodes[k]->setID(id); // the node in the old array gets the index of
192  // the same node in the new array
193  }
194  // the other nodes are not copied and get the index of the nodes they
195  // will have been collapsed with
196  else
197  {
198  nodes[k]->setID(nodes[id_map[k]]->getID());
199  }
200  }
201  return new_nodes;
202 }
203 
205  MeshLib::Element const* const element)
206 {
207  unsigned const nNodes(element->getNumberOfBaseNodes());
208  unsigned count(nNodes);
209 
210  for (unsigned i = 0; i < nNodes - 1; ++i)
211  {
212  for (unsigned j = i + 1; j < nNodes; ++j)
213  {
214  if (element->getNode(i)->getID() == element->getNode(j)->getID())
215  {
216  count--;
217  break;
218  }
219  }
220  }
221  return count;
222 }
223 
224 template <typename T>
225 void fillNodeProperty(std::vector<T>& new_prop,
226  std::vector<T> const& old_prop,
227  std::vector<size_t>
228  node_ids)
229 {
230  std::size_t const n_nodes = node_ids.size();
231  for (std::size_t i = 0; i < n_nodes; ++i)
232  {
233  if (node_ids[i] != i)
234  {
235  continue;
236  }
237  new_prop.push_back(old_prop[i]);
238  }
239 }
240 
241 template <typename T>
242 void fillElemProperty(std::vector<T>& new_prop,
243  std::vector<T> const& old_prop,
244  std::vector<size_t>
245  elem_ids)
246 {
247  std::transform(elem_ids.cbegin(), elem_ids.cend(),
248  std::back_inserter(new_prop),
249  [&](std::size_t const i) { return old_prop[i]; });
250 }
251 
253  MeshLib::Properties const& props,
254  std::vector<std::size_t> const& node_ids,
255  std::vector<std::size_t> const& elem_ids) const
256 {
257  auto const prop_names = props.getPropertyVectorNames();
258  MeshLib::Properties new_properties;
259 
260  for (auto name : prop_names)
261  {
262  if (props.existsPropertyVector<int>(name, MeshItemType::Node, 1))
263  {
264  auto p = props.getPropertyVector<int>(name, MeshItemType::Node, 1);
265  auto new_node_vec = new_properties.createNewPropertyVector<int>(
267  fillNodeProperty(*new_node_vec, *p, node_ids);
268  continue;
269  }
270  if (props.existsPropertyVector<float>(name, MeshItemType::Node, 1))
271  {
272  auto p =
273  props.getPropertyVector<float>(name, MeshItemType::Node, 1);
274  auto new_node_vec = new_properties.createNewPropertyVector<float>(
276  fillNodeProperty(*new_node_vec, *p, node_ids);
277  continue;
278  }
279  if (props.existsPropertyVector<double>(name, MeshItemType::Node, 1))
280  {
281  auto p =
282  props.getPropertyVector<double>(name, MeshItemType::Node, 1);
283  auto new_node_vec = new_properties.createNewPropertyVector<double>(
285  fillNodeProperty(*new_node_vec, *p, node_ids);
286  continue;
287  }
288  if (props.existsPropertyVector<int>(name, MeshItemType::Cell, 1))
289  {
290  auto p = props.getPropertyVector<int>(name, MeshItemType::Cell, 1);
291  auto new_cell_vec = new_properties.createNewPropertyVector<int>(
293  fillElemProperty(*new_cell_vec, *p, elem_ids);
294  continue;
295  }
296  if (props.existsPropertyVector<float>(name, MeshItemType::Cell, 1))
297  {
298  auto p =
299  props.getPropertyVector<float>(name, MeshItemType::Cell, 1);
300  auto new_cell_vec = new_properties.createNewPropertyVector<float>(
302  fillElemProperty(*new_cell_vec, *p, elem_ids);
303  continue;
304  }
305  if (props.existsPropertyVector<double>(name, MeshItemType::Cell, 1))
306  {
307  auto p =
308  props.getPropertyVector<double>(name, MeshItemType::Cell, 1);
309  auto new_cell_vec = new_properties.createNewPropertyVector<double>(
311  fillElemProperty(*new_cell_vec, *p, elem_ids);
312  continue;
313  }
314  WARN("PropertyVector {:s} not being converted.", name);
315  }
316  return new_properties;
317 }
318 
320  MeshLib::Element const* const element,
321  std::vector<MeshLib::Node*> const& nodes,
322  std::vector<MeshLib::Element*>& elements) const
323 {
324  if (element->getGeomType() == MeshElemType::QUAD)
325  {
326  return this->subdivideQuad(element, nodes, elements);
327  }
328  if (element->getGeomType() == MeshElemType::HEXAHEDRON)
329  {
330  return this->subdivideHex(element, nodes, elements);
331  }
332  if (element->getGeomType() == MeshElemType::PYRAMID)
333  {
334  return this->subdividePyramid(element, nodes, elements);
335  }
336  if (element->getGeomType() == MeshElemType::PRISM)
337  {
338  return this->subdividePrism(element, nodes, elements);
339  }
340  return 0;
341 }
342 
344  MeshLib::Element const* const element,
345  unsigned n_unique_nodes,
346  std::vector<MeshLib::Node*> const& nodes,
347  std::vector<MeshLib::Element*>& elements,
348  unsigned min_elem_dim) const
349 {
350  /***************
351  * TODO: modify neighbouring elements if one elements has been subdivided
352  ***************/
353  if (element->getGeomType() == MeshElemType::TRIANGLE && min_elem_dim == 1)
354  {
355  elements.push_back(this->constructLine(element, nodes));
356  return 1;
357  }
358  if ((element->getGeomType() == MeshElemType::QUAD) ||
359  (element->getGeomType() == MeshElemType::TETRAHEDRON))
360  {
361  if (n_unique_nodes == 3 && min_elem_dim < 3)
362  {
363  elements.push_back(this->constructTri(element, nodes));
364  }
365  else if (min_elem_dim == 1)
366  {
367  elements.push_back(this->constructLine(element, nodes));
368  }
369  return 1;
370  }
371  if (element->getGeomType() == MeshElemType::HEXAHEDRON)
372  {
373  return reduceHex(element, n_unique_nodes, nodes, elements,
374  min_elem_dim);
375  }
376  if (element->getGeomType() == MeshElemType::PYRAMID)
377  {
378  this->reducePyramid(element, n_unique_nodes, nodes, elements,
379  min_elem_dim);
380  return 1;
381  }
382  if (element->getGeomType() == MeshElemType::PRISM)
383  {
384  return reducePrism(element, n_unique_nodes, nodes, elements,
385  min_elem_dim);
386  }
387 
388  ERR("Unknown element type.");
389  return 0;
390 }
391 
393  MeshLib::Element const* const quad,
394  std::vector<MeshLib::Node*> const& nodes,
395  std::vector<MeshLib::Element*>& new_elements) const
396 {
397  std::array tri1_nodes{nodes[quad->getNode(0)->getID()],
398  nodes[quad->getNode(1)->getID()],
399  nodes[quad->getNode(2)->getID()]};
400  new_elements.push_back(new MeshLib::Tri(tri1_nodes));
401 
402  std::array tri2_nodes{nodes[quad->getNode(0)->getID()],
403  nodes[quad->getNode(2)->getID()],
404  nodes[quad->getNode(3)->getID()]};
405  new_elements.push_back(new MeshLib::Tri(tri2_nodes));
406 
407  return 2;
408 }
409 
411  MeshLib::Element const* const hex,
412  std::vector<MeshLib::Node*> const& nodes,
413  std::vector<MeshLib::Element*>& new_elements) const
414 {
415  std::array<Node*, 6> prism1_nodes;
416  prism1_nodes[0] = nodes[hex->getNode(0)->getID()];
417  prism1_nodes[1] = nodes[hex->getNode(2)->getID()];
418  prism1_nodes[2] = nodes[hex->getNode(1)->getID()];
419  prism1_nodes[3] = nodes[hex->getNode(4)->getID()];
420  prism1_nodes[4] = nodes[hex->getNode(6)->getID()];
421  prism1_nodes[5] = nodes[hex->getNode(5)->getID()];
422  auto* prism1(new MeshLib::Prism(prism1_nodes));
423  this->subdividePrism(prism1, nodes, new_elements);
424  delete prism1;
425 
426  std::array<Node*, 6> prism2_nodes;
427  prism2_nodes[0] = nodes[hex->getNode(4)->getID()];
428  prism2_nodes[1] = nodes[hex->getNode(6)->getID()];
429  prism2_nodes[2] = nodes[hex->getNode(7)->getID()];
430  prism2_nodes[3] = nodes[hex->getNode(0)->getID()];
431  prism2_nodes[4] = nodes[hex->getNode(2)->getID()];
432  prism2_nodes[5] = nodes[hex->getNode(3)->getID()];
433  auto* prism2(new MeshLib::Prism(prism2_nodes));
434  this->subdividePrism(prism2, nodes, new_elements);
435  delete prism2;
436 
437  return 6;
438 }
439 
441  MeshLib::Element const* const pyramid,
442  std::vector<MeshLib::Node*> const& nodes,
443  std::vector<MeshLib::Element*>& new_elements) const
444 {
445  auto addTetrahedron =
446  [&pyramid, &nodes, &new_elements](std::size_t id0, std::size_t id1,
447  std::size_t id2, std::size_t id3)
448  {
449  std::array<Node*, 4> tet_nodes;
450  tet_nodes[0] = nodes[pyramid->getNode(id0)->getID()];
451  tet_nodes[1] = nodes[pyramid->getNode(id1)->getID()];
452  tet_nodes[2] = nodes[pyramid->getNode(id2)->getID()];
453  tet_nodes[3] = nodes[pyramid->getNode(id3)->getID()];
454  new_elements.push_back(new MeshLib::Tet(tet_nodes));
455  };
456 
457  addTetrahedron(0, 1, 2, 4);
458 
459  addTetrahedron(0, 2, 3, 4);
460 
461  return 2;
462 }
463 
465  MeshLib::Element const* const prism,
466  std::vector<MeshLib::Node*> const& nodes,
467  std::vector<MeshLib::Element*>& new_elements) const
468 {
469  auto addTetrahedron =
470  [&prism, &nodes, &new_elements](std::size_t id0, std::size_t id1,
471  std::size_t id2, std::size_t id3)
472  {
473  std::array<Node*, 4> tet_nodes;
474  tet_nodes[0] = nodes[prism->getNode(id0)->getID()];
475  tet_nodes[1] = nodes[prism->getNode(id1)->getID()];
476  tet_nodes[2] = nodes[prism->getNode(id2)->getID()];
477  tet_nodes[3] = nodes[prism->getNode(id3)->getID()];
478  new_elements.push_back(new MeshLib::Tet(tet_nodes));
479  };
480 
481  addTetrahedron(0, 1, 2, 3);
482 
483  addTetrahedron(3, 2, 4, 5);
484 
485  addTetrahedron(2, 1, 3, 4);
486 
487  return 3;
488 }
489 
490 unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem,
491  unsigned n_unique_nodes,
492  std::vector<MeshLib::Node*> const& nodes,
493  std::vector<MeshLib::Element*>& new_elements,
494  unsigned min_elem_dim) const
495 {
496  // TODO?
497  // if two diametral nodes collapse, all kinds of bizarre (2D-)element
498  // combinations could be the result. this case is currently not implemented!
499 
500  if (n_unique_nodes == 7)
501  {
502  // reduce to prism + pyramid
503  for (unsigned i = 0; i < 7; ++i)
504  {
505  for (unsigned j = i + 1; j < 8; ++j)
506  {
507  if (org_elem->getNode(i)->getID() ==
508  org_elem->getNode(j)->getID())
509  {
510  const std::array<unsigned, 4> base_nodes(
511  this->lutHexCuttingQuadNodes(i, j));
512  std::array pyr_nodes{
513  nodes[org_elem->getNode(base_nodes[0])->getID()],
514  nodes[org_elem->getNode(base_nodes[1])->getID()],
515  nodes[org_elem->getNode(base_nodes[2])->getID()],
516  nodes[org_elem->getNode(base_nodes[3])->getID()],
517  nodes[org_elem->getNode(i)->getID()]};
518  new_elements.push_back(new MeshLib::Pyramid(pyr_nodes));
519 
520  if (i < 4 && j >= 4)
521  {
522  std::swap(i, j);
523  }
524  std::array prism_nodes{
525  nodes[org_elem->getNode(base_nodes[0])->getID()],
526  nodes[org_elem->getNode(base_nodes[3])->getID()],
527  nodes[org_elem->getNode(this->lutHexDiametralNode(j))
528  ->getID()],
529  nodes[org_elem->getNode(base_nodes[1])->getID()],
530  nodes[org_elem->getNode(base_nodes[2])->getID()],
531  nodes[org_elem->getNode(this->lutHexDiametralNode(i))
532  ->getID()]};
533  new_elements.push_back(new MeshLib::Prism(prism_nodes));
534  return 2;
535  }
536  }
537  }
538  }
539  else if (n_unique_nodes == 6)
540  {
541  // reduce to prism
542  for (unsigned i = 0; i < 6; ++i)
543  {
544  const MeshLib::Element* face(org_elem->getFace(i));
545  if (face->getNode(0)->getID() == face->getNode(1)->getID() &&
546  face->getNode(2)->getID() == face->getNode(3)->getID())
547  {
548  std::array prism_nodes{
549  nodes[org_elem
550  ->getNode(
552  *org_elem, face->getNode(0))))
553  ->getID()],
554  nodes[org_elem
555  ->getNode(
557  *org_elem, face->getNode(1))))
558  ->getID()],
559  nodes[org_elem
560  ->getNode(getNodeIDinElement(*org_elem,
561  face->getNode(2)))
562  ->getID()],
563  nodes[org_elem
564  ->getNode(
566  *org_elem, face->getNode(2))))
567  ->getID()],
568  nodes[org_elem
569  ->getNode(
571  *org_elem, face->getNode(3))))
572  ->getID()],
573  nodes[org_elem
574  ->getNode(getNodeIDinElement(*org_elem,
575  face->getNode(0)))
576  ->getID()]};
577  new_elements.push_back(new MeshLib::Prism(prism_nodes));
578  delete face;
579  return 1;
580  }
581  if (face->getNode(0)->getID() == face->getNode(3)->getID() &&
582  face->getNode(1)->getID() == face->getNode(2)->getID())
583  {
584  std::array prism_nodes{
585  nodes[org_elem
586  ->getNode(
588  *org_elem, face->getNode(0))))
589  ->getID()],
590  nodes[org_elem
591  ->getNode(
593  *org_elem, face->getNode(3))))
594  ->getID()],
595  nodes[org_elem
596  ->getNode(getNodeIDinElement(*org_elem,
597  face->getNode(2)))
598  ->getID()],
599  nodes[org_elem
600  ->getNode(
602  *org_elem, face->getNode(1))))
603  ->getID()],
604  nodes[org_elem
605  ->getNode(
607  *org_elem, face->getNode(2))))
608  ->getID()],
609  nodes[org_elem
610  ->getNode(getNodeIDinElement(*org_elem,
611  face->getNode(0)))
612  ->getID()]};
613  new_elements.push_back(new MeshLib::Prism(prism_nodes));
614  delete face;
615  return 1;
616  }
617  delete face;
618  }
619  // reduce to four tets -> divide into 2 prisms such that each has one
620  // collapsed node
621  for (unsigned i = 0; i < 7; ++i)
622  {
623  for (unsigned j = i + 1; j < 8; ++j)
624  {
625  if (org_elem->getNode(i)->getID() ==
626  org_elem->getNode(j)->getID())
627  {
628  for (unsigned k = i; k < 7; ++k)
629  {
630  for (unsigned l = k + 1; l < 8; ++l)
631  {
632  if (!(i == k && j == l) && org_elem->isEdge(i, j) &&
633  org_elem->isEdge(k, l) &&
634  org_elem->getNode(k)->getID() ==
635  org_elem->getNode(l)->getID())
636  {
637  const std::pair<unsigned, unsigned> back(
638  this->lutHexBackNodes(i, j, k, l));
639  if (back.first ==
640  std::numeric_limits<unsigned>::max() ||
641  back.second ==
642  std::numeric_limits<unsigned>::max())
643  {
644  ERR("Unexpected error during Hex "
645  "reduction");
646  return 0;
647  }
648 
649  std::array<unsigned, 4> cutting_plane(
650  this->lutHexCuttingQuadNodes(back.first,
651  back.second));
652  std::array pris1_nodes{
653  const_cast<MeshLib::Node*>(
654  org_elem->getNode(back.first)),
655  const_cast<MeshLib::Node*>(
656  org_elem->getNode(cutting_plane[0])),
657  const_cast<MeshLib::Node*>(
658  org_elem->getNode(cutting_plane[3])),
659  const_cast<MeshLib::Node*>(
660  org_elem->getNode(back.second)),
661  const_cast<MeshLib::Node*>(
662  org_elem->getNode(cutting_plane[1])),
663  const_cast<MeshLib::Node*>(
664  org_elem->getNode(cutting_plane[2]))};
665  auto* prism1(new MeshLib::Prism(pris1_nodes));
666  unsigned nNewElements = this->reducePrism(
667  prism1, 5, nodes, new_elements,
668  min_elem_dim);
669  delete prism1;
670 
671  std::array pris2_nodes{
672  const_cast<MeshLib::Node*>(
673  org_elem->getNode(
674  this->lutHexDiametralNode(
675  back.first))),
676  const_cast<MeshLib::Node*>(
677  org_elem->getNode(cutting_plane[0])),
678  const_cast<MeshLib::Node*>(
679  org_elem->getNode(cutting_plane[3])),
680  const_cast<MeshLib::Node*>(
681  org_elem->getNode(
682  this->lutHexDiametralNode(
683  back.second))),
684  const_cast<MeshLib::Node*>(
685  org_elem->getNode(cutting_plane[1])),
686  const_cast<MeshLib::Node*>(
687  org_elem->getNode(cutting_plane[2]))};
688  auto* prism2(new MeshLib::Prism(pris2_nodes));
689  nNewElements += this->reducePrism(
690  prism2, 5, nodes, new_elements,
691  min_elem_dim);
692  delete prism2;
693  return nNewElements;
694  }
695  }
696  }
697  }
698  }
699  }
700  }
701  else if (n_unique_nodes == 5)
702  {
703  MeshLib::Element* tet1(this->constructFourNodeElement(org_elem, nodes));
704  std::array<std::size_t, 4> first_four_nodes = {
705  {tet1->getNode(0)->getID(), tet1->getNode(1)->getID(),
706  tet1->getNode(2)->getID(), tet1->getNode(3)->getID()}};
707  unsigned fifth_node(
708  this->findPyramidTopNode(*org_elem, first_four_nodes));
709 
710  bool tet_changed(false);
711  if (tet1->getGeomType() == MeshElemType::QUAD)
712  {
713  delete tet1;
714  tet_changed = true;
715  auto** tet1_nodes = new MeshLib::Node*[4];
716  tet1_nodes[0] = nodes[first_four_nodes[0]];
717  tet1_nodes[1] = nodes[first_four_nodes[1]];
718  tet1_nodes[2] = nodes[first_four_nodes[2]];
719  tet1_nodes[3] = nodes[org_elem->getNode(fifth_node)->getID()];
720  new_elements.push_back(new MeshLib::Tet(tet1_nodes));
721  }
722  else
723  {
724  new_elements.push_back(tet1);
725  }
726 
727  std::array tet2_nodes = {(tet_changed) ? nodes[first_four_nodes[0]]
728  : nodes[first_four_nodes[1]],
729  nodes[first_four_nodes[2]],
730  nodes[first_four_nodes[3]],
731  nodes[org_elem->getNode(fifth_node)->getID()]};
732  new_elements.push_back(new MeshLib::Tet(tet2_nodes));
733  return 2;
734  }
735  else if (n_unique_nodes == 4)
736  {
737  MeshLib::Element* elem(
738  this->constructFourNodeElement(org_elem, nodes, min_elem_dim));
739  if (elem)
740  {
741  new_elements.push_back(elem);
742  return 1;
743  }
744  }
745  else if (n_unique_nodes == 3 && min_elem_dim < 3)
746  {
747  new_elements.push_back(this->constructTri(org_elem, nodes));
748  return 1;
749  }
750  else if (min_elem_dim == 1)
751  {
752  new_elements.push_back(this->constructLine(org_elem, nodes));
753  return 1;
754  }
755  return 0;
756 }
757 
759  unsigned n_unique_nodes,
760  std::vector<MeshLib::Node*> const& nodes,
761  std::vector<MeshLib::Element*>& new_elements,
762  unsigned min_elem_dim) const
763 {
764  if (n_unique_nodes == 4)
765  {
766  MeshLib::Element* elem(
767  this->constructFourNodeElement(org_elem, nodes, min_elem_dim));
768  if (elem)
769  {
770  new_elements.push_back(elem);
771  }
772  }
773  else if (n_unique_nodes == 3 && min_elem_dim < 3)
774  {
775  new_elements.push_back(this->constructTri(org_elem, nodes));
776  }
777  else if (n_unique_nodes == 2 && min_elem_dim == 1)
778  {
779  new_elements.push_back(this->constructLine(org_elem, nodes));
780  }
781 }
782 
783 unsigned MeshRevision::reducePrism(MeshLib::Element const* const org_elem,
784  unsigned n_unique_nodes,
785  std::vector<MeshLib::Node*> const& nodes,
786  std::vector<MeshLib::Element*>& new_elements,
787  unsigned min_elem_dim) const
788 {
789  auto addTetrahedron =
790  [&org_elem, &nodes, &new_elements](std::size_t id0, std::size_t id1,
791  std::size_t id2, std::size_t id3)
792  {
793  std::array tet_nodes{nodes[org_elem->getNode(id0)->getID()],
794  nodes[org_elem->getNode(id1)->getID()],
795  nodes[org_elem->getNode(id2)->getID()],
796  nodes[org_elem->getNode(id3)->getID()]};
797  new_elements.push_back(new MeshLib::Tet(tet_nodes));
798  };
799 
800  // TODO?
801  // In theory a node from the bottom triangle and a node from the top
802  // triangle that are not connected by an edge could collapse, resulting in a
803  // combination of tri and quad elements. This case is currently not tested.
804 
805  // if one of the non-triangle edges collapsed, elem can be reduced to a
806  // pyramid, otherwise it will be two tets
807  if (n_unique_nodes == 5)
808  {
809  for (unsigned i = 0; i < 5; ++i)
810  {
811  for (unsigned j = i + 1; j < 6; ++j)
812  {
813  if (i != j && org_elem->getNode(i)->getID() ==
814  org_elem->getNode(j)->getID())
815  {
816  // non triangle edge collapsed
817  if (i % 3 == j % 3)
818  {
819  addTetrahedron((i + 1) % 3, (i + 2) % 3, i,
820  (i + 1) % 3 + 3);
821  addTetrahedron((i + 1) % 3 + 3, (i + 2) % 3, i,
822  (i + 2) % 3 + 3);
823  return 2;
824  }
825 
826  // triangle edge collapsed
827  const unsigned i_offset = (i > 2) ? i - 3 : i + 3;
828  const unsigned j_offset = (i > 2) ? j - 3 : j + 3;
829  const unsigned k = lutPrismThirdNode(i, j);
830  if (k == std::numeric_limits<unsigned>::max())
831  {
832  ERR("Unexpected error during prism reduction.");
833  return 0;
834  }
835  const unsigned k_offset = (i > 2) ? k - 3 : k + 3;
836 
837  addTetrahedron(i_offset, j_offset, k_offset, i);
838 
839  const unsigned l =
840  (MathLib::isCoplanar(*org_elem->getNode(i_offset),
841  *org_elem->getNode(k_offset),
842  *org_elem->getNode(i),
843  *org_elem->getNode(k)))
844  ? j
845  : i;
846  const unsigned l_offset = (i > 2) ? l - 3 : l + 3;
847  addTetrahedron(l_offset, k_offset, i, k);
848  return 2;
849  }
850  }
851  }
852  }
853  else if (n_unique_nodes == 4)
854  {
855  MeshLib::Element* elem(
856  this->constructFourNodeElement(org_elem, nodes, min_elem_dim));
857  if (elem)
858  {
859  new_elements.push_back(elem);
860  }
861  }
862  else if (n_unique_nodes == 3 && min_elem_dim < 3)
863  {
864  new_elements.push_back(this->constructTri(org_elem, nodes));
865  }
866  else if (n_unique_nodes == 2 && min_elem_dim == 1)
867  {
868  new_elements.push_back(this->constructLine(org_elem, nodes));
869  }
870  return 1;
871 }
872 
874  MeshLib::Element const* const element,
875  const std::vector<MeshLib::Node*>& nodes) const
876 {
877  std::array<Node*, 2> line_nodes;
878  line_nodes[0] = nodes[element->getNode(0)->getID()];
879  line_nodes[1] = nullptr;
880  for (unsigned i = 1; i < element->getNumberOfBaseNodes(); ++i)
881  {
882  if (element->getNode(i)->getID() != element->getNode(0)->getID())
883  {
884  line_nodes[1] = nodes[element->getNode(i)->getID()];
885  break;
886  }
887  }
888  assert(line_nodes[1] != nullptr);
889  return new MeshLib::Line(line_nodes);
890 }
891 
893  MeshLib::Element const* const element,
894  const std::vector<MeshLib::Node*>& nodes) const
895 {
896  // TODO?
897  // In theory three unique nodes could also be reduced to two lines e.g. with
898  // a quad where two diametral nodes collapse. This case is currently not
899  // implemented!
900  std::array<Node*, 3> tri_nodes;
901  tri_nodes[0] = nodes[element->getNode(0)->getID()];
902  tri_nodes[2] = nullptr;
903  for (unsigned i = 1; i < element->getNumberOfBaseNodes(); ++i)
904  {
905  if (element->getNode(i)->getID() != tri_nodes[0]->getID())
906  {
907  tri_nodes[1] = nodes[element->getNode(i)->getID()];
908  for (unsigned j = i + 1; j < element->getNumberOfBaseNodes(); ++j)
909  {
910  if (element->getNode(j)->getID() != tri_nodes[1]->getID())
911  {
912  tri_nodes[2] = nodes[element->getNode(j)->getID()];
913  break;
914  }
915  }
916  if (tri_nodes[2])
917  {
918  break;
919  }
920  }
921  }
922  assert(tri_nodes[2] != nullptr);
923  return new MeshLib::Tri(tri_nodes);
924 }
925 
927  MeshLib::Element const* const element,
928  std::vector<MeshLib::Node*> const& nodes,
929  unsigned min_elem_dim) const
930 {
931  std::array<Node*, 4> new_nodes;
932  unsigned count(0);
933  new_nodes[count++] = nodes[element->getNode(0)->getID()];
934  for (unsigned i = 1; i < element->getNumberOfBaseNodes(); ++i)
935  {
936  if (count > 3)
937  {
938  break;
939  }
940  bool unique_node(true);
941  for (unsigned j = 0; j < i; ++j)
942  {
943  if (element->getNode(i)->getID() == element->getNode(j)->getID())
944  {
945  unique_node = false;
946  break;
947  }
948  }
949  if (unique_node)
950  {
951  new_nodes[count++] = nodes[element->getNode(i)->getID()];
952  };
953  }
954 
955  // test if quad or tet
956  const bool isQuad(MathLib::isCoplanar(*new_nodes[0], *new_nodes[1],
957  *new_nodes[2], *new_nodes[3]));
958  if (isQuad && min_elem_dim < 3)
959  {
960  MeshLib::Element* elem(new MeshLib::Quad(new_nodes));
961  for (unsigned i = 1; i < 3; ++i)
962  {
963  if (elem->validate().none())
964  {
965  return elem;
966  }
967 
968  // change node order if not convex
969  MeshLib::Node* tmp = new_nodes[i + 1];
970  new_nodes[i + 1] = new_nodes[i];
971  new_nodes[i] = tmp;
972  }
973  return elem;
974  }
975  if (!isQuad)
976  {
977  return new MeshLib::Tet(new_nodes);
978  }
979  // is quad but min elem dim == 3
980 
981  return nullptr;
982 }
983 
985  MeshLib::Element const& element,
986  std::array<std::size_t, 4> const& base_node_ids)
987 {
988  const std::size_t nNodes(element.getNumberOfBaseNodes());
989  for (std::size_t i = 0; i < nNodes; ++i)
990  {
991  bool top_node = true;
992  for (unsigned j = 0; j < 4; ++j)
993  {
994  if (element.getNode(i)->getID() == base_node_ids[j])
995  {
996  top_node = false;
997  }
998  }
999  if (top_node)
1000  {
1001  return i;
1002  }
1003  }
1004  return std::numeric_limits<unsigned>::max(); // should never be reached if
1005  // called correctly
1006 }
1007 
1009 {
1010  return _hex_diametral_nodes[id];
1011 }
1012 
1013 std::array<unsigned, 4> MeshRevision::lutHexCuttingQuadNodes(unsigned id1,
1014  unsigned id2)
1015 {
1016  std::array<unsigned, 4> nodes{};
1017  if (id1 == 0 && id2 == 1)
1018  {
1019  nodes[0] = 3;
1020  nodes[1] = 2;
1021  nodes[2] = 5;
1022  nodes[3] = 4;
1023  }
1024  else if (id1 == 1 && id2 == 2)
1025  {
1026  nodes[0] = 0;
1027  nodes[1] = 3;
1028  nodes[2] = 6;
1029  nodes[3] = 5;
1030  }
1031  else if (id1 == 2 && id2 == 3)
1032  {
1033  nodes[0] = 1;
1034  nodes[1] = 0;
1035  nodes[2] = 7;
1036  nodes[3] = 6;
1037  }
1038  else if (id1 == 3 && id2 == 0)
1039  {
1040  nodes[0] = 2;
1041  nodes[1] = 1;
1042  nodes[2] = 4;
1043  nodes[3] = 7;
1044  }
1045  else if (id1 == 4 && id2 == 5)
1046  {
1047  nodes[0] = 0;
1048  nodes[1] = 1;
1049  nodes[2] = 6;
1050  nodes[3] = 7;
1051  }
1052  else if (id1 == 5 && id2 == 6)
1053  {
1054  nodes[0] = 1;
1055  nodes[1] = 2;
1056  nodes[2] = 7;
1057  nodes[3] = 4;
1058  }
1059  else if (id1 == 6 && id2 == 7)
1060  {
1061  nodes[0] = 2;
1062  nodes[1] = 3;
1063  nodes[2] = 4;
1064  nodes[3] = 5;
1065  }
1066  else if (id1 == 7 && id2 == 4)
1067  {
1068  nodes[0] = 3;
1069  nodes[1] = 0;
1070  nodes[2] = 5;
1071  nodes[3] = 6;
1072  }
1073  else if (id1 == 0 && id2 == 4)
1074  {
1075  nodes[0] = 3;
1076  nodes[1] = 7;
1077  nodes[2] = 5;
1078  nodes[3] = 1;
1079  }
1080  else if (id1 == 1 && id2 == 5)
1081  {
1082  nodes[0] = 0;
1083  nodes[1] = 4;
1084  nodes[2] = 6;
1085  nodes[3] = 2;
1086  }
1087  else if (id1 == 2 && id2 == 6)
1088  {
1089  nodes[0] = 1;
1090  nodes[1] = 5;
1091  nodes[2] = 7;
1092  nodes[3] = 3;
1093  }
1094  else if (id1 == 3 && id2 == 7)
1095  {
1096  nodes[0] = 2;
1097  nodes[1] = 6;
1098  nodes[2] = 4;
1099  nodes[3] = 0;
1100  }
1101 
1102  else if (id1 == 1 && id2 == 0)
1103  {
1104  nodes[0] = 2;
1105  nodes[1] = 3;
1106  nodes[2] = 4;
1107  nodes[3] = 5;
1108  }
1109  else if (id1 == 2 && id2 == 1)
1110  {
1111  nodes[0] = 3;
1112  nodes[1] = 0;
1113  nodes[2] = 5;
1114  nodes[3] = 6;
1115  }
1116  else if (id1 == 3 && id2 == 2)
1117  {
1118  nodes[0] = 0;
1119  nodes[1] = 1;
1120  nodes[2] = 6;
1121  nodes[3] = 7;
1122  }
1123  else if (id1 == 0 && id2 == 3)
1124  {
1125  nodes[0] = 1;
1126  nodes[1] = 2;
1127  nodes[2] = 7;
1128  nodes[3] = 4;
1129  }
1130  else if (id1 == 5 && id2 == 4)
1131  {
1132  nodes[0] = 1;
1133  nodes[1] = 0;
1134  nodes[2] = 7;
1135  nodes[3] = 6;
1136  }
1137  else if (id1 == 6 && id2 == 5)
1138  {
1139  nodes[0] = 2;
1140  nodes[1] = 1;
1141  nodes[2] = 4;
1142  nodes[3] = 7;
1143  }
1144  else if (id1 == 7 && id2 == 6)
1145  {
1146  nodes[0] = 3;
1147  nodes[1] = 2;
1148  nodes[2] = 5;
1149  nodes[3] = 4;
1150  }
1151  else if (id1 == 4 && id2 == 7)
1152  {
1153  nodes[0] = 0;
1154  nodes[1] = 3;
1155  nodes[2] = 6;
1156  nodes[3] = 5;
1157  }
1158  else if (id1 == 4 && id2 == 0)
1159  {
1160  nodes[0] = 7;
1161  nodes[1] = 3;
1162  nodes[2] = 1;
1163  nodes[3] = 5;
1164  }
1165  else if (id1 == 5 && id2 == 1)
1166  {
1167  nodes[0] = 4;
1168  nodes[1] = 0;
1169  nodes[2] = 2;
1170  nodes[3] = 6;
1171  }
1172  else if (id1 == 6 && id2 == 2)
1173  {
1174  nodes[0] = 5;
1175  nodes[1] = 1;
1176  nodes[2] = 3;
1177  nodes[3] = 7;
1178  }
1179  else if (id1 == 7 && id2 == 3)
1180  {
1181  nodes[0] = 6;
1182  nodes[1] = 2;
1183  nodes[2] = 0;
1184  nodes[3] = 4;
1185  }
1186  return nodes;
1187 }
1188 
1189 std::pair<unsigned, unsigned> MeshRevision::lutHexBackNodes(unsigned i,
1190  unsigned j,
1191  unsigned k,
1192  unsigned l)
1193 {
1194  // collapsed edges are *not* connected
1195  std::pair<unsigned, unsigned> back(std::numeric_limits<unsigned>::max(),
1196  std::numeric_limits<unsigned>::max());
1197  if (lutHexDiametralNode(i) == k)
1198  {
1199  back.first = i;
1200  back.second = lutHexDiametralNode(l);
1201  }
1202  else if (lutHexDiametralNode(i) == l)
1203  {
1204  back.first = i;
1205  back.second = lutHexDiametralNode(k);
1206  }
1207  else if (lutHexDiametralNode(j) == k)
1208  {
1209  back.first = j;
1210  back.second = lutHexDiametralNode(l);
1211  }
1212  else if (lutHexDiametralNode(j) == l)
1213  {
1214  back.first = j;
1215  back.second = lutHexDiametralNode(k);
1216  }
1217  // collapsed edges *are* connected
1218  else if (i == k)
1219  {
1220  back.first = lutHexDiametralNode(l);
1221  back.second = j;
1222  }
1223  else if (i == l)
1224  {
1225  back.first = lutHexDiametralNode(k);
1226  back.second = j;
1227  }
1228  else if (j == k)
1229  {
1230  back.first = lutHexDiametralNode(l);
1231  back.second = i;
1232  }
1233  else if (j == l)
1234  {
1235  back.first = lutHexDiametralNode(k);
1236  back.second = i;
1237  }
1238  return back;
1239 }
1240 
1241 unsigned MeshRevision::lutPrismThirdNode(unsigned id1, unsigned id2)
1242 {
1243  if ((id1 == 0 && id2 == 1) || (id1 == 1 && id2 == 2))
1244  {
1245  return 2;
1246  }
1247  if ((id1 == 1 && id2 == 2) || (id1 == 2 && id2 == 1))
1248  {
1249  return 0;
1250  }
1251  if ((id1 == 0 && id2 == 2) || (id1 == 2 && id2 == 0))
1252  {
1253  return 1;
1254  }
1255  if ((id1 == 3 && id2 == 4) || (id1 == 4 && id2 == 3))
1256  {
1257  return 5;
1258  }
1259  if ((id1 == 4 && id2 == 5) || (id1 == 5 && id2 == 4))
1260  {
1261  return 3;
1262  }
1263  if ((id1 == 3 && id2 == 5) || (id1 == 5 && id2 == 3))
1264  {
1265  return 4;
1266  }
1267  return std::numeric_limits<unsigned>::max();
1268 }
1269 } // end namespace MeshLib
Definition of Duplicate functions.
Definition of the Grid class.
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 MeshRevision class.
Definition of the Mesh class.
Collects error flags for mesh elements.
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCube(P const &center, double half_len) const
Definition: Grid.h:252
std::size_t getID() const
Definition: Point3dWithID.h:62
virtual MeshElemType getGeomType() const =0
virtual const Node * getNode(unsigned idx) const =0
virtual ElementErrorCode validate() const =0
virtual unsigned getNumberOfBaseNodes() const =0
virtual const Element * getFace(unsigned i) const =0
Returns the i-th face of the element.
virtual bool isEdge(unsigned i, unsigned j) const =0
Returns true if these two indices form an edge and false otherwise.
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
MeshLib::Element * constructTri(MeshLib::Element const *const element, const std::vector< MeshLib::Node * > &nodes) const
Mesh & _mesh
The original mesh used for constructing the class.
Definition: MeshRevision.h:210
MeshLib::Properties copyProperties(MeshLib::Properties const &props, std::vector< std::size_t > const &node_ids, std::vector< std::size_t > const &elem_ids) const
unsigned getNumberOfCollapsableNodes(double eps=std::numeric_limits< double >::epsilon()) const
Returns the number of potentially collapsible nodes.
static unsigned lutHexDiametralNode(unsigned id)
unsigned subdividePrism(MeshLib::Element const *const prism, const std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > &new_elements) const
Subdivides a prism with nonplanar quad faces into two tets.
void reducePyramid(MeshLib::Element const *const org_elem, unsigned n_unique_nodes, const std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned min_elem_dim) const
static unsigned findPyramidTopNode(MeshLib::Element const &element, std::array< std::size_t, 4 > const &base_node_ids)
unsigned subdivideHex(MeshLib::Element const *const hex, const std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > &new_elements) const
Subdivides a Hex with nonplanar faces into tets.
static std::array< unsigned, 4 > lutHexCuttingQuadNodes(unsigned id1, unsigned id2)
unsigned subdivideQuad(MeshLib::Element const *const quad, const std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > &new_elements) const
Subdivides a nonplanar quad into two triangles.
static unsigned lutPrismThirdNode(unsigned id1, unsigned id2)
static std::pair< unsigned, unsigned > lutHexBackNodes(unsigned i, unsigned j, unsigned k, unsigned l)
MeshLib::Element * constructFourNodeElement(MeshLib::Element const *const element, const std::vector< MeshLib::Node * > &nodes, unsigned min_elem_dim=1) const
unsigned subdividePyramid(MeshLib::Element const *const pyramid, const std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > &new_elements) const
Subdivides a pyramid with a nonplanar base into two tets.
std::vector< MeshLib::Node * > constructNewNodesArray(const std::vector< std::size_t > &id_map) const
static unsigned getNumberOfUniqueNodes(MeshLib::Element const *const element)
std::vector< std::size_t > collapseNodeIndices(double eps) const
MeshLib::Mesh * simplifyMesh(const std::string &new_mesh_name, double eps, unsigned min_elem_dim=1) const
unsigned reduceHex(MeshLib::Element const *const org_elem, unsigned n_unique_nodes, const std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned min_elem_dim) const
std::size_t reduceElement(MeshLib::Element const *const element, unsigned n_unique_nodes, const std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > &elements, unsigned min_elem_dim) const
static const std::array< unsigned, 8 > _hex_diametral_nodes
Definition: MeshRevision.h:212
MeshRevision(MeshLib::Mesh &mesh)
unsigned reducePrism(MeshLib::Element const *const org_elem, unsigned n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned min_elem_dim) const
MeshLib::Element * constructLine(MeshLib::Element const *const element, const std::vector< MeshLib::Node * > &nodes) const
std::size_t subdivideElement(MeshLib::Element const *const element, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &elements) const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:95
void resetNodeIDs()
Resets the IDs of all mesh-nodes to their position in the node vector.
Definition: Mesh.cpp:136
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
Properties & getProperties()
Definition: Mesh.h:123
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:89
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::vector< std::string > getPropertyVectorNames() const
Definition: Properties.cpp:35
PropertyVector< T > const * getPropertyVector(std::string const &name) const
bool existsPropertyVector(std::string const &name) const
PropertyVector< T > * createNewPropertyVector(std::string const &name, MeshItemType mesh_item_type, std::size_t n_components=1)
void cleanupVectorElements(std::vector< T1 * > const &items, std::vector< T2 * > const &dependent_items)
Definition: Algorithm.h:300
bool isCoplanar(const MathLib::Point3d &a, const MathLib::Point3d &b, const MathLib::Point3d &c, const MathLib::Point3d &d)
Checks if the four given points are located on a plane.
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition: Point3d.h:48
TemplateElement< MeshLib::TetRule4 > Tet
Definition: Tet.h:25
TemplateElement< MeshLib::LineRule2 > Line
Definition: Line.h:25
TemplateElement< MeshLib::TriRule3 > Tri
Definition: Tri.h:26
unsigned getNodeIDinElement(Element const &element, const MeshLib::Node *node)
Returns the position of the given node in the node array of this element.
Definition: Element.cpp:212
void fillNodeProperty(std::vector< T > &new_prop, std::vector< T > const &old_prop, std::vector< size_t > node_ids)
Element * copyElement(Element const *const element, const std::vector< Node * > &nodes, std::vector< std::size_t > const *const id_map)
void fillElemProperty(std::vector< T > &new_prop, std::vector< T > const &old_prop, std::vector< size_t > elem_ids)