diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index ad2c06fa8ecb3977d4c957c851b7a4b440b23e0b..0e8375c15cd6d1a03b4c64fc313dca526367667c 100644 --- a/GeoLib/AnalyticalGeometry.cpp +++ b/GeoLib/AnalyticalGeometry.cpp @@ -359,17 +359,6 @@ double calcTriangleArea(MathLib::Point3d const& a, return 0.5 * w.getLength(); } -double calcTetrahedronVolume(MathLib::Point3d const& x1, - MathLib::Point3d const& x2, - MathLib::Point3d const& x3, - MathLib::Point3d const& x4) -{ - const MathLib::Vector3 ab(x1, x2); - const MathLib::Vector3 ac(x1, x3); - const MathLib::Vector3 ad(x1, x4); - return std::abs(MathLib::scalarTriple(ac, ad, ab)) / 6.0; -} - void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat) { // *** some frequently used terms *** diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h index 0c75ef137598cfd32f451686c8ee2122bfaec848..b7a495ff3c14740dc9c3e2d9d93f27c6766b7349 100644 --- a/GeoLib/AnalyticalGeometry.h +++ b/GeoLib/AnalyticalGeometry.h @@ -184,15 +184,6 @@ void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts); double calcTriangleArea(MathLib::Point3d const& a, MathLib::Point3d const& b, MathLib::Point3d const& c); -/** - * Calculates the volume of a tetrahedron. - * The formula is V=1/6*|a(b x c)| with a=x1->x2, b=x1->x3 and c=x1->x4. - */ -double calcTetrahedronVolume(MathLib::Point3d const& x1, - MathLib::Point3d const& x2, - MathLib::Point3d const& x3, - MathLib::Point3d const& x4); - /** * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c. * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c3fcff04588712aba2f73eb085f958e17f43e6f0 --- /dev/null +++ b/MathLib/GeometricBasics.cpp @@ -0,0 +1,28 @@ +/** + * + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "GeometricBasics.h" +#include "Point3d.h" +#include "Vector3.h" + +namespace MathLib +{ +double calcTetrahedronVolume(MathLib::Point3d const& a, + MathLib::Point3d const& b, + MathLib::Point3d const& c, + MathLib::Point3d const& d) +{ + const MathLib::Vector3 ab(a, b); + const MathLib::Vector3 ac(a, c); + const MathLib::Vector3 ad(a, d); + return std::abs(MathLib::scalarTriple(ac, ad, ab)) / 6.0; +} + +} // end namespace MathLib diff --git a/MathLib/GeometricBasics.h b/MathLib/GeometricBasics.h new file mode 100644 index 0000000000000000000000000000000000000000..c33fd7b6f596101f6b6444f1a0a431d994e2f4ae --- /dev/null +++ b/MathLib/GeometricBasics.h @@ -0,0 +1,33 @@ +/** + * + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef GEOMETRIC_BASICS_H_ +#define GEOMETRIC_BASICS_H_ + +#include <cstddef> + +namespace MathLib +{ + +template <typename T, std::size_t DIM> class TemplatePoint; +typedef MathLib::TemplatePoint<double,3> Point3d; + +/** + * Calculates the volume of a tetrahedron. + * The formula is V=1/6*|a(b x c)| with a=x1->x2, b=x1->x3 and c=x1->x4. + */ +double calcTetrahedronVolume(MathLib::Point3d const& x1, + MathLib::Point3d const& x2, + MathLib::Point3d const& x3, + MathLib::Point3d const& x4); + +} // end namespace MathLib + +#endif diff --git a/MeshLib/Elements/HexRule8.cpp b/MeshLib/Elements/HexRule8.cpp index 8847cfc43d4f6621ba7f3fa6414eab50a09ab10f..db669c01b53c10ac0ea4482d9d063b97db6e9a7f 100644 --- a/MeshLib/Elements/HexRule8.cpp +++ b/MeshLib/Elements/HexRule8.cpp @@ -14,6 +14,7 @@ #include "logog/include/logog.hpp" #include "GeoLib/AnalyticalGeometry.h" +#include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" #include "Quad.h" @@ -62,12 +63,12 @@ const Element* HexRule8::getFace(const Element* e, unsigned i) double HexRule8::computeVolume(Node const* const* _nodes) { - return GeoLib::calcTetrahedronVolume(*_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0]) - + GeoLib::calcTetrahedronVolume(*_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0]) - + GeoLib::calcTetrahedronVolume(*_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0]) - + GeoLib::calcTetrahedronVolume(*_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2]) - + GeoLib::calcTetrahedronVolume(*_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2]) - + GeoLib::calcTetrahedronVolume(*_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2]); + return MathLib::calcTetrahedronVolume(*_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0]) + + MathLib::calcTetrahedronVolume(*_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0]) + + MathLib::calcTetrahedronVolume(*_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0]) + + MathLib::calcTetrahedronVolume(*_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2]) + + MathLib::calcTetrahedronVolume(*_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2]) + + MathLib::calcTetrahedronVolume(*_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2]); } bool HexRule8::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps) diff --git a/MeshLib/Elements/PrismRule6.cpp b/MeshLib/Elements/PrismRule6.cpp index 8a7643c0afd359fecd29ba3d8f813525aabde782..8a1f72546f5d29e870035c50e69bb5af8802c8d2 100644 --- a/MeshLib/Elements/PrismRule6.cpp +++ b/MeshLib/Elements/PrismRule6.cpp @@ -12,6 +12,7 @@ #include "logog/include/logog.hpp" #include "GeoLib/AnalyticalGeometry.h" +#include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" #include "Quad.h" @@ -71,9 +72,9 @@ const Element* PrismRule6::getFace(const Element* e, unsigned i) double PrismRule6::computeVolume(Node const* const* _nodes) { - return GeoLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]) - + GeoLib::calcTetrahedronVolume(*_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3]) - + GeoLib::calcTetrahedronVolume(*_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3]); + return MathLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]) + + MathLib::calcTetrahedronVolume(*_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3]) + + MathLib::calcTetrahedronVolume(*_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3]); } bool PrismRule6::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps) diff --git a/MeshLib/Elements/PyramidRule5.cpp b/MeshLib/Elements/PyramidRule5.cpp index 1f50cac8fc2e889146f3eac55a694d1ce0726f3c..98b3edcfdace8225dc088e41996b199976dcac29 100644 --- a/MeshLib/Elements/PyramidRule5.cpp +++ b/MeshLib/Elements/PyramidRule5.cpp @@ -12,6 +12,7 @@ #include "logog/include/logog.hpp" #include "GeoLib/AnalyticalGeometry.h" +#include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" #include "Quad.h" @@ -70,8 +71,8 @@ const Element* PyramidRule5::getFace(const Element* e, unsigned i) double PyramidRule5::computeVolume(Node const* const* _nodes) { - return GeoLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4]) - + GeoLib::calcTetrahedronVolume(*_nodes[2], *_nodes[3], *_nodes[0], *_nodes[4]); + return MathLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4]) + + MathLib::calcTetrahedronVolume(*_nodes[2], *_nodes[3], *_nodes[0], *_nodes[4]); } bool PyramidRule5::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps) diff --git a/MeshLib/Elements/TetRule4.cpp b/MeshLib/Elements/TetRule4.cpp index 1d9ca6e4546abab2fadfc7881d8f28992dbe9dc8..daeee3683381dd24d33567be34a4b91cb56f7bba 100644 --- a/MeshLib/Elements/TetRule4.cpp +++ b/MeshLib/Elements/TetRule4.cpp @@ -14,6 +14,7 @@ #include "logog/include/logog.hpp" #include "GeoLib/AnalyticalGeometry.h" +#include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" #include "Tri.h" @@ -54,7 +55,7 @@ const Element* TetRule4::getFace(const Element* e, unsigned i) double TetRule4::computeVolume(Node const* const* _nodes) { - return GeoLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]); + return MathLib::calcTetrahedronVolume(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]); } bool TetRule4::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)