TestMeshGenerator.cpp
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#include <gtest/gtest.h>
#include <cmath>
#include <memory>
#include <numeric>
#include "MathLib/MathTools.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Elements/Hex.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshGenerators/MeshGenerator.h"
#include "MeshLib/Node.h"
using namespace MeshLib;
TEST(MeshLib, MeshGeneratorRegularHex)
{
const double L = 10.0;
const std::size_t n_subdivisions = 9;
const double dL = L / static_cast<double>(n_subdivisions);
std::unique_ptr<Mesh> msh(
MeshGenerator::generateRegularHexMesh(L, n_subdivisions));
ASSERT_EQ(std::pow(n_subdivisions, 3), msh->getNumberOfElements());
ASSERT_EQ(std::pow(n_subdivisions + 1, 3), msh->getNumberOfNodes());
// check nodes
const Node& node0 = *msh->getNode(0);
const Node& node1 = *msh->getNode(1);
const Node& node_n = *msh->getNode(msh->getNumberOfNodes() - 2);
const Node& node_n1 = *msh->getNode(msh->getNumberOfNodes() - 1);
ASSERT_DOUBLE_EQ(.0, node0[0]);
ASSERT_DOUBLE_EQ(.0, node0[1]);
ASSERT_DOUBLE_EQ(.0, node0[2]);
ASSERT_DOUBLE_EQ(dL, node1[0]);
ASSERT_DOUBLE_EQ(.0, node1[1]);
ASSERT_DOUBLE_EQ(.0, node1[2]);
ASSERT_DOUBLE_EQ(L - dL, node_n[0]);
ASSERT_DOUBLE_EQ(L, node_n[1]);
ASSERT_DOUBLE_EQ(L, node_n[2]);
ASSERT_DOUBLE_EQ(L, node_n1[0]);
ASSERT_DOUBLE_EQ(L, node_n1[1]);
ASSERT_DOUBLE_EQ(L, node_n1[2]);
// check elements
const Element& ele0 = *msh->getElement(0);
const Element& ele_n = *msh->getElement(msh->getNumberOfElements() - 1);
const std::size_t offset_y0 = (n_subdivisions + 1);
const std::size_t offset_z0 = (n_subdivisions + 1) * (n_subdivisions + 1);
ASSERT_EQ(0u, getNodeIndex(ele0, 0));
ASSERT_EQ(1u, getNodeIndex(ele0, 1));
ASSERT_EQ(offset_y0 + 1, getNodeIndex(ele0, 2));
ASSERT_EQ(offset_y0, getNodeIndex(ele0, 3));
ASSERT_EQ(offset_z0, getNodeIndex(ele0, 4));
ASSERT_EQ(offset_z0 + 1, getNodeIndex(ele0, 5));
ASSERT_EQ(offset_z0 + offset_y0 + 1, getNodeIndex(ele0, 6));
ASSERT_EQ(offset_z0 + offset_y0, getNodeIndex(ele0, 7));
const std::size_t offset_yn0 = (n_subdivisions + 1) * (n_subdivisions - 1);
const std::size_t offset_yn1 = (n_subdivisions + 1) * n_subdivisions;
const std::size_t offset_zn0 =
(n_subdivisions + 1) * (n_subdivisions + 1) * (n_subdivisions - 1);
const std::size_t offset_zn1 =
(n_subdivisions + 1) * (n_subdivisions + 1) * n_subdivisions;
ASSERT_EQ(offset_zn0 + offset_yn0 + n_subdivisions - 1,
getNodeIndex(ele_n, 0));
ASSERT_EQ(offset_zn0 + offset_yn0 + n_subdivisions, getNodeIndex(ele_n, 1));
ASSERT_EQ(offset_zn0 + offset_yn1 + n_subdivisions, getNodeIndex(ele_n, 2));
ASSERT_EQ(offset_zn0 + offset_yn1 + n_subdivisions - 1,
getNodeIndex(ele_n, 3));
ASSERT_EQ(offset_zn1 + offset_yn0 + n_subdivisions - 1,
getNodeIndex(ele_n, 4));
ASSERT_EQ(offset_zn1 + offset_yn0 + n_subdivisions, getNodeIndex(ele_n, 5));
ASSERT_EQ(offset_zn1 + offset_yn1 + n_subdivisions, getNodeIndex(ele_n, 6));
ASSERT_EQ(offset_zn1 + offset_yn1 + n_subdivisions - 1,
getNodeIndex(ele_n, 7));
std::unique_ptr<Mesh> msh2(MeshGenerator::generateRegularHexMesh(
n_subdivisions, n_subdivisions, n_subdivisions, L / n_subdivisions));
ASSERT_EQ(msh->getNumberOfNodes(), msh2->getNumberOfNodes());
ASSERT_DOUBLE_EQ(
0, MathLib::sqrDist(*(msh->getNode(msh->getNumberOfNodes() - 1)),
*(msh2->getNode(msh->getNumberOfNodes() - 1))));
unsigned n_x(10);
unsigned n_y(5);
unsigned n_z(2);
double delta(1.2);
std::unique_ptr<Mesh> hex_mesh(
MeshGenerator::generateRegularHexMesh(n_x, n_y, n_z, delta));
ASSERT_EQ(n_x * n_y * n_z, hex_mesh->getNumberOfElements());
ASSERT_EQ((n_x + 1) * (n_y + 1) * (n_z + 1), hex_mesh->getNumberOfNodes());
const MeshLib::Node* node(
hex_mesh->getNode(hex_mesh->getNumberOfNodes() - 1));
ASSERT_DOUBLE_EQ(n_x * delta, (*node)[0]);
ASSERT_DOUBLE_EQ(n_y * delta, (*node)[1]);
ASSERT_DOUBLE_EQ(n_z * delta, (*node)[2]);
}
TEST(MeshLib, MeshGeneratorRegularQuad)
{
unsigned n_x(10);
unsigned n_y(5);
double delta(1.2);
std::unique_ptr<Mesh> quad_mesh(
MeshGenerator::generateRegularQuadMesh(n_x, n_y, delta));
ASSERT_EQ(n_x * n_y, quad_mesh->getNumberOfElements());
ASSERT_EQ((n_x + 1) * (n_y + 1), quad_mesh->getNumberOfNodes());
const MeshLib::Node* node(
quad_mesh->getNode(quad_mesh->getNumberOfNodes() - 1));
ASSERT_DOUBLE_EQ(n_x * delta, (*node)[0]);
ASSERT_DOUBLE_EQ(n_y * delta, (*node)[1]);
ASSERT_DOUBLE_EQ(0, (*node)[2]);
const double L = 10.0;
const std::size_t n_subdivisions = 9;
std::unique_ptr<Mesh> quad_mesh2(
MeshGenerator::generateRegularQuadMesh(L, n_subdivisions));
ASSERT_EQ(n_subdivisions * n_subdivisions,
quad_mesh2->getNumberOfElements());
node = quad_mesh2->getNode(quad_mesh2->getNumberOfNodes() - 1);
ASSERT_DOUBLE_EQ(L, (*node)[0]);
ASSERT_DOUBLE_EQ(L, (*node)[1]);
}
TEST(MeshLib, MeshGeneratorRegularPrism)
{
double const l_x(10);
double const l_y(10);
double const l_z(8);
unsigned n_x(10);
unsigned n_y(5);
unsigned n_z(4);
std::unique_ptr<Mesh> mesh(
MeshGenerator::generateRegularPrismMesh(l_x, l_y, l_z, n_x, n_y, n_z));
ASSERT_EQ(2 * n_x * n_y * n_z, mesh->getNumberOfElements());
ASSERT_EQ((n_x + 1) * (n_y + 1) * (n_z + 1), mesh->getNumberOfNodes());
std::size_t count(0);
for (std::size_t k = 0; k < (n_z + 1); ++k)
{
for (std::size_t j = 0; j < (n_y + 1); ++j)
{
for (std::size_t i = 0; i < (n_x + 1); ++i)
{
const MeshLib::Node* node(mesh->getNode(count++));
ASSERT_DOUBLE_EQ(static_cast<double>(i), (*node)[0]);
ASSERT_DOUBLE_EQ(static_cast<double>(j * 2), (*node)[1]);
ASSERT_DOUBLE_EQ(static_cast<double>(k * 2), (*node)[2]);
}
}
}
}
TEST(MeshLib, MeshGeneratorRegularPyramid)
{
const double L = 10.0;
const std::size_t n_subdivisions = 9;
const double dL = L / static_cast<double>(n_subdivisions);
std::vector<std::unique_ptr<BaseLib::ISubdivision>> vec_div;
vec_div.reserve(3);
for (unsigned i = 0; i < 3; i++)
{
vec_div.emplace_back(
new BaseLib::UniformSubdivision(L, n_subdivisions));
}
std::unique_ptr<Mesh> msh(MeshGenerator::generateRegularPyramidMesh(
*vec_div[0], *vec_div[1], *vec_div[2]));
// check number of generates nodes and elements
ASSERT_EQ(6 * std::pow(n_subdivisions, 3), msh->getNumberOfElements());
ASSERT_EQ(std::pow(n_subdivisions + 1, 3) + std::pow(n_subdivisions, 3),
msh->getNumberOfNodes());
// check the positions of the diagonal edge nodes (lower left front and
// upper right back)
const Node& node0 = *msh->getNode(0);
const Node& node1 = *msh->getNode(1);
const std::size_t L_index = static_cast<std::size_t>(L);
const Node& node_n = *msh->getNode(L_index * L_index * L_index - 2);
const Node& node_n1 = *msh->getNode(L_index * L_index * L_index - 1);
ASSERT_DOUBLE_EQ(.0, node0[0]);
ASSERT_DOUBLE_EQ(.0, node0[1]);
ASSERT_DOUBLE_EQ(.0, node0[2]);
ASSERT_DOUBLE_EQ(dL, node1[0]);
ASSERT_DOUBLE_EQ(.0, node1[1]);
ASSERT_DOUBLE_EQ(.0, node1[2]);
ASSERT_DOUBLE_EQ(L - dL, node_n[0]);
ASSERT_DOUBLE_EQ(L, node_n[1]);
ASSERT_DOUBLE_EQ(L, node_n[2]);
ASSERT_DOUBLE_EQ(L, node_n1[0]);
ASSERT_DOUBLE_EQ(L, node_n1[1]);
ASSERT_DOUBLE_EQ(L, node_n1[2]);
// check if the domain volume equals the volume of the elements
long double const element_volumes = std::accumulate(
cbegin(msh->getElements()), cend(msh->getElements()), 0.,
[](double const volume, auto* const element)
{ return volume + element->computeVolume(); });
EXPECT_NEAR(L * L * L, element_volumes, 1e-10);
// test node order of the elements
for (auto const element : msh->getElements())
{
ASSERT_TRUE(element->testElementNodeOrder());
}
// validate the elements
for (auto const element : msh->getElements())
{
ASSERT_TRUE(element->validate().none());
}
}
