From f1d6e5ac6e52d8d826ccbef46316dffe7efe35ba Mon Sep 17 00:00:00 2001
From: geay <anthony.geay@cea.fr>
Date: Thu, 22 May 2014 14:39:14 +0200
Subject: [PATCH] some improvements and on the road of 3D manangement

---
 .../MEDCouplingCartesianAMRMesh.cxx           | 45 +++++-----
 src/MEDCoupling/MEDCouplingStructuredMesh.cxx | 89 ++++++++++---------
 src/MEDCoupling/MEDCouplingStructuredMesh.hxx |  1 -
 src/MEDCoupling_Swig/MEDCouplingBasicsTest.py |  4 +
 4 files changed, 70 insertions(+), 69 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx
index d97d294cf..b73244d43 100644
--- a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx
@@ -177,6 +177,7 @@ class InternalPatch : public RefCountObjectOnly
 {
 public:
   InternalPatch():_nb_of_true(0) { }
+  int getDimension() const { return (int)_part.size(); }
   double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
   int getNumberOfCells() const { return (int)_crit.size(); }
   void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
@@ -261,7 +262,7 @@ void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const Intern
     }
   int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
   int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
-  std::vector<double > ratio_inner(dimRatioInner);
+  std::vector<double> ratio_inner(dimRatioInner);
   double minRatio(1.e10);
   for(int i=0; i<dimRatioInner; i++)
     {
@@ -274,21 +275,17 @@ void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const Intern
   cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
 }
 
-void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, const int axisId, bool& cutFound, int& cutPlace)
+void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace)
 {
   cutPlace=-1; cutFound=false;
   int minCellDirection(bso.getMinCellDirection());
-
-  int sortedDims[2];
-  sortedDims[0]=axisId==0?1:0;
-  sortedDims[1]=axisId==0?0:1;
-  static const int dim(2);
+  const int dim(patchToBeSplit->getDimension());
   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
   for(int id=0;id<dim;id++)
     {
-      const std::vector<int>& signature(signatures[sortedDims[id]]);
+      const std::vector<int>& signature(signatures[id]);
       std::vector<int> hole;
-      std::vector<double> distance ;
+      std::vector<double> distance;
       int len((int)signature.size());
       for(int i=0;i<len;i++)
         if(signature[i]==0)
@@ -296,14 +293,14 @@ void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch
             hole.push_back(i);
       if (hole.size()>0)
         {
-          double center(((double)len/2.)+0.5);
+          double center(((double)len/2.));
           for(std::size_t i=0;i<hole.size();i++)
-            distance.push_back(fabs(hole[i]+1+0.5-center));//anouar ! why 0.5 ?
+            distance.push_back(fabs(hole[i]+1.-center));
 
           double distanceMin=*std::min_element(distance.begin(),distance.end());
           int posDistanceMin=std::find(distance.begin(),distance.end(),distanceMin)-distance.begin()-1;
-
           cutFound=true;
+          axisId=id;
           cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
           return ;
         }
@@ -312,14 +309,14 @@ void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch
 
 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
 {
-  cutFound=false; cutPlace=-1; axisId=-1;
+  cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
 
   const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
   int sign,minCellDirection(bso.getMinCellDirection());
-  static const int dim = 2 ;
+  const int dim(patchToBeSplit->getDimension());
 
-  std::vector<int> zeroCrossDims(2,-1);
-  std::vector<int> zeroCrossVals(2,-1);
+  std::vector<int> zeroCrossDims(dim,-1);
+  std::vector<int> zeroCrossVals(dim,-1);
   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
   for (int id=0;id<dim;id++)
     {
@@ -357,9 +354,9 @@ void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const Interna
             if (edge[i]==max_cross)
               max_cross_list.push_back(zero_cross[i]+1) ;
 
-          float center = (signature.size()/2.0)+0.5;
+          double center((signature.size()/2.0));
           for (unsigned int i=0;i<max_cross_list.size();i++)
-            distance.push_back(fabs(max_cross_list[i]+1+0.5-center));
+            distance.push_back(fabs(max_cross_list[i]+1-center));
 
           float distance_min=*min_element(distance.begin(),distance.end()) ;
           int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
@@ -403,7 +400,7 @@ void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const Interna
 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
 {
   cutFound=false;
-  if(patchToBeSplit->getEfficiency() <= bso.getEffeciencySnd())
+  if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd())
     {
       if(rangeOfAxisId>=2*bso.getMinCellDirection())
         {
@@ -469,19 +466,19 @@ void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL
       for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
         {
           //
-          int axisId,rangeOfAxisId;
+          int axisId,rangeOfAxisId,cutPlace;
           bool cutFound;
-          int cutPlace;
           MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
           if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
             { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
-          FindHole(bso,*it,axisId,cutFound,cutPlace);
+          FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
           if(cutFound)
             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
           FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
           if(cutFound)
             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
           TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
+          std::cerr << axisId << std::endl;
           if(cutFound)
             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
           listOfPatchesOK.push_back(DealWithNoCut(*it));
@@ -489,9 +486,7 @@ void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL
       listOfPatches=listOfPatchesTmp;
     }
   for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
-    {
-      addPatch((*it)->getConstPart(),factors);
-    }
+    addPatch((*it)->getConstPart(),factors);
 }
 
 void MEDCouplingCartesianAMRMesh::removePatch(int patchId)
diff --git a/src/MEDCoupling/MEDCouplingStructuredMesh.cxx b/src/MEDCoupling/MEDCouplingStructuredMesh.cxx
index 19c054b6e..8db11f5a8 100644
--- a/src/MEDCoupling/MEDCouplingStructuredMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingStructuredMesh.cxx
@@ -808,8 +808,53 @@ std::vector< std::vector<int> > MEDCouplingStructuredMesh::ComputeSignaturePerAx
           }
         break;
       }
+    case 3:
+      {
+        int nx(st[0]),ny(st[1]),nz(st[2]);
+        ret[0].resize(nx); ret[1].resize(ny); ret[2].resize(nz);
+        std::vector<int>& retX(ret[0]);
+        for(int i=0;i<nx;i++)
+          {
+            int cnt(0);
+            for(int k=0;k<nz;k++)
+              {
+                int offz(k*nx*ny+i);
+                for(int j=0;j<ny;j++)
+                  if(crit[offz+j*nx])
+                    cnt++;
+              }
+            retX[i]=cnt;
+          }
+        std::vector<int>& retY(ret[1]);
+        for(int j=0;j<ny;j++)
+          {
+            int cnt(0),offy(j*nx);
+            for(int k=0;k<nz;k++)
+              {
+                int offz(k*nx*ny+offy);
+                for(int i=0;i<nx;i++)
+                  if(crit[offz+i])
+                    cnt++;
+              }
+            retY[j]=cnt;
+          }
+        std::vector<int>& retZ(ret[2]);
+        for(int k=0;k<nz;k++)
+          {
+            int cnt(0),offz(k*nx*ny);
+            for(int j=0;j<ny;j++)
+              {
+                int offy(offz+j*nx);
+                for(int i=0;i<nx;i++)
+                  if(crit[offy+i])
+                    cnt++;
+              }
+            retZ[k]=cnt;
+          }
+        break;
+      }
     default:
-       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeSignatureOf : only dimensions 1, 2 are supported !");
+       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeSignatureOf : only dimensions 1, 2 and 3 are supported !");
   }
   return ret;
 }
@@ -982,48 +1027,6 @@ int MEDCouplingStructuredMesh::FindMinimalPartOf3D(const std::vector<int>& st, c
   return ret;
 }
 
-void MEDCouplingStructuredMesh::ExtractVecOfBool(const std::vector<int>& st, const std::vector<bool>& crit, const std::vector< std::pair<int,int> >& partCompactFormat, std::vector<bool>& reducedCrit)
-{
-  int nbt(DeduceNumberOfGivenRangeInCompactFrmt(partCompactFormat));
-  std::vector<int> dims(GetDimensionsFromCompactFrmt(partCompactFormat));
-  reducedCrit.resize(nbt);
-  switch((int)st.size())
-  {
-    case 1:
-      {
-        int nx(dims[0]);
-        int kk(partCompactFormat[0].first);
-        for(int i=0;i<nx;i++)
-          reducedCrit[i]=crit[kk+i];
-        break;
-      }
-    case 2:
-      {
-        int nx(dims[0]),ny(dims[1]);
-        int kk(partCompactFormat[0].first+partCompactFormat[1].first*nx),it(0);
-        for(int j=0;j<ny;j++,kk+=nx)
-          for(int i=0;i<nx;i++,it++)
-            reducedCrit[it]=crit[kk+i];
-        break;
-      }
-    case 3:
-      {
-        int nx(dims[0]),ny(dims[1]),nz(dims[2]);
-        int kk(partCompactFormat[0].first+partCompactFormat[1].first*nx+partCompactFormat[2].first*nx*ny),it(0);
-        for(int k=0;k<nz;k++,kk+=nx*ny)
-          for(int j=0;j<ny;j++)
-            {
-              int kk2(j*nx);
-              for(int i=0;i<nx;i++,it++)
-                reducedCrit[it]=crit[kk+kk2+i];
-            }
-        break;
-      }
-    default:
-      throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractVecOfBool : Only dimension 1, 2 and 3 are supported actually !");
-  }
-}
-
 /*!
  * This method computes given the nodal structure defined by [ \a nodeStBg , \a nodeStEnd ) the zipped form.
  * std::distance( \a nodeStBg, \a nodeStEnd ) is equal to the space dimension. The returned value is equal to
diff --git a/src/MEDCoupling/MEDCouplingStructuredMesh.hxx b/src/MEDCoupling/MEDCouplingStructuredMesh.hxx
index 1175ab37d..2cde046e1 100644
--- a/src/MEDCoupling/MEDCouplingStructuredMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingStructuredMesh.hxx
@@ -98,7 +98,6 @@ namespace ParaMEDMEM
     static int FindMinimalPartOf1D(const std::vector<int>& st, const std::vector<bool>& crit, std::vector< std::pair<int,int> >& partCompactFormat);
     static int FindMinimalPartOf2D(const std::vector<int>& st, const std::vector<bool>& crit, std::vector< std::pair<int,int> >& partCompactFormat);
     static int FindMinimalPartOf3D(const std::vector<int>& st, const std::vector<bool>& crit, std::vector< std::pair<int,int> >& partCompactFormat);
-    static void ExtractVecOfBool(const std::vector<int>& st, const std::vector<bool>& crit, const std::vector< std::pair<int,int> >& partCompactFormat, std::vector<bool>& reducedCrit);
   protected:
     static int ZipNodeStructure(const int *nodeStBg, const int *nodeStEnd, int zipNodeSt[3]);
   protected:
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index 530995096..3a0d610bb 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -15136,6 +15136,10 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         #
         self.assertEqual(MEDCouplingStructuredMesh.ChangeReferenceToGlobalOfCompactFrmt([(8,32),(4,17)],[(0,24),(2,12)]),[(8,32),(6,16)])
         self.assertEqual(MEDCouplingStructuredMesh.ChangeReferenceFromGlobalOfCompactFrmt([(8,32),(4,17)],[(8,32),(6,16)]),[(0,24),(2,12)])
+        self.assertTrue(amr.getImageMesh().isEqual(im,1e-12))
+        m=amr.getImageMesh().asSingleCell().build1SGTUnstructured()
+        self.assertTrue(m.getNodalConnectivity().isEqual(DataArrayInt([1,0,2,3])))
+        self.assertTrue(m.getCoords().isEqualWithoutConsideringStr(DataArrayDouble([(0,0),(2,0),(0,2),(2,2)]),1e-12))
         pass
 
     def setUp(self):
-- 
2.39.2