Extend IMesh condensation to 3D and 1D.

[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingIMesh.cxx

diff --git a/src/MEDCoupling/MEDCouplingIMesh.cxx b/src/MEDCoupling/MEDCouplingIMesh.cxx
index 05e5875353daffd85717b511ac5a88ba34ea0ae8..b7f8e78cb050676c9ec4d01b3be3be4038533130 100644 (file)
--- a/src/MEDCoupling/MEDCouplingIMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingIMesh.cxx
@@ -164,7 +164,7 @@ MEDCouplingCMesh *MEDCouplingIMesh::convertToCartesian() const
   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCMesh> ret(MEDCouplingCMesh::New());
   try
   { ret->copyTinyInfoFrom(this); }
-  catch(INTERP_KERNEL::Exception& e) { }
+  catch(INTERP_KERNEL::Exception& ) { }
   int spaceDim(getSpaceDimension());
   std::vector<std::string> infos(buildInfoOnComponents());
   for(int i=0;i<spaceDim;i++)
@@ -222,39 +222,87 @@ void MEDCouplingIMesh::CondenseFineToCoarse(DataArrayDouble *coarseDA, const std
       throw INTERP_KERNEL::Exception(oss.str().c_str());
     }
   int nbTuplesFine(fineDA->getNumberOfTuples());
-  if(nbTuplesFine%nbOfTuplesInCoarseExp!=0)
+  if(nbTuplesFine%nbOfTuplesInFineExp!=0)
     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : Invalid nb of tuples in fine DataArray regarding its structure !");
   int factN(nbTuplesFine/nbOfTuplesInFineExp);
   int fact(FindIntRoot(factN,meshDim));
   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids(BuildExplicitIdsFrom(coarseSt,fineLocInCoarse));
-  const int *idsPtr(ids->begin());
   double *outPtr(coarseDA->getPointer());
   const double *inPtr(fineDA->begin());
-  coarseDA->setPartOfValuesSimple3(0.,ids->begin(),ids->end(),0,nbCompo,1);
   //
+  std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
   switch(meshDim)
   {
+    case 1:
+      {
+        int offset(fineLocInCoarse[0].first);
+        for(int i=0;i<dims[0];i++)
+          {
+            double *loc(outPtr+(offset+i)*nbCompo);
+            for(int ifact=0;ifact<fact;ifact++,inPtr+=nbCompo)
+              {
+                if(ifact!=0)
+                  std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
+                else
+                  std::copy(inPtr,inPtr+nbCompo,loc);
+              }
+          }
+        break;
+      }
     case 2:
       {
-        int kk(0);
-        std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
-        for(int it=0;it<dims[1];it++)
+        int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first);
+        for(int j=0;j<dims[1];j++)
+          {
+            for(int jfact=0;jfact<fact;jfact++)
+              {
+                for(int i=0;i<dims[0];i++)
+                  {
+                    double *loc(outPtr+(kk+i)*nbCompo);
+                    for(int ifact=0;ifact<fact;ifact++,inPtr+=nbCompo)
+                      {
+                        if(jfact!=0 || ifact!=0)
+                          std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
+                        else
+                          std::copy(inPtr,inPtr+nbCompo,loc);
+                      }
+                  }
+              }
+            kk+=coarseSt[0];
+          }
+        break;
+      }
+    case 3:
+      {
+        int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first+coarseSt[0]*coarseSt[1]*fineLocInCoarse[2].first);
+        for(int k=0;k<dims[2];k++)
           {
-            for(int i=0;i<fact;i++)
+            for(int kfact=0;kfact<fact;kfact++)
               {
-                for(int j=0;j<dims[0];j++,inPtr+=fact)
+                for(int j=0;j<dims[1];j++)
                   {
-                    double *loc(outPtr+idsPtr[kk+j]*nbCompo);
-                    std::transform(inPtr,inPtr+fact,loc,loc,std::plus<double>());
+                    for(int jfact=0;jfact<fact;jfact++)
+                      {
+                        for(int i=0;i<dims[0];i++)
+                          {
+                            double *loc(outPtr+(kk+i+j*coarseSt[0])*nbCompo);
+                            for(int ifact=0;ifact<fact;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);
+                              }
+                          }
+                      }
                   }
               }
-            kk+=it;
+            kk+=coarseSt[0]*coarseSt[1];
           }
         break;
       }
     default:
-      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : only dimensions 2 supported !");
+      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : only dimensions 1, 2 and 3 supported !");
   }
 }
 
@@ -558,7 +606,7 @@ int MEDCouplingIMesh::getCellContainingPoint(const double *pos, double eps) cons
     {
       int nbOfCells(_structure[i]-1);
       double ref(pos[i]);
-      int tmp((ref-_origin[i])/_dxyz[i]);
+      int tmp((int)((ref-_origin[i])/_dxyz[i]));
       if(tmp>=0 && tmp<nbOfCells)
         {
           ret+=coeff*tmp;
@@ -833,9 +881,12 @@ int MEDCouplingIMesh::FindIntRoot(int val, int order)
   else//order==3
     {
       double retf(std::pow(val,0.3333333333333333));
-      int ret((int)round(retf));
-      if(ret*ret*ret!=val)
+      int ret((int)retf),ret2(ret+1);
+      if(ret*ret*ret!=val && ret2*ret2*ret2!=val)
         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect cublic root !");
-      return ret;
+      if(ret*ret*ret==val)
+        return ret;
+      else
+        return ret2;
     }
 }