From 1692e901cf5c9708a080048c327c458df2bcecdd Mon Sep 17 00:00:00 2001 From: vbd Date: Mon, 10 Sep 2007 13:05:51 +0000 Subject: [PATCH] staffan : * doc updates --- src/INTERP_KERNEL/IntersectorTetra.cxx | 2 +- src/INTERP_KERNEL/TransformedTriangle.cxx | 85 ++++++++++++++----- src/INTERP_KERNEL/TransformedTriangle.hxx | 8 +- .../TransformedTriangle_inline.hxx | 1 - .../TransformedTriangle_math.cxx | 26 ++++-- 5 files changed, 89 insertions(+), 33 deletions(-) diff --git a/src/INTERP_KERNEL/IntersectorTetra.cxx b/src/INTERP_KERNEL/IntersectorTetra.cxx index ac67374a3..5ccec2ab4 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.cxx +++ b/src/INTERP_KERNEL/IntersectorTetra.cxx @@ -225,7 +225,7 @@ namespace MEDMEM // reset if it is very small to keep the matrix sparse // is this a good idea? - if(epsilonEqual(totalVolume, 0.0, 1.0e-11)) + if(epsilonEqual(totalVolume, 0.0, 1.0e-14)) { totalVolume = 0.0; } diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 6986e509d..35528e788 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -11,12 +11,25 @@ #include "VectorUtils.hxx" - +/** + * Class representing a circular order of a set of points around their barycenter. + * It is used with the STL sort() algorithm to sort the point of the two polygons + * + */ class ProjectedCentralCircularSortOrder { -public: +public: + + /// Enumeration of different planes to project on when calculating order enum CoordType { XY, XZ, YZ }; + /** + * Constructor + * + * @param barycenter double[3] containing the barycenter of the points to be compared + * @param type plane to project on when comparing. The comparison will not work if all the points are in a plane perpendicular + * to the plane being projected on + */ ProjectedCentralCircularSortOrder(const double* barycenter, const CoordType type) : _aIdx((type == YZ) ? 2 : 0), _bIdx((type == XY) ? 1 : 2), @@ -25,6 +38,15 @@ public: { } + /** + * Comparison operator + * Compares the relative position between two points in their ordering around the barycenter. + * + * @param pt1 a double[3] representing a point + * @param pt2 a double[3] representing a point + * @return true if the angle of the difference vector between pt1 and the barycenter is greater than that + * of the difference vector between pt2 and the barycenter. + */ bool operator()(const double* pt1, const double* pt2) { // calculate angles with the axis @@ -35,19 +57,23 @@ public: } private: + /// indices of X, Y, Z coordinates in double[3] arrays to use : these depend on the projection plane const int _aIdx, _bIdx; + + /// values of projected coordinates for barycenter const double _a, _b; }; -class Vector3Cmp -{ -public: - bool operator()(double* const& pt1, double* const& pt2) - { - LOG(6, "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]))); - return (pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]); - } -}; + +//class Vector3Cmp +//{ +// public: +// bool operator()(double* const& pt1, double* const& pt2) +// { +// LOG(6, "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]))); +// return (pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]); +// } +// }; namespace INTERP_UTILS { @@ -98,7 +124,7 @@ namespace INTERP_UTILS } - /* + /** * Destructor * * Deallocates the memory used to store the points of the polygons. @@ -190,12 +216,13 @@ namespace INTERP_UTILS } //////////////////////////////////////////////////////////////////////////// - /// PRIVATE //////////////////////////////////////////////////////////////// + // PRIVATE ///////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// - /// High-level methods called directly by calculateIntersectionVolume() //////// + // High-level methods called directly by calculateIntersectionVolume() /////// //////////////////////////////////////////////////////////////////////////////////// + /** * Calculates the intersection polygons A and B, performing the intersection tests * and storing the corresponding points in the vectors _polygonA and _polygonB @@ -539,7 +566,6 @@ namespace INTERP_UTILS /** * Sorts the given intersection polygon in circular order around its barycenter. - * @pre the intersection polygons have been calculated with calculateIntersectionPolygons() * @post the vertices in _polygonA and _polygonB are sorted in circular order around their * respective barycenters @@ -571,11 +597,11 @@ namespace INTERP_UTILS // We keep the test here anyway, to avoid interdependency. // is triangle parallel to x == 0 ? - if(isTriangleInPlaneOfFacet(OZX)) + if(isTriangleParallelToFacet(OZX)) { type = SortOrder::YZ; } - else if(isTriangleInPlaneOfFacet(OYZ)) + else if(isTriangleParallelToFacet(OYZ)) { type = SortOrder::XZ; } @@ -637,7 +663,7 @@ namespace INTERP_UTILS //////////////////////////////////////////////////////////////////////////////////// - /// Detection of (very) degenerate cases //////////// + // Detection of (very) degenerate cases ///////////// //////////////////////////////////////////////////////////////////////////////////// /** @@ -646,7 +672,7 @@ namespace INTERP_UTILS * @param facet one of the facets of the tetrahedron * @return true if PQR lies in the plane of the facet, false if not */ - bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) + bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const { // coordinate to check @@ -663,12 +689,25 @@ namespace INTERP_UTILS return true; } - /* + /** + * Checks if the triangle is parallel to the given facet + * + * @param facet one of the facets of the unit tetrahedron + * @return true if triangle is parallel to facet, false if not + */ + bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const + { + // coordinate to check + const int coord = static_cast(facet); + return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]); + } + + /** * Determines whether the triangle is below the z-plane. * * @return true if the z-coordinate of the three corners of the triangle are all less than 0, false otherwise. */ - bool TransformedTriangle::isTriangleBelowTetraeder() + bool TransformedTriangle::isTriangleBelowTetraeder() const { for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) { @@ -681,11 +720,11 @@ namespace INTERP_UTILS return true; } - /* + /** * Prints the coordinates of the triangle to std::cout * */ - void TransformedTriangle::dumpCoords() + void TransformedTriangle::dumpCoords() const { std::cout << "Coords : "; for(int i = 0 ; i < 3; ++i) diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 391fc05c4..df91c099d 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -116,7 +116,7 @@ namespace INTERP_UTILS double calculateIntersectionVolume(); - void dumpCoords(); + void dumpCoords() const; private: @@ -136,9 +136,11 @@ namespace INTERP_UTILS /// Detection of degenerate triangles //////////// //////////////////////////////////////////////////////////////////////////////////// - bool isTriangleInPlaneOfFacet(const TetraFacet facet); + bool isTriangleInPlaneOfFacet(const TetraFacet facet) const; + + bool isTriangleParallelToFacet(const TetraFacet facet) const; - bool isTriangleBelowTetraeder(); + bool isTriangleBelowTetraeder() const; //////////////////////////////////////////////////////////////////////////////////// /// Intersection test methods and intersection point calculations //////// diff --git a/src/INTERP_KERNEL/TransformedTriangle_inline.hxx b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx index d3f2d80ab..8670fb65a 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_inline.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx @@ -13,7 +13,6 @@ * member variable array_triangleSurroundsEdgeCache. * */ - inline void TransformedTriangle::preCalculateTriangleSurroundsEdge() { for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index 6ed3a7477..44bfe1c11 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -16,11 +16,16 @@ namespace INTERP_UTILS { //////////////////////////////////////////////////////////////////////////////////// - /// Constants ///////////////// + // Tables ///////////////// //////////////////////////////////////////////////////////////////////////////////// + + /// Table with first coordinate (a) used to calculate double product c^pq_ab = p_a * q_b - p_b * q_a (index to be used : DoubleProduct) const int TransformedTriangle::DP_OFFSET_1[8] = {1, 2, 0, 2, 0, 1, 4, 1}; + + /// Table with second coordinate (a) used to calculate double product c^pq_ab = p_a * q_b - p_b * q_a (index to be used : DoubleProduct) const int TransformedTriangle::DP_OFFSET_2[8] = {2, 0, 1, 3, 3, 3, 0, 4}; + /// Coordinates used to calculate triple products by the expanding one of the three rows of the determinant (index to be used : 3*Corner + row) const int TransformedTriangle::COORDINATE_FOR_DETERMINANT_EXPANSION[12] = { // row 1, 2, 3 @@ -30,6 +35,7 @@ namespace INTERP_UTILS 0, 1, 3 // Z }; + /// Double products used to calculate triple products by expanding one of the three rows of the determinant (index to be used : 3*Corner + row) const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_DETERMINANT_EXPANSION[12] = { // row 1, 2, 3 @@ -39,11 +45,16 @@ namespace INTERP_UTILS C_YH, C_XH, C_XY // Z }; - //const double TransformedTriangle::MACH_EPS = 1.0e-15; + /// The machine epsilon, used in precision corrections const long double TransformedTriangle::MACH_EPS = std::numeric_limits::epsilon(); + + /// 4.0 * the machine epsilon, represents the precision of multiplication when performing corrections corrections ( f in Grandy ) const long double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS; + + /// Threshold for resetting double and triple products to zero; ( F / f in Grandy ) const long double TransformedTriangle::THRESHOLD_F = 20.0; + /// Threshold for what is considered a small enough angle to warrant correction of triple products by Grandy, [57] const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1; //////////////////////////////////////////////////////////////////////////////////// @@ -106,7 +117,6 @@ namespace INTERP_UTILS } - // -- (2) check that each double product statisfies Grandy, [47], else set to 0 for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1)) { @@ -152,7 +162,7 @@ namespace INTERP_UTILS /** * Checks if the double products for a given segment are consistent, as defined by - * Grandy, [46] + * Grandy, [46]. * * @param seg Segment for which to check consistency of double products * @return true if the double products are consistent, false if not @@ -186,6 +196,13 @@ namespace INTERP_UTILS } #ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx + + /** + * Sets the three double product associated with a given segment and a given corner to 0.0. + * + * @param seg a segment of the triangle + * @param corner a corner of the tetrahedron + */ void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner) { // set the three corresponding double products to 0.0 @@ -453,7 +470,6 @@ namespace INTERP_UTILS * @param project indicates whether or not to perform projection as inidicated in Grandy, p.446 * @return triple product associated with corner (see Grandy, [50]-[52]) */ - double TransformedTriangle::calcTByDevelopingRow(const TetraCorner corner, const int row, const bool project) const { -- 2.39.2