diff --git a/MeshToolsLib/MeshEditing/MeshRevision.cpp b/MeshToolsLib/MeshEditing/MeshRevision.cpp index fa46df8476c..315d4fb429f 100644 --- a/MeshToolsLib/MeshEditing/MeshRevision.cpp +++ b/MeshToolsLib/MeshEditing/MeshRevision.cpp @@ -22,8 +22,43 @@ #include "MathLib/GeometricBasics.h" #include "MeshLib/Elements/Elements.h" #include "MeshLib/Mesh.h" +#include "MeshLib/Properties.h" #include "MeshLib/Utils/DuplicateMeshComponents.h" +namespace MeshToolsLib +{ +/// Lookup-table for returning the third node of bottom or top triangle given +/// the other two. +unsigned lutPrismThirdNode(unsigned id1, unsigned id2) +{ + if ((id1 == 0 && id2 == 1) || (id1 == 1 && id2 == 0)) + { + return 2; + } + if ((id1 == 1 && id2 == 2) || (id1 == 2 && id2 == 1)) + { + return 0; + } + if ((id1 == 0 && id2 == 2) || (id1 == 2 && id2 == 0)) + { + return 1; + } + if ((id1 == 3 && id2 == 4) || (id1 == 4 && id2 == 3)) + { + return 5; + } + if ((id1 == 4 && id2 == 5) || (id1 == 5 && id2 == 4)) + { + return 3; + } + if ((id1 == 3 && id2 == 5) || (id1 == 5 && id2 == 3)) + { + return 4; + } + return std::numeric_limits::max(); +} +} // namespace MeshToolsLib + namespace { /// Subdivides a nonplanar quad into two triangles. @@ -238,432 +273,446 @@ MeshLib::Element* constructFourNodeElement( return nullptr; } -} // namespace - -namespace MeshToolsLib -{ -const std::array MeshRevision::_hex_diametral_nodes = { - {6, 7, 4, 5, 2, 3, 0, 1}}; - -MeshRevision::MeshRevision(MeshLib::Mesh& mesh) : _mesh(mesh) {} -unsigned MeshRevision::getNumberOfCollapsibleNodes(double eps) const +/// Reduces a pyramid element by removing collapsed nodes and constructing a new +/// elements from the remaining nodes. +void reducePyramid(MeshLib::Element const* const org_elem, + unsigned n_unique_nodes, + std::vector const& nodes, + std::vector& new_elements, + unsigned min_elem_dim) { - std::vector id_map(this->collapseNodeIndices(eps)); - std::size_t nNodes(id_map.size()); - unsigned count(0); - for (std::size_t i = 0; i < nNodes; ++i) + if (n_unique_nodes == 4) { - if (i != id_map[i]) + MeshLib::Element* elem( + constructFourNodeElement(org_elem, nodes, min_elem_dim)); + if (elem) { - count++; + new_elements.push_back(elem); } } - return count; -} - -MeshLib::Mesh* MeshRevision::simplifyMesh(const std::string& new_mesh_name, - double eps, - unsigned min_elem_dim) const -{ - if (this->_mesh.getNumberOfElements() == 0) - { - return nullptr; - } - - std::vector const& elements(this->_mesh.getElements()); - auto const node_ids = collapseNodeIndices(eps); - std::vector new_nodes = - this->constructNewNodesArray(node_ids); - std::vector new_elements; - std::vector element_ids; - - for (std::size_t k(0); k < elements.size(); ++k) + else if (n_unique_nodes == 3 && min_elem_dim < 3) { - MeshLib::Element const* const elem(elements[k]); - unsigned n_unique_nodes(this->getNumberOfUniqueNodes(elem)); - if (n_unique_nodes == elem->getNumberOfBaseNodes() && - elem->getDimension() >= min_elem_dim) - { - ElementErrorCode e(elem->validate()); - if (e[ElementErrorFlag::NonCoplanar]) - { - std::size_t const n_new_elements( - subdivideElement(elem, new_nodes, new_elements)); - if (n_new_elements == 0) - { - ERR("Element {:d} has unknown element type.", k); - _mesh.resetNodeIDs(); - BaseLib::cleanupVectorElements(new_nodes, new_elements); - return nullptr; - } - element_ids.insert(element_ids.end(), n_new_elements, k); - } - else - { - new_elements.push_back(MeshLib::copyElement(elem, new_nodes)); - element_ids.push_back(k); - } - } - else if (n_unique_nodes < elem->getNumberOfBaseNodes() && - n_unique_nodes > 1) - { - std::size_t const n_new_elements(reduceElement( - elem, n_unique_nodes, new_nodes, new_elements, min_elem_dim)); - element_ids.insert(element_ids.end(), n_new_elements, k); - } - else - { - ERR("Something is wrong, more unique nodes than actual nodes"); - } + new_elements.push_back(constructTri(org_elem, nodes)); } - - auto const& props = _mesh.getProperties(); - MeshLib::Properties const new_properties = - copyProperties(props, node_ids, element_ids); - - _mesh.resetNodeIDs(); - if (!new_elements.empty()) + else if (n_unique_nodes == 2 && min_elem_dim == 1) { - return new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements, - new_properties); + new_elements.push_back(constructLine(org_elem, nodes)); } - - BaseLib::cleanupVectorElements(new_nodes, new_elements); - return nullptr; } -std::vector MeshRevision::collapseNodeIndices(double eps) const +/// Reduces a prism element by removing collapsed nodes and constructing one or +/// two new elements from the remaining nodes. +/// @return The number of newly created elements +unsigned reducePrism(MeshLib::Element const* const org_elem, + unsigned n_unique_nodes, + std::vector const& nodes, + std::vector& new_elements, + unsigned min_elem_dim) { - const std::vector& nodes(_mesh.getNodes()); - const std::size_t nNodes(_mesh.getNumberOfNodes()); - std::vector id_map(nNodes); - const double half_eps(eps / 2.0); - const double sqr_eps(eps * eps); - std::iota(id_map.begin(), id_map.end(), 0); + auto addTetrahedron = + [&org_elem, &nodes, &new_elements](std::size_t id0, std::size_t id1, + std::size_t id2, std::size_t id3) + { + std::array tet_nodes{nodes[org_elem->getNode(id0)->getID()], + nodes[org_elem->getNode(id1)->getID()], + nodes[org_elem->getNode(id2)->getID()], + nodes[org_elem->getNode(id3)->getID()]}; + new_elements.push_back(new MeshLib::Tet(tet_nodes)); + }; - GeoLib::Grid const grid(nodes.begin(), nodes.end(), 64); + // TODO? + // In theory a node from the bottom triangle and a node from the top + // triangle that are not connected by an edge could collapse, resulting in a + // combination of tri and quad elements. This case is currently not tested. - for (std::size_t k = 0; k < nNodes; ++k) + // if one of the non-triangle edges collapsed, elem can be reduced to a + // pyramid, otherwise it will be two tets + if (n_unique_nodes == 5) { - MeshLib::Node const* const node(nodes[k]); - if (node->getID() != k) - { - continue; - } - std::vector const*> node_vectors( - grid.getPntVecsOfGridCellsIntersectingCube(*node, half_eps)); - - const std::size_t nVectors(node_vectors.size()); - for (std::size_t i = 0; i < nVectors; ++i) + for (unsigned i = 0; i < 5; ++i) { - const std::vector& cell_vector(*node_vectors[i]); - const std::size_t nGridCellNodes(cell_vector.size()); - for (std::size_t j = 0; j < nGridCellNodes; ++j) + for (unsigned j = i + 1; j < 6; ++j) { - MeshLib::Node const* const test_node(cell_vector[j]); - // are node indices already identical (i.e. nodes will be - // collapsed) - if (id_map[node->getID()] == id_map[test_node->getID()]) + if (i != j && org_elem->getNode(i)->getID() == + org_elem->getNode(j)->getID()) { - continue; - } + // non triangle edge collapsed + if (i % 3 == j % 3) + { + addTetrahedron((i + 1) % 3, (i + 2) % 3, i, + (i + 1) % 3 + 3); + addTetrahedron((i + 1) % 3 + 3, (i + 2) % 3, i, + (i + 2) % 3 + 3); + return 2; + } - // if test_node has already been collapsed to another node x, - // ignore it (if the current node would need to be collapsed - // with x it would already have happened when x was tested) - if (test_node->getID() != id_map[test_node->getID()]) - { - continue; - } + // triangle edge collapsed + const unsigned i_offset = (i > 2) ? i - 3 : i + 3; + const unsigned j_offset = (i > 2) ? j - 3 : j + 3; + const unsigned k = MeshToolsLib::lutPrismThirdNode(i, j); + if (k == std::numeric_limits::max()) + { + ERR("Unexpected error during prism reduction."); + return 0; + } + const unsigned k_offset = (i > 2) ? k - 3 : k + 3; - // calc distance - if (MathLib::sqrDist(*node, *test_node) < sqr_eps) - { - id_map[test_node->getID()] = node->getID(); + addTetrahedron(i_offset, j_offset, k_offset, i); + + const unsigned l = + (MathLib::isCoplanar(*org_elem->getNode(i_offset), + *org_elem->getNode(k_offset), + *org_elem->getNode(i), + *org_elem->getNode(k))) + ? j + : i; + const unsigned l_offset = (i > 2) ? l - 3 : l + 3; + addTetrahedron(l_offset, k_offset, i, k); + return 2; } } } } - return id_map; -} - -std::vector MeshRevision::constructNewNodesArray( - const std::vector& id_map) const -{ - const std::vector& nodes(_mesh.getNodes()); - const std::size_t nNodes(nodes.size()); - std::vector new_nodes; - new_nodes.reserve(nNodes); - for (std::size_t k = 0; k < nNodes; ++k) + else if (n_unique_nodes == 4) { - // all nodes that have not been collapsed with other nodes are copied - // into new array - if (nodes[k]->getID() == id_map[k]) - { - std::size_t const id(new_nodes.size()); - new_nodes.push_back(new MeshLib::Node( - (*nodes[k])[0], (*nodes[k])[1], (*nodes[k])[2], id)); - nodes[k]->setID(id); // the node in the old array gets the index of - // the same node in the new array - } - // the other nodes are not copied and get the index of the nodes they - // will have been collapsed with - else + MeshLib::Element* elem( + constructFourNodeElement(org_elem, nodes, min_elem_dim)); + if (elem) { - nodes[k]->setID(nodes[id_map[k]]->getID()); + new_elements.push_back(elem); } } - return new_nodes; + else if (n_unique_nodes == 3 && min_elem_dim < 3) + { + new_elements.push_back(constructTri(org_elem, nodes)); + } + else if (n_unique_nodes == 2 && min_elem_dim == 1) + { + new_elements.push_back(constructLine(org_elem, nodes)); + } + return 1; } -unsigned MeshRevision::getNumberOfUniqueNodes( - MeshLib::Element const* const element) +/// Lookup-table for returning four nodes connected to the two nodes (id1, id2) +/// forming an edge in a Hex. +std::array lutHexCuttingQuadNodes(unsigned id1, unsigned id2) { - unsigned const nNodes(element->getNumberOfBaseNodes()); - unsigned count(nNodes); - - for (unsigned i = 0; i < nNodes - 1; ++i) + std::array nodes{}; + if (id1 == 0 && id2 == 1) { - for (unsigned j = i + 1; j < nNodes; ++j) - { - if (element->getNode(i)->getID() == element->getNode(j)->getID()) - { - count--; - break; - } - } + nodes[0] = 3; + nodes[1] = 2; + nodes[2] = 5; + nodes[3] = 4; } - return count; -} - -template -void fillNodeProperty(std::vector& new_prop, - std::vector const& old_prop, - std::vector - node_ids) -{ - std::size_t const n_nodes = node_ids.size(); - for (std::size_t i = 0; i < n_nodes; ++i) + else if (id1 == 1 && id2 == 2) { - if (node_ids[i] != i) - { - continue; - } - new_prop.push_back(old_prop[i]); + nodes[0] = 0; + nodes[1] = 3; + nodes[2] = 6; + nodes[3] = 5; } -} - -template -void fillElemProperty(std::vector& new_prop, - std::vector const& old_prop, - std::vector - elem_ids) -{ - std::transform(elem_ids.cbegin(), elem_ids.cend(), - std::back_inserter(new_prop), - [&](std::size_t const i) { return old_prop[i]; }); -} - -MeshLib::Properties MeshRevision::copyProperties( - MeshLib::Properties const& props, - std::vector const& node_ids, - std::vector const& elem_ids) const -{ - auto const prop_names = props.getPropertyVectorNames(); - MeshLib::Properties new_properties; - - for (auto name : prop_names) + else if (id1 == 2 && id2 == 3) { - if (props.existsPropertyVector(name, MeshLib::MeshItemType::Node, - 1)) - { - auto p = props.getPropertyVector( - name, MeshLib::MeshItemType::Node, 1); - auto new_node_vec = new_properties.createNewPropertyVector( - name, MeshLib::MeshItemType::Node, 1); - fillNodeProperty(*new_node_vec, *p, node_ids); - continue; - } - if (props.existsPropertyVector(name, MeshLib::MeshItemType::Node, - 1)) - { - auto p = props.getPropertyVector( - name, MeshLib::MeshItemType::Node, 1); - auto new_node_vec = new_properties.createNewPropertyVector( - name, MeshLib::MeshItemType::Node, 1); - fillNodeProperty(*new_node_vec, *p, node_ids); - continue; - } - if (props.existsPropertyVector(name, - MeshLib::MeshItemType::Node, 1)) - { - auto p = props.getPropertyVector( - name, MeshLib::MeshItemType::Node, 1); - auto new_node_vec = new_properties.createNewPropertyVector( - name, MeshLib::MeshItemType::Node, 1); - fillNodeProperty(*new_node_vec, *p, node_ids); - continue; - } - if (props.existsPropertyVector(name, MeshLib::MeshItemType::Cell, - 1)) - { - auto p = props.getPropertyVector( - name, MeshLib::MeshItemType::Cell, 1); - auto new_cell_vec = new_properties.createNewPropertyVector( - name, MeshLib::MeshItemType::Cell, 1); - fillElemProperty(*new_cell_vec, *p, elem_ids); - continue; - } - if (props.existsPropertyVector(name, MeshLib::MeshItemType::Cell, - 1)) - { - auto p = props.getPropertyVector( - name, MeshLib::MeshItemType::Cell, 1); - auto new_cell_vec = new_properties.createNewPropertyVector( - name, MeshLib::MeshItemType::Cell, 1); - fillElemProperty(*new_cell_vec, *p, elem_ids); - continue; - } - if (props.existsPropertyVector(name, - MeshLib::MeshItemType::Cell, 1)) - { - auto p = props.getPropertyVector( - name, MeshLib::MeshItemType::Cell, 1); - auto new_cell_vec = new_properties.createNewPropertyVector( - name, MeshLib::MeshItemType::Cell, 1); - fillElemProperty(*new_cell_vec, *p, elem_ids); - continue; - } - WARN("PropertyVector {:s} not being converted.", name); + nodes[0] = 1; + nodes[1] = 0; + nodes[2] = 7; + nodes[3] = 6; } - return new_properties; -} - -std::size_t MeshRevision::subdivideElement( - MeshLib::Element const* const element, - std::vector const& nodes, - std::vector& elements) const -{ - if (element->getGeomType() == MeshLib::MeshElemType::QUAD) + else if (id1 == 3 && id2 == 0) { - return subdivideQuad(element, nodes, elements); + nodes[0] = 2; + nodes[1] = 1; + nodes[2] = 4; + nodes[3] = 7; } - if (element->getGeomType() == MeshLib::MeshElemType::HEXAHEDRON) + else if (id1 == 4 && id2 == 5) { - return subdivideHex(element, nodes, elements); + nodes[0] = 0; + nodes[1] = 1; + nodes[2] = 6; + nodes[3] = 7; } - if (element->getGeomType() == MeshLib::MeshElemType::PYRAMID) + else if (id1 == 5 && id2 == 6) { - return subdividePyramid(element, nodes, elements); + nodes[0] = 1; + nodes[1] = 2; + nodes[2] = 7; + nodes[3] = 4; } - if (element->getGeomType() == MeshLib::MeshElemType::PRISM) + else if (id1 == 6 && id2 == 7) { - return subdividePrism(element, nodes, elements); + nodes[0] = 2; + nodes[1] = 3; + nodes[2] = 4; + nodes[3] = 5; } - return 0; -} - -std::size_t MeshRevision::reduceElement( - MeshLib::Element const* const element, - unsigned n_unique_nodes, - std::vector const& nodes, - std::vector& elements, - unsigned min_elem_dim) const -{ - /*************** - * TODO: modify neighbouring elements if one elements has been subdivided - ***************/ - if (element->getGeomType() == MeshLib::MeshElemType::TRIANGLE && - min_elem_dim == 1) + else if (id1 == 7 && id2 == 4) { - elements.push_back(constructLine(element, nodes)); - return 1; + nodes[0] = 3; + nodes[1] = 0; + nodes[2] = 5; + nodes[3] = 6; } - if ((element->getGeomType() == MeshLib::MeshElemType::QUAD) || - (element->getGeomType() == MeshLib::MeshElemType::TETRAHEDRON)) + else if (id1 == 0 && id2 == 4) { - if (n_unique_nodes == 3 && min_elem_dim < 3) - { - elements.push_back(constructTri(element, nodes)); - } - else if (min_elem_dim == 1) - { - elements.push_back(constructLine(element, nodes)); - } - return 1; + nodes[0] = 3; + nodes[1] = 7; + nodes[2] = 5; + nodes[3] = 1; } - if (element->getGeomType() == MeshLib::MeshElemType::HEXAHEDRON) + else if (id1 == 1 && id2 == 5) { - return reduceHex(element, n_unique_nodes, nodes, elements, - min_elem_dim); + nodes[0] = 0; + nodes[1] = 4; + nodes[2] = 6; + nodes[3] = 2; } - if (element->getGeomType() == MeshLib::MeshElemType::PYRAMID) + else if (id1 == 2 && id2 == 6) { - this->reducePyramid(element, n_unique_nodes, nodes, elements, - min_elem_dim); - return 1; + nodes[0] = 1; + nodes[1] = 5; + nodes[2] = 7; + nodes[3] = 3; } - if (element->getGeomType() == MeshLib::MeshElemType::PRISM) + else if (id1 == 3 && id2 == 7) { - return reducePrism(element, n_unique_nodes, nodes, elements, - min_elem_dim); + nodes[0] = 2; + nodes[1] = 6; + nodes[2] = 4; + nodes[3] = 0; } - ERR("Unknown element type."); - return 0; -} - -unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, - unsigned n_unique_nodes, - std::vector const& nodes, - std::vector& new_elements, - unsigned min_elem_dim) const -{ - // TODO? - // if two diametral nodes collapse, all kinds of bizarre (2D-)element - // combinations could be the result. this case is currently not implemented! - - if (n_unique_nodes == 7) + else if (id1 == 1 && id2 == 0) { - // reduce to prism + pyramid - for (unsigned i = 0; i < 7; ++i) - { - for (unsigned j = i + 1; j < 8; ++j) - { - if (org_elem->getNode(i)->getID() == - org_elem->getNode(j)->getID()) - { - const std::array base_nodes( - this->lutHexCuttingQuadNodes(i, j)); - std::array pyr_nodes{ - nodes[org_elem->getNode(base_nodes[0])->getID()], - nodes[org_elem->getNode(base_nodes[1])->getID()], - nodes[org_elem->getNode(base_nodes[2])->getID()], - nodes[org_elem->getNode(base_nodes[3])->getID()], - nodes[org_elem->getNode(i)->getID()]}; - new_elements.push_back(new MeshLib::Pyramid(pyr_nodes)); - - if (i < 4 && j >= 4) - { - std::swap(i, j); - } - std::array prism_nodes{ - nodes[org_elem->getNode(base_nodes[0])->getID()], - nodes[org_elem->getNode(base_nodes[3])->getID()], - nodes[org_elem->getNode(this->lutHexDiametralNode(j)) - ->getID()], - nodes[org_elem->getNode(base_nodes[1])->getID()], - nodes[org_elem->getNode(base_nodes[2])->getID()], - nodes[org_elem->getNode(this->lutHexDiametralNode(i)) - ->getID()]}; - new_elements.push_back(new MeshLib::Prism(prism_nodes)); - return 2; - } - } - } + nodes[0] = 2; + nodes[1] = 3; + nodes[2] = 4; + nodes[3] = 5; } - else if (n_unique_nodes == 6) + else if (id1 == 2 && id2 == 1) + { + nodes[0] = 3; + nodes[1] = 0; + nodes[2] = 5; + nodes[3] = 6; + } + else if (id1 == 3 && id2 == 2) + { + nodes[0] = 0; + nodes[1] = 1; + nodes[2] = 6; + nodes[3] = 7; + } + else if (id1 == 0 && id2 == 3) + { + nodes[0] = 1; + nodes[1] = 2; + nodes[2] = 7; + nodes[3] = 4; + } + else if (id1 == 5 && id2 == 4) + { + nodes[0] = 1; + nodes[1] = 0; + nodes[2] = 7; + nodes[3] = 6; + } + else if (id1 == 6 && id2 == 5) + { + nodes[0] = 2; + nodes[1] = 1; + nodes[2] = 4; + nodes[3] = 7; + } + else if (id1 == 7 && id2 == 6) + { + nodes[0] = 3; + nodes[1] = 2; + nodes[2] = 5; + nodes[3] = 4; + } + else if (id1 == 4 && id2 == 7) + { + nodes[0] = 0; + nodes[1] = 3; + nodes[2] = 6; + nodes[3] = 5; + } + else if (id1 == 4 && id2 == 0) + { + nodes[0] = 7; + nodes[1] = 3; + nodes[2] = 1; + nodes[3] = 5; + } + else if (id1 == 5 && id2 == 1) + { + nodes[0] = 4; + nodes[1] = 0; + nodes[2] = 2; + nodes[3] = 6; + } + else if (id1 == 6 && id2 == 2) + { + nodes[0] = 5; + nodes[1] = 1; + nodes[2] = 3; + nodes[3] = 7; + } + else if (id1 == 7 && id2 == 3) + { + nodes[0] = 6; + nodes[1] = 2; + nodes[2] = 0; + nodes[3] = 4; + } + return nodes; +} + +/// Lookup-table for returning the diametral node id of the given node id in a +/// Hex. +unsigned lutHexDiametralNode(unsigned const id) +{ + constexpr std::array hex_diametral_nodes = { + {6, 7, 4, 5, 2, 3, 0, 1}}; + + return hex_diametral_nodes[id]; +} + +/// When a hex is subdivided into two prisms, this returns the nodes of the hex +/// edge that will serve as the back of one of the prisms. +std::pair lutHexBackNodes(unsigned const i, + unsigned const j, + unsigned const k, + unsigned const l) +{ + // collapsed edges are *not* connected + std::pair back(std::numeric_limits::max(), + std::numeric_limits::max()); + if (lutHexDiametralNode(i) == k) + { + back.first = i; + back.second = lutHexDiametralNode(l); + } + else if (lutHexDiametralNode(i) == l) + { + back.first = i; + back.second = lutHexDiametralNode(k); + } + else if (lutHexDiametralNode(j) == k) + { + back.first = j; + back.second = lutHexDiametralNode(l); + } + else if (lutHexDiametralNode(j) == l) + { + back.first = j; + back.second = lutHexDiametralNode(k); + } + // collapsed edges *are* connected + else if (i == k) + { + back.first = lutHexDiametralNode(l); + back.second = j; + } + else if (i == l) + { + back.first = lutHexDiametralNode(k); + back.second = j; + } + else if (j == k) + { + back.first = lutHexDiametralNode(l); + back.second = i; + } + else if (j == l) + { + back.first = lutHexDiametralNode(k); + back.second = i; + } + return back; +} + +// In an element with 5 unique nodes, return the node that will be the top of +// the resulting pyramid. +unsigned findPyramidTopNode(MeshLib::Element const& element, + std::array const& base_node_ids) +{ + const std::size_t nNodes(element.getNumberOfBaseNodes()); + for (std::size_t i = 0; i < nNodes; ++i) + { + bool top_node = true; + for (unsigned j = 0; j < 4; ++j) + { + if (element.getNode(i)->getID() == base_node_ids[j]) + { + top_node = false; + } + } + if (top_node) + { + return i; + } + } + return std::numeric_limits::max(); // should never be reached if + // called correctly +} + +/// Reduces a hexahedron element by removing collapsed nodes and constructing +/// one or more new elements from the remaining nodes. +/// @return The number of newly created elements +unsigned reduceHex(MeshLib::Element const* const org_elem, + unsigned n_unique_nodes, + std::vector const& nodes, + std::vector& new_elements, + unsigned min_elem_dim) +{ + // TODO? + // if two diametral nodes collapse, all kinds of bizarre (2D-)element + // combinations could be the result. this case is currently not implemented! + + if (n_unique_nodes == 7) + { + // reduce to prism + pyramid + for (unsigned i = 0; i < 7; ++i) + { + for (unsigned j = i + 1; j < 8; ++j) + { + if (org_elem->getNode(i)->getID() == + org_elem->getNode(j)->getID()) + { + const std::array base_nodes( + lutHexCuttingQuadNodes(i, j)); + std::array pyr_nodes{ + nodes[org_elem->getNode(base_nodes[0])->getID()], + nodes[org_elem->getNode(base_nodes[1])->getID()], + nodes[org_elem->getNode(base_nodes[2])->getID()], + nodes[org_elem->getNode(base_nodes[3])->getID()], + nodes[org_elem->getNode(i)->getID()]}; + new_elements.push_back(new MeshLib::Pyramid(pyr_nodes)); + + if (i < 4 && j >= 4) + { + std::swap(i, j); + } + std::array prism_nodes{ + nodes[org_elem->getNode(base_nodes[0])->getID()], + nodes[org_elem->getNode(base_nodes[3])->getID()], + nodes[org_elem->getNode(lutHexDiametralNode(j)) + ->getID()], + nodes[org_elem->getNode(base_nodes[1])->getID()], + nodes[org_elem->getNode(base_nodes[2])->getID()], + nodes[org_elem->getNode(lutHexDiametralNode(i)) + ->getID()]}; + new_elements.push_back(new MeshLib::Prism(prism_nodes)); + return 2; + } + } + } + } + else if (n_unique_nodes == 6) { // reduce to prism for (unsigned i = 0; i < 6; ++i) @@ -674,28 +723,24 @@ unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, { std::array prism_nodes{ nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(0)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(0)))) ->getID()], nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(1)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(1)))) ->getID()], nodes[org_elem ->getNode(getNodeIDinElement(*org_elem, face->getNode(2))) ->getID()], nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(2)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(2)))) ->getID()], nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(3)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(3)))) ->getID()], nodes[org_elem ->getNode(getNodeIDinElement(*org_elem, @@ -710,28 +755,24 @@ unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, { std::array prism_nodes{ nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(0)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(0)))) ->getID()], nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(3)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(3)))) ->getID()], nodes[org_elem ->getNode(getNodeIDinElement(*org_elem, face->getNode(2))) ->getID()], nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(1)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(1)))) ->getID()], nodes[org_elem - ->getNode( - this->lutHexDiametralNode(getNodeIDinElement( - *org_elem, face->getNode(2)))) + ->getNode(lutHexDiametralNode(getNodeIDinElement( + *org_elem, face->getNode(2)))) ->getID()], nodes[org_elem ->getNode(getNodeIDinElement(*org_elem, @@ -762,7 +803,7 @@ unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, org_elem->getNode(l)->getID()) { const std::pair back( - this->lutHexBackNodes(i, j, k, l)); + lutHexBackNodes(i, j, k, l)); if (back.first == std::numeric_limits::max() || back.second == @@ -774,8 +815,8 @@ unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, } std::array cutting_plane( - this->lutHexCuttingQuadNodes(back.first, - back.second)); + lutHexCuttingQuadNodes(back.first, + back.second)); std::array pris1_nodes{ const_cast( org_elem->getNode(back.first)), @@ -790,32 +831,30 @@ unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, const_cast( org_elem->getNode(cutting_plane[2]))}; auto* prism1(new MeshLib::Prism(pris1_nodes)); - unsigned nNewElements = this->reducePrism( - prism1, 5, nodes, new_elements, - min_elem_dim); + unsigned nNewElements = + reducePrism(prism1, 5, nodes, new_elements, + min_elem_dim); delete prism1; std::array pris2_nodes{ const_cast( org_elem->getNode( - this->lutHexDiametralNode( - back.first))), + lutHexDiametralNode(back.first))), const_cast( org_elem->getNode(cutting_plane[0])), const_cast( org_elem->getNode(cutting_plane[3])), const_cast( org_elem->getNode( - this->lutHexDiametralNode( - back.second))), + lutHexDiametralNode(back.second))), const_cast( org_elem->getNode(cutting_plane[1])), const_cast( org_elem->getNode(cutting_plane[2]))}; auto* prism2(new MeshLib::Prism(pris2_nodes)); - nNewElements += this->reducePrism( - prism2, 5, nodes, new_elements, - min_elem_dim); + nNewElements += + reducePrism(prism2, 5, nodes, new_elements, + min_elem_dim); delete prism2; return nNewElements; } @@ -831,8 +870,7 @@ unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, std::array first_four_nodes = { {tet1->getNode(0)->getID(), tet1->getNode(1)->getID(), tet1->getNode(2)->getID(), tet1->getNode(3)->getID()}}; - unsigned fifth_node( - this->findPyramidTopNode(*org_elem, first_four_nodes)); + unsigned fifth_node(findPyramidTopNode(*org_elem, first_four_nodes)); bool tet_changed(false); if (tet1->getGeomType() == MeshLib::MeshElemType::QUAD) @@ -882,404 +920,384 @@ unsigned MeshRevision::reduceHex(MeshLib::Element const* const org_elem, return 0; } -void MeshRevision::reducePyramid(MeshLib::Element const* const org_elem, - unsigned n_unique_nodes, - std::vector const& nodes, - std::vector& new_elements, - unsigned min_elem_dim) const +/// Subdivides an element if it has a face that is not coplanar +/// @param element the element that will be subdivided +/// @param nodes vector containing the nodes the elements originated by the +/// subdivision are based on +/// @param elements vector of MeshLib::Elements; the elements originated by the +/// subdivision will be inserted into elements +/// @return the number of elements originated by the subdivision +std::size_t subdivideElement(MeshLib::Element const* const element, + std::vector const& nodes, + std::vector& elements) { - if (n_unique_nodes == 4) + if (element->getGeomType() == MeshLib::MeshElemType::QUAD) { - MeshLib::Element* elem( - constructFourNodeElement(org_elem, nodes, min_elem_dim)); - if (elem) - { - new_elements.push_back(elem); - } + return subdivideQuad(element, nodes, elements); } - else if (n_unique_nodes == 3 && min_elem_dim < 3) + if (element->getGeomType() == MeshLib::MeshElemType::HEXAHEDRON) { - new_elements.push_back(constructTri(org_elem, nodes)); + return subdivideHex(element, nodes, elements); } - else if (n_unique_nodes == 2 && min_elem_dim == 1) + if (element->getGeomType() == MeshLib::MeshElemType::PYRAMID) { - new_elements.push_back(constructLine(org_elem, nodes)); + return subdividePyramid(element, nodes, elements); + } + if (element->getGeomType() == MeshLib::MeshElemType::PRISM) + { + return subdividePrism(element, nodes, elements); } + return 0; } -unsigned MeshRevision::reducePrism(MeshLib::Element const* const org_elem, - unsigned n_unique_nodes, - std::vector const& nodes, - std::vector& new_elements, - unsigned min_elem_dim) const +// Revises an element by removing collapsed nodes, using the nodes vector from +// the result mesh. +std::size_t reduceElement(MeshLib::Element const* const element, + unsigned n_unique_nodes, + std::vector const& nodes, + std::vector& elements, + unsigned min_elem_dim) { - auto addTetrahedron = - [&org_elem, &nodes, &new_elements](std::size_t id0, std::size_t id1, - std::size_t id2, std::size_t id3) + /*************** + * TODO: modify neighbouring elements if one elements has been subdivided + ***************/ + if (element->getGeomType() == MeshLib::MeshElemType::TRIANGLE && + min_elem_dim == 1) { - std::array tet_nodes{nodes[org_elem->getNode(id0)->getID()], - nodes[org_elem->getNode(id1)->getID()], - nodes[org_elem->getNode(id2)->getID()], - nodes[org_elem->getNode(id3)->getID()]}; - new_elements.push_back(new MeshLib::Tet(tet_nodes)); - }; - - // TODO? - // In theory a node from the bottom triangle and a node from the top - // triangle that are not connected by an edge could collapse, resulting in a - // combination of tri and quad elements. This case is currently not tested. - - // if one of the non-triangle edges collapsed, elem can be reduced to a - // pyramid, otherwise it will be two tets - if (n_unique_nodes == 5) + elements.push_back(constructLine(element, nodes)); + return 1; + } + if ((element->getGeomType() == MeshLib::MeshElemType::QUAD) || + (element->getGeomType() == MeshLib::MeshElemType::TETRAHEDRON)) { - for (unsigned i = 0; i < 5; ++i) + if (n_unique_nodes == 3 && min_elem_dim < 3) { - for (unsigned j = i + 1; j < 6; ++j) - { - if (i != j && org_elem->getNode(i)->getID() == - org_elem->getNode(j)->getID()) - { - // non triangle edge collapsed - if (i % 3 == j % 3) - { - addTetrahedron((i + 1) % 3, (i + 2) % 3, i, - (i + 1) % 3 + 3); - addTetrahedron((i + 1) % 3 + 3, (i + 2) % 3, i, - (i + 2) % 3 + 3); - return 2; - } - - // triangle edge collapsed - const unsigned i_offset = (i > 2) ? i - 3 : i + 3; - const unsigned j_offset = (i > 2) ? j - 3 : j + 3; - const unsigned k = lutPrismThirdNode(i, j); - if (k == std::numeric_limits::max()) - { - ERR("Unexpected error during prism reduction."); - return 0; - } - const unsigned k_offset = (i > 2) ? k - 3 : k + 3; - - addTetrahedron(i_offset, j_offset, k_offset, i); - - const unsigned l = - (MathLib::isCoplanar(*org_elem->getNode(i_offset), - *org_elem->getNode(k_offset), - *org_elem->getNode(i), - *org_elem->getNode(k))) - ? j - : i; - const unsigned l_offset = (i > 2) ? l - 3 : l + 3; - addTetrahedron(l_offset, k_offset, i, k); - return 2; - } - } + elements.push_back(constructTri(element, nodes)); } - } - else if (n_unique_nodes == 4) - { - MeshLib::Element* elem( - constructFourNodeElement(org_elem, nodes, min_elem_dim)); - if (elem) + else if (min_elem_dim == 1) { - new_elements.push_back(elem); + elements.push_back(constructLine(element, nodes)); } + return 1; } - else if (n_unique_nodes == 3 && min_elem_dim < 3) + if (element->getGeomType() == MeshLib::MeshElemType::HEXAHEDRON) { - new_elements.push_back(constructTri(org_elem, nodes)); + return reduceHex(element, n_unique_nodes, nodes, elements, + min_elem_dim); } - else if (n_unique_nodes == 2 && min_elem_dim == 1) + if (element->getGeomType() == MeshLib::MeshElemType::PYRAMID) { - new_elements.push_back(constructLine(org_elem, nodes)); + reducePyramid(element, n_unique_nodes, nodes, elements, min_elem_dim); + return 1; + } + if (element->getGeomType() == MeshLib::MeshElemType::PRISM) + { + return reducePrism(element, n_unique_nodes, nodes, elements, + min_elem_dim); } - return 1; -} -unsigned MeshRevision::findPyramidTopNode( - MeshLib::Element const& element, - std::array const& base_node_ids) + ERR("Unknown element type."); + return 0; +} +/// Calculates the number of unique nodes in an element (i.e. uncollapsed +/// nodes). +unsigned getNumberOfUniqueNodes(MeshLib::Element const* const element) { - const std::size_t nNodes(element.getNumberOfBaseNodes()); - for (std::size_t i = 0; i < nNodes; ++i) + unsigned const nNodes(element->getNumberOfBaseNodes()); + unsigned count(nNodes); + + for (unsigned i = 0; i < nNodes - 1; ++i) { - bool top_node = true; - for (unsigned j = 0; j < 4; ++j) + for (unsigned j = i + 1; j < nNodes; ++j) { - if (element.getNode(i)->getID() == base_node_ids[j]) + if (element->getNode(i)->getID() == element->getNode(j)->getID()) { - top_node = false; + count--; + break; } } - if (top_node) + } + return count; +} + +template +void fillNodeProperty(std::vector& new_prop, + std::vector const& old_prop, + std::vector + node_ids) +{ + std::size_t const n_nodes = node_ids.size(); + for (std::size_t i = 0; i < n_nodes; ++i) + { + if (node_ids[i] != i) { - return i; + continue; } + new_prop.push_back(old_prop[i]); } - return std::numeric_limits::max(); // should never be reached if - // called correctly } -unsigned MeshRevision::lutHexDiametralNode(unsigned id) +template +void fillElemProperty(std::vector& new_prop, + std::vector const& old_prop, + std::vector + elem_ids) { - return _hex_diametral_nodes[id]; + std::transform(elem_ids.cbegin(), elem_ids.cend(), + std::back_inserter(new_prop), + [&](std::size_t const i) { return old_prop[i]; }); } -std::array MeshRevision::lutHexCuttingQuadNodes(unsigned id1, - unsigned id2) +/// Copies all scalar arrays according to the restructured Node- and +/// Element-vectors after the mesh revision process (i.e. collapsed nodes, split +/// elements, etc.) +MeshLib::Properties copyProperties(MeshLib::Properties const& props, + std::vector const& node_ids, + std::vector const& elem_ids) { - std::array nodes{}; - if (id1 == 0 && id2 == 1) - { - nodes[0] = 3; - nodes[1] = 2; - nodes[2] = 5; - nodes[3] = 4; - } - else if (id1 == 1 && id2 == 2) - { - nodes[0] = 0; - nodes[1] = 3; - nodes[2] = 6; - nodes[3] = 5; - } - else if (id1 == 2 && id2 == 3) - { - nodes[0] = 1; - nodes[1] = 0; - nodes[2] = 7; - nodes[3] = 6; - } - else if (id1 == 3 && id2 == 0) - { - nodes[0] = 2; - nodes[1] = 1; - nodes[2] = 4; - nodes[3] = 7; - } - else if (id1 == 4 && id2 == 5) - { - nodes[0] = 0; - nodes[1] = 1; - nodes[2] = 6; - nodes[3] = 7; - } - else if (id1 == 5 && id2 == 6) + auto const prop_names = props.getPropertyVectorNames(); + MeshLib::Properties new_properties; + + for (auto name : prop_names) { - nodes[0] = 1; - nodes[1] = 2; - nodes[2] = 7; - nodes[3] = 4; + if (props.existsPropertyVector(name, MeshLib::MeshItemType::Node, + 1)) + { + auto p = props.getPropertyVector( + name, MeshLib::MeshItemType::Node, 1); + auto new_node_vec = new_properties.createNewPropertyVector( + name, MeshLib::MeshItemType::Node, 1); + fillNodeProperty(*new_node_vec, *p, node_ids); + continue; + } + if (props.existsPropertyVector(name, MeshLib::MeshItemType::Node, + 1)) + { + auto p = props.getPropertyVector( + name, MeshLib::MeshItemType::Node, 1); + auto new_node_vec = new_properties.createNewPropertyVector( + name, MeshLib::MeshItemType::Node, 1); + fillNodeProperty(*new_node_vec, *p, node_ids); + continue; + } + if (props.existsPropertyVector(name, + MeshLib::MeshItemType::Node, 1)) + { + auto p = props.getPropertyVector( + name, MeshLib::MeshItemType::Node, 1); + auto new_node_vec = new_properties.createNewPropertyVector( + name, MeshLib::MeshItemType::Node, 1); + fillNodeProperty(*new_node_vec, *p, node_ids); + continue; + } + if (props.existsPropertyVector(name, MeshLib::MeshItemType::Cell, + 1)) + { + auto p = props.getPropertyVector( + name, MeshLib::MeshItemType::Cell, 1); + auto new_cell_vec = new_properties.createNewPropertyVector( + name, MeshLib::MeshItemType::Cell, 1); + fillElemProperty(*new_cell_vec, *p, elem_ids); + continue; + } + if (props.existsPropertyVector(name, MeshLib::MeshItemType::Cell, + 1)) + { + auto p = props.getPropertyVector( + name, MeshLib::MeshItemType::Cell, 1); + auto new_cell_vec = new_properties.createNewPropertyVector( + name, MeshLib::MeshItemType::Cell, 1); + fillElemProperty(*new_cell_vec, *p, elem_ids); + continue; + } + if (props.existsPropertyVector(name, + MeshLib::MeshItemType::Cell, 1)) + { + auto p = props.getPropertyVector( + name, MeshLib::MeshItemType::Cell, 1); + auto new_cell_vec = new_properties.createNewPropertyVector( + name, MeshLib::MeshItemType::Cell, 1); + fillElemProperty(*new_cell_vec, *p, elem_ids); + continue; + } + WARN("PropertyVector {:s} not being converted.", name); } - else if (id1 == 6 && id2 == 7) + return new_properties; +} +} // namespace + +namespace MeshToolsLib +{ +MeshRevision::MeshRevision(MeshLib::Mesh& mesh) : _mesh(mesh) {} + +unsigned MeshRevision::getNumberOfCollapsibleNodes(double eps) const +{ + std::vector id_map(this->collapseNodeIndices(eps)); + std::size_t nNodes(id_map.size()); + unsigned count(0); + for (std::size_t i = 0; i < nNodes; ++i) { - nodes[0] = 2; - nodes[1] = 3; - nodes[2] = 4; - nodes[3] = 5; + if (i != id_map[i]) + { + count++; + } } - else if (id1 == 7 && id2 == 4) + return count; +} + +MeshLib::Mesh* MeshRevision::simplifyMesh(const std::string& new_mesh_name, + double eps, + unsigned min_elem_dim) const +{ + if (this->_mesh.getNumberOfElements() == 0) { - nodes[0] = 3; - nodes[1] = 0; - nodes[2] = 5; - nodes[3] = 6; + return nullptr; } - else if (id1 == 0 && id2 == 4) + + std::vector const& elements(this->_mesh.getElements()); + auto const node_ids = collapseNodeIndices(eps); + std::vector new_nodes = + this->constructNewNodesArray(node_ids); + std::vector new_elements; + std::vector element_ids; + + for (std::size_t k(0); k < elements.size(); ++k) { - nodes[0] = 3; - nodes[1] = 7; - nodes[2] = 5; - nodes[3] = 1; - } - else if (id1 == 1 && id2 == 5) - { - nodes[0] = 0; - nodes[1] = 4; - nodes[2] = 6; - nodes[3] = 2; - } - else if (id1 == 2 && id2 == 6) - { - nodes[0] = 1; - nodes[1] = 5; - nodes[2] = 7; - nodes[3] = 3; - } - else if (id1 == 3 && id2 == 7) - { - nodes[0] = 2; - nodes[1] = 6; - nodes[2] = 4; - nodes[3] = 0; + MeshLib::Element const* const elem(elements[k]); + unsigned const n_unique_nodes(getNumberOfUniqueNodes(elem)); + if (n_unique_nodes == elem->getNumberOfBaseNodes() && + elem->getDimension() >= min_elem_dim) + { + ElementErrorCode e(elem->validate()); + if (e[ElementErrorFlag::NonCoplanar]) + { + std::size_t const n_new_elements( + subdivideElement(elem, new_nodes, new_elements)); + if (n_new_elements == 0) + { + ERR("Element {:d} has unknown element type.", k); + _mesh.resetNodeIDs(); + BaseLib::cleanupVectorElements(new_nodes, new_elements); + return nullptr; + } + element_ids.insert(element_ids.end(), n_new_elements, k); + } + else + { + new_elements.push_back(MeshLib::copyElement(elem, new_nodes)); + element_ids.push_back(k); + } + } + else if (n_unique_nodes < elem->getNumberOfBaseNodes() && + n_unique_nodes > 1) + { + std::size_t const n_new_elements(reduceElement( + elem, n_unique_nodes, new_nodes, new_elements, min_elem_dim)); + element_ids.insert(element_ids.end(), n_new_elements, k); + } + else + { + ERR("Something is wrong, more unique nodes than actual nodes"); + } } - else if (id1 == 1 && id2 == 0) - { - nodes[0] = 2; - nodes[1] = 3; - nodes[2] = 4; - nodes[3] = 5; - } - else if (id1 == 2 && id2 == 1) - { - nodes[0] = 3; - nodes[1] = 0; - nodes[2] = 5; - nodes[3] = 6; - } - else if (id1 == 3 && id2 == 2) - { - nodes[0] = 0; - nodes[1] = 1; - nodes[2] = 6; - nodes[3] = 7; - } - else if (id1 == 0 && id2 == 3) - { - nodes[0] = 1; - nodes[1] = 2; - nodes[2] = 7; - nodes[3] = 4; - } - else if (id1 == 5 && id2 == 4) - { - nodes[0] = 1; - nodes[1] = 0; - nodes[2] = 7; - nodes[3] = 6; - } - else if (id1 == 6 && id2 == 5) - { - nodes[0] = 2; - nodes[1] = 1; - nodes[2] = 4; - nodes[3] = 7; - } - else if (id1 == 7 && id2 == 6) - { - nodes[0] = 3; - nodes[1] = 2; - nodes[2] = 5; - nodes[3] = 4; - } - else if (id1 == 4 && id2 == 7) - { - nodes[0] = 0; - nodes[1] = 3; - nodes[2] = 6; - nodes[3] = 5; - } - else if (id1 == 4 && id2 == 0) - { - nodes[0] = 7; - nodes[1] = 3; - nodes[2] = 1; - nodes[3] = 5; - } - else if (id1 == 5 && id2 == 1) - { - nodes[0] = 4; - nodes[1] = 0; - nodes[2] = 2; - nodes[3] = 6; - } - else if (id1 == 6 && id2 == 2) - { - nodes[0] = 5; - nodes[1] = 1; - nodes[2] = 3; - nodes[3] = 7; - } - else if (id1 == 7 && id2 == 3) + auto const& props = _mesh.getProperties(); + MeshLib::Properties const new_properties = + copyProperties(props, node_ids, element_ids); + + _mesh.resetNodeIDs(); + if (!new_elements.empty()) { - nodes[0] = 6; - nodes[1] = 2; - nodes[2] = 0; - nodes[3] = 4; + return new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements, + new_properties); } - return nodes; + + BaseLib::cleanupVectorElements(new_nodes, new_elements); + return nullptr; } -std::pair MeshRevision::lutHexBackNodes(unsigned i, - unsigned j, - unsigned k, - unsigned l) +std::vector MeshRevision::collapseNodeIndices(double eps) const { - // collapsed edges are *not* connected - std::pair back(std::numeric_limits::max(), - std::numeric_limits::max()); - if (lutHexDiametralNode(i) == k) - { - back.first = i; - back.second = lutHexDiametralNode(l); - } - else if (lutHexDiametralNode(i) == l) - { - back.first = i; - back.second = lutHexDiametralNode(k); - } - else if (lutHexDiametralNode(j) == k) - { - back.first = j; - back.second = lutHexDiametralNode(l); - } - else if (lutHexDiametralNode(j) == l) - { - back.first = j; - back.second = lutHexDiametralNode(k); - } - // collapsed edges *are* connected - else if (i == k) - { - back.first = lutHexDiametralNode(l); - back.second = j; - } - else if (i == l) - { - back.first = lutHexDiametralNode(k); - back.second = j; - } - else if (j == k) - { - back.first = lutHexDiametralNode(l); - back.second = i; - } - else if (j == l) + const std::vector& nodes(_mesh.getNodes()); + const std::size_t nNodes(_mesh.getNumberOfNodes()); + std::vector id_map(nNodes); + const double half_eps(eps / 2.0); + const double sqr_eps(eps * eps); + std::iota(id_map.begin(), id_map.end(), 0); + + GeoLib::Grid const grid(nodes.begin(), nodes.end(), 64); + + for (std::size_t k = 0; k < nNodes; ++k) { - back.first = lutHexDiametralNode(k); - back.second = i; + MeshLib::Node const* const node(nodes[k]); + if (node->getID() != k) + { + continue; + } + std::vector const*> node_vectors( + grid.getPntVecsOfGridCellsIntersectingCube(*node, half_eps)); + + const std::size_t nVectors(node_vectors.size()); + for (std::size_t i = 0; i < nVectors; ++i) + { + const std::vector& cell_vector(*node_vectors[i]); + const std::size_t nGridCellNodes(cell_vector.size()); + for (std::size_t j = 0; j < nGridCellNodes; ++j) + { + MeshLib::Node const* const test_node(cell_vector[j]); + // are node indices already identical (i.e. nodes will be + // collapsed) + if (id_map[node->getID()] == id_map[test_node->getID()]) + { + continue; + } + + // if test_node has already been collapsed to another node x, + // ignore it (if the current node would need to be collapsed + // with x it would already have happened when x was tested) + if (test_node->getID() != id_map[test_node->getID()]) + { + continue; + } + + // calc distance + if (MathLib::sqrDist(*node, *test_node) < sqr_eps) + { + id_map[test_node->getID()] = node->getID(); + } + } + } } - return back; + return id_map; } -unsigned MeshRevision::lutPrismThirdNode(unsigned id1, unsigned id2) +std::vector MeshRevision::constructNewNodesArray( + const std::vector& id_map) const { - if ((id1 == 0 && id2 == 1) || (id1 == 1 && id2 == 2)) - { - return 2; - } - if ((id1 == 1 && id2 == 2) || (id1 == 2 && id2 == 1)) - { - return 0; - } - if ((id1 == 0 && id2 == 2) || (id1 == 2 && id2 == 0)) - { - return 1; - } - if ((id1 == 3 && id2 == 4) || (id1 == 4 && id2 == 3)) - { - return 5; - } - if ((id1 == 4 && id2 == 5) || (id1 == 5 && id2 == 4)) - { - return 3; - } - if ((id1 == 3 && id2 == 5) || (id1 == 5 && id2 == 3)) + const std::vector& nodes(_mesh.getNodes()); + const std::size_t nNodes(nodes.size()); + std::vector new_nodes; + new_nodes.reserve(nNodes); + for (std::size_t k = 0; k < nNodes; ++k) { - return 4; + // all nodes that have not been collapsed with other nodes are copied + // into new array + if (nodes[k]->getID() == id_map[k]) + { + std::size_t const id(new_nodes.size()); + new_nodes.push_back(new MeshLib::Node( + (*nodes[k])[0], (*nodes[k])[1], (*nodes[k])[2], id)); + nodes[k]->setID(id); // the node in the old array gets the index of + // the same node in the new array + } + // the other nodes are not copied and get the index of the nodes they + // will have been collapsed with + else + { + nodes[k]->setID(nodes[id_map[k]]->getID()); + } } - return std::numeric_limits::max(); + return new_nodes; } + } // namespace MeshToolsLib diff --git a/MeshToolsLib/MeshEditing/MeshRevision.h b/MeshToolsLib/MeshEditing/MeshRevision.h index 66efb863d29..2c8eacdab54 100644 --- a/MeshToolsLib/MeshEditing/MeshRevision.h +++ b/MeshToolsLib/MeshEditing/MeshRevision.h @@ -19,14 +19,11 @@ #include #include -#include "MeshLib/Properties.h" - // forward declaration namespace MeshLib { class Mesh; class Node; -class Element; } // namespace MeshLib namespace MeshToolsLib @@ -46,8 +43,6 @@ class MeshRevision */ explicit MeshRevision(MeshLib::Mesh& mesh); - virtual ~MeshRevision() = default; - /// Returns the number of potentially collapsible nodes unsigned getNumberOfCollapsibleNodes( double eps = std::numeric_limits::epsilon()) const; @@ -76,102 +71,8 @@ class MeshRevision std::vector constructNewNodesArray( const std::vector& id_map) const; - /// Calculates the number of unique nodes in an element (i.e. uncollapsed - /// nodes) - static unsigned getNumberOfUniqueNodes( - MeshLib::Element const* const element); - - /** - * Copies all scalar arrays according to the restructured Node- and - * Element-vectors after the mesh revision process (i.e. collapsed nodes, - * split elements, etc.) - */ - MeshLib::Properties copyProperties( - MeshLib::Properties const& props, - std::vector const& node_ids, - std::vector const& elem_ids) const; - - /// Subdivides an element if it has a face that is not coplanar - /// @param element the element that will be subdivided - /// @param nodes vector containing the nodes the elements originated by the - /// subdivision are based on - /// @param elements vector of MeshLib::Elements; the elements originated by - /// the subdivision will be inserted into elements - /// @return the number of elements originated by the subdivision - std::size_t subdivideElement( - MeshLib::Element const* const element, - std::vector const& nodes, - std::vector& elements) const; - - // Revises an element by removing collapsed nodes, using the nodes vector - // from the result mesh. - std::size_t reduceElement(MeshLib::Element const* const element, - unsigned n_unique_nodes, - const std::vector& nodes, - std::vector& elements, - unsigned min_elem_dim) const; - - /// Cleans up all nodes and elements if something went wrong - void cleanUp(std::vector& nodes, - std::vector& new_elements) const; - - /** - * Reduces a hexahedron element by removing collapsed nodes and constructing - * one or more new elements from the remaining nodes. - * @return The number of newly created elements - */ - unsigned reduceHex(MeshLib::Element const* const org_elem, - unsigned n_unique_nodes, - const std::vector& nodes, - std::vector& new_elements, - unsigned min_elem_dim) const; - /// Reduces a pyramid element by removing collapsed nodes and constructing a - /// new elements from the remaining nodes. - void reducePyramid(MeshLib::Element const* const org_elem, - unsigned n_unique_nodes, - const std::vector& nodes, - std::vector& new_elements, - unsigned min_elem_dim) const; - /** - * Reduces a prism element by removing collapsed nodes and constructing one - * or two new elements from the remaining nodes. - * @return The number of newly created elements - */ - unsigned reducePrism(MeshLib::Element const* const org_elem, - unsigned n_unique_nodes, - std::vector const& nodes, - std::vector& new_elements, - unsigned min_elem_dim) const; - - // In an element with 5 unique nodes, return the node that will be the top - // of the resulting pyramid - static unsigned findPyramidTopNode( - MeshLib::Element const& element, - std::array const& base_node_ids); - - /// Lookup-table for returning the diametral node id of the given node id in - /// a Hex - static unsigned lutHexDiametralNode(unsigned id); - - /// Lookup-table for returning four nodes connected to the two nodes (id1, - /// id2) forming an edge in a Hex - static std::array lutHexCuttingQuadNodes(unsigned id1, - unsigned id2); - - /// When a hex is subdivided into two prisms, this returns the nodes of the - /// hex edge that will serve as the back of one of the prisms. - static std::pair lutHexBackNodes(unsigned i, unsigned j, - unsigned k, - unsigned l); - - /// Lookup-table for returning the third node of bottom or top triangle - /// given the other two - static unsigned lutPrismThirdNode(unsigned id1, unsigned id2); - /// The original mesh used for constructing the class MeshLib::Mesh& _mesh; - - static const std::array _hex_diametral_nodes; }; } // namespace MeshToolsLib diff --git a/Tests/MeshLib/TestMeshRevision.cpp b/Tests/MeshLib/TestMeshRevision.cpp index 625be93f5fd..3e7f90bdd5b 100644 --- a/Tests/MeshLib/TestMeshRevision.cpp +++ b/Tests/MeshLib/TestMeshRevision.cpp @@ -537,3 +537,41 @@ TEST(MeshEditing, Prism2Tri) delete result; } + +namespace MeshToolsLib +{ +extern unsigned lutPrismThirdNode(unsigned, unsigned); +} + +TEST(MeshEditing, PrismLutThirdNode) +{ + constexpr auto max = std::numeric_limits::max(); + // clang-format off + // Expected results in a matrix 6x6 to cover all possible combinations + constexpr std::array, 6> expected = + // 0 1 2 3 4 5 + {{{{max, 2 , 1 , max, max, max}}, // 0 + {{2 , max, 0 , max, max, max}}, // 1 + {{1 , 0 , max, max, max, max}}, // 2 + {{max, max, max, max, 5 , 4 }}, // 3 + {{max, max, max, 5 , max, 3 }}, // 4 + {{max, max, max, 4 , 3 , max}}}}; // 5 + // clang-format on + + for (unsigned i = 0; i < 6; ++i) + { + for (unsigned j = 0; j < 6; ++j) + { + EXPECT_EQ(MeshToolsLib::lutPrismThirdNode(i, j), expected[i][j]) + << "for (i,j) = " << i << ", " << j; + } + } + + // Test for outside values + EXPECT_EQ(MeshToolsLib::lutPrismThirdNode(-1, -1), max); + EXPECT_EQ(MeshToolsLib::lutPrismThirdNode(0, -1), max); + EXPECT_EQ(MeshToolsLib::lutPrismThirdNode(-1, 0), max); + EXPECT_EQ(MeshToolsLib::lutPrismThirdNode(5, 6), max); + EXPECT_EQ(MeshToolsLib::lutPrismThirdNode(6, 5), max); + EXPECT_EQ(MeshToolsLib::lutPrismThirdNode(6, 6), max); +}