Revision 7c388cd32b25673329004897491f44189c5272ef authored by Lars Bilke on 19 July 2021, 09:13:09 UTC, committed by Lars Bilke on 20 July 2021, 18:41:53 UTC
1 parent ca35d9c
Is2DMeshOnRotatedVerticalPlane.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
*
* 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)->getCoords();
auto const x2 = element->getNode(1)->getCoords();
auto const x3 = element->getNode(2)->getCoords();
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)->getCoords();
auto const x2 = element->getNode(1)->getCoords();
auto const x3 = element->getNode(2)->getCoords();
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
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...