Salome HOME
Some useful helpers into medcoupling python module
[tools/medcoupling.git] / src / INTERP_KERNEL / PlanarIntersector.txx
index 6f0c5aa5d0f7d8363446a06d42154dfe066f098d..5824d6923dfa381650dc35d39ef22abf60e3de0c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2014  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2016  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
@@ -32,7 +32,7 @@ namespace INTERP_KERNEL
   template<class MyMeshType, class MyMatrix>
   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),_max_distance_3Dsurf_intersect(md3DSurf),_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();
@@ -93,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)
@@ -344,7 +344,13 @@ namespace INTERP_KERNEL
     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++)
@@ -374,12 +380,12 @@ namespace INTERP_KERNEL
               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);
+            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++)
@@ -401,7 +407,7 @@ namespace INTERP_KERNEL
   }
   
   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