From e83dc6e4809ed462e945ad3c199ad415c9a8b4a3 Mon Sep 17 00:00:00 2001
From: geay <anthony.geay@cea.fr>
Date: Tue, 17 Jun 2014 17:25:00 +0200
Subject: [PATCH] Add 3D algorithms for Spread and Condense with ghosts.

---
 src/MEDCoupling/MEDCouplingIMesh.cxx          | 151 ++++++++++++------
 src/MEDCoupling/MEDCouplingIMesh.hxx          |   1 +
 src/MEDCoupling_Swig/MEDCouplingBasicsTest.py |  31 ++++
 3 files changed, 137 insertions(+), 46 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingIMesh.cxx b/src/MEDCoupling/MEDCouplingIMesh.cxx
index f5821c9e6..d26364a52 100644
--- a/src/MEDCoupling/MEDCouplingIMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingIMesh.cxx
@@ -484,8 +484,44 @@ void MEDCouplingIMesh::CondenseFineToCoarseGhost(const std::vector<int>& coarseS
           }
         break;
       }
+    case 3:
+      {
+        int nxwg(coarseSt[0]+2*ghostSize),nxywg((coarseSt[0]+2*ghostSize)*(coarseSt[1]+2*ghostSize));
+        int kk(fineLocInCoarse[0].first+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize)+nxywg*(fineLocInCoarse[2].first+ghostSize)),fact2(facts[2]),fact1(facts[1]),fact0(facts[0]);
+        inPtr+=(dims[0]*fact0+2*ghostSize)*(dims[1]*fact1+2*ghostSize)*ghostSize*nbCompo;
+        for(int k=0;k<dims[2];k++)
+          {
+            for(int kfact=0;kfact<fact2;kfact++)
+              {
+                inPtr+=ghostSize*(dims[0]*fact0+2*ghostSize)*nbCompo;
+                for(int j=0;j<dims[1];j++)
+                  {
+                    int kky(j*nxwg);
+                    for(int jfact=0;jfact<fact1;jfact++)
+                      {
+                        inPtr+=ghostSize*nbCompo;
+                        for(int i=0;i<dims[0];i++)
+                          {
+                            double *loc(outPtr+(kky+kk+i)*nbCompo);
+                            for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
+                              {
+                                if(kfact!=0 || jfact!=0 || ifact!=0)
+                                  std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
+                                else
+                                  std::copy(inPtr,inPtr+nbCompo,loc);
+                              }
+                          }
+                        inPtr+=ghostSize*nbCompo;
+                      }
+                  }
+                inPtr+=ghostSize*(dims[0]*fact0+2*ghostSize)*nbCompo;
+              }
+            kk+=nxywg;
+          }
+        break;
+      }
     default:
-      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : only dimensions 1, 2 supported !");
+      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : only dimensions 1, 2, 3 supported !");
   }
 }
 
@@ -636,11 +672,11 @@ void MEDCouplingIMesh::SpreadCoarseToFineGhost(const DataArrayDouble *coarseDA,
   double *outPtr(fineDA->getPointer());
   const double *inPtr(coarseDA->begin());
   //
-  std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
   switch(meshDim)
   {
     case 1:
       {
+        std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
         int offset(fineLocInCoarse[0].first+ghostSize-1),fact0(facts[0]);//offset is always >=0 thanks to the fact that ghostSize>=1 !
         for(int i=0;i<ghostSize;i++)
           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
@@ -658,53 +694,27 @@ void MEDCouplingIMesh::SpreadCoarseToFineGhost(const DataArrayDouble *coarseDA,
       }
     case 2:
       {
-        int nxwg(coarseSt[0]+2*ghostSize),fact0(facts[0]),fact1(facts[1]);
-        int kk(fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].first+ghostSize-1));//kk is always >=0 thanks to the fact that ghostSize>=1 !
-        for(int jg=0;jg<ghostSize;jg++)
-          {
-            for(int ig=0;ig<ghostSize;ig++)
-              outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
-            int kk0(kk+1);
-            for(int ig=0;ig<dims[0];ig++,kk0++)
-              for(int ifact=0;ifact<fact0;ifact++)
-                outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
-            for(int ik=0;ik<ghostSize;ik++)
-              outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
-          }
-        for(int j=0;j<dims[1];j++)
-          {
-            kk=fineLocInCoarse[0].first-1+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize+j);
-            for(int jfact=0;jfact<fact1;jfact++)
-              {
-                for(int ig=0;ig<ghostSize;ig++)
-                  outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
-                int kk0(kk+1);//1 not ghost. We make the hypothesis that factors is >= ghostlev
-                for(int i=0;i<dims[0];i++,kk0++)
-                  {
-                    const double *loc(inPtr+kk0*nbCompo);
-                    for(int ifact=0;ifact<fact0;ifact++)
-                      outPtr=std::copy(loc,loc+nbCompo,outPtr);
-                  }
-                for(int ig=0;ig<ghostSize;ig++)
-                  outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
-              }
-          }
-        kk=fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].second+ghostSize);
-        for(int jg=0;jg<ghostSize;jg++)
-          {
-            for(int ig=0;ig<ghostSize;ig++)
-              outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
-            int kk0(kk+1);
-            for(int ig=0;ig<dims[0];ig++,kk0++)
-              for(int ifact=0;ifact<fact0;ifact++)
-                outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
-            for(int ik=0;ik<ghostSize;ik++)
-              outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
-          }
+        SpreadCoarseToFineGhost2D(inPtr,outPtr,nbCompo,coarseSt,fineLocInCoarse,facts,ghostSize);
+        break;
+      }
+    case 3:
+      {
+        std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
+        int fact0(facts[0]),fact1(facts[1]),fact2(facts[2]);
+        int nxyWgCoarse((coarseSt[0]+2*ghostSize)*(coarseSt[1]+2*ghostSize)),nxyWgFine((dims[0]*fact0+2*ghostSize)*(dims[1]*fact1+2*ghostSize));
+        int offset((fineLocInCoarse[2].first+ghostSize-1)*nxyWgCoarse);//offset is always >=0 thanks to the fact that ghostSize>=1 !
+        for(int i=0;i<ghostSize;i++,outPtr+=nxyWgFine*nbCompo)
+          SpreadCoarseToFineGhost2D(inPtr+offset*nbCompo,outPtr,nbCompo,coarseSt,fineLocInCoarse,facts,ghostSize);
+        offset+=nxyWgCoarse;
+        for(int i=0;i<dims[2];i++,offset+=nxyWgCoarse)
+          for(int j=0;j<fact2;j++,outPtr+=nxyWgFine*nbCompo)
+            SpreadCoarseToFineGhost2D(inPtr+offset*nbCompo,outPtr,nbCompo,coarseSt,fineLocInCoarse,facts,ghostSize);
+        for(int i=0;i<ghostSize;i++,outPtr+=nxyWgFine*nbCompo)
+          SpreadCoarseToFineGhost2D(inPtr+offset*nbCompo,outPtr,nbCompo,coarseSt,fineLocInCoarse,facts,ghostSize);
         break;
       }
     default:
-      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : only dimensions 1, 2 supported !");
+      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : only dimensions 1, 2, 3 supported !");
   }
 }
 
@@ -1407,3 +1417,52 @@ int MEDCouplingIMesh::FindIntRoot(int val, int order)
         return ret2;
     }
 }
+
+void MEDCouplingIMesh::SpreadCoarseToFineGhost2D(const double *inPtr, double *outPtr, int nbCompo, const std::vector<int>& coarseSt, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts, int ghostSize)
+{
+  double *outPtrWork(outPtr);
+  std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
+  int nxwg(coarseSt[0]+2*ghostSize),fact0(facts[0]),fact1(facts[1]);
+  int kk(fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].first+ghostSize-1));//kk is always >=0 thanks to the fact that ghostSize>=1 !
+  for(int jg=0;jg<ghostSize;jg++)
+    {
+      for(int ig=0;ig<ghostSize;ig++)
+        outPtrWork=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtrWork);
+      int kk0(kk+1);
+      for(int ig=0;ig<dims[0];ig++,kk0++)
+        for(int ifact=0;ifact<fact0;ifact++)
+          outPtrWork=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtrWork);
+      for(int ik=0;ik<ghostSize;ik++)
+        outPtrWork=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtrWork);
+    }
+  for(int j=0;j<dims[1];j++)
+    {
+      kk=fineLocInCoarse[0].first-1+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize+j);
+      for(int jfact=0;jfact<fact1;jfact++)
+        {
+          for(int ig=0;ig<ghostSize;ig++)
+            outPtrWork=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtrWork);
+          int kk0(kk+1);//1 not ghost. We make the hypothesis that factors is >= ghostlev
+          for(int i=0;i<dims[0];i++,kk0++)
+            {
+              const double *loc(inPtr+kk0*nbCompo);
+              for(int ifact=0;ifact<fact0;ifact++)
+                outPtrWork=std::copy(loc,loc+nbCompo,outPtrWork);
+            }
+          for(int ig=0;ig<ghostSize;ig++)
+            outPtrWork=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtrWork);
+        }
+    }
+  kk=fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].second+ghostSize);
+  for(int jg=0;jg<ghostSize;jg++)
+    {
+      for(int ig=0;ig<ghostSize;ig++)
+        outPtrWork=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtrWork);
+      int kk0(kk+1);
+      for(int ig=0;ig<dims[0];ig++,kk0++)
+        for(int ifact=0;ifact<fact0;ifact++)
+          outPtrWork=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtrWork);
+      for(int ik=0;ik<ghostSize;ik++)
+        outPtrWork=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtrWork);
+    }
+}
diff --git a/src/MEDCoupling/MEDCouplingIMesh.hxx b/src/MEDCoupling/MEDCouplingIMesh.hxx
index 0fbc783d1..5a203700e 100644
--- a/src/MEDCoupling/MEDCouplingIMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingIMesh.hxx
@@ -110,6 +110,7 @@ namespace ParaMEDMEM
     void checkSpaceDimension() const;
     static void CheckSpaceDimension(int spaceDim);
     static int FindIntRoot(int val, int order);
+    static void SpreadCoarseToFineGhost2D(const double *inPtr, double *outPtr, int nbCompo, const std::vector<int>& coarseSt, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts, int ghostSize);
   private:
     int _space_dim;
     double _origin[3];
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index 68e88ad55..36a7125cd 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -15495,6 +15495,37 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertTrue(att5.getFieldOn(att5.getMyGodFather().getMeshAtPosition((1,3)),"YY").isEqualWithoutConsideringStr(exp25,1e-12))
         pass
 
+    def testSwig2AMR11(self):
+        """ Some tests in 3D with CondenseFineToCoarseGhost and SpreadCoarseToFineGhost"""
+        coarse=DataArrayDouble((6+4)*(7+4)*(5+4)) ; coarse.iota()
+        fine=DataArrayDouble((4*2+4)*(2*3+4)*(3*4+4))
+        MEDCouplingIMesh.SpreadCoarseToFineGhost(coarse,[6,7,5],fine,[(1,5),(2,4),(1,4)],[2,3,4],2)
+        exp0=DataArrayDouble([252.,252.,253.,253.,254.,254.,255.,255.,256.,256.,257.,257.,252.,252.,253.,253.,254.,254.,255.,255.,256.,256.,257.,257.,262.,262.,263.,263.,264.,264.,265.,265.,266.,266.,267.,267.,262.,262.,263.,263.,264.,264.,265.,265.,266.,266.,267.,267.,262.,262.,263.,263.,264.,264.,265.,265.,266.,266.,267.,267.,272.,272.,273.,273.,274.,274.,275.,275.,276.,276.,277.,277.,272.,272.,273.,273.,274.,274.,275.,275.,276.,276.,277.,277.,272.,272.,273.,273.,274.,274.,275.,275.,276.,276.,277.,277.,282.,282.,283.,283.,284.,284.,285.,285.,286.,286.,287.,287.,282.,282.,283.,283.,284.,284.,285.,285.,286.,286.,287.,287.])
+        exp1=DataArrayDouble([362.,362.,363.,363.,364.,364.,365.,365.,366.,366.,367.,367.,362.,362.,363.,363.,364.,364.,365.,365.,366.,366.,367.,367.,372.,372.,373.,373.,374.,374.,375.,375.,376.,376.,377.,377.,372.,372.,373.,373.,374.,374.,375.,375.,376.,376.,377.,377.,372.,372.,373.,373.,374.,374.,375.,375.,376.,376.,377.,377.,382.,382.,383.,383.,384.,384.,385.,385.,386.,386.,387.,387.,382.,382.,383.,383.,384.,384.,385.,385.,386.,386.,387.,387.,382.,382.,383.,383.,384.,384.,385.,385.,386.,386.,387.,387.,392.,392.,393.,393.,394.,394.,395.,395.,396.,396.,397.,397.,392.,392.,393.,393.,394.,394.,395.,395.,396.,396.,397.,397.])
+        exp2=DataArrayDouble([472.,472.,473.,473.,474.,474.,475.,475.,476.,476.,477.,477.,472.,472.,473.,473.,474.,474.,475.,475.,476.,476.,477.,477.,482.,482.,483.,483.,484.,484.,485.,485.,486.,486.,487.,487.,482.,482.,483.,483.,484.,484.,485.,485.,486.,486.,487.,487.,482.,482.,483.,483.,484.,484.,485.,485.,486.,486.,487.,487.,492.,492.,493.,493.,494.,494.,495.,495.,496.,496.,497.,497.,492.,492.,493.,493.,494.,494.,495.,495.,496.,496.,497.,497.,492.,492.,493.,493.,494.,494.,495.,495.,496.,496.,497.,497.,502.,502.,503.,503.,504.,504.,505.,505.,506.,506.,507.,507.,502.,502.,503.,503.,504.,504.,505.,505.,506.,506.,507.,507.])
+        exp3=DataArrayDouble([582.,582.,583.,583.,584.,584.,585.,585.,586.,586.,587.,587.,582.,582.,583.,583.,584.,584.,585.,585.,586.,586.,587.,587.,592.,592.,593.,593.,594.,594.,595.,595.,596.,596.,597.,597.,592.,592.,593.,593.,594.,594.,595.,595.,596.,596.,597.,597.,592.,592.,593.,593.,594.,594.,595.,595.,596.,596.,597.,597.,602.,602.,603.,603.,604.,604.,605.,605.,606.,606.,607.,607.,602.,602.,603.,603.,604.,604.,605.,605.,606.,606.,607.,607.,602.,602.,603.,603.,604.,604.,605.,605.,606.,606.,607.,607.,612.,612.,613.,613.,614.,614.,615.,615.,616.,616.,617.,617.,612.,612.,613.,613.,614.,614.,615.,615.,616.,616.,617.,617.])
+        exp4=DataArrayDouble([692.,692.,693.,693.,694.,694.,695.,695.,696.,696.,697.,697.,692.,692.,693.,693.,694.,694.,695.,695.,696.,696.,697.,697.,702.,702.,703.,703.,704.,704.,705.,705.,706.,706.,707.,707.,702.,702.,703.,703.,704.,704.,705.,705.,706.,706.,707.,707.,702.,702.,703.,703.,704.,704.,705.,705.,706.,706.,707.,707.,712.,712.,713.,713.,714.,714.,715.,715.,716.,716.,717.,717.,712.,712.,713.,713.,714.,714.,715.,715.,716.,716.,717.,717.,712.,712.,713.,713.,714.,714.,715.,715.,716.,716.,717.,717.,722.,722.,723.,723.,724.,724.,725.,725.,726.,726.,727.,727.,722.,722.,723.,723.,724.,724.,725.,725.,726.,726.,727.,727.])
+        exp=DataArrayDouble.Aggregate([exp0,exp0,exp1,exp1,exp1,exp1,exp2,exp2,exp2,exp2,exp3,exp3,exp3,exp3,exp4,exp4])
+        self.assertTrue(fine.isEqual(exp,1e-12))
+        fine.iota()
+        coarse.iota(0.5)
+        MEDCouplingIMesh.CondenseFineToCoarseGhost([6,7,5],fine,[(1,5),(2,4),(1,4)],[2,3,4],coarse,2)
+        amr=MEDCouplingCartesianAMRMesh("mesh",3,[7,8,6],[0.,0.,0.],[1.,1.,1.])
+        amr.addPatch([(1,5),(2,4),(1,4)],[2,3,4])
+        att=MEDCouplingAMRAttribute(amr,[("YY",1)],2)
+        att.alloc()
+        exp1=DataArrayDouble(990) ; exp1.iota(0.5)
+        ids=DataArrayInt([373,374,375,376,383,384,385,386,483,484,485,486,493,494,495,496,593,594,595,596,603,604,605,606])
+        vals=DataArrayDouble([11004.,11052.,11100.,11148.,11868.,11916.,11964.,12012.,22524.,22572.,22620.,22668.,23388.,23436.,23484.,23532.,34044.,34092.,34140.,34188.,34908.,34956.,35004.,35052.])
+        exp1[ids]=vals
+        self.assertTrue(coarse.isEqual(exp1,1e-12))
+        MEDCouplingStructuredMesh.MultiplyPartOf([10,11,9],[(3,7),(4,6),(3,6)],1/24.,coarse)
+        #
+        exp2=DataArrayDouble(990) ; exp2.iota(0.5)
+        exp2[ids]=vals/24.
+        self.assertTrue(coarse.isEqual(exp2,1e-12))
+        pass
+
     pass
 
 if __name__ == '__main__':
-- 
2.39.2