src/INTERP_KERNEL/CurveIntersector.txx

   1 // Copyright (C) 2007-2019  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 #ifndef __CURVEINTERSECTOR_TXX__
  21 #define __CURVEINTERSECTOR_TXX__
  22
  23 #include "CurveIntersector.hxx"
  24 #include "InterpolationUtils.hxx"
  25
  26 #include <limits>
  27
  28 namespace INTERP_KERNEL
  29 {
  30   template<class MyMeshType, class MyMatrix>
  31   CurveIntersector<MyMeshType,MyMatrix>
  32   ::CurveIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
  33                      double precision, double tolerance, double medianLine, int printLevel):
  34     _meshT(meshT),
  35     _meshS(meshS),
  36     _tolerance(tolerance),
  37     _precision(precision),
  38     _median_line(medianLine),
  39     _print_level(printLevel)
  40   {
  41     if ( SPACEDIM != 1 && SPACEDIM != 2 )
  42       throw Exception("CurveIntersector(): space dimension of mesh must be 1 or 2");
  43     if ( MESHDIM != 1 )
  44       throw Exception("CurveIntersector(): mesh dimension must be 1");
  45
  46     _connectT = meshT.getConnectivityPtr();
  47     _connectS = meshS.getConnectivityPtr();
  48     _connIndexT = meshT.getConnectivityIndexPtr();
  49     _connIndexS = meshS.getConnectivityIndexPtr();
  50     _coordsT = meshT.getCoordinatesPtr();
  51     _coordsS = meshS.getCoordinatesPtr();
  52   }
  53
  54   template<class MyMeshType, class MyMatrix>
  55   CurveIntersector<MyMeshType,MyMatrix>::~CurveIntersector()
  56   {
  57   }
  58
  59   //================================================================================
  60   /*!
  61     \brief creates the bounding boxes for all the cells of mesh \a mesh
  62
  63     \param mesh structure pointing to the mesh
  64     \param bbox vector containing the bounding boxes
  65   */
  66   //================================================================================
  67
  68   template<class MyMeshType, class MyMatrix>
  69   void CurveIntersector<MyMeshType,MyMatrix>::createBoundingBoxes (const MyMeshType&    mesh,
  70                                                                    std::vector<double>& bbox)
  71   {
  72     long nbelems = mesh.getNumberOfElements();
  73     bbox.resize(2*SPACEDIM* nbelems);
  74     const double* coords = mesh.getCoordinatesPtr();
  75     const ConnType* conn = mesh.getConnectivityPtr();
  76     const ConnType* conn_index = mesh.getConnectivityIndexPtr();
  77     int ibox=0;
  78     for(long icell=0; icell<nbelems; icell++)
  79       {
  80         int nb_nodes_per_elem = conn_index[icell+1]-conn_index[icell];
  81         //initializing bounding box limits
  82         for(int idim=0; idim<SPACEDIM; idim++)
  83           {
  84             bbox[2*SPACEDIM*ibox+2*idim]   =  std::numeric_limits<double>::max();
  85             bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
  86           }
  87         //updating the bounding box with each node of the element
  88         for (int j=0; j<nb_nodes_per_elem; j++)
  89           {
  90             const double* coord_node = coords +
  91               SPACEDIM*OTT<ConnType,numPol>
  92               ::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[icell]+j)]);
  93             for(int idim=0; idim<SPACEDIM; idim++)
  94               {
  95                 double x = *(coord_node+idim);
  96                 bbox[ibox*2*SPACEDIM + 2*idim]   =
  97                   ( bbox[ibox*2*SPACEDIM + 2*idim] < x ) ? bbox[ibox*2*SPACEDIM + 2*idim] : x;
  98                 bbox[ibox*2*SPACEDIM + 2*idim+1] =
  99                   ( bbox[ibox*2*SPACEDIM + 2*idim+1] > x ) ? bbox[ibox*2*SPACEDIM + 2*idim+1] : x;
 100               }
 101           }
 102         ibox++;
 103       }
 104   }
 105
 106   /*!
 107     Computes the bounding box of a given element. iP in numPol mode.
 108   */
 109   template<class MyMeshType, class MyMatrix>
 110   void CurveIntersector<MyMeshType,MyMatrix>::getElemBB (double*           bb,
 111                                                          const MyMeshType& mesh,
 112                                                          ConnType          iP,
 113                                                          ConnType          nb_nodes)
 114   {
 115     const double* coords = mesh.getCoordinatesPtr();
 116     const ConnType* conn_index = mesh.getConnectivityIndexPtr();
 117     const ConnType* conn = mesh.getConnectivityPtr();
 118     //initializing bounding box limits
 119     for(int idim=0; idim<SPACEDIM; idim++)
 120       {
 121         bb[2*idim  ] =  std::numeric_limits<double>::max();
 122         bb[2*idim+1] = -std::numeric_limits<double>::max();
 123       }
 124
 125     for (ConnType i=0; i<nb_nodes; i++)
 126       {
 127         //MN: iP= cell index, not node index, use of connectivity array ?
 128         const double* coord_node = coords +
 129           SPACEDIM*(OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[OTT<ConnType,numPol>::ind2C(iP)]+i)]));
 130         for(int idim=0; idim<SPACEDIM; idim++)
 131           {
 132             double x = *(coord_node+idim);
 133             bb[2*idim  ] = (x<bb[2*idim  ]) ? x : bb[2*idim  ];
 134             bb[2*idim+1] = (x>bb[2*idim+1]) ? x : bb[2*idim+1];
 135           }
 136       }
 137   }
 138
 139   /*!
 140    * \param [in] startOfSeg - input coming from intersectSegments or intersectSegmentsInternal
 141    * \param [in] endOfSeg - input coming from intersectSegments or intersectSegmentsInternal. Assume that endOfSeg>startOfSeg.
 142    * \param [in] pt - position of point that the method computes the bary coords for.
 143    */
 144   template<class MyMeshType, class MyMatrix>
 145   bool CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
 146   {
 147     double deno(endOfSeg-startOfSeg);
 148     startPos=(endOfSeg-pt)/deno;
 149     endPos=1.-startPos;
 150     return startPos>=0. && endPos>=0.;
 151   }
 152
 153   /*! Readjusts a set of bounding boxes so that they are extended
 154     in all dimensions for avoiding missing interesting intersections
 155
 156     \param bbox vector containing the bounding boxes
 157   */
 158   template<class MyMeshType, class MyMatrix>
 159   void CurveIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes (std::vector<double>& bbox,
 160                                                                    double adjustmentEpsAbs)
 161   {
 162     long size = bbox.size()/(2*SPACEDIM);
 163     for (int i=0; i<size; i++)
 164       {
 165         for(int idim=0; idim<SPACEDIM; idim++)
 166           {
 167             bbox[i*2*SPACEDIM+2*idim  ] -= adjustmentEpsAbs;
 168             bbox[i*2*SPACEDIM+2*idim+1] += adjustmentEpsAbs;
 169           }
 170       }
 171   }
 172
 173   /*!
 174    * @param icellT id in target mesh in format of MyMeshType.
 175    * @param coordsT output val that stores coordinates of the target cell
 176    * automatically resized to the right length.
 177    * @return true if segment is quadratic and in this case coordinates of medium node
 178    * are placed in the middle of coordsT
 179    */
 180   template<class MyMeshType, class MyMatrix>
 181   bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const
 182   {
 183     int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
 184     coordsT.resize(SPACEDIM*nbNodesT);
 185     for (ConnType iT=0; iT<nbNodesT; iT++)
 186       {
 187         for(int idim=0; idim<SPACEDIM; idim++)
 188           {
 189             coordsT[SPACEDIM*iT+idim] =
 190               _coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
 191           }
 192       }
 193     if ( nbNodesT > 2 )
 194       {
 195         for(int idim=0; idim<SPACEDIM; idim++)
 196           std::swap( coordsT[SPACEDIM*1+idim], coordsT[SPACEDIM*2+idim]);
 197         return true;
 198       }
 199     return false;
 200   }
 201
 202   template<class MyMeshType, class MyMatrix>
 203   typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const
 204   {
 205     int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
 206     if(nodeIdInCellT>=0 && nodeIdInCellT<nbNodesT)
 207       return OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+nodeIdInCellT)]);
 208     else
 209       throw Exception("getNodeIdOfTargetCellAt : error in nodeId in cell");
 210   }
 211
 212   /*!
 213    * @param icellS id in source mesh in format of MyMeshType.
 214    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
 215    * @return true if segment is quadratic and in this case coordinates of medium node
 216    * are placed in the middle of coordsS
 217    */
 218   template<class MyMeshType, class MyMatrix>
 219   bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS) const
 220   {
 221     int nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
 222     coordsS.resize(SPACEDIM*nbNodesS);
 223     for(ConnType iS=0; iS<nbNodesS; iS++)
 224       {
 225         for(int idim=0; idim<SPACEDIM; idim++)
 226           {
 227             coordsS[SPACEDIM*iS+idim] =
 228               _coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
 229           }
 230       }
 231     if ( nbNodesS > 2 )
 232       {
 233         for(int idim=0; idim<SPACEDIM; idim++)
 234           std::swap( coordsS[SPACEDIM*1+idim], coordsS[SPACEDIM*2+idim]);
 235         return true;
 236       }
 237     return false;
 238   }
 239
 240   template<class MyMeshType, class MyMatrix>
 241   typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(ConnType icellS, ConnType nodeIdInCellS) const
 242   {
 243     int nbNodesS(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]);
 244     if(nodeIdInCellS>=0 && nodeIdInCellS<nbNodesS)
 245       return OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+nodeIdInCellS)]);
 246     else
 247       throw Exception("getNodeIdOfSourceCellAt : error in nodeId in cell");
 248   }
 249
 250   /*!
 251    * \brief Return dual segments of given segment
 252    *  \param icell - given segment in C mode
 253    *  \param mesh - mesh
 254    *  \param segments - dual segments
 255    */
 256   template<class MyMeshType, class MyMatrix>
 257   void CurveIntersector<MyMeshType,MyMatrix>::getDualSegments(ConnType                   icell,
 258                                                               const MyMeshType&          mesh,
 259                                                               std::vector<TDualSegment>& segments)
 260   {
 261     // get coordinates of cell nodes
 262     int nbNodes;
 263     std::vector<double> ncoords;
 264     std::vector<int>    nodeIds;
 265     {
 266       const ConnType *connect   = mesh.getConnectivityPtr();
 267       const ConnType *connIndex = mesh.getConnectivityIndexPtr();
 268       const double *coords      = mesh.getCoordinatesPtr();
 269
 270       nbNodes = connIndex[icell+1] - connIndex[icell];
 271
 272       ncoords.resize(SPACEDIM*nbNodes);
 273       nodeIds.resize(nbNodes);
 274
 275       for(int i=0; i<nbNodes; i++)
 276         for(int idim=0; idim<SPACEDIM; idim++)
 277           {
 278             nodeIds[i] = connect[OTT<ConnType,numPol>::conn2C(connIndex[OTT<ConnType,numPol>::ind2C(icell)]+i)];
 279             ncoords[SPACEDIM*i+idim] = coords[SPACEDIM*OTT<ConnType,numPol>::coo2C(nodeIds[i])+idim];
 280           }
 281       if ( nbNodes > 2 ) // quadratic segment, put medium node in the middle
 282         {
 283           for(int idim=0; idim<SPACEDIM; idim++)
 284             std::swap( ncoords[SPACEDIM*1+idim], ncoords[SPACEDIM*2+idim]);
 285           std::swap( nodeIds[1], nodeIds[2] );
 286         }
 287     }
 288
 289     // fill segments
 290     segments.clear();
 291     segments.reserve( 2*nbNodes );
 292     for(int i=0; i<nbNodes-1; i++)
 293       {
 294         segments.push_back(TDualSegment());
 295         TDualSegment& seg1 = segments.back();
 296         segments.push_back(TDualSegment());
 297         TDualSegment& seg2 = segments.back();
 298
 299         seg1._nodeId = nodeIds[i];
 300         seg2._nodeId = nodeIds[i+1];
 301
 302         seg1._coords.resize( SPACEDIM * 2 );
 303         seg2._coords.resize( SPACEDIM * 2 );
 304
 305         for(int idim=0; idim<SPACEDIM; idim++)
 306           {
 307             double c1 = ncoords[SPACEDIM*i+idim];
 308             double c2 = ncoords[SPACEDIM*(i+1)+idim];
 309             double m = 0.5 * ( c1 + c2 );
 310             seg1._coords[ idim ] = c1;
 311             seg1._coords[ SPACEDIM + idim ] = m;
 312             seg2._coords[ idim ] = m;
 313             seg2._coords[ SPACEDIM + idim ] = c2;
 314           }
 315       }
 316   }
 317
 318   template<class MyMeshType, class MyMatrix>
 319   bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS,
 320                                                              double& xs0, double& xs1, double& xt0, double& xt1) const
 321   {
 322     xt0 = coordsT[0]; xt1 = coordsT[1];
 323     xs0 = coordsS[0]; xs1 = coordsS[1];
 324     if ( SPACEDIM == 2 )
 325       {
 326         // Pass 2D->1D
 327
 328         enum { X=0, Y };
 329
 330         // check if two segments overlap in 2D within tolerance
 331
 332         const double* t0 = coordsT;
 333         const double* t1 = coordsT + 2;
 334         double t01[2] = { t1[X]-t0[X], t1[Y]-t0[Y] }; // tgt segment direction
 335         double tSize = sqrt( t01[X]*t01[X] + t01[Y]*t01[Y] ); // tgt segment size
 336         if ( tSize < _precision )
 337           return false; // degenerated segment
 338         t01[X] /= tSize, t01[Y] /= tSize; // normalize t01
 339
 340         const double* s0 = coordsS;
 341         const double* s1 = coordsS + 2;
 342         double t0s0[2] = { s0[X]-t0[X], s0[Y]-t0[Y] };
 343         double t0s1[2] = { s1[X]-t0[X], s1[Y]-t0[Y] };
 344         double nt01_x_t0s0 = t0s0[X] * t01[Y] - t0s0[Y] * t01[X]; // t0s0 dot norm of t01
 345         double nt01_x_t0s1 = t0s1[X] * t01[Y] - t0s1[Y] * t01[X]; // t0s1 dot norm of t01
 346         double dist_ts0 = fabs( nt01_x_t0s0 ); // dist from tgt seg to s0
 347         double dist_ts1 = fabs( nt01_x_t0s1 ); // dist from tgt seg to s1
 348         bool s0_out_of_tol = ( dist_ts0 > _tolerance );
 349         bool s1_out_of_tol = ( dist_ts1 > _tolerance );
 350         if ( nt01_x_t0s0 * nt01_x_t0s1 > 0 && ( s0_out_of_tol || s1_out_of_tol ))
 351           return false; // tgt segment is to far from src segment
 352
 353         double S0[2] = { s0[X], s0[Y] };
 354         double S1[2] = { s1[X], s1[Y] };
 355         if ( s0_out_of_tol ) // put s0 within tolerance
 356           {
 357             double t = _tolerance * nt01_x_t0s0 / dist_ts0; // signed tolerance
 358             double r = ( nt01_x_t0s0 - t ) / ( nt01_x_t0s0 - nt01_x_t0s1 );
 359             S0[X] = s0[X] * ( 1.-r ) + s1[X] * r;
 360             S0[Y] = s0[Y] * ( 1.-r ) + s1[Y] * r;
 361           }
 362         if ( s1_out_of_tol ) // put s1 within tolerance
 363           {
 364             double t = _tolerance * nt01_x_t0s1 / dist_ts1; // signed tolerance
 365             double r = ( nt01_x_t0s1 - t ) / ( nt01_x_t0s1 - nt01_x_t0s0 );
 366             S1[X] = s1[X] * ( 1.-r ) + s0[X] * r;
 367             S1[Y] = s1[Y] * ( 1.-r ) + s0[Y] * r;
 368           }
 369
 370         // project tgt and src segments to median line
 371
 372         double s01[2] = { S1[X]-S0[X], S1[Y]-S0[Y] }; // src segment direction
 373         double sSize = sqrt( s01[X]*s01[X] + s01[Y]*s01[Y] ); // src segment size
 374         if ( sSize < _precision )
 375           return false; // degenerated segment
 376         s01[X] /= sSize, s01[Y] /= sSize; // normalize s01
 377
 378         // make t01 and s01 codirected
 379         double t01_x_s01 = t01[X] * s01[X] + t01[Y] * s01[Y]; // t01 dot s01
 380         if ( t01_x_s01 < 0 )
 381           s01[X] = -s01[X], s01[Y] = -s01[Y];
 382
 383         double medianDir[2] = {
 384           t01[X] * ( 1.-_median_line) + s01[X] * _median_line,
 385           t01[Y] * ( 1.-_median_line) + s01[Y] * _median_line
 386         };
 387         double medianSize = sqrt( medianDir[X]*medianDir[X] + medianDir[Y]*medianDir[Y] );
 388         if ( medianSize < std::numeric_limits<double>::min() )
 389           return false; // strange...
 390         medianDir[X] /= medianSize, medianDir[Y] /= medianSize;
 391
 392         xt0 = t0[X] * medianDir[X] + t0[Y] * medianDir[Y];
 393         xt1 = t1[X] * medianDir[X] + t1[Y] * medianDir[Y];
 394         xs0 = S0[X] * medianDir[X] + S0[Y] * medianDir[Y];
 395         xs1 = S1[X] * medianDir[X] + S1[Y] * medianDir[Y];
 396
 397       } // if ( SPACEDIM == 2 )
 398     return true;
 399   }
 400
 401   /*!
 402    * \brief Return length of intersection of two segments
 403    */
 404   template<class MyMeshType, class MyMatrix>
 405   double CurveIntersector<MyMeshType,MyMatrix>::intersectSegmentsInternal(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const
 406   {
 407     if(!projectionThis(coordsT,coordsS,xs0,xs1,xt0,xt1))
 408       return 0.;
 409
 410     if ( xt0 > xt1 ) std::swap( xt0, xt1 );
 411     if ( xs0 > xs1 ) std::swap( xs0, xs1 );
 412
 413     double x0 = std::max( xt0, xs0 );
 414     double x1 = std::min( xt1, xs1 );
 415     return ( x0 < x1 ) ? ( x1 - x0 ) : 0.;
 416   }
 417
 418   /*!
 419    * \brief Return length of intersection of two segments
 420    */
 421   template<class MyMeshType, class MyMatrix>
 422   double CurveIntersector<MyMeshType,MyMatrix>::intersectSegments(const double *coordsT, const double *coordsS) const
 423   {
 424     double xs0,xs1,xt0,xt1;
 425     return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1);
 426   }
 427
 428 }
 429
 430 #endif