diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index 0b8cd960aa89d4d24477ffb6d432251f132c0653..fec6cbbdcfef3d399004354642637a2f2572671d 100644 --- a/GeoLib/AnalyticalGeometry.cpp +++ b/GeoLib/AnalyticalGeometry.cpp @@ -339,4 +339,20 @@ double scalarTriple(GeoLib::Point const& u, GeoLib::Point const& v, GeoLib::Poin return result; } +bool dividedByPlane(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib::Point& c, const GeoLib::Point& d) +{ + double abc(0), abd(0); + for (unsigned x=0; x<3; ++x) + { + const unsigned y=(x+1)%3; + abc = (b[x] - a[x])*(c[y] - a[y]) - (b[y] - a[y])*(c[x] - a[x]); + abd = (b[x] - a[x])*(d[y] - a[y]) - (b[y] - a[y])*(d[x] - a[x]); + + if ((abc>0 && abd<0) || (abc<0 && abd>0)) + return true; + continue; + } + return false; +} + } // end namespace GeoLib diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h index a991293f5091cf8c7885d8b7c7c5d2734cb04353..9623083b768a644a16d7ee871ba84a0336eea0e4 100644 --- a/GeoLib/AnalyticalGeometry.h +++ b/GeoLib/AnalyticalGeometry.h @@ -135,9 +135,22 @@ bool lineSegmentIntersect (const GeoLib::Point& a, const GeoLib::Point& b, */ GeoLib::Point* triangleLineIntersection(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c, GeoLib::Point const& p, GeoLib::Point const& q); -// Calculates the scalar triple (u x v) . w +/// Calculates the scalar triple (u x v) . w double scalarTriple(GeoLib::Point const& u, GeoLib::Point const& v, GeoLib::Point const& w); +/** + * Checks if a und b can be placed on a plane such that c and d lie on different sides of said plane. + * (In 2D space this checks if c and d are on different sides of a line through a and b.) + * @param a first point on plane + * @param b second point on plane + * @param c first point to test + * @param d second point to test + * @return true, if such a plane can be found, false otherwise + */ + bool dividedByPlane(const GeoLib::Point& a, const GeoLib::Point& b, + const GeoLib::Point& c, const GeoLib::Point& d); + + } // end namespace GeoLib #endif /* ANALYTICAL_GEOMETRY_H_ */ diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h index a609ff86fb0452c1305d176f1ad61ef34bc65577..170b503077fec459ed22514acd429cb040236256 100644 --- a/MeshLib/Elements/Element.h +++ b/MeshLib/Elements/Element.h @@ -16,6 +16,7 @@ #define ELEMENT_H_ #include <vector> +#include <limits> #include "MeshEnums.h" #include "Mesh.h" @@ -132,6 +133,11 @@ public: */ unsigned getValue() const { return _value; } + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const = 0; + /** * Set the index value for external information. * @param value an unsigned value for linking with external information diff --git a/MeshLib/Elements/TemplateHex.h b/MeshLib/Elements/TemplateHex.h index b98f5b072dd41955ca88a2e18848fc38bcd5724a..007882a75331609157318935141d295973cd8746 100644 --- a/MeshLib/Elements/TemplateHex.h +++ b/MeshLib/Elements/TemplateHex.h @@ -98,6 +98,11 @@ public: /// Returns true if these two indices form an edge and false otherwise bool isEdge(unsigned i, unsigned j) const; + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const; + /** * Method clone is inherited from class Element. It makes a deep copy of the Hex instance. * @return an exact copy of the object diff --git a/MeshLib/Elements/TemplateHex.tpp b/MeshLib/Elements/TemplateHex.tpp index 92d2075a67676865486b78f5605196fdb9fe1d30..f2e7538d8c215055b5f412ff2e61cd7605ab1195 100644 --- a/MeshLib/Elements/TemplateHex.tpp +++ b/MeshLib/Elements/TemplateHex.tpp @@ -155,6 +155,20 @@ unsigned TemplateHex<NNODES,CELLHEXTYPE>::identifyFace(Node* nodes[3]) const return std::numeric_limits<unsigned>::max(); } +template <unsigned NNODES, CellType CELLHEXTYPE> +bool TemplateHex<NNODES,CELLHEXTYPE>::isValid() const +{ + for (unsigned i=0; i<6; ++i) + { + const MeshLib::Quad* quad (dynamic_cast<const MeshLib::Quad*>(this->getFace(i))); + const bool quad_is_valid (quad->isValid()); + delete quad; + if (!quad_is_valid) + return false; + } + return true; +} + template <unsigned NNODES, CellType CELLHEXTYPE> Element* TemplateHex<NNODES,CELLHEXTYPE>::reviseElement() const { diff --git a/MeshLib/Elements/TemplateLine.h b/MeshLib/Elements/TemplateLine.h index 0b462c1d45b873cf767fc7f4f30304b0cf7db1ff..73095891d860e9ac3fc1141b68c9197a2a9009d8 100644 --- a/MeshLib/Elements/TemplateLine.h +++ b/MeshLib/Elements/TemplateLine.h @@ -72,6 +72,11 @@ public: */ virtual CellType getCellType() const { return CELLLINETYPE; } + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const { return (this->_length) > std::numeric_limits<double>::epsilon(); } + /** * Method clone is inherited from class Element. It makes a deep copy of the TemplateLine instance. * @return an exact copy of the object diff --git a/MeshLib/Elements/TemplatePrism.h b/MeshLib/Elements/TemplatePrism.h index eb14f544dc2bb3b3ef98b23cb76d75f566c9bf3e..343a32e4fcf1fcafbfb0a00688b86246ad78c979 100644 --- a/MeshLib/Elements/TemplatePrism.h +++ b/MeshLib/Elements/TemplatePrism.h @@ -96,6 +96,11 @@ public: /// Returns true if these two indeces form an edge and false otherwise bool isEdge(unsigned i, unsigned j) const; + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const; + /** * Method clone is inherited from class Element. It makes a deep copy of the * Hex instance employing the copy constructor of class TemplatePrism. diff --git a/MeshLib/Elements/TemplatePrism.tpp b/MeshLib/Elements/TemplatePrism.tpp index 328e31a83146e7a1ffea47890ed2259cce37f0f2..dafee3a33ec5f11ccd2507f088f371c49857fd09 100644 --- a/MeshLib/Elements/TemplatePrism.tpp +++ b/MeshLib/Elements/TemplatePrism.tpp @@ -164,6 +164,20 @@ unsigned TemplatePrism<NNODES,CELLPRISMTYPE>::identifyFace(Node* nodes[3]) const return std::numeric_limits<unsigned>::max(); } +template <unsigned NNODES, CellType CELLPRISMTYPE> +bool TemplatePrism<NNODES,CELLPRISMTYPE>::isValid() const +{ + for (unsigned i=1; i<4; ++i) + { + const MeshLib::Quad* quad (dynamic_cast<const MeshLib::Quad*>(this->getFace(i))); + const bool quad_is_valid (quad->isValid()); + delete quad; + if (!quad_is_valid) + return false; + } + return true; +} + template <unsigned NNODES, CellType CELLPRISMTYPE> Element* TemplatePrism<NNODES,CELLPRISMTYPE>::reviseElement() const { diff --git a/MeshLib/Elements/TemplatePyramid.h b/MeshLib/Elements/TemplatePyramid.h index e38d6e8c764c6d85d63414891d84c97c67dbc3fb..802bab5d46067b64d0f84b2bb5d192ecdba10476 100644 --- a/MeshLib/Elements/TemplatePyramid.h +++ b/MeshLib/Elements/TemplatePyramid.h @@ -94,6 +94,11 @@ public: /// Returns true if these two indeces form an edge and false otherwise bool isEdge(unsigned i, unsigned j) const; + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const; + /** * Method clone is inherited from class Element. It makes a deep copy of the * TemplatePyramid instance employing the copy constructor of class TemplatePyramid. diff --git a/MeshLib/Elements/TemplatePyramid.tpp b/MeshLib/Elements/TemplatePyramid.tpp index 8efa44a025cb13bd3cc41fed003f96034c7cce4d..0c136901dfc5c4f488ec14dff2f3e95dee449281 100644 --- a/MeshLib/Elements/TemplatePyramid.tpp +++ b/MeshLib/Elements/TemplatePyramid.tpp @@ -166,6 +166,16 @@ unsigned TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::identifyFace(Node* nodes[3]) c return std::numeric_limits<unsigned>::max(); } +template <unsigned NNODES, CellType CELLPYRAMIDTYPE> +bool TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::isValid() const +{ + const MeshLib::Quad* base (dynamic_cast<const MeshLib::Quad*>(this->getFace(4))); + const bool base_is_valid (base->isValid()); + delete base; + return base_is_valid; +} + + template <unsigned NNODES, CellType CELLPYRAMIDTYPE> Element* TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::reviseElement() const { diff --git a/MeshLib/Elements/TemplateQuad.h b/MeshLib/Elements/TemplateQuad.h index d58db41d8bae6afeac27758522bb502e177e4452..0a083f06117ab7cd023d220698739a0a0737eb89 100644 --- a/MeshLib/Elements/TemplateQuad.h +++ b/MeshLib/Elements/TemplateQuad.h @@ -97,6 +97,10 @@ public: */ virtual bool isPntInside(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const; + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const; /** * Method clone is inherited from class Element. It makes a deep copy of the TemplateQuad instance. diff --git a/MeshLib/Elements/TemplateQuad.tpp b/MeshLib/Elements/TemplateQuad.tpp index f632083310db0d93c5e1f96e238ae98fe0073742..8f8fe18782cad454db41105b811b163ab246cda8 100644 --- a/MeshLib/Elements/TemplateQuad.tpp +++ b/MeshLib/Elements/TemplateQuad.tpp @@ -119,6 +119,32 @@ unsigned TemplateQuad<NNODES,CELLQUADTYPE>::identifyFace(Node* nodes[3]) const return std::numeric_limits<unsigned>::max(); } +template <unsigned NNODES, CellType CELLQUADTYPE> +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()) + { + // check if quad is convex + if (GeoLib::dividedByPlane(*_nodes[0], *_nodes[2], *_nodes[1], *_nodes[3]) && + GeoLib::dividedByPlane(*_nodes[1], *_nodes[3], *_nodes[0], *_nodes[2])) + return true; + } + return false; +} + template <unsigned NNODES, CellType CELLQUADTYPE> Element* TemplateQuad<NNODES,CELLQUADTYPE>::reviseElement() const { diff --git a/MeshLib/Elements/TemplateTet.h b/MeshLib/Elements/TemplateTet.h index 72bc27b896167183cf4193571cdf44224270d17b..9abc25724cea775328c3e2f3aa684201a9c10c27 100644 --- a/MeshLib/Elements/TemplateTet.h +++ b/MeshLib/Elements/TemplateTet.h @@ -93,6 +93,11 @@ public: /// Returns true if these two indeces form an edge and false otherwise bool isEdge(unsigned i, unsigned j) const; + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const { return (this->_volume > std::numeric_limits<double>::epsilon()); } + /** * Method clone is inherited from class Element. It makes a deep copy of the TemplateTet instance. * @return an exact copy of the object diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h index a6acc86c3d5fa5f1852f83351d90df555a306a3a..4541b9a62093403e6a68aa05b389a6e70baafa49 100644 --- a/MeshLib/Elements/TemplateTri.h +++ b/MeshLib/Elements/TemplateTri.h @@ -94,6 +94,12 @@ public: */ virtual bool isPntInside(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const; + /** + * Tests if the element is geometrically valid, i.e. convex with volume > 0. + */ + virtual bool isValid() const { return (this->_area > std::numeric_limits<double>::epsilon()); } + + /** * Method clone is inherited from class Element. It makes a deep copy of the TemplateTri instance. * @return an exact copy of the object diff --git a/scripts/cmake/Functions.cmake b/scripts/cmake/Functions.cmake index 27b48eeac7e83afd7b46bc087f9b19c9e7800248..8a5161e64c9462cbd41183430f31f10a63f963d8 100644 --- a/scripts/cmake/Functions.cmake +++ b/scripts/cmake/Functions.cmake @@ -32,7 +32,8 @@ MACRO(GET_SOURCE_FILES SOURCE_FILES) GET_CURRENT_SOURCE_SUBDIRECTORY(DIRECTORY) SOURCE_GROUP( "${DIRECTORY}${DIR}" FILES ${GET_SOURCE_FILES_HEADERS} - ${GET_SOURCE_FILES_SOURCES}) + ${GET_SOURCE_FILES_SOURCES} + ${GET_SOURCE_FILES_TEMPLATES}) ENDMACRO()