diff --git a/MeshLib/MeshGenerators/QuadraticMeshGenerator.cpp b/MeshLib/MeshGenerators/QuadraticMeshGenerator.cpp
index 5d6b68ccbe68d7fd59997eac57bba4439d519c74..27851860882e061581fef71be6689b7259ad6a9b 100644
--- a/MeshLib/MeshGenerators/QuadraticMeshGenerator.cpp
+++ b/MeshLib/MeshGenerators/QuadraticMeshGenerator.cpp
@@ -9,140 +9,133 @@
#include "QuadraticMeshGenerator.h"
-#include "BaseLib/makeVectorUnique.h"
-
-#include "MeshLib/Node.h"
#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Elements/Hex.h"
#include "MeshLib/Elements/Line.h"
-#include "MeshLib/Elements/Tri.h"
#include "MeshLib/Elements/Quad.h"
-#include "MeshLib/Elements/Hex.h"
+#include "MeshLib/Elements/Tri.h"
#include "MeshLib/MeshEditing/DuplicateMeshComponents.h"
-#include "MeshLib/Properties.h"
-#include "MeshLib/PropertyVector.h"
-
-namespace MeshLib
-{
-namespace
-{
+#include "MeshLib/Node.h"
-struct Edge
+/// Given an (linear) element divide all its edges by inserting a point in the
+/// middle and return a new element.
+template <typename QuadraticElement>
+std::unique_ptr<QuadraticElement> convertLinearToQuadratic(
+ MeshLib::Element const& e)
{
- Edge(std::size_t id1, std::size_t id2)
- : _id1(std::min(id1, id2)), _id2(std::max(id1, id2))
- {}
-
- bool operator==(Edge const& r) const
+ auto const n_all_nodes = QuadraticElement::n_all_nodes;
+ auto const n_base_nodes = QuadraticElement::n_base_nodes;
+ assert(n_base_nodes == e.getNumberOfBaseNodes());
+
+ // Copy base nodes of element to the quadratic element new nodes'.
+ std::array<MeshLib::Node*, n_all_nodes> nodes;
+ for (int i = 0; i < n_base_nodes; i++)
+ nodes[i] = const_cast<MeshLib::Node*>(e.getNode(i));
+
+ // For each edge create a middle node.
+ int const number_of_edges = e.getNumberOfEdges();
+ for (int i = 0; i < number_of_edges; i++)
{
- return (_id1 == r._id1 && _id2 == r._id2);
- }
-
- std::size_t _id1;
- std::size_t _id2;
+ auto const& a = *e.getEdgeNode(i, 0);
+ auto const& b = *e.getEdgeNode(i, 1);
- std::size_t _edge_id = 0;
-};
+ nodes[n_base_nodes + i] = new MeshLib::Node(
+ (a[0] + b[0]) / 2, (a[1] + b[1]) / 2, (a[2] + b[2]) / 2);
+ }
-bool operator< (Edge const& l, Edge const& r)
-{
- return (l._id1 != r._id1) ? (l._id1 < r._id1) : l._id2 < r._id2;
+ return std::make_unique<QuadraticElement>(nodes, e.getID());
}
-template <typename T_ELEMENT>
-T_ELEMENT* createQuadraticElement(
- MeshLib::Element const* e, std::vector<Edge> const& vec_edges,
- std::vector<MeshLib::Node*> const& vec_new_nodes,
- const std::size_t n_mesh_base_nodes)
+/// Return a new quadratic element corresponding to the linear element's type.
+std::unique_ptr<MeshLib::Element> createQuadraticElement(
+ MeshLib::Element const& e)
{
- auto const n_all_nodes = T_ELEMENT::n_all_nodes;
- auto const n_base_nodes = T_ELEMENT::n_base_nodes;
- auto** nodes = new MeshLib::Node*[n_all_nodes];
- for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
- nodes[i] =
- const_cast<MeshLib::Node*>(vec_new_nodes[e->getNode(i)->getID()]);
- for (unsigned i = 0; i < e->getNumberOfEdges(); i++)
+ if (e.getCellType() == MeshLib::CellType::LINE2)
{
- auto itr = std::find(
- vec_edges.begin(), vec_edges.end(),
- Edge(e->getEdgeNode(i, 0)->getID(), e->getEdgeNode(i, 1)->getID()));
- assert(itr != vec_edges.end());
- nodes[n_base_nodes + i] =
- vec_new_nodes[n_mesh_base_nodes + itr->_edge_id];
+ return convertLinearToQuadratic<MeshLib::Line3>(e);
+ }
+ if (e.getCellType() == MeshLib::CellType::TRI3)
+ {
+ return convertLinearToQuadratic<MeshLib::Tri6>(e);
+ }
+ if (e.getCellType() == MeshLib::CellType::QUAD4)
+ {
+ return convertLinearToQuadratic<MeshLib::Quad8>(e);
+ }
+ if (e.getCellType() == MeshLib::CellType::HEX8)
+ {
+ return convertLinearToQuadratic<MeshLib::Hex20>(e);
}
- return new T_ELEMENT(nodes);
-}
-} // no named namespace
+ OGS_FATAL("Mesh element type %s is not supported",
+ MeshLib::CellType2String(e.getCellType()).c_str());
+}
-std::unique_ptr<Mesh> createQuadraticOrderMesh(Mesh const& org_mesh)
+struct nodeByCoordinatesComparator
{
- std::vector<MeshLib::Node*> vec_new_nodes = MeshLib::copyNodeVector(org_mesh.getNodes());
-
- // collect edges
- std::vector<Edge> vec_edges;
- for (MeshLib::Element const* e : org_mesh.getElements())
- {
- for (unsigned i=0; i<e->getNumberOfEdges(); i++)
- {
- auto node0 = e->getEdgeNode(i, 0);
- auto node1 = e->getEdgeNode(i, 1);
- vec_edges.emplace_back(node0->getID(), node1->getID());
- }
- }
- BaseLib::makeVectorUnique(vec_edges);
- for (std::size_t i=0; i<vec_edges.size(); i++)
- vec_edges[i]._edge_id = i;
- INFO("Found %d edges in the mesh", vec_edges.size());
-
- // create mid-point nodes
- double coords[3];
- for (Edge const& edge : vec_edges)
+ bool operator()(MeshLib::Node* a, MeshLib::Node* b) const
{
- auto const& node0 = *vec_new_nodes[edge._id1];
- auto const& node1 = *vec_new_nodes[edge._id2];
- for (unsigned i=0; i<3; i++)
- coords[i] = (node0[i] + node1[i]) * 0.5;
- vec_new_nodes.push_back(new MeshLib::Node(coords));
+ return *a < *b;
}
+};
+
+namespace MeshLib
+{
+std::unique_ptr<Mesh> createQuadraticOrderMesh(Mesh const& linear_mesh)
+{
+ // Clone the linear mesh nodes.
+ auto quadratic_mesh_nodes = MeshLib::copyNodeVector(linear_mesh.getNodes());
- // create new elements with the quadratic nodes
- std::vector<MeshLib::Element*> vec_new_eles;
- for (MeshLib::Element const* e : org_mesh.getElements())
+ // Temporary container for unique quadratic nodes with O(log(n)) search.
+ std::set<MeshLib::Node*, nodeByCoordinatesComparator> unique_nodes;
+
+ // Create new elements with the quadratic nodes
+ std::vector<MeshLib::Element*> quadratic_elements;
+ auto const& linear_mesh_elements = linear_mesh.getElements();
+ for (MeshLib::Element const* e : linear_mesh_elements)
{
- if (e->getCellType() == MeshLib::CellType::LINE2)
- {
- vec_new_eles.push_back(createQuadraticElement<MeshLib::Line3>(
- e, vec_edges, vec_new_nodes, org_mesh.getNumberOfNodes()));
- }
- else if (e->getCellType() == MeshLib::CellType::TRI3)
- {
- vec_new_eles.push_back(createQuadraticElement<MeshLib::Tri6>(
- e, vec_edges, vec_new_nodes, org_mesh.getNumberOfNodes()));
- }
- else if (e->getCellType() == MeshLib::CellType::QUAD4)
- {
- vec_new_eles.push_back(createQuadraticElement<MeshLib::Quad8>(
- e, vec_edges, vec_new_nodes, org_mesh.getNumberOfNodes()));
- }
- else if (e->getCellType() == MeshLib::CellType::HEX8)
+ auto quadratic_element = createQuadraticElement(*e);
+
+ // Replace the base nodes with cloned linear nodes.
+ int const number_base_nodes = quadratic_element->getNumberOfBaseNodes();
+ for (int i = 0; i < number_base_nodes; ++i)
{
- vec_new_eles.push_back(createQuadraticElement<MeshLib::Hex20>(
- e, vec_edges, vec_new_nodes, org_mesh.getNumberOfNodes()));
+ quadratic_element->setNode(
+ i, quadratic_mesh_nodes[quadratic_element->getNodeIndex(i)]);
}
- else
+
+ // Make the new (middle-edge) nodes unique.
+ int const number_all_nodes = quadratic_element->getNumberOfNodes();
+ for (int i = number_base_nodes; i < number_all_nodes; ++i)
{
- OGS_FATAL("Mesh element type %s is not supported", MeshLib::CellType2String(e->getCellType()).c_str());
+ Node* original_node =
+ const_cast<Node*>(quadratic_element->getNode(i));
+
+ auto it = unique_nodes.insert(original_node);
+ if (!it.second) // same node was already inserted before, no
+ // insertion
+ {
+ // Replace the element's node with the unique node.
+ quadratic_element->setNode(i, *it.first);
+ // And delete the original node
+ delete original_node;
+ }
}
+
+ quadratic_elements.push_back(quadratic_element.release());
}
- std::unique_ptr<MeshLib::Mesh> new_mesh(
- new MeshLib::Mesh(org_mesh.getName(), vec_new_nodes, vec_new_eles,
- org_mesh.getProperties().excludeCopyProperties(
- std::vector<MeshLib::MeshItemType>(
- 1, MeshLib::MeshItemType::Node)),
- org_mesh.getNumberOfNodes()));
- return new_mesh;
+ // Add the unique quadratic nodes to the cloned linear nodes.
+ quadratic_mesh_nodes.reserve(linear_mesh.getNodes().size() +
+ unique_nodes.size());
+ std::copy(unique_nodes.begin(), unique_nodes.end(),
+ std::back_inserter(quadratic_mesh_nodes));
+
+ return std::make_unique<MeshLib::Mesh>(
+ linear_mesh.getName(), quadratic_mesh_nodes, quadratic_elements,
+ linear_mesh.getProperties().excludeCopyProperties(
+ std::vector<MeshLib::MeshItemType>(1, MeshLib::MeshItemType::Node)),
+ linear_mesh.getNumberOfNodes());
}
-} // namespace MeshLib
-
+} // namespace MeshLib
diff --git a/Tests/MeshLib/TestQuadraticMesh.cpp b/Tests/MeshLib/TestQuadraticMesh.cpp
index 2beafb3c85d88a77003a465cf3551574f2af74cd..0e79acb6440b5ecb6c4861215708e8d56418707e 100644
--- a/Tests/MeshLib/TestQuadraticMesh.cpp
+++ b/Tests/MeshLib/TestQuadraticMesh.cpp
@@ -82,29 +82,31 @@ TEST(MeshLib, QuadraticOrderMesh_Quad)
ASSERT_FALSE(mesh->isBaseNode(e->getNodeIndex(i)));
}
- auto isIn = [](MeshLib::Node const* node, std::initializer_list<std::size_t> list)
- {
- return list.end() != std::find(list.begin(), list.end(), node->getID());
- };
-
- for (MeshLib::Node const* node : mesh->getNodes())
- {
- if (isIn(node, {4}))
- {
- ASSERT_EQ(4u, node->getElements().size());
- ASSERT_EQ(21u, node->getConnectedNodes().size());
- }
- else if (isIn(node, {0, 2, 6, 8, 9, 10, 11, 13, 15, 18, 19, 20}))
- {
- ASSERT_EQ(1u, node->getElements().size());
- ASSERT_EQ(8u, node->getConnectedNodes().size());
- }
- else
- {
- ASSERT_EQ(2u, node->getElements().size());
- ASSERT_EQ(13u, node->getConnectedNodes().size());
- }
- }
+ auto const& mesh_nodes = mesh->getNodes();
+
+ // Count nodes shared by four elements and also connected to all 21 other
+ // nodes.
+ ASSERT_EQ(1, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+ [](Node* const n) {
+ return (n->getElements().size() == 4) &&
+ (n->getConnectedNodes().size() == 21);
+ }));
+
+ // Count nodes belonging to one element and also connected to all 8 other
+ // nodes of that corner element.
+ ASSERT_EQ(12, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+ [](Node* const n) {
+ return (n->getElements().size() == 1) &&
+ (n->getConnectedNodes().size() == 8);
+ }));
+
+ // Count nodes shared by two elements and also connected to the 13 other
+ // nodes of the two elements.
+ ASSERT_EQ(8, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+ [](Node* const n) {
+ return (n->getElements().size() == 2) &&
+ (n->getConnectedNodes().size() == 13);
+ }));
}
TEST(MeshLib, QuadraticOrderMesh_LineQuad)
@@ -159,32 +161,37 @@ TEST(MeshLib, QuadraticOrderMesh_LineQuad)
ASSERT_FALSE(mesh->isBaseNode(e->getNodeIndex(i)));
}
- auto isIn = [](MeshLib::Node const* node, std::initializer_list<std::size_t> list)
- {
- return list.end() != std::find(list.begin(), list.end(), node->getID());
- };
-
- for (MeshLib::Node const* node : mesh->getNodes())
- {
- if (isIn(node, {4}))
- {
- ASSERT_EQ(6u, node->getElements().size());
- ASSERT_EQ(21u, node->getConnectedNodes().size());
- }
- else if (isIn(node, {0, 2, 6, 8, 9, 10, 11, 13, 15, 18, 19, 20}))
- {
- ASSERT_EQ(1u, node->getElements().size());
- ASSERT_EQ(8u, node->getConnectedNodes().size());
- }
- else if (isIn(node, {3, 5, 14, 16}))
- {
- ASSERT_EQ(3u, node->getElements().size());
- ASSERT_EQ(13u, node->getConnectedNodes().size());
- }
- else
- {
- ASSERT_EQ(2u, node->getElements().size());
- ASSERT_EQ(13u, node->getConnectedNodes().size());
- }
- }
+ auto const& mesh_nodes = mesh->getNodes();
+
+ // Count nodes shared by six elements and also connected to all 21 other
+ // nodes.
+ ASSERT_EQ(1, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+ [](Node* const n) {
+ return (n->getElements().size() == 6) &&
+ (n->getConnectedNodes().size() == 21);
+ }));
+
+ // Count nodes belonging to one element and also connected to all 8 other
+ // nodes of that corner element.
+ ASSERT_EQ(12, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+ [](Node* const n) {
+ return (n->getElements().size() == 1) &&
+ (n->getConnectedNodes().size() == 8);
+ }));
+
+ // Count nodes shared by three elements (quads and the line) and also
+ // connected to the 13 other nodes of the two elements.
+ ASSERT_EQ(4, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+ [](Node* const n) {
+ return (n->getElements().size() == 3) &&
+ (n->getConnectedNodes().size() == 13);
+ }));
+
+ // Count nodes shared by two elements (quads) and also connected to the 13
+ // other nodes of the two elements.
+ ASSERT_EQ(4, std::count_if(mesh_nodes.begin(), mesh_nodes.end(),
+ [](Node* const n) {
+ return (n->getElements().size() == 2) &&
+ (n->getConnectedNodes().size() == 13);
+ }));
}