From 2e174c959335c1f55b4f0e86f9387635c10b4393 Mon Sep 17 00:00:00 2001 From: vbd Date: Mon, 10 Sep 2007 14:41:29 +0000 Subject: [PATCH] staffan : * doc updates --- src/INTERP_KERNEL/IntersectorTetra.cxx | 17 +- src/INTERP_KERNEL/IntersectorTetra.hxx | 25 +- src/INTERP_KERNEL/MeshUtils.hxx | 15 +- src/INTERP_KERNEL/TransformedTriangle.cxx | 711 +++++++++--------- src/INTERP_KERNEL/TransformedTriangle.hxx | 98 ++- .../TransformedTriangle_inline.hxx | 20 +- .../TransformedTriangle_intersect.cxx | 72 +- .../TransformedTriangle_math.cxx | 31 +- src/INTERP_KERNEL/VectorUtils.hxx | 25 +- 9 files changed, 503 insertions(+), 511 deletions(-) diff --git a/src/INTERP_KERNEL/IntersectorTetra.cxx b/src/INTERP_KERNEL/IntersectorTetra.cxx index 5ccec2ab4..33baf6eec 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.cxx +++ b/src/INTERP_KERNEL/IntersectorTetra.cxx @@ -17,12 +17,12 @@ using namespace INTERP_UTILS; namespace MEDMEM { - /* + /** * Constructor * * @param srcMesh mesh containing the source elements * @param targetMesh mesh containing the target elements - * + * @param targetCell global number of the target cell */ IntersectorTetra::IntersectorTetra(const MESH& srcMesh, const MESH& targetMesh, int targetCell) : _srcMesh(srcMesh), filtered(0) @@ -44,7 +44,7 @@ namespace MEDMEM } - /* + /** * Destructor * */ @@ -58,18 +58,19 @@ namespace MEDMEM } } - /* - * Calculates the volume of intersection of an element in the source mesh and an element - * in the target mesh. The method is based on the algorithm of Grandy. It first calculates the transformation - * that takes the target tetrahedron into the unit tetrahedron. After that, the + /** + * Calculates the volume of intersection of an element in the source mesh and the target element. + * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the * faces of the source element are triangulated and the calculated transformation is applied * to each triangle. The algorithm of Grandy, implemented in INTERP_UTILS::TransformedTriangle is used * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions * divided by the determinant of the transformation. * + * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated + * with triangulated faces to avoid having to recalculate these. + * * @pre The element in _targetMesh referenced by targetCell is of type MED_TETRA4. * @param srcCell global number of the source element (1 <= srcCell < # source cells) - * @param targetCell global number of the target element (1 <= targetCell < # target cells) - this element must be a tetrahedron */ double IntersectorTetra::intersectSourceCell(int element) { diff --git a/src/INTERP_KERNEL/IntersectorTetra.hxx b/src/INTERP_KERNEL/IntersectorTetra.hxx index acb7ed885..df639fad6 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.hxx +++ b/src/INTERP_KERNEL/IntersectorTetra.hxx @@ -25,9 +25,6 @@ namespace INTERP_UTILS TriangleFaceKey(int node1, int node2, int node3) { sort3Ints(_nodes, node1, node2, node3); - // assert(_nodes[0] < _nodes[1]); - //assert(_nodes[0] < _nodes[2]); - //assert(_nodes[1] < _nodes[2]); _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29; } @@ -62,7 +59,7 @@ namespace INTERP_UTILS { - /* + /** * Class calculating the volume of intersection of individual 3D elements. * */ @@ -80,6 +77,12 @@ namespace INTERP_UTILS private: + /// disallow copying + IntersectorTetra(const IntersectorTetra& t); + + /// disallow assignment + IntersectorTetra& operator=(const IntersectorTetra& t); + TetraAffineTransform* _t; hash_map< int, double*> _nodes; @@ -88,8 +91,6 @@ namespace INTERP_UTILS inline void checkIsOutside(const double* pt, bool* isOutside) const; inline void calculateNode(int globalNodeNum); inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); - // inline void sort3Ints(int* sorted, int node1, int node2, int node3); - //inline std::string createFaceKey(int node1, int node2, int node3); const MEDMEM::MESH& _srcMesh; @@ -162,18 +163,6 @@ namespace INTERP_UTILS } } } - -#if 0 - inline std::string IntersectorTetra::createFaceKey(int node1, int node2, int node3) - { - int sorted[3]; - sort3Ints(sorted, node1, node2, node3); - - std::stringstream sstr; - sstr << node1 << "-" << node2 << "-" << node3; - return sstr.str(); - } -#endif }; diff --git a/src/INTERP_KERNEL/MeshUtils.hxx b/src/INTERP_KERNEL/MeshUtils.hxx index ae8d81a97..9677a5444 100644 --- a/src/INTERP_KERNEL/MeshUtils.hxx +++ b/src/INTERP_KERNEL/MeshUtils.hxx @@ -23,11 +23,7 @@ namespace INTERP_UTILS */ inline const double* getCoordsOfNode(int node, int element, const MESH& mesh) { - assert(node >= 1); - assert(node <= mesh.getNumberOfNodes()); - const int nodeOffset = node - 1; - const int elemIdx = mesh.getConnectivityIndex(MED_NODAL, MED_CELL)[element - 1] - 1; - const int connIdx = mesh.getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_CELL, MED_ALL_ELEMENTS)[elemIdx + nodeOffset] - 1; + const int connIdx = getGlobalNumberOfNode(node, element, mesh); return &(mesh.getCoordinates(MED_FULL_INTERLACE)[3*connIdx]); } @@ -47,6 +43,15 @@ namespace INTERP_UTILS return static_cast(type) - 300; } + /** + * Returns the global number of the node of an element. + * (1 <= node <= #nodes of element) + * + * @param node the node for which the global number is sought + * @param element an element of the mesh + * @param mesh a mesh + * @return the node's global number (its coordinates in the coordinates array are at [3*globalNumber, 3*globalNumber + 2] + */ inline int getGlobalNumberOfNode(int node, int element, const MESH& mesh) { assert(node >= 1); diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index b46ed2c7c..f9e8059a5 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -10,83 +10,72 @@ #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 +namespace INTERP_UTILS { -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), - _a(barycenter[_aIdx]), - _b(barycenter[_bIdx]) - { - } /** - * Comparison operator. - * Compares the relative position between two points in their ordering around the barycenter. + * 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 * - * @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) + class ProjectedCentralCircularSortOrder { - // calculate angles with the axis - const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b); - const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b); + public: - return ang1 > ang2; - } - -private: - /// index corresponding to first coordinate of plane on which points are projected - const int _aIdx; + /// Enumeration of different planes to project on when calculating order + enum CoordType { XY, XZ, YZ }; - /// index corresponding to second coordinate of plane on which points are projected - const int _bIdx; + /** + * 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), + _a(barycenter[_aIdx]), + _b(barycenter[_bIdx]) + { + } - /// value of first projected coordinate of the barycenter - const double _a; - - /// value of second projected coordinate of the barycenter - const double _b; -}; + /** + * 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 + const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b); + const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b); + return ang1 > ang2; + } -//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]); -// } -// }; + private: + /// index corresponding to first coordinate of plane on which points are projected + const int _aIdx; + + /// index corresponding to second coordinate of plane on which points are projected + const int _bIdx; -namespace INTERP_UTILS -{ - //////////////////////////////////////////////////////////////////////////// - /// PUBLIC //////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////// + /// value of first projected coordinate of the barycenter + const double _a; + + /// value of second projected coordinate of the barycenter + const double _b; + }; + // ---------------------------------------------------------------------------------- + // TransformedTriangle PUBLIC + // ---------------------------------------------------------------------------------- + /** * Constructor * @@ -205,29 +194,25 @@ namespace INTERP_UTILS double volB = 0.0; // if triangle is not in h = 0 plane, calculate volume under B - if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) - { - LOG(2, "---- Treating polygon B ... "); + if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) + { + LOG(2, "---- Treating polygon B ... "); - calculatePolygonBarycenter(B, barycenter); - sortIntersectionPolygon(B, barycenter); - volB = calculateVolumeUnderPolygon(B, barycenter); - LOG(2, "Volume is " << sign * volB); - } + calculatePolygonBarycenter(B, barycenter); + sortIntersectionPolygon(B, barycenter); + volB = calculateVolumeUnderPolygon(B, barycenter); + LOG(2, "Volume is " << sign * volB); + } LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********"); return sign * (volA + volB); } - - // ////////////////////////////////////////////////////////////////////////// - // PRIVATE ///////////////////////////////////////////////////////////////// - // ////////////////////////////////////////////////////////////////////////// - - // ////////////////////////////////////////////////////////////////////////////////// - // High-level methods called directly by calculateIntersectionVolume() - // ////////////////////////////////////////////////////////////////////////////////// + + // ---------------------------------------------------------------------------------- + // TransformedTriangle PRIVATE + // ---------------------------------------------------------------------------------- /** * Calculates the intersection polygons A and B, performing the intersection tests @@ -367,377 +352,377 @@ namespace INTERP_UTILS #else - // -- segment intersections - for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1)) - { - - // segment - facet - for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1)) + // -- segment intersections + for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1)) { - if(testSegmentFacetIntersection(seg, facet)) + + // segment - facet + for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1)) { - double* ptA = new double[3]; - calcIntersectionPtSegmentFacet(seg, facet, ptA); - _polygonA.push_back(ptA); - LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); - if(facet == XYZ) + if(testSegmentFacetIntersection(seg, facet)) { - double* ptB = new double[3]; - copyVector3(ptA, ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); - } + double* ptA = new double[3]; + calcIntersectionPtSegmentFacet(seg, facet, ptA); + _polygonA.push_back(ptA); + LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); + if(facet == XYZ) + { + double* ptB = new double[3]; + copyVector3(ptA, ptB); + _polygonB.push_back(ptB); + LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); + } + } } - } - // segment - edge - for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) - { - if(testSegmentEdgeIntersection(seg, edge)) + // segment - edge + for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) { - double* ptA = new double[3]; - calcIntersectionPtSegmentEdge(seg, edge, ptA); - _polygonA.push_back(ptA); - LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A"); - if(edge >= XY) + if(testSegmentEdgeIntersection(seg, edge)) { - double* ptB = new double[3]; - copyVector3(ptA, ptB); - _polygonB.push_back(ptB); + double* ptA = new double[3]; + calcIntersectionPtSegmentEdge(seg, edge, ptA); + _polygonA.push_back(ptA); + LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A"); + if(edge >= XY) + { + double* ptB = new double[3]; + copyVector3(ptA, ptB); + _polygonB.push_back(ptB); + } } } - } - // segment - corner - for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) - { - if(testSegmentCornerIntersection(seg, corner)) + // segment - corner + for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) { - double* ptA = new double[3]; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptA); - _polygonA.push_back(ptA); - LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A"); - if(corner != O) + if(testSegmentCornerIntersection(seg, corner)) { - double* ptB = new double[3]; - _polygonB.push_back(ptB); - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B"); + double* ptA = new double[3]; + copyVector3(&COORDS_TET_CORNER[3 * corner], ptA); + _polygonA.push_back(ptA); + LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A"); + if(corner != O) + { + double* ptB = new double[3]; + _polygonB.push_back(ptB); + copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); + LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B"); + } } } - } #endif #ifdef OPTIMIZE - // segment - ray - for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) - { - if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner)) + // segment - ray + for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) { - double* ptB = new double[3]; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); + if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner)) + { + double* ptB = new double[3]; + copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); + _polygonB.push_back(ptB); + LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); + } } - } - // segment - halfstrip - for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1)) - { + // segment - halfstrip + for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1)) + { #if 0 - const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX - const bool doTest = - !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] && - !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]]; + const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX + const bool doTest = + !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] && + !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]]; - if(doTest && testSegmentHalfstripIntersection(seg, edge)) + if(doTest && testSegmentHalfstripIntersection(seg, edge)) #endif - if(testSegmentHalfstripIntersection(seg, edge)) - { - double* ptB = new double[3]; - calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); + if(testSegmentHalfstripIntersection(seg, edge)) + { + double* ptB = new double[3]; + calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); + _polygonB.push_back(ptB); + LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); + } } - } #else - // segment - ray - for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) - { - if(testSegmentRayIntersection(seg, corner)) + // segment - ray + for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) { - double* ptB = new double[3]; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); + if(testSegmentRayIntersection(seg, corner)) + { + double* ptB = new double[3]; + copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); + _polygonB.push_back(ptB); + LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); + } } - } - // segment - halfstrip - for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1)) - { - if(testSegmentHalfstripIntersection(seg, edge)) + // segment - halfstrip + for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1)) { - double* ptB = new double[3]; - calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); + if(testSegmentHalfstripIntersection(seg, edge)) + { + double* ptB = new double[3]; + calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); + _polygonB.push_back(ptB); + LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); + } } - } #endif - } + } - // inclusion tests - for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1)) - { - // { XYZ - inclusion only possible if in Tetrahedron? - // tetrahedron - if(testCornerInTetrahedron(corner)) + // inclusion tests + for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1)) { - double* ptA = new double[3]; - copyVector3(&_coords[5*corner], ptA); - _polygonA.push_back(ptA); - LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); - } + // { XYZ - inclusion only possible if in Tetrahedron? + // tetrahedron + if(testCornerInTetrahedron(corner)) + { + double* ptA = new double[3]; + copyVector3(&_coords[5*corner], ptA); + _polygonA.push_back(ptA); + LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); + } - // XYZ - plane - if(testCornerOnXYZFacet(corner)) - { - double* ptB = new double[3]; - copyVector3(&_coords[5*corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); - } + // XYZ - plane + if(testCornerOnXYZFacet(corner)) + { + double* ptB = new double[3]; + copyVector3(&_coords[5*corner], ptB); + _polygonB.push_back(ptB); + LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); + } + + // projection on XYZ - facet + if(testCornerAboveXYZFacet(corner)) + { + double* ptB = new double[3]; + copyVector3(&_coords[5*corner], ptB); + ptB[2] = 1 - ptB[0] - ptB[1]; + assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); + _polygonB.push_back(ptB); + LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); + } - // projection on XYZ - facet - if(testCornerAboveXYZFacet(corner)) - { - double* ptB = new double[3]; - copyVector3(&_coords[5*corner], ptB); - ptB[2] = 1 - ptB[0] - ptB[1]; - assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); - _polygonB.push_back(ptB); - LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); } } - } - - /** - * Calculates the barycenters of the given intersection polygon. - * - * @pre the intersection polygons have been calculated with calculateIntersectionPolygons() - * - * @param poly one of the two intersection polygons - * @param barycenter array of three doubles where barycenter is stored - * - */ - void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter) - { - LOG(3,"--- Calculating polygon barycenter"); + /** + * Calculates the barycenters of the given intersection polygon. + * + * @pre the intersection polygons have been calculated with calculateIntersectionPolygons() + * + * @param poly one of the two intersection polygons + * @param barycenter array of three doubles where barycenter is stored + * + */ + void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter) + { + LOG(3,"--- Calculating polygon barycenter"); - // get the polygon points - std::vector& polygon = (poly == A) ? _polygonA : _polygonB; + // get the polygon points + std::vector& polygon = (poly == A) ? _polygonA : _polygonB; - // calculate barycenter - const int m = polygon.size(); + // calculate barycenter + const int m = polygon.size(); - for(int j = 0 ; j < 3 ; ++j) - { - barycenter[j] = 0.0; - } + for(int j = 0 ; j < 3 ; ++j) + { + barycenter[j] = 0.0; + } - if(m != 0) - { - for(int i = 0 ; i < m ; ++i) + if(m != 0) { - const double* pt = polygon[i]; - for(int j = 0 ; j < 3 ; ++j) + for(int i = 0 ; i < m ; ++i) { - barycenter[j] += pt[j] / double(m); + const double* pt = polygon[i]; + for(int j = 0 ; j < 3 ; ++j) + { + barycenter[j] += pt[j] / double(m); + } } } + LOG(3,"Barycenter is " << vToStr(barycenter)); } - LOG(3,"Barycenter is " << vToStr(barycenter)); - } - /** - * 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 - * - * @param poly one of the two intersection polygons - * @param barycenter array of three doubles with the coordinates of the barycenter - * - */ - void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter) - { - LOG(3,"--- Sorting polygon ..."); - - using ::ProjectedCentralCircularSortOrder; - typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor - typedef SortOrder::CoordType CoordType; + /** + * 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 + * + * @param poly one of the two intersection polygons + * @param barycenter array of three doubles with the coordinates of the barycenter + * + */ + void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter) + { + LOG(3,"--- Sorting polygon ..."); - // get the polygon points - std::vector& polygon = (poly == A) ? _polygonA : _polygonB; + using ::ProjectedCentralCircularSortOrder; + typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor + typedef SortOrder::CoordType CoordType; - if(polygon.size() == 0) - return; + // get the polygon points + std::vector& polygon = (poly == A) ? _polygonA : _polygonB; - // determine type of sorting - CoordType type = SortOrder::XY; - if(poly == A) // B is on h = 0 plane -> ok - { - // NB : the following test is never true if we have eliminated the - // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume(). - // We keep the test here anyway, to avoid interdependency. + if(polygon.size() == 0) + return; - // is triangle parallel to x == 0 ? - if(isTriangleParallelToFacet(OZX)) + // determine type of sorting + CoordType type = SortOrder::XY; + if(poly == A) // B is on h = 0 plane -> ok { - type = SortOrder::YZ; - } - else if(isTriangleParallelToFacet(OYZ)) - { - type = SortOrder::XZ; + // NB : the following test is never true if we have eliminated the + // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume(). + // We keep the test here anyway, to avoid interdependency. + + // is triangle parallel to x == 0 ? + if(isTriangleParallelToFacet(OZX)) + { + type = SortOrder::YZ; + } + else if(isTriangleParallelToFacet(OYZ)) + { + type = SortOrder::XZ; + } } - } - // create order object - SortOrder order(barycenter, type); + // create order object + SortOrder order(barycenter, type); - // sort vector with this object - // NB : do not change place of first object, with respect to which the order - // is defined - sort((polygon.begin()), polygon.end(), order); + // sort vector with this object + // NB : do not change place of first object, with respect to which the order + // is defined + sort((polygon.begin()), polygon.end(), order); - LOG(3,"Sorted polygon is "); - for(int i = 0 ; i < polygon.size() ; ++i) - { - LOG(3,vToStr(polygon[i])); - } + LOG(3,"Sorted polygon is "); + for(int i = 0 ; i < polygon.size() ; ++i) + { + LOG(3,vToStr(polygon[i])); + } - } + } - /** - * Calculates the volume between the given polygon and the z = 0 plane. - * - * @pre the intersection polygones have been calculated with calculateIntersectionPolygons(), - * and they have been sorted in circular order with sortIntersectionPolygons(void) - * - * @param poly one of the two intersection polygons - * @param barycenter array of three doubles with the coordinates of the barycenter - * @return the volume between the polygon and the z = 0 plane - * - */ - double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter) - { - LOG(2,"--- Calculating volume under polygon"); + /** + * Calculates the volume between the given polygon and the z = 0 plane. + * + * @pre the intersection polygones have been calculated with calculateIntersectionPolygons(), + * and they have been sorted in circular order with sortIntersectionPolygons(void) + * + * @param poly one of the two intersection polygons + * @param barycenter array of three doubles with the coordinates of the barycenter + * @return the volume between the polygon and the z = 0 plane + * + */ + double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter) + { + LOG(2,"--- Calculating volume under polygon"); - // get the polygon points - std::vector& polygon = (poly == A) ? _polygonA : _polygonB; + // get the polygon points + std::vector& polygon = (poly == A) ? _polygonA : _polygonB; - double vol = 0.0; - const int m = polygon.size(); + double vol = 0.0; + const int m = polygon.size(); - for(int i = 0 ; i < m ; ++i) - { - const double* ptCurr = polygon[i]; // pt "i" - const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0) + for(int i = 0 ; i < m ; ++i) + { + const double* ptCurr = polygon[i]; // pt "i" + const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0) - const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2]; - const double factor2 = - ptCurr[0]*(ptNext[1] - barycenter[1]) - + ptNext[0]*(barycenter[1] - ptCurr[1]) - + barycenter[0]*(ptCurr[1] - ptNext[1]); - vol += (factor1 * factor2) / 6.0; - } + const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2]; + const double factor2 = + ptCurr[0]*(ptNext[1] - barycenter[1]) + + ptNext[0]*(barycenter[1] - ptCurr[1]) + + barycenter[0]*(ptCurr[1] - ptNext[1]); + vol += (factor1 * factor2) / 6.0; + } - LOG(2,"Abs. Volume is " << vol); - return vol; - } + LOG(2,"Abs. Volume is " << vol); + return vol; + } - //////////////////////////////////////////////////////////////////////////////////// - // Detection of (very) degenerate cases ///////////// - //////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////// + // Detection of (very) degenerate cases ///////////// + //////////////////////////////////////////////////////////////////////////////////// - /** - * Checks if the triangle lies in the plane of a given facet - * - * @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) const - { + /** + * Checks if the triangle lies in the plane of a given facet + * + * @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) const + { - // coordinate to check - const int coord = static_cast(facet); + // coordinate to check + const int coord = static_cast(facet); - for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) - { - if(_coords[5*c + coord] != 0.0) + for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) { - return false; + if(_coords[5*c + coord] != 0.0) + { + return false; + } } - } - return true; - } + 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]); - } + /** + * 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() const - { - for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) + /** + * 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() const { - // check z-coords for all points - if(_coords[5*c + 2] >= 0.0) + for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) { - return false; + // check z-coords for all points + if(_coords[5*c + 2] >= 0.0) + { + return false; + } } + return true; } - return true; - } - /** - * Prints the coordinates of the triangle to std::cout - * - */ - void TransformedTriangle::dumpCoords() const - { - std::cout << "Coords : "; - for(int i = 0 ; i < 3; ++i) + /** + * Prints the coordinates of the triangle to std::cout + * + */ + void TransformedTriangle::dumpCoords() const { - std::cout << vToStr(&_coords[5*i]) << ","; + std::cout << "Coords : "; + for(int i = 0 ; i < 3; ++i) + { + std::cout << vToStr(&_coords[5*i]) << ","; + } + std::cout << std::endl; } - std::cout << std::endl; - } -}; // NAMESPACE + }; // NAMESPACE diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 73c5a8ead..8152e894c 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -3,9 +3,13 @@ #include + + #ifdef OPTIMIZE +/// OPT_INLINE will be replaced by "inline" if OPTIMIZE is defined and else nothing #define OPT_INLINE inline #else +/// OPT_INLINE will be replaced by "inline" if OPTIMIZE is defined and else nothing #define OPT_INLINE #endif @@ -71,7 +75,9 @@ namespace INTERP_UTILS * When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it * is not known directly. It is then added to the polygon A and/or B as necessary. * - * + * OPTIMIZE : + * If OPTIMIZE is defined, a large number of methods will be prefixed with inline and some optimizations concerning the tests + * with zero double products will be used. */ class TransformedTriangle { @@ -119,6 +125,10 @@ namespace INTERP_UTILS private: + + // ---------------------------------------------------------------------------------- + // High-level methods called directly by calculateIntersectionVolume() + // ---------------------------------------------------------------------------------- void calculateIntersectionPolygons(); void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter); @@ -127,9 +137,9 @@ namespace INTERP_UTILS double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter); - //////////////////////////////////////////////////////////////////////////////////// - /// Detection of degenerate triangles //////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Detection of degenerate triangles + // ---------------------------------------------------------------------------------- bool isTriangleInPlaneOfFacet(const TetraFacet facet) const; @@ -137,9 +147,9 @@ namespace INTERP_UTILS bool isTriangleBelowTetraeder() const; - //////////////////////////////////////////////////////////////////////////////////// - /// Intersection test methods and intersection point calculations //////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Intersection test methods and intersection point calculations + // ---------------------------------------------------------------------------------- OPT_INLINE bool testSurfaceEdgeIntersection(const TetraEdge edge) const; @@ -169,9 +179,9 @@ namespace INTERP_UTILS OPT_INLINE bool testCornerAboveXYZFacet(const TriCorner corner) const; - //////////////////////////////////////////////////////////////////////////////////// - /// Utility methods used in intersection tests /////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Utility methods used in intersection tests + // ---------------------------------------------------------------------------------- bool testTriangleSurroundsEdge(const TetraEdge edge) const; @@ -187,14 +197,16 @@ namespace INTERP_UTILS bool testTriangleSurroundsRay(const TetraCorner corner) const; - //////////////////////////////////////////////////////////////////////////////////// - /// Double and triple product calculations /////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Double and triple product calculations + // ---------------------------------------------------------------------------------- - void preCalculateDoubleProducts(void); + bool areDoubleProductsConsistent(const TriSegment seg) const; + void preCalculateDoubleProducts(void); + OPT_INLINE void resetDoubleProducts(const TriSegment seg, const TetraCorner corner); double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const; @@ -211,48 +223,63 @@ namespace INTERP_UTILS double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const; - //////////////////////////////////////////////////////////////////////////////////// - /// Member variables /////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Member variables + // ---------------------------------------------------------------------------------- private: - // order : - // [ p_x, p_y, p_z, p_h, p_H, q_x, q_y, q_z, q_h, q_H, r_x, r_y, r_z, r_h, r_H ] + /// Array holding the coordinates of the triangle's three corners + /// order : + /// [ p_x, p_y, p_z, p_h, p_H, q_x, q_y, q_z, q_h, q_H, r_x, r_y, r_z, r_h, r_H ] double _coords[15]; - /// flags showing whether the double and triple products have been precalculated for this class - bool _isDoubleProductsCalculated, _isTripleProductsCalculated; + /// Flag showing whether the double products have been calculated yet + bool _isDoubleProductsCalculated; + + /// Flag showing whether the triple products have been calculated yet + bool _isTripleProductsCalculated; - /// array containing the 24 double products + /// Array containing the 24 double products. /// order : c^PQ_YZ, ... ,cPQ_10, ... c^QR_YZ, ... c^RP_YZ /// following order in enumeration DoubleProduct double _doubleProducts[24]; - /// array containing the 4 triple products + /// Array containing the 4 triple products. /// order : t_O, t_X, t_Y, t_Z double _tripleProducts[4]; - /// arrays holding the points in the two intersection polygons A and B + /// Vector holding the points of the intersection polygon A. /// these points are allocated in calculateIntersectionPolygons() and liberated in the destructor - std::vector _polygonA, _polygonB; + std::vector _polygonA; - /// vectors holding the coordinates of the barycenters of the polygons A and B - /// these points are calculated in calculatePolygonBarycenter - double _barycenterA[3], _barycenterB[3]; + /// Vector holding the points of the intersection polygon B. + /// These points are allocated in calculateIntersectionPolygons() and liberated in the destructor + std::vector _polygonB; + + /// Array holding the coordinates of the barycenter of the polygon A + /// This point is calculated in calculatePolygonBarycenter + double _barycenterA[3]; + + /// Array holding the coordinates of the barycenter of the polygon B + /// This point is calculated in calculatePolygonBarycenter + double _barycenterB[3]; - // used for debugging + /// Array of flags indicating which of the four triple products have been correctly calculated. + /// Used for asserts in debug mode bool _validTP[4]; #ifdef OPTIMIZE - void preCalculateTriangleSurroundsEdge(); + void preCalculateTriangleSurroundsEdge(); + + /// Array holding results of the test testTriangleSurroundsEdge() for all the edges. + /// These are calculated in preCalculateTriangleSurroundsEdge(). bool _triangleSurroundsEdgeCache[NO_TET_EDGE]; - bool _isOutsideTetra; #endif - //////////////////////////////////////////////////////////////////////////////////// - /// Constants ///////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Constants + // ---------------------------------------------------------------------------------- // offsets : 0 -> x, 1 -> y, 2 -> z, 3 -> h, 4 -> H // corresponds to order of double products in DoubleProduct @@ -309,9 +336,6 @@ namespace INTERP_UTILS // double products used in segment - ray test static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21]; - - inline bool isTriangleOutsideTetra(void) const; - }; diff --git a/src/INTERP_KERNEL/TransformedTriangle_inline.hxx b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx index 8670fb65a..2953118c3 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_inline.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx @@ -5,9 +5,10 @@ // It replaces those methods if OPTIMIZE is defined. // NB : most of these methods have documentation in their corresponding .cxx - file. -//////////////////////////////////////////////////////////////////////////////////////// -/// Optimization methods. These are only defined and used if OPTIMIZE is defined. -//////////////////////////////////////////////////////////////////////////////////////// +// ---------------------------------------------------------------------------------- +// Optimization methods. These are only defined and used if OPTIMIZE is defined. +// ----------------------------------------------------------------------------------- + /** * Calls TransformedTriangle::testTriangleSurroundsEdge for edges OX to ZX and stores the result in * member variable array_triangleSurroundsEdgeCache. @@ -22,9 +23,10 @@ inline void TransformedTriangle::preCalculateTriangleSurroundsEdge() } -//////////////////////////////////////////////////////////////////////////////////////// -/// TransformedTriangle_math.cxx /////// -//////////////////////////////////////////////////////////////////////////////////////// +// ---------------------------------------------------------------------------------- +// TransformedTriangle_math.cxx +// ---------------------------------------------------------------------------------- + inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner) { // set the three corresponding double products to 0.0 @@ -71,9 +73,9 @@ inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const Dou return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; } -//////////////////////////////////////////////////////////////////////////////////////// -/// TransformedTriangle_intersect.cxx /////// -//////////////////////////////////////////////////////////////////////////////////////// +// ---------------------------------------------------------------------------------- +// TransformedTriangle_intersect.cxx +// ---------------------------------------------------------------------------------- inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const { return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge); diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index 2605bb3c5..b198dea92 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -5,17 +5,15 @@ #include #include "VectorUtils.hxx" -#define SEG_RAY_TABLE 1 // seems correct - namespace INTERP_UTILS { - //////////////////////////////////////////////////////////////////////////////////// - /// Correspondance tables describing all the variations of formulas. ////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Correspondance tables describing all the variations of formulas. + // ---------------------------------------------------------------------------------- - // correspondance facet - double product - // Grandy, table IV + /// Correspondance between facets and double products. + /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] = { C_XH, C_XY, C_ZX, // OYZ @@ -24,7 +22,8 @@ namespace INTERP_UTILS C_XH, C_YH, C_ZH // XYZ }; - // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION + /// Signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION + /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] = { 1.0, 1.0, -1.0, @@ -33,7 +32,8 @@ namespace INTERP_UTILS 1.0, 1.0, 1.0 }; - // coordinates of corners of tetrahedron + /// Coordinates of corners of tetrahedron. + /// Use 3*Corner + coordinate as index const double TransformedTriangle::COORDS_TET_CORNER[12] = { 0.0, 0.0, 0.0, @@ -42,10 +42,10 @@ namespace INTERP_UTILS 0.0, 0.0, 1.0 }; - // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION - // for the calculation of the coordinates (x,y,z) of the intersection points - // for Segment-Facet and Segment-Edge intersections - // -1 indicates that the coordinate is 0 + /// Indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION + /// for the calculation of the coordinates (x,y,z) of the intersection points + /// for Segment-Facet and Segment-Edge intersections. + /// Use 3*facet + coordinate as index. -1 indicates that the coordinate is 0. const int TransformedTriangle::DP_INDEX[12] = { // x, y, z @@ -55,7 +55,9 @@ namespace INTERP_UTILS 9, 10, 11 // XYZ }; - // correspondance edge - corners + /// Correspondance edge - corners + /// Gives the two corners associated with each edge + /// Use 2*edge + {0, 1} as index const TransformedTriangle::TetraCorner TransformedTriangle::CORNERS_FOR_EDGE[12] = { O, X, // OX @@ -66,8 +68,8 @@ namespace INTERP_UTILS Z, X // ZX }; - // correspondance edge - facets - // facets shared by each edge + /// Correspondance edge - facets. + /// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] = { OXY, OZX, // OX @@ -78,7 +80,8 @@ namespace INTERP_UTILS OZX, XYZ // ZX }; - // edges meeting at a given corner + /// Correspondance corners - edges + /// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] = { OX, OY, OZ, // O @@ -87,7 +90,10 @@ namespace INTERP_UTILS OZ, ZX, YZ // Z }; - // NB : some uncertainty whether these last are correct + /// Double products to use in halfstrip intersection tests + /// Use 4*(offset_edge) + {0,1,2,3} as index. offset_edge = edge - 3 (so that XY -> 0, YZ -> 1, ZX -> 2) + /// Entries with offset 0 and 1 are for the first condition (positive product) + /// and those with offset 2 and 3 are for the second condition (negative product). const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_HALFSTRIP_INTERSECTION[12] = { C_10, C_01, C_ZH, C_10, // XY @@ -95,30 +101,20 @@ namespace INTERP_UTILS C_XY, C_10, C_YH, C_XY // ZX }; - // double products to use in segment-ray test - // dp 1 -> cond 1 - // dp 2-7 -> cond 3 -#if SEG_RAY_TABLE==1 + /// Double products to use in segment-ray test + /// Use 7*corner_offset + {0,1,2,3,4,5,6} as index. corner_offset = corner - 1 (so that X -> 0, Y-> 1, Z->2) + /// Entries with offset 0 are for first condition (zero double product) and the rest are for condition 3 (in the same + /// order as in the article) const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] = { C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01 // Z }; -#else - - const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] = - { - C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X - C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_YZ, // Y - C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_ZX // Z - }; -#endif - - //////////////////////////////////////////////////////////////////////////////////// - /// Intersection test methods and intersection point calculations //////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Intersection test methods and intersection point calculations + // ---------------------------------------------------------------------------------- #ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the given edge of the tetrahedron intersects the triangle PQR. (Grandy, eq [17]) @@ -691,9 +687,9 @@ namespace INTERP_UTILS #endif - //////////////////////////////////////////////////////////////////////////////////// - /// Utility methods used in intersection tests /////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ///////////////////////////////////////////////////////////////////////////////// + // Utility methods used in intersection tests /////////////// + // ///////////////////////////////////////////////////////////////////////////////// /** * Tests if the triangle PQR surrounds the axis on which the * given edge of the tetrahedron lies. diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index 44bfe1c11..7abb59412 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -9,20 +9,17 @@ #include "VectorUtils.hxx" -#undef SECOND_CORNER_RESET -#undef FIXED_DELTA 1.0e-15 - namespace INTERP_UTILS { - //////////////////////////////////////////////////////////////////////////////////// - // Tables ///////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // 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) + /// Table with second coordinate (b) 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) @@ -57,9 +54,9 @@ namespace INTERP_UTILS /// 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; - //////////////////////////////////////////////////////////////////////////////////// - /// Double and triple product calculations /////////////// - //////////////////////////////////////////////////////////////////////////////////// + // ---------------------------------------------------------------------------------- + // Double and triple product calculations + // ---------------------------------------------------------------------------------- /** * Pre-calculates all double products for this triangle, and stores @@ -103,15 +100,6 @@ namespace INTERP_UTILS const TetraCorner minCorner = distances.begin()->second; resetDoubleProducts(seg, minCorner); - // test : reset also for second corner if distance is small enough -#ifdef SECOND_CORNER_RESET - std::map::const_iterator iter = distances.begin(); - ++iter; - if(iter->first < 1.0e-12) - { - resetDoubleProducts(seg, iter->second); - } -#endif distances.clear(); } @@ -133,11 +121,8 @@ namespace INTERP_UTILS const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2]; const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; -#ifdef FIXED_DELTA - const double delta = FIXED_DELTA; -#else + const long double delta = MULT_PREC_F * ( std::abs(term1) + std::abs(term2) ); -#endif if( epsilonEqual(_doubleProducts[8*seg + dp], 0.0, THRESHOLD_F * delta)) { diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index f8793e77b..24ac03432 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -7,17 +7,22 @@ #include #include +/// Precision used for tests of 3D part of INTERP_KERNEL #define VOL_PREC 1.0e-6 + +/// Default relative tolerance in epsilonEqualRelative #define DEFAULT_REL_TOL 1.0e-6 + +/// Default absolute tolerance in epsilonEqual and epsilonEqualRelative #define DEFAULT_ABS_TOL 5.0e-12 namespace INTERP_UTILS { - /////////////////////////////////////////////////////////////////////// - // Math operations for vectors represented by double[3] - arrays //// - /////////////////////////////////////////////////////////////////////// + // ------------------------------------------------------------------- + // Math operations for vectors represented by double[3] - arrays + // ------------------------------------------------------------------- - /* + /** * Copies a double[3] vector from src to dest * * @param src source vector @@ -32,7 +37,7 @@ namespace INTERP_UTILS } } - /* + /** * Creates a string representation of a double[3] vector * * @param pt a 3-vector @@ -45,7 +50,7 @@ namespace INTERP_UTILS return ss.str(); } - /* + /** * Calculates the cross product of two double[3] - vectors. * * @param v1 vector v1 @@ -59,7 +64,7 @@ namespace INTERP_UTILS res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } - /* + /** * Calculates the dot product of two double[3] - vectors * * @param v1 vector v1 @@ -71,7 +76,7 @@ namespace INTERP_UTILS return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]; } - /* + /** * Calculates norm of a double[3] vector * * @param v a vector v @@ -82,7 +87,7 @@ namespace INTERP_UTILS return sqrt(dot(v,v)); } - /* + /** * Compares doubles using an absolute tolerance * This is suitable mainly for comparisons with 0.0 * @@ -97,7 +102,7 @@ namespace INTERP_UTILS // return std::abs(x - y) < errTol; } - /* + /** * Compares doubles using a relative tolerance * This is suitable mainly for comparing larger values to each other. Before performing the relative test, * an absolute test is performed to guard from problems when comparing to 0.0 -- 2.39.2