* 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
{
* 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
* 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
{
*
* @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)
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]);
for(vector<int>::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)
{
delete targetElems[i];
}
+ std::cout << "Intersector filtered out " << intersector->filtered << " elements" << std::endl;
delete intersector;
}
*
*/
Intersector3D::Intersector3D(const MESH& srcMesh, const MESH& targetMesh)
- : _srcMesh(srcMesh), _targetMesh(targetMesh)
+ : _srcMesh(srcMesh), _targetMesh(targetMesh), filtered(0)
{
}
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
* @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<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
{
-
+
// get cell model for the element
CELLMODEL cellModel(type);
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:
//? 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;
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
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<INTERP_UTILS::TransformedTriangle>& 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);
}
};
BBTree.H\
MeshUtils.hxx\
Intersector3D.hxx\
-Log.hxx
+Log.hxx\
+TransformedTriangle_inline.hxx
+
# Libraries targets
TetraAffineTransform.cxx\
Interpolation3D.cxx\
Intersector3D.cxx\
-#Log.cxx
-#Interpolation3DSurf.cxx\
-#Interpolation2D.cxx\
+Interpolation3DSurf.cxx\
+Interpolation2D.cxx
# Executables targets
BIN =
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
# 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
#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.
#LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lutil
# for gcov
#LDFLAGS+= -lgcov
-LDFLAGS+= -pg
+#LDFLAGS+= -pg
#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
LDFLAGSFORBIN+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace -lSALOMEBasics
endif
-LDFLAGSFORBIN+= -pg
+#LDFLAGSFORBIN+= -pg
LIBSFORBIN=$(BOOSTLIBS) $(MPI_LIBS)
/*
* 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)
#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
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<int, double>::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<int, double>& col = m[i];
+ for(map<int, double>::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<double>* f=mesh.getVolume(sup);
+ tab = f->getValue();
+ // delete sup;
+}
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 = " <<sum_row - sVol[i] );
+ }
+
+ // target elements
+ const double* tVol = new double[tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
+ getVolumes(tMesh, tVol);
+ for(int i = 0; i < tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
+ {
+ const double sum_col = sumCol(m, i);
+ if(!epsilonEqualRelative(sum_col, tVol[i], VOL_PREC))
+ {
+ LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
+ ok = false;
+ }
+ LOG(1, "diff = " <<sum_col - tVol[i] );
+ }
+ delete[] sVol;
+ delete[] tVol;
+ return ok;
+}
+
bool Interpolation3DTest::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
{
bool compatitable = true;
LOG(1, std::endl << "=== -> 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 );
}
{
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));
#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 );
// 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();
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;
void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const;
+
};
#endif
#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
// --- 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 );
// 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):
// -----------------------------------------------------------------------------------------------------
#include "VectorUtils.hxx"
-#undef MERGE_CALC
-#undef COORDINATE_CORRECTION 1.0e-15
-
class ProjectedCentralCircularSortOrder
* 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
_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();
}
+ /*
+ * Destructor
+ *
+ * Deallocates the memory used to store the points of the polygons.
+ * This memory is allocated in calculateIntersectionPolygons().
+ */
TransformedTriangle::~TransformedTriangle()
{
// delete elements of polygons
* 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
*
*/
}
- // 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
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 ... ");
}
LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
-
+
return sign * (volA + volB);
}
// 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);
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))
{
// 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)
*
* @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)
* 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)
{
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<double*> pAUnique;
- pAUnique.reserve(_polygonA.size());
- Vector3Cmp cmp;
- unique_copy(_polygonA.begin(), _polygonA.end(), back_inserter(pAUnique), cmp);
- //for(std::vector<double*>::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<double*>::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;
}
return true;
}
+ /*
+ * Prints the coordinates of the triangle to std::cout
+ *
+ */
void TransformedTriangle::dumpCoords()
{
std::cout << "Coords : ";
#include <vector>
-#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
* 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.
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 };
void dumpCoords();
-
private:
////////////////////////////////////////////////////////////////////////////////////
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 ////////
////////////////////////////////////////////////////////////////////////////////////
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;
+#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
};
}
-/**
- * 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);
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
{
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] =
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
#define SEG_RAY_TABLE 1 // seems correct
-
-
-
-
namespace INTERP_UTILS
{
////////////////////////////////////////////////////////////////////////////////////
- /// Constants /////////////////
+ /// Correspondance tables describing all the variations of formulas. //////////////
////////////////////////////////////////////////////////////////////////////////////
// correspondance facet - double product
C_ZH, C_ZX, C_YZ, // OXY
C_XH, C_YH, C_ZH // XYZ
};
+#if 0
+ template<TetraFacet facet, TriSegment seg>
+ inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection()
+ {
+ return NO_DP;
+ }
+
+ template<>
+ inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection<OYZ, PQ>
+ {
+ return C_XH;
+ }
+
+#define DEF_DP_FOR_SEG_FACET(FACET,SEG,DP) template<> inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection<FACET,SEG> { 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] =
{
};
#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])
}
}
-#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
{
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);
/**
* 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;
}
}
// 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);
/**
* 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
{
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;
}
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
{
*
* @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)
{
/**
* 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
{
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;
}
-#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
{
* (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
{
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];
const double H = _coords[5*corner + 4];
return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
-
}
#endif
* 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));
// 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 )
{
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
//? 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
{
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
//? 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
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;
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
{
//? 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
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))
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;
}
_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 );
}
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
};
}
#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)
{
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)];
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;
}
_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
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;
* @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
* @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
{
* @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
* @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
* @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)
{
* @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)
{
