getCellsContainingPoints is sensitive to 2D quadratic static cells.

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

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 9da7ee6b4167d9709e6c0dadbfd3df67e6f072d2..92bed51db1a58cfae044871c51b9383aa0e01126 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -16,7 +16,7 @@
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
-// Author : Anthony Geay (CEA/DEN)
+// Author : Anthony Geay (EDF R&D)
 
 #include "MEDCouplingFieldDiscretization.hxx"
 #include "MEDCouplingCMesh.hxx"
@@ -140,7 +140,7 @@ MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField
     case MEDCouplingFieldDiscretizationKriging::TYPE:
       return new MEDCouplingFieldDiscretizationKriging;
     default:
-      throw INTERP_KERNEL::Exception("Choosen discretization is not implemented yet.");
+      throw INTERP_KERNEL::Exception("Chosen discretization is not implemented yet.");
   }
 }
 
@@ -288,8 +288,7 @@ void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const
   if(!arr)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::integral : input array is NULL !");
   MCAuto<MEDCouplingFieldDouble> vol=getMeasureField(mesh,isWAbs);
-  int nbOfCompo=arr->getNumberOfComponents();
-  int nbOfElems=getNumberOfTuples(mesh);
+  std::size_t nbOfCompo(arr->getNumberOfComponents()),nbOfElems(getNumberOfTuples(mesh));
   if(nbOfElems!=arr->getNumberOfTuples())
     {
       std::ostringstream oss; oss << "MEDCouplingFieldDiscretization::integral : field is not correct ! number of tuples in array is " << arr->getNumberOfTuples();
@@ -297,10 +296,9 @@ void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const
       throw INTERP_KERNEL::Exception(oss.str().c_str());
     }
   std::fill(res,res+nbOfCompo,0.);
-  const double *arrPtr=arr->getConstPointer();
-  const double *volPtr=vol->getArray()->getConstPointer();
+  const double *arrPtr(arr->begin()),*volPtr(vol->getArray()->begin());
   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
-  for (int i=0;i<nbOfElems;i++)
+  for(std::size_t i=0;i<nbOfElems;i++)
     {
       std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,(double *)tmp,std::bind2nd(std::multiplies<double>(),volPtr[i]));
       std::transform((double *)tmp,(double *)tmp+nbOfCompo,res,res,std::plus<double>());
@@ -363,7 +361,7 @@ void MEDCouplingFieldDiscretization::finishUnserialization(const std::vector<dou
 
 /*!
  * This method is typically the first step of renumbering. The implementation is empty it is not a bug only gauss is impacted
- * virtualy by this method.
+ * virtually by this method.
  */
 void MEDCouplingFieldDiscretization::renumberCells(const int *old2NewBg, bool check)
 {
@@ -555,7 +553,7 @@ int MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode(con
               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)
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || (int)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());
@@ -799,7 +797,7 @@ int MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCod
               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)
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || (int)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());
@@ -862,7 +860,7 @@ void MEDCouplingFieldDiscretizationOnNodes::checkCoherencyBetween(const MEDCoupl
 {
   if(!mesh || !da)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::checkCoherencyBetween : NULL input mesh or DataArray !");
-  if(mesh->getNumberOfNodes()!=da->getNumberOfTuples())
+  if(mesh->getNumberOfNodes()!=(int)da->getNumberOfTuples())
     {
       std::ostringstream message;
       message << "Field on nodes invalid because there are " << mesh->getNumberOfNodes();
@@ -1143,7 +1141,7 @@ void MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(const MEDCoupl
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has no discretization per cell !");
   if(!mesh)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween : NULL input mesh or DataArray !");
-  int nbOfTuples=_discr_per_cell->getNumberOfTuples();
+  std::size_t 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 !");
 }
@@ -1185,7 +1183,7 @@ bool MEDCouplingFieldDiscretizationPerCell::isEqualWithoutConsideringStr(const M
 
 /*!
  * This method is typically the first step of renumbering. The impact on _discr_per_cell is necessary here.
- * virtualy by this method.
+ * virtually by this method.
  */
 void MEDCouplingFieldDiscretizationPerCell::renumberCells(const int *old2NewBg, bool check)
 {
@@ -1395,7 +1393,7 @@ int MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode(
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
   int nbOfSplit=(int)idsPerType.size();
   int nbOfTypes=(int)code.size()/3;
-  int ret=0;
+  std::size_t ret(0);
   for(int i=0;i<nbOfTypes;i++)
     {
       int nbOfEltInChunk=code[3*i+1];
@@ -1410,7 +1408,7 @@ int MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode(
               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)
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || (int)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());
@@ -1462,16 +1460,16 @@ DataArrayInt *MEDCouplingFieldDiscretizationGauss::getOffsetArr(const MEDCouplin
 {
   if(!mesh)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getOffsetArr : NULL input mesh !");
-  int nbOfTuples=mesh->getNumberOfCells();
+  std::size_t nbOfTuples(mesh->getNumberOfCells());
   MCAuto<DataArrayInt> ret=DataArrayInt::New();
   ret->alloc(nbOfTuples+1,1);
-  int *retPtr=ret->getPointer();
-  const int *start=_discr_per_cell->getConstPointer();
+  int *retPtr(ret->getPointer());
+  const int *start(_discr_per_cell->begin());
   if(_discr_per_cell->getNumberOfTuples()!=nbOfTuples)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getOffsetArr : mismatch between the mesh and the discretization ids array length !");
   int maxPossible=(int)_loc.size();
   retPtr[0]=0;
-  for(int i=0;i<nbOfTuples;i++,start++)
+  for(std::size_t i=0;i<nbOfTuples;i++,start++)
     {
       if(*start>=0 && *start<maxPossible)
         retPtr[i+1]=retPtr[i]+_loc[*start].getNumberOfGaussPt();
@@ -1541,15 +1539,14 @@ DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValue
     {
       INTERP_KERNEL::GaussCoords calculator;
       //
-      const MEDCouplingGaussLocalization& cli=_loc[locIds[i]];//curLocInfo
-      INTERP_KERNEL::NormalizedCellType typ=cli.getType();
-      const std::vector<double>& wg=cli.getWeights();
+      const MEDCouplingGaussLocalization& cli(_loc[locIds[i]]);//curLocInfo
+      INTERP_KERNEL::NormalizedCellType typ(cli.getType());
+      const std::vector<double>& wg(cli.getWeights());
       calculator.addGaussInfo(typ,INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension(),
           &cli.getGaussCoords()[0],(int)wg.size(),&cli.getRefCoords()[0],
           INTERP_KERNEL::CellModel::GetCellModel(typ).getNumberOfNodes());
       //
-      int nbt=parts2[i]->getNumberOfTuples();
-      for(const int *w=parts2[i]->getConstPointer();w!=parts2[i]->getConstPointer()+nbt;w++)
+      for(const int *w=parts2[i]->begin();w!=parts2[i]->end();w++)
         calculator.calculateCoords(cli.getType(),coords,spaceDim,conn+connI[*w]+1,valsToFill+spaceDim*(ptrOffsets[*w]));
     }
   ret->copyStringInfoFrom(*umesh->getCoords());
@@ -1666,7 +1663,7 @@ void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplin
     {
       if(dc[i]>=nbOfDesc)
         {
-          std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has an undefined gauss location ! Should never happend !";
+          std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has an undefined gauss location ! Should never happened !";
           throw INTERP_KERNEL::Exception(oss.str().c_str());
         }
       if(dc[i]<0)
@@ -1680,10 +1677,10 @@ void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplin
           throw INTERP_KERNEL::Exception(oss.str().c_str());
         }
     }
-  int nbOfTuples=getNumberOfTuples(mesh);
+  std::size_t nbOfTuples(getNumberOfTuples(mesh));
   if(nbOfTuples!=da->getNumberOfTuples())
     {
-      std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
+      std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " having " << da->getNumberOfTuples() << " !";
       throw INTERP_KERNEL::Exception(oss.str().c_str());
     }
 }
@@ -1745,7 +1742,7 @@ void MEDCouplingFieldDiscretizationGauss::getValueOn(const DataArrayDouble *arr,
 
 void MEDCouplingFieldDiscretizationGauss::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
 {
-  throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
+  throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applicable for Gauss points !");
 }
 
 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
@@ -1821,7 +1818,7 @@ DataArrayInt *MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCe
   if(!mesh)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null mesh !");
   MCAuto<DataArrayInt> nbOfNodesPerCell=buildNbOfGaussPointPerCellField();//check of _discr_per_cell not NULL pointer
-  int nbOfCells=mesh->getNumberOfCells();
+  std::size_t nbOfCells(mesh->getNumberOfCells());
   if(_discr_per_cell->getNumberOfTuples()!=nbOfCells)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : mismatch of nb of tuples of cell ids array and number of cells !");
   nbOfNodesPerCell->computeOffsetsFull();
@@ -2166,7 +2163,7 @@ int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCod
               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)
+          if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || (int)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());
@@ -2587,7 +2584,7 @@ double MEDCouplingFieldDiscretizationGaussNE::getIJK(const MEDCouplingMesh *mesh
 
 void MEDCouplingFieldDiscretizationGaussNE::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const
 {
-  int nbOfTuples=getNumberOfTuples(mesh);
+  std::size_t nbOfTuples(getNumberOfTuples(mesh));
   if(nbOfTuples!=da->getNumberOfTuples())
     {
       std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
@@ -2637,7 +2634,7 @@ void MEDCouplingFieldDiscretizationGaussNE::getValueOn(const DataArrayDouble *ar
 
 void MEDCouplingFieldDiscretizationGaussNE::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
 {
-  throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
+  throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applicable for Gauss points !");
 }
 
 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
@@ -2797,7 +2794,7 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::getValueOnMulti(const Da
 {
   if(!arr || !arr->isAllocated())
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::getValueOnMulti : input array is null or not allocated !");
-  int nbOfRows(getNumberOfMeshPlaces(mesh));
+  std::size_t nbOfRows(getNumberOfMeshPlaces(mesh));
   if(arr->getNumberOfTuples()!=nbOfRows)
     {
       std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationKriging::getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !";
@@ -3027,9 +3024,9 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec(const
  */
 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const
 {
-  int spaceDimension=arr->getNumberOfComponents();
+  std::size_t spaceDimension(arr->getNumberOfComponents());
   delta=spaceDimension+1;
-  int szOfMatrix=arr->getNumberOfTuples();
+  std::size_t szOfMatrix(arr->getNumberOfTuples());
   if(szOfMatrix*szOfMatrix!=matr->getNumberOfTuples())
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::performDrift : invalid size");
   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
@@ -3037,7 +3034,7 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataA
   const double *srcWork=matr->getConstPointer();
   const double *srcWork2=arr->getConstPointer();
   double *destWork=ret->getPointer();
-  for(int i=0;i<szOfMatrix;i++)
+  for(std::size_t i=0;i<szOfMatrix;i++)
     {
       destWork=std::copy(srcWork,srcWork+szOfMatrix,destWork);
       srcWork+=szOfMatrix;
@@ -3049,7 +3046,7 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataA
   std::fill(destWork,destWork+spaceDimension+1,0.); destWork+=spaceDimension+1;
   MCAuto<DataArrayDouble> arrNoI=arr->toNoInterlace();
   srcWork2=arrNoI->getConstPointer();
-  for(int i=0;i<spaceDimension;i++)
+  for(std::size_t i=0;i<spaceDimension;i++)
     {
       destWork=std::copy(srcWork2,srcWork2+szOfMatrix,destWork);
       srcWork2+=szOfMatrix;