/**
* \file
* \brief Tests for GeoLib::Surface::isPntInSfc()
*
* \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 <array>
#include <boost/math/constants/constants.hpp>
#include <memory>
#include <random>
#include <vector>
#include "GeoLib/AnalyticalGeometry.h"
#include "GeoLib/GEOObjects.h"
#include "GeoLib/Point.h"
#include "GeoLib/Surface.h"
#include "GeoLib/Triangle.h"
#include "MathLib/Point3d.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshGenerators/MeshGenerator.h"
#include "MeshLib/Node.h"
#include "MeshLib/convertMeshToGeo.h"
inline double constant(double /*unused*/, double /*unused*/)
{
return 0.0;
}
inline double coscos(double x, double y)
{
return std::cos(x) * std::cos(y);
}
inline MathLib::Point3d getEdgeMiddlePoint(MathLib::Point3d const& a,
MathLib::Point3d const& b)
{
return MathLib::Point3d(std::array<double, 3>(
{{(a[0] + b[0]) / 2, (a[1] + b[1]) / 2, (a[2] + b[2]) / 2}}));
}
inline std::tuple<MathLib::Point3d, MathLib::Point3d, MathLib::Point3d>
getEdgeMiddlePoints(GeoLib::Triangle const& tri)
{
return std::make_tuple(
getEdgeMiddlePoint(*tri.getPoint(0), *tri.getPoint(1)),
getEdgeMiddlePoint(*tri.getPoint(1), *tri.getPoint(2)),
getEdgeMiddlePoint(*tri.getPoint(2), *tri.getPoint(0)));
}
/// Computes rotation matrix according to the z,x',z'' convention
inline Eigen::Matrix3d getRotMat(double alpha, double beta, double gamma)
{
Eigen::Matrix3d rot_mat;
rot_mat(0, 0) =
cos(alpha) * cos(gamma) - sin(alpha) * cos(beta) * sin(gamma);
rot_mat(0, 1) =
sin(alpha) * cos(gamma) + cos(alpha) * cos(beta) * sin(gamma);
rot_mat(0, 2) = sin(beta) * sin(gamma);
rot_mat(1, 0) =
-cos(alpha) * sin(gamma) - sin(alpha) * cos(beta) * cos(gamma);
rot_mat(1, 1) =
-sin(alpha) * sin(gamma) + cos(alpha) * cos(beta) * cos(gamma);
rot_mat(1, 2) = sin(beta) * cos(gamma);
rot_mat(2, 0) = sin(alpha) * sin(beta);
rot_mat(2, 1) = -cos(alpha) * sin(beta);
rot_mat(2, 2) = cos(beta);
return rot_mat;
}
TEST(GeoLib, SurfaceIsPointInSurface)
{
std::vector<std::function<double(double, double)>> surface_functions;
surface_functions.emplace_back(constant);
surface_functions.emplace_back(coscos);
for (const auto& f : surface_functions)
{
std::random_device rd;
std::string name("Surface");
// generate ll and ur in random way
std::mt19937 random_engine_mt19937(rd());
std::normal_distribution<> normal_dist_ll(-10, 2);
std::normal_distribution<> normal_dist_ur(10, 2);
MathLib::Point3d ll(std::array<double, 3>(
{{normal_dist_ll(random_engine_mt19937),
normal_dist_ll(random_engine_mt19937), 0.0}}));
MathLib::Point3d ur(std::array<double, 3>(
{{normal_dist_ur(random_engine_mt19937),
normal_dist_ur(random_engine_mt19937), 0.0}}));
for (std::size_t k(0); k < 3; ++k)
{
if (ll[k] > ur[k])
{
std::swap(ll[k], ur[k]);
}
}
// random discretization of the domain
std::default_random_engine re(rd());
std::uniform_int_distribution<std::size_t> uniform_dist(2, 25);
std::array<std::size_t, 2> n_steps = {
{uniform_dist(re), uniform_dist(re)}};
std::unique_ptr<MeshLib::Mesh> sfc_mesh(
MeshLib::MeshGenerator::createSurfaceMesh(name, ll, ur, n_steps,
f));
// random rotation angles
std::normal_distribution<> normal_dist_angles(
0, boost::math::double_constants::two_pi);
std::array<double, 3> euler_angles = {
{normal_dist_angles(random_engine_mt19937),
normal_dist_angles(random_engine_mt19937),
normal_dist_angles(random_engine_mt19937)}};
Eigen::Matrix3d rot_mat(
getRotMat(euler_angles[0], euler_angles[1], euler_angles[2]));
std::vector<MeshLib::Node*> const& nodes(sfc_mesh->getNodes());
GeoLib::rotatePoints<MeshLib::Node>(rot_mat, nodes);
Eigen::Vector3d const normal({0, 0, 1.0});
Eigen::Vector3d const surface_normal(rot_mat * normal);
double const scaling(1e-6);
Eigen::Vector3d const displacement(scaling * surface_normal);
GeoLib::GEOObjects geometries;
MeshLib::convertMeshToGeo(*sfc_mesh, geometries);
std::vector<GeoLib::Surface*> const& sfcs(
*geometries.getSurfaceVec(name));
GeoLib::Surface const* const sfc(sfcs.front());
std::vector<GeoLib::Point*> const& pnts(*geometries.getPointVec(name));
double const eps(std::numeric_limits<double>::epsilon());
// test triangle edge point of the surface triangles
for (auto const p : pnts)
{
EXPECT_TRUE(sfc->isPntInSfc(*p, eps));
MathLib::Point3d q(*p);
for (std::size_t k(0); k < 3; ++k)
{
q[k] += displacement[k];
}
EXPECT_FALSE(sfc->isPntInSfc(q, eps));
}
// test edge middle points of the triangles
for (std::size_t k(0); k < sfc->getNumberOfTriangles(); ++k)
{
MathLib::Point3d p;
MathLib::Point3d q;
MathLib::Point3d r;
std::tie(p, q, r) = getEdgeMiddlePoints(*(*sfc)[k]);
EXPECT_TRUE(sfc->isPntInSfc(p, eps));
EXPECT_TRUE(sfc->isPntInSfc(q, eps));
EXPECT_TRUE(sfc->isPntInSfc(r, eps));
}
}
}