diff --git a/src/MeshLoader.cpp b/src/MeshLoader.cpp index fc4c20011d8e85c6c97b6a7cac314e78e7099f04..ec0e39b3210debdddd0a22a2568e977578f8034d 100644 --- a/src/MeshLoader.cpp +++ b/src/MeshLoader.cpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include "MeshLoader.h" @@ -220,3 +222,70 @@ std::vector> MeshLoader::NodesNeighbors() { return ans; } + +double GetAngleWorker(Node& a, Node& b, Node& c, Node& d) { + double tmp1x = (c.coord[1] - a.coord[1]) * (b.coord[2] - a.coord[2]) - + (b.coord[1] - a.coord[1]) * (c.coord[2] - a.coord[2]); + + double tmp1y = (b.coord[0] - a.coord[0]) * (c.coord[2] - a.coord[2]) - + (c.coord[0] - a.coord[0]) * (b.coord[2] - a.coord[2]); + + double tmp1z = (c.coord[0] - a.coord[0]) * (b.coord[1] - a.coord[1]) - + (b.coord[0] - a.coord[0]) * (c.coord[1] - a.coord[1]); + + double tmp1 = tmp1x * tmp1x + tmp1y * tmp1y + tmp1z * tmp1z; + + double tmp2x = (d.coord[1] - a.coord[1]) * (b.coord[2] - a.coord[2]) - + (b.coord[1] - a.coord[1]) * (d.coord[2] - a.coord[2]); + + double tmp2y = (b.coord[0] - a.coord[0]) * (d.coord[2] - a.coord[2]) - + (d.coord[0] - a.coord[0]) * (b.coord[2] - a.coord[2]); + + double tmp2z = (d.coord[0] - a.coord[0]) * (b.coord[1] - a.coord[1]) - + (b.coord[0] - a.coord[0]) * (d.coord[1] - a.coord[1]); + + double tmp2 = tmp2x * tmp2x + tmp2y * tmp2y + tmp2z * tmp2z; + + double chsl = tmp1 + tmp2 - ((d.coord[0] - c.coord[0]) * (d.coord[0] - c.coord[0]) + + (d.coord[1] - c.coord[1]) * (d.coord[1] - c.coord[1]) + + (d.coord[2] - c.coord[2]) * (d.coord[2] - c.coord[2])) * + ((b.coord[0] - a.coord[0]) * (b.coord[0] - a.coord[0]) + + (b.coord[1] - a.coord[1]) * (b.coord[1] - a.coord[1]) + + (b.coord[2] - a.coord[2]) * (b.coord[2] - a.coord[2])); + + double znm = 2 * std::sqrt(tmp1) * std::sqrt(tmp2); + + return acos(chsl / znm); +} + +double MeshLoader::GetAngle(const FiniteElement& elem) { + if (elem.nodes_id.size() < 4) + return 0.0; + Node a = GetNodeById(elem.nodes_id[0]); + Node b = GetNodeById(elem.nodes_id[1]); + Node c = GetNodeById(elem.nodes_id[2]); + Node d = GetNodeById(elem.nodes_id[3]); + + double alpha = GetAngleWorker(a, b, c, d); + double betha = GetAngleWorker(b, c, a, d); + double hama = GetAngleWorker(a, c, b, d); + + return std::max(std::max(alpha, betha), hama); +} + +void MeshLoader::SortByAngle(const std::string& filepath) { + std::vector elements = FiniteElements(); + if (elements.empty()) + return; + std::ofstream file(filepath + ".txt"); + std::ranges::sort(elements, std::less<>{}, [&](const FiniteElement& elem) { + return GetAngle(elem); + }); + std::ranges::copy( + elements | std::ranges::views::drop(9) | std::ranges::views::take(10) | + std::ranges::views::transform( + [](const FiniteElement& elem) { + return elem.nodes_id; + }), + std::ostream_iterator>(file)); +}