From: ageay <ageay>
Date: Fri, 5 Mar 2010 15:40:13 +0000 (+0000)
Subject: getElementContainingPoint-getValueOn P0
X-Git-Tag: V5_1_main_FINAL~171
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=dde8f52747f20c9fb67a781eb79c31f6c3678e12;p=tools%2Fmedcoupling.git

getElementContainingPoint-getValueOn P0
---

diff --git a/src/MEDCoupling/MEDCouplingCMesh.cxx b/src/MEDCoupling/MEDCouplingCMesh.cxx
index d1b6aec3c..4f4b5150f 100644
--- a/src/MEDCoupling/MEDCouplingCMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingCMesh.cxx
@@ -200,6 +200,12 @@ MEDCouplingFieldDouble *MEDCouplingCMesh::buildOrthogonalField() const
   return ret;
 }
 
+int MEDCouplingCMesh::getElementContainingPoint(const double *pos, double eps) const
+{
+  //not implemented yet !
+  return -1;
+}
+
 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
 {
   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
diff --git a/src/MEDCoupling/MEDCouplingCMesh.hxx b/src/MEDCoupling/MEDCouplingCMesh.hxx
index 4f4096a44..201a828d5 100644
--- a/src/MEDCoupling/MEDCouplingCMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingCMesh.hxx
@@ -48,6 +48,7 @@ namespace ParaMEDMEM
     MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
     MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
     MEDCouplingFieldDouble *buildOrthogonalField() const;
+    int getElementContainingPoint(const double *pos, double eps) const;
     void rotate(const double *center, const double *vector, double angle);
     void translate(const double *vector);
     MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
index a286dc69b..39e440fc2 100644
--- a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
@@ -132,6 +132,12 @@ MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::buildOrthogonalField() const
   throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::buildOrthogonalField not implemented yet !");
 }
 
+int MEDCouplingExtrudedMesh::getElementContainingPoint(const double *pos, double eps) const
+{
+  //not implemented yet
+  return -1;
+}
+
 MEDCouplingExtrudedMesh::~MEDCouplingExtrudedMesh()
 {
   if(_mesh2D)
diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
index 8480bcf2d..ee3abeeef 100644
--- a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
@@ -49,6 +49,7 @@ namespace ParaMEDMEM
     MEDCouplingFieldDouble *getMeasureField(bool) const;
     MEDCouplingFieldDouble *getMeasureFieldOnNode(bool) const;
     MEDCouplingFieldDouble *buildOrthogonalField() const;
+    int getElementContainingPoint(const double *pos, double eps) const;
     static int findCorrespCellByNodalConn(const std::vector<int>& nodalConnec,
                                           const int *revNodalPtr, const int *revNodalIndxPtr) throw(INTERP_KERNEL::Exception);
     void rotate(const double *center, const double *vector, double angle);
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index f3756291c..266ee4bdd 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -26,6 +26,8 @@
 
 using namespace ParaMEDMEM;
 
+const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12;
+
 const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
 
 const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
@@ -34,6 +36,10 @@ const char MEDCouplingFieldDiscretizationP1::REPR[]="P1";
 
 const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
 
+MEDCouplingFieldDiscretization::MEDCouplingFieldDiscretization():_precision(DFLT_PRECISION)
+{
+}
+
 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type)
 {
   switch(type)
@@ -108,6 +114,12 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(cons
   return mesh->getMeasureField(isAbs);
 }
 
+void MEDCouplingFieldDiscretizationP0::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
+{
+  int id=mesh->getElementContainingPoint(loc,_precision);
+  arr->getTuple(id,res);
+}
+
 /*!
  * Nothing to do. It's not a bug.
  */
@@ -183,6 +195,10 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(cons
   return mesh->getMeasureFieldOnNode(isAbs);
 }
 
+void MEDCouplingFieldDiscretizationP1::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
+{
+}
+
 void MEDCouplingFieldDiscretizationP1::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const
 {
   int oldNbOfElems=old2New->getNbOfElems();
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
index fdb2fd115..0d9b90c26 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
@@ -35,6 +35,8 @@ namespace ParaMEDMEM
   {
   public:
     static MEDCouplingFieldDiscretization *New(TypeOfField type);
+    double getPrecision() const { return _precision; }
+    void setPrecision(double val) { _precision=val; }
     static TypeOfField getTypeOfFieldFromStringRepr(const char *repr) throw(INTERP_KERNEL::Exception);
     virtual TypeOfField getEnum() const = 0;
     virtual bool isEqual(const MEDCouplingFieldDiscretization *other) const = 0;
@@ -45,8 +47,14 @@ namespace ParaMEDMEM
     virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0;
     virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0;
     virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const = 0;
+    virtual void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const = 0;
     virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0;
     virtual void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const = 0;
+  protected:
+    MEDCouplingFieldDiscretization();
+  protected:
+    double _precision;
+    static const double DFLT_PRECISION;
   };
 
   class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP0 : public MEDCouplingFieldDiscretization
@@ -61,6 +69,7 @@ namespace ParaMEDMEM
     void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
     void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
     MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const;
+    void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const;
     void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const;
     MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
   public:
@@ -80,6 +89,7 @@ namespace ParaMEDMEM
     void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
     void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
     MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const;
+    void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const;
     MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
     void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const;
     static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di);
diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
index dc2d27652..45d04a2e0 100644
--- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
@@ -165,12 +165,21 @@ void MEDCouplingFieldDouble::measureAccumulate(bool isWAbs, double *res) const
 
 void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception)
 {
-  _time_discr->checkNoTimePresence();
+  const DataArrayDouble *arr=_time_discr->getArray();
+  _type->getValueOn(arr,_mesh,spaceLoc,res);
 }
 
 void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception)
 {
-  _time_discr->checkTimePresence(time);
+  std::vector< const DataArrayDouble *> arrs=_time_discr->getArraysForTime(time);
+  std::vector<double> res2;
+  for(std::vector< const DataArrayDouble *>::const_iterator iter=arrs.begin();iter!=arrs.end();iter++)
+    {
+      int sz=res2.size();
+      res2.resize(sz+(*iter)->getNumberOfComponents());
+      _type->getValueOn(*iter,_mesh,spaceLoc,&res2[sz]);
+    }
+  _time_discr->getValueForTime(time,res2,res);
 }
 
 void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId)
diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx
index 3bbdec6e2..bd5491fb0 100644
--- a/src/MEDCoupling/MEDCouplingMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingMesh.hxx
@@ -56,6 +56,7 @@ namespace ParaMEDMEM
     virtual void getBoundingBox(double *bbox) const = 0;
     virtual MEDCouplingFieldDouble *getMeasureField(bool isAbs) const = 0;
     virtual MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const = 0;
+    virtual int getElementContainingPoint(const double *pos, double eps) const = 0;
     virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const;
     virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const;
     virtual MEDCouplingFieldDouble *buildOrthogonalField() const = 0;
diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx
index 0cac14c02..3a4a4738b 100644
--- a/src/MEDCoupling/MEDCouplingPointSet.cxx
+++ b/src/MEDCoupling/MEDCouplingPointSet.cxx
@@ -401,6 +401,13 @@ bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double*
  * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
  */
 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
+{
+  double *coords=_coords->getPointer();
+  int nbNodes=getNumberOfNodes();
+  rotate3DAlg(center,vect,angle,nbNodes,coords);
+}
+
+void MEDCouplingPointSet::rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
 {
   double sina=sin(angle);
   double cosa=cos(angle);
@@ -422,8 +429,6 @@ void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, dou
   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
   //rotation matrix computed.
-  double *coords=_coords->getPointer();
-  int nbNodes=getNumberOfNodes();
   double tmp[3];
   for(int i=0; i<nbNodes; i++)
     {
@@ -438,13 +443,18 @@ void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, dou
  * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
  */
 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
+{
+  double *coords=_coords->getPointer();
+  int nbNodes=getNumberOfNodes();
+  rotate2DAlg(center,angle,nbNodes,coords);
+}
+
+void MEDCouplingPointSet::rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
 {
   double cosa=cos(angle);
   double sina=sin(angle);
   double matrix[4];
   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
-  double *coords=_coords->getPointer();
-  int nbNodes=getNumberOfNodes();
   double tmp[2];
   for(int i=0; i<nbNodes; i++)
     {
diff --git a/src/MEDCoupling/MEDCouplingPointSet.hxx b/src/MEDCoupling/MEDCouplingPointSet.hxx
index 544580419..1236cebb1 100644
--- a/src/MEDCoupling/MEDCouplingPointSet.hxx
+++ b/src/MEDCoupling/MEDCouplingPointSet.hxx
@@ -55,6 +55,8 @@ namespace ParaMEDMEM
     void tryToShareSameCoords(MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception);
     static DataArrayDouble *mergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2);
     static MEDCouplingPointSet *buildInstanceFromMeshType(MEDCouplingMeshType type);
+    static void rotate2DAlg(const double *center, double angle, int nbNodes, double *coords);
+    static void rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords);
     virtual MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const = 0;
     virtual MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const = 0;
     virtual void findBoundaryNodes(std::vector<int>& nodes) const = 0;
diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx
index fa900e7a4..0f3337323 100644
--- a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx
@@ -395,7 +395,12 @@ void MEDCouplingNoTimeLabel::checkTimePresence(double time) const throw(INTERP_K
   throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
 }
 
-DataArrayDouble *MEDCouplingNoTimeLabel::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+std::vector< const DataArrayDouble *> MEDCouplingNoTimeLabel::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception)
+{
+  throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingNoTimeLabel::getValueForTime(double time, const std::vector<double>& vals, double *res) const
 {
   throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
 }
@@ -590,18 +595,23 @@ void MEDCouplingWithTimeStep::checkTimePresence(double time) const throw(INTERP_
     }
 }
 
-DataArrayDouble *MEDCouplingWithTimeStep::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+std::vector< const DataArrayDouble *> MEDCouplingWithTimeStep::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception)
 {
   if(std::fabs(time-_time)<=_time_tolerance)
     {
-      if(_array)
-        _array->incrRef();
-      return _array;
+      std::vector< const DataArrayDouble *> ret(1);
+      ret[0]=_array;
+      return ret;
     }
   else
     throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
 }
 
+void MEDCouplingWithTimeStep::getValueForTime(double time, const std::vector<double>& vals, double *res) const
+{
+  std::copy(vals.begin(),vals.end(),res);
+}
+
 void MEDCouplingWithTimeStep::getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception)
 {
   if(std::fabs(time-_time)<=_time_tolerance)
@@ -666,18 +676,23 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::performCpy(bool d
   return new MEDCouplingConstOnTimeInterval(*this,deepCpy);
 }
 
-DataArrayDouble *MEDCouplingConstOnTimeInterval::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+std::vector< const DataArrayDouble *> MEDCouplingConstOnTimeInterval::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception)
 {
   if(time>_start_time-_time_tolerance && time<_end_time+_time_tolerance)
     {
-      if(_array)
-        _array->incrRef();
-      return _array;
+      std::vector< const DataArrayDouble *> ret(1);
+      ret[0]=_array;
+      return ret;
     }
   else
     throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
 }
 
+void MEDCouplingConstOnTimeInterval::getValueForTime(double time, const std::vector<double>& vals, double *res) const
+{
+  std::copy(vals.begin(),vals.end(),res);
+}
+
 bool MEDCouplingConstOnTimeInterval::areCompatible(const MEDCouplingTimeDiscretization *other) const
 {
   if(!MEDCouplingTimeDiscretization::areCompatible(other))
diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx
index 158c6bb51..07aea3bfb 100644
--- a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx
+++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx
@@ -61,8 +61,8 @@ namespace ParaMEDMEM
     virtual void setArrays(const std::vector<DataArrayDouble *>& arrays, TimeLabel *owner) throw(INTERP_KERNEL::Exception);
     DataArrayDouble *getArray() const { return _array; }
     virtual DataArrayDouble *getEndArray() const { return _array; }
-    //! Warning contrary to getArray method this method returns an object to deal with.
-    virtual DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) = 0;
+    virtual std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception) = 0;
+    virtual void getValueForTime(double time, const std::vector<double>& vals, double *res) const = 0; 
     virtual void getArrays(std::vector<DataArrayDouble *>& arrays) const;
     virtual bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
     virtual bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
@@ -104,7 +104,8 @@ namespace ParaMEDMEM
     MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const;
     void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) { }
     void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception);
-    DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+    std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception);
+    void getValueForTime(double time, const std::vector<double>& vals, double *res) const;
     bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
     bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
     double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception);
@@ -144,6 +145,8 @@ namespace ParaMEDMEM
     double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; }
     double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; }
     DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+    std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception);
+    void getValueForTime(double time, const std::vector<double>& vals, double *res) const;
     void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception);
     void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception);
   public:
@@ -168,7 +171,8 @@ namespace ParaMEDMEM
     MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const;
     bool areCompatible(const MEDCouplingTimeDiscretization *other) const;
     bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const;
-    DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+    std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception);
+    void getValueForTime(double time, const std::vector<double>& vals, double *res) const;
     void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception);
     void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception);
     TypeOfTimeDiscretization getEnum() const { return DISCRETIZATION; }
diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 86ad3e08e..e76f933e5 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -21,6 +21,8 @@
 #include "CellModel.hxx"
 #include "VolSurfUser.txx"
 #include "InterpolationUtils.hxx"
+#include "PointLocatorAlgos.txx"
+#include "BBTree.txx"
 
 #include <sstream>
 #include <limits>
@@ -907,6 +909,140 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const
   return ret;
 }
 
+/*{
+    const double *coordsPtr=_coords->getConstPointer();
+    BBTree<SPACEDIM,int> myTree(&bbox[0],0,0,nbNodes,-prec);
+    double bb[2*SPACEDIM];
+    double prec2=prec*prec;
+    for(int i=0;i<nbNodes;i++)
+      {
+        if(std::find(c.begin(),c.end(),i)!=c.end())
+          continue;
+        for(int j=0;j<SPACEDIM;j++)
+          {
+            bb[2*j]=coordsPtr[SPACEDIM*i+j];
+            bb[2*j+1]=coordsPtr[SPACEDIM*i+j];
+          }
+        std::vector<int> intersectingElems;
+        myTree.getIntersectingElems(bb,intersectingElems);
+        if(intersectingElems.size()>1)
+          {
+            std::vector<int> commonNodes;
+            for(std::vector<int>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
+              if(*it!=i)
+                if(INTERP_KERNEL::distance2<SPACEDIM>(coordsPtr+SPACEDIM*i,coordsPtr+SPACEDIM*(*it))<prec2)
+                  commonNodes.push_back(*it);
+            if(!commonNodes.empty())
+              {
+                cI.push_back(cI.back()+commonNodes.size()+1);
+                c.push_back(i);
+                c.insert(c.end(),commonNodes.begin(),commonNodes.end());
+              }
+          }
+          }*/
+
+int MEDCouplingUMesh::getElementContainingPoint(const double *pos, double eps) const
+{
+  std::vector<int> elts;
+  getElementsContainingPoint(pos,eps,elts);
+  if(elts.empty())
+    return -1;
+  return elts.front();
+}
+
+void MEDCouplingUMesh::getElementsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const
+{
+  std::vector<int> eltsIndex;
+  getElementsContainingPoints(pos,1,eps,elts,eltsIndex);
+}
+
+namespace ParaMEDMEM
+{
+  template<const int SPACEDIMM>
+  class DummyClsMCUG
+  {
+  public:
+    static const int MY_SPACEDIM=SPACEDIMM;
+    static const int MY_MESHDIM=8;
+    typedef int MyConnType;
+    static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
+  };
+}
+
+template<int SPACEDIM>
+void MEDCouplingUMesh::getElementsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
+                                                      double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const
+{
+  std::vector<double> bbox;
+  eltsIndex.resize(nbOfPoints+1);
+  eltsIndex[0]=0;
+  elts.clear();
+  getBoundingBoxForBBTree(bbox);
+  int nbOfCells=getNumberOfCells();
+  const int *conn=_nodal_connec->getConstPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  double bb[2*SPACEDIM];
+  for(int j=0;j<SPACEDIM;j++)
+    { bb[2*j]=std::numeric_limits<double>::max(); bb[2*j+1]=-std::numeric_limits<double>::max(); }
+  BBTree<SPACEDIM,int> myTree(&bbox[0],0,0,nbOfCells,-eps);
+  for(int i=0;i<nbOfPoints;i++)
+    {
+      eltsIndex[i+1]=eltsIndex[i];
+      for(int j=0;j<SPACEDIM;j++)
+        {
+          bb[2*j]=std::min(bb[2*j],pos[SPACEDIM*i+j]);
+          bb[2*j+1]=std::max(bb[2*j+1],pos[SPACEDIM*i+j]);
+        }
+      std::vector<int> candidates;
+      myTree.getIntersectingElems(bb,candidates);
+      for(std::vector<int>::const_iterator iter=candidates.begin();iter!=candidates.end();iter++)
+        {
+          int sz=connI[(*iter)+1]-connI[*iter]-1;
+          if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCUG<SPACEDIM> >::isElementContainsPoint(pos+i*SPACEDIM,
+                                                                                               (INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]],
+                                                                                               coords,conn+connI[*iter]+1,sz,eps))
+            {
+              eltsIndex[i+1]++;
+              elts.push_back(*iter);
+            }
+        }
+    }
+}
+
+void MEDCouplingUMesh::getElementsContainingPoints(const double *pos, int nbOfPoints, double eps,
+                                                   std::vector<int>& elts, std::vector<int>& eltsIndex) const
+{
+  int spaceDim=getSpaceDimension();
+  int mDim=getMeshDimension();
+  if(spaceDim==3)
+    {
+      if(mDim==3)
+        {
+          const double *coords=_coords->getConstPointer();
+          getElementsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
+        }
+      /*else if(mDim==2)
+        {
+          
+        }*/
+      else
+        throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !");
+    }
+  else if(spaceDim==2)
+    {
+      if(mDim==2)
+        {
+          const double *coords=_coords->getConstPointer();
+          getElementsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
+        }
+      else
+        throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !");
+    }
+
+
+  
+}
+
 /*!
  * This method is only available for a mesh with meshDim==2 and spaceDim==2||spaceDim==3.
  * This method returns a vector 'cells' where all detected butterfly cells have been added to cells.
@@ -937,6 +1073,40 @@ void MEDCouplingUMesh::checkButterflyCells(std::vector<int>& cells) const
     }
 }
 
+/*!
+ * This method aggregate the bbox of each cell and put it into bbox parameter.
+ * @param bbox out parameter of size 2*spacedim*nbOfcells.
+ */
+void MEDCouplingUMesh::getBoundingBoxForBBTree(std::vector<double>& bbox) const
+{
+  int spaceDim=getSpaceDimension();
+  int nbOfCells=getNumberOfCells();
+  bbox.resize(2*nbOfCells*spaceDim);
+  for(int i=0;i<nbOfCells*spaceDim;i++)
+    {
+      bbox[2*i]=std::numeric_limits<double>::max();
+      bbox[2*i+1]=-std::numeric_limits<double>::max();
+    }
+  const double *coordsPtr=_coords->getConstPointer();
+  const int *conn=_nodal_connec->getConstPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  for(int i=0;i<nbOfCells;i++)
+    {
+      int offset=connI[i]+1;
+      int nbOfNodesForCell=connI[i+1]-offset;
+      for(int j=0;j<nbOfNodesForCell;j++)
+        {
+          int nodeId=conn[offset+j];
+          if(nodeId>=0)
+            for(int k=0;k<spaceDim;k++)
+              {
+                bbox[2*spaceDim*i+2*k]=std::min(bbox[2*spaceDim*i+2*k],coordsPtr[spaceDim*nodeId+k]);
+                bbox[2*spaceDim*i+2*k+1]=std::max(bbox[2*spaceDim*i+2*k+1],coordsPtr[spaceDim*nodeId+k]);
+              }
+        }
+    }
+}
+
 namespace ParaMEDMEMImpl
 {
   class ConnReader
diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx
index 9c036342f..e785cff54 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.hxx
@@ -72,7 +72,12 @@ namespace ParaMEDMEM
     MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
     MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
     MEDCouplingFieldDouble *buildOrthogonalField() const;
+    int getElementContainingPoint(const double *pos, double eps) const;
+    void getElementsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const;
+    void getElementsContainingPoints(const double *pos, int nbOfPoints, double eps,
+                                     std::vector<int>& elts, std::vector<int>& eltsIndex) const;
     void checkButterflyCells(std::vector<int>& cells) const;
+    void getBoundingBoxForBBTree(std::vector<double>& bbox) const;
     bool checkConsecutiveCellTypes() const;
     MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
     DataArrayDouble *getBarycenterAndOwner() const;
@@ -85,6 +90,9 @@ namespace ParaMEDMEM
     void checkFullyDefined() const throw(INTERP_KERNEL::Exception);
     //tools
     MEDCouplingUMesh *buildPartOfMySelfKeepCoords(const int *start, const int *end) const;
+    template<int SPACEDIM>
+    void getElementsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
+                                        double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const;
   private:
     //! this iterator stores current position in _nodal_connec array.
     mutable int _iterator;
diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
index 3a34d1f84..feca3dc3c 100644
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
@@ -312,3 +312,9 @@ DataArrayDouble *MEDCouplingUMeshDesc::getBarycenterAndOwner() const
   //not implemented yet.
   return 0;
 }
+
+int MEDCouplingUMeshDesc::getElementContainingPoint(const double *pos, double eps) const
+{
+  //not implemented yet.
+  return -1;
+}
diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
index ebd34115a..94876129b 100644
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
@@ -57,6 +57,7 @@ namespace ParaMEDMEM
     void renumberNodes(const int *newNodeNumbers, int newNbOfNodes);
     MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
     MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
+    int getElementContainingPoint(const double *pos, double eps) const;
     MEDCouplingFieldDouble *buildOrthogonalField() const;
     DataArrayInt *zipCoordsTraducer();
     MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
index c53720bed..0b8f9b70d 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
@@ -1473,6 +1473,102 @@ void MEDCouplingBasicsTest::testBuildOrthogonalField()
   targetMesh->decrRef();
 }
 
+void MEDCouplingBasicsTest::testGetElementsContainingPoint()
+{
+  MEDCouplingUMesh *targetMesh=build2DTargetMesh_1();
+  double pos[12]={0.,0.,0.4,0.4,0.,0.4,0.1,0.1,0.25,0.,0.65,0.};
+  std::vector<int> t1,t2;
+  //2D basic
+  targetMesh->getElementsContainingPoints(pos,6,1e-12,t1,t2);
+  CPPUNIT_ASSERT_EQUAL(6,(int)t1.size());
+  CPPUNIT_ASSERT_EQUAL(7,(int)t2.size());
+  const int expectedValues1[6]={0,4,3,0,1,2};
+  const int expectedValues2[7]={0,1,2,3,4,5,6};
+  CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues1));
+  CPPUNIT_ASSERT(std::equal(t2.begin(),t2.end(),expectedValues2));
+  //2D with no help of bounding box.
+  double center[2]={0.2,0.2};
+  MEDCouplingPointSet::rotate2DAlg(center,0.78539816339744830962,6,pos);
+  targetMesh->rotate(center,0,0.78539816339744830962);
+  t1.clear(); t2.clear();
+  targetMesh->getElementsContainingPoints(pos,6,1e-12,t1,t2);
+  CPPUNIT_ASSERT_EQUAL(6,(int)t1.size());
+  CPPUNIT_ASSERT_EQUAL(7,(int)t2.size());
+  CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues1));
+  CPPUNIT_ASSERT(std::equal(t2.begin(),t2.end(),expectedValues2));
+  //2D outside
+  const double pos1bis[2]={-0.3303300858899107,-0.11819805153394641};
+  CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos1bis,1e-12));
+  targetMesh->decrRef();
+  //test limits 2D
+  targetMesh=build2DTargetMesh_1();
+  const double pos2[2]={0.2,-0.05};
+  t1.clear();
+  targetMesh->getElementsContainingPoint(pos2,1e-12,t1);
+  CPPUNIT_ASSERT_EQUAL(2,(int)t1.size());
+  const int expectedValues3[2]={0,1};
+  CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues3));
+  const double pos3[2]={0.2,0.2};
+  t1.clear();
+  targetMesh->getElementsContainingPoint(pos3,1e-12,t1);
+  CPPUNIT_ASSERT_EQUAL(5,(int)t1.size());
+  const int expectedValues4[5]={0,1,2,3,4};
+  CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues4));
+  CPPUNIT_ASSERT_EQUAL(0,targetMesh->getElementContainingPoint(pos3,1e-12));
+  targetMesh->decrRef();
+  //3D
+  targetMesh=build3DTargetMesh_1();
+  const double pos4[3]={25.,25.,25.};
+  CPPUNIT_ASSERT_EQUAL(0,targetMesh->getElementContainingPoint(pos4,1e-12));
+  const double pos5[3]={50.,50.,50.};
+  t1.clear();
+  targetMesh->getElementsContainingPoint(pos5,1e-12,t1);
+  CPPUNIT_ASSERT_EQUAL(8,(int)t1.size());
+  const int expectedValues5[8]={0,1,2,3,4,5,6,7};
+  CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues5));
+  const double pos6[3]={0., 50., 0.};
+  t1.clear();
+  targetMesh->getElementsContainingPoint(pos6,1e-12,t1);
+  CPPUNIT_ASSERT_EQUAL(2,(int)t1.size());
+  const int expectedValues6[2]={0,2};
+  CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues6));
+  //3D outside
+  const double pos7[3]={-1.0,-1.0,0.};
+  CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos7,1e-12));
+  //3D outside 2
+  const double center2[3]={0.,0.,0.};
+  const double vec2[3]={0.,-1.,0.};
+  targetMesh->rotate(center2,vec2,0.78539816339744830962);
+  const double pos8[3]={-25,25.,12.};
+  CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos8,1e-12));
+  //
+  targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::testGetValueOn1()
+{
+  MEDCouplingUMesh *targetMesh=build2DTargetMesh_1();
+  MEDCouplingFieldDouble *fieldOnCells=MEDCouplingFieldDouble::New(ON_CELLS);
+  int nbOfCells=targetMesh->getNumberOfCells();
+  fieldOnCells->setMesh(targetMesh);
+  DataArrayDouble *array=DataArrayDouble::New();
+  array->alloc(nbOfCells,2);
+  fieldOnCells->setArray(array);
+  double *tmp=array->getPointer();
+  for(int i=0;i<nbOfCells;i++)
+    { tmp[2*i]=7.+(double)i; tmp[2*i+1]=17.+(double)i; }
+  array->decrRef();
+  //
+  const double pos1[2]={0.25,0.};
+  double res[2];
+  fieldOnCells->getValueOn(pos1,res);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(8.,res[0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(18.,res[1],1e-12);
+  //
+  fieldOnCells->decrRef();
+  targetMesh->decrRef();
+}
+
 void MEDCouplingBasicsTest::test2DInterpP0P0_1()
 {
   MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1();
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
index f2c054659..d1a50bfea 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
@@ -62,6 +62,8 @@ namespace ParaMEDMEM
     CPPUNIT_TEST( testMergeNodesOnField );
     CPPUNIT_TEST( testCheckConsecutiveCellTypes );
     CPPUNIT_TEST( testBuildOrthogonalField );
+    CPPUNIT_TEST( testGetElementsContainingPoint );
+    CPPUNIT_TEST( testGetValueOn1 );
     CPPUNIT_TEST( test2DInterpP0P0_1 );
     CPPUNIT_TEST( test2DInterpP0P0PL_1 );
     CPPUNIT_TEST( test2DInterpP0P0PL_2 );
@@ -143,6 +145,8 @@ namespace ParaMEDMEM
     void testMergeNodesOnField();
     void testCheckConsecutiveCellTypes();
     void testBuildOrthogonalField();
+    void testGetElementsContainingPoint();
+    void testGetValueOn1();
     void test2DInterpP0P0_1();
     void test2DInterpP0P0PL_1();
     void test2DInterpP0P0PL_2();