X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FPlanarIntersector.txx;h=f5a9ae7298aa01722bdae885b7c7f12f2e123afe;hb=b219559763498c4bd10c730cd3d2c62b1eed45db;hp=df6593a0192990e70ce58f8845e188220dfddd71;hpb=48782c06022ca2caa36f849cb5a29ea4fe2aaa83;p=tools%2Fmedcoupling.git

diff --git a/src/INTERP_KERNEL/PlanarIntersector.txx b/src/INTERP_KERNEL/PlanarIntersector.txx
index df6593a01..f5a9ae729 100644
--- a/src/INTERP_KERNEL/PlanarIntersector.txx
+++ b/src/INTERP_KERNEL/PlanarIntersector.txx
@@ -1,21 +1,22 @@
-//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2019  CEA/DEN, EDF R&D
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+// Author : Anthony Geay (CEA/DEN)
 #ifndef __PLANARINTERSECTOR_TXX__
 #define __PLANARINTERSECTOR_TXX__
 
@@ -29,9 +30,9 @@
 namespace INTERP_KERNEL
 {
   template<class MyMeshType, class MyMatrix>
-  PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel):
+  PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double minDot3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel):
     _meshT(meshT),_meshS(meshS),
-    _dim_caracteristic(dimCaracteristic),_precision(precision),_median_plane(medianPlane),
+    _dim_caracteristic(dimCaracteristic),_max_distance_3Dsurf_intersect(md3DSurf),_min_dot_btw_3Dsurf_intersect(minDot3DSurf),_precision(precision),_median_plane(medianPlane),
     _do_rotate(doRotate),_orientation(orientation),_print_level(printLevel)
   {
     _connectT=meshT.getConnectivityPtr();
@@ -69,7 +70,7 @@ namespace INTERP_KERNEL
     int ibox=0;
     for(long icell=0; icell<nbelems; icell++)
       {
-        int nb_nodes_per_elem =conn_index[icell+1]-conn_index[icell];
+        ConnType nb_nodes_per_elem =conn_index[icell+1]-conn_index[icell];
         //initializing bounding box limits
         for(int idim=0; idim<SPACEDIM; idim++)
           {
@@ -77,7 +78,7 @@ namespace INTERP_KERNEL
             bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
           }
         //updating the bounding box with each node of the element
-        for (int j=0; j<nb_nodes_per_elem; j++)
+        for (ConnType j=0; j<nb_nodes_per_elem; j++)
           {
             const double* coord_node=coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[icell]+j)]);
             for(int idim=0; idim<SPACEDIM; idim++)
@@ -92,7 +93,7 @@ namespace INTERP_KERNEL
   }
 
   /*!
-    Computes the bouding box of a given element. iP in numPol mode.
+    Computes the bounding box of a given element. iP in numPol mode.
   */
   template<class MyMeshType, class MyMatrix>
   void PlanarIntersector<MyMeshType,MyMatrix>::getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes)
@@ -127,12 +128,12 @@ namespace INTERP_KERNEL
     \param bbox vector containing the bounding boxes
   */
   template<class MyMeshType, class MyMatrix>
-  void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double Surf3DAdjustmentEps)
+  void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs)
   {
     /* We build the segment tree for locating possible matching intersections*/
   
-    long size = bbox.size()/(2*SPACEDIM);
-    for (int i=0; i<size; i++)
+    std::size_t size = bbox.size()/(2*SPACEDIM);
+    for (std::size_t i=0; i<size; i++)
       {
         double max=- std::numeric_limits<double>::max();
         for(int idim=0; idim<SPACEDIM; idim++)
@@ -142,8 +143,8 @@ namespace INTERP_KERNEL
           }
         for(int idim=0; idim<SPACEDIM; idim++)
           {            
-            bbox[i*2*SPACEDIM+2*idim  ] -= Surf3DAdjustmentEps*max;
-            bbox[i*2*SPACEDIM+2*idim+1] += Surf3DAdjustmentEps*max;
+            bbox[i*2*SPACEDIM+2*idim  ] -= surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
+            bbox[i*2*SPACEDIM+2*idim+1] += surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
           }
       }
   }
@@ -155,7 +156,7 @@ namespace INTERP_KERNEL
   template<class MyMeshType, class MyMatrix>
   void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT)
   {
-    int nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+    ConnType nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
     coordsT.resize(SPACEDIM*nbNodesT);
     for (ConnType iT=0; iT<nbNodesT; iT++)
       for(int idim=0; idim<SPACEDIM; idim++)
@@ -169,12 +170,48 @@ namespace INTERP_KERNEL
   template<class MyMeshType, class MyMatrix>
   void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS)
   {
-    int nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+    ConnType nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
     coordsS.resize(SPACEDIM*nbNodesS);
     for (ConnType iS=0; iS<nbNodesS; iS++)
       for(int idim=0; idim<SPACEDIM; idim++)
         coordsS[SPACEDIM*iS+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
   }
+
+  /*!
+   * @param icellT id in target mesh in format of MyMeshType.
+   * @param offset is a value in C format that indicates the number of circular permutation.
+   * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
+   */
+  template<class MyMeshType, class MyMatrix>
+  void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinatesPermute(ConnType icellT, ConnType offset, std::vector<double>& coordsT)
+  {
+    ConnType nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+    coordsT.resize(SPACEDIM*nbNodesT);
+    for (ConnType iTTmp=0; iTTmp<nbNodesT; iTTmp++)
+      {
+        ConnType iT=(iTTmp+offset)%nbNodesT;
+        for(int idim=0; idim<SPACEDIM; idim++)
+          coordsT[SPACEDIM*iTTmp+idim]=_coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
+      }
+  }
+
+  /*!
+   * @param icellS id in source mesh in format of MyMeshType.
+   * @param offset is a value in C format that indicates the number of circular permutation.
+   * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
+   */
+  template<class MyMeshType, class MyMatrix>
+  void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinatesPermute(ConnType icellS, ConnType offset, std::vector<double>& coordsS)
+  {
+    ConnType nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+    coordsS.resize(SPACEDIM*nbNodesS);
+    for (ConnType iSTmp=0; iSTmp<nbNodesS; iSTmp++)
+      {
+        ConnType iS=(iSTmp+offset)%nbNodesS;
+        for(int idim=0; idim<SPACEDIM; idim++)
+          coordsS[SPACEDIM*iSTmp+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
+      }
+  }
   
   /*!
    * @param icellT id in target mesh in format of MyMeshType.
@@ -188,58 +225,6 @@ namespace INTERP_KERNEL
   template<class MyMeshType, class MyMatrix>
   void PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation)
   {
-    /*double epsilon=_precision*_dim_caracteristic;
-      coordsT.resize(SPACEDIM*nbNodesT);
-      coordsS.resize(SPACEDIM*nbNodesS);
-      int nb_dist_NodesT=nbNodesT;
-      int nb_dist_NodesS=nbNodesS;
-      int i_last = nbNodesT - 1;
-      const double * Pi_last=_coordsT +_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+i_last)];
-
-      for (int iT=0; iT<nbNodesT; iT++)
-      {
-      const double * PiT = _coordsT + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)]);
-      if(distance2<SPACEDIM>(Pi_last, PiT)>epsilon)
-      {
-      for (int idim=0; idim<SPACEDIM; idim++)
-      coordsT[SPACEDIM*iT+idim]=PiT[idim];
-      i_last=iT; Pi_last = PiT;
-      }
-      else 
-      nb_dist_NodesT--;
-      }
-      coordsT.resize(nb_dist_NodesT*SPACEDIM);
-      i_last = nbNodesS - 1;
-      Pi_last=_coordsS + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+i_last)]);
-      for (int iS=0; iS<nbNodesS; iS++)
-      {
-      const double * PiS=_coordsS+SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)]);
-      if(distance2<SPACEDIM>(Pi_last, PiS)>epsilon)
-      {
-      for (int idim=0; idim<SPACEDIM; idim++)
-      coordsS[SPACEDIM*iS+idim]=PiS[idim];
-      i_last=iS; Pi_last = PiS;
-      }
-      else
-      nb_dist_NodesS--;
-      }
-      coordsS.resize(nb_dist_NodesS*SPACEDIM);
-      //project cells T and S on the median plane
-      // and rotate the median plane
-      if(SPACEDIM==3) 
-      orientation = projectionThis(&coordsT[0], &coordsS[0], nb_dist_NodesT, nb_dist_NodesS);
-
-      //DEBUG PRINTS
-      if(_print_level >= 3) 
-      {
-      std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
-      std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " <<  nbNodesT << std::endl;
-      for(int iT =0; iT< nbNodesT; iT++)
-      {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
-      std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " <<  nbNodesS << std::endl;
-      for(int iS =0; iS< nbNodesS; iS++)
-      {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
-      }*/
     coordsT.resize(SPACEDIM*nbNodesT);
     coordsS.resize(SPACEDIM*nbNodesS);
     for (int idim=0; idim<SPACEDIM; idim++)
@@ -266,16 +251,36 @@ namespace INTERP_KERNEL
           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
       }
   }
+  
+  /*!
+   * Filtering out zero surfaces and badly oriented surfaces
+   * _orientation = -1,0,1,2
+   * -1 : the intersection is taken into account if target and cells have different orientation
+   * 0 : the intersection is always taken into account
+   * 1 : the intersection is taken into account if target and cells have the same orientation
+   * 2 : the absolute value of intersection is always taken into account 
+   */
+  template<class MyMeshType, class MyMatrix>
+  double PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(double val) const
+  {
+    if(_orientation==0)
+      return val;
+    if(_orientation==2)
+      return fabs(val);
+    if (( val > 0.0 && _orientation==1) || ( val < 0.0 && _orientation==-1 ))
+      return _orientation*val;
+    return 0.;
+  }
 
   template<class MyMeshType, class MyMatrix>
-  int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB)
+  int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, ConnType nb_NodesA, ConnType nb_NodesB)
   {
-    return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_median_plane,_do_rotate);
+    return Projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_max_distance_3Dsurf_intersect,_min_dot_btw_3Dsurf_intersect,_median_plane,_do_rotate);
   }
 
   template<class MyMeshType, class MyMatrix>
-  int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B, 
-                                                         int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate)
+  int PlanarIntersector<MyMeshType,MyMatrix>::Projection(double *Coords_A, double *Coords_B, 
+                                                         ConnType nb_NodesA, ConnType nb_NodesB, double epsilon, double md3DSurf, double minDot3DSurf, double median_plane, bool do_rotate)
   {
     double normal_A[3]={0,0,0};
     double normal_B[3]={0,0,0};
@@ -284,11 +289,12 @@ namespace INTERP_KERNEL
     bool same_orientation;
 
     //Find the normal to cells A and B
-    int i_A1=1;
-    while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++;
-    int i_A2=i_A1+1;
+    ConnType i_A1(1);
+    while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon)
+      i_A1++;
+    ConnType i_A2(i_A1+1);
     crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
-    double normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
+    double normA(sqrt(dotprod<SPACEDIM>(normal_A,normal_A)));
     while(i_A2<nb_NodesA && normA < epsilon)
       {
         crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
@@ -296,11 +302,12 @@ namespace INTERP_KERNEL
         normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
 
       }
-    int i_B1=1;
-    while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) i_B1++;
-    int i_B2=i_B1+1;
+    ConnType i_B1(1);
+    while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon)
+      i_B1++;
+    ConnType i_B2(i_B1+1);
     crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
-    double normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
+    double normB(sqrt(dotprod<SPACEDIM>(normal_B,normal_B)));
     while(i_B2<nb_NodesB && normB < epsilon)
       {
         crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
@@ -308,66 +315,99 @@ namespace INTERP_KERNEL
         normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
       }
 
+    //fabien option
+    if(md3DSurf>0.)
+      {
+        double coords_GA[3];
+        for(int i=0;i<3;i++)
+          {
+            coords_GA[i]=0.;
+            for (int j=0;j<nb_NodesA;j++)
+              coords_GA[i]+=Coords_A[3*j+i];
+            coords_GA[i]/=(double)nb_NodesA;
+          }
+        double G1[3],G2[3],G3[3];
+        for(int i=0;i<3;i++)
+          {
+            G1[i]=Coords_B[i]-coords_GA[i];
+            G2[i]=Coords_B[i+3]-coords_GA[i];
+            G3[i]=Coords_B[i+6]-coords_GA[i];
+          }
+        double prodvect[3];
+        prodvect[0]=G1[1]*G2[2]-G1[2]*G2[1];
+        prodvect[1]=G1[2]*G2[0]-G1[0]*G2[2];
+        prodvect[2]=G1[0]*G2[1]-G1[1]*G2[0];
+        double prodscal=prodvect[0]*G3[0]+prodvect[1]*G3[1]+prodvect[2]*G3[2];
+        if(fabs(prodscal)>md3DSurf)
+          return 0;
+      }
     if(i_A2<nb_NodesA && i_B2<nb_NodesB)
       {
         //Build the normal of the median plane
-        same_orientation = dotprod<SPACEDIM>(normal_A,normal_B)>=0;
+
+        double dotProd(dotprod<SPACEDIM>(normal_A,normal_B)/(normA*normB));
+
+        if(fabs(dotProd)<minDot3DSurf)
+          return 0;
+
+        same_orientation=(dotProd>=0);
         
         if(!same_orientation)
-          for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1;
+          for(int idim =0; idim< SPACEDIM; idim++)
+            normal_A[idim] *=-1;
         
-        double normB= sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
+        double normBB(sqrt(dotprod<SPACEDIM>(normal_B,normal_B)));
         
         for(int idim =0; idim< SPACEDIM; idim++)
-          linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normB;
+          linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normBB;
         double norm= sqrt(dotprod<SPACEDIM>(linear_comb,linear_comb));
 
         //Necessarily: norm>epsilon, no need to check
-        for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim]/=norm;
+        for(int idim =0; idim< SPACEDIM; idim++)
+          linear_comb[idim]/=norm;
         
         //Project the nodes of A and B on the median plane
-        for(int i_A=0; i_A<nb_NodesA; i_A++)
+        for(ConnType i_A=0; i_A<nb_NodesA; i_A++)
           {
             proj = dotprod<SPACEDIM>(&Coords_A[SPACEDIM*i_A],linear_comb);
             for(int idim =0; idim< SPACEDIM; idim++)
               Coords_A[SPACEDIM*i_A+idim] -=  proj*linear_comb[idim];
           }
-        for(int i_B=0; i_B<nb_NodesB; i_B++)
+        for(ConnType i_B=0; i_B<nb_NodesB; i_B++)
           {
             proj = dotprod<SPACEDIM>(Coords_B+SPACEDIM*i_B,linear_comb);
             for(int idim =0; idim< SPACEDIM; idim++)
               Coords_B[SPACEDIM*i_B+idim] -=  proj*linear_comb[idim];
           }
         
-        //Buid the matrix sending  A into the Oxy plane and apply it to A and B  
+        //Build the matrix sending  A into the Oxy plane and apply it to A and B  
         if(do_rotate)
           {
             TranslationRotationMatrix rotation;
             //rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
-            rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
-            for (int i=0; i<nb_NodesA; i++)    rotation.transform_vector(Coords_A+SPACEDIM*i);
-            for (int i=0; i<nb_NodesB; i++)    rotation.transform_vector(Coords_B+SPACEDIM*i);
+            Rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
+            for (int i=0; i<nb_NodesA; i++)
+              rotation.transform_vector(Coords_A+SPACEDIM*i);
+            for (int i=0; i<nb_NodesB; i++)
+              rotation.transform_vector(Coords_B+SPACEDIM*i);
           }
-        if(same_orientation)
-          return 1;
-        else return -1;
+        return same_orientation?1:-1;
       }
     else
       {
-        std::cout << " Maille dégénérée " << "epsilon = " << epsilon << std::endl;
+        std::cout << " Degenerated cell " << "epsilon = " << epsilon << std::endl;
         std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl;
         std::cout << " distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1])= " <<  distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1]) << std::endl;
         std::cout << "abs(normal_A) = " << fabs(normal_A[0]) << " ; " <<fabs( normal_A[1]) << " ; " << fabs(normal_A[2]) << std::endl;
         std::cout << " i_B1= " << i_B1 << " i_B2= " << i_B2 << std::endl; 
         std::cout << " distance2<SPACEDIM>(&Coords_B[0],&Coords_B[i_B1])= " <<  distance2<SPACEDIM>(Coords_B,Coords_B+i_B1) << std::endl;
         std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl;
-
         return 1;
       }
   }
   
   template<class MyMeshType, class MyMatrix>
-  void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
+  void PlanarIntersector<MyMeshType,MyMatrix>::Rotate3DTriangle(double* PP1, double*PP2, double*PP3,
                                                                 TranslationRotationMatrix& rotation_matrix)
   {
     //initializes