Salome HOME
projects
/
tools
/
medcoupling.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Add test for .mesh file format
[tools/medcoupling.git]
/
src
/
INTERP_KERNEL
/
TransformedTriangleIntersect.cxx
diff --git
a/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx
b/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx
index fa4c8a53581e8b5498cd5766f531a38ea0964b78..36727ef2624ebfce5ed18c859c807574e8966583 100644
(file)
--- a/
src/INTERP_KERNEL/TransformedTriangleIntersect.cxx
+++ b/
src/INTERP_KERNEL/TransformedTriangleIntersect.cxx
@@
-1,4
+1,4
@@
-// Copyright (C) 2007-20
15 CEA/DEN, EDF R&D
+// Copyright (C) 2007-20
24 CEA, EDF
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@
-28,10
+28,10
@@
namespace INTERP_KERNEL
{
// ----------------------------------------------------------------------------------
{
// ----------------------------------------------------------------------------------
- // Correspond
a
nce tables describing all the variations of formulas.
+ // Correspond
e
nce tables describing all the variations of formulas.
// ----------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------
- /// \brief Correspond
a
nce between facets and double products.
+ /// \brief Correspond
e
nce between facets and double products.
///
/// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] =
///
/// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] =
@@
-78,7
+78,7
@@
namespace INTERP_KERNEL
9, 10, 11 // XYZ
};
9, 10, 11 // XYZ
};
- /// \brief Correspond
a
nce edge - corners.
+ /// \brief Correspond
e
nce edge - corners.
///
/// Gives the two corners associated with each edge
/// Use 2*edge + {0, 1} as index
///
/// Gives the two corners associated with each edge
/// Use 2*edge + {0, 1} as index
@@
-92,7
+92,7
@@
namespace INTERP_KERNEL
Z, X // ZX
};
Z, X // ZX
};
- /// \brief Correspond
a
nce edge - facets.
+ /// \brief Correspond
e
nce edge - facets.
///
/// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index
const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
///
/// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index
const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
@@
-105,7
+105,7
@@
namespace INTERP_KERNEL
OZX, XYZ // ZX
};
OZX, XYZ // ZX
};
- /// \brief Correspond
a
nce corners - edges.
+ /// \brief Correspond
e
nce corners - edges.
///
/// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index
const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
///
/// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index
const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
@@
-168,13
+168,12
@@
namespace INTERP_KERNEL
const double alpha = tA / (tA - tB);
// calculate point
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)
{
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>();
#if 0
pt[i] = (1 - alpha) * getCoordinateForTetCorner<corners[0], i>() +
alpha * getCoordinateForTetCorner<corners[0], i>();
@@
-234,7
+233,7
@@
namespace INTERP_KERNEL
/**
* Tests if the given segment of the triangle intersects the given edge of the tetrahedron (Grandy, eq. [20]
/**
* 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
* @param seg segment of the triangle
* @param edge edge of tetrahedron
* @return true if the segment intersects the edge
@@
-249,13
+248,13
@@
namespace INTERP_KERNEL
for(int i = 0 ; i < 2 ; ++i)
{
facet[i] = FACET_FOR_EDGE[2*edge + i];
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];
// 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;
if(dp1 == DoubleProduct( edge ))
{
idx1 = 2;
@@
-266,7
+265,7
@@
namespace INTERP_KERNEL
idx2 = 2;
dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
}
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 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);
@@
-343,28
+342,24
@@
namespace INTERP_KERNEL
for(int j = 0 ; j < 2; ++j)
{
const int dpIdx = DP_INDEX[3*facets[j] + i];
for(int j = 0 ; j < 2; ++j)
{
const int dpIdx = DP_INDEX[3*facets[j] + i];
- const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
- const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
- c[j] = dpIdx < 0.0 ? 0.0 : sign * calcStableC(seg, dp);
+ if (dpIdx != -1)
+ {
+ const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
+ const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
+ c[j] = sign * calcStableC(seg, dp);
+ }
+ else
+ c[j] = 0.0;
}
}
-
- // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2)
pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
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.
}
}
/**
* 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
*
* @param seg segment of the triangle
* @param corner corner of the tetrahedron
@@
-372,8
+367,6
@@
namespace INTERP_KERNEL
*/
bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
{
*/
bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
{
-
-
// facets meeting at a given corner
static const TetraFacet FACETS_FOR_CORNER[12] =
{
// facets meeting at a given corner
static const TetraFacet FACETS_FOR_CORNER[12] =
{
@@
-388,9
+381,7
@@
namespace INTERP_KERNEL
{
const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
if(testSegmentIntersectsFacet(seg, facet))
{
const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
if(testSegmentIntersectsFacet(seg, facet))
- {
- return true;
- }
+ return true;
}
return false;
}
return false;
@@
-432,11
+423,10
@@
namespace INTERP_KERNEL
};
const TetraFacet facet = FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
};
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);
// special case : facet H = 0
const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
- LOG(4, "Halfstrip tests (" << s
eg << ", " << edge
<< ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
+ LOG(4, "Halfstrip tests (" << s
trTriS(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);
LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] );
return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
@@
-497,7
+487,6
@@
namespace INTERP_KERNEL
/**
* 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])
/**
* 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)
*
* @param seg segment of the triangle PQR
* @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z)
@@
-573,6
+562,7
@@
namespace INTERP_KERNEL
// if two or more c-values are zero we disallow x-edge intersection
// Grandy, p.446
// 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 )
const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
if(numZeros >= 2 )