]> SALOME platform Git repositories - modules/med.git/commitdiff
Salome HOME
PointLocator implementation for P1P1 for 1D and 2D curve.
authorAnthony Geay <anthony.geay@edf.fr>
Tue, 21 Jul 2015 07:28:25 +0000 (09:28 +0200)
committerAnthony Geay <anthony.geay@edf.fr>
Tue, 21 Jul 2015 07:28:25 +0000 (09:28 +0200)
src/INTERP_KERNEL/CurveIntersector.hxx
src/INTERP_KERNEL/CurveIntersector.txx
src/INTERP_KERNEL/CurveIntersectorP1P1PL.hxx [new file with mode: 0644]
src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx [new file with mode: 0644]
src/INTERP_KERNEL/InterpolationCurve.txx
src/MEDCoupling_Swig/MEDCouplingRemapperTest.py

index f2f846e1a5fadff954dbf480c7786156c8b4fcb7..e1c3e69cec3076eb620fbc6c91abeb889d87335d 100644 (file)
@@ -41,12 +41,16 @@ namespace INTERP_KERNEL
     void createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox);
     void adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEpsAbs);
     static void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes);
-
+    static bool ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos);
   protected :
-    bool getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT);
-    bool getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS);
-    double intersectSegments(double *Coords_T, double *Coords_S);
-
+    bool projectionThis(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const;
+    bool getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const;
+    typename MyMeshType::MyConnType getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const;
+    bool getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS) const;
+    typename MyMeshType::MyConnType getNodeIdOfSourceCellAt(ConnType icellT, ConnType nodeIdInCellT) const;
+    double intersectSegments(const double *coordsT, const double *coordsS) const;
+    double intersectSegmentsInternal(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const;
+    
     struct TDualSegment
     {
       std::vector<double> _coords;
index 3dc7ff60ad2e4a85d6aea6dd5f31a9151491b98a..c311605f7e7c3067567cf4ef7e5bad05a946e572 100644 (file)
@@ -103,12 +103,9 @@ namespace INTERP_KERNEL
       }
   }
 
-  //================================================================================
   /*!
     Computes the bouding box of a given element. iP in numPol mode.
   */
-  //================================================================================
-
   template<class MyMeshType, class MyMatrix>
   void CurveIntersector<MyMeshType,MyMatrix>::getElemBB (double*           bb,
                                                          const MyMeshType& mesh,
@@ -138,15 +135,26 @@ namespace INTERP_KERNEL
           }
       }
   }
+  
+  /*!
+   * \param [in] startOfSeg - input coming from intersectSegments or intersectSegmentsInternal
+   * \param [in] endOfSeg - input coming from intersectSegments or intersectSegmentsInternal. Assume that endOfSeg>startOfSeg.
+   * \param [in] pt - position of point that the method computes the bary coords for.
+   */
+  template<class MyMeshType, class MyMatrix>
+  bool CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
+  {
+    double deno(endOfSeg-startOfSeg);
+    startPos=(endOfSeg-pt)/deno;
+    endPos=1.-startPos;
+    return startPos>=0. && endPos>=0.;
+  }
 
-  //================================================================================
   /*! Readjusts a set of bounding boxes so that they are extended
     in all dimensions for avoiding missing interesting intersections
 
     \param bbox vector containing the bounding boxes
   */
-  //================================================================================
-
   template<class MyMeshType, class MyMatrix>
   void CurveIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes (std::vector<double>& bbox,
                                                                    double adjustmentEpsAbs)
@@ -162,7 +170,6 @@ namespace INTERP_KERNEL
       }
   }
 
-  //================================================================================
   /*!
    * @param icellT id in target mesh in format of MyMeshType.
    * @param coordsT output val that stores coordinates of the target cell
@@ -170,14 +177,10 @@ namespace INTERP_KERNEL
    * @return true if segment is quadratic and in this case coordinates of medium node
    * are placed in the middle of coordsT
    */
-  //================================================================================
-
   template<class MyMeshType, class MyMatrix>
-  bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates
-  (ConnType icellT, std::vector<double>& coordsT)
+  bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const
   {
-    int nbNodesT = _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] -
-      _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+    int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
     coordsT.resize(SPACEDIM*nbNodesT);
     for (ConnType iT=0; iT<nbNodesT; iT++)
       {
@@ -195,22 +198,27 @@ namespace INTERP_KERNEL
       }
     return false;
   }
+  
+  template<class MyMeshType, class MyMatrix>
+  typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const
+  {
+    int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
+    if(nodeIdInCellT>=0 && nodeIdInCellT<nbNodesT)
+      return OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+nodeIdInCellT)]);
+    else
+      throw Exception("getNodeIdOfTargetCellAt : error in nodeId in cell");
+  }
 
-  //================================================================================
   /*!
    * @param icellS id in source mesh in format of MyMeshType.
    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
    * @return true if segment is quadratic and in this case coordinates of medium node
    * are placed in the middle of coordsS
    */
-  //================================================================================
-
   template<class MyMeshType, class MyMatrix>
-  bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates
-  (ConnType icellS, std::vector<double>& coordsS)
+  bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS) const
   {
-    int nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] -
-      _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+    int nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
     coordsS.resize(SPACEDIM*nbNodesS);
     for(ConnType iS=0; iS<nbNodesS; iS++)
       {
@@ -229,15 +237,22 @@ namespace INTERP_KERNEL
     return false;
   }
 
-  //================================================================================
+  template<class MyMeshType, class MyMatrix>
+  typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(ConnType icellS, ConnType nodeIdInCellS) const
+  {
+    int nbNodesS(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]);
+    if(nodeIdInCellS>=0 && nodeIdInCellS<nbNodesS)
+      return OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+nodeIdInCellS)]);
+    else
+      throw Exception("getNodeIdOfSourceCellAt : error in nodeId in cell");
+  }
+
   /*!
    * \brief Return dual segments of given segment
    *  \param icell - given segment in C mode
    *  \param mesh - mesh
    *  \param segments - dual segments
    */
-  //================================================================================
-
   template<class MyMeshType, class MyMatrix>
   void CurveIntersector<MyMeshType,MyMatrix>::getDualSegments(ConnType                   icell,
                                                               const MyMeshType&          mesh,
@@ -300,18 +315,12 @@ namespace INTERP_KERNEL
       }
   }
 
-  //================================================================================
-  /*!
-   * \brief Return length of intersection of two segments
-   */
-  //================================================================================
-
   template<class MyMeshType, class MyMatrix>
-  double CurveIntersector<MyMeshType,MyMatrix>::intersectSegments(double *Coords_T,
-                                                                  double *Coords_S)
+  bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS,
+                                                             double& xs0, double& xs1, double& xt0, double& xt1) const
   {
-    double xt0 = Coords_T[0], xt1 = Coords_T[1];
-    double xs0 = Coords_S[0], xs1 = Coords_S[1];
+    xt0 = coordsT[0]; xt1 = coordsT[1];
+    xs0 = coordsS[0]; xs1 = coordsS[1];
     if ( SPACEDIM == 2 )
       {
         // Pass 2D->1D
@@ -320,15 +329,16 @@ namespace INTERP_KERNEL
 
         // check if two segments overlap in 2D within tolerance
 
-        double* t0 = Coords_T;
-        double* t1 = Coords_T + 2;
+        const double* t0 = coordsT;
+        const double* t1 = coordsT + 2;
         double t01[2] = { t1[X]-t0[X], t1[Y]-t0[Y] }; // tgt segment direction
         double tSize = sqrt( t01[X]*t01[X] + t01[Y]*t01[Y] ); // tgt segment size
-        if ( tSize < _precision ) return 0; // degenerated segment
+        if ( tSize < _precision )
+          return false; // degenerated segment
         t01[X] /= tSize, t01[Y] /= tSize; // normalize t01
 
-        double* s0 = Coords_S;
-        double* s1 = Coords_S + 2;
+        const double* s0 = coordsS;
+        const double* s1 = coordsS + 2;
         double t0s0[2] = { s0[X]-t0[X], s0[Y]-t0[Y] };
         double t0s1[2] = { s1[X]-t0[X], s1[Y]-t0[Y] };
         double nt01_x_t0s0 = t0s0[X] * t01[Y] - t0s0[Y] * t01[X]; // t0s0 dot norm of t01
@@ -338,7 +348,7 @@ namespace INTERP_KERNEL
         bool s0_out_of_tol = ( dist_ts0 > _tolerance );
         bool s1_out_of_tol = ( dist_ts1 > _tolerance );
         if ( nt01_x_t0s0 * nt01_x_t0s1 > 0 && ( s0_out_of_tol || s1_out_of_tol ))
-          return 0; // tgt segment is to far from src segment
+          return false; // tgt segment is to far from src segment
 
         double S0[2] = { s0[X], s0[Y] };
         double S1[2] = { s1[X], s1[Y] };
@@ -361,7 +371,8 @@ namespace INTERP_KERNEL
 
         double s01[2] = { S1[X]-S0[X], S1[Y]-S0[Y] }; // src segment direction
         double sSize = sqrt( s01[X]*s01[X] + s01[Y]*s01[Y] ); // src segment size
-        if ( sSize < _precision ) return 0; // degenerated segment
+        if ( sSize < _precision )
+          return false; // degenerated segment
         s01[X] /= sSize, s01[Y] /= sSize; // normalize s01
 
         // make t01 and s01 codirected
@@ -375,7 +386,7 @@ namespace INTERP_KERNEL
         };
         double medianSize = sqrt( medianDir[X]*medianDir[X] + medianDir[Y]*medianDir[Y] );
         if ( medianSize < std::numeric_limits<double>::min() )
-          return 0; // strange...
+          return false; // strange...
         medianDir[X] /= medianSize, medianDir[Y] /= medianSize;
 
         xt0 = t0[X] * medianDir[X] + t0[Y] * medianDir[Y];
@@ -384,7 +395,18 @@ namespace INTERP_KERNEL
         xs1 = S1[X] * medianDir[X] + S1[Y] * medianDir[Y];
 
       } // if ( SPACEDIM == 2 )
-
+    return true;
+  }
+  
+  /*!
+   * \brief Return length of intersection of two segments
+   */
+  template<class MyMeshType, class MyMatrix>
+  double CurveIntersector<MyMeshType,MyMatrix>::intersectSegmentsInternal(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const
+  {
+    if(!projectionThis(coordsT,coordsS,xs0,xs1,xt0,xt1))
+      return 0.;
+    
     if ( xt0 > xt1 ) std::swap( xt0, xt1 );
     if ( xs0 > xs1 ) std::swap( xs0, xs1 );
 
@@ -392,6 +414,16 @@ namespace INTERP_KERNEL
     double x1 = std::min( xt1, xs1 );
     return ( x0 < x1 ) ? ( x1 - x0 ) : 0.;
   }
+  
+  /*!
+   * \brief Return length of intersection of two segments
+   */
+  template<class MyMeshType, class MyMatrix>
+  double CurveIntersector<MyMeshType,MyMatrix>::intersectSegments(const double *coordsT, const double *coordsS) const
+  {
+    double xs0,xs1,xt0,xt1;
+    return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1);
+  }
 
 }
 
diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P1PL.hxx b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.hxx
new file mode 100644 (file)
index 0000000..d9bab05
--- /dev/null
@@ -0,0 +1,50 @@
+// Copyright (C) 2007-2015  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
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// Author : Anthony Geay (EDF R&D)
+
+#ifndef __CURVEINTERSECTORP1P1PL_HXX__
+#define __CURVEINTERSECTORP1P1PL_HXX__
+
+#include "CurveIntersector.hxx"
+
+namespace INTERP_KERNEL
+{
+  template<class MyMeshType, class MyMatrix>
+  class CurveIntersectorP1P1PL : public CurveIntersector<MyMeshType,MyMatrix>
+  {
+  public:
+    static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
+    static const int MESHDIM=MyMeshType::MY_MESHDIM;
+    typedef typename MyMeshType::MyConnType ConnType;
+    static const NumberingPolicy numPol=MyMeshType::My_numPol;
+
+    CurveIntersectorP1P1PL(const MyMeshType& meshT, const MyMeshType& meshS,
+                           double  precision, double tolerance,
+                           double medianLine, int printLevel);
+  public:
+    void intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res);
+    int getNumberOfRowsOfResMatrix() const;
+    int getNumberOfColsOfResMatrix() const;
+  private:
+    static void AppendValueInMatrix(MyMatrix& res, ConnType nodeIdT, ConnType nodeIdS0, double val0, ConnType nodeIdS1, double val1);
+    static void AppendValueInMatrix2(typename MyMatrix::value_type& resRow, ConnType nodeIdS0, double val0);
+  };
+}
+
+#endif
diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx
new file mode 100644 (file)
index 0000000..02b09eb
--- /dev/null
@@ -0,0 +1,112 @@
+// Copyright (C) 2007-2015  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
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// Author : Anthony Geay (EDF R&D)
+
+#ifndef __CURVEINTERSECTORP1P1PL_TXX__
+#define __CURVEINTERSECTORP1P1PL_TXX__
+
+#include "CurveIntersectorP1P1PL.hxx"
+#include "CurveIntersector.txx"
+
+#include <cassert>
+
+namespace INTERP_KERNEL
+{
+  template<class MyMeshType, class MyMatrix>
+  CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::CurveIntersectorP1P1PL(const MyMeshType& meshT, const MyMeshType& meshS, double precision, double tolerance, double medianLine, int printLevel):CurveIntersector<MyMeshType,MyMatrix>(meshT, meshS, precision, tolerance, medianLine, printLevel)
+  {
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  int CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfRowsOfResMatrix() const
+  {
+    return CurveIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfNodes();
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  int CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfColsOfResMatrix() const
+  {
+    return CurveIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  void CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::AppendValueInMatrix2(typename MyMatrix::value_type& resRow, ConnType nodeIdS0, double val0)
+  {
+    typename MyMatrix::value_type::const_iterator iterRes(resRow.find(OTT<ConnType,numPol>::indFC(nodeIdS0)));    
+    if(iterRes==resRow.end())
+      {
+        resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(nodeIdS0),val0));
+      }
+    else
+      {
+        double val((*iterRes).second+val0);
+        resRow.erase(OTT<ConnType,numPol>::indFC(nodeIdS0));
+        resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(nodeIdS0),val));
+      }
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  void CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::AppendValueInMatrix(MyMatrix& res, ConnType nodeIdT, ConnType nodeIdS0, double val0, ConnType nodeIdS1, double val1)
+  {
+    typename MyMatrix::value_type& resRow(res[nodeIdT]);
+    AppendValueInMatrix2(resRow,nodeIdS0,val0);
+    AppendValueInMatrix2(resRow,nodeIdS1,val1);
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  void CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res)
+  {
+    std::vector<double> coordsT;
+    if(CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(icellT,coordsT))
+      throw INTERP_KERNEL::Exception("Invalid target cell detected for meshdim==1. Only SEG2 supported !");
+    assert(coordsT.size()/SPACEDIM==2);
+    for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
+      {
+        std::vector<double> coordsS;
+        if(CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(*iter,coordsS))
+          throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==1. Only SEG2 supported !");
+        assert(coordsS.size()/SPACEDIM==2);
+        double xs0,xs1,xt0,xt1;
+        double lgth(CurveIntersector<MyMeshType,MyMatrix>::intersectSegmentsInternal(&coordsT[0],&coordsS[0],xs0,xs1,xt0,xt1));
+        ConnType nodeIdS0(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(*iter,0));
+        ConnType nodeIdS1(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(*iter,1));
+        if(lgth>0.)
+          {
+            double a,b;
+            // for first
+            ConnType nodeIdT0(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(icellT,0));
+            if(CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(xs0,xs1,xt0,a,b))
+              {
+                a*=lgth; b*=lgth;
+                AppendValueInMatrix(res,nodeIdT0,nodeIdS0,a,nodeIdS1,b);
+              }
+            //
+            ConnType nodeIdT1(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(icellT,1));
+            typename MyMatrix::value_type& resRow1=res[nodeIdT1];
+            if(CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(xs0,xs1,xt1,a,b))
+              {
+                a*=lgth; b*=lgth;
+                AppendValueInMatrix(res,nodeIdT1,nodeIdS0,a,nodeIdS1,b);
+              }
+          }
+      }
+  }
+}
+
+#endif
index 35e78917290425591dfe5e7bd5ec0ff435b28513..239370a61ee5b85ea6f1f4e303dbc92826c17f66 100644 (file)
@@ -26,6 +26,7 @@
 #include "CurveIntersectorP1P0.txx"
 #include "CurveIntersectorP0P1.txx"
 #include "CurveIntersectorP1P1.txx"
+#include "CurveIntersectorP1P1PL.txx"
 #include "BBTree.txx"
 
 #include <time.h>
@@ -90,39 +91,78 @@ namespace INTERP_KERNEL
     CurveIntersector<MyMeshType,MatrixType>* intersector=0;
     if(method=="P0P0")
       {
-        intersector = new CurveIntersectorP0P0<MyMeshType,MatrixType>
-          (myMeshT, myMeshS,
-           InterpolationOptions::getPrecision(),
-           InterpolationOptions::getBoundingBoxAdjustmentAbs(),
-           InterpolationOptions::getMedianPlane(),
-           InterpolationOptions::getPrintLevel());
+        switch (InterpolationOptions::getIntersectionType())
+          {
+            case Triangulation:
+              {
+                intersector = new CurveIntersectorP0P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
+                                                                              InterpolationOptions::getPrecision(),
+                                                                              InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+                                                                              InterpolationOptions::getMedianPlane(),
+                                                                              InterpolationOptions::getPrintLevel());
+                break;
+              }
+            default:
+              throw INTERP_KERNEL::Exception("For P0P0 in 1D or 2D curve only Triangulation supported for the moment !");
+          }
       }
     else if(method=="P0P1")
       {
-        intersector = new CurveIntersectorP0P1<MyMeshType,MatrixType>
-          (myMeshT, myMeshS,
-           InterpolationOptions::getPrecision(),
-           InterpolationOptions::getBoundingBoxAdjustmentAbs(),
-           InterpolationOptions::getMedianPlane(),
-           InterpolationOptions::getPrintLevel());
+        switch (InterpolationOptions::getIntersectionType())
+          {
+            case Triangulation:
+              {
+                intersector = new CurveIntersectorP0P1<MyMeshType,MatrixType>(myMeshT, myMeshS,
+                                                                              InterpolationOptions::getPrecision(),
+                                                                              InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+                                                                              InterpolationOptions::getMedianPlane(),
+                                                                              InterpolationOptions::getPrintLevel());
+                break;
+              }
+            default:
+              throw INTERP_KERNEL::Exception("For P0P1 in 1D or 2D curve only Triangulation supported for the moment !");
+          }
       }
     else if(method=="P1P0")
       {
-        intersector = new CurveIntersectorP1P0<MyMeshType,MatrixType>
-          (myMeshT, myMeshS,
-           InterpolationOptions::getPrecision(),
-           InterpolationOptions::getBoundingBoxAdjustmentAbs(),
-           InterpolationOptions::getMedianPlane(),
-           InterpolationOptions::getPrintLevel());
+        switch (InterpolationOptions::getIntersectionType())
+          {
+            case Triangulation:
+              {
+                intersector = new CurveIntersectorP1P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
+                                                                              InterpolationOptions::getPrecision(),
+                                                                              InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+                                                                              InterpolationOptions::getMedianPlane(),
+                                                                              InterpolationOptions::getPrintLevel());
+                break;
+              }
+            default:
+              throw INTERP_KERNEL::Exception("For P1P0 in 1D or 2D curve only Triangulation supported for the moment !");
+          }
       }
     else if(method=="P1P1")
       {
-        intersector = new CurveIntersectorP1P1<MyMeshType,MatrixType>
-          (myMeshT, myMeshS,
-           InterpolationOptions::getPrecision(),
-           InterpolationOptions::getBoundingBoxAdjustmentAbs(),
-           InterpolationOptions::getMedianPlane(),
-           InterpolationOptions::getPrintLevel());
+        switch (InterpolationOptions::getIntersectionType())
+          {
+          case Triangulation:
+            intersector = new CurveIntersectorP1P1<MyMeshType,MatrixType>
+              (myMeshT, myMeshS,
+               InterpolationOptions::getPrecision(),
+               InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+               InterpolationOptions::getMedianPlane(),
+               InterpolationOptions::getPrintLevel());
+            break;
+          case PointLocator:
+            intersector = new CurveIntersectorP1P1PL<MyMeshType,MatrixType>
+              (myMeshT, myMeshS,
+               InterpolationOptions::getPrecision(),
+               InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+               InterpolationOptions::getMedianPlane(),
+               InterpolationOptions::getPrintLevel());
+            break;
+          default:
+            throw INTERP_KERNEL::Exception("For P1P1 in 1D or 2D curve only Triangulation and PointLocator supported !");
+          }
       }
     else
       throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
index 469f924631338e936bd45c2964362a555f23959a..ba5ba6428de3f48e3b263fce123994e59c474a8d 100644 (file)
@@ -786,6 +786,74 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertTrue(coarse.isEqual(trgField.getArray(),1e-12))
         pass
     
+    @unittest.skipUnless(MEDCouplingHasNumPyBindings() and MEDCouplingHasSciPyBindings(),"requires numpy AND scipy")
+    def test1DPointLocator1(self):
+        """This test focuses on PointLocator for P1P1 in 1D and 2DCurve."""
+        from numpy import array
+        from scipy.sparse import diags,csr_matrix,identity
+        ## basic case 1D
+        arrS=DataArrayInt.Range(0,11,1).convertToDblArr()
+        arrT=DataArrayDouble([0.1,1.7,5.5,9.6])
+        mS=MEDCouplingCMesh() ; mS.setCoords(arrS)
+        mT=MEDCouplingCMesh() ; mT.setCoords(arrT)
+        rem=MEDCouplingRemapper()
+        rem.setIntersectionType(PointLocator)
+        self.assertEqual(rem.prepare(mS.buildUnstructured(),mT.buildUnstructured(),"P1P1"),1)
+        m=rem.getCrudeCSRMatrix()
+        rowSum=m.sum(axis=1)
+        m=diags(array(1/rowSum.transpose()),[0])*m
+        # expected matrix
+        row=array([0,0,1,1,2,2,3,3])
+        col=array([0,1,1,2,5,6,9,10])
+        data=array([0.9,0.1,0.3,0.7,0.5,0.5,0.4,0.6])
+        mExp0=csr_matrix((data,(row,col)),shape=(4,11))
+        # compute diff and check
+        diff=abs(m-mExp0)
+        self.assertAlmostEqual(diff.max(),0.,14)
+        ## full specific case 1D where target=source
+        rem=MEDCouplingRemapper()
+        rem.setIntersectionType(PointLocator)
+        self.assertEqual(rem.prepare(mS.buildUnstructured(),mS.buildUnstructured(),"P1P1"),1)
+        m=rem.getCrudeCSRMatrix()
+        rowSum=m.sum(axis=1)
+        m=diags(array(1/rowSum.transpose()),[0])*m
+        # expected matrix
+        mExp1=identity(11)
+        diff=abs(m-mExp1)
+        self.assertAlmostEqual(diff.max(),0.,14)
+        ## case where some points in target are not in source
+        arrT=DataArrayDouble([-0.2,0.1,1.7,5.5,10.3])
+        mT=MEDCouplingCMesh() ; mT.setCoords(arrT)
+        mT=mT.buildUnstructured()
+        rem=MEDCouplingRemapper()
+        rem.setIntersectionType(PointLocator)
+        self.assertEqual(rem.prepare(mS.buildUnstructured(),mT,"P1P1"),1)
+        m=rem.getCrudeCSRMatrix()
+        row=array([1,1,2,2,3,3])
+        col=array([0,1,1,2,5,6])
+        data=array([0.9,0.1,0.3,0.7,0.5,0.5])
+        mExp2=csr_matrix((data,(row,col)),shape=(5,11))
+        diff=abs(m-mExp2)
+        self.assertAlmostEqual(diff.max(),0.,14)
+        ## basic case 2D Curve
+        arrS=DataArrayInt.Range(0,11,1).convertToDblArr()
+        arrT=DataArrayDouble([0.1,1.7,5.5,9.6])
+        mS=MEDCouplingCMesh() ; mS.setCoords(arrS)
+        mT=MEDCouplingCMesh() ; mT.setCoords(arrT)
+        mS=mS.buildUnstructured() ; mS.changeSpaceDimension(2)
+        mT=mT.buildUnstructured() ; mT.changeSpaceDimension(2)
+        mS.rotate([-1.,-1.],1.2)
+        mT.rotate([-1.,-1.],1.2)
+        rem=MEDCouplingRemapper()
+        rem.setIntersectionType(PointLocator)
+        self.assertEqual(rem.prepare(mS,mT,"P1P1"),1)
+        m=rem.getCrudeCSRMatrix()
+        rowSum=m.sum(axis=1)
+        m=diags(array(1/rowSum.transpose()),[0])*m
+        diff=abs(m-mExp0)
+        self.assertAlmostEqual(diff.max(),0.,14)
+        pass
+    
     def build2DSourceMesh_1(self):
         sourceCoords=[-0.3,-0.3, 0.7,-0.3, -0.3,0.7, 0.7,0.7]
         sourceConn=[0,3,1,0,2,3]