From 69b57141e40aaea3563b30ccc8acd034f0e76d9e Mon Sep 17 00:00:00 2001 From: vbd Date: Mon, 10 Sep 2007 10:46:09 +0000 Subject: [PATCH] staffan : * added class IntersectorTetra which treats intersection between one tetrahedron and source elements in a more effective manner, by caching intermediary results --- src/INTERP_KERNEL/IntersectorTetra.cxx | 240 ++++++++++++++++++ src/INTERP_KERNEL/IntersectorTetra.hxx | 181 +++++++++++++ src/INTERP_KERNEL/Makefile.in | 9 +- src/INTERP_KERNEL/MeshUtils.hxx | 11 + .../Test/Interpolation3DTest.cxx | 7 +- .../Test/Interpolation3DTest.hxx | 22 +- src/INTERP_KERNEL/Test/Makefile.in | 38 +-- src/INTERP_KERNEL/Test/PerfTest.cxx | 2 +- src/INTERP_KERNEL/TransformedTriangle.hxx | 4 + .../TransformedTriangle_intersect.cxx | 22 +- .../TransformedTriangle_math.cxx | 2 +- 11 files changed, 486 insertions(+), 52 deletions(-) create mode 100644 src/INTERP_KERNEL/IntersectorTetra.cxx create mode 100644 src/INTERP_KERNEL/IntersectorTetra.hxx diff --git a/src/INTERP_KERNEL/IntersectorTetra.cxx b/src/INTERP_KERNEL/IntersectorTetra.cxx new file mode 100644 index 000000000..ac67374a3 --- /dev/null +++ b/src/INTERP_KERNEL/IntersectorTetra.cxx @@ -0,0 +1,240 @@ +#include "IntersectorTetra.hxx" + +#include "TetraAffineTransform.hxx" +#include "TransformedTriangle.hxx" +#include "MeshUtils.hxx" +#include "VectorUtils.hxx" +#include "Log.hxx" + +#include +#include +#include +#include + +using namespace MEDMEM; +using namespace MED_EN; +using namespace INTERP_UTILS; + +namespace MEDMEM +{ + /* + * Constructor + * + * @param srcMesh mesh containing the source elements + * @param targetMesh mesh containing the target elements + * + */ + IntersectorTetra::IntersectorTetra(const MESH& srcMesh, const MESH& targetMesh, int targetCell) + : _srcMesh(srcMesh), filtered(0) + { + const medGeometryElement targetType = targetMesh.getElementType(MED_CELL, targetCell); + + // maybe we should do something more civilized here + assert(targetType == MED_TETRA4); + + // get array of points of target tetraeder + const double* tetraCorners[4]; + for(int i = 0 ; i < 4 ; ++i) + { + tetraCorners[i] = getCoordsOfNode(i + 1, targetCell, targetMesh); + } + + // create AffineTransform from tetrahedron + _t = new TetraAffineTransform( tetraCorners ); + + } + + /* + * Destructor + * + */ + IntersectorTetra::~IntersectorTetra() + { + delete _t; + + for(hash_map< int, double*>::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter) + { + delete[] iter->second; + } + } + + /* + * Calculates the volume of intersection of an element in the source mesh and an element + * in the target mesh. The method is based on the algorithm of Grandy. It first calculates the transformation + * that takes the target tetrahedron into the unit tetrahedron. After that, the + * faces of the source element are triangulated and the calculated transformation is applied + * to each triangle. The algorithm of Grandy, implemented in INTERP_UTILS::TransformedTriangle is used + * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions + * divided by the determinant of the transformation. + * + * @pre The element in _targetMesh referenced by targetCell is of type MED_TETRA4. + * @param srcCell global number of the source element (1 <= srcCell < # source cells) + * @param targetCell global number of the target element (1 <= targetCell < # target cells) - this element must be a tetrahedron + */ + double IntersectorTetra::intersectSourceCell(int element) + { + + //{ could be done on outside + // check if we have planar tetra element + if(_t->determinant() == 0.0) + { + // tetra is planar + LOG(2, "Planar tetra -- volume 0"); + return 0.0; + } + + // get type of cell + const medGeometryElement type = _srcMesh.getElementType(MED_CELL, element); + + // get cell model for the element + const CELLMODEL cellModel(type); + + assert(cellModel.getDimension() == 3); + assert(element >= 1); + assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS)); + + // halfspace filtering + bool isOutside[8] = {true, true, true, true, true, true, true, true}; + bool isTargetOutside = false; + + // calculate the coordinates of the nodes + for(int i = 1; i <= cellModel.getNumberOfNodes() ; ++i) + { + // we could store mapping local -> global numbers too, but not sure it is worth it + const int globalNodeNum = getGlobalNumberOfNode(i, element, _srcMesh); + if(_nodes.find(globalNodeNum) == _nodes.end()) + { + calculateNode(globalNodeNum); + } + + checkIsOutside(_nodes[globalNodeNum], isOutside); + + // local caching of globalNodeNum + // not sure this is efficient + // globalNodeNumbers[i] = globalNodeNum; + } + + // halfspace filtering check + // NB : might not be beneficial for caching of triangles + for(int i = 0; i < 8; ++i) + { + if(isOutside[i]) + { + isTargetOutside = true; + } + } + + double totalVolume = 0.0; + + if(!isTargetOutside) + { + for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i) + { + const medGeometryElement faceType = cellModel.getConstituentType(1, i); + const CELLMODEL faceModel(faceType); + + assert(faceModel.getDimension() == 2); + + int faceNodes[faceModel.getNumberOfNodes()]; + + // get the nodes of the face + for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) + { + const int locNodeNum = cellModel.getNodeConstituent(1, i, j); + assert(locNodeNum >= 1); + assert(locNodeNum <= cellModel.getNumberOfNodes()); + + faceNodes[j-1] = getGlobalNumberOfNode(locNodeNum, element, _srcMesh);//globalNodeNumbers[locNodeNum]; + } + + switch(faceType) + { + case MED_TRIA3: + { + // create the face key + TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); + + // calculate the triangle if needed + if(_volumes.find(key) == _volumes.end()) + { + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); + calculateVolume(tri, key); + totalVolume += _volumes[key]; + } else { + // count negative as face has reversed orientation + totalVolume -= _volumes[key]; + } + } + break; + + case MED_QUAD4: + + // simple triangulation of faces along a diagonal : + // + // 2 ------ 3 + // | / | + // | / | + // | / | + // | / | + // | / | + // | / | + // 1 ------ 4 + // + //? not sure if this always works + { + // calculate the triangles if needed + + // local nodes 1, 2, 3 + TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); + if(_volumes.find(key1) == _volumes.end()) + { + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); + calculateVolume(tri, key1); + totalVolume += _volumes[key1]; + } else { + // count negative as face has reversed orientation + totalVolume -= _volumes[key1]; + } + + // local nodes 1, 3, 4 + TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]); + if(_volumes.find(key2) == _volumes.end()) + { + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]); + calculateVolume(tri, key2); + totalVolume += _volumes[key2]; + } + else + { + // count negative as face has reversed orientation + totalVolume -= _volumes[key2]; + } + } + break; + + default: + std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl; + assert(false); + } + } + } + else + { + ++filtered; + } + + // reset if it is very small to keep the matrix sparse + // is this a good idea? + if(epsilonEqual(totalVolume, 0.0, 1.0e-11)) + { + totalVolume = 0.0; + } + + LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant()); + + // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation + // that should be used (which is equivalent to dividing by the determinant) + return std::abs(1.0 / _t->determinant() * totalVolume) ; + } + +}; diff --git a/src/INTERP_KERNEL/IntersectorTetra.hxx b/src/INTERP_KERNEL/IntersectorTetra.hxx new file mode 100644 index 000000000..acb7ed885 --- /dev/null +++ b/src/INTERP_KERNEL/IntersectorTetra.hxx @@ -0,0 +1,181 @@ +#ifndef __INTERSECTOR_TETRA_HXX__ +#define __INTERSECTOR_TETRA_HXX__ + +#include "MEDMEM_define.hxx" +#include "MEDMEM_Mesh.hxx" + +#include "Intersector.hxx" +#include +#include +#include + +#include "TetraAffineTransform.hxx" +#include "TransformedTriangle.hxx" + +using __gnu_cxx::hash_map; + +namespace INTERP_UTILS +{ + class TriangleFaceKey + { + public: + int _nodes[3]; + int _hashVal; + + TriangleFaceKey(int node1, int node2, int node3) + { + sort3Ints(_nodes, node1, node2, node3); + // assert(_nodes[0] < _nodes[1]); + //assert(_nodes[0] < _nodes[2]); + //assert(_nodes[1] < _nodes[2]); + _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29; + } + + bool operator==(const TriangleFaceKey& rhs) const + { + return _nodes[0] == rhs._nodes[0] && _nodes[1] == rhs._nodes[1] && _nodes[2] == rhs._nodes[2]; + } + + int hashVal() const + { + return _hashVal; + } + + inline void sort3Ints(int* sorted, int node1, int node2, int node3); + }; +} + +namespace __gnu_cxx +{ + template<> + class hash + { + public: + int operator()(const INTERP_UTILS::TriangleFaceKey& key) const + { + return key.hashVal(); + } + }; +}; + +namespace INTERP_UTILS +{ + + + /* + * Class calculating the volume of intersection of individual 3D elements. + * + */ + class IntersectorTetra + { + + public : + IntersectorTetra(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, int targetCell); + + ~IntersectorTetra(); + + double intersectSourceCell(int srcCell); + + mutable int filtered; + + private: + + TetraAffineTransform* _t; + + hash_map< int, double*> _nodes; + hash_map< TriangleFaceKey, double > _volumes; + + inline void checkIsOutside(const double* pt, bool* isOutside) const; + inline void calculateNode(int globalNodeNum); + inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); + // inline void sort3Ints(int* sorted, int node1, int node2, int node3); + //inline std::string createFaceKey(int node1, int node2, int node3); + + const MEDMEM::MESH& _srcMesh; + + }; + + inline void IntersectorTetra::checkIsOutside(const double* pt, bool* isOutside) const + { + isOutside[0] = isOutside[0] && (pt[0] <= 0.0); + isOutside[1] = isOutside[1] && (pt[0] >= 1.0); + isOutside[2] = isOutside[2] && (pt[1] <= 0.0); + isOutside[3] = isOutside[3] && (pt[1] >= 1.0); + isOutside[4] = isOutside[4] && (pt[2] <= 0.0); + isOutside[5] = isOutside[5] && (pt[2] >= 1.0); + isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] <= 0.0); + isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0); + } + + inline void IntersectorTetra::calculateNode(int globalNodeNum) + { + assert(globalNodeNum >= 0); + assert(globalNodeNum < _srcMesh.getNumberOfNodes()); + const double* node = &(_srcMesh.getCoordinates(MED_EN::MED_FULL_INTERLACE)[3*globalNodeNum]); + double* transformedNode = new double[3]; + assert(transformedNode != 0); + + _t->apply(transformedNode, node); + _nodes[globalNodeNum] = transformedNode; + } + + inline void IntersectorTetra::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key) + { + const double vol = tri.calculateIntersectionVolume(); + _volumes.insert(make_pair(key, vol)); + } + + inline void TriangleFaceKey::sort3Ints(int* sorted, int node1, int node2, int node3) + { + if(node1 < node2) + { + if(node1 < node3) + { + // node 1 is min + sorted[0] = node1; + sorted[1] = node2 < node3 ? node2 : node3; + sorted[2] = node2 < node3 ? node3 : node2; + } + else + { + // node3 , node1, node2 + sorted[0] = node3; + sorted[1] = node1; + sorted[2] = node2; + } + } + else // x2 < x1 + { + if(node2 < node3) + { + // node 2 is min + sorted[0] = node2; + sorted[1] = node1 < node3 ? node1 : node3; + sorted[2] = node1 < node3 ? node3 : node1; + } + else + { + // node 3, node 2, node 1 + sorted[0] = node3; + sorted[1] = node2; + sorted[2] = node1; + } + } + } + +#if 0 + inline std::string IntersectorTetra::createFaceKey(int node1, int node2, int node3) + { + int sorted[3]; + sort3Ints(sorted, node1, node2, node3); + + std::stringstream sstr; + sstr << node1 << "-" << node2 << "-" << node3; + return sstr.str(); + } +#endif + +}; + + +#endif diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index 306be5718..7550eb131 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -54,7 +54,8 @@ BBTree.H\ MeshUtils.hxx\ Intersector3D.hxx\ Log.hxx\ -TransformedTriangle_inline.hxx +TransformedTriangle_inline.hxx\ +IntersectorTetra.hxx # Libraries targets @@ -73,7 +74,8 @@ TetraAffineTransform.cxx\ Interpolation3D.cxx\ Intersector3D.cxx\ Interpolation3DSurf.cxx\ -Interpolation2D.cxx +Interpolation2D.cxx\ +IntersectorTetra.cxx # Executables targets BIN = @@ -90,8 +92,7 @@ CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) CPPFLAGS+=$(BOOST_CPPFLAGS) # optimization -#CXXFLAGS+= -O1 -DOPTIMIZE -#CPPFLAGS+= -O1 -DOPTIMIZE + CXXFLAGS+= -O2 -DOPTIMIZE_FILTER -DOPTIMIZE CPPFLAGS+= -O2 -DOPTIMIZE_FILTER -DOPTIMIZE diff --git a/src/INTERP_KERNEL/MeshUtils.hxx b/src/INTERP_KERNEL/MeshUtils.hxx index 498ea0aa5..ae8d81a97 100644 --- a/src/INTERP_KERNEL/MeshUtils.hxx +++ b/src/INTERP_KERNEL/MeshUtils.hxx @@ -4,6 +4,8 @@ #include "MEDMEM_define.hxx" #include +#include + using namespace MEDMEM; using namespace MED_EN; @@ -45,6 +47,15 @@ namespace INTERP_UTILS return static_cast(type) - 300; } + inline int getGlobalNumberOfNode(int node, int element, const MESH& mesh) + { + assert(node >= 1); + assert(node <= mesh.getNumberOfNodes()); + const int nodeOffset = node - 1; + const int elemIdx = mesh.getConnectivityIndex(MED_NODAL, MED_CELL)[element - 1] - 1; + return mesh.getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_CELL, MED_ALL_ELEMENTS)[elemIdx + nodeOffset] - 1; + } + }; diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx index cdb5dcc23..6429f864f 100644 --- a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx +++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx @@ -259,7 +259,12 @@ void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const ch m = interpolator->interpol_maillages(sMesh, tMesh); - testVolumes(m, sMesh, tMesh); + // if reflexive, check volumes + if(strcmp(mesh1path,mesh2path) == 0) + { + const bool row_and_col_sums_ok = testVolumes(m, sMesh, tMesh); + CPPUNIT_ASSERT_EQUAL_MESSAGE("Row or column sums incorrect", true, row_and_col_sums_ok); + } LOG(1, "Intersection calculation done. " << std::endl ); diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx index 767488265..b52bb16c4 100644 --- a/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx +++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx @@ -36,24 +36,26 @@ class Interpolation3DTest : public CppUnit::TestFixture // multi - element CPPUNIT_TEST( tetraComplexIncluded ); -#endif + CPPUNIT_TEST( dividedUnitTetraSimplerReflexive ); CPPUNIT_TEST( dividedUnitTetraReflexive ); - //#if 0 - //CPPUNIT_TEST( nudgedDividedUnitTetra ); - //CPPUNIT_TEST( nudgedDividedUnitTetraSimpler ); - //CPPUNIT_TEST( dividedGenTetra ); + + CPPUNIT_TEST( nudgedDividedUnitTetra ); + CPPUNIT_TEST( nudgedDividedUnitTetraSimpler ); + CPPUNIT_TEST( dividedGenTetra ); CPPUNIT_TEST( boxReflexive ); CPPUNIT_TEST( boxReflexiveModerate ); - //CPPUNIT_TEST( tetraBoxes ); - //CPPUNIT_TEST( moderateBoxes ); - //CPPUNIT_TEST( moderateBoxesSmaller ); + CPPUNIT_TEST( tetraBoxes ); + CPPUNIT_TEST( moderateBoxes ); +#endif + CPPUNIT_TEST( moderateBoxesSmaller ); CPPUNIT_TEST( moderateBoxSmallReflexive ); +#if 0 CPPUNIT_TEST( moderateBoxEvenSmallerReflexive ); CPPUNIT_TEST( tinyBoxReflexive ); - //CPPUNIT_TEST( simpleHexaBox ); - //#endif + CPPUNIT_TEST( simpleHexaBox ); +#endif CPPUNIT_TEST_SUITE_END(); diff --git a/src/INTERP_KERNEL/Test/Makefile.in b/src/INTERP_KERNEL/Test/Makefile.in index 94e48888a..c3fdbfae9 100644 --- a/src/INTERP_KERNEL/Test/Makefile.in +++ b/src/INTERP_KERNEL/Test/Makefile.in @@ -34,7 +34,11 @@ VPATH=.:@srcdir@:@top_srcdir@/idl @COMMENCE@ # header files -EXPORT_HEADERS = TestingUtils.hxx +EXPORT_HEADERS = CppUnitTest.hxx \ + TransformedTriangleTest.hxx \ + TransformedTriangleIntersectTest.hxx \ + Interpolation3DTest.hxx \ + TestingUtils.hxx # Libraries targets @@ -46,9 +50,10 @@ LIB_CLIENT_IDL = # Executables targets -BIN = PerfTest +BIN = PerfTest TestInterpKernel -BIN_SRC = +BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx \ +Interpolation3DTest.cxx BIN_CLIENT_IDL = @@ -60,20 +65,20 @@ CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/sr CPPFLAGS+=$(BOOST_CPPFLAGS) $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/src/INTERP_KERNEL #include CppUnit -#CXXFLAGS+= -I/usr/include/cppunit -#CPPFLAGS+= -I/usr/include/cppunit +CXXFLAGS+= -I/usr/include/cppunit +CPPFLAGS+= -I/usr/include/cppunit # for log -CXXFLAGS+= -DLOG_LEVEL=0 -CPPFLAGS+= -DLOG_LEVEL=0 +CXXFLAGS+= -DLOG_LEVEL=3 #-DOPTIMIZE -O2 +CPPFLAGS+= -DLOG_LEVEL=3 #-DOPTIMIZE -O2 # for gcov #CXXFLAGS+=-fprofile-arcs -ftest-coverage #CPPFLAGS+=-fprofile-arcs -ftest-coverage #for gprof -CXXFLAGS+=-pg -CPPFLAGS+=-pg +#CXXFLAGS+=-pg +#CPPFLAGS+=-pg #LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. @@ -82,7 +87,7 @@ LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil # for gcov #LDFLAGS+= -lgcov -LDFLAGS+= -pg +#LDFLAGS+= -pg #LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. @@ -99,15 +104,12 @@ endif LIBS = @LIBS@ @CPPUNIT_LIBS@ -LDFLAGSFORBIN += $(LDFLAGS) -lm $(MED2_LIBS) $(HDF5_LIBS) \ - -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome -lmed_V2_1 -lmedmem \ - -linterpkernel -#LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \ -# -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \ -# -lcppunit -linterpkernel +LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \ + -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \ + -lcppunit -linterpkernel -LDFLAGSFORBIN += -pg +#LDFLAGSFORBIN += -pg -UNIT_TEST_PROG = #PerfTest +UNIT_TEST_PROG = PerfTest TestInterpKernel @CONCLUDE@ diff --git a/src/INTERP_KERNEL/Test/PerfTest.cxx b/src/INTERP_KERNEL/Test/PerfTest.cxx index 958152527..f72c0e883 100644 --- a/src/INTERP_KERNEL/Test/PerfTest.cxx +++ b/src/INTERP_KERNEL/Test/PerfTest.cxx @@ -18,7 +18,7 @@ int main(int argc, char** argv) calcIntersectionMatrix(mesh1path.c_str(), mesh1.c_str(), mesh2path.c_str(), mesh2.c_str(), m); - //dumpIntersectionMatrix(m); + // dumpIntersectionMatrix(m); return 0; diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 97033e7a5..a313814c2 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -21,6 +21,8 @@ class TransformedTriangleTest; class TransformedTriangleIntersectTest; + + namespace INTERP_UTILS { @@ -320,4 +322,6 @@ namespace INTERP_UTILS }; +#include "TransformedTriangle_tables.hxx" + #endif diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index de37ecc49..5ab7a5f81 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -23,23 +23,6 @@ namespace INTERP_UTILS C_ZH, C_ZX, C_YZ, // OXY C_XH, C_YH, C_ZH // XYZ }; -#if 0 - template - inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection() - { - return NO_DP; - } - - template<> - inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection - { - return C_XH; - } - -#define DEF_DP_FOR_SEG_FACET(FACET,SEG,DP) template<> inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection { return DP; } - DEF_DP_FOR_SEG_FACET() - -#endif // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] = @@ -181,8 +164,13 @@ namespace INTERP_UTILS 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]; +#if 0 + pt[i] = (1 - alpha) * getCoordinateForTetCorner() + + alpha * getCoordinateForTetCorner(); +#endif LOG(6, pt[i] ); assert(pt[i] >= 0.0); assert(pt[i] <= 1.0); diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index 3e3bbcf63..3ef665293 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -321,7 +321,7 @@ namespace INTERP_UTILS _isTripleProductsCalculated = true; } - /* + /** * Calculates the angle between an edge of the tetrahedron and the triangle * * @param edge edge of the tetrahedron -- 2.39.2