* LOG_LEVEL (which can be defined at compile-time for each file by passing option -DLOG_LEVEL=x to gcc)
* than the message is logged.
*
- *
- *
*/
+
+//#define LOG_LEVEL 4
/// define LOG_LEVEL here if it is not already defined
#ifndef LOG_LEVEL
#define LOG_LEVEL 0
-
-
-
-
-
-
-
-
#endif
* which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
* the corresponding array element will be true if and only if it was true when the method was called and pt lies in the halfspace.
*
- * @param pt double[3] containing the coordiantes of a transformed point
+ * @param pt double[3] containing the coordinates of a transformed point
* @param isOutside bool[8] which indicate the results of earlier checks.
*/
template<class MyMeshType>
// halfspace filtering check
// NB : might not be beneficial for caching of triangles
for(int i = 0; i < 8; ++i)
- {
- if(isOutside[i])
- {
- isTargetOutside = true;
- }
- }
+ if(isOutside[i])
+ isTargetOutside = true;
double totalVolume = 0.0;
TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
double vol = tri.calculateIntersectionVolume();
+ LOG(1, "ii = " << ii << " Volume=" << vol)
totalVolume += vol;
}
TransformedTriangle::~TransformedTriangle()
{
// delete elements of polygons
- for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
- {
- delete[] *it;
- }
- for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
- {
- delete[] *it;
- }
+ for(auto& it: _polygonA)
+ delete[] it;
+ for(auto& it: _polygonB)
+ delete[] it;
}
/**
if(_polygonA.size() > 2)
{
LOG(2, "---- Treating polygon A ... ");
+#if LOG_LEVEL > 0
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: _polygonA)
+ LOG(3,vToStr(pt));
+#endif
calculatePolygonBarycenter(A, barycenter);
sortIntersectionPolygon(A, barycenter);
volA = calculateVolumeUnderPolygon(A, barycenter);
double volB = 0.0;
// if triangle is not in h = 0 plane, calculate volume under B
- if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
+ if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) // second condition avoids double counting in case triangle fully included in h=0 facet
{
LOG(2, "---- Treating polygon B ... ");
-
+#if LOG_LEVEL > 0
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: _polygonB)
+ LOG(3,vToStr(pt));
+#endif
calculatePolygonBarycenter(B, barycenter);
sortIntersectionPolygon(B, barycenter);
volB = calculateVolumeUnderPolygon(B, barycenter);
LOG(2, "Volume is " << sign * volB);
}
- LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
+#if LOG_LEVEL >= 2
+ LOG(2, "############ Triangle :")
+ dumpCoords();
+ LOG(2, "vol A = " << volA);
+ LOG(2, "vol B = " << volB);
+ LOG(2, "TOTAL = " << sign*(volA+volB));
+#endif
return _volume = sign * (volA + volB);
}
// ----------------------------------------------------------------------------------
- // TransformedTriangle PRIVATE
+ // TransformedTriangle PROTECTED
// ----------------------------------------------------------------------------------
/**
*/
void TransformedTriangle::calculateIntersectionAndProjectionPolygons()
{
+#if LOG_LEVEL >= 3
+ std::cout << " @@@@@@@@ COORDS @@@@@@ " << std::endl;
+ dumpCoords();
+#endif
+
assert(_polygonA.size() == 0);
assert(_polygonB.size() == 0);
// avoid reallocations in push_back() by pre-allocating enough memory
double* ptA = new double[3];
calcIntersectionPtSurfaceEdge(edge, ptA);
_polygonA.push_back(ptA);
- LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
+ LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptA) << " added to A ");
if(edge >= XY)
{
double* ptB = new double[3];
copyVector3(ptA, ptB);
_polygonB.push_back(ptB);
- LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B ");
+ LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptB) << " added to B ");
}
-
+
}
}
double* ptB = new double[3];
copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
_polygonB.push_back(ptB);
- LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
+ LOG(3,"Surface-ray (corner " << strTC(corner) << "): " << vToStr(ptB) << " added to B");
}
}
bool isZero[NO_DP];
- // check beforehand which double-products are zero
+ // check beforehand which double-products are zero.
for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
- {
- isZero[dp] = (calcStableC(seg, dp) == 0.0);
- }
+ isZero[dp] = (calcStableC(seg, dp) == 0.0);
// segment - facet
for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
{
// is this test worth it?
const bool doTest =
- !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
- !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
- !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
+ !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
+ !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
+ !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
if(doTest && testSegmentFacetIntersection(seg, facet))
{
double* ptA = new double[3];
calcIntersectionPtSegmentFacet(seg, facet, ptA);
_polygonA.push_back(ptA);
- LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
+ LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
if(facet == XYZ)
{
double* ptB = new double[3];
copyVector3(ptA, ptB);
_polygonB.push_back(ptB);
- LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
+ LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
}
}
double* ptA = new double[3];
calcIntersectionPtSegmentEdge(seg, edge, ptA);
_polygonA.push_back(ptA);
- LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
+ LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
if(edge >= XY)
{
double* ptB = new double[3];
copyVector3(ptA, ptB);
_polygonB.push_back(ptB);
+ LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to B");
}
}
}
-
+
// segment - corner
for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
{
const bool doTest =
- isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
- isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
- isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
+ isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
+ isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
+ isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
if(doTest && testSegmentCornerIntersection(seg, corner))
{
double* ptA = new double[3];
copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
_polygonA.push_back(ptA);
- LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
+ LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
if(corner != O)
{
double* ptB = new double[3];
_polygonB.push_back(ptB);
copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
+ LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
}
}
}
- // segment - ray
- for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
+ // segment - ray
+ for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
+ {
+ if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
{
- if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
- {
- double* ptB = new double[3];
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
- }
+ double* ptB = new double[3];
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-ray (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
}
-
- // segment - halfstrip
- for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
- {
+ }
+
+ // segment - halfstrip
+ for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
+ {
#if 0
- const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
- const bool doTest =
- !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
- !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
-
+ const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
+ const bool doTest =
+ !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
+ !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
+
- if(doTest && testSegmentHalfstripIntersection(seg, edge))
+ if(doTest && testSegmentHalfstripIntersection(seg, edge))
#endif
- if(testSegmentHalfstripIntersection(seg, edge))
- {
- double* ptB = new double[3];
- calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
- }
- }
+ if(testSegmentHalfstripIntersection(seg, edge))
+ {
+ double* ptB = new double[3];
+ calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
+ }
+ }
}
- // inclusion tests
- for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
+ // inclusion tests
+ for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
+ {
+ // { XYZ - inclusion only possible if in Tetrahedron?
+ // tetrahedron
+ if(testCornerInTetrahedron(corner))
{
- // { XYZ - inclusion only possible if in Tetrahedron?
- // tetrahedron
- if(testCornerInTetrahedron(corner))
- {
- double* ptA = new double[3];
- copyVector3(&_coords[5*corner], ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
- }
-
- // XYZ - plane
- if(testCornerOnXYZFacet(corner))
- {
- double* ptB = new double[3];
- copyVector3(&_coords[5*corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
- }
+ double* ptA = new double[3];
+ copyVector3(&_coords[5*corner], ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Inclusion tetrahedron (corner " << strTriC(corner) << "): " << vToStr(ptA) << " added to A");
+ }
- // projection on XYZ - facet
- if(testCornerAboveXYZFacet(corner))
- {
- double* ptB = new double[3];
- copyVector3(&_coords[5*corner], ptB);
- ptB[2] = 1 - ptB[0] - ptB[1];
- assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
- _polygonB.push_back(ptB);
- LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
- }
+ // XYZ - plane
+ if(testCornerOnXYZFacet(corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&_coords[5*corner], ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Inclusion XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B");
+ }
+ // projection on XYZ - facet
+ if(testCornerAboveXYZFacet(corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&_coords[5*corner], ptB);
+ ptB[2] = 1 - ptB[0] - ptB[1];
+ assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
+ _polygonB.push_back(ptB);
+ LOG(3,"Projection XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B");
}
}
+ }
+
/**
* Calculates the intersection polygon A, performing the intersection tests
* and storing the corresponding point in the vector _polygonA.
bool isZero[NO_DP];
// check beforehand which double-products are zero
+ // Test for "== 0.0" here is OK since doubleProduct has been fixed for rounding to zero already.
for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
- {
- isZero[dp] = (calcStableC(seg, dp) == 0.0);
- }
+ isZero[dp] = (calcStableC(seg, dp) == 0.0);
// segment - facet
for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
*/
bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
{
-
// coordinate to check
const int coord = static_cast<int>(facet);
for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
- {
- if(_coords[5*c + coord] != 0.0)
- {
- return false;
- }
- }
-
+ if(_coords[5*c + coord] != 0.0)
+ return false;
+
return true;
}
double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
if(epsilonEqual(sign, 0.))
- {
- sign = 0.;
- }
+ sign = 0.;
return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0;
}
bool TransformedTriangle::isTriangleBelowTetraeder() const
{
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)
- {
- return false;
- }
- }
+ // check z-coords for all points
+ if(_coords[5*c + 2] >= 0.0)
+ return false;
+
return true;
}
{
std::cout << "Coords : ";
for(int i = 0 ; i < 3; ++i)
- {
- std::cout << vToStr(&_coords[5*i]) << ",";
- }
+ std::cout << vToStr(&_coords[5*i]) << ",";
+
std::cout << std::endl;
}
-
+
+
} // NAMESPACE
};
- // include definitions of inline methods
+ 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
+ static const DoubleProduct DOUBLE_PRODUCTS[12] =
+ {
+ C_YZ, C_ZX, C_XY, // O
+ C_YZ, C_ZH, C_YH, // X
+ C_ZX, C_ZH, C_XH, // Y
+ C_XY, C_YH, C_XH // Z
+ };
+
+ for(int i = 0 ; i < 3 ; ++i) {
+ const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i];
+
+ LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner);
+ _doubleProducts[8*seg + dp] = 0.0;
+ };
+ }
+
+ inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const
+ {
+ return _doubleProducts[8*seg + dp];
+ }
+
+ inline double TransformedTriangle::calcStableT(const TetraCorner corner) const
+ {
+ // assert(_isTripleProductsCalculated);
+ // assert(_validTP[corner]);
+ return _tripleProducts[corner];
+ }
+
+ inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const
+ {
+
+ // find the points of the triangle
+ // 0 -> P, 1 -> Q, 2 -> R
+ const int pt1 = seg;
+ const int pt2 = (seg + 1) % 3;
+
+ // find offsets
+ const int off1 = DP_OFFSET_1[dp];
+ const int off2 = DP_OFFSET_2[dp];
+
+ return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
+ }
+
+ // ----------------------------------------------------------------------------------
+ // TransformedTriangle_intersect.cxx
+ // ----------------------------------------------------------------------------------
+ inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const
+ {
+ return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge);
+ }
+
+ inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const
+ {
+ return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet);
+ }
+
+ inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const
+ {
+ return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner );
+ }
+
+ inline bool TransformedTriangle::testCornerInTetrahedron(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] // z
+ };
+
+ for(int i = 0 ; i < 4 ; ++i)
+ {
+ if(pt[i] < 0.0 || pt[i] > 1.0)
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const
+ {
+#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;
+ }
+
+ 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];
+
+ return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
+
+ }
+
+ inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const
+ {
+
+ // assert(edge < H01);
+
+ // correspondence 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) && !epsilonEqual(t1 - t2, 0.0); // tuleap26461
+ }
+
+ 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 correspondence facet a = 0 <=> offset for coordinate a in _coords
+ // and also correspondence 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 );
+
+ }
-#include "TransformedTriangleInline.hxx"
}
+++ /dev/null
-// Copyright (C) 2007-2024 CEA, EDF
-//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-//
-
-#ifndef __TRANSFORMEDTRIANGLEINLINE_HXX__
-#define __TRANSFORMEDTRIANGLEINLINE_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
- static const DoubleProduct DOUBLE_PRODUCTS[12] =
- {
- C_YZ, C_ZX, C_XY, // O
- C_YZ, C_ZH, C_YH, // X
- C_ZX, C_ZH, C_XH, // Y
- C_XY, C_YH, C_XH // Z
- };
-
- for(int i = 0 ; i < 3 ; ++i) {
- const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i];
-
- LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner);
- _doubleProducts[8*seg + dp] = 0.0;
- };
-}
-
-inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const
-{
- return _doubleProducts[8*seg + dp];
-}
-
-inline double TransformedTriangle::calcStableT(const TetraCorner corner) const
-{
- // assert(_isTripleProductsCalculated);
- // assert(_validTP[corner]);
- return _tripleProducts[corner];
-}
-
-inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const
-{
-
- // find the points of the triangle
- // 0 -> P, 1 -> Q, 2 -> R
- const int pt1 = seg;
- const int pt2 = (seg + 1) % 3;
-
- // find offsets
- const int off1 = DP_OFFSET_1[dp];
- const int off2 = DP_OFFSET_2[dp];
-
- return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
-}
-
-// ----------------------------------------------------------------------------------
-// TransformedTriangle_intersect.cxx
-// ----------------------------------------------------------------------------------
-inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const
-{
- return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge);
-}
-
-inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const
-{
- return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet);
-}
-
-inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const
-{
- return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner );
-}
-
-inline bool TransformedTriangle::testCornerInTetrahedron(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] // z
- };
-
- for(int i = 0 ; i < 4 ; ++i)
- {
- if(pt[i] < 0.0 || pt[i] > 1.0)
- {
- return false;
- }
- }
- return true;
-}
-
-inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const
-{
-#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;
-}
-
-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];
-
- return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
-
-}
-
-inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const
-{
-
- // assert(edge < H01);
-
- // correspondence 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) && !epsilonEqual(t1 - t2, 0.0); // tuleap26461
-}
-
-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 correspondence facet a = 0 <=> offset for coordinate a in _coords
- // and also correspondence 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
const double alpha = tA / (tA - tB);
// calculate point
- LOG(4, "corner A = " << corners[0] << " corner B = " << corners[1] );
+ LOG(4, "corner A = " << strTC(corners[0]) << " corner B = " << strTC(corners[1]) );
LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha );
for(int i = 0; i < 3; ++i)
{
- pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] +
- alpha * COORDS_TET_CORNER[3*corners[1] + i];
+ pt[i] = (1 - alpha)*COORDS_TET_CORNER[3*corners[0] + i] + alpha*COORDS_TET_CORNER[3*corners[1] + i];
#if 0
pt[i] = (1 - alpha) * getCoordinateForTetCorner<corners[0], i>() +
alpha * getCoordinateForTetCorner<corners[0], i>();
/**
* 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
* @return true if the segment intersects the edge
for(int i = 0 ; i < 2 ; ++i)
{
facet[i] = FACET_FOR_EDGE[2*edge + i];
-
+
// find the two c-values -> the two for the other edges of the facet
int idx1 = 0 ;
int idx2 = 1;
DoubleProduct dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
DoubleProduct dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
-
+
if(dp1 == DoubleProduct( edge ))
{
idx1 = 2;
idx2 = 2;
dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
}
-
+
const double c1 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1]*calcStableC(seg, dp1);
const double c2 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2]*calcStableC(seg, dp2);
const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
c[j] = dpIdx < 0.0 ? 0.0 : sign * calcStableC(seg, dp);
}
-
- // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2)
pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
-
- // strange bug with -O2 enabled : assertion fails when we don't have the following
- // trace - line
- //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]). If OPTIMIZE is defined, the double products that should be zero are not verified.
+ * (Grandy, eq. [21]).
*
* @param seg segment of the triangle
* @param corner corner of the tetrahedron
*/
bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
{
-
-
// facets meeting at a given corner
static const TetraFacet FACETS_FOR_CORNER[12] =
{
{
const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
if(testSegmentIntersectsFacet(seg, facet))
- {
- return true;
- }
+ return true;
}
return false;
};
const TetraFacet facet = FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
-
// special case : facet H = 0
const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
- LOG(4, "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
+ LOG(4, "Halfstrip tests (" << strTriS(seg) << ", " << strTE(edge) << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] );
return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
/**
* 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)
// if two or more c-values are zero we disallow x-edge intersection
// Grandy, p.446
+ // test for == 0.0 is OK here, if you look at the way double products were pre-comuted:
const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
if(numZeros >= 2 )
#include <algorithm>
#include <sstream>
+#include <iomanip>
#include <numeric>
#include <string>
#include <cmath>
inline const std::string vToStr(const double* pt)
{
std::stringstream ss(std::ios::out);
- ss << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]";
+ ss << std::setprecision(16) << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]";
return ss.str();
}
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] );
+ LOG(2, "diff = " <<sum_row - sVol[i] );
}
// target elements
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] );
+ LOG(2, "diff = " <<sum_col - tVol[i] );
}
delete[] sVol;
delete[] tVol;
IntersectionMatrix matrix1;
calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
-
#if LOG_LEVEL >= 2
dumpIntersectionMatrix(matrix1);
#endif
void MeshTestToolkit<SPACEDIM,MESHDIM>::intersectMeshes(const char* mesh1, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
{
const std::string path1 = std::string(mesh1) + std::string(".med");
- std::cout << "here :" << path1 << std::endl;
const std::string path2 = std::string(mesh2) + std::string(".med");
+ std::cout << "here :" << path1 << " with " << path2 << std::endl;
intersectMeshes(path1.c_str(), mesh1, path2.c_str(), mesh2, correctVol, prec, doubleTest);
}
// std::pair<int, int> eff = countNumberOfMatrixEntries(m);
// LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = "
// << double(eff.first) / double(numTargetElems * numSrcElems));
- LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = "
- << double(eff.second) / double(eff.first));
+// LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = "
+// << double(eff.second) / double(eff.first));
LOG(1, "Intersection calculation done. " << std::endl );
#ifndef MEDCOUPLING_MICROMED
// These test suites need MEDLoader to load some test files:
CPPUNIT_TEST_SUITE_REGISTRATION( InterpolationOptionsTest );
+CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementTetraTests );
CPPUNIT_TEST_SUITE_REGISTRATION( HexaTests );
CPPUNIT_TEST_SUITE_REGISTRATION( MultiElement2DTests );
CPPUNIT_TEST_SUITE_REGISTRATION( MultiElementTetraTests );
-CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementTetraTests );
CPPUNIT_TEST_SUITE_REGISTRATION( ThreeDSurfProjectionTest );
#endif
// --- generic Main program from KERNEL_SRC/src/Basics/Test
