Generalization of unstructured grid supported by the remapper.

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

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 1c99da1c2b1224b733b40056103d396f41649f54..a927a731c068e65c8c6dafb02d2c7831e300dc77 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -479,6 +479,41 @@ int MEDCouplingFieldDiscretizationP0::getNumberOfTuples(const MEDCouplingMesh *m
   return mesh->getNumberOfCells();
 }
 
+/*!
+ * mesh is not used here. It is not a bug !
+ */
+int MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
+{
+  if(code.size()%3!=0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
+  int nbOfSplit=(int)idsPerType.size();
+  int nbOfTypes=(int)code.size()/3;
+  int ret=0;
+  for(int i=0;i<nbOfTypes;i++)
+    {
+      int nbOfEltInChunk=code[3*i+1];
+      if(nbOfEltInChunk<0)
+        throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
+      int pos=code[3*i+2];
+      if(pos!=-1)
+        {
+          if(pos<0 || pos>=nbOfSplit)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+          const DataArrayInt *ids(idsPerType[pos]);
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+        }
+      ret+=nbOfEltInChunk;
+    }
+  return ret;
+}
+
 int MEDCouplingFieldDiscretizationP0::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
 {
   if(!mesh)
@@ -497,7 +532,7 @@ DataArrayInt *MEDCouplingFieldDiscretizationP0::getOffsetArr(const MEDCouplingMe
   return ret;
 }
 
-void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
+void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
                                                              const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 {
   if(!mesh)
@@ -505,13 +540,13 @@ void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMe
   const int *array=old2NewBg;
   if(check)
     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
-  for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
+  for(std::vector<DataArray *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
     {
       if(*it)
         (*it)->renumberInPlace(array);
     }
   if(check)
-    delete [] array;
+    free(const_cast<int *>(array));
 }
 
 DataArrayDouble *MEDCouplingFieldDiscretizationP0::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
@@ -543,10 +578,10 @@ void MEDCouplingFieldDiscretizationP0::checkCompatibilityWithNature(NatureOfFiel
 {
 }
 
-void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
+void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
 {
-  if(!mesh)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::checkCoherencyBetween : NULL input mesh !");
+  if(!mesh || !da)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::checkCoherencyBetween : NULL input mesh or DataArray !");
   if(mesh->getNumberOfCells()!=da->getNumberOfTuples())
     {
       std::ostringstream message;
@@ -684,6 +719,41 @@ int MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuples(const MEDCouplingMe
   return mesh->getNumberOfNodes();
 }
 
+/*!
+ * mesh is not used here. It is not a bug !
+ */
+int MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
+{
+  if(code.size()%3!=0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
+  int nbOfSplit=(int)idsPerType.size();
+  int nbOfTypes=(int)code.size()/3;
+  int ret=0;
+  for(int i=0;i<nbOfTypes;i++)
+    {
+      int nbOfEltInChunk=code[3*i+1];
+      if(nbOfEltInChunk<0)
+        throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
+      int pos=code[3*i+2];
+      if(pos!=-1)
+        {
+          if(pos<0 || pos>=nbOfSplit)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+          const DataArrayInt *ids(idsPerType[pos]);
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+        }
+      ret+=nbOfEltInChunk;
+    }
+  return ret;
+}
+
 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
 {
   if(!mesh)
@@ -694,7 +764,7 @@ int MEDCouplingFieldDiscretizationOnNodes::getNumberOfMeshPlaces(const MEDCoupli
 /*!
  * Nothing to do here.
  */
-void MEDCouplingFieldDiscretizationOnNodes::renumberArraysForCell(const MEDCouplingMesh *, const std::vector<DataArrayDouble *>& arrays,
+void MEDCouplingFieldDiscretizationOnNodes::renumberArraysForCell(const MEDCouplingMesh *, const std::vector<DataArray *>& arrays,
                                                                   const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 {
 }
@@ -732,10 +802,10 @@ void MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds(c
   trueTupleRestriction=ret2.retn();
 }
 
-void MEDCouplingFieldDiscretizationOnNodes::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
+void MEDCouplingFieldDiscretizationOnNodes::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
 {
-  if(!mesh)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::checkCoherencyBetween : NULL input mesh !");
+  if(!mesh || !da)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::checkCoherencyBetween : NULL input mesh or DataArray !");
   if(mesh->getNumberOfNodes()!=da->getNumberOfTuples())
     {
       std::ostringstream message;
@@ -789,7 +859,7 @@ MEDCouplingMesh *MEDCouplingFieldDiscretizationOnNodes::buildSubMeshDataRange(co
 
 /*!
  * This method returns a tuple ids selection from cell ids selection [start;end).
- * This method is called by MEDCouplingFieldDiscretizationP0::buildSubMeshData to return parameter \b di.
+ * This method is called by MEDCouplingFieldDiscretizationOnNodes::buildSubMeshData to return parameter \b di.
  * Here for P1 only nodes fetched by submesh of mesh[startCellIds:endCellIds) is returned !
  *
  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
@@ -1005,12 +1075,12 @@ std::size_t MEDCouplingFieldDiscretizationPerCell::getHeapMemorySize() const
   return ret;
 }
 
-void MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
+void MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
 {
   if(!_discr_per_cell)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has no discretization per cell !");
   if(!mesh)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween : NULL input mesh !");
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween : NULL input mesh or DataArray !");
   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
   if(nbOfTuples!=mesh->getNumberOfCells())
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has a discretization per cell but it's not matching the underlying mesh !");
@@ -1067,7 +1137,7 @@ void MEDCouplingFieldDiscretizationPerCell::renumberCells(const int *old2NewBg,
   _discr_per_cell=dpc;
   //
   if(check)
-    delete [] const_cast<int *>(array);
+    free(const_cast<int *>(array));
 }
 
 void MEDCouplingFieldDiscretizationPerCell::buildDiscrPerCellIfNecessary(const MEDCouplingMesh *mesh)
@@ -1248,6 +1318,47 @@ const char *MEDCouplingFieldDiscretizationGauss::getRepr() const
   return REPR;
 }
 
+/*!
+ * mesh is not used here. It is not a bug !
+ */
+int MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
+{
+  if(!_discr_per_cell || !_discr_per_cell->isAllocated() || _discr_per_cell->getNumberOfComponents()!=1)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode");
+  if(code.size()%3!=0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
+  int nbOfSplit=(int)idsPerType.size();
+  int nbOfTypes=(int)code.size()/3;
+  int ret=0;
+  for(int i=0;i<nbOfTypes;i++)
+    {
+      int nbOfEltInChunk=code[3*i+1];
+      if(nbOfEltInChunk<0)
+        throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
+      int pos=code[3*i+2];
+      if(pos!=-1)
+        {
+          if(pos<0 || pos>=nbOfSplit)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+          const DataArrayInt *ids(idsPerType[pos]);
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+        }
+      ret+=nbOfEltInChunk;
+    }
+  if(ret!=_discr_per_cell->getNumberOfTuples())
+    {
+      std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : input code points to " << ret << " cells whereas discretization percell array lgth is " <<  _discr_per_cell->getNumberOfTuples() << " !";
+    }
+  return getNumberOfTuples(0);//0 is not an error ! It is to be sure that input mesh is not used
+}
+
 int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh *) const throw(INTERP_KERNEL::Exception)
 {
   int ret=0;
@@ -1255,8 +1366,17 @@ int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh
     throw INTERP_KERNEL::Exception("Discretization is not initialized!");
   const int *dcPtr=_discr_per_cell->getConstPointer();
   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
+  int maxSz=(int)_loc.size();
   for(const int *w=dcPtr;w!=dcPtr+nbOfTuples;w++)
-    ret+=_loc[*w].getNumberOfGaussPt();
+    {
+      if(*w>=0 && *w<maxSz)
+        ret+=_loc[*w].getNumberOfGaussPt();
+      else
+        {
+          std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuples : At cell #" << std::distance(dcPtr,w) << " localization id is " << *w << " should be in [0," << maxSz << ") !";
+          throw INTERP_KERNEL::Exception(oss.str().c_str());
+        }
+    }
   return ret;
 }
 
@@ -1297,7 +1417,7 @@ DataArrayInt *MEDCouplingFieldDiscretizationGauss::getOffsetArr(const MEDCouplin
   return ret.retn();
 }
 
-void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
+void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
                                                                 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 {
   if(!mesh)
@@ -1321,12 +1441,12 @@ void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplin
         array2[j]=array3[array[i]]+k;
     }
   delete [] array3;
-  for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
+  for(std::vector<DataArray *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
     if(*it)
       (*it)->renumberInPlace(array2);
   delete [] array2;
   if(check)
-    delete [] const_cast<int*>(array);
+    free(const_cast<int*>(array));
 }
 
 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
@@ -1460,10 +1580,10 @@ double MEDCouplingFieldDiscretizationGauss::getIJK(const MEDCouplingMesh *mesh,
   return da->getIJ(offset+nodeIdInCell,compoId);
 }
 
-void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
+void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
 {
-  if(!mesh)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween : NULL input mesh !");
+  if(!mesh || !da)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween : NULL input mesh or DataArray !");
   MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(mesh,da);
   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
     (*iter).checkCoherency();
@@ -1925,6 +2045,44 @@ bool MEDCouplingFieldDiscretizationGaussNE::isEqualIfNotWhy(const MEDCouplingFie
   return ret;
 }
 
+int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
+{
+  if(code.size()%3!=0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
+  int nbOfSplit=(int)idsPerType.size();
+  int nbOfTypes=(int)code.size()/3;
+  int ret=0;
+  for(int i=0;i<nbOfTypes;i++)
+    {
+      int nbOfEltInChunk=code[3*i+1];
+      if(nbOfEltInChunk<0)
+        throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
+      int pos=code[3*i+2];
+      if(pos!=-1)
+        {
+          if(pos<0 || pos>=nbOfSplit)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+          const DataArrayInt *ids(idsPerType[pos]);
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+        }
+      ret+=nbOfEltInChunk;
+    }
+  if(!mesh)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : NULL input mesh !");
+  if(ret!=mesh->getNumberOfCells())
+    {
+      std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : input code points to " << ret << " number of cells should be " <<  mesh->getNumberOfCells() << " !";
+    }
+  return getNumberOfTuples(mesh);
+}
+
 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuples(const MEDCouplingMesh *mesh) const throw(INTERP_KERNEL::Exception)
 {
   if(!mesh)
@@ -1969,7 +2127,7 @@ DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::getOffsetArr(const MEDCoupl
   return ret;
 }
 
-void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
+void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
                                                                   const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 {
   if(!mesh)
@@ -1997,12 +2155,12 @@ void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCoupl
         array2[j]=array3[array[i]]+k;
     }
   delete [] array3;
-  for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
+  for(std::vector<DataArray *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
     if(*it)
       (*it)->renumberInPlace(array2);
   delete [] array2;
   if(check)
-    delete [] const_cast<int *>(array);
+    free(const_cast<int *>(array));
 }
 
 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
@@ -2209,7 +2367,7 @@ double MEDCouplingFieldDiscretizationGaussNE::getIJK(const MEDCouplingMesh *mesh
   return da->getIJ(offset+nodeIdInCell,compoId);
 }
 
-void MEDCouplingFieldDiscretizationGaussNE::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
+void MEDCouplingFieldDiscretizationGaussNE::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
 {
   int nbOfTuples=getNumberOfTuples(mesh);
   if(nbOfTuples!=da->getNumberOfTuples())