From 0a68c72a271c88d7d9f558de71acf3940bd1dd4f Mon Sep 17 00:00:00 2001
From: Anthony Geay <anthony.geay@edf.fr>
Date: Mon, 6 Feb 2017 08:52:20 +0100
Subject: [PATCH] Management of 3D surf voronoi for Gauss

---
 src/MEDCoupling/MEDCouplingFieldDouble.cxx    | 15 +++----
 .../MEDCouplingGaussLocalization.cxx          |  8 ++--
 .../MEDCouplingBasicsTest5.py                 | 41 +++++++++++++++++++
 3 files changed, 51 insertions(+), 13 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
index 96e977352..2b51bd607 100644
--- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
@@ -2212,16 +2212,13 @@ MCAuto<MEDCouplingFieldDouble> MEDCouplingFieldDouble::voronoize(double eps) con
   checkConsistencyLight();
   const MEDCouplingMesh *mesh(getMesh());
   INTERP_KERNEL::AutoCppPtr<Voronizer> vor;
-  if(mesh->getSpaceDimension()==2 && mesh->getSpaceDimension()==2)
-    {
-      vor=new Voronizer2D;
-    }
-  else if(mesh->getSpaceDimension()==3 && mesh->getSpaceDimension()==3)
-    {
-      vor=new Voronizer3D;
-    }
+  int meshDim(mesh->getMeshDimension()),spaceDim(mesh->getSpaceDimension());
+  if(meshDim==2 && (spaceDim==2 || spaceDim==3))
+    vor=new Voronizer2D;
+  else if(meshDim==3 && spaceDim==3)
+    vor=new Voronizer3D;
   else
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::voronoize : only 2D and 3D are supported for the moment !");
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::voronoize : only 2D, 3D surf, and 3D are supported for the moment !");
   return voronoizeGen(vor,eps);
 }
 
diff --git a/src/MEDCoupling/MEDCouplingGaussLocalization.cxx b/src/MEDCoupling/MEDCouplingGaussLocalization.cxx
index f75cd82a8..6c69fa0c6 100644
--- a/src/MEDCoupling/MEDCouplingGaussLocalization.cxx
+++ b/src/MEDCoupling/MEDCouplingGaussLocalization.cxx
@@ -218,9 +218,9 @@ MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::localizePtsInRefCooForEach
   //
   int nbPts(ptsInRefCoo->getNumberOfTuples());
   INTERP_KERNEL::NormalizedCellType typ(getType());
-  int dim(INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension());
+  int dim(INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension()),outDim(mesh->getSpaceDimension());
   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
-  ret->alloc(nbPts*nbCells,dim);
+  ret->alloc(nbPts*nbCells,outDim);
   double *retPtr(ret->getPointer());
   if(dim!=ptsInRefCoo->getNumberOfComponents())
     throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : number of components of input coo is not equal to dim of element !");
@@ -228,8 +228,8 @@ MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::localizePtsInRefCooForEach
   INTERP_KERNEL::GaussCoords calculator;
   calculator.addGaussInfo(typ,dim, ptsInRefCoo->begin(),nbPts,&_ref_coord[0],getNumberOfPtsInRefCell());
   //
-  for(int i=0;i<nbCells;i++,retPtr+=nbPts*dim)
-    calculator.calculateCoords(getType(),coords,dim,conn+connI[i]+1,retPtr);
+  for(int i=0;i<nbCells;i++,retPtr+=nbPts*outDim)
+    calculator.calculateCoords(getType(),coords,outDim,conn+connI[i]+1,retPtr);
   return ret;
 }
 
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py
index 3f5bdda99..3f3af10cb 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py
@@ -4720,6 +4720,47 @@ class MEDCouplingBasicsTest5(unittest.TestCase):
         #
         self.assertEqual(fieldOnCell.getMesh().buildDescendingConnectivity()[0].getNumberOfCells(),2*7+22)# last little test to reduce chance of errors. For humans there 21 but last tiny edge is split into 2 subedges due to alg
         pass
+
+    def testVoronoi3DSurf_1(self):
+        tmp=MEDCouplingCMesh("mesh")
+        arr=DataArrayDouble(5) ; arr.iota()
+        tmp.setCoords(arr,arr)
+        tmp=tmp.build1SGTUnstructured()
+        conn=tmp.getNodalConnectivity()
+        conn.rearrange(4)
+        conn.reversePerTuple()
+        conn.circularPermutationPerTuple(2)
+        conn.rearrange(1)
+        coo=tmp.getCoords().deepCopy()
+        coo.circularPermutationPerTuple(2) ; coo*=0.1
+        coo.reverse()
+        coo2=DataArrayDouble(len(tmp.getCoords())*tmp.getSpaceDimension()) ; coo2.iota() ; coo2.rearrange(tmp.getSpaceDimension())
+        coo2*=0.14
+        coo2.circularPermutationPerTuple(2)
+        tmp.getCoords()[:]+=coo2*coo
+        #
+        tmp.changeSpaceDimension(3,0.)    # force 3D surf
+        tmp.rotate([0,0,0],[1,0,0],pi/3)  # force 3D surf
+        #
+        field=MEDCouplingFieldDouble(ON_GAUSS_PT)
+        field.setName("MyFieldPG") ; field.setMesh(tmp)
+        field.setGaussLocalizationOnType(NORM_QUAD4,[-1.,-1.,1.,-1.,1.,1.,-1.,1.],[0.8,-0.8, 0.8,0.8, -0.8,0.8, -0.8,-0.8, 0.,0., 0.2,0.2, 0.1,0.3],[0.1,0.1,0.1,0.1,0.1,0.1,0.4])
+        arr=DataArrayDouble(field.getNumberOfTuplesExpected()) ; arr.iota() 
+        field.setArray(arr)
+        field.checkConsistencyLight()
+        #####
+        fieldOnCell=field.voronoize(1e-12);
+        fieldOnCell.checkConsistencyLight()
+        self.assertEqual(field.getMesh().getSpaceDimension(),fieldOnCell.getMesh().getSpaceDimension())
+        self.assertTrue(fieldOnCell.getArray().isEqual(field.getArray(),1e-12))
+        meaRef=field.getMesh().getMeasureField(True).getArray()
+        mea=fieldOnCell.getMesh().getMeasureField(True).getArray()
+        self.assertEqual(field.getDiscretization().getNbOfGaussLocalization(),1)
+        self.assertEqual(field.getDiscretization().getGaussLocalization(0).getNumberOfGaussPt(),7)
+        mea.rearrange(7)
+        mea2=mea.sumPerTuple()
+        self.assertTrue(mea2.isEqual(meaRef,1e-12))
+        pass
     
     pass
 
-- 
2.39.2