Dealing polyhedra for 3D P0P0 non conservative interpolation.

diff --git a/src/INTERP_KERNEL/Intersector3D.txx b/src/INTERP_KERNEL/Intersector3D.txx
index 501bb861d59675a4948e880b70c460358aabb2f1..ace541a5b09c2723eb7918f4e0db5fd2b96a3aca 100644 (file)
--- a/src/INTERP_KERNEL/Intersector3D.txx
+++ b/src/INTERP_KERNEL/Intersector3D.txx
@@ -36,14 +36,14 @@ namespace INTERP_KERNEL
   template<class MyMeshType, class MyMatrix>
   void Intersector3D<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const
   {
-    const ConnType *myConectT=_target_mesh.getConnectivityPtr();
-    const ConnType *myConIndexT=_target_mesh.getConnectivityIndexPtr();
-    const double *myCoordsT=_target_mesh.getCoordinatesPtr();
-    int nbNodesT=myConIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-myConIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+    int nbNodesT=_target_mesh.getNumberOfNodesOfElement(icellT);
     coordsT.resize(SPACEDIM*nbNodesT);
+    std::vector<double>::iterator iter=coordsT.begin();
     for (ConnType iT=0; iT<nbNodesT; iT++)
-      for(int idim=0; idim<SPACEDIM; idim++)
-        coordsT[SPACEDIM*iT+idim]=myCoordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(myConectT[OTT<ConnType,numPol>::conn2C(myConIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
+      {
+        const double *coordsCur=getCoordsOfNode(iT,icellT,_target_mesh);
+        iter=std::copy(coordsCur,coordsCur+SPACEDIM,iter);
+      }
   }
 
   /*!
@@ -52,14 +52,14 @@ namespace INTERP_KERNEL
   template<class MyMeshType, class MyMatrix>
   void Intersector3D<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS) const
   {
-    const ConnType *myConectS=_src_mesh.getConnectivityPtr();
-    const ConnType *myConIndexS=_src_mesh.getConnectivityIndexPtr();
-    const double *myCoordsS=_src_mesh.getCoordinatesPtr();
-    int nbNodesS=myConIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-myConIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+    int nbNodesS=_src_mesh.getNumberOfNodesOfElement(icellS);
     coordsS.resize(SPACEDIM*nbNodesS);
+    std::vector<double>::iterator iter=coordsS.begin();
     for (ConnType iS=0; iS<nbNodesS; iS++)
-      for(int idim=0; idim<SPACEDIM; idim++)
-        coordsS[SPACEDIM*iS+idim]=myCoordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(myConectS[OTT<ConnType,numPol>::conn2C(myConIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
+      {
+        const double *coordsCur=getCoordsOfNode(iS,icellS,_src_mesh);
+        iter=std::copy(coordsCur,coordsCur+SPACEDIM,iter);
+      }
   }
 
   /*!

diff --git a/src/INTERP_KERNEL/MeshUtils.hxx b/src/INTERP_KERNEL/MeshUtils.hxx
index a39634ce27456c03daf0431aaa5e0fa2f24f0a2b..9025a75218748806b4b8580815fb00a6b6419d60 100644 (file)
--- a/src/INTERP_KERNEL/MeshUtils.hxx
+++ b/src/INTERP_KERNEL/MeshUtils.hxx
@@ -29,7 +29,7 @@ namespace INTERP_KERNEL
    * @param      node       the node for which the global number is sought (ALWAYS in C mode)
    * @param      element    an element of the mesh (in numPol policy)
    * @param      mesh       a mesh
-   * @return    the node's global number so that (its coordinates in the coordinates array are at [SPACEDIM*globalNumber, SPACEDIM*globalNumber + 2]
+   * @return    the node's global number so that (its coordinates in the coordinates array are at [SPACEDIM*globalNumber, SPACEDIM*globalNumber + SPACEDIM]
    */
   template<class MyMeshType>
   inline typename MyMeshType::MyConnType getGlobalNumberOfNode(typename MyMeshType::MyConnType node, typename MyMeshType::MyConnType element, const MyMeshType& mesh)
@@ -37,7 +37,20 @@ namespace INTERP_KERNEL
     typedef typename MyMeshType::MyConnType ConnType;
     const NumberingPolicy numPol=MyMeshType::My_numPol;
     const ConnType elemIdx=OTT<ConnType,numPol>::conn2C(mesh.getConnectivityIndexPtr()[OTT<ConnType,numPol>::ind2C(element)]);
-    return OTT<ConnType,numPol>::coo2C(mesh.getConnectivityPtr()[elemIdx + node]);
+    if(mesh.getTypeOfElement(element)!=INTERP_KERNEL::NORM_POLYHED)
+      return OTT<ConnType,numPol>::coo2C(mesh.getConnectivityPtr()[elemIdx + node]);
+    else
+      {
+        const ConnType *startNodalConnOfElem=mesh.getConnectivityPtr()+elemIdx;
+        ConnType ptr=0,ret=0;
+        while(startNodalConnOfElem[ret]==-1 || ptr!=node)
+          {
+            ret++;
+            if(startNodalConnOfElem[ret]!=-1)
+              ptr++;
+          }
+        return OTT<ConnType,numPol>::coo2C(startNodalConnOfElem[ret]);
+      }
   }
 
   /**
@@ -53,7 +66,8 @@ namespace INTERP_KERNEL
   {
     typedef typename MyMeshType::MyConnType ConnType;
     const ConnType connIdx = getGlobalNumberOfNode(node, element, mesh);
-    return mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*connIdx;
+    const double *ret=mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*connIdx;
+    return ret;
   }
 
   /**

diff --git a/src/INTERP_KERNEL/PointLocator3DIntersectorP0P0.txx b/src/INTERP_KERNEL/PointLocator3DIntersectorP0P0.txx
index 5de6bae9c41f34f122933eb3a7df2c938f9e0735..cc9600a7dccab43664fcafc43da7fdfa66dc2ad1 100644 (file)
--- a/src/INTERP_KERNEL/PointLocator3DIntersectorP0P0.txx
+++ b/src/INTERP_KERNEL/PointLocator3DIntersectorP0P0.txx
@@ -66,7 +66,7 @@ namespace INTERP_KERNEL
         const CellModel& cmTypeS=CellModel::getCellModel(tS);
         std::vector<ConnType> connOfCurCellS;
         Intersector3DP0P0<MyMeshType,MyMatrix>::getConnOfSourceCell(OTT<ConnType,numPol>::indFC(*iterCellS),connOfCurCellS);
-        if(PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(bary,&connOfCurCellS[0],coordsS,cmTypeS,_precision))
+        if(PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(bary,&connOfCurCellS[0],connOfCurCellS.size(),coordsS,cmTypeS,_precision))
           {
             resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(*iterCellS),1));
           }

diff --git a/src/INTERP_KERNEL/PointLocator3DIntersectorP0P1.txx b/src/INTERP_KERNEL/PointLocator3DIntersectorP0P1.txx
index c59f8a962100cacea5d226cbec1b3a3f9979aafc..33f529edc86238238ef8e19b3c4e6303b63e1704 100644 (file)
--- a/src/INTERP_KERNEL/PointLocator3DIntersectorP0P1.txx
+++ b/src/INTERP_KERNEL/PointLocator3DIntersectorP0P1.txx
@@ -64,7 +64,7 @@ namespace INTERP_KERNEL
         Intersector3DP0P1<MyMeshType,MyMatrix>::getConnOfSourceCell(OTT<ConnType,numPol>::indFC(*iterCellS),connOfCurCellS);
         for(int nodeIdT=0;nodeIdT<nbNodesT;nodeIdT++)
           {
-            if(PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(&coordsTarget[nodeIdT*SPACEDIM],&connOfCurCellS[0],coordsS,cmTypeS,_precision))
+            if(PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(&coordsTarget[nodeIdT*SPACEDIM],&connOfCurCellS[0],connOfCurCellS.size(),coordsS,cmTypeS,_precision))
               {
                 ConnType curNodeTInCmode=OTT<ConnType,numPol>::coo2C(startOfCellNodeConnT[nodeIdT]);
                 typename MyMatrix::value_type& resRow=res[curNodeTInCmode];

diff --git a/src/INTERP_KERNEL/PointLocator3DIntersectorP1P0.txx b/src/INTERP_KERNEL/PointLocator3DIntersectorP1P0.txx
index e666f3a2923147ebb4b9880a5ed68cb6c82d3dc5..bc16de9a27da2aec0cfb609c9417993adbba42eb 100644 (file)
--- a/src/INTERP_KERNEL/PointLocator3DIntersectorP1P0.txx
+++ b/src/INTERP_KERNEL/PointLocator3DIntersectorP1P0.txx
@@ -65,7 +65,7 @@ namespace INTERP_KERNEL
         const CellModel& cmTypeS=CellModel::getCellModel(tS);
         std::vector<ConnType> connOfCurCellS;
         Intersector3DP1P0<MyMeshType,MyMatrix>::getConnOfSourceCell(OTT<ConnType,numPol>::indFC(*iterCellS),connOfCurCellS);
-        if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(baryT,&connOfCurCellS[0],coordsS,cmTypeS,_precision) )
+        if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(baryT,&connOfCurCellS[0],connOfCurCellS.size(),coordsS,cmTypeS,_precision) )
           {
             double resLoc[4];
             std::vector<double> srcCell;

diff --git a/src/INTERP_KERNEL/PointLocator3DIntersectorP1P1.txx b/src/INTERP_KERNEL/PointLocator3DIntersectorP1P1.txx
index 8b29fc250d0f6fac0ad61a492719a6f7e4e19089..9670cffaa77d93d7bf9fabc248ad3c07e76fd9eb 100644 (file)
--- a/src/INTERP_KERNEL/PointLocator3DIntersectorP1P1.txx
+++ b/src/INTERP_KERNEL/PointLocator3DIntersectorP1P1.txx
@@ -67,7 +67,7 @@ namespace INTERP_KERNEL
             typename MyMatrix::value_type& resRow=res[OTT<ConnType,numPol>::ind2C(startOfCellNodeConnT[nodeIdT])];
             std::vector<ConnType> connOfCurCellS;
             Intersector3DP1P1<MyMeshType,MyMatrix>::getConnOfSourceCell(OTT<ConnType,numPol>::indFC(*iterCellS),connOfCurCellS);
-            if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(&CoordsT[nodeIdT*SPACEDIM],&connOfCurCellS[0],coordsS,cmTypeS,_precision) )
+            if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(&CoordsT[nodeIdT*SPACEDIM],&connOfCurCellS[0],connOfCurCellS.size(),coordsS,cmTypeS,_precision) )
               {
                 double resLoc[4];
                 std::vector<double> localCoordsS;

diff --git a/src/INTERP_KERNEL/PointLocatorAlgos.txx b/src/INTERP_KERNEL/PointLocatorAlgos.txx
index 1db1d3c7fa0ff5bcb45c00d93c0f02c105de641a..5379ce1caa9d2d884da199817be05bb378e373f4 100644 (file)
--- a/src/INTERP_KERNEL/PointLocatorAlgos.txx
+++ b/src/INTERP_KERNEL/PointLocatorAlgos.txx
@@ -131,21 +131,22 @@ namespace INTERP_KERNEL
       return ret;
     }
 
-    static bool isElementContainsPointAlg3D(const double *ptToTest, const int *conn_elem, const double *coords, const CellModel& cmType, double eps)
+    static bool isElementContainsPointAlg3D(const double *ptToTest, const typename MyMeshType::MyConnType *conn_elem, int conn_elem_sz, const double *coords, const CellModel& cmType, double eps)
     {
       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
       typedef typename MyMeshType::MyConnType ConnType;
       const NumberingPolicy numPol=MyMeshType::My_numPol;
       
-      int nbfaces = cmType.getNumberOfSons();
-      int* sign = new int[nbfaces];
+      int nbfaces = cmType.getNumberOfSons2(conn_elem,conn_elem_sz);
+      int *sign = new int[nbfaces];
+      int *connOfSon = new int[conn_elem_sz];
       for (int iface=0; iface<nbfaces; iface++)
         {
-          const unsigned* connface=cmType.getNodesConstituentTheSon(iface);
-          const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[0]]));
-          const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[1]]));
-          const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[2]]));
-                                                        
+          NormalizedCellType typeOfSon;
+          cmType.fillSonCellNodalConnectivity2(iface,conn_elem,conn_elem_sz,connOfSon,typeOfSon);
+          const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[0]));
+          const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[1]));
+          const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[2]));                                                       
           double Vol=triple_product(AA,BB,CC,ptToTest);
           if (Vol<-eps)
             sign[iface]=-1;
@@ -156,10 +157,11 @@ namespace INTERP_KERNEL
         }
       bool ret=decideFromSign(sign, nbfaces);
       delete [] sign;
+      delete [] connOfSon;
       return ret;
     }
 
-    static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords, const int *conn_elem)
+    static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords, const typename MyMeshType::MyConnType *conn_elem, int conn_elem_sz)
     {
       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
       typedef typename MyMeshType::MyConnType ConnType;
@@ -183,7 +185,7 @@ namespace INTERP_KERNEL
                         
       if (SPACEDIM==3)
         {
-          return isElementContainsPointAlg3D(ptToTest,conn_elem,coords,cmType,1e-12);
+          return isElementContainsPointAlg3D(ptToTest,conn_elem,conn_elem_sz,coords,cmType,1e-12);
         }
     }
         
@@ -198,8 +200,9 @@ namespace INTERP_KERNEL
       const ConnType* conn=_mesh.getConnectivityPtr();
       const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
       const ConnType* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
+      int conn_elem_sz=conn_index[i+1]-conn_index[i];
       NormalizedCellType type=_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
-      return isElementContainsPoint(x,type,coords,conn_elem);
+      return isElementContainsPoint(x,type,coords,conn_elem,conn_elem_sz);
     }
                 
     static bool decideFromSign(const int* sign, int nbelem)

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 84178f580c4ad2163a3aad17193e18263f70d6f0..01790c40854f2615a1f3c67e3df3dd97ea810235 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -572,8 +572,12 @@ INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(int cellId) co
 
 int MEDCouplingUMesh::getNumberOfNodesInCell(int cellId) const
 {
-  int *ptI=_nodal_connec_index->getPointer();
-  return ptI[cellId+1]-ptI[cellId]-1;
+  const int *ptI=_nodal_connec_index->getConstPointer();
+  const int *pt=_nodal_connec->getConstPointer();
+  if(pt[ptI[cellId]]!=INTERP_KERNEL::NORM_POLYHED)
+    return ptI[cellId+1]-ptI[cellId]-1;
+  else
+    return std::count_if(pt+ptI[cellId]+1,pt+ptI[cellId+1],std::bind2nd(std::not_equal_to<int>(),-1));
 }
 
 void MEDCouplingUMesh::setConnectivity(DataArrayInt *conn, DataArrayInt *connIndex, bool isComputingTypes)

diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
index 4d43338063cdda7d9f48a926a1f80acf58cac0a4..6aaea0edf2c7ac854aee85c5720e673108669e8b 100644 (file)
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
@@ -1652,6 +1652,142 @@ void MEDCouplingBasicsTest::test3DInterpP0P0PL_1()
   targetMesh->decrRef();
 }
 
+void MEDCouplingBasicsTest::test3DInterpP0P0PL_2()
+{
+  MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1();
+  MEDCouplingUMesh *targetMesh=build3DTargetMesh_1();
+  std::vector<int> cellsIds(targetMesh->getNumberOfCells());
+  for(int i=0;i<targetMesh->getNumberOfCells();i++)
+    cellsIds[i]=i;
+  targetMesh->convertToPolyTypes(cellsIds);
+  //
+  MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
+  MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
+  INTERP_KERNEL::Interpolation3D myInterpolator;
+  vector<map<int,double> > res;
+  myInterpolator.setPrecision(1e-12);
+  myInterpolator.setIntersectionType(INTERP_KERNEL::PointLocator);
+  myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+  CPPUNIT_ASSERT_EQUAL(8,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[1][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[1][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[2][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[2][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[3][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[3][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[4][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[4][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[5][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[5][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[6][9],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][3],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][4],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][5],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][11],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(21.,sumAll(res),1e-12);
+  //clean up
+  sourceMesh->decrRef();
+  targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test3DInterpP0P0PL_3()
+{
+  MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1();
+  MEDCouplingUMesh *targetMesh=build3DTargetMesh_1();
+  std::vector<int> cellsIds(sourceMesh->getNumberOfCells());
+  for(int i=0;i<sourceMesh->getNumberOfCells();i++)
+    cellsIds[i]=i;
+  sourceMesh->convertToPolyTypes(cellsIds);
+  //
+  MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
+  MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
+  INTERP_KERNEL::Interpolation3D myInterpolator;
+  vector<map<int,double> > res;
+  myInterpolator.setPrecision(1e-12);
+  myInterpolator.setIntersectionType(INTERP_KERNEL::PointLocator);
+  myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+  CPPUNIT_ASSERT_EQUAL(8,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[1][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[1][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[2][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[2][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[3][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[3][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[4][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[4][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[5][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[5][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[6][9],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][3],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][4],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][5],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][11],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(21.,sumAll(res),1e-12);
+  //clean up
+  sourceMesh->decrRef();
+  targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test3DInterpP0P0PL_4()
+{
+  MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1();
+  MEDCouplingUMesh *targetMesh=build3DTargetMesh_1();
+  std::vector<int> cellsIds(sourceMesh->getNumberOfCells());
+  for(int i=0;i<sourceMesh->getNumberOfCells();i++)
+    cellsIds[i]=i;
+  sourceMesh->convertToPolyTypes(cellsIds);
+  cellsIds.resize(targetMesh->getNumberOfCells());
+  for(int j=0;j<targetMesh->getNumberOfCells();j++)
+    cellsIds[j]=j;
+  targetMesh->convertToPolyTypes(cellsIds);
+  //
+  MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
+  MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
+  INTERP_KERNEL::Interpolation3D myInterpolator;
+  vector<map<int,double> > res;
+  myInterpolator.setPrecision(1e-12);
+  myInterpolator.setIntersectionType(INTERP_KERNEL::PointLocator);
+  myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+  CPPUNIT_ASSERT_EQUAL(8,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[0][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[1][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[1][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[2][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[2][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[3][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[3][8],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[4][6],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[4][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[5][7],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[5][10],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[6][9],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][3],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][4],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][5],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,res[7][11],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(21.,sumAll(res),1e-12);
+  //clean up
+  sourceMesh->decrRef();
+  targetMesh->decrRef();
+}
+
 void MEDCouplingBasicsTest::test3DInterpP0P1_1()
 {
   MEDCouplingUMesh *sourceMesh=build3DTargetMesh_1();

diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
index 7ce21d7da72b133e5687d7cd1b347ab2fbe8a804..a63117dc339e6f35f0f6c6b049902748473217ba 100644 (file)
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
@@ -70,6 +70,9 @@ namespace ParaMEDMEM
     CPPUNIT_TEST( testInterpolationCC );
     CPPUNIT_TEST( test3DInterpP0P0_1 );
     CPPUNIT_TEST( test3DInterpP0P0PL_1 );
+    CPPUNIT_TEST( test3DInterpP0P0PL_2 );
+    CPPUNIT_TEST( test3DInterpP0P0PL_3 );
+    CPPUNIT_TEST( test3DInterpP0P0PL_4 );
     CPPUNIT_TEST( test3DInterpP0P1_1 );
     CPPUNIT_TEST( test3DInterpP0P1PL_1 );
     CPPUNIT_TEST( test3DInterpP1P0_1 );
@@ -125,6 +128,9 @@ namespace ParaMEDMEM
     void test3DSurfInterpP0P0_3();
     void test3DInterpP0P0_1();
     void test3DInterpP0P0PL_1();
+    void test3DInterpP0P0PL_2();
+    void test3DInterpP0P0PL_3();
+    void test3DInterpP0P0PL_4();
     void test3DInterpP0P1_1();
     void test3DInterpP0P1PL_1();
     void test3DInterpP1P0_1();