src/INTERP_KERNEL/CurveIntersector.txx

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