From 7cdf24ff04373e30ba5a159bee7796a289a1ffcc Mon Sep 17 00:00:00 2001 From: ageay Date: Thu, 1 Aug 2013 13:57:46 +0000 Subject: [PATCH] Start debugging 3D interpolation error on OCTA12 in target mesh --- src/INTERP_KERNEL/SplitterTetra.hxx | 27 ++---------- src/INTERP_KERNEL/SplitterTetra.txx | 9 +++- src/INTERP_KERNEL/TetraAffineTransform.cxx | 13 +++--- src/INTERP_KERNEL/TetraAffineTransform.hxx | 2 +- .../UnitTetraIntersectionBaryTest.cxx | 11 +++-- .../Test/MEDCouplingRemapperTest.cxx | 44 +++++++++++++++++++ .../Test/MEDCouplingRemapperTest.hxx | 3 ++ 7 files changed, 71 insertions(+), 38 deletions(-) diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index 2f3c5f0b8..1c8b1e4b9 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -152,7 +152,7 @@ namespace INTERP_KERNEL */ TriangleFaceKey(int node1, int node2, int node3) { - sort3Ints(_nodes, node1, node2, node3); + Sort3Ints(_nodes, node1, node2, node3); _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29; } @@ -201,7 +201,7 @@ namespace INTERP_KERNEL return _hashVal; } - inline void sort3Ints(int* sorted, int node1, int node2, int node3); + inline static void Sort3Ints(int* sorted, int node1, int node2, int node3); private: /// global numbers of the three nodes, sorted in ascending order @@ -219,7 +219,7 @@ namespace INTERP_KERNEL * @param x2 second integer * @param x3 third integer */ - inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3) + inline void TriangleFaceKey::Sort3Ints(int* sorted, int x1, int x2, int x3) { if(x1 < x2) { @@ -316,8 +316,6 @@ namespace INTERP_KERNEL void clearVolumesCache(); private: - // member functions - inline void createAffineTransform(const double** corners); inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const; inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const; inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum); @@ -347,11 +345,7 @@ namespace INTERP_KERNEL HashMap< int , double* > _nodes; /// HashMap relating triangular faces to calculated volume contributions, used for caching - HashMap< TriangleFaceKey, double -// #ifdef WIN32 -// , hash_compare -// #endif - > _volumes; + HashMap< TriangleFaceKey, double > _volumes; /// reference to the source mesh const MyMeshType& _src_mesh; @@ -366,19 +360,6 @@ namespace INTERP_KERNEL static const double SPARSE_TRUNCATION_LIMIT; }; - /** - * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors - * - * @param corners double*[4] array containing pointers to four double[3] arrays with the - * coordinates of the corners of the tetrahedron - */ - template - inline void SplitterTetra::createAffineTransform(const double** corners) - { - // create AffineTransform from tetrahedron - _t = new TetraAffineTransform( corners ); - } - /** * Function used to filter out elements by checking if they belong to one of the halfspaces * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 6364c79d9..9312b5fa6 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -59,6 +59,13 @@ namespace INTERP_KERNEL } /** + * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a + * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish + * if a partial computation of a triangle has already been performed (to increase performance). + * + * The \a srcMesh can contain polyhedron cells. + * + * * Constructor creating object from the four corners of the tetrahedron. * * @param srcMesh mesh containing the source elements @@ -75,7 +82,7 @@ namespace INTERP_KERNEL _coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2]; _coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2]; // create the affine transform - createAffineTransform(tetraCorners); + _t=new TetraAffineTransform(_coords); } /** diff --git a/src/INTERP_KERNEL/TetraAffineTransform.cxx b/src/INTERP_KERNEL/TetraAffineTransform.cxx index 87690657c..a71ec1342 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.cxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.cxx @@ -38,13 +38,12 @@ namespace INTERP_KERNEL * with corners specified in pts. If the tetrahedron is degenerate or almost degenerate, * construction succeeds, but the determinant of the transform is set to 0. * - * @param pts a 4x3 matrix containing 4 points (pts[0], ..., pts[3]) of 3 coordinates each + * @param pts a 4x3 matrix containing 4 points (P0X,P0Y,P0Z,P1X,P1Y,P1Z,...) of 3 coordinates each */ - TetraAffineTransform::TetraAffineTransform(const double** pts) + TetraAffineTransform::TetraAffineTransform(const double *pts) { LOG(2,"Creating transform from tetraeder : "); - LOG(2, vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3])); // three last points -> linear transform for(int i = 0; i < 3 ; ++i) @@ -52,13 +51,13 @@ namespace INTERP_KERNEL for(int j = 0 ; j < 3 ; ++j) { // NB we insert columns, not rows - _linear_transform[3*j + i] = (pts[i+1])[j] - (pts[0])[j]; + _linear_transform[3*j + i] = pts[(i+1)*3+j] - pts[j]; } } // remember _linear_transform for the reverse transformation memcpy( _back_linear_transform, _linear_transform, 9*sizeof(double)); - memcpy( _back_translation, pts[0], 3*sizeof(double)); + memcpy( _back_translation, pts, 3*sizeof(double)); calculateDeterminant(); @@ -81,7 +80,7 @@ namespace INTERP_KERNEL // and O is the position vector of the point that is mapped onto the origin for(int i = 0 ; i < 3 ; ++i) { - _translation[i] = -(_linear_transform[3*i]*(pts[0])[0] + _linear_transform[3*i+1]*(pts[0])[1] + _linear_transform[3*i+2]*(pts[0])[2]) ; + _translation[i] = -(_linear_transform[3*i]*pts[0] + _linear_transform[3*i+1]*pts[1] + _linear_transform[3*i+2]*pts[2]) ; } // precalculate determinant (again after inversion of transform) @@ -95,7 +94,7 @@ namespace INTERP_KERNEL for(int i = 0; i < 4 ; ++i) { double v[3]; - apply(v, pts[i]); + apply(v, pts+3*i); LOG(4, vToStr(v)) for(int j = 0; j < 3; ++j) { diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx index de1e8ee6c..258085df2 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.hxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx @@ -35,7 +35,7 @@ namespace INTERP_KERNEL { public: - TetraAffineTransform(const double** pts); + TetraAffineTransform(const double *pts); void apply(double* destPt, const double* srcPt) const; diff --git a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx index 582174532..6cda7ee9c 100644 --- a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx +++ b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx @@ -301,15 +301,14 @@ namespace INTERP_TEST void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply() { - double nodes[4][3] = { {-4.0, 9.0, 3.0 }, - {11.0, 0.0, 2.0 }, - { 0.0, 0.0, 0.0 }, - { 2.0, 1.0,10.0 }}; + double nodes[12] = { -4.0, 9.0, 3.0, + 11.0, 0.0, 2.0, + 0.0, 0.0, 0.0, + 2.0, 1.0,10.0 }; // double pSrc[3] = { -4.0, 9.0, 3.0 }; double pSrc[3] = { 40., -20., 100. }; double pDest[] = {1,1,1}; - const double* n[4] = { &nodes[0][0], &nodes[1][0], &nodes[2][0], &nodes[3][0] }; - TetraAffineTransform a(&n[0]); + TetraAffineTransform a(nodes); a.apply( pDest, pSrc ); a.reverseApply( pDest, pDest ); CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[0], pDest[0], 1e-12); diff --git a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx index 22bbd6b8b..f62ee6832 100644 --- a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx @@ -1272,3 +1272,47 @@ void MEDCouplingRemapperTest::testPartialTransfer1() sourceMesh->decrRef(); targetMesh->decrRef(); } + +void MEDCouplingRemapperTest::testBugNonRegression1() +{ + // source + DataArrayDouble *coordsSrc(DataArrayDouble::New()); + const double coordsSrcData[18]={-6.25,3.6084391824351605,264.85199999999998,-6.25,3.6084391824351605,289.05200000000002,-6.2499999999999991,-3.6084391824351618,264.85199999999998,-6.2499999999999991,-3.6084391824351618,289.05200000000002,-1.7763568394002505e-15,4.4408920985006262e-15,264.85199999999998,-1.7763568394002505e-15,4.4408920985006262e-15,289.05200000000002}; + coordsSrc->useArray(coordsSrcData,false,CPP_DEALLOC,6,3); + DataArrayInt *connSrc(DataArrayInt::New()),*connISrc(DataArrayInt::New()); + const int connSrcData[7]={16,2,0,4,3,1,5}; + connSrc->useArray(connSrcData,false,CPP_DEALLOC,7,1); + const int connISrcData[2]={0,7}; + connISrc->useArray(connISrcData,false,CPP_DEALLOC,2,1); + MEDCouplingUMesh *srcMesh(MEDCouplingUMesh::New("source",3)); + srcMesh->setCoords(coordsSrc); + srcMesh->setConnectivity(connSrc,connISrc,true); + coordsSrc->decrRef(); connSrc->decrRef(); connISrc->decrRef(); + // target + DataArrayDouble *coordsTrg(DataArrayDouble::New()); +const double coordsTrgData[36]={-2,1.1547005383792521,264.85199999999998,-2,0.57735026918962618,264.85199999999998,-2.5,0.2886751345948132,264.85199999999998,-2.5,1.443375672974065,264.85199999999998,-3.0000000000000004,1.1547005383792526,264.85199999999998,-3.0000000000000004,0.57735026918962662,264.85199999999998,-2,1.1547005383792521,289.05200000000002,-2,0.57735026918962618,289.05200000000002,-2.5,0.2886751345948132,289.05200000000002,-2.5,1.443375672974065,289.05200000000002,-3.0000000000000004,1.1547005383792526,289.05200000000002,-3.0000000000000004,0.57735026918962662,289.05200000000002}; + coordsTrg->useArray(coordsTrgData,false,CPP_DEALLOC,12,3); + DataArrayInt *connTrg=DataArrayInt::New(); + const int connTrgData[44]={31,0,1,2,5,4,3,-1,7,6,9,10,11,8,-1,3,9,6,0,-1,4,10,9,3,-1,5,11,10,4,-1,2,8,11,5,-1,1,7,8,2,-1,0,6,7,1}; + connTrg->useArray(connTrgData,false,CPP_DEALLOC,44,1); + DataArrayInt *connITrg=DataArrayInt::New(); + const int connITrgData[2]={0,44}; + connITrg->useArray(connITrgData,false,CPP_DEALLOC,2,1); + MEDCouplingUMesh *trgMesh=MEDCouplingUMesh::New("target",3); + trgMesh->setCoords(coordsTrg); + trgMesh->setConnectivity(connTrg,connITrg,true); + coordsTrg->decrRef(); connTrg->decrRef(); connITrg->decrRef(); + // Go ! + const double valExpected(20.957814771583468); + MEDCouplingRemapper remapper; + remapper.setPrecision(1e-12); + remapper.setIntersectionType(INTERP_KERNEL::Triangulation); + CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(srcMesh,trgMesh,"P0P0")); + std::vector > matrx(remapper.getCrudeMatrix()); + CPPUNIT_ASSERT_EQUAL(1,(int)matrx.size()); + CPPUNIT_ASSERT_EQUAL(1,(int)matrx[0].size()); + std::cerr << std::endl << matrx[0][0] << std::endl; + // + srcMesh->decrRef(); trgMesh->decrRef(); +} + diff --git a/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx b/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx index c7d408967..d3a5df100 100644 --- a/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx @@ -43,6 +43,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( testExtruded2 ); CPPUNIT_TEST( testPrepareEx1 ); CPPUNIT_TEST( testPartialTransfer1 ); + CPPUNIT_TEST( testBugNonRegression1 ); CPPUNIT_TEST_SUITE_END(); public: void test2DInterpP0P0_1(); @@ -55,6 +56,8 @@ namespace ParaMEDMEM void testExtruded2(); void testPrepareEx1(); void testPartialTransfer1(); + // + void testBugNonRegression1(); private: static MEDCouplingUMesh *build1DTargetMesh_2(); static MEDCouplingUMesh *build2DTargetMesh_3(); -- 2.39.2