From 4f3fb5a026df7d4be4fa63578f6b9426ce6bc0e6 Mon Sep 17 00:00:00 2001
From: Anthony Geay <anthony.geay@edf.fr>
Date: Thu, 6 Oct 2016 08:55:15 +0200
Subject: [PATCH] Correct bug on P0P0 interpolation of Extruded<->Extruded
 case. Now if Z lev are different each other it works...

---
 src/MEDCoupling/MEDCouplingRemapper.cxx       | 10 +++--
 .../MEDCouplingRemapperTest.py                | 45 +++++++++++++++++++
 2 files changed, 52 insertions(+), 3 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx
index a8c8dae00..276802f81 100644
--- a/src/MEDCoupling/MEDCouplingRemapper.cxx
+++ b/src/MEDCoupling/MEDCouplingRemapper.cxx
@@ -620,9 +620,11 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
   const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
   if(methC!="P0P0")
     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
-  MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
-  MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
-  INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
+  MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
+  MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
+  MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
+  MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
+  INTERP_KERNEL::Interpolation2D interpolation2D(*this);
   std::vector<std::map<int,double> > matrix2D;
   int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
   MEDCouplingUMesh *s1D,*t1D;
@@ -632,6 +634,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
   std::vector<std::map<int,double> > matrix1D;
   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
+  if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
+    interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
   int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
   s1D->decrRef();
   t1D->decrRef();
diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
index 061af89b7..75d137660 100644
--- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
@@ -1019,6 +1019,51 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertEqual(rem.getCrudeMatrix(),[{0:1.},{1:1.}])
         pass
     
+    def testExtrudedOnDiffZLev1(self):
+        """Non regression bug : This test is base on P0P0 ExtrudedExtruded. This test checks that if the input meshes are not based on a same plane // OXY the interpolation works"""
+        arrX=DataArrayDouble([0,1]) ; arrY=DataArrayDouble([0,1]) ; arrZ=DataArrayDouble([0,1,2])
+        src=MEDCouplingCMesh() ; src.setCoords(arrX,arrY,arrZ)
+        arrX=DataArrayDouble([0.5,1.5]) ; arrY=DataArrayDouble([0.5,1.5]) ; arrZ=DataArrayDouble([0.5,2])
+        trg=MEDCouplingCMesh() ; trg.setCoords(arrX,arrY,arrZ)
+        #
+        src=MEDCouplingMappedExtrudedMesh(src) ; trg=MEDCouplingMappedExtrudedMesh(trg)
+        pt1=src.getMesh2D().getCoords().getHiddenCppPointer() ; pt2=trg.getMesh2D().getCoords().getHiddenCppPointer()
+        #
+        rem=MEDCouplingRemapper()
+        rem.prepare(src,trg,"P0P0")
+        self.checkMatrix(rem.getCrudeMatrix(),[{0:0.125,1:0.25}],src.getNumberOfCells(),1e-12)
+        #
+        self.assertEqual(src.getMesh2D().getSpaceDimension(),3)
+        self.assertEqual(trg.getMesh2D().getSpaceDimension(),3)
+        self.assertEqual(src.getMesh2D().getCoords().getHiddenCppPointer(),pt1)
+        self.assertEqual(trg.getMesh2D().getCoords().getHiddenCppPointer(),pt2)
+        #
+        rem2=MEDCouplingRemapper()
+        rem2.setIntersectionType(Geometric2D)
+        rem2.prepare(src,trg,"P0P0")
+        self.checkMatrix(rem2.getCrudeMatrix(),[{0:0.125,1:0.25}],src.getNumberOfCells(),1e-12)
+        pass
+
+    def checkMatrix(self,mat1,mat2,nbCols,eps):
+        self.assertEqual(len(mat1),len(mat2))
+        for i in xrange(len(mat1)):
+            self.assertTrue(max(mat2[i].keys())<nbCols)
+            self.assertTrue(max(mat1[i].keys())<nbCols)
+            self.assertTrue(min(mat2[i].keys())>=0)
+            self.assertTrue(min(mat1[i].keys())>=0)
+            s1=set(mat1[i].keys()) ; s2=set(mat2[i].keys())
+            for elt in s1.intersection(s2):
+                self.assertTrue(abs(mat1[i][elt]-mat2[i][elt])<eps)
+                pass
+            for elt in s1.difference(s2):
+                self.assertTrue(abs(mat1[i][elt])<eps)
+                pass
+            for elt in s2.difference(s1):
+                self.assertTrue(abs(mat2[i][elt])<eps)
+                pass
+            pass
+        pass
+    
     def build2DSourceMesh_1(self):
         sourceCoords=[-0.3,-0.3, 0.7,-0.3, -0.3,0.7, 0.7,0.7]
         sourceConn=[0,3,1,0,2,3]
-- 
2.39.2