]> SALOME platform Git repositories - modules/med.git/commitdiff
Salome HOME
Start debugging 3D interpolation error on OCTA12 in target mesh B4PolyhIntersect
authorageay <ageay>
Thu, 1 Aug 2013 13:57:46 +0000 (13:57 +0000)
committerageay <ageay>
Thu, 1 Aug 2013 13:57:46 +0000 (13:57 +0000)
src/INTERP_KERNEL/SplitterTetra.hxx
src/INTERP_KERNEL/SplitterTetra.txx
src/INTERP_KERNEL/TetraAffineTransform.cxx
src/INTERP_KERNEL/TetraAffineTransform.hxx
src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx
src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx
src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx

index 2f3c5f0b8ac0d4ce475b7d91869638655cae1905..1c8b1e4b9d583507fdc45f618db0ee70a3c84fe7 100644 (file)
@@ -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<TriangleFaceKey,TriangleFaceKeyComparator> 
-// #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<class MyMeshType>
-  inline void SplitterTetra<MyMeshType>::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
index 6364c79d9e96ae2f5d4c7bd54a76514bd7348a52..9312b5fa6eb8d31dfc81259d5e712d3eab7825d2 100644 (file)
@@ -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);
   }
   
   /**
index 87690657c8026de81c94e815d04ac71d2eb054d4..a71ec134233c63107e60e6ba185ec267f665822c 100644 (file)
@@ -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)
             {
index de1e8ee6c8312a813ce17124b7f1359413f868eb..258085df28d35ada4dbc663ac03ec9fc365c4150 100644 (file)
@@ -35,7 +35,7 @@ namespace INTERP_KERNEL
   {
 
   public:
-    TetraAffineTransform(const double** pts);
+    TetraAffineTransform(const double *pts);
 
     void apply(double* destPt, const double* srcPt) const;
 
index 582174532a034d558dff3a573172b1c9f1970d88..6cda7ee9cd290dfbe659c03019381b11a952325e 100644 (file)
@@ -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);
index 22bbd6b8bb702d6a6630abf1d6ac97f4fad7f4bb..f62ee6832fdf45017802ed2ee6907a72efaa709e 100644 (file)
@@ -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<std::map<int,double> > 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();
+}
+
index c7d40896720d1c5107924669342dcc697b74223e..d3a5df1000fd659ff8f2689df1b7bc8e74963a27 100644 (file)
@@ -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();