Skip to content
Snippets Groups Projects
Commit ae36e829 authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

Merge branch 'MathLibMeshLibCodeImprovements' into 'master' 

Rewrite AngleSkewMetric and code improvements using Eigen

See merge request ogs/ogs!3327
parents b5f0dd26 9b54dd39
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MathLib/MathTools.cpp
+ 12
− 10
View file @ ae36e829
......@@ -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,
......
MathLib/MathTools.h
+ 1
− 1
View file @ ae36e829
......@@ -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,
......
MeshLib/MeshQuality/AngleSkewMetric.cpp
+ 103
− 137
View file @ ae36e829
......@@ -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;
}
}
}
......
MeshLib/MeshQuality/AngleSkewMetric.h
+ 0
− 14
View file @ ae36e829
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment