Salome HOME
On the road of the refinement for AMR mesh driven by a vector of bool.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingIMesh.cxx
index 6ef737b8ed5562ad42ef79d67a65af1972a0f1dd..142fac4b839f76d57871e907d253fd690aa290cb 100644 (file)
@@ -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++)
@@ -176,6 +176,251 @@ MEDCouplingCMesh *MEDCouplingIMesh::convertToCartesian() const
   return ret.retn();
 }
 
+/*!
+ * This method refines \a this uniformaly along all of its dimensions. In case of success the space covered by \a this will remain
+ * the same before the invocation except that the number of cells will be multiplied by \a factor ^ this->getMeshDimension().
+ * The origin of \a this will be not touched only spacing and node structure will be changed.
+ * This method can be useful for AMR users.
+ */
+void MEDCouplingIMesh::refineWithFactor(const std::vector<int>& factors)
+{
+  if((int)factors.size()!=_space_dim)
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::refineWithFactor : refinement factors must have size equal to spaceDim !");
+  checkCoherency();
+  std::vector<int> structure(_structure,_structure+3);
+  std::vector<double> dxyz(_dxyz,_dxyz+3);
+  for(int i=0;i<_space_dim;i++)
+    {
+      if(factors[i]<=0)
+        {
+          std::ostringstream oss; oss << "MEDCouplingIMesh::refineWithFactor : factor for axis #" << i << " (" << factors[i] << ")is invalid ! Must be > 0 !";
+          throw INTERP_KERNEL::Exception(oss.str().c_str());
+        }
+      int factAbs(std::abs(factors[i]));
+      double fact2(1./(double)factors[i]);
+      structure[i]=(_structure[i]-1)*factAbs+1;
+      dxyz[i]=fact2*_dxyz[i];
+    }
+  std::copy(structure.begin(),structure.end(),_structure);
+  std::copy(dxyz.begin(),dxyz.end(),_dxyz);
+  declareAsNew();
+}
+
+/*!
+ * This static method is useful to condense field on cells of a MEDCouplingIMesh instance coming from a refinement ( MEDCouplingIMesh::refineWithFactor for example)
+ * to a coarse MEDCouplingIMesh instance. So this method can be seen as a specialization in P0P0 conservative interpolation non overlaping from fine image mesh
+ * to a coarse image mesh. Only tuples ( deduced from \a fineLocInCoarse ) of \a coarseDA will be modified. Other tuples of \a coarseDA will be let unchanged.
+ *
+ * \param [in,out] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
+ * \param [in] coarseSt The cell structure of coarse mesh.
+ * \param [in] fineDA The DataArray containing the cell field on uniformly refined mesh
+ * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
+ * \param [in] facts The refinement coefficient per axis.
+ * \sa SpreadCoarseToFine
+ */
+void MEDCouplingIMesh::CondenseFineToCoarse(DataArrayDouble *coarseDA, const std::vector<int>& coarseSt, const DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts)
+{
+  if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the parameters 1 or 3 are NULL or not allocated !");
+  int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
+  int nbCompo(fineDA->getNumberOfComponents());
+  if(coarseDA->getNumberOfComponents()!=nbCompo)
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the number of components of fine DA and coarse one mismatches !");
+  if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the size of fineLocInCoarse (4th param) and facts (5th param) must be equal to the sier of coarseSt (2nd param) !");
+  if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
+    {
+      std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
+  int nbTuplesFine(fineDA->getNumberOfTuples());
+  if(nbTuplesFine%nbOfTuplesInFineExp!=0)
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : Invalid nb of tuples in fine DataArray regarding its structure !");
+  int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
+  if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
+    {
+      std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
+  // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
+  double *outPtr(coarseDA->getPointer());
+  const double *inPtr(fineDA->begin());
+  //
+  std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
+  switch(meshDim)
+  {
+    case 1:
+      {
+        int offset(fineLocInCoarse[0].first),fact0(facts[0]);
+        for(int i=0;i<dims[0];i++)
+          {
+            double *loc(outPtr+(offset+i)*nbCompo);
+            for(int ifact=0;ifact<fact0;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(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact1(facts[1]),fact0(facts[0]);
+        for(int j=0;j<dims[1];j++)
+          {
+            for(int jfact=0;jfact<fact1;jfact++)
+              {
+                for(int i=0;i<dims[0];i++)
+                  {
+                    double *loc(outPtr+(kk+i)*nbCompo);
+                    for(int ifact=0;ifact<fact0;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),fact2(facts[2]),fact1(facts[1]),fact0(facts[0]);
+        for(int k=0;k<dims[2];k++)
+          {
+            for(int kfact=0;kfact<fact2;kfact++)
+              {
+                for(int j=0;j<dims[1];j++)
+                  {
+                    for(int jfact=0;jfact<fact1;jfact++)
+                      {
+                        for(int i=0;i<dims[0];i++)
+                          {
+                            double *loc(outPtr+(kk+i+j*coarseSt[0])*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);
+                              }
+                          }
+                      }
+                  }
+              }
+            kk+=coarseSt[0]*coarseSt[1];
+          }
+        break;
+      }
+    default:
+      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : only dimensions 1, 2 and 3 supported !");
+  }
+}
+
+/*!
+ * This method spreads the values of coarse data \a coarseDA into \a fineDA.
+ *
+ * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
+ * \param [in] coarseSt The cell structure of coarse mesh.
+ * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
+ * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
+ * \param [in] facts The refinement coefficient per axis.
+ * \sa CondenseFineToCoarse
+ */
+void MEDCouplingIMesh::SpreadCoarseToFine(const DataArrayDouble *coarseDA, const std::vector<int>& coarseSt, DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts)
+{
+  if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the parameters 1 or 3 are NULL or not allocated !");
+  int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
+  int nbCompo(fineDA->getNumberOfComponents());
+  if(coarseDA->getNumberOfComponents()!=nbCompo)
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the number of components of fine DA and coarse one mismatches !");
+  if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the size of fineLocInCoarse (4th param) and facts (5th param) must be equal to the sier of coarseSt (2nd param) !");
+  if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
+    {
+      std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
+  int nbTuplesFine(fineDA->getNumberOfTuples());
+  if(nbTuplesFine%nbOfTuplesInFineExp!=0)
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : Invalid nb of tuples in fine DataArray regarding its structure !");
+  int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
+  if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
+    {
+      std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
+  // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
+  double *outPtr(fineDA->getPointer());
+  const double *inPtr(coarseDA->begin());
+  //
+  std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
+  switch(meshDim)
+  {
+    case 1:
+      {
+        int offset(fineLocInCoarse[0].first),fact0(facts[0]);
+        for(int i=0;i<dims[0];i++)
+          {
+            const double *loc(inPtr+(offset+i)*nbCompo);
+            for(int ifact=0;ifact<fact0;ifact++)
+              outPtr=std::copy(loc,loc+nbCompo,outPtr);
+          }
+        break;
+      }
+    case 2:
+      {
+        int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact0(facts[0]),fact1(facts[1]);
+        for(int j=0;j<dims[1];j++)
+          {
+            for(int jfact=0;jfact<fact1;jfact++)
+              {
+                for(int i=0;i<dims[0];i++)
+                  {
+                    const double *loc(inPtr+(kk+i)*nbCompo);
+                    for(int ifact=0;ifact<fact0;ifact++)
+                      outPtr=std::copy(loc,loc+nbCompo,outPtr);
+                  }
+              }
+            kk+=coarseSt[0];
+          }
+        break;
+      }
+    case 3:
+      {
+        int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first+coarseSt[0]*coarseSt[1]*fineLocInCoarse[2].first),fact0(facts[0]),fact1(facts[2]),fact2(facts[2]);
+        for(int k=0;k<dims[2];k++)
+          {
+            for(int kfact=0;kfact<fact2;kfact++)
+              {
+                for(int j=0;j<dims[1];j++)
+                  {
+                    for(int jfact=0;jfact<fact1;jfact++)
+                      {
+                        for(int i=0;i<dims[0];i++)
+                          {
+                            const double *loc(inPtr+(kk+i+j*coarseSt[0])*nbCompo);
+                            for(int ifact=0;ifact<fact0;ifact++)
+                              outPtr=std::copy(loc,loc+nbCompo,outPtr);
+                          }
+                      }
+                  }
+              }
+            kk+=coarseSt[0]*coarseSt[1];
+          }
+        break;
+      }
+    default:
+      throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : only dimensions 1, 2 and 3 supported !");
+  }
+}
+
 void MEDCouplingIMesh::setSpaceDimension(int spaceDim)
 {
   if(spaceDim==_space_dim)
@@ -327,30 +572,6 @@ void MEDCouplingIMesh::checkCoherency2(double eps) const
   checkCoherency1(eps);
 }
 
-void MEDCouplingIMesh::getSplitCellValues(int *res) const
-{
-  int meshDim(getMeshDimension());
-  for(int l=0;l<meshDim;l++)
-    {
-      int val=1;
-      for(int p=0;p<meshDim-l-1;p++)
-        val*=_structure[p]-1;
-      res[meshDim-l-1]=val;
-    }
-}
-
-void MEDCouplingIMesh::getSplitNodeValues(int *res) const
-{
-  int spaceDim(getSpaceDimension());
-  for(int l=0;l<spaceDim;l++)
-    {
-      int val=1;
-      for(int p=0;p<spaceDim-l-1;p++)
-        val*=_structure[p];
-      res[spaceDim-l-1]=val;
-    }
-}
-
 void MEDCouplingIMesh::getNodeGridStructure(int *res) const
 {
   checkSpaceDimension();
@@ -500,7 +721,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;
@@ -670,15 +891,15 @@ void MEDCouplingIMesh::unserialization(const std::vector<double>& tinyInfoD, con
 void MEDCouplingIMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
 {
   checkCoherency();
-  std::ostringstream extent;
+  std::ostringstream extent,origin,spacing;
   for(int i=0;i<3;i++)
     {
       if(i<_space_dim)
-        { extent << "0 " <<  _structure[i]-1 << " "; }
+        { extent << "0 " <<  _structure[i]-1 << " "; origin << _origin[i] << " "; spacing << _dxyz[i] << " "; }
       else
-        { extent << "0 0 "; }
+        { extent << "0 0 "; origin << "0 "; spacing << "0 "; }
     }
-  ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
+  ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\" Origin=\"" << origin.str() << "\" Spacing=\"" << spacing.str() << "\">\n";
   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
   ofs << "      <PointData>\n" << pointData << std::endl;
   ofs << "      </PointData>\n";
@@ -753,3 +974,34 @@ void MEDCouplingIMesh::CheckSpaceDimension(int spaceDim)
     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CheckSpaceDimension : input spaceDim must be in [0,1,2,3] !");
 }
 
+int MEDCouplingIMesh::FindIntRoot(int val, int order)
+{
+  if(order==0)
+    return 1;
+  if(val<0)
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : input val is < 0 ! Not possible to compute a root !");
+  if(order==1)
+    return val;
+  if(order!=2 && order!=3)
+    throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the order available are 0,1,2 or 3 !");
+  double valf((double)val);
+  if(order==2)
+    {
+      double retf(sqrt(valf));
+      int ret((int)retf);
+      if(ret*ret!=val)
+        throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect square root !");
+      return ret;
+    }
+  else//order==3
+    {
+      double retf(std::pow(val,0.3333333333333333));
+      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 !");
+      if(ret*ret*ret==val)
+        return ret;
+      else
+        return ret2;
+    }
+}