Salome HOME
getCellsContainingPoints is sensitive to 2D quadratic static cells.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingFieldDouble.cxx
index 504451a71f6e3cae108062e06b116d97d4459686..6e6d5dc21c5c55c1b0e7802e9846071fa86940a4 100644 (file)
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
-// Author : Anthony Geay (CEA/DEN)
+// Author : Anthony Geay (EDF R&D)
 
 #include "MEDCouplingFieldDouble.hxx"
 #include "MEDCouplingFieldTemplate.hxx"
 #include "MEDCouplingFieldT.txx"
 #include "MEDCouplingFieldInt.hxx"
+#include "MEDCouplingFieldFloat.hxx"
 #include "MEDCouplingUMesh.hxx"
 #include "MEDCouplingTimeDiscretization.hxx"
 #include "MEDCouplingFieldDiscretization.hxx"
 #include "MCAuto.txx"
 #include "MEDCouplingVoronoi.hxx"
 #include "MEDCouplingNatureOfField.hxx"
+#include "MEDCouplingMemArray.txx"
 
 #include "InterpKernelAutoPtr.hxx"
 #include "InterpKernelGaussCoords.hxx"
@@ -174,9 +176,9 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::deepCopy() const
  * \endif
  * \sa clone()
  */
-MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildNewTimeReprFromThis(TypeOfTimeDiscretization td, bool deepCopy) const
+MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildNewTimeReprFromThis(TypeOfTimeDiscretization td, bool deepCpy) const
 {
-  MEDCouplingTimeDiscretization *tdo=timeDiscr()->buildNewTimeReprFromThis(td,deepCopy);
+  MEDCouplingTimeDiscretization *tdo=timeDiscr()->buildNewTimeReprFromThis(td,deepCpy);
   MCAuto<MEDCouplingFieldDiscretization> disc;
   if(_type)
     disc=_type->clone();
@@ -188,8 +190,8 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildNewTimeReprFromThis(TypeOfT
 }
 
 /*!
- * This method converts a field on nodes (\a this) to a cell field (returned field). The convertion is a \b non \b conservative remapping !
- * This method is useful only for users that need a fast convertion from node to cell spatial discretization. The algorithm applied is simply to attach
+ * This method converts a field on nodes (\a this) to a cell field (returned field). The conversion is a \b non \b conservative remapping !
+ * This method is useful only for users that need a fast conversion from node to cell spatial discretization. The algorithm applied is simply to attach
  * to each cell the average of values on nodes constituting this cell.
  *
  * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
@@ -244,8 +246,8 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::nodeToCellDiscretization() const
 }
 
 /*!
- * This method converts a field on cell (\a this) to a node field (returned field). The convertion is a \b non \b conservative remapping !
- * This method is useful only for users that need a fast convertion from cell to node spatial discretization. The algorithm applied is simply to attach
+ * This method converts a field on cell (\a this) to a node field (returned field). The conversion is a \b non \b conservative remapping !
+ * This method is useful only for users that need a fast conversion from cell to node spatial discretization. The algorithm applied is simply to attach
  * to each node the average of values on cell sharing this node. If \a this lies on a mesh having orphan nodes the values applied on them will be NaN (division by 0.).
  *
  * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
@@ -367,76 +369,6 @@ bool MEDCouplingFieldDouble::areCompatibleForMeld(const MEDCouplingFieldDouble *
   return true;
 }
 
-/*!
- * Permutes values of \a this field according to a given permutation array for cells
- * renumbering. The underlying mesh is deeply copied and its cells are also permuted. 
- * The number of cells remains the same; for that the permutation array \a old2NewBg
- * should not contain equal ids.
- * ** Warning, this method modifies the mesh aggreagated by \a this (by performing a deep copy ) **.
- *
- *  \param [in] old2NewBg - the permutation array in "Old to New" mode. Its length is
- *         to be equal to \a this->getMesh()->getNumberOfCells().
- *  \param [in] check - if \c true, \a old2NewBg is transformed to a new permutation
- *         array, so that its maximal cell id to correspond to (be less than) the number
- *         of cells in mesh. This new array is then used for the renumbering. If \a 
- *         check == \c false, \a old2NewBg is used as is, that is less secure as validity 
- *         of ids in \a old2NewBg is not checked.
- *  \throw If the mesh is not set.
- *  \throw If the spatial discretization of \a this field is NULL.
- *  \throw If \a check == \c true and \a old2NewBg contains equal ids.
- *  \throw If mesh nature does not allow renumbering (e.g. structured mesh).
- * 
- *  \if ENABLE_EXAMPLES
- *  \ref cpp_mcfielddouble_renumberCells "Here is a C++ example".<br>
- *  \ref  py_mcfielddouble_renumberCells "Here is a Python example".
- *  \endif
- */
-void MEDCouplingFieldDouble::renumberCells(const int *old2NewBg, bool check)
-{
-  renumberCellsWithoutMesh(old2NewBg,check);
-  MCAuto<MEDCouplingMesh> m=_mesh->deepCopy();
-  m->renumberCells(old2NewBg,check);
-  setMesh(m);
-  updateTime();
-}
-
-/*!
- * Permutes values of \a this field according to a given permutation array for cells
- * renumbering. The underlying mesh is \b not permuted. 
- * The number of cells remains the same; for that the permutation array \a old2NewBg
- * should not contain equal ids.
- * This method performs a part of job of renumberCells(). The reasonable use of this
- * method is only for multi-field instances lying on the same mesh to avoid a
- * systematic duplication and renumbering of _mesh attribute. 
- * \warning Use this method with a lot of care!
- *  \param [in] old2NewBg - the permutation array in "Old to New" mode. Its length is
- *         to be equal to \a this->getMesh()->getNumberOfCells().
- *  \param [in] check - if \c true, \a old2NewBg is transformed to a new permutation
- *         array, so that its maximal cell id to correspond to (be less than) the number
- *         of cells in mesh. This new array is then used for the renumbering. If \a 
- *         check == \c false, \a old2NewBg is used as is, that is less secure as validity 
- *         of ids in \a old2NewBg is not checked.
- *  \throw If the mesh is not set.
- *  \throw If the spatial discretization of \a this field is NULL.
- *  \throw If \a check == \c true and \a old2NewBg contains equal ids.
- *  \throw If mesh nature does not allow renumbering (e.g. structured mesh).
- */
-void MEDCouplingFieldDouble::renumberCellsWithoutMesh(const int *old2NewBg, bool check)
-{
-  if(!_mesh)
-    throw INTERP_KERNEL::Exception("Expecting a defined mesh to be able to operate a renumbering !");
-  if(_type.isNull())
-    throw INTERP_KERNEL::Exception("Expecting a spatial discretization to be able to operate a renumbering !");
-  //
-  _type->renumberCells(old2NewBg,check);
-  std::vector<DataArrayDouble *> arrays;
-  timeDiscr()->getArrays(arrays);
-  std::vector<DataArray *> arrays2(arrays.size()); std::copy(arrays.begin(),arrays.end(),arrays2.begin());
-  _type->renumberArraysForCell(_mesh,arrays2,old2NewBg,check);
-  //
-  updateTime();
-}
-
 /*!
  * Permutes values of \a this field according to a given permutation array for node
  * renumbering. The underlying mesh is deeply copied and its nodes are also permuted. 
@@ -449,7 +381,7 @@ void MEDCouplingFieldDouble::renumberCellsWithoutMesh(const int *old2NewBg, bool
  *  \throw If the spatial discretization of \a this field is NULL.
  *  \throw If \a check == \c true and \a old2NewBg contains equal ids.
  *  \throw If mesh nature does not allow renumbering (e.g. structured mesh).
- *  \throw If values at merged nodes deffer more than \a eps.
+ *  \throw If values at merged nodes differ more than \a eps.
  * 
  *  \if ENABLE_EXAMPLES
  *  \ref cpp_mcfielddouble_renumberNodes "Here is a C++ example".<br>
@@ -489,7 +421,7 @@ void MEDCouplingFieldDouble::renumberNodes(const int *old2NewBg, double eps)
  *         the values differ more than \a eps, an exception is thrown.
  *  \throw If the mesh is not set.
  *  \throw If the spatial discretization of \a this field is NULL.
- *  \throw If values at merged nodes deffer more than \a eps.
+ *  \throw If values at merged nodes differ more than \a eps.
  */
 void MEDCouplingFieldDouble::renumberNodesWithoutMesh(const int *old2NewBg, int newNbOfNodes, double eps)
 {
@@ -522,21 +454,32 @@ DataArrayInt *MEDCouplingFieldDouble::findIdsInRange(double vmin, double vmax) c
   return getArray()->findIdsInRange(vmin,vmax);
 }
 
-MEDCouplingFieldInt *MEDCouplingFieldDouble::convertToIntField() const
+template<class U>
+typename Traits<U>::FieldType *ConvertToUField(const MEDCouplingFieldDouble *self)
 {
-  MCAuto<MEDCouplingFieldTemplate> tmp(MEDCouplingFieldTemplate::New(*this));
+  MCAuto<MEDCouplingFieldTemplate> tmp(MEDCouplingFieldTemplate::New(*self));
   int t1,t2;
-  double t0(getTime(t1,t2));
-  MCAuto<MEDCouplingFieldInt> ret(MEDCouplingFieldInt::New(*tmp,getTimeDiscretization()));
+  double t0(self->getTime(t1,t2));
+  MCAuto<typename Traits<U>::FieldType > ret(Traits<U>::FieldType::New(*tmp,self->getTimeDiscretization()));
   ret->setTime(t0,t1,t2);
-  if(getArray())
+  if(self->getArray())
     {
-      MCAuto<DataArrayInt> arr(getArray()->convertToIntArr());
+      MCAuto<typename Traits<U>::ArrayType> arr(self->getArray()->convertToOtherTypeOfArr<U>());
       ret->setArray(arr);
     }
   return ret.retn();
 }
 
+MEDCouplingFieldInt *MEDCouplingFieldDouble::convertToIntField() const
+{
+  return ConvertToUField<int>(this);
+}
+
+MEDCouplingFieldFloat *MEDCouplingFieldDouble::convertToFloatField() const
+{
+  return ConvertToUField<float>(this);
+}
+
 MEDCouplingFieldDouble::MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td):MEDCouplingFieldT<double>(type,MEDCouplingTimeDiscretization::New(td))
 {
 }
@@ -548,7 +491,7 @@ MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldTemplate& f
 {
 }
 
-MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCopy):MEDCouplingFieldT<double>(other,deepCopy)
+MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCpy):MEDCouplingFieldT<double>(other,deepCpy)
 {
 }
 
@@ -616,7 +559,7 @@ double MEDCouplingFieldDouble::getMaxValue() const
 /*!
  * Returns the maximal value and all its locations within \a this scalar field.
  * Only the first of available data arrays is checked.
- *  \param [out] tupleIds - a new instance of DataArrayInt containg indices of
+ *  \param [out] tupleIds - a new instance of DataArrayInt containing indices of
  *               tuples holding the maximal value. The caller is to delete it using
  *               decrRef() as it is no more needed.
  *  \return double - the maximal value among all values of the first array of \a this filed.
@@ -680,7 +623,7 @@ double MEDCouplingFieldDouble::getMinValue() const
 /*!
  * Returns the minimal value and all its locations within \a this scalar field.
  * Only the first of available data arrays is checked.
- *  \param [out] tupleIds - a new instance of DataArrayInt containg indices of
+ *  \param [out] tupleIds - a new instance of DataArrayInt containing indices of
  *               tuples holding the minimal value. The caller is to delete it using
  *               decrRef() as it is no more needed.
  *  \return double - the minimal value among all values of the first array of \a this filed.
@@ -740,19 +683,6 @@ double MEDCouplingFieldDouble::norm2() const
   return getArray()->norm2();
 }
 
-/*!
- * This method returns the max norm of \a this field.
- * \f[
- * \max_{0 \leq i < nbOfEntity}{abs(val[i])}
- * \f]
- *  \throw If the data array is not set.
- */
-double MEDCouplingFieldDouble::normMax() const
-{
-  if(getArray()==0)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::normMax : no default array defined !");
-  return getArray()->normMax();
-}
 
 /*!
  * Computes the weighted average of values of each component of \a this field, the weights being the
@@ -899,6 +829,48 @@ void MEDCouplingFieldDouble::normL2(double *res) const
   _type->normL2(_mesh,getArray(),res);
 }
 
+/*!
+ * Returns the \c infinite norm of values of a given component of \a this field:
+* \f[
+ * \max_{0 \leq i < nbOfEntity}{abs(val[i])}
+ * \f]
+ *  \param [in] compId - an index of the component of interest.
+ *  \throw If \a compId is not valid.
+           A valid range is ( 0 <= \a compId < \a this->getNumberOfComponents() ).
+ *  \throw If the data array is not set.
+ */
+double MEDCouplingFieldDouble::normMax(int compId) const
+{
+  if(getArray()==0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::normMax : no default array defined !");
+  int nbComps=getArray()->getNumberOfComponents();
+  if(compId<0 || compId>=nbComps)
+    {
+      std::ostringstream oss; oss << "MEDCouplingFieldDouble::normMax : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
+      throw INTERP_KERNEL::Exception(oss.str());
+    }
+  INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
+  getArray()->normMaxPerComponent(res);
+  return res[compId];
+}
+
+/*!
+ * Returns the \c infinite norm of values of each component of \a this field:
+ * \f[
+ * \max_{0 \leq i < nbOfEntity}{abs(val[i])}
+ * \f]
+ *  \param [out] res - pointer to an array of result values, of size at least \a
+ *         this->getNumberOfComponents(), that is to be allocated by the caller.
+ *  \throw If the data array is not set.
+ *
+ */
+void MEDCouplingFieldDouble::normMax(double *res) const
+{
+  if(getArray()==0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::normMax : no default array defined !");
+  getArray()->normMaxPerComponent(res);
+}
+
 /*!
  * Computes a sum of values of a given component of \a this field multiplied by
  * values returned by buildMeasureField().
@@ -1340,7 +1312,7 @@ void MEDCouplingFieldDouble::applyFunc(int nbOfComp, double val)
  *   - "2*x + z"               produces (5.,5.,5.,5.)
  *   - "2*x + 0*y + z"         produces (9.,9.,9.,9.)
  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,4.)
- *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,4.)
+ *   - "2*y*IVec + z*KVec + x" produces (7.,1.,8.,1.)
  *
  *  \param [in] nbOfComp - the number of components for \a this field to have.
  *  \param [in] func - the function used to compute values of \a this field.
@@ -1380,7 +1352,7 @@ void MEDCouplingFieldDouble::applyFunc(int nbOfComp, const std::string& func)
  *   - "2*x + z"               produces (5.,5.,5.,5.)
  *   - "2*x + 0*y + z"         produces (9.,9.,9.,9.)
  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,4.)
- *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,4.)
+ *   - "2*y*IVec + z*KVec + x" produces (7.,1.,8.,1.)
  *
  *  \param [in] nbOfComp - the number of components for \a this field to have.
  *  \param [in] func - the function used to compute values of \a this field.
@@ -1644,7 +1616,7 @@ double MEDCouplingFieldDouble::getIJK(int cellId, int nodeIdInCell, int compoId)
  *  \throw If \a other == NULL.
  *  \throw If any of the meshes is not well defined.
  *  \throw If the two meshes do not match.
- *  \throw If field values at merged nodes (if any) deffer more than \a eps.
+ *  \throw If field values at merged nodes (if any) differ more than \a eps.
  *
  *  \if ENABLE_EXAMPLES
  *  \ref cpp_mcfielddouble_changeUnderlyingMesh "Here is a C++ example".<br>
@@ -1692,7 +1664,7 @@ void MEDCouplingFieldDouble::changeUnderlyingMesh(const MEDCouplingMesh *other,
  *  \throw If any of the meshes is not set or is not well defined.
  *  \throw If the two meshes do not match.
  *  \throw If the two fields are not coherent for merge.
- *  \throw If field values at merged nodes (if any) deffer more than \a eps.
+ *  \throw If field values at merged nodes (if any) differ more than \a eps.
  *
  *  \if ENABLE_EXAMPLES
  *  \ref cpp_mcfielddouble_substractInPlaceDM "Here is a C++ example".<br>
@@ -1725,7 +1697,7 @@ void MEDCouplingFieldDouble::substractInPlaceDM(const MEDCouplingFieldDouble *f,
  *  \throw If the mesh is not well defined.
  *  \throw If the spatial discretization of \a this field is NULL.
  *  \throw If the data array is not set.
- *  \throw If field values at merged nodes (if any) deffer more than \a epsOnVals.
+ *  \throw If field values at merged nodes (if any) differ more than \a epsOnVals.
  */
 bool MEDCouplingFieldDouble::mergeNodes(double eps, double epsOnVals)
 {
@@ -1764,7 +1736,7 @@ bool MEDCouplingFieldDouble::mergeNodes(double eps, double epsOnVals)
  *  \throw If the mesh is not well defined.
  *  \throw If the spatial discretization of \a this field is NULL.
  *  \throw If the data array is not set.
- *  \throw If field values at merged nodes (if any) deffer more than \a epsOnVals.
+ *  \throw If field values at merged nodes (if any) differ more than \a epsOnVals.
  */
 bool MEDCouplingFieldDouble::mergeNodesCenter(double eps, double epsOnVals)
 {
@@ -1801,7 +1773,7 @@ bool MEDCouplingFieldDouble::mergeNodesCenter(double eps, double epsOnVals)
  *  \throw If the mesh is not well defined.
  *  \throw If the spatial discretization of \a this field is NULL.
  *  \throw If the data array is not set.
- *  \throw If field values at merged nodes (if any) deffer more than \a epsOnVals.
+ *  \throw If field values at merged nodes (if any) differ more than \a epsOnVals.
  */
 bool MEDCouplingFieldDouble::zipCoords(double epsOnVals)
 {
@@ -1842,7 +1814,7 @@ bool MEDCouplingFieldDouble::zipCoords(double epsOnVals)
  *  \throw If the mesh is not well defined.
  *  \throw If the spatial discretization of \a this field is NULL.
  *  \throw If the data array is not set.
- *  \throw If field values at merged cells (if any) deffer more than \a epsOnVals.
+ *  \throw If field values at merged cells (if any) differ more than \a epsOnVals.
  */
 bool MEDCouplingFieldDouble::zipConnectivity(int compType, double epsOnVals)
 {
@@ -1852,7 +1824,7 @@ bool MEDCouplingFieldDouble::zipConnectivity(int compType, double epsOnVals)
   if(_type.isNull())
     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform zipConnectivity !");
   MCAuto<MEDCouplingUMesh> meshC2((MEDCouplingUMesh *)meshC->deepCopy());
-  int oldNbOfCells=meshC2->getNumberOfCells();
+  std::size_t oldNbOfCells(meshC2->getNumberOfCells());
   MCAuto<DataArrayInt> arr=meshC2->zipConnectivityTraducer(compType);
   if(meshC2->getNumberOfCells()!=oldNbOfCells)
     {
@@ -2016,7 +1988,6 @@ MCAuto<MEDCouplingFieldDouble> MEDCouplingFieldDouble::convertQuadraticCellsToLi
         if(!disc2)
           throw INTERP_KERNEL::Exception("convertQuadraticCellsToLinear : Not a ON_GAUSS_PT field");
         std::set<INTERP_KERNEL::NormalizedCellType> gt2(umesh->getAllGeoTypes());
-        const DataArrayDouble *arr(getArray());
         std::vector< MCAuto<DataArrayInt> > cellIdsV;
         std::vector< MCAuto<MEDCouplingUMesh> > meshesV;
         std::vector< MEDCouplingGaussLocalization > glV;
@@ -2071,7 +2042,7 @@ MCAuto<MEDCouplingFieldDouble> MEDCouplingFieldDouble::convertQuadraticCellsToLi
 
 /*!
  * This is expected to be a 3 components vector field on nodes (if not an exception will be thrown). \a this is also expected to lie on a MEDCouplingPointSet mesh.
- * Finaly \a this is also expected to be consistent.
+ * Finally \a this is also expected to be consistent.
  * In these conditions this method returns a newly created field (to be dealed by the caller).
  * The returned field will also 3 compo vector field be on nodes lying on the same mesh than \a this.
  * 
@@ -2980,12 +2951,12 @@ MCAuto<MEDCouplingFieldDouble> MEDCouplingFieldDouble::voronoizeGen(const Voroni
         ptsInReal=gl.localizePtsInRefCooForEachCell(vorCellsForCurDisc->getCoords(),subMesh);
       }
       int nbPtsPerCell(vorCellsForCurDisc->getNumberOfNodes());
-      for(std::size_t i=0;i<ids.size();i++)
+      for(std::size_t j=0;j<ids.size();j++)
         {
           MCAuto<MEDCouplingUMesh> elt(vorCellsForCurDisc->clone(false));
-          MCAuto<DataArrayDouble> coo(ptsInReal->selectByTupleIdSafeSlice(i*nbPtsPerCell,(i+1)*nbPtsPerCell,1));
-          elt->setCoords(coo);
-          cells[ids[i]]=elt;
+          MCAuto<DataArrayDouble> coo4(ptsInReal->selectByTupleIdSafeSlice(j*nbPtsPerCell,(j+1)*nbPtsPerCell,1));
+          elt->setCoords(coo4);
+          cells[ids[j]]=elt;
         }
     }
   std::vector< const MEDCouplingUMesh * > cellsPtr(VecAutoToVecOfCstPt(cells));