18#include <range/v3/algorithm/copy.hpp>
19#include <range/v3/algorithm/transform.hpp>
20#include <range/v3/range/conversion.hpp>
21#include <range/v3/view/transform.hpp>
38 if ((id1 == 0 && id2 == 1) || (id1 == 1 && id2 == 0))
42 if ((id1 == 1 && id2 == 2) || (id1 == 2 && id2 == 1))
46 if ((id1 == 0 && id2 == 2) || (id1 == 2 && id2 == 0))
50 if ((id1 == 3 && id2 == 4) || (id1 == 4 && id2 == 3))
54 if ((id1 == 4 && id2 == 5) || (id1 == 5 && id2 == 4))
58 if ((id1 == 3 && id2 == 5) || (id1 == 5 && id2 == 3))
62 return std::numeric_limits<unsigned>::max();
68template <
typename ElementType>
70 std::span<MeshLib::Node* const>
const element_nodes,
71 std::vector<MeshLib::Node*>
const& nodes,
72 std::array<std::size_t, ElementType::n_all_nodes>
const local_ids)
75 auto lookup_in = [](
auto const& values)
77 return ranges::views::transform([&values](std::size_t
const n)
78 {
return values[n]; });
81 std::array<MeshLib::Node*, ElementType::n_all_nodes> new_nodes;
82 ranges::copy(local_ids | lookup_in(element_nodes) | ids | lookup_in(nodes),
85 return std::make_unique<ElementType>(new_nodes);
90 std::vector<MeshLib::Node*>
const& nodes,
91 std::vector<MeshLib::Element*>& new_elements)
93 std::array<std::size_t, 3>
const tri1_node_ids{0, 1, 2};
94 new_elements.push_back(
95 createElement<MeshLib::Tri>(quad->
nodes(), nodes, tri1_node_ids)
98 std::array<std::size_t, 3>
const tri2_node_ids{0, 2, 3};
99 new_elements.push_back(
100 createElement<MeshLib::Tri>(quad->
nodes(), nodes, tri2_node_ids)
108 std::vector<MeshLib::Node*>
const& nodes,
109 std::vector<MeshLib::Element*>& new_elements)
111 auto addTetrahedron =
112 [&prism, &nodes, &new_elements](std::array<std::size_t, 4>
const ids)
114 new_elements.push_back(
118 addTetrahedron({0, 1, 2, 3});
119 addTetrahedron({3, 2, 4, 5});
120 addTetrahedron({2, 1, 3, 4});
127 std::vector<MeshLib::Node*>
const& nodes,
128 std::vector<MeshLib::Element*>& new_elements)
143 std::vector<MeshLib::Node*>
const& nodes,
144 std::vector<MeshLib::Element*>& new_elements)
146 auto addTetrahedron =
147 [&pyramid, &nodes, &new_elements](std::array<std::size_t, 4>
const ids)
149 new_elements.push_back(
154 addTetrahedron({0, 1, 2, 4});
155 addTetrahedron({0, 2, 3, 4});
163 const std::vector<MeshLib::Node*>& nodes)
165 std::array<std::size_t, 2> line_node_ids = {0, 0};
170 line_node_ids[1] = i;
174 assert(line_node_ids[1] != 0);
182 const std::vector<MeshLib::Node*>& nodes)
188 std::array<MeshLib::Node*, 3> tri_nodes;
190 tri_nodes[2] =
nullptr;
193 if (element->
getNode(i)->
getID() != tri_nodes[0]->getID())
198 if (element->
getNode(j)->
getID() != tri_nodes[1]->getID())
210 assert(tri_nodes[2] !=
nullptr);
218 std::vector<MeshLib::Node*>
const& nodes,
219 unsigned const min_elem_dim = 1)
221 std::array<MeshLib::Node*, 4> new_nodes;
223 new_nodes[count++] = nodes[element->
getNode(0)->
getID()];
230 bool unique_node(
true);
231 for (
unsigned j = 0; j < i; ++j)
241 new_nodes[count++] = nodes[element->
getNode(i)->
getID()];
247 *new_nodes[2], *new_nodes[3]));
248 if (isQuad && min_elem_dim < 3)
251 for (
unsigned i = 1; i < 3; ++i)
260 new_nodes[i + 1] = new_nodes[i];
277 unsigned const n_unique_nodes,
278 std::vector<MeshLib::Node*>
const& nodes,
279 std::vector<MeshLib::Element*>& new_elements,
280 unsigned const min_elem_dim)
282 if (n_unique_nodes == 4)
288 new_elements.push_back(elem);
291 else if (n_unique_nodes == 3 && min_elem_dim < 3)
295 else if (n_unique_nodes == 2 && min_elem_dim == 1)
305 unsigned const n_unique_nodes,
306 std::vector<MeshLib::Node*>
const& nodes,
307 std::vector<MeshLib::Element*>& new_elements,
308 unsigned const min_elem_dim)
310 auto addTetrahedron =
311 [&org_elem, &nodes, &new_elements](std::array<std::size_t, 4>
const ids)
313 new_elements.push_back(
325 if (n_unique_nodes == 5)
327 for (
unsigned i = 0; i < 5; ++i)
329 for (
unsigned j = i + 1; j < 6; ++j)
338 {(i + 1) % 3, (i + 2) % 3, i, (i + 1) % 3 + 3});
340 {(i + 1) % 3 + 3, (i + 2) % 3, i, (i + 2) % 3 + 3});
345 const unsigned i_offset = (i > 2) ? i - 3 : i + 3;
346 const unsigned j_offset = (i > 2) ? j - 3 : j + 3;
348 if (k == std::numeric_limits<unsigned>::max())
350 ERR(
"Unexpected error during prism reduction.");
353 const unsigned k_offset = (i > 2) ? k - 3 : k + 3;
355 addTetrahedron({i_offset, j_offset, k_offset, i});
364 const unsigned l_offset = (i > 2) ? l - 3 : l + 3;
365 addTetrahedron({l_offset, k_offset, i, k});
371 else if (n_unique_nodes == 4)
377 new_elements.push_back(elem);
380 else if (n_unique_nodes == 3 && min_elem_dim < 3)
384 else if (n_unique_nodes == 2 && min_elem_dim == 1)
395 if (id1 == 0 && id2 == 1)
399 if (id1 == 1 && id2 == 2)
403 if (id1 == 2 && id2 == 3)
407 if (id1 == 3 && id2 == 0)
411 if (id1 == 4 && id2 == 5)
415 if (id1 == 5 && id2 == 6)
419 if (id1 == 6 && id2 == 7)
423 if (id1 == 7 && id2 == 4)
427 if (id1 == 0 && id2 == 4)
431 if (id1 == 1 && id2 == 5)
435 if (id1 == 2 && id2 == 6)
439 if (id1 == 3 && id2 == 7)
443 if (id1 == 1 && id2 == 0)
447 if (id1 == 2 && id2 == 1)
451 if (id1 == 3 && id2 == 2)
455 if (id1 == 0 && id2 == 3)
459 if (id1 == 5 && id2 == 4)
463 if (id1 == 6 && id2 == 5)
467 if (id1 == 7 && id2 == 6)
471 if (id1 == 4 && id2 == 7)
475 if (id1 == 4 && id2 == 0)
479 if (id1 == 5 && id2 == 1)
483 if (id1 == 6 && id2 == 2)
487 if (id1 == 7 && id2 == 3)
493 "lutHexCuttingQuadNodes() for nodes {} and {} does not have a valid "
502 constexpr std::array<unsigned, 8> hex_diametral_node_ids = {
503 {6, 7, 4, 5, 2, 3, 0, 1}};
505 return hex_diametral_node_ids[id];
551 return {std::numeric_limits<unsigned>::max(),
552 std::numeric_limits<unsigned>::max()};
558 std::array<std::size_t, 4>
const& base_node_ids)
561 for (std::size_t i = 0; i < nNodes; ++i)
563 bool top_node =
true;
564 for (
unsigned j = 0; j < 4; ++j)
576 return std::numeric_limits<unsigned>::max();
584 unsigned const n_unique_nodes,
585 std::vector<MeshLib::Node*>
const& nodes,
586 std::vector<MeshLib::Element*>& new_elements,
587 unsigned const min_elem_dim)
593 if (n_unique_nodes == 7)
596 for (
unsigned i = 0; i < 7; ++i)
598 for (
unsigned j = i + 1; j < 8; ++j)
603 const std::array<unsigned, 4> base_node_ids(
605 std::array<std::size_t, 5>
const pyr_node_ids = {
606 base_node_ids[0], base_node_ids[1], base_node_ids[2],
607 base_node_ids[3], i};
608 new_elements.push_back(
617 std::array<std::size_t, 6>
const prism_node_ids{
618 base_node_ids[0], base_node_ids[3],
621 new_elements.push_back(
630 else if (n_unique_nodes == 6)
633 for (
unsigned i = 0; i < 6; ++i)
639 std::array<std::size_t, 6>
const prism_node_ids{
641 getNodeIDinElement(*org_elem, face->
getNode(0))),
643 getNodeIDinElement(*org_elem, face->
getNode(1))),
644 getNodeIDinElement(*org_elem, face->
getNode(2)),
646 getNodeIDinElement(*org_elem, face->
getNode(2))),
648 getNodeIDinElement(*org_elem, face->
getNode(3))),
649 getNodeIDinElement(*org_elem, face->
getNode(0))};
651 new_elements.push_back(
661 std::array<std::size_t, 6>
const prism_node_ids{
663 getNodeIDinElement(*org_elem, face->
getNode(0))),
665 getNodeIDinElement(*org_elem, face->
getNode(3))),
666 getNodeIDinElement(*org_elem, face->
getNode(2)),
668 getNodeIDinElement(*org_elem, face->
getNode(1))),
670 getNodeIDinElement(*org_elem, face->
getNode(2))),
671 getNodeIDinElement(*org_elem, face->
getNode(0))};
672 new_elements.push_back(
683 for (
unsigned i = 0; i < 7; ++i)
685 for (
unsigned j = i + 1; j < 8; ++j)
690 for (
unsigned k = i; k < 7; ++k)
692 for (
unsigned l = k + 1; l < 8; ++l)
694 if (!(i == k && j == l) && org_elem->
isEdge(i, j) &&
699 const std::pair<unsigned, unsigned> back(
702 std::numeric_limits<unsigned>::max() ||
704 std::numeric_limits<unsigned>::max())
706 ERR(
"Unexpected error during Hex "
711 std::array<unsigned, 4>
const cutting_plane(
714 std::array<std::size_t, 6>
const pris1_node_ids{
715 back.first, cutting_plane[0],
716 cutting_plane[3], back.second,
717 cutting_plane[1], cutting_plane[2]};
719 org_elem->
nodes(), nodes, pris1_node_ids);
720 unsigned nNewElements =
722 new_elements, min_elem_dim);
724 std::array<std::size_t, 6>
const pris2_node_ids{
732 org_elem->
nodes(), nodes, pris2_node_ids);
735 new_elements, min_elem_dim);
744 else if (n_unique_nodes == 5)
747 std::array<std::size_t, 4>
const first_four_node_ids = {
750 unsigned const fifth_node =
753 bool tet_changed(
false);
758 std::array
const tet1_nodes = {
759 nodes[first_four_node_ids[0]], nodes[first_four_node_ids[1]],
760 nodes[first_four_node_ids[2]],
766 new_elements.push_back(tet1);
769 std::array
const tet2_nodes = {
770 (tet_changed) ? nodes[first_four_node_ids[0]]
771 : nodes[first_four_node_ids[1]],
772 nodes[first_four_node_ids[2]], nodes[first_four_node_ids[3]],
777 else if (n_unique_nodes == 4)
783 new_elements.push_back(elem);
787 else if (n_unique_nodes == 3 && min_elem_dim < 3)
792 else if (min_elem_dim == 1)
808 std::vector<MeshLib::Node*>
const& nodes,
809 std::vector<MeshLib::Element*>& elements)
833 unsigned const n_unique_nodes,
834 std::vector<MeshLib::Node*>
const& nodes,
835 std::vector<MeshLib::Element*>& elements,
836 unsigned const min_elem_dim)
850 if (n_unique_nodes == 3 && min_elem_dim < 3)
854 else if (min_elem_dim == 1)
862 return reduceHex(element, n_unique_nodes, nodes, elements,
867 reducePyramid(element, n_unique_nodes, nodes, elements, min_elem_dim);
872 return reducePrism(element, n_unique_nodes, nodes, elements,
876 ERR(
"Unknown element type.");
884 unsigned count(nNodes);
886 for (
unsigned i = 0; i < nNodes - 1; ++i)
888 for (
unsigned j = i + 1; j < nNodes; ++j)
903 std::vector<size_t>
const& node_ids)
905 std::size_t
const n_nodes = node_ids.size();
906 for (std::size_t i = 0; i < n_nodes; ++i)
908 if (node_ids[i] != i)
919 std::vector<size_t>
const& elem_ids)
921 assert(new_prop.
size() == old_prop.
size());
922 assert(new_prop.
size() == elem_ids.size());
923 ranges::transform(elem_ids, new_prop.
begin(),
924 [&](std::size_t
const i) { return old_prop[i]; });
931 std::vector<std::size_t>
const& node_ids,
932 std::vector<std::size_t>
const& elem_ids)
937 for (
auto name : prop_names)
999 WARN(
"PropertyVector {:s} not being converted.", name);
1001 return new_properties;
1012 std::size_t
const nNodes = id_map.size();
1014 for (std::size_t i = 0; i < nNodes; ++i)
1026 unsigned const min_elem_dim)
const
1035 std::vector<MeshLib::Node*> new_nodes =
1037 std::vector<MeshLib::Element*> new_elements;
1038 std::vector<std::size_t> element_ids;
1040 for (std::size_t k(0); k < elements.size(); ++k)
1043 unsigned const n_unique_nodes(getNumberOfUniqueNodes(elem));
1050 std::size_t
const n_new_elements(
1051 subdivideElement(elem, new_nodes, new_elements));
1052 if (n_new_elements == 0)
1054 ERR(
"Element {:d} has unknown element type.", k);
1059 element_ids.insert(element_ids.end(), n_new_elements, k);
1064 element_ids.push_back(k);
1067 else if (n_unique_nodes < elem->getNumberOfBaseNodes() &&
1070 std::size_t
const n_new_elements(reduceElement(
1071 elem, n_unique_nodes, new_nodes, new_elements, min_elem_dim));
1072 element_ids.insert(element_ids.end(), n_new_elements, k);
1076 ERR(
"Something is wrong, more unique nodes than actual nodes");
1082 copyProperties(props, node_ids, element_ids);
1085 if (!new_elements.empty())
1087 return new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements,
1097 double const eps)
const
1101 std::vector<std::size_t> id_map(nNodes);
1102 const double half_eps(eps / 2.0);
1103 const double sqr_eps(eps * eps);
1104 std::iota(id_map.begin(), id_map.end(), 0);
1108 for (std::size_t k = 0; k < nNodes; ++k)
1111 if (node->
getID() != k)
1115 std::vector<std::vector<MeshLib::Node*>
const*>
const node_vectors(
1118 const std::size_t nVectors(node_vectors.size());
1119 for (std::size_t i = 0; i < nVectors; ++i)
1121 const std::vector<MeshLib::Node*>& cell_vector(*node_vectors[i]);
1122 const std::size_t nGridCellNodes(cell_vector.size());
1123 for (std::size_t j = 0; j < nGridCellNodes; ++j)
1128 if (id_map[node->
getID()] == id_map[test_node->
getID()])
1136 if (test_node->
getID() != id_map[test_node->
getID()])
1153 const std::vector<std::size_t>& id_map)
const
1156 const std::size_t nNodes(nodes.size());
1157 std::vector<MeshLib::Node*> new_nodes;
1158 new_nodes.reserve(nNodes);
1159 for (std::size_t k = 0; k < nNodes; ++k)
1163 if (nodes[k]->getID() == id_map[k])
1165 std::size_t
const id(new_nodes.size());
1167 (*nodes[k])[0], (*nodes[k])[1], (*nodes[k])[2],
id));
1168 nodes[k]->setID(
id);
1175 nodes[k]->setID(nodes[id_map[k]]->getID());
Definition of Duplicate functions.
Definition of the Grid class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the class Properties that implements a container of properties.
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 ¢er, double half_len) const
std::size_t getID() const
virtual MeshElemType getGeomType() const =0
virtual ElementErrorCode validate() const =0
virtual const Element * getFace(unsigned i) const =0
Returns the i-th face of the element.
virtual unsigned getNumberOfBaseNodes() const =0
virtual const Node * getNode(unsigned idx) const =0
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.
constexpr std::span< Node *const > nodes() const
Span of element's nodes, their pointers actually.
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
void resetNodeIDs()
Resets the IDs of all mesh-nodes to their position in the node vector.
Properties & getProperties()
std::size_t getNumberOfNodes() const
Get the number of nodes.
std::size_t getNumberOfElements() const
Get the number of elements.
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
std::vector< std::string > getPropertyVectorNames() const
bool existsPropertyVector(std::string_view name) const
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
PropertyVector< T > const * getPropertyVector(std::string_view name) const
constexpr std::size_t size() const
constexpr void push_back(const PROP_VAL_TYPE &value)
constexpr PROP_VAL_TYPE * begin()
void cleanupVectorElements(std::vector< T * > &items)
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)
MeshLib specific, lazy, non-owning, non-mutating, composable range views.
TemplateElement< MeshLib::TetRule4 > Tet
TemplateElement< MeshLib::TriRule3 > Tri
Element * copyElement(Element const *const element, const std::vector< Node * > &nodes, std::vector< std::size_t > const *const id_map)
std::size_t reduceElement(MeshLib::Element const *const element, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &elements, unsigned const min_elem_dim)
unsigned lutHexDiametralNode(unsigned const id)
MeshLib::Element * constructTri(MeshLib::Element const *const element, const std::vector< MeshLib::Node * > &nodes)
MeshLib::Element * constructLine(MeshLib::Element const *const element, const std::vector< MeshLib::Node * > &nodes)
unsigned reduceHex(MeshLib::Element const *const org_elem, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned const min_elem_dim)
unsigned reducePrism(MeshLib::Element const *const org_elem, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned const min_elem_dim)
void fillNodeProperty(MeshLib::PropertyVector< T > &new_prop, MeshLib::PropertyVector< T > const &old_prop, std::vector< size_t > const &node_ids)
unsigned getNumberOfUniqueNodes(MeshLib::Element const *const element)
unsigned subdivideHex(MeshLib::Element const *const hex, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a Hex with nonplanar faces into tets.
std::unique_ptr< MeshLib::Element > createElement(std::span< MeshLib::Node *const > const element_nodes, std::vector< MeshLib::Node * > const &nodes, std::array< std::size_t, ElementType::n_all_nodes > const local_ids)
std::size_t subdivideElement(MeshLib::Element const *const element, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &elements)
void reducePyramid(MeshLib::Element const *const org_elem, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned const min_elem_dim)
MeshLib::Properties copyProperties(MeshLib::Properties const &props, std::vector< std::size_t > const &node_ids, std::vector< std::size_t > const &elem_ids)
unsigned subdividePrism(MeshLib::Element const *const prism, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a prism with nonplanar quad faces into two tets.
std::array< unsigned, 4 > lutHexCuttingQuadNodes(unsigned id1, unsigned id2)
std::pair< unsigned, unsigned > lutHexBackNodes(unsigned const i, unsigned const j, unsigned const k, unsigned const l)
unsigned subdividePyramid(MeshLib::Element const *const pyramid, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a pyramid with a nonplanar base into two tets.
void fillElemProperty(MeshLib::PropertyVector< T > &new_prop, MeshLib::PropertyVector< T > const &old_prop, std::vector< size_t > const &elem_ids)
unsigned findPyramidTopNode(MeshLib::Element const &element, std::array< std::size_t, 4 > const &base_node_ids)
MeshLib::Element * constructFourNodeElement(MeshLib::Element const *const element, std::vector< MeshLib::Node * > const &nodes, unsigned const min_elem_dim=1)
unsigned subdivideQuad(MeshLib::Element const *const quad, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a nonplanar quad into two triangles.