]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/INTERP_KERNEL/CurveIntersector.txx
Salome HOME
Copyright update 2021
[tools/medcoupling.git] / 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 }