]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/INTERP_KERNEL/CurveIntersector.txx
Salome HOME
Update copyrights 2014.
[tools/medcoupling.git] / 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