GitHub PullRequestReviewCommentEvent at 2014-03-12 08:09:31
Repository: ogs
Actor: rinkk
Payload: {"comment":{"url":"https://api.github.com/repos/ufz/ogs/pulls/comments/10508440","id":10508440,"diff_hunk":"@@ -0,0 +1,812 @@\n+/**\n+ * \\file   MeshRevision.cpp\n+ * \\author Karsten Rink\n+ * \\date   2014-02-14\n+ * \\brief  Implementation of the MeshRevision class. \n+ *\n+ * \\copyright\n+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)\n+ *            Distributed under a Modified BSD License.\n+ *              See accompanying file LICENSE.txt or\n+ *              http://www.opengeosys.org/project/license\n+ *\n+ */\n+\n+// ThirdParty/logog\n+#include \"logog/include/logog.hpp\"\n+\n+#include \"MeshRevision.h\"\n+\n+// GeoLib\n+#include \"Grid.h\"\n+\n+// MeshLib\n+#include \"Mesh.h\"\n+#include \"Node.h\"\n+#include \"Elements/Element.h\"\n+#include \"Elements/Line.h\"\n+#include \"Elements/Tri.h\"\n+#include \"Elements/Quad.h\"\n+#include \"Elements/Tet.h\"\n+#include \"Elements/Hex.h\"\n+#include \"Elements/Pyramid.h\"\n+#include \"Elements/Prism.h\"\n+\n+\n+namespace MeshLib {\n+\n+MeshRevision::MeshRevision(const MeshLib::Mesh &mesh) :\n+\t_mesh (mesh)\n+{}\n+\n+\n+MeshLib::Mesh* MeshRevision::collapseNodes(const std::string &new_mesh_name, double eps)\n+{\n+\tstd::vector<MeshLib::Node*> new_nodes = this->constructNewNodesArray(this->collapseNodeIndeces(eps));\n+\tstd::vector<MeshLib::Element*> new_elements;\n+\tnew_elements.reserve(this->_mesh.getNElements());\n+\tconst std::vector<MeshLib::Element*> elements (this->_mesh.getElements());\n+\tfor (auto elem = elements.begin(); elem != elements.end(); ++elem)\n+\t\tnew_elements.push_back(copyElement(*elem, new_nodes));\n+\treturn new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements);\n+}\n+\n+unsigned MeshRevision::getNCollapsableNodes(double eps) const\n+{\n+\tstd::vector<std::size_t> id_map (this->collapseNodeIndeces(eps));\n+\tstd::size_t nNodes (id_map.size());\n+\tunsigned count (0);\n+\tfor (std::size_t i=0; i<nNodes; ++i)\n+\t\tif (i != id_map[i])\n+\t\t\tcount++;\n+\treturn count;\n+}\n+\n+MeshLib::Mesh* MeshRevision::simplifyMesh(const std::string &new_mesh_name, double eps, unsigned min_elem_dim)\n+{\n+\tif (this->_mesh.getNElements() == 0)\n+\t\treturn nullptr;\n+\n+\tstd::vector<MeshLib::Node*> new_nodes = this->constructNewNodesArray(this->collapseNodeIndeces(eps));\n+\tstd::vector<MeshLib::Element*> new_elements;\n+\n+\tconst std::vector<MeshLib::Element*> elements (this->_mesh.getElements());\n+\tfor (auto elem = elements.begin(); elem != elements.end(); ++elem)\n+\t{\n+\t\tunsigned n_unique_nodes (this->getNUniqueNodes(*elem));\n+\t\tif (n_unique_nodes == (*elem)->getNNodes() && (*elem)->getDimension() >= min_elem_dim)\n+\t\t{\n+\t\t\tElementErrorCode e ((*elem)->validate());\n+\t\t\tif (e[ElementErrorFlag::NonCoplanar])\n+\t\t\t\tthis->subdivideElement(*elem, new_nodes, new_elements);\n+\t\t\telse\n+\t\t\t\tnew_elements.push_back(copyElement(*elem, new_nodes));\n+\t\t}\n+\t\telse if (n_unique_nodes < (*elem)->getNNodes() && n_unique_nodes>1)\n+\t\t\treduceElement(*elem, n_unique_nodes, new_nodes, new_elements, min_elem_dim);\n+\t\telse\n+\t\t\tstd::cout << \"Error: Something is wrong, more unique nodes than actual nodes\" << std::endl;\n+\n+\t}\n+\tif (!new_elements.empty())\n+\t\treturn new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements);\n+\n+\t// clean up\n+\tfor (auto node = new_nodes.begin(); node != new_nodes.end(); ++node)\n+\t\tdelete *node;\n+\treturn nullptr;\n+}\n+\n+std::vector<std::size_t> MeshRevision::collapseNodeIndeces(double eps) const\n+{\n+\tconst std::vector<MeshLib::Node*> nodes (_mesh.getNodes());\n+\tconst std::size_t nNodes (_mesh.getNNodes());\n+\tstd::vector<std::size_t> id_map(nNodes);\n+\tfor (std::size_t k=0; k<nNodes; ++k)\n+\t\tid_map[k]=k;\n+\n+\tGeoLib::Grid<MeshLib::Node>* grid(new GeoLib::Grid<MeshLib::Node>(nodes.begin(), nodes.end(), 64));\n+\n+\tfor (size_t k=0; k<nNodes; ++k)\n+\t{\n+\t\tstd::vector<std::vector<MeshLib::Node*> const*> node_vectors;\n+\t\tMeshLib::Node const*const node(nodes[k]);\n+\t\tgrid->getPntVecsOfGridCellsIntersectingCube(node->getCoords(), eps, node_vectors);\n+\n+\t\tconst size_t nVectors (node_vectors.size());\n+\t\tfor (size_t i=0; i<nVectors; ++i)\n+\t\t{\n+\t\t\tstd::vector<MeshLib::Node*> cell_vector (*node_vectors[i]);\n+\t\t\tconst size_t nGridCellNodes (cell_vector.size());\n+\t\t\tfor (size_t j=0; j<nGridCellNodes; ++j)\n+\t\t\t{\n+\t\t\t\tMeshLib::Node* test_node(cell_vector[j]);\n+\t\t\t\t// are node indices already identical (i.e. nodes will be collapsed)\n+\t\t\t\tif (id_map[node->getID()] == id_map[test_node->getID()])\n+\t\t\t\t\tcontinue;\n+\n+\t\t\t\t// if test_node has already been collapsed to another node x, ignore it\n+\t\t\t\t// (if the current node would need to be collapsed with x it would already have happened when x was tested)\n+\t\t\t\tif (test_node->getID() != id_map[test_node->getID()])\n+\t\t\t\t\tcontinue;\n+\n+\t\t\t\t// calc distance\n+\t\t\t\tif (MathLib::sqrDist(node->getCoords(), test_node->getCoords()) < eps)\n+\t\t\t\t\tid_map[test_node->getID()] = node->getID();\n+\t\t\t}\n+\t\t}\n+\t}\n+\treturn id_map;\n+}\n+\n+std::vector<MeshLib::Node*> MeshRevision::constructNewNodesArray(const std::vector<std::size_t> &id_map)\n+{\n+\tstd::vector<MeshLib::Node*> nodes (_mesh.getNodes());\n+\tconst std::size_t nNodes (nodes.size());\n+\tstd::vector<MeshLib::Node*> new_nodes;\n+\tfor (std::size_t k=0; k<nNodes; ++k)\n+\t{\n+\t\t// all nodes that have not been collapsed with other nodes are copied into new array\n+\t\tif (nodes[k]->getID() == id_map[k])\n+\t\t{\n+\t\t\tunsigned id = static_cast<unsigned>(new_nodes.size());\n+\t\t\tnew_nodes.push_back (new MeshLib::Node( (*nodes[k])[0], (*nodes[k])[1], (*nodes[k])[2], id) );\n+\t\t\tnodes[k]->setID(id); // the node in the old array gets the index of the same node in the new array\n+\t\t}\n+\t\t// the other nodes are not copied and get the index of the nodes they will have been collapsed with\n+\t\telse\n+\t\t\tnodes[k]->setID(nodes[id_map[k]]->getID());\n+\t}\n+\treturn new_nodes;\n+}\n+\n+unsigned MeshRevision::getNUniqueNodes(MeshLib::Element const*const element) const\n+{\n+\tunsigned nNodes (element->getNNodes());\n+\tunsigned count (nNodes);\n+\n+\tfor (unsigned i=0; i<nNodes-1; ++i)\n+\t\tfor (unsigned j=i+1; j<nNodes; ++j)\n+\t\t\tif (element->getNode(i)->getID()==element->getNode(j)->getID())\n+\t\t\t{\n+\t\t\t\tcount--;\n+\t\t\t\tbreak;\n+\t\t\t}\n+\treturn count;\n+}\n+\n+MeshLib::Element* MeshRevision::copyElement(MeshLib::Element const*const element, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Element* new_elem (nullptr);\n+\tif (element->getGeomType() == MeshElemType::LINE)\n+\t\treturn this->copyLine(element, nodes);\n+\telse if (element->getGeomType() == MeshElemType::TRIANGLE)\n+\t\treturn this->copyTri(element, nodes);\n+\telse if (element->getGeomType() == MeshElemType::QUAD)\n+\t\treturn this->copyQuad(element, nodes);\n+\telse if (element->getGeomType() == MeshElemType::TETRAHEDRON)\n+\t\treturn this->copyTet(element, nodes);\n+\telse if (element->getGeomType() == MeshElemType::HEXAHEDRON)\n+\t\treturn this->copyHex(element, nodes);\n+\telse if (element->getGeomType() == MeshElemType::PYRAMID)\n+\t\treturn this->copyPyramid(element, nodes);\n+\telse if (element->getGeomType() == MeshElemType::PRISM)\n+\t\treturn this->copyPrism(element, nodes);\n+\n+\tERR (\"Error: Unknown element type.\");\n+\treturn nullptr;\n+}\n+\n+void MeshRevision::subdivideElement(MeshLib::Element const*const element, const std::vector<MeshLib::Node*> &nodes, std::vector<MeshLib::Element*> &elements) const\n+{\n+\tif (element->getGeomType() == MeshElemType::QUAD)\n+\t\tthis->subdivideQuad(element, nodes, elements);\n+\telse if (element->getGeomType() == MeshElemType::HEXAHEDRON)\n+\t\tthis->subdivideHex(element, nodes, elements);\n+\telse if (element->getGeomType() == MeshElemType::PYRAMID)\n+\t\tthis->subdividePyramid(element, nodes, elements);\n+\telse if (element->getGeomType() == MeshElemType::PRISM)\n+\t\tthis->subdividePrism(element, nodes, elements);\n+\telse\n+\t\tERR (\"Error: Unknown element type.\");\n+}\n+\n+void MeshRevision::reduceElement(MeshLib::Element const*const element, \n+\t\t\t\t\t\t\t\t unsigned n_unique_nodes, \n+\t\t\t\t\t\t\t\t const std::vector<MeshLib::Node*> &nodes,\n+\t\t\t\t\t\t\t\t std::vector<MeshLib::Element*> &elements,\n+\t\t\t\t\t\t\t\t unsigned min_elem_dim) const\n+{\n+\t/***************\n+\t * TODO: modify neighbouring elements if one elements has been subdivided\n+\t ***************/\n+\tMeshLib::Element* new_elem (nullptr);\n+\tif (element->getGeomType() == MeshElemType::TRIANGLE && min_elem_dim == 1)\n+\t\telements.push_back (this->constructLine(element, nodes));\n+\telse if ((element->getGeomType() == MeshElemType::QUAD) || \n+\t\t     (element->getGeomType() == MeshElemType::TETRAHEDRON))\n+\t{\n+\t\tif (n_unique_nodes == 3 && min_elem_dim < 3)\n+\t\t\telements.push_back (this->constructTri(element, nodes));\n+\t\telse if (min_elem_dim == 1)\n+\t\t\telements.push_back (this->constructLine(element, nodes));\n+\t}\n+\telse if (element->getGeomType() == MeshElemType::HEXAHEDRON)\n+\t\tthis->reduceHex(element, n_unique_nodes, nodes, elements, min_elem_dim);\n+\telse if (element->getGeomType() == MeshElemType::PYRAMID)\n+\t\tthis->reducePyramid(element, n_unique_nodes, nodes, elements, min_elem_dim);\n+\telse if (element->getGeomType() == MeshElemType::PRISM)\n+\t\tthis->reducePrism(element, n_unique_nodes, nodes, elements, min_elem_dim);\n+\telse\n+\t\tERR (\"Error: Unknown element type.\");\n+\treturn;\n+}\n+\n+MeshLib::Element* MeshRevision::copyLine(MeshLib::Element const*const org_elem, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[2];\n+\tnew_nodes[0] = nodes[org_elem->getNode(0)->getID()];\n+\tnew_nodes[1] = nodes[org_elem->getNode(1)->getID()];\n+\treturn new MeshLib::Line(new_nodes, org_elem->getValue());\n+}\n+\n+MeshLib::Element* MeshRevision::copyTri(MeshLib::Element const*const org_elem, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[3];\n+\tfor (unsigned i=0; i<3; ++i)\n+\t\tnew_nodes[i] = nodes[org_elem->getNode(i)->getID()];\n+\treturn new MeshLib::Tri(new_nodes, org_elem->getValue());\n+}\n+\n+MeshLib::Element* MeshRevision::copyQuad(MeshLib::Element const*const org_elem, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[4];\n+\tfor (unsigned i=0; i<4; ++i)\n+\t\tnew_nodes[i] = nodes[org_elem->getNode(i)->getID()];\n+\treturn new MeshLib::Quad(new_nodes, org_elem->getValue());\n+}\n+\n+MeshLib::Element* MeshRevision::copyTet(MeshLib::Element const*const org_elem, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[4];\n+\tfor (unsigned i=0; i<4; ++i)\n+\t\tnew_nodes[i] = nodes[org_elem->getNode(i)->getID()];\n+\treturn new MeshLib::Tet(new_nodes, org_elem->getValue());\n+}\n+\n+MeshLib::Element* MeshRevision::copyHex(MeshLib::Element const*const org_elem, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[8];\n+\tfor (unsigned i=0; i<8; ++i)\n+\t\tnew_nodes[i] = nodes[org_elem->getNode(i)->getID()];\n+\treturn new MeshLib::Hex(new_nodes, org_elem->getValue());\n+}\n+\n+MeshLib::Element* MeshRevision::copyPyramid(MeshLib::Element const*const org_elem, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[5];\n+\tfor (unsigned i=0; i<5; ++i)\n+\t\tnew_nodes[i] = nodes[org_elem->getNode(i)->getID()];\n+\treturn new MeshLib::Pyramid(new_nodes, org_elem->getValue());\n+}\n+\n+MeshLib::Element* MeshRevision::copyPrism(MeshLib::Element const*const org_elem, const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[6];\n+\tfor (unsigned i=0; i<6; ++i)\n+\t\tnew_nodes[i] = nodes[org_elem->getNode(i)->getID()];\n+\treturn new MeshLib::Prism(new_nodes, org_elem->getValue());\n+}\n+\n+unsigned MeshRevision::subdivideQuad(MeshLib::Element const*const quad, const std::vector<MeshLib::Node*> &nodes, std::vector<MeshLib::Element*> &new_elements) const\n+{\n+\tMeshLib::Node** tri1_nodes = new MeshLib::Node*[3];\n+\ttri1_nodes[0] = nodes[quad->getNode(0)->getID()];\n+\ttri1_nodes[1] = nodes[quad->getNode(1)->getID()];\n+\ttri1_nodes[2] = nodes[quad->getNode(2)->getID()];\n+\tnew_elements.push_back(new MeshLib::Tri(tri1_nodes, quad->getValue()));\n+\n+\tMeshLib::Node** tri2_nodes = new MeshLib::Node*[3];\n+\ttri2_nodes[0] = nodes[quad->getNode(0)->getID()];\n+\ttri2_nodes[1] = nodes[quad->getNode(2)->getID()];\n+\ttri2_nodes[2] = nodes[quad->getNode(3)->getID()];\n+\tnew_elements.push_back(new MeshLib::Tri(tri2_nodes, quad->getValue()));\n+\n+\treturn 2;\n+}\n+\n+unsigned MeshRevision::subdivideHex(MeshLib::Element const*const hex, const std::vector<MeshLib::Node*> &nodes, std::vector<MeshLib::Element*> &new_elements) const\n+{\n+\tMeshLib::Node** prism1_nodes = new MeshLib::Node*[6];\n+\tprism1_nodes[0] = nodes[hex->getNode(0)->getID()];\n+\tprism1_nodes[1] = nodes[hex->getNode(2)->getID()];\n+\tprism1_nodes[2] = nodes[hex->getNode(1)->getID()];\n+\tprism1_nodes[3] = nodes[hex->getNode(4)->getID()];\n+\tprism1_nodes[4] = nodes[hex->getNode(6)->getID()];\n+\tprism1_nodes[5] = nodes[hex->getNode(5)->getID()];\n+\tMeshLib::Prism* prism1 (new MeshLib::Prism(prism1_nodes, hex->getValue()));\n+\tthis->subdividePrism(prism1, nodes, new_elements);\n+\tdelete prism1;\n+\n+\tMeshLib::Node** prism2_nodes = new MeshLib::Node*[6];\n+\tprism2_nodes[0] = nodes[hex->getNode(4)->getID()];\n+\tprism2_nodes[1] = nodes[hex->getNode(6)->getID()];\n+\tprism2_nodes[2] = nodes[hex->getNode(7)->getID()];\n+\tprism2_nodes[3] = nodes[hex->getNode(0)->getID()];\n+\tprism2_nodes[4] = nodes[hex->getNode(2)->getID()];\n+\tprism2_nodes[5] = nodes[hex->getNode(3)->getID()];\n+\tMeshLib::Prism* prism2 (new MeshLib::Prism(prism2_nodes, hex->getValue()));\n+\tthis->subdividePrism(prism2, nodes, new_elements);\n+\tdelete prism2;\n+\n+\treturn 6;\n+}\n+\n+unsigned MeshRevision::subdividePyramid(MeshLib::Element const*const pyramid, const std::vector<MeshLib::Node*> &nodes, std::vector<MeshLib::Element*> &new_elements) const\n+{\n+\tMeshLib::Node** tet1_nodes = new MeshLib::Node*[4];\n+\ttet1_nodes[0] = nodes[pyramid->getNode(0)->getID()];\n+\ttet1_nodes[1] = nodes[pyramid->getNode(1)->getID()];\n+\ttet1_nodes[2] = nodes[pyramid->getNode(2)->getID()];\n+\ttet1_nodes[3] = nodes[pyramid->getNode(4)->getID()];\n+\tnew_elements.push_back(new MeshLib::Tet(tet1_nodes, pyramid->getValue()));\n+\n+\tMeshLib::Node** tet2_nodes = new MeshLib::Node*[4];\n+\ttet2_nodes[0] = nodes[pyramid->getNode(0)->getID()];\n+\ttet2_nodes[1] = nodes[pyramid->getNode(2)->getID()];\n+\ttet2_nodes[2] = nodes[pyramid->getNode(3)->getID()];\n+\ttet2_nodes[3] = nodes[pyramid->getNode(4)->getID()];\n+\tnew_elements.push_back(new MeshLib::Tet(tet2_nodes, pyramid->getValue()));\n+\n+\treturn 2;\n+}\n+\n+unsigned MeshRevision::subdividePrism(MeshLib::Element const*const prism, const std::vector<MeshLib::Node*> &nodes, std::vector<MeshLib::Element*> &new_elements) const\n+{\n+\tMeshLib::Node** tet1_nodes = new MeshLib::Node*[4];\n+\ttet1_nodes[0] = nodes[prism->getNode(0)->getID()];\n+\ttet1_nodes[1] = nodes[prism->getNode(1)->getID()];\n+\ttet1_nodes[2] = nodes[prism->getNode(2)->getID()];\n+\ttet1_nodes[3] = nodes[prism->getNode(3)->getID()];\n+\tnew_elements.push_back(new MeshLib::Tet(tet1_nodes, prism->getValue()));\n+\n+\tMeshLib::Node** tet2_nodes = new MeshLib::Node*[4];\n+\ttet2_nodes[0] = nodes[prism->getNode(3)->getID()];\n+\ttet2_nodes[1] = nodes[prism->getNode(2)->getID()];\n+\ttet2_nodes[2] = nodes[prism->getNode(4)->getID()];\n+\ttet2_nodes[3] = nodes[prism->getNode(5)->getID()];\n+\tnew_elements.push_back(new MeshLib::Tet(tet2_nodes, prism->getValue()));\n+\n+\tMeshLib::Node** tet3_nodes = new MeshLib::Node*[4];\n+\ttet3_nodes[0] = nodes[prism->getNode(2)->getID()];\n+\ttet3_nodes[1] = nodes[prism->getNode(1)->getID()];\n+\ttet3_nodes[2] = nodes[prism->getNode(3)->getID()];\n+\ttet3_nodes[3] = nodes[prism->getNode(4)->getID()];\n+\tnew_elements.push_back(new MeshLib::Tet(tet3_nodes, prism->getValue()));\n+\n+\treturn 3;\n+}\n+\n+unsigned MeshRevision::reduceHex(MeshLib::Element const*const org_elem, \n+\t\t\t\t\t\t\t unsigned n_unique_nodes, \n+\t\t\t\t\t\t\t const std::vector<MeshLib::Node*> &nodes,\n+\t\t\t\t\t\t\t std::vector<MeshLib::Element*> &new_elements,\n+\t\t\t\t\t\t\t unsigned min_elem_dim) const\n+{\n+\t// TODO?\n+\t// if two diametral nodes collapse, all kinds of bizarre (2D-)element combinations could be the result.\n+\t// this case is currently not implemented!\n+\n+\tif (n_unique_nodes == 7)\n+\t{\n+\t\t// reduce to prism + pyramid\n+\t\tfor (unsigned i=0; i<7; ++i)\n+\t\t\tfor (unsigned j=i+1; j<8; ++j)\n+\t\t\t\tif (org_elem->getNode(i)->getID() == org_elem->getNode(j)->getID())\n+\t\t\t\t{\n+\t\t\t\t\tconst std::array<unsigned,4> base_nodes (this->lutHexCuttingQuadNodes(i,j));\n+\t\t\t\t\tMeshLib::Node** pyr_nodes = new MeshLib::Node*[5];\n+\t\t\t\t\tpyr_nodes[0] = nodes[org_elem->getNode(base_nodes[0])->getID()];\n+\t\t\t\t\tpyr_nodes[1] = nodes[org_elem->getNode(base_nodes[1])->getID()];\n+\t\t\t\t\tpyr_nodes[2] = nodes[org_elem->getNode(base_nodes[2])->getID()];\n+\t\t\t\t\tpyr_nodes[3] = nodes[org_elem->getNode(base_nodes[3])->getID()];\n+\t\t\t\t\tpyr_nodes[4] = nodes[org_elem->getNode(i)->getID()];\n+\t\t\t\t\tnew_elements.push_back (new MeshLib::Pyramid(pyr_nodes, org_elem->getValue()));\n+\n+\t\t\t\t\tif (i<4 && j>=4) std::swap(i,j);\n+\t\t\t\t\tMeshLib::Node** prism_nodes = new MeshLib::Node*[6];\n+\t\t\t\t\tprism_nodes[0] = nodes[org_elem->getNode(base_nodes[0])->getID()];\n+\t\t\t\t\tprism_nodes[1] = nodes[org_elem->getNode(base_nodes[3])->getID()];\n+\t\t\t\t\tprism_nodes[2] = nodes[org_elem->getNode(this->lutHexDiametralNode(j))->getID()];\n+\t\t\t\t\tprism_nodes[3] = nodes[org_elem->getNode(base_nodes[1])->getID()];\n+\t\t\t\t\tprism_nodes[4] = nodes[org_elem->getNode(base_nodes[2])->getID()];\n+\t\t\t\t\tprism_nodes[5] = nodes[org_elem->getNode(this->lutHexDiametralNode(i))->getID()];\n+\t\t\t\t\tnew_elements.push_back (new MeshLib::Prism(prism_nodes, org_elem->getValue()));\n+\t\t\t\t\treturn 2;\n+\t\t\t\t}\n+\t}\n+\telse if (n_unique_nodes == 6)\n+\t{\n+\t\t// reduce to prism\n+\t\tfor (unsigned i=0; i<6; ++i)\n+\t\t{\n+\t\t\tconst MeshLib::Element* face (org_elem->getFace(i));\n+\t\t\tif (face->getNode(0)->getID() == face->getNode(1)->getID() && face->getNode(2)->getID() == face->getNode(3)->getID())\n+\t\t\t{\n+\t\t\t\tMeshLib::Node** prism_nodes = new MeshLib::Node*[6];\n+\t\t\t\tprism_nodes[0] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(0))))->getID()];\n+\t\t\t\tprism_nodes[1] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(1))))->getID()];\n+\t\t\t\tprism_nodes[2] = nodes[org_elem->getNode(org_elem->getNodeIDinElement(face->getNode(2)))->getID()];\n+\t\t\t\tprism_nodes[3] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(2))))->getID()];\n+\t\t\t\tprism_nodes[4] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(3))))->getID()];\n+\t\t\t\tprism_nodes[5] = nodes[org_elem->getNode(org_elem->getNodeIDinElement(face->getNode(0)))->getID()];\n+\t\t\t\tnew_elements.push_back (new MeshLib::Prism(prism_nodes, org_elem->getValue()));\n+\t\t\t\treturn 1;\n+\t\t\t}\n+\t\t\tif (face->getNode(0)->getID() == face->getNode(3)->getID() && face->getNode(1)->getID() == face->getNode(2)->getID())\n+\t\t\t{\n+\t\t\t\tMeshLib::Node** prism_nodes = new MeshLib::Node*[6];\n+\t\t\t\tprism_nodes[0] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(0))))->getID()];\n+\t\t\t\tprism_nodes[1] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(3))))->getID()];\n+\t\t\t\tprism_nodes[2] = nodes[org_elem->getNode(org_elem->getNodeIDinElement(face->getNode(2)))->getID()];\n+\t\t\t\tprism_nodes[3] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(1))))->getID()];\n+\t\t\t\tprism_nodes[4] = nodes[org_elem->getNode(this->lutHexDiametralNode(org_elem->getNodeIDinElement(face->getNode(2))))->getID()];\n+\t\t\t\tprism_nodes[5] = nodes[org_elem->getNode(org_elem->getNodeIDinElement(face->getNode(0)))->getID()];\n+\t\t\t\treturn 1;\n+\t\t\t}\n+\t\t\tdelete face;\n+\t\t}\n+\t\t// reduce to four tets -> divide into 2 prisms such that each has one collapsed node\n+\t\tfor (unsigned i=0; i<7; ++i)\n+\t\t\tfor (unsigned j=i+1; j<8; ++j)\n+\t\t\tif (org_elem->getNode(i)->getID() == org_elem->getNode(j)->getID())\n+\t\t\t{\n+\t\t\t\tfor (unsigned k=i; k<7; ++k)\n+\t\t\t\t\tfor (unsigned l=k+1; l<8; ++l)\n+\t\t\t\t\t\tif (!(i==k && j==l) && org_elem->isEdge(i,j) && org_elem->isEdge(k,l) && \n+\t\t\t\t\t\t\torg_elem->getNode(k)->getID() == org_elem->getNode(l)->getID())\n+\t\t\t\t\t\t{\n+\t\t\t\t\t\t\tconst std::pair<unsigned, unsigned> back (this->lutHexBackNodes(i,j,k,l));\n+\t\t\t\t\t\t\tif (back.first == std::numeric_limits<unsigned>::max() || back.second == std::numeric_limits<unsigned>::max())\n+\t\t\t\t\t\t\t{\n+\t\t\t\t\t\t\t\tERR (\"Unexpected error during Hex reduction\");\n+\t\t\t\t\t\t\t\treturn 0;\n+\t\t\t\t\t\t\t}\n+\n+\t\t\t\t\t\t\tstd::array<unsigned, 4> cutting_plane (this->lutHexCuttingQuadNodes(back.first, back.second));\n+\t\t\t\t\t\t\tMeshLib::Node** pris1_nodes = new MeshLib::Node*[6];\n+\t\t\t\t\t\t\tpris1_nodes[0] = const_cast<MeshLib::Node*>(org_elem->getNode(back.first));\n+\t\t\t\t\t\t\tpris1_nodes[1] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[0])); \n+\t\t\t\t\t\t\tpris1_nodes[2] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[3]));\n+\t\t\t\t\t\t\tpris1_nodes[3] = const_cast<MeshLib::Node*>(org_elem->getNode(back.second)); \n+\t\t\t\t\t\t\tpris1_nodes[4] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[1])); \n+\t\t\t\t\t\t\tpris1_nodes[5] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[2]));\n+\t\t\t\t\t\t\tMeshLib::Prism* prism1 (new MeshLib::Prism(pris1_nodes, org_elem->getValue()));\n+\t\t\t\t\t\t\tunsigned nNewElements = this->reducePrism(prism1, 5, nodes, new_elements, min_elem_dim);\n+\t\t\t\t\t\t\tdelete prism1;\n+\n+\t\t\t\t\t\t\tMeshLib::Node** pris2_nodes = new MeshLib::Node*[6];\n+\t\t\t\t\t\t\tpris2_nodes[0] = const_cast<MeshLib::Node*>(org_elem->getNode(this->lutHexDiametralNode(back.first)));\n+\t\t\t\t\t\t\tpris2_nodes[1] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[0]));\n+\t\t\t\t\t\t\tpris2_nodes[2] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[3]));\n+\t\t\t\t\t\t\tpris2_nodes[3] = const_cast<MeshLib::Node*>(org_elem->getNode(this->lutHexDiametralNode(back.second)));\n+\t\t\t\t\t\t\tpris2_nodes[4] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[1]));\n+\t\t\t\t\t\t\tpris2_nodes[5] = const_cast<MeshLib::Node*>(org_elem->getNode(cutting_plane[2]));\n+\t\t\t\t\t\t\tMeshLib::Prism* prism2 (new MeshLib::Prism(pris2_nodes, org_elem->getValue()));\n+\t\t\t\t\t\t\tnNewElements += this->reducePrism(prism2, 5, nodes, new_elements, min_elem_dim);\n+\t\t\t\t\t\t\tdelete prism2;\n+\t\t\t\t\t\t\treturn nNewElements;\n+\t\t\t\t\t\t}\n+\t\t\t}\n+\t}\n+\telse if (n_unique_nodes == 5)\n+\t{\n+\t\tMeshLib::Element* tet1 (this->constructFourNodeElement(org_elem, nodes));\n+\t\tstd::array<unsigned, 4> first_four_nodes = { tet1->getNode(0)->getID(), tet1->getNode(1)->getID(), tet1->getNode(2)->getID(), tet1->getNode(3)->getID() };\n+\t\tunsigned fifth_node (this->findPyramidTopNode(*org_elem, first_four_nodes));\n+\n+\t\tbool tet_changed (false);\n+\t\tif (tet1->getGeomType() == MeshElemType::QUAD)\n+\t\t{\n+\t\t\tdelete tet1;\n+\t\t\ttet_changed =true;\n+\t\t\tMeshLib::Node** tet1_nodes = new MeshLib::Node*[4];\n+\t\t\ttet1_nodes[0] = nodes[first_four_nodes[0]];\n+\t\t\ttet1_nodes[1] = nodes[first_four_nodes[1]];\n+\t\t\ttet1_nodes[2] = nodes[first_four_nodes[2]];\n+\t\t\ttet1_nodes[3] = nodes[org_elem->getNode(fifth_node)->getID()];\n+\t\t\tnew_elements.push_back(new MeshLib::Tet(tet1_nodes, org_elem->getValue()));\n+\t\t}\n+\t\telse \n+\t\t\tnew_elements.push_back(tet1);\n+\n+\t\tMeshLib::Node** tet2_nodes = new MeshLib::Node*[4];\n+\t\ttet2_nodes[0] = (tet_changed) ? nodes[first_four_nodes[0]] : nodes[first_four_nodes[1]];\n+\t\ttet2_nodes[1] = nodes[first_four_nodes[2]];\n+\t\ttet2_nodes[2] = nodes[first_four_nodes[3]];\n+\t\ttet2_nodes[3] = nodes[org_elem->getNode(fifth_node)->getID()];\n+\t\tnew_elements.push_back(new MeshLib::Tet(tet2_nodes, org_elem->getValue()));\n+\t\treturn 2;\n+\t}\n+\telse if (n_unique_nodes == 4)\n+\t{\n+\t\tMeshLib::Element* elem (this->constructFourNodeElement(org_elem, nodes, min_elem_dim));\n+\t\tif (elem)\n+\t\t\tnew_elements.push_back (elem);\n+\t}\n+\telse if (n_unique_nodes == 3 && min_elem_dim < 3)\n+\t\tnew_elements.push_back (this->constructTri(org_elem, nodes));\n+\telse if (min_elem_dim == 1)\n+\t\tnew_elements.push_back (this->constructLine(org_elem, nodes));\n+\treturn 1;\n+}\n+\n+void MeshRevision::reducePyramid(MeshLib::Element const*const org_elem, \n+\t\t\t\t\t\t\t\t unsigned n_unique_nodes, \n+\t\t\t\t\t\t\t\t const std::vector<MeshLib::Node*> &nodes, \n+\t\t\t\t\t\t\t\t std::vector<MeshLib::Element*> &new_elements,\n+\t\t\t\t\t\t\t\t unsigned min_elem_dim) const\n+{\n+\tif (n_unique_nodes == 4)\n+\t{\n+\t\tMeshLib::Element* elem (this->constructFourNodeElement(org_elem, nodes, min_elem_dim));\n+\t\tif (elem)\n+\t\t\tnew_elements.push_back (elem);\n+\t}\n+\telse if (n_unique_nodes == 3 && min_elem_dim < 3)\n+\t\tnew_elements.push_back (this->constructTri(org_elem, nodes));\n+\telse if (n_unique_nodes == 2 && min_elem_dim == 1)\n+\t\tnew_elements.push_back (this->constructLine(org_elem, nodes));\n+\treturn;\n+}\n+\n+unsigned MeshRevision::reducePrism(MeshLib::Element const*const org_elem, \n+\t\t\t\t\t\t\t   unsigned n_unique_nodes, \n+\t\t\t\t\t\t\t   const std::vector<MeshLib::Node*> &nodes, \n+\t\t\t\t\t\t\t   std::vector<MeshLib::Element*> &new_elements,\n+\t\t\t\t\t\t\t   unsigned min_elem_dim) const\n+{\n+\t// TODO?\n+\t// In theory a node from the bottom triangle and a node from the top triangle that are not connected by an edge\n+\t// could collapse, resulting in a combination of tri and quad elements. This case is currently not tested.\n+\n+\t// if one of the non-triangle edges collapsed, elem can be reduced to a pyramid, otherwise it will be two tets\n+\tif (n_unique_nodes == 5)\n+\t{\n+\t\tfor (unsigned i=0; i<5; ++i)\n+\t\t\tfor (unsigned j=i+1; j<6; ++j)\n+\t\t\t\tif (i!=j && org_elem->getNode(i)->getID() == org_elem->getNode(j)->getID())\n+\t\t\t\t{\n+\t\t\t\t\t// non triangle edge collapsed\n+\t\t\t\t\tif (i%3 == j%3)\n+\t\t\t\t\t{\n+\t\t\t\t\t\tMeshLib::Node** tet1_nodes = new MeshLib::Node*[4];\n+\t\t\t\t\t\ttet1_nodes[0] = nodes[org_elem->getNode((i+1)%3)->getID()];\n+\t\t\t\t\t\ttet1_nodes[1] = nodes[org_elem->getNode((i+2)%3)->getID()];\n+\t\t\t\t\t\ttet1_nodes[2] = nodes[org_elem->getNode(i)->getID()];\n+\t\t\t\t\t\ttet1_nodes[3] = nodes[org_elem->getNode((i+1)%3+3)->getID()];\n+\t\t\t\t\t\tnew_elements.push_back (new MeshLib::Tet(tet1_nodes, org_elem->getValue()));\n+\n+\t\t\t\t\t\tMeshLib::Node** tet2_nodes = new MeshLib::Node*[4];\n+\t\t\t\t\t\ttet2_nodes[0] = nodes[org_elem->getNode((i+1)%3+3)->getID()];\n+\t\t\t\t\t\ttet2_nodes[1] = nodes[org_elem->getNode((i+2)%3)->getID()];\n+\t\t\t\t\t\ttet2_nodes[2] = nodes[org_elem->getNode(i)->getID()];\t\t\t\t\t\n+\t\t\t\t\t\ttet2_nodes[3] = nodes[org_elem->getNode((i+2)%3+3)->getID()];\n+\t\t\t\t\t\tnew_elements.push_back (new MeshLib::Tet(tet2_nodes, org_elem->getValue()));\n+\t\t\t\t\t\treturn 2;\n+\t\t\t\t\t}\n+\n+\t\t\t\t\t// triangle edge collapsed\n+\t\t\t\t\tconst int offset = (i>2) ? -3 : 3;\n+\t\t\t\t\tconst unsigned k = this->lutPrismThirdNode(i,j);\n+\t\t\t\t\tif (k == std::numeric_limits<unsigned>::max())\n+\t\t\t\t\t{\n+\t\t\t\t\t\tERR (\"Unexpected error during prism reduction.\")\n+\t\t\t\t\t\treturn 0;\n+\t\t\t\t\t}\n+\t\t\t\t\tMeshLib::Node** tet1_nodes = new MeshLib::Node*[4];\n+\t\t\t\t\ttet1_nodes[0] = nodes[org_elem->getNode(i+offset)->getID()];\n+\t\t\t\t\ttet1_nodes[1] = nodes[org_elem->getNode(j+offset)->getID()];\n+\t\t\t\t\ttet1_nodes[2] = nodes[org_elem->getNode(k+offset)->getID()];\n+\t\t\t\t\ttet1_nodes[3] = nodes[org_elem->getNode(i)->getID()];\n+\t\t\t\t\tnew_elements.push_back (new MeshLib::Tet(tet1_nodes, org_elem->getValue()));\n+\n+\t\t\t\t\tunsigned l = (GeoLib::isCoplanar(*org_elem->getNode(i+offset), *org_elem->getNode(k+offset), *org_elem->getNode(i), *org_elem->getNode(k))) ? j : i;\n+\t\t\t\t\tMeshLib::Node** tet2_nodes = new MeshLib::Node*[4];\n+\t\t\t\t\ttet2_nodes[0] = nodes[org_elem->getNode(l+offset)->getID()];\n+\t\t\t\t\ttet2_nodes[1] = nodes[org_elem->getNode(k+offset)->getID()];\n+\t\t\t\t\ttet2_nodes[2] = nodes[org_elem->getNode(i)->getID()];\t\t\t\t\t\n+\t\t\t\t\ttet2_nodes[3] = nodes[org_elem->getNode(k)->getID()];\n+\t\t\t\t\tnew_elements.push_back (new MeshLib::Tet(tet2_nodes, org_elem->getValue()));\n+\t\t\t\t\treturn 2;\n+\t\t\t\t}\n+\t}\n+\telse if (n_unique_nodes == 4)\n+\t{\n+\t\tMeshLib::Element* elem (this->constructFourNodeElement(org_elem, nodes, min_elem_dim));\n+\t\tif (elem)\n+\t\t\tnew_elements.push_back (elem);\n+\t}\n+\telse if (n_unique_nodes == 3 && min_elem_dim < 3)\n+\t\tnew_elements.push_back (this->constructTri(org_elem, nodes));\n+\telse if (n_unique_nodes == 2 && min_elem_dim == 1)\n+\t\tnew_elements.push_back (this->constructLine(org_elem, nodes));\n+\treturn 1;\n+}\n+\n+MeshLib::Element* MeshRevision::constructLine(MeshLib::Element const*const element, \n+\t\t\t\t\t\t\t\t\t\t\t  const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** line_nodes = new MeshLib::Node*[2];\n+\tline_nodes[0] = nodes[element->getNode(0)->getID()];\n+\tfor (unsigned i=1; i<element->getNNodes(); ++i)\n+\t{\n+\t\tif (element->getNode(i)->getID() != element->getNode(0)->getID())\n+\t\t{\n+\t\t\tline_nodes[1] = nodes[element->getNode(i)->getID()];\n+\t\t\tbreak;\n+\t\t}\n+\t}\n+\treturn new MeshLib::Line(line_nodes, element->getValue());\n+}\n+\n+MeshLib::Element* MeshRevision::constructTri(MeshLib::Element const*const element, \n+\t\t\t\t\t\t\t\t\t\t\t const std::vector<MeshLib::Node*> &nodes) const\n+{\n+\tMeshLib::Node** tri_nodes = new MeshLib::Node*[3];\n+\ttri_nodes[0] = nodes[element->getNode(0)->getID()];\n+\tfor (unsigned i=1; i<element->getNNodes(); ++i)\n+\t\tif (element->getNode(i)->getID() != tri_nodes[0]->getID())\n+\t\t{\n+\t\t\ttri_nodes[1] = nodes[element->getNode(i)->getID()];\n+\t\t\tfor (unsigned j=i+1; j<element->getNNodes(); ++j)\n+\t\t\t\tif (element->getNode(j)->getID() != tri_nodes[1]->getID())\n+\t\t\t\t{\n+\t\t\t\t\ttri_nodes[2] = nodes[element->getNode(j)->getID()];\n+\t\t\t\t\t// TODO? \n+\t\t\t\t\t// In theory three unique nodes could also be reduced to two lines e.g. with  \n+\t\t\t\t\t// a quad where two diametral nodes collapse. This case is currently not implemented!\n+\t\t\t\t\treturn new MeshLib::Tri(tri_nodes, element->getValue());\n+\t\t\t\t}\n+\t\t}\n+\t// this should never be reached\n+\treturn nullptr;\n+}\n+\n+MeshLib::Element* MeshRevision::constructFourNodeElement(MeshLib::Element const*const element, const std::vector<MeshLib::Node*> &nodes, unsigned min_elem_dim) const\n+{\n+\tMeshLib::Node** new_nodes = new MeshLib::Node*[4];\n+\tunsigned count(0);\n+\tnew_nodes[count++] = nodes[element->getNode(0)->getID()];\n+\tfor (unsigned i=1; i<element->getNNodes(); ++i)\n+\t{\n+\t\tif (count>3)\n+\t\t\tbreak;\n+\t\tbool unique_node (true);\n+\t\tfor (unsigned j=0; j<i; ++j)\n+\t\t{\n+\t\t\tif (element->getNode(i)->getID() == element->getNode(j)->getID())\n+\t\t\t{\n+\t\t\t\tunique_node = false;\n+\t\t\t\tbreak;\n+\t\t\t}\n+\t\t}\n+\t\tif (unique_node)\n+\t\t\tnew_nodes[count++] = nodes[element->getNode(i)->getID()];; \n+\t}\n+\t\n+\t// test if quad or tet\n+\tconst bool isQuad (GeoLib::isCoplanar(*new_nodes[0], *new_nodes[1], *new_nodes[2], *new_nodes[3]));\n+\tif (isQuad && min_elem_dim < 3)\n+\t{\n+\t\tMeshLib::Element* elem (new MeshLib::Quad(new_nodes, element->getValue()));\n+\t\tfor (unsigned i=1; i<3; ++i)\n+\t\t{\n+\t\t\tif (elem->validate().none())\n+\t\t\t\treturn elem;\n+\t\t\telse\n+\t\t\t{\n+\t\t\t\t// change node order if not convex\n+\t\t\t\tMeshLib::Node* tmp = new_nodes[i+1];\n+\t\t\t\tnew_nodes[i+1] = new_nodes[i];\n+\t\t\t\tnew_nodes[i] = tmp;\n+\t\t\t}\n+\t\t}\n+\t\treturn elem;\n+\t}\n+\telse if (!isQuad)\n+\t\treturn new MeshLib::Tet(new_nodes, element->getValue());\n+\telse // is quad but min elem dim == 3\n+\t\treturn nullptr;\n+}\n+\n+unsigned MeshRevision::findPyramidTopNode(const MeshLib::Element &element, const std::array<unsigned,4> &base_node_ids) const\n+{\n+\tconst std::size_t nNodes (element.getNNodes());\n+\tfor (std::size_t i=0; i<nNodes; ++i)\n+\t{\n+\t\tbool top_node=true;\n+\t\tfor (unsigned j=0; j<4; ++j)\n+\t\t\tif (element.getNode(i)->getID() == base_node_ids[j])\n+\t\t\t\ttop_node=false;\n+\t\tif (top_node)\n+\t\t\treturn i;\n+\t}\n+\treturn std::numeric_limits<unsigned>::max(); // should never be reached if called correctly\n+}\n+\n+unsigned MeshRevision::lutHexDiametralNode(unsigned id) const\n+{\n+\tswitch (id)\n+\t{\n+\t\tcase 0: return 6;","path":"MeshLib/MeshEditing/MeshRevision.cpp","position":766,"original_position":742,"commit_id":"929cde23390c868087ab8355ebb8067f8f74efac","original_commit_id":"7047b7dda1b39b3cc7b08a85d646aa42b877a21a","user":{"login":"rinkk","id":798445,"avatar_url":"https://gravatar.com/avatar/4a3a2e3b8f5a29db71fc45f7fa47e7df?d=https%3A%2F%2Fidenticons.github.com%2Fe5743c1d0e506418ae0b86fa09134651.png&r=x","gravatar_id":"4a3a2e3b8f5a29db71fc45f7fa47e7df","url":"https://api.github.com/users/rinkk","html_url":"https://github.com/rinkk","followers_url":"https://api.github.com/users/rinkk/followers","following_url":"https://api.github.com/users/rinkk/following{/other_user}","gists_url":"https://api.github.com/users/rinkk/gists{/gist_id}","starred_url":"https://api.github.com/users/rinkk/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/rinkk/subscriptions","organizations_url":"https://api.github.com/users/rinkk/orgs","repos_url":"https://api.github.com/users/rinkk/repos","events_url":"https://api.github.com/users/rinkk/events{/privacy}","received_events_url":"https://api.github.com/users/rinkk/received_events","type":"User","site_admin":false},"body":"good idea.","created_at":"2014-03-12T08:09:31Z","updated_at":"2014-03-12T08:09:31Z","html_url":"https://github.com/ufz/ogs/pull/344#discussion_r10508440","pull_request_url":"https://api.github.com/repos/ufz/ogs/pulls/344","_links":{"self":{"href":"https://api.github.com/repos/ufz/ogs/pulls/comments/10508440"},"html":{"href":"https://github.com/ufz/ogs/pull/344#discussion_r10508440"},"pull_request":{"href":"https://api.github.com/repos/ufz/ogs/pulls/344"}}}}