17 #include <boost/math/constants/constants.hpp>
23 using namespace boost::math::double_constants;
29 template <
unsigned long N>
31 std::array<MeshLib::Node, N>
const& nodes)
33 double min_angle(two_pi);
34 double max_angle(0.0);
36 for (decltype(N) i = 0; i < N; ++i)
40 min_angle = std::min(angle, min_angle);
41 max_angle = std::max(angle, max_angle);
43 return {min_angle, max_angle};
51 return std::max((max_angle - third_pi) / two_thirds_pi,
52 (third_pi - min_angle) / third_pi);
61 return std::max((max_angle - half_pi) / half_pi,
62 (half_pi - min_angle) / half_pi);
67 std::array<double, 4> min;
68 std::array<double, 4> max;
69 for (
auto face_number = 0; face_number < 4; ++face_number)
71 std::unique_ptr<Element const> face{elem.
getFace(face_number)};
72 std::array
const nodes = {*face->
getNode(0), *face->getNode(1),
74 std::tie(min[face_number], max[face_number]) =
getMinMaxAngle(nodes);
77 double const min_angle = *std::min_element(min.begin(), min.end());
78 double const max_angle = *std::max_element(max.begin(), max.end());
80 return std::max((max_angle - third_pi) / two_thirds_pi,
81 (third_pi - min_angle) / third_pi);
86 std::array<double, 6> min;
87 std::array<double, 6> max;
88 for (
auto face_number = 0; face_number < 6; ++face_number)
90 std::unique_ptr<Element const> face{elem.
getFace(face_number)};
91 std::array
const nodes = {*face->
getNode(0), *face->getNode(1),
92 *face->getNode(2), *face->getNode(3)};
93 std::tie(min[face_number], max[face_number]) =
getMinMaxAngle(nodes);
96 double const min_angle = *std::min_element(min.begin(), min.end());
97 double const max_angle = *std::max_element(max.begin(), max.end());
99 return std::max((max_angle - half_pi) / half_pi,
100 (half_pi - min_angle) / half_pi);
106 std::unique_ptr<Element const> f0{elem.
getFace(0)};
107 std::array
const nodes_f0 = {*f0->
getNode(0), *f0->getNode(1),
109 auto const& [min_angle_tri0, max_angle_tri0] =
getMinMaxAngle(nodes_f0);
112 std::unique_ptr<Element const> f4{elem.
getFace(4)};
113 std::array
const nodes_f4 = {*f4->
getNode(0), *f4->getNode(1),
115 auto const& [min_angle_tri1, max_angle_tri1] =
getMinMaxAngle(nodes_f4);
117 auto const min_angle_tri = std::min(min_angle_tri0, min_angle_tri1);
118 auto const max_angle_tri = std::max(max_angle_tri0, max_angle_tri1);
120 double const tri_criterion(
121 std::max((max_angle_tri - third_pi) / two_thirds_pi,
122 (third_pi - min_angle_tri) / third_pi));
124 std::array<double, 3> min;
125 std::array<double, 3> max;
126 for (
int i = 1; i < 4; ++i)
128 std::unique_ptr<Element const> f{elem.
getFace(i)};
129 std::array
const nodes = {*f->
getNode(0), *f->getNode(1),
130 *f->getNode(2), *f->getNode(3)};
134 double const min_angle_quad = *std::min_element(min.begin(), min.end());
135 double const max_angle_quad = *std::max_element(max.begin(), max.end());
137 double const quad_criterion(std::max((max_angle_quad - half_pi) / half_pi,
138 (half_pi - min_angle_quad) / half_pi));
140 return std::min(tri_criterion, quad_criterion);
145 void AngleSkewMetric::calculateQuality()
147 for (
auto const e : _mesh.getElements())
149 switch (e->getGeomType())
151 case MeshElemType::LINE:
152 _element_quality_metric[e->getID()] = -1.0;
154 case MeshElemType::TRIANGLE:
157 case MeshElemType::QUAD:
158 _element_quality_metric[e->getID()] =
checkQuad(*e);
160 case MeshElemType::TETRAHEDRON:
163 case MeshElemType::HEXAHEDRON:
166 case MeshElemType::PRISM:
167 _element_quality_metric[e->getID()] =
checkPrism(*e);
Definition of the AngleSkewMetric class.
Definition of the Node class.
virtual const Node * getNode(unsigned idx) const =0
virtual const Element * getFace(unsigned i) const =0
Returns the i-th face of the element.
double getAngle(Point3d const &p0, Point3d const &p1, Point3d const &p2)
double checkQuad(Element const &elem)
double checkHexahedron(Element const &elem)
double checkTetrahedron(Element const &elem)
double checkPrism(Element const &elem)
std::tuple< double, double > getMinMaxAngle(std::array< MeshLib::Node, N > const &nodes)
double checkTriangle(Element const &elem)