From: vbd Date: Thu, 2 Aug 2007 14:40:40 +0000 (+0000) Subject: try to put my new TransformedTriangle_intersect.cxx file in the repos X-Git-Tag: trio_trio_coupling~88 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=7e7774fbfecb4fc3e069bb8623d21b406b0694f6;p=tools%2Fmedcoupling.git try to put my new TransformedTriangle_intersect.cxx file in the repos --- diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index 69f55d794..049b4c9a6 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -78,6 +78,8 @@ namespace MEDMEM const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); + std::cout << "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements " << std::endl; + // create empty maps for all source elements matrix.resize(numSrcElems); @@ -136,11 +138,11 @@ namespace MEDMEM { RegionNode* currNode = nodes.top(); nodes.pop(); - // std::cout << "Popping node " << std::endl; + std::cout << "Popping node " << std::endl; if(currNode->getSrcRegion().getNumberOfElements() == 1) { - // std::cout << " - One element" << std::endl; + std::cout << " - One element" << std::endl; // volume calculation MeshElement* srcElement = *(currNode->getSrcRegion().getBeginElements()); @@ -155,7 +157,7 @@ namespace MEDMEM iter != currNode->getTargetRegion().getEndElements() ; ++iter) { const double vol = calculateIntersectionVolume(*srcElement, **iter); - if(vol != 0.0) + // if(vol != 0.0) { const int targetIdx = (*iter)->getIndex(); @@ -167,7 +169,7 @@ namespace MEDMEM else // recursion { - // std::cout << " - Recursion" << std::endl; + std::cout << " - Recursion" << std::endl; RegionNode* leftNode = new RegionNode(); RegionNode* rightNode = new RegionNode(); @@ -178,19 +180,22 @@ namespace MEDMEM currNode->getSrcRegion().split(leftNode->getSrcRegion(), rightNode->getSrcRegion(), axis); + std::cout << "After split, left src region has " << leftNode->getSrcRegion().getNumberOfElements() << + " elements and right src region has " << rightNode->getSrcRegion().getNumberOfElements() << " elements" << std::endl; + // ugly hack to avoid problem with enum which does not start at 0 // I guess I ought to implement ++ for it instead ... // Anyway, it basically chooses the next axis, circually axis = (axis != BoundingBox::ZMAX) ? static_cast(axis + 1) : BoundingBox::XMAX; // add target elements of curr node that overlap the two new nodes - // std::cout << " -- Adding target elements" << std::endl; + std::cout << " -- Adding target elements" << std::endl; int numLeftElements = 0; int numRightElements = 0; for(vector::const_iterator iter = currNode->getTargetRegion().getBeginElements() ; iter != currNode->getTargetRegion().getEndElements() ; ++iter) { - // std::cout << " --- New target node" << std::endl; + //std::cout << " --- New target node" << std::endl; if(!leftNode->getSrcRegion().isDisjointWithElementBoundingBox(**iter)) { @@ -204,8 +209,11 @@ namespace MEDMEM ++numRightElements; } + } + std::cout << "Left target region has " << numLeftElements << " elements and right target region has " << numRightElements << " elements" << std::endl; + if(numLeftElements != 0) { nodes.push(leftNode); @@ -228,7 +236,7 @@ namespace MEDMEM } delete currNode; - // std::cout << "Next iteration. Nodes left : " << nodes.size() << std::endl; + std::cout << "Next iteration. Nodes left : " << nodes.size() << std::endl; } @@ -286,7 +294,9 @@ namespace MEDMEM // create AffineTransform TetraAffineTransform T( tetraCorners ); - // T.dump(); + std::cout << "Transform : " << std::endl; + T.dump(); + std::cout << std::endl; // triangulate source element faces (assumed tetraeder for the time being) // do nothing @@ -296,10 +306,10 @@ namespace MEDMEM double volume = 0.0; // std::cout << "num triangles = " << triangles.size() << std::endl; - + int i = 0; for(vector::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter) { - // std::cout << std::endl << "= > Triangle " << ++i << std::endl; + std::cout << std::endl << "= > Triangle " << ++i << std::endl; iter->dumpCoords(); volume += iter->calculateIntersectionVolume(); } diff --git a/src/INTERP_KERNEL/Interpolation3D.hxx b/src/INTERP_KERNEL/Interpolation3D.hxx index 40246802e..eaf5ece8b 100644 --- a/src/INTERP_KERNEL/Interpolation3D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D.hxx @@ -14,17 +14,14 @@ // typedefs typedef std::vector< std::map< int, double > > IntersectionMatrix; - -#ifdef TESTING_INTERP_KERNEL -class Interpolation3DTest; -#endif - namespace INTERP_UTILS { class MeshElement; class MeshRegion; }; +class Interpolation3DTest; + namespace MEDMEM { /** @@ -34,9 +31,8 @@ namespace MEDMEM class Interpolation3D : public Interpolation { -#ifdef TESTING_INTERP_KERNEL friend class ::Interpolation3DTest; -#endif + public : diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index 17d330ba9..394aa7cc2 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -39,11 +39,17 @@ EXPORT_PYSCRIPTS = \ EXPORT_HEADERS = \ +Interpolation.hxx\ +Interpolation3D.hxx\ +RegionNode.hxx\ +MeshRegion.hxx\ +MeshElement.hxx\ +BoundingBox.hxx\ +TetraAffineTransform.hxx\ +TranslationRotationMatrix.hxx\ Interpolation2D.hxx\ Interpolation3DSurf.hxx\ -Interpolation.hxx\ InterpolationUtils.hxx\ -TranslationRotationMatrix.hxx\ BBTree.H # Libraries targets @@ -54,8 +60,14 @@ LIB_SRC = \ TransformedTriangle.cxx\ TransformedTriangle_intersect.cxx\ TransformedTriangle_math.cxx\ -Interpolation3DSurf.cxx\ -Interpolation2D.cxx\ +MeshElement.cxx\ +MeshRegion.cxx\ +RegionNode.cxx\ +BoundingBox.cxx\ +TetraAffineTransform.cxx\ +Interpolation3D.cxx +#Interpolation3DSurf.cxx\ +#Interpolation2D.cxx\ # Executables targets BIN = @@ -63,7 +75,7 @@ BIN_SRC = BIN_SERVER_IDL = BIN_CLIENT_IDL = -TEST_PROGS = make_plane test_INTERPOL_2D +TEST_PROGS = #make_plane test_INTERPOL_2D LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome @@ -71,17 +83,17 @@ LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) CPPFLAGS+=$(BOOST_CPPFLAGS) -# enable testing -CXXFLAGS+= -DTESTING_INTERP_KERNEL -CPPFLAGS+= -DTESTING_INTERP_KERNEL + #LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil +#LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lutil #LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. LDFLAGSFORBIN+= -lm $(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 -lmedmem $(BOOST_LIBS) -lutil +#LDFLAGSFORBIN+= -lm $(HDF5_LIBS) -lmedmem $(BOOST_LIBS) -lutil ifeq ($(MED_WITH_KERNEL),yes) CPPFLAGS+= ${KERNEL_CXXFLAGS} diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index 1c4d22c54..620e36d86 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -99,8 +99,8 @@ namespace INTERP_UTILS const double* MeshElement::getCoordsOfNode(int node) const { const int nodeOffset = node - 1; - const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index]; - return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + nodeOffset]); + const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1; + return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + 3*nodeOffset]); } /** @@ -111,34 +111,41 @@ namespace INTERP_UTILS */ void MeshElement::triangulate(std::vector& 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]; - - double transformedPts[9]; - + 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]; - const double* pt = &(_mesh->getCoordinates(MED_FULL_INTERLACE)[connNodeIdx - 1]); + _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])); diff --git a/src/INTERP_KERNEL/MeshRegion.cxx b/src/INTERP_KERNEL/MeshRegion.cxx index 6a58263fd..72fd8388d 100644 --- a/src/INTERP_KERNEL/MeshRegion.cxx +++ b/src/INTERP_KERNEL/MeshRegion.cxx @@ -51,6 +51,14 @@ namespace INTERP_UTILS _box = new BoundingBox(pts, numNodes); + } else { + const int numNodes = element->getNumberNodes(); + + for(int i = 1 ; i <= numNodes ; ++i) + { + const double* pt = element->getCoordsOfNode(i); + _box->updateWithPoint(pt); + } } } diff --git a/src/INTERP_KERNEL/Test/Makefile.in b/src/INTERP_KERNEL/Test/Makefile.in index 46c821400..5e08e1d13 100644 --- a/src/INTERP_KERNEL/Test/Makefile.in +++ b/src/INTERP_KERNEL/Test/Makefile.in @@ -36,7 +36,8 @@ VPATH=.:@srcdir@:@top_srcdir@/idl # header files EXPORT_HEADERS = CppUnitTest.hxx \ TransformedTriangleTest.hxx \ - TransformedTriangleIntersectTest.hxx + TransformedTriangleIntersectTest.hxx \ + Interpolation3DTest.hxx # Libraries targets @@ -50,7 +51,8 @@ LIB_CLIENT_IDL = BIN = TestInterpKernel -BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx +BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx \ +Interpolation3DTest.cxx BIN_CLIENT_IDL = @@ -64,13 +66,10 @@ CPPFLAGS+=$(BOOST_CPPFLAGS) $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/sr CXXFLAGS+= -I/usr/include/cppunit CPPFLAGS+= -I/usr/include/cppunit -# enable testing -CXXFLAGS+= -DTESTING_INTERP_KERNEL -CPPFLAGS+= -DTESTING_INTERP_KERNEL - #LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil +#LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lmedmem -lutil -lmed #LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. @@ -89,7 +88,9 @@ LIBS = @LIBS@ @CPPUNIT_LIBS@ LDFLAGSFORBIN += $(LDFLAGS) -lm $(MED2_LIBS) $(HDF5_LIBS) \ -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome -lmed_V2_1 -lmedmem \ -lcppunit -linterpkernel - +#LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \ +# -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \ +# -lcppunit -linterpkernel UNIT_TEST_PROG = TestInterpKernel @CONCLUDE@ diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx index 777dd06c2..d8e5e6c2e 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.hxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx @@ -4,22 +4,22 @@ #include #include -#ifdef TESTING_INTERP_KERNEL -class TetraAffineTransformTest; -#endif +#include "VectorUtils.hxx" namespace INTERP_UTILS { class TetraAffineTransform { -#ifdef TESTING_INTERP_KERNEL - friend class ::TetraAffineTransformTest; -#endif + public: TetraAffineTransform(const double** pts) { + + std::cout << "Creating transform from tetraeder : " << std::endl; + std::cout << vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]) << ", " << std::endl; + #if 0 do { #endif @@ -61,6 +61,20 @@ namespace INTERP_UTILS // precalculate determinant (again after inversion of transform) calculateDeterminant(); + + 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; + 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)); + } + } + + std::cout << " ok" << std::endl; } void apply(double* destPt, const double* srcPt) const @@ -111,13 +125,14 @@ namespace INTERP_UTILS std::cout << "A = " << std::endl << "["; for(int i = 0; i < 3; ++i) { - cout << _linearTransform[3*i] << ", " << _linearTransform[3*i + 1] << ", " << _linearTransform[3*i + 1]; - if(i !=2 ) cout << endl; + cout << _linearTransform[3*i] << ", " << _linearTransform[3*i + 1] << ", " << _linearTransform[3*i + 2]; + if(i != 2 ) cout << endl; } cout << "]" << endl; cout << "b = " << "[" << _translation[0] << ", " << _translation[1] << ", " << _translation[2] << "]" << endl; } + private: void invertLinearTransform() diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 56e57a945..d2807adc7 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -155,6 +155,7 @@ namespace INTERP_UTILS if(isTriangleBelowTetraeder()) { + std::cout << std::endl << "Triangle is below tetraeder - V = 0.0" << std::endl << std::endl ; return 0.0; } @@ -171,6 +172,7 @@ namespace INTERP_UTILS if(sign == 0.0) { + std::cout << std::endl << "Triangle is perpendicular to z-plane - V = 0.0" << std::endl << std::endl; return 0.0; } @@ -369,10 +371,23 @@ namespace INTERP_UTILS { double* ptB = new double[3]; copyVector3(&_coords[5*corner], ptB); - _polygonA.push_back(ptB); + _polygonB.push_back(ptB); std::cout << "Inclusion XYZ-plane : " << vToStr(ptB) << " added to B" << std::endl; } + + // projection on XYZ - facet + if(testCornerAboveXYZFacet(corner)) + { + double* ptB = new double[3]; + copyVector3(&_coords[5*corner], ptB); + ptB[2] = 1 - ptB[0] - ptB[1]; + assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); + _polygonB.push_back(ptB); + std::cout << "Projection XYZ-plane : " << vToStr(ptB) << " added to B" << std::endl; + } + } + } /** diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 89e990907..79a828403 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -3,11 +3,11 @@ #include -#ifdef TESTING_INTERP_KERNEL + class TransformedTriangleTest; class TransformedTriangleIntersectTest; class TransformedTriangleCalcVolumeTest; -#endif + namespace INTERP_UTILS { @@ -26,11 +26,11 @@ namespace INTERP_UTILS public: -#ifdef TESTING_INTERP_KERNEL + friend class ::TransformedTriangleTest; friend class ::TransformedTriangleIntersectTest; friend class ::TransformedTriangleCalcVolumeTest; -#endif + /** * Enumerations representing the different geometric elements of the unit tetrahedron @@ -118,7 +118,8 @@ namespace INTERP_UTILS bool testCornerInTetrahedron(const TriCorner corner) const; bool testCornerOnXYZFacet(const TriCorner corner) const; - + + bool testCornerAboveXYZFacet(const TriCorner corner) const; //////////////////////////////////////////////////////////////////////////////////// /// Utility methods used in intersection tests /////////////// @@ -132,7 +133,7 @@ namespace INTERP_UTILS bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const; - bool testSegmentIntersectsH(const TriSegment seg); + bool testSegmentIntersectsHPlane(const TriSegment seg) const; bool testSurfaceAboveCorner(const TetraCorner corner) const; @@ -210,6 +211,14 @@ namespace INTERP_UTILS // coordinates of corners of tetrahedron static const double COORDS_TET_CORNER[12]; + + // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION + // for the calculation of the coordinates (x,y,z) of the intersection points + // for Segment-Facet and Segment-Edge intersections + static const int TransformedTriangle::DP_INDEX[12]; + + // correspondance edge - corners + static const TetraCorner CORNERS_FOR_EDGE[12]; }; diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index 2916da10b..0980b3b20 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -3,6 +3,7 @@ #include #include #include +#include namespace INTERP_UTILS { @@ -29,12 +30,22 @@ namespace INTERP_UTILS C_YH, C_XH, C_XY // Z }; - const long double TransformedTriangle::MACH_EPS = 1.0e-15; - const long double TransformedTriangle::MULT_PREC_F = 4.0*TransformedTriangle::MACH_EPS; + //const double TransformedTriangle::MACH_EPS = 1.0e-15; + const long double TransformedTriangle::MACH_EPS = std::numeric_limits::epsilon(); + const long double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS; const long double TransformedTriangle::THRESHOLD_F = 20.0; const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1; + // coordinates of corner ptTetCorner + /* const double TransformedTriangle::COORDS_TET_CORNER[12] = + { + 0.0, 0.0, 0.0, // O + 1.0, 0.0, 0.0, // X + 0.0, 1.0, 0.0, // Y + 0.0, 0.0, 1.0 // Z + }; + */ //////////////////////////////////////////////////////////////////////////////////// /// Double and triple product calculations /////////////// //////////////////////////////////////////////////////////////////////////////////// @@ -104,9 +115,9 @@ namespace INTERP_UTILS // coordinates of corner const double ptTetCorner[3] = { - corner == TT::X ? 1.0 : 0.0, - corner == TT::Y ? 1.0 : 0.0, - corner == TT::Z ? 1.0 : 0.0 + COORDS_TET_CORNER[3*corner ], + COORDS_TET_CORNER[3*corner + 1], + COORDS_TET_CORNER[3*corner + 2] }; // dist^2 = ( PQ x CP )^2 / |PQ|^2 where C is the corner point @@ -114,7 +125,9 @@ namespace INTERP_UTILS // difference vectors const double diffPQ[3] = { ptQ[0] - ptP[0], ptQ[1] - ptP[1], ptQ[2] - ptP[2] }; const double diffCornerP[3] = { ptP[0] - ptTetCorner[0], ptP[1] - ptTetCorner[1], ptP[2] - ptTetCorner[2] }; - + + //{ use functions in VectorUtils for this + // cross product of difference vectors const double cross[3] = { @@ -148,8 +161,7 @@ namespace INTERP_UTILS for(int i = 0 ; i < 3 ; ++i) { DoubleProduct dp = DOUBLE_PRODUCTS[3*min_corner + i]; - // std::cout << std::endl << "in code inconsistent dp :" << dp << std::endl; - + // std::cout << std::endl << "inconsistent dp :" << dp << std::endl; _doubleProducts[8*seg + dp] = 0.0; } @@ -172,18 +184,19 @@ namespace INTERP_UTILS const int off1 = DP_OFFSET_1[dp]; const int off2 = DP_OFFSET_2[dp]; - const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2]; - const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; - const double absTerm1 = (term1 > 0.0 ? term1 : -term1); - const double absTerm2 = (term2 > 0.0 ? term2 : -term2); + const long double term1 = static_cast(_coords[5*pt1 + off1]) * static_cast(_coords[5*pt2 + off2]); + const long double term2 = static_cast(_coords[5*pt1 + off2]) * static_cast(_coords[5*pt2 + off1]); - const double delta = MULT_PREC_F*(absTerm1 + absTerm2); - const double absDoubleProduct = _doubleProducts[8*seg + dp] > 0.0 ? _doubleProducts[8*seg + dp] : -_doubleProducts[8*seg + dp]; + const long double delta = MULT_PREC_F * ( std::abs(term1) + std::abs(term2) ); - if(absDoubleProduct < THRESHOLD_F*delta) + if( std::abs(static_cast(_doubleProducts[8*seg + dp])) < THRESHOLD_F*delta ) { - //std::cout << "Double product " << 8*seg+dp << " = " << absDoubleProduct << " is imprecise, reset to 0.0" << std::endl; + 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; + } _doubleProducts[8*seg + dp] = 0.0; + } } } @@ -202,8 +215,10 @@ namespace INTERP_UTILS if(_isTripleProductsCalculated) return; + 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; bool isGoodRowFound = false; // find edge / row to use @@ -227,7 +242,7 @@ namespace INTERP_UTILS // find normal to PQR - cross PQ and PR const double pq[3] = { - _coords[5*Q] - _coords[5*P], + _coords[5*Q] - _coords[5*P], _coords[5*Q + 1] - _coords[5*P + 1], _coords[5*Q + 2] - _coords[5*P + 2] }; @@ -288,12 +303,16 @@ namespace INTERP_UTILS { _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false); } + _validTP[corner] = true; } else { // this value will not be used // we set it to whatever + // std::cout << "Triple product not calculated for corner " << corner << std::endl; _tripleProducts[corner] = -3.14159265; + _validTP[corner] = false; + } } @@ -332,6 +351,7 @@ namespace INTERP_UTILS double TransformedTriangle::calcStableT(const TetraCorner corner) const { assert(_isTripleProductsCalculated); + // assert(_validTP[corner]); return _tripleProducts[corner]; } diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index e0dc18162..af5159671 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -3,12 +3,13 @@ #include #include -#include +#include +#include namespace INTERP_UTILS { - void copyVector3(const double* src, double* dest) + inline void copyVector3(const double* src, double* dest) { for(int i = 0 ; i < 3 ; ++i) { @@ -16,31 +17,31 @@ namespace INTERP_UTILS } } - const std::string vToStr(const double* pt) + inline const std::string vToStr(const double* pt) { std::stringstream ss(std::ios::out); ss << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]"; return ss.str(); } - void cross(const double* v1, const double* v2,double* res) + inline void cross(const double* v1, const double* v2,double* res) { res[0] = v1[1]*v2[2] - v1[2]*v2[1]; res[1] = v1[2]*v2[0] - v1[0]*v2[2]; res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } - double dot(const double* v1, const double* v2) + inline double dot(const double* v1, const double* v2) { return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]; } - double norm(const double* v) + inline double norm(const double* v) { return sqrt(dot(v,v)); } - double angleBetweenVectors(const double* v1, const double* v2, const double* n) + inline double angleBetweenVectors(const double* v1, const double* v2, const double* n) { const double denominator = dot(v1, v2); double v3[3]; @@ -51,6 +52,10 @@ namespace INTERP_UTILS } + inline bool epsilonEqual(const double x, const double y, const double errTol = 1.0e-12) + { + return std::abs(x - y) < errTol; + } };