From 97c91c8bde35f8df4ffd5154f9e3d1684f343cf5 Mon Sep 17 00:00:00 2001
From: Karsten Rink <karsten.rink@ufz.de>
Date: Thu, 20 Feb 2014 10:13:12 +0100
Subject: [PATCH] moved the 4-points-on-a-plane-check into AnalyticalGeometry
---
GeoLib/AnalyticalGeometry.cpp | 14 ++++++++++++++
GeoLib/AnalyticalGeometry.h | 4 ++++
MeshLib/Elements/TemplateQuad.tpp | 12 +-----------
3 files changed, 19 insertions(+), 11 deletions(-)
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index fec6cbbdcfe..16f9877c11d 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -355,4 +355,18 @@ bool dividedByPlane(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib
return false;
}
+bool pointsOnAPlane(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib::Point& c, const GeoLib::Point& d)
+{
+ const GeoLib::Point AB(b[0]-a[0], b[1]-a[1], b[2]-a[2]);
+ const GeoLib::Point AC(c[0]-a[0], c[1]-a[1], c[2]-a[2]);
+ const GeoLib::Point AD(d[0]-a[0], d[1]-a[1], d[2]-a[2]);
+
+ double squared_scalar_triple = pow(GeoLib::scalarTriple(AC, AD, AB), 2);
+ double normalisation_factor = (AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]) *
+ (AC[0]*AC[0]+AC[1]*AC[1]+AC[2]*AC[2]) *
+ (AD[0]*AD[0]+AD[1]*AD[1]+AD[2]*AD[2]);
+
+ return (squared_scalar_triple/normalisation_factor < std::numeric_limits<double>::epsilon());
+}
+
} // end namespace GeoLib
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 9623083b768..7ff3f05b271 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -150,6 +150,10 @@ double scalarTriple(GeoLib::Point const& u, GeoLib::Point const& v, GeoLib::Poin
bool dividedByPlane(const GeoLib::Point& a, const GeoLib::Point& b,
const GeoLib::Point& c, const GeoLib::Point& d);
+ /// Checks if the four given points are located on a plane.
+ bool pointsOnAPlane(const GeoLib::Point& a, const GeoLib::Point& b,
+ const GeoLib::Point& c, const GeoLib::Point& d);
+
} // end namespace GeoLib
diff --git a/MeshLib/Elements/TemplateQuad.tpp b/MeshLib/Elements/TemplateQuad.tpp
index 8f8fe18782c..70cd0aa073b 100644
--- a/MeshLib/Elements/TemplateQuad.tpp
+++ b/MeshLib/Elements/TemplateQuad.tpp
@@ -125,17 +125,7 @@ bool TemplateQuad<NNODES,CELLQUADTYPE>::isValid() const
if (this->_area <= std::numeric_limits<double>::epsilon())
return false;
- const GeoLib::Point AB((*_nodes[1])[0]-(*_nodes[0])[0], (*_nodes[1])[1]-(*_nodes[0])[1], (*_nodes[1])[2]-(*_nodes[0])[2]);
- const GeoLib::Point AC((*_nodes[2])[0]-(*_nodes[0])[0], (*_nodes[2])[1]-(*_nodes[0])[1], (*_nodes[2])[2]-(*_nodes[0])[2]);
- const GeoLib::Point AD((*_nodes[3])[0]-(*_nodes[0])[0], (*_nodes[3])[1]-(*_nodes[0])[1], (*_nodes[3])[2]-(*_nodes[0])[2]);
-
- // check if all points lie on the same plane
- double squared_scalar_triple = pow(GeoLib::scalarTriple(AC, AD, AB), 2);
- double normalisation_factor = (AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]) *
- (AC[0]*AC[0]+AC[1]*AC[1]+AC[2]*AC[2]) *
- (AD[0]*AD[0]+AD[1]*AD[1]+AD[2]*AD[2]);
-
- if (squared_scalar_triple/normalisation_factor < std::numeric_limits<double>::epsilon())
+ if (GeoLib::pointsOnAPlane(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]))
{
// check if quad is convex
if (GeoLib::dividedByPlane(*_nodes[0], *_nodes[2], *_nodes[1], *_nodes[3]) &&
--
GitLab