#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),
{
}
+ /**
+ * 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
}
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
{
}
- /*
+ /**
* Destructor
*
* Deallocates the memory used to store the points of the polygons.
}
////////////////////////////////////////////////////////////////////////////
- /// 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
/**
* 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
// 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;
}
////////////////////////////////////////////////////////////////////////////////////
- /// Detection of (very) degenerate cases ////////////
+ // Detection of (very) degenerate cases /////////////
////////////////////////////////////////////////////////////////////////////////////
/**
* @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
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<int>(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))
{
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)
{
////////////////////////////////////////////////////////////////////////////////////
- /// 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
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
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<double>::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;
////////////////////////////////////////////////////////////////////////////////////
}
-
// -- (2) check that each double product statisfies Grandy, [47], else set to 0
for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1))
{
/**
* 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
}
#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
* @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
{
