Improve perf of P1P0,P0P1,P1P1

diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx
index 292800556345481bb54efcef0fe0f62b6f674f46..9b9aebbe2ef3246b7e4063a7635a39bcb340a208 100644 (file)
--- a/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx
+++ b/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx
@@ -78,9 +78,8 @@ namespace INTERP_KERNEL
   void PolyhedronIntersectorP0P1<MyMeshType,MyMatrix>::intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res)
   {
     SplitterTetra<MyMeshType>* subTetras[24];
-    int nbOfNodesT=Intersector3D<MyMeshType,MyMatrix>::_target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell));
     releaseArrays();
-    _split.splitTargetCell(targetCell,nbOfNodesT,_tetra);
+    _split.splitTargetCell2(targetCell,_tetra);
     for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
       {
         for(typename std::vector<SplitterTetra<MyMeshType>*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
@@ -92,7 +91,13 @@ namespace INTERP_KERNEL
                 double volume = tmp->intersectSourceCell(*iterCellS);
                 if(volume!=0.)
                   {
-                    typename MyMatrix::value_type& resRow=res[tmp->getId(0)];
+                    int targetNodeId(tmp->getId(0));
+                    if(targetNodeId<0)
+                      {
+                        std::ostringstream oss; oss << "PolyhedronIntersectorP0P1::intersectCells : On target cell #" <<  targetCell << " the splitting into tetra4 leads to the creation of an additional point that interacts with source cell Id #" << *iterCellS << " !";
+                        throw INTERP_KERNEL::Exception(oss.str().c_str());
+                      }
+                    typename MyMatrix::value_type& resRow=res[targetNodeId];
                     typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(*iterCellS));
                     if(iterRes==resRow.end())
                       resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(*iterCellS),volume));

diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx
index ccfa1d2563b19ff8dda4c2b8c8dadd0b459d2971..59baa8c68ac48fcdf7e738e6b3e23919689bc864 100644 (file)
--- a/src/INTERP_KERNEL/SplitterTetra.hxx
+++ b/src/INTERP_KERNEL/SplitterTetra.hxx
@@ -363,7 +363,7 @@ namespace INTERP_KERNEL
     
     SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
 
-    SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12]);
+    SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const int *conn = 0);
 
     ~SplitterTetra();
 
@@ -388,8 +388,8 @@ namespace INTERP_KERNEL
     void clearVolumesCache();
 
   private:
-    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 static void CheckIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL);
+    inline static void CheckIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL);
     inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
     inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
@@ -442,7 +442,7 @@ namespace INTERP_KERNEL
    * @param isOutside bool[8] which indicate the results of earlier checks. 
    */
   template<class MyMeshType>
-  inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
+  inline void SplitterTetra<MyMeshType>::CheckIsOutside(const double* pt, bool* isOutside, const double errTol)
   {
     isOutside[0] = isOutside[0] && (pt[0] < errTol);
     isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
@@ -455,7 +455,7 @@ namespace INTERP_KERNEL
   }
   
   template<class MyMeshType>
-  inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
+  inline void SplitterTetra<MyMeshType>::CheckIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol)
   {
     isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
     isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));

diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx
index 4bc794ae6d45e58c16be404d3f68759c9c2e597f..0cb18c057644164d1a7ec29d530f950958ddc270 100644 (file)
--- a/src/INTERP_KERNEL/SplitterTetra.txx
+++ b/src/INTERP_KERNEL/SplitterTetra.txx
@@ -99,9 +99,12 @@ namespace INTERP_KERNEL
    * \param [in] tetraCorners  array 4*3 doubles containing corners of input tetrahedron (P0X,P0Y,P0Y,P1X,P1Y,P1Z,P2X,P2Y,P2Z,P3X,P3Y,P3Z).
    */
   template<class MyMeshType>
-  SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12]): _t(0),_src_mesh(srcMesh)
+  SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const int *conn): _t(0),_src_mesh(srcMesh)
   {
-    _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3;
+    if(!conn)
+      { _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3; }
+    else
+      { _conn[0]=conn[0]; _conn[1]=conn[1]; _conn[2]=conn[2]; _conn[3]=conn[3]; }
     _coords[0]=tetraCorners[0]; _coords[1]=tetraCorners[1]; _coords[2]=tetraCorners[2]; _coords[3]=tetraCorners[3]; _coords[4]=tetraCorners[4]; _coords[5]=tetraCorners[5];
     _coords[6]=tetraCorners[6]; _coords[7]=tetraCorners[7]; _coords[8]=tetraCorners[8]; _coords[9]=tetraCorners[9]; _coords[10]=tetraCorners[10]; _coords[11]=tetraCorners[11];
     // create the affine transform
@@ -207,8 +210,7 @@ namespace INTERP_KERNEL
             //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
             calculateNode(globalNodeNum);
           }
-
-        checkIsOutside(_nodes[globalNodeNum], isOutside);       
+        CheckIsOutside(_nodes[globalNodeNum], isOutside);       
       }
 
     // halfspace filtering check
@@ -622,8 +624,8 @@ namespace INTERP_KERNEL
             calculateNode2(globalNodeNum, polyCoords[i]);
           }
 
-        checkIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside, precision);
-        checkIsOutside(_nodes[globalNodeNum], isOutside, precision);
+        CheckIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside, precision);
+        CheckIsOutside(_nodes[globalNodeNum], isOutside, precision);
       }
 
     // halfspace filtering check
@@ -848,7 +850,7 @@ namespace INTERP_KERNEL
     for(int i = 0;i<(int)nbOfNodes4Type;++i)
     {
       _t->apply(nodes[i], tetraCorners[i]);
-      checkIsOutside(nodes[i], isOutside);
+      CheckIsOutside(nodes[i], isOutside);
     }
 
     // halfspace filtering check
@@ -938,11 +940,13 @@ namespace INTERP_KERNEL
     SplitIntoTetras(_splitting_pol,gt,cellConn,refConn+_target_mesh.getConnectivityIndexPtr()[targetCell+1],coords,tetrasNodalConn,addCoords);
     std::size_t nbTetras(tetrasNodalConn.size()/4); tetra.resize(nbTetras);
     double tmp[12];
+    int tmp2[4];
     for(std::size_t i=0;i<nbTetras;i++)
       {
         for(int j=0;j<4;j++)
           {
             int cellId(tetrasNodalConn[4*i+j]);
+            tmp2[j]=cellId;
             if(cellId>=0)
               {
                 tmp[j*3+0]=coords[3*cellId+0];
@@ -956,7 +960,7 @@ namespace INTERP_KERNEL
                 tmp[j*3+2]=addCoords[3*(-cellId-1)+2];
               }
           }
-        tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp);
+        tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp,tmp2);
       }
   }