src/INTERP_KERNEL/CurveIntersector.txx

   1 // Copyright (C) 2007-2014  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   /*!
 108     Computes the bouding box of a given element. iP in numPol mode.
 109   */
 110   //================================================================================
 111
 112   template<class MyMeshType, class MyMatrix>
 113   void CurveIntersector<MyMeshType,MyMatrix>::getElemBB (double*           bb,
 114                                                          const MyMeshType& mesh,
 115                                                          ConnType          iP,
 116                                                          ConnType          nb_nodes)
 117   {
 118     const double* coords = mesh.getCoordinatesPtr();
 119     const ConnType* conn_index = mesh.getConnectivityIndexPtr();
 120     const ConnType* conn = mesh.getConnectivityPtr();
 121     //initializing bounding box limits
 122     for(int idim=0; idim<SPACEDIM; idim++)
 123       {
 124         bb[2*idim  ] =  std::numeric_limits<double>::max();
 125         bb[2*idim+1] = -std::numeric_limits<double>::max();
 126       }
 127
 128     for (ConnType i=0; i<nb_nodes; i++)
 129       {
 130         //MN: iP= cell index, not node index, use of connectivity array ?
 131         const double* coord_node = coords +
 132           SPACEDIM*(OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[OTT<ConnType,numPol>::ind2C(iP)]+i)]));
 133         for(int idim=0; idim<SPACEDIM; idim++)
 134           {
 135             double x = *(coord_node+idim);
 136             bb[2*idim  ] = (x<bb[2*idim  ]) ? x : bb[2*idim  ];
 137             bb[2*idim+1] = (x>bb[2*idim+1]) ? x : bb[2*idim+1];
 138           }
 139       }
 140   }
 141
 142   //================================================================================
 143   /*! Readjusts a set of bounding boxes so that they are extended
 144     in all dimensions for avoiding missing interesting intersections
 145
 146     \param bbox vector containing the bounding boxes
 147   */
 148   //================================================================================
 149
 150   template<class MyMeshType, class MyMatrix>
 151   void CurveIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes (std::vector<double>& bbox,
 152                                                                    double adjustmentEpsAbs)
 153   {
 154     long size = bbox.size()/(2*SPACEDIM);
 155     for (int i=0; i<size; i++)
 156       {
 157         for(int idim=0; idim<SPACEDIM; idim++)
 158           {
 159             bbox[i*2*SPACEDIM+2*idim  ] -= adjustmentEpsAbs;
 160             bbox[i*2*SPACEDIM+2*idim+1] += adjustmentEpsAbs;
 161           }
 162       }
 163   }
 164
 165   //================================================================================
 166   /*!
 167    * @param icellT id in target mesh in format of MyMeshType.
 168    * @param coordsT output val that stores coordinates of the target cell
 169    * automatically resized to the right length.
 170    * @return true if segment is quadratic and in this case coordinates of medium node
 171    * are placed in the middle of coordsT
 172    */
 173   //================================================================================
 174
 175   template<class MyMeshType, class MyMatrix>
 176   bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates
 177   (ConnType icellT, std::vector<double>& coordsT)
 178   {
 179     int nbNodesT = _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] -
 180       _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
 181     coordsT.resize(SPACEDIM*nbNodesT);
 182     for (ConnType iT=0; iT<nbNodesT; iT++)
 183       {
 184         for(int idim=0; idim<SPACEDIM; idim++)
 185           {
 186             coordsT[SPACEDIM*iT+idim] =
 187               _coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
 188           }
 189       }
 190     if ( nbNodesT > 2 )
 191       {
 192         for(int idim=0; idim<SPACEDIM; idim++)
 193           std::swap( coordsT[SPACEDIM*1+idim], coordsT[SPACEDIM*2+idim]);
 194         return true;
 195       }
 196     return false;
 197   }
 198
 199   //================================================================================
 200   /*!
 201    * @param icellS id in source mesh in format of MyMeshType.
 202    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
 203    * @return true if segment is quadratic and in this case coordinates of medium node
 204    * are placed in the middle of coordsS
 205    */
 206   //================================================================================
 207
 208   template<class MyMeshType, class MyMatrix>
 209   bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates
 210   (ConnType icellS, std::vector<double>& coordsS)
 211   {
 212     int nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] -
 213       _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
 214     coordsS.resize(SPACEDIM*nbNodesS);
 215     for(ConnType iS=0; iS<nbNodesS; iS++)
 216       {
 217         for(int idim=0; idim<SPACEDIM; idim++)
 218           {
 219             coordsS[SPACEDIM*iS+idim] =
 220               _coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
 221           }
 222       }
 223     if ( nbNodesS > 2 )
 224       {
 225         for(int idim=0; idim<SPACEDIM; idim++)
 226           std::swap( coordsS[SPACEDIM*1+idim], coordsS[SPACEDIM*2+idim]);
 227         return true;
 228       }
 229     return false;
 230   }
 231
 232   //================================================================================
 233   /*!
 234    * \brief Return dual segments of given segment
 235    *  \param icell - given segment in C mode
 236    *  \param mesh - mesh
 237    *  \param segments - dual segments
 238    */
 239   //================================================================================
 240
 241   template<class MyMeshType, class MyMatrix>
 242   void CurveIntersector<MyMeshType,MyMatrix>::getDualSegments(ConnType                   icell,
 243                                                               const MyMeshType&          mesh,
 244                                                               std::vector<TDualSegment>& segments)
 245   {
 246     // get coordinates of cell nodes
 247     int nbNodes;
 248     std::vector<double> ncoords;
 249     std::vector<int>    nodeIds;
 250     {
 251       const ConnType *connect   = mesh.getConnectivityPtr();
 252       const ConnType *connIndex = mesh.getConnectivityIndexPtr();
 253       const double *coords      = mesh.getCoordinatesPtr();
 254
 255       nbNodes = connIndex[icell+1] - connIndex[icell];
 256
 257       ncoords.resize(SPACEDIM*nbNodes);
 258       nodeIds.resize(nbNodes);
 259
 260       for(int i=0; i<nbNodes; i++)
 261         for(int idim=0; idim<SPACEDIM; idim++)
 262           {
 263             nodeIds[i] = connect[OTT<ConnType,numPol>::conn2C(connIndex[OTT<ConnType,numPol>::ind2C(icell)]+i)];
 264             ncoords[SPACEDIM*i+idim] = coords[SPACEDIM*OTT<ConnType,numPol>::coo2C(nodeIds[i])+idim];
 265           }
 266       if ( nbNodes > 2 ) // quadratic segment, put medium node in the middle
 267         {
 268           for(int idim=0; idim<SPACEDIM; idim++)
 269             std::swap( ncoords[SPACEDIM*1+idim], ncoords[SPACEDIM*2+idim]);
 270           std::swap( nodeIds[1], nodeIds[2] );
 271         }
 272     }
 273
 274     // fill segments
 275     segments.clear();
 276     segments.reserve( 2*nbNodes );
 277     for(int i=0; i<nbNodes-1; i++)
 278       {
 279         segments.push_back(TDualSegment());
 280         TDualSegment& seg1 = segments.back();
 281         segments.push_back(TDualSegment());
 282         TDualSegment& seg2 = segments.back();
 283
 284         seg1._nodeId = nodeIds[i];
 285         seg2._nodeId = nodeIds[i+1];
 286
 287         seg1._coords.resize( SPACEDIM * 2 );
 288         seg2._coords.resize( SPACEDIM * 2 );
 289
 290         for(int idim=0; idim<SPACEDIM; idim++)
 291           {
 292             double c1 = ncoords[SPACEDIM*i+idim];
 293             double c2 = ncoords[SPACEDIM*(i+1)+idim];
 294             double m = 0.5 * ( c1 + c2 );
 295             seg1._coords[ idim ] = c1;
 296             seg1._coords[ SPACEDIM + idim ] = m;
 297             seg2._coords[ idim ] = m;
 298             seg2._coords[ SPACEDIM + idim ] = c2;
 299           }
 300       }
 301   }
 302
 303   //================================================================================
 304   /*!
 305    * \brief Return length of intersection of two segments
 306    */
 307   //================================================================================
 308
 309   template<class MyMeshType, class MyMatrix>
 310   double CurveIntersector<MyMeshType,MyMatrix>::intersectSegments(double *Coords_T,
 311                                                                   double *Coords_S)
 312   {
 313     double xt0 = Coords_T[0], xt1 = Coords_T[1];
 314     double xs0 = Coords_S[0], xs1 = Coords_S[1];
 315     if ( SPACEDIM == 2 )
 316       {
 317         // Pass 2D->1D
 318
 319         enum { X=0, Y };
 320
 321         // check if two segments overlap in 2D within tolerance
 322
 323         double* t0 = Coords_T;
 324         double* t1 = Coords_T + 2;
 325         double t01[2] = { t1[X]-t0[X], t1[Y]-t0[Y] }; // tgt segment direction
 326         double tSize = sqrt( t01[X]*t01[X] + t01[Y]*t01[Y] ); // tgt segment size
 327         if ( tSize < _precision ) return 0; // degenerated segment
 328         t01[X] /= tSize, t01[Y] /= tSize; // normalize t01
 329
 330         double* s0 = Coords_S;
 331         double* s1 = Coords_S + 2;
 332         double t0s0[2] = { s0[X]-t0[X], s0[Y]-t0[Y] };
 333         double t0s1[2] = { s1[X]-t0[X], s1[Y]-t0[Y] };
 334         double nt01_x_t0s0 = t0s0[X] * t01[Y] - t0s0[Y] * t01[X]; // t0s0 dot norm of t01
 335         double nt01_x_t0s1 = t0s1[X] * t01[Y] - t0s1[Y] * t01[X]; // t0s1 dot norm of t01
 336         double dist_ts0 = fabs( nt01_x_t0s0 ); // dist from tgt seg to s0
 337         double dist_ts1 = fabs( nt01_x_t0s1 ); // dist from tgt seg to s1
 338         bool s0_out_of_tol = ( dist_ts0 > _tolerance );
 339         bool s1_out_of_tol = ( dist_ts1 > _tolerance );
 340         if ( nt01_x_t0s0 * nt01_x_t0s1 > 0 && ( s0_out_of_tol || s1_out_of_tol ))
 341           return 0; // tgt segment is to far from src segment
 342
 343         double S0[2] = { s0[X], s0[Y] };
 344         double S1[2] = { s1[X], s1[Y] };
 345         if ( s0_out_of_tol ) // put s0 within tolerance
 346           {
 347             double t = _tolerance * nt01_x_t0s0 / dist_ts0; // signed tolerance
 348             double r = ( nt01_x_t0s0 - t ) / ( nt01_x_t0s0 - nt01_x_t0s1 );
 349             S0[X] = s0[X] * ( 1.-r ) + s1[X] * r;
 350             S0[Y] = s0[Y] * ( 1.-r ) + s1[Y] * r;
 351           }
 352         if ( s1_out_of_tol ) // put s1 within tolerance
 353           {
 354             double t = _tolerance * nt01_x_t0s1 / dist_ts1; // signed tolerance
 355             double r = ( nt01_x_t0s1 - t ) / ( nt01_x_t0s1 - nt01_x_t0s0 );
 356             S1[X] = s1[X] * ( 1.-r ) + s0[X] * r;
 357             S1[Y] = s1[Y] * ( 1.-r ) + s0[Y] * r;
 358           }
 359
 360         // project tgt and src segments to median line
 361
 362         double s01[2] = { S1[X]-S0[X], S1[Y]-S0[Y] }; // src segment direction
 363         double sSize = sqrt( s01[X]*s01[X] + s01[Y]*s01[Y] ); // src segment size
 364         if ( sSize < _precision ) return 0; // degenerated segment
 365         s01[X] /= sSize, s01[Y] /= sSize; // normalize s01
 366
 367         // make t01 and s01 codirected
 368         double t01_x_s01 = t01[X] * s01[X] + t01[Y] * s01[Y]; // t01 dot s01
 369         if ( t01_x_s01 < 0 )
 370           s01[X] = -s01[X], s01[Y] = -s01[Y];
 371
 372         double medianDir[2] = {
 373           t01[X] * ( 1.-_median_line) + s01[X] * _median_line,
 374           t01[Y] * ( 1.-_median_line) + s01[Y] * _median_line
 375         };
 376         double medianSize = sqrt( medianDir[X]*medianDir[X] + medianDir[Y]*medianDir[Y] );
 377         if ( medianSize < std::numeric_limits<double>::min() )
 378           return 0; // strange...
 379         medianDir[X] /= medianSize, medianDir[Y] /= medianSize;
 380
 381         xt0 = t0[X] * medianDir[X] + t0[Y] * medianDir[Y];
 382         xt1 = t1[X] * medianDir[X] + t1[Y] * medianDir[Y];
 383         xs0 = S0[X] * medianDir[X] + S0[Y] * medianDir[Y];
 384         xs1 = S1[X] * medianDir[X] + S1[Y] * medianDir[Y];
 385
 386       } // if ( SPACEDIM == 2 )
 387
 388     if ( xt0 > xt1 ) std::swap( xt0, xt1 );
 389     if ( xs0 > xs1 ) std::swap( xs0, xs1 );
 390
 391     double x0 = std::max( xt0, xs0 );
 392     double x1 = std::min( xt1, xs1 );
 393     return ( x0 < x1 ) ? ( x1 - x0 ) : 0.;
 394   }
 395
 396 }
 397
 398 #endif