diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 469c3ab0c033277b4801001d51c403e2bb701f45..5d7017408e0e27d0471fc812588c52c4f8b96509 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -20,12 +20,10 @@ namespace MathLib
double calcProjPntToLineAndDists(Point3d const& pp, Point3d const& pa,
Point3d const& pb, double& lambda, double& d0)
{
- auto const a =
- Eigen::Map<Eigen::Vector3d>(const_cast<double*>(pa.getCoords()));
- auto const b =
- Eigen::Map<Eigen::Vector3d>(const_cast<double*>(pb.getCoords()));
- auto const p =
- Eigen::Map<Eigen::Vector3d>(const_cast<double*>(pp.getCoords()));
+ auto const a = Eigen::Map<Eigen::Vector3d const>(pa.getCoords());
+ auto const b = Eigen::Map<Eigen::Vector3d const>(pb.getCoords());
+ auto const p = Eigen::Map<Eigen::Vector3d const>(pp.getCoords());
+
// g(lambda) = a + lambda v, v = b-a
Eigen::Vector3d const v = b - a;
@@ -42,13 +40,17 @@ double calcProjPntToLineAndDists(Point3d const& pp, Point3d const& pa,
return (p - proj_pnt).norm();
}
-double getAngle (const double p0[3], const double p1[3], const double p2[3])
+double getAngle(Point3d const& p0, Point3d const& p1, Point3d const& p2)
{
- const double v0[3] = {p0[0]-p1[0], p0[1]-p1[1], p0[2]-p1[2]};
- const double v1[3] = {p2[0]-p1[0], p2[1]-p1[1], p2[2]-p1[2]};
+ auto const a = Eigen::Map<Eigen::Vector3d const>(p0.getCoords());
+ auto const b = Eigen::Map<Eigen::Vector3d const>(p1.getCoords());
+ auto const c = Eigen::Map<Eigen::Vector3d const>(p2.getCoords());
+ Eigen::Vector3d const v0 = a - b;
+ Eigen::Vector3d const v1 = c - b;
// apply Cauchy Schwarz inequality
- return std::acos (scalarProduct<double,3> (v0,v1) / (std::sqrt(scalarProduct<double,3>(v0,v0)) * sqrt(scalarProduct<double,3>(v1,v1))));
+ return std::acos(
+ (v0.transpose() * v1 / (v0.norm() * v1.norm()))(0, 0));
}
double scalarTriple(Eigen::Vector3d const& u, Eigen::Vector3d const& v,
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 6018091a454ec1fa979adb89e82e32935504f0eb..699d3ff7225ed3ea8c2368d25c29b90ad54060d2 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -90,7 +90,7 @@ double calcProjPntToLineAndDists(MathLib::Point3d const& pp,
* @param p2 end point of edge 1
* @return the angle between the edges
*/
-double getAngle (const double p0[3], const double p1[3], const double p2[3]);
+double getAngle(Point3d const& p0, Point3d const& p1, Point3d const& p2);
/// Calculates the scalar triple (u x v) . w
double scalarTriple(Eigen::Vector3d const& u, Eigen::Vector3d const& v,
diff --git a/MeshLib/MeshQuality/AngleSkewMetric.cpp b/MeshLib/MeshQuality/AngleSkewMetric.cpp
index d5f927d99f9f21fbb1a6d6121ea96a49ea6f00da..ac8e9d75abced59108b51d17a0767999ff3f9b60 100644
--- a/MeshLib/MeshQuality/AngleSkewMetric.cpp
+++ b/MeshLib/MeshQuality/AngleSkewMetric.cpp
@@ -25,190 +25,156 @@ using namespace boost::math::double_constants;
namespace MeshLib
{
-AngleSkewMetric::AngleSkewMetric(Mesh const& mesh) :
- ElementQualityMetric(mesh)
-{}
-void AngleSkewMetric::calculateQuality ()
+namespace
{
- const std::vector<MeshLib::Element*>& elements(_mesh.getElements());
- const std::size_t nElements (_mesh.getNumberOfElements());
+template <unsigned long N>
+std::tuple<double, double> getMinMaxAngle(
+ std::array<MeshLib::Node, N> const& nodes)
+{
+ double min_angle(two_pi);
+ double max_angle(0.0);
- for (std::size_t k(0); k < nElements; k++)
+ for (decltype(N) i = 0; i < N; ++i)
{
- Element const& elem (*elements[k]);
- switch (elem.getGeomType())
- {
- case MeshElemType::LINE:
- _element_quality_metric[k] = -1.0;
- break;
- case MeshElemType::TRIANGLE:
- _element_quality_metric[k] = checkTriangle (elem);
- break;
- case MeshElemType::QUAD:
- _element_quality_metric[k] = checkQuad (elem);
- break;
- case MeshElemType::TETRAHEDRON:
- _element_quality_metric[k] = checkTetrahedron (elem);
- break;
- case MeshElemType::HEXAHEDRON:
- _element_quality_metric[k] = checkHexahedron (elem);
- break;
- case MeshElemType::PRISM:
- _element_quality_metric[k] = checkPrism (elem);
- break;
- default:
- break;
- }
+ const double angle(MathLib::getAngle(nodes[i], nodes[(i + 1) % N],
+ nodes[(i + 2) % N]));
+ min_angle = std::min(angle, min_angle);
+ max_angle = std::max(angle, max_angle);
}
+ return {min_angle, max_angle};
}
-double AngleSkewMetric::checkTriangle (Element const& elem) const
+double checkTriangle(Element const& elem)
{
- double const* const node0 (elem.getNode(0)->getCoords());
- double const* const node1 (elem.getNode(1)->getCoords());
- double const* const node2 (elem.getNode(2)->getCoords());
-
- double min_angle(two_pi);
- double max_angle(0.0);
- getMinMaxAngleFromTriangle (node0, node1, node2, min_angle, max_angle);
-
+ std::array const nodes = {*elem.getNode(0), *elem.getNode(1),
+ *elem.getNode(2)};
+ auto const& [min_angle, max_angle] = getMinMaxAngle(nodes);
return std::max((max_angle - third_pi) / two_thirds_pi,
(third_pi - min_angle) / third_pi);
}
-double AngleSkewMetric::checkQuad (Element const& elem) const
+double checkQuad(Element const& elem)
{
- double const* const node0 (elem.getNode(0)->getCoords());
- double const* const node1 (elem.getNode(1)->getCoords());
- double const* const node2 (elem.getNode(2)->getCoords());
- double const* const node3 (elem.getNode(3)->getCoords());
-
- double min_angle (two_pi);
- double max_angle (0.0);
-
- getMinMaxAngleFromQuad (node0, node1, node2, node3, min_angle, max_angle);
+ std::array const nodes = {*elem.getNode(0), *elem.getNode(1),
+ *elem.getNode(2), *elem.getNode(3)};
+ auto const& [min_angle, max_angle] = getMinMaxAngle(nodes);
return std::max((max_angle - half_pi) / half_pi,
(half_pi - min_angle) / half_pi);
}
-double AngleSkewMetric::checkTetrahedron (Element const& elem) const
+double checkTetrahedron(Element const& elem)
{
- double const* const node0 (elem.getNode(0)->getCoords());
- double const* const node1 (elem.getNode(1)->getCoords());
- double const* const node2 (elem.getNode(2)->getCoords());
- double const* const node3 (elem.getNode(3)->getCoords());
-
- double min_angle (two_pi);
- double max_angle (0.0);
-
- // first triangle (0,1,2)
- getMinMaxAngleFromTriangle(node0, node1, node2, min_angle, max_angle);
- // second triangle (0,1,3)
- getMinMaxAngleFromTriangle(node0, node1, node3, min_angle, max_angle);
- // third triangle (0,2,3)
- getMinMaxAngleFromTriangle(node0, node2, node3, min_angle, max_angle);
- // fourth triangle (1,2,3)
- getMinMaxAngleFromTriangle(node1, node2, node3, min_angle, max_angle);
+ std::array<double, 4> min;
+ std::array<double, 4> max;
+ for (auto face_number = 0; face_number < 4; ++face_number)
+ {
+ auto const& face = *elem.getFace(face_number);
+ std::array const nodes = {*face.getNode(0), *face.getNode(1),
+ *face.getNode(2)};
+ std::tie(min[face_number], max[face_number]) = getMinMaxAngle(nodes);
+ }
+
+ double const min_angle = *std::min_element(min.begin(), min.end());
+ double const max_angle = *std::max_element(max.begin(), max.end());
return std::max((max_angle - third_pi) / two_thirds_pi,
(third_pi - min_angle) / third_pi);
}
-double AngleSkewMetric::checkHexahedron (Element const& elem) const
+double checkHexahedron(Element const& elem)
{
- double const* const node0 (elem.getNode(0)->getCoords());
- double const* const node1 (elem.getNode(1)->getCoords());
- double const* const node2 (elem.getNode(2)->getCoords());
- double const* const node3 (elem.getNode(3)->getCoords());
- double const* const node4 (elem.getNode(4)->getCoords());
- double const* const node5 (elem.getNode(5)->getCoords());
- double const* const node6 (elem.getNode(6)->getCoords());
- double const* const node7 (elem.getNode(7)->getCoords());
-
- double min_angle (two_pi);
- double max_angle (0.0);
-
- // first surface (0,1,2,3)
- getMinMaxAngleFromQuad (node0, node1, node2, node3, min_angle, max_angle);
- // second surface (0,3,7,4)
- getMinMaxAngleFromQuad (node0, node3, node7, node4, min_angle, max_angle);
- // third surface (4,5,6,7)
- getMinMaxAngleFromQuad (node4, node5, node6, node7, min_angle, max_angle);
- // fourth surface (5,1,2,6)
- getMinMaxAngleFromQuad (node5, node1, node2, node6, min_angle, max_angle);
- // fifth surface (5,1,0,4)
- getMinMaxAngleFromQuad (node5, node1, node0, node4, min_angle, max_angle);
- // sixth surface (6,2,3,7)
- getMinMaxAngleFromQuad (node6, node2, node3, node7, min_angle, max_angle);
+ std::array<double, 6> min;
+ std::array<double, 6> max;
+ for (auto face_number = 0; face_number < 6; ++face_number)
+ {
+ auto const& face = *elem.getFace(face_number);
+ std::array const nodes = {*face.getNode(0), *face.getNode(1),
+ *face.getNode(2), *face.getNode(3)};
+ std::tie(min[face_number], max[face_number]) = getMinMaxAngle(nodes);
+ }
+
+ double const min_angle = *std::min_element(min.begin(), min.end());
+ double const max_angle = *std::max_element(max.begin(), max.end());
return std::max((max_angle - half_pi) / half_pi,
(half_pi - min_angle) / half_pi);
}
-double AngleSkewMetric::checkPrism (Element const& elem) const
+double checkPrism(Element const& elem)
{
- double const* const node0 (elem.getNode(0)->getCoords());
- double const* const node1 (elem.getNode(1)->getCoords());
- double const* const node2 (elem.getNode(2)->getCoords());
- double const* const node3 (elem.getNode(3)->getCoords());
- double const* const node4 (elem.getNode(4)->getCoords());
- double const* const node5 (elem.getNode(5)->getCoords());
+ // face 0: triangle (0,1,2)
+ auto const& f0 = *elem.getFace(0);
+ std::array const nodes_f0 = {*f0.getNode(0), *f0.getNode(1),
+ *f0.getNode(2)};
+ auto const& [min_angle_tri0, max_angle_tri0] = getMinMaxAngle(nodes_f0);
- double min_angle_tri (two_pi);
- double max_angle_tri (0.0);
+ // face 4: triangle (3,4,5)
+ auto const& f4 = *elem.getFace(4);
+ std::array const nodes_f4 = {*f4.getNode(0), *f4.getNode(1),
+ *f4.getNode(2)};
+ auto const& [min_angle_tri1, max_angle_tri1] = getMinMaxAngle(nodes_f4);
- // first triangle (0,1,2)
- getMinMaxAngleFromTriangle (node0, node1, node2, min_angle_tri, max_angle_tri);
- // second surface (3,4,5)
- getMinMaxAngleFromTriangle (node3, node4, node5, min_angle_tri, max_angle_tri);
+ auto const min_angle_tri = std::min(min_angle_tri0, min_angle_tri1);
+ auto const max_angle_tri = std::max(max_angle_tri0, max_angle_tri1);
double const tri_criterion(
std::max((max_angle_tri - third_pi) / two_thirds_pi,
(third_pi - min_angle_tri) / third_pi));
- double min_angle_quad (two_pi);
- double max_angle_quad (0.0);
- // surface (0,3,4,1)
- getMinMaxAngleFromQuad (node0, node3, node4, node1, min_angle_quad, max_angle_quad);
- // surface (2,5,3,0)
- getMinMaxAngleFromQuad (node2, node5, node3, node0, min_angle_quad, max_angle_quad);
- // surface (1,2,5,4)
- getMinMaxAngleFromQuad (node1, node2, node5, node4, min_angle_quad, max_angle_quad);
+ std::array<double, 3> min;
+ std::array<double, 3> max;
+ for (int i = 1; i < 4; ++i)
+ {
+ auto const& f = *elem.getFace(i);
+ std::array const nodes = {*f.getNode(0), *f.getNode(1), *f.getNode(2),
+ *f.getNode(3)};
+ std::tie(min[i - 1], max[i - 1]) = getMinMaxAngle(nodes);
+ }
+
+ double const min_angle_quad = *std::min_element(min.begin(), min.end());
+ double const max_angle_quad = *std::max_element(max.begin(), max.end());
double const quad_criterion(std::max((max_angle_quad - half_pi) / half_pi,
(half_pi - min_angle_quad) / half_pi));
- return std::min (tri_criterion, quad_criterion);
+ return std::min(tri_criterion, quad_criterion);
}
-void AngleSkewMetric::getMinMaxAngleFromQuad (
- double const* const n0, double const* const n1,
- double const* const n2, double const* const n3,
- double &min_angle, double &max_angle) const
-{
- const double* nodes[4] = {n0, n1, n2, n3};
- for (unsigned i=0; i<4; ++i)
- {
- const double angle (MathLib::getAngle (nodes[i], nodes[(i+1)%4], nodes[(i+2)%4]));
- min_angle = std::min(angle, min_angle);
- max_angle = std::max(angle, max_angle);
- }
-}
+} // end unnamed namespace
+
+AngleSkewMetric::AngleSkewMetric(Mesh const& mesh) :
+ ElementQualityMetric(mesh)
+{}
-void AngleSkewMetric::getMinMaxAngleFromTriangle(double const* const n0,
- double const* const n1,
- double const* const n2,
- double &min_angle,
- double &max_angle) const
+void AngleSkewMetric::calculateQuality()
{
- const double* nodes[3] = {n0, n1, n2};
- for (unsigned i=0; i<3; ++i)
+ for (auto const e : _mesh.getElements())
{
- const double angle (MathLib::getAngle (nodes[i], nodes[(i+1)%3], nodes[(i+2)%3]));
- min_angle = std::min(angle, min_angle);
- max_angle = std::max(angle, max_angle);
+ switch (e->getGeomType())
+ {
+ case MeshElemType::LINE:
+ _element_quality_metric[e->getID()] = -1.0;
+ break;
+ case MeshElemType::TRIANGLE:
+ _element_quality_metric[e->getID()] = checkTriangle(*e);
+ break;
+ case MeshElemType::QUAD:
+ _element_quality_metric[e->getID()] = checkQuad(*e);
+ break;
+ case MeshElemType::TETRAHEDRON:
+ _element_quality_metric[e->getID()] = checkTetrahedron(*e);
+ break;
+ case MeshElemType::HEXAHEDRON:
+ _element_quality_metric[e->getID()] = checkHexahedron(*e);
+ break;
+ case MeshElemType::PRISM:
+ _element_quality_metric[e->getID()] = checkPrism(*e);
+ break;
+ default:
+ break;
+ }
}
}
diff --git a/MeshLib/MeshQuality/AngleSkewMetric.h b/MeshLib/MeshQuality/AngleSkewMetric.h
index c2765b2cf25b2075feda80f9fd771c21f1c9bf72..b93fd9bcc91fa95bd4c8a47464d00f35c0bbe9f1 100644
--- a/MeshLib/MeshQuality/AngleSkewMetric.h
+++ b/MeshLib/MeshQuality/AngleSkewMetric.h
@@ -28,19 +28,5 @@ public:
explicit AngleSkewMetric(Mesh const& mesh);
void calculateQuality() override;
-
-private:
- double checkTriangle(Element const& elem) const;
- double checkQuad(Element const& elem) const;
- double checkTetrahedron(Element const& elem) const;
- double checkHexahedron(Element const& elem) const;
- double checkPrism (Element const& elem) const;
- void getMinMaxAngleFromQuad(double const* const n0,
- double const* const n1, double const* const n2,
- double const* const n3, double &min_angle,
- double &max_angle) const;
- void getMinMaxAngleFromTriangle(double const* const n0,
- double const* const n1, double const* const n2,
- double &min_angle, double &max_angle) const;
};
} // namespace MeshLib