}
+ // to be removed
assert(faceType == MED_TRIA3);
// create transformed triangles from face
case MED_TRIA3:
triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6]));
break;
+
+ // add other cases here to treat hexahedra, pyramides, etc
default:
std::cout << "Only elements with triangular faces are supported at the moment." << std::endl;
}
-
-#if 0
- void MeshElement::triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
- {
- std::cout << "Triangulating element " << getIndex() << std::endl;
- switch(_type)
- {
- case MED_TETRA4 :
- // triangles == faces
- const int beginFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index];
- const int endFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index + 1];
-
- std::cout << "elements has faces at indices " << beginFaceIdx << " to " << endFaceIdx << std::endl;
- assert(endFaceIdx - beginFaceIdx == 4);
-
- for(int i = beginFaceIdx ; i < endFaceIdx ; ++i) // loop over faces of element
- {
- const int faceIdx = _mesh->getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS)[i - 1];
- const int beginNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx - 1];
- const int endNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx];
- std::cout << "Face " << faceIdx << " with nodes in [" << beginNodeIdx << "," << endNodeIdx << "[" << std::endl;
- assert(endNodeIdx - beginNodeIdx == 3);
-
- double transformedPts[9];
-
- for(int j = 0 ; j < 3 ; ++j) // loop over nodes of face
- {
- // { optimise here using getCoordinatesForNode ?
- // could maybe use the connNodeIdx directly and only transform each node once
- // instead of once per face
- const int connNodeIdx =
- _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_FACE, MED_ALL_ELEMENTS)[beginNodeIdx + j - 1] - 1;
- const double* pt = &(_mesh->getCoordinates(MED_FULL_INTERLACE)[3*connNodeIdx]);
-
- //const double* pt = getCoordsOfNode(j + 1);
-
- // transform
- T.apply(&transformedPts[3*j], pt);
- std::cout << "Transforming : " << vToStr(pt) << " -> " << vToStr(&transformedPts[3*j]) << std::endl;
- }
-
- triangles.push_back(TransformedTriangle(&transformedPts[0], &transformedPts[3], &transformedPts[6]));
- }
-
- break;
- default:
- break;
- }
- }
-#endif
-
int MeshElement::getIndex() const
{
return _index + 1;
#include <fstream>
#include <cassert>
#include <cmath>
+#include "VectorUtils.hxx"
namespace INTERP_UTILS
{
const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
pt[i] = -( sign * calcStableC(seg, dp) ) / s;
- assert(pt[i] > 0.0);
- assert(pt[i] < 1.0);
+ // std::cout << "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] << std::endl;
+ assert(pt[i] >= 0.0);
+ assert(pt[i] <= 1.0);
}
}
// products 3 and 4 for each edge -> third condition
// NB : some uncertainty whether these last are correct
static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] =
+ {
+ C_10, C_01, C_ZH, C_10, // XY
+ C_01, C_XY, C_XH, C_01, // YZ
+ C_XY, C_10, C_YH, C_XY // ZX
+ };
+ /*static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] =
{
C_10, C_01, C_ZH, C_10, // XY
C_01, C_XY, C_XH, C_ZX, // YZ
C_XY, C_10, C_YH, C_XY // ZX
};
-
+ */
/*static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] =
{
C_10, C_01, C_ZH, C_XY, // XY
assert(pt[i] >= 0.0);
assert(pt[i] <= 1.0);
}
- assert(std::abs(pt[0] + pt[1] + pt[2] - 1.0) < 1.0e-9);
+ assert(epsilonEqual(pt[0] + pt[1] + pt[2] - 1.0, 0.0, 1.0e-9));
}
/**
bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const
{
assert(corner == X || corner == Y || corner == Z);
+ // std::cout << "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner << std::endl;
// readjust index since O is not used
const int cornerIdx = static_cast<int>(corner) - 1;
// 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] =
+ {
+ 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
+ };*/
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
};
-
+ */
// facets to use
//? not sure this is correct
- static const TetraFacet FACET_SEGMENT_RAY_INTERSECTION[3] =
+ static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
{
OZX, // X
- OXY, // Y
- OYZ, // Z
+ OYZ, // Y
+ OZX, // Z
};
const DoubleProduct dp0 = DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx];
+ const double cVal0 = calcStableC(seg, dp0);
//? epsilon-equality here?
- if(dp0 == 0.0) // cond. 1
+ // cond. 1
+ if(cVal0 != 0.0)
{
+ // std::cout << "SR fails at cond 1 cVal0 = " << cVal0 << std::endl;
+ return false;
+ }
- if(testSegmentIntersectsFacet(seg, FACET_SEGMENT_RAY_INTERSECTION[cornerIdx])) // cond. 2.1
- {
- const double H1 = _coords[5*seg + 4];
- const double H2 = _coords[5*( (seg + 1) % 3) + 4];
-
- // S_H -> cond. 2.2
- // std::cout << "H1 = " << H1 << ", H2= " << H2 << std::endl;
- if(H1*H2 <= 0.0 && H1 != H2) // should equality be in "epsilon-sense" ?
- {
+ // cond 2
+ const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
+ const bool cond22 = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
- const double cVals[6] =
- {
- calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]),
- calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]),
- calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]),
- calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]),
- calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]),
- calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]),
- };
-
- // cond. 3
- return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
- }
- }
+ if(!(cond21 || cond22))
+ {
+ // std::cout << "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 << std::endl;
+ return false;
}
- return false;
+ // cond 3
+ const double cVals[6] =
+ {
+ calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]),
+ calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]),
+ calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]),
+ calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]),
+ calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]),
+ calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]),
+ };
+
+ // cond. 3
+ if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0)
+ {
+ std::cout << "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] << std::endl;
+ }
+ return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
+
}
/**
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;
}
// 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;
//? 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, 2, false) * normal ) >= 0.0;
- return ( calcStableT(corner) * normal ) >= 0.0;
+ return ( calcTByDevelopingRow(corner, 2, false) * normal ) >= 0.0;
+ //return ( calcStableT(corner) * normal ) >= 0.0;
}
/**
#include <cassert>
#include <cmath>
#include <limits>
+#include "VectorUtils.hxx"
namespace INTERP_UTILS
{
//{ use functions in VectorUtils for this
// cross product of difference vectors
- const double cross[3] =
- {
- diffPQ[1]*diffCornerP[2] - diffPQ[2]*diffCornerP[1],
- diffPQ[2]*diffCornerP[0] - diffPQ[0]*diffCornerP[2],
- diffPQ[0]*diffCornerP[1] - diffPQ[1]*diffCornerP[0],
- };
-
+
+ double crossProd[3];
+ cross(diffPQ, diffCornerP, crossProd);
- const double cross_squared = cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2];
- const double norm_diffPQ_squared = diffPQ[0]*diffPQ[0] + diffPQ[1]*diffPQ[1] + diffPQ[2]*diffPQ[2];
+ const double cross_squared = dot(crossProd, crossProd);
+ const double norm_diffPQ_squared = dot(diffPQ, diffPQ);
const double dist = cross_squared / norm_diffPQ_squared;
// update minimum (initializs with first corner)
double TransformedTriangle::calcStableT(const TetraCorner corner) const
{
assert(_isTripleProductsCalculated);
- // assert(_validTP[corner]);
+ assert(_validTP[corner]);
return _tripleProducts[corner];
}
const double cPQ = calcStableC(PQ, dp);
// coordinate to use for projection (Grandy, [57]) with edges
- // OX, OY, OZ, YZ, ZX, XY in order :
+ // OX, OY, OZ, XY, YZ, ZX in order :
// (y, z, x, h, h, h)
// for the first three we could also use {2, 0, 1}
static const int PROJECTION_COORDS[6] = { 1, 2, 0, 3, 3, 3 } ;
const double p_term = _coords[5*P + offset] * cQR * (1.0 - alpha * coordValues[0] * cQR) ;
const double q_term = _coords[5*Q + offset] * cRP * (1.0 - alpha * coordValues[1] * cRP) ;
const double r_term = _coords[5*R + offset] * cPQ * (1.0 - alpha * coordValues[2] * cPQ) ;
-
- // NB : using plus also for the middle term compensates for a double product
- // which is inversely ordered
- // std::cout << "Triple product for corner " << corner << ", row " << row << " = " << sign*( p_term + q_term + r_term ) << std::endl;
- return sign*( p_term + q_term + r_term );
+
+ // check that we are not so close to zero that numerical errors could take over,
+ // otherwise reset to zero (cf Grandy, p. 446)
+ const long double delta = MULT_PREC_F * (std::abs(p_term) + std::abs(q_term) + std::abs(r_term));
+ if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * delta) )
+ {
+ return 0.0;
+ }
+ else
+ {
+ // NB : using plus also for the middle term compensates for a double product
+ // which is inversely ordered
+ // std::cout << "Triple product for corner " << corner << ", row " << row << " = " << sign*( p_term + q_term + r_term ) << std::endl;
+ return sign*( p_term + q_term + r_term );
+ }
}
