/**
* \file
* \copyright
* Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
* Created on April 22, 2021, 11:52 AM
*/
#include "Is2DMeshOnRotatedVerticalPlane.h"
#include <algorithm>
#include <limits>
#include "BaseLib/Error.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/Node.h"
namespace MeshLib
{
bool is2DMeshOnRotatedVerticalPlane(Mesh const& mesh)
{
auto const mesh_dimension = mesh.getDimension();
if (mesh_dimension != 2)
{
OGS_FATAL(
"A 2D mesh is required for this computation but the provided mesh, "
"mesh {:s}, has {:d}D elements.",
mesh.getName(), mesh_dimension);
}
auto const& elements = mesh.getElements();
bool const has_inclined_element =
std::any_of(elements.begin(), elements.end(),
[&mesh_dimension](auto const& element)
{ return element->space_dimension_ != mesh_dimension; });
if (!has_inclined_element)
{
return false;
}
const bool is_rotated_around_y_axis = std::all_of(
elements.cbegin(), elements.cend(),
[](auto const& element)
{
// 3 nodes are enough to make up a plane.
auto const x1 = element->getNode(0)->data();
auto const x2 = element->getNode(1)->data();
auto const x3 = element->getNode(2)->data();
double const a0 = x2[0] - x1[0];
double const a2 = x2[2] - x1[2];
double const b0 = x3[0] - x1[0];
double const b2 = x3[2] - x1[2];
double const e_n_1 = -a0 * b2 + a2 * b0;
return std::fabs(e_n_1) < std::numeric_limits<double>::epsilon();
});
const bool is_rotated_around_z_axis = std::all_of(
elements.cbegin(), elements.cend(),
[](auto const& element)
{
// 3 nodes are enough to make up a plane.
auto const x1 = element->getNode(0)->data();
auto const x2 = element->getNode(1)->data();
auto const x3 = element->getNode(2)->data();
double const a0 = x2[0] - x1[0];
double const a1 = x2[1] - x1[1];
double const b0 = x3[0] - x1[0];
double const b1 = x3[1] - x1[1];
double const e_n_2 = a0 * b1 - a1 * b0;
return std::fabs(e_n_2) < std::numeric_limits<double>::epsilon();
});
if (!(is_rotated_around_y_axis || is_rotated_around_z_axis))
{
OGS_FATAL(
"2D Mesh {:s} is on an inclined plane, which is neither a vertical "
"nor horizontal plane that is required for the present "
"computation.",
mesh.getName());
}
return true;
}
}; // namespace MeshLib