From 3d0c29faabb6e3fd5c667ab7ccd34db65d7071f5 Mon Sep 17 00:00:00 2001 From: vbd Date: Mon, 3 Sep 2007 11:58:07 +0000 Subject: [PATCH] staffan : fixed bug with OPTIMIZE extended documentation added test per column / row for reflexive mesh tests made more methods inline --- src/INTERP_KERNEL/BoundingBox.cxx | 6 +- src/INTERP_KERNEL/Interpolation3D.cxx | 13 +- src/INTERP_KERNEL/Intersector3D.cxx | 104 ++++--- src/INTERP_KERNEL/Intersector3D.hxx | 53 ++-- src/INTERP_KERNEL/Makefile.in | 23 +- src/INTERP_KERNEL/MeshElement.cxx | 2 +- .../Test/Interpolation3DTest.cxx | 85 +++++- .../Test/Interpolation3DTest.hxx | 30 +- src/INTERP_KERNEL/Test/Makefile.in | 4 +- src/INTERP_KERNEL/Test/TestInterpKernel.cxx | 4 +- .../Test/TransformedTriangleIntersectTest.cxx | 9 + src/INTERP_KERNEL/TransformedTriangle.cxx | 108 ++----- src/INTERP_KERNEL/TransformedTriangle.hxx | 25 +- .../TransformedTriangle_inline.hxx | 281 +++++++++++------- .../TransformedTriangle_intersect.cxx | 160 +++++----- .../TransformedTriangle_math.cxx | 48 ++- src/INTERP_KERNEL/VectorUtils.hxx | 4 +- 17 files changed, 564 insertions(+), 395 deletions(-) diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index 2a4e8d921..d5f3feabb 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -74,7 +74,7 @@ namespace INTERP_UTILS * Determines if the intersection with a given box is empty * * @param box BoundingBox with which intersection is tested - * @returns true if intersection between boxes is empty, false if not + * @return true if intersection between boxes is empty, false if not */ bool BoundingBox::isDisjointWith(const BoundingBox& box) const { @@ -119,7 +119,7 @@ namespace INTERP_UTILS * Gets a coordinate of the box * * @param coord coordinate to set - * @returns value of coordinate + * @return value of coordinate * */ double BoundingBox::getCoordinate(const BoxCoord coord) const @@ -164,7 +164,7 @@ namespace INTERP_UTILS * Checks if the box is valid, which it is if its minimum coordinates are * smaller than its maximum coordinates in all directions. * - * @returns true if the box is valid, false if not + * @return true if the box is valid, false if not */ bool BoundingBox::isValid() const { diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index 318d58c29..71572583e 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -47,7 +47,7 @@ namespace MEDMEM * * @param srcMesh 3-dimensional source mesh * @param targetMesh 3-dimesional target mesh, containing only tetraedra - * @returns vector containing for each element i of the source mesh, a map giving for each element j + * @return vector containing for each element i of the source mesh, a map giving for each element j * of the target mesh which i intersects, the volume of the intersection */ IntersectionMatrix Interpolation3D::interpol_maillages(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) @@ -185,7 +185,7 @@ namespace MEDMEM const int targetIdx = (*iter)->getIndex(); const double vol = intersector->intersectCells(srcIdx, targetIdx); //if(!epsilonEqual(vol, 0.0)) - if(vol != 0.0) + //if(vol != 0.0) { volumes->insert(make_pair(targetIdx, vol)); LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]); @@ -311,20 +311,22 @@ namespace MEDMEM for(vector::const_iterator iter = intersectElems.begin() ; iter != intersectElems.end() ; ++iter) { -#if 0 + const int srcIdx = *iter + 1; const double vol = intersector->intersectCells(srcIdx, targetIdx); - if(!epsilonEqual(vol, 0.0)) + // if(!epsilonEqual(vol, 0.0)) { matrix[srcIdx - 1].insert(make_pair(targetIdx, vol)); } -#endif + } } #endif + + // free allocated memory for(int i = 0 ; i < numSrcElems ; ++i) { @@ -335,6 +337,7 @@ namespace MEDMEM delete targetElems[i]; } + std::cout << "Intersector filtered out " << intersector->filtered << " elements" << std::endl; delete intersector; } diff --git a/src/INTERP_KERNEL/Intersector3D.cxx b/src/INTERP_KERNEL/Intersector3D.cxx index a41ed2ac7..9fd419442 100644 --- a/src/INTERP_KERNEL/Intersector3D.cxx +++ b/src/INTERP_KERNEL/Intersector3D.cxx @@ -23,7 +23,7 @@ namespace MEDMEM * */ Intersector3D::Intersector3D(const MESH& srcMesh, const MESH& targetMesh) - : _srcMesh(srcMesh), _targetMesh(targetMesh) + : _srcMesh(srcMesh), _targetMesh(targetMesh), filtered(0) { } @@ -97,6 +97,13 @@ namespace MEDMEM volume += iter->calculateIntersectionVolume(); } + // reset if it is very small to keep the matrix sparse + // is this a good idea? + if(epsilonEqual(volume, 0.0, 1.0e-11)) + { + volume = 0.0; + } + LOG(2, "Volume = " << volume << ", det= " << T.determinant()); // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation @@ -112,10 +119,10 @@ namespace MEDMEM * @param triangles vector in which the transformed triangles are stored * @param T affine transform that is applied to the nodes of the triangles */ -#ifdef OPTIMIZE_ // optimized version +#ifdef OPTIMIZE_FILTER // optimized version void Intersector3D::triangulate(const medGeometryElement type, const int element, std::vector& triangles, const TetraAffineTransform& T) const { - + // get cell model for the element CELLMODEL cellModel(type); @@ -128,56 +135,61 @@ namespace MEDMEM triangles.reserve(2 * cellModel.getNumberOfConstituents(1)); // loop over faces - const int numFaces = cellModel.getNumberOfConstituents(1); - double* transformedNodeList[numFaces]; - medGeometryElement faceTypes[numFaces]; - bool isOutsideTetra = true; + double** transformedNodes = new double*[cellModel.getNumberOfNodes()]; + bool isOutside[8] = {true, true, true, true, true, true, true, true}; + bool isTargetOutside = false; - for(int i = 1 ; i <= numFaces ; ++i) + // calculate the coordinates of the nodes + for(int i = 1; i <= cellModel.getNumberOfNodes() ; ++i) { - faceTypes[i - 1] = cellModel.getConstituentType(1, i); - CELLMODEL faceModel(faceTypes[i-1]); + const double* node = getCoordsOfNode(i, element, _srcMesh); + transformedNodes[i-1] = new double[3]; + assert(transformedNodes[i-1] != 0); - assert(faceModel.getDimension() == 2); - // assert(faceModel.getNumberOfNodes() == 3); - - transformedNodeList[i - 1] = new double[3 * faceModel.getNumberOfNodes()]; - - // loop over nodes of face - for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) - { - // offset of node from cellIdx - int localNodeNumber = cellModel.getNodeConstituent(1, i, j); - - assert(localNodeNumber >= 1); - assert(localNodeNumber <= cellModel.getNumberOfNodes()); - - const double* node = getCoordsOfNode(localNodeNumber, element, _srcMesh); - - // transform - //{ not totally efficient since we transform each node once per face - T.apply(&transformedNodeList[i - 1][3*(j-1)], node); - //T.apply(&transformedNodes[3*(j-1)], node); + T.apply(transformedNodes[i - 1], node); - LOG(4, "Node " << localNodeNumber << " = " << vToStr(node) << " transformed to " - << vToStr(&transformedNodeList[i - 1][3*(j-1)])); + // LOG(4, "Node " << i << " = " << vToStr(node) << " transformed to " << vToStr(&transformedNodes[i - 1]); + checkIsOutside(transformedNodes[i - 1], isOutside); + } + // std:cout << "here1" << std::endl; + // check if we need to calculate intersection volume + for(int i = 0; i < 8; ++i) + { + if(isOutside[i]) + { + isTargetOutside = true; } - - isOutsideTetra = isOutsideTetra && isTriangleOutsideTetra(transformedNodeList[i - 1]); } + // std:cout << "here2" << std::endl; - if(!isOutsideTetra) + if(!isTargetOutside) { - for(int i = 1 ; i <= numFaces ; ++i) + for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i) { + const medGeometryElement faceType = cellModel.getConstituentType(1, i); + CELLMODEL faceModel(faceType); + + assert(faceModel.getDimension() == 2); + int nodes[faceModel.getNumberOfNodes()]; + + // get the nodes of the face + for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) + { + nodes[j-1] = cellModel.getNodeConstituent(1, i, j) - 1; + + assert(nodes[j-1] +1 >= 0); + assert(nodes[j-1] +1 <= cellModel.getNumberOfNodes()); + } + // std::cout << "here3 : " << nodes[0] << "," << nodes[1] << ", " << nodes[2] << std::endl; // create transformed triangles from face - switch(faceTypes[i - 1]) + switch(faceType) { case MED_TRIA3: // simple take the triangle as it is - triangles.push_back(TransformedTriangle(&transformedNodeList[i - 1][0], &transformedNodeList[i - 1][3], &transformedNodeList[i - 1][6])); + //std::cout << "here3.1 : " << transformedNodes[nodes[0]] << "," << transformedNodes[nodes[1]] << ", " << transformedNodes[nodes[2]] << std::endl; + triangles.push_back(TransformedTriangle(transformedNodes[nodes[0]], transformedNodes[nodes[1]], transformedNodes[nodes[2]])); break; case MED_QUAD4: @@ -195,10 +207,10 @@ namespace MEDMEM //? not sure if this always works // local nodes 1, 2, 3 - triangles.push_back(TransformedTriangle(&transformedNodeList[i - 1][0], &transformedNodeList[i - 1][3], &transformedNodeList[i - 1][6])); + triangles.push_back(TransformedTriangle(transformedNodes[nodes[0]], transformedNodes[nodes[1]], transformedNodes[nodes[2]])); // local nodes 1, 3, 4 - triangles.push_back(TransformedTriangle(&transformedNodeList[i - 1][0], &transformedNodeList[i - 1][6], &transformedNodeList[i - 1][9])); + triangles.push_back(TransformedTriangle(transformedNodes[nodes[0]], transformedNodes[nodes[2]], transformedNodes[nodes[3]])); break; @@ -206,12 +218,20 @@ namespace MEDMEM std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl; assert(false); } + // std:cout << "here4" << std::endl; } } - for(int i = 1 ; i <= numFaces ; ++i) + else + { + ++filtered; + } + + for(int i = 0 ; i < cellModel.getNumberOfNodes() ; ++i) { - delete[] transformedNodeList[i-1]; + delete[] transformedNodes[i]; } + delete[] transformedNodes; + // std:cout << "here5" << std::endl; } #else // un-optimized version diff --git a/src/INTERP_KERNEL/Intersector3D.hxx b/src/INTERP_KERNEL/Intersector3D.hxx index 8642d2cf8..d9b52d37f 100644 --- a/src/INTERP_KERNEL/Intersector3D.hxx +++ b/src/INTERP_KERNEL/Intersector3D.hxx @@ -30,37 +30,50 @@ namespace MEDMEM virtual double intersectCells(int srcCell, int targetCell); - + mutable int filtered; private: - inline bool isTriangleOutsideTetra(double* points) const; + //inline bool isTriangleOutsideTetra(double* points) const; + inline void checkIsOutside(const double* pt, bool* isOutside) const; void triangulate(const MED_EN::medGeometryElement type, const int element, std::vector& triangles, const INTERP_UTILS::TetraAffineTransform& T) const; const MESH& _srcMesh; const MESH& _targetMesh; - + }; - inline bool Intersector3D::isTriangleOutsideTetra(double* points) const +// inline bool Intersector3D::isTriangleOutsideTetra(double* points) const +// { +// const bool leftX = points[0] <= 0.0 && points[3] <= 0.0 && points[6] <= 0.0; +// const bool rightX = points[0] >= 1.0 && points[3] >= 1.0 && points[6] >= 1.0; +// const bool leftY = points[1] <= 0.0 && points[4] <= 0.0 && points[7] <= 0.0; +// const bool rightY = points[1] >= 1.0 && points[4] >= 1.0 && points[7] >= 1.0; +// const bool leftZ = points[2] <= 0.0 && points[5] <= 0.0 && points[8] <= 0.0; +// const bool rightZ = points[2] >= 1.0 && points[5] >= 1.0 && points[8] >= 1.0; +// const double h[3] = +// { +// 1 - points[0] - points[1] - points[2], +// 1 - points[3] - points[4] - points[5], +// 1 - points[6] - points[7] - points[8] +// }; +// const bool leftH = h[0] <= 0.0 && h[1] <= 0.0 && h[2] <= 0.0; +// const bool rightH = h[0] >= 1.0 && h[1] >= 1.0 && h[2] >= 1.0; + +// return leftX || rightX || leftY || rightY || leftZ || rightZ || leftH || rightH; +// } + + inline void Intersector3D::checkIsOutside(const double* pt, bool* isOutside) const { - const bool leftX = points[0] <= 0.0 && points[3] <= 0.0 && points[6] <= 0.0; - const bool rightX = points[0] >= 1.0 && points[3] >= 1.0 && points[6] >= 1.0; - const bool leftY = points[1] <= 0.0 && points[4] <= 0.0 && points[7] <= 0.0; - const bool rightY = points[1] >= 1.0 && points[4] >= 1.0 && points[7] >= 1.0; - const bool leftZ = points[2] <= 0.0 && points[5] <= 0.0 && points[8] <= 0.0; - const bool rightZ = points[2] >= 1.0 && points[5] >= 1.0 && points[8] >= 1.0; - const double h[3] = - { - 1 - points[0] - points[1] - points[2], - 1 - points[3] - points[4] - points[5], - 1 - points[6] - points[7] - points[8] - }; - const bool leftH = h[0] <= 0.0 && h[1] <= 0.0 && h[2] <= 0.0; - const bool rightH = h[0] >= 1.0 && h[1] >= 1.0 && h[2] >= 1.0; - - return leftX || rightX || leftY || rightY || leftZ || rightZ || leftH || rightH; + isOutside[0] = isOutside[0] && (pt[0] <= 0.0); + isOutside[1] = isOutside[1] && (pt[0] >= 1.0); + isOutside[2] = isOutside[2] && (pt[1] <= 0.0); + isOutside[3] = isOutside[3] && (pt[1] >= 1.0); + isOutside[4] = isOutside[4] && (pt[2] <= 0.0); + isOutside[5] = isOutside[5] && (pt[2] >= 1.0); + isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] <= 0.0); + isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0); } }; diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index dab294301..306be5718 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -53,7 +53,9 @@ InterpolationUtils.hxx\ BBTree.H\ MeshUtils.hxx\ Intersector3D.hxx\ -Log.hxx +Log.hxx\ +TransformedTriangle_inline.hxx + # Libraries targets @@ -70,9 +72,8 @@ BoundingBox.cxx\ TetraAffineTransform.cxx\ Interpolation3D.cxx\ Intersector3D.cxx\ -#Log.cxx -#Interpolation3DSurf.cxx\ -#Interpolation2D.cxx\ +Interpolation3DSurf.cxx\ +Interpolation2D.cxx # Executables targets BIN = @@ -80,7 +81,7 @@ BIN_SRC = BIN_SERVER_IDL = BIN_CLIENT_IDL = -TEST_PROGS = #make_plane test_INTERPOL_2D +TEST_PROGS = testUnitTetra LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome @@ -91,8 +92,8 @@ CPPFLAGS+=$(BOOST_CPPFLAGS) # optimization #CXXFLAGS+= -O1 -DOPTIMIZE #CPPFLAGS+= -O1 -DOPTIMIZE -CXXFLAGS+= -DOPTIMIZE -O2 -CPPFLAGS+= -DOPTIMIZE -O2 +CXXFLAGS+= -O2 -DOPTIMIZE_FILTER -DOPTIMIZE +CPPFLAGS+= -O2 -DOPTIMIZE_FILTER -DOPTIMIZE # for log CXXFLAGS+= -DLOG_LEVEL=0 @@ -103,8 +104,8 @@ CPPFLAGS+= -DLOG_LEVEL=0 #CPPFLAGS+=-fprofile-arcs -ftest-coverage #for gprof -CXXFLAGS+=-pg -CPPFLAGS+=-pg +#CXXFLAGS+=-pg +#CPPFLAGS+=-pg #LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. @@ -112,7 +113,7 @@ LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil #LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lutil # for gcov #LDFLAGS+= -lgcov -LDFLAGS+= -pg +#LDFLAGS+= -pg #LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. @@ -127,7 +128,7 @@ ifeq ($(MED_WITH_KERNEL),yes) LDFLAGSFORBIN+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace -lSALOMEBasics endif -LDFLAGSFORBIN+= -pg +#LDFLAGSFORBIN+= -pg LIBSFORBIN=$(BOOSTLIBS) $(MPI_LIBS) diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index 3069662c7..a12c5365d 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -95,7 +95,7 @@ namespace INTERP_UTILS /* * Comparison operator based on the bounding boxes of the elements * - * @returns true if the coordinate _coord of the bounding box of elem1 is + * @return true if the coordinate _coord of the bounding box of elem1 is * strictly smaller than that of the bounding box of elem2 */ bool ElementBBoxOrder::operator()( MeshElement* elem1, MeshElement* elem2) diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx index cc365c224..cdb5dcc23 100644 --- a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx +++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx @@ -9,6 +9,9 @@ #include "VectorUtils.hxx" +#include "MEDMEM_Field.hxx" +#include "MEDMEM_Support.hxx" + // levels : // 1 - titles and volume results // 2 - symmetry / diagonal results and intersection matrix output @@ -23,6 +26,41 @@ using namespace MEDMEM; using namespace std; using namespace INTERP_UTILS; +using namespace MED_EN; + +double Interpolation3DTest::sumRow(const IntersectionMatrix& m, int i) const +{ + double vol = 0.0; + for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter) + { + if(iter->count(i) != 0.0) + { + map::const_iterator iter2 = iter->find(i); + vol += iter2->second; + } + } + return vol; +} + +double Interpolation3DTest::sumCol(const IntersectionMatrix& m, int i) const +{ + double vol = 0.0; + const std::map& col = m[i]; + for(map::const_iterator iter = col.begin() ; iter != col.end() ; ++iter) + { + vol += std::abs(iter->second); + } + return vol; +} + + +void Interpolation3DTest::getVolumes(MESH& mesh, const double*& tab) const +{ + SUPPORT *sup=new SUPPORT(&mesh,"dummy",MED_CELL); + FIELD* f=mesh.getVolume(sup); + tab = f->getValue(); + // delete sup; +} double Interpolation3DTest::sumVolume(const IntersectionMatrix& m) const { @@ -45,6 +83,43 @@ double Interpolation3DTest::sumVolume(const IntersectionMatrix& m) const return vol; } +bool Interpolation3DTest::testVolumes(const IntersectionMatrix& m, MESH& sMesh, MESH& tMesh) const +{ + bool ok = true; + + // source elements + const double* sVol = new double[sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)]; + getVolumes(sMesh, sVol); + + for(int i = 0; i < sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i) + { + const double sum_row = sumRow(m, i+1); + if(!epsilonEqualRelative(sum_row, sVol[i], VOL_PREC)) + { + LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row ); + ok = false; + } + LOG(1, "diff = " < intersecting src = " << mesh1 << ", target = " << mesh2 ); LOG(5, "Loading " << mesh1 << " from " << mesh1path); - const MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1); - + MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1); + LOG(5, "Loading " << mesh2 << " from " << mesh2path); - const MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2); + MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2); m = interpolator->interpol_maillages(sMesh, tMesh); + testVolumes(m, sMesh, tMesh); + LOG(1, "Intersection calculation done. " << std::endl ); } @@ -213,7 +290,7 @@ void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mes { LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol ); CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1)); - + if(isTestReflexive) { CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1)); diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx index 91a76c6b7..767488265 100644 --- a/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx +++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx @@ -7,15 +7,17 @@ #define ERR_TOL 1.0e-8 using MEDMEM::Interpolation3D; +class MEDMEM::MESH; class Interpolation3DTest : public CppUnit::TestFixture { // single - element CPPUNIT_TEST_SUITE( Interpolation3DTest ); - +#if 0 CPPUNIT_TEST( tetraReflexiveUnit ); CPPUNIT_TEST( tetraReflexiveGeneral ); + CPPUNIT_TEST( tetraNudgedSimpler ); CPPUNIT_TEST( tetraNudged ); CPPUNIT_TEST( tetraCorner ); @@ -34,21 +36,24 @@ class Interpolation3DTest : public CppUnit::TestFixture // multi - element CPPUNIT_TEST( tetraComplexIncluded ); +#endif CPPUNIT_TEST( dividedUnitTetraSimplerReflexive ); CPPUNIT_TEST( dividedUnitTetraReflexive ); - CPPUNIT_TEST( nudgedDividedUnitTetra ); - CPPUNIT_TEST( nudgedDividedUnitTetraSimpler ); - CPPUNIT_TEST( dividedGenTetra ); + //#if 0 + //CPPUNIT_TEST( nudgedDividedUnitTetra ); + //CPPUNIT_TEST( nudgedDividedUnitTetraSimpler ); + //CPPUNIT_TEST( dividedGenTetra ); CPPUNIT_TEST( boxReflexive ); CPPUNIT_TEST( boxReflexiveModerate ); - CPPUNIT_TEST( tetraBoxes ); - CPPUNIT_TEST( moderateBoxes ); - CPPUNIT_TEST( moderateBoxesSmaller ); + //CPPUNIT_TEST( tetraBoxes ); + //CPPUNIT_TEST( moderateBoxes ); + //CPPUNIT_TEST( moderateBoxesSmaller ); CPPUNIT_TEST( moderateBoxSmallReflexive ); CPPUNIT_TEST( moderateBoxEvenSmallerReflexive ); CPPUNIT_TEST( tinyBoxReflexive ); - CPPUNIT_TEST( simpleHexaBox ); + //CPPUNIT_TEST( simpleHexaBox ); + //#endif CPPUNIT_TEST_SUITE_END(); @@ -220,6 +225,14 @@ private: Interpolation3D* interpolator; + double sumRow(const IntersectionMatrix& m, int i) const; + + double sumCol(const IntersectionMatrix& m, int i) const; + + void getVolumes( MEDMEM::MESH& mesh,const double*& tab) const; + + bool testVolumes(const IntersectionMatrix& m, MEDMEM::MESH& sMesh, MEDMEM::MESH& tMesh) const; + double sumVolume(const IntersectionMatrix& m) const; bool areCompatitable( const IntersectionMatrix& m1, const IntersectionMatrix& m2) const; @@ -235,6 +248,7 @@ private: void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const; + }; #endif diff --git a/src/INTERP_KERNEL/Test/Makefile.in b/src/INTERP_KERNEL/Test/Makefile.in index f9d22c043..94e48888a 100644 --- a/src/INTERP_KERNEL/Test/Makefile.in +++ b/src/INTERP_KERNEL/Test/Makefile.in @@ -64,8 +64,8 @@ CPPFLAGS+=$(BOOST_CPPFLAGS) $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/sr #CPPFLAGS+= -I/usr/include/cppunit # for log -CXXFLAGS+= -DLOG_LEVEL=1 -CPPFLAGS+= -DLOG_LEVEL=1 +CXXFLAGS+= -DLOG_LEVEL=0 +CPPFLAGS+= -DLOG_LEVEL=0 # for gcov #CXXFLAGS+=-fprofile-arcs -ftest-coverage diff --git a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx index 39793a3e3..1d3948c66 100644 --- a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx +++ b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx @@ -25,8 +25,8 @@ // --- Registers the fixture into the 'registry' //CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTestMultiElement ); -//CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTest ); -CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleIntersectTest ); +CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTest ); +//CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleIntersectTest ); //CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleTest ); //CPPUNIT_TEST_SUITE_REGISTRATION( TestBogusClass ); diff --git a/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx b/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx index e3b4e7ef3..7ac217231 100644 --- a/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx +++ b/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx @@ -17,6 +17,15 @@ // all these cases is therefore of low priority. //////////////////////////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////////////////////////// +/// ! NB ! +/// These tests will not pass if the library has been compiled with the -DOPTIMIZED flag. This is because +/// the tests for zero double products in testSegmentEdgeIntersection, testSegmentCornerIntersection and +/// testSegmentRayIntersection have in these cases been moved to the calling method +/// TransformedTriangle::calculateIntersectionPolygons(). This set of tests is meant to assure the correct +/// functioning of the lower level, and thus calls these methods directly. +//////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////// // Intersections tested (number indicates first triangle which contains the intersection): // ----------------------------------------------------------------------------------------------------- diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 3cea450fd..6986e509d 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -10,9 +10,6 @@ #include "VectorUtils.hxx" -#undef MERGE_CALC -#undef COORDINATE_CORRECTION 1.0e-15 - class ProjectedCentralCircularSortOrder @@ -62,7 +59,6 @@ namespace INTERP_UTILS * Constructor * * The coordinates are copied to the internal member variables - * (slower but safer - investigate if this is necessary) * * @param p array of three doubles containing coordinates of P * @param q array of three doubles containing coordinates of Q @@ -91,15 +87,6 @@ namespace INTERP_UTILS _coords[5*Q + 4] = 1 - q[0] - q[1]; _coords[5*R + 4] = 1 - r[0] - r[1]; -#ifdef COORDINATE_CORRECTION - - // correction of small (imprecise) coordinate values - for(int i = 0 ; i < 15 ; ++i) - { - _coords[i] = epsilonEqual(_coords[i], 0.0, COORDINATE_CORRECTION) ? 0.0 : _coords[i]; - } -#endif - // initialise rest of data preCalculateDoubleProducts(); @@ -111,6 +98,12 @@ namespace INTERP_UTILS } + /* + * Destructor + * + * Deallocates the memory used to store the points of the polygons. + * This memory is allocated in calculateIntersectionPolygons(). + */ TransformedTriangle::~TransformedTriangle() { // delete elements of polygons @@ -128,7 +121,7 @@ namespace INTERP_UTILS * Calculates the volume of intersection between the triangle and the * unit tetrahedron. * - * @returns volume of intersection of this triangle with unit tetrahedron, + * @return volume of intersection of this triangle with unit tetrahedron, * as described in Grandy * */ @@ -159,23 +152,12 @@ namespace INTERP_UTILS } - // normalize + // normalize sign sign = sign > 0.0 ? 1.0 : -1.0; LOG(2, "-- Calculating intersection polygons ... "); calculateIntersectionPolygons(); -#ifdef MERGE_CALC - const bool mergeCalculation = isPolygonAOnHFacet(); - if(mergeCalculation) - { - // move points in B to A to avoid missing points - // NB : need to remove elements from B in order to handle deletion properly - _polygonA.insert(_polygonA.end(), _polygonB.begin(), _polygonB.end()); - _polygonB.clear(); - } -#endif - double barycenter[3]; // calculate volume under A @@ -191,11 +173,7 @@ namespace INTERP_UTILS double volB = 0.0; // if triangle is not in h = 0 plane, calculate volume under B -#ifdef MERGE_CALC - if((!mergeCalculation) && _polygonB.size() > 2) -#else - if(!isTriangleInPlaneOfFacet(XYZ) && _polygonB.size() > 2) -#endif + if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) { LOG(2, "---- Treating polygon B ... "); @@ -206,7 +184,7 @@ namespace INTERP_UTILS } LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********"); - + return sign * (volA + volB); } @@ -423,7 +401,7 @@ namespace INTERP_UTILS // segment - ray for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) { - if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*corner]] && testSegmentRayIntersection(seg, corner)) + if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner)) { double* ptB = new double[3]; copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); @@ -588,6 +566,10 @@ namespace INTERP_UTILS 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. + // is triangle parallel to x == 0 ? if(isTriangleInPlaneOfFacet(OZX)) { @@ -606,8 +588,6 @@ namespace INTERP_UTILS // NB : do not change place of first object, with respect to which the order // is defined sort((polygon.begin()), polygon.end(), order); - //stable_sort((++polygon.begin()), polygon.end(), order); - LOG(3,"Sorted polygon is "); for(int i = 0 ; i < polygon.size() ; ++i) @@ -625,7 +605,7 @@ namespace INTERP_UTILS * * @param poly one of the two intersection polygons * @param barycenter array of three doubles with the coordinates of the barycenter - * @returns the volume between the polygon and the z = 0 plane + * @return the volume between the polygon and the z = 0 plane * */ double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter) @@ -664,7 +644,7 @@ namespace INTERP_UTILS * Checks if the triangle lies in the plane of a given facet * * @param facet one of the facets of the tetrahedron - * @returns true if PQR lies in the plane of the facet, false if not + * @return true if PQR lies in the plane of the facet, false if not */ bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) { @@ -674,60 +654,26 @@ namespace INTERP_UTILS for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) { - // ? should have epsilon-equality here? -#ifdef EPS_TESTING - if(!epsilonEqualRelative(_coords[5*c + coord], 0.0, TEST_EPS * _coords[5*c + coord])) -#else - if(_coords[5*c + coord] != 0.0) -#endif - { - return false; - } - } - - return true; - } - - bool TransformedTriangle::isPolygonAOnHFacet() const - { - // need to have vector of unique points in order to determine the "real" number of - // of points in the polygon, to avoid problems when polygon A has less than 3 points - - using ::Vector3Cmp; - std::vector pAUnique; - pAUnique.reserve(_polygonA.size()); - Vector3Cmp cmp; - unique_copy(_polygonA.begin(), _polygonA.end(), back_inserter(pAUnique), cmp); - //for(std::vector::const_iterator iter = pAUnique.begin() ; iter != pAUnique.end() ; ++iter) - //std::cout << "next : " << vToStr(*iter) << std::endl; - - // std::cout << "paunique size = " << pAUnique.size() << std::endl; - if(pAUnique.size() < 3) - { - return false; - } - for(std::vector::const_iterator iter = _polygonA.begin() ; iter != _polygonA.end() ; ++iter) - { - const double* pt = *iter; - const double h = 1.0 - pt[0] - pt[1] - pt[2]; - //// std::cout << "h = " << h << std::endl; - - //if(h != 0.0) - if(!epsilonEqual(h, 0.0)) + if(_coords[5*c + coord] != 0.0) { return false; } } - // std::cout << "Polygon A is on h = 0 facet" << std::endl; + return true; } + /* + * 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() { for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) { // check z-coords for all points - if(_coords[5*c + 2] > 0.0) + if(_coords[5*c + 2] >= 0.0) { return false; } @@ -735,6 +681,10 @@ namespace INTERP_UTILS return true; } + /* + * Prints the coordinates of the triangle to std::cout + * + */ void TransformedTriangle::dumpCoords() { std::cout << "Coords : "; diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index f506bcc1d..97033e7a5 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -3,12 +3,6 @@ #include -#undef EPS_TESTING // does not give good results - -#ifdef EPS_TESTING -#define TEST_EPS 1.0e-14 -#endif - #ifdef OPTIMIZE #define OPT_INLINE inline #else @@ -60,7 +54,7 @@ namespace INTERP_UTILS * calculateIntersectionPolygons() : * This method goes through all the possible ways in which the triangle can intersect the tetrahedron and tests for these * types of intersections in accordance with the formulas described in Grandy. These tests are implemented in the test* - methods. - * The formulas in the article are stated for one case each only, while the calculation must take into account all cases. + * The formulas in the article are stated for one case each only, while the calculation must take into account all cases. * To this end, a number of tables, implemented as static const arrays of different types, are used. The tables * mainly contain values of the different enumeration types described at the beginning of the class interface. For example, * the formula Grandy gives for the segment-halfstrip intersection tests ([30]) is for use with the halfstrip above the zx edge. @@ -79,14 +73,13 @@ namespace INTERP_UTILS public: - friend class ::TransformedTriangleTest; friend class ::TransformedTriangleIntersectTest; - /** * Enumerations representing the different geometric elements of the unit tetrahedron - * and the triangle. + * and the triangle. The end element, NO_* gives the number of elements in the enumeration + * and can be used as end element in loops. */ /// Corners of tetrahedron enum TetraCorner { O = 0, X, Y, Z, NO_TET_CORNER }; @@ -117,7 +110,6 @@ namespace INTERP_UTILS void dumpCoords(); - private: //////////////////////////////////////////////////////////////////////////////////// @@ -133,16 +125,13 @@ namespace INTERP_UTILS double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter); //////////////////////////////////////////////////////////////////////////////////// - /// Detection of (very) degenerate cases //////////// + /// Detection of degenerate triangles //////////// //////////////////////////////////////////////////////////////////////////////////// bool isTriangleInPlaneOfFacet(const TetraFacet facet); bool isTriangleBelowTetraeder(); - bool isPolygonAOnHFacet() const; - - //////////////////////////////////////////////////////////////////////////////////// /// Intersection test methods and intersection point calculations //////// //////////////////////////////////////////////////////////////////////////////////// @@ -181,11 +170,11 @@ namespace INTERP_UTILS bool testTriangleSurroundsEdge(const TetraEdge edge) const; - bool testEdgeIntersectsTriangle(const TetraEdge edge) const; + OPT_INLINE bool testEdgeIntersectsTriangle(const TetraEdge edge) const; - bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const; + OPT_INLINE bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const; - bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const; + OPT_INLINE bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const; bool testSegmentIntersectsHPlane(const TriSegment seg) const; diff --git a/src/INTERP_KERNEL/TransformedTriangle_inline.hxx b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx index 9b805b529..63c37d026 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_inline.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx @@ -1,3 +1,25 @@ +#ifndef _TRANSFORMED_TRIANGLE_INLINE_HXX +#define _TRANSFORMED_TRIANGLE_INLINE_HXX + +// This file contains inline versions of some of the methods in the TransformedTriangle*.cxx files. +// 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. +//////////////////////////////////////////////////////////////////////////////////////// +inline void TransformedTriangle::preCalculateTriangleSurroundsEdge() +{ + for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) + { + _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge); + } +} + + +//////////////////////////////////////////////////////////////////////////////////////// +/// TransformedTriangle_math.cxx /////// +//////////////////////////////////////////////////////////////////////////////////////// inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner) { // set the three corresponding double products to 0.0 @@ -17,32 +39,11 @@ inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const }; } -/** - * Returns the stable double product c_{xy}^{pq} - * - * @pre The stable double products have been calculated with preCalculateDoubleProducts. - * @param seg segment of triangle - * @param dp double product sought - * - * @returns stabilised double product c_{xy}^{pq} - * - */ inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const { return _doubleProducts[8*seg + dp]; } -/** - * Returns the stable triple product t_X for a given corner - * The triple product gives the signed volume of the tetrahedron between - * this corner and the triangle PQR. These triple products have been calculated - * in a way to avoid problems with cancellation. - * - * @pre double products have already been calculated - * @pre triple products have already been calculated - * @param corner corner for which the triple product is calculated - * @returns triple product associated with corner (see Grandy, eqs. [50]-[52]) - */ inline double TransformedTriangle::calcStableT(const TetraCorner corner) const { assert(_isTripleProductsCalculated); @@ -50,17 +51,6 @@ inline double TransformedTriangle::calcStableT(const TetraCorner corner) const return _tripleProducts[corner]; } -/** - * Calculates the given double product c_{xy}^{pq} = x_p*y_q - y_p*x_q for a - * a segment PQ of the triangle. This method does not compensate for - * precision errors. - * - * @param seg segment of triangle - * @param dp double product sought - * - * @returns double product c_{xy}^{pq} - * - */ inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const { @@ -76,57 +66,24 @@ 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]; } -inline void TransformedTriangle::preCalculateTriangleSurroundsEdge() -{ - for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) - { - _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge); - } -} - -/** - * Tests if the given edge of the tetrahedron intersects the triangle PQR. (Grandy, eq [17]) - * - * @param edge edge of tetrahedron - * @returns true if PQR intersects the edge, and the edge is not in the plane of the triangle. - */ +//////////////////////////////////////////////////////////////////////////////////////// +/// TransformedTriangle_intersect.cxx /////// +//////////////////////////////////////////////////////////////////////////////////////// inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const { return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge); } -/** - * Tests if the given segment of the triangle intersects the given facet of the tetrahedron. - * (Grandy, eq. [19]) - * - * @param seg segment of the triangle - * @param facet facet of the tetrahedron - * @returns true if the segment intersects the facet - */ inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const { return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet); } -/** - * Tests if the triangle PQR intersects the ray pointing in the upwards z - direction - * from the given corner of the tetrahedron. (Grandy eq. [29]) - * - * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z) - * @returns true if the upwards ray from the corner intersects the triangle. - */ inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const { return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner ); } -/** - * Tests if the given corner of triangle PQR lies in the interior of the unit tetrahedron - * (0 <= x_p, y_p, z_p, h_p <= 1) - * - * @param corner corner of the triangle PQR - * @returns true if the corner lies in the interior of the unit tetrahedron. - */ inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const { const double pt[4] = @@ -144,48 +101,162 @@ inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) return false; } } - return true; + return true; } -/** - * Tests if the given corner of triangle PQR lies on the facet h = 0 (the XYZ facet) - * (0 <= x_p, y_p, z_p <= 1 && h_p = 0) - * - * @param corner corner of the triangle PQR - * @returns true if the corner lies on the facet h = 0 - */ inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const - { - const double pt[4] = - { - _coords[5*corner], // x - _coords[5*corner + 1], // y - _coords[5*corner + 2], // z - _coords[5*corner + 3] // h - }; +{ +#if 0 + const double pt[4] = + { + _coords[5*corner], // x + _coords[5*corner + 1], // y + _coords[5*corner + 2], // z + _coords[5*corner + 3] // h + }; +#endif + const double* pt = &_coords[5*corner]; - if(pt[3] != 0.0) - { - return false; - } - - for(int i = 0 ; i < 3 ; ++i) - { - if(pt[i] < 0.0 || pt[i] > 1.0) - { - return false; - } - } - return true; - } + if(pt[3] != 0.0) + { + return false; + } + + for(int i = 0 ; i < 3 ; ++i) + { + if(pt[i] < 0.0 || pt[i] > 1.0) + { + return false; + } + } + return true; +} inline bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const - { - const double x = _coords[5*corner]; - const double y = _coords[5*corner + 1]; - const double h = _coords[5*corner + 3]; - const double H = _coords[5*corner + 4]; +{ + const double x = _coords[5*corner]; + const double y = _coords[5*corner + 1]; + const double h = _coords[5*corner + 3]; + const double H = _coords[5*corner + 4]; - return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0; + return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0; - } +} + +inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const +{ + + assert(edge < H01); + + // correspondance edge - triple products + // for edges OX, ..., ZX (Grandy, table III) + static const TetraCorner TRIPLE_PRODUCTS[12] = + { + X, O, // OX + Y, O, // OY + Z, O, // OZ + X, Y, // XY + Y, Z, // YZ + Z, X, // ZX + }; + + // Grandy, [16] + const double t1 = calcStableT(TRIPLE_PRODUCTS[2*edge]); + const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]); + + //? should equality with zero use epsilon? + LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 ); + return (t1*t2 <= 0.0) && (t1 - t2 != 0.0); +} + +inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const +{ +#if 0 + const double signs[3] = + { + SIGN_FOR_SEG_FACET_INTERSECTION[3*facet], + SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 1], + SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 2] + }; +#endif + + const double* signs = &SIGN_FOR_SEG_FACET_INTERSECTION[3*facet]; + const double c1 = signs[0]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet]); + const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]); + const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]); + + return (c1*c3 > 0.0) && (c2*c3 > 0.0); +} + +inline bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const +{ + // use correspondance facet a = 0 <=> offset for coordinate a in _coords + // and also correspondance segment AB => corner A + const double coord1 = _coords[5*seg + facet]; + const double coord2 = _coords[5*( (seg + 1) % 3) + facet]; + + //? should we use epsilon-equality here in second test? + LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); + + return (coord1*coord2 <= 0.0) && (coord1 != coord2); +} + +inline bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const +{ + // get the H - coordinates + const double coord1 = _coords[5*seg + 4]; + const double coord2 = _coords[5*( (seg + 1) % 3) + 4]; + //? should we use epsilon-equality here in second test? + LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); + + return (coord1*coord2 <= 0.0) && (coord1 != coord2); +} + +inline bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) const +{ + // ? There seems to be an error in Grandy -> it should be C_XY instead of C_YZ in [28]. + // ? I haven't really figured out why, but it seems to work. + const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY); + + LOG(6, "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) ); + LOG(6, "] - stable : " << calcStableT(corner) ); + + //? we don't care here if the triple product is "invalid", that is, the triangle does not surround one of the + // edges going out from the corner (Grandy [53]) + if(!_validTP[corner]) + { + return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0; + } + else + { + return ( calcStableT(corner) * normal ) >= 0.0; + } +} + +inline bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const +{ + assert(corner == X || corner == Y || corner == Z); + + // double products to use for the possible corners + static const DoubleProduct DP_FOR_RAY_INTERSECTION[4] = + { + DoubleProduct(0), // O - only here to fill out and make indices match + C_10, // X + C_01, // Y + C_XY // Z + }; + + const DoubleProduct dp = DP_FOR_RAY_INTERSECTION[corner]; + + const double cPQ = calcStableC(PQ, dp); + const double cQR = calcStableC(QR, dp); + const double cRP = calcStableC(RP, dp); + + //? NB here we have no correction for precision - is this good? + // Our authority Grandy says nothing + LOG(5, "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" ); + + return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 ); + +} +#endif diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index 6a1aed56e..de37ecc49 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -7,15 +7,11 @@ #define SEG_RAY_TABLE 1 // seems correct - - - - namespace INTERP_UTILS { //////////////////////////////////////////////////////////////////////////////////// - /// Constants ///////////////// + /// Correspondance tables describing all the variations of formulas. ////////////// //////////////////////////////////////////////////////////////////////////////////// // correspondance facet - double product @@ -27,6 +23,24 @@ namespace INTERP_UTILS C_ZH, C_ZX, C_YZ, // OXY C_XH, C_YH, C_ZH // XYZ }; +#if 0 + template + inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection() + { + return NO_DP; + } + + template<> + inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection + { + return C_XH; + } + +#define DEF_DP_FOR_SEG_FACET(FACET,SEG,DP) template<> inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection { return DP; } + DEF_DP_FOR_SEG_FACET() + +#endif + // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] = { @@ -118,24 +132,23 @@ namespace INTERP_UTILS }; #endif - - - + //////////////////////////////////////////////////////////////////////////////////// /// Intersection test methods and intersection point calculations //////// //////////////////////////////////////////////////////////////////////////////////// -#ifndef OPTIMIZE +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the given edge of the tetrahedron intersects the triangle PQR. (Grandy, eq [17]) * * @param edge edge of tetrahedron - * @returns true if PQR intersects the edge, and the edge is not in the plane of the triangle. + * @return true if PQR intersects the edge, and the edge is not in the plane of the triangle. */ bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const { return testTriangleSurroundsEdge(edge) && testEdgeIntersectsTriangle(edge); } #endif + /** * Calculates the point of intersection between the given edge of the tetrahedron and the * triangle PQR. (Grandy, eq [22]) @@ -176,14 +189,14 @@ namespace INTERP_UTILS } } -#ifndef OPTIMIZE +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the given segment of the triangle intersects the given facet of the tetrahedron. * (Grandy, eq. [19]) * * @param seg segment of the triangle * @param facet facet of the tetrahedron - * @returns true if the segment intersects the facet + * @return true if the segment intersects the facet */ bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const { @@ -228,6 +241,7 @@ namespace INTERP_UTILS const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx]; const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx]; pt[i] = -( sign * calcStableC(seg, dp) ) / s; + LOG(4, "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] ); LOG(4, "c(" << seg << ", " << dp << ") = " << sign * calcStableC(seg, dp) ); assert(pt[i] >= 0.0); @@ -239,22 +253,17 @@ namespace INTERP_UTILS /** * Tests if the given segment of the triangle intersects the given edge of the tetrahedron (Grandy, eq. [20] - * + * If the OPTIMIZE is defined, it does not do the test the double product that should be zero. * @param seg segment of the triangle * @param edge edge of tetrahedron - * @returns true if the segment intersects the edge + * @return true if the segment intersects the edge */ bool TransformedTriangle::testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const { #ifndef OPTIMIZE // in this case, we have already checked if the double product is zero -#ifdef EPS_TESTING - if(!epsilonEqual(calcStableC(seg,DoubleProduct( edge )), 0.0, TEST_EPS)) - -#else if(calcStableC(seg,DoubleProduct( edge )) != 0.0) -#endif { return false; } @@ -368,7 +377,9 @@ namespace INTERP_UTILS } // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2) + pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator; + //std::cout << "pt[i] = " << pt[i] << std::endl; assert(pt[i] >= 0.0); // check we are in tetraeder assert(pt[i] <= 1.0); @@ -378,11 +389,11 @@ namespace INTERP_UTILS /** * Tests if the given segment of the triangle intersects the given corner of the tetrahedron. - * (Grandy, eq. [21]) + * (Grandy, eq. [21]). If OPTIMIZE is defined, the double products that should be zero are not verified. * * @param seg segment of the triangle * @param corner corner of the tetrahedron - * @returns true if the segment intersects the corner + * @return true if the segment intersects the corner */ bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const { @@ -404,11 +415,8 @@ namespace INTERP_UTILS const TetraEdge edge = EDGES_FOR_CORNER[3*corner + i]; const DoubleProduct dp = DoubleProduct( edge ); const double c = calcStableC(seg, dp); -#ifdef EPS_TESTING - if(!epsilonEqual(c, 0.0, TEST_EPS)) -#else + if(c != 0.0) -#endif { return false; } @@ -428,13 +436,13 @@ namespace INTERP_UTILS return false; } -#ifndef OPTIMIZE +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the triangle PQR intersects the ray pointing in the upwards z - direction * from the given corner of the tetrahedron. (Grandy eq. [29]) * * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z) - * @returns true if the upwards ray from the corner intersects the triangle. + * @return true if the upwards ray from the corner intersects the triangle. */ bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const { @@ -448,7 +456,7 @@ namespace INTERP_UTILS * * @param seg segment of the triangle * @param edge edge of the h = 0 plane of the tetrahedron (XY, YZ, ZX) - * @returns true if the upwards ray from the corner intersects the triangle. + * @return true if the upwards ray from the corner intersects the triangle. */ bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edge) { @@ -548,10 +556,11 @@ namespace INTERP_UTILS /** * Tests if the given segment of triangle PQR intersects the ray pointing * in the upwards z - direction from the given corner of the tetrahedron. (Grandy eq. [29]) + * If OPTIMIZE is defined, the double product that should be zero is not verified. * * @param seg segment of the triangle PQR * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z) - * @returns true if the upwards ray from the corner intersects the segment. + * @return true if the upwards ray from the corner intersects the segment. */ bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const { @@ -570,18 +579,14 @@ namespace INTERP_UTILS OZX, // Z }; - +#ifndef OPTIMIZE // in this case we have already checked that the double product is zero const DoubleProduct dp0 = DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx]; const double cVal0 = calcStableC(seg, dp0); -#ifndef OPTIMIZE // in this case we have already checked that the double product is zero + //? epsilon-equality here? // cond. 1 -#ifdef EPS_TESTING - if(epsilonEqual(cVal0, 0.0, TEST_EPS)) -#else if(cVal0 != 0.0) -#endif { LOG(4, "SR fails at cond 1 cVal0 = " << cVal0 ); return false; @@ -619,13 +624,13 @@ namespace INTERP_UTILS } -#ifndef OPTIMIZE +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the given corner of triangle PQR lies in the interior of the unit tetrahedron * (0 <= x_p, y_p, z_p, h_p <= 1) * * @param corner corner of the triangle PQR - * @returns true if the corner lies in the interior of the unit tetrahedron. + * @return true if the corner lies in the interior of the unit tetrahedron. */ bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const { @@ -653,7 +658,7 @@ namespace INTERP_UTILS * (0 <= x_p, y_p, z_p <= 1 && h_p = 0) * * @param corner corner of the triangle PQR - * @returns true if the corner lies on the facet h = 0 + * @return true if the corner lies on the facet h = 0 */ bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const { @@ -680,6 +685,12 @@ namespace INTERP_UTILS return true; } + /* + * Tests if the given corner of the triangle lies above the XYZ-facet of the tetrahedron. + * + * @param corner corner of the triangle + * @return true if the corner lies above the XYZ facet, false if not. + */ bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const { const double x = _coords[5*corner]; @@ -688,7 +699,6 @@ namespace INTERP_UTILS const double H = _coords[5*corner + 4]; return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0; - } #endif @@ -701,14 +711,13 @@ namespace INTERP_UTILS * given edge of the tetrahedron lies. * * @param edge edge of tetrahedron - * @returns true if PQR surrounds edge, false if not (see Grandy, eq. [53]) + * @return true if PQR surrounds edge, false if not (see Grandy, eq. [53]) */ bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const { // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III // so we can use the edge directly - // optimization : we could use _doubleProducts directly here const double cPQ = calcStableC(PQ, DoubleProduct(edge)); const double cQR = calcStableC(QR, DoubleProduct(edge)); const double cRP = calcStableC(RP, DoubleProduct(edge)); @@ -717,11 +726,7 @@ namespace INTERP_UTILS // if two or more c-values are zero we disallow x-edge intersection // Grandy, p.446 -#ifdef EPS_TESTING - const int numZeros = (epsilonEqual(cPQ,0.0) ? 1 : 0) + (epsilonEqual(cQR,0.0) ? 1 : 0) + (epsilonEqual(cRP, 0.0) ? 1 : 0); -#else const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0); -#endif if(numZeros >= 2 ) { @@ -731,11 +736,12 @@ namespace INTERP_UTILS return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2; } +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the corners of the given edge lie on different sides of the triangle PQR. * * @param edge edge of the tetrahedron - * @returns true if the corners of the given edge lie on different sides of the triangle PQR + * @return true if the corners of the given edge lie on different sides of the triangle PQR * or if one (but not both) lies in the plane of the triangle. */ bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const @@ -761,20 +767,18 @@ namespace INTERP_UTILS //? should equality with zero use epsilon? LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 ); -#ifdef EPS_TESTING - return (t1*t2 <= 0.0) && (!epsilonEqualRelative(t1, t2, TEST_EPS * std::max(t1, t2))); -#else return (t1*t2 <= 0.0) && (t1 - t2 != 0.0); -#endif } +#endif +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the given facet of the tetrahedron surrounds the axis on which the * given segment of the triangle lies. * * @param seg segment of the triangle * @param facet facet of the tetrahedron - * @returns true if the facet surrounds the segment, false if not (see Grandy, eq. [18]) + * @return true if the facet surrounds the segment, false if not (see Grandy, eq. [18]) */ bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const { @@ -785,30 +789,23 @@ namespace INTERP_UTILS SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 1], SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 2] }; - + const double c1 = signs[0]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet]); const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]); const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]); -#ifdef EPS_TSTING - if(epsilonEqual(c1, 0.0, TEST_EPS) || epsilonEqual(c2, 0.0, TEST_EPS) || epsilonEqual(c3, 0.0, TEST_EPS)) - { - return false; - } - else -#endif - { - return (c1*c3 > 0.0) && (c2*c3 > 0.0); - } + return (c1*c3 > 0.0) && (c2*c3 > 0.0); } +#endif +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the corners of the given segment lie on different sides of the given facet. * * @param seg segment of the triangle * @param facet facet of the tetrahedron - * @returns true if the corners of the given segment lie on different sides of the facet + * @return true if the corners of the given segment lie on different sides of the facet * or if one (but not both) lies in the plane of the facet. (see Grandy, eq. [18]) */ bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const @@ -820,13 +817,20 @@ namespace INTERP_UTILS //? should we use epsilon-equality here in second test? LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); -#ifdef EPS_TESTING - return (coord1*coord2 <= 0.0) && epsilonEqualRelative(coord1,coord2, TEST_EPS * std::max(coord1, coord2)); -#else + return (coord1*coord2 <= 0.0) && (coord1 != coord2); -#endif } +#endif +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx + /* + * Tests if the H-coordinates (1 - x - y) for the two end of a segment of the triangle + * lie on different sides of the H = 0 plane. + * + * @param seg segment of the Triangle + * @return true if the two points of the triangle lie on different sides of the H = 0 plane : + * that is, if their H-coordinates have different signs + */ bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const { // get the H - coordinates @@ -834,32 +838,31 @@ namespace INTERP_UTILS const double coord2 = _coords[5*( (seg + 1) % 3) + 4]; //? should we use epsilon-equality here in second test? LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); -#ifdef EPS_TESTING - return (coord1*coord2 <= 0.0) && epsilonEqualRelative(coord1,coord2, TEST_EPS * std::max(coord1, coord2)); -#else + return (coord1*coord2 <= 0.0) && (coord1 != coord2); -#endif } +#endif +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the triangle PQR lies above a given corner in the z-direction (implying that the * ray pointing upward in the z-direction from the corner can intersect the triangle) * (Grandy eq.[28]) * * @param corner corner of the tetrahedron - * @returns true if the triangle lies above the corner in the z-direction + * @return true if the triangle lies above the corner in the z-direction */ bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) const { - //? is it always YZ here ? - //? changed to XY ! + // ? There seems to be an error in Grandy -> it should be C_XY instead of C_YZ in [28]. + // ? I haven't really figured out why, but it seems to work. const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY); + LOG(6, "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) ); LOG(6, "] - stable : " << calcStableT(corner) ); //? we don't care here if the triple product is "invalid", that is, the triangle does not surround one of the // edges going out from the corner (Grandy [53]) - //return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0; if(!_validTP[corner]) { return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0; @@ -869,13 +872,15 @@ namespace INTERP_UTILS return ( calcStableT(corner) * normal ) >= 0.0; } } +#endif +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx /** * Tests if the triangle PQR surrounds the ray pointing upwards in the z-direction * from the given corner. * * @param corner corner on face XYZ of tetrahedron - * @returns true if PQR surrounds ray, false if not (see Grandy, eq. [18]) + * @return true if PQR surrounds ray, false if not (see Grandy, eq. [18]) */ bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const { @@ -899,8 +904,9 @@ namespace INTERP_UTILS //? NB here we have no correction for precision - is this good? // Our authority Grandy says nothing LOG(5, "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" ); + return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 ); } - +#endif }; // NAMESPACE diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index 49de4a361..3e3bbcf63 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -101,14 +101,13 @@ namespace INTERP_UTILS resetDoubleProducts(seg, iter->second); } #endif + distances.clear(); } - distances.clear(); } // -- (2) check that each double product statisfies Grandy, [47], else set to 0 - for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1)) { for(DoubleProduct dp = C_YZ ; dp <= C_10 ; dp = DoubleProduct(dp + 1)) @@ -132,11 +131,16 @@ namespace INTERP_UTILS if( epsilonEqual(_doubleProducts[8*seg + dp], 0.0, THRESHOLD_F * delta)) { + // debug output +#if LOG_LEVEL >= 5 if(_doubleProducts[8*seg + dp] != 0.0) { LOG(5, "Double product for (seg,dp) = (" << seg << ", " << dp << ") = " ); LOG(5, std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" ); } +#endif + + _doubleProducts[8*seg + dp] = 0.0; } @@ -146,12 +150,18 @@ namespace INTERP_UTILS _isDoubleProductsCalculated = true; } + /* + * Checks if the double products for a given segment are consistent, as defined by + * Grandy, [46] + * + * @param seg Segment for which to check consistency of double products + * @return true if the double products are consistent, false if not + */ bool TransformedTriangle::areDoubleProductsConsistent(const TriSegment seg) const { const double term1 = _doubleProducts[8*seg + C_YZ] * _doubleProducts[8*seg + C_XH]; const double term2 = _doubleProducts[8*seg + C_ZX] * _doubleProducts[8*seg + C_YH]; const double term3 = _doubleProducts[8*seg + C_XY] * _doubleProducts[8*seg + C_ZH]; - LOG(2, "for seg " << seg << " consistency " << term1 + term2 + term3 ); LOG(2, "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 ); @@ -173,11 +183,9 @@ namespace INTERP_UTILS } return !((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 ); - // return (num_zero == 3) || (num_neg != 0 && num_neg != 3 && num_pos != 0 && num_pos != 3); - // return epsilonEqual(term1 + term2 + term3, 0.0); } -#ifndef OPTIMIZE +#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner) { // set the three corresponding double products to 0.0 @@ -197,6 +205,14 @@ namespace INTERP_UTILS }; } #endif OPTIMIZE + + /* + * Calculate the shortest distance between a tetrahedron corner and a triangle segment. + * + * @param corner corner of the tetrahedron + * @param seg segment of the triangle + * @return shortest distance from the corner to the segment + */ double TransformedTriangle::calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const { // NB uses fact that TriSegment <=> TriCorner that is first point of segment (PQ <=> P) @@ -254,7 +270,6 @@ namespace INTERP_UTILS { LOG(6, "- Triple product for corner " << corner ); - for(int row = 1 ; row < 4 ; ++row) { const DoubleProduct dp = DP_FOR_DETERMINANT_EXPANSION[3*corner + (row - 1)]; @@ -283,12 +298,10 @@ namespace INTERP_UTILS if(minAngle < TRIPLE_PRODUCT_ANGLE_THRESHOLD) { _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, true); - //_tripleProducts[corner] = calcTByDevelopingRow(corner, 1, false); } else { _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false); - //_tripleProducts[corner] = calcTByDevelopingRow(corner, 1, false); } _validTP[corner] = true; } @@ -308,6 +321,12 @@ namespace INTERP_UTILS _isTripleProductsCalculated = true; } + /* + * Calculates the angle between an edge of the tetrahedron and the triangle + * + * @param edge edge of the tetrahedron + * @return angle between triangle and edge + */ double TransformedTriangle::calculateAngleEdgeTriangle(const TetraEdge edge) const { // find normal to PQR - cross PQ and PR @@ -350,9 +369,6 @@ namespace INTERP_UTILS const double lenNormal = norm(normal); const double lenEdgeVec = norm(edgeVec); const double dotProd = dot(normal, edgeVec); - // return asin( dotProd / ( lenNormal * lenEdgeVec ) ); - //# assert(dotProd / ( lenNormal * lenEdgeVec ) + 1.0 >= 0.0); - //# assert(dotProd / ( lenNormal * lenEdgeVec ) - 1.0 <= 0.0); //? is this more stable? -> no subtraction // return asin( dotProd / ( lenNormal * lenEdgeVec ) ) + 3.141592625358979 / 2.0; @@ -368,7 +384,7 @@ namespace INTERP_UTILS * @param seg segment of triangle * @param dp double product sought * - * @returns stabilised double product c_{xy}^{pq} + * @return stabilised double product c_{xy}^{pq} * */ double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const @@ -386,7 +402,7 @@ namespace INTERP_UTILS * @pre double products have already been calculated * @pre triple products have already been calculated * @param corner corner for which the triple product is calculated - * @returns triple product associated with corner (see Grandy, eqs. [50]-[52]) + * @return triple product associated with corner (see Grandy, eqs. [50]-[52]) */ double TransformedTriangle::calcStableT(const TetraCorner corner) const { @@ -404,7 +420,7 @@ namespace INTERP_UTILS * @param seg segment of triangle * @param dp double product sought * - * @returns double product c_{xy}^{pq} + * @return double product c_{xy}^{pq} * */ double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const @@ -435,7 +451,7 @@ namespace INTERP_UTILS * @param corner corner for which the triple product is calculated * @param row row (1 <= row <= 3) used to calculate the determinant * @param project indicates whether or not to perform projection as inidicated in Grandy, p.446 - * @returns triple product associated with corner (see Grandy, [50]-[52]) + * @return triple product associated with corner (see Grandy, [50]-[52]) */ double TransformedTriangle::calcTByDevelopingRow(const TetraCorner corner, const int row, const bool project) const diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index c969b9af0..f8793e77b 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -89,7 +89,7 @@ namespace INTERP_UTILS * @param x first value * @param y second value * @param errTol maximum allowed absolute difference that is to be treated as equality - * @returns true if |x - y| < errTol, false otherwise + * @return true if |x - y| < errTol, false otherwise */ inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL) { @@ -106,7 +106,7 @@ namespace INTERP_UTILS * @param y second value * @param relTol maximum allowed relative difference that is to be treated as equality * @param absTol maximum allowed absolute difference that is to be treated as equality - * @returns true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise + * @return true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise */ inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL) { -- 2.39.2