T.dump();
std::cout << std::endl;
+ if(T.determinant() == 0.0)
+ {
+ // tetra is planar
+ std::cout << "Planar tetra -- volume 0" << std::endl;
+ return 0.0;
+ }
+
// triangulate source element faces (assumed tetraeder for the time being)
// do nothing
vector<TransformedTriangle> triangles;
#include "TetraAffineTransform.hxx"
#include "TransformedTriangle.hxx"
+#include "MEDMEM_CellModel.hxx"
namespace INTERP_UTILS
{
* @param T affine transform that is applied to the nodes of the triangles
*/
void MeshElement::triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
+ {
+ // get cell model for the element
+ CELLMODEL cellModel(_type);
+
+ assert(cellModel.getDimension() == 3);
+
+ // start index in connectivity array for cell
+ // const int cellIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1;
+
+ // assert(cellIdx >= 0);
+ //assert(cellIdx < _mesh->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
+
+ // loop over faces
+ for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i)
+ {
+ medGeometryElement faceType = cellModel.getConstituentType(1, i);
+ CELLMODEL faceModel(faceType);
+
+ assert(faceModel.getDimension() == 2);
+
+ double transformedNodes[3 * faceModel.getNumberOfNodes()];
+ const double* coords = _mesh->getCoordinates(MED_FULL_INTERLACE);
+
+ // 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);
+
+ // transform
+ //{ not totally efficient since we transform each node once per face
+ T.apply(&transformedNodes[3*(j-1)], node);
+
+ }
+
+ assert(faceType == MED_TRIA3);
+
+ // create triangles from face
+ switch(faceType)
+ {
+ case MED_TRIA3:
+ triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6]));
+ break;
+
+ 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)
break;
}
}
-
+#endif
int MeshElement::getIndex() const
{
return _index + 1;
_linearTransform[3*j + i] = (pts[i+1])[j] - (pts[0])[j];
}
}
-#if 0
+
calculateDeterminant();
+
+ std::cout << "determinant before inverse = " << _determinant << std::endl;
+
+ // check that tetra is non-planar -> determinant is not zero
+ // otherwise set _determinant to zero to signal caller that transformation did not work
+ if(epsilonEqual(_determinant, 0.0))
+ {
+ _determinant = 0.0;
+ return;
+ }
+
+#if 0
+
//assert(_determinant != 0.0);
if(_determinant < 0.0)
// precalculate determinant (again after inversion of transform)
calculateDeterminant();
+ // self-check
+ std::cout << "transform determinant is " << _determinant << std::endl;
std::cout << "*Self-check : Applying transformation to original points ... ";
for(int i = 0; i < 4 ; ++i)
{
double v[3];
apply(v, pts[i]);
- // std::cout << vToStr(v) << std::endl;
+ std::cout << vToStr(v) << std::endl;
for(int j = 0; j < 3; ++j)
{
assert(epsilonEqual(v[j], (3*i+j == 3 || 3*i+j == 7 || 3*i+j == 11 ) ? 1.0 : 0.0));
// calculate point
for(int i = 0; i < 3; ++i)
{
- // std::cout << "alpha= " << alpha << std::endl;
+ 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;
- assert(pt[i] >= 0.0);
- assert(pt[i] <= 1.0);
+ std::cout << pt[i] << std::endl;
+ //assert(pt[i] >= 0.0);
+ //assert(pt[i] <= 1.0);
}
}
const double cQR = calcStableC(QR, DoubleProduct(edge));
const double cRP = calcStableC(RP, DoubleProduct(edge));
- // std::cout << "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl;
+ std::cout << "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl;
// if two or more c-values are zero we disallow x-edge intersection
// Grandy, p.446
X, O, // OX
Y, O, // OY
Z, O, // OZ
+ X, Y, // XY
Y, Z, // YZ
Z, X, // ZX
- X, Y // XY
};
// Grandy, [16]
const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]);
//? should equality with zero use epsilon?
+ std::cout << "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 << std::endl;
return (t1*t2 <= 0.0) && (t1 - t2 != 0.0);
}
//? 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) << "] - 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;
}
/**
if(_isTripleProductsCalculated)
return;
- std::cout << "Precalculating triple products" << std::endl;
+ // std::cout << "Precalculating triple products" << std::endl;
for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1))
{
- std::cout << "- Triple product for corner " << corner << std::endl;
+ // std::cout << "- Triple product for corner " << corner << std::endl;
bool isGoodRowFound = false;
// find edge / row to use
{
// this value will not be used
// we set it to whatever
- // std::cout << "Triple product not calculated for corner " << corner << std::endl;
+ std::cout << "Triple product not calculated for corner " << corner << std::endl;
_tripleProducts[corner] = -3.14159265;
_validTP[corner] = false;
