src/INTERP_KERNEL/CurveIntersector.txx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #pragma once
  22
  23 #include "CurveIntersector.hxx"
  24 #include "InterpolationUtils.hxx"
  25 #include "PointLocatorAlgos.txx"
  26
  27 #include <limits>
  28
  29 namespace INTERP_KERNEL
  30 {
  31   template<class MyMeshType, class MyMatrix>
  32   CurveIntersector<MyMeshType,MyMatrix>
  33   ::CurveIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
  34                      double precision, double tolerance, double medianLine, int printLevel):
  35     _meshT(meshT),
  36     _meshS(meshS),
  37     _tolerance(tolerance),
  38     _precision(precision),
  39     _median_line(medianLine),
  40     _print_level(printLevel)
  41   {
  42     if ( SPACEDIM != 1 && SPACEDIM != 2 )
  43       throw Exception("CurveIntersector(): space dimension of mesh must be 1 or 2");
  44     if ( MESHDIM != 1 )
  45       throw Exception("CurveIntersector(): mesh dimension must be 1");
  46
  47     _connectT = meshT.getConnectivityPtr();
  48     _connectS = meshS.getConnectivityPtr();
  49     _connIndexT = meshT.getConnectivityIndexPtr();
  50     _connIndexS = meshS.getConnectivityIndexPtr();
  51     _coordsT = meshT.getCoordinatesPtr();
  52     _coordsS = meshS.getCoordinatesPtr();
  53   }
  54
  55   template<class MyMeshType, class MyMatrix>
  56   CurveIntersector<MyMeshType,MyMatrix>::~CurveIntersector()
  57   {
  58   }
  59
  60   //================================================================================
  61   /*!
  62     \brief creates the bounding boxes for all the cells of mesh \a mesh
  63
  64     \param mesh structure pointing to the mesh
  65     \param bbox vector containing the bounding boxes
  66   */
  67   //================================================================================
  68
  69   template<class MyMeshType, class MyMatrix>
  70   void CurveIntersector<MyMeshType,MyMatrix>::createBoundingBoxes (const MyMeshType&    mesh,
  71                                                                    std::vector<double>& bbox)
  72   {
  73     long nbelems = mesh.getNumberOfElements();
  74     bbox.resize(2*SPACEDIM* nbelems);
  75     const double* coords = mesh.getCoordinatesPtr();
  76     const ConnType* conn = mesh.getConnectivityPtr();
  77     const ConnType* conn_index = mesh.getConnectivityIndexPtr();
  78     int ibox=0;
  79     for(long icell=0; icell<nbelems; icell++)
  80       {
  81         ConnType nb_nodes_per_elem = conn_index[icell+1]-conn_index[icell];
  82         //initializing bounding box limits
  83         for(int idim=0; idim<SPACEDIM; idim++)
  84           {
  85             bbox[2*SPACEDIM*ibox+2*idim]   =  std::numeric_limits<double>::max();
  86             bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
  87           }
  88         //updating the bounding box with each node of the element
  89         for (ConnType j=0; j<nb_nodes_per_elem; j++)
  90           {
  91             const double* coord_node = coords +
  92               SPACEDIM*OTT<ConnType,numPol>
  93               ::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[icell]+j)]);
  94             for(int idim=0; idim<SPACEDIM; idim++)
  95               {
  96                 double x = *(coord_node+idim);
  97                 bbox[ibox*2*SPACEDIM + 2*idim]   =
  98                   ( bbox[ibox*2*SPACEDIM + 2*idim] < x ) ? bbox[ibox*2*SPACEDIM + 2*idim] : x;
  99                 bbox[ibox*2*SPACEDIM + 2*idim+1] =
 100                   ( bbox[ibox*2*SPACEDIM + 2*idim+1] > x ) ? bbox[ibox*2*SPACEDIM + 2*idim+1] : x;
 101               }
 102           }
 103         ibox++;
 104       }
 105   }
 106
 107   /*!
 108     Computes the bounding box of a given element. iP in numPol mode.
 109   */
 110   template<class MyMeshType, class MyMatrix>
 111   void CurveIntersector<MyMeshType,MyMatrix>::getElemBB (double*           bb,
 112                                                          const MyMeshType& mesh,
 113                                                          ConnType          iP,
 114                                                          ConnType          nb_nodes)
 115   {
 116     const double* coords = mesh.getCoordinatesPtr();
 117     const ConnType* conn_index = mesh.getConnectivityIndexPtr();
 118     const ConnType* conn = mesh.getConnectivityPtr();
 119     //initializing bounding box limits
 120     for(int idim=0; idim<SPACEDIM; idim++)
 121       {
 122         bb[2*idim  ] =  std::numeric_limits<double>::max();
 123         bb[2*idim+1] = -std::numeric_limits<double>::max();
 124       }
 125
 126     for (ConnType i=0; i<nb_nodes; i++)
 127       {
 128         //MN: iP= cell index, not node index, use of connectivity array ?
 129         const double* coord_node = coords +
 130           SPACEDIM*(OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[OTT<ConnType,numPol>::ind2C(iP)]+i)]));
 131         for(int idim=0; idim<SPACEDIM; idim++)
 132           {
 133             double x = *(coord_node+idim);
 134             bb[2*idim  ] = (x<bb[2*idim  ]) ? x : bb[2*idim  ];
 135             bb[2*idim+1] = (x>bb[2*idim+1]) ? x : bb[2*idim+1];
 136           }
 137       }
 138   }
 139
 140   /*!
 141    * \param [in] startOfSeg - input coming from intersectSegments or intersectSegmentsInternal
 142    * \param [in] endOfSeg - input coming from intersectSegments or intersectSegmentsInternal. NO Assume about sort
 143    * \param [in] pt - position of point that the method computes the bary coords for.
 144    */
 145   template<class MyMeshType, class MyMatrix>
 146   void CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
 147   {
 148     double deno(endOfSeg-startOfSeg);
 149     startPos = (endOfSeg-pt)/deno;
 150     startPos = std::max(startPos,0.); startPos = std::min(startPos,1.);
 151     endPos=1.-startPos;
 152   }
 153
 154   /*! Readjusts a set of bounding boxes so that they are extended
 155     in all dimensions for avoiding missing interesting intersections
 156
 157     \param bbox vector containing the bounding boxes
 158   */
 159   template<class MyMeshType, class MyMatrix>
 160   void CurveIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes (std::vector<double>& bbox,
 161                                                                    double adjustmentEpsAbs)
 162   {
 163     std::size_t size = bbox.size()/(2*SPACEDIM);
 164     for (std::size_t i=0; i<size; i++)
 165       {
 166         for(int idim=0; idim<SPACEDIM; idim++)
 167           {
 168             bbox[i*2*SPACEDIM+2*idim  ] -= adjustmentEpsAbs;
 169             bbox[i*2*SPACEDIM+2*idim+1] += adjustmentEpsAbs;
 170           }
 171       }
 172   }
 173
 174   /*!
 175    * @param icellT id in target mesh in format of MyMeshType.
 176    * @param coordsT output val that stores coordinates of the target cell
 177    * automatically resized to the right length.
 178    * @return true if segment is quadratic and in this case coordinates of medium node
 179    * are placed in the middle of coordsT
 180    */
 181   template<class MyMeshType, class MyMatrix>
 182   bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const
 183   {
 184     ConnType nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
 185     coordsT.resize(SPACEDIM*nbNodesT);
 186     for (ConnType iT=0; iT<nbNodesT; iT++)
 187       {
 188         for(int idim=0; idim<SPACEDIM; idim++)
 189           {
 190             coordsT[SPACEDIM*iT+idim] =
 191               _coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
 192           }
 193       }
 194     if ( nbNodesT > 2 )
 195       {
 196         for(int idim=0; idim<SPACEDIM; idim++)
 197           std::swap( coordsT[SPACEDIM*1+idim], coordsT[SPACEDIM*2+idim]);
 198         return true;
 199       }
 200     return false;
 201   }
 202
 203   template<class MyMeshType, class MyMatrix>
 204   typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const
 205   {
 206     ConnType nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
 207     if(nodeIdInCellT>=0 && nodeIdInCellT<nbNodesT)
 208       return OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+nodeIdInCellT)]);
 209     else
 210       throw Exception("getNodeIdOfTargetCellAt : error in nodeId in cell");
 211   }
 212
 213   /*!
 214    * @param icellS id in source mesh in format of MyMeshType.
 215    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
 216    * @return true if segment is quadratic and in this case coordinates of medium node
 217    * are placed in the middle of coordsS
 218    */
 219   template<class MyMeshType, class MyMatrix>
 220   bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS) const
 221   {
 222     ConnType nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
 223     coordsS.resize(SPACEDIM*nbNodesS);
 224     for(ConnType iS=0; iS<nbNodesS; iS++)
 225       {
 226         for(int idim=0; idim<SPACEDIM; idim++)
 227           {
 228             coordsS[SPACEDIM*iS+idim] =
 229               _coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
 230           }
 231       }
 232     if ( nbNodesS > 2 )
 233       {
 234         for(int idim=0; idim<SPACEDIM; idim++)
 235           std::swap( coordsS[SPACEDIM*1+idim], coordsS[SPACEDIM*2+idim]);
 236         return true;
 237       }
 238     return false;
 239   }
 240
 241   template<class MyMeshType, class MyMatrix>
 242   typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(ConnType icellS, ConnType nodeIdInCellS) const
 243   {
 244     ConnType nbNodesS(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]);
 245     if(nodeIdInCellS>=0 && nodeIdInCellS<nbNodesS)
 246       return OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+nodeIdInCellS)]);
 247     else
 248       throw Exception("getNodeIdOfSourceCellAt : error in nodeId in cell");
 249   }
 250
 251   /*!
 252    * \brief Return dual segments of given segment
 253    *  \param icell - given segment in C mode
 254    *  \param mesh - mesh
 255    *  \param segments - dual segments
 256    */
 257   template<class MyMeshType, class MyMatrix>
 258   void CurveIntersector<MyMeshType,MyMatrix>::getDualSegments(ConnType                   icell,
 259                                                               const MyMeshType&          mesh,
 260                                                               std::vector<TDualSegment>& segments)
 261   {
 262     // get coordinates of cell nodes
 263     ConnType nbNodes;
 264     std::vector<double>   ncoords;
 265     std::vector<ConnType> nodeIds;
 266     {
 267       const ConnType *connect   = mesh.getConnectivityPtr();
 268       const ConnType *connIndex = mesh.getConnectivityIndexPtr();
 269       const double *coords      = mesh.getCoordinatesPtr();
 270
 271       nbNodes = connIndex[icell+1] - connIndex[icell];
 272
 273       ncoords.resize(SPACEDIM*nbNodes);
 274       nodeIds.resize(nbNodes);
 275
 276       for(ConnType i=0; i<nbNodes; i++)
 277         for(int idim=0; idim<SPACEDIM; idim++)
 278           {
 279             nodeIds[i] = connect[OTT<ConnType,numPol>::conn2C(connIndex[OTT<ConnType,numPol>::ind2C(icell)]+i)];
 280             ncoords[SPACEDIM*i+idim] = coords[SPACEDIM*OTT<ConnType,numPol>::coo2C(nodeIds[i])+idim];
 281           }
 282       if ( nbNodes > 2 ) // quadratic segment, put medium node in the middle
 283         {
 284           for(int idim=0; idim<SPACEDIM; idim++)
 285             std::swap( ncoords[SPACEDIM*1+idim], ncoords[SPACEDIM*2+idim]);
 286           std::swap( nodeIds[1], nodeIds[2] );
 287         }
 288     }
 289
 290     // fill segments
 291     segments.clear();
 292     segments.reserve( 2*nbNodes );
 293     for(ConnType i=0; i<nbNodes-1; i++)
 294       {
 295         segments.push_back(TDualSegment());
 296         TDualSegment& seg1 = segments.back();
 297         segments.push_back(TDualSegment());
 298         TDualSegment& seg2 = segments.back();
 299
 300         seg1._nodeId = nodeIds[i];
 301         seg2._nodeId = nodeIds[i+1];
 302
 303         seg1._coords.resize( SPACEDIM * 2 );
 304         seg2._coords.resize( SPACEDIM * 2 );
 305
 306         for(int idim=0; idim<SPACEDIM; idim++)
 307           {
 308             double c1 = ncoords[SPACEDIM*i+idim];
 309             double c2 = ncoords[SPACEDIM*(i+1)+idim];
 310             double m = 0.5 * ( c1 + c2 );
 311             seg1._coords[ idim ] = c1;
 312             seg1._coords[ SPACEDIM + idim ] = m;
 313             seg2._coords[ idim ] = m;
 314             seg2._coords[ SPACEDIM + idim ] = c2;
 315           }
 316       }
 317   }
 318
 319   template<class MyMeshType, class MyMatrix>
 320   bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const
 321   {
 322     enum { X=0, Y };
 323     switch(SPACEDIM)
 324       {
 325       case 1:
 326         {
 327           xt  = coordsT[0];
 328           xs0 = coordsS[0]; xs1 = coordsS[1];
 329           return true;
 330         }
 331       case 2:
 332         {
 333           const double *s0(coordsS),*s1(coordsS + 2);
 334           double s01[2] = { s1[X]-s0[X], s1[Y]-s0[Y] }; // src segment direction
 335           double sSize = sqrt( s01[X]*s01[X] + s01[Y]*s01[Y] ); // src segment size
 336           if( sSize < this->_precision )
 337             return false;
 338           s01[X] /= sSize; s01[Y] /= sSize; // normalize s01
 339           double t[2] = { coordsT[X]-s0[X], coordsT[Y]-s0[Y] };
 340           xs0 = 0. ; xs1 = sSize; xt = s01[X]*t[X] + s01[Y]*t[Y];
 341           double proj_t_on_s[2] = { s0[X]+xt*s01[X], s0[Y]+xt*s01[Y] };
 342           double dist_t_s_vect[2] = { coordsT[X]-proj_t_on_s[X], coordsT[Y]-proj_t_on_s[Y] };
 343           double dist_t_s = sqrt( dist_t_s_vect[X]*dist_t_s_vect[X]+dist_t_s_vect[Y]*dist_t_s_vect[Y] );
 344           return dist_t_s < this->_precision;
 345         }
 346       default:
 347         throw Exception("CurveIntersector::projectionThis : space dimension of mesh must be 1 or 2");
 348       }
 349   }
 350
 351   template<class MyMeshType, class MyMatrix>
 352   bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS,
 353                                                              double& xs0, double& xs1, double& xt0, double& xt1) const
 354   {
 355     xt0 = coordsT[0]; xt1 = coordsT[1];
 356     xs0 = coordsS[0]; xs1 = coordsS[1];
 357     if ( SPACEDIM == 2 )
 358       {
 359         // Pass 2D->1D
 360
 361         enum { X=0, Y };
 362
 363         // check if two segments overlap in 2D within tolerance
 364
 365         const double* t0 = coordsT;
 366         const double* t1 = coordsT + 2;
 367         double t01[2] = { t1[X]-t0[X], t1[Y]-t0[Y] }; // tgt segment direction
 368         double tSize = sqrt( t01[X]*t01[X] + t01[Y]*t01[Y] ); // tgt segment size
 369         if ( tSize < _precision )
 370           return false; // degenerated segment
 371         t01[X] /= tSize, t01[Y] /= tSize; // normalize t01
 372
 373         const double* s0 = coordsS;
 374         const double* s1 = coordsS + 2;
 375         double t0s0[2] = { s0[X]-t0[X], s0[Y]-t0[Y] };
 376         double t0s1[2] = { s1[X]-t0[X], s1[Y]-t0[Y] };
 377         double nt01_x_t0s0 = t0s0[X] * t01[Y] - t0s0[Y] * t01[X]; // t0s0 dot norm of t01
 378         double nt01_x_t0s1 = t0s1[X] * t01[Y] - t0s1[Y] * t01[X]; // t0s1 dot norm of t01
 379         double dist_ts0 = fabs( nt01_x_t0s0 ); // dist from tgt seg to s0
 380         double dist_ts1 = fabs( nt01_x_t0s1 ); // dist from tgt seg to s1
 381         bool s0_out_of_tol = ( dist_ts0 > _tolerance );
 382         bool s1_out_of_tol = ( dist_ts1 > _tolerance );
 383         if ( nt01_x_t0s0 * nt01_x_t0s1 > 0 && ( s0_out_of_tol || s1_out_of_tol ))
 384           return false; // tgt segment is to far from src segment
 385
 386         double S0[2] = { s0[X], s0[Y] };
 387         double S1[2] = { s1[X], s1[Y] };
 388         if ( s0_out_of_tol ) // put s0 within tolerance
 389           {
 390             double t = _tolerance * nt01_x_t0s0 / dist_ts0; // signed tolerance
 391             double r = ( nt01_x_t0s0 - t ) / ( nt01_x_t0s0 - nt01_x_t0s1 );
 392             S0[X] = s0[X] * ( 1.-r ) + s1[X] * r;
 393             S0[Y] = s0[Y] * ( 1.-r ) + s1[Y] * r;
 394           }
 395         if ( s1_out_of_tol ) // put s1 within tolerance
 396           {
 397             double t = _tolerance * nt01_x_t0s1 / dist_ts1; // signed tolerance
 398             double r = ( nt01_x_t0s1 - t ) / ( nt01_x_t0s1 - nt01_x_t0s0 );
 399             S1[X] = s1[X] * ( 1.-r ) + s0[X] * r;
 400             S1[Y] = s1[Y] * ( 1.-r ) + s0[Y] * r;
 401           }
 402
 403         // project tgt and src segments to median line
 404
 405         double s01[2] = { S1[X]-S0[X], S1[Y]-S0[Y] }; // src segment direction
 406         double sSize = sqrt( s01[X]*s01[X] + s01[Y]*s01[Y] ); // src segment size
 407         if ( sSize < _precision )
 408           return false; // degenerated segment
 409         s01[X] /= sSize, s01[Y] /= sSize; // normalize s01
 410
 411         // make t01 and s01 codirected
 412         double t01_x_s01 = t01[X] * s01[X] + t01[Y] * s01[Y]; // t01 dot s01
 413         if ( t01_x_s01 < 0 )
 414           s01[X] = -s01[X], s01[Y] = -s01[Y];
 415
 416         double medianDir[2] = {
 417           t01[X] * ( 1.-_median_line) + s01[X] * _median_line,
 418           t01[Y] * ( 1.-_median_line) + s01[Y] * _median_line
 419         };
 420         double medianSize = sqrt( medianDir[X]*medianDir[X] + medianDir[Y]*medianDir[Y] );
 421         if ( medianSize < std::numeric_limits<double>::min() )
 422           return false; // strange...
 423         medianDir[X] /= medianSize, medianDir[Y] /= medianSize;
 424
 425         xt0 = t0[X] * medianDir[X] + t0[Y] * medianDir[Y];
 426         xt1 = t1[X] * medianDir[X] + t1[Y] * medianDir[Y];
 427         xs0 = S0[X] * medianDir[X] + S0[Y] * medianDir[Y];
 428         xs1 = S1[X] * medianDir[X] + S1[Y] * medianDir[Y];
 429
 430       } // if ( SPACEDIM == 2 )
 431     return true;
 432   }
 433
 434   /*!
 435    * \brief Return length of intersection of two segments
 436    */
 437   template<class MyMeshType, class MyMatrix>
 438   double CurveIntersector<MyMeshType,MyMatrix>::intersectSegmentsInternal(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const
 439   {
 440     if(!projectionThis(coordsT,coordsS,xs0,xs1,xt0,xt1))
 441       return 0.;
 442
 443     if ( xt0 > xt1 ) std::swap( xt0, xt1 );
 444     if ( xs0 > xs1 ) std::swap( xs0, xs1 );
 445
 446     double x0 = std::max( xt0, xs0 );
 447     double x1 = std::min( xt1, xs1 );
 448     return ( x0 < x1 ) ? ( x1 - x0 ) : 0.;
 449   }
 450
 451   template<class MyMeshType>
 452   class DummyMyMeshType1D
 453   {
 454   public:
 455     static const int MY_SPACEDIM=1;
 456     static const int MY_MESHDIM=8;
 457     typedef mcIdType MyConnType;
 458     static const INTERP_KERNEL::NumberingPolicy My_numPol=MyMeshType::My_numPol;
 459     // begin
 460     // useless, but for windows compilation ...
 461     const double *getCoordinatesPtr() const { return nullptr; }
 462     const MyConnType *getConnectivityPtr() const { return nullptr; }
 463     const MyConnType *getConnectivityIndexPtr() const { return nullptr; }
 464     INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType) const { return (INTERP_KERNEL::NormalizedCellType)0; }
 465     // end
 466   };
 467
 468   /*!
 469    * This method determines if a target point ( \a coordsT ) is in source seg2 contained in \a coordsS. To do so _precision attribute is used.
 470    * If target point is in, \a xs0, \a xs1 and \a xt are set to 1D referential for a further barycentric computation.
 471    * This method deals with SPACEDIM == 2 (see projectionThis).
 472    */
 473   template<class MyMeshType, class MyMatrix>
 474   bool CurveIntersector<MyMeshType,MyMatrix>::isPtIncludedInSeg(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const
 475   {
 476     if(!projectionThis(coordsT,coordsS,xs0,xs1,xt))
 477       return false;
 478     constexpr ConnType TAB[2]={0,1};
 479     const double coordsS_1D[2]={xs0,xs1};
 480     const double *coordsT_1D(&xt);
 481     return PointLocatorAlgos<DummyMyMeshType1D<MyMeshType>>::isElementContainsPoint(coordsT_1D,NORM_SEG2,coordsS_1D,TAB,2,this->_precision);
 482   }
 483
 484   /*!
 485    * \brief Return length of intersection of two segments
 486    */
 487   template<class MyMeshType, class MyMatrix>
 488   double CurveIntersector<MyMeshType,MyMatrix>::intersectSegments(const double *coordsT, const double *coordsS) const
 489   {
 490     double xs0,xs1,xt0,xt1;
 491     return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1);
 492   }
 493 }