// minimum coordinate of one is greater than the maximum coordinate of the other
//? stable version
- /* const double tol = 1.0e-2;
+ /*const double tol = 1.0e-2*_coords[c];
if(_coords[c] > otherMaxCoord + tol
|| _coords[c + 3] < otherMinCoord - tol)
*/
+
if(_coords[c] > otherMaxCoord
|| _coords[c + 3] < otherMinCoord)
- {
+
+ {
return true;
}
+
}
return false;
}
#include "RegionNode.hxx"
#include "TetraAffineTransform.hxx"
#include "TransformedTriangle.hxx"
+#include "VectorUtils.hxx"
using namespace INTERP_UTILS;
using namespace MEDMEM;
iter != currNode->getTargetRegion().getEndElements() ; ++iter)
{
const double vol = calculateIntersectionVolume(*srcElement, **iter);
- if(vol != 0.0)
+ if(!epsilonEqual(vol, 0.0, 1.0e-10))
+ // if(vol != 0.0)
{
const int targetIdx = (*iter)->getIndex();
// (a), (b) and (c) not yet implemented
// (d) : without fine-level filtering (a) - (c) for the time being
-
- // std::cout << "Source : ";
- // srcElement.dumpCoords();
- // std::cout << "Target : ";
- // targetElement.dumpCoords();
+ // std::cout << "------------------" << std::endl;
+ // std::cout << "Source : ";
+ srcElement.dumpCoords();
+ // std::cout << "Target : ";
+ targetElement.dumpCoords();
// get array of points of target tetraeder
const double* tetraCorners[4];
TetraAffineTransform T( tetraCorners );
// check if we have planar tetra element
- if(T.determinant() == 0.0)
+ if(epsilonEqual(T.determinant(), 0.0, 1.0e-16))
{
// tetra is planar
// std::cout << "Planar tetra -- volume 0" << std::endl;
for(vector<TransformedTriangle>::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter)
{
// std::cout << std::endl << "= > Triangle " << ++i << std::endl;
- // iter->dumpCoords();
+ iter->dumpCoords();
volume += iter->calculateIntersectionVolume();
}
//? fault in article, Grandy, [8] : it is the determinant of the inverse transformation that
// should be used
+
return std::abs(1.0 / T.determinant() * volume) ;
}
assert(faceType == MED_TRIA3);
// create transformed triangles from face
+
switch(faceType)
{
case MED_TRIA3:
+ // std::cout << std::endl<< "** Adding triangle " << i << std::endl;
triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6]));
break;
void MeshElement::dumpCoords() const
{
- std::cout << "Element " << _index + 1 << " has nodes " << std::endl;
+ // std::cout << "Element " << _index + 1 << " has nodes " << std::endl;
for(int i = 1 ; i <= getNumberNodes() ; ++i)
{
- std::cout << vToStr(getCoordsOfNode(i)) << ", ";
+ // std::cout << vToStr(getCoordsOfNode(i)) << ", ";
}
- std::cout << std::endl;
+ // std::cout << std::endl;
}
#include "VectorUtils.hxx"
+#undef NULL_COORD_CORRECTION
+
namespace INTERP_UTILS
{
// and O is the position vector of the point that is mapped onto the origin
for(int i = 0 ; i < 3 ; ++i)
{
- _translation[i] = -_linearTransform[3*i]*(pts[0])[0] - _linearTransform[3*i+1]*(pts[0])[1] - _linearTransform[3*i+2]*(pts[0])[2] ;
+ _translation[i] = -(_linearTransform[3*i]*(pts[0])[0] + _linearTransform[3*i+1]*(pts[0])[1] + _linearTransform[3*i+2]*(pts[0])[2]) ;
}
// precalculate determinant (again after inversion of transform)
// translation
dest[i] += _translation[i];
+
+#ifdef NULL_COORD_CORRECTION
+ if(epsilonEqual(dest[i], 0.0, 1.0e-15))
+ {
+ dest[i] = 0.0;
+ }
+#endif
+
}
if(selfAllocation)
{
//{ we copy the matrix for the lu-factorization
// maybe inefficient
+
double lu[9];
for(int i = 0 ; i < 9; ++i)
{
// calculate LU factorization
int idx[3];
- factorizeLU(lu, idx);
+ // double parity = 1.0;
+ factorizeLU(lu, idx/*, &parity*/);
+
+ //_determinant = parity / (lu[3*idx[0]] * lu[3*idx[1] + 1] * lu[3*idx[2] + 2]);
+ //std::cout << "lu-determinant = " << _determinant << std::endl;
// calculate inverse by forward and backward substitution
// store in _linearTransform
/////////////////////////////////////////////////
/// Auxiliary methods for inverse calculation ///
/////////////////////////////////////////////////
- void factorizeLU(double* lu, int* idx) const
+ void factorizeLU(double* lu, int* idx/*, double* parity*/) const
{
// 3 x 3 LU factorization
// initialise idx
int tmp = idx[k];
idx[k] = idx[row];
idx[row] = tmp;
+ /* if(k != row)
+ *parity = -(*parity);
+ */
// calculate row
for(int j = k + 1 ; j < 3 ; ++j)
lu[3*idx[j] + k] /= lu[3*idx[k] + k];
for(int s = k + 1; s < 3 ; ++s)
{
- // case s = k will always become zero, and then replaced by
+ // case s = k will always become zero, and then be replaced by
// the l-value
// there's no need to calculate this explicitly
#include <cassert>
#include <algorithm>
#include <functional>
+#include <iterator>
#include "VectorUtils.hxx"
#include <math.h>
+#include <vector>
+
+#undef MERGE_CALC
+#undef COORDINATE_CORRECTION 1.0e-15
class CircularSortOrder
{
const double _a, _b;
};
+class Vector3Cmp
+{
+ public:
+ bool operator()(double* const& pt1, double* const& pt2)
+ {
+ // std::cout << "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2])) << std::endl;
+ return (pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]);
+ }
+};
+
namespace INTERP_UTILS
{
////////////////////////////////////////////////////////////////////////////
}
// h coordinate
+
_coords[5*P + 3] = 1 - p[0] - p[1] - p[2];
_coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
_coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
-
+
+ // order of substractions give different results ?
+ /*
+ _coords[5*P + 3] = 1 - (p[0] - (p[1] - p[2]));
+ _coords[5*Q + 3] = 1 - (q[0] - (q[1] - q[2]));
+ std::cout << "old = " << 1 - q[0] - q[1] - q[2] << " calculated = " << 1 - (q[0] - (q[1] - q[2])) << " stored : " << _coords[5*Q + 3] << std::endl;
+ _coords[5*R + 3] = 1 - (r[0] -(r[1] - r[2]));
+ */
+
+
// H coordinate
_coords[5*P + 4] = 1 - p[0] - p[1];
_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();
preCalculateTripleProducts();
double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
+
if(sign == 0.0)
{
// std::cout << std::endl << "Triangle is perpendicular to z-plane - V = 0.0" << std::endl << std::endl;
return 0.0;
}
+
// normalize
sign = sign > 0.0 ? 1.0 : -1.0;
// std::cout << std::endl << "-- Calculating intersection polygons ... " << std::endl;
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
{
// std::cout << std::endl << "-- Treating polygon B ... " << std::endl;
+ // std::cout << _coords[5*P + 3] << ", " << _coords[5*Q + 3] << ", " << _coords[5*R+ 3] << std::endl;
+
calculatePolygonBarycenter(B, barycenter);
sortIntersectionPolygon(B, barycenter);
volB = calculateVolumeUnderPolygon(B, barycenter);
vol += (factor1 * factor2) / 6.0;
}
- // // std::cout << "Abs. Volume is " << vol << std::endl;
+ // std::cout << "Abs. Volume is " << vol << std::endl;
return vol;
}
for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
{
// ? should have epsilon-equality here?
+ //if(!epsilonEqual(_coords[5*c + coord], 0.0, 1.0e-15))
if(_coords[5*c + coord] != 0.0)
{
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))
+ {
+ return false;
+ }
+ }
+ // std::cout << "Polygon A is on h = 0 facet" << std::endl;
+ return true;
+ }
bool TransformedTriangle::isTriangleBelowTetraeder()
{
void TransformedTriangle::dumpCoords()
{
- std::cout << "Coords : ";
+ // 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;
+ // std::cout << std::endl;
}
}; // NAMESPACE
bool isTriangleBelowTetraeder();
+ bool isPolygonAOnHFacet() const;
+
+
////////////////////////////////////////////////////////////////////////////////////
/// Intersection test methods and intersection point calculations ////////
////////////////////////////////////////////////////////////////////////////////////
void preCalculateDoubleProducts(void);
+ void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
+
double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const;
void preCalculateTripleProducts(void);
#include <cmath>
#include "VectorUtils.hxx"
+#define SEG_RAY_TABLE 1 // seems correct
+#undef TEST_EPS 0.0//1.0e-14
+
namespace INTERP_UTILS
{
// : (x,y,z)* = (1-alpha)*A + alpha*B where
// alpha = t_A / (t_A - t_B)
-
-
const TetraCorner corners[2] =
{
CORNERS_FOR_EDGE[2*edge],
const double alpha = tA / (tA - tB);
// calculate point
+ // std::cout << "corner A = " << corners[0] << " corner B = " << corners[1] << std::endl;
+ // std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl;
for(int i = 0; i < 3; ++i)
{
- // std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl;
pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] +
alpha * COORDS_TET_CORNER[3*corners[1] + i];
- // std::cout << pt[i] << std::endl;
+ // // std::cout << pt[i] << std::endl;
assert(pt[i] >= 0.0);
assert(pt[i] <= 1.0);
}
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;
- // std::cout << "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] << std::endl;
+ // std::cout << "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] << std::endl;
+ // std::cout << "c(" << seg << ", " << dp << ") = " << sign * calcStableC(seg, dp) << std::endl;
assert(pt[i] >= 0.0);
assert(pt[i] <= 1.0);
}
OZX, XYZ // ZX
};
- if(calcStableC(seg,DoubleProduct( edge )) != 0.0)
+ if(calcStableC(seg,DoubleProduct( edge )) != 0.0)
+ // if(!epsilonEqual(calcStableC(seg,DoubleProduct( edge )), 0.0, TEST_EPS))
{
return false;
}
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;
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);
- isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0;
+ //isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0;
+ if(c1*c2 > 0.0)
+ {
+ isFacetCondVerified = true;
+ }
}
+
if(!isFacetCondVerified)
{
return false;
const DoubleProduct dp = DoubleProduct( edge );
const double c = calcStableC(seg, dp);
if(c != 0.0)
+ // if(!epsilonEqual(c, 0.0, TEST_EPS))
{
return false;
}
// dp 1 -> cond 1
// dp 2-7 -> cond 3
//? NB : last two rows are not completely understood and may contain errors
- /*static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] =
+#if SEG_RAY_TABLE==1
+ static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] =
{
C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01 // Z
- };*/
+ };
+#else
+
static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] =
{
C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_YZ, // Y
C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_ZX // Z
};
- /*static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] =
- {
- C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
- C_01, C_XH, C_ZH, C_XY, C_10, C_XH, C_ZX, // Y
- C_XY, C_YH, C_XH, C_10, C_01, C_ZH, C_YZ // Z
- };
- */
+#endif
+
// facets to use
//? not sure this is correct
static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
// if two or more c-values are zero we disallow x-edge intersection
// Grandy, p.446
- const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
+ const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
+ //? const int numZeros = (epsilonEqual(cPQ,0.0) ? 1 : 0) + (epsilonEqual(cQR,0.0) ? 1 : 0) + (epsilonEqual(cRP, 0.0) ? 1 : 0);
if(numZeros >= 2 )
{
X, O, // OX
Y, O, // OY
Z, O, // OZ
- X, Y, // XY
+ X, Y, // XY
Y, Z, // YZ
Z, X, // ZX
};
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);
+ if(false)
+ //if(epsilonEqual(c1, 0.0, TEST_EPS) || epsilonEqual(c2, 0.0, TEST_EPS) || epsilonEqual(c3, 0.0, TEST_EPS))
+ {
+ return false;
+ }
+ else
+ {
+ return (c1*c3 > 0.0) && (c2*c3 > 0.0);
+ }
+
}
/**
//? should we use epsilon-equality here in second test?
// std::cout << "coord1 : " << coord1 << " coord2 : " << coord2 << std::endl;
+ //return (coord1*coord2 <= 0.0) && epsilonEqual(coord1,coord2);
return (coord1*coord2 <= 0.0) && (coord1 != coord2);
}
//? is it always YZ here ?
//? changed to XY !
const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY);
- // std::cout << "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) << "] - stable : " << calcStableT(corner) << std::endl;
+ //std::cout << "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) << std::endl;
+ //std::cout << "] - stable : " << calcStableT(corner) << std::endl;
+
//? 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;
- //return ( calcStableT(corner) * normal ) >= 0.0;
+ //return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0;
+ if(!_validTP[corner])
+ {
+ return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0;
+ }
+ else
+ {
+ return ( calcStableT(corner) * normal ) >= 0.0;
+ }
}
/**
#include "VectorUtils.hxx"
+#undef SECOND_CORNER_RESET
+#undef FIXED_DELTA 1.0e-15
+
+
namespace INTERP_UTILS
{
const double term3 = _doubleProducts[8*seg + C_XY] * _doubleProducts[8*seg + C_ZH];
// std::cout << std::endl;
- //std::cout << "for seg " << seg << " consistency " << term1 + term2 + term3 << std::endl;
- //std::cout << "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 << std::endl;
+ // std::cout << "for seg " << seg << " consistency " << term1 + term2 + term3 << std::endl;
+ // std::cout << "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 << std::endl;
// if(term1 + term2 + term3 != 0.0)
const int num_zero = (term1 == 0.0 ? 1 : 0) + (term2 == 0.0 ? 1 : 0) + (term3 == 0.0 ? 1 : 0);
// * two terms zero
// * all terms positive
// * all terms negative
- if((num_zero == 1 && num_neg != 1) || num_zero == 2 || num_neg == 0 || num_neg == 3 )
+ if((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 )
{
std::map<double, TetraCorner> distances;
//std::cout << "inconsistent! "<< std::endl;
// first element -> minimum distance
const TetraCorner minCorner = distances.begin()->second;
- // set the three corresponding double products to 0.0
- static const DoubleProduct DOUBLE_PRODUCTS[12] =
+ resetDoubleProducts(seg, minCorner);
+ // test : reset also for second corner if distance is small enough
+#ifdef SECOND_CORNER_RESET
+ std::map<double, TetraCorner>::const_iterator iter = distances.begin();
+ ++iter;
+ if(iter->first < 1.0e-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) {
- DoubleProduct dp = DOUBLE_PRODUCTS[3*minCorner + i];
-
- // std::cout << std::endl << "inconsistent dp :" << dp << std::endl;
- _doubleProducts[8*seg + dp] = 0.0;
- }
-
+ resetDoubleProducts(seg, iter->second);
+ }
+#endif
}
+
}
for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1))
{
- for(DoubleProduct dp = C_YZ ; dp <= C_10 ; dp = DoubleProduct(dp + 1))
+ for(DoubleProduct dp = C_YZ ; dp <= C_10 ; dp = DoubleProduct(dp + 1))
{
// find the points of the triangle
// 0 -> P, 1 -> Q, 2 -> R
const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2];
const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
-
+#ifdef FIXED_DELTA
+ const double delta = FIXED_DELTA;
+#else
const long double delta = MULT_PREC_F * ( std::abs(term1) + std::abs(term2) );
+#endif
- if( std::abs(_doubleProducts[8*seg + dp]) < THRESHOLD_F*delta )
+ if( epsilonEqual(_doubleProducts[8*seg + dp], 0.0, THRESHOLD_F * delta))
{
if(_doubleProducts[8*seg + dp] != 0.0)
{
- // std::cout << "Double product for (seg,dp) = (" << seg << ", " << dp << ") = "
- // << std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" << std::endl;
+ // std::cout << "Double product for (seg,dp) = (" << seg << ", " << dp << ") = "
+ // << std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" << std::endl;
}
_doubleProducts[8*seg + dp] = 0.0;
_isDoubleProductsCalculated = true;
}
+ 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];
+
+ // std::cout << std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner << std::endl;
+ _doubleProducts[8*seg + dp] = 0.0;
+ };
+ }
+
double TransformedTriangle::calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const
{
// NB uses fact that TriSegment <=> TriCorner that is first point of segment (PQ <=> P)
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, minRow, false);
+ //_tripleProducts[corner] = calcTByDevelopingRow(corner, 1, false);
}
_validTP[corner] = true;
}
1.0, 0.0, 0.0, // OX
0.0, 1.0, 0.0, // OY
0.0, 0.0, 1.0, // OZ
- -1.0,-1.0, 1.0, // XY
+ -1.0, 1.0, 0.0, // XY
0.0,-1.0, 1.0, // YZ
1.0, 0.0,-1.0 // ZX
};
EDGE_VECTORS[3*edge + 1],
EDGE_VECTORS[3*edge + 2],
};
-
+
+ //return angleBetweenVectors(normal, edgeVec);
+
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;
return 3.14159262535 - acos( dotProd / ( lenNormal * lenEdgeVec ) );
+
}
/**
const double sumDPProd = coordDPProd[0] + coordDPProd[1] + coordDPProd[2];
const double sumDPProdSq = dot(coordDPProd, coordDPProd);
- alpha = sumDPProd / sumDPProdSq;
+ // alpha = sumDPProd / sumDPProdSq;
+ alpha = (sumDPProdSq != 0.0) ? sumDPProd / sumDPProdSq : 0.0;
}
const double cQRbar = cQR * (1.0 - alpha * coordValues[0] * cQR);
const double cPQbar = cPQ * (1.0 - alpha * coordValues[2] * cPQ);
// check sign of barred products - should not change
- assert(cQRbar * cQR >= 0.0);
- assert(cRPbar * cRP >= 0.0);
- assert(cPQbar * cPQ >= 0.0);
+ // assert(cQRbar * cQR >= 0.0);
+ //assert(cRPbar * cRP >= 0.0);
+ //assert(cPQbar * cPQ >= 0.0);
const double p_term = _coords[5*P + offset] * cQRbar;
const double q_term = _coords[5*Q + offset] * cRPbar;
// check that we are not so close to zero that numerical errors could take over,
// otherwise reset to zero (cf Grandy, p. 446)
+#ifdef FIXED_DELTA
+ const double delta = FIXED_DELTA;
+#else
const long double delta = MULT_PREC_F * (std::abs(p_term) + std::abs(q_term) + std::abs(r_term));
-
+#endif
+
if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * delta) )
{
+ // std::cout << "Reset imprecise triple product for corner " << corner << " to zero" << std::endl;
return 0.0;
}
else
